-<i> This program grew out of a student project by Abner Salgado at Texas A&M University. Most of the work for this program is by him.
+<i>
+This program grew out of a student project by Abner Salgado at Texas A&M University. Most of the work for
+this program is by him.
</i>
<a name="Intro"></a>
<h1> Introduction </h1>
-The purpose of this program is to show how to effectively solve the incompressible time-dependent Navier-Stokes equations.
-In previous tutorial programs (see for instance ****) we have seen how to solve the time-independent Stokes equations using a Schur-complement approach. However
\ No newline at end of file
+<a name="Motivation"></a>
+<h2> Motivation </h2>
+The purpose of this program is to show how to effectively solve the incompressible time-dependent
+Navier-Stokes equations. These equations describe the flow of a viscous incompressible fluid and they read
+@f{align*}
+ u_t + u \cdot \nabla u - \nu \Delta u + \nabla p = f, \\
+ \nabla \cdot u = 0,
+@f}
+supplemented by the boundary condition $u=0$ Here $u$ represents the velocity of the flow and $p$ the pressure.
+
+In previous tutorial programs (see for instance step-20 and step-22) we have seen
+how to solve the time-independent Stokes equations using a Schur-complement approach. For the
+time-dependent case, after time discretization, we would arrive at a system like
+@f{align*}
+ \frac1\tau u^k - \nu \Delta u^k + \nabla p^k = F^k, \\
+ \nabla \cdot u^k = 0,
+@f}
+where $\tau$ is the time-step. Although the structure of this system is similar to the Stokes system and thus
+it could be solved using a Schur-complement approach, it turns out that the condition number of the
+Schur-complement is proportional to $\tau^{-2}$, which makes the system very difficult to solve.
+
+<a name="Projection"></a>
+<h2> Projection Methods </h2>
+
+As we explained above, we need to come up with a different approach to solve the time-dependent Navier-Stokes
+equations. The difficulty in their solution comes from the fact that the velocity and the pressure are coupled
+through the constraint
+@f[
+ \nabla \cdot u = 0.
+@f]
+Projection methods aim at decoupling this constraint from the diffusion (Laplace) operator.
+
+Let us shortly describe how the projection methods look like in a semi-discrete setting. The objective is to
+obtain a sequence of velocities $\{u^k\}$ and pressures $\{p^k\}$. We will also obtain a sequence $\{\phi^k\}.$
+Suppose that from the initial conditions, and an application of a first order method we have found
+$(u^0,p^0,\phi^0=0)$ and $(u^1,p^1,\phi^1=p^1-p^0)$. Then the projection method consists of three steps:
+
+<ul>
+ <li> <b>Step 0</b>: Extrapolation. Define:
+ @f[
+ u^\star = 2u^k - u^{k-1}, \quad p^\sharp = p^k + \frac43 \phi^k - \frac13 \phi^{k-1}.
+ @f]
+ <li> <b>Step 1</b>: Diffusion step. We find $u^{k+1}$ that solves
+ @f[
+ \frac1{2\tau}\left( 3u^{k+1} - 4u^k + u^{k-1} \right)
+ + u^\star \cdot\nabla u^{k+1} + \frac12 \left( \nabla \cdot u^\star \right) u^{k+1}
+ -\nu \Delta u^{k+1} + \nabla p^\sharp
+ = f^{k+1},
+ \quad
+ u|_{\partial\Omega} = 0.
+ @f]
+ <li> <b>Step 2</b>: Projection. Find $\phi^{k+1}$ that solves
+ @f[
+ \Delta \phi^{k+1} = \frac3{2\tau} \nabla \cdot u^{k+1},
+ \quad
+ \partial_n \phi^{k+1} = 0.
+ @f]
+ <li> <b>Step 3</b>: Pressure correction. Here we have two options:
+ <ul>
+ <li> <i>Incremental Method in Standard Form</i>. The pressure is updated by:
+ @f[
+ p^{k+1} = p^k + \phi^{k+1}.
+ @f]
+ <li> <i>Incremental Method in Rotational Form</i>. In this case
+ @f[
+ p^{k+1} = p^k + \phi^{k+1} - \nu \nabla \cdot u^{k+1}.
+ @f]
+ </ul>
+</ul>
+
+Without going into details, which can be found in the corresponding literature, let us remark a few things about the
+projection methods that we have just described:
+<ul>
+ <li> The advection term $u\cdot\nabla u$ is replaced by its <i>skew symmetric form</i>
+ @f[
+ u \cdot \nabla u + \frac12 \left( \nabla\cdot u \right) u.
+ @f]
+ This is consistent with the equation and it is needed to guarantee unconditional stability of the
+ time-stepping scheme. Moreover, to linearize the term we use the second order extrapolation $u^\star$ of
+ $u^{k+1}$.
+ <li> The projection step is a realization of the Helmholtz decomposition
+ @f[
+ L^2(\Omega)^d = H \oplus \nabla H^1(\Omega),
+ @f]
+ where
+ @f[
+ H = \left\{ v \in L^2(\Omega)^d:\ \nabla\cdot v =0 \ v\cdot n|_{\partial\Omega} = 0 \right\}.
+ @f]
+ Indeed, if we use this decomposition on $u^{k+1}$ we obtain
+ @f[
+ u^{k+1} = v^{k+1} + \nabla \left( \frac{2\tau}{3} \phi^{k+1} \right),
+ @f]
+ with $v^{k+1}\in H$. Taking the divergence of this equation we arrive the projection equation.
+ <li> The most accurate form of the methods is the rotational one. We included the standard because it was
+ so easy to implement that it is worth having. Moreover, it is the one (in my experience) that should be
+ used if, for instance, the viscosity $\nu$ is variable.
+</ul>
+
+<a name ="testcase"></a>
+<h2> The Test Case </h2>
+
+The test case that we are going to carry out consist of the flow around a square obstacle. The geometry is
+the following
+
+@image html geometry.gif
+
+with $H=4.1$.
+
+We impose no-slip boundary conditions on both the top and bottom walls and the obstacle. On the left side we
+have the inflow boundary condition
+@f[
+ U = 4 U_m y (H-y)/H^2, \quad V = 0,
+@f]
+with $U_m = 1.5$, which corresponds to the Poiseuille flow for this configuration. Finally, on the right
+vertical wall we impose the condition that the vertical component of the velocity should be zero.
+The Reynolds number is $Re=100$ and the final time $T=10$.
+
+
+<a name="results"></a>
<h1>Results</h1>
+
+We run the code with the following <code>parameter-file.prm</code>, which can be found in the
+same directory as the source:
+@verbatim
+ # First a global definition
+ # the type of method we want to use
+ set Method_Form = rotational
+
+ subsection Physical data
+ # In this subsection we declare the physical data
+ # The initial and final time, and the Reynolds number
+ set initial_time = 0.
+ set final_time = 10.
+ set Reynolds = 100
+ end
+
+ subsection Time step data
+ # In this subsection we declare the data that is to be used for time discretization,
+ # i.e. the time step dt
+ set dt = 5e-3
+ end
+
+ subsection Space discretization
+ # In this subsection we declare the data that is relevant to the space discretization
+ # we set the number of global refines the triangulation must have
+ # and the degree k of the pair Q_(k+1)--Q_k of velocity--pressure finite element spaces
+ set n_of_refines = 3
+ set pressure_fe_degree = 1
+ end
+
+ subsection Data solve velocity
+ # In this section we declare the parameters that are going to control the solution process
+ # for the velocity.
+ set max_iterations = 1000 # maximal number of iterations that GMRES must make
+ set eps = 1e-6 # stopping criterion
+ set Krylov_size = 30 # size of the Krylov subspace to be used in GMRES
+ set off_diagonals = 60 # number of off diagonals that ILU must compute
+ set diag_strength = 0.01 # diagonal strengthening value
+ set update_prec = 10 # this number indicates how often the preconditioner must be updated
+ end
+
+ #The output frequency
+ set output = 5
+
+ #Finally we set the verbosity level
+ set verbose = false
+@endverbatim
+
+Since the <code>verbose</code> parameter is set to <code>false</code>, we do not get any kind of output.
+If we we were to set it to <code>true</code> we would get information on what the program is doing and
+how many steps each iterative process had to make to converge, etc.
+
+A plot of the obtained results for $t=0,1,5,10$ is the following
+
+@image html visit0001.jpeg
+
+@image html visit0002.jpeg
+
+@image html visit0003.jpeg
+
+@image html visit0004.jpeg
+
+The contour lines correspond to the pressure, the color to the voriticity and the arrows to the velocity.
+Here we only show a zoom into a region near the obstacle which is, after all, where the interesting things happen.
+
+<a name="extensions"></a>
+<h2> Possible Extensions </h2>
+
+This program can be extended in the following directions:
+<ul>
+ <li> Adaptive mesh refinement. As we have seen, we computed everything on a single fixed mesh.
+ Using adaptive mesh refinement can lead to increased accuracy while not significantly increasing the
+ computational time.
+ <li> Adaptive time-stepping. Although there is no theory about projection methods with variable time-step,
+ practice shows that they perform very well. This could be implemented too.
+ <li> Variable density incompressible flows. There are projection-like methods for the case of incompressible
+ flows with variable density.
+</ul>