New generated Stokes operators
This issue shall document and discuss the introduction of the Stokes composite operators that are built from the new generated operators. Also point out if I made a mistake and feel free to correct this description.
List of Stokes operators
It seems(?) that we are particularly interested in three types of Stokes operators, particularly differing in the momentum terms.
Let
K = \begin{bmatrix}
A & B^T \\
B & 0
\end{bmatrix}
be the discrete Stokes operator for a \mathbb{P}_2-\mathbb{P}_1
Taylor-Hood discretization.
The operators differ in the way the A
block is defined. Let's summarize.
Proposed name (prefix) | Strong form of momentum term ("velocity part") | Weak form that defines A
|
---|---|---|
P2P1StokesConstant |
- \Delta u |
\displaystyle\int_\Omega \nabla u : \nabla v |
P2P1StokesEpsilon |
\displaystyle- \nabla \cdot \Big(2 \mu \varepsilon(u)\Big) |
\displaystyle\int_\Omega 2 \mu \Big(\varepsilon(u) : \varepsilon(v)\Big) |
P2P1StokesFull |
\displaystyle - \nabla \cdot \Big( 2 \mu \varepsilon(u) \Big) - \frac{2}{\mathrm{dim}} \nabla (\nabla \cdot u) |
\displaystyle\int_\Omega 2 \mu \Big(\varepsilon(u) : \varepsilon(v)\Big) - \frac{2}{\mathrm{dim}} (\nabla \cdot u) \cdot (\nabla \cdot v) |
where \varepsilon(w) = \frac{1}{2} (\nabla w + (\nabla w)^T)
.
\mu = \mu(x)
will eventually be a scalar finite element function (in \mathbb{P}_2
) that is supplied by the user.
Possible file structure
I propose the following structure - open to suggestions.
hyteg/src/
- hyteg/
- terraneo/
- ...
- hyteg_operators/ # submodule with generated operators
- hyteg_operators_composites/ # new module
- stokes/
- viscousblock/
- P2ViscousBlockLaplaceOperator.*pp
- P2ViscousBlockEpsilonOperator.*pp
- P2ViscousBlockFullOperator.*pp
- divergence/
- P2VectortoP1DivergenceOperator.*pp
- gradient/
- P1toP2VectorGradientOperator.*pp
- P2P1StokesConstantOperator.*pp
- P2P1StokesEpsilonOperator.*pp
- P2P1StokesFullOperator.*pp
- ...
Blending
For each relevant blending map, additional operators will be added to the respective files. For instance
// P2P1StokesEpsilonOperator.hpp
class P2P1StokesEpsilonOperator { ... };
class P2P1StokesEpsilonIcosahedralShellMapOperator { ... };
Some notes
- I am not sure yet which quadrature rules to choose for the blending cases. For now I chose rules that integrate degree 3 polynomials exactly for blending, and degree 2 polynomials without blending (the latter should be sufficient to maintain convergence order(?)). But we may want to discuss this.
- The new composites shall only use/depend on the newly generated operators and never introduce any of the existing constant stencil ops, elementwise ops etc. If something is missing, we should adapt the HFG.
- Lumped and diagonal operators are already in the making. Those are needed for all sorts of preconditioners.
- I currently only consider Taylor-Hood-style operators. Stabilized
\mathbb{P}_1-\mathbb{P}_1
operators should be straightforward to generate, too. - In the future, we want to enable the HFG to generate the vector-operators directly such that they do not have to be composed manually. This will likely remove the need for all operators in the subdirectories
viscousblock
,gradient
,divergence
, etc. But for starters, we will probably go with the composed versions.
Checklist
-
generate scalar block operators -
compile composites for viscous blocks (!710 (merged)) -
compile div and grad blocks (!710 (merged)) -
compile Stokes operators (!710 (merged)) -
tests (!710 (merged)) -
compile/implement diagonal operators to use preconditioners (moved to #250)