From: kronbichler Date: Wed, 31 Aug 2011 13:32:11 +0000 (+0000) Subject: Fix stabilization for alpha=2: need to weight by entropy then. Also, the entropy... X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=7f936abc03b7eb830402db2bfe0413d1b8aa72b7;p=dealii-svn.git Fix stabilization for alpha=2: need to weight by entropy then. Also, the entropy needs to be computed against the average temperature. git-svn-id: https://svn.dealii.org/trunk@24229 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/examples/step-32/step-32.cc b/deal.II/examples/step-32/step-32.cc index f525d4737a..c96c0db93b 100644 --- a/deal.II/examples/step-32/step-32.cc +++ b/deal.II/examples/step-32/step-32.cc @@ -899,6 +899,7 @@ class BoussinesqFlowProblem void project_temperature_field (); double get_maximal_velocity () const; double get_cfl_number () const; + double get_entropy_variation (const double average_temperature) const; std::pair get_extrapolated_temperature_range () const; void solve (); void output_results (); @@ -917,6 +918,8 @@ class BoussinesqFlowProblem const std::vector > &old_old_strain_rates, const double global_u_infty, const double global_T_variation, + const double average_temperature, + const double global_entropy_variatiion, const double cell_diameter) const; public: @@ -1033,6 +1036,7 @@ class BoussinesqFlowProblem void local_assemble_temperature_rhs (const std::pair global_T_range, const double global_max_velocity, + const double global_entropy_variation, const typename DoFHandler::active_cell_iterator &cell, Assembly::Scratch::TemperatureRHS &scratch, Assembly::CopyData::TemperatureRHS &data); @@ -1426,6 +1430,87 @@ double BoussinesqFlowProblem::get_cfl_number () const +template +double +BoussinesqFlowProblem::get_entropy_variation (const double average_temperature) const +{ + // only do this if we really need entropy + // variation + if (parameters.stabilization_alpha != 2) + return 1.; + + // record maximal entropy on Gauss quadrature + // points + const QGauss quadrature_formula (parameters.temperature_degree+1); + const unsigned int n_q_points = quadrature_formula.size(); + + FEValues fe_values (temperature_fe, quadrature_formula, + update_values | update_JxW_values); + std::vector old_temperature_values(n_q_points); + std::vector old_old_temperature_values(n_q_points); + + double min_entropy = std::numeric_limits::max(), + max_entropy = -std::numeric_limits::max(), + area = 0, + entropy_integrated = 0; + + typename DoFHandler::active_cell_iterator + cell = temperature_dof_handler.begin_active(), + endc = temperature_dof_handler.end(); + for (; cell!=endc; ++cell) + if (cell->subdomain_id() == + Utilities::System::get_this_mpi_process(MPI_COMM_WORLD)) + { + fe_values.reinit (cell); + fe_values.get_function_values (old_temperature_solution, + old_temperature_values); + fe_values.get_function_values (old_old_temperature_solution, + old_old_temperature_values); + for (unsigned int q=0; q &old_temperature, const std::vector > &old_old_strain_rates, const double global_u_infty, const double global_T_variation, + const double average_temperature, + const double global_entropy_variation, const double cell_diameter) const { if (global_u_infty == 0) @@ -1577,30 +1664,36 @@ compute_viscosity (const std::vector &old_temperature, 2 * EquationData::eta * strain_rate * strain_rate) / (EquationData::density(T) * EquationData::specific_heat)); - const double residual - = std::abs((dT_dt + u_grad_T - kappa_Delta_T - gamma) * - std::pow((old_temperature[q]+old_old_temperature[q]) / 2, - parameters.stabilization_alpha-1.)); + double residual + = std::abs(dT_dt + u_grad_T - kappa_Delta_T - gamma); + if (parameters.stabilization_alpha == 2) + residual *= std::abs(T - average_temperature); max_residual = std::max (residual, max_residual); max_velocity = std::max (std::sqrt (u*u), max_velocity); } + const double max_viscosity = (parameters.stabilization_beta * + max_velocity * cell_diameter); if (timestep_number == 0) - return parameters.stabilization_beta * max_velocity * cell_diameter; + return max_viscosity; else { Assert (old_time_step > 0, ExcInternalError()); - const double global_scaling = parameters.stabilization_c_R * - global_u_infty * global_T_variation * - std::pow(global_Omega_diameter, parameters.stabilization_alpha - 2.); + double entropy_viscosity; + if (parameters.stabilization_alpha == 2) + entropy_viscosity = (parameters.stabilization_c_R * + cell_diameter * cell_diameter * + max_residual / + global_entropy_variation); + else + entropy_viscosity = (parameters.stabilization_c_R * + cell_diameter * global_Omega_diameter * + max_velocity * max_residual / + (global_u_infty * global_T_variation)); - return (parameters.stabilization_beta * - max_velocity * - std::min (cell_diameter, - std::pow(cell_diameter,parameters.stabilization_alpha) * max_residual / - global_scaling)); + return std::min (max_viscosity, entropy_viscosity); } } @@ -2399,7 +2492,7 @@ BoussinesqFlowProblem::build_stokes_preconditioner () Amg_data.elliptic = true; Amg_data.higher_order_elements = true; Amg_data.smoother_sweeps = 2; -// Amg_data.aggregation_threshold = 0.02; + Amg_data.aggregation_threshold = 0.02; Mp_preconditioner->initialize (stokes_preconditioner_matrix.block(1,1)); Amg_preconditioner->initialize (stokes_preconditioner_matrix.block(0,0), @@ -2725,6 +2818,7 @@ template void BoussinesqFlowProblem:: local_assemble_temperature_rhs (const std::pair global_T_range, const double global_max_velocity, + const double global_entropy_variation, const typename DoFHandler::active_cell_iterator &cell, Assembly::Scratch::TemperatureRHS &scratch, Assembly::CopyData::TemperatureRHS &data) @@ -2786,6 +2880,8 @@ local_assemble_temperature_rhs (const std::pair global_T_range, scratch.old_old_strain_rates, global_max_velocity, global_T_range.second - global_T_range.first, + 0.5 * (global_T_range.second + global_T_range.first), + global_entropy_variation, cell->diameter()); for (unsigned int q=0; q::assemble_temperature_system (const double maxim const std::pair global_T_range = get_extrapolated_temperature_range(); + // use midpoint between maximum and minimum + // temperature for definition of average + // temperature in entropy viscosity. Could + // also use the integral average, but the + // results are not very sensitive to this + // choice. + const double average_temperature = 0.5 * (global_T_range.first + + global_T_range.second); + const double global_entropy_variation = + get_entropy_variation (average_temperature); + typedef FilteredIterator::active_cell_iterator> SubdomainFilter; @@ -2960,6 +3067,7 @@ void BoussinesqFlowProblem::assemble_temperature_system (const double maxim this, global_T_range, maximal_velocity, + global_entropy_variation, std_cxx1x::_1, std_cxx1x::_2, std_cxx1x::_3), @@ -3056,7 +3164,7 @@ void BoussinesqFlowProblem::solve () // succeeds in 30 steps or less. unsigned int n_iterations = 0; double reduction = 0; - const double solver_tolerance = 1e-8 * stokes_rhs.l2_norm(); + const double solver_tolerance = 1e-7 * stokes_rhs.l2_norm(); SolverControl solver_control (30, solver_tolerance); try diff --git a/deal.II/examples/step-32/step-32.prm b/deal.II/examples/step-32/step-32.prm index 40b4d05057..b89786da2c 100644 --- a/deal.II/examples/step-32/step-32.prm +++ b/deal.II/examples/step-32/step-32.prm @@ -48,7 +48,7 @@ subsection Stabilization parameters set beta = 0.078 # The c_R factor in the entropy viscosity stabilization. - set c_R = 0.11 + set c_R = 0.5 end