From 54831fc156aff5c6e149a526f1d907186679214d Mon Sep 17 00:00:00 2001 From: bangerth Date: Thu, 24 Apr 2008 00:32:41 +0000 Subject: [PATCH] Add a bit to the text. git-svn-id: https://svn.dealii.org/trunk@16008 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/examples/step-31/doc/intro.dox | 125 ++++++++++--------------- 1 file changed, 48 insertions(+), 77 deletions(-) diff --git a/deal.II/examples/step-31/doc/intro.dox b/deal.II/examples/step-31/doc/intro.dox index dbb770c0de..733c64d08f 100644 --- a/deal.II/examples/step-31/doc/intro.dox +++ b/deal.II/examples/step-31/doc/intro.dox @@ -25,12 +25,18 @@ Boussinesq equations that read as follows: - \nabla \cdot \kappa \nabla T &=& \gamma. @f} +These equations fall into the class of vector-valued problems (a +toplevel overview of this topic can be found in the @ref vector_valued module). Here, ${\mathbf u}$ is the velocity field, $p$ the pressure, and $T$ the temperature of the fluid. $\varepsilon ({\mathbf u}) = \frac 12 [(\nabla{\mathbf u}) + (\nabla {\mathbf u})^T]$ is the symmetric gradient of the velocity. As can be seen, velocity and pressure solve a Stokes equation describing the motion of an incompressible -fluid. The forcing term of the fluid motion is the buoyancy of the +fluid, an equation we have previously considered in @ref step_22 "step-22"; we +will draw extensively on the experience we have gained in that program, in +particular with regard to efficient linear Stokes solvers. + +The forcing term of the fluid motion is the buoyancy of the fluid, expressed as the product of the Rayleigh number $\mathrm{Ra}$, the temperature $T$ and the gravity vector ${\mathbf g}$. (A possibly more intuitive formulation would use $\mathrm{Ra} \; (T-\bar T) @@ -47,15 +53,22 @@ temperature differences by moving around, the third equation states how the fluid motion affects the temperature field: it is an advection diffusion equation, i.e. the temperature is attached to the fluid particles and advected along in the flow field, with an additional -diffusion (heat conduction) term. - -In these equations, $\eta$ and $\kappa$ denote the viscosity and -diffusivity coefficients. In the more general case and in many -physical applications, they may and $\eta$ often will depend on the -temperature, but we will neglect this dependence for the purpose of -this tutorial program. The term $\gamma$ on the right hand side -denotes the heat sources and may be a spatially and temporally varying -function. +diffusion (heat conduction) term. In many applications, the diffusion +coefficient is fairly small, and the temperature equation is in fact +transport, not diffusion dominated and therefore in character more hyperbolic +than elliptic; we will have to take this into account when developing a stable +discretization. + +In the equations above, the term $\gamma$ on the right hand side denotes the +heat sources and may be a spatially and temporally varying function. $\eta$ +and $\kappa$ denote the viscosity and diffusivity coefficients. In the more +general case, they may and $\eta$ often will depend on the temperature, but we +will neglect this dependence for the purpose of this tutorial program even +though it is an important factor in physical applications: Most materials +become more fluid as the get hotter (i.e. $\eta$ decreases with $T$); +sometimes, as in the case of rock minerals at temperatures close to their +melting point, $\eta$ may change by orders of magnitude over the typical range +of temperatures. $\mathrm{Ra}$, called the Rayleigh @@ -87,26 +100,24 @@ velocities in absolute terms. For example, the Rayleigh number in the earth mantle has a Rayleigh number larger than $10^6$. Yet the velocities are small: the material is in fact solid rock but it is so hot and under pressure that it can flow very slowly, on the order of -at most a few centimeters per year. Nevertheless, the can lead to +at most a few centimeters per year. Nevertheless, this can lead to mixing over time scales of many million years, a time scale much shorter than for the same amount of heat to be distributed by thermal -conductivity. +conductivity and a time scale of relevance to affect the evolution of the +earth's interior and surface structure. -The equations covered here fall into the class of vector-valued problems. A -toplevel overview of this topic can be found in the @ref vector_valued module.

%Boundary and initial conditions

-Since the Boussinesq equations are derived under the assumption that -inertia of the fluid's motion does not play a role, the flow field is -at each time entirely determined by buoyancy difference at that time, -not by the flow field at previous times. This is reflected by the fact -that the first two equations above are the steady state Stokes -equation. Consequently, we do not need initial conditions for either -velocities or pressure. On the other hand, the temperature field does -satisfy an equation with a time derivative, so we need initial -conditions for $T$. +Since the Boussinesq equations are derived under the assumption that inertia +of the fluid's motion does not play a role, the flow field is at each time +entirely determined by buoyancy difference at that time, not by the flow field +at previous times. This is reflected by the fact that the first two equations +above are the steady state Stokes equation that do not contain a time +derivative. Consequently, we do not need initial conditions for either +velocities or pressure. On the other hand, the temperature field does satisfy +an equation with a time derivative, so we need initial conditions for $T$. As for boundary conditions: if $\kappa>0$ then the temperature satisfies a second order differential equation that requires @@ -138,16 +149,13 @@ where now we have a Stokes system -\nabla \cdot \eta \varepsilon ({\mathbf u}) + \nabla p &=& f, \\ \nabla\cdot \mathbf u &=& 0, @f} -where $\nabla \cdot \varepsilon (\cdot)$ is an operator similar to the +where $\nabla \cdot \eta \varepsilon (\cdot)$ is an operator similar to the Laplacian $\Delta$ applied to a vector field. Given the similarity to what we have done in @ref step_21 "step-21", it may not come as a surprise that we choose a similar approach, although we will have to make adjustments for the change in operator -in the top-left corner of the differential operator. As it will turn -out, it is actually a lot easier to solve the Stokes problem than it -is to solve the mixed Laplacian, and we will make full use of this -fact! +in the top-left corner of the differential operator.

Time stepping

@@ -187,7 +195,9 @@ field. Like solving the mixed Laplace equations, solving the Stokes equations requires us to choose particular pairs of finite elements for -velocities and pressure variables. In the current program, we use the +velocities and pressure variables. Because this has already been discussed in +@ref step_22 "step-22", we only cover this topic briefly: +Here, we use the stable pair $Q_{p+1}^d \times Q_p, p\ge 1$. These are continuous elements, so we can form the weak form of the Stokes equation without problem by integrating by parts and substituting continuous functions @@ -262,53 +272,14 @@ operator splitting where we first solve the Stokes system for the velocities and pressures using the old temperature field, and then solve for the new temperature field using the just computed velocity field. +The methods to do so are not new: Solving the linear equations coming from the +Stokes system has been discussed in great detail in @ref step_22 "step-22". In +particular, in the results section of that program, we have discussed a number +of alternative linear solver strategies. We will use the best alternative +identified there, using a GMRES solver preconditioned by a block matrix +involving the Schur complement. +XXXXXX more specifically (state matrix, preconditioner) XXXXXXXX - - - -@code - { - Vector xx(dof_handler.n_dofs()); - std::vector p_component (dof_handler.n_dofs()); - std::vector p_boundary (dof_handler.n_dofs()); - std::vector component_mask (dim+2, false); - component_mask[dim] = true; - DoFTools::extract_dofs (dof_handler, component_mask, p_component); - DoFTools::extract_boundary_dofs (dof_handler, component_mask, p_boundary); - - QGauss quadrature_formula(degree+2); - FEValues fe_values (fe, quadrature_formula, - update_values | update_gradients | update_JxW_values); - std::vector > > solution_grads(fe_values.n_quadrature_points, - std::vector > (dim+2)); - std::vector local_dof_indices (fe.dofs_per_cell); - - typename DoFHandler::active_cell_iterator - cell = dof_handler.begin_active(), - endc = dof_handler.end(); - for (; cell!=endc; ++cell) - { - fe_values.reinit (cell); - fe_values.get_function_gradients (solution, solution_grads); - - cell->get_dof_indices (local_dof_indices); - - for (unsigned int q=0; q