From 1ea5fbee3e653c24889bbb958676349baf22d90a Mon Sep 17 00:00:00 2001 From: Wolfgang Bangerth Date: Tue, 8 Feb 2011 14:41:04 +0000 Subject: [PATCH] Use a slightly different preconditioner: the diagonal part of the eps:eps form, instead of the (already diagonal) grad:grad form. This doesn't make a difference here, but cut the number of inner iterations in Jennifer's (nonlinear) testcase using BiCGStab by approximately 1/4 (e.g. from 40 to 30). git-svn-id: https://svn.dealii.org/trunk@23303 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/examples/step-32/step-32.cc | 67 +++++++++++++++-------------- 1 file changed, 35 insertions(+), 32 deletions(-) diff --git a/deal.II/examples/step-32/step-32.cc b/deal.II/examples/step-32/step-32.cc index 6d07266426..1e58451b5a 100644 --- a/deal.II/examples/step-32/step-32.cc +++ b/deal.II/examples/step-32/step-32.cc @@ -133,7 +133,7 @@ namespace EquationData // interpolate the following values with a physically realistic model: // const double g0 = 10.7; /* m / s^2 */ // const double g1 = 9.81; /* m / s^2 */ - + const double r = p.norm(); return -(1.245e-6 * r + 7.714e13/r/r) * p / p.norm(); } @@ -362,8 +362,8 @@ namespace Assembly FEValues stokes_fe_values; - std::vector > grad_phi_u; - std::vector phi_p; + std::vector > grads_phi_u; + std::vector phi_p; }; template @@ -374,7 +374,7 @@ namespace Assembly : stokes_fe_values (stokes_fe, stokes_quadrature, update_flags), - grad_phi_u (stokes_fe.dofs_per_cell), + grads_phi_u (stokes_fe.dofs_per_cell), phi_p (stokes_fe.dofs_per_cell) {} @@ -387,7 +387,7 @@ namespace Assembly stokes_fe_values (scratch.stokes_fe_values.get_fe(), scratch.stokes_fe_values.get_quadrature(), scratch.stokes_fe_values.get_update_flags()), - grad_phi_u (scratch.grad_phi_u), + grads_phi_u (scratch.grads_phi_u), phi_p (scratch.phi_p) {} @@ -521,7 +521,7 @@ namespace Assembly std::vector > old_strain_rates; std::vector > old_old_strain_rates; - + std::vector old_temperature_values; std::vector old_old_temperature_values; std::vector > old_temperature_grads; @@ -551,7 +551,7 @@ namespace Assembly old_old_velocity_values (quadrature.size()), old_strain_rates (quadrature.size()), old_old_strain_rates (quadrature.size()), - + old_temperature_values (quadrature.size()), old_old_temperature_values(quadrature.size()), old_temperature_grads(quadrature.size()), @@ -849,10 +849,10 @@ class BoussinesqFlowProblem struct Parameters { Parameters (); - + static void declare_parameters (ParameterHandler &prm); void parse_parameters (ParameterHandler &prm); - + double end_time; unsigned int initial_global_refinement; @@ -860,12 +860,12 @@ class BoussinesqFlowProblem bool generate_graphical_output; unsigned int graphical_output_interval; - + double stabilization_alpha; double stabilization_c_R; double stabilization_beta; }; - + ConditionalOStream pcout; Parameters parameters; @@ -975,7 +975,7 @@ BoussinesqFlowProblem::Parameters::Parameters () stabilization_c_R (0.11), stabilization_beta (0.078) {} - + template @@ -1013,7 +1013,7 @@ declare_parameters (ParameterHandler &prm) prm.declare_entry ("c_R", "0.11", Patterns::Double (0), "The c_R factor in the entropy viscosity " - "stabilization."); + "stabilization."); prm.declare_entry ("beta", "0.078", Patterns::Double (0), "The beta factor in the artificial viscosity " @@ -1036,7 +1036,7 @@ parse_parameters (ParameterHandler &prm) generate_graphical_output = prm.get_bool ("Generate graphical output"); graphical_output_interval = prm.get_integer ("Time steps between graphical output"); - + prm.enter_subsection ("Stabilization parameters"); { stabilization_alpha = prm.get_double ("alpha"); @@ -1046,7 +1046,7 @@ parse_parameters (ParameterHandler &prm) prm.leave_subsection (); } - + // @sect4{BoussinesqFlowProblem::BoussinesqFlowProblem} @@ -1303,7 +1303,7 @@ BoussinesqFlowProblem::get_extrapolated_temperature_range () const } } } - + double min_temperature, max_temperature; #ifdef DEAL_II_COMPILER_SUPPORTS_MPI MPI_Allreduce (&max_local_temperature, &max_temperature, 1, MPI_DOUBLE, @@ -1355,7 +1355,7 @@ compute_viscosity (const std::vector &old_temperature, const SymmetricTensor<2,dim> strain_rate = (old_strain_rates[q] + old_old_strain_rates[q]) / 2; - + const double T = (old_temperature[q] + old_old_temperature[q]) / 2; const double dT_dt = (old_temperature[q] - old_old_temperature[q]) / old_time_step; @@ -1369,7 +1369,7 @@ compute_viscosity (const std::vector &old_temperature, = ((EquationData::radiogenic_heating * EquationData::density(T) + 2 * EquationData::eta * strain_rate * strain_rate) / - (EquationData::density(T) * EquationData::specific_heat)); + (EquationData::density(T) * EquationData::specific_heat)); const double residual = std::abs((dT_dt + u_grad_T - kappa_Delta_T - gamma) * @@ -2014,15 +2014,18 @@ local_assemble_stokes_preconditioner (const typename DoFHandler::active_cel { for (unsigned int k=0; k global_T_range, + 2 * EquationData::eta * extrapolated_strain_rate * extrapolated_strain_rate) / (EquationData::density(old_Ts) * EquationData::specific_heat)); - + for (unsigned int i=0; i::Postprocessor : public DataPostprocessor { public: Postprocessor (const unsigned int partition); - + virtual void compute_derived_quantities_vector (const std::vector > &uh, @@ -3050,7 +3053,7 @@ compute_derived_quantities_vector (const std::vector > // pressure computed_quantities[q](dim) = uh[q](dim) * EquationData::pressure_scaling; - + // temperature computed_quantities[q](dim+1) = uh[q](dim+1); @@ -3155,7 +3158,7 @@ void BoussinesqFlowProblem::output_results () TrilinosWrappers::MPI::Vector joint_solution; joint_solution.reinit (joint_dof_handler.locally_owned_dofs(), MPI_COMM_WORLD); - + { //double minimal_pressure = stokes_solution.block(1)(0); //for (unsigned int i=0; i::run (const std::string parameter_filename) << 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); } - + prm.read_input (parameter_file); parameters.parse_parameters (prm); } - - + + GridGenerator::hyper_shell (triangulation, Point(), EquationData::R0, @@ -3550,7 +3553,7 @@ void BoussinesqFlowProblem::run (const std::string parameter_filename) { stokes_solution.sadd (1.+time_step/old_time_step, -time_step/old_time_step, old_old_stokes_solution); - temperature_solution.sadd (1.+time_step/old_time_step, + temperature_solution.sadd (1.+time_step/old_time_step, -time_step/old_time_step, old_old_temperature_solution); } @@ -3576,7 +3579,7 @@ int main (int argc, char *argv[]) parameter_filename = argv[1]; else parameter_filename = "step-32.prm"; - + const int dim = 3; BoussinesqFlowProblem flow_problem; flow_problem.run (parameter_filename); -- 2.39.5