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.
+where $u$ represents the velocity of the flow and $p$ the pressure. This system of equations is supplemented by
+the initial condition
+@f[
+ u |_{t=0} = u_0,
+@f]
+with $u_0$ sufficiently smooth and solenoidal, and suitable boundary conditions. For instance, an admissible boundary
+condition, is
+@f[
+ u|_{\partial\Omega} = u_b.
+@f]
+It is possible to prescribe other boundary conditions as well. In the test case that we solve here the boundary
+is partitioned into two disjoint subsets $\partial\Omega = \Gamma_1 \cup \Gamma_2$ and we have
+@f[
+ u|_{\Gamma_1} = u_b, \qquad u\times n|_{\Gamma_2} = 0.
+@f]
+where $n$ is the outer unit normal. The boundary condition on $\Gamma_2$ is used to model outflow conditions.
In previous tutorial programs (see for instance @ref step_20 "step-20" and
@ref step_22 "step-22") we have seen
@f]
where
@f[
- H = \left\{ v \in L^2(\Omega)^d:\ \nabla\cdot v =0 \ v\cdot n|_{\partial\Omega} = 0 \right\}.
+ 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[
@f]
with $v^{k+1}\in H$. Taking the divergence of this equation we arrive at the projection equation.
<li> The more accurate of the two variants outlined above is the rotational
- one. That said, in the program below we use the standard form because it is
- much easier to implement. Moreover, in the author's experience, it is the one that should be
- used if, for instance, the viscosity $\nu$ is variable.
+ one. However, the program below implements both variants. Moreover, in the author's experience,
+ the standard form is the one that should be used if, for instance, the viscosity $\nu$ is variable.
</ul>
+<a name = "implementation"></a>
+<h3> Implementation </h3>
+
+Our implementation of the projection methods follows <i>verbatim</i> the description given above. We must note,
+however, that as opposed to most other problems that have several solution components, we do not use
+vector-valued finite elements. Instead, we use separate finite elements for each component of the velocity
+and the pressure, respectively, and use different <code>DoFHandler</code>'s for those as well. The main
+reason for doing this is that, as we see from the description of the scheme, the <code>dim</code> components
+of the velocity and the pressure are decoupled. As a consequence, the equations for all the velocity components
+look all the same, have the same system matrix, and can be solved in parallel. Obviously, this approach
+has also its disadvantages. For instance, we need to keep several <code>DoFHandler</code>s and iterators
+synchronized when assembling matrices and right hand sides; obtaining quantities that are inherent to
+vector-valued functions (i.e. divergences) becomes a little awkward, and others.
+
<a name ="testcase"></a>
-<h3> The testcase </h3>
+<h3> The Testcase </h3>
The testcase that we use for this program consists of the flow around a square obstacle. The geometry is
as follows: