From: Liang Zhao Date: Mon, 8 Aug 2016 14:56:34 +0000 (-0400) Subject: Introduce step-57 X-Git-Tag: v8.5.0-rc1~227^2~1 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=ef96bf86db38155af85bbe161219097f22f077e6;p=dealii.git Introduce step-57 --- diff --git a/doc/doxygen/tutorial/tutorial.h.in b/doc/doxygen/tutorial/tutorial.h.in index 20eb41c08b..9cd6ab19cf 100644 --- a/doc/doxygen/tutorial/tutorial.h.in +++ b/doc/doxygen/tutorial/tutorial.h.in @@ -419,6 +419,12 @@ * step-56 * Geometric Multigrid for Stokes. * + * + * + * step-57 + * Incompressible, stationary Navier Stokes equations. + * + * * * * @@ -503,7 +509,8 @@ * step-43, * step-44, * step-55, - * step-56 + * step-56, + * step-57 * * * @@ -602,7 +609,8 @@ * step-32, * step-33, * step-42, - * step-43 + * step-43, + * step-57 * * * @@ -740,7 +748,8 @@ * step-32, * step-43, * step-55, - * step-56 + * step-56, + * step-57 * * * @@ -757,7 +766,8 @@ * step-33, * step-41, * step-42, - * step-44 + * step-44, + * step-57 * * * @@ -859,7 +869,8 @@ * step-35, * step-46, * step-55, - * step-56 + * step-56, + * step-57 * * * @@ -942,7 +953,8 @@ * step-32, * step-35, * step-55, - * step-56 + * step-56, + * step-57 * * * diff --git a/doc/news/changes/major/20170117LiangZhaoTimoHeister b/doc/news/changes/major/20170117LiangZhaoTimoHeister new file mode 100644 index 0000000000..ae5d789d54 --- /dev/null +++ b/doc/news/changes/major/20170117LiangZhaoTimoHeister @@ -0,0 +1,4 @@ +New: The tutorial step-57 shows how to solve the stationary Navier-Stokes +equations using Newton's method. +
+(Liang Zhao, Timo Heister, 2017/01/17) diff --git a/examples/step-57/CMakeLists.txt b/examples/step-57/CMakeLists.txt new file mode 100644 index 0000000000..4cca8afbdd --- /dev/null +++ b/examples/step-57/CMakeLists.txt @@ -0,0 +1,51 @@ +## +# CMake script for the step-57 tutorial program: +## + +# Set the name of the project and target: +SET(TARGET "step-57") + +# Declare all source files the target consists of. Here, this is only +# the one step-X.cc file, but as you expand your project you may wish +# to add other source files as well. If your project becomes much larger, +# you may want to either replace the following statement by something like +# FILE(GLOB_RECURSE TARGET_SRC "source/*.cc") +# FILE(GLOB_RECURSE TARGET_INC "include/*.h") +# SET(TARGET_SRC ${TARGET_SRC} ${TARGET_INC}) +# or switch altogether to the large project CMakeLists.txt file discussed +# in the "CMake in user projects" page accessible from the "User info" +# page of the documentation. +SET(TARGET_SRC + ${TARGET}.cc + ) + +# Usually, you will not need to modify anything beyond this point... + +CMAKE_MINIMUM_REQUIRED(VERSION 2.8.8) + +FIND_PACKAGE(deal.II 8.5.0 QUIET + HINTS ${deal.II_DIR} ${DEAL_II_DIR} ../ ../../ $ENV{DEAL_II_DIR} + ) +IF(NOT ${deal.II_FOUND}) + MESSAGE(FATAL_ERROR "\n" + "*** Could not locate a (sufficiently recent) version of deal.II. ***\n\n" + "You may want to either pass a flag -DDEAL_II_DIR=/path/to/deal.II to cmake\n" + "or set an environment variable \"DEAL_II_DIR\" that contains this path." + ) +ENDIF() + +# +# Are all dependencies fulfilled? +# +IF(NOT DEAL_II_WITH_UMFPACK) + MESSAGE(FATAL_ERROR " +Error! The deal.II library found at ${DEAL_II_PATH} was not configured with + DEAL_II_WITH_UMFPACK = ON +One or all of these are OFF in your installation but are required for this tutorial step." + ) +ENDIF() + + +DEAL_II_INITIALIZE_CACHED_VARIABLES() +PROJECT(${TARGET}) +DEAL_II_INVOKE_AUTOPILOT() diff --git a/examples/step-57/doc/builds-on b/examples/step-57/doc/builds-on new file mode 100644 index 0000000000..f8c8247063 --- /dev/null +++ b/examples/step-57/doc/builds-on @@ -0,0 +1 @@ +step-15 step-22 diff --git a/examples/step-57/doc/gnuplot.gpl b/examples/step-57/doc/gnuplot.gpl new file mode 100644 index 0000000000..f41fd83c04 --- /dev/null +++ b/examples/step-57/doc/gnuplot.gpl @@ -0,0 +1,65 @@ +# Gnuplot script to plot velocity profiles together +# with reference data. + +set term svg size 800,900 dynamic +set style fill solid 0.25 border + +set size 1.0,1.0 +set tmargin at screen 0.95 +set bmargin at screen 0.10 +set rmargin at screen 0.95 +set lmargin at screen 0.10 + +set style line 1 lt 1 lc rgb '#0c0887' # blue +set style line 2 lt 1 lc rgb '#4b03a1' # purple-blue +set style line 3 lt 1 lc rgb '#7d03a8' # purple +set style line 4 lt 1 lc rgb '#a82296' # purple +set style line 5 lt 1 lc rgb '#cb4679' # magenta +set style line 6 lt 1 lc rgb '#e56b5d' # red +set style line 7 lt 1 lc rgb '#f89441' # orange +set style line 8 lt 1 lc rgb '#fdc328' # orange +set style line 9 lt 1 lc rgb '#f0f921' # yellow +set style line 10 lt 1 lc rgb '#00a000' # dark green +set style line 11 lt 1 lc rgb '#d00000' # dark red + + + +set output "step-57.compare-Re400.svg" +set xlabel "horizonal velocity" +set ylabel "y coordinate" +set xrange [-0.5: 1.0] +set title "horizontal velocity at x=0.5, Re=400" +set key right bottom +plot "400-line-4.txt" using 2:1 w l title "Numerical Solution", \ +"ref_2d_ghia_u.txt" using 3:1 ls 10 pt 4 w p title "Reference [Ghia 82]" + + + + +set output "step-57.compare-Re7500.svg" +set xlabel "horizonal velocity" +set ylabel "y coordinate" +set xrange [-0.5: 1.0] +set title "horizontal velocity at x=0.5, Re=7500" +set key right bottom +plot "7500-line-4.txt" using 2:1 w l title "Numerical Solution", \ +"ref_2d_ghia_u.txt" using 7:1 ls 10 pt 4 w p title "Reference [Ghia 82]", \ +"ref_2d_erturk_u.txt" using 5:1 ls 11 pt 6 w p title "Reference [Erturk]" + + + +set output "step-57.converge-Re7500.svg" +set xlabel "horizonal velocity" +set ylabel "y coordinate" +set xrange [-0.5: 1.0] +set title "horizontal velocity at x=0.5, Re=7500" +set key right bottom +plot "7500-line-0.txt" using 2:1 w l ls 1 title "Refinement 0", \ +"7500-line-1.txt" using 2:1 w l ls 2 title "Refinement 1", \ +"7500-line-2.txt" using 2:1 w l ls 3 title "Refinement 2", \ +"7500-line-3.txt" using 2:1 w l ls 4 title "Refinement 3", \ +"7500-line-4.txt" using 2:1 w l ls 5 title "Refinement 4", \ +"ref_2d_ghia_u.txt" using 7:1 ls 10 pt 4 w p title "Reference [Ghia 82]", \ +"ref_2d_erturk_u.txt" using 5:1 ls 11 pt 6 w p title "Reference [Erturk]" + + diff --git a/examples/step-57/doc/intro.dox b/examples/step-57/doc/intro.dox new file mode 100644 index 0000000000..be97c6f08d --- /dev/null +++ b/examples/step-57/doc/intro.dox @@ -0,0 +1,334 @@ +
+ +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. + + + +

Introduction

+ +

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 +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 + +@f{eqnarray*} +- \nu \Delta\textbf{u} + (\textbf{u} \cdot \nabla)\textbf{u} + \nabla p &=& \textbf{f}\\ +- \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 +\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

+ +Moving the right-hand side terms to the left, a nonlinear function is created as + +@f{eqnarray*} +F(\mathbf{u}, p) = \left( + \begin{array}{c} + - \nu \Delta\mathbf{u} + (\mathbf{u} \cdot \nabla)\mathbf{u} + \nabla p - \mathbf{f} \\ + - \nabla \cdot \mathbf{u} \\ + \end{array} + \right). +@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 +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} +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 +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 +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}), +@f} +where $x^{k+1}=x^{k}+\delta x^{k}$. + +Then we can evaluate the update term by solving the system +@f{eqnarray*} +\nabla F(\textbf{x}^{k}) \delta \textbf{x}^{k} = -F(\textbf{x}^{k}). +@f} +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 + +@f{eqnarray*} + & &\nabla F(\mathbf{u}^{k}, p^{k}) (\delta \mathbf{u}^{k}, \delta p^{k}) \\ + \\ + &=& \lim_{\epsilon \to 0} \frac{1}{\epsilon} (F(\mathbf{u}^{k}+\epsilon \delta \mathbf{u}^{k}, p^{k}+\epsilon\nabla\delta p^{k}) - (F(\mathbf{u}^{k}, p^{k}))\\ + \\ + &=& \lim_{\epsilon \to 0} \frac{1}{\epsilon} \left( + \begin{array}{c} + - \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{array} + \right)\\ + \\ + &=& \left( + \begin{array}{c} + - \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{array} + \right). +@f} + +Therefore, we arrive at the linearized system: + +@f{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} = \mathbf{g}, \\ +- \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 +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. + +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

+ +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. + +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 NSE 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 + +@f{eqnarray*} +- \nu_{1} \Delta\textbf{u} + \nabla p &=& \textbf{f}\\ +- \nabla \cdot \textbf{u} &=& 0 +@f} + +to get the initial guess for + +@f{eqnarray*} +- \nu_{1} \Delta\textbf{u} + (\textbf{u} \cdot \nabla)\textbf{u} + \nabla p &=& \textbf{f},\\ +- \nabla \cdot \textbf{u} &=& 0, +@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. + +Then the solution to + +@f{eqnarray*} +- \nu_{i} \Delta\textbf{u} + (\textbf{u} \cdot \nabla)\textbf{u} + \nabla p &=& \textbf{f},\\ +- \nabla \cdot \textbf{u} &=& 0. +@f} + +acts as the initial guess for + +@f{eqnarray*} +- \nu_{i+1} \Delta\textbf{u} + (\textbf{u} \cdot \nabla)\textbf{u} + \nabla p &=& \textbf{f},\\ +- \nabla \cdot \textbf{u} &=& 0. +@f} + +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 +@f{eqnarray*} + \left( + \begin{array}{cc} + A & B^{T} \\ + B & 0 \\ + \end{array} + \right) + \left( + \begin{array}{c} + U \\ + P \\ + \end{array} + \right) + = + \left( + \begin{array}{c} + F \\ + 0 \\ + \end{array} + \right). +@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. +Instead of solving the above system, we can solve the equivalent system + +@f{eqnarray*} + \left( + \begin{array}{cc} + A + \gamma B^TW^{-1}B & B^{T} \\ + B & 0 \\ + \end{array} + \right) + \left( + \begin{array}{c} + U \\ + P \\ + \end{array} + \right) + = + \left( + \begin{array}{c} + F \\ + 0 \\ + \end{array} + \right) +@f} +with a parameter $\gamma$ and an invertible matrix W. Here +$\gamma B^TW^{-1}B & B^{T}$ is the Augmented Lagrangian term and + 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 \\ + 0 & \tilde{S} \end{array}\right)^{-1}, +@f} + +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 + +@f{eqnarray*} +\tilde{S}^{-1} \approx -(\nu+\gamma)M_p^{-1}. +@f} +See [1] for details. + +We decompose $P^{-1}$ as + +@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) +\left(\begin{array}{cc} I & 0 \\ 0 & \tilde{S}^{-1} \end{array}\right). +@f} + +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 +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 [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 +$f=0$. The boundary condition is + +@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}. +@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 +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 +guess. Therefore, the nonlinear error is always brought below the tolerance of +Newton's iteration and the discretization error is reduced with each mesh +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, +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 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 +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. + +

Reference

+
    + +
  1. An Augmented Lagrangian-Based Approach to the Oseen Problem, M. Benzi and M. Olshanskii, SIAM J. SCI. COMPUT. 2006 +
  2. Efficient augmented Lagrangian-type preconditioning for the Oseen problem using Grad-Div stabilization, Timo Heister and Gerd Rapin +
  3. http://www.cfd-online.com/Wiki/Lid-driven_cavity_problem +
  4. High-Re solution for incompressible flow using the Navier-Stokes Equations and a Multigrid Method, U. Ghia, K. N. Ghia, and C. T. Shin +
  5. Numerical solutions of 2-D steady incompressible driven cavity flow at high Reynolds numbers, E. Erturk, T.C. Corke and C. Gokcol +
  6. Implicit Weighted ENO Schemes for the Three-Dimensional Incompressible Navier-Stokes Equations, Yang et al, 1998 +
  7. The 2D lid-driven cavity problem revisited, C. Bruneau and M. Saad, 2006 +
  8. https://en.wikipedia.org/wiki/Reynolds_number +
diff --git a/examples/step-57/doc/kind b/examples/step-57/doc/kind new file mode 100644 index 0000000000..e62f4e7222 --- /dev/null +++ b/examples/step-57/doc/kind @@ -0,0 +1 @@ +fluids diff --git a/examples/step-57/doc/results.dox b/examples/step-57/doc/results.dox new file mode 100644 index 0000000000..f1475d646d --- /dev/null +++ b/examples/step-57/doc/results.dox @@ -0,0 +1,430 @@ +

Results

+ +Now we use the method we discussed above to solve Navier Stokes equations with +viscosity 1/400 and 1/7500. + +

Test case 1: Re=400

+ +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. + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +
 Mesh0Mesh1Mesh2Mesh3Mesh4
Newton's iter Residual FGMRES Residual FGMRES Residual FGMRES Residual FGMRES Residual FGMRES
1 7.40396e-33 1.05562e-3 3 4.94796e-4 3 2.5624e-4 2 1.26733e-4 2
2 3.86766e-3 4 1.3549e-5 3 1.41981e-6 3 1.29108e-6 4 6.14794e-7 4
3 1.60421e-34 1.24836e-9 3 9.11557e-11 3 3.35933e-11 3 5.86734e-11 2
4 9.26748e-4 4 2.75537e-14 4 1.39986e-14 5 2.18864e-14 5 3.38787e-14 5
5 1.34601e-54                
6 2.5235e-8 5                
7 1.38899e-12 4                
8 4.68224e-15 4                
+ +The following figures show the sequence of the 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. + + + + + + + + + + + +
+ + + + + +
+ + + +
+ +This picture is the graphical streamline result of lid-driven cavity with Re=400. + + +Then the solution is compared with a reference solution +from [4] and the reference solution data can be found in the file "ref_2d_ghia_u.txt". + + + +

Test case 2: Re=7500

+ +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 +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 +in further iterations. + + + +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 +iterations are executed for solving this test case. + + + +Also we show the residual from each step of Newton's iteration on every +mesh. The quadratic convergence is shown clearly in the table. + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +
 Mesh0Mesh1Mesh2Mesh3Mesh4
Newton's iter Residual FGMRES Residual FGMRES Residual FGMRES Residual FGMRES Residual FGMRES
1 1.89223e-6 6 4.2506e-3 3 1.42993e-3 3 4.87932e-4 2 1.89981e-04 2
2 3.16439e-98 1.3732e-3 7 4.15062e-4 7 9.11191e-5 8 1.35553e-58
3 1.7628e-149 2.19455e-4 6 1.78805e-5 6 5.26782e-7 7 9.37391e-9 7
4     8.82693e-6 6 6.82096e-9 7 2.27696e-11 8 1.25899e-139
5     1.29739e-77 1.25167e-13 9 1.76128e-14 10    
6     4.43518e-117            
7     6.42323e-15 9           
+The sequence of generated grids looks like this: + + + + + + + + + + +
+ + + + + +
+ + + +
+We compare our solution with reference solution from [5]. + +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 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 +meshes. + + + + + +

Possibilities for extensions

+ +

Compare to other solvers

+ +It is easy to compare the currently implemented linear solver to just using +UMFPACK for the whole linear system. You need to remove the nullspace +containing the constant pressures and it is done in step-56. More interesting +is the comparison to other state of the art preconditioners like PCD. It turns +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 +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. + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +
RefDoFsIterative: Total/s (Setup/s)Direct: Total/s (Setup/s)
595390.10 (0.06)0.13 (0.12)
6375070.58 (0.37)1.03 (0.97)
71487393.59 (2.73)7.78 (7.53)
859238729.17 (24.94)(>4GB RAM)
+ + +

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 +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/doc/tooltip b/examples/step-57/doc/tooltip new file mode 100644 index 0000000000..25b3ae43f8 --- /dev/null +++ b/examples/step-57/doc/tooltip @@ -0,0 +1 @@ +The Navier Stokes equations via Newton's iteration diff --git a/examples/step-57/ref_2d_erturk_u.txt b/examples/step-57/ref_2d_erturk_u.txt new file mode 100644 index 0000000000..3af17f6ead --- /dev/null +++ b/examples/step-57/ref_2d_erturk_u.txt @@ -0,0 +1,25 @@ +# x-velocity along vertical line with x=0.5 from Erturk,Corke 2005 +# y u_x@Re=1000 2500 5000 7500 10'000 12'500 15'000 17'500 20'000 21'000 +1.000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 +0.990 0.8486 0.7704 0.6866 0.6300 0.5891 0.5587 0.5358 0.5183 0.5048 0.5003 +0.980 0.7065 0.5924 0.5159 0.4907 0.4837 0.4833 0.4850 0.4871 0.4889 0.4895 +0.970 0.5917 0.4971 0.4749 0.4817 0.4891 0.4941 0.4969 0.4982 0.4985 0.4983 +0.960 0.5102 0.4607 0.4739 0.4860 0.4917 0.4937 0.4937 0.4925 0.4906 0.4897 +0.950 0.4582 0.4506 0.4738 0.4824 0.4843 0.4833 0.4811 0.4784 0.4754 0.4742 +0.940 0.4276 0.4470 0.4683 0.4723 0.4711 0.4684 0.4653 0.4622 0.4592 0.4580 +0.930 0.4101 0.4424 0.4582 0.4585 0.4556 0.4523 0.4492 0.4463 0.4436 0.4425 +0.920 0.3993 0.4353 0.4452 0.4431 0.4398 0.4366 0.4338 0.4312 0.4287 0.4277 +0.910 0.3913 0.4256 0.4307 0.4275 0.4243 0.4216 0.4190 0.4166 0.4142 0.4132 +0.900 0.3838 0.4141 0.4155 0.4123 0.4095 0.4070 0.4047 0.4024 0.4001 0.3992 +0.500 -0.0620 -0.0403 -0.0319 -0.0287 -0.0268 -0.0256 -0.0247 -0.0240 -0.0234 -0.0232 +0.200 -0.3756 -0.3228 -0.3100 -0.3038 -0.2998 -0.2967 -0.2942 -0.2920 -0.2899 -0.2892 +0.180 -0.3869 -0.3439 -0.3285 -0.3222 -0.3179 -0.3146 -0.3119 -0.3096 -0.3074 -0.3066 +0.160 -0.3854 -0.3688 -0.3467 -0.3406 -0.3361 -0.3326 -0.3297 -0.3271 -0.3248 -0.3239 +0.140 -0.3690 -0.3965 -0.3652 -0.3587 -0.3543 -0.3506 -0.3474 -0.3446 -0.3422 -0.3412 +0.120 -0.3381 -0.4200 -0.3876 -0.3766 -0.3721 -0.3685 -0.3652 -0.3622 -0.3595 -0.3585 +0.100 -0.2960 -0.4250 -0.4168 -0.3978 -0.3899 -0.3859 -0.3827 -0.3797 -0.3769 -0.3758 +0.080 -0.2472 -0.3979 -0.4419 -0.4284 -0.4142 -0.4054 -0.4001 -0.3965 -0.3936 -0.3925 +0.060 -0.1951 -0.3372 -0.4272 -0.4491 -0.4469 -0.4380 -0.4286 -0.4206 -0.4143 -0.4121 +0.040 -0.1392 -0.2547 -0.3480 -0.3980 -0.4259 -0.4407 -0.4474 -0.4490 -0.4475 -0.4463 +0.020 -0.0757 -0.1517 -0.2223 -0.2633 -0.2907 -0.3113 -0.3278 -0.3412 -0.3523 -0.3562 +0.000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 diff --git a/examples/step-57/ref_2d_ghia_u.txt b/examples/step-57/ref_2d_ghia_u.txt new file mode 100644 index 0000000000..5e27d978cf --- /dev/null +++ b/examples/step-57/ref_2d_ghia_u.txt @@ -0,0 +1,19 @@ +# x-velocity along vertical line with x=0.5 from Ghia 1982 +# y u_x@Re=100 400 1000 3200 5000 7500 10000 +1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 +0.9766 0.8412 0.7584 0.6593 0.5324 0.4822 0.4724 0.4722 +0.9688 0.7887 0.6844 0.5749 0.4830 0.4612 0.4705 0.4778 +0.9609 0.7372 0.6176 0.5112 0.4655 0.4599 0.4732 0.4807 +0.9531 0.6872 0.5589 0.4660 0.4610 0.4604 0.4717 0.4780 +0.8516 0.2315 0.2909 0.3330 0.3468 0.3356 0.3423 0.3464 +0.7344 0.0033 0.1626 0.1872 0.1979 0.2009 0.2059 0.2067 +0.6172 -0.1364 0.0214 0.0570 0.0716 0.0818 0.0834 0.0834 +0.5000 -0.2058 -0.1148 -0.0608 -0.0427 -0.0304 -0.0380 -0.0311 +0.4531 -0.2109 -0.1712 -0.1065 -0.8664 -0.0740 -0.0750 -0.0754 +0.2813 -0.1566 -0.3273 -0.2781 -0.2443 -0.2286 -0.2318 -0.2319 +0.1719 -0.1015 -0.2430 -0.3829 -0.3432 -0.3305 -0.3239 -0.3271 +0.1016 -0.0643 -0.1461 -0.2973 -0.4193 -0.4044 -0.3832 -0.3800 +0.0703 -0.0478 -0.1034 -0.2222 -0.3783 -0.4364 -0.4303 -0.4166 +0.0625 -0.0419 -0.0927 -0.2020 -0.3534 -0.4290 -0.4359 -0.4254 +0.0547 -0.0372 -0.0819 -0.1811 -0.3241 -0.4117 -0.4315 -0.4274 +0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 diff --git a/examples/step-57/step-57.cc b/examples/step-57/step-57.cc new file mode 100644 index 0000000000..468df2e0f2 --- /dev/null +++ b/examples/step-57/step-57.cc @@ -0,0 +1,898 @@ +/* --------------------------------------------------------------------- + * + * Copyright (C) 2008 - 2016 by the deal.II authors + * + * This file is part of the deal.II library. + * + * The deal.II library is free software; you can use it, redistribute + * it, and/or modify it under the terms of the GNU Lesser General + * Public License as published by the Free Software Foundation; either + * version 2.1 of the License, or (at your option) any later version. + * The full text of the license can be found in the file LICENSE at + * the top level of the deal.II distribution. + * + * --------------------------------------------------------------------- + * + * Author: Liang Zhao and Timo Heister, Clemson University, 2016 + */ + +// @sect3{Include files} + +// As usual, we start by including some well-known files: +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include + +#include +#include +#include + +#include +#include +#include +#include + +// To transfer solutions between meshes, this file is included: +#include + +// This file includes UMFPACK: the direct solver: +#include + +// And the one for ILU preconditioner: +#include + + +#include +#include +#include + +namespace Step57 +{ + using namespace dealii; + + // @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 ConstraintMatrix for Dirichlet boundary + // conditions at the initial step and a zero ConstraintMatrix for the Newton + // is defined by 1/Re which has been discussed in the introduction. + + template + class StationaryNavierStokes + { + public: + StationaryNavierStokes(const unsigned int degree); + void run(const unsigned int refinement); + + private: + void setup_dofs(); + void initialize_system(); + void assemble_system(const bool initial_step, + const bool assemble_matrix, + const bool assemble_rhs); + void assemble_matrix(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; + double gamma; + const unsigned int degree; + std::vector dofs_per_block; + + Triangulation triangulation; + FESystem fe; + DoFHandler dof_handler; + + ConstraintMatrix zero_constraints; + ConstraintMatrix nonzero_constraints; + + BlockSparsityPattern sparsity_pattern; + BlockSparseMatrix system_matrix; + SparseMatrix pressure_mass_matrix; + + BlockVector present_solution; + BlockVector newton_update; + BlockVector system_rhs; + 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 ConstraintMatrix 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: + template + class BoundaryValues : public Function + { + public: + BoundaryValues() : Function(dim+1) {} + virtual double value(const Point &p, + const unsigned int component) const; + + virtual void vector_value(const Point &p, + Vector &values) const; + }; + + template + double BoundaryValues::value(const Point &p, + const unsigned int component) const + { + Assert (component < this->n_components, + ExcIndexRange (component, 0, this->n_components)); + if (component == 0 && std::abs(p[dim-1]-1.0) < 1e-10) + return 1.0; + + return 0; + } + + template + void BoundaryValues::vector_value ( const Point &p, + Vector &values ) const + { + for (unsigned int c = 0; c < this->n_components; ++c) + values(c) = BoundaryValues::value (p, c); + } + + + // @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. + // + template + class BlockSchurPreconditioner : public Subscriptor + { + public: + BlockSchurPreconditioner (double gamma, + double viscosity, + const BlockSparseMatrix &S, + const SparseMatrix &P, + const PreconditionerMp &Mppreconditioner); + + void vmult (BlockVector &dst, + const BlockVector &src) const; + + private: + const double gamma; + const double viscosity; + const BlockSparseMatrix &stokes_matrix; + const SparseMatrix &pressure_mass_matrix; + const PreconditionerMp &mp_preconditioner; + 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. + + template + BlockSchurPreconditioner:: + BlockSchurPreconditioner (double gamma, + double viscosity, + const BlockSparseMatrix &S, + const SparseMatrix &P, + const PreconditionerMp &Mppreconditioner) + : + gamma (gamma), + viscosity (viscosity), + stokes_matrix (S), + pressure_mass_matrix (P), + mp_preconditioner (Mppreconditioner) + { + A_inverse.initialize(stokes_matrix.block(0,0)); + } + + template + void + BlockSchurPreconditioner:: + vmult (BlockVector &dst, + const BlockVector &src) const + { + Vector utmp(src.block(0)); + + { + SolverControl solver_control(1000, 1e-6 * src.block(1).l2_norm()); + SolverCG<> cg (solver_control); + + dst.block(1) = 0.0; + cg.solve(pressure_mass_matrix, + dst.block(1), src.block(1), + mp_preconditioner); + dst.block(1) *= -(viscosity+gamma); + } + + { + stokes_matrix.block(0,1).vmult(utmp, dst.block(1)); + utmp *= -1.0; + utmp += src.block(0); + } + + A_inverse.vmult (dst.block(0), utmp); + } + + // @sect3{StationaryNavierStokes class implementation} + // @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. + // + + template + StationaryNavierStokes::StationaryNavierStokes(const unsigned int degree) + : + viscosity(1.0/7500.0), + gamma(1.0), + degree(degree), + triangulation(Triangulation::maximum_smoothing), + fe(FE_Q(degree+1), dim, + FE_Q(degree), 1), + dof_handler(triangulation) + {} + + // @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() + { + system_matrix.clear(); + pressure_mass_matrix.clear(); + + // The first step is to associate DoFs with a given mesh. + dof_handler.distribute_dofs (fe); + + // 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); + + dofs_per_block.resize (2); + DoFTools::count_dofs_per_block (dof_handler, dofs_per_block, block_component); + unsigned int dof_u = dofs_per_block[0]; + unsigned int dof_p = dofs_per_block[1]; + + // In Newton's scheme, we first apply the boundary condition on the solution + // obtained from the initial step. To make sure the boundary conditions 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(); + + DoFTools::make_hanging_node_constraints(dof_handler, nonzero_constraints); + VectorTools::interpolate_boundary_values(dof_handler, + 0, + BoundaryValues(), + nonzero_constraints, + fe.component_mask(velocities)); + } + nonzero_constraints.close(); + + { + zero_constraints.clear(); + + DoFTools::make_hanging_node_constraints(dof_handler, zero_constraints); + VectorTools::interpolate_boundary_values(dof_handler, + 0, + ZeroFunction(dim+1), + zero_constraints, + fe.component_mask(velocities)); + } + zero_constraints.close(); + + std::cout << " Number of active cells: " + << triangulation.n_active_cells() + << std::endl + << " Number of degrees of freedom: " + << dof_handler.n_dofs() + << " (" << dof_u << '+' << dof_p << ')' + << std::endl; + } + + // @sect4{StationaryNavierStokes::initialize_system} + // On each mesh the sparsity pattern and the size of the linear system + // are different. This function initializes them after mesh refinement. + + template + void StationaryNavierStokes::initialize_system() + { + { + BlockDynamicSparsityPattern dsp (dofs_per_block, dofs_per_block); + DoFTools::make_sparsity_pattern (dof_handler, dsp, nonzero_constraints); + sparsity_pattern.copy_from (dsp); + } + + system_matrix.reinit (sparsity_pattern); + + present_solution.reinit (dofs_per_block); + newton_update.reinit (dofs_per_block); + system_rhs.reinit (dofs_per_block); + } + + // @sect4{StationaryNavierStokes::assemble_system} + + // This function builds the system matrix and right hand side that we + // actually work on. "initial_step" is given for applying different + // constraints (nonzero for the initial step and zero for the others). The + // other two flags are to determine whether to assemble the system matrix + // or the right hand side vector, respectively. + + template + void StationaryNavierStokes::assemble_system(const bool initial_step, + const bool assemble_matrix, + const bool assemble_rhs) + { + if (assemble_matrix) + system_matrix = 0; + + if (assemble_rhs) + system_rhs = 0; + + QGauss quadrature_formula(degree+2); + + FEValues fe_values (fe, + quadrature_formula, + update_values | + update_quadrature_points | + update_JxW_values | + update_gradients ); + + const unsigned int dofs_per_cell = fe.dofs_per_cell; + const unsigned int n_q_points = quadrature_formula.size(); + + const FEValuesExtractors::Vector velocities (0); + const FEValuesExtractors::Scalar pressure (dim); + + FullMatrix local_matrix (dofs_per_cell, dofs_per_cell); + Vector local_rhs (dofs_per_cell); + + std::vector local_dof_indices (dofs_per_cell); + + // For the linearized system, we create temporary storage for present velocity + // and gradient, and present pressure. In practice, they are + // all obtained through their shape functions at quadrature points. + + std::vector > present_velocity_values (n_q_points); + std::vector > present_velocity_gradients (n_q_points); + std::vector present_pressure_values (n_q_points); + + std::vector div_phi_u (dofs_per_cell); + std::vector > phi_u (dofs_per_cell); + std::vector > grad_phi_u (dofs_per_cell); + std::vector phi_p (dofs_per_cell); + + typename DoFHandler::active_cell_iterator + cell = dof_handler.begin_active(), + endc = dof_handler.end(); + + for (; cell!=endc; ++cell) + { + fe_values.reinit(cell); + + local_matrix = 0; + local_rhs = 0; + + fe_values[velocities].get_function_values(evaluation_point, + present_velocity_values); + + fe_values[velocities].get_function_gradients(evaluation_point, + present_velocity_gradients); + + fe_values[pressure].get_function_values(evaluation_point, + present_pressure_values); + + // The assembly is similar to step-22. An additional term with gamma as a coefficient + // is the Augmented Lagrangian (AL), which is assembled via grad-div stabilization. + // As we discussed in the introduction, the bottom right block of the system matrix should be + // zero. Since the pressure mass matrix is used while creating the preconditioner, + // we assemble it here and then move it into a separate SparseMatrix at the end (same as in step-22). + + for (unsigned int q=0; qget_dof_indices (local_dof_indices); + + const ConstraintMatrix &constraints_used = initial_step ? nonzero_constraints : zero_constraints; + // Finally we move pressure mass matrix into a separate matrix: + + if (assemble_matrix) + { + constraints_used.distribute_local_to_global(local_matrix, + local_dof_indices, + system_matrix); + } + + if (assemble_rhs) + { + constraints_used.distribute_local_to_global(local_rhs, + local_dof_indices, + system_rhs); + } + + } + + if (assemble_matrix) + { + pressure_mass_matrix.reinit(sparsity_pattern.block(1,1)); + pressure_mass_matrix.copy_from(system_matrix.block(1,1)); + system_matrix.block(1,1) = 0; + } + } + + template + void StationaryNavierStokes::assemble_matrix(const bool initial_step) + { + assemble_system(initial_step, true, false); + } + + template + void StationaryNavierStokes::assemble_rhs(const bool initial_step) + { + assemble_system(initial_step, false, true); + } + // @sect4{StationaryNavierStokes::solve} + // In this function, we use FGMRES together with the block preconditioner, + // which is defined at the beginning of the program, to solve the linear + // system. What we obtain at this step is the solution vector. If this is + // the initial step, the solution vector gives us an initial guess for the + // Navier Stokes equations. For the initial step, nonzero constraints are + // applied in order to make sure boundary conditions are satisfied. In the + // following steps, we will solve for the Newton update so zero + // constraints are used. + template + void StationaryNavierStokes::solve (const bool initial_step) + { + const ConstraintMatrix &constraints_used = initial_step ? nonzero_constraints : zero_constraints; + + SolverControl solver_control (system_matrix.m(),1e-4*system_rhs.l2_norm(), true); + SolverFGMRES > gmres(solver_control); + + SparseILU pmass_preconditioner; + pmass_preconditioner.initialize (pressure_mass_matrix, + SparseILU::AdditionalData()); + + const BlockSchurPreconditioner > + preconditioner (gamma, + viscosity, + system_matrix, + pressure_mass_matrix, + pmass_preconditioner); + + gmres.solve (system_matrix, + newton_update, + system_rhs, + preconditioner); + std::cout << " ****FGMRES steps: " << solver_control.last_step() << std::endl; + + constraints_used.distribute(newton_update); + } + + // @sect4{StationaryNavierStokes::refine_mesh} + // + // After finding a good initial guess on the coarse mesh, we hope to + // decrease the error through refining the mesh. Here we do adaptive + // refinement similar to step-15 except that we use the Kelly estimator on + // the velocity only. We also need to transfer the current solution to the + // next mesh using the SolutionTransfer class. + template + void StationaryNavierStokes::refine_mesh() + { + Vector estimated_error_per_cell (triangulation.n_active_cells()); + FEValuesExtractors::Vector velocity(0); + KellyErrorEstimator::estimate (dof_handler, + QGauss(degree+1), + typename FunctionMap::type(), + present_solution, + estimated_error_per_cell, + fe.component_mask(velocity)); + + GridRefinement::refine_and_coarsen_fixed_number (triangulation, + estimated_error_per_cell, + 0.3, 0.0); + + triangulation.prepare_coarsening_and_refinement(); + SolutionTransfer > solution_transfer(dof_handler); + solution_transfer.prepare_for_coarsening_and_refinement(present_solution); + 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. + setup_dofs(); + + BlockVector tmp (dofs_per_block); + + // Transfer solution from coarse to fine mesh and apply boundary value + // constraints to the new transfered solution. Note that present_solution + // is still a vector corresponding to the old mesh. + solution_transfer.interpolate(present_solution, tmp); + nonzero_constraints.distribute(tmp); + + // Finally set up matrix and vectors and set the present_solution to the + // interpolated data. + initialize_system(); + present_solution = tmp; + } + + // @sect4{StationaryNavierStokes::newton_iteration} + // + // 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. + + template + void StationaryNavierStokes::newton_iteration(const double tolerance, + const unsigned int max_iteration, + const unsigned int max_refinement, + const bool is_initial_step, + const bool output_result) + { + double current_res; + double last_res; + bool first_step = is_initial_step; + + for (unsigned int refinement = 0; refinement < max_refinement+1; ++refinement) + { + unsigned int outer_iteration = 0; + last_res = 1.0; + current_res = 1.0; + std::cout << "*****************************************" << std::endl; + std::cout << "************ refinement = " << refinement << " ************ " << std::endl; + std::cout << "viscosity= " << viscosity << std::endl; + std::cout << "*****************************************" << std::endl; + + while ((first_step || (current_res > tolerance)) && outer_iteration < max_iteration) + { + if (first_step) + { + setup_dofs(); + initialize_system(); + evaluation_point = present_solution; + assemble_matrix(first_step); + assemble_rhs(first_step); + solve(first_step); + present_solution = newton_update; + nonzero_constraints.distribute(present_solution); + first_step = false; + evaluation_point = present_solution; + assemble_rhs(first_step); + current_res = system_rhs.l2_norm(); + std::cout << "******************************" << std::endl; + std::cout << " The residual of initial guess is " << current_res << std::endl; + std::cout << " Initialization complete! " << std::endl; + last_res = current_res; + } + + else + { + evaluation_point = present_solution; + assemble_matrix(first_step); + if (outer_iteration == 0) + assemble_rhs(first_step); + 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. + + for (double alpha = 1.0; alpha > 1e-5; alpha *= 0.5) + { + evaluation_point = present_solution; + evaluation_point.add(alpha, newton_update); + nonzero_constraints.distribute(evaluation_point); + assemble_rhs(first_step); + current_res = system_rhs.l2_norm(); + std::cout << " alpha = " << std::setw(6) << alpha << std::setw(0) + << " res = " << current_res << std::endl; + if (current_res < last_res) + break; + } + { + present_solution = evaluation_point; + std::cout << " ---- Iteration " << outer_iteration << " residual: " << current_res << std::endl; + last_res = current_res; + } + + } + ++outer_iteration; + + if (output_result) + { + output_results (max_iteration*refinement+outer_iteration); + + if (current_res <= tolerance) + process_solution(refinement); + } + } + + if (refinement < max_refinement) + { + refine_mesh(); + std::cout << "*****************************************" << std::endl + << " Do refinement ------ " << std::endl; + } + } + + + } + + // @sect4{StationaryNavierStokes::compute_initial_guess} + // + // This function will provide us with an initial guess by using a + // continuation method as we discussed in the introduction. The Reynolds + // number is increased step-by-step until we reach the target value. By + // experiment, the solution to Stokes is good enough to be the initial guess + // of NSE with Reynolds number 1000 so we start there. To make sure the + // solution from previous problem is close enough to the next one, the step + // size must be small enough. + template + void StationaryNavierStokes::compute_initial_guess(double step_size) + { + const double target_Re = 1.0/viscosity; + + bool is_initial_step = true; + + for (double Re=1000.0; Re < target_Re; Re = std::min(Re+step_size, target_Re)) + { + viscosity = 1.0/Re; + std::cout << "*****************************************" << std::endl; + std::cout << " Searching for initial guess with Re = " << Re << std::endl; + std::cout << "*****************************************" << std::endl; + + newton_iteration(1e-12, 50, 0, is_initial_step, false); + is_initial_step = false; + } + } + + // @sect4{StationaryNavierStokes::output_results} + // This function is the same as in step-22. + template + void StationaryNavierStokes::output_results (const unsigned int output_index) const + { + std::vector solution_names (dim, "velocity"); + solution_names.push_back ("pressure"); + + std::vector + data_component_interpretation + (dim, DataComponentInterpretation::component_is_part_of_vector); + data_component_interpretation + .push_back (DataComponentInterpretation::component_is_scalar); + DataOut data_out; + data_out.attach_dof_handler (dof_handler); + data_out.add_data_vector (present_solution, solution_names, + DataOut::type_dof_data, + data_component_interpretation); + data_out.build_patches (); + + std::ostringstream filename; + filename << 1.0/viscosity + << "-solution-" + << Utilities::int_to_string (output_index, 4) + << ".vtk"; + + std::ofstream output (filename.str().c_str()); + data_out.write_vtk (output); + } + + // @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 + // can be compared with data from the literature. + template + void StationaryNavierStokes::process_solution(unsigned int refinement) + { + std::ostringstream filename; + filename << (1.0/viscosity) << "-line-" << refinement << ".txt"; + + std::ofstream f (filename.str().c_str()); + f << "# y u_x u_y" << std::endl; + + Point p; + p(0)= 0.5; + p(1)= 0.5; + + f << std::scientific; + + for (unsigned int i=0; i<=100; ++i) + { + p(dim-1) = i/100.0; + + Vector tmp_vector(dim+1); + VectorTools::point_value(dof_handler, present_solution, p, tmp_vector); + f << p(dim-1); + + for (int j=0; j + void StationaryNavierStokes::run(const unsigned int refinement) + { + + GridGenerator::hyper_cube(triangulation); + triangulation.refine_global(5); + + const double Re = 1.0/viscosity; + + // 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, + // transfer solutions, and repeat. + if (Re > 1000.0) + { + std::cout << " Searching for initial guess ... " << std::endl; + const double step_size = 2000.0; + compute_initial_guess(step_size); + std::cout << "*****************************************" << std::endl + << " Initial guess obtained " << std::endl + << " * " << std::endl + << " * " << std::endl + << " * " << std::endl + << " * " << std::endl + << "*****************************************" << std::endl; + + std::cout << " Computing solution with target Re = " << Re << std::endl; + viscosity = 1.0/Re; + newton_iteration(1e-12, 50, refinement, false, true); + } + else + { + // When the viscosity is larger than 1/1000, the solution to Stokes + // equations is good enough as an initial guess. If so, we do not need + // to search for the initial guess using a continuation + // method. Newton's iteration can be started directly. + + newton_iteration(1e-12, 50, refinement, true, true); + } + } +} + +int main() +{ + try + { + using namespace dealii; + using namespace Step57; + + StationaryNavierStokes<2> flow(/* degree = */1); + flow.run(4); + } + catch (std::exception &exc) + { + std::cerr << std::endl << std::endl + << "----------------------------------------------------" + << std::endl; + std::cerr << "Exception on processing: " << std::endl + << exc.what() << std::endl + << "Aborting!" << std::endl + << "----------------------------------------------------" + << std::endl; + return 1; + } + catch (...) + { + std::cerr << std::endl << std::endl + << "----------------------------------------------------" + << std::endl; + std::cerr << "Unknown exception!" << std::endl + << "Aborting!" << std::endl + << "----------------------------------------------------" + << std::endl; + return 1; + } + return 0; +}