From 66e90392f76b3f65f706698d3cbf26ad5bd7af96 Mon Sep 17 00:00:00 2001 From: bangerth Date: Wed, 3 Sep 2008 12:58:32 +0000 Subject: [PATCH] Remove all comments and other documentation copied from step-31. Only new stuff should be commented on. git-svn-id: https://svn.dealii.org/trunk@16729 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/examples/step-32/doc/intro.dox | 898 ----------------------- deal.II/examples/step-32/doc/results.dox | 272 ------- deal.II/examples/step-32/step-31.cc | 584 +-------------- 3 files changed, 2 insertions(+), 1752 deletions(-) diff --git a/deal.II/examples/step-32/doc/intro.dox b/deal.II/examples/step-32/doc/intro.dox index 33ce876a6f..acaa3cb3b7 100644 --- a/deal.II/examples/step-32/doc/intro.dox +++ b/deal.II/examples/step-32/doc/intro.dox @@ -14,901 +14,3 @@ California Institute of Technology.

Introduction

- -

The Boussinesq equations

- -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, u is the velocity field, p the pressure, and T -the temperature of the fluid. $\varepsilon ({\mathbf u}) = \frac 12 -[(\nabla{\mathbf u}) + (\nabla {\mathbf u})^T]$ is the symmetric -gradient of the velocity. As can be seen, velocity and pressure -solve a Stokes equation describing the motion of an incompressible -fluid, an equation we have previously considered in @ref step_22 "step-22"; we -will draw extensively on the experience we have gained in that program, in -particular with regard to efficient linear Stokes solvers. - -The forcing term of the fluid motion is the buoyancy of the -fluid, expressed as the product of the Rayleigh number $\mathrm{Ra}$, -the temperature T and the gravity vector g. (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 T); -sometimes, as in the case of rock minerals at temperatures close to their -melting point, $\eta$ may change by orders of magnitude over the typical range -of temperatures. - -$\mathrm{Ra}$, called the Rayleigh -number, 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 u 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 T 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 Navier-Stokes -equations hard to solve when the Reynolds number -$\mathrm{Re}$ 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. - - - -

%Boundary and initial conditions

- -Since the Boussinesq equations are derived under the assumption that inertia -of the fluid's motion does not play a role, the flow field is at each time -entirely determined by buoyancy difference at that time, not by the flow field -at previous times. This is reflected by the fact that the first two equations -above are the steady state Stokes equation that do not contain a time -derivative. Consequently, we do not need initial conditions for either -velocities or pressure. On the other hand, the temperature field does satisfy -an equation with a time derivative, so we need initial conditions for T. - -As for boundary conditions: if $\kappa>0$ then the temperature -satisfies a second order differential equation that requires -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 u=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. - - -

Solution approach

- -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 u and p 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. - - -

Time stepping

- -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 n 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 Backward -Differentiation Formula scheme of order 2 (BDF-2 in short) 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 -k 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 -Courant-Friedrichs-Lewy (CFL) condition 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 results section, but for the moment of importance is that -this CFL condition means that the time step -size k 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. - - -

Weak form and space discretization for the Stokes part

- -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. - - -

Stabilization, weak form and space discretization for the temperature equation

- -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 T 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 results section) 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 T — 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 T. 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 h, -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 n (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 k 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 k: -@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 k 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. - - - -

Linear solvers

- -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. - - -
Linear solvers for the Stokes problem
- -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 S 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 S 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 A 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 dim 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 x-components of the velocity are numbered, then -the y-components, and then the z-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 A 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 A 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 u=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 A. 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: - - - - -
Linear solvers for the temperature equation
- -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. - - - -

Implementation details

- -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 everything -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 downsides 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 update_hessians 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 all 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: - -In essence this means that all bottlenecks apart from the algebraic -multigrid have been removed. diff --git a/deal.II/examples/step-32/doc/results.dox b/deal.II/examples/step-32/doc/results.dox index 6a2f055b1a..a79c9a63f4 100644 --- a/deal.II/examples/step-32/doc/results.dox +++ b/deal.II/examples/step-32/doc/results.dox @@ -1,275 +1,3 @@

Results

- -

Numerical experiments to determine optimal parameters

- -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: - -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. - - -

Choosing ck and β

- -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. - - -
Results for Q1 elements
- -Here is what we get for -$\beta=0.01, \beta=0.1$, and $\beta=0.5$, different choices of $c_k$, and -bilinear elements (temperature_degree=1) in 2d: - - - - - - - - - - - -
- @image html "step-33.timestep.q1.beta=0.01.png" "" width=4cm - - @image html "step-33.timestep.q1.beta=0.03.png" "" width=4cm -
- @image html "step-33.timestep.q1.beta=0.1.png" "" width=4cm - - @image html "step-33.timestep.q1.beta=0.5.png" "" width=4cm -
- -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. - - -
Results for Q2 elements
- -One can repeat the same sequence of experiments for higher order -elements as well. Here are the graphs for bi-quadratic shape functions -(temperature_degree=2) for the temperature, while we -retain the $Q_2/Q_1$ stable Taylor-Hood element for the Stokes system: - - - - - - - - - - -
- @image html "step-33.timestep.q2.beta=0.01.png" "" width=4cm - - @image html "step-33.timestep.q2.beta=0.03.png" "" width=4cm -
- @image html "step-33.timestep.q2.beta=0.1.png" "" width=4cm -
- -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 grid point away (which for $Q_1$ elements -is $h_K$, but for $Q_2$ elements is $h_K/2$). - - -
Results for 3d
- -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 diameter 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. - - - - -
Conclusions
- -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 temperature_degree). 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. - - -

Possible extensions

- -Parallelization -> step-33 - diff --git a/deal.II/examples/step-32/step-31.cc b/deal.II/examples/step-32/step-31.cc index 2936194255..7e8da2e2d7 100644 --- a/deal.II/examples/step-32/step-31.cc +++ b/deal.II/examples/step-32/step-31.cc @@ -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 p and temperature - // T. - - // Secondly, we set an initial - // condition for all problem variables, - // i.e., for u, p and T, - // so the function has dim+2 - // 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 class TemperatureRightHandSide : public Function { @@ -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 InverseMatrix class template} - - // This class is an interface to - // calculate the action of an - // "inverted" matrix on a vector - // (using the vmult - // 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 - // Preconditioner. template 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 P - // 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 - // vmult 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 BlockSchurPreconditioner : public Subscriptor { @@ -434,26 +273,6 @@ namespace LinearSolvers tmp (stokes_matrix->block(1,1).row_map) {} - - // This is the vmult - // 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 A. - // 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 void BlockSchurPreconditioner::vmult ( TrilinosWrappers::BlockVector &dst, @@ -469,20 +288,6 @@ namespace LinearSolvers // @sect3{The BoussinesqFlowProblem 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 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 BoussinesqFlowProblem::BoussinesqFlowProblem () : @@ -754,46 +546,6 @@ compute_viscosity(const std::vector &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 setup_matrices - // 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 - // no_normal_flux_boundaries - // sets the no flux b.c. to those - // boundaries with boundary indicator - // zero. template void BoussinesqFlowProblem::setup_dofs () { @@ -845,67 +597,6 @@ void BoussinesqFlowProblem::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 collect_sizes() - // 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 - // make_stokes_sparsity_pattern - // 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::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::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::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 @@ -1126,11 +781,6 @@ BoussinesqFlowProblem::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 (new TrilinosWrappers::PreconditionSSOR( stokes_preconditioner_matrix.block(1,1),1.2)); @@ -1143,73 +793,6 @@ BoussinesqFlowProblem::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 void BoussinesqFlowProblem::assemble_stokes_system () { @@ -1252,38 +835,6 @@ void BoussinesqFlowProblem::assemble_stokes_system () std::vector 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 old_solution_values - // 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 pressure_boundary_values; std::vector boundary_values (n_face_q_points); @@ -1299,19 +850,6 @@ void BoussinesqFlowProblem::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::active_cell_iterator cell = stokes_dof_handler.begin_active(), endc = stokes_dof_handler.end(); @@ -1332,31 +870,6 @@ void BoussinesqFlowProblem::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 - // rebuild_matrices - // flag. for (unsigned int k=0; k::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::faces_per_cell; ++face_no) @@ -1422,16 +921,6 @@ void BoussinesqFlowProblem::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 - // local_dof_indices. - // 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::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 void BoussinesqFlowProblem::assemble_temperature_matrix () { @@ -1506,15 +973,6 @@ void BoussinesqFlowProblem::assemble_temperature_matrix () std::vector phi_T (dofs_per_cell); std::vector > 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::active_cell_iterator cell = temperature_dof_handler.begin_active(), endc = temperature_dof_handler.end(); @@ -1599,18 +1057,6 @@ void BoussinesqFlowProblem::assemble_temperature_system () std::vector 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 > present_stokes_values (n_q_points, Vector(dim+1)); @@ -1634,15 +1080,6 @@ void BoussinesqFlowProblem::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::active_cell_iterator cell = temperature_dof_handler.begin_active(), endc = temperature_dof_handler.end(); @@ -1773,15 +1210,8 @@ template void BoussinesqFlowProblem::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 mp_inverse (stokes_preconditioner_matrix.block(1,1), *Mp_preconditioner); @@ -1790,8 +1220,6 @@ void BoussinesqFlowProblem::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::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::solve () temperature_rhs, preconditioner); - // produce a consistent temperature field temperature_constraints.distribute (temperature_solution); std::cout << " " -- 2.39.5