From: bangerth Date: Thu, 14 Aug 2008 23:03:11 +0000 (+0000) Subject: Commit a version in which artificial diffusion appears to work. This version also... X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=6b92377cc60c674f1586a7d0fa96a3f75c9cb369;p=dealii-svn.git Commit a version in which artificial diffusion appears to work. This version also contains a function that writes the artificial diffusion parameter into the output file. git-svn-id: https://svn.dealii.org/trunk@16557 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 7a072b7a59..25c9c23216 100644 --- a/deal.II/examples/step-31/step-31.cc +++ b/deal.II/examples/step-31/step-31.cc @@ -325,6 +325,7 @@ class BoussinesqFlowProblem double get_maximal_velocity () const; double get_maximal_temperature () const; void solve (); + void test (Vector &) const; void output_results () const; void refine_mesh (); @@ -1706,20 +1707,24 @@ double compute_viscosity( 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(); + const double beta = 0.1; + const double alpha = 1; + + if (global_u_infty == 0) + return 5e-3 * cell_diameter; + + const 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); + double max_residual = 0; + double max_velocity = 0; for (unsigned int q=0; q < n_q_points; ++q) { @@ -1727,33 +1732,35 @@ double compute_viscosity( / old_time_step; double u_grad_T = 0.; for (unsigned int d=0; d::assemble_rhs_T () 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. + const double global_Omega_diameter = GridTools::diameter (triangulation); // Now, let's start the loop // over all cells in the @@ -1863,6 +1870,8 @@ void BoussinesqFlowProblem::assemble_rhs_T () for (; cell!=endc; ++cell) { local_rhs = 0; + local_matrix = 0; + fe_values.reinit (cell); fe_values.get_function_values (solution, present_solution_values); @@ -1879,7 +1888,6 @@ void BoussinesqFlowProblem::assemble_rhs_T () gamma_values, dim+1); // build matrix contributions - local_matrix = 0; // define diffusion. take the // maximum of what we really @@ -1887,16 +1895,14 @@ void BoussinesqFlowProblem::assemble_rhs_T () // of diffusion (determined // impirically) to keep the // scheme stable - const double kappa = std::max (5e-3 * cell->diameter(), - 1e-5); - - 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); + const double kappa = 1e-6; + const double nu = + 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, global_u_infty, global_T_infty, + global_Omega_diameter, cell->diameter(), + old_time_step); for (unsigned int q=0; q::assemble_rhs_T () phi_T[i] - time_step * - artificial_diffusion * + nu * ((1+time_step/old_time_step) * old_grad_T - time_step / old_time_step * old_old_grad_T) * @@ -1973,7 +1979,7 @@ void BoussinesqFlowProblem::assemble_rhs_T () present_u * old_grad_T * phi_T[i] - time_step * - artificial_diffusion * + nu * old_grad_T * grad_phi_T[i] + time_step * @@ -1998,6 +2004,83 @@ void BoussinesqFlowProblem::assemble_rhs_T () +template +void +BoussinesqFlowProblem::test (Vector &viscosity) const +{ + QGauss 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 n_q_points = quadrature_formula.size(); + + std::vector > present_solution_values (n_q_points, + Vector(dim+2)); + std::vector > old_solution_values(n_q_points, + Vector(dim+2)); + std::vector > old_old_solution_values(n_q_points, + Vector(dim+2)); + std::vector > > + old_solution_grads(n_q_points, + std::vector >(dim+2)); + 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)); + RightHandSide right_hand_side; + std::vector gamma_values (n_q_points); + + const FEValuesExtractors::Scalar temperature (dim+1); + + const double global_u_infty = get_maximal_velocity(); + const double global_T_infty = get_maximal_temperature(); + const double global_Omega_diameter = GridTools::diameter (triangulation); + + typename DoFHandler::active_cell_iterator + cell = dof_handler.begin_active(), + endc = dof_handler.end(); + + for (unsigned int index=0; cell!=endc; ++cell, ++index) + { + fe_values.reinit (cell); + + fe_values.get_function_values (solution, present_solution_values); + fe_values.get_function_values (old_solution, old_solution_values); + fe_values.get_function_values (old_old_solution, old_old_solution_values); + + 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); + + const double kappa = 1e-6; + right_hand_side.value_list (fe_values.get_quadrature_points(), + gamma_values, dim+1); + + viscosity(index) = + (timestep_number == 0 + ? + 5e-3 * cell->diameter() + : + 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, global_u_infty, global_T_infty, + global_Omega_diameter, cell->diameter(), + old_time_step)); + } +} + + + // @sect4{BoussinesqFlowProblem::solve} template @@ -2116,6 +2199,9 @@ void BoussinesqFlowProblem::output_results () const solution_names.push_back ("p"); solution_names.push_back ("T"); + Vector viscosity (triangulation.n_active_cells()); + test (viscosity); + DataOut data_out; data_out.attach_dof_handler (dof_handler); @@ -2131,6 +2217,8 @@ void BoussinesqFlowProblem::output_results () const DataOut::type_dof_data, data_component_interpretation); + data_out.add_data_vector (viscosity, "viscosity"); + data_out.build_patches (); std::ostringstream filename;