diff --git a/include/interfaces/navier_stokes-chorin-temam.h b/include/interfaces/navier_stokes-chorin-temam.h new file mode 100644 index 00000000..cef819ce --- /dev/null +++ b/include/interfaces/navier_stokes-chorin-temam.h @@ -0,0 +1,324 @@ +/*! \addtogroup equations + * @{ + */ + +/** + * Navier Stokes Equation using Chorin-Temam projection method + */ + +#ifndef _pidomus_navier_stokes_h_ +#define _pidomus_navier_stokes_h_ + +#include "pde_system_interface.h" + +#include + +#include + +#include +#include + +#include +#include +#include + +//////////////////////////////////////////////////////////////////////////////// +/// Navier Stokes interface: + +template +class NavierStokes + : + public PDESystemInterface, LAC> +{ + +public: + ~NavierStokes () {}; + NavierStokes (); + + void declare_parameters (ParameterHandler &prm); + + template + void + energies_and_residuals( + const typename DoFHandler::active_cell_iterator &cell, + FEValuesCache &scratch, + std::vector &energies, + std::vector> &residuals, + bool compute_only_system_terms) const; + + void + compute_system_operators( + const std::vector>, + LinearOperator &, + LinearOperator &, + LinearOperator &) const; + + void + set_matrix_couplings(std::vector &couplings) const; + +private: + + /** + * AMG preconditioner for velocity-velocity matrix. + */ + mutable ParsedJacobiPreconditioner AMG_u; + + /** + * AMG preconditioner for velocity-velocity matrix. + */ + mutable ParsedJacobiPreconditioner AMG_v; + + /** + * AMG preconditioner for velocity-velocity matrix. + */ + mutable ParsedAMGPreconditioner AMG_p; + + // PHYSICAL PARAMETERS: + //////////////////////////////////////////// + + /** + * Density + */ + double rho; + + /** + * Viscosity + */ + double nu; + + double initial_delta_t; + /** + * Solver tolerance for CG + */ + double CG_solver_tolerance; + + /** + * Solver tolerance for GMRES + */ + double GMRES_solver_tolerance; +}; + +template +NavierStokes:: +NavierStokes() + : + PDESystemInterface, LAC>( + "Navier Stokes Interface", + dim+dim+1, + 1, + "FESystem[FE_Q(2)^d-FE_Q(1)-FE_Q(2)^d]", + (dim==3)?"v,v,v,p,u,u,u":"v,v,p,u,u", + "1,0,1"), + AMG_u("Prec for u", 1.4), + AMG_v("Prec for v", 1.4), + AMG_p("AMG for p") +{ + this->init(); +} + +template +void NavierStokes:: +set_matrix_couplings(std::vector &couplings) const +{ + couplings[0] = "1,1,1;1,1,1;1,1,1"; +} + +template +void NavierStokes:: +declare_parameters (ParameterHandler &prm) +{ + PDESystemInterface,LAC>:: + declare_parameters(prm); + this->add_parameter(prm, &rho, + "rho [kg m^3]", "1.0", + Patterns::Double(0.0), + "Density"); + this->add_parameter(prm, &nu, + "nu [Pa s]", "1.0", + Patterns::Double(0.0), + "Viscosity"); + this->add_parameter(prm, &initial_delta_t, + "initial delta t", "1e-1", + Patterns::Double(0.0), + "Initial Delta t"); + this->add_parameter(prm, &CG_solver_tolerance, + "CG Solver tolerance", "1e-8", + Patterns::Double(0.0)); + this->add_parameter(prm, &GMRES_solver_tolerance, + "GMRES Solver tolerance", "1e-8", + Patterns::Double(0.0)); +} + +template +template +void +NavierStokes:: +energies_and_residuals(const typename DoFHandler::active_cell_iterator &cell, + FEValuesCache &fe_cache, + std::vector &, + std::vector > &residual, + bool compute_only_system_terms) const +{ + const FEValuesExtractors::Vector aux_velocity(0); + const FEValuesExtractors::Scalar pressure(dim); + const FEValuesExtractors::Vector velocity(dim+1); + + ResidualType et = 0; + double dummy = 0.0; + + this->reinit (et, cell, fe_cache); + + // Velocity: + auto &us = fe_cache.get_values("solution", "u", velocity, et); + auto &ues = fe_cache.get_values("explicit_solution", "ue", velocity, dummy); + auto &grad_ues = fe_cache.get_gradients("explicit_solution", "grad_u", velocity, dummy); + auto &div_us = fe_cache.get_divergences("solution", "div_u", velocity, et); + auto &sym_grad_ues = fe_cache.get_symmetric_gradients("explicit_solution", "sym_grad_ue", velocity, dummy); + + auto &vs = fe_cache.get_values("solution", "v", aux_velocity, et); + auto &ves = fe_cache.get_values("explicit_solution", "ve", aux_velocity, dummy); + auto &sym_grad_vs = fe_cache.get_symmetric_gradients("solution", "sym_grad_v", aux_velocity, et); + auto &grad_vs = fe_cache.get_gradients("solution", "grad_v", aux_velocity, et); + auto &div_vs = fe_cache.get_divergences("solution", "div_v", aux_velocity, et); + + // Pressure: + auto &ps = fe_cache.get_values("solution", "p", pressure, et); + auto &pes = fe_cache.get_values("explicit_solution", "pe", pressure, dummy); + auto &grad_ps = fe_cache.get_gradients("solution", "grad_p", pressure,et); + auto &grad_pes = fe_cache.get_gradients("explicit_solution", "grad_pe", pressure,dummy); + + + const unsigned int n_quad_points = us.size(); + auto &JxW = fe_cache.get_JxW_values(); + const auto delta_t = this->get_timestep() != this->get_timestep() ? initial_delta_t : this->get_timestep() ; + auto &fev = fe_cache.get_current_fe_values(); + + for (unsigned int quad=0; quad &grad_dp = grad_dps[quad]; + + // Pressure: + const ResidualType &p = ps[quad]; + const double &pe = pes[quad]; + const Tensor<1, dim, ResidualType> &grad_p = grad_ps[quad]; + const Tensor<1, dim, double> &grad_pe = grad_pes[quad]; + + // Velocity: + const Tensor<1, dim, ResidualType> &u = us[quad]; + const Tensor<1, dim, double> &ue = ues[quad]; + const Tensor<2, dim, double> &grad_ue = grad_ues[quad]; + const ResidualType &div_u = div_us[quad]; + + const Tensor<1, dim, ResidualType> &v = vs[quad]; + const Tensor<1, dim, double> vd = SacadoTools::to_double(v); + const Tensor<2, dim, ResidualType> &grad_v = grad_vs[quad]; + const Tensor<2, dim, ResidualType> &sym_grad_v = sym_grad_vs[quad]; + const Tensor<2, dim, double> &sym_grad_ue = sym_grad_ues[quad]; + const Tensor<1, dim, double> &ve = ves[quad]; + const ResidualType &div_v = div_vs[quad]; + + for (unsigned int i=0; i +void +NavierStokes::compute_system_operators( + const std::vector> matrices, + LinearOperator &system_op, + LinearOperator &prec_op, + LinearOperator &prec_op_finer) const +{ + const unsigned int num_blocks = 3; + typedef LATrilinos::VectorType::BlockType BVEC; + typedef LATrilinos::VectorType VEC; + + static ReductionControl solver_control(matrices[0]->m(), CG_solver_tolerance); + static SolverCG solver_CG(solver_control); + + const DoFHandler &dh = this->get_dof_handler(); + const ParsedFiniteElement fe = this->pfe; + AMG_v.initialize_preconditioner( matrices[0]->block(0,0)); //, fe, dh); + AMG_p.initialize_preconditioner( matrices[0]->block(1,1), fe, dh); + AMG_u.initialize_preconditioner( matrices[0]->block(2,2)); //, fe, dh); + + //////////////////////////////////////////////////////////////////////////// + // SYSTEM MATRIX: + //////////////////////////////////////////////////////////////////////////// + std::array, num_blocks>, num_blocks> S; + for (unsigned int i = 0; i(matrices[0]->block(i,j) ); + system_op = BlockLinearOperator< VEC >(S); + + + //////////////////////////////////////////////////////////////////////////// + // PRECONDITIONER MATRIX: + //////////////////////////////////////////////////////////////////////////// + auto Av = linear_operator(matrices[0]->block(0,0)); + auto Av_inv = inverse_operator(Av, solver_CG, AMG_v); + + auto Ap = linear_operator(matrices[0]->block(1,1)); + auto Ap_inv = inverse_operator(Ap, solver_CG, AMG_p); + + auto Au = linear_operator(matrices[0]->block(2,2)); + auto Au_inv = inverse_operator(Au, solver_CG, AMG_u); + + // Preconditioner + ////////////////////////////// + + BlockLinearOperator diag_inv + = block_diagonal_operator( + { + { + Av_inv, Ap_inv, Au_inv + } + } + ); + prec_op = diag_inv; + + /* prec_op = block_forward_substitution( */ + /* BlockLinearOperator< VEC >(S), */ + /* diag_inv); */ +} + +#endif + +/*! @} */ diff --git a/source/pidomus.cc b/source/pidomus.cc index 0d6d530b..80771c51 100644 --- a/source/pidomus.cc +++ b/source/pidomus.cc @@ -405,7 +405,7 @@ piDoMUS::solve_jacobian_system(const typename LAC::VectorTyp if (we_are_parallel == false && use_direct_solver == true) { - + std::cout << "DIRETTO!" << std::endl; SparseDirectUMFPACK inverse; inverse.factorize((sMAT &) *matrices[0]); inverse.vmult((sVEC &)dst, (sVEC &)src); diff --git a/tests/navier_stokes-chorin-temam_00.cc b/tests/navier_stokes-chorin-temam_00.cc new file mode 100644 index 00000000..0f012b63 --- /dev/null +++ b/tests/navier_stokes-chorin-temam_00.cc @@ -0,0 +1,43 @@ +#include "pidomus.h" +#include "interfaces/navier_stokes-chorin-temam.h" +#include "tests.h" + +/** + * Test: Navier Stokes interface. + * Method: Direct + * Problem: Non time depending Navier Stokes Equations + * Exact solution: + * \f[ + * u=\big( 2*(x^2)*y, -2*x*(y^2) \big) + * \textrm{ and }p=xy; + * \f] + */ + +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); + + const int dim = 2; + + NavierStokes energy; + piDoMUS navier_stokes ("",energy); + ParameterAcceptor::initialize( + SOURCE_DIR "/parameters/navier_stokes_chorin-temam_00.prm", + "used_parameters.prm"); + + navier_stokes.run (); + + auto sol = navier_stokes.get_solution(); + for (unsigned int i = 0 ; i energy; + piDoMUS<2,2,LATrilinos> navier_stokes ("",energy); + ParameterAcceptor::initialize( + SOURCE_DIR "/parameters/navier_stokes_chorin-temam_01.prm", + "used_parameters.prm"); + + navier_stokes.run (); + + auto sol = navier_stokes.get_solution(); + for (unsigned int i = 0 ; i set Adaptive refinement = false set Enable finer preconditioner = false - set Initial global refinement = 4 + set Initial global refinement = 5 set Jacobian solver tolerance = 1e-6 set Max iterations = 50 set Max iterations finer prec. = 0 - set Number of cycles = 3 + set Number of cycles = 1 set Overwrite Newton's iterations = true set Print some useful informations about processes = true set Refine mesh during transient = false diff --git a/tests/parameters/navier_stokes_chorin-temam_00.prm b/tests/parameters/navier_stokes_chorin-temam_00.prm new file mode 100644 index 00000000..3040a3e2 --- /dev/null +++ b/tests/parameters/navier_stokes_chorin-temam_00.prm @@ -0,0 +1,86 @@ +# Listing of Parameters +# --------------------- +subsection Exact solution + set Function expression = t*cos(x); t*y*sin(x); -1*y*sin(x);\ + t*cos(x); t*y*sin(x) + set Variable names = x,y,t +end + +subsection Forcing terms + set IDs and component masks = 0=v + set IDs and expressions = 0 = \ + -t^2*sin(2*x)/2 + 0.1*t*cos(x) + sin(y) + cos(x);\ + t^2*y + x*cos(y) + y*sin(x);\ + 0; 0; 0 + set Known component names = v,v,p,u,u + set Used constants = d=0.01 +end + +subsection Dirichlet boundary conditions + set IDs and component masks = 0 = v;p + set IDs and expressions = 0 = t*cos(x); t*y*sin(x); -1*y*sin(x);\ + t*cos(x); t*y*sin(x) + set Known component names = v,v,p,u,u + set Used constants = k=0.01 +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 IMEX Parameters + set Final time = 1e0 + set Initial time = 0.0 + set Step size = 1e-1 +end + + +subsection Navier Stokes Interface + set Block of differential components = 1,0,1 + set Blocking of the finite element = v,v,p,u,u + set CG Solver tolerance = 1e-12 + set GMRES Solver tolerance = 1e-12 + set Finite element space = FESystem[FE_Q(2)^d-FE_Q(1)-FE_Q(2)^d] + set nu [Pa s] = 0.1 + set rho [kg m^3] = 1.0 + set initial delta t = 1e-3 +end + +subsection Output Parameters + set Files to save in run directory = + set Incremental run prefix = + set Output format = vtu + set Output partitioning = true + set Problem base name = Navier_Stokes_chorin-temam +end + + +subsection piDoMUS<2, 2, LADealII> + set Adaptive refinement = false + set Enable finer preconditioner = false + set Initial global refinement = 3 + set Number of cycles = 1 + set Jacobian solver tolerance = 1e-12 + set Use direct solver if available = true +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 = v,v,p,u,u + set Solution names for latex = v,v,p,u,u + 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 = Linfty; AddUp;L2; L2, Linfty, H1; AddUp + set Rate key = + end +end diff --git a/tests/parameters/navier_stokes_chorin-temam_01.prm b/tests/parameters/navier_stokes_chorin-temam_01.prm new file mode 100644 index 00000000..a65bc953 --- /dev/null +++ b/tests/parameters/navier_stokes_chorin-temam_01.prm @@ -0,0 +1,168 @@ +# Parameter file generated with +# D2K_GIT_BRANCH= master +# D2K_GIT_SHORTREV= 9558556 +# DEAL_II_GIT_BRANCH= spherical-manifold +# DEAL_II_GIT_SHORTREV= 8095e57 +subsection AMG for p + set Aggregation threshold = 0.1 + set Elliptic = true +end +subsection Dirichlet boundary conditions + set IDs and component masks = 0 = u;p % 1 = u % 2 = u % 3 = u + set IDs and expressions = 0 = t*cos(x); t*y*sin(x); -1*y*sin(x);t*cos(x); t*y*sin(x) %1 = t*cos(x); t*y*sin(x); -1*y*sin(x);t*cos(x); t*y*sin(x) %2 = t*cos(x); t*y*sin(x); -1*y*sin(x);t*cos(x); t*y*sin(x) %3 = t*cos(x); t*y*sin(x); -1*y*sin(x);t*cos(x); t*y*sin(x) + set Known component names = v,v,p,u,u + set Used constants = k=0.01 +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 = v,v,p,u,u + set Solution names for latex = v,v,p,u,u + 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 = Linfty; AddUp;L2; L2, Linfty, H1; AddUp + set Rate key = + end +end +subsection Exact solution + set Function constants = + set Function expression = t*cos(x); t*y*sin(x); -1*y*sin(x);t*cos(x); t*y*sin(x) + set Variable names = x,y,t +end +subsection Forcing terms + set IDs and component masks = 0=v + set IDs and expressions = 0 = -t^2*sin(2*x)/2 + 0.1*t*cos(x) + sin(y) + cos(x);t^2*y + x*cos(y) + y*sin(x);0; 0; 0 + set Known component names = v,v,p,u,u + set Used constants = d=0.01 +end +subsection IDA Solver Parameters + set Absolute error tolerance = 1e-4 + set Final time = 1 + set Ignore algebraic terms for error computations = false + set Initial condition Newton max iterations = 5 + set Initial condition Newton parameter = 0.33 + set Initial condition type = use_y_diff + set Initial condition type after restart = use_y_dot + set Initial step size = 1e-4 + set Initial time = 0 + set Maximum number of nonlinear iterations = 10 + set Maximum order of BDF = 5 + set Min step size = 5e-5 + set Relative error tolerance = 1e-3 + set Seconds between each output = 1e-1 + set Show output of time steps = true + set Use local tolerances = false +end +subsection IMEX Parameters + set Absolute error tolerance = 1e-8 + set Final time = 1e0 + set Initial time = 0.0 + set Intervals between outputs = 1 + set Maximum number of inner nonlinear iterations = 3 + set Maximum number of outer nonlinear iterations = 5 + set Method used = fixed_alpha + set Newton relaxation parameter = 1 + set Number of elements in backtracking sequence = 5 + set Print useful informations = false + set Relative error tolerance = 1e-6 + set Step size = 1e-4 + set Update continuously Jacobian = false + set Use the KINSOL solver = true +end +subsection Initial solution + set Function constants = + set Function expression = 0; 0; 0; 0; 0 + set Variable names = x,y,t +end +subsection Initial solution_dot + set Function constants = + set Function expression = 0; 0; 0; 0; 0 + set Variable names = x,y,t +end +subsection KINSOL for IMEX + set Level of verbosity of the KINSOL solver = 0 + set Maximum number of iteration before Jacobian update = 200 + set Maximum number of iterations = 200 + set Step tolerance = 1e-8 + set Strategy = global_newton + set Tolerance for residuals = 1e-6 + set Use internal KINSOL direct solver = false +end +subsection Navier Stokes Interface + set Block of differential components = 1,0,1 + set Blocking of the finite element = v,v,p,u,u + set CG Solver tolerance = 1e-8 + set Finite element space = FESystem[FE_Q(2)^d-FE_Q(1)-FE_Q(2)^d] + set GMRES Solver tolerance = 1e-12 + set initial delta t = 1e-4 + set nu [Pa s] = 0.1 + set rho [kg m^3] = 1.0 +end +subsection Neumann boundary conditions + set IDs and component masks = + set IDs and expressions = + set Known component names = v,v,p,u,u + set Used constants = +end +subsection Output Parameters + set Files to save in run directory = + set Incremental run prefix = + set Output format = vtu + set Output partitioning = true + set Problem base name = Navier_Stokes_chorin-temam + set Solution names = u + set Subdivisions = 1 +end +subsection Prec for u + set Min Diagonal = 0.000000 + set Number of sweeps = 2 + set Omega = 1.400000 +end +subsection Prec for v + set Min Diagonal = 0.000000 + set Number of sweeps = 2 + set Omega = 1.400000 +end +subsection Refinement + set Bottom fraction = 0.100000 + set Maximum number of cells (if available) = 0 + set Order (optimize) = 2 + set Refinement strategy = fraction + set Top fraction = 0.300000 +end +subsection Time derivative of Dirichlet boundary conditions + set IDs and component masks = + set IDs and expressions = + set Known component names = v,v,p,u,u + set Used constants = +end +subsection Zero average constraints + set Known component names = v,v,p,u,u + set Zero average on boundary = + set Zero average on whole domain = +end +subsection piDoMUS<2, 2, LATrilinos> + set Adaptive refinement = true + set Enable finer preconditioner = false + set Initial global refinement = 6 + set Jacobian solver tolerance = 1e-6 + set Number of cycles = 1 + set Overwrite Newton's iterations = true + set Print some useful informations about processes = true + set Refine mesh during transient = false + set Threshold for solver's restart = 1e-2 + set Time stepper = imex + set Use direct solver if available = false +end