From 3611e2dfd118ac79994b687da21d93ca17fad4b4 Mon Sep 17 00:00:00 2001 From: David Wells Date: Sat, 27 Oct 2018 13:20:17 -0400 Subject: [PATCH] step-57: Reword some sentences for clarity. --- examples/step-57/doc/intro.dox | 51 ++++++++------- examples/step-57/doc/results.dox | 2 +- examples/step-57/step-57.cc | 106 ++++++++++++++----------------- 3 files changed, 74 insertions(+), 85 deletions(-) diff --git a/examples/step-57/doc/intro.dox b/examples/step-57/doc/intro.dox index 24b3681bce..d16ed6a45b 100644 --- a/examples/step-57/doc/intro.dox +++ b/examples/step-57/doc/intro.dox @@ -16,7 +16,7 @@ Award No. EAR-0949446 and The University of California-Davis.

Navier Stokes Equations

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 @@ -25,11 +25,11 @@ satisfying @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 @@ -38,7 +38,7 @@ solve is necessary.

Linearization of Navier-Stokes Equations

-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) = @@ -48,32 +48,31 @@ 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} @@ -117,22 +116,21 @@ Therefore, we arrive at the linearized system: + \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: @@ -152,8 +150,9 @@ Now, Newton's iteration can be used to solve for the update terms:

Finding an Initial Guess

-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, @@ -298,8 +297,8 @@ 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. -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*} @@ -319,9 +318,9 @@ refinement. 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 diff --git a/examples/step-57/doc/results.dox b/examples/step-57/doc/results.dox index 98c6f78791..ad70d60b3e 100644 --- a/examples/step-57/doc/results.dox +++ b/examples/step-57/doc/results.dox @@ -326,7 +326,7 @@ 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 decreases as we perform Newton iterations, 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 diff --git a/examples/step-57/step-57.cc b/examples/step-57/step-57.cc index 34029fb3a3..112783b0e3 100644 --- a/examples/step-57/step-57.cc +++ b/examples/step-57/step-57.cc @@ -71,17 +71,14 @@ namespace Step57 // @sect3{The NavierStokesProblem 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 class StationaryNavierStokes { @@ -138,19 +135,16 @@ namespace Step57 BlockVector 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 dim+1. 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 class BoundaryValues : public Function { @@ -176,20 +170,18 @@ namespace Step57 // @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 BlockSchurPreconditioner : public Subscriptor { @@ -211,10 +203,10 @@ namespace Step57 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 BlockSchurPreconditioner::BlockSchurPreconditioner( @@ -264,9 +256,7 @@ namespace Step57 // @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. - // - + // gamma. template StationaryNavierStokes::StationaryNavierStokes(const unsigned int degree) : viscosity(1.0 / 7500.0) @@ -612,14 +602,14 @@ namespace Step57 // // 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 is_initial_step tells 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 + // determines whether or not graphical output should be produced. template void StationaryNavierStokes::newton_iteration( const double tolerance, @@ -666,11 +656,11 @@ namespace Step57 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 + // alpha 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; -- 2.39.5