]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Minor edits to the introduction of step-70.
authorWolfgang Bangerth <bangerth@colostate.edu>
Fri, 15 May 2020 19:27:15 +0000 (13:27 -0600)
committerMatthias Maier <tamiko@43-1.org>
Sat, 16 May 2020 17:42:41 +0000 (12:42 -0500)
examples/step-70/doc/intro.dox

index 78d4f5c574ab780158ce2e163acc543b96ace323..9f094b584b97b26a502abfa07210a35666e7387d 100644 (file)
@@ -24,7 +24,7 @@ dependent) difference between the two domains, namely:
 $\Omega\setminus\Omega^{\text{imp}}$.
 
 For rotating impellers, the use of Arbitrary Lagrangian Eulerian formulations
-(in which the fluid domain is smoothly deformed to follow the deformations
+(in which the fluid domain -- along with the mesh! -- is smoothly deformed to follow the deformations
 of the immersed solid) is not possible, unless only small times (i.e.,
 small fluid domain deformations) are considered. If one wants to track the
 evolution of the flow across multiple rotations of the impellers, the resulting
@@ -36,17 +36,20 @@ grid (that may or may not be locally refined in time to better capture the solid
 motion) is coupled with a rotating, independent, grid.
 
 In order to maintain the same notations used in step-60, we use $\Omega$ to
-denote the domain in $R^{\text{spacedim}}$ representing the container of both
-the fluid and the impeller, and we use $\Gamma$ in $R^{\text{dim}}$ to denote
+denote the domain in ${\mathbb R}^{\text{spacedim}}$ representing the container of both
+the fluid and the impeller, and we use $\Gamma$ in ${\mathbb R}^{\text{dim}}$ to denote
 either the full impeller (when its `spacedim` measure is non-negligible, i.e.,
 when we can represent it as a grid of dimension `dim` equal to `spacedim`),
 a co-dimension one representation of a thin impeller, or just the boundary of
 the full impeller.
 
 The domain $\Gamma$ is embedded in $\Omega$ ($\Gamma \subseteq \Omega$) and it
-is non-matching. We solve a partial differential equation on $\Omega$,
+is non-matching: It does not, in general, align with any of the
+features of the volume mesh. We solve a partial differential equation on $\Omega$,
 enforcing some conditions on the solution of the problem on the embedded
-domain $\Gamma$ by some penalization techniques.
+domain $\Gamma$ by some penalization techniques. In the current case,
+the condition is that the velocity of the fluid at points on $\Gamma$
+equal the velocity of the solid impeller at that point.
 
 The technique we describe here is presented in the literature using one of many
 names: the <b>immersed finite element method</b> and the <b>fictitious boundary
@@ -73,7 +76,7 @@ step-60).
 There are two interesting scenarios that occur when one wants to enforce
 conditions on the embedded domain $\Gamma$:
 
-- the geometrical dimension `dim` of the embedded domain $\Gamma$ is the same of
+- The geometrical dimension `dim` of the embedded domain $\Gamma$ is the same of
 the domain $\Omega$ (`spacedim`), that is, the spacedim-dimensional measure of
 $\Gamma$ is not zero. In this case, the imposition of the Dirichlet boundary
 boundary condition on $\Gamma$ is done through a volumetric penalization. If the
@@ -82,17 +85,20 @@ to as $\mathcal{L}^2$ penalization whereas if the penalization depends
 on both the velocity and its gradient, it is an $\mathcal{H}^1$ penalization.
 The case of the $\mathcal{L}^2$ penalization is very similar to a Darcy-type
 approach. Both $\mathcal{L}^2$ and $\mathcal{H}^1$ penalizations have been
-analyzed extensively (ex: see Angot 1999).
+analyzed extensively (see, for example, Angot 1999).
 
-- the embedded domain $\Gamma$ has an intrinsic dimension `dim` which is smaller
+- The embedded domain $\Gamma$ has an intrinsic dimension `dim` which is smaller
 than that of $\Omega$ (`spacedim`), thus its spacedim-dimensional measure is
 zero; for example it is a curve embedded in a two dimensional domain, or a
-surface embedded in a three-dimensional domain. In this case, the boundary
-condition is imposed weekly on $\Gamma$ by applying the Nitsche method (see
+surface embedded in a three-dimensional domain. This is of course
+physically impossible, but one may consider very thin sheets of metal
+moving in a fluid as essentially lower-dimensional if the thickness of
+the sheet is negligible. In this case, the boundary
+condition is imposed weakly on $\Gamma$ by applying the <a href="https://en.wikipedia.org/wiki/Joachim_Nitsche">Nitsche</a> method (see
 Freund, 1995).
 
 Both approaches have very similar requirements and result in highly
-similar formulations, therefore we treat them almost in the same way.
+similar formulations. Thus, we treat them almost in the same way.
 
 In this tutorial program we are not interested in further details on $\Gamma$:
 we assume that the dimension of the embedded domain (`dim`) is always smaller by
@@ -109,6 +115,9 @@ regular function $g$ on $\Gamma$, find the solution $(\textbf{u},p)$ to
   \textbf{u} &=& 0 \text{ on } \partial\Omega.
 @f}
 
+This equation, which we have normalized by scaling the time units in
+such a way that the viscosity has a numerical value of 1, describes
+slow, viscous flow such as honey or lava.
 The main goal of this tutorial is to show how to impose the velocity field
 condition $\mathbf{u} = \mathbf{g}$ on a non-matching $\Gamma$ in a weak way,
 using a penalization method. A more extensive discussion of the Stokes
@@ -157,17 +166,19 @@ we obtain the following variational problem:
  - (q, \textrm{div}\; \textbf{u})_{\Omega}&=& 0
 @f}
 
-where $(\cdot, \cdot)_{\Omega}$ represents the $L^2$ scalar product.
+where $(\cdot, \cdot)_{\Omega}$ represents the $L^2$ scalar
+product. This is the same variational form used in step-22.
 
 This variational formulation does not take into account the embedded domain.
-Contrarily to step-60, we do not enforce strongly the constraints of
+Contrary to step-60, we do not enforce strongly the constraints of
 $\textbf{u}$ on $\Gamma$, but enforce them weakly via a penalization term.
 
 The analysis of this weak imposition of the boundary condition depends on the
-spacedim-dimensional measure of $\Gamma$ as ether positive (if `dim` is equal
+spacedim-dimensional measure of $\Gamma$ as either positive (if `dim` is equal
 to `spacedim`) or zero (if `dim` is smaller than `spacedim`). We discuss both
 scenarios.
 
+
 <h4>Co-dimension one case</h4>
 
 In this case, we assume that $\Gamma$ is the boundary of the actual impeller,
@@ -177,24 +188,25 @@ considering a weak imposition of the Dirichlet boundary condition on $\Gamma$,
 following the Nitsche method. This is achieved by using the following modified formulation
 on the fluid domain, where no strong conditions on the test functions on $\Gamma$ are imposed:
 
-@f{eqnarray*}
-(\nabla \textbf{v}, \nabla \textbf{u})_{\Omega\setminus\Omega^{\text{imp}}} - (\textrm{div}\;  \textbf{v}, p)_{\Omega\setminus\Omega^{\text{imp}}}
+@f{multline*}{
+(\nabla \textbf{v}, \nabla \textbf{u})_{\Omega\setminus\Omega^{\text{imp}}} - (\textrm{div}\;  \textbf{v}, p)_{\Omega\setminus\Omega^{\text{imp}}}
   - (q, \textrm{div}\; \textbf{u})_{\Omega\setminus\Omega^{\text{imp}}} \\
-  - &(\textbf{v},\nabla \textbf{u} \cdot \textbf{n})_{\Gamma}
+  - (\textbf{v},\nabla \textbf{u} \cdot \textbf{n})_{\Gamma}
   + (\textbf{v}\cdot \textbf{n},p)_{\Gamma} \\
- -  &(\nabla\textbf{v}\cdot \textbf{n},\textbf{u})_{\Gamma}
- + (q, \textbf{u} \cdot n)_{\Gamma}
+ -  (\nabla\textbf{v}\cdot \textbf{n},\textbf{u})_{\Gamma}
+ + (q, \textbf{u} \cdot \textbf{n})_{\Gamma}
  + \beta (\textbf{v},\textbf{u})_{\Gamma} \\
-=  - &(\nabla\textbf{v}\cdot \textbf{n},\textbf{g})_{\Gamma} + (q, \textbf{g} \cdot n)_{\Gamma}
+=  - (\nabla\textbf{v}\cdot \textbf{n},\textbf{g})_{\Gamma} + (q, \textbf{g} \cdot \textbf{n})_{\Gamma}
  + \beta (\textbf{v},\textbf{g})_{\Gamma}.
 @f}
 
-The integrals over $\Gamma$ are face integrals. It can be shown (see  Freund,
+The integrals over $\Gamma$ are lower-dimensional integrals. It can be shown (see  Freund,
 1995) that there exists a positive constant
 $C_1$ so that if $\beta > C_1$, the weak imposition of the boundary will
 be consistent and stable. The first two additional integrals on $\Gamma$ (the
 second line in the equation above) appear naturally after integrating by parts,
-when one does not assume that $\mathbf{v}$ is zero on $\Gamma$.
+when one does not assume that $\mathbf{v}$ is zero on
+$\Gamma$.
 
 The third line in the equation above contains two terms that are added to ensure
 consistency of the weak form, and a stabilization term, that is there to enforce
@@ -209,28 +221,34 @@ consistent.
 
 We note that an alternative (non-symmetric) formulation can be used :
 
-@f{eqnarray*}
-(\nabla \textbf{v}, \nabla \textbf{u})_{\Omega\setminus\Omega^{\text{imp}}} - & (\textrm{div}\;  \textbf{v}, p)_{\Omega\setminus\Omega^{\text{imp}}}
+@f{multline*}{
+(\nabla \textbf{v}, \nabla \textbf{u})_{\Omega\setminus\Omega^{\text{imp}}} -  (\textrm{div}\;  \textbf{v}, p)_{\Omega\setminus\Omega^{\text{imp}}}
   - (q, \textrm{div}\; \textbf{u})_{\Omega\setminus\Omega^{\text{imp}}} \\
-  &-(\textbf{v},\nabla \textbf{u} \cdot \textbf{n})_{\Gamma}
+  -(\textbf{v},\nabla \textbf{u} \cdot \textbf{n})_{\Gamma}
   + (\textbf{v}\cdot \textbf{n},p)_{\Gamma} \\
-   &+(\nabla\textbf{v}\cdot \textbf{n},\textbf{u})_{\Gamma}
- - (q, \textbf{u} \cdot n)_{\Gamma}
+   +(\nabla\textbf{v}\cdot \textbf{n},\textbf{u})_{\Gamma}
+ - (q, \textbf{u} \cdot \textbf{n})_{\Gamma}
  + \beta (\textbf{v},\textbf{u})_{\Gamma} \\
-=   &(\nabla\textbf{v}\cdot \textbf{n},\textbf{g})_{\Gamma} - (q, \textbf{g} \cdot n)_{\Gamma}
- + \beta (\textbf{v},\textbf{g})_{\Gamma},
+=   (\nabla\textbf{v}\cdot \textbf{n},\textbf{g})_{\Gamma} - (q, \textbf{g} \cdot \textbf{n})_{\Gamma}
+ + \beta (\textbf{v},\textbf{g})_{\Gamma}.
 @f}
-in which case the stability and consistency conditions become $\beta > 0$. In
+Note the different sign of the first terms on the third and fourth lines.
+In this case, the stability and consistency conditions become $\beta > 0$. In
 the symmetric case, the value of $\beta$ is dependent on $h$, and it is in
 general chosen such that $\beta = C h^{-1} $ with $h$
 a measure of size of the face being integrated and $C$ a constant such that
-$1 \leq C \leq 10$.
+$1 \leq C \leq 10$. This is as one usually does with the Nitsche
+penalty method to enforcing Dirichlet boundary conditions.
 
+The non-symmetric approach, on the other hand, is related to how one
+enforced continuity for the non-symmetric interior penalty method for
+discontinuous Galerkin methods (the "NIPG" method @cite Riviere1999).
 Even if the non-symmetric case seems advantageous w.r.t.
 possible choices of stabilization parameters, we opt for the symmetric
 discretization, since in this case it can be shown that the dual problem is
 also consistent, leading to a solution where not only the energy norm of the
-solution converges with the correct order, but also its $L^2$ norm.
+solution converges with the correct order, but also its $L^2$
+norm. Furthermore, the resulting matrix remains symmetric.
 
 The above formulation works under the assumption that the domain is discretized
 exactly. However, if the deformation of the impeller is a rigid body
@@ -242,19 +260,24 @@ $\Gamma$, using the same penalization technique, and testing with test
 functions $\mathbf{v}$ which are globally continuous over $\Omega$.
 
 This results in the following (intermediate) formulation:
-@f{eqnarray*}
-(\nabla \textbf{v}, \nabla \textbf{u})_{\Omega} - (\textrm{div}\;  \textbf{v}, p)_{\Omega}
+@f{multline*}{
+(\nabla \textbf{v}, \nabla \textbf{u})_{\Omega} - (\textrm{div}\;  \textbf{v}, p)_{\Omega}
   - (q, \textrm{div}\; \textbf{u})_{\Omega} \\
-  - &(\textbf{v},  \lbrack \nabla \textbf{u} \rbrack \cdot \textbf{n})_{\Gamma}
+  - (\textbf{v},  \lbrack \nabla \textbf{u} \rbrack \cdot \textbf{n})_{\Gamma}
   + (\textbf{v}\cdot \textbf{n},\lbrack p \rbrack )_{\Gamma} \\
- -  &(\lbrack \nabla\textbf{v} \rbrack \cdot \textbf{n},\textbf{u})_{\Gamma}
+ -  (\lbrack \nabla\textbf{v} \rbrack \cdot \textbf{n},\textbf{u})_{\Gamma}
  + (\lbrack q \rbrack, \textbf{u} \cdot n)_{\Gamma}
  + 2\beta (\textbf{v},\textbf{u})_{\Gamma} \\
-=  - &(\lbrack \nabla\textbf{v}\rbrack\cdot \textbf{n},\textbf{g})_{\Gamma} + (\lbrack q\rbrack, \textbf{g} \cdot n)_{\Gamma}
+=  - (\lbrack \nabla\textbf{v}\rbrack\cdot \textbf{n},\textbf{g})_{\Gamma} + (\lbrack q\rbrack, \textbf{g} \cdot n)_{\Gamma}
  + 2\beta (\textbf{v},\textbf{g})_{\Gamma},
 @f}
 where the jump terms, denoted with $\lbrack \cdot \rbrack$, are computed with
-respect to a fixed orientation of the normal vector $\textbf{n}$.
+respect to a fixed orientation of the normal vector $\textbf{n}$. The
+factor of 2 appears in front of $\beta$ since we see every part of
+$\Gamma$ twice, once from within the fluid and once from within the
+obstacle moving around in it. (For all of the other integrals over
+$\Gamma$, we visit each part of $\Gamma$ twice, but with opposite
+signs, and consequently get the jump terms.)
 
 Here we notice that, unlike in discontinuous Galerkin methods, the test
 and trial functions are continuous across $\Gamma$. Moreover, if $\Gamma$ is
@@ -265,7 +288,7 @@ number of points, all the contributions on $\Gamma$, with the exception of
 the stabilization ones, can be neglected from the formulation, resulting in
 the following final form of the variational formulation:
 
-@f{eqnarray*}
+@f{multline*}{
 (\nabla \textbf{v}, \nabla \textbf{u})_{\Omega} - & (\textrm{div}\;  \textbf{v}, p)_{\Omega}
   - (q, \textrm{div}\; \textbf{u})_{\Omega}  + 2\beta (\textbf{v},\textbf{u})_{\Gamma} \\
 =  2\beta (\textbf{v},\textbf{g})_{\Gamma}.
@@ -283,11 +306,13 @@ method. Furthermore, as in step-60, we still need to integrate over the
 non-matching embedded grid in order to construct the boundary term necessary
 to impose the boundary condition over $\Gamma$.
 
+
 <h4>Co-dimension zero case</h4>
 
-In this case, $\Gamma$ has the same dimension, but is embedded into $\Omega$.
+In this case, $\Gamma$ has the same dimension, but is embedded into
+$\Omega$. We can think of this as a thick object moving around in the fluid.
 In the case of $\mathcal{L}^2$ penalization, the additional penalization
-term, can be interpreted as a Darcy term within $\Gamma$, resulting in:
+term can be interpreted as a Darcy term within $\Gamma$, resulting in:
 
 @f{eqnarray*}
 (\nabla \textbf{v}, \nabla \textbf{u})_{\Omega} - & (\textrm{div}\;  \textbf{v}, p)_{\Omega}
@@ -295,7 +320,8 @@ term, can be interpreted as a Darcy term within $\Gamma$, resulting in:
 =  \beta (\textbf{v},\textbf{g})_{\Gamma}.
 @f}
 
-The $\mathcal{L}^2$ penalization consists in adding a volumetric term that
+Here, integrals over $\Gamma$ are simply integrals over a part of the volume.
+The $\mathcal{L}^2$ penalization thus consists in adding a volumetric term that
 constrains the velocity of the fluid to adhere to the velocity of the rigid body
 within $\Gamma$. Also in this case, $\beta$ must be chosen sufficiently large
 in order to ensure that the Dirichlet boundary condition in $\Gamma$ is
@@ -319,12 +345,14 @@ Notice that the $L^2$ penalization (`dim` equal to `spacedim`) and the Nitsche
 penalization (`dim` equal to `spacedim-1`) result in the exact same numerical
 implementation, thanks to the dimension independent capabilities of deal.II.
 
-<h4>Representation of $\Omega$ and $\Gamma$</h4>
+
+<h4>Representation of Ω and Γ</h4>
 
 In this tutorial, both the embedded grid $\Gamma$ and the embedding
 grid are described using a parallel::distributed::Triangulation. These two
-triangulations can be built from the GridGenerator namespace or by reading
-a mesh file produced with another application (e.g. GMSH). This is slightly
+triangulations can be built from functions in the GridGenerator namespace or by reading
+a mesh file produced with another application (e.g. GMSH, see the
+discussion in step-49). This is slightly
 more general than what was previously done in step-60.
 
 The addition of the immersed boundary method, whether
@@ -349,7 +377,7 @@ and compute the integral on $\hat K$ using a quadrature formula. For example:
 
 Computing this sum is non-trivial because we have to evaluate $(v_j \circ F_{K})
 (\hat x_i)$. In general, if $\Gamma$ and $\Omega$ are not aligned, the point
-$F_{K}(\hat x_i)$ is completely arbitrary with respect to $\Omega$, and unless
+$y_i = F_{K}(\hat x_i)$ is completely arbitrary with respect to $\Omega$, and unless
 we figure out a way to interpolate all basis functions of $V_h(\Omega)$ on an
 arbitrary point on $\Omega$, we cannot compute the integral needed.
 
@@ -359,19 +387,25 @@ taken (as shown in the picture below):
 
 - For a given cell $K$ in $\Gamma$ compute the real point $y_i \dealcoloneq F_{K} (\hat
 x_i)$, where $x_i$ is one of the quadrature points used for the integral on $K
-\subseteq \Gamma$.
+\subseteq \Gamma$. This is the easy part:
+FEValues::quadrature_point() gives us the real-space locations of all
+quadrature points.
 
 - Find the cell of $\Omega$ in which $y_i$ lies. We shall call this element $T$.
 
-- To evaluate the basis function use the inverse of the mapping $G_T$ that
-transforms the reference element $\hat T$ into the element $T$: $v_j(y_i) = \hat
-v_j \circ G^{-1}_{T} (y_i)$.
+- Find the reference coordinates within $T$ of $y_i$. For this, we
+need the inverse of the mapping $G_T$ that
+transforms the reference element $\hat T$ into the element $T$: $\hat y_i = G^{-1}_{T} (y_i)$.
+
+- Evaluate the basis function $v_j$ of the $\Omega$ mesh at this
+  point $\hat y_i$. This is, again, relatively simple using FEValues.
+
 
 <p align="center"> <img
   src="https://www.dealii.org/images/steps/developer/step-60.C_interpolation.png"
   alt=""> </p>
 
-In step-60, the three steps above were computed by calling, in turn,
+In step-60, the second through fourth steps above were computed by calling, in turn,
 
 - GridTools::find_active_cell_around_point(), followed by
 
@@ -396,13 +430,13 @@ in this step.
 When looking at the formulation for both the penalization approach ($\mathcal{L}^2$ or $\mathcal{H}^1$)
 and the Nitsche method, we can come to the conclusion that we only require limited information
 related to $\Gamma$ on its quadrature points that is:
-- Their reference location w.r.t. to the embedding fluid cell
-- The weights times the jacobian of the transformation, i.e JxW.
+- Their reference location with regard to the embedding fluid cell;
+- The weights times the Jacobian of the transformation, i.e. `JxW`.
 
 Since both of these properties are point-properties (or point-vectors) that are
 attached to the solid material, they can be idealized as a set of disconnected
 infinitesimally small particles, which carry the required information with the
-movement of the solid. Deal.II already possesses the ability to distribute and
+movement of the solid. deal.II already possesses the ability to distribute and
 store such a set of particles in large-scale parallel computations in form of
 the ParticleHandler class (for details on the implementation see Gassmöller et
 al., 2018), and we will make use of this functionality in this tutorial.
@@ -413,7 +447,7 @@ Thus, the approach taken in this step is as follows:
 - Call the Particles::ParticleHandler::insert_global_particles() function,
   to distribute the particles across processors, *following the solid
   triangulation*
-- Attach the quadrature weight as property to each Particles::Particle object.
+- Attach the quadrature weight as property to each Particles::Particle object.
 
 This structure is relatively expensive to generate, but must only be generated
 once per simulation. Once the Particles::ParticleHandler is generated and the
@@ -428,7 +462,10 @@ Since the Particles::ParticleHandler can manage the exchange of particles from
 one processor to the other, the embedded
 triangulation can be moved or deformed by displacing the particles.
 The only constraint associated with this displacement is that particles should
-be displaced by a distance that is no larger than the size of a cell.
+be displaced by a distance that is no larger than the size of one
+cell. That's because that is the limit to which
+Particles::ParticleHandler can track which cell a particle that leaves
+its current cell now resides in.
 
 Once the entire problem (the Stokes problem and the immersed boundary
 imposition) is assembled,
@@ -463,27 +500,29 @@ demonstrating the time-reversibility of the flow.
 
 <h3>References</h3>
 
-- Freund, J., Stenberg, R. (1995). "On weakly imposed boundary conditions for
+<ul>
+<li> Freund, J., Stenberg, R. (1995). "On weakly imposed boundary conditions for
   second order problems". Proceedings of the Ninth International Conference on
   Finite Elements in Fluids. 327-336.
 
-- Angot, Philippe, Charles-Henri Bruneau and Pierre Fabrie. 1999. "A penalization
+<li> Angot, Philippe, Charles-Henri Bruneau and Pierre Fabrie. 1999. "A penalization
   method to take into account obstacles in incompressible viscous flows."
   Numerische Mathematik 81.4 : 497-520.
 
-- Glowinski, R., T.-W. Pan, T.I. Hesla, and D.D. Joseph. 1999. “A Distributed
+<li> Glowinski, R., T.-W. Pan, T.I. Hesla, and D.D. Joseph. 1999. “A Distributed
   Lagrange Multiplier/fictitious Domain Method for Particulate Flows.”
   International Journal of Multiphase Flow 25 (5). Pergamon: 755–94.
 
-- Boffi, D., L. Gastaldi, L. Heltai, and C.S. Peskin. 2008. “On the
+<li> Boffi, D., L. Gastaldi, L. Heltai, and C.S. Peskin. 2008. “On the
   Hyper-Elastic Formulation of the Immersed Boundary Method.” Computer Methods
   in Applied Mechanics and Engineering 197 (25–28).
 
-- Heltai, L., and F. Costanzo. 2012. “Variational Implementation of Immersed
+<li> Heltai, L., and F. Costanzo. 2012. “Variational Implementation of Immersed
   Finite Element Methods.” Computer Methods in Applied Mechanics and Engineering
   229–232.
 
-- Gassmöller, R., H. Lokavarapu, E. Heien, E. G. Puckett, and W. Bangerth.
+<li> Gassmöller, R., H. Lokavarapu, E. Heien, E. G. Puckett, and W. Bangerth.
   2018.  "Flexible and Scalable Particle‐in‐Cell Methods With Adaptive Mesh
   Refinement for Geodynamic Computations." Geochemistry, Geophysics, Geosystems
   19(9).  3596-3604.
+</ul>

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.