test

Author: ESeNonFossiIo

include/interfaces/ALE_navier_stokes.h

@@ -59,24 +59,30 @@ class ALENavierStokes
 
   void connect_to_signals() const
   {
-    // first of all we get the struct Signals from pidomus
     auto &signals = this->get_signals();
-    signals.end_run.connect(
-      [&]()
+    auto &pcout = this->get_pcout();
+
+    signals.end_cycle.connect([&]()
     {
       clean();
     });
   }
 
   void
-  clean() const;
+  clean() const
+  {
+    mapping = nullptr;
+  };
 
   const Mapping<dim,spacedim>
   &get_default_mapping() const;
 
   const Mapping<dim,spacedim>
-  &get_fe_mapping() const;
-  
+  &get_fe_mapping() const
+  {
+    return StaticMappingQ1<dim,spacedim>::mapping;
+  };
+
 private:
 
 // Physical parameter
@@ -87,7 +93,7 @@ class ALENavierStokes
   bool AMG_u_use_inverse_operator;
   bool AMG_d_use_inverse_operator;
   bool AMG_Ap_use_inverse_operator;
-  
+
   /**
    * AMG preconditioner for velocity-velocity matrix.
    */
@@ -102,7 +108,7 @@ class ALENavierStokes
    * AMG preconditioner for the pressure stifness matrix.
    */
   mutable ParsedAMGPreconditioner AMG_Ap;
-  
+
   /**
    * Jacobi preconditioner for the pressure mass matrix.
    */
@@ -143,21 +149,6 @@ set_matrix_couplings(std::vector<std::string> &couplings) const
   couplings[2] = "0,0,0; 0,0,0; 0,0,1"; // TODO: remove unused couplings
 }
 
-template <int dim, int spacedim, typename LAC>
-void ALENavierStokes<dim,spacedim,LAC>::
-clean() const
-{
-  mapping = nullptr;
-}
-
-template <int dim, int spacedim, typename LAC>
-const Mapping<dim,spacedim> &
-ALENavierStokes<dim,spacedim,LAC>::
-get_fe_mapping() const
-{
-  return StaticMappingQ1<dim,spacedim>::mapping;
-}
-
 template <int dim, int spacedim, typename LAC>
 const Mapping<dim,spacedim> &
 ALENavierStokes<dim,spacedim,LAC>::
@@ -169,6 +160,7 @@ get_default_mapping() const
 
   if (!mapping)
     {
+      std::cout << "update mapping" << std::endl << std::flush;
       const unsigned int degree = this->get_fe().degree;
       const DoFHandler<dim,spacedim> &dh = this->get_dof_handler();
       mapping = SP(new MappingQEulerian< dim, typename LAC::VectorType, spacedim >(degree, dh, solution));
@@ -209,10 +201,10 @@ declare_parameters (ParameterHandler &prm)
                       Patterns::Bool(),
                       "Enable the use of inverse operator for AMG u");
 
-                      this->add_parameter(prm, &AMG_Ap_use_inverse_operator,
-                                          "AMG Ap - use inverse operator", "false",
-                                          Patterns::Bool(),
-                                          "Enable the use of inverse operator for AMG Ap");
+  this->add_parameter(prm, &AMG_Ap_use_inverse_operator,
+                      "AMG Ap - use inverse operator", "false",
+                      Patterns::Bool(),
+                      "Enable the use of inverse operator for AMG Ap");
 }
 
 template <int dim, int spacedim, typename LAC>
@@ -256,7 +248,7 @@ energies_and_residuals(const typename DoFHandler<dim,spacedim>::active_cell_iter
 
 // pressure:
   auto &ps = fe_cache.get_values( "solution", "p", pressure, et);
-auto &grad_ps = fe_cache.get_gradients( "solution", "grad_p", pressure, et);
+  auto &grad_ps = fe_cache.get_gradients( "solution", "grad_p", pressure, et);
 
 // Jacobian:
   auto &JxW = fe_cache.get_JxW_values();
@@ -267,7 +259,7 @@ auto &grad_ps = fe_cache.get_gradients( "solution", "grad_p", pressure, et);
 
   ResidualType res_u = 0;
   ResidualType res_d = 0;
-  
+
   for (unsigned int quad=0; quad<n_quad_points; ++quad)
     {
 // variables:
@@ -290,7 +282,7 @@ auto &grad_ps = fe_cache.get_gradients( "solution", "grad_p", pressure, et);
 // pressure:
       const ResidualType &p = ps[quad];
       const Tensor<1, dim, ResidualType> &grad_p = grad_ps[quad];
-      
+
 // jacobian of ALE transformation:
       const double J_ale = SacadoTools::to_double(determinant(F));
 
@@ -316,26 +308,26 @@ auto &grad_ps = fe_cache.get_gradients( "solution", "grad_p", pressure, et);
 // pressure:
           auto p_test = fev[pressure].value(i,quad);
           auto grad_p_test = fev[pressure].gradient(i,quad);
-          
+
           // Preconditioner:
           residual[1][i] +=
             (
-              p * p_test //* J_ale
+              p * p_test * J_ale
             )*JxW[quad];
 
           residual[2][i] +=
             (
-              grad_p * grad_p_test //* J_ale
+              grad_p * grad_p_test * J_ale
             )*JxW[quad];
-              
+
           // Navier Stokes:
           res_u =
 // time derivative term
             rho*scalar_product( u_dot  , u_test )
 //
             + scalar_product( grad_u * ( F_inv * ( u_old - d_dot ) )  , u_test )
 //
-            + scalar_product(sigma , grad_u_test * F_inv )
+            + scalar_product(sigma * Ft_inv, grad_u_test )
 // divergence free constriant
             - div_u  * p_test;
 
@@ -369,6 +361,7 @@ ALENavierStokes<dim,spacedim,LAC>::compute_system_operators(
 
   AMG_d.initialize_preconditioner<dim, spacedim>( matrices[0]->block(0,0), fe, dh);
   AMG_u.initialize_preconditioner<dim, spacedim>( matrices[0]->block(1,1), fe, dh);
+  AMG_Ap.initialize_preconditioner<dim, spacedim>( matrices[2]->block(2,2), fe, dh);
   jac_M.initialize_preconditioner<>(matrices[1]->block(2,2));
 
 ////////////////////////////////////////////////////////////////////////////
@@ -432,20 +425,20 @@ ALENavierStokes<dim,spacedim,LAC>::compute_system_operators(
     }
 
   auto Ap = linear_operator< TrilinosWrappers::MPI::Vector >( matrices[2]->block(2,2) );
-    LinearOperator<BVEC> Ap_inv;
-    if (AMG_Ap_use_inverse_operator)
-      {
-        Ap_inv = inverse_operator(Ap,
-                                  solver_GMRES,
-                                  jac_M);
-      }
-    else
-      {
-        Ap_inv = linear_operator<BVEC>(matrices[2]->block(2,2),
-                                       AMG_Ap);
-      }
+  LinearOperator<BVEC> Ap_inv;
+  if (AMG_Ap_use_inverse_operator)
+    {
+      Ap_inv = inverse_operator(Ap,
+                                solver_GMRES,
+                                jac_M);
+    }
+  else
+    {
+      Ap_inv = linear_operator<BVEC>(matrices[2]->block(2,2),
+                                     AMG_Ap);
+    }
   auto alpha = this->get_alpha();
-  auto Schur_inv = nu * Mp_inv ;//+ (rho * alpha * Ap_inv);
+  auto Schur_inv = 1/nu * Mp_inv;// + rho * alpha * Ap_inv;
 
   if (AMG_d_use_inverse_operator)
     {

include/pidomus_signals.h

@@ -94,6 +94,7 @@ struct Signals
   boost::signals2::signal<void ()> begin_refine_and_transfer_solutions;
   boost::signals2::signal<void ()> begin_assemble_matrices;
   boost::signals2::signal<void ()> begin_solver_should_restart;
+  boost::signals2::signal<void ()> begin_cycle;
   boost::signals2::signal<void ()> begin_run;
 
   boost::signals2::signal<void ()> end_make_grid_fe;
@@ -105,6 +106,7 @@ struct Signals
   boost::signals2::signal<void ()> end_refine_and_transfer_solutions;
   boost::signals2::signal<void ()> end_assemble_matrices;
   boost::signals2::signal<void ()> end_solver_should_restart;
+  boost::signals2::signal<void ()> end_cycle;
   boost::signals2::signal<void ()> end_run;
 
 };

source/pidomus.cc

@@ -126,6 +126,7 @@ void piDoMUS<dim, spacedim, LAC>::run ()
 
   for (current_cycle = 0; current_cycle < n_cycles; ++current_cycle)
     {
+      signals.begin_cycle();
       if (current_cycle == 0)
         {
           make_grid_fe();
@@ -145,6 +146,7 @@ void piDoMUS<dim, spacedim, LAC>::run ()
           euler.start_ode(solution, solution_dot);
         }
       eh.error_from_exact(interface.get_error_mapping(), *dof_handler, locally_relevant_solution, exact_solution);
+      signals.end_cycle();
     }
 
   eh.output_table(pcout);