<h4>Variational versus collocation-type discretizations</h4>
-Following Step-9, Step-12, and Step-33, at this point it might look tempting
-to base a discretization of Euler's equations on a (semi-discrete) variational
-formulation:
+Following Step-9, Step-12, Step-33, and Step-67, at this point it might look
+tempting to base a discretization of Euler's equations on a (semi-discrete)
+variational formulation:
@f{align*}
(\partial_t\mathbf{u}_{h},\textbf{v}_h)_{L^2(\Omega)}
- ( \mathbb{f}(\mathbf{u}_{h}) ,\text{grad} \, \textbf{v}_{h})_{L^2(\Omega)}
have proven to be some of the best approaches for simulations in the subsonic
shockless regime and similarly benign situations.
-However, in the transonic and supersonic regime, and shock-hydrodynamics
+<!-- In particular, tutorial Step-67 focuses on Euler's equation of gas
+dynamics in the subsonic regime using dG techniques. -->
+
+However, in the transonic and supersonic regimes, and shock-hydrodynamics
applications the use of variational schemes might be questionable. In fact,
at the time of this writing, most shock-hydrodynamics codes are still
firmly grounded on finite volume methods. The main reason for failure of
to compute the block-matrices (required in order to advance the solution)
with relative ease. However, the interactions that have to be computed in
the context of time-explicit collocation-type schemes (such as finite
-differences and/or the scheme presented in this tutorial) can be are
+differences and/or the scheme presented in this tutorial) can be
better described as <i>interactions between nodes</i> (not between DOFs).
In addition, in our case we do not solve a linear equation in order to
advance the solution. This leaves very little reason to use vector-valued
(\mathbf{U}_i^{n},\mathbf{U}_j^{n}, \textbf{n}_{ij}),
\lambda_{\text{max}} (\mathbf{U}_j^{n}, \mathbf{U}_i^{n},
\textbf{n}_{ji}) \} \|\mathbf{c}_{ij}\|$ if $i \not = j$ is the so
- called <i>graph viscosity</i>. The graph viscosity serves as a
+ called <i>graph-viscosity</i>. The graph-viscosity serves as a
stabilization term, it is somewhat the discrete counterpart of
$\epsilon \Delta \mathbf{u}$ that appears in the notion of viscosity
solution described above. We will base our construction of $d_{ij}$ on
update. The core principle remains unchanged, though: we do not loop over
cells but rather over all edges of the sparsity graph.
-@note It is not uncommon to encounter such fully algebraic schemes (i.e.
+@note It is not uncommon to encounter such fully-algebraic schemes (i.e.
no bilinear forms, no cell loops, and no quadrature) outside of the finite
element community in the wider CFD community. There is a rich history of
application of this kind of schemes, also called <i>edge-based</i> or
@cite Rainald2008 for a historical overview). However, it is important to
highlight that the algebraic structure of the scheme (presented in this
tutorial) and the node-loops are not just a performance gimmick. Actually, the
-structure of this scheme was borne out of theoretical necessity: the proof of
+structure of this scheme was born out of theoretical necessity: the proof of
pointwise stability of the scheme hinges on the specific algebraic structure of
the scheme. In addition, it is not possible to compute the algebraic
viscosities $d_{ij}$ using cell-loops since they depend nonlinearly on
@f{align*}
\mathbf{m}_i \dealcoloneq \mathbf{m}_i - (\widehat{\boldsymbol{\nu}}_i
\cdot \mathbf{m}_i) \widehat{\boldsymbol{\nu}}_i \ \
- \text{for all }\mathbf{x}_i \in \partial\Omega^r
+ \text{where} \ \
+ \widehat{\boldsymbol{\nu}}_i \dealcoloneq
+ \frac{\int_{\partial\Omega} \phi_i \widehat{\boldsymbol{\nu}} \,
+ \, \mathrm{d}\mathbf{s}}{\big|\int_{\partial\Omega} \phi_i
+ \widehat{\boldsymbol{\nu}} \, \mathrm{d}\mathbf{s}\big|}
+ \ \ \text{for all }\mathbf{x}_i \in \partial\Omega^r
+ \ \ \ \ \boldsymbol{(1)}
@f}
- that removes the normal component of $\mathbf{m}$. Here the definition of
- nodal normal $\widehat{\boldsymbol{\nu}}_i$ is very much arbitrary (there is
- no unique definition) but it should be consistent upon refinement with the
- underlying geometry.
+ that removes the normal component of $\mathbf{m}$. This is a somewhat
+ naive idea that preserves a few fundamental properties of the PDE as we
+ explain below.
This is approach is usually called "explicit treatment of boundary conditions".
The well seasoned finite element person might find this approach questionable.
@f{align*}
\int_{\Omega} \rho(\mathbf{x},t_2) \, \mathrm{d}\mathbf{x} =
\int_{\Omega} \rho(\mathbf{x},t_1) \, \mathrm{d}\mathbf{x} \ , \ \
-\int_{\Omega} E(\mathbf{x},t_2) \, \mathrm{d}\mathbf{x} =
-\int_{\Omega} E(\mathbf{x},t_1) \, \mathrm{d}\mathbf{x} \ , \ \
-\int_{\Omega} \mathbf{m}(\mathbf{x},t_2) \, \mathrm{d}\mathbf{x} =
-\int_{\Omega} \mathbf{m}(\mathbf{x},t_1) \, \mathrm{d}\mathbf{x}
+\int_{\Omega} \mathbf{m}(\mathbf{x},t_2) \, \mathrm{d}\mathbf{x}
+ \int_{t_1}^{t_2} \! \int_{\partial\Omega} p \boldsymbol{\nu} \,
-\mathrm{d}\mathbf{s} \mathrm{d}t\, .
+\mathrm{d}\mathbf{s} \mathrm{d}t =
+\int_{\Omega} \mathbf{m}(\mathbf{x},t_1) \,
+\mathrm{d}\mathbf{x} \ , \ \
+\int_{\Omega} E(\mathbf{x},t_2) \, \mathrm{d}\mathbf{x} =
+\int_{\Omega} E(\mathbf{x},t_1) \, \mathrm{d}\mathbf{x} \ \ \ \
+\boldsymbol{(2)}
@f}
-Note that momentum is not a conserved quantity (interaction with walls leads to
-momentum gain/loss). Even though we will not use reflecting boundary conditions
-in the entirety of the domain, we would like to know that our implementation of
+Note that momentum is NOT a conserved quantity (interaction with walls leads to
+momentum gain/loss): however $\mathbf{m}$ has to satisfy a momentum balance.
+Even though we will not use reflecting boundary conditions in the entirety of
+the domain, we would like to know that our implementation of reflecting
boundary conditions is consistent with the conservation properties mentioned
-above. In order to guarantee such conservation property it is necessary to
-modify the values of the vectors $\mathbf{c}_{ij}$ as follows
+above. In particular, if we use the projection $\boldsymbol{(1)}$ in the
+entirety of the domain the following discrete mass-balance can be guaranteed:
@f{align*}
-\mathbf{c}_{ij} \dealcoloneq \mathbf{c}_{ij} + \int_{\partial \Omega^r}
-(\boldsymbol{\nu}_j - \boldsymbol{\nu}(\mathbf{s})) \phi_i \phi_j \,
-\mathrm{d}\mathbf{s}
-\ \ \text{whenever both }\mathbf{x}_i\text{ and }\mathbf{x}_j \text{ lie on the
-boundary}
+\sum_{i \in \mathcal{V}} m_i \rho_i^{n+1} =
+\sum_{i \in \mathcal{V}} m_i \rho_i^{n} \ , \ \
+\sum_{i \in \mathcal{V}} m_i \mathbf{m}_i^{n+1}
++ \tau_n \int_{\partial\Omega} \Big(\sum_{i \in \mathcal{V}} p_i^{n} \phi_i\Big)
+\widehat{\boldsymbol{\nu}} \mathrm{d}\mathbf{s} =
+\sum_{i \in \mathcal{V}} m_i \mathbf{m}_i^{n} \ , \ \
+\sum_{i \in \mathcal{V}} m_i E_i^{n+1} = \sum_{i \in \mathcal{V}} m_i
+E_i^{n} \ \ \ \
+\boldsymbol{(3)}
@f}
-This consistent modification of the $\mathbf{c}_{ij}$ is a direct consequence
-of simple integration by parts arguments, see page 12 of @cite GuermondEtAl2018
-for more details.
+where $p_i$ is the pressure at the nodes that lie at the boundary. Clearly
+$\boldsymbol{(3)}$ is the discrete counterpart of $\boldsymbol{(2)}$. The
+proof of identity $\boldsymbol{(3)}$ is omitted, but we briefly mention that
+it hinges on the definition of the <i>nodal normal</i>
+$\widehat{\boldsymbol{\nu}}_i$ provided in $\boldsymbol{(1)}$. We also note that
+this enforcement of reflecting boundary conditions is different from the one
+originally advanced in @cite GuermondEtAl2018.
projects / dealii.git / commitdiff