From a302a022251895a59c15709ce07ec2566338e2b3 Mon Sep 17 00:00:00 2001 From: kronbichler Date: Wed, 3 Dec 2008 13:48:45 +0000 Subject: [PATCH] Simplified the temperature rhs assembly. Use IC preconditioner for temperature system instead of SSOR. Could skip two parameters in artificial viscosity call. git-svn-id: https://svn.dealii.org/trunk@17826 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/examples/step-31/step-31.cc | 148 ++++++++++++---------------- 1 file changed, 65 insertions(+), 83 deletions(-) diff --git a/deal.II/examples/step-31/step-31.cc b/deal.II/examples/step-31/step-31.cc index 7a232a5328..f9e7107a6d 100644 --- a/deal.II/examples/step-31/step-31.cc +++ b/deal.II/examples/step-31/step-31.cc @@ -592,7 +592,6 @@ class BoussinesqFlowProblem void output_results () const; void refine_mesh (const unsigned int max_grid_level); - static double compute_viscosity(const std::vector &old_temperature, const std::vector &old_old_temperature, @@ -605,12 +604,11 @@ class BoussinesqFlowProblem const std::vector &gamma_values, const double global_u_infty, const double global_T_variation, - const double global_Omega_diameter, - const double cell_diameter, - const double old_time_step); + const double cell_diameter); Triangulation triangulation; + double global_Omega_diameter; const unsigned int stokes_degree; FESystem stokes_fe; @@ -977,9 +975,7 @@ compute_viscosity (const std::vector &old_temperature, const std::vector &gamma_values, const double global_u_infty, const double global_T_variation, - const double global_Omega_diameter, - const double cell_diameter, - const double old_time_step) + const double cell_diameter) { const double beta = 0.015 * dim; const double alpha = 1; @@ -1506,19 +1502,23 @@ BoussinesqFlowProblem::build_stokes_preconditioner () // matrix we want it to apply to. // // Finally, we also initialize the - // preconditioner for the inversion of the - // pressure mass matrix. This matrix is - // symmetric and well-behaved, so we can - // chose a simple preconditioner. We stick - // with an incomple Cholesky (IC) + // preconditioner for the inversion of + // the pressure mass matrix. This matrix + // is symmetric and well-behaved, so we + // can chose a simple preconditioner. We + // stick with an incomple Cholesky (IC) // factorization preconditioner, which is - // designed for symmetric matrices. We wrap - // the preconditioners into a - // boost::shared_ptr pointer, - // which makes it easier to recreate the - // preconditioner next time around since we - // do not have to care about destroying the - // previously used object. + // designed for symmetric matrices. We + // could have also chosen an SSOR + // preconditioner with relaxation factor + // around 1.2, but IC is cheaper for our + // example. We wrap the preconditioners + // into a boost::shared_ptr + // pointer, which makes it easier to + // recreate the preconditioner next time + // around since we do not have to care + // about destroying the previously used + // object. amg_data.elliptic = true; amg_data.higher_order_elements = true; amg_data.aggregation_threshold = 5e-2; @@ -2007,7 +2007,6 @@ void BoussinesqFlowProblem::assemble_temperature_system () const double global_u_infty = get_maximal_velocity(); const std::pair global_T_range = get_extrapolated_temperature_range(); - const double global_Omega_diameter = GridTools::diameter (triangulation); // Now, let's start the loop over all // cells in the triangulation. Again, we @@ -2093,8 +2092,7 @@ void BoussinesqFlowProblem::assemble_temperature_system () gamma_values, global_u_infty, global_T_range.second - global_T_range.first, - global_Omega_diameter, cell->diameter(), - old_time_step); + cell->diameter()); for (unsigned int q=0; q::assemble_temperature_system () phi_T[k] = temperature_fe_values.shape_value (k, q); } - const double old_T = old_temperature_values[q]; - const double old_old_T = old_old_temperature_values[q]; - - const Tensor<1,dim> old_grad_T = old_temperature_grads[q]; - const Tensor<1,dim> old_old_grad_T = old_old_temperature_grads[q]; - + const double old_Ts = use_bdf2_scheme ? + (old_temperature_values[q] * + (time_step + old_time_step) / old_time_step + - + old_old_temperature_values[q] * + (time_step * time_step) / + (old_time_step * (time_step + old_time_step))) + : + old_temperature_values[q]; + + const Tensor<1,dim> ext_grad_T = use_bdf2_scheme ? + (old_temperature_grads[q] * + (1+time_step/old_time_step) + - + old_old_temperature_grads[q] * + time_step / old_time_step) + : + old_temperature_grads[q]; Tensor<1,dim> extrapolated_u; for (unsigned int d=0; d::assemble_temperature_system () extrapolated_u[d] = old_stokes_values[q](d); } - if (use_bdf2_scheme == true) - { - for (unsigned int i=0; iget_dof_indices (local_dof_indices); temperature_constraints.distribute_local_to_global (local_rhs, local_dof_indices, @@ -2294,12 +2271,15 @@ void BoussinesqFlowProblem::solve () // moderately ill conditioned even for // small mesh sizes and we get a // reasonably good preconditioner by - // simple means, for example SSOR with a - // relaxation parameter of 1.2. As a - // solver, we choose the conjugate - // gradient method CG. As before, we tell - // the solver to use Trilinos vectors via - // the template argument + // simple means, for example with an + // incomplete Cholesky decomposition + // preconditioner (IC) as we also use for + // preconditioning the pressure mass + // matrix solver. As a solver, we choose + // the conjugate gradient method CG. As + // before, we tell the solver to use + // Trilinos vectors via the template + // argument // TrilinosWrappers::Vector. // Finally, we solve, distribute the // hanging node constraints and write out @@ -2323,8 +2303,8 @@ void BoussinesqFlowProblem::solve () 1e-8*temperature_rhs.l2_norm()); SolverCG cg (solver_control); - TrilinosWrappers::PreconditionSSOR preconditioner; - preconditioner.initialize (temperature_matrix, 1.2); + TrilinosWrappers::PreconditionIC preconditioner; + preconditioner.initialize (temperature_matrix); cg.solve (temperature_matrix, temperature_solution, temperature_rhs, preconditioner); @@ -2760,6 +2740,8 @@ void BoussinesqFlowProblem::run () GridGenerator::hyper_cube (triangulation); + global_Omega_diameter = GridTools::diameter (triangulation); + triangulation.refine_global (initial_refinement); setup_dofs(); -- 2.39.5