From: bangerth Date: Mon, 11 Aug 2008 23:02:04 +0000 (+0000) Subject: Document the stabilization approach we are going to take. X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=489cba2ce2b65d12ef3388739910b928bfe5ab19;p=dealii-svn.git Document the stabilization approach we are going to take. git-svn-id: https://svn.dealii.org/trunk@16511 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/examples/step-31/doc/intro.dox b/deal.II/examples/step-31/doc/intro.dox index a898a0fa67..921046597c 100644 --- a/deal.II/examples/step-31/doc/intro.dox +++ b/deal.II/examples/step-31/doc/intro.dox @@ -176,7 +176,7 @@ in the top-left corner of the differential operator. 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 time lag scheme: first solve the Stokes equations for velocity and +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 @@ -187,25 +187,43 @@ system \\ \nabla \cdot {\mathbf u}^n &=& 0, @f} -and then +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 +term $\frac{\frac 32 T^{n}-2T^{n-1}+\frac 12 T^{n-2}}{k_n}$ where +$k_n$ is the current time step. @f{eqnarray*} - T^n + \frac 32 T^n - k_n\nabla \cdot \kappa \nabla T^n &=& - T^{n-1} - + - k_n\gamma + 2 T^{n-1} + - + \frac 12 T^{n-2} - - k_n{\mathbf u}^n \cdot \nabla T^{n-1}, + k_n{\mathbf u}^n \cdot \nabla T^{n-1} + + + k_n\gamma. @f} -where $k_n$ is the time step size. Note how the temperature equation is +Note how the temperature equation is solved semi-explicitly: diffusion is treated implicitly whereas advection is treated explicitly using the just-computed velocity -field. +field. We will have to modify the equation some more (by adding +stabilization terms) later on, but the general principle remains the +same. 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

+

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 @@ -245,36 +263,233 @@ that we implement is (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 we do with the temperature -equation. By default, not all discretizations of this equation are -equally stable unless we either do something like upwinding, -stabilization, or all of this. XXXXXXXXXX TALK ABOUT CHOICE OF FE -AND DG/CG XXXXX +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*} - (\tau_h,T^n_h) + \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) 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) = + \|\mathbf{u}\|_{L^\infty(\Omega)} \|T\|_{L^\infty(\Omega)} + |\textrm{diam}(\Omega)|^{\alpha-2}$. +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 also show that it produces better convergence orders +for the linear transport equation than for example streamline +diffusion. For $\alpha=2$, no theoretical results are available +currently, but numerical tests indicate that it produces results that +are much better than for $\alpha=1$. + +A more practical question is how to introduce this artificial +diffusion into the equations we would like to solve. To this end note +that the equations now have the form +@f{eqnarray*} + \frac{\partial T}{\partial t} + - k_n(\nabla \tau_h, \kappa \nabla T^n_h) + {\mathbf u} \cdot \nabla T + - + \nabla \cdot (\kappa+\nu(T)) \nabla T &=& \gamma, +@f} +and are consequently 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. Using the BDF-2 scheme introduced above, +this yields +@f{eqnarray*} + \frac 32 T^n + - + k_n\nabla \cdot \kappa \nabla T^n &=& - (\tau_h,T^{n-1}_h+k_n\gamma) + 2 T^{n-1} - - k_n (\tau_h, {\mathbf u}^n \cdot \nabla T^{n-1}_h), + \frac 12 T^{n-2} + \\ + && + + + k_n\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_n{\mathbf u}^n \cdot \nabla (2T^{n-1}-T^{n-2}) + \\ + && + + + k_n\gamma. @f} -for all discrete test functions $\tau_h$. One typically integrated the -last term on the right hand side by parts on each cell, yielding the -following taking into account that the velocity field is assumed to be -divergence-free: +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_n$ 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 where we calculate the +viscosity $\nu_\alpha$ using the temperature field at time $t_{n-\frac +32}$ using $\frac 12 (T^{n-1}+T^{n-2})$; the third line contains the +advection term, and the fourth the sources. Note that for +artificial diffusion and advection operate on an extrapolated +temperature at 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}$. 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. + +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*} - (\tau_h,T^n_h) + \frac 32 (\tau_h,T_h^n) + - k_n(\nabla \tau_h, \kappa \nabla T^n_h) + k_n (\nabla \tau_h, \kappa \nabla T_h^n) &=& - (\tau_h,T^{n-1}_h+k_n\gamma) + \left(\tau_h,2 T_h^{n-1} + - + \frac 12 T_h^{n-2} + - + k_n{\mathbf u}_h^n \cdot \nabla (2T_h^{n-1}-T_h^{n-2}) + - k_n \sum_K \left\{ ({\mathbf u}^n \cdot \nabla \tau_h, T^{n-1}_h)_K - - (\tau_h, (\mathbf n \cdot {\mathbf u}^n) T^{n-1}_h)_{\partial K} \right\}. + k_n\gamma \right) + \\ + && + - + k_n \left(\nabla \tau_h, + \nu_\alpha\left(\frac 12 T_h^{n-1}+\frac 12 T_h^{n-2}\right) + \ \nabla (2T_h^{n-1}-T_h^{n-2}) + \right) +@f} +for all discrete test functions $\tau_h$. This then results in a +matrix equation of form +@f{eqnarray*} + \left(\frac 32 M+k_n A\right) T_h^n = F, @f} -XXXXXXXXXX WHAT DO I DO WITH THE DIV U TERM? XXXXX +which given the structure of matrix on the left (the sum of two +positive definite matrices) is easily solved using the Conjugate +Gradient method.