From: bangerth Date: Wed, 13 Aug 2008 21:27:16 +0000 (+0000) Subject: Use a Q(degree) element for the temperature as well. This allows us to get rid of... X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=47c8a3d41e02fc59c028ad65d31b2d81bfcbed45;p=dealii-svn.git Use a Q(degree) element for the temperature as well. This allows us to get rid of all the face terms in computing the right hand side of the temperature equation. On the other hand, in its current (unstabilized) version, the scheme is not stable for the transport equation, so we need to add a minimal amount of diffusion and ensure that the time step is small enough. The former restriction will be lifted once we get to proper nonlinear stabilization schemes. git-svn-id: https://svn.dealii.org/trunk@16538 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/examples/step-31/step-31.cc b/deal.II/examples/step-31/step-31.cc index a27a3446e6..26e07b896c 100644 --- a/deal.II/examples/step-31/step-31.cc +++ b/deal.II/examples/step-31/step-31.cc @@ -46,7 +46,6 @@ #include #include -#include #include #include #include @@ -874,7 +873,7 @@ BoussinesqFlowProblem::BoussinesqFlowProblem (const unsigned int degree) degree (degree), fe (FE_Q(degree+1), dim, FE_Q(degree), 1, - FE_DGQ(degree-1), 1), + FE_Q(degree), 1), dof_handler (triangulation), time_step (0), rebuild_matrices (true), @@ -1438,7 +1437,7 @@ void BoussinesqFlowProblem::assemble_system () const double Rayleigh_number = 10; std::vector > phi_u (dofs_per_cell); - std::vector > phi_grads_u (dofs_per_cell); + std::vector > grads_phi_u (dofs_per_cell); std::vector div_phi_u (dofs_per_cell); std::vector phi_p (dofs_per_cell); std::vector phi_T (dofs_per_cell); @@ -1506,7 +1505,7 @@ void BoussinesqFlowProblem::assemble_system () phi_u[k] = fe_values[velocities].value (k,q); if (rebuild_matrices) { - phi_grads_u[k] = fe_values[velocities].symmetric_gradient(k,q); + grads_phi_u[k] = fe_values[velocities].symmetric_gradient(k,q); div_phi_u[k] = fe_values[velocities].divergence (k, q); phi_p[k] = fe_values[pressure].value (k, q); phi_T[k] = fe_values[temperature].value (k, q); @@ -1516,15 +1515,29 @@ void BoussinesqFlowProblem::assemble_system () for (unsigned int i=0; i phi_i_u = fe_values[velocities].value (i, q); + // define viscosity and + // diffusion. for the + // latter, take the + // maximum of what we + // really want and the + // minimal amount of + // diffusion + // (determined + // impirically) to keep + // the scheme stable + const double eta = 1, + kappa = std::max (5e-4 * cell->diameter(), + 1e-6); + if (rebuild_matrices) for (unsigned int j=0; j gravity = ( (dim == 2) ? (Point (0,1)) : @@ -1857,22 +1870,10 @@ void BoussinesqFlowProblem::assemble_rhs_T () fe_values.JxW(q); } - -//TODO: unify the code that actually does the assembly down below for (unsigned int face_no=0; face_no::faces_per_cell; ++face_no) - if (cell->at_boundary(face_no) - || - ((cell->neighbor(face_no)->has_children() == false) - && - (cell->neighbor(face_no)->level() == cell->level()))) + if (cell->at_boundary(face_no)) { - // cell either at - // boundary or with a - // neighbor that has the - // same refinement level - // and is not further - // refined fe_face_values.reinit (cell, face_no); fe_face_values.get_function_values (old_solution, @@ -1880,25 +1881,9 @@ void BoussinesqFlowProblem::assemble_rhs_T () fe_face_values.get_function_values (solution, present_solution_values_face); - if (cell->at_boundary(face_no)) - temperature_boundary_values - .value_list (fe_face_values.get_quadrature_points(), - neighbor_temperature); - else - { - const typename DoFHandler::active_cell_iterator - neighbor = cell->neighbor(face_no); - - fe_face_values_neighbor.reinit (neighbor, - cell->neighbor_of_neighbor(face_no)); - - fe_face_values_neighbor - .get_function_values (old_solution, - old_solution_values_face_neighbor); - - for (unsigned int q=0; q::assemble_rhs_T () fe_face_values.JxW(q); } } - else - if (cell->neighbor(face_no)->has_children()) - { - // neighbor is further - // refined. loop over - // all sub faces - for (unsigned int subface_no=0; - subface_no::max_children_per_face; - ++subface_no) - { - fe_subface_values.reinit (cell, face_no, subface_no); - - fe_subface_values.get_function_values (old_solution, - old_solution_values_face); - fe_subface_values.get_function_values (solution, - present_solution_values_face); - - const typename DoFHandler::active_cell_iterator - neighbor = cell->neighbor_child_on_subface (face_no, subface_no); - - fe_face_values_neighbor.reinit (neighbor, - cell->neighbor_of_neighbor(face_no)); - - fe_face_values_neighbor - .get_function_values (old_solution, - old_solution_values_face_neighbor); - - for (unsigned int q=0; q present_u_face; - for (unsigned int d=0; d= 0); - - for (unsigned int i=0; i::active_cell_iterator - neighbor = cell->neighbor (face_no); - - const std::pair faceno_subfaceno= - cell->neighbor_of_coarser_neighbor(face_no); - const unsigned int neighbor_face_no = faceno_subfaceno.first, - neighbor_subface_no = faceno_subfaceno.second; - - Assert (neighbor->neighbor_child_on_subface (neighbor_face_no, - neighbor_subface_no) - == cell, - ExcInternalError()); - - fe_subface_values_neighbor.reinit (neighbor, - neighbor_face_no, - neighbor_subface_no); - - fe_subface_values_neighbor - .get_function_values (old_solution, - old_solution_values_face_neighbor); - - for (unsigned int q=0; q present_u_face; - for (unsigned int d=0; d= 0); - - for (unsigned int i=0; iget_dof_indices (local_dof_indices); for (unsigned int i=0; i::solve () } // for DGQ1 needs to be /15 time_step = GridTools::minimal_cell_diameter(triangulation) / - std::max (get_maximal_velocity(), .05) / 2; + std::max (get_maximal_velocity(), .05) / 4; assemble_rhs_T (); { @@ -2125,8 +1997,8 @@ void BoussinesqFlowProblem::solve () SolverControl solver_control (system_matrix.block(2,2).m(), 1e-8*system_rhs.block(2).l2_norm()); SolverCG<> cg (solver_control); - PreconditionJacobi<> preconditioner; - preconditioner.initialize (system_matrix.block(2,2)); + PreconditionSSOR<> preconditioner; + preconditioner.initialize (system_matrix.block(2,2), 1.2); try { @@ -2291,7 +2163,7 @@ void BoussinesqFlowProblem::run () GridGenerator::hyper_cube (triangulation); - triangulation.refine_global (6); + triangulation.refine_global (6); break; }