From 3d153ef88281d711716eb217e71f3c48f56ec129 Mon Sep 17 00:00:00 2001 From: kronbichler Date: Thu, 14 Aug 2008 17:58:01 +0000 Subject: [PATCH] Added calculation of artificial viscosity. git-svn-id: https://svn.dealii.org/trunk@16553 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/examples/step-31/step-31.cc | 143 ++++++++++++++++++++++++++-- 1 file changed, 133 insertions(+), 10 deletions(-) diff --git a/deal.II/examples/step-31/step-31.cc b/deal.II/examples/step-31/step-31.cc index 6f60dfec62..ff1431057b 100644 --- a/deal.II/examples/step-31/step-31.cc +++ b/deal.II/examples/step-31/step-31.cc @@ -62,7 +62,6 @@ // This is Trilinos #include #include -#include #include #include #include @@ -394,6 +393,7 @@ class BoussinesqFlowProblem void assemble_system (); void assemble_rhs_T (); double get_maximal_velocity () const; + double get_maximal_temperature () const; void solve (); void output_results () const; void refine_mesh (); @@ -1736,6 +1736,68 @@ void BoussinesqFlowProblem::assemble_system () +template +double compute_viscosity( + const std::vector > present_solution, + const std::vector > old_solution, + const std::vector > old_old_solution, + const std::vector > > old_solution_grads, + const std::vector > > old_old_solution_grads, + const std::vector > > old_solution_hessians, + const std::vector > > old_old_solution_hessians, + const std::vector gamma_values, + const double kappa, + const double beta, + const double global_u_infty, + const double global_T_infty, + const double global_Omega_diameter, + const double cell_diameter, + const double alpha, + const double old_time_step + ) +{ + unsigned int n_q_points = old_solution.size(); + + // Stage 1: calculate residual + std::vector residual (n_q_points); + std::vector velocity_norms (n_q_points); + + for (unsigned int q=0; q < n_q_points; ++q) + { + const double dT_dt = (old_solution[q](dim+1) - old_old_solution[q](dim+1)) + / old_time_step; + double u_grad_T = 0.; + for (unsigned int d=0; d::assemble_rhs_T () QGauss face_quadrature_formula(degree+2); FEValues fe_values (fe, quadrature_formula, update_values | update_gradients | + update_hessians | update_quadrature_points | update_JxW_values); const unsigned int dofs_per_cell = fe.dofs_per_cell; const unsigned int n_q_points = quadrature_formula.size(); - const unsigned int n_face_q_points = face_quadrature_formula.size(); Vector local_rhs (dofs_per_cell); FullMatrix local_matrix (dofs_per_cell, dofs_per_cell); @@ -1806,19 +1868,31 @@ void BoussinesqFlowProblem::assemble_rhs_T () std::vector > > old_old_solution_grads(n_q_points, std::vector >(dim+2)); + std::vector > > old_solution_hessians( + n_q_points, + std::vector >(dim+2)); + std::vector > > old_old_solution_hessians( + n_q_points, + std::vector >(dim+2)); TemperatureBoundaryValues temperature_boundary_values; + RightHandSide right_hand_side; + std::vector gamma_values (n_q_points); const FEValuesExtractors::Scalar temperature (dim+1); std::vector phi_T (dofs_per_cell); std::vector > grad_phi_T (dofs_per_cell); + + const double global_u_infty = get_maximal_velocity(); + const double global_T_infty = get_maximal_temperature(); + const double global_Omega_diameter = 2.; // to be modified later. // Now, let's start the loop // over all cells in the // triangulation. The first // actions within the loop - // are, as usual, the evaluation + // are, 0as usual, the evaluation // of the FE basis functions // and the old and present // solution at the quadrature @@ -1839,6 +1913,11 @@ void BoussinesqFlowProblem::assemble_rhs_T () fe_values.get_function_gradients (old_solution, old_solution_grads); fe_values.get_function_gradients (old_old_solution, old_old_solution_grads); + fe_values.get_function_hessians (old_solution, old_solution_hessians); + fe_values.get_function_hessians (old_old_solution, old_old_solution_hessians); + + right_hand_side.value_list (fe_values.get_quadrature_points(), + gamma_values, dim+1); // build matrix contributions local_matrix = 0; @@ -1852,8 +1931,13 @@ void BoussinesqFlowProblem::assemble_rhs_T () const double kappa = std::max (5e-3 * cell->diameter(), 1e-5); - const double artificial_diffusion = 0; - + const double artificial_diffusion = + compute_viscosity (present_solution_values, old_solution_values, + old_old_solution_values, old_solution_grads, old_old_solution_grads, + old_solution_hessians, old_old_solution_hessians, gamma_values, + kappa, /* beta = */ 1., global_u_infty, global_T_infty, + global_Omega_diameter, cell->diameter(), /* alpha = */ 1., + old_time_step); for (unsigned int q=0; q::assemble_rhs_T () phi_T[k] = fe_values[temperature].value (k, q); } - const Point p = fe_values.quadrature_point(q); - const double gamma = RightHandSide().value (p, dim+1); - const double old_T = old_solution_values[q](dim+1); const double old_old_T = old_old_solution_values[q](dim+1); @@ -1912,7 +1993,7 @@ void BoussinesqFlowProblem::assemble_rhs_T () grad_phi_T[i] + time_step * - gamma * phi_T[i]) + gamma_values[q] * phi_T[i]) * fe_values.JxW(q); } @@ -1937,7 +2018,7 @@ void BoussinesqFlowProblem::assemble_rhs_T () old_grad_T * grad_phi_T[i] + time_step * - gamma * phi_T[i]) + gamma_values[q] * phi_T[i]) * fe_values.JxW(q); } @@ -2189,6 +2270,48 @@ double BoussinesqFlowProblem::get_maximal_velocity () const + + // @sect4{BoussinesqFlowProblem::get_maximal_velocity} +template +double BoussinesqFlowProblem::get_maximal_temperature () const +{ + QGauss quadrature_formula(degree+2); + const unsigned int n_q_points + = quadrature_formula.size(); + + FEValues fe_values (fe, quadrature_formula, + update_values); + std::vector > old_solution_values(n_q_points, + Vector(dim+2)); + std::vector > old_old_solution_values(n_q_points, + Vector(dim+2)); + double max_temperature = 0; + + typename DoFHandler::active_cell_iterator + cell = dof_handler.begin_active(), + endc = dof_handler.end(); + for (; cell!=endc; ++cell) + { + fe_values.reinit (cell); + fe_values.get_function_values (old_solution, old_solution_values); + fe_values.get_function_values (old_old_solution, old_old_solution_values); + + for (unsigned int q=0; q void BoussinesqFlowProblem::run () -- 2.39.5