<h3> Navier Stokes Equations </h3>
In this tutorial we show how to solve the incompressible Navier
-Stokes equations (NSE) by Newton's method. The flow we consider here
+Stokes equations (NSE) with Newton's method. The flow we consider here
is assumed to be steady. In a domain $\Omega \subset
\mathbb{R}^{d}$, $d=2,3$, with a piecewise smooth boundary
$\partial \Omega$, and a given force field $\textbf{f}$, we seek
@f{eqnarray*}
- \nu \Delta\textbf{u} + (\textbf{u} \cdot \nabla)\textbf{u} + \nabla p &=& \textbf{f}\\
-- \nabla \cdot \textbf{u} &=& 0
+- \nabla \cdot \textbf{u} &=& 0.
@f}
-Different from the Stokes equations as discussed in step-22, the NSE are a
-nonlinear system because of the convective term $(\textbf{u} \cdot
+Unlike the Stokes equations as discussed in step-22, the NSE are a
+nonlinear system of equations because of the convective term $(\textbf{u} \cdot
\nabla)\textbf{u}$. The first step of computing a numerical solution
is to linearize the system and this will be done using Newton's method. A
time-dependent problem is discussed in step-35, where the system is linearized
<h3> Linearization of Navier-Stokes Equations </h3>
-Moving the right-hand side terms to the left, a nonlinear function is created as
+We define a nonlinear function whose root is a solution to the NSE by
@f{eqnarray*}
F(\mathbf{u}, p) =
\end{pmatrix}.
@f}
-$F(\textbf{u}, p)$ is a nonlinear function whose root is
-the solution to the NSE. Assuming the initial guess is good enough to
+Assuming the initial guess is good enough to
guarantee the convergence of Newton's iteration and denoting
-$\textbf{x} = (\textbf{u}, p)$, Newton's iteration on a vector field
+$\textbf{x} = (\textbf{u}, p)$, Newton's iteration on a vector function
can be defined as
@f{eqnarray*}
\textbf{x}^{k+1} = \textbf{x}^{k} - (\nabla F(\textbf{x}^{k}))^{-1} F(\textbf{x}^{k}),
@f}
where
-$\textbf{x}^{k+1}$ is the approximate solution in step k+1,
-$\textbf{x}^{k}$ represents the solution from the last step, and $\nabla
+$\textbf{x}^{k+1}$ is the approximate solution in step $k+1$,
+$\textbf{x}^{k}$ represents the solution from the previous step, and $\nabla
F(\textbf{x}^{k})$ is the Jacobian matrix evaluated at
$\textbf{x}^{k}$.
A similar iteration can be found in step-15.
-From Newton's iteration formula, we can observe that the new
+The Newton iteration formula implies the new
solution is obtained by adding an update term to the old solution. Instead
of evaluating the Jacobian matrix and taking its inverse, we consider
the update term as a whole, that is
@f{eqnarray*}
-\delta \textbf{x}^{k} = - (\nabla F(\textbf{x}^{k}))^{-1} F(\textbf{x}^{k}),
+ \delta \textbf{x}^{k} = - (\nabla F(\textbf{x}^{k}))^{-1} F(\textbf{x}^{k}),
@f}
where $\textbf{x}^{k+1}=\textbf{x}^{k}+\delta \textbf{x}^{k}$.
-Then we can evaluate the update term by solving the system
+We can find the update term by solving the system
@f{eqnarray*}
\nabla F(\textbf{x}^{k}) \delta \textbf{x}^{k} = -F(\textbf{x}^{k}).
@f}
+ \mathbf{u}^{k} \cdot \nabla \delta \mathbf{u}^{k}
+ \delta \mathbf{u}^{k} \cdot \nabla \mathbf{u}^{k}
+ \nabla \delta p^{k}
- = \mathbf{g}, \\
+ = -F(\mathbf{x}^k), \\
-\nabla \cdot\delta \mathbf{u}^{k}
= \nabla \cdot \mathbf{u}^{k},
@f}
-where $\textbf{g} =\textbf{f} + \nu \Delta \textbf{u}^k - (\textbf{u}^k
-\cdot \nabla)\textbf{u}^k -\nabla p^k$ and $\textbf{u}^k$ and $p^k$ are the solutions from the
+where $\textbf{u}^k$ and $p^k$ are the solutions from the
previous iteration. Additionally, the
right hand side of the second equation is not zero since the discrete
solution is not exactly divergence free (divergence free for the continuous
solution). The right hand side here acts as a correction which leads the
discrete solution of the velocity to be divergence free along Newton's
iteration. In this linear system, the only unknowns are the
-update terms $\delta \textbf{u}^{k}$ and $\delta p^{k}$, and we can use a similar strategy
-to the one used in step-22. The weak form is
-derived like it is done in step-22.
+update terms $\delta \textbf{u}^{k}$ and $\delta p^{k}$, and we can use a
+similar strategy to the one used in step-22 (and derive the weak form in the
+same way).
Now, Newton's iteration can be used to solve for the update terms:
<h3> Finding an Initial Guess </h3>
-Getting Newton's method to converge, the initial guess needs to be close
-enough to the solution, so it is crucial to find a good starting value.
+The initial guess needs to be close enough to the solution for Newton's method
+to converge; hence, finding a good starting value is crucial to the nonlinear
+solver.
When the viscosity $\nu$ is large, a good initial guess can be obtained
by solving the Stokes equation with viscosity $\nu$. While problem dependent,
<h3> Test Case </h3>
-Here we use the lid driven cavity flow as our test case, see [3] for details.
-The computational domain is the unit square and the right-hand side
+We use the lid driven cavity flow as our test case; see [3] for details.
+The computational domain is the unit square and the right-hand side is
$f=0$. The boundary condition is
@f{eqnarray*}
Inside the loop, we involve three solvers: one for $\tilde{A}^{-1}$,
one for $M_p^{-1}$ and one for $Gx=b$. The first two
solvers are invoked in the preconditioner and the outer solver gives us
-the update term. Overall convergence is controlled by the nonlinear residual
-and Newton's method does not have to require an exact Jacobian, so for the outer
-linear solver we employ FGMRES with a relative tolerance of only 1e-4. In fact,
+the update term. Overall convergence is controlled by the nonlinear residual;
+as Newton's method does not require an exact Jacobian, we employ FGMRES with a
+relative tolerance of only 1e-4 for the outer linear solver. In fact,
we use the truncated Newton solve for this system.
As described in step-22, the inner linear solves are also not required
to be done very accurately. Here we use CG with a relative
// @sect3{The <code>NavierStokesProblem</code> class template}
- // This is the main function and its member functions.
- // As explained in the introduction, what we obtain at each step is the
- // Newton's update, so we define two variables: the present solution
- // and the update. Additionally, the evaluation point is
- // for temporarily holding Newton update in line search. A sparse matrix
- // for the pressure mass matrix is created for the operator of a block Schur
- // complement preconditioner. We use one AffineConstraints object for
- // Dirichlet boundary conditions at the initial step and a zero
- // AffineConstraints object for the Newton is defined by 1/Re which has been
- // discussed in the introduction.
-
+ // This class manages the matrices and vectors described in the
+ // introduction: in particular, we store a BlockVector for the current
+ // solution, current Newton update, and the line search update. We also
+ // store two AffineConstraints objects: one which enforces the Dirichlet
+ // boundary conditions and one that sets all boundary values to zero. The
+ // first constrains the solution vector while the second constraints the
+ // updates (i.e., we never update boundary values, so we force the relevant
+ // update vector values to be zero).
template <int dim>
class StationaryNavierStokes
{
BlockVector<double> evaluation_point;
};
- // @sect3{Boundary values and right hand side}
- // In this problem we set the velocity along the upper surface of the cavity
- // to be one and zero on the other three walls. The right hand side function
- // is zero so we do not need to set the right hand side function in this
- // tutorial. The number of components of the boundary function is dim+1.
- // In practice, the boundary values are
- // applied to our solution through an AffineConstraints object which is
- // obtained by using VectorTools::interpolate_boundary_values. The components
- // of boundary value functions are required to be chosen according to the
- // finite element space. Therefore we have to define the boundary value of
- // pressure even though we actually do not need it.
-
- // The following function represents the boundary values:
+ // @sect3{Boundary values and right hand side} In this problem we set the
+ // velocity along the upper surface of the cavity to be one and zero on the
+ // other three walls. The right hand side function is zero so we do not need
+ // to set the right hand side function in this tutorial. The number of
+ // components of the boundary function is <code>dim+1</code>. We will
+ // ultimately use VectorTools::interpolate_boundary_values to set boundary
+ // values, which requires the boundary value functions to have the same
+ // number of components as the solution, even if all are not used. Put
+ // another way: to make this function happy we define boundary values for
+ // the pressure even though we will never actually use them.
template <int dim>
class BoundaryValues : public Function<dim>
{
// @sect3{BlockSchurPreconditioner for Navier Stokes equations}
//
- // The block
- // Schur complement preconditioner is defined in this part. As discussed in
- // the introduction, the preconditioner in Krylov iterative methods is
- // implemented as a matrix-vector product operator. In practice, the Schur
- // complement preconditioner is decomposed as a product of three matrices (as
- // presented in the first section). The $\tilde{A}^{-1}$ in the first factor
- // involves a solve for the linear system $\tilde{A}x=b$. Here we solve
- // this system via a direct solver for simplicity. The computation involved
- // in the second factor is a simple matrix-vector multiplication. The Schur
- // complement $\tilde{S}$ can be well approximated by the pressure mass
- // matrix and its inverse can be obtained through an inexact solver. Because
- // the pressure mass matrix is symmetric and positive definite, we can use
- // CG to solve the corresponding linear system.
- //
+ // As discussed in the introduction, the preconditioner in Krylov iterative
+ // methods is implemented as a matrix-vector product operator. In practice,
+ // the Schur complement preconditioner is decomposed as a product of three
+ // matrices (as presented in the first section). The $\tilde{A}^{-1}$ in the
+ // first factor involves a solve for the linear system $\tilde{A}x=b$. Here
+ // we solve this system via a direct solver for simplicity. The computation
+ // involved in the second factor is a simple matrix-vector
+ // multiplication. The Schur complement $\tilde{S}$ can be well approximated
+ // by the pressure mass matrix and its inverse can be obtained through an
+ // inexact solver. Because the pressure mass matrix is symmetric and
+ // positive definite, we can use CG to solve the corresponding linear
+ // system.
template <class PreconditionerMp>
class BlockSchurPreconditioner : public Subscriptor
{
SparseDirectUMFPACK A_inverse;
};
- // We can notice that the initialization of the inverse of the matrix at (0,0)
- // corner is completed in the constructor. If so, every application of the
- // preconditioner then no longer requires the computation of the matrix
- // factors.
+ // We can notice that the initialization of the inverse of the matrix at the
+ // top left corner is completed in the constructor. If so, every application
+ // of the preconditioner then no longer requires the computation of the
+ // matrix factors.
template <class PreconditionerMp>
BlockSchurPreconditioner<PreconditionerMp>::BlockSchurPreconditioner(
// @sect4{StationaryNavierStokes::StationaryNavierStokes}
// The constructor of this class looks very similar to the one in step-22. The
// only difference is the viscosity and the Augmented Lagrangian coefficient
- // gamma.
- //
-
+ // <code>gamma</code>.
template <int dim>
StationaryNavierStokes<dim>::StationaryNavierStokes(const unsigned int degree)
: viscosity(1.0 / 7500.0)
//
// This function implements the Newton iteration with given tolerance, maximum
// number of iterations, and the number of mesh refinements to do.
- // "is_initial_step" is the flag to tell us whether "setup_system" is
- // necessary, and which part, system matrix or right hand side vector, should
- // be assembled. If we do a line search, the right hand side is already
- // assembled while checking the residual norm in the last iteration.
- // Therefore, we just need to assemble the system matrix at the current
- // iteration. The last argument "output_result" is whether output should be
- // produced.
-
+ //
+ // The argument <code>is_initial_step</code> tells us whether
+ // <code>setup_system</code> is necessary, and which part, system matrix or
+ // right hand side vector, should be assembled. If we do a line search, the
+ // right hand side is already assembled while checking the residual norm in
+ // the last iteration. Therefore, we just need to assemble the system matrix
+ // at the current iteration. The last argument <code>output_result</code>
+ // determines whether or not graphical output should be produced.
template <int dim>
void StationaryNavierStokes<dim>::newton_iteration(
const double tolerance,
solve(first_step);
// To make sure our solution is getting close to the exact
- // solution, we let the solution be updated with a weight alpha
- // such that the new residual is smaller than the one of last
- // step, which is done in the following loop. Also the line
- // search method can be located in step-15.
-
+ // solution, we let the solution be updated with a weight
+ // <code>alpha</code> such that the new residual is smaller
+ // than the one of last step, which is done in the following
+ // loop. This is the same line search algorithm used in
+ // step-15.
for (double alpha = 1.0; alpha > 1e-5; alpha *= 0.5)
{
evaluation_point = present_solution;
projects / dealii.git / commitdiff