diff --git a/include/base_interface.h b/include/base_interface.h index e69908c0..01c795cc 100644 --- a/include/base_interface.h +++ b/include/base_interface.h @@ -247,6 +247,7 @@ class BaseInterface : public ParameterAcceptor, public SimulatorAccess &prec_op, LinearOperator &prec_op_finer) const; + /** * Compute linear operators needed by the problem. When using * deal.II vector and matrix types, this function is empty, since a diff --git a/include/interfaces/ALE_navier_stokes.h b/include/interfaces/ALE_navier_stokes.h index e4f53a28..9b41bdac 100644 --- a/include/interfaces/ALE_navier_stokes.h +++ b/include/interfaces/ALE_navier_stokes.h @@ -12,12 +12,15 @@ #include "pde_system_interface.h" +#include + #include #include #include #include #include +#include //////////////////////////////////////////////////////////////////////////////// /// ALE Navier Stokes interface: @@ -54,6 +57,32 @@ class ALENavierStokes void set_matrix_couplings(std::vector &couplings) const; + void connect_to_signals() const + { + auto &signals = this->get_signals(); + auto &pcout = this->get_pcout(); + + signals.end_cycle.connect([&]() + { + clean(); + }); + } + + void + clean() const + { + mapping = nullptr; + }; + + const Mapping + &get_default_mapping() const; + + const Mapping + &get_fe_mapping() const + { + return StaticMappingQ1::mapping; + }; + private: // Physical parameter @@ -62,22 +91,39 @@ class ALENavierStokes bool Mp_use_inverse_operator; bool AMG_u_use_inverse_operator; + bool AMG_v_use_inverse_operator; bool AMG_d_use_inverse_operator; + bool AMG_Ap_use_inverse_operator; /** * AMG preconditioner for velocity-velocity matrix. */ mutable ParsedAMGPreconditioner AMG_u; + /** + * AMG preconditioner for velocity-velocity matrix. + */ + mutable ParsedAMGPreconditioner AMG_v; + /** * AMG preconditioner for the pressure stifness matrix. */ mutable ParsedAMGPreconditioner AMG_d; + /** + * AMG preconditioner for the pressure stifness matrix. + */ + mutable ParsedAMGPreconditioner AMG_Ap; + /** * Jacobi preconditioner for the pressure mass matrix. */ mutable ParsedJacobiPreconditioner jac_M; + + /** + * Mapping.. + */ + mutable shared_ptr > mapping; }; @@ -87,18 +133,49 @@ ALENavierStokes() : PDESystemInterface, LAC>( "ALE Navier Stokes Interface", - dim+dim+1, - 2, - "FESystem[FE_Q(2)^d-FE_Q(2)^d-FE_Q(1)]", - "d,d,u,u,p", - "1,1,0"), + dim+dim+1+dim, + 3, + "FESystem[FE_Q(2)^d-FE_Q(2)^d-FE_Q(1)-FE_Q(2)^d]", + dim == 2 ? "d,d,u,u,p,v,v" : "d,d,d,u,u,u,p,v,v,v", + "1,1,0,1"), AMG_u("AMG for u"), + AMG_v("AMG for v"), AMG_d("AMG for d"), + AMG_Ap("AMG for Ap"), jac_M("Jacobi for M") { this->init(); } +template +void ALENavierStokes:: +set_matrix_couplings(std::vector &couplings) const +{ + couplings[0] = "1,1,1,1; 1,1,1,1; 1,1,1,1; 1,1,1,1"; // TODO: remove unused couplings + couplings[1] = "1,1,1,1; 1,1,1,1; 1,1,1,1; 1,1,1,1"; // TODO: remove unused couplings + couplings[2] = "1,1,1,1; 1,1,1,1; 1,1,1,1; 1,1,1,1"; // TODO: remove unused couplings +} + +template +const Mapping & +ALENavierStokes:: +get_default_mapping() const +{ + const typename LAC::VectorType &solution = this->get_locally_relevant_solution(); + if (solution.size() == 0) + return StaticMappingQ1::mapping; + + if (!mapping) + { + // std::cout << "update mapping" << std::endl << std::flush; + const unsigned int degree = this->get_fe().degree; + const DoFHandler &dh = this->get_dof_handler(); + mapping = SP(new MappingQEulerian< dim, typename LAC::VectorType, spacedim >(degree, dh, solution)); + } + + return *mapping; +} + template void ALENavierStokes:: declare_parameters (ParameterHandler &prm) @@ -130,6 +207,16 @@ declare_parameters (ParameterHandler &prm) "AMG u - use inverse operator", "false", Patterns::Bool(), "Enable the use of inverse operator for AMG u"); + + this->add_parameter(prm, &AMG_v_use_inverse_operator, + "AMG v - use inverse operator", "false", + Patterns::Bool(), + "Enable the use of inverse operator for AMG v"); + + 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 @@ -137,14 +224,6 @@ void ALENavierStokes:: parse_parameters_call_back () {} -template -void ALENavierStokes:: -set_matrix_couplings(std::vector &couplings) const -{ - couplings[0] = "1,1,1; 1,1,1; 1,1,0"; // TODO: Select only not null entries - couplings[1] = "0,0,0; 0,0,0; 0,0,1"; -} - template template void @@ -158,90 +237,34 @@ energies_and_residuals(const typename DoFHandler::active_cell_iter const FEValuesExtractors::Vector displacement(0); const FEValuesExtractors::Vector velocity(dim); const FEValuesExtractors::Scalar pressure(2*dim); + const FEValuesExtractors::Vector velocity_aux(2*dim+1); ResidualType et = 0; double dummy = 0.0; - double h = cell->diameter(); - - for (unsigned int face=0; face < GeometryInfo::faces_per_cell; ++face) - { - unsigned int face_id = cell->face(face)->boundary_id(); - if (cell->face(face)->at_boundary()) - { - auto nitsche = this->get_dirichlet_bcs(); - if (nitsche.acts_on_id(face_id)) - { - bool check = false; - for (unsigned int i = 0; ireinit(et, cell, face, fe_cache); -// Displacement: - auto &d_ = fe_cache.get_values("solution", "d", displacement, et); - auto &d_dot_ = fe_cache.get_values( "solution_dot", "d_dot", displacement, et); - -// Velocity: -// auto &grad_u_ = fe_cache.get_gradients("solution", "u", velocity, et); - auto &u_ = fe_cache.get_values("solution", "grad_u", velocity, et); - -// Pressure: -// auto &grad_p_ = fe_cache.get_gradients("solution", "p", pressure, et); - - auto &fev = fe_cache.get_current_fe_values(); - auto &q_points = fe_cache.get_quadrature_points(); - auto &JxW = fe_cache.get_JxW_values(); - - for (unsigned int q=0; q &d = d_[q]; - const Tensor<1, dim, ResidualType> &d_dot = d_dot_[q]; - -// Velocity: - const Tensor<1, dim, ResidualType> &u = u_[q]; - - for (unsigned int i=0; iat_boundary - } - } - } // end loop over faces - this->reinit (et, cell, fe_cache); // displacement: -// auto &ds = fe_cache.get_values( "solution", "d", displacement, et); auto &grad_ds = fe_cache.get_gradients( "solution", "grad_d", displacement, et); -// auto &div_ds = fe_cache.get_divergences( "solution", "div_d", displacement, et); auto &Fs = fe_cache.get_deformation_gradients( "solution", "Fd", displacement, et); auto &ds_dot = fe_cache.get_values( "solution_dot", "d_dot", displacement, et); + auto &vs = fe_cache.get_values( "solution", "v", velocity_aux, et); + auto &div_vs = fe_cache.get_divergences( "solution", "div_v", velocity_aux, et); + // velocity: auto &us = fe_cache.get_values( "solution", "u", velocity, et); auto &grad_us = fe_cache.get_gradients( "solution", "grad_u", velocity, et); auto &div_us = fe_cache.get_divergences( "solution", "div_u", velocity, et); - auto &sym_grad_us = fe_cache.get_symmetric_gradients( "solution", "u", velocity, et); + auto &sym_grad_us = fe_cache.get_symmetric_gradients( "solution", "sym_grad_u", velocity, et); auto &us_dot = fe_cache.get_values( "solution_dot", "u_dot", velocity, et); // Previous time step solution: - auto &u_olds = fe_cache.get_values("explicit_solution","u",velocity,dummy); - // auto &ds_dot_old = fe_cache.get_values("explicit_solution","d_dot",displacement,dummy); - + auto &u_olds = fe_cache.get_values("explicit_solution","uo",velocity,dummy); // pressure: auto &ps = fe_cache.get_values( "solution", "p", pressure, et); + auto &grad_ps = fe_cache.get_gradients( "solution", "grad_p", pressure, et); // Jacobian: auto &JxW = fe_cache.get_JxW_values(); @@ -250,43 +273,61 @@ energies_and_residuals(const typename DoFHandler::active_cell_iter auto &fev = fe_cache.get_current_fe_values(); + ResidualType res_u = 0; + ResidualType res_d = 0; + for (unsigned int quad=0; quad &v = vs[quad]; + const ResidualType &div_v= div_vs[quad]; + // velocity: const ResidualType &div_u = div_us[quad]; + const Tensor<1, dim, ResidualType> &u = us[quad]; const Tensor<1, dim, ResidualType> &u_dot = us_dot[quad]; const Tensor<2, dim, ResidualType> &grad_u = grad_us[quad]; const Tensor <2, dim, ResidualType> &sym_grad_u = sym_grad_us[quad]; + + // displacement const Tensor<1, dim, ResidualType> &d_dot = ds_dot[quad]; const Tensor<2, dim, ResidualType> &grad_d = grad_ds[quad]; const Tensor <2, dim, ResidualType> &F = Fs[quad]; - ResidualType J = determinant(F); - const Tensor <2, dim, ResidualType> &F_inv = invert(F); - const Tensor <2, dim, ResidualType> &Ft_inv = transpose(F_inv); + const Tensor <2, dim, ResidualType> F_inv = invert(F); + const Tensor <2, dim, ResidualType> Ft_inv = transpose(F_inv); + // const Tensor <2, dim, ResidualType> &F_sac = Fs[quad]; + // Tensor <2, dim, double> F; + // for (unsigned int i = 0; i< dim; ++i) + // for (unsigned int j = 0; j< dim; ++j) + // F[i][j] = SacadoTools::to_double(F_sac[i][j]); + // const Tensor <2, dim, double> F_inv = invert(F); + // const Tensor <2, dim, double> Ft_inv = transpose(F_inv); + + // jacobian of ALE transformation: + const ResidualType J_ale = determinant(F); + // const double J_ale = SacadoTools::to_double(determinant(F)); // Previous time step solution: const Tensor<1, dim, ResidualType> &u_old = u_olds[quad]; // pressure: const ResidualType &p = ps[quad]; + const Tensor<1, dim, ResidualType> &grad_p = grad_ps[quad]; -// others: - auto J_ale = J; // jacobian of ALE transformation -// auto div_u_ale = (J_ale * (F_inv * u) ); - Tensor <2, dim, ResidualType> Id; + Tensor <2, dim, ResidualType> p_x_Id; for (unsigned int i = 0; i sigma = - - Id + my_rho * ( nu* sym_grad_u * F_inv + ( Ft_inv * transpose(sym_grad_u) ) ) ; + - p_x_Id + my_rho * ( sym_grad_u * F_inv + ( Ft_inv * transpose(sym_grad_u) ) ) ; for (unsigned int i=0; i::active_cell_iter // pressure: auto p_test = fev[pressure].value(i,quad); -// auto q = fev[pressure].value(i,quad); + auto grad_p_test = fev[pressure].gradient(i,quad); + // Preconditioner: residual[1][i] += ( - (1./nu)*p*p_test + p * p_test * J_ale )*JxW[quad]; - residual[0][i] += + residual[2][i] += ( - // time derivative term - rho*scalar_product( u_dot * J_ale , u_test ) -// - + scalar_product( grad_u * ( F_inv * ( u_old - d_dot ) ) * J_ale , u_test ) + grad_p * grad_p_test * J_ale + )*JxW[quad]; + + // Navier Stokes: + res_u = J_ale * ( +// time derivative term + rho * scalar_product( + u_dot + grad_u * ( F_inv * ( u_old - d_dot ) ), u_test ) // - + scalar_product( J_ale * sigma * Ft_inv, grad_u_test ) + + scalar_product(sigma * Ft_inv, grad_u_test )) // divergence free constriant - - div_u * p_test -// pressure term - - p * div_u_test -// Impose armonicity of d and v=d_dot - + scalar_product( grad_d , grad_d_test ) - )*JxW[quad]; + - div_v * p_test + + (v - J_ale * F_inv * u_old) * v_test; + + // ALE: + res_d = +// Impose armonicity of d + scalar_product( grad_d , grad_d_test ); + residual[0][i] += (res_u + res_d)*JxW[quad]; } } @@ -345,31 +393,35 @@ ALENavierStokes::compute_system_operators( AMG_d.initialize_preconditioner( matrices[0]->block(0,0), fe, dh); AMG_u.initialize_preconditioner( matrices[0]->block(1,1), fe, dh); + AMG_v.initialize_preconditioner( matrices[0]->block(3,3), fe, dh); + AMG_Ap.initialize_preconditioner( matrices[2]->block(2,2), fe, dh); jac_M.initialize_preconditioner<>(matrices[1]->block(2,2)); //////////////////////////////////////////////////////////////////////////// // SYSTEM MATRIX: - std::array, 3>, 3> S; - for (unsigned int i = 0; i<3; ++i) - for (unsigned int j = 0; j<3; ++j) + std::array, 4>, 4> S; + for (unsigned int i = 0; i<4; ++i) + for (unsigned int j = 0; j<4; ++j) S[i][j] = linear_operator< BVEC >(matrices[0]->block(i,j) ); + + S[2][2] = null_operator< TrilinosWrappers::MPI::Vector >(S[2][2]); system_op = BlockLinearOperator< VEC >(S); //////////////////////////////////////////////////////////////////////////// // PRECONDITIONER MATRIX: - std::array, 3>, 3> P; - for (unsigned int i = 0; i<3; ++i) - for (unsigned int j = 0; j<3; ++j) + std::array, 4>, 4> P; + for (unsigned int i = 0; i<4; ++i) + for (unsigned int j = 0; j<4; ++j) P[i][j] = linear_operator( matrices[0]->block(i,j) ); static ReductionControl solver_control_pre(5000, 1e-8); static SolverCG solver_CG(solver_control_pre); static SolverGMRES solver_GMRES(solver_control_pre); - for (unsigned int i = 0; i<3; ++i) - for (unsigned int j = 0; j<3; ++j) + for (unsigned int i = 0; i<4; ++i) + for (unsigned int j = 0; j<4; ++j) if (i!=j) P[i][j] = null_operator< TrilinosWrappers::MPI::Vector >(P[i][j]); @@ -405,7 +457,21 @@ ALENavierStokes::compute_system_operators( jac_M); } - auto Schur_inv = nu * Mp_inv; + auto Ap = linear_operator< TrilinosWrappers::MPI::Vector >( matrices[2]->block(2,2) ); + LinearOperator Ap_inv; + if (AMG_Ap_use_inverse_operator) + { + Ap_inv = inverse_operator(Ap, + solver_GMRES, + jac_M); + } + else + { + Ap_inv = linear_operator(matrices[2]->block(2,2), + AMG_Ap); + } + auto alpha = this->get_alpha(); + auto Schur_inv = 1/nu * Mp_inv;// + rho * alpha * Ap_inv; if (AMG_d_use_inverse_operator) { @@ -419,6 +485,18 @@ ALENavierStokes::compute_system_operators( AMG_d); } + if (AMG_v_use_inverse_operator) + { + P[3][3] = inverse_operator( S[3][3], + solver_CG, + AMG_v); + } + else + { + P[3][3] = linear_operator(matrices[0]->block(3,3), + AMG_v); + } + P[1][1] = A_inv; P[1][2] = A_inv * Bt * Schur_inv; P[2][1] = null_operator(B); diff --git a/include/pidomus.h b/include/pidomus.h index 6e2b68c8..68e782c4 100644 --- a/include/pidomus.h +++ b/include/pidomus.h @@ -43,6 +43,7 @@ #include #include #include + #include #include #include diff --git a/include/pidomus_signals.h b/include/pidomus_signals.h index c33a553f..d1b1032c 100644 --- a/include/pidomus_signals.h +++ b/include/pidomus_signals.h @@ -94,6 +94,8 @@ struct Signals boost::signals2::signal begin_refine_and_transfer_solutions; boost::signals2::signal begin_assemble_matrices; boost::signals2::signal begin_solver_should_restart; + boost::signals2::signal begin_cycle; + boost::signals2::signal begin_run; boost::signals2::signal end_make_grid_fe; boost::signals2::signal end_setup_dofs; @@ -104,6 +106,8 @@ struct Signals boost::signals2::signal end_refine_and_transfer_solutions; boost::signals2::signal end_assemble_matrices; boost::signals2::signal end_solver_should_restart; + boost::signals2::signal end_cycle; + boost::signals2::signal end_run; }; diff --git a/source/pidomus.cc b/source/pidomus.cc index e3ab1b84..61fe44a7 100644 --- a/source/pidomus.cc +++ b/source/pidomus.cc @@ -193,12 +193,14 @@ piDoMUS::piDoMUS (const std::string &name, template void piDoMUS::run () { + signals.begin_run(); interface.set_stepper(time_stepper); interface.connect_to_signals(); for (current_cycle = 0; current_cycle < n_cycles; ++current_cycle) { + signals.begin_cycle(); if (current_cycle == 0) { make_grid_fe(); @@ -235,10 +237,11 @@ void piDoMUS::run () imex.solve_dae(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); - + signals.end_run(); } diff --git a/tests/ALE_navier_stokes_00.cc b/tests/ALE_navier_stokes_00.cc new file mode 100644 index 00000000..c92502f9 --- /dev/null +++ b/tests/ALE_navier_stokes_00.cc @@ -0,0 +1,35 @@ +#include "pidomus.h" +#include "interfaces/ALE_navier_stokes.h" +#include "tests.h" + +/** + * Test: ALE Navier Stokes interface. + */ + +using namespace dealii; +int main (int argc, char *argv[]) +{ + + Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv, + numbers::invalid_unsigned_int); + + initlog(); + deallog.depth_file(1); + deallog.threshold_double(1.0e-3); + + ALENavierStokes<2,2, LADealII> energy; + piDoMUS<2,2, LADealII> navier_stokes ("",energy); + ParameterAcceptor::initialize( + SOURCE_DIR "/parameters/ALE_navier_stokes_00.prm", + "used_parameters.prm"); + + navier_stokes.run (); + + auto sol = navier_stokes.get_solution(); + for (unsigned int i = 0 ; i energy; + piDoMUS<2,2> navier_stokes ("",energy); + ParameterAcceptor::initialize( + SOURCE_DIR "/parameters/ALE_navier_stokes_01.prm", + "used_parameters.prm"); + + navier_stokes.run (); + + auto sol = navier_stokes.get_solution(); + for (unsigned int i = 0 ; i energy; + piDoMUS<2,2> navier_stokes ("",energy); + ParameterAcceptor::initialize( + SOURCE_DIR "/parameters/ALE_navier_stokes_02.prm", + "used_parameters.prm"); + + navier_stokes.run (); + + auto sol = navier_stokes.get_solution(); + for (unsigned int i = 0 ; i 1 = 0,0 + set Optional Point 2 = 1,1 +end +subsection Error Tables + set Compute error = true + set Error file format = tex + set Error precision = 3 + set Output error tables = true + set Solution names = d,d,u,u,p,v,v + set Solution names for latex = d,d,u,u,p,v,v + set Table names = error + set Write error files = false + subsection Table 0 + set Add convergence rates = true + set Extra terms = cells,dofs + set Latex table caption = error + set List of error norms to compute = L2; AddUp; L2, Linfty, H1; AddUp; L2; L2, Linfty, H1; AddUp; + set Rate key = + end +end +subsection IMEX Parameters + set Absolute error tolerance = 1.0e-8 + set Final time = 0.4 + set Initial time = 0.0 + set Print useful informations = false + set Relative error tolerance = 1.0e-6 + set Step size = 2.0e-1 + set Use the KINSOL solver = true +end +subsection KINSOL for IMEX + set Level of verbosity of the KINSOL solver = 0 + set Maximum number of iteration before Jacobian update = 10 + set Maximum number of iterations = 200 + set Strategy = global_newton + set Use internal KINSOL direct solver = false +end +subsection piDoMUS<2, 2, LADealII> + set Adaptive refinement = false + set Initial global refinement = 5 + set Jacobian solver tolerance = 1e-6 + set Number of cycles = 1 + set Time stepper = euler + set Use direct solver if available = true +end \ No newline at end of file diff --git a/tests/parameters/ALE_navier_stokes_01.prm b/tests/parameters/ALE_navier_stokes_01.prm new file mode 100644 index 00000000..4e23de5e --- /dev/null +++ b/tests/parameters/ALE_navier_stokes_01.prm @@ -0,0 +1,93 @@ +# Problem +################################################################################ +subsection Exact solution + set Function constants = + set Function expression = 0;0;cos(x); y*sin(x); x*sin(y); 0;0 + set Variable names = x,y,t +end + +subsection Dirichlet boundary conditions + set IDs and component masks = 0=u;p;d + set IDs and expressions = 0=0;0;cos(x); y*sin(x); x*sin(y);0;0 + set Known component names = d,d,u,u,p,v,v + set Used constants = +end +subsection Forcing terms + set IDs and component masks = 0=u + set IDs and expressions = 0=0;0;\ + -sin(x)*cos(x) + sin(y) + 0.1*cos(x);\ + x*cos(y) + y*sin(x)*sin(x) + y*cos(x)*cos(x);\ + 0;0;0 + set Known component names = d,d,u,u,p,v,v + set Used constants = k=0.0 +end + +# Solver +################################################################################ +subsection ALE Navier Stokes Interface + set AMG d - use inverse operator = true + set AMG v - use inverse operator = true + set AMG u - use inverse operator = true + set Invert Mp using inverse operator = true + set nu [Pa s] = 0.1 + set rho [Kg m^-d] = 1.0 +end +subsection AMG for d + set Aggregation threshold = 1 +end +subsection AMG for v + set Aggregation threshold = 1 +end +subsection AMG for u + set Aggregation threshold = 1 +end +subsection Jacobi for M + set Omega = 1.4 +end +subsection Domain + set Colorize = false + set Grid to generate = rectangle + set Optional Point 1 = 0,0 + set Optional Point 2 = 1,1 +end +subsection Error Tables + set Compute error = true + set Error file format = tex + set Error precision = 3 + set Output error tables = true + set Solution names = d,d,u,u,p,v,v + set Solution names for latex = d,d,u,u,p,v,v + set Table names = error + set Write error files = false + subsection Table 0 + set Add convergence rates = true + set Extra terms = cells,dofs + set Latex table caption = error + set List of error norms to compute = L2; AddUp; L2, Linfty, H1; AddUp; L2; L2, Linfty, H1; AddUp; + set Rate key = + end +end +subsection IMEX Parameters + set Absolute error tolerance = 1.0e-8 + set Final time = 0.4 + set Initial time = 0.0 + set Print useful informations = false + set Relative error tolerance = 1.0e-6 + set Step size = 2.0e-1 + set Use the KINSOL solver = true +end +subsection KINSOL for IMEX + set Level of verbosity of the KINSOL solver = 0 + set Maximum number of iteration before Jacobian update = 1 + set Maximum number of iterations = 200 + set Strategy = global_newton + set Use internal KINSOL direct solver = false +end +subsection piDoMUS<2, 2, LATrilinos> + set Adaptive refinement = false + set Initial global refinement = 3 + set Jacobian solver tolerance = 1e-6 + set Number of cycles = 1 + set Time stepper = euler + set Use direct solver if available = false +end \ No newline at end of file diff --git a/tests/parameters/ALE_navier_stokes_02.prm b/tests/parameters/ALE_navier_stokes_02.prm new file mode 100644 index 00000000..895ad64b --- /dev/null +++ b/tests/parameters/ALE_navier_stokes_02.prm @@ -0,0 +1,106 @@ +# Problem +################################################################################ +subsection Exact solution + set Function constants = + set Function expression = 0;0;cos(x); y*sin(x); x*sin(y);0;0 + set Variable names = x,y,t +end + +subsection Dirichlet boundary conditions + set IDs and component masks = 0=u;p;d + set IDs and expressions = 0=sin(x*t);0;cos(x); y*sin(x); x*sin(y);0;0 + set Known component names = d,d,u,u,p,v,v + set Used constants = +end +subsection Forcing terms + set IDs and component masks = 0=u + set IDs and expressions = 0=0;0;\ + -sin(x)*cos(x) + sin(y) + 0.1*cos(x);\ + x*cos(y) + y*sin(x)*sin(x) + y*cos(x)*cos(x);\ + 0;0;0 + set Known component names = d,d,u,u,p,v,v + set Used constants = k=0.0 +end + +# Solver +################################################################################ +subsection ALE Navier Stokes Interface + set AMG d - use inverse operator = false + set AMG u - use inverse operator = false + set AMG v - use inverse operator = false + set Invert Mp using inverse operator = true + set nu [Pa s] = 10 + set rho [Kg m^-d] = 1.0 +end +subsection AMG for Ap + set Aggregation threshold = 5 +end +subsection AMG for d + set Aggregation threshold = 10 +end +subsection AMG for u + set Aggregation threshold = 10 +end +subsection AMG for v + set Aggregation threshold = 10 +end +subsection Jacobi for M + set Omega = 1.4 +end +subsection Domain + set Colorize = false + set Grid to generate = rectangle + set Optional Point 1 = 0,0 + set Optional Point 2 = 1,1 +end +subsection Error Tables + set Compute error = true + set Error file format = tex + set Error precision = 3 + set Output error tables = true + set Solution names = d,d,u,u,p,v,v + set Solution names for latex = d,d,u,u,p,v,v + set Table names = error + set Write error files = false + subsection Table 0 + set Add convergence rates = true + set Extra terms = cells,dofs + set Latex table caption = error + set List of error norms to compute = L2; AddUp; L2, Linfty, H1; AddUp; L2; L2, Linfty, H1; AddUp; + set Rate key = + end +end +subsection IMEX Parameters + set Absolute error tolerance = 1.0e-8 + set Final time = 0.4 + set Initial time = 0.0 + set Print useful informations = false + set Relative error tolerance = 1.0e-6 + set Step size = 1.0e-3 + set Use the KINSOL solver = false +end +subsection KINSOL for IMEX + set Level of verbosity of the KINSOL solver = 0 + set Maximum number of iteration before Jacobian update = 20 + set Maximum number of iterations = 200 + set Strategy = global_newton + set Use internal KINSOL direct solver = false +end +subsection piDoMUS<2, 2, LATrilinos> + set Adaptive refinement = false + set Initial global refinement = 6 + set Jacobian solver tolerance = 1e-4 + set Number of cycles = 1 + set Time stepper = euler + set Use direct solver if available = true +end + +subsection Output Parameters + set Files to save in run directory = + set Incremental run prefix = + set Output format = vtu + set Output partitioning = false + set Problem base name = solution + set Solution names = u + set Subdivisions = 1 +end \ No newline at end of file diff --git a/tests/parameters/navier_stokes_04.prm b/tests/parameters/navier_stokes_04.prm index b19ec65e..8625b104 100644 --- a/tests/parameters/navier_stokes_04.prm +++ b/tests/parameters/navier_stokes_04.prm @@ -78,7 +78,7 @@ end subsection Forcing terms set IDs and component masks = 0=u set IDs and expressions = 0=\ - 2*k*cos(x)*cos(y)-sin(x)*sin(y)^2*cos(x)-sin(x)*cos(x)*cos(y)^2+cos(x)*cos(y); \ + 2*k*cos(x)*cos(y)-sin(x)*sin(y)*sin(y)*cos(x)-sin(x)*cos(x)*cos(y)^2+cos(x)*cos(y); \ 2*k*sin(x)*sin(y)+sin(x)^2*sin(y)*cos(y)-sin(x)*sin(y)+ sin(y)*cos(x)^2*cos(y); 0 set Known component names = u,u,p set Used constants = k=1.0 diff --git a/tests/tests_to_check/ode_argument_trilinos_01.cc b/tests/tests_to_check/ode_argument_trilinos_01.cc index 66ac240e..f357cc81 100644 --- a/tests/tests_to_check/ode_argument_trilinos_01.cc +++ b/tests/tests_to_check/ode_argument_trilinos_01.cc @@ -18,7 +18,7 @@ #include -#include + #include #include