\nabla \cdot {\mathbf u} &=& 0.
@f}
The right hand side of the first equation appears a wee bit
-unmotivated. Here's how things should really be. On the right side, we
+unmotivated. Here's how things should really be. We
need the external forces that act on the fluid, which we assume are
given by gravity only. In the current case, we assume that the fluid
does expand slightly for the purposes of this gravity force, but not
- A final reason to do this is discussed in the results section and
concerns possible extensions to the model we use here. It has to do
- with the fact that while the temperature equation we use here does not
+ with the fact that the temperature equation (see below) we use here does not
include a term that contains the pressure. It should, however:
rock, like gas, heats up as you compress it. Consequently,
material that rises up cools adiabatically, and cold material that
sinks down heats adiabatically. We discuss this further below.
@note There is, however, a downside to this procedure. In the earth,
-the dynamic pressure is several order of magnitude smaller than the
+the dynamic pressure is several orders of magnitude smaller than the
total pressure. If we use the equations above and solve all variables
to, say, 4 digits of accuracy, then we may be able to get the velocity
and the total pressure right, but we will have no accuracy at all if
could store all elements with even indices on one processor and all odd ones
on another. That's not very efficient, of course, but it's
possible. Furthermore, the elements of these partitionings do not
- necessarily by mutually exclusive. This is important because when
+ necessarily be mutually exclusive. This is important because when
postprocessing solutions, we need access to all locally relevant or at least
the locally active degrees of freedom (see the module on @ref distributed
for a definition, as well as the discussion in step-40). Which elements the
<li>The right hand side of the temperature equation contains the rate of
%internal heating $\gamma$. The earth does heat naturally through several mechanisms:
radioactive decay, chemical separation (heavier elements sink to the bottom,
- lighter ones rise to the top; the countercurrents dissipate emergy equal to
+ lighter ones rise to the top; the countercurrents dissipate energy equal to
the loss of potential energy by this separation process); heat release
by crystallization of liquid metal as the solid inner core of the earth
grows; and heat dissipation from viscous friction as the fluid moves.
const double year_in_seconds = 60*60*24*365.2425;
-//TODO: document in intro.dox
- // scale not by R1-R0, but by a
- // typical length scale, say 10km,
- // of variation ("plume
- // diameter"). this choice also
- // roughly equilibrates the sizes
- // of the velocity and pressure
- // components of the solution
- // vectors
const double pressure_scaling = eta / 10000;
template <int dim>
Tensor<1,dim> gravity_vector (const Point<dim> &p)
{
-// interpolate the following values with a physically realistic model:
+// Interpolate the following values with a physically realistic model:
// const double g0 = 10.7; /* m / s^2 */
// const double g1 = 9.81; /* m / s^2 */
template <int dim>
double adiabatic_pressure (const Point<dim> &p)
{
- // the static, adiabatic pressure
+ // The static, adiabatic pressure
// satisfies
// dP/dr = -g rho
- // assuming a constant density,
+ // Assuming a constant density,
// we can integrate the pressure
// equation in depth to get that
// the adiabatic pressure equals
// $P(r) = rho_0 \int_r^{R_1} g(r) dr$
//
- // using the model for the
+ // Using the model for the
// gravity vector above, this
// yields the following formula:
const double r = p.norm();
template <int dim>
double adiabatic_temperature (const Point<dim> &p)
{
- // the static, adiabatic
+ // The static, adiabatic
// temperature satisfies
- // $dT/dr = -T alpha/c_P g$
+ // $dT/dr = -T \alpha/c_P g$
- // let's assume constant gravity,
+ // Let's assume constant gravity,
// then we get by integration
const double r = p.norm();
// @sect3{Linear solvers and preconditioners}
+// TODO (MK): update
+
// In comparison to step-31, we did one
// change in the linear algebra of the
// problem: We exchange the
// struct is equipped with a constructor
// that create an FEValues object for a
// @ref FiniteElement "finite element", a
-// @ref Quadrature "quadrature formula"
+// @ref Quadrature "quadrature formula", the
+// @ref Mapping "mapping" that describes the
+// interpolation of curved boundaries,
// and some
// @ref UpdateFlags "update flags".
// Moreover, we manually
// frequently have to query velocities and
// temperatures at arbitrary quadrature
// points; consequently, rather than
-// "localizing" a vector whenever we need a
-// localized vector, we solve linear
+// importing ghost information of a vector
+// whenever we need access to degrees of freedom
+// that are relevant locally but owned by
+// another processor,
+// we solve linear
// systems in %parallel but then immediately
-// localize the solution for further
+// initialize a vector including ghost entries
+// the solution for further
// processing. The various
// <code>*_solution</code> vectors are
// therefore filled immediately after
// @sect4{The BoussinesqFlowProblem helper functions}
//
-// Except two small details, this
-// function is the very same as in
+// Except two small details, the function
+// to compute the global maximum of
+// the velocity is the very same as in
// step-31. The first detail is
// actually common to all functions
// that implement loop over all cells
// each processor calculate the
// maximum among its cells, and then
// do a global communication
-// operation called
-// <code>MaxAll</code> that searches
+// operation
+// <code>Utilities::MPI::max</code> that searches
// for the maximum value among all
// the maximum values of the
// individual processors. MPI
// will do the right thing even if we
// work without MPI and on a single
// machine only. The call to
-// <code>MaxAll</code> needs three
+// <code>Utilities::MPI::max</code> needs two
// arguments, namely the local
-// maximum (input), a field for the
-// global maximum (output), and an
-// integer value one that says that
-// we only work on one double.
+// maximum (input) and the MPI communicator,
+// which is MPI_COMM_WORLD in this example.
template <int dim>
double BoussinesqFlowProblem<dim>::get_maximal_velocity () const
{
// Similar function to before, but we now
-// compute the cfl number, i.e., maximal
+// compute the CFL number, i.e., maximal
// velocity on a cell divided by the cell
-// diameter
+// diameter. This number is necessary to
+// determine the time step size, as we use
+// a semi-explicit time stepping scheme for
+// the temperature equation (see step-31 for
+// a discussion).
template <int dim>
double BoussinesqFlowProblem<dim>::get_cfl_number () const
{
+ // Next comes the computation of the global entropy variation
+ // $\|E(T)-avg(E)\|_\infty$. This is needed for the evaluation of the
+ // stabilization in the temperature equation as explained in the
+ // introduction. The entropy variation is actually only needed if we use
+ // $\alpha=2$ as a power in the residual computation. The infinity norm is
+ // computed by the maxima over quadrature points, as usual in discrete
+ // computations.
template <int dim>
double
BoussinesqFlowProblem<dim>::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<dim> quadrature_formula (parameters.temperature_degree+1);
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);
+ // Note how we preset the data fields for
+ // the minimum and maximum values: The
+ // minimum is initialized with a bigger
+ // and the maximum with a smaller number
+ // than one that is going to
+ // appear. These numbers will be
+ // overwritten in the cell loop or in the
+ // communication step at the latest.
double min_entropy = std::numeric_limits<double>::max(),
max_entropy = -std::numeric_limits<double>::max(),
area = 0,
}
}
- // do MPI data exchange: we need to sum over
+ // MPI data exchange: we need to sum over
// the two integrals (area,
// entropy_integrated), and get the extrema
// for maximum and minimum. combine
- // MPI_Allreduce for two values since that is
+ // MPI_Allreduce for two values since it is
// an expensive operation
const double local_for_sum[2] = { entropy_integrated, area },
local_for_max[2] = { -min_entropy, max_entropy };
-// Again, this is only a slightly
+// Again, the evaluation of the extrapolated
+// temperature range is only a slightly
// modified version of the respective
// function in step-31. What is new is
// that each processor works on its
std::vector<double> old_temperature_values(n_q_points);
std::vector<double> old_old_temperature_values(n_q_points);
- // This presets the minimum with a bigger
- // and the maximum with a smaller number
- // than one that is going to appear. Will
- // be overwritten in the cell loop or in
- // the communication step at the
- // latest.
double min_local_temperature = std::numeric_limits<double>::max(),
max_local_temperature = -std::numeric_limits<double>::max();
projects / dealii-svn.git / commitdiff