This tutorial presents a first-order scheme for solving compressible
Euler's equations that is based on three ingredients: a
-<i>collocation</i>-type discretization of Euler's equations in context of
-finite elements; a graph-viscosity stabilization based on a
+<i>collocation</i>-type discretization of Euler's equations in the context
+of finite elements; a graph-viscosity stabilization based on a
<i>guaranteed</i> upper bound of the local wave speed; and explicit
time-stepping. As such, the ideas and techniques presented in this tutorial
step are drastically different from those used in Step-33, which focuses on
the use of automatic differentiation. From a programming perspective this
tutorial will focus on a number of techniques found in large-scale
-computations: hybrid thread-MPI parallelization; efficient
-local numbering of degrees of freedom; concurrent post-processing and
-write-out of results using worker threads; as well as checkpointing and
-restart.
+computations: hybrid thread-MPI parallelization; efficient local numbering
+of degrees of freedom; concurrent post-processing and write-out of results
+using worker threads; as well as checkpointing and restart.
It should be noted that first-order schemes in the context of hyperbolic
conservation laws require prohibitively many degrees of freedom to resolve
The compressible Euler's equations of gas dynamics are written in
conservative form as follows:
-
@f{align}
\mathbf{u}_t + \text{div} \, \mathbb{f}(\mathbf{u}) = \boldsymbol{0} ,
@f}
-
where $\mathbf{u}(\textbf{x},t):\mathbb{R}^{d} \times \mathbb{R}
\rightarrow \mathbb{R}^{d+2}$, and $\mathbb{f}(\mathbf{u}):\mathbb{R}^{d+2}
\rightarrow \mathbb{R}^{(d+2) \times d}$, and $d \geq 1$ is the space
denotes the density, $\textbf{m} \in \mathbb{R}^d$ is the momentum, and $E
\in \mathbb{R}^+$ is the total energy of the system. The flux of the system
$\mathbb{f}(\mathbf{u})$ is defined as
-
@f{align*}
\mathbb{f}(\textbf{u})
=
\tfrac{\textbf{m}^\top}{\rho} (E + p)
\end{bmatrix},
@f}
-
where $\mathbb{I} \in \mathbb{R}^{d \times d}$ is the identity matrix and
$\otimes$ denotes the tensor product. Here, we have introduced the pressure
$p$ that, in general, is defined by an closed-form equation of state.
\tfrac{|\textbf{m}|_{\ell^2}^2}{2\,\rho}
\Big),
@f}
-
where the factor $\gamma \in (1,5/3]$ denotes the
<a href="https://en.wikipedia.org/wiki/Heat_capacity_ratio">ratio of
specific heats</a>, and $|\cdot|_{\ell^2}$ denotes the Euclidian norm.
solution is formally defined as $\mathbf{u} := \lim_{\epsilon \rightarrow
0^+} \mathbf{u}^{\epsilon}$ where $\mathbf{u}^{\epsilon}$ is the solution
of the parabolic regularization
-
@f{align}
\mathbf{u}_t^{\epsilon} + \text{div} \, \mathbb{f}(\mathbf{u}^{\epsilon})
= {\epsilon} \Delta \mathbf{u}^{\epsilon}.
@f}
-
Such solutions, which are understood as the solution recovered in the
zero-viscosity limit, are often refered to as <i>viscosity solutions</i>.
Global existence and uniqueness of such solutions is a widely open issue.
However, we know at least that if such viscosity solutions exists they have
to satisfy the constraint $\textbf{u}(\mathbf{x},t) \in \mathcal{B}$ for
all $\mathbf{x} \in \Omega$ and $t \geq 0$ where
-
@f{align}
\mathcal{B} = \big\{ \textbf{u} =
[\rho, \textbf{m},E]^{\top} \in \mathbb{R}^{d+2} \, \big |
s(\mathbf{u}) \geq \min_{x \in \Omega} s(\mathbf{u}_0(\mathbf{x}))
\big\}.
@f}
-
Here, $s(\mathbf{u})$ denotes the specific entropy
-
@f{align}
s(\mathbf{u}) = \ln \Big(\frac{p(\mathbf{u})}{\rho^{\gamma}}\Big).
@f}
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:
-
@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)}
+ s_h(\mathbf{u}_{h},\textbf{v}_h)_{L^2(\Omega)} = \boldsymbol{0}
\quad\forall \textbf{v}_h \in \mathbb{V}_h.
@f}
-
Here, $\mathbb{V}_h$ is an appropriate finite element space, and
$s_h(\cdot,\cdot)_{L^2(\Omega)}$ is some linear stabilization method
(possibly complemented with some ad-hoc shock-capturing technique, see for
from an analysis perspective, variational discretizations are conceived
to provide some notion of global (integral) stabiliy, meaning an
estimate of the form
-
@f{align*}
|\!|\!| \mathbf{u}_{h}(t) |\!|\!| \leq |\!|\!| \mathbf{u}_{h}(0) |\!|\!|
@f}
-
holds true, where $|\!|\!| \cdot |\!|\!| $ could represent the
$L^2(\Omega)$-norm or, more generally, some discrete (possibly mesh
dependent) energy-norm. Variational discretizations of hyperbolic
In this tutorial step we therefore depart from variational schemes. We will
present a completely algebraic formulation (with the flavor of a
collocation-type scheme) that preserves constraints pointwise, i.e.,
-
@f{align*}
\textbf{u}_h(\mathbf{x}_i,t) \in \mathcal{B}
\;\text{at every node}\;\mathbf{x}_i\;\text{of the mesh}.
@f}
-
Contrary to finite difference/volume schemes, the scheme implemented in
this step maximizes the use of finite element software infrastructure,
works in any mesh, in any space dimension, and is theoretically guaranteed
a scalar-valued shape function $\phi_i$.
With this notation at hand we can define the scheme as:
-
@f{align*}
m_i \frac{\mathbf{U}_i^{n+1} - \mathbf{U}_i^{n}}{\tau}
+ \sum_{j \in \mathcal{I}(i)} \mathbb{f}(\mathbf{U}_j^{n})\cdot
\mathbf{c}_{ij} - d_{ij} \mathbf{U}_j^{n} = \boldsymbol{0} \, ,
@f}
-
-Where
+where
- $m_i := \int_{\Omega} \phi_i \, \mathrm{d}\mathbf{x}$
- $\tau$ is the time step size
- $\mathbf{c}_{ij} := \int_{\Omega} \nabla\phi_j\phi_i \,
forward strategy for computing the solution $\textbf{U}^{n+1}$ at a new
time $t_{n+1} = t_n + \tau_n$ given a known state $\textbf{U}^{n}$ at time
$t_n$:
-
@f{align*}
&\textbf{for } i \in \mathcal{V} \\
&\ \ \ \ \{\mathbf{c}_{ij}\}_{j \in \mathcal{I}(i)} \leftarrow
dimensions, first-order polynomial space $\mathbb{Q}^1$, and that
$\mathbf{x}_i$ is an interior node (i.e. $\mathbf{x}_i$ is not on the boundary
of the domain ) then: $\{\textbf{U}_j^n\}_{j \in \mathcal{I}(i)}$ should contain
-nine state-vectors (i.e. all the states in the patch/macro element associated to
-the shape function $\phi_i$). This is one of the major differences with the
-usual cell-based loop where the gather functionality (encoded in
-FEValuesBase<dim, spacedim>.get_function_values() in the case of deal.ii) only
-collects values for the local cell (just a subset of the patch).
+nine state vector elements (i.e. all the states in the patch/macro element
+associated to the shape function $\phi_i$). This is one of the major
+differences with the usual cell-based loop where the gather functionality
+(encoded in FEValuesBase<dim, spacedim>.get_function_values() in the case
+of deal.ii) only collects values for the local cell (just a subset of the
+patch).
The actual implementation will deviate from above code in one key aspect:
the time-step size $\tau$ has to be chosen subject to a CFL condition
// indices <code>(i,l)</code> of the entry in order to retrieve its
// actual value. We should expect <code>gather</code> to be slightly
// more expensive than <code>gather_get_entry</code>. The use of
- // <code>gather</code> will be limited to the task of computing the
- // algebraic viscosity $d_{ij}$ in the particular case that when
+ // <code>gather</code> will be limited to the task of computing the
+ // algebraic viscosity $d_{ij}$ in the particular case that when
// both $i$ and $j$ lie at the boundary.
//
- // @note The reader should be aware that accessing an arbitrary
- // <code>(i,l)</code> entry of a matrix (say for instance Trilinos or PETSc
- // matrices) is in general unacceptably expensive. Here is where we might
- // want to keep an eye on complexity: we want this operation to have
- // constant complexity, which is the case of the current implementation
+ // @note The reader should be aware that accessing an arbitrary
+ // <code>(i,l)</code> entry of a matrix (say for instance Trilinos or PETSc
+ // matrices) is in general unacceptably expensive. Here is where we might
+ // want to keep an eye on complexity: we want this operation to have
+ // constant complexity, which is the case of the current implementation
// using deal.ii matrices.
template <typename T1, std::size_t k, typename T2, typename T3>
// <code>gather</code> (second interface): this second function
// signature having two input arguments will be used to gather the
- // state at a node <code>i</code> and return it as a
+ // state at a node <code>i</code> and return it as a
// <code>Tensor<1,problem_dimension></code> for our convenience.
template <typename T1, std::size_t k, typename T2>
}
// <code>scatter</code>: this function has three input arguments, the
- // first one is meant to be a "global object" (say a locally owned or
+ // first one is meant to be a "global object" (say a locally owned or
// locally relevant vector), the second argument which could be a
// <code>Tensor<1,problem_dimension></code>, and the last argument
// which represents a index of the global object. This function will be
// here $T$ denotes elements,
// $\text{supp}(\phi_i)$ the support of the shape function $\phi_i$,
// $F$ are faces of the element $T$, and $\mathbf{x}_{q,F}$
- // are quadrature points on such face. Note that this formula for
- // $\widehat{\boldsymbol{\nu}}_i$ is nothing else than some form of
- // weighted averaging. Other more sophisticated definitions for $\nu_i$
+ // are quadrature points on such face. Note that this formula for
+ // $\widehat{\boldsymbol{\nu}}_i$ is nothing else than some form of
+ // weighted averaging. Other more sophisticated definitions for $\nu_i$
// are possible but none of them have much influence in theory or practice.
template <int dim>
// In this section we describe the implementation of the class members of
// the <code>ProblemDescription</code> class. Most of the code here is
- // specific for compressible Euler's equations with an ideal gas law.
+ // specific for compressible Euler's equations with an ideal gas law.
// If we wanted to re-purpose Step-69 for a different conservation law
// (say for: instance the shallow water equation) most of the
// implementation of this class would have to change. But most of the other
// approximation for the intermediate pressure $p^*$, see for instance
// Equation (4.46), page 128 in @cite Toro2009.
//
- // The estimate returned by <code>lambda_max_two_rarefaction</code>
- // is guaranteed to be an upper bound, it is in general quite sharp, and
- // overall sufficient for our purposes. However, for some specific situations
- // (in particular when one of states is close to vacuum conditions) such
- // an estimate will be overly pessimistic. That's why we used a second
- // estimate to avoid this degeneracy that will be invoked by a call to the
- // function <code>lambda_max_expansion</code>. The most important function
- // here is <code>compute_lambda_max</code> which takes the minimum between
+ // The estimate returned by <code>lambda_max_two_rarefaction</code>
+ // is guaranteed to be an upper bound, it is in general quite sharp, and
+ // overall sufficient for our purposes. However, for some specific situations
+ // (in particular when one of states is close to vacuum conditions) such
+ // an estimate will be overly pessimistic. That's why we used a second
+ // estimate to avoid this degeneracy that will be invoked by a call to the
+ // function <code>lambda_max_expansion</code>. The most important function
+ // here is <code>compute_lambda_max</code> which takes the minimum between
// the estimates returned by <code>lambda_max_two_rarefaction</code> and
// <code>lambda_max_expansion</code>.
//
// primitive state $[\rho, u, p, a]$ and a given pressure $p^\ast$
// @cite GuermondPopov2016 Eqn. (3.7):
// @f{align*}
- // \lambda^- = u - a\,\sqrt{1 + \frac{\gamma+1}{2\gamma}
+ // \lambda^- = u - a\,\sqrt{1 + \frac{\gamma+1}{2\gamma}
// \left(\frac{p^\ast-p}{p}\right)_+}
// @f}
- // Here, the $(\cdot)_{+}$ denotes the positive part of the given
+ // Here, the $(\cdot)_{+}$ denotes the positive part of the given
// argument.
DEAL_II_ALWAYS_INLINE inline double
// Analougously @cite GuermondPopov2016 Eqn. (3.8):
// @f{align*}
- // \lambda^+ = u + a\,\sqrt{1 + \frac{\gamma+1}{2\gamma}
+ // \lambda^+ = u + a\,\sqrt{1 + \frac{\gamma+1}{2\gamma}
// \left(\frac{p^\ast-p}{p}\right)_+}
// @f}
}
} // namespace
- // The following is the main function that we are going to call in order to
+ // The following is the main function that we are going to call in order to
// compute $\lambda_{\text{max}} (\mathbf{U}_i^{n},\mathbf{U}_j^{n},
// \textbf{n}_{ij})$. We simply compute both maximal wavespeed estimates
// and return the minimum.
static constexpr auto gamma = ProblemDescription<dim>::gamma;
// The following lambda function translates a given primitive 1d state
- // (density $\rho$, velocity $u$, and pressure $p$) into a
- // conserved n-dimensional state (density $\rho$, momentum
+ // (density $\rho$, velocity $u$, and pressure $p$) into a
+ // conserved n-dimensional state (density $\rho$, momentum
// $\mathbf{m}$, and total energy $E$). Note that we
// <a href="https://en.cppreference.com/w/cpp/language/lambda">capture</a>
// the <code>this</code> pointer and thus access to
// \mathbf{U}_j^{n}, \textbf{n}_{ij}) = \lambda_{\text{max}}
// (\mathbf{U}_j^{n}, \mathbf{U}_i^{n}, \textbf{n}_{ji})$ do not
// necessarily hold true. The only mathematically safe solution for this
- // dilemma is to compute both of them $d_{ij}$ and $d_{ji}$ and
+ // dilemma is to compute both of them $d_{ij}$ and $d_{ji}$ and
// take the maximum.
//
// Overall, the computation of $d_{ij}$ is quite expensive. In
// order to save some computing time we exploit the fact that the viscosity
- // matrix has to be symmetric (as mentioned above): we only compute
+ // matrix has to be symmetric (as mentioned above): we only compute
// the upper-triangular entries of $d_{ij}$ and copy the
// corresponding entries to the lower-triangular counterpart.
//
{
const auto j = jt->column();
- // We only compute $d_{ij}$ if $j < i$ (upper triangular
+ // We only compute $d_{ij}$ if $j < i$ (upper triangular
// entries) and later copy the values over to $d_{ji}$.
if (j >= i)
continue;
double d = norm * lambda_max;
// If both support points happen to be at the boundary we
- // have to compute $d_{ji}$ as well and then take
+ // have to compute $d_{ji}$ as well and then take
// $max(d_{ij},d_{ji})$:
if (boundary_normal_map.count(i) != 0 &&
boundary_normal_map.count(j) != 0)
projects / dealii.git / commitdiff