The active flux scheme on Cartesian grids and its low Mach number limit

TL;DR

This paper introduces an active flux finite volume scheme for 2D acoustic equations on Cartesian grids, demonstrating stability, accuracy in discontinuous solutions, and low Mach number compliance without additional fixes.

Contribution

It extends the active flux scheme to two-dimensional Cartesian grids for acoustic equations, showing stability, accuracy, and low Mach number compliance.

Findings

Stable simulation of discontinuous solutions.

Low Mach number compliance without fixes.

Effective implementation on 2D Cartesian grids.

Abstract

Finite volume schemes for hyperbolic conservation laws require a numerical intercell flux. In one spatial dimension the numerical flux can be successfully obtained by solving (exactly or approximately) Riemann problems that are introduced at cell interfaces. This is more challenging in multiple spatial dimensions. The active flux scheme is a finite volume scheme that considers continuous reconstructions instead. The intercell flux is obtained using additional degrees of freedom distributed along the cell boundary. For their time evolution an exact evolution operator is employed, which naturally ensures the correct direction of information propagation and provides stability. This paper presents an implementation of active flux for the acoustic equations on two-dimensional Cartesian grids and demonstrates its ability to simulate discontinuous solutions with an explicit time stepping in a stable manner. Additionally, it is shown that the active flux scheme for linear acoustics is low Mach number compliant without the need for any fix.

Figures

Equations

$${\partial }_{t}q+\nabla \cdot \mathbf{f}(q)$$

$$q$$

$$\mathbf{f}=(f^x,f^y)f^x,f^y$$

$$\begin{array}{rl}{\partial }_t\rho +\nabla \cdot (\rho \mathbf{v})& =0\\ {\partial }_t(\rho \mathbf{v})+\nabla \cdot (\rho \mathbf{v}\otimes \mathbf{v}+p1)& =0\end{array}$$

$$\begin{array}{rl}{\partial }_{t}p+c\nabla \cdot \mathbf{v}& =0,\\ {\partial }_t\mathbf{v}+c\nabla p& =0.\end{array}$$

$$p(0,\mathbf{x})$$

$${\partial }_t(\nabla \times \mathbf{v})$$

$$\begin{array}{rl}{\partial }_{t}p+\frac{c}{ϵ}\nabla \cdot \mathbf{v}& =0,\\ {\partial }_t\mathbf{v}+\frac{c}{ϵ}\nabla p& =0.\end{array}$$

$${\partial }_t^2\mathbf{v}-c^2\nabla (\nabla \cdot \mathbf{v})$$

$$\begin{array}{rl}M\left[f\right](\mathbf{x},r)& =\frac{1}{4\pi }{\oint }_{S^2}\mathrm{d}\mathbf{y}f(\mathbf{x}+r\mathbf{y})=\frac{1}{4\pi }\int {\text{displaylimits}}_0^{2\pi }\mathrm{d}\phi \int {\text{displaylimits}}_0^{\pi }\mathrm{d}\vartheta \sin \vartheta f\left(\mathbf{x}+r\cdot \mathbf{n}\right)\end{array}$$

$$\mathbf{n}=\left(\begin{array}{c}\sin \vartheta \cos \phi \\ \sin \vartheta \sin \phi \\ \cos \vartheta \end{array}\right)$$

$$\displaystyle M^{\scalebox{0.5}{$2D$}}\left[f\right](x,y,r)=\frac{1}{2\pi r}\int\displaylimits_{0}^{2\pi}\mathrm{d}\varphi\int\displaylimits_{0}^{r}\mathrm{d}sf\left(x+s\cos\varphi,y+s\sin\varphi\right)\frac{s}{\sqrt{r^{2}-s^{2}}}$$

$$\displaystyle M^{\scalebox{0.5}{$2D$}}\left[x^{p}y^{q}\right](0,0,r)$$

$$\int {\text{displaylimits}}_0^r\mathrm{d}s\frac{s^m}{\sqrt{r^2-s^2}}$$

$$p(t,\mathbf{x})$$

$$\mathbf{v}(t,\mathbf{x})$$

$$p(t,\mathbf{x})=$$

$$\displaystyle\begin{split}\mathbf{v}(t,\mathbf{x})=&\,\mathbf{v}_{0}(\mathbf{x})-\frac{1}{ct}\partial_{r}\Big{(}r^{2}M\left[p_{0}\mathbf{n}\right](\mathbf{x},r)\Big{)}\Big{|}_{r=ct}\\ &+\int\limits_{0}^{ct}\frac{1}{r}\partial_{r}\left(\frac{1}{r}\partial_{r}\left(r^{3}M\left[\left(\mathbf{v}_{0}\cdot\mathbf{n}\right)\mathbf{n}\right](\mathbf{x},r)\right)-rM\left[\mathbf{v}_{0}\right](\mathbf{x},r)\right)\mathrm{d}r\end{split}$$

$${\partial }_t{\int }_{\mathcal{C}}\mathrm{d}\mathbf{x}q+{\int }_{\partial \mathcal{C}}\mathrm{d}\mathbf{x}\mathbf{n}\cdot \mathbf{f}(q)$$

$${\partial }_t{q}_{\mathcal{C}}+\frac{1}{|\mathcal{C}|}\sum _{e\subset \partial \mathcal{C}}|e|{\bar{f}}_e$$

$$q_{\text{recon}}(\mathbf{x})=q_{\mathcal{C}}\text{if}\mathbf{x}\in \mathcal{C}$$

$$\frac{1}{|\mathcal{C}|}{\int }_{\mathcal{C}}\mathrm{d}\mathbf{x}{q}_{\text{recon}}(\mathbf{x})$$

$$q_{\text{recon}}(\mathbf{p})$$

$${\mathcal{C}}_{ij}=\left[\left(i-\frac{1}{2}\right)\Delta x,\left(i+\frac{1}{2}\right)\Delta x\right]\times \left[\left(j-\frac{1}{2}\right)\Delta y,\left(j+\frac{1}{2}\right)\Delta y\right]\subset {\mathbb{R}}^2$$

$$q_{\text{recon},ij}:\left[-\frac{\Delta x}{2},\frac{\Delta x}{2}\right]\times \left[-\frac{\Delta y}{2},\frac{\Delta y}{2}\right]\to {\mathbb{R}}^n$$

$$q_{\text{recon}}(\mathbf{x})=q_{\text{recon},ij}(\mathbf{x}-{\mathbf{x}}_{ij})\mathbf{x}\in {\mathcal{C}}_{ij}$$

$$\mathbf{\xi }$$

$$a_{00}+a_{10}\xi +a_{20}{\xi }^2+a_{01}\eta +a_{11}\xi \eta +a_{21}{\xi }^2\eta +a_{02}{\eta }^2+a_{12}\xi {\eta }^2+a_{22}{\xi }^2{\eta }^2$$

$$q_{\text{recon},ij}(\mathbf{\xi })=\sum _{m=1}^9{c}_{m,ij}{b}_m(\mathbf{\xi })$$

$$b_m(\mathbf{\xi })$$

$$b_9({\mathbf{\xi }}_m)$$

$$\frac{1}{|\mathcal{C}|}{\int }_{\mathcal{C}}\mathrm{d}\mathbf{x}{q}_{\text{recon}}(\mathbf{\xi }(\mathbf{x}))$$

$$c_{m,ij}=q_{m,ij}\forall m\in \{1,\dots ,8\}$$

Full text

Keywords: Active flux; acoustic equations; hyperbolic conservation laws; low Mach number limit; vorticity preserving.

Mathematics Subject Classification: 35L65; 35L45; 65M08.

The active flux scheme on Cartesian grids and its low Mach number limit111We thank Yifan Bai of the University of Michigan for providing the results shown in Figure 15. WB acknowledges support of the German National Academic Foundation and the support of the German Academic Exchange Service (DAAD) with funds from the German Federal Ministry of Education and Research (BMBF) and the European Union (FP7-PEOPLE-2013-COFUND – grant agreement no. 605728).

Wasilij Barsukow222Institute of Mathematics, Zurich University, 8057 Zurich, Switzerland and Institute of Mathematics, Wuerzburg University, Emil-Fischer-Strasse 40, 97074 Wuerzburg, Germany, wasilij.barsukow@math.uzh.ch                   Jonathan Hohm333Institute of Mathematics, Wuerzburg University, Emil-Fischer-Strasse 40, 97074 Wuerzburg, Germany

Christian Klingenberg444Institute of Mathematics, Wuerzburg University, Emil-Fischer-Strasse 40, 97074 Wuerzburg, Germany                   Philip L. Roe555Department of Aerospace Engineering, University of Michigan, 1320 Beal Avenue, Ann Arbor, MI 48109, USA.

1 Introduction

Conservation laws, such as the Euler equations of ideal hydrodynamics, express conservation of density, momentum and energy of a compressible fluid. Numerical methods need to reflect the conservation property, as by the Lax-Wendroff theorem only conservative numerical schemes are able to converge to a (weak) solution of the equations. Cell-based methods consider a partition of space into computational cells. A natural class of numerical schemes for conservation laws are finite volume schemes. They interpret the discrete degree of freedom in a computational cell as the average of the dependent variable. Numerical fluxes assigned to the cell boundaries automatically ensure conservativity of the method.

In one spatial dimension the numerical flux can be successfully obtained from a so-called Riemann solver. The dependent variables are reconstructed in a piecewise constant manner allowing for jumps at the locations of cell interfaces. The Riemann problems that arise at cell interfaces are then solved (see e.g. [God59]). It is often possible to solve the Riemann problem exactly, but an approximate solution can lead to a similarly accurate scheme while requiring less computations. Therefore approximate Riemann solvers are very popular (e.g. [Bou04, Tor09]). Among others, relaxation (e.g. [JX95]) can be a successful strategy to construct them.

In multiple spatial dimensions one faces additional challenges. First of all, even approximate Riemann solvers are complicated (see e.g. [Zhe12, BK17] for exact solutions to multi-dimensional Riemann problems and e.g. [Col90, Bal12] for approximate solvers). This explains why dimensionally split methods are so popular: they consider the Riemann problems arising at the edges of the computational cell independently. The multi-dimensional problem is thus replaced by several one-dimensional ones.

Additionally, the solutions to Euler equations in multiple spatial dimensions exhibit a number of phenomena absent in one dimension. Particularly prominent features are vortices, that strongly dominate realistic multi-dimensional flows. In the limit of low Mach numbers, the solutions of compressible Euler equations tend to those of incompressible Euler equations. This is a multi-dimensional feature as well, because incompressible flows are trivial in one spatial dimension. Many finite volume methods suffer from artefacts when applied e.g. to the regime of low Mach number (see e.g. [GV99]). While Riemann solvers seem to yield very good results in one spatial dimension it is unclear whether this is true in multiple spatial dimensions, or whether numerical methods based on different concepts are more adequate. It has been, for instance, noticed that the failure in the regime of low Mach number is observed even when using exact multi-dimensional Riemann solvers ([GM04, BK17]).

It thus seems relevant to investigate alternative ways how multi-dimensional numerical methods can be constructed. Recently, such a method has been proposed by [ER13], as an extension of a method from [VL79] for linear advection. It has been given the name active flux. It involves a continuous reconstruction, together with pointwise degrees of freedom along the cell boundary. It thus does not require the solution of any Riemann problem. Continuous reconstructions are found not to stand in the way of computing discontinuous solutions. Indeed, even schemes that employ reconstructions with jumps across cell interfaces are mostly unable to resolve a shock wave sharply. They resolve a shock wave by a succession of jumps. Instead, a method with a continuous reconstruction resolves the discontinuity by a continuous function with a sharp gradient (see Figure 1). The details of the method are discussed in section 3 below.

This paper aims at contributing to a deeper understanding of the active flux method. First, a very general formulation of the method is given, that is independent of the grid or the equations (section 3). It thus highlights the most important features that distinguish the method from conventional finite volume schemes. Active flux in multiple spatial dimensions has so far been introduced for triangular computational cells. For certain applications triangles might be too complicated a mesh. This paper for the first time presents an implementation on two-dimensional Cartesian grids (section 4).

One of the objectives behind the development of the active flux method is to improve multi-dimensional simulations by including more multi-dimensional information in the scheme. Following [ER13], in this paper the active flux method is applied to the equations of linear acoustics. Linear acoustics is obtained when the Euler equations are linearized around a constant background with no velocity. In multiple spatial dimensions, linear acoustics cannot be reduced to some kind of multi-dimensional advection. This makes it a valuable system in order to investigate the performance of numerical methods in multiple spatial dimensions (as e.g. also in [DOR10]).

The active flux scheme pursues a strategy alternative to discontinuous reconstructions. The fact that the method does not use a Riemann solver raises the question whether maybe it is better suited for the low Mach number limit than Riemann solvers. Linear acoustics exhibits a low Mach number limit similar to that of the Euler equations and there exists an efficient methodology to study the behaviour of numerical schemes in this limit on Cartesian grids (introduced in [Bar19, Bar17]). Here, this methodology is extended to take into account the additional degrees of freedom along the cell boundary and is applied to the active flux scheme in section 6. It is shown that indeed, active flux is able to resolve the low Mach number limit of linear acoustics.

The presence of point values along the boundary makes it easy to evaluate the fluxes necessary for the update of the cell average. These point values, however, also require a time evolution. Here, following [ER13], the exact solution operator of linear acoustics is applied onto reconstructed initial data. However, this represents only one possible choice of an evolution operator for the point values and other evolutions are possible. The presentation of the general philosophy of active flux in section 3 therefore does not insist on a particular way of evolving the point values. This becomes particularly important for nonlinear equations without an available exact evolution operator. Recent work on approximate evolution operators for nonlinear equations is [HKS19].

2 Acoustic equations

This paper considers numerical methods for $n\times n$ systems of conservation laws in $d=2$ spatial dimensions:

[TABLE]

Only vectors with $d$ components are set in boldface symbols. Additionally, indices are never denoting a derivative.

The most prominent example of a system of conservation laws are the Euler equations

[TABLE]

with velocity $\mathbf{v}=(u,v)^{T}$, pressure $p$ and density $\rho$.

As a stepping stone towards a detailed understanding of numerical methods for the Euler equations in multiple spatial dimensions the acoustic equations are studied first. They are obtained as a linearization of the Euler equations around the state of constant density and pressure and vanishing velocity (see also e.g. [BK17]). Their symmetrized version reads

[TABLE]

with $c>0$ the speed of sound of the background state. This system is strongly hyperbolic in two spatial dimensions with eigenvalues $\{\pm c,0\}$. It has to be augmented by initial data

[TABLE]

and possibly boundary conditions.

The curl of the velocity is called vorticity. In the context of Euler equations, there is great interest in methods that are able to compute reliably flows that contain vorticity, so that the appearance of vortical structures is caused by physics and affected as little as possible by numerical error.

For the acoustic equations vorticity is stationary:

[TABLE]

Thus, for acoustics, an interesting class of methods are those that keep a discrete version of the vorticity exactly stationary (vorticity preserving), see e.g. [MR01, JT06, MT09, Bar19].

The low Mach number limit for the equations of linear acoustics is the limit $\epsilon\to 0$ of

[TABLE]

It is obtained by analogy from the low Mach number limit of the Euler equations. There, the pressure gradient $\nabla p$ is rescaled as $\frac{\nabla p}{\epsilon^{2}}$. The same is done here for acoustics, which after symmetrization yields (2.8). The interpretation of the Mach number $\epsilon$ as the ratio between advective and acoustic speeds in the setting of the Euler equations, however, is not applicable to the acoustic equations. Nevertheless this limit is useful in the design and understanding of low Mach number compliant schemes.

The low Mach number limit in (2.8) is the same as the limit $t\to\infty$ of the initial value problem (2.5)–(2.6). With appropriate boundary conditions, the long time solution is governed by the stationary states (for more details see [Bar19]).

Note that system (2.5) can be reduced to a wave equation for $p$, and a vector wave equation

[TABLE]

for $\mathbf{v}$. This latter, however, cannot generally be reduced to several scalar wave equations for the components of $\mathbf{v}$, because $\nabla(\nabla\cdot\mathbf{v})-\Delta\mathbf{v}=\nabla\times(\nabla\times\mathbf{v})$.

The numerical method will make use of the exact solution to the initial value problem (2.5)–(2.6) for non-differentiable initial data. The required (distributional) solution operator for initial data of such low regularity has been derived in detail in [BK17]. The solution formulae contain spherical means, that also appear in the solution of the scalar wave equation ([Joh78]). It is helpful to first consider three spatial dimensions:

Definition 1** (Spherical mean).**

The spherical mean $M\left[f\right](\mathbf{x},r)$ of an integrable function $f$ that depends on $\mathbf{x}=(x,y,z)\in\mathbb{R}^{3}$ is given by

[TABLE]

with the outward normal vector given by

[TABLE]

This paper concentrates on the two-dimensional case, though. If $f$ does not depend on $z$, then, with $s:=r\sin\vartheta$, it is possible to rewrite equation 2.10 as

[TABLE]

In both cases, if $f$ is polynomial, the spherical means can be evaluated analytically. Indeed, by shifting the point of integration one is left with

[TABLE]

Both integrals can easily be evaluated analytically because the angular integration bounds here are multiples of $\frac{\pi}{2}$. The second integral is

[TABLE]

In [ER13] for the initial value problem (2.5)–(2.6) the following solution formula appears:

[TABLE]

These formulas are derived from the classical Poisson formula for the scalar wave equation in second-order form [Whi11], which was used in [AGH00, Hag15] to derive high-order finite-difference methods for the scalar problem. For the acoustic equations in first-order form, and for the linearized Euler equations, (2.14) was applied to the pressure, and a similar equation to the velocity components, and employed in [ER13] to create an active flux scheme, with vorticity treated as a source term. The full extension to vortical flows (2.14)–(2.15) appears in [FR15], together with investigations of the active flux method on unstructured grids, including boundary conditions and nonlinearity. Finite-volume schemes based on the full extension are reported in [FG18].

For the numerical method it will be necessary to consider initial data that are continuous, but possess discontinuous first derivatives, such that terms like $\nabla\cdot\nabla p_{0}$ require clarification. The interpretation of this formula in the sense of distributions has been achieved in [BK17]. At the same time it has been shown that the formula can be rewritten as follows:

Theorem 1** (Solution operator).**

The solution to the initial value problem (2.5)–(2.6) is given by

[TABLE]

The derivatives are to be interpreted in the sense of distributions, if necessary.

For the proof of the theorem, and more details on its distributional version, see [BK17]. Below it is shown that for the reconstructions considered here all the derivatives in (2.16)–(2.17) exist in the strong sense and it is not necessary to consider anything distributionally. When trying to use formula (2.14)–(2.15), one would be forced to interpret them in the sense of distributions in order to compute the second derivatives correctly. Nothing seems to be gained when using (2.16)–(2.17), because they also contain second derivatives. However this time, the derivatives are all expressed with respect to $r$ and the initial data in the particular setup of this paper will turn out to have continuous derivatives with respect to $r$ (section 4.3). Thus when using formulae (2.16)–(2.17) it is possible to avoid dealing with distributional solutions, which is a substantial advantage.

It should be noted, that the exact solution may be expressed in various representations, which differ by the way they can be used. In this respect the above formulae are very different from those obtained, for example, in [Ost97, LMMW04] using bicharacteristics. There, analytical relations are derived, which connect the solution at time $t>0$ with the data at initial time via a so called mantle integral. This latter involves the solution at all intermediate times. Therefore, the knowledge of the data at initial time does not allow to immediately compute the solution at a later time exactly. With the above formulae (2.14)–(2.15) or (2.16)–(2.17), on other hand, the solution at time $t$ is given as a functional of the data at initial time only and thus can be determined immediately.

3 General structure of the active flux method

This section outlines the idea of the active flux method on general polygonal (unstructured) meshes and without specifying the hyperbolic system of PDEs that is to be solved. Several references to implementations of active flux for particular equations that appear in the literature are given in section 3.3. In section 4 finally, Cartesian meshes with rectangular cells are introduced, on which the active flux method is used to solve the acoustic equations (2.5).

The active flux method is an extension of the finite volume method. Therefore finite volume methods are first reviewed in section 3.1, before the active flux method is introduced in section 3.2.

3.1 Finite volume methods

For the finite volume method the discrete degree of freedom $q_{\mathcal{C}}$ is interpreted as the average of $q$ over the polygonal numerical cell $\mathcal{C}\subset\mathbb{R}^{d}$. For Cartesian grids, also the notation $q_{\mathcal{C}_{ij}}=:q_{ij}$ is standard. However, our presentation of the method for the moment does not depend on the nature of the grid. In order to construct the method one first integrates the conservation law (2.1) over the cell and applies Gauss’ law. This gives rise to fluxes through the cell boundary:

[TABLE]

Here $\mathbf{n}$ is the outward normal on the cell boundary $\partial\mathcal{C}$ of cell $\mathcal{C}$.

The finite volume method now replaces the exact expression (3.1) by a (yet undetermined) numerical approximation $\bar{f}_{e}$ of the flux through an edge $e\subset\mathcal{\partial}C$:

[TABLE]

Here $|e|$ is the length of the edge $e$ and $|\mathcal{C}|$ the area of the cell $\mathcal{C}$.

One way to obtain a numerical flux is by using a Riemann solver: One considers the initial value problem (2.1) with initial data $q_{\text{recon}}:\mathbb{R}^{d}\to\mathbb{R}^{n}$ that are piecewise constant

[TABLE]

This step is called reconstruction. The resulting Riemann problems across the cell boundaries are solved exactly or approximately in time. The numerical flux is then taken to be the time average of the flux of this solution through the cell boundary.

Alternatively, one might consider an initial value problem given by more complicated reconstructions inside the cells, which is referred to as generalized Riemann problem (see e.g. [GMP84, BLFR89, BAF03, TT02a, TT02b] and many others). This leads to schemes that are higher order discretizations of the original PDE. In this case the reconstruction $q_{\text{recon}}$ in cell $\mathcal{C}$ does not only depend on $q_{\mathcal{C}}$, but also on the values in neighbouring cells. Conservation still requires the average $\frac{1}{|\mathcal{C}|}\int_{\mathcal{C}}\mathrm{d}\mathbf{x}\,q_{\text{recon}}(\mathbf{x})$ of the reconstruction over any cell $\mathcal{C}$ to match the discrete degree of freedom $q_{\mathcal{C}}$ in this cell.

3.2 Active flux method

Now the active flux scheme is introduced. Its presentation in this section again is given without specifying the type of grid. In section 3.3 a brief overview of an implementation on one-dimensional grids and on triangles is given, before section 4 presents the novel implementation on Cartesian grids.

The active flux method differs from usual finite volume methods in how the numerical flux is obtained. Additional degrees of freedom666Recall that throughout the paper indices never denote a derivative. $q_{\mathbf{p}}$, $\mathbf{p}\in\partial\mathcal{C}$ are introduced, which are given the interpretation of point values and are distributed at a finite number of locations $\mathbf{p}$ along the cell boundary $\partial\mathcal{C}$. This immediately allows to use them in a quadrature formula in order to approximate the flux through the cell boundary. The update formula for the cell average $q_{\mathcal{C}}$ remains equation (3.2) just as for usual finite volume methods.

The additional degrees of freedom require an update procedure which shall be detailed next. Here, a reconstruction $q_{\text{recon}}:\mathbb{R}^{d}\to\mathbb{R}^{n}$ is used as initial data. The resulting initial value problem (2.1) is then solved (exactly or approximately) at the location $\mathbf{p}\in\partial\mathcal{C}$ of the point value at the cell boundary. In the implementation of active flux for the acoustic equations below, an exact solution of the initial value problem is used (equations (2.16)–(2.17)).

The reconstruction $q_{\text{recon}}$ has to fulfill both

[TABLE]

and

[TABLE]

The latter condition implies that the reconstruction is continuous across the cell boundary at least in all points $\mathbf{p}$. In fact, in this paper the reconstructions will be continuous across all of $\partial\mathcal{C}$. The increased number of conditions makes the appearance of high order interpolation polynomials natural, which make the scheme high order in space. High order temporal accuracy is obtained by using one or several substeps during the time integration (as discussed in section 4.4).

To summarize, the active flux method promotes the solution along the cell boundary (and thus in a sense the cell boundary flux) to the status of an independent degree of freedom (hence the name active flux). In usual finite volume methods on the other hand the cell boundary flux is a derived quantity. Another difference to usual finite volume methods is that the reconstruction is used to evolve the point values along the cell boundary (which are then used to compute the flux), rather than to obtain a flux directly by e.g. solving a Riemann problem. The algorithmic structure of the active flux method is shown in Figure 2 and in algorithm 1.

3.3 The standard active flux method

The active flux method as formulated in section 3.2 is very general. It is not restricted to a particular grid, nor to a particular form of the conservation law (2.1). At the same time it leaves a lot of freedom concerning its design, as neither the precise locations of the point values along the cell boundary nor the reconstruction have been specified in section 3.2.

Two particular designs appear in the literature and shall be briefly discussed here for the sake of an overview. The first goes back to [VL79] and considers a one-dimensional grid. The cell boundaries are points, and an additional degree of freedom is located at each one. The reconstruction is piecewise parabolic, which is the lowest degree polynomial that fulfills conditions (3.4) and (3.5) in this setting. The reconstruction is continuous across the cell boundaries due to (3.5), but has discontinuous derivative there. [VL79] discusses the application to linear advection; Burgers’ equation has been considered in [ER11a, Roe17], nonlinear hyperbolic systems in [ER11b].

An extension to two-dimensional triangular grids has been introduced in [ER13, Eym13]. Inspired by P2 interpolations the point values are located at vertices and edge midpoints of the triangles. (3.4) and (3.5) are 7 equations in this case. The reconstruction is taken to be a subset of the biparabolic polynomials (see [ER13], Table 2 for details). In [ER13] the method is used to numerically solve both linear advection and linear acoustics. Extensions to nonlinear problems are discussed in [Fan17, Mae17].

In both cases the reconstructions are piecewise parabolic, thus yielding third order accuracy in space. In order to reach the same accuracy for the temporal discretization, the point values along the cell boundary are evolved with half time steps from the same data. The quadrature formula for the numerical flux then is a space-time Simpson rule.

4 Active flux on a Cartesian mesh for the acoustic equations

In this chapter, an implementation of the active flux scheme on rectangular meshes is considered for the first time for the acoustic equations. Acoustic equations in multiple spatial dimensions are an important stepping stone in the development of the active flux method for the Euler equations. At the same time Cartesian grids are easy to implement and are used in a variety of applications.

A Cartesian mesh consists of rectangular cells $\mathcal{C}_{ij}$ with width $\Delta x$ and height $\Delta y$

[TABLE]

indexed by $(i,j)\in\mathbb{Z}^{2}$. Therefore $|\mathcal{C}|=\Delta x\Delta y$. The cell average $q_{\mathcal{C}_{ij}}$ is denoted by $q_{ij}$.

The algorithm consists of the four steps mentioned in section 3.2. The distribution of the point values along the cell boundary is presented in section 4.1. The reconstruction on rectangular cells is discussed in section 4.2. It is used as initial data in order to advance the point values at the cell boundaries forward in time. The corresponding solution operator for the equations of linear acoustics is presented in section 4.3. A quadrature rule is then applied to the point values along the cell boundaries in order to determine the numerical flux (section 4.4), and the fluxes are used to update the cell average.

4.1 Distribution of point values along the cell boundary

In this and the next section, the reconstruction is adapted to a rectangular cell. As on triangular grids mentioned in section 3.3, the point values are taken to be located at the corners of the cell (node values) and at edge centers (edge values). These locations are indexed as in Figure 4. Each cell $\mathcal{C}_{ij}$ has eight point values along the cell boundary: four node and four edge values. These degrees of freedom are correspondingly denoted by $q_{m,ij}$, $m=1,...,8$. They are shared by the adjacent cells, i.e. for example $q_{4,ij}=q_{8,i+1,j}$.

For usual finite volume methods, the cell average is the only degree of freedom per cell. Now additionally to the cell average there are 8 degrees of freedom distributed as point values along the cell boundary: There are four node values, each shared by four cells. Similarly, there are four edge values, which are shared by two cells each. Therefore per cell one is left with

• •

1 cell average

• •

1 node value

• •

2 edge values (horizontal and vertical)

In Figure 4 these degrees of freedom are shown in different colors. Indeed, each of them forms a lattice with spacings $\Delta x$ and $\Delta y$ in $x$ and $y$ direction, respectively. When solving linear acoustics with the active flux method, there are thus 12 variables per cell, i.e. 4 degrees of freedom per cell with 3 variables per degree of freedom.

4.2 Reconstruction

The reconstruction is subject to the constraints (3.4) and (3.5) from section 3: it has to be exact in the eight point values along the cell boundary and the average of the reconstruction over the cell has to yield the cell average. These are nine conditions. It is helpful to define the reconstruction $q_{\text{recon},ij}$ in any cell to refer to the cell midpoint as $\mathbf{x}=0$. Therefore denote by

[TABLE]

the reconstruction of any quantity $q$ in cell $\mathcal{C}_{ij}$ and by

[TABLE]

the reconstruction on the entire grid as introduced at the end of section 3.1.

In this coordinate frame, the locations of the point values at the cell boundary are denoted by $\mathbf{x}_{m}$, $m=1,\ldots,8$. E.g. $\mathbf{x}_{4}=(\Delta x/2,0)$.

To simplify notation it is also helpful to define a mapping to a reference cell. Consider coordinates $\mathbf{x}=(x,y)$ aligned with the cell edges and centered at a given cell. They can be transformed to reference coordinates $\boldsymbol{\xi}=(\xi,\eta)$ of a square cell with dimensions $\Delta x=\Delta y=2$. The edges of the square are located at $\xi=-1$, $\xi=1$, $\eta=-1$ and $\eta=1$. The transformation between the coordinates is

[TABLE]

$\boldsymbol{\xi}_{m}$ are the locations of the point values at the cell boundaries in the reference cell. They correspond to $\mathbf{x}_{m}=J\boldsymbol{\xi}_{m}$. The reconstruction $q_{\text{recon},ij}(\mathbf{x}(\boldsymbol{\xi}))$ will be denoted by the same symbol $q_{\text{recon},ij}(\boldsymbol{\xi})$ whenever there is no confusion possible.

A straightforward derivation of the reconstruction polynomial is to solve the linear system that arises by inserting a general biquadratic polynomial

[TABLE]

into the 9 equations that arise from (3.4) and (3.5). The polynomial consists of all possible combinations of $\xi$ and $\eta$ each up to second degree. It has $9$ free parameters. Hence, the solution of the linear system is unique.

This reconstruction can be obtained in a more elegant way, reminiscent of Lagrange polynomials. This representation of the interpolation will also be useful later. In every cell, one seeks a reconstruction in the form

[TABLE]

with basis functions $b_{m}(\boldsymbol{\xi})$ and coefficients $c_{m,ij}\in\mathbb{R}^{n}$. The expressions for $b_{m}$ are the same for every cell while the coefficients $c_{m,ij}$ vary, as they depend on the average value in the cell and the point values at the cell boundaries.

The reconstruction can be organized nicely by choosing the basis functions at most biquadratic with

[TABLE]

This implies that

[TABLE]

The value of $c_{9,ij}$ is to be determined after the basis functions have been obtained explicitly.

Theorem 2** (Interpolation basis).**

The following polynomials fulfill (4.7)–(4.8)

[TABLE]

Proof.

The result is verified by explicit calculation. ∎

Note: These products are polynomials that vanish along straight lines in $\xi$–$\eta$-plane, as shown in Figures 6 and 6 for $b_{1}$ and $b_{2}$.

Additionally, using (4.10) and computing (4.9) one obtains

[TABLE]

(4.6) becomes a biquadratic polynomial upon inserting the basis functions $b_{1}$, …$b_{9}$. With the coefficients (4.10) and (4.20) it must be equation (4.5) because the interpolation polynomial in this case is unique.

Theorem 3** (Continuous reconstruction).**

The biparabolic reconstruction (4.6) with basis functions from Theorem 2 and coefficients from (4.10) and (4.20) is continuous across cell interfaces.

Proof.

From (4.5) it is obvious that the reconstruction is parabolic along the coordinate axes. Consider an edge common to two cells. The two reconstructions on both sides of the edge both reduce to quadratic functions along the edge. Moreover, they agree in three points (at the two vertices and at the edge midpoint). This information specifies the quadratic function uniquely, and the two quadratic parabolas must agree.∎

Note that the derivatives perpendicular to the common cell edge in general are different. An example of a reconstruction for an arbitrary cell is sketched in Figure 8.

4.3 Updating the point values at cell boundaries

In order to advance in time the point values located at cell boundaries, one solves the initial value problem (2.1) with initial data

[TABLE]

at the location of the degree of freedom.

The exact solution of the acoustic equations on a rectangular mesh is given by equations (2.16)–(2.17). The reconstruction $q_{\text{recon}}$ is defined piecewise according to (4.3). The evolution of the point values at the cell boundaries therefore involves the reconstruction in several adjacent cells.

For example, the integration for the spherical mean of an edge value is split up into integrations over two hemispheres: The integration bounds for a vertical ($x=\mathrm{const}$) edge are $\varphi\in\{\frac{\pi}{2},\frac{3\pi}{2}\}$ and $\varphi\in\{0,\pi\}$ for a horizontal ($y=\mathrm{const}$) edge. It is convenient to define the partial spherical mean

[TABLE]

and $M_{0}^{2\pi}\equiv M$. Thus for a vertical edge, the spherical mean is computed as

[TABLE]

Analogously the integration for the spherical mean of a node value is divided into four different integrals with integration bounds $\varphi\in\{0,\frac{\pi}{2},\pi,\frac{3\pi}{2}\}$.

In a practical implementation the computational cost can be reduced by precomputing (analytically)

[TABLE]

for all relevant values of $\ell,m$, at all relevant locations $\mathbf{x}_{m}$, $m=1,\ldots,8$ and for all relevant angular domains. The spherical means are then obtained as linear combinations of these precomputed values. As equation (2.13) shows, they are polynomials in $r$, and this even allows to compute the derivatives and integrals with respect to $r$ in (2.16)–(2.17) analytically. All the operations like spherical means and the exact evolution operator thus can be precomputed in such a way that the final formula is merely a linear combination of the values in neighbouring cells. Thus, the evolution operator can be programmed in a way that makes its evaluation inexpensive as compared to the evaluation of an approximate Riemann solver for a standard finite volume scheme.

4.4 Time integration and numerical flux

In order to advance the cell average in time the fluxes through the cell boundaries are required. The strategy proposed here follows closely that of [ER13, Eym13] on triangular grids.

Consider notation of Figure 9, where an arbitrary edge is depicted and the edge value is denoted by M and the two node values by L and R.

The biparabolic reconstruction implies formally third order accuracy of the evolution of point values along the cell boundary. It is natural to strive for the same order of accuracy in other parts of the scheme, thus making it overall third order.

For the approximation of the average flux through a cell edge $e$

[TABLE]

two node values and one edge value are available. As the solution is thought of as continuous across cell boundaries, these values can be immediately used in the flux function $\mathbf{f}(q)$. Then Simpson’s rule as a quadrature formula leads to a scheme of third order.

In order to achieve the same accuracy in time, it is necessary to introduce an intermediate time step $n+\frac{1}{2}$. The point values at cell boundaries are computed at two time steps $n+\frac{1}{2}$ and $n+1$. Both go back to the initial data at time step $n$, see the green an blue regions in Figure 8. All fluxes are obtained by applying the flux function of the PDE to the point values along the boundary of the cell. Thus $f^{n}$ only depends on $q^{n}$, $f^{n+\frac{1}{2}}$ only on $q^{n+\frac{1}{2}}$, $f^{n+1}$ only on $q^{n+1}$, which, themselves, have been computed from the values $q^{n}$ during the update of the point values.

The space-time Simpson’s rule approximating (4.31) is then given by

[TABLE]

The numerical flux for a rectangular cell is visually presented in Figure 9. The cell average for $t^{n+1}$ is calculated with the finite volume method (3.2).

As the time integration amounts to an explicit scheme, there is a CFL-type time step restriction. All spherical means involved in the evolution of point values at cell boundaries have to remain inside the cell. The strongest constraint comes from the edge midpoint: the radius of the corresponding spherical mean has to be smaller than

[TABLE]

Thus the largest possible time step of the active flux method for the acoustic equations on a rectangular mesh is

[TABLE]

Experimentally, the active flux method was observed to be stable with time steps very close to this bound, suggesting that the bound is sharp.

As can be seen from Figure 11, the update of a cell average thus involves the reconstruction inside the cell and the reconstruction in all its 8 neighbours. The numerical fluxes and the cell average for an arbitrary cell are sketched in Figure 11.

5 Numerical results

In this section numerical results of the active flux scheme for the acoustic equations are presented. Given analytical expressions for the initial data, the point values at cell boundaries are initialized by evaluating these expressions at the corresponding locations. The initial cell average is computed with a quadrature formula following Simpson’s rule.

5.1 Stationary vortex

The first setup is the stationary vortex ($r=\sqrt{x^{2}+y^{2}}$)

[TABLE]

with $\mathbf{n}_{\varphi}=(-\sin(\varphi),\cos(\varphi))^{T}$. The simulation has been performed on a $50\times 50$ Cartesian mesh with $\Delta x=\Delta y=0.03$. Here, $p_{0}=0$ and the CFL number is $0.45$.

Contour plots of the cell averages of the absolute value of the velocity $\left|\mathbf{v}\right|=\sqrt{u^{2}+v^{2}}$ for this setup are shown in Figures13–13 for the times $t=0$ and $t=100$.

Additionally, the radial plot of the cell averages of the absolute value of the velocity $\left|\mathbf{u}\right|$ at $t=100$ and the exact solution are shown in Figure 14. The radial plots in Figure 14 suggest that the presented method is able to preserve stationary states to very high accuracy. This is not just due to the high order of the scheme. As is shown in the next section, the active flux method on rectangular grids is stationarity preserving. This is a property introduced in [Bar19] which means that the stationary states of the scheme discretize all the analytic stationary states.

Additionally, Figure 15 shows the same setup on a grid with a similar number of triangular elements. Average values for every element are shown.

5.2 Acoustic waves

A second test case is that of smooth acoustic waves propagating obliquely to the grid directions:

[TABLE]

Here $\delta=0.1$, $w=0.5\cos\alpha$ and $\alpha=\arctan\frac{1}{2}$ is used. The initial setup on a $50\times 50$ grid covering $[0,1]^{2}$ with periodic boundaries is shown in Figure 16 (upper left). This problem can be transformed into a one-dimensional setup along $\mathbf{b}$ which allows to derive an exact solution by the method of characteristics. Figure 16 (upper right) shows the error of the numerical solution at the midpoints of the vertical edges on $M\times M$ grids for different values of $M$. For this instationary setup a CFL number of 0.45 has been used. One observes a third order convergence of the method.

5.3 Spherical shock tube

The last test case is that of a radial shock tube. The initial data are chosen

[TABLE]

The simulation is performed on a $100\times 100$ grid. The results are shown in Figure 17. One observes that despite a continuous reconstruction the method is able to compute discontinuous solutions. The radial scatter plot demonstrates the symmetry of the solution. The over- and undershoots are due to the high order of the scheme and no limiting being employed.

6 Low Mach number compliance and vorticity preservation

Low Mach number compliance of a numerical scheme for linear acoustics has been found to be related to further concepts in [Bar19]. In particular, it has been shown that stationarity preserving schemes for the acoustic equations are just those that comply with the limit of low Mach number.

Preserving stationary states is a challenge to numerical methods. Whereas in one spatial dimension the stationary states of (2.1) are characterized by just an ODE, the multi-dimensional stationary states themselves fulfill a PDE. In particular, stationary states of (2.5) are characterized by a constant pressure and a divergencefree velocity. It is impossible to require a numerical method to keep stationary any sampling a divergencefree vector field. Thus, when treating stationary states numerically in multiple spatial dimensions it is necessary to content oneself with less.

[Bar19] made the following observation: even though the stationary states of (2.5) are all divergencefree vector fields, many schemes only keep trivial states stationary, such as the constant state $q=\mathrm{const}$. The stationary states of such numerical schemes thus are not a discretization of all stationary states of the PDE. In [Bar19] therefore, stationarity preserving schemes were introduced. They possess discrete stationary states that discretize all the analytic stationary states. In view of what has been said above, such schemes of course do not keep stationary any sampling of the analytic stationary states. However, for example in the case of linear acoustics, they keep stationary a particular discretization of the divergence. Initial data, sampled from a divergencefree vector field, evolve in time towards one of the discrete stationary states, which then persists for an arbitrarily long time without being diffused away. For a thorough presentation of these concepts the reader is referred to [Bar19, Bar17, Bar18].

It turns out that for the acoustic equations stationarity preservation is equivalent to the scheme being vorticity preserving, i.e. possessing a discrete vorticity that remains stationary. Therefore a numerical scheme that is not vorticity preserving also fails to discretize all the stationary states. It introduces so much diffusion that all but the most trivial stationary states are decaying in time. Numerical results of section 5 suggest that this is not the case for the active flux scheme on Cartesian grids.

In this chapter, it is shown that the active flux method on a rectangular mesh is stationarity preserving when used for numerically solving the two-dimensional acoustic equations. The framework presented in [Bar19] has to be generalized in order to take into account the point values at cell boundaries that are present in the active flux scheme. The proof involves applying the discrete Fourier transform to the numerical scheme, which is discussed in section 6.1. Stationarity preservation of the active flux scheme is proven in Corollary 1 in section 6.2. The discrete vorticity that is kept stationary by the scheme is discussed there as well.

6.1 The active flux scheme in Fourier space

Stationarity preservation can be most efficiently studied in Fourier space. A prerequisite is that the scheme under consideration is linear. For the active flux scheme this is established in the following

Theorem 4**.**

The active flux scheme for the acoustic equations on Cartesian grids is linear in the discrete degrees of freedom.

Proof.

The reconstruction 4.6 is linear in the cell average and the point values at cell boundaries. The spherical mean (equation 2.10) is a linear functional. Thus the update of the point values at cell boundaries from time step $n$ to the time step $n+\frac{1}{2}$ or $n+1$ is a linear function of the quantities at time step $n$. The numerical flux obtained by the Simpson’s rule is a linear functional of the involved point values at cell boundaries. Finally, the finite volume method used to update the cell average is linear in the fluxes. ∎

As the acoustic equations are a linear problem, and thus the active flux scheme is linear as well, it is possible to consider it in Fourier space. This is analogous to the procedure in [Bar19], and switching to Fourier space allows to prove stationarity preservation easily. In the following, the nomenclature from Figure 4 in section 4.1 is used to differentiate between the different kinds of degrees of freedom. The additional point values at cell boundaries that are present in the active flux scheme can be easily incorporated in the framework. The discrete Fourier transform with respect to the spatial variable can be applied to each of the lattices $q_{ij}$, $q^{\scalebox{0.5}{N}}_{ij}$, $q^{\scalebox{0.5}{EV}}_{ij}$, $q^{\scalebox{0.5}{EH}}_{ij}$ defined in section 4.1: Any function $q_{ij}(t)$ is expressed as

[TABLE]

with the wave vector $\mathbf{k}=(k_{x},k_{y})$ and the imaginary unit $\mathbb{i}$. The hat in $\hat{q}$ denotes the discrete spatial Fourier transform of $q$. The overall prefactor $\exp(\mathbb{i}\Delta xk_{x}i+\mathbb{i}\Delta yk_{y}j)$ will appear in all terms.

For linear acoustics solved with the active flux method, the Fourier transforms shall be ordered as follows (see also Figure 4)

[TABLE]

The translation operators that convey shifts by one cell in $x$-direction and $y$-direction are defined as

[TABLE]

such that Equation 6.1 becomes

[TABLE]

The Fourier transform of any two-dimensional linear finite difference formula corresponds to a Laurent polynomial in $t_{x}$ and $t_{y}$. For example:

[TABLE]

A more complicated finite difference formula can be, for example, expressed as

[TABLE]

Also the reconstruction $q_{\text{recon},ij}(\mathbf{x})$ in cell $\mathcal{C}_{ij}$ is a linear function of the degrees of freedom in $\mathcal{C}_{ij}$ which allows to apply the discrete Fourier transform to it:

[TABLE]

Note that the spatial variable $\mathbf{x}$ is unaffected, because the origin $\mathbf{x}=0$ has been chosen to be the cell midpoint of the corresponding cell, and the discrete Fourier transform only acts on $i,j$. Then upon the discrete Fourier transform one can rewrite the reconstruction in any other cell $\mathcal{C}_{i+s,j+p}$ as

[TABLE]

On the other hand, the evolution operator (2.16)–(2.17) only acts on $\mathbf{x}$. It is thus possible to construct a $12\times 12$ matrix $\mathcal{A}$ which describes the evolution of the Fourier modes $\hat{Q}$ introduced in (6.2):

[TABLE]

Theorem 5**.**

The active flux scheme for linear acoustics in two spatial dimensions is stationarity preserving iff \det\Big{(}\mathcal{A}(\Delta t;t_{x},t_{y})-\mathbb{1}_{12\times 12}\Big{)}=0 $\forall\Delta t,t_{x},t_{y}$.

Proof.

This is Theorem 2.11 from [Bar19]. ∎

The matrix $\mathcal{A}$ can be constructed explicitly. However, this involves computations of considerable length. Therefore in the following section a simpler strategy is applied. Its main ingredient is the observation that the active flux scheme uses the exact evolution operator for the point values at cell boundaries.

6.2 Stationarity preserving reconstruction

Clearly, at continuous level, if the initial data for (2.5) fulfill

[TABLE]

then they remain stationary. (This statement needs to be understood in the sense of distributions, if necessary – see [BK17] for more details.) Therefore one path towards understanding the stationary states of the active flux scheme is to study under which conditions the reconstruction fulfills

[TABLE]

Assume this to be the case for some choice of the discrete values on the grid. Then at least the point values at cell boundaries will remain precisely stationary over one time step. This shall be studied now in more detail; the question of whether the cell averages remain stationary is postponed and taken up at the end.

Theorem 5 states that a scheme is stationarity preserving (and thus vorticity preserving) if there is a zero eigenvalue. The corresponding eigenvector can be seen as a non-zero choice of

[TABLE]

that implies that all the degrees of freedom remain stationary. Therefore

Theorem 6**.**

The active flux scheme for linear acoustics is stationarity preserving, if there is a non-zero choice of $\hat{Q}$ that implies (6.11).

Note: As such a non-zero $\hat{Q}$ spans an eigenspace, any multiple of it implies (6.11) as well.

Observe the tremendous simplification from Theorem 5 to Theorem 6. Whereas showing \det\Big{(}\mathcal{A}(t_{x},t_{y})-\mathbb{1}_{12\times 12}\Big{)}=0 requires writing down explicitly all the spherical means (like in (4.30)) and evolution operators, the statement of Theorem 6 refers to a property of the initial data only. This is only possible because the evolution operators are exact and one thus can express in simple words which data they keep stationary.

Theorem 7**.**

If $\hat{Q}$ is such that

[TABLE]

[TABLE]

then $p_{\text{recon}}=0$ and $\mathrm{div\,}\mathbf{v}_{\text{recon}}=0$ uniformly.

Proof.

Consider the biparabolic reconstruction of equation (4.6) and Theorem 2 together with equation (4.20). Applying the discrete Fourier transform one can rewrite

[TABLE]

$p_{\text{recon}}(\mathbf{x})=0$ is thus clear. Moreover, one easily can compute $\mathrm{div\,}\mathbf{v}_{\text{recon}}(\mathbf{x})$ by differentiating the basis functions $b_{m}$ as given in theorem 2 and verify $\mathrm{div\,}\mathbf{v}_{\text{recon}}=0$. The computation is lengthy but uneventful, and is thus omitted. ∎

Note: From (6.17), (6.14) it follows that

[TABLE]

and from (6.15)–(6.16)

[TABLE]

These equations are discretizations of $\mathrm{div\,}\mathbf{v}=0$.

Note 2: Having thus found one eigenvector that belongs to the eigenvalue 0 of $\mathcal{A}(t_{x},t_{y})-\mathbb{1}_{12\times 12}$ it is not a priori clear whether there are more. This can be checked by actually computing the kernel of $\mathcal{A}(t_{x},t_{y})-\mathbb{1}_{12\times 12}$. Due to the extreme length of the expressions, the one-dimensionality of the kernel could so far only be verified using mathematica.

Corollary 1**.**

The active flux scheme for linear acoustics in two spatial dimensions is stationarity preserving.

Proof.

The proof consists of two parts. First, stationarity of the average values has to be checked. Second, the multi-step integration procedure (section 4.4) has to be taken into account.

1. i)

Assume the initial data to fulfill (6.14)–(6.17). By the above theorem this implies stationarity of the point values at cell boundaries and thus the three brackets in (4.32) are all equal. The values of $p$ are all zero, and thus the fluxes of $u$ and $v$ are zero. The change of the cell average of the pressure is

[TABLE]

Equations (6.20), (6.22) imply that (6.24) vanishes identically. Thus the cell averages are stationary. 2. ii)

Assume again that the initial data fulfill (6.14)–(6.17). Then at time step $n+\frac{1}{2}$ and at time step $n+1$ they are equal to those at time step $n$. This is because they both are computed from the same initial data at time step $n$. The update of the cell average happens only at time step $n+1$. Therefore stationarity preservation is valid independently of how many steps are used for the integration in time.

This completes the proof.∎

Stationarity preservation has several consequences: The eigenvector $\hat{Q}$ of theorem 7 is the Fourier transform of those data that the active flux scheme keeps exactly stationary. In [Bar19] it is shown that many numerical schemes add so much diffusion that only trivial (e.g. constant) stationary states are stationary at the numerical level as well. The active flux scheme on the other hand is stationarity preserving, i.e. it keeps stationary a discretization of all the stationary states of the PDE. By inverting the Fourier transform one obtains from (6.20)–(6.23) the following discrete relations that characterize the discrete stationary states:

[TABLE]

Here the following notation has been used:

[TABLE]

The notation is combined for different directions, e.g.

[TABLE]

By [Bar19] a stationarity preserving scheme is vorticity preserving. This means that there exists a discretization of $\nabla\times\mathbf{v}$ that is kept stationary by the scheme (even if the solution itself is not stationary). This is a discrete counterpart to Equation (2.7). The discrete Fourier transform of the discrete vorticity is given by the left eigenvector $\Omega$ belonging to eigenvalue zero. As $\mathcal{A}$ depends on $\Delta t$, one is facing the slightly surprising situation that the discrete vorticity might depend on $\Delta t$. Indeed, this might have occurred for the discrete stationary states already, but turns out not to be the case. Additionally, the amount of computations made it impossible to determine the explicit shape of the discrete vorticity. However, its existence is clear by the results of [Bar19]. From the dimension of the kernel it is also clear that there is only one discrete vorticity that is invariant. Only this particular discrete vorticity will be kept stationary while other discretizations will undergo some evolution.

That the discrete vorticity depends on $\Delta t$ can be shown as follows: Consider an expansion of $\mathcal{A}$ in powers of $\Delta t$:

[TABLE]

Independently of $\Delta t$, the vector $\hat{Q}$ from theorem 7 is a right eigenvector of $\mathcal{A}-\mathbb{1}$ corresponding to an eigenvalue 0. Therefore one concludes that $\mathcal{A}^{(1)}\hat{Q}=\mathcal{A}^{(2)}\hat{Q}=\mathcal{A}^{(3)}\hat{Q}=0$. It is possible to compute the left eigenvector $\Omega^{(1)}$ of $\mathcal{A}^{(1)}$:

[TABLE]

This would correspond to a discrete vorticity

[TABLE]

One can verify, however, that $\Omega^{(1)}\mathcal{A}^{(2)}\neq 0$ and thus $\Omega\neq\Omega^{(1)}$.

Finally, low Mach compliance is another consequence of stationarity preservation as discussed in [Bar19]. Von Neumann stability of the method implies that no Fourier mode is growing. Thus the long time solution is governed by the discrete stationary states, and stationarity preserving schemes discretize all of the analytic stationary states (for more details see [Bar19]). As $\epsilon$ goes to zero, one thus observes that the numerical solution obtained with the active flux scheme is a discretization of the limit solution.

In order for finite volume methods to be stable under explicit time integration upwinding is introduced. Many numerical schemes then strongly diffuse stationary states, and equivalently are unable to resolve the low Mach number limit. One remedy is to use implicit time integration without upwinding (central discretizations) (e.g. [CDK12] and many others). If fully explicit time integration is required (e.g. if both low Mach number regions and shocks are present in the flow), then many numerical schemes are made low Mach compliant by introducing some kind of fix (e.g. [DJOR16, BEK*+*17]). For the active flux scheme presented here such a fix is not needed: the exact evolution operators at the boundary use the correct directions of information propagation while not spoiling the low Mach number limit.

7 Summary

The active flux scheme is a finite volume method with additional degrees of freedom located on the cell boundary. Introduced in [VL79] for one-dimensional linear advection, it has since been extended to triangular grids ([ER13]) and other systems of hyperbolic PDEs. This paper presents an implementation for two-dimensional Cartesian grids. This underlines the viewpoint that the active flux is a concept that can be used with considerable flexibility with respect to the computational grid.

Cartesian grids have several advantages. Apart from the ease of implementation and low computational cost, the limit of low Mach numbers can be studied very efficiently for linear acoustics on Cartesian grids ([Bar19]). Low Mach compliance of a numerical scheme for linear acoustics is linked to it being stationarity preserving, i.e. to be able to keep stationary a discretization of all the analytic stationary states. Here, the framework of [Bar19] is used in order to show that the active flux possesses this property, and thus is low Mach compliant. The stationary states of the active flux scheme are obtained explicitly. They are closely linked to the solutions of the active flux scheme in the limit of low Mach numbers.

Finally, as stationarity preservation is equivalent to vorticity preservation, it is thus possible to show that the active flux scheme for linear acoustics is vorticity preserving. This can often be checked easily if it is known which discretization of the vorticity remains stationary. The approach via stationarity preservation, on the other hand, allows to check vorticity preservation without having to know the discrete vorticity in advance. We expect analogous results to hold true for the active flux scheme on three-dimensional Cartesian grids as well.

Several numerical examples show the good performance of the active flux scheme in practice. In particular they show that active flux seems to be stable with the maximum CFL condition. To study this theoretically is subject of future work. As the scheme is of high order limiting is a further aspect of future investigation.

Low Mach compliance is shown for the scheme endowed with an exact evolution operator for the point values at cell boundaries. This is different from Riemann solver based schemes, where even an exact Riemann solver does not in general lead to a low Mach compliant scheme. This, therefore, is clearly an advantage of the active flux scheme. However, the exact solution operator is generally unavailable for nonlinear systems of hyperbolic PDEs. Therefore approximate evolutions will be necessary. With the results of this paper in mind, future work shall focus on choosing them such that low Mach compliance is retained.