]> https://gitweb.dealii.org/ - code-gallery.git/commitdiff
Time-dependent Navier-Stokes
authorJie Cheng <chengjiehust@gmail.com>
Mon, 13 Nov 2017 05:57:57 +0000 (00:57 -0500)
committerJie Cheng <chengjiehust@gmail.com>
Fri, 27 Apr 2018 04:29:19 +0000 (00:29 -0400)
time_dependent_navier_stokes/CMakeLists.txt [new file with mode: 0644]
time_dependent_navier_stokes/Readme.md [new file with mode: 0644]
time_dependent_navier_stokes/doc/Re100.mp4 [new file with mode: 0644]
time_dependent_navier_stokes/doc/Re100.png [new file with mode: 0644]
time_dependent_navier_stokes/doc/author [new file with mode: 0644]
time_dependent_navier_stokes/doc/builds-on [new file with mode: 0644]
time_dependent_navier_stokes/doc/dependencies [new file with mode: 0644]
time_dependent_navier_stokes/doc/entry-name [new file with mode: 0644]
time_dependent_navier_stokes/doc/tooltip [new file with mode: 0644]
time_dependent_navier_stokes/time_dependent_navier_stokes.cc [new file with mode: 0644]

diff --git a/time_dependent_navier_stokes/CMakeLists.txt b/time_dependent_navier_stokes/CMakeLists.txt
new file mode 100644 (file)
index 0000000..8eac2e3
--- /dev/null
@@ -0,0 +1,53 @@
+##
+#  CMake script for the time-dependent Navier-Stokes solver:
+##
+
+# Set the name of the project and target:
+SET(TARGET "fluid")
+
+# 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
+  time_dependent_navier_stokes.cc
+  )
+
+# Usually, you will not need to modify anything beyond this point...
+
+CMAKE_MINIMUM_REQUIRED(VERSION 2.8.12)
+
+FIND_PACKAGE(deal.II 9.0.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) # keep in one line
+  MESSAGE(FATAL_ERROR "
+Error! This tutorial requires a deal.II library that was configured with the following options:
+    DEAL_II_WITH_UMFPACK = ON
+However, the deal.II library found at ${DEAL_II_PATH} was configured with these options
+    DEAL_II_WITH_UMFPACK = ${DEAL_II_WITH_UMFPACK}
+which conflict with the requirements."
+    )
+ENDIF()
+
+
+DEAL_II_INITIALIZE_CACHED_VARIABLES()
+PROJECT(${TARGET})
+DEAL_II_INVOKE_AUTOPILOT()
diff --git a/time_dependent_navier_stokes/Readme.md b/time_dependent_navier_stokes/Readme.md
new file mode 100644 (file)
index 0000000..989029b
--- /dev/null
@@ -0,0 +1,138 @@
+Time-dependent Navier-Stokes
+------------------------------------------
+
+### General description of the problem ###
+
+We solve the time-dependent Navier-Stokes equations with implicit-explicit (IMEX) scheme.
+Here is the equations we want to solve:
+@f{eqnarray*}
+{\bold{u}}_{,t} - \nu {\nabla}^2\bold{u} + (\bold{u}\cdot\nabla)\bold{u} + \nabla p = \bold{f}
+@f}
+
+The idea is as follows: we use backward Euler time for time discretization. The diffusion term
+is treated implicitly and the convection term is treated explicitly. Let $(u, p)$ denote the
+velocity and pressure, respectively and $(v, q)$ denote the corresponding test functions, we
+end up with the following linear system:
+@f{eqnarray*}
+m(u^{n+1}, v) + \Delta{t}\cdot a((u^{n+1}, p^{n+1}), (v, q))=m(u^n, v)-\Delta{t}c(u^n;u^n, v)
+@f}
+
+where $a((u, p), (v, q))$ is the bilinear form of the diffusion term:
+@f{eqnarray*}
+a((u, p), (v, q)) = \int_\Omega \nu\nabla{u}\nabla{v}-p\nabla\cdot v-q\nabla\cdot ud\Omega
+@f}
+
+$m(u, v)$ is the mass matrix:
+@f{eqnarray*}
+m(u, v) = \int_{\Omega} u \cdot v d\Omega
+@f}
+
+and $c(u;u, v)$ is the convection term:
+@f{eqnarray*}
+c(u;u, v) = \int_{\Omega} (u \cdot \nabla u) \cdot v d\Omega
+@f}
+
+Substracting $m(u^n, v) + \Delta{t}a((u^n, p^n), (v, q))$ from both sides of the equation,
+we have the incremental form:
+@f{eqnarray*}
+m(\Delta{u}, v) + \Delta{t}\cdot a((\Delta{u}, \Delta{p}), (v, q)) = \Delta{t}(-a(u^n, p^n), (q, v)) - \Delta{t}c(u^n;u^n, v)
+@f}
+
+
+The system we want to solve can be written in matrix 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}
+
+#### Grad-Div stablization ####
+
+Similar to step-57, we add $\gamma B^T M_p^{-1} B$ to the upper left block of the system,
+thus the system becomes:
+
+@f{eqnarray*}
+    \left(
+      \begin{array}{cc}
+        \tilde{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}
+where $\tilde{A} = A + \gamma B^T M_p^{-1} B$.
+
+Detailed explaination of the Grad-Div stablization can be found in [1].
+
+#### Block preconditioner ####
+
+The block preconditioner is pretty much the same as in step-22, except for two additional terms,
+namely the inertial term (mass matrix) and the Grad-Div term.
+
+The block preconditioner can be written as:
+@f{eqnarray*}
+    P^{-1}
+    =
+    \left(
+      \begin{array}{cc}
+        {\tilde{A}}^{-1} & 0 \\
+        {\tilde{S}}^{-1}B{\tilde{A}}^{-1} & -{\tilde{S}}^{-1} \\
+      \end{array}
+    \right)
+@f}
+where ${\tilde{S}}$ is the Schur complement of ${\tilde{A}}$, which can be decomposed 
+into the Schur complements of the mass matrix, diffusion matrix, and the Grad-Div term:
+@f{eqnarray*}
+    {\tilde{S}}^{-1}
+    \approx 
+    {S_{mass}}^{-1} + {S_{diff}}^{-1} + {S_{Grad-Div}}^{-1}
+    \approx 
+    [B^T (diag M)^{-1} B]^{-1} + \Delta{t}(\nu + \gamma)M_p^{-1}
+@f}
+
+For more information about preconditioning incompressible Navier-Stokes equations, please refer
+to [1] and [2].
+
+#### Test case ####
+We test the code with a classical benchmark case, flow past a cylinder, with Reynold's
+number 100. The geometry setup of the case can be found on
+[this webpage](http://www.featflow.de/en/benchmarks/cfdbenchmarking/flow/dfg_benchmark2_re100.html).
+
+### Acknowledgements ###
+Thank Wolfgang Bangerth, Timo Heister and Martin Kronbichler for their helpful discussions
+on my numerical formulation and implementation.
+
+------------------------------------------
+### References ###
+[1] Timo Heister. A massively parallel finite element framework with application to incompressible flows. Doctoral dissertation, University of Gottingen, 2011.
+
+[2] M. Kronbichler, A. Diagne and H. Holmgren. A fast massively parallel two-phase flow solver for microfluidic chip simulation, International Journal of High Performance Computing Applications, 2016.
diff --git a/time_dependent_navier_stokes/doc/Re100.mp4 b/time_dependent_navier_stokes/doc/Re100.mp4
new file mode 100644 (file)
index 0000000..b050e0f
Binary files /dev/null and b/time_dependent_navier_stokes/doc/Re100.mp4 differ
diff --git a/time_dependent_navier_stokes/doc/Re100.png b/time_dependent_navier_stokes/doc/Re100.png
new file mode 100644 (file)
index 0000000..77027e3
Binary files /dev/null and b/time_dependent_navier_stokes/doc/Re100.png differ
diff --git a/time_dependent_navier_stokes/doc/author b/time_dependent_navier_stokes/doc/author
new file mode 100644 (file)
index 0000000..2995ade
--- /dev/null
@@ -0,0 +1 @@
+Jie Cheng <chengjiehust@gmail.com>
diff --git a/time_dependent_navier_stokes/doc/builds-on b/time_dependent_navier_stokes/doc/builds-on
new file mode 100644 (file)
index 0000000..9818e36
--- /dev/null
@@ -0,0 +1 @@
+step-20 step-22 step-57
diff --git a/time_dependent_navier_stokes/doc/dependencies b/time_dependent_navier_stokes/doc/dependencies
new file mode 100644 (file)
index 0000000..3f4cba3
--- /dev/null
@@ -0,0 +1 @@
+DEAL_II_WITH_CXX11
diff --git a/time_dependent_navier_stokes/doc/entry-name b/time_dependent_navier_stokes/doc/entry-name
new file mode 100644 (file)
index 0000000..03d7b99
--- /dev/null
@@ -0,0 +1 @@
+Time-dependent Navier-Stokes
diff --git a/time_dependent_navier_stokes/doc/tooltip b/time_dependent_navier_stokes/doc/tooltip
new file mode 100644 (file)
index 0000000..ca55f62
--- /dev/null
@@ -0,0 +1 @@
+Solving time-dependent incompressible Navier-Stokes problem with Grad-Div stablization using IMEX scheme.
diff --git a/time_dependent_navier_stokes/time_dependent_navier_stokes.cc b/time_dependent_navier_stokes/time_dependent_navier_stokes.cc
new file mode 100644 (file)
index 0000000..85f0976
--- /dev/null
@@ -0,0 +1,1008 @@
+#include <deal.II/base/function.h>
+#include <deal.II/base/logstream.h>
+#include <deal.II/base/quadrature_lib.h>
+#include <deal.II/base/tensor.h>
+#include <deal.II/base/utilities.h>
+#include <deal.II/base/timer.h>
+
+#include <deal.II/lac/block_sparse_matrix.h>
+#include <deal.II/lac/block_vector.h>
+#include <deal.II/lac/constraint_matrix.h>
+#include <deal.II/lac/dynamic_sparsity_pattern.h>
+#include <deal.II/lac/full_matrix.h>
+#include <deal.II/lac/precondition.h>
+#include <deal.II/lac/solver_cg.h>
+#include <deal.II/lac/solver_gmres.h>
+#include <deal.II/lac/sparse_direct.h>
+#include <deal.II/lac/sparse_ilu.h>
+
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/grid_out.h>
+#include <deal.II/grid/grid_refinement.h>
+#include <deal.II/grid/grid_tools.h>
+#include <deal.II/grid/manifold_lib.h>
+#include <deal.II/grid/tria.h>
+#include <deal.II/grid/tria_accessor.h>
+#include <deal.II/grid/tria_iterator.h>
+
+#include <deal.II/dofs/dof_accessor.h>
+#include <deal.II/dofs/dof_handler.h>
+#include <deal.II/dofs/dof_renumbering.h>
+#include <deal.II/dofs/dof_tools.h>
+
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/fe_system.h>
+#include <deal.II/fe/fe_values.h>
+
+#include <deal.II/numerics/data_out.h>
+#include <deal.II/numerics/error_estimator.h>
+#include <deal.II/numerics/matrix_tools.h>
+#include <deal.II/numerics/vector_tools.h>
+#include <deal.II/numerics/solution_transfer.h>
+
+#include <deal.II/physics/elasticity/standard_tensors.h>
+
+#include <fstream>
+#include <iostream>
+#include <sstream>
+
+namespace fluid
+{
+  using namespace dealii;
+
+  // @sect3{Create the triangulation}
+
+  // The code to create triangulation is copied from Martin Kronbichler's code
+  // (https://github.com/kronbichler/adaflo/blob/master/tests/flow_past_cylinder.cc)
+  // with very few modifications.
+  // Helper function used in both 2d and 3d:
+  void create_triangulation_2d(Triangulation<2> &tria, bool compute_in_2d = true)
+  {
+    SphericalManifold<2> boundary(Point<2>(0.5, 0.2));
+    Triangulation<2> left, middle, right, tmp, tmp2;
+    GridGenerator::subdivided_hyper_rectangle(
+      left,
+      std::vector<unsigned int>({3U, 4U}),
+      Point<2>(),
+      Point<2>(0.3, 0.41),
+      false);
+    GridGenerator::subdivided_hyper_rectangle(
+      right,
+      std::vector<unsigned int>({18U, 4U}),
+      Point<2>(0.7, 0),
+      Point<2>(2.5, 0.41),
+      false);
+
+    // Create middle part first as a hyper shell.
+    GridGenerator::hyper_shell(middle, Point<2>(0.5, 0.2), 0.05, 0.2, 4, true);
+    middle.set_manifold(0, boundary);
+    middle.refine_global(1);
+
+    // Then move the vertices to the points where we want them to be to create a
+    // slightly asymmetric cube with a hole
+    for (Triangulation<2>::cell_iterator cell = middle.begin();
+         cell != middle.end();
+         ++cell)
+      for (unsigned int v = 0; v < GeometryInfo<2>::vertices_per_cell; ++v)
+        {
+          Point<2> &vertex = cell->vertex(v);
+          if (std::abs(vertex[0] - 0.7) < 1e-10 &&
+              std::abs(vertex[1] - 0.2) < 1e-10)
+            vertex = Point<2>(0.7, 0.205);
+          else if (std::abs(vertex[0] - 0.6) < 1e-10 &&
+                   std::abs(vertex[1] - 0.3) < 1e-10)
+            vertex = Point<2>(0.7, 0.41);
+          else if (std::abs(vertex[0] - 0.6) < 1e-10 &&
+                   std::abs(vertex[1] - 0.1) < 1e-10)
+            vertex = Point<2>(0.7, 0);
+          else if (std::abs(vertex[0] - 0.5) < 1e-10 &&
+                   std::abs(vertex[1] - 0.4) < 1e-10)
+            vertex = Point<2>(0.5, 0.41);
+          else if (std::abs(vertex[0] - 0.5) < 1e-10 &&
+                   std::abs(vertex[1] - 0.0) < 1e-10)
+            vertex = Point<2>(0.5, 0.0);
+          else if (std::abs(vertex[0] - 0.4) < 1e-10 &&
+                   std::abs(vertex[1] - 0.3) < 1e-10)
+            vertex = Point<2>(0.3, 0.41);
+          else if (std::abs(vertex[0] - 0.4) < 1e-10 &&
+                   std::abs(vertex[1] - 0.1) < 1e-10)
+            vertex = Point<2>(0.3, 0);
+          else if (std::abs(vertex[0] - 0.3) < 1e-10 &&
+                   std::abs(vertex[1] - 0.2) < 1e-10)
+            vertex = Point<2>(0.3, 0.205);
+          else if (std::abs(vertex[0] - 0.56379) < 1e-4 &&
+                   std::abs(vertex[1] - 0.13621) < 1e-4)
+            vertex = Point<2>(0.59, 0.11);
+          else if (std::abs(vertex[0] - 0.56379) < 1e-4 &&
+                   std::abs(vertex[1] - 0.26379) < 1e-4)
+            vertex = Point<2>(0.59, 0.29);
+          else if (std::abs(vertex[0] - 0.43621) < 1e-4 &&
+                   std::abs(vertex[1] - 0.13621) < 1e-4)
+            vertex = Point<2>(0.41, 0.11);
+          else if (std::abs(vertex[0] - 0.43621) < 1e-4 &&
+                   std::abs(vertex[1] - 0.26379) < 1e-4)
+            vertex = Point<2>(0.41, 0.29);
+        }
+
+    // Refine once to create the same level of refinement as in the
+    // neighboring domains:
+    middle.refine_global(1);
+
+    // Must copy the triangulation because we cannot merge triangulations with
+    // refinement:
+    GridGenerator::flatten_triangulation(middle, tmp2);
+
+    // Left domain is requred in 3d only.
+    if (compute_in_2d)
+    {
+      GridGenerator::merge_triangulations(tmp2, right, tria);
+    }
+    else
+    {
+      GridGenerator::merge_triangulations(left, tmp2, tmp);
+      GridGenerator::merge_triangulations(tmp, right, tria);
+    }
+  }
+
+  // Create 2D triangulation:
+  void create_triangulation(Triangulation<2> &tria)
+  {
+    create_triangulation_2d(tria);
+    // Set the cylinder boundary to 1, the right boundary (outflow) to 2, the rest to 0.
+    for (Triangulation<2>::active_cell_iterator cell = tria.begin();
+         cell != tria.end();
+         ++cell)
+    {
+      for (unsigned int f = 0; f < GeometryInfo<2>::faces_per_cell; ++f)
+      {
+        if (cell->face(f)->at_boundary())
+          {
+            if (std::abs(cell->face(f)->center()[0] - 2.5) < 1e-12)
+            {
+              cell->face(f)->set_all_boundary_ids(2);
+            }
+            else if (Point<2>(0.5, 0.2).distance(cell->face(f)->center()) <= 0.05)
+              {
+                cell->face(f)->set_all_manifold_ids(10);
+                cell->face(f)->set_all_boundary_ids(1);
+              }
+            else
+            {
+              cell->face(f)->set_all_boundary_ids(0);
+            }
+          }
+      }
+    }
+  }
+
+  // Create 3D triangulation:
+  void create_triangulation(Triangulation<3> &tria)
+  {
+    Triangulation<2> tria_2d;
+    create_triangulation_2d(tria_2d, false);
+    GridGenerator::extrude_triangulation(tria_2d, 5, 0.41, tria);
+    // Set the cylinder boundary to 1, the right boundary (outflow) to 2, the rest to 0.
+    for (Triangulation<3>::active_cell_iterator cell = tria.begin();
+        cell != tria.end(); ++cell)
+    {
+      for (unsigned int f = 0; f<GeometryInfo<3>::faces_per_cell; ++f)
+      {
+        if (cell->face(f)->at_boundary())
+        {
+          if (std::abs(cell->face(f)->center()[0]-2.5) < 1e-12)
+          {
+            cell->face(f)->set_all_boundary_ids(2);
+          }
+          else if (Point<3>(0.5, 0.2, cell->face(f)->center()[2]).distance
+            (cell->face(f)->center()) <= 0.05)
+          {
+            cell->face(f)->set_all_manifold_ids(10);
+            cell->face(f)->set_all_boundary_ids(1);
+          }
+          else
+          {
+            cell->face(f)->set_all_boundary_ids(0);
+          }
+        }
+      }
+    }
+  }
+
+  // @sect3{Time stepping}
+  class Time
+  {
+  public:
+    Time(const double time_end, const double delta_t)
+      : timestep(0), time_current(0.0), time_end(time_end), delta_t(delta_t)
+    {
+    }
+    virtual ~Time() {}
+    double current() const { return time_current; }
+    double end() const { return time_end; }
+    double get_delta_t() const { return delta_t; }
+    unsigned int get_timestep() const { return timestep; }
+    void increment()
+    {
+      time_current += delta_t;
+      ++timestep;
+    }
+
+  private:
+    unsigned int timestep;
+    double time_current;
+    const double time_end;
+    const double delta_t;
+  };
+
+  // @sect3{Boundary values}
+
+  // Dirichlet boundary conditions for the velocity inlet and walls
+  template <int dim>
+  class BoundaryValues : public Function<dim>
+  {
+  public:
+    BoundaryValues() : Function<dim>(dim + 1) {}
+    virtual double value(const Point<dim> &p,
+                         const unsigned int component) const;
+
+    virtual void vector_value(const Point<dim> &p,
+                              Vector<double> &values) const;
+  };
+
+  template <int dim>
+  double BoundaryValues<dim>::value(const Point<dim> &p,
+                                    const unsigned int component) const
+  {
+    Assert(component < this->n_components,
+           ExcIndexRange(component, 0, this->n_components));
+    if (component == 0 && std::abs(p[0] - 0.3) < 1e-10)
+      {
+        double U = 1.5;
+        double y = p[1];
+        return 4 * U * y * (0.41 - y) / (0.41 * 0.41);
+      }
+    return 0;
+  }
+
+  template <int dim>
+  void BoundaryValues<dim>::vector_value(const Point<dim> &p,
+                                         Vector<double> &values) const
+  {
+    for (unsigned int c = 0; c < this->n_components; ++c)
+      values(c) = BoundaryValues<dim>::value(p, c);
+  }
+
+  // @sect3{Preconditioners}
+
+  // The LHS of the system matrix is the same as Stokes equation for IMEX scheme.
+  // A block preconditioner as in step-22 is used here.
+  // @sect4{Inner preconditioner}
+
+  // Adapted from step-22, used to solve for ${\tilde{A}}^{-1}$
+  template <int dim>
+  struct InnerPreconditioner;
+
+  template <>
+  struct InnerPreconditioner<2>
+  {
+    typedef SparseDirectUMFPACK type;
+  };
+
+  template <>
+  struct InnerPreconditioner<3>
+  {
+    typedef SparseILU<double> type;
+  };
+
+  // @sect4{Inverse matrix}
+
+  // This is used for ${\tilde{S}}^{-1}$ and ${\tilde{A}}^{-1}$, which are symmetric so we use CG
+  // solver inside
+  template <class MatrixType, class PreconditionerType>
+  class InverseMatrix : public Subscriptor
+  {
+  public:
+    InverseMatrix(const MatrixType &m,
+                  const PreconditionerType &preconditioner);
+    void vmult(Vector<double> &dst, const Vector<double> &src) const;
+
+  private:
+    const SmartPointer<const MatrixType> matrix;
+    const SmartPointer<const PreconditionerType> preconditioner;
+  };
+
+  template <class MatrixType, class PreconditionerType>
+  InverseMatrix<MatrixType, PreconditionerType>::InverseMatrix(
+    const MatrixType &m, const PreconditionerType &preconditioner)
+    : matrix(&m), preconditioner(&preconditioner)
+  {
+  }
+
+  template <class MatrixType, class PreconditionerType>
+  void InverseMatrix<MatrixType, PreconditionerType>::vmult(
+    Vector<double> &dst, const Vector<double> &src) const
+  {
+    SolverControl solver_control(src.size(), 1e-6 * src.l2_norm());
+    SolverCG<> cg(solver_control);
+    dst = 0;
+    cg.solve(*matrix, dst, src, *preconditioner);
+  }
+
+  // @sect4{Approximate Schur complement of mass matrix}
+
+  // The Schur complement of mass matrix is written as $S_M = BM^{-1}B^T$
+  // Similar to step-20, we use $B(diag(M))^{-1}B^T$ to approximate it.
+  class ApproximateMassSchur : public Subscriptor
+  {
+  public:
+    ApproximateMassSchur(const BlockSparseMatrix<double> &M);
+    void vmult(Vector<double> &dst, const Vector<double> &src) const;
+
+  private:
+    const SmartPointer<const BlockSparseMatrix<double>> mass_matrix;
+    mutable Vector<double> tmp1, tmp2;
+  };
+
+  ApproximateMassSchur::ApproximateMassSchur(
+    const BlockSparseMatrix<double> &M)
+    : mass_matrix(&M), tmp1(M.block(0, 0).m()), tmp2(M.block(0, 0).m())
+  {
+  }
+
+  void ApproximateMassSchur::vmult(Vector<double> &dst,
+                                         const Vector<double> &src) const
+  {
+    mass_matrix->block(0, 1).vmult(tmp1, src);
+    mass_matrix->block(0, 0).precondition_Jacobi(tmp2, tmp1);
+    mass_matrix->block(1, 0).vmult(dst, tmp2);
+  }
+
+  // @sect4{The inverse matrix of the system Schur complement}
+
+  // The inverse of the total Schur complement is the sum of the inverse of
+  // diffusion, Grad-Div term, and mass Schur complements. Note that the first
+  // two components add up to $\Delta{t}(\nu + \gamma)M_p^{-1}$ as introduced in step-57,
+  // in which the additional $\Delta{t}$ comes from the time discretization,
+  // and the last component is obtained by wrapping a <code>InverseMatrix<\code>
+  // around <code>ApproximateMassSchur<\code>.
+  template <class PreconditionerSm, class PreconditionerMp>
+  class SchurComplementInverse : public Subscriptor
+  {
+  public:
+    SchurComplementInverse(
+      double gamma, double viscosity, double dt,
+      const InverseMatrix<ApproximateMassSchur, PreconditionerSm> &Sm_inv,
+      const InverseMatrix<SparseMatrix<double>, PreconditionerMp> &Mp_inv);
+    void vmult(Vector<double> &dst, const Vector<double> &src) const;
+  private:
+    const double gamma;
+    const double viscosity;
+    const double dt;
+    const SmartPointer<const InverseMatrix<ApproximateMassSchur,
+      PreconditionerSm>> Sm_inverse;
+    const SmartPointer<const InverseMatrix<SparseMatrix<double>,
+      PreconditionerMp>> Mp_inverse;
+  };
+
+  template <class PreconditionerSm, class PreconditionerMp>
+  SchurComplementInverse<PreconditionerSm, PreconditionerMp>::SchurComplementInverse(
+    double gamma, double viscosity, double dt,
+    const InverseMatrix<ApproximateMassSchur, PreconditionerSm> &Sm_inv,
+    const InverseMatrix<SparseMatrix<double>, PreconditionerMp> &Mp_inv) :
+    gamma(gamma), viscosity(viscosity), dt(dt), Sm_inverse(&Sm_inv), Mp_inverse(&Mp_inv)
+  {
+  }
+
+  template <class PreconditionerSm, class PreconditionerMp>
+  void SchurComplementInverse<PreconditionerSm, PreconditionerMp>::vmult(
+    Vector<double> &dst, const Vector<double> &src) const
+  {
+    Vector<double> tmp(src.size());
+    Sm_inverse->vmult(dst, src);
+    Mp_inverse->vmult(tmp, src);
+    tmp *= (viscosity + gamma) * dt;
+    dst += tmp;
+  }
+
+  // @sect4{The block Schur preconditioner}
+
+  // The block Schur preconditioner has the same form as in step-22, which is written as
+  // $P^{-1} = [\tilde{A}}^{-1}, 0; {\tilde{S}}^{-1}B{\tilde{A}}^{-1}, -{\tilde{S}}^{-1}]$
+  // Note that ${\tilde{A}}^{-1}$ has contributions from the diffusion, Grad-Div and mass terms.
+  // This class has three template arguments: PreconditionerA is needed for ${\tilde{A}}^{-1}$,
+  // PreconditionerSm and PreconditionerMp are used in the inverse of the Schur complement
+  // of $\tilde{A}$, namely ${\tilde{S}}^{-1}$.
+  template <class PreconditionerA, class PreconditionerSm, class PreconditionerMp>
+  class BlockSchurPreconditioner : public Subscriptor
+  {
+  public:
+    BlockSchurPreconditioner(
+      const BlockSparseMatrix<double> &system_m,
+      const InverseMatrix<SparseMatrix<double>, PreconditionerA> &A_inv,
+      const SchurComplementInverse<PreconditionerSm, PreconditionerMp> &S_inv);
+    void vmult(BlockVector<double> &dst, const BlockVector<double> &src) const;
+
+  private:
+    const SmartPointer<const BlockSparseMatrix<double>> system_matrix;
+    const SmartPointer<
+      const InverseMatrix<SparseMatrix<double>, PreconditionerA>> A_inverse;
+    const SmartPointer<
+      const SchurComplementInverse<PreconditionerSm, PreconditionerMp>> S_inverse;
+    mutable Vector<double> tmp;
+  };
+
+  template <class PreconditionerA, class PreconditionerSm, class PreconditionerMp>
+  BlockSchurPreconditioner<PreconditionerA, PreconditionerSm, PreconditionerMp>::
+    BlockSchurPreconditioner(
+      const BlockSparseMatrix<double> &system_m,
+      const InverseMatrix<SparseMatrix<double>, PreconditionerA> &A_inv,
+      const SchurComplementInverse<PreconditionerSm, PreconditionerMp> &S_inv)
+    : system_matrix(&system_m), A_inverse(&A_inv), S_inverse(&S_inv),
+      tmp(system_matrix->block(1, 1).m())
+  {
+  }
+
+  template <class PreconditionerA, class PreconditionerSm, class PreconditionerMp>
+  void BlockSchurPreconditioner<PreconditionerA, PreconditionerSm, PreconditionerMp>::vmult(
+    BlockVector<double> &dst, const BlockVector<double> &src) const
+  {
+    A_inverse->vmult(dst.block(0), src.block(0));
+    system_matrix->block(1, 0).residual(tmp, dst.block(0), src.block(1));
+    tmp *= -1;
+    S_inverse->vmult(dst.block(1), tmp);
+  }
+
+  // @sect3{The time-dependent Navier-Stokes class template}
+  template <int dim>
+  class NavierStokes
+  {
+  public:
+    NavierStokes(const unsigned int degree);
+    void run();
+
+  private:
+    void setup();
+    void assemble(bool assemble_lhs);
+
+    std::pair<unsigned int, double> solve_linear_system(bool update_preconditioner);
+    void output_results(const unsigned int index) const;
+    void process_solution(std::ofstream& out) const;
+    const ConstraintMatrix &get_constraints() const;
+
+    double viscosity;
+    double gamma;
+    const unsigned int degree;
+    std::vector<types::global_dof_index> dofs_per_block;
+
+    Triangulation<dim> triangulation;
+    FESystem<dim> fe;
+    DoFHandler<dim> dof_handler;
+    QGauss<dim> quadrature_formula;
+    QGauss<dim-1> face_quadrature_formula;
+
+    ConstraintMatrix zero_constraints;
+    ConstraintMatrix nonzero_constraints;
+
+    BlockSparsityPattern sparsity_pattern;
+    BlockSparseMatrix<double> system_matrix;
+    // We need both velocity mass and pressure mass, so we use a block sparse matrix to store it.
+    BlockSparseMatrix<double> mass_matrix;
+
+    BlockVector<double> solution;
+    BlockVector<double> solution_increment;
+    BlockVector<double> system_rhs;
+
+    Time time;
+    mutable TimerOutput timer;
+
+    // We use shared pointers for all the preconditioning-related stuff
+    std::shared_ptr<ApproximateMassSchur> approximate_Sm;
+    std::shared_ptr<PreconditionIdentity> preconditioner_Sm;
+    std::shared_ptr<InverseMatrix<ApproximateMassSchur, PreconditionIdentity>> Sm_inverse;
+
+    std::shared_ptr<SparseILU<double>> preconditioner_Mp;
+    std::shared_ptr<InverseMatrix<SparseMatrix<double>, SparseILU<double>>> Mp_inverse;
+
+    std::shared_ptr<SchurComplementInverse<PreconditionIdentity,
+      SparseILU<double>>> S_inverse;
+
+    std::shared_ptr<typename InnerPreconditioner<dim>::type> preconditioner_A;
+    std::shared_ptr<InverseMatrix<SparseMatrix<double>, 
+      typename InnerPreconditioner<dim>::type>> A_inverse;
+
+    std::shared_ptr<BlockSchurPreconditioner
+      <typename InnerPreconditioner<dim>::type, PreconditionIdentity, SparseILU<double>>> preconditioner;
+  };
+
+  // @sect4{NavierStokes::NavierStokes}
+  template <int dim>
+  NavierStokes<dim>::NavierStokes(const unsigned int degree)
+    : viscosity(0.001),
+      gamma(1),
+      degree(degree),
+      triangulation(Triangulation<dim>::maximum_smoothing),
+      fe(FE_Q<dim>(degree + 1), dim, FE_Q<dim>(degree), 1),
+      dof_handler(triangulation),
+      quadrature_formula(degree+2),
+      face_quadrature_formula(degree+2),
+      time(1e-2, 1e-3),
+      timer(std::cout, TimerOutput::summary, TimerOutput::wall_times)
+  {
+  }
+
+  // @sect4{NavierStokes::setup}
+  template <int dim>
+  void NavierStokes<dim>::setup()
+  {
+    timer.enter_subsection("Setup system");
+    dof_handler.distribute_dofs(fe);
+    DoFRenumbering::Cuthill_McKee(dof_handler);
+
+    // 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
+    std::vector<unsigned int> 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];
+
+    // The Dirichlet boundary condition is applied to boundaries 0 and 1.
+    FEValuesExtractors::Vector velocities(0);
+    {
+      nonzero_constraints.clear();
+
+      DoFTools::make_hanging_node_constraints(dof_handler, nonzero_constraints);
+      VectorTools::interpolate_boundary_values(dof_handler,
+                                               0,
+                                               BoundaryValues<dim>(),
+                                               nonzero_constraints,
+                                               fe.component_mask(velocities));
+      VectorTools::interpolate_boundary_values(dof_handler,
+                                               1,
+                                               BoundaryValues<dim>(),
+                                               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,
+        Functions::ZeroFunction<dim>(dim + 1),
+        zero_constraints,
+        fe.component_mask(velocities));
+      VectorTools::interpolate_boundary_values(
+        dof_handler,
+        1,
+        Functions::ZeroFunction<dim>(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 vertices: " << triangulation.n_vertices()
+              << std::endl
+              << "   Number of degrees of freedom: " << dof_handler.n_dofs()
+              << " (" << dof_u << '+' << dof_p << ')' << std::endl;
+
+    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);
+    mass_matrix.reinit(sparsity_pattern);
+
+    solution.reinit(dofs_per_block);
+    solution_increment.reinit(dofs_per_block);
+    system_rhs.reinit(dofs_per_block);
+
+    timer.leave_subsection();
+  }
+
+  // @sect4{NavierStokes::setup}
+
+  // A helper function to determine which constrint to use based on the current timestep
+  template <int dim>
+  const ConstraintMatrix &NavierStokes<dim>::get_constraints() const
+  {
+    return time.get_timestep() == 0 ? nonzero_constraints : zero_constraints;
+  }
+
+  // @sect4{NavierStokes::assemble}
+
+  // Note that we only need to assemble the LHS for twice: once with the nonzero constraint
+  // and once for zero constraint. But we must assemble the RHS at every time step.
+  template <int dim>
+  void NavierStokes<dim>::assemble(bool assemble_lhs)
+  {
+    timer.enter_subsection("Assemble system");
+    if (assemble_lhs)
+      {
+        system_matrix = 0;
+        mass_matrix = 0;
+      }
+
+    system_rhs = 0;
+
+    FEValues<dim> 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<double> local_matrix(dofs_per_cell, dofs_per_cell);
+    FullMatrix<double> local_mass_matrix(dofs_per_cell, dofs_per_cell);
+    Vector<double> local_rhs(dofs_per_cell);
+
+    std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
+
+    std::vector<Tensor<1, dim>> current_velocity_values(n_q_points);
+    std::vector<Tensor<2, dim>> current_velocity_gradients(n_q_points);
+    std::vector<double> current_velocity_divergences(n_q_points);
+    std::vector<double> current_pressure_values(n_q_points);
+
+    std::vector<double> div_phi_u(dofs_per_cell);
+    std::vector<Tensor<1, dim>> phi_u(dofs_per_cell);
+    std::vector<Tensor<2, dim>> grad_phi_u(dofs_per_cell);
+    std::vector<double> phi_p(dofs_per_cell);
+
+    typename DoFHandler<dim>::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;
+        local_mass_matrix = 0;
+
+        fe_values[velocities].get_function_values(solution,
+                                                  current_velocity_values);
+
+        fe_values[velocities].get_function_gradients(
+          solution, current_velocity_gradients);
+
+        fe_values[velocities].get_function_divergences(
+          solution, current_velocity_divergences);
+
+        fe_values[pressure].get_function_values(solution,
+                                                current_pressure_values);
+
+        for (unsigned int q = 0; q < n_q_points; ++q)
+          {
+            for (unsigned int k = 0; k < dofs_per_cell; ++k)
+              {
+                div_phi_u[k] = fe_values[velocities].divergence(k, q);
+                grad_phi_u[k] = fe_values[velocities].gradient(k, q);
+                phi_u[k] = fe_values[velocities].value(k, q);
+                phi_p[k] = fe_values[pressure].value(k, q);
+              }
+
+            for (unsigned int i = 0; i < dofs_per_cell; ++i)
+              {
+                if (assemble_lhs)
+                  {
+                    for (unsigned int j = 0; j < dofs_per_cell; ++j)
+                      {
+                        // $LHS = a((u, p), (v, q))*dt + m(u, v)
+                        //      = ((grad_v, nu*grad_u) - (div_v, p) - (q, div_u))*dt +
+                        //        m(u, v)$ plus Grad-Div term.
+                        local_matrix(i, j) +=
+                          ((viscosity *
+                             scalar_product(grad_phi_u[j], grad_phi_u[i]) -
+                           div_phi_u[i] * phi_p[j] - phi_p[i] * div_phi_u[j] +
+                           gamma*div_phi_u[j]*div_phi_u[i]) *
+                            time.get_delta_t() +
+                           phi_u[i] * phi_u[j]) *
+                          fe_values.JxW(q);
+                        // Besides the velocity and pressure mass matrices, we also
+                        // assemble $B^T$ and $B$ into the block mass matrix for convenience
+                        // because we need to use them to compute the Schur complement.
+                        // As a result $M = [M_u, B^T; B, M_p]$.
+                        local_mass_matrix(i, j) +=
+                          (phi_u[i] * phi_u[j] + phi_p[i] * phi_p[j] -
+                           div_phi_u[i] * phi_p[j] - phi_p[i] * div_phi_u[j]) *
+                          fe_values.JxW(q);
+                      }
+                  }
+                // $RHS = - dt*[ a((u_prev, p_prev), (v, q)) + c(u_prev; u_prev, v)]$
+                // plus Grad-Div term.
+                local_rhs(i) -=
+                  (viscosity * scalar_product(current_velocity_gradients[q],
+                                              grad_phi_u[i]) -
+                    current_velocity_divergences[q] * phi_p[i] -
+                    current_pressure_values[q] * div_phi_u[i] +
+                    current_velocity_gradients[q] * current_velocity_values[q] *
+                      phi_u[i] +
+                    gamma * current_velocity_divergences[q] * div_phi_u[i]) *
+                  fe_values.JxW(q) * time.get_delta_t();
+              }
+          }
+
+        cell->get_dof_indices(local_dof_indices);
+
+        const ConstraintMatrix &constraints_used = get_constraints();
+
+        if (assemble_lhs)
+          {
+            constraints_used.distribute_local_to_global(local_matrix,
+                                                        local_rhs,
+                                                        local_dof_indices,
+                                                        system_matrix,
+                                                        system_rhs);
+            constraints_used.distribute_local_to_global(local_mass_matrix,
+                                                        local_dof_indices,
+                                                        mass_matrix);
+          }
+        else
+          {
+            constraints_used.distribute_local_to_global(
+              local_rhs, local_dof_indices, system_rhs);
+          }
+      }
+    timer.leave_subsection();
+  }
+
+  // @sect4{NavierStokes::solve_linear_system}
+
+  // Only updates the preconditioners when we assemble the LHS of the system.
+  template <int dim>
+  std::pair<unsigned int, double> NavierStokes<dim>::solve_linear_system(bool update_precondition)
+  {
+    const ConstraintMatrix &constraints_used = get_constraints();
+
+    if (update_precondition)
+    {
+      timer.enter_subsection("Precondition linear system");
+
+      preconditioner.reset();
+      A_inverse.reset();
+      preconditioner_A.reset();
+      S_inverse.reset();
+      Mp_inverse.reset();
+      preconditioner_Mp.reset();
+      Sm_inverse.reset();
+      preconditioner_Sm.reset();
+      approximate_Sm.reset();
+
+      approximate_Sm.reset(new ApproximateMassSchur(mass_matrix));
+      preconditioner_Sm.reset(new PreconditionIdentity());
+      Sm_inverse.reset(new InverseMatrix<ApproximateMassSchur, PreconditionIdentity>
+        (*approximate_Sm, *preconditioner_Sm));
+      preconditioner_Mp.reset(new SparseILU<double>());
+      preconditioner_Mp->initialize(mass_matrix.block(1,1));
+      Mp_inverse.reset(new InverseMatrix<SparseMatrix<double>, SparseILU<double>>
+        (mass_matrix.block(1,1), *preconditioner_Mp)); 
+      S_inverse.reset(new SchurComplementInverse<PreconditionIdentity,
+        SparseILU<double>>(gamma, viscosity, time.get_delta_t(), *Sm_inverse, *Mp_inverse));
+      preconditioner_A.reset(new typename InnerPreconditioner<dim>::type());
+      preconditioner_A->initialize(system_matrix.block(0,0),
+        typename InnerPreconditioner<dim>::type::AdditionalData());
+
+      A_inverse.reset(new InverseMatrix<SparseMatrix<double>,
+        typename InnerPreconditioner<dim>::type>(system_matrix.block(0,0), *preconditioner_A)); 
+      preconditioner.reset(new BlockSchurPreconditioner<
+        typename InnerPreconditioner<dim>::type, PreconditionIdentity,
+        SparseILU<double>>(system_matrix, *A_inverse, *S_inverse));
+
+      timer.leave_subsection();
+    }
+
+    // Solve with GMRES solver.
+    timer.enter_subsection("Solve linear system");
+    SolverControl solver_control(system_matrix.m(),
+                                 1e-8 * system_rhs.l2_norm());
+    GrowingVectorMemory<BlockVector<double>> vector_memory;
+    SolverGMRES<BlockVector<double>>::AdditionalData gmres_data;
+    gmres_data.max_n_tmp_vectors = 100;
+    SolverGMRES<BlockVector<double>> gmres(
+      solver_control, vector_memory, gmres_data);
+    gmres.solve(system_matrix, solution_increment, system_rhs, *preconditioner);
+
+    constraints_used.distribute(solution_increment);
+    timer.leave_subsection();
+
+    return {solver_control.last_step(), solver_control.last_value()};
+  }
+
+  // @sect4{NavierStokes::run}
+
+  template <int dim>
+  void NavierStokes<dim>::run()
+  {
+    create_triangulation(triangulation);
+    triangulation.refine_global(2);
+    setup();
+
+    std::ofstream out("grid.eps");
+    GridOut grid_out;
+    grid_out.write_eps(triangulation, out);
+
+    std::ofstream out2("force.txt");
+    out2 << std::setw(13) << std::left << "Time/s"
+      << std::setw(13) << std::left << " Drag" << std::setw(13)
+      << std::left << " Lift" << std::endl;
+
+    // In IMEX scheme we do not need to implement the Newton's method, what we need
+    // to do at every time step is simple:
+    // 1. Solve for the solution increment; 2. Update the solution.
+    output_results(time.get_timestep());
+    while (time.current() <= time.end())
+      {
+        std::cout << "*****************************************" << std::endl;
+        std::cout << "Time = " << time.current() << std::endl;
+
+        assemble(time.get_timestep() < 2);
+
+        auto state = solve_linear_system(time.get_timestep() < 2);
+        solution.add(1.0, solution_increment);
+
+        // solution is distributed using nonzero_constraints all the time
+        nonzero_constraints.distribute(solution);
+        solution_increment = 0;
+
+        std::cout << " FGMRES steps = " << state.first 
+          << " residual = " << std::setw(6) << state.second << std::endl;
+
+        time.increment();
+
+        if (time.get_timestep() % 1 == 0)
+        {
+          output_results(time.get_timestep());
+          process_solution(out2);
+        }
+      }
+
+    out2.close();
+  }
+
+  // @sect4{NavierStokes::output_result}
+
+  template <int dim>
+  void NavierStokes<dim>::output_results(const unsigned int output_index) const
+  {
+    timer.enter_subsection("Output");
+    std::cout << " Writing results..." << std::endl;
+    std::vector<std::string> solution_names(dim, "velocity");
+    solution_names.push_back("pressure");
+
+    std::vector<DataComponentInterpretation::DataComponentInterpretation>
+      data_component_interpretation(
+        dim, DataComponentInterpretation::component_is_part_of_vector);
+    data_component_interpretation.push_back(
+      DataComponentInterpretation::component_is_scalar);
+    DataOut<dim> data_out;
+    data_out.attach_dof_handler(dof_handler);
+    data_out.add_data_vector(solution,
+                             solution_names,
+                             DataOut<dim>::type_dof_data,
+                             data_component_interpretation);
+    data_out.build_patches();
+
+    std::ostringstream filename;
+    filename << "Re100-"
+             << Utilities::int_to_string(output_index, 6) << ".vtu";
+
+    std::ofstream output(filename.str().c_str());
+    data_out.write_vtu(output);
+    timer.leave_subsection();
+  }
+
+  // @sect4{NavierStokes::process_solution}
+
+  // This function is used to calculate the drag and lift coefficients on the cylinder.
+  // We first calculate the traction of the fluid, which is nothing but the product of the
+  // stress tensor and the normal of the cylindrical surface, and then integrate it along
+  // the cylindrical surface and negate it.
+  template <int dim>
+  void NavierStokes<dim>::process_solution(std::ofstream& out) const
+  {
+    timer.enter_subsection("Process solution");
+    
+    Tensor<1, dim> force;
+
+    FEFaceValues<dim> fe_face_values(fe,
+                            face_quadrature_formula,
+                            update_values | update_quadrature_points |
+                            update_JxW_values | update_normal_vectors |
+                            update_gradients);
+
+    const unsigned int n_q_points = face_quadrature_formula.size();
+
+    const FEValuesExtractors::Vector velocities(0);
+    const FEValuesExtractors::Scalar pressure(dim);
+
+    std::vector<double> p(n_q_points);
+    std::vector<SymmetricTensor<2, dim>> grad_sym_v(n_q_points);
+
+    for (auto cell = dof_handler.begin_active(); cell != dof_handler.end(); ++cell)
+    {
+      for (unsigned int f = 0; f < GeometryInfo<2>::faces_per_cell; ++f)
+      {
+        if (cell->face(f)->at_boundary() && cell->face(f)->boundary_id() == 1)
+        {
+          fe_face_values.reinit(cell, f);
+          fe_face_values[pressure].get_function_values(solution, p);
+          fe_face_values[velocities].get_function_symmetric_gradients(solution, grad_sym_v);
+          for (unsigned int q = 0; q < n_q_points; ++q)
+          {
+            const Tensor<1, dim> &N = fe_face_values.normal_vector(q);
+            SymmetricTensor<2, dim> stress = -p[q]*Physics::Elasticity::StandardTensors<dim>::I
+              + viscosity*grad_sym_v[q];
+            force -= stress*N*fe_face_values.JxW(q);
+          }
+        }
+      }
+    }
+
+    double drag_coef = 2*force[0]/(0.1);
+    double lift_coef = 2*force[dim-1]/(0.1);
+
+    out.precision(6);
+    out.width(12);
+     
+    out << std::scientific << std::left << 
+      time.current() << " " << drag_coef << " " << lift_coef << std::endl;
+
+    timer.leave_subsection();
+  }
+}
+
+// @sect3{main function}
+
+int main()
+{
+  try
+    {
+      using namespace dealii;
+      using namespace fluid;
+
+      NavierStokes<2> flow(/* degree = */ 1);
+      flow.run();
+    }
+  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;
+}

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.