]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Remove all comments and other documentation copied from step-31. Only new stuff shoul...
authorbangerth <bangerth@0785d39b-7218-0410-832d-ea1e28bc413d>
Wed, 3 Sep 2008 12:58:32 +0000 (12:58 +0000)
committerbangerth <bangerth@0785d39b-7218-0410-832d-ea1e28bc413d>
Wed, 3 Sep 2008 12:58:32 +0000 (12:58 +0000)
git-svn-id: https://svn.dealii.org/trunk@16729 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/examples/step-32/doc/intro.dox
deal.II/examples/step-32/doc/results.dox
deal.II/examples/step-32/step-31.cc

index 33ce876a6fddcb7d0b50e669c9fc04a5bc9ca06e..acaa3cb3b7034d34a52acec825ad13e670152d5f 100644 (file)
@@ -14,901 +14,3 @@ California Institute of Technology.
 
 <a name="Intro"></a>
 <h1>Introduction</h1>
-
-<h3>The Boussinesq equations</h3>
-
-This program deals with an interesting physical problem: how does a
-fluid (i.e. a liquid or gas) behave if it experiences differences in
-buoyancy caused by temperature differences? It is clear that those
-parts of the fluid that are hotter (and therefore lighter) are going
-to rise up and those that are cooler (and denser) are going to sink
-down with gravity.
-
-In cases where the fluid moves slowly enough such that inertia effects
-can be neglected, the equations that describe such behavior are the
-Boussinesq equations that read as follows:
-@f{eqnarray*}
-  -\nabla \cdot \eta \varepsilon ({\mathbf u}) + \nabla p &=& 
-  \mathrm{Ra} \; T \mathbf{g},
-  \\
-  \nabla \cdot {\mathbf u} &=& 0,
-  \\
-  \frac{\partial T}{\partial t}
-  +
-  {\mathbf u} \cdot \nabla T
-  -
-  \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, <b>u</b> is the velocity field, <i>p</i> the pressure, and <i>T</i>
-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, 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 <i>T</i> and the gravity vector <b>g</b>. (A possibly
-more intuitive formulation would use $\mathrm{Ra} \; (T-\bar T)
-\mathbf{g}$ as right hand side where $\bar T$ is the average
-temperature, and the right hand side then describes the forces due to
-local deviations from the average density; this formulation is
-entirely equivalent if the gravity vector results from a gravity
-potential $\phi$, i.e. $\mathbf{g}=-\nabla\phi$, and yields the exact
-same solution except for the pressure which will now be $p+\mathrm{Ra}
-\;\bar T \phi$.)
-
-While the first two equations describe how the fluid reacts to
-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 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, which we assume
-constant for this tutorial program. The more general case when $\eta$ depends on
-the temperature is an important factor in physical applications: Most materials
-become more fluid as they get hotter (i.e., $\eta$ decreases with <i>T</i>);
-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 <a
-href="http://en.wikipedia.org/wiki/Rayleigh_number">Rayleigh
-number</a>, is a dimensionless number that describes the ratio of heat
-transport due to convection induced by buoyancy changes from
-temperature differences, and of heat transport due to thermal
-diffusion. A small Rayleigh number implies that buoyancy is not strong
-relative to viscosity and fluid motion <b>u</b> is slow enough so
-that heat diffusion $\kappa\Delta T$ is the dominant heat transport
-term. On the other hand, a fluid with a high Rayleigh number will show
-vigorous convection that dominates heat conduction. 
-
-For most fluids for which we are interested in computing thermal
-convection, the Rayleigh number is very large, often $10^6$ or
-larger. From the structure of the equations, we see that this will
-lead to large pressure differences and large velocities. Consequently,
-the convection term in the convection-diffusion equation for <i>T</i> will
-also be very large and an accurate solution of this equation will
-require us to choose small time steps. Problems with large Rayleigh
-numbers are therefore hard to solve numerically for similar reasons
-that make solving the <a
-href="http://en.wikipedia.org/wiki/Navier-stokes_equations">Navier-Stokes
-equations</a> hard to solve when the <a
-href="http://en.wikipedia.org/wiki/Reynolds_number">Reynolds number
-$\mathrm{Re}$</a> is large.
-
-Note that a large Rayleigh number does not necessarily involve large
-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, 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 and a time scale of relevance to affect the evolution of the
-earth's interior and surface structure.
-
-
-
-<h3>%Boundary and initial conditions</h3>
-
-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 <i>T</i>.
-
-As for boundary conditions: if $\kappa>0$ then the temperature
-satisfies a second order differential equation that requires
-boundary data all around the boundary for all times. These can either be a
-prescribed boundary temperature $T|_{\partial\Omega}=T_b$ (Dirichlet boundary
-conditions), or a prescribed thermal flux $\mathbf{n}\cdot\kappa\nabla
-T|_{\partial\Omega}=\phi$; in this program, we will use an insulated boundary
-condition, i.e. prescribe no thermal flux: $\phi=0$. 
-
-Similarly, the velocity field requires us to pose boundary conditions. These
-may be no-slip no-flux conditions <b>u</b>=0 on $\partial\Omega$ if the fluid
-sticks to the boundary, or no normal flux conditions $\mathbf n \cdot \mathbf
-u = 0$ if the fluid can flow along but not across the boundary, or any number
-of other conditions that are physically reasonable. In this program, we will
-use no normal flux conditions.
-
-
-<h3>Solution approach</h3>
-
-Like the equations solved in @ref step_21 "step-21", we here have a
-system of differential-algebraic equations (DAE): with respect to the time
-variable, only the temperature equation is a differential equation
-whereas the Stokes system for <b>u</b> and <i>p</i> has no
-time-derivatives and is therefore of the sort of an algebraic
-constraint that has to hold at each time instant. The main difference
-to @ref step_21 "step-21" is that the algebraic constraint there was a
-mixed Laplace system of the form
-@f{eqnarray*}
-  \mathbf u + {\mathbf K}\lambda \nabla p &=& 0, \\
-  \nabla\cdot \mathbf u &=& f,
-@f}
-where now we have a Stokes system
-@f{eqnarray*}
-  -\nabla \cdot \eta \varepsilon ({\mathbf u}) + \nabla p &=& f, \\
-  \nabla\cdot \mathbf u &=& 0,
-@f}
-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. 
-
-
-<h4>Time stepping</h4>
-
-The structure of the problem as a DAE allows us to use the same
-strategy as we have already used in @ref step_21 "step-21", i.e. we
-use a time lag scheme: first solve the Stokes equations for velocity and
-pressure using the temperature field from the previous time step, then
-with the new velocities update the temperature field for the current
-time step. In other words, in time step <i>n</i> we first solve the Stokes
-system
-@f{eqnarray*}
-  -\nabla \cdot \eta \varepsilon ({\mathbf u}^n) + \nabla p^n &=& 
-  \mathrm{Ra} \; T^{n-1} \mathbf{g},
-  \\
-  \nabla \cdot {\mathbf u}^n &=& 0,
-@f}
-and then the temperature equation with the so-computed velocity field 
-${\mathbf u}^n$. In contrast to @ref step_21 "step-21", we'll use a
-higher order time stepping scheme here, namely the <a 
-href="http://en.wikipedia.org/wiki/Backward_differentiation_formula">Backward
-Differentiation Formula scheme of order 2 (BDF-2 in short)</a> that
-replaces the time derivative $\frac{\partial T}{\partial t}$ by the (one-sided)
-difference quotient $\frac{\frac 32 T^{n}-2T^{n-1}+\frac 12 T^{n-2}}{k}$ with
-<i>k</i> the time step size.
-
-This gives the discretized-in-time temperature equation
-@f{eqnarray*}
-  \frac 32 T^n
-  -
-  k\nabla \cdot \kappa \nabla T^n 
-  &=& 
-  2 T^{n-1}
-  -
-  \frac 12 T^{n-2}
-  -
-  k{\mathbf u}^n \cdot \nabla (2T^{n-1}-T^{n-2})
-  +
-  k\gamma.
-@f}
-Note how the temperature equation is
-solved semi-explicitly: diffusion is treated implicitly whereas
-advection is treated explicitly using the just-computed velocity
-field but only previously computed temperature fields. The
-temperature terms appearing in the advection term are forward
-projected to the current time:
-$T^n \approx T^{n-1} + k_n
-\frac{\partial T}{\partial t} \approx T^{n-1} + k_n
-\frac{T^{n-1}-T^{n-2}}{k_n} = 2T^{n-1}-T^{n-2}$. We need this projection
-for maintaining the order of accuracy of the BDF-2 scheme. In other words, the
-temperature fields we use in the explicit right hand side are first
-order approximations of the current temperature field &mdash; not
-quite an explicit time stepping scheme, but by character not too far
-away either.
-
-The introduction of the temperature extrapolation limits the time step
-by a <a href="http://en.wikipedia.org/wiki/Courant–Friedrichs–Lewy_condition">
-Courant-Friedrichs-Lewy (CFL) condition</a> just like it was in 
-@ref step_21 "step-21". (We wouldn't have had that stability condition if
-we treated the advection term implicitly since the BDF-2 scheme is A-stable, 
-at the price that we needed to build a new temperature matrix at each time 
-step.) We will discuss the exact choice of time step in the <a
-href="#Results">results section</a>, but for the moment of importance is that 
-this CFL condition means that the time step
-size <i>k</i> may change from time step to time step, and that we have to
-modify the above formula slightly. If $k_n,k_{n-1}$ are the time steps
-sizes of the current and previous time step, then we use the
-approximations 
-$\frac{\partial T}{\partial t} \approx
- \frac 1{k_n}
- \left(
-       \frac{2k_n+k_{n-1}}{k_n+k_{n-1}} T^{n} 
-       -
-       \frac{k_n+k_{n-1}}{k_{n-1}}T^{n-1} 
-       +
-       \frac{k_n^2}{k_{n-1}(k_n+k_{n-1})} T^{n-2}
- \right)$
-and
-$T^n \approx 
-   T^{n-1} + k_n \frac{\partial T}{\partial t} 
-   \approx
-   T^{n-1} + k_n
-   \frac{T^{n-1}-T^{n-2}}{k_{n-1}}
-   =
-   \left(1+\frac{k_n}{k_{n-1}}\right)T^{n-1}-\frac{k_n}{k_{n-1}}T^{n-2}$,
-and above equation is generalized as follows:
-@f{eqnarray*}
-  \frac{2k_n+k_{n-1}}{k_n+k_{n-1}} T^n
-  -
-  k_n\nabla \cdot \kappa \nabla T^n 
-  &=& 
-  \frac{k_n+k_{n-1}}{k_{n-1}} T^{n-1}
-  -
-  \frac{k_n^2}{k_{n-1}(k_n+k_{n-1})} T^{n-2}
-  -
-  k_n{\mathbf u}^n \cdot \nabla \left[
-    \left(1+\frac{k_n}{k_{n-1}}\right)T^{n-1}-\frac{k_n}{k_{n-1}}T^{n-2}
-  \right]
-  +
-  k_n\gamma.
-@f}
-That's not an easy to read equation, but will provide us with the
-desired higher order accuracy. As a consistency check, it is easy to
-verify that it reduces to the same equation as above if $k_n=k_{n-1}$.
-
-As a final remark we note that the choice of a higher order time
-stepping scheme of course forces us to keep more time steps in memory;
-in particular, we here will need to have $T^{n-2}$ around, a vector
-that we could previously discard. This seems like a nuisance that we
-were able to avoid previously by using only a first order time
-stepping scheme, but as we will see below when discussing the topic of
-stabilization, we will need this vector anyway and so keeping it
-around for time discretization is essentially for free and gives us
-the opportunity to use a higher order scheme.
-
-
-<h4>Weak form and space discretization for the Stokes part</h4>
-
-Like solving the mixed Laplace equations, solving the Stokes equations
-requires us to choose particular pairs of finite elements for
-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
-by their discrete counterparts:
-@f{eqnarray*}
-  (\nabla {\mathbf v}_h, \eta \varepsilon ({\mathbf u}^n_h))
-  -
-  (\nabla \cdot {\mathbf v}_h, p^n_h) 
-  &=& 
-  ({\mathbf v}_h, \mathrm{Ra} \; T^{n-1}_h \mathbf{g}),
-  \\
-  (q_h, \nabla \cdot {\mathbf u}^n_h) &=& 0,
-@f}
-for all test functions $\mathbf v_h, q_h$. The first term of the first
-equation is considered as the inner product between tensors, i.e.
-$(\nabla {\mathbf v}_h, \eta \varepsilon ({\mathbf u}^n_h))_\Omega
- = \int_\Omega \sum_{i,j=1}^d [\nabla {\mathbf v}_h]_{ij}
-           \eta [\varepsilon ({\mathbf u}^n_h)]_{ij}\, dx$.
-Because the second tensor in this product is symmetric, the
-anti-symmetric component of $\nabla {\mathbf v}_h$ plays no role and
-it leads to the entirely same form if we use the symmetric gradient of
-$\mathbf v_h$ instead. Consequently, the formulation we consider and
-that we implement is
-@f{eqnarray*}
-  (\varepsilon({\mathbf v}_h), \eta \varepsilon ({\mathbf u}^n_h))
-  -
-  (\nabla \cdot {\mathbf v}_h, p^n_h) 
-  &=& 
-  ({\mathbf v}_h, \mathrm{Ra} \; T^{n-1}_h \mathbf{g}),
-  \\
-  (q_h, \nabla \cdot {\mathbf u}^n_h) &=& 0.
-@f}
-
-This is exactly the same as what we already discussed in 
-@ref step_22 "step-22" and there is not much more to say about this here.
-
-
-<h4>Stabilization, weak form and space discretization for the temperature equation</h4>
-
-The more interesting question is what to do with the temperature
-advection-diffusion equation. By default, not all discretizations of
-this equation are equally stable unless we either do something like
-upwinding, stabilization, or all of this. One way to achieve this is
-to use discontinuous elements (i.e. the FE_DGQ class that we used, for
-example, in the discretization of the transport equation in 
-@ref step_12 "step-12", or in discretizing the pressure in 
-@ref step_20 "step-20" and @ref step_21 "step-21") and to define a
-flux at the interface between cells that takes into account
-upwinding. If we had a pure advection problem this would probably be
-the simplest way to go. However, here we have some diffusion as well,
-and the discretization of the Laplace operator with discontinuous
-elements is cumbersome because of the significant number of additional
-terms that need to be integrated on each face between
-cells. Discontinuous elements also have the drawback that the use of
-numerical fluxes introduces an additional numerical diffusion that
-acts everywhere, whereas we would really like to minimize the effect
-of numerical diffusion to a minimum and only apply it where it is
-necessary to stabilize the scheme.
-
-A better alternative is therefore to add some nonlinear viscosity to
-the model. Essentially, what this does is to transform the temperature
-equation from the form 
-@f{eqnarray*}
-  \frac{\partial T}{\partial t}
-  +
-  {\mathbf u} \cdot \nabla T
-  -
-  \nabla \cdot \kappa \nabla T &=& \gamma
-@f}
-to something like
-@f{eqnarray*}
-  \frac{\partial T}{\partial t}
-  +
-  {\mathbf u} \cdot \nabla T
-  -
-  \nabla \cdot (\kappa+\nu(T)) \nabla T &=& \gamma,
-@f}
-where $\nu(T)$ is an addition viscosity (diffusion) term that only
-acts in the vicinity of shocks and other discontinuities. $\nu(T)$ is
-chosen in such a way that if <i>T</i> satisfies the original equations, the
-additional viscosity is zero.
-
-To achieve this, the literature contains a number of approaches. We
-will here follow one developed by Guermond and Popov that builds on a
-suitably defined residual and a limiting procedure for the additional
-viscosity. To this end, let us define a residual $R_\alpha(T)$ as follows:
-@f{eqnarray*}
-  R_\alpha(T)
-  =
-  \left(
-  \frac{\partial T}{\partial t}
-  +
-  {\mathbf u} \cdot \nabla T
-  -
-  \nabla \cdot \kappa \nabla T - \gamma
-  \right)
-  T^{\alpha-1}
-@f}
-where we will later choose the stabilization exponent $\alpha$ from
-within the range $[1,2]$. Note that $R_\alpha(T)$ will be zero if $T$
-satisfies the temperature equation, since then the term in parentheses
-will be zero. Multiplying terms out, we get the following, entirely
-equivalent form:
-@f{eqnarray*}
-  R_\alpha(T)
-  =
-  \frac 1\alpha
-  \frac{\partial (T^\alpha)}{\partial t}
-  +
-  \frac 1\alpha
-  {\mathbf u} \cdot \nabla (T^\alpha)
-  -
-  \frac 1\alpha
-  \nabla \cdot \kappa \nabla (T^\alpha)
-  +
-  \kappa(\alpha-1)
-  T^{\alpha-2} |\nabla T|^\alpha
-  - 
-  \gamma
-  T^{\alpha-1}
-@f}
-
-With this residual, we can now define the artificial viscosity as
-a piecewise constant function defined on each cell $K$ with diameter
-$h_K$ separately as
-follows:
-@f{eqnarray*}
-  \nu_\alpha(T)|_K
-  =
-  \beta
-  \|\mathbf{u}\|_{L^\infty(K)}
-  \min\left\{
-    h_K,
-    h_K^\alpha
-    \frac{\|R_\alpha(T)\|_{L^\infty(K)}}{c(\mathbf{u},T)}
-  \right\}
-@f}
-
-Here, $\beta$ is a stabilization constant (a dimensional analysis
-reveals that it is unitless and therefore independent of scaling; we will
-discuss its choice in the <a href="#Results">results section</a>) 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) = 
- c_R\ \|\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) and $c_R$ is a dimensionless
-constant. 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
-the stabilization term that injects artificial diffusion will there be
-of size $h_K^{\alpha+1}$ &mdash; i.e. rather small, just as we hope it to
-be when no additional diffusion is necessary. On the other hand, if we
-are on or close to a discontinuity of the temperature field, then the
-residual will be large; the minimum operation in the definition of
-$\nu_\alpha(T)$ will then ensure that the stabilization has size $h_K$
-&mdash; the optimal amount of artificial viscosity to ensure stability of
-the scheme.
-
-It is certainly a good questions whether this scheme really works?
-Computations by Guermond and Popov have shown that this form of
-stabilization actually performs much better than most of the other
-stabilization schemes that are around (for example streamline
-diffusion, to name only the simplest one). Furthermore, for $\alpha\in
-[1,2)$ they can even prove that it produces better convergence orders
-for the linear transport equation than for example streamline
-diffusion. For $\alpha=2$, no theoretical results are currently
-available, but numerical tests indicate that the results
-are considerably better than for $\alpha=1$.
-
-A more practical question is how to introduce this artificial
-diffusion into the equations we would like to solve. Note that the
-numerical viscosity $\nu(T)$ is temperature-dependent, so the equation
-we want to solve is nonlinear in <i>T</i> &mdash; not what one desires from a
-simple method to stabilize an equation, and even less so if we realize
-that $\nu(T)$ is non-differentiable in <i>T</i>. However, there is no
-reason to despair: we still have to discretize in time and we can
-treat the term explicitly.
-
-In the definition of the stabilization parameter, we approximate the time
-derivative by $\frac{\partial T}{\partial t} \approx 
-\frac{T^{n-1}-T^{n-2}}{k^{n-1}}$. This approximation makes only use
-of available time data and this is the reason why we need to store data of two 
-previous time steps (which enabled us to use the BDF-2 scheme without
-additional storage cost). We could now simply evaluate the rest of the
-terms at $t_{n-1}$, but then the discrete residual would be nothing else than
-a backward Euler approximation, which is only first order accurate. So, in
-case of smooth solutions, the residual would be still of the order <i>h</i>, 
-despite the second order time accuracy in the outer BDF-2 scheme and the 
-spatial FE discretization. This is certainly not what we want to have
-(in fact, we desired to have small residuals in regions where the solution
-behaves nicely), so a bit more care is needed. The key to this problem
-is to observe that the first derivative as we constructed it is actually 
-centered at $t_{n-\frac{3}{2}}$. We get the desired second order accurate
-residual calculation if we evaluate all spatial terms at $t_{n-\frac{3}{2}}$
-by using the approximation $\frac 12 T^{n-1}+\frac 12 T^{n-2}$, which means
-that we calculate the nonlinear viscosity as a function of this
-intermediate temperature, $\nu_\alpha = 
-\nu_\alpha\left(\frac 12 T^{n-1}+\frac 12 T^{n-2}\right)$. Note that this
-evaluation of the residual is nothing else than a Crank-Nicholson scheme,
-so we can be sure that now everything is alright. One might wonder whether
-it is a problem that the numerical viscosity now is not evaluated at 
-time <i>n</i> (as opposed to the rest of the equation). However, this offset
-is uncritical: For smooth solutions, $\nu_\alpha$ will vary continuously,
-so the error in time offset is <i>k</i> times smaller than the nonlinear
-viscosity itself, i.e., it is a small higher order contribution that is 
-left out. That's fine because the term itself is already at the level of
-discretization error in smooth regions.
-
-Using the BDF-2 scheme introduced above,
-this yields for the simpler case of uniform time steps of size <i>k</i>:
-@f{eqnarray*}
-  \frac 32 T^n
-  -
-  k\nabla \cdot \kappa \nabla T^n 
-  &=& 
-  2 T^{n-1}
-  -
-  \frac 12 T^{n-2}
-  \\
-  &&
-  +
-  k\nabla \cdot 
-  \left[
-    \nu_\alpha\left(\frac 12 T^{n-1}+\frac 12 T^{n-2}\right) 
-    \ \nabla (2T^{n-1}-T^{n-2})
-  \right]
-  \\
-  &&
-  -
-  k{\mathbf u}^n \cdot \nabla (2T^{n-1}-T^{n-2})
-  \\
-  &&
-  +
-  k\gamma.
-@f}
-On the left side of this equation remains the term from the time
-derivative and the original (physical) diffusion which we treat
-implicitly (this is actually a nice term: the matrices that result
-from the left hand side are the mass matrix and a multiple of the
-Laplace matrix &mdash; both are positive definite and if the time step
-size <i>k</i> is small, the sum is simple to invert). On the right hand
-side, the terms in the first line result from the time derivative; in
-the second line is the artificial diffusion at time $t_{n-\frac
-32}$; the third line contains the
-advection term, and the fourth the sources. Note that the
-artificial diffusion operates on the extrapolated
-temperature at the current time in the same way as we have discussed
-the advection works in the section on time stepping.
-
-The form for non-uniform time steps that we will have to use in
-reality is a bit more complicated (which is why we showed the simpler
-form above first) and reads:
-@f{eqnarray*}
-  \frac{2k_n+k_{n-1}}{k_n+k_{n-1}} T^n
-  -
-  k_n\nabla \cdot \kappa \nabla T^n 
-  &=& 
-  \frac{k_n+k_{n-1}}{k_{n-1}} T^{n-1}
-  -
-  \frac{k_n^2}{k_{n-1}(k_n+k_{n-1})} T^{n-2}
-  \\
-  &&
-  +
-  k_n\nabla \cdot 
-  \left[
-    \nu_\alpha\left(\frac 12 T^{n-1}+\frac 12 T^{n-2}\right) 
-    \ \nabla  \left[
-    \left(1+\frac{k_n}{k_{n-1}}\right)T^{n-1}-\frac{k_n}{k_{n-1}}T^{n-2}
-  \right]
-  \right]
-  \\
-  &&
-  -
-  k_n{\mathbf u}^n \cdot \nabla \left[
-    \left(1+\frac{k_n}{k_{n-1}}\right)T^{n-1}-\frac{k_n}{k_{n-1}}T^{n-2}
-  \right]
-  \\
-  &&
-  +
-  k_n\gamma.
-@f}
-
-After settling all these issues, the weak form follows naturally from
-the strong form shown in the last equation, and we immediately arrive
-at the weak form of the discretized equations:
-@f{eqnarray*}
-  \frac{2k_n+k_{n-1}}{k_n+k_{n-1}} (\tau_h,T_h^n)
-  +
-  k_n (\nabla \tau_h, \kappa \nabla T_h^n)
-  &=& 
-  \biggl(\tau_h, 
-  \frac{k_n+k_{n-1}}{k_{n-1}} T_h^{n-1}
-  -
-  \frac{k_n^2}{k_{n-1}(k_n+k_{n-1})} T_h^{n-2}
-  \\
-  &&\qquad\qquad
-  -
-  k_n{\mathbf u}_h^n \cdot \nabla \left[
-    \left(1+\frac{k_n}{k_{n-1}}\right)T^{n-1}-\frac{k_n}{k_{n-1}}T^{n-2}
-  \right]
-  +
-  k_n\gamma \biggr)
-  \\
-  &&
-  -
-  k_n \left(\nabla \tau_h, 
-    \nu_\alpha\left(\frac 12 T_h^{n-1}+\frac 12 T_h^{n-2}\right) 
-    \ \nabla \left[
-    \left(1+\frac{k_n}{k_{n-1}}\right)T^{n-1}-\frac{k_n}{k_{n-1}}T^{n-2}
-  \right]
-  \right)
-@f}
-for all discrete test functions $\tau_h$. Here, the diffusion term has been
-integrated by parts, and we have used that we will impose no thermal flux,
-$\mathbf{n}\cdot\kappa\nabla T|_{\partial\Omega}=0$.
-
-This then results in a
-matrix equation of form
-@f{eqnarray*}
-  \left( \frac{2k_n+k_{n-1}}{k_n+k_{n-1}} M+k_n A_T\right) T_h^n = F(U_h^n,T_h^{n-1},T_h^{n-2}),
-@f}
-which given the structure of matrix on the left (the sum of two
-positive definite matrices) is easily solved using the Conjugate
-Gradient method.
-
-
-
-<h4>Linear solvers</h4>
-
-As explained above, our approach to solving the joint system for
-velocities/pressure on the one hand and temperature on the other is to use an
-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. 
-
-
-<h5>Linear solvers for the Stokes problem</h5>
-
-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 that turned out to be more
-efficient than the original approach. The best alternative
-identified there we to use a GMRES solver preconditioned by a block
-matrix involving the Schur complement. Specifically, the Stokes
-operator leads to a block structured matrix
-@f{eqnarray*}
-  \left(\begin{array}{cc}
-    A & B^T \\ B & 0
-  \end{array}\right)
-@f}
-and as discussed there a good preconditioner is
-@f{eqnarray*}
-  P^{-1}
-  = 
-  \left(\begin{array}{cc}
-    A^{-1} & 0 \\ S^{-1} B A^{-1} & -S^{-1}
-  \end{array}\right)
-@f}
-where <i>S</i> is the Schur complement of the Stokes operator
-$S=B^TA^{-1}B$. Of course, this preconditioner is not useful because we
-can't form the various inverses of matrices, but we can use the
-following as a preconditioner:
-@f{eqnarray*}
-  \tilde P^{-1}
-  = 
-  \left(\begin{array}{cc}
-    \tilde A^{-1} & 0 \\ \tilde S^{-1} B \tilde A^{-1} & -\tilde S^{-1}
-  \end{array}\right)
-@f}
-where $\tilde A^{-1},\tilde S^{-1}$ are approximations to the inverse
-matrices. In particular, it turned out that <i>S</i> is spectrally
-equivalent to the mass matrix and consequently replacing $\tilde
-S^{-1}$ by a CG solver applied to the mass matrix on the pressure
-space was a good choice.
-
-It was more complicated to come up with a good replacement $\tilde
-A^{-1}$, which corresponds to the discretized symmetric Laplacian of
-the vector-valued velocity field, i.e. 
-$A_{ij} = (\varepsilon {\mathbf v}_i, \eta \varepsilon ({\mathbf
-v}_j))$.
-In @ref step_22 "step-22" we used a sparse LU decomposition (using the
-SparseDirectUMFPACK class) of <i>A</i> for $\tilde A^{-1}$ &mdash; the
-perfect preconditioner &mdash; in 2d, but for 3d memory and compute
-time is not usually sufficient to actually compute this decomposition;
-consequently, we only use an incomplete LU decomposition (ILU, using
-the SparseILU class) in 3d.
-
-For this program, we would like to go a bit further. To this end, note
-that the symmetrized bilinear form on vector fields, 
-$(\varepsilon {\mathbf v}_i, \eta \varepsilon ({\mathbf v}_j))$
-is not too far away from the nonsymmetrized version,
-$(\nabla {\mathbf v}_i, \eta \nabla {\mathbf v}_j)
-= \sum_{k,l=1}^d 
-  (\partial_k ({\mathbf v}_i)_l, \eta \partial_k ({\mathbf v}_j)_l)
-$. The latter,
-however, has the advantage that the <code>dim</code> vector components
-of the test functions are not coupled (well, almost, see below),
-i.e. the resulting matrix is block-diagonal: one block for each vector
-component, and each of these blocks is equal to the Laplace matrix for
-this vector component. So assuming we order degrees of freedom in such
-a way that first all <i>x</i>-components of the velocity are numbered, then
-the <i>y</i>-components, and then the <i>z</i>-components, then the matrix
-$\hat A$ that is associated with this slightly different bilinear form has
-the form
-@f{eqnarray*}
-  \hat A =
-  \left(\begin{array}{ccc}
-    A_s & 0 & 0 \\ 0 & A_s & 0 \\ 0 & 0 & A_s
-  \end{array}\right)
-@f}
-where $A_s$ is a Laplace matrix of size equal to the number of shape functions
-associated with each component of the vector-valued velocity. With this
-matrix, one could be tempted to define our preconditioner for the
-velocity matrix <i>A</i> as follows:
-@f{eqnarray*}
-  \tilde A^{-1} =
-  \left(\begin{array}{ccc}
-    \tilde A_s^{-1} & 0 & 0 \\
-    0 & \tilde A_s^{-1} & 0 \\ 
-    0 & 0 & \tilde A_s^{-1}
-  \end{array}\right),
-@f}
-where $\tilde A_s^{-1}$ is a preconditioner for the Laplace matrix &mdash;
-something where we know very well how to build good preconditioners! 
-
-In reality, the story is not quite as simple: To make the matrix
-$\tilde A$ definite, we need to make the individual blocks $\tilde
-A_s$ definite by applying boundary conditions. One can try to do so by
-applying Dirichlet boundary conditions all around the boundary, and
-then the so-defined preconditioner $\tilde A^{-1}$ turns out to be a
-good preconditioner for <i>A</i> if the latter matrix results from a Stokes
-problem where we also have Dirichlet boundary conditions on the
-velocity components all around the domain, i.e. if we enforce <b>u</b>=0. 
-
-Unfortunately, this "if" is an "if and only if": in the program below
-we will want to use no-flux boundary conditions of the form $\mathbf u
-\cdot \mathbf n = 0$ (i.e. flow parallel to the boundary is allowed,
-but no flux through the boundary). In this case, it turns out that the
-block diagonal matrix defined above is not a good preconditioner
-because it neglects the coupling of components at the boundary. A
-better way to do things is therefore if we build the matrix $\hat A$
-as the vector Laplace matrix $\hat A_{ij} = (\nabla {\mathbf v}_i,
-\eta \nabla {\mathbf v}_j)$ and then apply the same boundary condition
-as we applied to <i>A</i>. If this is Dirichlet boundary conditions all
-around the domain, the $\hat A$ will decouple to three diagonal blocks
-as above, and if the boundary conditions are of the form $\mathbf u
-\cdot \mathbf n = 0$ then this will introduce a coupling of degrees of
-freedom at the boundary but only there. This, in fact, turns out to be
-a much better preconditioner than the one introduced above, and has
-almost all the benefits of what we hoped to get.
-
-
-To sum this whole story up, we can observe:
-<ul>
-  <li> Compared to building a preconditioner from the original matrix <i>A</i>
-  resulting from the symmetric gradient as we did in @ref step_22 "step-22",
-  we have to expect that the preconditioner based on the Laplace bilinear form
-  performs worse since it does not take into account the coupling between
-  vector components.
-
-  <li>On the other hand, preconditioners for the Laplace matrix are typically
-  more mature and perform better than ones for vector problems. For example,
-  at the time of this writing, Algebraic Multigrid (AMG) algorithms are very
-  well developed for scalar problems, but not so for vector problems.
-
-  <li>In building this preconditioner, we will have to build up the
-  matrix $\hat A$ and its preconditioner. While this means that we
-  have to store an additional matrix we didn't need before, the
-  preconditioner $\tilde A_s^{-1}$ is likely going to need much less
-  memory than storing a preconditioner for the coupled matrix
-  <i>A</i>. This is because the matrix $A_s$ has only a third of the
-  entries per row for all rows corresponding to interior degrees of
-  freedom, and contains coupling between vector components only on
-  those parts of the boundary where the boundary conditions introduce
-  such a coupling. Storing the matrix is therefore comparatively
-  cheap, and we can expect that computing and storing the
-  preconditioner $\tilde A_s$ will also be much cheaper compared to
-  doing so for the fully coupled matrix.
-</ul>
-
-
-
-<h5>Linear solvers for the temperature equation</h5>
-
-This is the easy part: The matrix for the temperature equation has the form
-$\alpha M + \beta A$, where $M,A$ are mass and stiffness matrices on the
-temperature space, and $\alpha,\beta$ are constants related the time stepping
-scheme and the current and previous time step. This being the sum of a
-symmetric positive definite and a symmetric positive semidefinite matrix, the
-result is also symmetric positive definite. Furthermore, $\frac\beta\alpha$ is
-a number proportional to the time step, and so becomes small whenever the mesh
-is fine, damping the effect of the then ill-conditioned stiffness matrix.
-
-As a consequence, inverting this matrix with the Conjugate Gradient algorithm,
-using a simple preconditioner, is trivial and very cheap compared to inverting
-the Stokes matrix.
-
-
-
-<h3>Implementation details</h3>
-
-One of the things worth explaining up front about the program below is the use
-of two different DoFHandler objects. If one looks at the structure of the
-equations above and the scheme for their solution, one realizes that there is
-little commonality that keeps the Stokes part and the temperature part
-together. In all previous tutorial programs in which we have discussed @ref
-vector_valued "vector-valued problems" we have always only used a single
-finite element with several vector components, and a single DoFHandler object.
-Sometimes, we have substructured the resulting matrix into blocks to
-facilitate particular solver schemes; this was, for example, the case in the
-@ref step_22 "step-22" program for the Stokes equations upon which the current
-program is based.
-
-We could of course do the same here. The linear system that we would get would
-look like this:
-@f{eqnarray*}
-  \left(\begin{array}{ccc}
-    A & B^T & 0 \\ B & 0 &0 \\ C & 0 & K
-  \end{array}\right)
-  \left(\begin{array}{ccc}
-    U^n \\ P^n \\ T^n
-  \end{array}\right)
-  =
-  \left(\begin{array}{ccc}
-    F_U(T^{n-1}) \\ 0 \\ F_T(U^n,T^{n-1},T^{n-1})
-  \end{array}\right).
-@f}
-The problem with this is: We never use the whole matrix at the same time. In
-fact, it never really exists at the same time: As explained above, $K$ and
-$F_T$ depend on the already computed solution $U^n$, in the first case through
-the time step (that depends on $U^n$ because it has to satisfy a CFL
-condition). So we can only assemble it once we've already solved the top left
-$2\times 2$ block Stokes system, and once we've moved on to the temperature
-equation we don't need the Stokes part any more. Furthermore, we don't
-actually build the matrix $C$: Because by the time we get to the temperature
-equation we already know $U^n$, and because we have to assemble the right hand
-side $F_T$ at this time anyway, we simply move the term $CU^n$ to the right
-hand side and assemble it along with all the other terms there. What this
-means is that there does not remain a part of the matrix where temperature
-variables and Stokes variables couple, and so a global enumeration of all
-degrees of freedom is no longer important: It is enough if we have an
-enumeration of all Stokes degrees of freedom, and of all temperature degrees
-of freedom independently.
-
-In essence, there is consequently not much use in putting <i>everything</i>
-into a block matrix (though there are of course the same good reasons to do so
-for the $2\times 2$ Stokes part), or, for that matter, in putting everything
-into the same DoFHandler object.
-
-But are there <i>downsides</i> to doing so? These exist, though they may not
-be obvious at first. The main problem is that if we need to create one global
-finite element that contains velocity, pressure, and temperature shape
-functions, and use this to initialize the DoFHandler. But we also use this
-finite element object to initialize all FEValues or FEFaceValues objects that
-we use. This may not appear to be that big a deal, but imagine what happens
-when, for example, we evaluate the residual
-$
-  R_\alpha(T)
-  =
-  \left(
-  \frac{\partial T}{\partial t}
-  +
-  {\mathbf u} \cdot \nabla T
-  -
-  \nabla \cdot \kappa \nabla T - \gamma
-  \right)
-  T^{\alpha-1}
-$
-that we need to compute the artificial viscosity $\nu_\alpha(T)|_K$.  For
-this, we need the Laplacian of the temperature, which we compute using the
-tensor of second derivatives (Hessians) of the shape functions (we have to
-give the <code>update_hessians</code> flag to the FEValues object for
-this). Now, if we have a finite that contains the shape functions for
-velocities, pressures, and temperatures, that means that we have to compute
-the Hessians of <i>all</i> shape functions, including the many higher order
-shape functions for the velocities. That's a lot of computations that we don't
-need, and indeed if one were to do that (as we had in an early version of the
-program), assembling the right hand side took about a quarter of the overall
-compute time.
-
-So what we will do is to use two different finite element objects, one for the
-Stokes components and one for the temperatures. With this come two different
-DoFHandlers, two sparsity patterns and two matrices for the Stokes and
-temperature parts, etc. And whenever we have to assemble something that
-contains both temperature and Stokes shape functions (in particular the right
-hand sides of Stokes and temperature equations), then we use two FEValues
-objects initialized with two cell iterators that we walk in parallel through
-the two DoFHandler objects associated with the same Triangulation object; for
-these two FEValues objects, we use of course the same quadrature objects so
-that we can iterate over the same set of quadrature points, but each FEValues
-object will get update flags only according to what it actually needs to
-compute. In particular, when we compute the residual as above, we only ask for
-the values of the Stokes shape functions, but also the Hessians of the
-temperature shape functions &mdash; much cheaper indeed, and as it turns out:
-assembling the right hand side of the temperature equation is now a component
-of the program that is hardly measurable.
-
-With these changes, timing the program yields that only the following
-operations are relevant for the overall run time: 
-<ul>
-  <li>Solving the Stokes system: 72% of the run time.
-  <li>Assembling the Stokes preconditioner and computing the algebraic
-      multigrid hierarchy using the Trilinos ML package: 11% of the
-      run time.
-  <li>The function <code>BoussinesqFlowProblem::setup_dofs</code>: 7%
-      of overall run time.
-  <li>Assembling the Stokes and temperature right hand side vectors as
-      well as assembling the matrices: 7%.
-</ul>
-In essence this means that all bottlenecks apart from the algebraic
-multigrid have been removed.
index 6a2f055b1adbf360a96391516bdb3e1875a3ef5f..a79c9a63f464305f450510dd70cf84e21f967593 100644 (file)
@@ -1,275 +1,3 @@
 <h1>Results</h1>
 
 
-
-<h3> Numerical experiments to determine optimal parameters </h3>
-
-The program as is has three parameters that we don't have much of a
-theoretical handle on how to choose in an optimal way. These are:
-<ul>
-  <li>The time step must satisfy a CFL condition
-      $k\le \min_K \frac{c_kh_K}{\|\mathbf{u}\|_{L^\infty(K)}}$. Here, $c_k$ is 
-      dimensionless, but what is the right value?
-  <li>In the computation of the artificial viscosity,
-@f{eqnarray*}
-  \nu_\alpha(T)|_K
-  =
-  \beta
-  \|\mathbf{u}\|_{L^\infty(K)}
-  \min\left\{
-    h_K,
-    h_K^\alpha
-    \frac{\|R_\alpha(T)\|_{L^\infty(K)}}{c(\mathbf{u},T)}
-  \right\},
-@f}
-      with $c(\mathbf{u},T) = 
-      c_R\ \|\mathbf{u}\|_{L^\infty(\Omega)} \ \mathrm{var}(T)
-      \ |\mathrm{diam}(\Omega)|^{\alpha-2}$.
-      Here, the choice of the dimensionless numbers $\beta,c_R$ is of
-      interest.
-</ul>
-In all of these cases, we will have to expect that the correct choice of each
-value depends on that of the others, and most likely also on the space
-dimension and polynomial degree of the finite element used for the
-temperature. Below we'll discuss a few numerical experiments to choose
-constants. 
-
-
-<h4> Choosing <i>c<sub>k</sub></i> and &beta; </h4>
-
-These two constants are definitely linked in some way. The reason is easy to
-see: In the case of a pure advection problem, 
-$\frac{\partial T}{\partial t} + \mathbf{u}\cdot\nabla T = \gamma$, any
-explicit scheme has to satisfy a CFL condition of the form 
-$k\le \min_K \frac{c_k^a h_K}{\|\mathbf{u}\|_{L^\infty(K)}}$. On the other hand,
-for a pure diffusion problem, 
-$\frac{\partial T}{\partial t} + \nu \Delta T = \gamma$,
-explicit schemes need to satisfy a condition
-$k\le \min_K \frac{c_k^d h_K^2}{\nu}$. So given the form of $\nu$ above, an
-advection diffusion problem like the one we have to solve here will result in
-a condition of the form
-$
-k\le \min_K \min \left\{
-  \frac{c_k^a h_K}{\|\mathbf{u}\|_{L^\infty(K)}},
-  \frac{c_k^d h_K^2}{\beta \|mathbf{u}\|_{L^\infty(K)} h_K}\right\}
-  =
-  \min_K \left( \min \left\{
-  c_k^a,
-  \frac{c_k^d}{\beta}\right\}
-  \frac{h_K}{\|\mathbf{u}\|_{L^\infty(K)}} \right)
-$.
-It follows that we have to face the fact that we might want to choose $\beta$
-larger to improve the stability of the numerical scheme (by increasing the
-amount of artificial diffusion), but we have to pay a price in the form of
-smaller, and consequently more time steps. In practice, one would therefore
-like to choose $\beta$ as small as possible to keep the transport problem
-sufficiently stabilized while at the same time trying to choose the time step
-as large as possible to reduce the overall amount of work.
-
-The find the right balance, the only way is to do a few computational
-experiments. Here's what we did: We modified the program slightly to allow
-less mesh refinement (so we don't always have to wait that long) and to choose 
-$
-  \nu(T)|_K
-  =
-  \beta
-  \|\mathbf{u}\|_{L^\infty(K)} h_K
-$ to eliminate the effect of of the constant $c_R$. We then run the program
-for different values $c_k,\beta$ and observe maximal and minimal temperatures
-in the domain. What we expect to see is this: If we choose the time step too
-big (i.e. choose a $c_k$ bigger than theoretically allowed) then we will get
-exponential growth of the temperature. If we choose $\beta$ too small, then
-the transport stabilization becomes insufficient and the solution will show
-significant oscillations but not exponential growth. 
-
-
-<h5>Results for Q<sub>1</sub> elements</h5>
-
-Here is what we get for
-$\beta=0.01, \beta=0.1$, and $\beta=0.5$, different choices of $c_k$, and
-bilinear elements (<code>temperature_degree=1</code>) in 2d:
-
-<table align="center" border="1" cellspacing="3" cellpadding="3">
-  <tr>
-    <td>
-        @image html "step-33.timestep.q1.beta=0.01.png" "" width=4cm
-    </td>
-    <td>
-        @image html "step-33.timestep.q1.beta=0.03.png" "" width=4cm
-    </td>
-  </tr>
-
-  <tr>
-    <td>
-        @image html "step-33.timestep.q1.beta=0.1.png" "" width=4cm
-    </td>
-    <td>
-        @image html "step-33.timestep.q1.beta=0.5.png" "" width=4cm
-    </td>
-  </tr>
-</table>
-
-The way to interpret these graphs goes like this: for $\beta=0.01$ and
-$c_k=\frac 12,\frac 14$, we see exponential growth or at least large
-variations, but if we choose 
-$k=\frac 18\frac{h_K}{\|\mathbf{u}\|_{L^\infty(K)}}$ 
-or smaller, then the scheme is 
-stable though a bit wobbly. For more artificial diffusion, we can choose 
-$k=\frac 14\frac{h_K}{\|\mathbf{u}\|_{L^\infty(K)}}$ 
-or smaller for $\beta=0.03$,
-$k=\frac 13\frac{h_K}{\|\mathbf{u}\|_{L^\infty(K)}}$ 
-or smaller for $\beta=0.1$, and again need
-$k=\frac 1{15}\frac{h_K}{\|\mathbf{u}\|_{L^\infty(K)}}$ 
-for $\beta=0.5$ (this time because much diffusion requires a small time
-step). 
-
-So how to choose? If we were simply interested in a large time step, then we
-would go with $\beta=0.1$ and
-$k=\frac 13\frac{h_K}{\|\mathbf{u}\|_{L^\infty(K)}}$.
-On the other hand, we're also interested in accuracy and here it may be of
-interest to actually investigate what these curves show. To this end note that
-we start with a zero temperature and that our sources are positive &mdash; so
-we would intuitively expect that the temperature can never drop below
-zero. But it does, a consequence of Gibb's phenomenon when using continuous
-elements to approximate a discontinuous solution. We can therefore see that
-choosing $\beta$ too small is bad: too little artificial diffusion leads to
-over- and undershoots that aren't diffused away. On the other hand, for large
-$\beta$, the minimum temperature drops below zero at the beginning but then
-quickly diffuses back to zero. 
-
-On the other hand, let's also look at the maximum temperature. Watching the
-movie of the solution, we see that initially the fluid is at rest. The source
-keeps heating the same volume of fluid whose temperature increases linearly at
-the beginning until its buoyancy is able to move it upwards. The hottest part
-of the fluid is therefore transported away from the solution and fluid taking
-its place is heated for only a short time before being moved out of the source
-region, therefore remaining cooler than the initial bubble. If $\kappa=0$
-(in the program it is nonzero but very small) then the hottest part of the
-fluid should be advected along with the flow with its temperature
-constant. That's what we can see in the graphs with the smallest $\beta$: Once
-the maximum temperature is reached, it hardly changes any more. On the other
-hand, the larger the artificial diffusion, the more the hot spot is
-diffused. Note that for this criterion, the time step size does not play a
-significant role.
-
-So to sum up, likely the best choice would appear to be $\beta=0.03$
-and $k=\frac 14\frac{h_K}{\|\mathbf{u}\|_{L^\infty(K)}}$. The curve is
-a bit wobbly, but overall pictures looks pretty reasonable with the
-exception of some over and undershoots close to the start time due to
-Gibb's phenomenon.
-
-
-<h5>Results for Q<sub>2</sub> elements</h5>
-
-One can repeat the same sequence of experiments for higher order
-elements as well. Here are the graphs for bi-quadratic shape functions
-(<code>temperature_degree=2</code>) for the temperature, while we
-retain the $Q_2/Q_1$ stable Taylor-Hood element for the Stokes system:
-
-<table align="center" border="1" cellspacing="3" cellpadding="3">
-  <tr>
-    <td>
-        @image html "step-33.timestep.q2.beta=0.01.png" "" width=4cm
-    </td>
-    <td>
-        @image html "step-33.timestep.q2.beta=0.03.png" "" width=4cm
-    </td>
-  </tr>
-
-  <tr>
-    <td>
-        @image html "step-33.timestep.q2.beta=0.1.png" "" width=4cm
-    </td>
-  </tr>
-</table>
-
-Again, small values of $\beta$ lead to less diffusion but we have to
-choose the time step very small to keep things under control. Too
-large values of $\beta$ make for more diffusion, but again require
-small time steps. The best value would appear to be $\beta=0.03$, as
-for the $Q_1$ element, and the we have to choose 
-$k=\frac 18\frac{h_K}{\|\mathbf{u}\|_{L^\infty(K)}}$ &mdash; exactly
-half the size for the $Q_1$ element, a fact that may not be surprising
-if we state the CFL condition as the requirement that the time step be
-small enough so that the distance transport advects in each time step
-is no longer than one <i>grid point</i> away (which for $Q_1$ elements
-is $h_K$, but for $Q_2$ elements is $h_K/2$).
-
-
-<h5>Results for 3d</h5>
-
-One can repeat these experiments in 3d and find the optimal time step
-for each value of $\beta$ and find the best value of $\beta$. What one
-finds is that for the same $\beta$ already used in 2d, the time steps
-needs to be a bit small, by around a factor of 1.2 or so. This is
-easily explained: the time step restriction is 
-$k=\min_K \frac{ch_K}{\|\mathbf{u}\|_{L^\infty(K)}}$ where $h_K$ is
-the <i>diameter</i> of the cell. However, what is really needed is the
-distance between mesh points, which is $\frac{h_K}{\sqrt{d}}$. So a
-more appropriate form would be
-$k=\min_K \frac{ch_K}{\|\mathbf{u}\|_{L^\infty(K)}\sqrt{d}}$.
-
-The second find is that one needs to choose $\beta$ slightly bigger
-(about $\beta=0.05$ or so). This then again reduces the time step we
-can take.
-
-
-
-
-<h5>Conclusions</h5>
-
-Concluding, $\beta=0.03$ appears to be a good choice for the
-stabilization parameter in 2d, and $\beta=0.05$ in 3d. In a dimension
-independent way, we can model this as $\beta=0.015d$. As we have seen
-in the sections above, in 2d
-$k=\frac 14 \frac 1{q_T}\frac{h_K}{\|\mathbf{u}\|_{L^\infty(K)}}$
-is an appropriate time step, where $q_T$ is the polynomial degree of
-the temperature shape functions (in the program, this corresponds to
-the variable <code>temperature_degree</code>). To reconcile this with
-the findings in 3d for the same $\beta$, we could write this as
-$k=\frac 1{2\sqrt{2}\sqrt{d}} \frac
-1{q_T}\frac{h_K}{\|\mathbf{u}\|_{L^\infty(K)}}$
-but this doesn't take into account that we also have to increase
-$\beta$ in 3d. The final form that takes all these factors in reads as
-follows:
-@f{eqnarray*}
-  k =
-  \frac 1{2\sqrt{2}} \frac 1{\sqrt{d}}
-  \frac 2d
-  \frac 1{q_T}
-  \frac{h_K}{\|\mathbf{u}\|_{L^\infty(K)}}
-  =
-  \frac 1{d\sqrt{2}\sqrt{d}}
-  \frac 1{q_T}
-  \frac{h_K}{\|\mathbf{u}\|_{L^\infty(K)}}.
-@f}
-In the first form (in the center of the equation), $\frac
-1{2\sqrt{2}}$ is a universal constant, $\frac 1{\sqrt{d}}$
-is the factor that accounts for the difference between cell diameter
-and grid point separation,
-$\frac 2d$ accounts for the increase in $\beta$ with space dimension,
-$\frac 1{q_T}$ accounts for the distance between grid points for
-higher order elements, and $\frac{h_K}{\|\mathbf{u}\|_{L^\infty(K)}}$
-for the local speed of transport relative to the cell size. This is
-the formula that we use in the program.
-
-As for the question of whether to use $Q_1$ or $Q_2$ elements for the
-temperature, the following considerations may be useful: First,
-solving the temperature equation is hardly a factor in the overall
-scheme since almost the entire compute time goes into solving the
-Stokes system in each time step. Higher order elements for the
-temperature equation are therefore not a significant drawback. On the
-other hand, if one compares the size of the over- and undershoots the
-solution produces due to the discontinuous source description, one
-notices that for the choice of $\beta$ and $k$ as above, the $Q_1$
-solution dips down to around $-0.47$, whereas the $Q_2$ solution only
-goes to $-0.13$ (remember that the exact solution should never become
-negative at all. This means that the $Q_2$ solution is significantly
-more accurate; the program therefore uses these higher order elements,
-despite the penalty we pay in terms of smaller time steps.
-
-
-<h3> Possible extensions </h3>
-
-Parallelization -> step-33
-
index 29361942555e9d9ff8f686ffde30fa30e2fe0247..7e8da2e2d76de23a856ca23a7878061e9f4e812c 100644 (file)
@@ -3,7 +3,7 @@
 
 /*    $Id$       */
 /*                                                                */
-/*    Copyright (C) 2007, 2008 by the deal.II authors */
+/*    Copyright (C) 2008 by the deal.II authors */
 /*                                                                */
 /*    This file is subject to QPL and may not be  distributed     */
 /*    without copyright and license information. Please refer     */
@@ -69,32 +69,6 @@ using namespace dealii;
 
                                 // @sect3{Equation data}
 
-                                // Again, the next stage in the program
-                                // is the definition of the equation 
-                                // data, that is, the various
-                                // boundary conditions, the right hand
-                                // side and the initial condition (remember
-                                // that we're about to solve a time-
-                                // dependent system). The basic strategy
-                                // for this definition is the same as in
-                                // step-22. Regarding the details, though,
-                                // there are some differences.
-
-                                // The first
-                                // thing is that we don't set any boundary
-                                // conditions on the velocity, as is
-                                // explained in the introduction. So
-                                // what is left are two conditions for
-                                // pressure <i>p</i> and temperature
-                                // <i>T</i>.
-
-                                // Secondly, we set an initial
-                                // condition for all problem variables,
-                                // i.e., for <b>u</b>, <i>p</i> and <i>T</i>,
-                                // so the function has <i>dim+2</i>
-                                // components.
-                                // In this case, we choose a very simple
-                                // test case, where everything is zero.
 
                                 // @sect4{Boundary values}
 namespace EquationData
@@ -162,18 +136,6 @@ namespace EquationData
 
 
                                   // @sect4{Right hand side}
-                                  // 
-                                  // The last definition of this kind
-                                  // is the one for the right hand
-                                  // side function. Again, the content
-                                  // of the function is very
-                                  // basic and zero in most of the
-                                  // components, except for a source
-                                  // of temperature in some isolated
-                                  // regions near the bottom of the
-                                  // computational domain, as is explained
-                                  // in the problem description in the
-                                  // introduction.
   template <int dim>
   class TemperatureRightHandSide : public Function<dim>
   {
@@ -226,34 +188,8 @@ namespace EquationData
 
                                   // @sect3{Linear solvers and preconditioners}
 
-                                  // This section introduces some
-                                  // objects that are used for the
-                                  // solution of the linear equations of
-                                  // Stokes system that we need to
-                                  // solve in each time step. The basic
-                                  // structure is still the same as
-                                  // in step-20, where Schur complement
-                                  // based preconditioners and solvers
-                                  // have been introduced, with the 
-                                  // actual interface taken from step-22.
 namespace LinearSolvers
 {
-
-                                  // @sect4{The <code>InverseMatrix</code> class template}
-
-                                  // This class is an interface to
-                                  // calculate the action of an
-                                  // "inverted" matrix on a vector
-                                  // (using the <code>vmult</code>
-                                  // operation)
-                                  // in the same way as the corresponding
-                                  // function in step-22: when the
-                                  // product of an object of this class
-                                  // is requested, we solve a linear
-                                  // equation system with that matrix
-                                  // using the CG method, accelerated
-                                  // by a preconditioner of (templated) class
-                                  // <code>Preconditioner</code>.
   template <class Matrix, class Preconditioner>
   class InverseMatrix : public Subscriptor
   {
@@ -301,103 +237,6 @@ namespace LinearSolvers
       }
   }
 
-                                  // @sect4{Schur complement preconditioner}
-
-                                  // This is the implementation
-                                  // of the Schur complement
-                                  // preconditioner as described
-                                  // in the section on improved
-                                  // solvers in step-22.
-                                  // 
-                                  // The basic 
-                                  // concept of the preconditioner is 
-                                  // different to the solution 
-                                  // strategy used in step-20 and 
-                                  // step-22. There, the Schur
-                                  // complement was used for a 
-                                  // two-stage solution of the linear
-                                  // system. Recall that the process
-                                  // in the Schur complement solver is
-                                  // a Gaussian elimination of
-                                  // a 2x2 block matrix, where each
-                                  // block is solved iteratively. 
-                                  // Here, the idea is to let 
-                                  // an iterative solver act on the
-                                  // whole system, and to use 
-                                  // a Schur complement for 
-                                  // preconditioning. As usual when
-                                  // dealing with preconditioners, we
-                                  // don't intend to exacly set up a 
-                                  // Schur complement, but rather use
-                                  // a good approximation to the
-                                  // Schur complement for the purpose of
-                                  // preconditioning.
-                                  // 
-                                  // So the question is how we can
-                                  // obtain a good preconditioner.
-                                  // Let's have a look at the 
-                                  // preconditioner matrix <i>P</i>
-                                  // acting on the block system, built
-                                  // as
-                                  // @f{eqnarray*}
-                                  //   P^{-1}
-                                  //   = 
-                                  //   \left(\begin{array}{cc}
-                                  //     A^{-1} & 0 \\ S^{-1} B A^{-1} & -S^{-1}
-                                  //   \end{array}\right)
-                                  // @f}
-                                  // using the Schur complement 
-                                  // $S = B A^{-1} B^T$. If we apply
-                                  // this matrix in the solution of 
-                                  // a linear system, convergence of
-                                  // an iterative Krylov-based solver
-                                  // will be governed by the matrix
-                                  // @f{eqnarray*}
-                                  //   P^{-1}\left(\begin{array}{cc}
-                                  //     A & B^T \\ B & 0
-                                  //   \end{array}\right) 
-                                  //  = 
-                                  //   \left(\begin{array}{cc}
-                                  //     I & A^{-1} B^T \\ 0 & 0
-                                  //   \end{array}\right),
-                                  // @f}
-                                  // which turns out to be very simple.
-                                  // A GMRES solver based on exact
-                                  // matrices would converge in two
-                                  // iterations, since there are
-                                  // only two distinct eigenvalues.
-                                  // Such a preconditioner for the
-                                  // blocked Stokes system has been 
-                                  // proposed by Silvester and Wathen,
-                                  // Fast iterative solution of 
-                                  // stabilised Stokes systems part II. 
-                                  // Using general block preconditioners.
-                                  // (SIAM J. Numer. Anal., 31 (1994),
-                                  // pp. 1352-1367).
-                                  // 
-                                  // The deal.II users who have already
-                                  // gone through the step-20 and step-22 
-                                  // tutorials can certainly imagine
-                                  // how we're going to implement this.
-                                  // We replace the inverse matrices
-                                  // in $P^{-1}$ using the InverseMatrix
-                                  // class, and the inverse Schur 
-                                  // complement will be approximated
-                                  // by the pressure mass matrix $M_p$.
-                                  // Having this in mind, we define a
-                                  // preconditioner class with a 
-                                  // <code>vmult</code> functionality,
-                                  // which is all we need for the
-                                  // interaction with the usual solver
-                                  // functions further below in the
-                                  // program code.
-                                  // 
-                                  // First the declarations. These
-                                  // are similar to the definition of
-                                  // the Schur complement in step-20,
-                                  // with the difference that we need
-                                  // some more preconditioners in
-                                  // the constructor.
   template <class PreconditionerA, class PreconditionerMp>
   class BlockSchurPreconditioner : public Subscriptor
   {
@@ -434,26 +273,6 @@ namespace LinearSolvers
                  tmp                     (stokes_matrix->block(1,1).row_map)
   {}
 
-
-                                  // This is the <code>vmult</code>
-                                  // function. We implement
-                                  // the action of $P^{-1}$ as described
-                                  // above in three successive steps.
-                                  // The first step multiplies
-                                  // the velocity vector by a 
-                                  // preconditioner of the matrix <i>A</i>.
-                                  // The resuling velocity vector
-                                  // is then multiplied by $B$ and
-                                  // subtracted from the pressure.
-                                  // This second step only acts on 
-                                  // the pressure vector and is 
-                                  // accomplished by the command
-                                  // SparseMatrix::residual. Next, 
-                                  // we change the sign in the 
-                                  // temporary pressure vector and
-                                  // finally multiply by the pressure
-                                  // mass matrix to get the final
-                                  // pressure vector.
   template <class PreconditionerA, class PreconditionerMp>
   void BlockSchurPreconditioner<PreconditionerA, PreconditionerMp>::vmult (
     TrilinosWrappers::BlockVector       &dst,
@@ -469,20 +288,6 @@ namespace LinearSolvers
 
 
                                 // @sect3{The <code>BoussinesqFlowProblem</code> class template}
-
-                                // The definition of this class is
-                                // mainly based on the step-22 tutorial
-                                // program. Most of the data types are
-                                // the same as there. However, we
-                                // deal with a time-dependent system now,
-                                // and there is temperature to take care
-                                // of as well, so we need some additional
-                                // function and variable declarations.
-                                // Furthermore, we have a slightly more
-                                // sophisticated solver we are going to
-                                // use, so there is a second pointer
-                                // to a sparse ILU for a pressure
-                                // mass matrix as well.
 template <int dim>
 class BoussinesqFlowProblem
 {
@@ -569,19 +374,6 @@ class BoussinesqFlowProblem
                                 // @sect3{BoussinesqFlowProblem class implementation}
 
                                 // @sect4{BoussinesqFlowProblem::BoussinesqFlowProblem}
-                                // 
-                                // The constructor of this class is
-                                // an extension of the constructor
-                                // in step-22. We need to include 
-                                // the temperature in the definition
-                                // of the finite element. As discussed
-                                // in the introduction, we are going 
-                                // to use discontinuous elements 
-                                // of one degree less than for pressure
-                                // there. Moreover, we initialize
-                                // the time stepping as well as the
-                                // options for the matrix assembly 
-                                // and preconditioning.
 template <int dim>
 BoussinesqFlowProblem<dim>::BoussinesqFlowProblem ()
                 :
@@ -754,46 +546,6 @@ compute_viscosity(const std::vector<double>          &old_temperature,
 
 
                                 // @sect4{BoussinesqFlowProblem::setup_dofs}
-                                // 
-                                // This function does the same as
-                                // in most other tutorial programs. 
-                                // As a slight difference, the 
-                                // program is called with a 
-                                // parameter <code>setup_matrices</code>
-                                // that decides whether to 
-                                // recreate the sparsity pattern
-                                // and the associated stiffness
-                                // matrix.
-                                // 
-                                // The body starts by assigning dofs on 
-                                // basis of the chosen finite element,
-                                // and then renumbers the dofs 
-                                // first using the Cuthill_McKee
-                                // algorithm (to generate a good
-                                // quality ILU during the linear
-                                // solution process) and then group
-                                // components of velocity, pressure
-                                // and temperature together. This 
-                                // happens in complete analogy to
-                                // step-22.
-                                // 
-                                // We then proceed with the generation
-                                // of the hanging node constraints
-                                // that arise from adaptive grid
-                                // refinement. Next we impose
-                                // the no-flux boundary conditions
-                                // $\vec{u}\cdot \vec{n}=0$ by adding
-                                // a respective constraint to the
-                                // hanging node constraints
-                                // matrix. The second parameter in 
-                                // the function describes the first 
-                                // of the velocity components
-                                // in the total dof vector, which is 
-                                // zero here. The parameter 
-                                // <code>no_normal_flux_boundaries</code>
-                                // sets the no flux b.c. to those
-                                // boundaries with boundary indicator
-                                // zero.
 template <int dim>
 void BoussinesqFlowProblem<dim>::setup_dofs ()
 {
@@ -845,67 +597,6 @@ void BoussinesqFlowProblem<dim>::setup_dofs ()
             << std::endl
             << std::endl;
 
-
-  
-                                // The next step is to 
-                                // create the sparsity 
-                                // pattern for the system matrix 
-                                // based on the Boussinesq 
-                                // system. As in step-22, 
-                                // we choose to create the
-                                // pattern not as in the
-                                // first tutorial programs,
-                                // but by using the blocked
-                                // version of 
-                                // CompressedSetSparsityPattern.
-                                // The reason for doing this 
-                                // is mainly a memory issue,
-                                // that is, the basic procedures
-                                // consume too much memory
-                                // when used in three spatial
-                                // dimensions as we intend
-                                // to do for this program.
-                                // 
-                                // So, in case we need
-                                // to recreate the matrices,
-                                // we first release the
-                                // stiffness matrix from the
-                                // sparsity pattern and then
-                                // set up an object of the 
-                                // BlockCompressedSetSparsityPattern
-                                // consisting of three blocks. 
-                                // Each of these blocks is
-                                // initialized with the
-                                // respective number of 
-                                // degrees of freedom. 
-                                // Once the blocks are 
-                                // created, the overall size
-                                // of the sparsity pattern
-                                // is initiated by invoking 
-                                // the <code>collect_sizes()</code>
-                                // command, and then the
-                                // sparsity pattern can be
-                                // filled with information.
-                                // Then, the hanging
-                                // node constraints are applied
-                                // to the temporary sparsity
-                                // pattern, which is finally
-                                // then completed and copied
-                                // into the general sparsity
-                                // pattern structure.
-  
-                                // Observe that we use a 
-                                // coupling argument for 
-                                // telling the function
-                                // <code>make_stokes_sparsity_pattern</code>
-                                // which components actually
-                                // will hold data and which 
-                                // we're going to neglect.
-                                // 
-                                // After these actions, we 
-                                // need to reassign the 
-                                // system matrix structure to
-                                // the sparsity pattern.
   stokes_partitioner.clear();
   {
     Epetra_Map map_u(n_u, 0, trilinos_communicator);
@@ -927,12 +618,6 @@ void BoussinesqFlowProblem<dim>::setup_dofs ()
 
     Table<2,DoFTools::Coupling> coupling (dim+1, dim+1);
 
-                                    // build the sparsity
-                                    // pattern. note that all dim
-                                    // velocities couple with each
-                                    // other and with the pressures,
-                                    // but that there is no
-                                    // pressure-pressure coupling:
     for (unsigned int c=0; c<dim+1; ++c)
       for (unsigned int d=0; d<dim+1; ++d)
        if (! ((c==dim) && (d==dim)))
@@ -1004,15 +689,6 @@ void BoussinesqFlowProblem<dim>::setup_dofs ()
                                         temperature_sparsity_pattern);
   }
 
-                                  // As last action in this function,
-                                  // we need to set the vectors
-                                  // for the solution, the old 
-                                  // solution (required for 
-                                  // time stepping) and the system
-                                  // right hand side to the 
-                                  // three-block structure given
-                                  // by velocity, pressure and
-                                  // temperature.
   stokes_solution.reinit (stokes_partitioner);
   stokes_rhs.reinit (stokes_partitioner);
 
@@ -1092,27 +768,6 @@ BoussinesqFlowProblem<dim>::build_stokes_preconditioner ()
   
   std::cout << "   Rebuilding Stokes preconditioner..." << std::flush;
       
-
-                                  // This last step of the assembly
-                                  // function sets up the preconditioners
-                                  // used for the solution of the
-                                  // system. We are going to use an
-                                  // ILU preconditioner for the
-                                  // velocity block (to be used
-                                  // by BlockSchurPreconditioner class)
-                                  // as well as an ILU preconditioner
-                                  // for the inversion of the 
-                                  // pressure mass matrix. Recall that
-                                  // the velocity-velocity block sits
-                                  // at position (0,0) in the 
-                                  // global system matrix, and
-                                  // the pressure mass matrix in
-                                  // (1,1). The 
-                                  // storage of these objects is
-                                  // as in step-22, that is, we
-                                  // include them using a 
-                                  // shared pointer structure from the
-                                  // boost library.
   assemble_stokes_preconditioner ();
       
   Amg_preconditioner = boost::shared_ptr<TrilinosWrappers::PreconditionAMG>
@@ -1126,11 +781,6 @@ BoussinesqFlowProblem<dim>::build_stokes_preconditioner ()
   Amg_preconditioner->initialize(stokes_preconditioner_matrix.block(0,0),
                                 true, true, null_space, false);
 
-                                  // TODO: we could throw away the (0,0)
-                                  // block here since things have been
-                                  // copied over to Trilinos. we need to
-                                  // keep the (1,1) block, though
-      
   Mp_preconditioner = boost::shared_ptr<TrilinosWrappers::PreconditionSSOR>
        (new TrilinosWrappers::PreconditionSSOR(
           stokes_preconditioner_matrix.block(1,1),1.2));
@@ -1143,73 +793,6 @@ BoussinesqFlowProblem<dim>::build_stokes_preconditioner ()
 
 
                                 // @sect4{BoussinesqFlowProblem::assemble_stokes_system}
-                                // 
-                                // The assembly of the Boussinesq 
-                                // system is acutally a two-step
-                                // procedure. One is to create
-                                // the Stokes system matrix and
-                                // right hand side for the 
-                                // velocity-pressure system as
-                                // well as the mass matrix for
-                                // temperature, and
-                                // the second is to create the
-                                // rhight hand side for the temperature
-                                // dofs. The reason for doing this
-                                // in two steps is simply that 
-                                // the time stepping we have chosen
-                                // needs the result from the Stokes
-                                // system at the current time step
-                                // for building the right hand
-                                // side of the temperature equation.
-                                // 
-                                // This function does the 
-                                // first of these two tasks.
-                                // There are two different situations
-                                // for calling this function. The
-                                // first one is when we reset the
-                                // mesh, and both the matrix and
-                                // the right hand side have to
-                                // be generated. The second situation
-                                // only sets up the right hand
-                                // side. The reason for having 
-                                // two different accesses is that
-                                // the matrix of the Stokes system
-                                // does not change in time unless
-                                // the mesh is changed, so we can
-                                // save a considerable amount of
-                                // work by doing the full assembly
-                                // only when it is needed.
-                                // 
-                                // Regarding the technical details
-                                // of implementation, not much has
-                                // changed from step-22. We reset
-                                // matrix and vector, create 
-                                // a quadrature formula on the 
-                                // cells and one on cell faces
-                                // (for implementing Neumann 
-                                // boundary conditions). Then,
-                                // we create a respective
-                                // FEValues object for both the 
-                                // cell and the face integration.
-                                // For the the update flags of
-                                // the first, we perform the
-                                // calculations of basis function
-                                // derivatives only in
-                                // case of a full assembly, since
-                                // they are not needed otherwise,
-                                // which makes the call of
-                                // the FEValues::reinit function
-                                // further down in the program 
-                                // more efficient.
-                                // 
-                                // The declarations proceed 
-                                // with some shortcuts for 
-                                // array sizes, the creation of
-                                // the local matrix and right 
-                                // hand side as well as the
-                                // vector for the indices of
-                                // the local dofs compared to
-                                // the global system.
 template <int dim>
 void BoussinesqFlowProblem<dim>::assemble_stokes_system ()
 {
@@ -1252,38 +835,6 @@ void BoussinesqFlowProblem<dim>::assemble_stokes_system ()
 
   std::vector<unsigned int> local_dof_indices (dofs_per_cell);
 
-                                  // These few declarations provide
-                                  // the structures for the evaluation
-                                  // of inhomogeneous Neumann boundary
-                                  // conditions from the function
-                                  // declaration made above.
-                                  // The vector <code>old_solution_values</code>
-                                  // evaluates the solution 
-                                  // at the old time level, since
-                                  // the temperature from the
-                                  // old time level enters the 
-                                  // Stokes system as a source
-                                  // term in the momentum equation.
-                                  // 
-                                  // Then, we create a variable
-                                  // to hold the Rayleigh number,
-                                  // the measure of buoyancy.
-                                  // 
-                                  // The set of vectors we create
-                                  // next hold the evaluations of
-                                  // the basis functions that will
-                                  // be used for creating the
-                                  // matrices. This gives faster
-                                  // access to that data, which
-                                  // increases the performance
-                                  // of the assembly. See step-22 
-                                  // for details.
-                                  // 
-                                  // The last few declarations 
-                                  // are used to extract the 
-                                  // individual blocks (velocity,
-                                  // pressure, temperature) from
-                                  // the total FE system.
   const EquationData::PressureBoundaryValues<dim> pressure_boundary_values;
   std::vector<double>               boundary_values (n_face_q_points);
 
@@ -1299,19 +850,6 @@ void BoussinesqFlowProblem<dim>::assemble_stokes_system ()
   const FEValuesExtractors::Vector velocities (0);
   const FEValuesExtractors::Scalar pressure (dim);
 
-                                  // Now start the loop over
-                                  // all cells in the problem.
-                                  // The first commands are all
-                                  // very familiar, doing the
-                                  // evaluations of the element
-                                  // basis functions, resetting
-                                  // the local arrays and 
-                                  // getting the values of the
-                                  // old solution at the
-                                  // quadrature point. Then we
-                                  // are ready to loop over
-                                  // the quadrature points 
-                                  // on the cell.
   typename DoFHandler<dim>::active_cell_iterator
     cell = stokes_dof_handler.begin_active(),
     endc = stokes_dof_handler.end();
@@ -1332,31 +870,6 @@ void BoussinesqFlowProblem<dim>::assemble_stokes_system ()
        {
          const double old_temperature = old_temperature_values[q];
 
-                                          // Extract the basis relevant
-                                          // terms in the inner products
-                                          // once in advance as shown
-                                          // in step-22 in order to 
-                                          // accelerate assembly.
-                                          // 
-                                          // Once this is done, we 
-                                          // start the loop over the
-                                          // rows and columns of the
-                                          // local matrix and feed
-                                          // the matrix with the relevant
-                                          // products. The right hand
-                                          // side is filled with the 
-                                          // forcing term driven by
-                                          // temperature in direction
-                                          // of gravity (which is 
-                                          // vertical in our example).
-                                          // Note that the right hand 
-                                          // side term is always generated,
-                                          // whereas the matrix 
-                                          // contributions are only
-                                          // updated when it is 
-                                          // requested by the
-                                          // <code>rebuild_matrices</code>
-                                          // flag.
          for (unsigned int k=0; k<dofs_per_cell; ++k)
            {
              phi_u[k] = stokes_fe_values[velocities].value (k,q);
@@ -1384,20 +897,6 @@ void BoussinesqFlowProblem<dim>::assemble_stokes_system ()
                             gravity * phi_u[i] * old_temperature)*
                            stokes_fe_values.JxW(q);
        }
-
-
-                                      // Next follows the assembly 
-                                      // of the face terms, result
-                                      // from Neumann boundary 
-                                      // conditions. Since these
-                                      // terms only enter the right
-                                      // hand side vector and not
-                                      // the matrix, there is no
-                                      // substantial benefit from
-                                      // extracting the data 
-                                      // before using it, so 
-                                      // we remain in the lines 
-                                      // of step-20 at this point.
       for (unsigned int face_no=0;
           face_no<GeometryInfo<dim>::faces_per_cell;
           ++face_no)
@@ -1422,16 +921,6 @@ void BoussinesqFlowProblem<dim>::assemble_stokes_system ()
                }
          }      
 
-                                      // The last step in the loop 
-                                      // over all cells is to
-                                      // enter the local contributions
-                                      // into the global matrix and 
-                                      // vector structures to the
-                                      // positions specified in 
-                                      // <code>local_dof_indices</code>.
-                                      // Again, we only add the 
-                                      // matrix data when it is 
-                                      // requested.
       cell->get_dof_indices (local_dof_indices);
 
       if (rebuild_stokes_matrix == true)
@@ -1457,28 +946,6 @@ void BoussinesqFlowProblem<dim>::assemble_stokes_system ()
 
 
                                 // @sect4{BoussinesqFlowProblem::assemble_temperature_system}
-                                // 
-                                // This function does the second
-                                // part of the assembly work, the
-                                // creation of the velocity-dependent
-                                // right hand side of the
-                                // temperature equation. The 
-                                // declarations in this function
-                                // are pretty much the same as the
-                                // ones used in the other 
-                                // assembly routine, except that we
-                                // restrict ourselves to vectors
-                                // this time. Though, we need to
-                                // perform more face integrals 
-                                // at this point, induced by the
-                                // use of discontinuous elements for 
-                                // the temperature (just
-                                // as it was in the first DG 
-                                // example in step-12) in combination
-                                // with adaptive grid refinement
-                                // and subfaces. The update 
-                                // flags at face level are the 
-                                // same as in step-12.
 template <int dim>
 void BoussinesqFlowProblem<dim>::assemble_temperature_matrix ()
 {
@@ -1506,15 +973,6 @@ void BoussinesqFlowProblem<dim>::assemble_temperature_matrix ()
   std::vector<double>                  phi_T       (dofs_per_cell);
   std::vector<Tensor<1,dim> >          grad_phi_T  (dofs_per_cell);
 
-                                  // Now, let's start the loop
-                                  // over all cells in the
-                                  // triangulation. The first
-                                  // actions within the loop
-                                  // are, 0as usual, the evaluation
-                                  // of the FE basis functions 
-                                  // and the old and present
-                                  // solution at the quadrature 
-                                  // points.
   typename DoFHandler<dim>::active_cell_iterator
     cell = temperature_dof_handler.begin_active(),
     endc = temperature_dof_handler.end();
@@ -1599,18 +1057,6 @@ void BoussinesqFlowProblem<dim>::assemble_temperature_system ()
 
   std::vector<unsigned int> local_dof_indices (dofs_per_cell);
 
-                                  // Here comes the declaration
-                                  // of vectors to hold the old
-                                  // and present solution values
-                                  // and gradients
-                                  // for both the cell as well as faces
-                                  // to the cell. Next comes the
-                                  // declaration of an object
-                                  // to hold the temperature 
-                                  // boundary values and a
-                                  // well-known extractor for
-                                  // accessing the temperature
-                                  // part of the FE system.
   std::vector<Vector<double> > present_stokes_values (n_q_points, 
                                                      Vector<double>(dim+1));
 
@@ -1634,15 +1080,6 @@ void BoussinesqFlowProblem<dim>::assemble_temperature_system ()
     global_T_range = get_extrapolated_temperature_range();
   const double global_Omega_diameter = GridTools::diameter (triangulation);
 
-                                  // Now, let's start the loop
-                                  // over all cells in the
-                                  // triangulation. The first
-                                  // actions within the loop
-                                  // are, 0as usual, the evaluation
-                                  // of the FE basis functions 
-                                  // and the old and present
-                                  // solution at the quadrature 
-                                  // points.
   typename DoFHandler<dim>::active_cell_iterator
     cell = temperature_dof_handler.begin_active(),
     endc = temperature_dof_handler.end();
@@ -1773,15 +1210,8 @@ template <int dim>
 void BoussinesqFlowProblem<dim>::solve ()
 {
   std::cout << "   Solving..." << std::endl;
-  
-                                  // Use the BlockMatrixArray structure
-                                  // for extracting only the upper left
-                                  // 2x2 blocks from the matrix that will
-                                  // be used for the solution of the
-                                  // blocked system.
+
   {
-                                    // Set up inverse matrix for
-                                    // pressure mass matrix
     LinearSolvers::InverseMatrix<TrilinosWrappers::SparseMatrix,
                                 TrilinosWrappers::PreconditionSSOR>
       mp_inverse (stokes_preconditioner_matrix.block(1,1), *Mp_preconditioner);
@@ -1790,8 +1220,6 @@ void BoussinesqFlowProblem<dim>::solve ()
                                             TrilinosWrappers::PreconditionSSOR>
       preconditioner (stokes_matrix, mp_inverse, *Amg_preconditioner);
 
-                                    // Set up GMRES solver and
-                                    // solve.
     SolverControl solver_control (stokes_matrix.m(),
                                  1e-6*stokes_rhs.l2_norm());
 
@@ -1806,13 +1234,6 @@ void BoussinesqFlowProblem<dim>::solve ()
               << " GMRES iterations for Stokes subsystem."
               << std::endl;
              
-                                    // Produce a constistent solution
-                                    // field (we can't do this on the 'up'
-                                    // vector since it does not have the
-                                    // temperature component, but
-                                    // hanging_node_constraints has
-                                    // constraints also for the
-                                    // temperature vector)
     stokes_constraints.distribute (stokes_solution);
   }
 
@@ -1838,7 +1259,6 @@ void BoussinesqFlowProblem<dim>::solve ()
              temperature_rhs,
              preconditioner);
 
-                                    // produce a consistent temperature field
     temperature_constraints.distribute (temperature_solution);
 
     std::cout << "   "

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.