$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
void assemble_temperature_system ();
void assemble_temperature_matrix ();
double get_maximal_velocity () const;
- double get_maximal_temperature () const;
+ std::pair<double,double> get_extrapolated_temperature_range () const;
void solve ();
void output_results () const;
void refine_mesh (const unsigned int max_grid_level);
const std::vector<double> &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
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 *
std::vector<Tensor<1,dim> > 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<double,double>
+ global_T_range = get_extrapolated_temperature_range();
const double global_Omega_diameter = GridTools::diameter (triangulation);
// Now, let's start the loop
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);
- // @sect4{BoussinesqFlowProblem::get_maximal_velocity}
+ // @sect4{BoussinesqFlowProblem::get_extrapolated_temperature_range}
template <int dim>
-double BoussinesqFlowProblem<dim>::get_maximal_temperature () const
+std::pair<double,double>
+BoussinesqFlowProblem<dim>::get_extrapolated_temperature_range () const
{
QGauss<dim> quadrature_formula(temperature_degree+2);
const unsigned int n_q_points = quadrature_formula.size();
std::vector<double> old_temperature_values(n_q_points);
std::vector<double> 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<dim>::active_cell_iterator
cell = temperature_dof_handler.begin_active(),
for (unsigned int q=0; q<n_q_points; ++q)
{
- double temperature =
+ const double temperature =
(1. + time_step/old_time_step) * old_temperature_values[q]-
time_step/old_time_step * old_old_temperature_values[q];
- max_temperature = std::max (max_temperature,
- temperature);
+ min_temperature = std::min (min_temperature, temperature);
+ max_temperature = std::max (max_temperature, temperature);
}
}
- return max_temperature;
+ return std::make_pair(min_temperature, max_temperature);
}