From ceda0772e80459fe71b9d4aae24b5b9cef3562c3 Mon Sep 17 00:00:00 2001 From: Wolfgang Bangerth Date: Wed, 9 Mar 2011 12:43:38 +0000 Subject: [PATCH] Make the discretization run time parameters and allow for locally conservative discretizations. Adjust a couple things when the input file is missing to avoid strange exceptions. git-svn-id: https://svn.dealii.org/trunk@23471 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/examples/step-32/step-32.cc | 100 +++++++++++++++++++-------- deal.II/examples/step-32/step-32.prm | 34 ++++++--- 2 files changed, 96 insertions(+), 38 deletions(-) diff --git a/deal.II/examples/step-32/step-32.cc b/deal.II/examples/step-32/step-32.cc index 473e8c26df..1dff0e7722 100644 --- a/deal.II/examples/step-32/step-32.cc +++ b/deal.II/examples/step-32/step-32.cc @@ -59,6 +59,7 @@ #include #include +#include #include #include #include @@ -928,9 +929,14 @@ class BoussinesqFlowProblem unsigned int adaptive_refinement_interval; - double stabilization_alpha; - double stabilization_c_R; - double stabilization_beta; + double stabilization_alpha; + double stabilization_c_R; + double stabilization_beta; + + unsigned int stokes_velocity_degree; + bool use_locally_conservative_discretization; + + unsigned int temperature_degree; }; @@ -942,7 +948,6 @@ class BoussinesqFlowProblem const MappingQ mapping; - const unsigned int stokes_degree; const FESystem stokes_fe; DoFHandler stokes_dof_handler; @@ -956,7 +961,6 @@ class BoussinesqFlowProblem TrilinosWrappers::MPI::BlockVector stokes_rhs; - const unsigned int temperature_degree; FE_Q temperature_fe; DoFHandler temperature_dof_handler; ConstraintMatrix temperature_constraints; @@ -1044,7 +1048,10 @@ BoussinesqFlowProblem::Parameters::Parameters () adaptive_refinement_interval (10), stabilization_alpha (2), stabilization_c_R (0.11), - stabilization_beta (0.078) + stabilization_beta (0.078), + stokes_velocity_degree (2), + use_locally_conservative_discretization (true), + temperature_degree (2) {} @@ -1096,6 +1103,25 @@ declare_parameters (ParameterHandler &prm) "and 0.078 for 3d."); } prm.leave_subsection (); + + prm.enter_subsection ("Discretization"); + { + prm.declare_entry ("Stokes velocity polynomial degree", "2", + Patterns::Integer (1), + "The polynomial degree to use for the velocity variables " + "in the Stokes system."); + prm.declare_entry ("Temperature polynomial degree", "2", + Patterns::Integer (1), + "The polynomial degree to use for the temperature variable."); + prm.declare_entry ("Use locally conservative discretization", "true", + Patterns::Bool (), + "Whether to use a Stokes discretization that is locally " + "conservative at the expense of a larger number of degrees " + "of freedom, or to go with a cheaper discretization " + "that does not locally conserve mass (although it is " + "globally conservative."); + } + prm.leave_subsection (); } @@ -1121,6 +1147,15 @@ parse_parameters (ParameterHandler &prm) stabilization_beta = prm.get_double ("beta"); } prm.leave_subsection (); + + prm.enter_subsection ("Discretization"); + { + stokes_velocity_degree = prm.get_integer ("Stokes velocity polynomial degree"); + temperature_degree = prm.get_integer ("Temperature polynomial degree"); + use_locally_conservative_discretization + = prm.get_bool ("Use locally conservative discretization"); + } + prm.leave_subsection (); } @@ -1184,14 +1219,20 @@ BoussinesqFlowProblem::BoussinesqFlowProblem () mapping (4), - stokes_degree (1), - stokes_fe (FE_Q(stokes_degree+1), dim, - FE_Q(stokes_degree), 1), + stokes_fe (FE_Q(parameters.stokes_velocity_degree), + dim, + (parameters.use_locally_conservative_discretization + ? + static_cast &> + (FE_DGP(parameters.stokes_velocity_degree-1)) + : + static_cast &> + (FE_Q(parameters.stokes_velocity_degree-1))), + 1), stokes_dof_handler (triangulation), - temperature_degree (2), - temperature_fe (temperature_degree), + temperature_fe (parameters.temperature_degree), temperature_dof_handler (triangulation), time_step (0), @@ -1260,7 +1301,7 @@ template double BoussinesqFlowProblem::get_maximal_velocity () const { const QIterated quadrature_formula (QTrapez<1>(), - stokes_degree+1); + parameters.stokes_velocity_degree); const unsigned int n_q_points = quadrature_formula.size(); FEValues fe_values (mapping, stokes_fe, quadrature_formula, update_values); @@ -1313,7 +1354,7 @@ std::pair BoussinesqFlowProblem::get_extrapolated_temperature_range () const { const QIterated quadrature_formula (QTrapez<1>(), - temperature_degree); + parameters.temperature_degree); const unsigned int n_q_points = quadrature_formula.size(); FEValues fe_values (mapping, temperature_fe, quadrature_formula, @@ -1568,7 +1609,7 @@ void BoussinesqFlowProblem::project_temperature_field () { assemble_temperature_matrix (); - QGauss quadrature(temperature_degree+2); + QGauss quadrature(parameters.temperature_degree+2); UpdateFlags update_flags = UpdateFlags(update_values | update_quadrature_points | update_JxW_values); @@ -2187,7 +2228,7 @@ BoussinesqFlowProblem::assemble_stokes_preconditioner () { stokes_preconditioner_matrix = 0; - const QGauss quadrature_formula(stokes_degree+2); + const QGauss quadrature_formula(parameters.stokes_velocity_degree+1); typedef FilteredIterator::active_cell_iterator> @@ -2394,7 +2435,7 @@ void BoussinesqFlowProblem::assemble_stokes_system () stokes_rhs=0; - const QGauss quadrature_formula(stokes_degree+2); + const QGauss quadrature_formula(parameters.stokes_velocity_degree+1); typedef FilteredIterator::active_cell_iterator> @@ -2521,7 +2562,7 @@ void BoussinesqFlowProblem::assemble_temperature_matrix () temperature_mass_matrix = 0; temperature_stiffness_matrix = 0; - const QGauss quadrature_formula(temperature_degree+2); + const QGauss quadrature_formula(parameters.temperature_degree+2); typedef FilteredIterator::active_cell_iterator> @@ -2802,7 +2843,7 @@ void BoussinesqFlowProblem::assemble_temperature_system (const double maxim temperature_rhs = 0; - const QGauss quadrature_formula(temperature_degree+2); + const QGauss quadrature_formula(parameters.temperature_degree+2); const std::pair global_T_range = get_extrapolated_temperature_range(); @@ -2964,7 +3005,7 @@ void BoussinesqFlowProblem::solve () // size in 3d double scaling = (dim==3)?0.25:1.0; double local_time_step = scaling/(1.6*dim*std::sqrt(1.*dim)) / - temperature_degree * + parameters.temperature_degree * GridTools::minimal_cell_diameter(triangulation) / std::max(1e-10,maximal_velocity); @@ -3261,7 +3302,7 @@ void BoussinesqFlowProblem::output_results () //norm of gradient { double my_cells_error[2] = {0, 0}; - QGauss<1> q_base(stokes_degree+1); + QGauss<1> q_base(parameters.stokes_velocity_degree); QIterated err_quadrature(q_base, 2); const unsigned int n_q_points = err_quadrature.size(); @@ -3523,7 +3564,7 @@ void BoussinesqFlowProblem::refine_mesh (const unsigned int max_grid_level) Vector estimated_error_per_cell (triangulation.n_active_cells()); KellyErrorEstimator::estimate (temperature_dof_handler, - QGauss(temperature_degree+1), + QGauss(parameters.temperature_degree+1), typename FunctionMap::type(), temperature_solution, estimated_error_per_cell, @@ -3639,17 +3680,18 @@ void BoussinesqFlowProblem::run (const std::string parameter_filename) { parameter_file.close (); - std::cerr << "*** Input parameter file <" - << parameter_filename << "> not found. Creating a" - << std::endl - << "*** template file of the same name." - << std::endl; + std::ostringstream message; + message << "Input parameter file <" + << parameter_filename << "> not found. Creating a" + << std::endl + << "template file of the same name." + << std::endl; std::ofstream parameter_out (parameter_filename.c_str()); prm.print_parameters (parameter_out, ParameterHandler::Text); - std::exit (1); + AssertThrow (false, ExcMessage (message.str().c_str())); } const bool success = prm.read_input (parameter_file); @@ -3765,12 +3807,12 @@ void BoussinesqFlowProblem::run (const std::string parameter_filename) // This is copied verbatim from step-31: int main (int argc, char *argv[]) { + Utilities::System::MPI_InitFinalize mpi_initialization(argc, argv); + try { deallog.depth_console (0); - Utilities::System::MPI_InitFinalize mpi_initialization(argc, argv); - std::string parameter_filename; if (argc>=2) parameter_filename = argv[1]; diff --git a/deal.II/examples/step-32/step-32.prm b/deal.II/examples/step-32/step-32.prm index b8777bc30a..40b4d05057 100644 --- a/deal.II/examples/step-32/step-32.prm +++ b/deal.II/examples/step-32/step-32.prm @@ -3,36 +3,52 @@ # The end time of the simulation in years. set End time = 1e8 +# Whether graphical output is to be generated or not. You may not want to get +# graphical output if the number of processors is large. +set Generate graphical output = false + # The number of adaptive refinement steps performed after initial global # refinement. set Initial adaptive refinement = 2 # The number of global refinement steps performed on the initial coarse mesh, # before the problem is first solved there. -set Initial global refinement = 1 +set Initial global refinement = 2 + +# The number of time steps between each generation of graphical output files. +set Time steps between graphical output = 50 # The number of time steps after which the mesh is to be adapted based on # computed error indicators. set Time steps between mesh refinement = 10 -# Whether graphical output is to be generated or not. You may not want to get -# graphical output if the number of processors is large. -set Generate graphical output = false -# The number of time steps between each generation of graphical output files. -set Time steps between graphical output = 50 +subsection Discretization + # The polynomial degree to use for the velocity variables in the Stokes + # system. + set Stokes velocity polynomial degree = 2 + + # The polynomial degree to use for the temperature variable. + set Temperature polynomial degree = 2 + + # Whether to use a Stokes discretization that is locally conservative at the + # expense of a larger number of degrees of freedom, or to go with a cheaper + # discretization that does not locally conserve mass (although it is + # globally conservative. + set Use locally conservative discretization = true +end subsection Stabilization parameters # The exponent in the entropy viscosity stabilization. set alpha = 2 - # The c_R factor in the entropy viscosity stabilization. - set c_R = 0.11 - # The beta factor in the artificial viscosity stabilization. An appropriate # value for 2d is 0.052 and 0.078 for 3d. set beta = 0.078 + + # The c_R factor in the entropy viscosity stabilization. + set c_R = 0.11 end -- 2.39.5