From 8ce49b150f22108d32b181746b80c6d552dd7506 Mon Sep 17 00:00:00 2001 From: David Wells Date: Fri, 26 Oct 2018 22:14:59 -0400 Subject: [PATCH] step-57: Improve typography. 1. Remove trailing whitespace. 2. Fix some minor spelling mistakes. 3. Consistently use MathJax for equations. --- examples/step-57/doc/intro.dox | 93 ++++++++++++++++---------------- examples/step-57/doc/results.dox | 59 ++++++++++---------- examples/step-57/step-57.cc | 37 +++++++------ 3 files changed, 99 insertions(+), 90 deletions(-) diff --git a/examples/step-57/doc/intro.dox b/examples/step-57/doc/intro.dox index 55c31b2010..612aff39c3 100644 --- a/examples/step-57/doc/intro.dox +++ b/examples/step-57/doc/intro.dox @@ -1,7 +1,7 @@
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 @@ -33,7 +33,7 @@ nonlinear system 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 +using the solution from the last time step and no nonlinear solve is necessary.

Linearization of Navier-Stokes Equations

@@ -49,13 +49,13 @@ F(\mathbf{u}, p) = \left( \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, @@ -72,7 +72,7 @@ the update term as a whole, that is @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*} @@ -112,10 +112,10 @@ Therefore, we arrive at the linearized system: 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 @@ -124,14 +124,17 @@ derived like it is done in step-22. 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. -
    +
  5. Initialization: Initial guess $u_0$ and $p_0$, tolerance $\tau$;
    6. +
  6. Linear solve to compute update term $\delta\textbf{u}^{k}$ and + $\delta p^k$;
    7. +
  7. Update the approximation: + $\textbf{u}^{k+1} = \textbf{u}^{k} + \delta\textbf{u}^{k}$ and + $p^{k+1} = p^{k} + \delta p^{k}$;
    8. +
  8. 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

@@ -144,8 +147,8 @@ 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 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} - @@ -168,7 +171,7 @@ to get the initial guess for @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 @@ -184,11 +187,11 @@ acts as the initial guess for - \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. -

The Solver and Preconditioner

+

The %Solver and Preconditioner

At each step of Newton's iteration, the problem results in solving a saddle point systems of the form @@ -215,7 +218,7 @@ 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*} @@ -239,16 +242,15 @@ Instead of solving the above system, we can solve the equivalent system \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} @@ -264,7 +266,7 @@ See [1] for details. 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) @@ -274,7 +276,7 @@ P^{-1} = 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 @@ -283,19 +285,19 @@ M_p^{-1}B)_{ij}$, as explained in [2].

Test Case

-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 @@ -310,21 +312,22 @@ 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, 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. -

Reference

+

References

  1. An Augmented Lagrangian-Based Approach to the Oseen Problem, M. Benzi and M. Olshanskii, SIAM J. SCI. COMPUT. 2006 diff --git a/examples/step-57/doc/results.dox b/examples/step-57/doc/results.dox index 42ac92872e..98c6f78791 100644 --- a/examples/step-57/doc/results.dox +++ b/examples/step-57/doc/results.dox @@ -1,18 +1,19 @@

    Results

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

    Test case 1: Re=400

    +

    Test case 1: Low Reynolds Number

    -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. - + @@ -48,7 +49,7 @@ every mesh is shown. The data in the table shows that Newton's iteration converg - + @@ -61,7 +62,7 @@ every mesh is shown. The data in the table shows that Newton's iteration converg - + @@ -74,7 +75,7 @@ every mesh is shown. The data in the table shows that Newton's iteration converg - + @@ -87,7 +88,7 @@ every mesh is shown. The data in the table shows that Newton's iteration converg - + @@ -118,9 +119,10 @@ every mesh is shown. The data in the table shows that Newton's iteration converg 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.
    Re=400$\mathrm{Re}=400$ Mesh0 Mesh1 Mesh2
    2 7.0849e-045.0000e+005 9.9458e-04 5 1.1409e-04
    3 1.9980e-055.0000e+005 4.5007e-05 5 2.9020e-08
    4 2.3165e-096.0000e+006 1.6891e-07 5 1.2338e-14
    5 1.2585e-137.0000e+007 1.4520e-11 6 1.9044e-13
    @@ -144,7 +146,8 @@ the coarse mesh is interpolated to the fine mesh to be used as an initial guess.
    -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$. Then the solution is compared with a reference solution @@ -152,22 +155,22 @@ from [4] and the reference solution data can be found in the file "ref_2d_ghia_u -

    Test case 2: Re=7500

    +

    Test case 2: High Reynolds Number

    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. -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. @@ -177,7 +180,7 @@ mesh. The quadratic convergence is clearly visible in the table. - + @@ -323,11 +326,11 @@ The following picture presents the graphical result. 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. @@ -346,14 +349,14 @@ out that the preconditioner here is very competitive, as can be seen in the 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. -
    Re=7500$\mathrm{Re}=7500$ Mesh0 Mesh1 Mesh2
    +
    - + @@ -385,17 +388,17 @@ consumes less memory. This will be even more pronounced in 3d.
    RefRefinement Cycle DoFs Iterative: Total/s (Setup/s) Direct: Total/s (Setup/s)
    -

    3D computations

    +

    3d computations

    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. - +

    Parallelization

    -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. diff --git a/examples/step-57/step-57.cc b/examples/step-57/step-57.cc index a8d77db990..e2255fe977 100644 --- a/examples/step-57/step-57.cc +++ b/examples/step-57/step-57.cc @@ -96,19 +96,29 @@ namespace Step57 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; @@ -287,7 +297,6 @@ namespace Step57 // @sect4{StationaryNavierStokes::setup_dofs} // This function initializes the DoFHandler enumerating the degrees of freedom // and constraints on the current mesh. - template void StationaryNavierStokes::setup_dofs() { @@ -300,7 +309,6 @@ namespace Step57 // 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 block_component(dim + 1, 0); block_component[dim] = 1; DoFRenumbering::component_wise(dof_handler, block_component); @@ -317,7 +325,6 @@ namespace Step57 // 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(); @@ -351,9 +358,8 @@ namespace Step57 } // @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 void StationaryNavierStokes::initialize_system() { @@ -375,10 +381,9 @@ namespace Step57 // 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 void StationaryNavierStokes::assemble(const bool initial_step, const bool assemble_matrix) @@ -513,7 +518,6 @@ namespace Step57 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)); @@ -554,12 +558,12 @@ namespace Step57 const AffineConstraints &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> gmres(solver_control); - SparseILU pmass_preconditioner; + SolverFGMRES> gmres(solver_control); + SparseILU pmass_preconditioner; pmass_preconditioner.initialize(pressure_mass_matrix, SparseILU::AdditionalData()); @@ -608,8 +612,8 @@ namespace Step57 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 tmp, whose size is + // according with the solution on the new mesh. setup_dofs(); BlockVector tmp(dofs_per_block); @@ -799,7 +803,7 @@ namespace Step57 // @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 void StationaryNavierStokes::process_solution(unsigned int refinement) @@ -828,7 +832,6 @@ namespace Step57 } } - // @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 @@ -841,11 +844,11 @@ namespace Step57 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) { -- 2.39.5