Reference documentation for deal.II version GIT fac6d58a2d 2022-05-23 17:55:02+00:00
The step-69 tutorial program

This tutorial depends on step-33, step-40.

1. Introduction
2. The commented program
1. Results
2. The plain program
This program was contributed by Matthias Maier (Texas A&M University), and Ignacio Tomas (Sandia National Laboratories $$^{\!\dagger}$$).

$$^\dagger$$Sandia National Laboratories is a multimission laboratory managed and operated by National Technology & Engineering Solutions of Sandia, LLC, a wholly owned subsidiary of Honeywell International Inc., for the U.S. Department of Energy's National Nuclear Security Administration under contract DE-NA0003525. This document describes objective technical results and analysis. Any subjective views or opinions that might be expressed in the paper do not necessarily represent the views of the U.S. Department of Energy or the United States Government.

Note
This tutorial step implements a first-order accurate guaranteed maximum wavespeed method based on a first-order graph viscosity for solving Euler's equations of gas dynamics [GuermondPopov2016]. As such it is presented primarily for educational purposes. For actual research computations you might want to consider exploring a corresponding high-performance implementation of a second-order accurate scheme that uses convex limiting techniques, and strong stability-preserving (SSP) time integration, see [GuermondEtAl2018] (website).
If you use this program as a basis for your own work, please consider citing it in your list of references. The initial version of this work was contributed to the deal.II project by the authors listed in the following citation:

# Introduction

This tutorial presents a first-order scheme for solving compressible Euler's equations that is based on three ingredients: a collocation-type discretization of Euler's equations in the context of finite elements; a graph-viscosity stabilization based on a guaranteed 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.

It should be noted that first-order schemes in the context of hyperbolic conservation laws require prohibitively many degrees of freedom to resolve certain key features of the simulated fluid, and thus, typically only serve as elementary building blocks in higher-order schemes [GuermondEtAl2018]. However, we hope that the reader still finds the tutorial step to be a good starting point (in particular with respect to the programming techniques) before jumping into full research codes such as the second-order scheme discussed in [GuermondEtAl2018].

### Euler's equations of gas dynamics

The compressible Euler's equations of gas dynamics are written in conservative form as follows:

\begin{align} \mathbf{u}_t + \text{div} \, \mathbb{f}(\mathbf{u}) = \boldsymbol{0} , \end{align}

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 dimension. We say that $$\mathbf{u} \in \mathbb{R}^{d+2}$$ is the state and $$\mathbb{f}(\mathbf{u}) \in \mathbb{R}^{(d+2) \times d}$$ is the flux of the system. In the case of Euler's equations the state is given by $$\textbf{u} = [\rho, \textbf{m},E]^{\top}$$: where $$\rho \in \mathbb{R}^+$$ 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

\begin{align*} \mathbb{f}(\textbf{u}) = \begin{bmatrix} \textbf{m}^\top \\ \rho^{-1} \textbf{m} \otimes \textbf{m} + \mathbb{I} p\\ \tfrac{\textbf{m}^\top}{\rho} (E + p) \end{bmatrix}, \end{align*}

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 a closed-form equation of state. In this tutorial we limit the discussion to the class of polytropic ideal gases for which the pressure is given by

\begin{align*} p = p(\textbf{u}) := (\gamma -1) \Big(E - \tfrac{|\textbf{m}|^2}{2\,\rho} \Big), \end{align*}

where the factor $$\gamma \in (1,5/3]$$ denotes the ratio of specific heats.

#### Solution theory

Hyperbolic conservation laws, such as

\begin{align*} \mathbf{u}_t + \text{div} \, \mathbb{f}(\mathbf{u}) = \boldsymbol{0}, \end{align*}

pose a significant challenge with respect to solution theory. An evident observation is that rewriting the equation in variational form and testing with the solution itself does not lead to an energy estimate because the pairing $$\langle \text{div} \, \mathbb{f}(\mathbf{u}), \mathbf{u}\rangle$$ (understood as the $$L^2(\Omega)$$ inner product or duality pairing) is not guaranteed to be non-negative. Notions such as energy-stability or $$L^2(\Omega)$$-stability are (in general) meaningless in this context.

Historically, the most fruitful step taken in order to deepen the understanding of hyperbolic conservation laws was to assume that the 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

\begin{align} \mathbf{u}_t^{\epsilon} + \text{div} \, \mathbb{f}(\mathbf{u}^{\epsilon}) - {\epsilon} \Delta \mathbf{u}^{\epsilon} = 0. \end{align}

Such solutions, which are understood as the solution recovered in the zero-viscosity limit, are often referred to as viscosity solutions. (This is, because physically $$\epsilon$$ can be understood as related to the viscosity of the fluid, i.e., a quantity that indicates the amount of friction neighboring gas particles moving at different speeds exert on each other. The Euler equations themselves are derived under the assumption of no friction, but can physically be expected to describe the limiting case of vanishing friction or viscosity.) Global existence and uniqueness of such solutions is an 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

\begin{align} \mathcal{B} = \big\{ \textbf{u} = [\rho, \textbf{m},E]^{\top} \in \mathbb{R}^{d+2} \, \big | \ \rho > 0 \, , \ \ E - \tfrac{|\textbf{m}|^2}{2 \rho} > 0 \, , \ s(\mathbf{u}) \geq \min_{x \in \Omega} s(\mathbf{u}_0(\mathbf{x})) \big\}. \end{align}

Here, $$s(\mathbf{u})$$ denotes the specific entropy

\begin{align} s(\mathbf{u}) = \ln \Big(\frac{p(\mathbf{u})}{\rho^{\gamma}}\Big). \end{align}

We will refer to $$\mathcal{B}$$ as the invariant set of Euler's equations. In other words, a state $$\mathbf{u}(\mathbf{x},t)\in\mathcal{B}$$ obeys positivity of the density, positivity of the internal energy, and a local minimum principle on the specific entropy. This condition is a simplified version of a class of pointwise stability constraints satisfied by the exact (viscosity) solution. By pointwise we mean that the constraint has to be satisfied at every point of the domain, not just in an averaged (integral, or high order moments) sense.

In context of a numerical approximation, a violation of such a constraint has dire consequences: it almost surely leads to catastrophic failure of the numerical scheme, loss of hyperbolicity, and overall, loss of well-posedness of the (discrete) problem. It would also mean that we have computed something that can not be interpreted physically. (For example, what are we to make of a computed solution with a negative density?) In the following we will formulate a scheme that ensures that the discrete approximation of $$\mathbf{u}(\mathbf{x},t)$$ remains in $$\mathcal{B}$$.

#### Variational versus collocation-type discretizations

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:

\begin{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. \end{align*}

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 instance Chapter 5 of [GuermondErn2004] and references therein). Most time-dependent discretization approaches described in the deal.II tutorials are based on such a (semi-discrete) variational approach. Fundamentally, from an analysis perspective, variational discretizations are conceived to provide some notion of global (integral) stability, meaning an estimate of the form

\begin{align*} |\!|\!| \mathbf{u}_{h}(t) |\!|\!| \leq |\!|\!| \mathbf{u}_{h}(0) |\!|\!| \end{align*}

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 conservation laws have been very popular since the mid eighties, in particular combined with SUPG-type stabilization and/or upwinding techniques (see the early work of [Brooks1982] and [Johnson1986]). They 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 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 variational schemes in such extreme regimes is the lack of pointwise stability. This stems from the fact that a priori bounds on integrated quantities (e.g. integrals of moments) have in general no implications on pointwise properties of the solution. While some of these problems might be alleviated by the (perpetual) chase of the right shock capturing scheme, finite difference-like and finite volume schemes still have an edge in many regards.

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.,

\begin{align*} \textbf{u}_h(\mathbf{x}_i,t) \in \mathcal{B} \;\text{at every node}\;\mathbf{x}_i\;\text{of the mesh}. \end{align*}

Contrary to finite difference/volume schemes, the scheme implemented in this step maximizes the use of finite element software infrastructure, works on any mesh, in any space dimension, and is theoretically guaranteed to always work, all the time, no exception. This illustrates that deal.II can be used far beyond the context of variational schemes in Hilbert spaces and that a large number of classes, modules and namespaces from deal.II can be adapted for such a purpose.

### Description of the scheme

Let $$\mathbb{V}_h$$ be scalar-valued finite dimensional space spanned by a basis $$\{\phi_i\}_{i \in \mathcal{V}}$$ where: $$\phi_i:\Omega \rightarrow \mathbb{R}$$ and $$\mathcal{V}$$ is the set of all indices (nonnegative integers) identifying each scalar Degree of Freedom (DOF) in the mesh. Therefore a scalar finite element functional $$u_h \in \mathbb{V}_h$$ can be written as $$u_h = \sum_{i \in \mathcal{V}} U_i \phi_i$$ with $$U_i \in \mathbb{R}$$. We introduce the notation for vector-valued approximation spaces $$\pmb{\mathbb{V}}_h := \{\mathbb{V}_h\}^{d+2}$$. Let $$\mathbf{u}_h \in \pmb{\mathbb{V}}_h$$, then it can be written as $$\mathbf{u}_h = \sum_{i \in \mathcal{V}} \mathbf{U}_i \phi_i$$ where $$\mathbf{U}_i \in \mathbb{R}^{d+2}$$ and $$\phi_i$$ is a scalar-valued shape function.

Note
We purposely refrain from using vector-valued finite element spaces in our notation. Vector-valued finite element spaces are natural for variational formulations of PDE systems (e.g. Navier-Stokes). In such context, the interactions that have to be computed describe interactions between DOFs: with proper renumbering of the vector-valued DoFHandler (i.e. initialized with an FESystem) it is possible 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 better described as interactions between nodes (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 finite element spaces both in theory and/or practice.

We will use the usual Lagrange finite elements: let $$\{\mathbf{x}_i\}_{i \in \mathcal{V}}$$ denote the set of all support points (see this glossary entry), where $$\mathbf{x}_i \in \mathbb{R}^d$$. Then each index $$i \in \mathcal{V}$$ uniquely identifies a support point $$\mathbf{x}_i$$, as well as a scalar-valued shape function $$\phi_i$$. With this notation at hand we can define the (explicit time stepping) scheme as:

\begin{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} - \sum_{j \in \mathcal{I}(i)} d_{ij} \mathbf{U}_j^{n} = \boldsymbol{0} \, , \end{align*}

where

• $$m_i \dealcoloneq \int_{\Omega} \phi_i \, \mathrm{d}\mathbf{x}$$ is the lumped mass matrix
• $$\tau$$ is the time step size
• $$\mathbf{c}_{ij} \dealcoloneq \int_{\Omega} \nabla\phi_j\phi_i \, \mathrm{d}\mathbf{x}$$ (note that $$\mathbf{c}_{ij}\in \mathbb{R}^d$$) is a vector-valued matrix that was used to approximate the divergence of the flux in a weak sense.
• $$\mathcal{I}(i) \dealcoloneq \{j \in \mathcal{V} \ | \ \mathbf{c}_{ij} \not \equiv \boldsymbol{0}\} \cup \{i\}$$ is the adjacency list containing all degrees of freedom coupling to the index $$i$$. In other words $$\mathcal{I}(i)$$ contains all nonzero column indices for row index i. $$\mathcal{I}(i)$$ will also be called a "stencil".
• $$\mathbb{f}(\mathbf{U}_j^{n})$$ is the flux $$\mathbb{f}$$ of the hyperbolic system evaluated for the state $$\mathbf{U}_j^{n}$$ associated with support point $$\mathbf{x}_j$$.
• $$d_{ij} \dealcoloneq \max \{ \lambda_{\text{max}} (\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 graph viscosity. 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 an estimate of the maximal local wavespeed $$\lambda_{\text{max}}$$ that will be explained in detail in a moment.
• the diagonal entries of the viscosity matrix are defined as $$d_{ii} = - \sum_{j \in \mathcal{I}(i)\backslash \{i\}} d_{ij}$$.
• $$\textbf{n}_{ij} = \frac{\mathbf{c}_{ij}}{ \|\mathbf{c}_{ij}\| }$$ is a normalization of the $$\textbf{c}_{ij}$$ matrix that enters the approximate Riemann solver with which we compute an the approximations $$\lambda_{\text{max}}$$ on the local wavespeed. (This will be explained further down below).

The definition of $$\lambda_{\text{max}} (\mathbf{U},\mathbf{V}, \textbf{n})$$ is far from trivial and we will postpone the precise definition in order to focus first on some algorithmic and implementation questions. We note that

• $$m_i$$ and $$\mathbf{c}_{ij}$$ do not evolve in time (provided we keep the discretization fixed). It thus makes sense to assemble these matrices/vectors once in a so called offline computation and reuse them in every time step. They are part of what we are going to call off-line data.
• At every time step we have to evaluate $$\mathbb{f}(\mathbf{U}_j^{n})$$ and $$d_{ij} \dealcoloneq \max \{ \lambda_{\text{max}} (\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}\|$$, which will constitute the bulk of the computational cost.

Consider the following pseudo-code, illustrating a possible straight 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$$:

\begin{align*} &\textbf{for } i \in \mathcal{V} \\ &\ \ \ \ \{\mathbf{c}_{ij}\}_{j \in \mathcal{I}(i)} \leftarrow \mathtt{gather\_cij\_vectors} (\textbf{c}, \mathcal{I}(i)) \\ &\ \ \ \ \{\textbf{U}_j^n\}_{j \in \mathcal{I}(i)} \leftarrow \mathtt{gather\_state\_vectors} (\textbf{U}^n, \mathcal{I}(i)) \\ &\ \ \ \ \ \textbf{U}_i^{n+1} \leftarrow \mathbf{U}_i^{n} \\ &\ \ \ \ \textbf{for } j \in \mathcal{I}(i)\backslash\{i\} \\ &\ \ \ \ \ \ \ \ \texttt{compute } d_{ij} \\ &\ \ \ \ \ \ \ \ \texttt{compute } \mathbb{f}(\mathbf{U}_j^{n}) \\ &\ \ \ \ \ \ \ \ \textbf{U}_i^{n+1} \leftarrow \textbf{U}_i^{n+1} - \frac{\tau_n}{m_i} \mathbb{f}(\mathbf{U}_j^{n})\cdot \mathbf{c}_{ij} + d_{ij} \mathbf{U}_j^{n} \\ &\ \ \ \ \textbf{end} \\ &\ \ \ \ \mathtt{scatter\_updated\_state} (\textbf{U}_i^{n+1}) \\ &\textbf{end} \end{align*}

We note here that:

• This "assembly" does not require any form of quadrature or cell-loops.
• Here $$\textbf{c}$$ and $$\textbf{U}^n$$ are a global matrix and a global vector containing all the vectors $$\mathbf{c}_{ij}$$ and all the states $$\mathbf{U}_j^n$$ respectively.
• $$\mathtt{gather\_cij\_vectors}$$, $$\mathtt{gather\_state\_vectors}$$, and $$\mathtt{scatter\_updated\_state}$$ are hypothetical implementations that either collect (from) or write (into) global matrices and vectors.
• If we assume a Cartesian mesh in two space 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 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

\begin{align*} \tau_n = c_{\text{cfl}}\,\min_{ i\in\mathcal{V}}\left(\frac{m_i}{-2\,d_{ii}^{n}}\right), \end{align*}

where $$0<c_{\text{cfl}}\le1$$ is a chosen constant. This will require to compute all $$d_{ij}$$ in a separate step prior to actually performing above 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. 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 edge-based or graph-based finite element schemes (see for instance [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 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 information that spans more than one cell (superposition does not hold: adding contributions from separate cells does not lead to the right result).

### Stable boundary conditions and conservation properties.

In the example considered in this tutorial step we use three different types of boundary conditions: essential-like boundary conditions (we prescribe a state at the left boundary of our domain), outflow boundary conditions (also called "do-nothing" boundary conditions) at the right boundary of the domain, and "reflecting" boundary conditions $$\mathbf{m} \cdot \boldsymbol{\nu} = 0$$ (also called "slip" boundary conditions) at the top, bottom, and surface of the obstacle. We will not discuss much about essential and "do-nothing" boundary conditions since their implementation is relatively easy and the reader will be able to pick-up the implementation directly from the (documented) source code. In this portion of the introduction we will focus only on the "reflecting" boundary conditions which are somewhat more tricky.

Note
At the time of this writing (early 2020) it is not unreasonable to say that both analysis and implementation of stable boundary conditions for hyperbolic systems of conservation laws is an open issue. For the case of variational formulations, stable boundary conditions are those leading to a well-posed (coercive) bilinear form. But for general hyperbolic systems of conservation laws (and for the algebraic formulation used in this tutorial) coercivity has no applicability and/or meaning as a notion of stability. In this tutorial step we will use preservation of the invariant set as our main notion of stability which (at the very least) guarantees well-posedness of the discrete problem.

For the case of the reflecting boundary conditions we will proceed as follows:

• For every time step advance in time satisfying no boundary condition at all.
• Let $$\partial\Omega^r$$ be the portion of the boundary where we want to enforce reflecting boundary conditions. At the end of the time step we enforce reflecting boundary conditions strongly in a post-processing step where we execute the projection

\begin{align*} \mathbf{m}_i \dealcoloneq \mathbf{m}_i - (\widehat{\boldsymbol{\nu}}_i \cdot \mathbf{m}_i) \widehat{\boldsymbol{\nu}}_i \ \ \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)} \end{align*}

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. No doubt, when solving parabolic, or elliptic equations, we typically enforce essential (Dirichlet-like) boundary conditions by making them part of the approximation space $$\mathbb{V}$$, and treat natural (e.g. Neumann) boundary conditions as part of the variational formulation. We also know that explicit treatment of boundary conditions (in the context of parabolic PDEs) almost surely leads to catastrophic consequences. However, in the context of nonlinear hyperbolic equations we have that:

• It is relatively easy to prove that (for the case of reflecting boundary conditions) explicit treatment of boundary conditions is not only conservative but also guarantees preservation of the property $$\mathbf{U}_i \in \mathcal{B}$$ for all $$i \in \mathcal{V}$$ (well-posedness). This is perhaps the most important reason to use explicit enforcement of boundary conditions.
• To the best of our knowledge: we are not aware of any mathematical result proving that it is possible to guarantee the property $$\mathbf{U}_i \in \mathcal{B}$$ for all $$i \in \mathcal{V}$$ when using either direct enforcement of boundary conditions into the approximation space, or weak enforcement using the Nitsche penalty method (which is for example widely used in discontinuous Galerkin schemes). In addition, some of these traditional ideas lead to quite restrictive time step constraints.
• There is enough numerical evidence suggesting that explicit treatment of Dirichlet-like boundary conditions is stable under CFL conditions and does not introduce any loss in accuracy.

If $$\mathbf{u}_t + \text{div} \, \mathbb{f}(\mathbf{u}) = \boldsymbol{0}$$ represents Euler's equation with reflecting boundary conditions on the entirety of the boundary (i.e. $$\partial\Omega^r \equiv \partial\Omega$$) and we integrate in space and time $$\int_{\Omega}\int_{t_1}^{t_2}$$ we would obtain

\begin{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} \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 = \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)} \end{align*}

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 particular, if we use the projection $$\boldsymbol{(1)}$$ in the entirety of the domain the following discrete mass-balance can be guaranteed:

\begin{align*} \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)} \end{align*}

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 nodal normal $$\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 [GuermondEtAl2018].

# The commented program

### Include files

The set of include files is quite standard. The most intriguing part is the fact that we will rely solely on deal.II data structures for MPI parallelization, in particular parallel::distributed::Triangulation and LinearAlgebra::distributed::Vector included through distributed/tria.h and lac/la_parallel_vector.h. Instead of a Trilinos, or PETSc specific matrix class, we will use a non-distributed SparseMatrix (lac/sparse_matrix.h) to store the local part of the $$\mathbf{c}_{ij}$$, $$\mathbf{n}_{ij}$$ and $$d_{ij}$$ matrices.

In addition to above deal.II specific includes, we also include four boost headers. The first two are for binary archives that we will use for implementing a check-pointing and restart mechanism.

#include <boost/archive/binary_iarchive.hpp>
#include <boost/archive/binary_oarchive.hpp>

The last two header files are for creating custom iterator ranges over integer intervals.

#include <boost/range/iterator_range.hpp>

For std::isnan, std::isinf, std::ifstream, std::async, and std::future

#include <cmath>
#include <fstream>
#include <future>

### Class template declarations

We begin our actual implementation by declaring all classes with their data structures and methods upfront. In contrast to previous example steps we use a more fine-grained encapsulation of concepts, data structures, and parameters into individual classes. A single class thus usually centers around either a single data structure (such as the Triangulation) in the Discretization class, or a single method (such as the make_one_step() function of the TimeStepping class). We typically declare parameter variables and scratch data object private and make methods and data structures used by other classes public.

Note
A cleaner approach would be to guard access to all data structures by getter/setter functions. For the sake of brevity, we refrain from that approach, though.

We also note that the vast majority of classes is derived from ParameterAcceptor. This facilitates the population of all the global parameters into a single (global) ParameterHandler. More explanations about the use of inheritance from ParameterAcceptor as a global subscription mechanism can be found in step-60.

namespace Step69
{
using namespace dealii;

We start with defining a number of types::boundary_id constants used throughout the tutorial step. This allows us to refer to boundary types by a mnemonic (such as do_nothing) rather than a numerical value.

namespace Boundaries
{
constexpr types::boundary_id do_nothing = 0;
constexpr types::boundary_id free_slip = 1;
constexpr types::boundary_id dirichlet = 2;
} // namespace Boundaries
unsigned int boundary_id
Definition: types.h:129

#### The Discretization class

The class Discretization contains all data structures concerning the mesh (triangulation) and discretization (mapping, finite element, quadrature) of the problem. As mentioned, we use the ParameterAcceptor class to automatically populate problem-specific parameters, such as the geometry information (length, etc.) or the refinement level (refinement) from a parameter file. This requires us to split the initialization of data structures into two functions: We initialize everything that does not depend on parameters in the constructor, and defer the creation of the mesh to the setup() method that can be called once all parameters are read in via ParameterAcceptor::initialize().

template <int dim>
class Discretization : public ParameterAcceptor
{
public:
Discretization(const MPI_Comm mpi_communicator,
TimerOutput & computing_timer,
const std::string &subsection = "Discretization");
void setup();
const MPI_Comm mpi_communicator;
const MappingQ<dim> mapping;
const FE_Q<dim> finite_element;
unsigned int refinement;
private:
TimerOutput &computing_timer;
double length;
double height;
double disk_position;
double disk_diameter;
};
Definition: fe_q.h:549
const ::parallel::distributed::Triangulation< dim, spacedim > * triangulation

#### The OfflineData class

The class OfflineData contains pretty much all components of the discretization that do not evolve in time, in particular, the DoFHandler, SparsityPattern, boundary maps, the lumped mass matrix, $$\mathbf{c}_{ij}$$ and $$\mathbf{n}_{ij}$$ matrices. Here, the term offline refers to the fact that all the class members of OfflineData have well-defined values independent of the current time step. This means that they can be initialized ahead of time (at time step zero) and are not meant to be modified at any later time step. For instance, the sparsity pattern should not change as we advance in time (we are not doing any form of adaptivity in space). Similarly, the entries of the lumped mass matrix should not be modified as we advance in time either.

We also compute and store a boundary_normal_map that contains a map from a global index of type types::global_dof_index of a boundary degree of freedom to a tuple consisting of a normal vector, the boundary id, and the position associated with the degree of freedom. We have to compute and store this geometric information in this class because we won't have access to geometric (or cell-based) information later on in the algebraic loops over the sparsity pattern.

Note
Even though this class currently does not have any parameters that could be read in from a parameter file we nevertheless derive from ParameterAcceptor and follow the same idiom of providing a setup() (and assemble()) method as for the class Discretization.
template <int dim>
class OfflineData : public ParameterAcceptor
{
public:
using BoundaryNormalMap =
std::tuple<Tensor<1, dim>, types::boundary_id, Point<dim>>>;
OfflineData(const MPI_Comm mpi_communicator,
TimerOutput & computing_timer,
const Discretization<dim> &discretization,
const std::string & subsection = "OfflineData");
void setup();
void assemble();
DoFHandler<dim> dof_handler;
std::shared_ptr<const Utilities::MPI::Partitioner> partitioner;
unsigned int n_locally_owned;
unsigned int n_locally_relevant;
SparsityPattern sparsity_pattern;
BoundaryNormalMap boundary_normal_map;
SparseMatrix<double> lumped_mass_matrix;
std::array<SparseMatrix<double>, dim> cij_matrix;
std::array<SparseMatrix<double>, dim> nij_matrix;
SparseMatrix<double> norm_matrix;
private:
const MPI_Comm mpi_communicator;
TimerOutput & computing_timer;
};
Definition: point.h:111
void assemble(const MeshWorker::DoFInfoBox< dim, DOFINFO > &dinfo, A *assembler)
Definition: loop.h:71
unsigned int global_dof_index
Definition: types.h:76

#### The ProblemDescription class

The member functions of this class are utility functions and data structures specific to Euler's equations:

• The type alias state_type is used for the states $$\mathbf{U}_i^n$$
• The type alias flux_type is used for the fluxes $$\mathbb{f}(\mathbf{U}_j^n)$$.
• The momentum function extracts $$\textbf{m}$$ out of the state vector $$[\rho,\textbf{m},E]$$ and stores it in a Tensor<1, dim>.
• The internal_energy function computes $$E - \frac{|\textbf{m}|^2}{2\rho}$$ from a given state vector $$[\rho,\textbf{m},E]$$.

The purpose of the class members component_names, pressure, and speed_of_sound is evident from their names. We also provide a function compute_lambda_max(), that computes the wave speed estimate mentioned above, $$\lambda_{max}(\mathbf{U},\mathbf{V},\mathbf{n})$$, which is used in the computation of the $$d_{ij}$$ matrix.

Note
The DEAL_II_ALWAYS_INLINE macro expands to a (compiler specific) pragma that ensures that the corresponding function defined in this class is always inlined, i.e., the function body is put in place for every invocation of the function, and no call (and code indirection) is generated. This is stronger than the inline keyword, which is more or less a (mild) suggestion to the compiler that the programmer thinks it would be beneficial to inline the function. DEAL_II_ALWAYS_INLINE should only be used rarely and with caution in situations such as this one, where we actually know (due to benchmarking) that inlining the function in question improves performance.

Finally, we observe that this is the only class in this tutorial step that is tied to a particular "physics" or "hyperbolic conservation law" (in this case Euler's equations). All the other classes are primarily "discretization" classes, very much agnostic of the particular physics being solved.

template <int dim>
class ProblemDescription
{
public:
static constexpr unsigned int problem_dimension = 2 + dim;
using state_type = Tensor<1, problem_dimension>;
const static std::array<std::string, problem_dimension> component_names;
static constexpr double gamma = 7. / 5.;
momentum(const state_type &U);
static DEAL_II_ALWAYS_INLINE inline double
internal_energy(const state_type &U);
static DEAL_II_ALWAYS_INLINE inline double pressure(const state_type &U);
static DEAL_II_ALWAYS_INLINE inline double
speed_of_sound(const state_type &U);
static DEAL_II_ALWAYS_INLINE inline flux_type flux(const state_type &U);
static DEAL_II_ALWAYS_INLINE inline double
compute_lambda_max(const state_type & U_i,
const state_type & U_j,
const Tensor<1, dim> &n_ij);
};
Definition: tensor.h:503
#define DEAL_II_ALWAYS_INLINE
Definition: config.h:102
static const char U
long double gamma(const unsigned int n)

#### The InitialValues class

The class InitialValues's only public data attribute is a std::function initial_state that computes the initial state of a given point and time. This function is used for populating the initial flow field as well as setting Dirichlet boundary conditions (at inflow boundaries) explicitly in every time step.

For the purpose of this example step we simply implement a homogeneous uniform flow field for which the direction and a 1D primitive state (density, velocity, pressure) are read from the parameter file.

It would be desirable to initialize the class in a single shot: initialize/set the parameters and define the class members that depend on these default parameters. However, since we do not know the actual values for the parameters, this would be sort of meaningless and unsafe in general (we would like to have mechanisms to check the consistency of the input parameters). Instead of defining another setup() method to be called (by-hand) after the call to ParameterAcceptor::initialize() we provide an "implementation" for the class member parse_parameters_call_back() which is automatically called when invoking ParameterAcceptor::initialize() for every class that inherits from ParameterAceptor.

template <int dim>
class InitialValues : public ParameterAcceptor
{
public:
using state_type = typename ProblemDescription<dim>::state_type;
InitialValues(const std::string &subsection = "InitialValues");
std::function<state_type(const Point<dim> &point, double t)> initial_state;
private:
Point< spacedim > point(const gp_Pnt &p, const double tolerance=1e-10)
Definition: utilities.cc:185

We declare a private callback function that will be wired up to the ParameterAcceptor::parse_parameters_call_back signal.

void parse_parameters_callback();
Tensor<1, dim> initial_direction;
Tensor<1, 3> initial_1d_state;
};

#### The TimeStepping class

With the OfflineData and ProblemDescription classes at hand we can now implement the explicit time-stepping scheme that was introduced in the discussion above. The main method of the TimeStepping class is make_one_step(vector_type &U, double t) that takes a reference to a state vector U and a time point t (as input arguments) computes the updated solution, stores it in the vector temp, swaps its contents with the vector U, and returns the chosen step-size $$\tau$$.

The other important method is prepare() which primarily sets the proper partition and sparsity pattern for the temporary vector temp and the matrix dij_matrix respectively.

template <int dim>
{
public:
static constexpr unsigned int problem_dimension =
ProblemDescription<dim>::problem_dimension;
using state_type = typename ProblemDescription<dim>::state_type;
using flux_type = typename ProblemDescription<dim>::flux_type;
using vector_type =
std::array<LinearAlgebra::distributed::Vector<double>, problem_dimension>;
TimeStepping(const MPI_Comm mpi_communicator,
TimerOutput & computing_timer,
const OfflineData<dim> & offline_data,
const InitialValues<dim> &initial_values,
const std::string & subsection = "TimeStepping");
void prepare();
double make_one_step(vector_type &U, double t);
private:
const MPI_Comm mpi_communicator;
TimerOutput & computing_timer;
vector_type temporary_vector;
double cfl_update;
};

#### The SchlierenPostprocessor class

At its core, the Schlieren class implements the class member compute_schlieren(). The main purpose of this class member is to compute an auxiliary finite element field schlieren, that is defined at each node by

$\text{schlieren}[i] = e^{\beta \frac{ |\nabla r_i| - \min_j |\nabla r_j| }{\max_j |\nabla r_j| - \min_j |\nabla r_j| } },$

where $$r$$ can in principle be any scalar quantity. In practice though, the density is a natural candidate, viz. $$r \dealcoloneq \rho$$. Schlieren postprocessing is a standard method for enhancing the contrast of a visualization inspired by actual experimental X-ray and shadowgraphy techniques of visualization. (See step-67 for another example where we create a Schlieren plot.)

template <int dim>
class SchlierenPostprocessor : public ParameterAcceptor
{
public:
static constexpr unsigned int problem_dimension =
ProblemDescription<dim>::problem_dimension;
using state_type = typename ProblemDescription<dim>::state_type;
using vector_type =
std::array<LinearAlgebra::distributed::Vector<double>, problem_dimension>;
SchlierenPostprocessor(
const MPI_Comm mpi_communicator,
TimerOutput & computing_timer,
const OfflineData<dim> &offline_data,
const std::string & subsection = "SchlierenPostprocessor");
void prepare();
void compute_schlieren(const vector_type &U);
private:
const MPI_Comm mpi_communicator;
TimerOutput & computing_timer;
unsigned int schlieren_index;
double schlieren_beta;
};

#### The MainLoop class

Now, all that is left to do is to chain the methods implemented in the TimeStepping, InitialValues, and SchlierenPostprocessor classes together. We do this in a separate class MainLoop that contains an object of every class and again reads in a number of parameters with the help of the ParameterAcceptor class.

template <int dim>
class MainLoop : public ParameterAcceptor
{
public:
using vector_type = typename TimeStepping<dim>::vector_type;
MainLoop(const MPI_Comm mpi_communnicator);
void run();
private:
vector_type interpolate_initial_values(const double t = 0);
void checkpoint(const vector_type &U,
const std::string &name,
double t,
unsigned int cycle);
void output(const vector_type &U,
const std::string &name,
double t,
unsigned int cycle);
const MPI_Comm mpi_communicator;
std::ostringstream timer_output;
TimerOutput computing_timer;
std::string base_name;
double t_final;
double output_granularity;
bool asynchronous_writeback;
bool resume;
Discretization<dim> discretization;
OfflineData<dim> offline_data;
InitialValues<dim> initial_values;
TimeStepping<dim> time_stepping;
SchlierenPostprocessor<dim> schlieren_postprocessor;
vector_type output_vector;
};
void run(const Iterator &begin, const typename identity< Iterator >::type &end, Worker worker, Copier copier, const ScratchData &sample_scratch_data, const CopyData &sample_copy_data, const unsigned int queue_length, const unsigned int chunk_size)
Definition: work_stream.h:474

### Implementation

#### Grid generation, setup of data structures

The first major task at hand is the typical triplet of grid generation, setup of data structures, and assembly. A notable novelty in this example step is the use of the ParameterAcceptor class that we use to populate parameter values: we first initialize the ParameterAcceptor class by calling its constructor with a string subsection denoting the correct subsection in the parameter file. Then, in the constructor body every parameter value is initialized to a sensible default value and registered with the ParameterAcceptor class with a call to ParameterAcceptor::add_parameter().

template <int dim>
Discretization<dim>::Discretization(const MPI_Comm mpi_communicator,
TimerOutput & computing_timer,
const std::string &subsection)
: ParameterAcceptor(subsection)
, mpi_communicator(mpi_communicator)
, triangulation(mpi_communicator)
, mapping(1)
, finite_element(1)
, computing_timer(computing_timer)
{
length = 4.;
add_parameter("length", length, "Length of computational domain");
height = 2.;
add_parameter("height", height, "Height of computational domain");
disk_position = 0.6;
disk_position,
"x position of immersed disk center point");
disk_diameter = 0.5;
disk_diameter,
"Diameter of immersed disk");
refinement = 5;
refinement,
"Number of refinement steps of the geometry");
}

Note that in the previous constructor we only passed the MPI communicator to the triangulation but we still have not initialized the underlying geometry/mesh. As mentioned earlier, we have to postpone this task to the setup() function that gets called after the ParameterAcceptor::initialize() function has populated all parameter variables with the final values read from the parameter file.

The setup() function is the last class member that has to be implemented. It creates the actual triangulation that is a benchmark configuration consisting of a channel with a disk obstacle, see [GuermondEtAl2018]. We construct the geometry by modifying the mesh generated by GridGenerator::hyper_cube_with_cylindrical_hole(). We refer to step-49, step-53, and step-54 for an overview how to create advanced meshes. We first create 4 temporary (non distributed) coarse triangulations that we stitch together with the GridGenerator::merge_triangulation() function. We center the disk at $$(0,0)$$ with a diameter of disk_diameter. The lower left corner of the channel has coordinates (-disk_position, -height/2) and the upper right corner has (length-disk_position, height/2).

template <int dim>
void Discretization<dim>::setup()
{
TimerOutput::Scope scope(computing_timer, "discretization - setup");
triangulation.clear();
Triangulation<dim> tria1, tria2, tria3, tria4, tria5, tria6;
tria1, disk_diameter / 2., disk_diameter, 0.5, 1, false);
tria2,
{2, 1},
Point<2>(-disk_diameter, disk_diameter),
Point<2>(disk_diameter, height / 2.));
tria3,
{2, 1},
Point<2>(-disk_diameter, -disk_diameter),
Point<2>(disk_diameter, -height / 2.));
tria4,
{6, 2},
Point<2>(disk_diameter, -disk_diameter),
Point<2>(length - disk_position, disk_diameter));
tria5,
{6, 1},
Point<2>(disk_diameter, disk_diameter),
Point<2>(length - disk_position, height / 2.));
tria6,
{6, 1},
Point<2>(disk_diameter, -height / 2.),
Point<2>(length - disk_position, -disk_diameter));
{&tria1, &tria2, &tria3, &tria4, &tria5, &tria6},
1.e-12,
true);
void merge_triangulations(const Triangulation< dim, spacedim > &triangulation_1, const Triangulation< dim, spacedim > &triangulation_2, Triangulation< dim, spacedim > &result, const double duplicated_vertex_tolerance=1.0e-12, const bool copy_manifold_ids=false)
void subdivided_hyper_rectangle(Triangulation< dim, spacedim > &tria, const std::vector< unsigned int > &repetitions, const Point< dim > &p1, const Point< dim > &p2, const bool colorize=false)
void hyper_cube_with_cylindrical_hole(Triangulation< dim > &triangulation, const double inner_radius=.25, const double outer_radius=.5, const double L=.5, const unsigned int repetitions=1, const bool colorize=false)

We have to fix up the left edge that is currently located at $$x=-$$disk_diameter and has to be shifted to $$x=-$$disk_position. As a last step the boundary has to be colorized with Boundaries::do_nothing on the right, dirichlet on the left and free_slip on the upper and lower outer boundaries and the obstacle.

for (const auto &cell : triangulation.active_cell_iterators())
{
for (const auto v : cell->vertex_indices())
{
if (cell->vertex(v)[0] <= -disk_diameter + 1.e-6)
cell->vertex(v)[0] = -disk_position;
}
}
for (const auto &cell : triangulation.active_cell_iterators())
{
for (const auto f : cell->face_indices())
{
const auto face = cell->face(f);
if (face->at_boundary())
{
const auto center = face->center();
if (center[0] > length - disk_position - 1.e-6)
face->set_boundary_id(Boundaries::do_nothing);
else if (center[0] < -disk_position + 1.e-6)
face->set_boundary_id(Boundaries::dirichlet);
else
face->set_boundary_id(Boundaries::free_slip);
}
}
}
}
Point< 3 > center

#### Assembly of offline matrices

Not much is done in the constructor of OfflineData other than initializing the corresponding class members in the initialization list.

template <int dim>
OfflineData<dim>::OfflineData(const MPI_Comm mpi_communicator,
TimerOutput & computing_timer,
const Discretization<dim> &discretization,
const std::string & subsection)
: ParameterAcceptor(subsection)
, dof_handler(discretization.triangulation)
, mpi_communicator(mpi_communicator)
, computing_timer(computing_timer)
, discretization(&discretization)
{}

Now we can initialize the DoFHandler, extract the IndexSet objects for locally owned and locally relevant DOFs, and initialize a Utilities::MPI::Partitioner object that is needed for distributed vectors.

template <int dim>
void OfflineData<dim>::setup()
{
IndexSet locally_owned;
IndexSet locally_relevant;
{
TimerOutput::Scope scope(computing_timer,
"offline_data - distribute dofs");
dof_handler.distribute_dofs(discretization->finite_element);
locally_owned = dof_handler.locally_owned_dofs();
n_locally_owned = locally_owned.n_elements();
locally_relevant = DoFTools::extract_locally_relevant_dofs(dof_handler);
n_locally_relevant = locally_relevant.n_elements();
partitioner =
std::make_shared<Utilities::MPI::Partitioner>(locally_owned,
locally_relevant,
mpi_communicator);
}
size_type n_elements() const
Definition: index_set.h:1834
IndexSet extract_locally_relevant_dofs(const DoFHandler< dim, spacedim > &dof_handler)
Definition: dof_tools.cc:1144

#### Translation to local index ranges

We are now in a position to create the sparsity pattern for our matrices. There are quite a few peculiarities that need a detailed explanation. We avoid using a distributed matrix class (as for example provided by Trilinos or PETSc) and instead rely on deal.II's own SparseMatrix object to store the local part of all matrices. This design decision is motivated by the fact that (a) we actually never perform a matrix-vector multiplication, and (b) we can always assemble the local part of a matrix exclusively on a given MPI rank. Instead, we will compute nonlinear updates while iterating over (the local part) of a connectivity stencil; a task for which deal.II's own SparsityPattern is specifically optimized for.

This design consideration has a caveat, though. What makes the deal.II SparseMatrix class fast is the compressed row storage (CSR) used in the SparsityPattern (see Sparsity patterns). This, unfortunately, does not play nicely with a global distributed index range because a sparsity pattern with CSR cannot contain "holes" in the index range. The distributed matrices offered by deal.II avoid this by translating from a global index range into a contiguous local index range. But this is precisely the type of index manipulation we want to avoid in our iteration over the stencil because it creates a measurable overhead.

The Utilities::MPI::Partitioner class already implements the translation from a global index range to a contiguous local (per MPI rank) index range: we don't have to reinvent the wheel. We just need to use that translation capability (once and only once) in order to create a "local" sparsity pattern for the contiguous index range $$[0,$$n_locally_relevant $$)$$. That capability can be invoked by the Utilities::MPI::Partitioner::global_to_local() function. Once the sparsity pattern is created using local indices, all that is left to do is to ensure that (when implementing our scatter and gather auxiliary functions) we always access elements of a distributed vector by a call to LinearAlgebra::distributed::Vector::local_element(). This way we avoid index translations altogether and operate exclusively with local indices.

{
computing_timer,
"offline_data - create sparsity pattern and set up matrices");

We have to create the "local" sparsity pattern by hand. We therefore loop over all locally owned and ghosted cells (see GlossArtificialCell) and extract the (global) dof_indices associated with the cell DOFs and renumber them using partitioner->global_to_local(index).

Note
In the case of a locally owned dof, such renumbering consist of applying a shift (i.e. we subtract an offset) such that now they will become a number in the integer interval $$[0,$$n_locally_owned $$)$$. However, in the case of a ghosted dof (i.e. not locally owned) the situation is quite different, since the global indices associated with ghosted DOFs will not be (in general) a contiguous set of integers.
DynamicSparsityPattern dsp(n_locally_relevant, n_locally_relevant);
const auto dofs_per_cell =
discretization->finite_element.n_dofs_per_cell();
std::vector<types::global_dof_index> dof_indices(dofs_per_cell);
for (const auto &cell : dof_handler.active_cell_iterators())
{
if (cell->is_artificial())
continue;
/* We transform the set of global dof indices on the cell to the
* corresponding "local" index range on the MPI process: */
cell->get_dof_indices(dof_indices);
std::transform(dof_indices.begin(),
dof_indices.end(),
dof_indices.begin(),
return partitioner->global_to_local(index);
});
/* And simply add, for each dof, a coupling to all other "local"
* dofs on the cell: */
for (const auto dof : dof_indices)
}
sparsity_pattern.copy_from(dsp);
lumped_mass_matrix.reinit(sparsity_pattern);
norm_matrix.reinit(sparsity_pattern);
for (auto &matrix : cij_matrix)
matrix.reinit(sparsity_pattern);
for (auto &matrix : nij_matrix)
matrix.reinit(sparsity_pattern);
}
}
@ matrix
Contents is actually a matrix.
void transform(const InputIterator &begin_in, const InputIterator &end_in, OutputIterator out, const Predicate &predicate, const unsigned int grainsize)
Definition: parallel.h:135

This concludes the setup of the DoFHandler and SparseMatrix objects. Next, we have to assemble various matrices. We define a number of helper functions and data structures in an anonymous namespace.

namespace
{

CopyData class that will be used to assemble the offline data matrices using WorkStream. It acts as a container: it is just a struct where WorkStream stores the local cell contributions. Note that it also contains a class member local_boundary_normal_map used to store the local contributions required to compute the normals at the boundary.

template <int dim>
struct CopyData
{
bool is_artificial;
std::vector<types::global_dof_index> local_dof_indices;
typename OfflineData<dim>::BoundaryNormalMap local_boundary_normal_map;
FullMatrix<double> cell_lumped_mass_matrix;
std::array<FullMatrix<double>, dim> cell_cij_matrix;
};

Next we introduce a number of helper functions that are all concerned about reading and writing matrix and vector entries. They are mainly motivated by providing slightly more efficient code and syntactic sugar for otherwise somewhat tedious code.

The first function we introduce, get_entry(), will be used to read the value stored at the entry pointed by a SparsityPattern iterator it of matrix. The function works around a small deficiency in the SparseMatrix interface: The SparsityPattern is concerned with all index operations of the sparse matrix stored in CRS format. As such the iterator already knows the global index of the corresponding matrix entry in the low-level vector stored in the SparseMatrix object. Due to the lack of an interface in the SparseMatrix for accessing the element directly with a SparsityPattern iterator, we unfortunately have to create a temporary SparseMatrix iterator. We simply hide this in the get_entry() function.

template <typename IteratorType>
get_entry(const SparseMatrix<double> &matrix, const IteratorType &it)
{
const SparseMatrix<double>::const_iterator matrix_iterator(
&matrix, it->global_index());
return matrix_iterator->value();
}

The set_entry() helper is the inverse operation of get_value(): Given an iterator and a value, it sets the entry pointed to by the iterator in the matrix.

template <typename IteratorType>
const IteratorType & it,
{
it->global_index());
matrix_iterator->value() = value;
}

gather_get_entry(): we note that $$\mathbf{c}_{ij} \in \mathbb{R}^d$$. If $$d=2$$ then $$\mathbf{c}_{ij} = [\mathbf{c}_{ij}^1,\mathbf{c}_{ij}^2]^\top$$. Which basically implies that we need one matrix per space dimension to store the $$\mathbf{c}_{ij}$$ vectors. Similar observation follows for the matrix $$\mathbf{n}_{ij}$$. The purpose of gather_get_entry() is to retrieve those entries and store them into a Tensor<1, dim> for our convenience.

template <std::size_t k, typename IteratorType>
gather_get_entry(const std::array<SparseMatrix<double>, k> &c_ij,
const IteratorType it)
{
Tensor<1, k> result;
for (unsigned int j = 0; j < k; ++j)
result[j] = get_entry(c_ij[j], it);
return result;
}

gather() (first interface): this first function signature, having three input arguments, will be used to retrieve the individual components (i,l) of a matrix. The functionality of gather_get_entry() and gather() is very much the same, but their context is different: the function gather() does not rely on an iterator (that actually knows the value pointed to) but rather on the indices (i,l) of the entry in order to retrieve its actual value. We should expect gather() to be slightly more expensive than gather_get_entry(). The use of gather() 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 (i,l) 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 <std::size_t k>
gather(const std::array<SparseMatrix<double>, k> &n_ij,
const unsigned int i,
const unsigned int j)
{
Tensor<1, k> result;
for (unsigned int l = 0; l < k; ++l)
result[l] = n_ij[l](i, j);
return result;
}
Tensor< 2, dim, Number > l(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
void gather(VectorizedArray< Number, width > &out, const std::array< Number *, width > &ptrs, const unsigned int offset)

gather() (second interface): this second function signature having two input arguments will be used to gather the state at a node i and return it as a Tensor<1,problem_dimension> for our convenience.

template <std::size_t k>
const unsigned int i)
{
Tensor<1, k> result;
for (unsigned int j = 0; j < k; ++j)
result[j] = U[j].local_element(i);
return result;
}

scatter(): this function has three input arguments, the first one is meant to be a "global object" (say a locally owned or locally relevant vector), the second argument which could be a Tensor<1,problem_dimension>, and the last argument which represents a index of the global object. This function will be primarily used to write the updated nodal values, stored as Tensor<1,problem_dimension>, into the global objects.

template <std::size_t k, int k2>
const Tensor<1, k2> & tensor,
const unsigned int i)
{
static_assert(k == k2,
"The dimensions of the input arguments must agree");
for (unsigned int j = 0; j < k; ++j)
U[j].local_element(i) = tensor[j];
}
} // namespace

We are now in a position to assemble all matrices stored in OfflineData: the lumped mass entries $$m_i$$, the vector-valued matrices $$\mathbf{c}_{ij}$$ and $$\mathbf{n}_{ij} = \frac{\mathbf{c}_{ij}}{|\mathbf{c}_{ij}|}$$, and the boundary normals $$\boldsymbol{\nu}_i$$.

In order to exploit thread parallelization we use the WorkStream approach detailed in the Parallel computing with multiple processors accessing shared memory. As customary this requires definition of

• Scratch data (i.e. input info required to carry out computations): in this case it is scratch_data.
• The worker: in our case this is the local_assemble_system() function that actually computes the local (i.e. current cell) contributions from the scratch data.
• A copy data: a struct that contains all the local assembly contributions, in this case CopyData<dim>().
• A copy data routine: in this case it is copy_local_to_global() in charge of actually coping these local contributions into the global objects (matrices and/or vectors)

Most of the following lines are spent in the definition of the worker local_assemble_system() and the copy data routine copy_local_to_global(). There is not much to say about the WorkStream framework since the vast majority of ideas are reasonably well-documented in step-9, step-13 and step-32 among others.

Finally, assuming that $$\mathbf{x}_i$$ is a support point at the boundary, the (nodal) normals are defined as:

\begin{align*} \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|} \end{align*}

We will compute the numerator of this expression first and store it in OfflineData<dim>::BoundaryNormalMap. We will normalize these vectors in a posterior loop.

template <int dim>
{
lumped_mass_matrix = 0.;
norm_matrix = 0.;
for (auto &matrix : cij_matrix)
matrix = 0.;
for (auto &matrix : nij_matrix)
matrix = 0.;
unsigned int dofs_per_cell =
discretization->finite_element.n_dofs_per_cell();

What follows is the initialization of the scratch data required by WorkStream

discretization->mapping,
discretization->finite_element,
{
computing_timer,
"offline_data - assemble lumped mass matrix, and c_ij");
const auto local_assemble_system =
[&](const typename DoFHandler<dim>::cell_iterator &cell,
CopyData<dim> & copy) {
copy.is_artificial = cell->is_artificial();
if (copy.is_artificial)
return;
copy.local_boundary_normal_map.clear();
copy.cell_lumped_mass_matrix.reinit(dofs_per_cell, dofs_per_cell);
for (auto &matrix : copy.cell_cij_matrix)
matrix.reinit(dofs_per_cell, dofs_per_cell);
const auto &fe_values = scratch.reinit(cell);
copy.local_dof_indices.resize(dofs_per_cell);
cell->get_dof_indices(copy.local_dof_indices);
std::transform(copy.local_dof_indices.begin(),
copy.local_dof_indices.end(),
copy.local_dof_indices.begin(),
return partitioner->global_to_local(index);
});
void reinit(const Triangulation< dim, spacedim > &tria)
@ update_values
Shape function values.
@ update_normal_vectors
Normal vectors.
@ update_JxW_values
typename ActiveSelector::cell_iterator cell_iterator
Definition: dof_handler.h:466
void copy(const T *begin, const T *end, U *dest)

We compute the local contributions for the lumped mass matrix entries $$m_i$$ and and vectors $$c_{ij}$$ in the usual fashion:

for (unsigned int q_point = 0; q_point < n_q_points; ++q_point)
{
const auto JxW = fe_values.JxW(q_point);
for (unsigned int j = 0; j < dofs_per_cell; ++j)
{
const auto value_JxW =
fe_values.shape_value(j, q_point) * JxW;
copy.cell_lumped_mass_matrix(j, j) += value_JxW;
for (unsigned int i = 0; i < dofs_per_cell; ++i)
{
const auto value = fe_values.shape_value(i, q_point);
for (unsigned int d = 0; d < dim; ++d)
copy.cell_cij_matrix[d](i, j) += value * grad_JxW[d];
} /* i */
} /* j */
} /* q */
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)

Now we have to compute the boundary normals. Note that the following loop does not do much unless the element has faces on the boundary of the domain.

for (const auto f : cell->face_indices())
{
const auto face = cell->face(f);
const auto id = face->boundary_id();
if (!face->at_boundary())
continue;
const auto &fe_face_values = scratch.reinit(cell, f);
const unsigned int n_face_q_points =
for (unsigned int j = 0; j < dofs_per_cell; ++j)
{
if (!discretization->finite_element.has_support_on_face(j, f))
continue;

Note that "normal" will only represent the contributions from one of the faces in the support of the shape function phi_j. So we cannot normalize this local contribution right here, we have to take it "as is", store it and pass it to the copy data routine. The proper normalization requires an additional loop on nodes. This is done in the copy function below.

if (id == Boundaries::free_slip)
{
for (unsigned int q = 0; q < n_face_q_points; ++q)
normal += fe_face_values.normal_vector(q) *
fe_face_values.shape_value(j, q);
}
const auto index = copy.local_dof_indices[j];
Point<dim> position;
for (const auto v : cell->vertex_indices())
if (cell->vertex_dof_index(v, 0) ==
partitioner->local_to_global(index))
{
position = cell->vertex(v);
break;
}
const auto old_id =
std::get<1>(copy.local_boundary_normal_map[index]);
copy.local_boundary_normal_map[index] =
std::make_tuple(normal, std::max(old_id, id), position);
}
}
};

Last, we provide a copy_local_to_global function as required for the WorkStream

const auto copy_local_to_global = [&](const CopyData<dim> &copy) {
if (copy.is_artificial)
return;
for (const auto &it : copy.local_boundary_normal_map)
{
std::get<0>(boundary_normal_map[it.first]) +=
std::get<0>(it.second);
std::get<1>(boundary_normal_map[it.first]) =
std::max(std::get<1>(boundary_normal_map[it.first]),
std::get<1>(it.second));
std::get<2>(boundary_normal_map[it.first]) = std::get<2>(it.second);
}
copy.cell_lumped_mass_matrix);
for (int k = 0; k < dim; ++k)
{
}
};
WorkStream::run(dof_handler.begin_active(),
dof_handler.end(),
local_assemble_system,
copy_local_to_global,
scratch_data,
CopyData<dim>());
}
void run(const std::vector< std::vector< Iterator >> &colored_iterators, Worker worker, Copier copier, const ScratchData &sample_scratch_data, const CopyData &sample_copy_data, const unsigned int queue_length=2 *MultithreadInfo::n_threads(), const unsigned int chunk_size=8)
Definition: work_stream.h:1275

At this point in time we are done with the computation of $$m_i$$ and $$\mathbf{c}_{ij}$$, but so far the matrix nij_matrix contains just a copy of the matrix cij_matrix. That's not what we really want: we have to normalize its entries. In addition, we have not filled the entries of the matrix norm_matrix and the vectors stored in the map OfflineData<dim>::BoundaryNormalMap are not normalized.

In principle, this is just offline data, it doesn't make much sense to over-optimize their computation, since their cost will get amortized over the many time steps that we are going to use. However, computing/storing the entries of the matrix norm_matrix and the normalization of nij_matrix are perfect to illustrate thread-parallel node-loops:

• we want to visit every node $$i$$ in the mesh/sparsity graph,
• and for every such node we want to visit to every $$j$$ such that $$\mathbf{c}_{ij} \not \equiv 0$$.

From an algebraic point of view, this is equivalent to: visiting every row in the matrix and for each one of these rows execute a loop on the columns. Node-loops is a core theme of this tutorial step (see the pseudo-code in the introduction) that will repeat over and over again. That's why this is the right time to introduce them.

We have the thread parallelization capability parallel::apply_to_subranges() that is somehow more general than the WorkStream framework. In particular, parallel::apply_to_subranges() can be used for our node-loops. This functionality requires four input arguments which we explain in detail (for the specific case of our thread-parallel node loops):

• The iterator indices.begin() points to a row index.
• The iterator indices.end() points to a numerically higher row index.
• The function on_subranges(i1,i2) (where i1 and i2 define a sub-range within the range spanned by the end and begin iterators defined in the two previous bullets) applies an operation to every iterator in such subrange. We may as well call on_subranges the "worker".
• Grainsize: minimum number of iterators (in this case representing rows) processed by each thread. We decided for a minimum of 4096 rows.

A minor caveat here is that the iterators indices.begin() and indices.end() supplied to parallel::apply_to_subranges() have to be random access iterators: internally, parallel::apply_to_subranges() will break the range defined by the indices.begin() and indices.end() iterators into subranges (we want to be able to read any entry in those subranges with constant complexity). In order to provide such iterators we resort to std_cxx20::ranges::iota_view.

The bulk of the following piece of code is spent defining the "worker" on_subranges: i.e. the operation applied at each row of the sub-range. Given a fixed row_index we want to visit every column/entry in such row. In order to execute such columns-loops we use std::for_each from the standard library, where:

• sparsity_pattern.begin(row_index) gives us an iterator starting at the first column of the row,
• sparsity_pattern.end(row_index) is an iterator pointing at the last column of the row,
• the last argument required by std::for_each is the operation applied at each nonzero entry (a lambda expression in this case) of such row.

We note that, parallel::apply_to_subranges() will operate on disjoint sets of rows (the subranges) and our goal is to write into these rows. Because of the simple nature of the operations we want to carry out (computation and storage of normals, and normalization of the $$\mathbf{c}_{ij}$$ of entries) threads cannot conflict attempting to write the same entry (we do not need a scheduler).

{
TimerOutput::Scope scope(computing_timer,
"offline_data - compute |c_ij|, and n_ij");
0, n_locally_relevant);
const auto on_subranges =
[&](const auto i1, const auto i2) {
for (const auto row_index :
*i2))
{
boost::integer_range< IncrementableType > iota_view
Definition: iota_view.h:46

First column-loop: we compute and store the entries of the matrix norm_matrix and write normalized entries into the matrix nij_matrix:

std::for_each(
sparsity_pattern.begin(row_index),
sparsity_pattern.end(row_index),
[&](const ::SparsityPatternIterators::Accessor &jt) {
const auto c_ij = gather_get_entry(cij_matrix, &jt);
const double norm = c_ij.norm();
set_entry(norm_matrix, &jt, norm);
for (unsigned int j = 0; j < dim; ++j)
set_entry(nij_matrix[j], &jt, c_ij[j] / norm);
});
}
};
indices.end(),
on_subranges,
4096);
void apply_to_subranges(const RangeType &begin, const typename identity< RangeType >::type &end, const Function &f, const unsigned int grainsize)
Definition: parallel.h:382

Finally, we normalize the vectors stored in OfflineData<dim>::BoundaryNormalMap. This operation has not been thread parallelized as it would neither illustrate any important concept nor lead to any noticeable speed gain.

for (auto &it : boundary_normal_map)
{
auto &normal = std::get<0>(it.second);
normal /= (normal.norm() + std::numeric_limits<double>::epsilon());
}
}
}
numbers::NumberTraits< Number >::real_type norm() const
SymmetricTensor< 2, dim, Number > epsilon(const Tensor< 2, dim, Number > &Grad_u)

At this point we are very much done with anything related to offline data.

#### Equation of state and approximate Riemann solver

In this section we describe the implementation of the class members of the ProblemDescription class. Most of the code here is specific to the 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 classes (in particular those defining loop structures) would remain unchanged.

We start by implementing a number of small member functions for computing momentum, internal_energy, pressure, speed_of_sound, and the flux f of the system. The functionality of each one of these functions is self-explanatory from their names.

template <int dim>
ProblemDescription<dim>::momentum(const state_type &U)
{
std::copy_n(&U[1], dim, &result[0]);
return result;
}
template <int dim>
DEAL_II_ALWAYS_INLINE inline double
ProblemDescription<dim>::internal_energy(const state_type &U)
{
const double &rho = U[0];
const auto m = momentum(U);
const double &E = U[dim + 1];
return E - 0.5 * m.norm_square() / rho;
}
template <int dim>
DEAL_II_ALWAYS_INLINE inline double
ProblemDescription<dim>::pressure(const state_type &U)
{
return (gamma - 1.) * internal_energy(U);
}
template <int dim>
DEAL_II_ALWAYS_INLINE inline double
ProblemDescription<dim>::speed_of_sound(const state_type &U)
{
const double &rho = U[0];
const double p = pressure(U);
return std::sqrt(gamma * p / rho);
}
template <int dim>
DEAL_II_ALWAYS_INLINE inline typename ProblemDescription<dim>::flux_type
ProblemDescription<dim>::flux(const state_type &U)
{
const double &rho = U[0];
const auto m = momentum(U);
const auto p = pressure(U);
const double &E = U[dim + 1];
flux_type result;
result[0] = m;
for (unsigned int i = 0; i < dim; ++i)
{
result[1 + i] = m * m[i] / rho;
result[1 + i][i] += p;
}
result[dim + 1] = m / rho * (E + p);
return result;
}
static constexpr double E
Definition: numbers.h:208
::VectorizedArray< Number, width > sqrt(const ::VectorizedArray< Number, width > &)

Now we discuss the computation of $$\lambda_{\text{max}} (\mathbf{U}_i^{n},\mathbf{U}_j^{n}, \textbf{n}_{ij})$$. The analysis and derivation of sharp upper-bounds of maximum wavespeeds of Riemann problems is a very technical endeavor and we cannot include an advanced discussion about it in this tutorial. In this portion of the documentation we will limit ourselves to sketch the main functionality of our implementation functions and point to specific academic references in order to help the (interested) reader trace the source (and proper mathematical justification) of these ideas.

In general, obtaining a sharp guaranteed upper-bound on the maximum wavespeed requires solving a quite expensive scalar nonlinear problem. This is typically done with an iterative solver. In order to simplify the presentation in this example step we decided not to include such an iterative scheme. Instead, we will just use an initial guess as a guess for an upper bound on the maximum wavespeed. More precisely, equations (2.11) (3.7), (3.8) and (4.3) of [GuermondPopov2016b] are enough to define a guaranteed upper bound on the maximum wavespeed. This estimate is returned by a call to the function lambda_max_two_rarefaction(). At its core the construction of such an upper bound uses the so-called two-rarefaction approximation for the intermediate pressure $$p^*$$, see for instance Equation (4.46), page 128 in [Toro2009].

The estimate returned by lambda_max_two_rarefaction() 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 lambda_max_expansion(). The most important function here is compute_lambda_max() which takes the minimum between the estimates returned by lambda_max_two_rarefaction() and lambda_max_expansion().

We start again by defining a couple of helper functions:

The first function takes a state U and a unit vector n_ij and computes the projected 1D state in direction of the unit vector.

namespace
{
template <int dim>
DEAL_II_ALWAYS_INLINE inline std::array<double, 4> riemann_data_from_state(
const typename ProblemDescription<dim>::state_type U,
const Tensor<1, dim> & n_ij)
{
Tensor<1, 3> projected_U;
projected_U[0] = U[0];

For this, we have to change the momentum to $$\textbf{m}\cdot n_{ij}$$ and have to subtract the kinetic energy of the perpendicular part from the total energy:

const auto m = ProblemDescription<dim>::momentum(U);
projected_U[1] = n_ij * m;
const auto perpendicular_m = m - projected_U[1] * n_ij;
projected_U[2] = U[1 + dim] - 0.5 * perpendicular_m.norm_square() / U[0];

We return the 1D state in primitive variables instead of conserved quantities. The return array consists of density $$\rho$$, velocity $$u$$, pressure $$p$$ and local speed of sound $$a$$:

return {{projected_U[0],
projected_U[1] / projected_U[0],
ProblemDescription<1>::pressure(projected_U),
ProblemDescription<1>::speed_of_sound(projected_U)}};
}

At this point we also define two small functions that return the positive and negative part of a double.

DEAL_II_ALWAYS_INLINE inline double positive_part(const double number)
{
return std::max(number, 0.);
}
DEAL_II_ALWAYS_INLINE inline double negative_part(const double number)
{
return -std::min(number, 0.);
}

Next, we need two local wavenumbers that are defined in terms of a primitive state $$[\rho, u, p, a]$$ and a given pressure $$p^\ast$$ [GuermondPopov2016] Eqn. (3.7):

\begin{align*} \lambda^- = u - a\,\sqrt{1 + \frac{\gamma+1}{2\gamma} \left(\frac{p^\ast-p}{p}\right)_+} \end{align*}

Here, the $$(\cdot)_{+}$$ denotes the positive part of the given argument.

DEAL_II_ALWAYS_INLINE inline double
lambda1_minus(const std::array<double, 4> &riemann_data,
const double p_star)
{
/* Implements formula (3.7) in Guermond-Popov-2016 */
constexpr double gamma = ProblemDescription<1>::gamma;
const auto u = riemann_data[1];
const auto p = riemann_data[2];
const auto a = riemann_data[3];
const double factor = (gamma + 1.0) / 2.0 / gamma;
const double tmp = positive_part((p_star - p) / p);
return u - a * std::sqrt(1.0 + factor * tmp);
}

Analougously [GuermondPopov2016] Eqn. (3.8):

\begin{align*} \lambda^+ = u + a\,\sqrt{1 + \frac{\gamma+1}{2\gamma} \left(\frac{p^\ast-p}{p}\right)_+} \end{align*}

DEAL_II_ALWAYS_INLINE inline double
lambda3_plus(const std::array<double, 4> &riemann_data, const double p_star)
{
/* Implements formula (3.8) in Guermond-Popov-2016 */
constexpr double gamma = ProblemDescription<1>::gamma;
const auto u = riemann_data[1];
const auto p = riemann_data[2];
const auto a = riemann_data[3];
const double factor = (gamma + 1.0) / 2.0 / gamma;
const double tmp = positive_part((p_star - p) / p);
return u + a * std::sqrt(1.0 + factor * tmp);
}

All that is left to do is to compute the maximum of $$\lambda^-$$ and $$\lambda^+$$ computed from the left and right primitive state ([GuermondPopov2016] Eqn. (2.11)), where $$p^\ast$$ is given by [GuermondPopov2016] Eqn (4.3):

DEAL_II_ALWAYS_INLINE inline double
lambda_max_two_rarefaction(const std::array<double, 4> &riemann_data_i,
const std::array<double, 4> &riemann_data_j)
{
constexpr double gamma = ProblemDescription<1>::gamma;
const auto u_i = riemann_data_i[1];
const auto p_i = riemann_data_i[2];
const auto a_i = riemann_data_i[3];
const auto u_j = riemann_data_j[1];
const auto p_j = riemann_data_j[2];
const auto a_j = riemann_data_j[3];
const double numerator = a_i + a_j - (gamma - 1.) / 2. * (u_j - u_i);
const double denominator =
a_i * std::pow(p_i / p_j, -1. * (gamma - 1.) / 2. / gamma) + a_j * 1.;
/* Formula (4.3) in Guermond-Popov-2016 */
const double p_star =
p_j * std::pow(numerator / denominator, 2. * gamma / (gamma - 1));
const double lambda1 = lambda1_minus(riemann_data_i, p_star);
const double lambda3 = lambda3_plus(riemann_data_j, p_star);
/* Formula (2.11) in Guermond-Popov-2016 */
return std::max(positive_part(lambda3), negative_part(lambda1));
}
::VectorizedArray< Number, width > pow(const ::VectorizedArray< Number, width > &, const Number p)

We compute the second upper bound of the maximal wavespeed that is, in general, not as sharp as the two-rarefaction estimate. But it will save the day in the context of near vacuum conditions when the two-rarefaction approximation might attain extreme values:

\begin{align*} \lambda_{\text{exp}} = \max(u_i,u_j) + 5. \max(a_i, a_j). \end{align*}

Note
The constant 5.0 multiplying the maximum of the sound speeds is neither an ad-hoc constant, nor a tuning parameter. It defines an upper bound for any $$\gamma \in [0,5/3]$$. Do not play with it!
DEAL_II_ALWAYS_INLINE inline double
lambda_max_expansion(const std::array<double, 4> &riemann_data_i,
const std::array<double, 4> &riemann_data_j)
{
const auto u_i = riemann_data_i[1];
const auto a_i = riemann_data_i[3];
const auto u_j = riemann_data_j[1];
const auto a_j = riemann_data_j[3];
return std::max(std::abs(u_i), std::abs(u_j)) + 5. * std::max(a_i, a_j);
}
} // namespace

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.

template <int dim>
DEAL_II_ALWAYS_INLINE inline double
ProblemDescription<dim>::compute_lambda_max(const state_type & U_i,
const state_type & U_j,
const Tensor<1, dim> &n_ij)
{
const auto riemann_data_i = riemann_data_from_state(U_i, n_ij);
const auto riemann_data_j = riemann_data_from_state(U_j, n_ij);
const double lambda_1 =
lambda_max_two_rarefaction(riemann_data_i, riemann_data_j);
const double lambda_2 =
lambda_max_expansion(riemann_data_i, riemann_data_j);
return std::min(lambda_1, lambda_2);
}

We conclude this section by defining static arrays component_names that contain strings describing the component names of our state vector. We have template specializations for dimensions one, two and three, that are used later in DataOut for naming the corresponding components:

template <>
const std::array<std::string, 3> ProblemDescription<1>::component_names{
{"rho", "m", "E"}};
template <>
const std::array<std::string, 4> ProblemDescription<2>::component_names{
{"rho", "m_1", "m_2", "E"}};
template <>
const std::array<std::string, 5> ProblemDescription<3>::component_names{
{"rho", "m_1", "m_2", "m_3", "E"}};

#### Initial values

The last preparatory step, before we discuss the implementation of the forward Euler scheme, is to briefly implement the InitialValues class.

In the constructor we initialize all parameters with default values, declare all parameters for the ParameterAcceptor class and connect the parse_parameters_call_back slot to the respective signal.

The parse_parameters_call_back slot will be invoked from ParameterAceptor after the call to ParameterAcceptor::initialize(). In that regard, its use is appropriate for situations where the parameters have to be postprocessed (in some sense) or some consistency condition between the parameters has to be checked.

template <int dim>
InitialValues<dim>::InitialValues(const std::string &subsection)
: ParameterAcceptor(subsection)
{
/* We wire up the slot InitialValues<dim>::parse_parameters_callback to
* the ParameterAcceptor::parse_parameters_call_back signal: */
std::bind(&InitialValues<dim>::parse_parameters_callback, this));
initial_direction[0] = 1.;
initial_direction,
"Initial direction of the uniform flow field");
initial_1d_state[0] = ProblemDescription<dim>::gamma;
initial_1d_state[1] = 3.;
initial_1d_state[2] = 1.;
initial_1d_state,
"Initial 1d state (rho, u, p) of the uniform flow field");
}
boost::signals2::signal< void()> parse_parameters_call_back

So far the constructor of InitialValues has defined default values for the two private members initial_direction and initial_1d_state and added them to the parameter list. But we have not defined an implementation of the only public member that we really care about, which is initial_state() (the function that we are going to call to actually evaluate the initial solution at the mesh nodes). At the top of the function, we have to ensure that the provided initial direction is not the zero vector.

Note
As commented, we could have avoided using the method parse_parameters_call_back  and defined a class member setup() in order to define the implementation of initial_state(). But for illustrative purposes we want to document a different way here and use the call back signal from ParameterAcceptor.
template <int dim>
void InitialValues<dim>::parse_parameters_callback()
{
AssertThrow(initial_direction.norm() != 0.,
"Initial shock front direction is set to the zero vector."));
initial_direction /= initial_direction.norm();
static ::ExceptionBase & ExcMessage(std::string arg1)
#define AssertThrow(cond, exc)
Definition: exceptions.h:1583

Next, we implement the initial_state function object with a lambda function computing a uniform flow field. For this we have to translate a given primitive 1d state (density $$\rho$$, velocity $$u$$, and pressure $$p$$) into a conserved n-dimensional state (density $$\rho$$, momentum $$\mathbf{m}$$, and total energy $$E$$).

initial_state = [this](const Point<dim> & /*point*/, double /*t*/) {
const double rho = initial_1d_state[0];
const double u = initial_1d_state[1];
const double p = initial_1d_state[2];
static constexpr double gamma = ProblemDescription<dim>::gamma;
state_type state;
state[0] = rho;
for (unsigned int i = 0; i < dim; ++i)
state[1 + i] = rho * u * initial_direction[i];
state[dim + 1] = p / (gamma - 1.) + 0.5 * rho * u * u;
return state;
};
}

#### The Forward Euler step

The constructor of the TimeStepping class does not contain any surprising code:

template <int dim>
const MPI_Comm mpi_communicator,
TimerOutput & computing_timer,
const OfflineData<dim> & offline_data,
const InitialValues<dim> &initial_values,
const std::string & subsection /*= "TimeStepping"*/)
: ParameterAcceptor(subsection)
, mpi_communicator(mpi_communicator)
, computing_timer(computing_timer)
, offline_data(&offline_data)
, initial_values(&initial_values)
{
cfl_update = 0.80;
cfl_update,
"Relative CFL constant used for update");
}

In the class member prepare() we initialize the temporary vector temp and the matrix dij_matrix. The vector temp will be used to store the solution update temporarily before its contents is swapped with the old vector.

template <int dim>
void TimeStepping<dim>::prepare()
{
TimerOutput::Scope scope(computing_timer,
"time_stepping - prepare scratch space");
for (auto &it : temporary_vector)
it.reinit(offline_data->partitioner);
dij_matrix.reinit(offline_data->sparsity_pattern);
}

It is now time to implement the forward Euler step. Given a (writable reference) to the old state U at time $$t$$ we update the state U in place and return the chosen time-step size. We first declare a number of read-only references to various different variables and data structures. We do this is mainly to have shorter variable names (e.g., sparsity instead of offline_data->sparsity_pattern).

template <int dim>
double TimeStepping<dim>::make_one_step(vector_type &U, double t)
{
const auto &n_locally_owned = offline_data->n_locally_owned;
const auto &n_locally_relevant = offline_data->n_locally_relevant;
indices_owned(0, n_locally_owned);
indices_relevant(0, n_locally_relevant);
const auto &sparsity = offline_data->sparsity_pattern;
const auto &lumped_mass_matrix = offline_data->lumped_mass_matrix;
const auto &norm_matrix = offline_data->norm_matrix;
const auto &nij_matrix = offline_data->nij_matrix;
const auto &cij_matrix = offline_data->cij_matrix;
const auto &boundary_normal_map = offline_data->boundary_normal_map;

Step 1: Computing the $$d_{ij}$$ graph viscosity matrix.

It is important to highlight that the viscosity matrix has to be symmetric, i.e., $$d_{ij} = d_{ji}$$. In this regard we note here that $$\int_{\Omega} \nabla \phi_j \phi_i \, \mathrm{d}\mathbf{x}= - \int_{\Omega} \nabla \phi_i \phi_j \, \mathrm{d}\mathbf{x}$$ (or equivalently $$\mathbf{c}_{ij} = - \mathbf{c}_{ji}$$) provided either $$\mathbf{x}_i$$ or $$\mathbf{x}_j$$ is a support point located away from the boundary. In this case we can check that $$\lambda_{\text{max}} (\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})$$ by construction, which guarantees the property $$d_{ij} = d_{ji}$$.

However, if both support points $$\mathbf{x}_i$$ or $$\mathbf{x}_j$$ happen to lie on the boundary, then, the equalities $$\mathbf{c}_{ij} = - \mathbf{c}_{ji}$$ and $$\lambda_{\text{max}} (\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})$$ 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 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 the upper-triangular entries of $$d_{ij}$$ and copy the corresponding entries to the lower-triangular counterpart.

We use again parallel::apply_to_subranges() for thread-parallel for loops. Pretty much all the ideas for parallel traversal that we introduced when discussing the assembly of the matrix norm_matrix and the normalization of nij_matrix above are used here again.

We define again a "worker" function on_subranges that computes the viscosity $$d_{ij}$$ for a subrange [i1, i2) of column indices:

{
TimerOutput::Scope scope(computing_timer,
"time_stepping - 1 compute d_ij");
const auto on_subranges =
[&](const auto i1, const auto i2) {
for (const auto i :
*i2))
{
const auto U_i = gather(U, i);

For a given column index i we iterate over the columns of the sparsity pattern from sparsity.begin(i) to sparsity.end(i):

for (auto jt = sparsity.begin(i); jt != sparsity.end(i); ++jt)
{
const auto j = jt->column();

We only compute $$d_{ij}$$ if $$j < i$$ (upper triangular entries) and later copy the values over to $$d_{ji}$$.

if (j >= i)
continue;
const auto U_j = gather(U, j);
const auto n_ij = gather_get_entry(nij_matrix, jt);
const double norm = get_entry(norm_matrix, jt);
const auto lambda_max =
ProblemDescription<dim>::compute_lambda_max(U_i, U_j, n_ij);
double d = norm * lambda_max;
double norm(const FEValuesBase< dim > &fe, const ArrayView< const std::vector< Tensor< 1, dim >>> &Du)
Definition: divergence.h:472

If both support points happen to be at the boundary we have to compute $$d_{ji}$$ as well and then take $$\max(d_{ij},d_{ji})$$. After this we can finally set the upper triangular and lower triangular entries.

if (boundary_normal_map.count(i) != 0 &&
boundary_normal_map.count(j) != 0)
{
const auto n_ji = gather(nij_matrix, j, i);
const auto lambda_max_2 =
ProblemDescription<dim>::compute_lambda_max(U_j,
U_i,
n_ji);
const double norm_2 = norm_matrix(j, i);
d = std::max(d, norm_2 * lambda_max_2);
}
set_entry(dij_matrix, jt, d);
dij_matrix(j, i) = d;
}
}
};
parallel::apply_to_subranges(indices_relevant.begin(),
indices_relevant.end(),
on_subranges,
4096);
}

Step 2: Compute diagonal entries $$d_{ii}$$ and $$\tau_{\text{max}}$$.

So far we have computed all off-diagonal entries of the matrix dij_matrix. We still have to fill its diagonal entries defined as $$d_{ii}^n = - \sum_{j \in \mathcal{I}(i)\backslash \{i\}} d_{ij}^n$$. We use again parallel::apply_to_subranges() for this purpose. While computing the $$d_{ii}$$s we also determine the largest admissible time-step, which is defined as

$\tau_n \dealcoloneq c_{\text{cfl}}\,\min_{i\in\mathcal{V}} \left(\frac{m_i}{-2\,d_{ii}^{n}}\right) \, .$

Note that the operation $$\min_{i \in \mathcal{V}}$$ is intrinsically global, it operates on all nodes: first we have to take the minimum over all threads (of a given node) and then we have to take the minimum over all MPI processes. In the current implementation:

• We store tau_max (per node) as std::atomic<double>. The internal implementation of std::atomic will take care of guarding any possible race condition when more than one thread attempts to read and/or write tau_max at the same time.
• In order to take the minimum over all MPI process we use the utility function Utilities::MPI::min.
std::atomic<double> tau_max{std::numeric_limits<double>::infinity()};
{
TimerOutput::Scope scope(computing_timer,
"time_stepping - 2 compute d_ii, and tau_max");

on_subranges() will be executed on every thread individually. The variable tau_max_on_subrange is thus stored thread locally.

const auto on_subranges =
[&](const auto i1, const auto i2) {
double tau_max_on_subrange = std::numeric_limits<double>::infinity();
for (const auto i :
*i2))
{
double d_sum = 0.;
for (auto jt = sparsity.begin(i); jt != sparsity.end(i); ++jt)
{
const auto j = jt->column();
if (j == i)
continue;
d_sum -= get_entry(dij_matrix, jt);
}

We store the negative sum of the d_ij entries at the diagonal position

dij_matrix.diag_element(i) = d_sum;

and compute the maximal local time-step size tau:

const double mass = lumped_mass_matrix.diag_element(i);
const double tau = cfl_update * mass / (-2. * d_sum);
tau_max_on_subrange = std::min(tau_max_on_subrange, tau);
}

tau_max_on_subrange contains the largest possible time-step size computed for the (thread local) subrange. At this point we have to synchronize the value over all threads. This is were we use the std::atomic<double> compare exchange update mechanism:

while (current_tau_max > tau_max_on_subrange &&
!tau_max.compare_exchange_weak(current_tau_max,
tau_max_on_subrange))
;
};
parallel::apply_to_subranges(indices_relevant.begin(),
indices_relevant.end(),
on_subranges,
4096);

After all threads have finished we can simply synchronize the value over all MPI processes:

T min(const T &t, const MPI_Comm &mpi_communicator)

This is a good point to verify that the computed tau_max is indeed a valid floating point number.

"I'm sorry, Dave. I'm afraid I can't do that. - We crashed."));
}

Step 3: Perform update.

At this point, we have computed all viscosity coefficients $$d_{ij}$$ and we know the maximal admissible time-step size $$\tau_{\text{max}}$$. This means we can now compute the update:

$\mathbf{U}_i^{n+1} = \mathbf{U}_i^{n} - \frac{\tau_{\text{max}} }{m_i} \sum_{j \in \mathcal{I}(i)} (\mathbb{f}(\mathbf{U}_j^{n}) - \mathbb{f}(\mathbf{U}_i^{n})) \cdot \mathbf{c}_{ij} - d_{ij} (\mathbf{U}_j^{n} - \mathbf{U}_i^{n})$

This update formula is slightly different from what was discussed in the introduction (in the pseudo-code). However, it can be shown that both equations are algebraically equivalent (they will produce the same numerical values). We favor this second formula since it has natural cancellation properties that might help avoid numerical artifacts.

{
TimerOutput::Scope scope(computing_timer,
"time_stepping - 3 perform update");
const auto on_subranges =
[&](const auto i1, const auto i2) {
for (const auto i : boost::make_iterator_range(i1, i2))
{
Assert(i < n_locally_owned, ExcInternalError());
const auto U_i = gather(U, i);
const auto f_i = ProblemDescription<dim>::flux(U_i);
const double m_i = lumped_mass_matrix.diag_element(i);
auto U_i_new = U_i;
for (auto jt = sparsity.begin(i); jt != sparsity.end(i); ++jt)
{
const auto j = jt->column();
const auto U_j = gather(U, j);
const auto f_j = ProblemDescription<dim>::flux(U_j);
const auto c_ij = gather_get_entry(cij_matrix, jt);
const auto d_ij = get_entry(dij_matrix, jt);
for (unsigned int k = 0; k < problem_dimension; ++k)
{
U_i_new[k] +=
tau_max / m_i *
(-(f_j[k] - f_i[k]) * c_ij + d_ij * (U_j[k] - U_i[k]));
}
}
scatter(temporary_vector, U_i_new, i);
}
};
parallel::apply_to_subranges(indices_owned.begin(),
indices_owned.end(),
on_subranges,
4096);
}
static ::ExceptionBase & ExcInternalError()
#define Assert(cond, exc)
Definition: exceptions.h:1473
IteratorRange< BaseIterator > make_iterator_range(const BaseIterator &begin, const typename identity< BaseIterator >::type &end)

Step 4: Fix up boundary states.

As a last step in the Forward Euler method, we have to fix up all boundary states. As discussed in the intro we

• advance in time satisfying no boundary condition at all,
• at the end of the time step enforce boundary conditions strongly in a post-processing step.

Here, we compute the correction

$\mathbf{m}_i \dealcoloneq \mathbf{m}_i - (\boldsymbol{\nu}_i \cdot \mathbf{m}_i) \boldsymbol{\nu}_i,$

which removes the normal component of $$\mathbf{m}$$.

{
TimerOutput::Scope scope(computing_timer,
"time_stepping - 4 fix boundary states");
for (auto it : boundary_normal_map)
{
const auto i = it.first;

We only iterate over the locally owned subset:

if (i >= n_locally_owned)
continue;
const auto &normal = std::get<0>(it.second);
const auto &id = std::get<1>(it.second);
const auto &position = std::get<2>(it.second);
auto U_i = gather(temporary_vector, i);

On free slip boundaries we remove the normal component of the momentum:

if (id == Boundaries::free_slip)
{
auto m = ProblemDescription<dim>::momentum(U_i);
m -= (m * normal) * normal;
for (unsigned int k = 0; k < dim; ++k)
U_i[k + 1] = m[k];
}

On Dirichlet boundaries we enforce initial conditions strongly:

else if (id == Boundaries::dirichlet)
{
U_i = initial_values->initial_state(position, t + tau_max);
}
scatter(temporary_vector, U_i, i);
}
}

Step 5: We now update the ghost layer over all MPI ranks, swap the temporary vector with the solution vector U (that will get returned by reference) and return the chosen time-step size $$\tau_{\text{max}}$$:

for (auto &it : temporary_vector)
it.update_ghost_values();
U.swap(temporary_vector);
return tau_max;
}

#### Schlieren postprocessing

At various intervals we will output the current state U of the solution together with a so-called Schlieren plot. The constructor of the SchlierenPostprocessor class again contains no surprises. We simply supply default values to and register two parameters:

• schlieren_beta: is an ad-hoc positive amplification factor in order to enhance the contrast in the visualization. Its actual value is a matter of taste.
• schlieren_index: is an integer indicating which component of the state $$[\rho, \mathbf{m},E]$$ are we going to use in order to generate the visualization.
template <int dim>
SchlierenPostprocessor<dim>::SchlierenPostprocessor(
const MPI_Comm mpi_communicator,
TimerOutput & computing_timer,
const OfflineData<dim> &offline_data,
const std::string & subsection /*= "SchlierenPostprocessor"*/)
: ParameterAcceptor(subsection)
, mpi_communicator(mpi_communicator)
, computing_timer(computing_timer)
, offline_data(&offline_data)
{
schlieren_beta = 10.;
schlieren_beta,
"Beta factor used in Schlieren-type postprocessor");
schlieren_index = 0;
schlieren_index,
"Use the corresponding component of the state vector for the "
"schlieren plot");
}

Again, the prepare() function initializes two temporary the vectors (r and schlieren).

template <int dim>
void SchlierenPostprocessor<dim>::prepare()
{
TimerOutput::Scope scope(computing_timer,
"schlieren_postprocessor - prepare scratch space");
r.reinit(offline_data->n_locally_relevant);
schlieren.reinit(offline_data->partitioner);
}

We now discuss the implementation of the class member SchlierenPostprocessor<dim>::compute_schlieren(), which basically takes a component of the state vector U and computes the Schlieren indicator for such component (the formula of the Schlieren indicator can be found just before the declaration of the class SchlierenPostprocessor). We start by noting that this formula requires the "nodal gradients" $$\nabla r_j$$. However, nodal values of gradients are not defined for $$\mathcal{C}^0$$ finite element functions. More generally, pointwise values of gradients are not defined for $$W^{1,p}(\Omega)$$ functions. The simplest technique we can use to recover gradients at nodes is weighted-averaging i.e.

$\nabla r_j \dealcoloneq \frac{1}{\int_{S_i} \omega_i(\mathbf{x}) \, \mathrm{d}\mathbf{x}} \int_{S_i} r_h(\mathbf{x}) \omega_i(\mathbf{x}) \, \mathrm{d}\mathbf{x} \ \ \ \ \ \mathbf{(*)}$

where $$S_i$$ is the support of the shape function $$\phi_i$$, and $$\omega_i(\mathbf{x})$$ is the weight. The weight could be any positive function such as $$\omega_i(\mathbf{x}) \equiv 1$$ (that would allow us to recover the usual notion of mean value). But as usual, the goal is to reuse the off-line data as much as possible. In this sense, the most natural choice of weight is $$\omega_i = \phi_i$$. Inserting this choice of weight and the expansion $$r_h(\mathbf{x}) = \sum_{j \in \mathcal{V}} r_j \phi_j(\mathbf{x})$$ into $$\mathbf{(*)}$$ we get :

$\nabla r_j \dealcoloneq \frac{1}{m_i} \sum_{j \in \mathcal{I}(i)} r_j \mathbf{c}_{ij} \ \ \ \ \ \mathbf{(**)} \, .$

Using this last formula we can recover averaged nodal gradients without resorting to any form of quadrature. This idea aligns quite well with the whole spirit of edge-based schemes (or algebraic schemes) where we want to operate on matrices and vectors as directly as it could be possible avoiding by all means assembly of bilinear forms, cell-loops, quadrature, or any other intermediate construct/operation between the input arguments (the state from the previous time-step) and the actual matrices and vectors required to compute the update.

The second thing to note is that we have to compute global minimum and maximum $$\max_j |\nabla r_j|$$ and $$\min_j |\nabla r_j|$$. Following the same ideas used to compute the time step size in the class member TimeStepping<dim>::step() we define $$\max_j |\nabla r_j|$$ and $$\min_j |\nabla r_j|$$ as atomic doubles in order to resolve any conflicts between threads. As usual, we use Utilities::MPI::max() and Utilities::MPI::min() to find the global maximum/minimum among all MPI processes.

Finally, it is not possible to compute the Schlieren indicator in a single loop over all nodes. The entire operation requires two loops over nodes:

• The first loop computes $$|\nabla r_i|$$ for all $$i \in \mathcal{V}$$ in the mesh, and the bounds $$\max_j |\nabla r_j|$$ and $$\min_j |\nabla r_j|$$.
• The second loop finally computes the Schlieren indicator using the formula

$\text{schlieren}[i] = e^{\beta \frac{ |\nabla r_i| - \min_j |\nabla r_j| }{\max_j |\nabla r_j| - \min_j |\nabla r_j| } } \, .$

This means that we will have to define two workers on_subranges for each one of these stages.

template <int dim>
void SchlierenPostprocessor<dim>::compute_schlieren(const vector_type &U)
{
computing_timer, "schlieren_postprocessor - compute schlieren plot");
const auto &sparsity = offline_data->sparsity_pattern;
const auto &lumped_mass_matrix = offline_data->lumped_mass_matrix;
const auto &cij_matrix = offline_data->cij_matrix;
const auto &boundary_normal_map = offline_data->boundary_normal_map;
const auto &n_locally_owned = offline_data->n_locally_owned;
const auto indices =
n_locally_owned);

We define the r_i_max and r_i_min in the current MPI process as atomic doubles in order to avoid race conditions between threads:

std::atomic<double> r_i_max{0.};
std::atomic<double> r_i_min{std::numeric_limits<double>::infinity()};

First loop: compute the averaged gradient at each node and the global maxima and minima of the gradients.

{
const auto on_subranges =
[&](const auto i1, const auto i2) {
double r_i_max_on_subrange = 0.;
double r_i_min_on_subrange = std::numeric_limits<double>::infinity();
for (const auto i : boost::make_iterator_range(i1, i2))
{
Assert(i < n_locally_owned, ExcInternalError());
for (auto jt = sparsity.begin(i); jt != sparsity.end(i); ++jt)
{
const auto j = jt->column();
if (i == j)
continue;
const auto U_js = U[schlieren_index].local_element(j);
const auto c_ij = gather_get_entry(cij_matrix, jt);
r_i += c_ij * U_js;
}

We fix up the gradient r_i at free slip boundaries similarly to how we fixed up boundary states in the forward Euler step. This avoids sharp, artificial gradients in the Schlieren plot at free slip boundaries and is a purely cosmetic choice.

const auto bnm_it = boundary_normal_map.find(i);
if (bnm_it != boundary_normal_map.end())
{
const auto &normal = std::get<0>(bnm_it->second);
const auto &id = std::get<1>(bnm_it->second);
if (id == Boundaries::free_slip)
r_i -= 1. * (r_i * normal) * normal;
else
r_i = 0.;
}

We remind the reader that we are not interested in the nodal gradients per se. We only want their norms in order to compute the Schlieren indicator (weighted with the lumped mass matrix $$m_i$$):

const double m_i = lumped_mass_matrix.diag_element(i);
r[i] = r_i.norm() / m_i;
r_i_max_on_subrange = std::max(r_i_max_on_subrange, r[i]);
r_i_min_on_subrange = std::min(r_i_min_on_subrange, r[i]);
}

We compare the current_r_i_max and current_r_i_min (in the current subrange) with r_i_max and r_i_min (for the current MPI process) and update them if necessary:

while (current_r_i_max < r_i_max_on_subrange &&
!r_i_max.compare_exchange_weak(current_r_i_max,
r_i_max_on_subrange))
;
while (current_r_i_min > r_i_min_on_subrange &&
!r_i_min.compare_exchange_weak(current_r_i_min,
r_i_min_on_subrange))
;
};
indices.end(),
on_subranges,
4096);
}

And synchronize r_i_max and r_i_min over all MPI processes.

T max(const T &t, const MPI_Comm &mpi_communicator)

Second loop: we now have the vector r and the scalars r_i_max and r_i_min at our disposal. We are thus in a position to actually compute the Schlieren indicator.

{
const auto on_subranges =
[&](const auto i1, const auto i2) {
for (const auto i : boost::make_iterator_range(i1, i2))
{
Assert(i < n_locally_owned, ExcInternalError());
schlieren.local_element(i) =
1. - std::exp(-schlieren_beta * (r[i] - r_i_min) /
(r_i_max - r_i_min));
}
};
indices.end(),
on_subranges,
4096);
}
::VectorizedArray< Number, width > exp(const ::VectorizedArray< Number, width > &)

And finally, exchange ghost elements.

schlieren.update_ghost_values();
}

#### The main loop

With all classes implemented it is time to create an instance of Discretization<dim>, OfflineData<dim>, InitialValues<dim>, TimeStepping<dim>, and SchlierenPostprocessor<dim>, and run the forward Euler step in a loop.

In the constructor of MainLoop<dim> we now initialize an instance of all classes, and declare a number of parameters controlling output. Most notable, we declare a boolean parameter resume that will control whether the program attempts to restart from an interrupted computation, or not.

template <int dim>
MainLoop<dim>::MainLoop(const MPI_Comm mpi_communicator)
: ParameterAcceptor("A - MainLoop")
, mpi_communicator(mpi_communicator)
, computing_timer(mpi_communicator,
timer_output,
TimerOutput::never,
TimerOutput::cpu_and_wall_times)
, pcout(std::cout, Utilities::MPI::this_mpi_process(mpi_communicator) == 0)
, discretization(mpi_communicator, computing_timer, "B - Discretization")
, offline_data(mpi_communicator,
computing_timer,
discretization,
"C - OfflineData")
, initial_values("D - InitialValues")
, time_stepping(mpi_communicator,
computing_timer,
offline_data,
initial_values,
"E - TimeStepping")
, schlieren_postprocessor(mpi_communicator,
computing_timer,
offline_data,
"F - SchlierenPostprocessor")
{
base_name = "test";
add_parameter("basename", base_name, "Base name for all output files");
t_final = 4.;
output_granularity = 0.02;
output_granularity,
"time interval for output");
asynchronous_writeback = true;
asynchronous_writeback,
"Write out solution in a background thread performing IO");
resume = false;
add_parameter("resume", resume, "Resume an interrupted computation.");
}
unsigned int this_mpi_process(const MPI_Comm &mpi_communicator)
Definition: mpi.cc:152

We start by implementing a helper function print_head() in an anonymous namespace that is used to output messages in the terminal with some nice formatting.

namespace
{
const std::string & secondary = "")
{
std::string((35 - header_size) / 2, ' ');
const auto secondary_size = secondary.size();
std::string((34 - secondary_size) / 2, ' ') + secondary +
std::string((35 - secondary_size) / 2, ' ');
/* clang-format off */
pcout << std::endl;
pcout << " ####################################################" << std::endl;
pcout << " ######### #########" << std::endl;
pcout << " #########" << padded_secondary << "#########" << std::endl;
pcout << " ######### #########" << std::endl;
pcout << " ####################################################" << std::endl;
pcout << std::endl;
/* clang-format on */
}
} // namespace

With print_head in place it is now time to implement the MainLoop<dim>::run() that contains the main loop of our program.

template <int dim>
{

We start by reading in parameters and initializing all objects. We note here that the call to ParameterAcceptor::initialize reads in all parameters from the parameter file (whose name is given as a string argument). ParameterAcceptor handles a global ParameterHandler that is initialized with subsections and parameter declarations for all class instances that are derived from ParameterAceptor. The call to initialize enters the subsection for each each derived class, and sets all variables that were added using ParameterAcceptor::add_parameter()

pcout << "Reading parameters and allocating objects... " << std::flush;
pcout << "done" << std::endl;
static void initialize(const std::string &filename="", const std::string &output_filename="", const ParameterHandler::OutputStyle output_style_for_output_filename=ParameterHandler::Short, ParameterHandler &prm=ParameterAcceptor::prm, const ParameterHandler::OutputStyle output_style_for_filename=ParameterHandler::DefaultStyle)

Next we create the triangulation, assemble all matrices, set up scratch space, and initialize the DataOut<dim> object. All of these operations are pretty standard and discussed in detail in the Discretization and OfflineData classes.

{
discretization.setup();
if (resume)
else
discretization.triangulation.refine_global(discretization.refinement);
pcout << "Number of active cells: "
<< discretization.triangulation.n_global_active_cells()
<< std::endl;
offline_data.setup();
offline_data.assemble();
pcout << "Number of degrees of freedom: "
<< offline_data.dof_handler.n_dofs() << std::endl;
time_stepping.prepare();
schlieren_postprocessor.prepare();
}

We will store the current time and state in the variable t and vector U:

double t = 0.;
unsigned int output_cycle = 0;
vector_type U;
for (auto &it : U)
it.reinit(offline_data.partitioner);

##### Resume

By default the boolean resume is set to false, i.e. the following code snippet is not run. However, if resume==true we indicate that we have indeed an interrupted computation and the program shall restart by reading in an old state consisting of t, output_cycle, and the state vector U from checkpoint files.

A this point we have already read in the stored refinement history of our parallel distributed mesh. What is missing are the actual state vector U, the time and output cycle. We use the SolutionTransfer class in combination with the distributed::Triangulation::load() / distributed::Triangulation::save() mechanism to read in the state vector. A separate boost::archive is used to retrieve t and output_cycle. The checkpoint files will be created in the output() routine discussed below.

if (resume)
{
parallel::distributed::
SolutionTransfer<dim, LinearAlgebra::distributed::Vector<double>>
solution_transfer(offline_data.dof_handler);
std::vector<LinearAlgebra::distributed::Vector<double> *> vectors;
std::transform(U.begin(),
U.end(),
std::back_inserter(vectors),
[](auto &it) { return &it; });
solution_transfer.deserialize(vectors);
for (auto &it : U)
it.update_ghost_values();
std::ios::binary);
boost::archive::binary_iarchive ia(file);
ia >> t >> output_cycle;
}
else
{
U = interpolate_initial_values();
}

With either the initial state set up, or an interrupted state restored it is time to enter the main loop:

output(U, base_name, t, output_cycle++);
for (unsigned int cycle = 1; t < t_final; ++cycle)
{

We first print an informative status message

std::ostringstream secondary;
head << "Cycle " << Utilities::int_to_string(cycle, 6) << " ("
<< std::fixed << std::setprecision(1) << t / t_final * 100
<< "%)";
secondary << "at time t = " << std::setprecision(8) << std::fixed << t;
std::string int_to_string(const unsigned int value, const unsigned int digits=numbers::invalid_unsigned_int)
Definition: utilities.cc:473

and then perform a single forward Euler step. Note that the state vector U is updated in place and that time_stepping.make_one_step() returns the chosen step size.

t += time_stepping.make_one_step(U, t);

Post processing, generating output and writing out the current state is a CPU and IO intensive task that we cannot afford to do every time step - in particular with explicit time stepping. We thus only schedule output by calling the output() function if we are past a threshold set by output_granularity.

if (t > output_cycle * output_granularity)
{
checkpoint(U, base_name, t, output_cycle);
output(U, base_name, t, output_cycle);
++output_cycle;
}
}

We wait for any remaining background output thread to finish before printing a summary and exiting.

computing_timer.print_summary();
pcout << timer_output.str() << std::endl;
}

The interpolate_initial_values takes an initial time "t" as input argument and populates a state vector U with the help of the InitialValues<dim>::initial_state object.

template <int dim>
typename MainLoop<dim>::vector_type
MainLoop<dim>::interpolate_initial_values(const double t)
{
pcout << "MainLoop<dim>::interpolate_initial_values(t = " << t << ')'
<< std::endl;
TimerOutput::Scope scope(computing_timer,
"main_loop - setup scratch space");
vector_type U;
for (auto &it : U)
it.reinit(offline_data.partitioner);
constexpr auto problem_dimension =
ProblemDescription<dim>::problem_dimension;

The function signature of InitialValues<dim>::initial_state is not quite right for VectorTools::interpolate(). We work around this issue by, first, creating a lambda function that for a given position x returns just the value of the ith component. This lambda in turn is converted to a Function with the help of the ScalarFunctionFromFunctionObject wrapper.

for (unsigned int i = 0; i < problem_dimension; ++i)
VectorTools::interpolate(offline_data.dof_handler,
[&](const Point<dim> &x) {
return initial_values.initial_state(x, t)[i];
}),
U[i]);
for (auto &it : U)
it.update_ghost_values();
return U;
}
void interpolate(const Mapping< dim, spacedim > &mapping, const DoFHandler< dim, spacedim > &dof, const Function< spacedim, typename VectorType::value_type > &function, VectorType &vec, const ComponentMask &component_mask=ComponentMask())

##### Output and checkpointing

We checkpoint the current state by doing the precise inverse operation to what we discussed for the resume logic:

template <int dim>
void MainLoop<dim>::checkpoint(const typename MainLoop<dim>::vector_type &U,
const std::string &name,
const double t,
const unsigned int cycle)
{
parallel::distributed::
SolutionTransfer<dim, LinearAlgebra::distributed::Vector<double>>
solution_transfer(offline_data.dof_handler);
std::vector<const LinearAlgebra::distributed::Vector<double> *> vectors;
std::transform(U.begin(),
U.end(),
std::back_inserter(vectors),
[](auto &it) { return &it; });
solution_transfer.prepare_for_serialization(vectors);
discretization.triangulation.save(name + "-checkpoint.mesh");
if (::Utilities::MPI::this_mpi_process(mpi_communicator) == 0)
{
boost::archive::binary_oarchive oa(file);
oa << t << cycle;
}
}

Writing out the final vtk files is quite an IO intensive task that can stall the main loop for a while. In order to avoid this we use an asynchronous IO strategy by creating a background thread that will perform IO while the main loop is allowed to continue. In order for this to work we have to be mindful of two things:

• Before running the output_worker thread, we have to create a copy of the state vector U. We store it in the vector output_vector.
• We have to avoid any MPI communication in the background thread, otherwise the program might deadlock. This implies that we have to run the postprocessing outside of the worker thread.
template <int dim>
void MainLoop<dim>::output(const typename MainLoop<dim>::vector_type &U,
const std::string & name,
const double t,
const unsigned int cycle)
{
pcout << "MainLoop<dim>::output(t = " << t << ')' << std::endl;

If the asynchronous writeback option is set we launch a background thread performing all the slow IO to disc. In that case we have to make sure that the background thread actually finished running. If not, we have to wait to for it to finish. We launch said background thread with std::async() that returns a std::future object. This std::future object contains the return value of the function, which is in our case simply void.

{
TimerOutput::Scope timer(computing_timer, "main_loop - stalled output");
}
constexpr auto problem_dimension =
ProblemDescription<dim>::problem_dimension;

At this point we make a copy of the state vector, run the schlieren postprocessor, and run DataOut<dim>::build_patches() The actual output code is standard: We create a DataOut instance, attach all data vectors we want to output and call DataOut<dim>::build_patches(). There is one twist, however. In order to perform asynchronous IO on a background thread we create the DataOut<dim> object as a shared pointer that we pass on to the worker thread to ensure that once we exit this function and the worker thread finishes the DataOut<dim> object gets destroyed again.

for (unsigned int i = 0; i < problem_dimension; ++i)
{
output_vector[i] = U[i];
output_vector[i].update_ghost_values();
}
schlieren_postprocessor.compute_schlieren(output_vector);
auto data_out = std::make_shared<DataOut<dim>>();
data_out->attach_dof_handler(offline_data.dof_handler);
const auto &component_names = ProblemDescription<dim>::component_names;
for (unsigned int i = 0; i < problem_dimension; ++i)
"schlieren_plot");
data_out->build_patches(discretization.mapping,
discretization.finite_element.degree - 1);

Next we create a lambda function for the background thread. We capture the this pointer as well as most of the arguments of the output function by value so that we have access to them inside the lambda function.

const auto output_worker = [this, name, t, cycle, data_out]() {
cycle,
true,
data_out->set_flags(flags);
data_out->write_vtu_with_pvtu_record(
"", name + "-solution", cycle, mpi_communicator, 6);
};

If the asynchronous writeback option is set we launch a new background thread with the help of std::async function. The function returns a std::future object that we can use to query the status of the background thread. At this point we can return from the output() function and resume with the time stepping in the main loop - the thread will run in the background.

if (asynchronous_writeback)
{
}
else
{
output_worker();
}
}
} // namespace Step69

And finally, the main function.

int main(int argc, char *argv[])
{
try
{
constexpr int dim = 2;
using namespace dealii;
using namespace Step69;
Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv);
MPI_Comm mpi_communicator(MPI_COMM_WORLD);
MainLoop<dim> main_loop(mpi_communicator);
main_loop.run();
}
catch (std::exception &exc)
{
std::cerr << std::endl
<< std::endl
<< "----------------------------------------------------"
<< std::endl;
std::cerr << "Exception on processing: " << std::endl
<< exc.what() << std::endl
<< "Aborting!" << std::endl
<< "----------------------------------------------------"
<< std::endl;
return 1;
}
catch (...)
{
std::cerr << std::endl
<< std::endl
<< "----------------------------------------------------"
<< std::endl;
std::cerr << "Unknown exception!" << std::endl
<< "Aborting!" << std::endl
<< "----------------------------------------------------"
<< std::endl;
return 1;
};
}

# Results

Running the program with default parameters in release mode takes about 1 minute on a 4 core machine (with hyperthreading):

# mpirun -np 4 ./step-69 | tee output
Reading parameters and allocating objects... done

####################################################
#########                                  #########
#########       create triangulation       #########
#########                                  #########
####################################################

Number of active cells:       36864

####################################################
#########                                  #########
#########       compute offline data       #########
#########                                  #########
####################################################

Number of degrees of freedom: 37376

####################################################
#########                                  #########
#########         set up time step         #########
#########                                  #########
####################################################

####################################################
#########                                  #########
#########    interpolate initial values    #########
#########                                  #########
#########                                  #########
####################################################

TimeLoop<dim>::interpolate_initial_values(t = 0)
TimeLoop<dim>::output(t = 0, checkpoint = 0)

####################################################
#########                                  #########
#########         enter main loop          #########
#########                                  #########
#########                                  #########
####################################################

####################################################
#########                                  #########
#########      Cycle  000001  (0.0%)       #########
#########      at time t = 0.00000000      #########
#########                                  #########
####################################################

[...]

####################################################
#########                                  #########
#########     Cycle  007553  (100.0%)      #########
#########      at time t = 3.99984036      #########
#########                                  #########
####################################################

TimeLoop<dim>::output(t = 4.00038, checkpoint = 1)

+------------------------------------------------------------------------+------------+------------+
| Total CPU time elapsed since start                                     |       357s |            |
|                                                                        |            |            |
| Section                                                    | no. calls |  CPU time  | % of total |
+------------------------------------------------------------+-----------+------------+------------+
| discretization - setup                                     |         1 |     0.113s |         0% |
| offline_data - assemble lumped mass matrix, and c_ij       |         1 |     0.167s |         0% |
| offline_data - compute |c_ij|, and n_ij                    |         1 |   0.00255s |         0% |
| offline_data - create sparsity pattern and set up matrices |         1 |    0.0224s |         0% |
| offline_data - distribute dofs                             |         1 |    0.0617s |         0% |
| offline_data - fix slip boundary c_ij                      |         1 |    0.0329s |         0% |
| schlieren_postprocessor - compute schlieren plot           |       201 |     0.811s |      0.23% |
| schlieren_postprocessor - prepare scratch space            |         1 |   7.6e-05s |         0% |
| time_loop - setup scratch space                            |         1 |     0.127s |         0% |
| time_loop - stalled output                                 |       200 |  0.000685s |         0% |
| time_step - 1 compute d_ij                                 |      7553 |       240s |        67% |
| time_step - 2 compute d_ii, and tau_max                    |      7553 |      11.5s |       3.2% |
| time_step - 3 perform update                               |      7553 |       101s |        28% |
| time_step - 4 fix boundary states                          |      7553 |     0.724s |       0.2% |
| time_step - prepare scratch space                          |         1 |   0.00245s |         0% |
+------------------------------------------------------------+-----------+------------+------------+


One thing that becomes evident is the fact that the program spends two thirds of the execution time computing the graph viscosity d_ij and about a third of the execution time in performing the update, where computing the flux $$f(U)$$ is the expensive operation. The preset default resolution is about 37k gridpoints, which amounts to about 148k spatial degrees of freedom in 2D. An animated schlieren plot of the solution looks as follows:

It is evident that 37k gridpoints for the first-order method is nowhere near the resolution needed to resolve any flow features. For comparison, here is a "reference" computation with a second-order method and about 9.5M gridpoints (github project page):

So, we give the first-order method a second chance and run it with about 2.4M gridpoints on a small compute server:

# mpirun -np 16 ./step-69 | tee output

[...]

####################################################
#########                                  #########
#########     Cycle  070216  (100.0%)      #########
#########      at time t = 3.99999231      #########
#########                                  #########
####################################################

TimeLoop<dim>::output(t = 4.00006, checkpoint = 1)

[...]

+------------------------------------------------------------------------+------------+------------+
| Total wallclock time elapsed since start                               |  6.75e+03s |            |
|                                                                        |            |            |
| Section                                                    | no. calls |  wall time | % of total |
+------------------------------------------------------------+-----------+------------+------------+
| discretization - setup                                     |         1 |      1.97s |         0% |
| offline_data - assemble lumped mass matrix, and c_ij       |         1 |      1.19s |         0% |
| offline_data - compute |c_ij|, and n_ij                    |         1 |    0.0172s |         0% |
| offline_data - create sparsity pattern and set up matrices |         1 |     0.413s |         0% |
| offline_data - distribute dofs                             |         1 |      1.05s |         0% |
| offline_data - fix slip boundary c_ij                      |         1 |     0.252s |         0% |
| schlieren_postprocessor - compute schlieren plot           |       201 |      1.82s |         0% |
| schlieren_postprocessor - prepare scratch space            |         1 |  0.000497s |         0% |
| time_loop - setup scratch space                            |         1 |      1.45s |         0% |
| time_loop - stalled output                                 |       200 |   0.00342s |         0% |
| time_step - 1 compute d_ij                                 |     70216 |  4.38e+03s |        65% |
| time_step - 2 compute d_ii, and tau_max                    |     70216 |       419s |       6.2% |
| time_step - 3 perform update                               |     70216 |  1.87e+03s |        28% |
| time_step - 4 fix boundary states                          |     70216 |        24s |      0.36% |
| time_step - prepare scratch space                          |         1 |    0.0227s |         0% |
+------------------------------------------------------------+-----------+------------+------------+


And with the following result:

That's substantially better, although of course at the price of having run the code for roughly 2 hours on 16 cores.

### Possibilities for extensions

The program showcased here is really only first-order accurate, as discussed above. The pictures above illustrate how much diffusion that introduces and how far the solution is from one that actually resolves the features we care about.

This can be fixed, but it would exceed what a tutorial is about. Nevertheless, it is worth showing what one can achieve by adding a second-order scheme. For example, here is a video computed with the following research code that shows (with a different color scheme) a 2d simulation that corresponds to the cases shown above:

This simulation was done with 38 million degrees of freedom (continuous $$Q_1$$ finite elements) per component of the solution vector. The exquisite detail of the solution is remarkable for these kinds of simulations, including in the sub-sonic region behind the obstacle.

One can also with relative ease further extend this to the 3d case:

Solving this becomes expensive, however: The simulation was done with 1,817 million degrees of freedom (continuous $$Q_1$$ finite elements) per component (for a total of 9.09 billion spatial degrees of freedom) and ran on 30,720 MPI ranks. The code achieved an average througput of 969M grid points per second (0.04M gridpoints per second per CPU). The front and back wall show a "Schlieren plot": the magnitude of the gradient of the density on an exponential scale from white (low) to black (high). All other cutplanes and the surface of the obstacle show the magnitude of the vorticity on a white (low) - yellow (medium) - red (high) scale. (The scales of the individual cutplanes have been adjusted for a nicer visualization.)

# The plain program

/* ---------------------------------------------------------------------
*
* Copyright (C) 2019 - 2022 by the deal.II authors
*
* This file is part of the deal.II library.
*
* The deal.II library is free software; you can use it, redistribute
* it, and/or modify it under the terms of the GNU Lesser General
* version 2.1 of the License, or (at your option) any later version.
* The full text of the license can be found in the file LICENSE.md at
* the top level directory of deal.II.
*
* ---------------------------------------------------------------------
*
* Authors: Matthias Maier, Texas A&M University;
* Ignacio Tomas, Texas A&M University, Sandia National Laboratories
*
* Sandia National Laboratories is a multimission laboratory managed and
* operated by National Technology & Engineering Solutions of Sandia, LLC, a
* wholly owned subsidiary of Honeywell International Inc., for the U.S.
* Department of Energy's National Nuclear Security Administration under
* contract DE-NA0003525. This document describes objective technical results
* and analysis. Any subjective views or opinions that might be expressed in
* the paper do not necessarily represent the views of the U.S. Department of
* Energy or the United States Government.
*/
#include <deal.II/fe/fe.h>
#include <deal.II/lac/sparse_matrix.templates.h>
#include <boost/archive/binary_iarchive.hpp>
#include <boost/archive/binary_oarchive.hpp>
#include <boost/range/iterator_range.hpp>
#include <cmath>
#include <fstream>
#include <future>
namespace Step69
{
using namespace dealii;
namespace Boundaries
{
constexpr types::boundary_id do_nothing = 0;
constexpr types::boundary_id free_slip = 1;
constexpr types::boundary_id dirichlet = 2;
} // namespace Boundaries
template <int dim>
class Discretization : public ParameterAcceptor
{
public:
Discretization(const MPI_Comm mpi_communicator,
TimerOutput & computing_timer,
const std::string &subsection = "Discretization");
void setup();
const MPI_Comm mpi_communicator;
const MappingQ<dim> mapping;
const FE_Q<dim> finite_element;
unsigned int refinement;
private:
TimerOutput &computing_timer;
double length;
double height;
double disk_position;
double disk_diameter;
};
template <int dim>
class OfflineData : public ParameterAcceptor
{
public:
using BoundaryNormalMap =
std::tuple<Tensor<1, dim>, types::boundary_id, Point<dim>>>;
OfflineData(const MPI_Comm mpi_communicator,
TimerOutput & computing_timer,
const Discretization<dim> &discretization,
const std::string & subsection = "OfflineData");
void setup();
void assemble();
DoFHandler<dim> dof_handler;
std::shared_ptr<const Utilities::MPI::Partitioner> partitioner;
unsigned int n_locally_owned;
unsigned int n_locally_relevant;
SparsityPattern sparsity_pattern;
BoundaryNormalMap boundary_normal_map;
SparseMatrix<double> lumped_mass_matrix;
std::array<SparseMatrix<double>, dim> cij_matrix;
std::array<SparseMatrix<double>, dim> nij_matrix;
SparseMatrix<double> norm_matrix;
private:
const MPI_Comm mpi_communicator;
TimerOutput & computing_timer;
};
template <int dim>
class ProblemDescription
{
public:
static constexpr unsigned int problem_dimension = 2 + dim;
using state_type = Tensor<1, problem_dimension>;
const static std::array<std::string, problem_dimension> component_names;
static constexpr double gamma = 7. / 5.;
momentum(const state_type &U);
static DEAL_II_ALWAYS_INLINE inline double
internal_energy(const state_type &U);
static DEAL_II_ALWAYS_INLINE inline double pressure(const state_type &U);
static DEAL_II_ALWAYS_INLINE inline double
speed_of_sound(const state_type &U);
static DEAL_II_ALWAYS_INLINE inline flux_type flux(const state_type &U);
static DEAL_II_ALWAYS_INLINE inline double
compute_lambda_max(const state_type & U_i,
const state_type & U_j,
const Tensor<1, dim> &n_ij);
};
template <int dim>
class InitialValues : public ParameterAcceptor
{
public:
using state_type = typename ProblemDescription<dim>::state_type;
InitialValues(const std::string &subsection = "InitialValues");
std::function<state_type(const Point<dim> &point, double t)> initial_state;
private:
void parse_parameters_callback();
Tensor<1, dim> initial_direction;
Tensor<1, 3> initial_1d_state;
};
template <int dim>
{
public:
static constexpr unsigned int problem_dimension =
ProblemDescription<dim>::problem_dimension;
using state_type = typename ProblemDescription<dim>::state_type;
using flux_type = typename ProblemDescription<dim>::flux_type;
using vector_type =
std::array<LinearAlgebra::distributed::Vector<double>, problem_dimension>;
TimeStepping(const MPI_Comm mpi_communicator,
TimerOutput & computing_timer,
const OfflineData<dim> & offline_data,
const InitialValues<dim> &initial_values,
const std::string & subsection = "TimeStepping");
void prepare();
double make_one_step(vector_type &U, double t);
private:
const MPI_Comm mpi_communicator;
TimerOutput & computing_timer;
vector_type temporary_vector;
double cfl_update;
};
template <int dim>
class SchlierenPostprocessor : public ParameterAcceptor
{
public:
static constexpr unsigned int problem_dimension =
ProblemDescription<dim>::problem_dimension;
using state_type = typename ProblemDescription<dim>::state_type;
using vector_type =
std::array<LinearAlgebra::distributed::Vector<double>, problem_dimension>;
SchlierenPostprocessor(
const MPI_Comm mpi_communicator,
TimerOutput & computing_timer,
const OfflineData<dim> &offline_data,
const std::string & subsection = "SchlierenPostprocessor");
void prepare();
void compute_schlieren(const vector_type &U);
private:
const MPI_Comm mpi_communicator;
TimerOutput & computing_timer;
unsigned int schlieren_index;
double schlieren_beta;
};
template <int dim>
class MainLoop : public ParameterAcceptor
{
public:
using vector_type = typename TimeStepping<dim>::vector_type;
MainLoop(const MPI_Comm mpi_communnicator);
void run();
private:
vector_type interpolate_initial_values(const double t = 0);
void checkpoint(const vector_type &U,
const std::string &name,
double t,
unsigned int cycle);
void output(const vector_type &U,
const std::string &name,
double t,
unsigned int cycle);
const MPI_Comm mpi_communicator;
std::ostringstream timer_output;
TimerOutput computing_timer;
std::string base_name;
double t_final;
double output_granularity;
bool asynchronous_writeback;
bool resume;
Discretization<dim> discretization;
OfflineData<dim> offline_data;
InitialValues<dim> initial_values;
TimeStepping<dim> time_stepping;
SchlierenPostprocessor<dim> schlieren_postprocessor;
vector_type output_vector;
};
template <int dim>
Discretization<dim>::Discretization(const MPI_Comm mpi_communicator,
TimerOutput & computing_timer,
const std::string &subsection)
: ParameterAcceptor(subsection)
, mpi_communicator(mpi_communicator)
, triangulation(mpi_communicator)
, mapping(1)
, finite_element(1)
, computing_timer(computing_timer)
{
length = 4.;
add_parameter("length", length, "Length of computational domain");
height = 2.;
add_parameter("height", height, "Height of computational domain");
disk_position = 0.6;
disk_position,
"x position of immersed disk center point");
disk_diameter = 0.5;
disk_diameter,
"Diameter of immersed disk");
refinement = 5;
refinement,
"Number of refinement steps of the geometry");
}
template <int dim>
void Discretization<dim>::setup()
{
TimerOutput::Scope scope(computing_timer, "discretization - setup");
triangulation.clear();
Triangulation<dim> tria1, tria2, tria3, tria4, tria5, tria6;
tria1, disk_diameter / 2., disk_diameter, 0.5, 1, false);
tria2,
{2, 1},
Point<2>(-disk_diameter, disk_diameter),
Point<2>(disk_diameter, height / 2.));
tria3,
{2, 1},
Point<2>(-disk_diameter, -disk_diameter),
Point<2>(disk_diameter, -height / 2.));
tria4,
{6, 2},
Point<2>(disk_diameter, -disk_diameter),
Point<2>(length - disk_position, disk_diameter));
tria5,
{6, 1},
Point<2>(disk_diameter, disk_diameter),
Point<2>(length - disk_position, height / 2.));
tria6,
{6, 1},
Point<2>(disk_diameter, -height / 2.),
Point<2>(length - disk_position, -disk_diameter));
{&tria1, &tria2, &tria3, &tria4, &tria5, &tria6},
1.e-12,
true);
for (const auto &cell : triangulation.active_cell_iterators())
{
for (const auto v : cell->vertex_indices())
{
if (cell->vertex(v)[0] <= -disk_diameter + 1.e-6)
cell->vertex(v)[0] = -disk_position;
}
}
for (const auto &cell : triangulation.active_cell_iterators())
{
for (const auto f : cell->face_indices())
{
const auto face = cell->face(f);
if (face->at_boundary())
{
const auto center = face->center();
if (center[0] > length - disk_position - 1.e-6)
face->set_boundary_id(Boundaries::do_nothing);
else if (center[0] < -disk_position + 1.e-6)
face->set_boundary_id(Boundaries::dirichlet);
else
face->set_boundary_id(Boundaries::free_slip);
}
}
}
}
template <int dim>
OfflineData<dim>::OfflineData(const MPI_Comm mpi_communicator,
TimerOutput & computing_timer,
const Discretization<dim> &discretization,
const std::string & subsection)
: ParameterAcceptor(subsection)
, dof_handler(discretization.triangulation)
, mpi_communicator(mpi_communicator)
, computing_timer(computing_timer)
, discretization(&discretization)
{}
template <int dim>
void OfflineData<dim>::setup()
{
IndexSet locally_owned;
IndexSet locally_relevant;
{
TimerOutput::Scope scope(computing_timer,
"offline_data - distribute dofs");
dof_handler.distribute_dofs(discretization->finite_element);
locally_owned = dof_handler.locally_owned_dofs();
n_locally_owned = locally_owned.n_elements();
locally_relevant = DoFTools::extract_locally_relevant_dofs(dof_handler);
n_locally_relevant = locally_relevant.n_elements();
partitioner =
std::make_shared<Utilities::MPI::Partitioner>(locally_owned,
locally_relevant,
mpi_communicator);
}
{
computing_timer,
"offline_data - create sparsity pattern and set up matrices");
DynamicSparsityPattern dsp(n_locally_relevant, n_locally_relevant);
const auto dofs_per_cell =
discretization->finite_element.n_dofs_per_cell();
std::vector<types::global_dof_index> dof_indices(dofs_per_cell);
for (const auto &cell : dof_handler.active_cell_iterators())
{
if (cell->is_artificial())
continue;
/* We transform the set of global dof indices on the cell to the
* corresponding "local" index range on the MPI process: */
cell->get_dof_indices(dof_indices);
std::transform(dof_indices.begin(),
dof_indices.end(),
dof_indices.begin(),
return partitioner->global_to_local(index);
});
/* And simply add, for each dof, a coupling to all other "local"
* dofs on the cell: */
for (const auto dof : dof_indices)
}
sparsity_pattern.copy_from(dsp);
lumped_mass_matrix.reinit(sparsity_pattern);
norm_matrix.reinit(sparsity_pattern);
for (auto &matrix : cij_matrix)
matrix.reinit(sparsity_pattern);
for (auto &matrix : nij_matrix)
matrix.reinit(sparsity_pattern);
}
}
namespace
{
template <int dim>
struct CopyData
{
bool is_artificial;
std::vector<types::global_dof_index> local_dof_indices;
typename OfflineData<dim>::BoundaryNormalMap local_boundary_normal_map;
FullMatrix<double> cell_lumped_mass_matrix;
std::array<FullMatrix<double>, dim> cell_cij_matrix;
};
template <typename IteratorType>
get_entry(const SparseMatrix<double> &matrix, const IteratorType &it)
{
const SparseMatrix<double>::const_iterator matrix_iterator(
&matrix, it->global_index());
return matrix_iterator->value();
}
template <typename IteratorType>
const IteratorType & it,
{
it->global_index());
matrix_iterator->value() = value;
}
template <std::size_t k, typename IteratorType>
gather_get_entry(const std::array<SparseMatrix<double>, k> &c_ij,
const IteratorType it)
{
Tensor<1, k> result;
for (unsigned int j = 0; j < k; ++j)
result[j] = get_entry(c_ij[j], it);
return result;
}
template <std::size_t k>
gather(const std::array<SparseMatrix<double>, k> &n_ij,
const unsigned int i,
const unsigned int j)
{
Tensor<1, k> result;
for (unsigned int l = 0; l < k; ++l)
result[l] = n_ij[l](i, j);
return result;
}
template <std::size_t k>
const unsigned int i)
{
Tensor<1, k> result;
for (unsigned int j = 0; j < k; ++j)
result[j] = U[j].local_element(i);
return result;
}
template <std::size_t k, int k2>
const Tensor<1, k2> & tensor,
const unsigned int i)
{
static_assert(k == k2,
"The dimensions of the input arguments must agree");
for (unsigned int j = 0; j < k; ++j)
U[j].local_element(i) = tensor[j];
}
} // namespace
template <int dim>
{
lumped_mass_matrix = 0.;
norm_matrix = 0.;
for (auto &matrix : cij_matrix)
matrix = 0.;
for (auto &matrix : nij_matrix)
matrix = 0.;
unsigned int dofs_per_cell =
discretization->finite_element.n_dofs_per_cell();
discretization->mapping,
discretization->finite_element,
{
computing_timer,
"offline_data - assemble lumped mass matrix, and c_ij");
const auto local_assemble_system =
[&](const typename DoFHandler<dim>::cell_iterator &cell,
CopyData<dim> & copy) {
copy.is_artificial = cell->is_artificial();
if (copy.is_artificial)
return;
copy.local_boundary_normal_map.clear();
copy.cell_lumped_mass_matrix.reinit(dofs_per_cell, dofs_per_cell);
for (auto &matrix : copy.cell_cij_matrix)
matrix.reinit(dofs_per_cell, dofs_per_cell);
const auto &fe_values = scratch.reinit(cell);
copy.local_dof_indices.resize(dofs_per_cell);
cell->get_dof_indices(copy.local_dof_indices);
std::transform(copy.local_dof_indices.begin(),
copy.local_dof_indices.end(),
copy.local_dof_indices.begin(),
return partitioner->global_to_local(index);
});
for (unsigned int q_point = 0; q_point < n_q_points; ++q_point)
{
const auto JxW = fe_values.JxW(q_point);
for (unsigned int j = 0; j < dofs_per_cell; ++j)
{
const auto value_JxW =
fe_values.shape_value(j, q_point) * JxW;
copy.cell_lumped_mass_matrix(j, j) += value_JxW;
for (unsigned int i = 0; i < dofs_per_cell; ++i)
{
const auto value = fe_values.shape_value(i, q_point);
for (unsigned int d = 0; d < dim; ++d)
copy.cell_cij_matrix[d](i, j) += value * grad_JxW[d];
} /* i */
} /* j */
} /* q */
for (const auto f : cell->face_indices())
{
const auto face = cell->face(f);
const auto id = face->boundary_id();
if (!face->at_boundary())
continue;
const auto &fe_face_values = scratch.reinit(cell, f);
const unsigned int n_face_q_points =
for (unsigned int j = 0; j < dofs_per_cell; ++j)
{
if (!discretization->finite_element.has_support_on_face(j, f))
continue;
if (id == Boundaries::free_slip)
{
for (unsigned int q = 0; q < n_face_q_points; ++q)
normal += fe_face_values.normal_vector(q) *
fe_face_values.shape_value(j, q);
}
const auto index = copy.local_dof_indices[j];
Point<dim> position;
for (const auto v : cell->vertex_indices())
if (cell->vertex_dof_index(v, 0) ==
partitioner->local_to_global(index))
{
position = cell->vertex(v);
break;
}
const auto old_id =
std::get<1>(copy.local_boundary_normal_map[index]);
copy.local_boundary_normal_map[index] =
std::make_tuple(normal, std::max(old_id, id), position);
}
}
};
const auto copy_local_to_global = [&](const CopyData<dim> &copy) {
if (copy.is_artificial)
return;
for (const auto &it : copy.local_boundary_normal_map)
{
std::get<0>(boundary_normal_map[it.first]) +=
std::get<0>(it.second);
std::get<1>(boundary_normal_map[it.first]) =
std::max(std::get<1>(boundary_normal_map[it.first]),
std::get<1>(it.second));
std::get<2>(boundary_normal_map[it.first]) = std::get<2>(it.second);
}
copy.cell_lumped_mass_matrix);
for (int k = 0; k < dim; ++k)
{
}
};
WorkStream::run(dof_handler.begin_active(),
dof_handler.end(),
local_assemble_system,
copy_local_to_global,
scratch_data,
CopyData<dim>());
}
{
TimerOutput::Scope scope(computing_timer,
"offline_data - compute |c_ij|, and n_ij");
0, n_locally_relevant);
const auto on_subranges =
[&](const auto i1, const auto i2) {
for (const auto row_index :
*i2))
{
std::for_each(
sparsity_pattern.begin(row_index),
sparsity_pattern.end(row_index),
[&](const ::SparsityPatternIterators::Accessor &jt) {
const auto c_ij = gather_get_entry(cij_matrix, &jt);
const double norm = c_ij.norm();
set_entry(norm_matrix, &jt, norm);
for (unsigned int j = 0; j < dim; ++j)
set_entry(nij_matrix[j], &jt, c_ij[j] / norm);
});
}
};
indices.end(),
on_subranges,
4096);
for (auto &it : boundary_normal_map)
{
auto &normal = std::get<0>(it.second);
normal /= (normal.norm() + std::numeric_limits<double>::epsilon());
}
}
}
template <int dim>
ProblemDescription<dim>::momentum(const state_type &U)
{
std::copy_n(&U[1], dim, &result[0]);
return result;
}
template <int dim>
DEAL_II_ALWAYS_INLINE inline double
ProblemDescription<dim>::internal_energy(const state_type &U)
{
const double &rho = U[0];
const auto m = momentum(U);
const double &E = U[dim + 1];
return E - 0.5 * m.norm_square() / rho;
}
template <int dim>
DEAL_II_ALWAYS_INLINE inline double
ProblemDescription<dim>::pressure(const state_type &U)
{
return (gamma - 1.) * internal_energy(U);
}
template <int dim>
DEAL_II_ALWAYS_INLINE inline double
ProblemDescription<dim>::speed_of_sound(const state_type &U)
{
const double &rho = U[0];
const double p = pressure(U);
return std::sqrt(gamma * p / rho);
}
template <int dim>
DEAL_II_ALWAYS_INLINE inline typename ProblemDescription<dim>::flux_type
ProblemDescription<dim>::flux(const state_type &U)
{
const double &rho = U[0];
const auto m = momentum(U);
const auto p = pressure(U);
const double &E = U[dim + 1];
flux_type result;
result[0] = m;
for (unsigned int i = 0; i < dim; ++i)
{
result[1 + i] = m * m[i] / rho;
result[1 + i][i] += p;
}
result[dim + 1] = m / rho * (E + p);
return result;
}
namespace
{
template <int dim>
DEAL_II_ALWAYS_INLINE inline std::array<double, 4> riemann_data_from_state(
const typename ProblemDescription<dim>::state_type U,
const Tensor<1, dim> & n_ij)
{
Tensor<1, 3> projected_U;
projected_U[0] = U[0];
const auto m = ProblemDescription<dim>::momentum(U);
projected_U[1] = n_ij * m;
const auto perpendicular_m = m - projected_U[1] * n_ij;
projected_U[2] = U[1 + dim] - 0.5 * perpendicular_m.norm_square() / U[0];
return {{projected_U[0],
projected_U[1] / projected_U[0],
ProblemDescription<1>::pressure(projected_U),
ProblemDescription<1>::speed_of_sound(projected_U)}};
}
DEAL_II_ALWAYS_INLINE inline double positive_part(const double number)
{
return std::max(number, 0.);
}
DEAL_II_ALWAYS_INLINE inline double negative_part(const double number)
{
return -std::min(number, 0.);
}
DEAL_II_ALWAYS_INLINE inline double
lambda1_minus(const std::array<double, 4> &riemann_data,
const double p_star)
{
/* Implements formula (3.7) in Guermond-Popov-2016 */
constexpr double gamma = ProblemDescription<1>::gamma;
const auto u = riemann_data[1];
const auto p = riemann_data[2];
const auto a = riemann_data[3];
const double factor = (gamma + 1.0) / 2.0 / gamma;
const double tmp = positive_part((p_star - p) / p);
return u - a * std::sqrt(1.0 + factor * tmp);
}
DEAL_II_ALWAYS_INLINE inline double
lambda3_plus(const std::array<double, 4> &riemann_data, const double p_star)
{
/* Implements formula (3.8) in Guermond-Popov-2016 */
constexpr double gamma = ProblemDescription<1>::gamma;
const auto u = riemann_data[1];
const auto p = riemann_data[2];
const auto a = riemann_data[3];
const double factor = (gamma + 1.0) / 2.0 / gamma;
const double tmp = positive_part((p_star - p) / p);
return u + a * std::sqrt(1.0 + factor * tmp);
}
DEAL_II_ALWAYS_INLINE inline double
lambda_max_two_rarefaction(const std::array<double, 4> &riemann_data_i,
const std::array<double, 4> &riemann_data_j)
{
constexpr double gamma = ProblemDescription<1>::gamma;
const auto u_i = riemann_data_i[1];
const auto p_i = riemann_data_i[2];
const auto a_i = riemann_data_i[3];
const auto u_j = riemann_data_j[1];
const auto p_j = riemann_data_j[2];
const auto a_j = riemann_data_j[3];
const double numerator = a_i + a_j - (gamma - 1.) / 2. * (u_j - u_i);
const double denominator =
a_i * std::pow(p_i / p_j, -1. * (gamma - 1.) / 2. / gamma) + a_j * 1.;
/* Formula (4.3) in Guermond-Popov-2016 */
const double p_star =
p_j * std::pow(numerator / denominator, 2. * gamma / (gamma - 1));
const double lambda1 = lambda1_minus(riemann_data_i, p_star);
const double lambda3 = lambda3_plus(riemann_data_j, p_star);
/* Formula (2.11) in Guermond-Popov-2016 */
return std::max(positive_part(lambda3), negative_part(lambda1));
}
DEAL_II_ALWAYS_INLINE inline double
lambda_max_expansion(const std::array<double, 4> &riemann_data_i,
const std::array<double, 4> &riemann_data_j)
{
const auto u_i = riemann_data_i[1];
const auto a_i = riemann_data_i[3];
const auto u_j = riemann_data_j[1];
const auto a_j = riemann_data_j[3];
return std::max(std::abs(u_i), std::abs(u_j)) + 5. * std::max(a_i, a_j);
}
} // namespace
template <int dim>
DEAL_II_ALWAYS_INLINE inline double
ProblemDescription<dim>::compute_lambda_max(const state_type & U_i,
const state_type & U_j,
const Tensor<1, dim> &n_ij)
{
const auto riemann_data_i = riemann_data_from_state(U_i, n_ij);
const auto riemann_data_j = riemann_data_from_state(U_j, n_ij);
const double lambda_1 =
lambda_max_two_rarefaction(riemann_data_i, riemann_data_j);
const double lambda_2 =
lambda_max_expansion(riemann_data_i, riemann_data_j);
return std::min(lambda_1, lambda_2);
}
template <>
const std::array<std::string, 3> ProblemDescription<1>::component_names{
{"rho", "m", "E"}};
template <>
const std::array<std::string, 4> ProblemDescription<2>::component_names{
{"rho", "m_1", "m_2", "E"}};
template <>
const std::array<std::string, 5> ProblemDescription<3>::component_names{
{"rho", "m_1", "m_2", "m_3", "E"}};
template <int dim>
InitialValues<dim>::InitialValues(const std::string &subsection)
: ParameterAcceptor(subsection)
{
/* We wire up the slot InitialValues<dim>::parse_parameters_callback to
the ParameterAcceptor::parse_parameters_call_back signal: */
std::bind(&InitialValues<dim>::parse_parameters_callback, this));
initial_direction[0] = 1.;
initial_direction,
"Initial direction of the uniform flow field");
initial_1d_state[0] = ProblemDescription<dim>::gamma;
initial_1d_state[1] = 3.;
initial_1d_state[2] = 1.;
initial_1d_state,
"Initial 1d state (rho, u, p) of the uniform flow field");
}
template <int dim>
void InitialValues<dim>::parse_parameters_callback()
{
AssertThrow(initial_direction.norm() != 0.,
"Initial shock front direction is set to the zero vector."));
initial_direction /= initial_direction.norm();
initial_state = [this](const Point<dim> & /*point*/, double /*t*/) {
const double rho = initial_1d_state[0];
const double u = initial_1d_state[1];
const double p = initial_1d_state[2];
static constexpr double gamma = ProblemDescription<dim>::gamma;
state_type state;
state[0] = rho;
for (unsigned int i = 0; i < dim; ++i)
state[1 + i] = rho * u * initial_direction[i];
state[dim + 1] = p / (gamma - 1.) + 0.5 * rho * u * u;
return state;
};
}
template <int dim>
const MPI_Comm mpi_communicator,
TimerOutput & computing_timer,
const OfflineData<dim> & offline_data,
const InitialValues<dim> &initial_values,
const std::string & subsection /*= "TimeStepping"*/)
: ParameterAcceptor(subsection)
, mpi_communicator(mpi_communicator)
, computing_timer(computing_timer)
, offline_data(&offline_data)
, initial_values(&initial_values)
{
cfl_update = 0.80;
cfl_update,
"Relative CFL constant used for update");
}
template <int dim>
void TimeStepping<dim>::prepare()
{
TimerOutput::Scope scope(computing_timer,
"time_stepping - prepare scratch space");
for (auto &it : temporary_vector)
it.reinit(offline_data->partitioner);
dij_matrix.reinit(offline_data->sparsity_pattern);
}
template <int dim>
double TimeStepping<dim>::make_one_step(vector_type &U, double t)
{
const auto &n_locally_owned = offline_data->n_locally_owned;
const auto &n_locally_relevant = offline_data->n_locally_relevant;
indices_owned(0, n_locally_owned);
indices_relevant(0, n_locally_relevant);
const auto &sparsity = offline_data->sparsity_pattern;
const auto &lumped_mass_matrix = offline_data->lumped_mass_matrix;
const auto &norm_matrix = offline_data->norm_matrix;
const auto &nij_matrix = offline_data->nij_matrix;
const auto &cij_matrix = offline_data->cij_matrix;
const auto &boundary_normal_map = offline_data->boundary_normal_map;
{
TimerOutput::Scope scope(computing_timer,
"time_stepping - 1 compute d_ij");
const auto on_subranges =
[&](const auto i1, const auto i2) {
for (const auto i :
*i2))
{
const auto U_i = gather(U, i);
for (auto jt = sparsity.begin(i); jt != sparsity.end(i); ++jt)
{
const auto j = jt->column();
if (j >= i)
continue;
const auto U_j = gather(U, j);
const auto n_ij = gather_get_entry(nij_matrix, jt);
const double norm = get_entry(norm_matrix, jt);
const auto lambda_max =
ProblemDescription<dim>::compute_lambda_max(U_i, U_j, n_ij);
double d = norm * lambda_max;
if (boundary_normal_map.count(i) != 0 &&
boundary_normal_map.count(j) != 0)
{
const auto n_ji = gather(nij_matrix, j, i);
const auto lambda_max_2 =
ProblemDescription<dim>::compute_lambda_max(U_j,
U_i,
n_ji);
const double norm_2 = norm_matrix(j, i);
d = std::max(d, norm_2 * lambda_max_2);
}
set_entry(dij_matrix, jt, d);
dij_matrix(j, i) = d;
}
}
};
parallel::apply_to_subranges(indices_relevant.begin(),
indices_relevant.end(),
on_subranges,
4096);
}
std::atomic<double> tau_max{std::numeric_limits<double>::infinity()};
{
TimerOutput::Scope scope(computing_timer,
"time_stepping - 2 compute d_ii, and tau_max");
const auto on_subranges =
[&](const auto i1, const auto i2) {
double tau_max_on_subrange = std::numeric_limits<double>::infinity();
for (const auto i :
*i2))
{
double d_sum = 0.;
for (auto jt = sparsity.begin(i); jt != sparsity.end(i); ++jt)
{
const auto j = jt->column();
if (j == i)
continue;
d_sum -= get_entry(dij_matrix, jt);
}
dij_matrix.diag_element(i) = d_sum;
const double mass = lumped_mass_matrix.diag_element(i);
const double tau = cfl_update * mass / (-2. * d_sum);
tau_max_on_subrange = std::min(tau_max_on_subrange, tau);
}
while (current_tau_max > tau_max_on_subrange &&
!tau_max.compare_exchange_weak(current_tau_max,
tau_max_on_subrange))
;
};
parallel::apply_to_subranges(indices_relevant.begin(),
indices_relevant.end(),
on_subranges,
4096);
"I'm sorry, Dave. I'm afraid I can't do that. - We crashed."));
}
{
TimerOutput::Scope scope(computing_timer,
"time_stepping - 3 perform update");
const auto on_subranges =
[&](const auto i1, const auto i2) {
for (const auto i : boost::make_iterator_range(i1, i2))
{
Assert(i < n_locally_owned, ExcInternalError());
const auto U_i = gather(U, i);
const auto f_i = ProblemDescription<dim>::flux(U_i);
const double m_i = lumped_mass_matrix.diag_element(i);
auto U_i_new = U_i;
for (auto jt = sparsity.begin(i); jt != sparsity.end(i); ++jt)
{
const auto j = jt->column();
const auto U_j = gather(U, j);
const auto f_j = ProblemDescription<dim>::flux(U_j);
const auto c_ij = gather_get_entry(cij_matrix, jt);
const auto d_ij = get_entry(dij_matrix, jt);
for (unsigned int k = 0; k < problem_dimension; ++k)
{
U_i_new[k] +=
tau_max / m_i *
(-(f_j[k] - f_i[k]) * c_ij + d_ij * (U_j[k] - U_i[k]));
}
}
scatter(temporary_vector, U_i_new, i);
}
};
parallel::apply_to_subranges(indices_owned.begin(),
indices_owned.end(),
on_subranges,
4096);
}
{
TimerOutput::Scope scope(computing_timer,
"time_stepping - 4 fix boundary states");
for (auto it : boundary_normal_map)
{
const auto i = it.first;
if (i >= n_locally_owned)
continue;
const auto &normal = std::get<0>(it.second);
const auto &id = std::get<1>(it.second);
const auto &position = std::get<2>(it.second);
auto U_i = gather(temporary_vector, i);
if (id == Boundaries::free_slip)
{
auto m = ProblemDescription<dim>::momentum(U_i);
m -= (m * normal) * normal;
for (unsigned int k = 0; k < dim; ++k)
U_i[k + 1] = m[k];
}
else if (id == Boundaries::dirichlet)
{
U_i = initial_values->initial_state(position, t + tau_max);
}
scatter(temporary_vector, U_i, i);
}
}
for (auto &it : temporary_vector)
it.update_ghost_values();
U.swap(temporary_vector);
return tau_max;
}
template <int dim>
SchlierenPostprocessor<dim>::SchlierenPostprocessor(
const MPI_Comm mpi_communicator,
TimerOutput & computing_timer,
const OfflineData<dim> &offline_data,
const std::string & subsection /*= "SchlierenPostprocessor"*/)
: ParameterAcceptor(subsection)
, mpi_communicator(mpi_communicator)
, computing_timer(computing_timer)
, offline_data(&offline_data)
{
schlieren_beta = 10.;
schlieren_beta,
"Beta factor used in Schlieren-type postprocessor");
schlieren_index = 0;
schlieren_index,
"Use the corresponding component of the state vector for the "
"schlieren plot");
}
template <int dim>
void SchlierenPostprocessor<dim>::prepare()
{
TimerOutput::Scope scope(computing_timer,
"schlieren_postprocessor - prepare scratch space");
r.reinit(offline_data->n_locally_relevant);
schlieren.reinit(offline_data->partitioner);
}
template <int dim>
void SchlierenPostprocessor<dim>::compute_schlieren(const vector_type &U)
{
computing_timer, "schlieren_postprocessor - compute schlieren plot");
const auto &sparsity = offline_data->sparsity_pattern;
const auto &lumped_mass_matrix = offline_data->lumped_mass_matrix;
const auto &cij_matrix = offline_data->cij_matrix;
const auto &boundary_normal_map = offline_data->boundary_normal_map;
const auto &n_locally_owned = offline_data->n_locally_owned;
const auto indices =
n_locally_owned);
std::atomic<double> r_i_max{0.};
std::atomic<double> r_i_min{std::numeric_limits<double>::infinity()};
{
const auto on_subranges =
[&](const auto i1, const auto i2) {
double r_i_max_on_subrange = 0.;
double r_i_min_on_subrange = std::numeric_limits<double>::infinity();
for (const auto i : boost::make_iterator_range(i1, i2))
{
Assert(i < n_locally_owned, ExcInternalError());
for (auto jt = sparsity.begin(i); jt != sparsity.end(i); ++jt)
{
const auto j = jt->column();
if (i == j)
continue;
const auto U_js = U[schlieren_index].local_element(j);
const auto c_ij = gather_get_entry(cij_matrix, jt);
r_i += c_ij * U_js;
}
const auto bnm_it = boundary_normal_map.find(i);
if (bnm_it != boundary_normal_map.end())
{
const auto &normal = std::get<0>(bnm_it->second);
const auto &id = std::get<1>(bnm_it->second);
if (id == Boundaries::free_slip)
r_i -= 1. * (r_i * normal) * normal;
else
r_i = 0.;
}
const double m_i = lumped_mass_matrix.diag_element(i);
r[i] = r_i.norm() / m_i;
r_i_max_on_subrange = std::max(r_i_max_on_subrange, r[i]);
r_i_min_on_subrange = std::min(r_i_min_on_subrange, r[i]);
}
while (current_r_i_max < r_i_max_on_subrange &&
!r_i_max.compare_exchange_weak(current_r_i_max,
r_i_max_on_subrange))
;
while (current_r_i_min > r_i_min_on_subrange &&
!r_i_min.compare_exchange_weak(current_r_i_min,
r_i_min_on_subrange))
;
};
indices.end(),
on_subranges,
4096);
}
{
const auto on_subranges =
[&](const auto i1, const auto i2) {
for (const auto i : boost::make_iterator_range(i1, i2))
{
Assert(i < n_locally_owned, ExcInternalError());
schlieren.local_element(i) =
1. - std::exp(-schlieren_beta * (r[i] - r_i_min) /
(r_i_max - r_i_min));
}
};
indices.end(),
on_subranges,
4096);
}
schlieren.update_ghost_values();
}
template <int dim>
MainLoop<dim>::MainLoop(const MPI_Comm mpi_communicator)
: ParameterAcceptor("A - MainLoop")
, mpi_communicator(mpi_communicator)
, computing_timer(mpi_communicator,
timer_output,
TimerOutput::never,
TimerOutput::cpu_and_wall_times)
, pcout(std::cout, Utilities::MPI::this_mpi_process(mpi_communicator) == 0)
, discretization(mpi_communicator, computing_timer, "B - Discretization")
, offline_data(mpi_communicator,
computing_timer,
discretization,
"C - OfflineData")
, initial_values("D - InitialValues")
, time_stepping(mpi_communicator,
computing_timer,
offline_data,
initial_values,
"E - TimeStepping")
, schlieren_postprocessor(mpi_communicator,
computing_timer,
offline_data,
"F - SchlierenPostprocessor")
{
base_name = "test";
add_parameter("basename", base_name, "Base name for all output files");
t_final = 4.;
output_granularity = 0.02;
output_granularity,
"time interval for output");
asynchronous_writeback = true;
asynchronous_writeback,
"Write out solution in a background thread performing IO");
resume = false;
add_parameter("resume", resume, "Resume an interrupted computation.");
}
namespace
{
const std::string & secondary = "")
{
std::string((35 - header_size) / 2, ' ');
const auto secondary_size = secondary.size();
std::string((34 - secondary_size) / 2, ' ') + secondary +
std::string((35 - secondary_size) / 2, ' ');
/* clang-format off */
pcout << std::endl;
pcout << " ####################################################" << std::endl;
pcout << " ######### #########" << std::endl;
pcout << " #########" << padded_secondary << "#########" << std::endl;
pcout << " ######### #########" << std::endl;
pcout << " ####################################################" << std::endl;
pcout << std::endl;
/* clang-format on */
}
} // namespace
template <int dim>
{
pcout << "Reading parameters and allocating objects... " << std::flush;
pcout << "done" << std::endl;
{
discretization.setup();
if (resume)
else
discretization.triangulation.refine_global(discretization.refinement);
pcout << "Number of active cells: "
<< discretization.triangulation.n_global_active_cells()
<< std::endl;
offline_data.setup();
offline_data.assemble();
pcout << "Number of degrees of freedom: "
<< offline_data.dof_handler.n_dofs() << std::endl;
time_stepping.prepare();
schlieren_postprocessor.prepare();
}
double t = 0.;
unsigned int output_cycle = 0;
vector_type U;
for (auto &it : U)
it.reinit(offline_data.partitioner);
if (resume)
{
parallel::distributed::
SolutionTransfer<dim, LinearAlgebra::distributed::Vector<double>>
solution_transfer(offline_data.dof_handler);
std::vector<LinearAlgebra::distributed::Vector<double> *> vectors;
std::transform(U.begin(),
U.end(),
std::back_inserter(vectors),
[](auto &it) { return &it; });
solution_transfer.deserialize(vectors);
for (auto &it : U)
it.update_ghost_values();
std::ios::binary);
boost::archive::binary_iarchive ia(file);
ia >> t >> output_cycle;
}
else
{
U = interpolate_initial_values();
}
output(U, base_name, t, output_cycle++);
for (unsigned int cycle = 1; t < t_final; ++cycle)
{
std::ostringstream secondary;
head << "Cycle " << Utilities::int_to_string(cycle, 6) << " ("
<< std::fixed << std::setprecision(1) << t / t_final * 100
<< "%)";
secondary << "at time t = " << std::setprecision(8) << std::fixed << t;
t += time_stepping.make_one_step(U, t);
if (t > output_cycle * output_granularity)
{
checkpoint(U, base_name, t, output_cycle);
output(U, base_name, t, output_cycle);
++output_cycle;
}
}
computing_timer.print_summary();
pcout << timer_output.str() << std::endl;
}
template <int dim>
typename MainLoop<dim>::vector_type
MainLoop<dim>::interpolate_initial_values(const double t)
{
pcout << "MainLoop<dim>::interpolate_initial_values(t = " << t << ')'
<< std::endl;
TimerOutput::Scope scope(computing_timer,
"main_loop - setup scratch space");
vector_type U;
for (auto &it : U)
it.reinit(offline_data.partitioner);
constexpr auto problem_dimension =
ProblemDescription<dim>::problem_dimension;
for (unsigned int i = 0; i < problem_dimension; ++i)
VectorTools::interpolate(offline_data.dof_handler,
[&](const Point<dim> &x) {
return initial_values.initial_state(x, t)[i];
}),
U[i]);
for (auto &it : U)
it.update_ghost_values();
return U;
}
template <int dim>
void MainLoop<dim>::checkpoint(const typename MainLoop<dim>::vector_type &U,
const std::string &name,
const double t,
const unsigned int cycle)
{
parallel::distributed::
SolutionTransfer<dim, LinearAlgebra::distributed::Vector<double>>
solution_transfer(offline_data.dof_handler);
std::vector<const LinearAlgebra::distributed::Vector<double> *> vectors;
std::transform(U.begin(),
U.end(),
std::back_inserter(vectors),
[](auto &it) { return &it; });
solution_transfer.prepare_for_serialization(vectors);
discretization.triangulation.save(name + "-checkpoint.mesh");
if (::Utilities::MPI::this_mpi_process(mpi_communicator) == 0)
{
boost::archive::binary_oarchive oa(file);
oa << t << cycle;
}
}
template <int dim>
void MainLoop<dim>::output(const typename MainLoop<dim>::vector_type &U,
const std::string & name,
const double t,
const unsigned int cycle)
{
pcout << "MainLoop<dim>::output(t = " << t << ')' << std::endl;
{
TimerOutput::Scope timer(computing_timer, "main_loop - stalled output");
}
constexpr auto problem_dimension =
ProblemDescription<dim>::problem_dimension;
for (unsigned int i = 0; i < problem_dimension; ++i)
{
output_vector[i] = U[i];
output_vector[i].update_ghost_values();
}
schlieren_postprocessor.compute_schlieren(output_vector);
auto data_out = std::make_shared<DataOut<dim>>();
data_out->attach_dof_handler(offline_data.dof_handler);
const auto &component_names = ProblemDescription<dim>::component_names;
for (unsigned int i = 0; i < problem_dimension; ++i)
"schlieren_plot");
data_out->build_patches(discretization.mapping,
discretization.finite_element.degree - 1);
const auto output_worker = [this, name, t, cycle, data_out]() {
cycle,
true,
data_out->set_flags(flags);
data_out->write_vtu_with_pvtu_record(
"", name + "-solution", cycle, mpi_communicator, 6);
};
if (asynchronous_writeback)
{
}
else
{
output_worker();
}
}
} // namespace Step69
int main(int argc, char *argv[])
{
try
{
constexpr int dim = 2;
using namespace dealii;
using namespace Step69;
Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv);
MPI_Comm mpi_communicator(MPI_COMM_WORLD);
MainLoop<dim> main_loop(mpi_communicator);
main_loop.run();
}
catch (std::exception &exc)
{
std::cerr << std::endl
<< std::endl
<< "----------------------------------------------------"
<< std::endl;
std::cerr << "Exception on processing: " << std::endl
<< exc.what() << std::endl
<< "Aborting!" << std::endl
<< "----------------------------------------------------"
<< std::endl;
return 1;
}
catch (...)
{
std::cerr << std::endl
<< std::endl
<< "----------------------------------------------------"
<< std::endl;
std::cerr << "Unknown exception!" << std::endl
<< "Aborting!" << std::endl
<< "----------------------------------------------------"
<< std::endl;
return 1;
};
}