<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
+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: we first solve the temperature equation (using an extrapolated
+velocity field), and then insert the new temperature solution into the right
+hand side of the velocity equation. The way we implement this in our code
+looks at things from a slightly different perspective, though. We first
+solve the Stokes equations for velocity and pressure using the temperature
+field from the previous time step, which means that we get the velocity for
+the previous time step. In other words, we first solve the Stokes system for
+time step <i>n-1</i> as
@f{eqnarray*}
- -\nabla \cdot \eta \varepsilon ({\mathbf u}^n) + \nabla p^n &=&
+ -\nabla \cdot \eta \varepsilon ({\mathbf u}^{n-1}) + \nabla p^{n-1} &=&
\mathrm{Ra} \; T^{n-1} \mathbf{g},
\\
- \nabla \cdot {\mathbf u}^n &=& 0,
+ \nabla \cdot {\mathbf u}^{n-1} &=& 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.
+and then the temperature equation with an extrapolated velocity field to
+time <i>n</i>.
-This gives the discretized-in-time temperature equation
+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
-
-
\frac 12 T^{n-2}
-
- k{\mathbf u}^n \cdot \nabla (2T^{n-1}-T^{n-2})
+ k(2{\mathbf u}^{n-1} - {\mathbf u}^{n-2} ) \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
+Note how the temperature equation is solved semi-explicitly: diffusion is
+treated implicitly whereas advection is treated explicitly using an
+extrapolation (or forward-projection) of temperature and velocity, including
+the just-computed velocity ${\mathbf u}^{n-1}$. The forward-projection to
+the current time level <i>n</i> is derived from a Taylor expansion, $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 second 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{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{\mathbf u}^{*,n} \cdot \nabla T^{*,n}
+
- k_n\gamma.
+ 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}$.
+
+where ${(\cdot)}^{*,n} = \left(1+\frac{k_n}{k_{n-1}}\right)(\cdot)^{n-1} -
+\frac{k_n}{k_{n-1}}(\cdot)^{n-2}$ denotes the extrapolation of velocity
+<b>u</b> and temperature <i>T</i> to time level <i>n</i>, using the values
+at the two previous time steps. 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;
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 {\mathbf v}_h, \eta \varepsilon ({\mathbf u}^{n-1}_h))
-
- (\nabla \cdot {\mathbf v}_h, p^n_h)
+ (\nabla \cdot {\mathbf v}_h, p^{n-1}_h)
&=&
({\mathbf v}_h, \mathrm{Ra} \; T^{n-1}_h \mathbf{g}),
\\
- (q_h, \nabla \cdot {\mathbf u}^n_h) &=& 0,
+ (q_h, \nabla \cdot {\mathbf u}^{n-1}_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
+$(\nabla {\mathbf v}_h, \eta \varepsilon ({\mathbf u}^{n-1}_h))_\Omega
= \int_\Omega \sum_{i,j=1}^d [\nabla {\mathbf v}_h]_{ij}
- \eta [\varepsilon ({\mathbf u}^n_h)]_{ij}\, dx$.
+ \eta [\varepsilon ({\mathbf u}^{n-1}_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))
+ (\varepsilon({\mathbf v}_h), \eta \varepsilon ({\mathbf u}^{n-1}_h))
-
- (\nabla \cdot {\mathbf v}_h, p^n_h)
+ (\nabla \cdot {\mathbf v}_h, p^{n-1}_h)
&=&
({\mathbf v}_h, \mathrm{Ra} \; T^{n-1}_h \mathbf{g}),
\\
- (q_h, \nabla \cdot {\mathbf u}^n_h) &=& 0.
+ (q_h, \nabla \cdot {\mathbf u}^{n-1}_h) &=& 0.
@f}
This is exactly the same as what we already discussed in
\\
&&
-
- k{\mathbf u}^n \cdot \nabla (2T^{n-1}-T^{n-2})
+ k(2{\mathbf u}^{n-1}-{\mathbf u}^{n-2}) \cdot \nabla (2T^{n-1}-T^{n-2})
\\
&&
+
\\
&&
-
- k_n{\mathbf u}^n \cdot \nabla \left[
+ k_n
+ \left[
+ \left(1+\frac{k_n}{k_{n-1}}\right){\mathbf u}^{n-1} -
+ \frac{k_n}{k_{n-1}}{\mathbf u}^{n-2}
+ \right]
+ \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]
\\
-
\frac{k_n^2}{k_{n-1}(k_n+k_{n-1})} T_h^{n-2}
\\
- &&\qquad\qquad
+ &&\qquad
-
- k_n{\mathbf u}_h^n \cdot \nabla \left[
+ k_n
+ \left[
+ \left(1+\frac{k_n}{k_{n-1}}\right){\mathbf u}^{n-1} -
+ \frac{k_n}{k_{n-1}}{\mathbf u}^{n-2}
+ \right]
+ \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]
+
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}),
+ \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-1}, U_h^{n-2},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
A & B^T & 0 \\ B & 0 &0 \\ C & 0 & K
\end{array}\right)
\left(\begin{array}{ccc}
- U^n \\ P^n \\ T^n
+ U^{n-1} \\ P^{n-1} \\ 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})
+ F_U(T^{n-1}) \\ 0 \\ F_T(U^{n-1},U^{n-2},T^{n-1},T^{n-2})
\end{array}\right).
@f}
The problem with this is: We never use the whole matrix at the same time. In
// @sect4{Schur complement preconditioner}
- // This is the implementation of
- // the Schur complement
- // preconditioner as described in
- // detail in the introduction. As
- // opposed to step-20 and step-22,
- // we solve the block system
- // all-at-once using GMRES, and use
- // the Schur complement of the
- // block structured matrix to build
- // a good preconditioner instead.
+ // This is the implementation of the
+ // Schur complement preconditioner as
+ // described in detail in the
+ // introduction. As opposed to step-20
+ // and step-22, we solve the block system
+ // all-at-once using GMRES, and use the
+ // Schur complement of the block
+ // structured matrix to build a good
+ // preconditioner instead.
//
// Let's have a look at the ideal
// preconditioner matrix
- // $P=\left(\begin{array}{cc} A & 0 \\ B &
- // -S \end{array}\right)$
- // described in the introduction. If
- // we apply this matrix in the
- // solution of a linear system,
- // convergence of an iterative
- // GMRES solver will be
- // governed by the matrix
+ // $P=\left(\begin{array}{cc} A & 0 \\ B
+ // & -S \end{array}\right)$ described in
+ // the introduction. If we apply this
+ // matrix in the solution of a linear
+ // system, convergence of an iterative
+ // GMRES 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}
- //
+ // \end{array}\right),
+ // @f}
// which indeed is very simple. A
// GMRES solver based on exact
// matrices would converge in two
// SIAM J. Numer. Anal., 31 (1994),
// pp. 1352-1367).
//
- // Replacing <i>P</i> by
- // $\tilde{P}$ does not change the
- // situation dramatically. The
- // product $P^{-1} A$ will still be
- // close to a matrix with
- // eigenvalues 0 and 1, which lets
- // us hope to be able to get a
- // number of GMRES iterations that
- // does not depend on the problem
- // size.
+ // Replacing <i>P</i> by $\tilde{P}$ does
+ // not change the situation
+ // dramatically. The product $P^{-1} A$
+ // will still be close to a matrix with
+ // eigenvalues 0 and 1, which lets us
+ // hope to be able to get a number of
+ // GMRES iterations that does not depend
+ // on the problem size.
//
- // The deal.II users who have already gone
- // through the step-20 and step-22
+ // 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 exact inverse matrices in
// $P^{-1}$ by some approximate inverses
- // built from the InverseMatrix class, and
- // the inverse Schur complement will be
- // approximated by the pressure mass matrix
- // $M_p$ (weighted by $\eta^{-1}$ as
- // mentioned in the introduction). As
+ // built from the InverseMatrix class,
+ // and the inverse Schur complement will
+ // be approximated by the pressure mass
+ // matrix $M_p$ (weighted by $\eta^{-1}$
+ // as mentioned in the introduction). As
// pointed out in the results section of
// step-22, we can replace the exact
// inverse of <i>A</i> by just the
- // application of a preconditioner, in this
- // case on a vector Laplace matrix as was
- // explained in the introduction. This does
- // increase the number of (outer) GMRES
- // iterations, but is still significantly
- // cheaper than an exact inverse, which
- // would require between 20 and 35 CG
+ // application of a preconditioner, in
+ // this case on a vector Laplace matrix
+ // as was explained in the
+ // introduction. This does increase the
+ // number of (outer) GMRES iterations,
+ // but is still significantly cheaper
+ // than an exact inverse, which would
+ // require between 20 and 35 CG
// iterations for <em>each</em> outer
// solver step (using the AMG
// preconditioner).
//
- // Having the above explanations 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
+ // Having the above explanations 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 and that the
- // matrices we use here are built
- // upon Trilinos:
+ // 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 and
+ // that the matrices we use here are
+ // built upon Trilinos:
template <class PreconditionerA, class PreconditionerMp>
class BlockSchurPreconditioner : public Subscriptor
{
const std::vector<Tensor<1,dim> > &old_old_temperature_grads,
const std::vector<Tensor<2,dim> > &old_temperature_hessians,
const std::vector<Tensor<2,dim> > &old_old_temperature_hessians,
- const std::vector<Vector<double> > &present_stokes_values,
+ const std::vector<Vector<double> > &old_stokes_values,
+ const std::vector<Vector<double> > &old_old_stokes_values,
const std::vector<double> &gamma_values,
const double global_u_infty,
const double global_T_variation,
TrilinosWrappers::BlockSparseMatrix stokes_preconditioner_matrix;
TrilinosWrappers::BlockVector stokes_solution;
+ TrilinosWrappers::BlockVector old_stokes_solution;
TrilinosWrappers::BlockVector stokes_rhs;
// The last of the tool functions computes
// the artificial viscosity parameter
// $\nu|_K$ on a cell $K$ as a function of
- // the extrapolated temperature, its gradient
- // and Hessian (second derivatives), the
- // velocity, the right hand side $\gamma$ all
- // on the quadrature points of the current
- // cell, and various other parameters as
- // described in detail in the introduction.
+ // the extrapolated temperature, its
+ // gradient and Hessian (second
+ // derivatives), the velocity, the right
+ // hand side $\gamma$ all on the quadrature
+ // points of the current cell, and various
+ // other parameters as described in detail
+ // in the introduction.
//
- // There are some universal constants
- // worth mentioning here. First, we
- // need to fix $\beta$; we choose
- // $\beta=0.015\cdot dim$, a choice
- // discussed in detail in the results
- // section of this tutorial
- // program. The second is the
- // exponent $\alpha$; $\alpha=1$
- // appears to work fine for the
- // current program, even though some
- // additional benefit might be
+ // There are some universal constants worth
+ // mentioning here. First, we need to fix
+ // $\beta$; we choose $\beta=0.015\cdot
+ // dim$, a choice discussed in detail in
+ // the results section of this tutorial
+ // program. The second is the exponent
+ // $\alpha$; $\alpha=1$ appears to work
+ // fine for the current program, even
+ // though some additional benefit might be
// expected from chosing $\alpha =
- // 2$. Finally, there is one thing
- // that requires special casing: In
- // the first time step, the velocity
- // equals zero, and the formula for
- // $\nu|_K$ is not defined. In that
- // case, we return $\nu|_K=5\cdot
- // 10^3 \cdot h_K$, a choice
- // admittedly more motivated by
- // heuristics than anything else (it
- // is in the same order of magnitude,
- // however, as the value returned for
- // most cells on the second time
- // step).
+ // 2$. Finally, there is one thing that
+ // requires special casing: In the first
+ // time step, the velocity equals zero, and
+ // the formula for $\nu|_K$ is not
+ // defined. In that case, we return
+ // $\nu|_K=5\cdot 10^3 \cdot h_K$, a choice
+ // admittedly more motivated by heuristics
+ // than anything else (it is in the same
+ // order of magnitude, however, as the
+ // value returned for most cells on the
+ // second time step).
//
// The rest of the function should be
- // mostly obvious based on the
- // material discussed in the
- // introduction:
+ // mostly obvious based on the material
+ // discussed in the introduction:
template <int dim>
double
BoussinesqFlowProblem<dim>::
const std::vector<Tensor<1,dim> > &old_old_temperature_grads,
const std::vector<Tensor<2,dim> > &old_temperature_hessians,
const std::vector<Tensor<2,dim> > &old_old_temperature_hessians,
- const std::vector<Vector<double> > &present_stokes_values,
+ const std::vector<Vector<double> > &old_stokes_values,
+ const std::vector<Vector<double> > &old_old_stokes_values,
const std::vector<double> &gamma_values,
const double global_u_infty,
const double global_T_variation,
{
Tensor<1,dim> u;
for (unsigned int d=0; d<dim; ++d)
- u[d] = present_stokes_values[q](d);
+ u[d] = (old_stokes_values[q](d) + old_old_stokes_values[q](d)) / 2;
const double dT_dt = (old_temperature[q] - old_old_temperature[q])
/ old_time_step;
}
// Lastly, we set the vectors for the
- // solution $\mathbf u$ and $T^k$, the old
- // solutions $T^{k-1}$ and $T^{k-2}$
- // (required for time stepping) and the
- // system right hand sides to their correct
- // sizes and block structure:
+ // Stokes solutions $\mathbf u^{n-1}$ and
+ // $\mathbf u^{n-2}$, as well as for the
+ // temperatures $T^{n}$, $T^{n-1}$ and
+ // $T^{n-2}$ (required for time stepping)
+ // and all the system right hand sides to
+ // their correct sizes and block
+ // structure:
stokes_solution.reinit (stokes_block_sizes);
+ old_stokes_solution.reinit (stokes_block_sizes);
stokes_rhs.reinit (stokes_block_sizes);
temperature_solution.reinit (temperature_dof_handler.n_dofs());
std::vector<unsigned int> local_dof_indices (dofs_per_cell);
- // Next comes the declaration of vectors to
- // hold the old and present solution values
- // and gradients at quadrature points of
- // the current cell. We also declarate an
- // object to hold the temperature right
- // hande side values
- // (<code>gamma_values</code>), and we
- // again use shortcuts for the temperature
- // basis functions. Eventually, we need to
- // find the maximum of velocity,
- // temperature and the diameter of the
- // computational domain which will be used
- // for the definition of the stabilization
+ // Next comes the declaration of vectors
+ // to hold the old and older solution
+ // values (as a notation for time levels
+ // <i>n-1</i> and <i>n-2</i>,
+ // respectively) and gradients at
+ // quadrature points of the current
+ // cell. We also declarate an object to
+ // hold the temperature right hande side
+ // values (<code>gamma_values</code>),
+ // and we again use shortcuts for the
+ // temperature basis
+ // functions. Eventually, we need to find
+ // the maximum of velocity, temperature
+ // and the diameter of the computational
+ // domain which will be used for the
+ // definition of the stabilization
// parameter.
- std::vector<Vector<double> > present_stokes_values (n_q_points,
+ std::vector<Vector<double> > old_stokes_values (n_q_points,
+ Vector<double>(dim+1));
+ std::vector<Vector<double> > old_old_stokes_values (n_q_points,
Vector<double>(dim+1));
std::vector<double> old_temperature_values (n_q_points);
std::vector<double> old_old_temperature_values(n_q_points);
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. Again, we need two
- // cell iterators that walk in parallel
- // through the cells of the two involved
- // DoFHandler objects for the Stokes and
- // temperature part. Within the loop, we
- // first set the local rhs to zero, and
- // then get the values and derivatives of
- // the old solution functions (and the
- // current velocity) at the quadrature
- // points, since they are going to be
- // needed for the definition of the
- // stabilization parameters and as
- // coefficients in the equation,
- // respectively.
+ // Now, let's start the loop over all
+ // cells in the triangulation. Again, we
+ // need two cell iterators that walk in
+ // parallel through the cells of the two
+ // involved DoFHandler objects for the
+ // Stokes and temperature part. Within
+ // the loop, we first set the local rhs
+ // to zero, and then get the values and
+ // derivatives of the old solution
+ // functions at the quadrature points,
+ // since they are going to be needed for
+ // the definition of the stabilization
+ // parameters and as coefficients in the
+ // equation, respectively.
typename DoFHandler<dim>::active_cell_iterator
cell = temperature_dof_handler.begin_active(),
endc = temperature_dof_handler.end();
gamma_values);
stokes_fe_values.get_function_values (stokes_solution,
- present_stokes_values);
+ old_stokes_values);
+ stokes_fe_values.get_function_values (old_stokes_solution,
+ old_old_stokes_values);
// Next, we calculate the
// artificial viscosity for
old_old_temperature_grads,
old_temperature_hessians,
old_old_temperature_hessians,
- present_stokes_values,
+ old_stokes_values,
+ old_old_stokes_values,
gamma_values,
global_u_infty,
global_T_range.second - global_T_range.first,
const Tensor<1,dim> old_old_grad_T = old_old_temperature_grads[q];
- Tensor<1,dim> present_u;
+ Tensor<1,dim> extrapolated_u;
for (unsigned int d=0; d<dim; ++d)
- present_u[d] = present_stokes_values[q](d);
+ {
+ if (use_bdf2_scheme == true)
+ extrapolated_u[d] =
+ old_stokes_values[q](d) * (1+time_step/old_time_step) -
+ old_old_stokes_values[q](d) * time_step/old_time_step;
+ else
+ extrapolated_u[d] = old_stokes_values[q](d);
+ }
if (use_bdf2_scheme == true)
{
old_old_T * phi_T[i]
-
time_step *
- present_u *
+ extrapolated_u *
((1+time_step/old_time_step) * old_grad_T
-
time_step / old_time_step * old_old_grad_T) *
local_rhs(i) += (old_T * phi_T[i]
-
time_step *
- present_u * old_grad_T * phi_T[i]
+ extrapolated_u * old_grad_T * phi_T[i]
-
time_step *
nu *
// @sect4{BoussinesqFlowProblem::solve}
//
- // This function solves the linear
- // systems of equations. Following to
- // the introduction, we start with
- // the Stokes system, where we need
- // to generate our block Schur
- // preconditioner. Since all the
- // relevant actions are implemented
- // in the class
+ // This function solves the linear systems
+ // of equations. Following to the
+ // introduction, we start with the Stokes
+ // system, where we need to generate our
+ // block Schur preconditioner. Since all
+ // the relevant actions are implemented in
+ // the class
// <code>BlockSchurPreconditioner</code>,
- // all we have to do is to
- // initialize the class
- // appropriately. What we need to
+ // all we have to do is to initialize the
+ // class appropriately. What we need to
// pass down is an
- // <code>InverseMatrix</code> object
- // for the pressure mass matrix,
- // which we set up using the
- // respective class together with
- // the IC preconditioner we already
- // generated, and the AMG
- // preconditioner for the
- // velocity-velocity matrix. Note
- // that both
- // <code>Mp_preconditioner</code> and
- // <code>Amg_preconditioner</code> are
- // only pointers, so we use
- // <code>*</code> to pass down the
- // actual preconditioner objects.
+ // <code>InverseMatrix</code> object for
+ // the pressure mass matrix, which we set
+ // up using the respective class together
+ // with the IC preconditioner we already
+ // generated, and the AMG preconditioner
+ // for the velocity-velocity matrix. Note
+ // that both <code>Mp_preconditioner</code>
+ // and <code>Amg_preconditioner</code> are
+ // only pointers, so we use <code>*</code>
+ // to pass down the actual preconditioner
+ // objects.
//
- // Once the preconditioner is
- // ready, we create a GMRES solver
- // for the block system. Since we
- // are working with Trilinos data
- // structures, we have to set the
- // respective template argument in
- // the solver. GMRES needs to
- // internally store temporary
- // vectors for each iteration (see
- // the discussion in the
- // results section of step-22)
- // – the more vectors it can
- // use, the better it will
+ // Once the preconditioner is ready, we
+ // create a GMRES solver for the block
+ // system. Since we are working with
+ // Trilinos data structures, we have to set
+ // the respective template argument in the
+ // solver. GMRES needs to internally store
+ // temporary vectors for each iteration
+ // (see the discussion in the results
+ // section of step-22) – the more
+ // vectors it can use, the better it will
// generally perform. To keep memory
- // demands in check, we
- // set the number of vectors to
- // 100. This means that up to 100
- // solver iterations, every
- // temporary vector can be
- // stored. If the solver needs to
- // iterate more often to get the
- // specified tolerance, it will
- // work on a reduced set of vectors
- // by restarting at every 100
- // iterations.
+ // demands in check, we set the number of
+ // vectors to 100. This means that up to
+ // 100 solver iterations, every temporary
+ // vector can be stored. If the solver
+ // needs to iterate more often to get the
+ // specified tolerance, it will work on a
+ // reduced set of vectors by restarting at
+ // every 100 iterations.
//
// With this all set up, we solve the system
// and distribute the constraints in the
cell->clear_refine_flag ();
// Before we can apply the mesh
- // refinement, we have to prepare
- // the solution vectors that should
- // be transfered to the new grid
- // (we will lose the old grid once
- // we have done the
- // refinement). What we definetely
+ // refinement, we have to prepare the
+ // solution vectors that should be
+ // transfered to the new grid (we will
+ // lose the old grid once we have done
+ // the refinement). What we definetely
// need are the current and the old
// temperature (BDF-2 time stepping
- // requires two old
- // solutions). Since the
- // SolutionTransfer objects only
- // support to transfer one object
- // per dof handler, we need to
- // collect the two temperature
- // solutions in one data
- // structure. Moreover, we choose
- // to transfer the Stokes solution,
- // too. The reason for doing so is
- // that the Stokes solution will
- // not change dramatically from
- // step to step, so we get a good
- // initial guess for the linear
- // solver when we reuse old data,
- // which reduces the number of
- // needed solver iterations. Next,
- // we initialize the
- // SolutionTransfer objects, by
- // attaching them to the old dof
- // handler. With this at place, we
- // can prepare the triangulation
- // and the data vectors for
+ // requires two old solutions). Since the
+ // SolutionTransfer objects only support
+ // to transfer one object per dof
+ // handler, we need to collect the two
+ // temperature solutions in one data
+ // structure. Moreover, we choose to
+ // transfer the Stokes solution, too. The
+ // reason for doing so is that the Stokes
+ // solution will not change dramatically
+ // from step to step, so we get a good
+ // initial guess for the linear solver
+ // when we reuse old data, which reduces
+ // the number of needed solver
+ // iterations. Next, we initialize the
+ // SolutionTransfer objects, by attaching
+ // them to the old dof handler. With this
+ // at place, we can prepare the
+ // triangulation and the data vectors for
// refinement (in this order).
std::vector<TrilinosWrappers::Vector> x_temperature (2);
x_temperature[0].reinit (temperature_solution);
// change in case we've remeshed
// before), and then do the
// solve. The solution is then
- // written to screen. Before going
- // on with the next time step, we
- // have to check whether we should
- // first finish the pre-refinement
- // steps or if we should remesh
- // (every fifth time step),
- // refining up to a level that is
- // consistent with initial
+ // written to screen. Before going on
+ // with the next time step, we have
+ // to check whether we should first
+ // finish the pre-refinement steps or
+ // if we should remesh (every fifth
+ // time step), refining up to a level
+ // that is consistent with initial
// refinement and pre-refinement
// steps. Last in the loop is to
// advance the solutions, i.e. to
- // copy the temperature solution to
- // the next "older" time level.
+ // copy the solutions to the next
+ // "older" time level.
assemble_stokes_system ();
build_stokes_preconditioner ();
assemble_temperature_matrix ();
time += time_step;
++timestep_number;
+ old_stokes_solution = stokes_solution;
old_old_temperature_solution = old_temperature_solution;
- old_temperature_solution = temperature_solution;
+ old_temperature_solution = temperature_solution;
}
- // Do all the above until we arrive
- // at time 100.
+ // Do all the above until we arrive at
+ // time 100.
while (time <= 100);
}
// @sect3{The <code>main</code> function}
//
- // The main function looks almost
- // the same as in all other
- // programs. The only difference is
- // that Trilinos wants to get the
- // arguments from calling the
- // function (argc and argv) in
- // order to correctly set up the
- // MPI system in case we use those
- // compilers (even though this
- // program is only meant to be run
- // in serial).
+ // The main function looks almost the same
+ // as in all other programs. The only
+ // difference is that Trilinos wants to get
+ // the arguments from calling the function
+ // (argc and argv) in order to correctly
+ // set up the MPI system in case those
+ // compilers are in use (even though this
+ // program is only meant to be run on one
+ // processor).
int main (int argc, char *argv[])
{
try