<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 — 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}$ — 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$
-— 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> — 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 — 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}$ — the
-perfect preconditioner — 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 —
-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 — 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.
<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 β </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 — 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)}}$ — 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
-
/* $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 */
// @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
// @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>
{
// @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
{
}
}
- // @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
{
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,
// @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
{
// @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 ()
:
// @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 ()
{
<< 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);
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)))
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);
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>
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));
// @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 ()
{
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);
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();
{
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);
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)
}
}
- // 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)
// @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 ()
{
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();
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));
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();
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);
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());
<< " 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);
}
temperature_rhs,
preconditioner);
- // produce a consistent temperature field
temperature_constraints.distribute (temperature_solution);
std::cout << " "