1. Remove trailing whitespace.
2. Fix some minor spelling mistakes.
3. Consistently use MathJax for equations.
<br>
<i>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
\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
+using the solution from the last time step and no nonlinear
solve is necessary.
<h3> Linearization of Navier-Stokes Equations </h3>
\right).
@f}
-$F(\textbf{u}, p)$ is a nonlinear function whose root is
+$F(\textbf{u}, p)$ is a nonlinear function whose root is
the solution to the NSE. 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
can be defined as
-@f{eqnarray*} \textbf{x}^{k+1} = \textbf{x}^{k} -
-(\nabla F(\textbf{x}^{k}))^{-1} F(\textbf{x}^{k}),
+@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,
@f{eqnarray*}
\delta \textbf{x}^{k} = - (\nabla F(\textbf{x}^{k}))^{-1} F(\textbf{x}^{k}),
@f}
-where $x^{k+1}=x^{k}+\delta x^{k}$.
+where $\textbf{x}^{k+1}=\textbf{x}^{k}+\delta \textbf{x}^{k}$.
Then we can evaluate the update term by solving the system
@f{eqnarray*}
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
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
+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
Now, Newton's iteration can be used to solve for the update terms:
<ol>
- <li> Initialization: Initial guess $u_0$ and $p_0$, tolerance $\tau$;
- <li> Linear solve to compute update term $\delta\textbf{u}^{k}$ and $\delta p^k$;
- <li> Update the approximation: $\textbf{u}^{k+1} = \textbf{u}^{k} + \delta\textbf{u}^{k}$ and $p^{k+1} = p^{k} + \delta p^{k}$;
- <li> Check residual norm: $E^{k+1} = \|F(\mathbf{u}^{k+1}, p^{k+1})\|$:
- <ul>
- <li>If $E^{k+1} \leq \tau$, STOP.
- <li>If $E^{k+1} > \tau$, back to step 2.
- </ul>
+ <li>Initialization: Initial guess $u_0$ and $p_0$, tolerance $\tau$;</li>
+ <li>Linear solve to compute update term $\delta\textbf{u}^{k}$ and
+ $\delta p^k$;</li>
+ <li>Update the approximation:
+ $\textbf{u}^{k+1} = \textbf{u}^{k} + \delta\textbf{u}^{k}$ and
+ $p^{k+1} = p^{k} + \delta p^{k}$;</li>
+ <li>Check residual norm: $E^{k+1} = \|F(\mathbf{u}^{k+1}, p^{k+1})\|$:
+ <ul>
+ <li>If $E^{k+1} \leq \tau$, STOP.</li>
+ <li>If $E^{k+1} > \tau$, back to step 2.</li>
+ </ul></li>
</ol>
<h3> Finding an Initial Guess </h3>
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 NSE with
+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} -
@f}
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.
+can be used as an initial guess for the NSE in Newton's iteration.
Then the solution to
- \nabla \cdot \textbf{u} &=& 0.
@f}
-This process is repeated with a sequence of viscosities, $\{\nu_i\}$ that is
+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.
-<h3>The Solver and Preconditioner </h3>
+<h3>The %Solver and Preconditioner </h3>
At each step of Newton's iteration, the problem results in solving a
saddle point systems of the form
@f}
This system matrix has the same block structure as the one in step-22. However,
-the matrix $A$ at (1, 1) corner is not symmetric because of the nonlinear term.
+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
@f{eqnarray*}
\end{array}
\right)
@f}
-with a parameter $\gamma$ and an invertible matrix W. Here
-$\gamma B^TW^{-1}B$ is the Augmented Lagrangian term and
- see [1] for details.
+with a parameter $\gamma$ and an invertible matrix $W$. Here
+$\gamma B^TW^{-1}B$ is the Augmented Lagrangian term; see [1] 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
@f{eqnarray*}
-P^{-1} = \left(\begin{array}{cc} \tilde{A} & B^T \\
+P^{-1} = \left(\begin{array}{cc} \tilde{A} & B^T \\
0 & \tilde{S} \end{array}\right)^{-1},
@f}
We decompose $P^{-1}$ as
-@f{eqnarray*}
+@f{eqnarray*}
P^{-1} =
\left(\begin{array}{cc} \tilde{A}^{-1} & 0 \\ 0 & I \end{array}\right)
\left(\begin{array}{cc} I & -B^T \\ 0 & I \end{array}\right)
Here two inexact solvers will be needed for $\tilde{A}^{-1}$ and
$\tilde{S}^{-1}$, respectively (see [1]). 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
+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
<h3> Test Case </h3>
-Here we use the lid driven cavity flow as our test case, see [3] for details.
+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
$f=0$. The boundary condition is
-@f{eqnarray*}
-(u(x, y), v(x,y)) &=& (1,0)
+@f{eqnarray*}
+(u(x, y), v(x,y)) &=& (1,0)
\qquad\qquad \textrm{if}\ y = 1 \\
(u(x, y), v(x,y)) &=& (0,0)
- \qquad\qquad \textrm{else}.
+ \qquad\qquad \textrm{otherwise}.
@f}
When solving this problem, the error consists of the nonlinear error (from
-Newton's iteration) and the discretization error (depending on mesh size). The
+Newton's iteration) and the discretization error (dependent on mesh size). The
nonlinear part decreases with each Newton iteration and the discretization error
reduces with mesh refinement. In this example, the solution from the coarse
mesh is transferred to successively finer meshes and used as an initial
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,
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
+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 Re=400 and Re=7500 are from [4] and [5],
-respectively, where "Re" represents the Reynold number and can be located
-at [8]. Here the viscosity is defined by 1/Re.
-Even though we can still find a solution for Re=10000 and the
+The cavity reference values for $\mathrm{Re}=400$ and $\mathrm{Re}=7500$ are
+from [4] and [5], respectively, where $\mathrm{Re}$ is the Reynolds number and
+can be located at [8]. 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
-Re=7500. This is because the solution is no longer stationary starting around
-Re=8000 but instead becomes periodic, see [7] for details.
+$\mathrm{Re}=7500$. This is because the solution is no longer stationary
+starting around $\mathrm{Re}=8000$ but instead becomes periodic, see [7] for
+details.
-<h3> Reference </h3>
+<h3> References </h3>
<ol>
<li> An Augmented Lagrangian-Based Approach to the Oseen Problem, M. Benzi and M. Olshanskii, SIAM J. SCI. COMPUT. 2006
<h1>Results</h1>
Now we use the method we discussed above to solve Navier Stokes equations with
-viscosity 1/400 and 1/7500.
+viscosity $1/400$ and $1/7500$.
-<h3> Test case 1: Re=400 </h3>
+<h3> Test case 1: Low Reynolds Number </h3>
-In the first test case the viscosity is set to be 1/400. As we discussed in the
+In the first test case the viscosity is set to be $1/400$. As we discussed in the
introduction, the initial guess is the solution to the corresponding Stokes
problem. In the following table, the residuals at each Newton's iteration on
-every mesh is shown. The data in the table shows that Newton's iteration converges quadratically.
+every mesh is shown. The data in the table shows that Newton's iteration
+converges quadratically.
<table align="center" class="doxtable">
<tr>
- <th>Re=400</th>
+ <th>$\mathrm{Re}=400$</th>
<th colspan="2">Mesh0</th>
<th colspan="2">Mesh1</th>
<th colspan="2">Mesh2</th>
<tr>
<td>2</td>
<td>7.0849e-04</td>
- <td>5.0000e+00</td>
+ <td>5</td>
<td>9.9458e-04</td>
<td>5</td>
<td>1.1409e-04</td>
<tr>
<td>3</td>
<td>1.9980e-05</td>
- <td>5.0000e+00</td>
+ <td>5</td>
<td>4.5007e-05</td>
<td>5</td>
<td>2.9020e-08</td>
<tr>
<td>4</td>
<td>2.3165e-09</td>
- <td>6.0000e+00</td>
+ <td>6</td>
<td>1.6891e-07</td>
<td>5</td>
<td>1.2338e-14</td>
<tr>
<td>5</td>
<td>1.2585e-13</td>
- <td>7.0000e+00</td>
+ <td>7</td>
<td>1.4520e-11</td>
<td>6</td>
<td>1.9044e-13</td>
The following figures show the sequence of generated grids. For the case
-of Re=400, the initial guess is obtained by solving Stokes on an $8 \times 8$
-mesh, and the mesh is refined adaptively. Between meshes, the solution from
-the coarse mesh is interpolated to the fine mesh to be used as an initial guess.
+of $\mathrm{Re}=400$, the initial guess is obtained by solving Stokes on an
+$8 \times 8$ mesh, and the mesh is refined adaptively. Between meshes, the
+solution from the coarse mesh is interpolated to the fine mesh to be used as an
+initial guess.
<TABLE ALIGN="center">
<tr>
</tr>
</table>
-This picture is the graphical streamline result of lid-driven cavity with Re=400.
+This picture is the graphical streamline result of lid-driven cavity with
+$\mathrm{Re}=400$.
<img src="https://www.dealii.org/images/steps/developer/step-57.Re400_Streamline.png" alt="">
Then the solution is compared with a reference solution
<img src="https://www.dealii.org/images/steps/developer/step-57.compare-Re400.svg" style="width:50%" alt="">
-<h3> Test case 2: Re=7500 </h3>
+<h3> Test case 2: High Reynolds Number </h3>
Newton's iteration requires a good initial guess. However, the nonlinear term
-dominates when the Reynold number is large, so that the solution to the Stokes
+dominates when the Reynolds number is large, so that the solution to the Stokes
equations may be far away from the exact solution. If the Stokes solution acts
as the initial guess, the convergence will be lost. The following picture
-shows that the nonlinear iteration gets stuck and the residual no longer decreases
+shows that the nonlinear iteration gets stuck and the residual no longer decreases
in further iterations.
<img src="https://www.dealii.org/images/steps/developer/step-57.Re7500_loss_convergence.svg" style="width:50%" alt="">
-The initial guess, therefore, has to be obtained via a continuation method
+The initial guess, therefore, has to be obtained via a continuation method
which has been discussed in the introduction. Here the step size in the continuation method, that is $|\nu_{i}-\nu_{i+1}|$, is 2000 and the initial
mesh is of size $32 \times 32$. After obtaining an initial guess, the mesh is
refined as in the previous test case. The following picture shows that at each
-refinement Newton's iteration has quadratic convergence. 52 steps of Newton's
+refinement Newton's iteration has quadratic convergence. 52 steps of Newton's
iterations are executed for solving this test case.
<img src="https://www.dealii.org/images/steps/developer/step-57.Re7500_get_convergence.svg" style="width:50%" alt="">
<table align="center" class="doxtable">
<tr>
- <th>Re=7500</th>
+ <th>$\mathrm{Re}=7500$</th>
<th colspan="2">Mesh0</th>
<th colspan="2">Mesh1</th>
<th colspan="2">Mesh2</th>
<img src="https://www.dealii.org/images/steps/developer/step-57.Re7500_Streamline.png" alt="">
Furthermore, the error consists of the nonlinear error,
-which is decreasing as Newton's iteration going on, and the discretization error,
+which is decreasing as Newton's iteration going on, and the discretization error,
which depends on the mesh size. That is why we have to refine the
mesh and repeat Newton's iteration on the next finer mesh. From the table above, we can
see that the residual (nonlinear error) is below 1e-12 on each mesh, but the
-following picture shows us the difference between solutions on subsequently finer
+following picture shows us the difference between solutions on subsequently finer
meshes.
<img src="https://www.dealii.org/images/steps/developer/step-57.converge-Re7500.svg" style="width:50%" alt="">
paper [2].
The following table shows the timing results between our iterative approach
- (FGMRES) compared to a direct solver (UMFPACK) for the whole system
+(FGMRES) compared to a direct solver (UMFPACK) for the whole system
with viscosity set to 1/400. Even though we use the same direct solver for
the velocity block in the iterative solver, it is considerably faster and
-consumes less memory. This will be even more pronounced in 3d.
+consumes less memory. This will be even more pronounced in 3d.
-<table align="center" border="1">
+<table align="center" class="doxtable">
<tr>
- <th>Ref</th>
+ <th>Refinement Cycle</th>
<th>DoFs</th>
<th>Iterative: Total/s (Setup/s)</th>
<th>Direct: Total/s (Setup/s)</th>
</table>
-<h4>3D computations</h4>
+<h4>3d computations</h4>
The code is set up to also run in 3d. Of course the reference values are
different, see [6] for example. High resolution computations are not doable
with this example as is, because a direct solver for the velocity block does
not work well in 3d. Rather, a parallel solver based on algebraic or geometric
multigrid is needed. See below.
-
+
<h4>Parallelization</h4>
-For larger computations, especially in 3D, it is necessary to implement MPI
+For larger computations, especially in 3d, it is necessary to implement MPI
parallel solvers and preconditioners. A good starting point would be step-55,
which uses algebraic multigrid for the velocity block for the Stokes
equations.
private:
void setup_dofs();
+
void initialize_system();
+
void assemble(const bool initial_step, const bool assemble_matrix);
+
void assemble_system(const bool initial_step);
+
void assemble_rhs(const bool initial_step);
+
void solve(const bool initial_step);
+
void refine_mesh();
+
void process_solution(unsigned int refinement);
+
void output_results(const unsigned int refinement_cycle) const;
+
void newton_iteration(const double tolerance,
const unsigned int max_iteration,
const unsigned int n_refinements,
const bool is_initial_step,
const bool output_result);
+
void compute_initial_guess(double step_size);
double viscosity;
// @sect4{StationaryNavierStokes::setup_dofs}
// This function initializes the DoFHandler enumerating the degrees of freedom
// and constraints on the current mesh.
-
template <int dim>
void StationaryNavierStokes<dim>::setup_dofs()
{
// We renumber the components to have all velocity DoFs come before
// the pressure DoFs to be able to split the solution vector in two blocks
// which are separately accessed in the block preconditioner.
- //
std::vector<unsigned int> block_component(dim + 1, 0);
block_component[dim] = 1;
DoFRenumbering::component_wise(dof_handler, block_component);
// remain satisfied during Newton's iteration, zero boundary conditions are
// used for the update $\delta u^k$. Therefore we set up two different
// constraint objects.
-
FEValuesExtractors::Vector velocities(0);
{
nonzero_constraints.clear();
}
// @sect4{StationaryNavierStokes::initialize_system}
- // On each mesh the sparsity pattern and the size of the linear system
+ // On each mesh the SparsityPattern and the size of the linear system
// are different. This function initializes them after mesh refinement.
-
template <int dim>
void StationaryNavierStokes<dim>::initialize_system()
{
// This function builds the system matrix and right hand side that we
// currently work on. The @p initial_step argument is used to determine
// which set of constraints we apply (nonzero for the initial step and zero
- // for the others). The @p assemble_matrix flag determines whether to
+ // for the others). The @p assemble_matrix argument determines whether to
// assemble the whole system or only the right hand side vector,
// respectively.
-
template <int dim>
void StationaryNavierStokes<dim>::assemble(const bool initial_step,
const bool assemble_matrix)
if (assemble_matrix)
{
// Finally we move pressure mass matrix into a separate matrix:
-
pressure_mass_matrix.reinit(sparsity_pattern.block(1, 1));
pressure_mass_matrix.copy_from(system_matrix.block(1, 1));
const AffineConstraints<double> &constraints_used =
initial_step ? nonzero_constraints : zero_constraints;
- SolverControl solver_control(system_matrix.m(),
+ SolverControl solver_control(system_matrix.m(),
1e-4 * system_rhs.l2_norm(),
true);
- SolverFGMRES<BlockVector<double>> gmres(solver_control);
- SparseILU<double> pmass_preconditioner;
+ SolverFGMRES<BlockVector<double>> gmres(solver_control);
+ SparseILU<double> pmass_preconditioner;
pmass_preconditioner.initialize(pressure_mass_matrix,
SparseILU<double>::AdditionalData());
triangulation.execute_coarsening_and_refinement();
// First the DoFHandler is set up and constraints are generated. Then we
- // create a temporary vector "tmp", whose size is according with the
- // solution on the new mesh.
+ // create a temporary BlockVector <code>tmp</code>, whose size is
+ // according with the solution on the new mesh.
setup_dofs();
BlockVector<double> tmp(dofs_per_block);
// @sect4{StationaryNavierStokes::process_solution}
// In our test case, we do not know the analytical solution. This function
- // outputs the velocity components along x=0.5 and y from 0 to 1 so they
+ // outputs the velocity components along $x=0.5$ and $0 \leq y \leq 1$ so they
// can be compared with data from the literature.
template <int dim>
void StationaryNavierStokes<dim>::process_solution(unsigned int refinement)
}
}
-
// @sect4{StationaryNavierStokes::run}
// This is the last step of this program. In this part, we generate the grid
// and run the other functions respectively. The max refinement can be set by
const double Re = 1.0 / viscosity;
- // If the viscosity is smaller than 1/1000, we have to first search for an
+ // If the viscosity is smaller than $1/1000$, we have to first search for an
// initial guess via a continuation method. What we should notice is the
// search is always on the initial mesh, that is the $8 \times 8$ mesh in
// this program. After that, we just do the same as we did when viscosity
- // is larger than 1/1000: run Newton's iteration, refine the mesh,
+ // is larger than $1/1000$: run Newton's iteration, refine the mesh,
// transfer solutions, and repeat.
if (Re > 1000.0)
{
projects / dealii.git / commitdiff