 Reference documentation for deal.II version Git b1a5775265 2021-07-23 17:28:58 -0400
The step-57 tutorial program

This tutorial depends on step-15, step-22.

1. Introduction
2. The commented program
1. Results
2. The plain program

This program was contributed by Liang Zhao and Timo Heister.

This material is based upon work partially supported by National Science Foundation grant DMS1522191 and the Computational Infrastructure in Geodynamics initiative (CIG), through the National Science Foundation under Award No. EAR-0949446 and The University of California-Davis.

Note
If you use this program as a basis for your own work, please consider citing it in your list of references. The initial version of this work was contributed to the deal.II project by the authors listed in the following citation: # Introduction

### Navier Stokes Equations

In this tutorial we show how to solve the incompressible Navier 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 a velocity field $$\textbf{u}$$ and a pressure field $$\textbf{p}$$ satisfying

\begin{eqnarray*} - \nu \Delta\textbf{u} + (\textbf{u} \cdot \nabla)\textbf{u} + \nabla p &=& \textbf{f}\\ - \nabla \cdot \textbf{u} &=& 0. \end{eqnarray*}

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 using the solution from the last time step and no nonlinear solve is necessary.

### Linearization of Navier-Stokes Equations

We define a nonlinear function whose root is a solution to the NSE by

\begin{eqnarray*} F(\mathbf{u}, p) = \begin{pmatrix} - \nu \Delta\mathbf{u} + (\mathbf{u} \cdot \nabla)\mathbf{u} + \nabla p - \mathbf{f} \\ - \nabla \cdot \mathbf{u} \end{pmatrix}. \end{eqnarray*}

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 function can be defined as

\begin{eqnarray*} \textbf{x}^{k+1} = \textbf{x}^{k} - (\nabla F(\textbf{x}^{k}))^{-1} F(\textbf{x}^{k}), \end{eqnarray*}

where $$\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.

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

\begin{eqnarray*} \delta \textbf{x}^{k} = - (\nabla F(\textbf{x}^{k}))^{-1} F(\textbf{x}^{k}), \end{eqnarray*}

where $$\textbf{x}^{k+1}=\textbf{x}^{k}+\delta \textbf{x}^{k}$$.

We can find the update term by solving the system

\begin{eqnarray*} \nabla F(\textbf{x}^{k}) \delta \textbf{x}^{k} = -F(\textbf{x}^{k}). \end{eqnarray*}

Here, the left of the previous equation represents the directional gradient of $$F(\textbf{x})$$ along $$\delta \textbf{x}^{k}$$ at $$\textbf{x}^{k}$$. By definition, the directional gradient is given by

\begin{eqnarray*} & &\nabla F(\mathbf{u}^{k}, p^{k}) (\delta \mathbf{u}^{k}, \delta p^{k}) \\ \\ &=& \lim_{\epsilon \to 0} \frac{1}{\epsilon} \left( F(\mathbf{u}^{k} + \epsilon \delta \mathbf{u}^{k}, p^{k} + \epsilon \nabla \delta p^{k}) - F(\mathbf{u}^{k}, p^{k}) \right)\\ \\ &=& \lim_{\epsilon \to 0} \frac{1}{\epsilon} \begin{pmatrix} - \epsilon \nu \Delta \delta \mathbf{u}^{k} + \epsilon \mathbf{u}^{k} \cdot \nabla \delta \mathbf{u}^{k} + \epsilon \delta \mathbf{u}^{k} \cdot \nabla \mathbf{u}^{k} + \epsilon^{2} \delta \mathbf{u}^{k} \cdot \nabla \delta \mathbf{u}^{k} + \epsilon \nabla \delta p^{k}\\ - \epsilon \nabla \cdot \delta \mathbf{u}^{k} \end{pmatrix} \\ \\ &=& \begin{pmatrix} - \nu \Delta \delta \mathbf{u}^{k} + \mathbf{u}^{k} \cdot \nabla \delta \mathbf{u}^{k} + \delta \mathbf{u}^{k} \cdot \nabla \mathbf{u}^{k} + \nabla \delta p^{k}\\ - \nabla \cdot \delta \mathbf{u}^{k} \end{pmatrix}. \end{eqnarray*}

Therefore, we arrive at the linearized system:

\begin{eqnarray*} -\nu \Delta \delta \mathbf{u}^{k} + \mathbf{u}^{k} \cdot \nabla \delta \mathbf{u}^{k} + \delta \mathbf{u}^{k} \cdot \nabla \mathbf{u}^{k} + \nabla \delta p^{k} = -F(\mathbf{x}^k), \\ -\nabla \cdot\delta \mathbf{u}^{k} = \nabla \cdot \mathbf{u}^{k}, \end{eqnarray*}

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 (and derive the weak form in the same way).

Now, Newton's iteration can be used to solve for the update terms:

1. Initialization: Initial guess $$u_0$$ and $$p_0$$, tolerance $$\tau$$;
2. Linear solve to compute update term $$\delta\textbf{u}^{k}$$ and $$\delta p^k$$;
3. Update the approximation: $$\textbf{u}^{k+1} = \textbf{u}^{k} + \delta\textbf{u}^{k}$$ and $$p^{k+1} = p^{k} + \delta p^{k}$$;
4. Check residual norm: $$E^{k+1} = \|F(\mathbf{u}^{k+1}, p^{k+1})\|$$:
• If $$E^{k+1} \leq \tau$$, STOP.
• If $$E^{k+1} > \tau$$, back to step 2.

### Finding an Initial Guess

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, this works for $$\nu \geq 1/400$$ for the test problem considered here.

However, the convective term $$(\mathbf{u}\cdot\nabla)\mathbf{u}$$ will be dominant if the viscosity is small, like $$1/7500$$ in test case 2. In this situation, we use a continuation method to set up a series of auxiliary NSEs with viscosity approaching the one in the target NSE. Correspondingly, we create a sequence $$\{\nu_{i}\}$$ with $$\nu_{n}= \nu$$, and accept that the solutions to two NSE with viscosity $$\nu_{i}$$ and $$\nu_{i+1}$$ are close if $$|\nu_{i} - \nu_{i+1}|$$ is small. Then we use the solution to the NSE with viscosity $$\nu_{i}$$ as the initial guess of the NSE with $$\nu_{i+1}$$. This can be thought of as a staircase from the Stokes equations to the NSE we want to solve.

That is, we first solve a Stokes problem

\begin{eqnarray*} -\nu_{1} \Delta \textbf{u} + \nabla p &=& \textbf{f}\\ -\nabla \cdot \textbf{u} &=& 0 \end{eqnarray*}

to get the initial guess for

\begin{eqnarray*} -\nu_{1} \Delta \textbf{u} + (\textbf{u} \cdot \nabla)\textbf{u} + \nabla p &=& \textbf{f},\\ -\nabla \cdot \textbf{u} &=& 0, \end{eqnarray*}

which also acts as the initial guess of the continuation method. Here $$\nu_{1}$$ is relatively large so that the solution to the Stokes problem with viscosity $$\nu_{1}$$ can be used as an initial guess for the NSE in Newton's iteration.

Then the solution to

\begin{eqnarray*} -\nu_{i} \Delta \textbf{u} + (\textbf{u} \cdot \nabla)\textbf{u} + \nabla p &=& \textbf{f},\\ -\nabla \cdot \textbf{u} &=& 0. \end{eqnarray*}

acts as the initial guess for

\begin{eqnarray*} -\nu_{i+1} \Delta \textbf{u} + (\textbf{u} \cdot \nabla)\textbf{u} + \nabla p &=& \textbf{f},\\ -\nabla \cdot \textbf{u} &=& 0. \end{eqnarray*}

This process is repeated with a sequence of viscosities $$\{\nu_i\}$$ that is determined experimentally so that the final solution can used as a starting guess for the Newton iteration.

### The Solver and Preconditioner

At each step of Newton's iteration, the problem results in solving a saddle point systems of the form

\begin{eqnarray*} \begin{pmatrix} A & B^{T} \\ B & 0 \end{pmatrix} \begin{pmatrix} U \\ P \end{pmatrix} = \begin{pmatrix} F \\ 0 \end{pmatrix}. \end{eqnarray*}

This system matrix has the same block structure as the one in step-22. However, the matrix $$A$$ at the top left corner is not symmetric because of the nonlinear term. Instead of solving the above system, we can solve the equivalent system

\begin{eqnarray*} \begin{pmatrix} A + \gamma B^TW^{-1}B & B^{T} \\ B & 0 \end{pmatrix} \begin{pmatrix} U \\ P \end{pmatrix} = \begin{pmatrix} F \\ 0 \end{pmatrix} \end{eqnarray*}

with a parameter $$\gamma$$ and an invertible matrix $$W$$. Here $$\gamma B^TW^{-1}B$$ is the Augmented Lagrangian term; see  for details.

Denoting the system matrix of the new system by $$G$$ and the right-hand side by $$b$$, we solve it iteratively with right preconditioning $$P^{-1}$$ as $$GP^{-1}y = b$$, where

\begin{eqnarray*} P^{-1} = \begin{pmatrix} \tilde{A} & B^T \\ 0 & \tilde{S} \end{pmatrix}^{-1} \end{eqnarray*}

with $$\tilde{A} = A + \gamma B^TW^{-1}B$$ and $$\tilde{S}$$ is the corresponding Schur complement $$\tilde{S} = B^T \tilde{A}^{-1} B$$. We let $$W = M_p$$ where $$M_p$$ is the pressure mass matrix, then $$\tilde{S}^{-1}$$ can be approximated by

\begin{eqnarray*} \tilde{S}^{-1} \approx -(\nu+\gamma)M_p^{-1}. \end{eqnarray*}

See  for details.

We decompose $$P^{-1}$$ as

\begin{eqnarray*} P^{-1} = \begin{pmatrix} \tilde{A}^{-1} & 0 \\ 0 & I \end{pmatrix} \begin{pmatrix} I & -B^T \\ 0 & I \end{pmatrix} \begin{pmatrix} I & 0 \\ 0 & \tilde{S}^{-1} \end{pmatrix}. \end{eqnarray*}

Here two inexact solvers will be needed for $$\tilde{A}^{-1}$$ and $$\tilde{S}^{-1}$$, respectively (see ). Since the pressure mass matrix is symmetric and positive definite, CG with ILU as a preconditioner is appropriate to use for $$\tilde{S}^{-1}$$. For simplicity, we use the direct solver UMFPACK for $$\tilde{A}^{-1}$$. The last ingredient is a sparse matrix-vector product with $$B^T$$. Instead of computing the matrix product in the augmented Lagrangian term in $$\tilde{A}$$, we assemble Grad-Div stabilization $$(\nabla \cdot \phi _{i}, \nabla \cdot \phi _{j}) \approx (B^T M_p^{-1}B)_{ij}$$, as explained in .

### Test Case

We use the lid driven cavity flow as our test case; see  for details. The computational domain is the unit square and the right-hand side is $$f=0$$. The boundary condition is

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; 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 tolerance of 1e-6 for the pressure mass matrix. As expected, we still see convergence of the nonlinear residual down to 1e-14. Also, we use a simple line search algorithm for globalization of the Newton method.
The cavity reference values for $$\mathrm{Re}=400$$ and $$\mathrm{Re}=7500$$ are from  and , respectively, where $$\mathrm{Re}$$ is the Reynolds number and can be located at . Here the viscosity is defined by $$1/\mathrm{Re}$$. Even though we can still find a solution for $$\mathrm{Re}=10000$$ and the references contain results for comparison, we limit our discussion here to $$\mathrm{Re}=7500$$. This is because the solution is no longer stationary starting around $$\mathrm{Re}=8000$$ but instead becomes periodic, see  for details.