examples/navier-stokes.md

@@ -10,6 +10,7 @@ rate-of-strain tensor $`\mathbf{D}(\mathbf{u})=\frac{1}{2}(\nabla\mathbf{u}+\nab
 and external volume force $`\mathbf{f}`$, w.r.t. to boundary conditions $`\mathbf{u}=\mathbf{g}`$ on $`\partial\Omega`$.
 For a coarse overview, see also [Wikipedia](https://en.wikipedia.org/wiki/Navier%E2%80%93Stokes_equations).
 
+
 Discretization
 --------------
 
@@ -34,9 +35,9 @@ with $`0=t^0<t^1<\ldots<t^N=T`$, timestep $`\tau=t^{k+1} - t^k`$, and given init
 
 
 ### Finite-Element formulation
-Let $`\Omega`$ be partitioned into quasi-uniform elements by a conforming triangulation $`\mathcal{T}_h`$. The functions $`\mathbf{u}`$ and $`p`$ are approximated by functions from the discrete functionspace $`\mathbb{T}_h := [\mathbb{V}_{2,\,g_i}]^d\times \mathbb{V}_1`$ with $`\mathbb{V}_{p} := \{v\in C^0(\Omega)\,:\, v|_T\in\mathbb{P}_p(T),\, \forall T\in\mathcal{T}_h\}`$ and $`\mathbb{V}_{p,\,g} := \{v\in \mathbb{V}_p\,:\,\operatorname{tr}_{\partial\Omega}=g\}`$, where $`\mathbb{P}_p`$ is the space of polynomials of order at most $`p`$. The space $`\mathbb{T}_h`$ will be denoted as Taylor-Hood space.
+Let $`\Omega`$ be partitioned into quasi-uniform elements by a conforming triangulation $`\mathcal{T}_h`$. The functions $`\mathbf{u}`$ and $`p`$ are approximated by functions from the discrete functionspace $`\mathbb{T}_h := [\mathbb{V}_{2,\,g_i}]^d\times \mathbb{V}_1`$ with $`\mathbb{V}_{p} := \{v\in C^0(\Omega)\,:\, v|_T\in\mathbb{P}_p(T),\, \forall T\in\mathcal{T}_h\}`$ and $`\mathbb{V}_{p,\,g} := \{v\in \mathbb{V}_p\,:\,\operatorname{tr}_{\partial\Omega} v=g\}`$, where $`\mathbb{P}_p`$ is the space of polynomials of order at most $`p`$. The space $`\mathbb{T}_h`$ will be denoted as Taylor-Hood space.
 
-Let `Grid grid` be the Duen grid type representing the triangulation $`\mathcal{T}_h`$ of the domain $`\Omega`$, then
+Let `Grid grid` be the Dune grid type representing the triangulation $`\mathcal{T}_h`$ of the domain $`\Omega`$, then
 the basis of the Taylor-Hood space can be written as a composition of Lagrange bases:
 
 ```c++
@@ -45,13 +46,94 @@ auto basis = makeBasis(grid.leafGridView(),
     composite(
       power<Grid::dimensionworld>(lagrange<2>(), flatInterleaved()),
       lagrange<1>(),
-      flatLexicographic()
-    )
+      flatLexicographic() )
   );
 ```
-At this point, any Dirichlet boundary conditions is ignored in the definition of the bases. We
+At this point, any Dirichlet boundary conditions are ignored in the definition of the bases. We
 build the Taylor-Hood basis as a composition of a product of $`d`$ $`\mathbb{V}_2`$ bases for
-the velocity components and a $`\mathbb{V}_1`$ basis for the pressure.
+the velocity components and a $`\mathbb{V}_1`$ basis for the pressure. We have to use
+`composite(...)` to combine the different types for velocity and pressure space, while we
+can use `power<d>(...)` for the `d` velocity component spaces.
 
-For efficiency, we use flat blocking strategies `flatInterleaved` and `flatLexicographic`
+For efficiency, we use flat indexing strategies `flatInterleaved` and `flatLexicographic`
 to build a global continuous numbering of the basis functions.
+
+
+Implementation
+--------------
+
+We split the implementation into three parts, 1. the time derivative, 2. the stokes
+operator and 3. any external forces. For addressing the different components of
+the Taylor-Hood basis, we introduce the treepaths `auto _v = Dune::Indices::_0` and
+`auto _p = Dune::Indices::_1` for velocity and pressure, respectively.
+
+### Problem framework
+We have an instationary problem consistenting of a sequence of stationary equations
+for each timestep, thus we combine a `ProblemInstat` with a `Problemstat`.
+
+```c++
+ProblemStat prob("stokes", grid, basis);
+prob.initialize(INIT_ALL);
+
+ProblemInstat probInstat("stokes", prob);
+probInstat.initialize(INIT_UH_OLD);
+```
+For convenience, we use class template argument deduction here.
+
+### Time derivative
+For simplicity of this example, we use a backward Euler time stepping and a linearization
+of the nonlinear advection term:
+
+```c++
+// define a constant fluid density
+double density = 1.0;
+Parameters::get("stokes->density", density);
+
+// a reference to 1/tau
+auto invTau = std::ref(probInstat.invTau());
+
+// <1/tau * u, v>
+auto opTime = makeOperator(tag::testvec_trialvec{},
+  density * invTau);
+prob.addMatrixOperator(opTime, _v, _v);
+
+// <1/tau * u^old, v>
+auto opTimeOld = makeOperator(tag::testvec{},
+  density * invTau * prob.solution(_v));
+prob.addVectorOperator(opTimeOld, _v);
+
+for (int i = 0; i < Grid::dimensionworld; ++i) {
+  // <(u^old * nabla)u_i, v_i>
+  auto opNonlin = makeOperator(tag::test_gradtrial{},
+    density * prob.solution(_v));
+  prob.addMatrixOperator(opNonlin, treepath(_v,i), treepath(_v,i));
+}
+```
+
+### Stokes operator
+The stokes part of the Navier-Stokes equation is a common pattern that emerges
+in several fluid equations.
+```math
+\int_\Omega \nu\nabla\mathbf{u}:\nabla\mathbf{v} + p\,\nabla\cdot\mathbf{v} + q\,\nabla\cdot\mathbf{u}\,\text{d}\mathbf{x}
+```
+for $`\mathbf{u},\mathbf{v}\in V`$ and $`p,q\in Q`$.
+
+```c++
+// define a constant fluid viscosity
+double viscosity = 1.0;
+Parameters::get("stokes->viscosity", viscosity);
+
+for (int i = 0; i < Grid::dimensionworld; ++i) {
+  // <viscosity*grad(u_i), grad(v_i)>
+  auto opL = makeOperator(tag::gradtest_gradtrial{}, viscosity);
+  prob.addMatrixOperator(opL, treepath(_v,i), treepath(_v,i));
+}
+
+// <d_i(v_i), p>
+auto opP = makeOperator(tag::divtestvec_trial{}, 1.0);
+prob.addMatrixOperator(opP, _v, _p);
+
+// <q, d_i(u_i)>
+auto opDiv = makeOperator(tag::test_divtrialvec{}, 1.0);
+prob.addMatrixOperator(opDiv, _p, _v);
+```
\ No newline at end of file