From: bangerth Date: Wed, 20 Aug 2008 11:48:58 +0000 (+0000) Subject: Use dynamic range of temperature to define viscosity, not maximum norm. X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=5aa2689e6eb87868ccb13ef8267fe1820fa9fd2c;p=dealii-svn.git Use dynamic range of temperature to define viscosity, not maximum norm. git-svn-id: https://svn.dealii.org/trunk@16603 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/examples/step-31/doc/intro.dox b/deal.II/examples/step-31/doc/intro.dox index 9ed7afb2bf..0e5d079072 100644 --- a/deal.II/examples/step-31/doc/intro.dox +++ b/deal.II/examples/step-31/doc/intro.dox @@ -438,8 +438,11 @@ reveals that it is unitless and therefore independent of scaling) and $c(\mathbf{u},T)$ is a normalization constant that must have units $\frac{m^{\alpha-1}K^\alpha}{s}$. We will choose it as $c(\mathbf{u},T) = - \|\mathbf{u}\|_{L^\infty(\Omega)} \|T\|_{L^\infty(\Omega)} - |\textrm{diam}(\Omega)|^{\alpha-2}$. + \|\mathbf{u}\|_{L^\infty(\Omega)} \ \mathrm{var}(T) + \ |\mathrm{diam}(\Omega)|^{\alpha-2}$, +where $\mathrm{var}(T)=\max_\Omega T - \min_\Omega T$ is the range of present +temperature values (remember that buoyancy is driven by temperature +variations, not the absolute temperature). To understand why this method works consider this: If on a particular cell $K$ the temperature field is smooth, then we expect the residual to be small there (in fact to be on the order of ${\cal O}(h_K)$) and diff --git a/deal.II/examples/step-31/step-31.cc b/deal.II/examples/step-31/step-31.cc index bd8e6d2883..df5e563308 100644 --- a/deal.II/examples/step-31/step-31.cc +++ b/deal.II/examples/step-31/step-31.cc @@ -325,7 +325,7 @@ class BoussinesqFlowProblem void assemble_temperature_system (); void assemble_temperature_matrix (); double get_maximal_velocity () const; - double get_maximal_temperature () const; + std::pair get_extrapolated_temperature_range () const; void solve (); void output_results () const; void refine_mesh (const unsigned int max_grid_level); @@ -1699,7 +1699,7 @@ double compute_viscosity( const std::vector &gamma_values, const double kappa, const double global_u_infty, - const double global_T_infty, + const double global_T_variation, const double global_Omega_diameter, const double cell_diameter, const double old_time_step @@ -1741,7 +1741,7 @@ double compute_viscosity( max_velocity = std::max (std::sqrt (u*u), max_velocity); } - const double global_scaling = global_u_infty * global_T_infty / + const double global_scaling = global_u_infty * global_T_variation / std::pow(global_Omega_diameter, alpha - 2.); return (beta * @@ -1932,7 +1932,8 @@ void BoussinesqFlowProblem::assemble_temperature_system () 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 std::pair + global_T_range = get_extrapolated_temperature_range(); const double global_Omega_diameter = GridTools::diameter (triangulation); // Now, let's start the loop @@ -1989,7 +1990,8 @@ void BoussinesqFlowProblem::assemble_temperature_system () old_old_temperature_hessians, present_stokes_values, gamma_values, - kappa, global_u_infty, global_T_infty, + kappa, global_u_infty, + global_T_range.second - global_T_range.first, global_Omega_diameter, cell->diameter(), old_time_step); @@ -2339,9 +2341,10 @@ double BoussinesqFlowProblem::get_maximal_velocity () const - // @sect4{BoussinesqFlowProblem::get_maximal_velocity} + // @sect4{BoussinesqFlowProblem::get_extrapolated_temperature_range} template -double BoussinesqFlowProblem::get_maximal_temperature () const +std::pair +BoussinesqFlowProblem::get_extrapolated_temperature_range () const { QGauss quadrature_formula(temperature_degree+2); const unsigned int n_q_points = quadrature_formula.size(); @@ -2351,7 +2354,12 @@ double BoussinesqFlowProblem::get_maximal_temperature () const std::vector old_temperature_values(n_q_points); std::vector old_old_temperature_values(n_q_points); - double max_temperature = 0; + double min_temperature = (1. + time_step/old_time_step) * + old_temperature_solution.linfty_norm() + + + time_step/old_time_step * + old_old_temperature_solution.linfty_norm(), + max_temperature = -min_temperature; typename DoFHandler::active_cell_iterator cell = temperature_dof_handler.begin_active(), @@ -2364,16 +2372,16 @@ double BoussinesqFlowProblem::get_maximal_temperature () const for (unsigned int q=0; q