1 Introduction

Abstract

A new Active Flux method for the multi-dimensional Euler equations is based on an additive operator splitting into acoustics and advection. The acoustic operator is solved in a locally linearized manner by using the exact evolution operator. The nonlinear advection operator is solved at third order accuracy using a new approximate evolution operator. To simplify the splitting, the new method uses primitive variables for the point values and for the reconstruction. In order to handle discontinuous solutions, a blended bound preserving limiting is used, that combines a priori and a posteriori approaches. The resulting method is able to resolve multi-dimensional Riemann problems as well as low Mach number flow, and has a large domain of stability.

Keywords: Active Flux, Euler equations

Mathematics Subject Classification (2010): 65M08, 65M70, 76M12, 35L45

Wasilij Barsukow¹¹1Bordeaux Institute of Mathematics, Bordeaux University and CNRS/UMR5251, Talence, 33405 France, wasilij.barsukow@math.u-bordeaux.fr

1 Introduction

The Active Flux method, originally introduced for 1-d linear advection in [vL77], has received much attention since the pioneering works [ER11b, ER13]. As degrees of freedom, it employs averages and point values at cell interfaces. These are independent, i.e. the Active Flux method consists of update equations for both the averages and the point values. The former is easy to obtain for conservation laws: since the point values are located at cell interfaces, the flux can be directly evaluated; no Riemann solvers are needed and one recognizes the continuous nature of the spatial approximation. While the update of the averages is “exact”, i.e. as accurate as the point values, it cannot incorporate upwinding, since the latter always implies some kind of additional, artificial diffusion.

It is the point values which need to include the upwinding necessary for stability. The initial method from [vL77] uses characteristic tracing and a reconstruction, whose value at the foot of the characteristic is taken as the point value at the next time step. This procedure naturally includes upwinding, of course, and yields a one-stage method stable up to $\text{CFL}=1$ . Generalizations of this update procedure appeared under the name of (approximate) evolution operators in e.g. [ER13, Fan17, BHKR19, Bar21a, CHLM24].

It was shown in [Bar21a] that one needs to go one order of accuracy beyond local linearization to obtain a third-order Active Flux method. For Burgers’ equation $\partial_{t}u+u\partial_{x}u=0$ , instead of the local linearization $x\mapsto x-u(x)t$ , for example, [Roe17] suggests to use

\displaystyle x\mapsto x-\frac{u(x)t}{1+t\partial_{x}u(x)}\simeq x-u(x)t+t^{2}% u(x)\partial_{x}u(x)+\mathcal{O}(t^{3})

(1)

while [Bar21a] proposes to iterate

\displaystyle x\mapsto x-u\Big{(}x-u(x)t\Big{)}t\simeq x-u(x)t+t^{2}u(x)% \partial_{x}u(x)+\mathcal{O}(t^{3})

(2)

In [Bar21a], a sufficiently accurate approximate evolution operator was achieved for any hyperbolic system of conservation laws in one spatial dimension. It has been shown that it is by no means sufficient to merely iterate the linearization a few times to achieve third order of accuracy – this only works for scalar conservation laws since their characteristics are straight even in the nonlinear case (see Section 3.4 for further explanations). An alternative, ADER-inspired approach is [HKS19].

Concerning multiple spatial dimensions, early effort ([Fan17, BHKR19]) focused on the equations of linear acoustics and the corresponding evolution operator. The exact solution for a general class of data (in particular those not differentiable) was obtained in [BK22], initially destined to study genuinely multi-dimensional Godunov methods. The exact solution was used to demonstrate that even the complete solution of the multi-dimensional Riemann problem does not render a Godunov method stationarity or vorticity preserving. In [BHKR19] the exact acoustic evolution operator was used to construct an Active Flux method on two-dimensional Cartesian grids. Contrary to the Godunov method, the Active Flux method did turn out to be stationarity and vorticity preserving.

To construct an evolution operator for nonlinear systems of conservation laws in multi-d is significantly more difficult, since third order accuracy requires going beyond local linearization. At this point, in [Abg22, AB23], a semi-discrete version of Active Flux appeared. There, the point value is evolved according to an ODE in time, with the space derivative replaced by a finite difference formula (a non-standard one, since it involves both point values and averages in some neighbourhood of the degree of freedom). The time discretization is provided by a standard Runge-Kutta method according to the method of lines. The advantage of this approach is its immediate applicability to nonlinear systems in multi-d (e.g. [AL24, ABK25, DBK25]), with the structure preserving properties still in place ([BKKL24]). A slight disadvantage are the reduced low CFL numbers (0.2 and less), which [Roe21] connects to the usage of a multistage Runge-Kutta method. In fact, particular choices of less compact finite difference formulas do yield stability all the way up to $\text{CFL}=1$ for RK3 ([AB23]), but they are difficult to generalize in multi-d. However, undeniably, a Runge-Kutta method adds its own numerical diffusion.

The continuous approximation lets the semi-discrete Active Flux method appear different from, say, the Discontinuous Galerkin method, but the absence of the evolution operator also moves it away from van Leer’s and Roe’s original ideas. In the context of implicit methods for linear advection it was shown in [BB24] that one-stage implicit Active Flux methods obtained using an evolution operator were significantly better than semi-discrete Active Flux methods integrated with standard (implicit) integrators, at the same (or less) cost. This allows to conjecture that it might be worth striving for an Active Flux method for the Euler equations in multi-d that is based on evolution operators.

To the author’s knowledge, two such evolution-operator based Active Flux methods for multi-d Euler equations are available in the literature. [RMF18] suggests to start out from a Lax-Wendroff/Cauchy-Kowalevskaya expansion. The terms $\mathcal{O}(t)$ involve locally linearized acoustics and advection separately, while terms $\mathcal{O}(t^{2})$ involve deviations from a pure local linearization and interactions between these operators. The treatment of these extra terms remains difficult. The different operators of the Euler equations were treated in [Fan17, Mae17], yet a clear prescription how the problematic terms are to be discretized seems still to be missing, as is an investigation what impact particular choices would have in practice.

The recent work [CHLM24] employs second-order approximations to the evolution operators of linear advection and acoustics, but does not split them. This requires treating all the ways how a tilted characteristic cone (subject to “wind”) can cut through the cells in the vicinity of a point value. Here, the high-order terms that would be required for a third-order method are not included.

In view of such difficulties there seems to exist a certain interest in accepting second-order accuracy in time in certain ingredients of the numerical method for the sake of computational efficiency. The practical value of a time integration whose order of accuracy is not quite matching that of space discretization might depend on the application in question and on the willingness to spend effort on implementation and to accept a possible increase in computational effort. It is, for example common practice to use e.g. RK3 in conjunction with arbitrarily high-order methods (e.g. [JS96, Vil19, AB23]). Formally second-order methods that employ space reconstruction of degree up to 6 can bring significant improvement even without the corresponding increase in order of accuracy in time ([LAB⁺24]).

This paper proposes to use

–

a simple operator split between acoustics and advection, unsurpassed in computational efficiency, yet second-order; with
–

a fully third-order approximate evolution operator for the (nonlinear) advection in multi-d, and
–

a locally linearized, yet otherwise exact evolution operator for acoustics.

The advantage of the exact evolution operator is the inclusion of all the multi-dimensional information. Moreover, since the linearized method reduces to [BHKR19], its von Neumann stability limit is the same (contrary to [CHLM24]). The present paper gives a large amount of detail concerning the efficient implementation of the exact evolution operator. A bound-preserving limiting strategy combining one similar to that in [DBK25] with an a posteriori limiting for the point values is presented as well.

The result is a method for the Euler equations, possibly having some similarity with [RMF18], yet also clear differences, e.g. in the choice of a non-conservative advective operator and the usage of primitive variables for the reconstruction and the point values (cell averages are, naturally, in conserved variables). The simultaneous usage of primitive and conservative formulations has been advertised in [AL24] for the semi-discrete Active Flux, and already in [Bar21a] primitive variables at point values were in use for the fully discrete Active Flux method. The difference of the present work to [CHLM24] is the usage of the exact evolution operator, with more accuracy at essentially the same computational cost due to efficient implementation, and the algorithmic simplification that is brought about by the operator splitting. The purpose of this paper therefore is to present this simple one-stage method for the multi-d Euler equations. Beside that, the implementation details concerning the exact evolution operator for linear acoustics, at that time not included in [BHKR19], hopefully can now serve as a reference for future use.

The paper’s organization is as follows: Section 2 introduces the efficient implementation of the evolution operator for linear acoustics, complemented by Appendix A detailing the computation of spherical means. The method for the Euler equations and its properties are described in Section 3, followed by the numerical results of Section 4. The Appendices B and C contain pseudocode implementations of the acoustic evolution operator.

In the following, computations often are performed in three space dimensions, and two-dimensional results are obtained by assuming that the functions involved do not depend on $z$ . Vectors are set in boldface with, in particular, $\mathbf{x}=(x,y,z)^{\text{T}}$ . The notation $\mathbf{a}\cdot\mathbf{b}$ for $\mathbf{a},\mathbf{b}\in\mathbb{R}^{n}$ denotes the Euclidean scalar product of $\mathbf{a}$ and $\mathbf{b}$ ; $S^{2}:=\{\mathbf{x}\in\mathbb{R}^{3}:|\mathbf{x}|=1\}$ is the unit sphere in three space dimensions.

2 Evolution operator for linear acoustics

The exact solution of linear acoustics

	$\displaystyle\partial_{t}\mathbf{v}+c\nabla p$	$\displaystyle=0$	$\displaystyle\mathbf{v}$	$\displaystyle\colon\mathbb{R}^{+}_{0}\times\mathbb{R}^{3}\to\mathbb{R}^{3}$		(3)
	$\displaystyle\partial_{t}p+c\nabla\cdot\mathbf{v}$	$\displaystyle=0$	$\displaystyle p$	$\displaystyle\colon\mathbb{R}^{+}_{0}\times\mathbb{R}^{3}\to\mathbb{R}$		(4)

on $\mathbb{R}^{3}$ endowed with initial data $(\mathbf{v}_{0},p_{0})$ is

$\displaystyle p(t,\mathbf{x})$	$\displaystyle=\partial_{r}\left(rM[p_{0}](\mathbf{x},r)\right)-\frac{1}{r}% \partial_{r}\left(r^{2}M[\mathbf{n}\cdot\mathbf{v}_{0}](\mathbf{x},r)\right)$	(5)
$\displaystyle\mathbf{v}(t,\mathbf{x})$	$\displaystyle=\frac{2}{3}\mathbf{v}_{0}(\mathbf{x})-\frac{1}{r}\partial_{r}% \left(r^{2}M[p_{0}\mathbf{n}](\mathbf{x},r)\right)+\partial_{r}\left(rM[(% \mathbf{v}_{0}\cdot\mathbf{n})\mathbf{n}](\mathbf{x},r)\right)$
	$\displaystyle-M\left[\mathbf{v}_{0}-3(\mathbf{v}_{0}\cdot\mathbf{n})\mathbf{n}% \right](\mathbf{x},r)-\int_{0}^{ct}\mathrm{d}r\frac{1}{r}M\left[\mathbf{v}_{0}% -3(\mathbf{v}_{0}\cdot\mathbf{n})\mathbf{n}\right](\mathbf{x},r)$	(6)

where all formulae are understood to be evaluated at $r=ct$ after performing all the differentiations (see [BK22] for the derivation) and $M[f]\colon\mathbb{R}^{3}\times\mathbb{R}^{+}_{0}\to\mathbb{R}$ is the spherical mean of $f\colon\mathbb{R}^{3}\to\mathbb{R}$ , defined as

\displaystyle M[f](\mathbf{x},r)

\displaystyle:=\frac{1}{4\pi}\int_{S^{2}}\mathrm{d}\mathbf{y}\,f(\mathbf{x}+r% \mathbf{y})

(7)

For the definition of spherical means such as $M[\mathbf{n}\cdot\mathbf{v}_{0}]$ , see Section A first. The data $(\mathbf{v}_{0},p_{0})$ needs generally to be twice continuously differentiable for the solution formula from [Fan17] to be applicable, while (5)–(6) has the advantage that it only requires regularity with respect to the distance variable $r$ . Thus, it remains valid even if the data are discontinuous, as long as the discontinuities are straight lines through $\mathbf{x}$ , and this is precisely what is needed in the numerical setting. The formulation (5)–(6) thus avoids having to deal with distributions.

This formula is valid in three spatial dimensions, but one can of course restrict oneself to initial data that do not depend on certain variables, or initial velocity fields with certain components zero. In the following, $u$ and $v$ will be used to denote the $x$ and $y$ components of $\mathbf{v}$ .

The degrees of freedom of Active Flux are located at edges and nodes, with data being piecewise polynomial. It thus makes sense to divide the solution up in contributions of half-planes/quadrants. Consider, more generally, initial data in some wedge $W=[\varphi_{\text{min}},\varphi_{\text{max}})\times[0,2\pi]$ centered at $\mathbf{x}$ and denote by $p^{W}(t,\mathbf{x}),\mathbf{v}^{W}(t,\mathbf{x})$ the solution obtained by taking in (5)–(6) everywhere spherical averages $M^{W}$ over the wedge instead of $M$ . Then, if $\{W^{i}\}_{i}$ is a finite partition of $[0,2\pi)\times[0,2\pi]$ , then

	$\displaystyle M[f](x,r)$	$\displaystyle=\frac{1}{4\pi}\int_{\bigcup_{i}W_{i}}\mathrm{d}\mathbf{y}f(% \mathbf{x}+r\mathbf{y})=\frac{1}{4\pi}\sum_{i}2(\Delta\varphi)_{i}\frac{1}{2(% \Delta\varphi)_{i}}\int_{W_{i}}\mathrm{d}\mathbf{y}f(\mathbf{x}+r\mathbf{y})$		(8)
		$\displaystyle=\frac{1}{2\pi}\sum_{i}(\Delta\varphi)_{i}M^{W_{i}}[f](x,r)$		(9)

Thus,

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

\displaystyle=\frac{1}{2\pi}\sum_{i}|W_{i}|p^{W_{i}}(t,\mathbf{x})

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

\displaystyle=\frac{1}{2\pi}\sum_{i}|W_{i}|\mathbf{v}^{W_{i}}(t,\mathbf{x})

(10)

The ability to combine solutions across wedges as a simple sum is a consequence of the solution formula above. Finally, when applying the solution formula for every degree of freedom, $\mathbf{x}$ is fixed and the result is simply a polynomial in $t$ with easily precomputable coefficients. The next sections state these polynomials; by linearity one can consider initial data in different variables separately. For the sake of giving the reader a quicker access to the Section 3, which goes on to develop a method for the Euler equations, detailed computations are found in the Appendix A.

2.1 Initial data in $p$

First, consider initial data in $p$ only. The reconstructions in use are always polynomial; hence assume that $p_{0}(x,y)=x^{a}y^{b}$ in some wedge $W$ . Then

$\displaystyle p^{W}(t,\mathbf{x})$	$\displaystyle=\partial_{r}\left(rM[p_{0}](\mathbf{x},r)\right)\Big{\|}_{r=ct}$	(11)
	$\displaystyle\overset{\eqref{eq:sphericalmeanscalarfunction}}{=}\sum_{k=0}^{a}% \sum_{\ell=0}^{b}{a\choose k}{b\choose\ell}x^{k}y^{\ell}(a-k+b-\ell+1)(ct)^{a-% k+b-\ell}\frac{\mu_{a-k,b-\ell}^{W}}{2\Delta\varphi}\eta_{a-k+b-\ell+1}$	(12)
	$\displaystyle\overset{\mathbf{x}=0}{=}(a+b+1)(ct)^{a+b}\frac{\mu_{ab}^{W}}{2% \Delta\varphi}\eta_{a+b+1}$	(13)
$\displaystyle\mathbf{v}^{W}(t,\mathbf{x})$	$\displaystyle=-\frac{1}{r}\partial_{r}\left(r^{2}M[p_{0}\mathbf{n}](\mathbf{x}% ,r)\right)\Big{\|}_{r=ct}$	(14)
	$\displaystyle\overset{\eqref{eq:sphericalmeanscalarfunction}}{=}-\frac{1}{2% \Delta\varphi}\sum_{k=0}^{a}\sum_{\ell=0}^{b}{a\choose k}{b\choose\ell}x^{k}y^% {\ell}(a-k+b-\ell+2)(ct)^{a-k+b-\ell}$	(15)
	$\displaystyle\phantom{mmmmmmmmmmmmmmmmmmmmmmmm}\times\left(\begin{array}[]{c}% \mu_{a-k+1,b-\ell}^{W}\\ \mu_{a-k,b-\ell+1}^{W}\\ \end{array}\right)\eta_{a-k+b-\ell+2}$	(18)
	$\displaystyle\overset{\mathbf{x}=0}{=}-\frac{1}{2\Delta\varphi}(a+b+2)(ct)^{a+% b}\left(\begin{array}[]{c}\mu_{a+1,b}^{W}\\ \mu_{a,b+1}^{W}\\ \end{array}\right)\eta_{a+b+2}$	(21)

2.2 Initial data in $u$

Consider now $\mathbf{v}_{0}=(x^{a}y^{b},0,0)$ . Then, the solution reads

$\displaystyle p^{W}(t,\mathbf{x})$	$\displaystyle=-\frac{1}{r}\partial_{r}\left(r^{2}M[\mathbf{n}\cdot\mathbf{v}_{% 0}](\mathbf{x},r)\right)\Big{\|}_{r=ct}$	(22)
	$\displaystyle\overset{\eqref{eq:sphericalmeanvectorxspr}}{=}-\frac{1}{2\Delta% \varphi}\sum_{k=0}^{a}\sum_{\ell=0}^{b}{a\choose k}{b\choose\ell}x^{k}y^{\ell}% (a-k+b-\ell+2)(ct)^{a-k+b-\ell}\mu_{a-k+1,b-\ell}^{W}\eta_{a-k+b-\ell+2}$	(23)
	$\displaystyle\overset{\mathbf{x}=0}{=}-\frac{1}{2\Delta\varphi}(a+b+2)(ct)^{a+% b}\mu_{a+1,b}^{W}\eta_{a+b+2}$	(24)
$\displaystyle\mathbf{v}^{W}(t,\mathbf{x})$	$\displaystyle=\frac{2}{3}\mathbf{v}_{0}(\mathbf{x})+\partial_{r}\left(rM[(% \mathbf{v}_{0}\cdot\mathbf{n})\mathbf{n}](\mathbf{x},r)\right)$
	$\displaystyle-M\left[\mathbf{v}_{0}-3(\mathbf{v}_{0}\cdot\mathbf{n})\mathbf{n}% \right](\mathbf{x},r)-\int_{0}^{ct}\mathrm{d}r\frac{1}{r}M\left[\mathbf{v}_{0}% -3(\mathbf{v}_{0}\cdot\mathbf{n})\mathbf{n}\right](\mathbf{x},r)\Big{\|}_{r=ct}$	(25)
	$\displaystyle\overset{\eqref{eq:sphericalmeanscalarfunction},\eqref{eq:% sphericalmeannnx}}{=}\frac{2}{3}\left(\begin{array}[]{c}x^{a}y^{b}\\ 0\\ \end{array}\right)$	(28)
	$\displaystyle+\frac{1}{2\Delta\varphi}\sum_{k=0}^{a}\sum_{\ell=0}^{b}{a\choose k% }{b\choose\ell}x^{k}y^{\ell}(a-k+b-\ell+1)(ct)^{a-k+b-\ell}\left(\begin{array}% []{c}\mu_{a-k+2,b-\ell}^{W}\\ \mu_{a-k+1,b-\ell+1}^{W}\\ \end{array}\right)\eta_{a-k+b-\ell+3}$	(31)
	$\displaystyle-\frac{1}{2\Delta\varphi}\sum_{k=0}^{a}\sum_{\ell=0}^{b}{a\choose k% }{b\choose\ell}x^{k}y^{\ell}\left(\begin{array}[]{c}\mu_{a-k,b-\ell}^{W}\eta_{% a-k+b-\ell+1}-3\mu_{a-k+2,b-\ell}^{W}\eta_{a-k+b-\ell+3}\\ -3\mu_{a-k+1,b-\ell+1}^{W}\eta_{a-k+b-\ell+3}\\ \end{array}\right)$	(34)
	$\displaystyle\phantom{mmmmmmmmmmmmmmm}\times\left((ct)^{a-k+b-\ell}+\int_{0}^{% ct}\mathrm{d}r\,r^{a-k+b-\ell-1}\right)$

The integral $\int_{0}^{ct}\mathrm{d}r\,r^{a-k+b-\ell-1}$ is

\displaystyle\frac{(ct)^{a-k+b-\ell}}{a-k+b-\ell}\qquad\text{if }a-k+b-\ell\neq 0

(35)

The condition is violated when $k=a$ and $\ell=b$ . One can show that for continuous initial data, the term $\int_{0}^{ct}\mathrm{d}r\frac{1}{r}M\left[\mathbf{v}_{0}-3(\mathbf{v}_{0}\cdot% \mathbf{n})\mathbf{n}\right](\mathbf{x},r)$ is finite and in this case one can ignore the contribution from the case $a-k+b-\ell=0$ . In general, in particular for discontinuous data, a (logarithmic) singularity is part of the solution (see [BK22]). Thus,

$\displaystyle\mathbf{v}^{W}(t,\mathbf{x})$	$\displaystyle=\frac{2}{3}\left(\begin{array}[]{c}x^{a}y^{b}\\ 0\\ \end{array}\right)$	(38)
	$\displaystyle+\frac{1}{2\Delta\varphi}\sum_{k=0}^{a}\sum_{\ell=0}^{b}{a\choose k% }{b\choose\ell}x^{k}y^{\ell}(ct)^{a-k+b-\ell}$
	$\displaystyle\phantom{mmm}\times\left[(a-k+b-\ell+1)\left(\begin{array}[]{c}% \mu_{a-k+2,b-\ell}^{W}\\ \mu_{a-k+1,b-\ell+1}^{W}\\ \end{array}\right)\eta_{a-k+b-\ell+3}\right.$	(41)
	$\displaystyle\phantom{mmmmm}\left.-\left(\begin{array}[]{c}\mu_{a-k,b-\ell}^{W% }\eta_{a-k+b-\ell+1}-3\mu_{a-k+2,b-\ell}^{W}\eta_{a-k+b-\ell+3}\\ -3\mu_{a-k+1,b-\ell+1}^{W}\eta_{a-k+b-\ell+3}\\ \end{array}\right)\left(1+\frac{1}{a-k+b-\ell}\right)\right]$	(44)

where it is understood that the fraction in the last term is only to be inculded if it is defined.

If $\mathbf{x}=0$ , the expression simplifies:

	$\displaystyle\mathbf{v}^{W}(t,0)$	$\displaystyle=\frac{2}{3}\left(\begin{array}[]{c}x^{a}y^{b}\\ 0\\ \end{array}\right)+\frac{1}{2\Delta\varphi}(ct)^{a+b}$		(47)
		$\displaystyle\times\left[(a+b+1)\left(\begin{array}[]{c}\mu_{a+2,b}^{W}\\ \mu_{a+1,b+1}^{W}\\ \end{array}\right)\eta_{a+b+3}-\left(\begin{array}[]{c}\mu_{a,b}^{W}\eta_{a+b+% 1}-3\mu_{a+2,b}^{W}\eta_{a+b+3}\\ -3\mu_{a+1,b+1}^{W}\eta_{a+b+3}\\ \end{array}\right)\left(1+\frac{1}{a+b}\right)\right]$		(52)

where, again, the last fraction is only included if it is defined.

3 An Active Flux method for the Euler equations

3.1 General remarks

On a Cartesian grid with cells

\displaystyle C_{ij}=\left[x_{i-\frac{1}{2}},x_{i+\frac{1}{2}}\right]\times% \left[y_{j-\frac{1}{2}},y_{j+\frac{1}{2}}\right]

(53)

the degrees of freedom of Active Flux are (see also Figure 1) cell averages

\displaystyle\bar{q}_{ij}(t)\simeq\frac{1}{|C_{ij}|}\int_{C_{ij}}\mathrm{d}% \mathbf{x}\,q(t,\mathbf{x})

(54)

and point values at nodes

\displaystyle q_{i+\frac{1}{2},j+\frac{1}{2}}(t)\simeq q(t,x_{i+\frac{1}{2}},y% _{j+\frac{1}{2}})

(55)

and at edge midpoints

\displaystyle q_{i+\frac{1}{2},j}(t)

\displaystyle\simeq q(t,x_{i+\frac{1}{2}},y_{j})

\displaystyle q_{i,j+\frac{1}{2}}(t)

\displaystyle\simeq q(t,x_{i},y_{j+\frac{1}{2}})

(56)

Refer to caption — Figure 1: Degrees of freedom of Active Flux: Averages and point value at nodes and edge midpoints.

Time is discretized at instants $t^{n}$ , $n\in\mathbb{N}_{0}$ , with $\Delta t$ a generic name for the difference $t^{n+1}-t^{n}$ . The values of the degrees of freedom at instant $t^{n}$ are denoted by a superscript $n$ .

A reconstruction in each cell serves as initial data for the evolution operator. The classical reconstruction in Active Flux agrees pointwise with the given point values, and its average equals the given cell average. With 9 conditions, a biparabolic reconstruction in

\displaystyle\mathrm{span}(1,x,x^{2},y,xy,xy^{2},y^{2}xy^{2},y^{2}y^{2})

(57)

is natural (as used in [BHKR19, HKS19], where the polynomial is given explicitly). Since the restriction of a biparabolic function on one of the (axis-aligned) edges is a parabola, and since it is uniquely determined by the three values present along each edge, the reconstructions form globally continuous initial data. As is described later, a biparabolic reconstruction is still used here, but computed slightly differently.

3.2 Update of the average

Recall that the update of the average $\bar{q}_{ij}:=\frac{1}{|C_{ij}|}\int\mathrm{d}\mathbf{x}\,q(t,\mathbf{x})$ is obtained immediately by integrating the conservation law

	$\displaystyle\partial_{t}q+\nabla\cdot\mathbf{f}(q)$	$\displaystyle=0$	$\displaystyle q$	$\displaystyle\colon\mathbb{R}^{+}_{0}\times\mathbb{R}^{2}\to\mathbb{R}^{m}$		(58)
	$\displaystyle\mathbf{f}$	$\displaystyle=(f^{x},f^{y})\quad f^{x},f^{y}\colon\mathbb{R}^{m}\to\mathbb{R}^% {m}$				(59)

over the cell $C_{ij}$ and over the time step $[t^{n},t^{n+1}]$ and replacing the integrals by space-time Simpson-rule quadratures (see e.g. [ER13]):

\displaystyle\frac{\bar{q}_{ij}^{n+1}-\bar{q}_{ij}^{n}}{\Delta t}

\displaystyle=-\frac{1}{\Delta x}\left(\hat{f}^{x}_{i+\frac{1}{2},j}-\hat{f}^{% x}_{i-\frac{1}{2},j}\right)-\frac{1}{\Delta y}\left(\hat{f}^{y}_{i,j+\frac{1}{% 2}}-\hat{f}^{y}_{i,j-\frac{1}{2}}\right)\simeq-\int_{\partial C_{ij}}\mathrm{d% }\mathbf{x}\,\mathbf{n}\cdot\mathbf{f}(q)

(60)

with

$\displaystyle\hat{f}^{x}_{i+\frac{1}{2},j}$	$\displaystyle:=\frac{1}{6}\left(\frac{f^{x}(q^{n}_{i+\frac{1}{2},j-\frac{1}{2}% })}{6}+\frac{2f^{x}(q^{n}_{i+\frac{1}{2},j})}{3}+\frac{f^{x}(q^{n}_{i+\frac{1}% {2},j+\frac{1}{2}})}{6}\right)$	(61)
	$\displaystyle\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!+\frac{2}{3}\left(\frac{f^{x}(q^{% n+\frac{1}{2}}_{i+\frac{1}{2},j-\frac{1}{2}})}{6}+\frac{2f^{x}(q^{n+\frac{1}{2% }}_{i+\frac{1}{2},j})}{3}+\frac{f^{x}(q^{n+\frac{1}{2}}_{i+\frac{1}{2},j+\frac% {1}{2}})}{6}\right)+\frac{1}{6}\left(\frac{f^{x}(q^{n+1}_{i+\frac{1}{2},j-% \frac{1}{2}})}{6}+\frac{2f^{x}(q^{n+1}_{i+\frac{1}{2},j})}{3}+\frac{f^{x}(q^{n% +1}_{i+\frac{1}{2},j+\frac{1}{2}})}{6}\right)$
	$\displaystyle\simeq\frac{1}{\Delta y}\int_{y_{j-\frac{1}{2}}}^{y_{j+\frac{1}{2% }}}\mathrm{d}y\frac{1}{\Delta t}\int_{t^{n}}^{t^{n+1}}\mathrm{d}t\,f^{x}(q(x_{% i+\frac{1}{2}},y,t))$	(62)
$\displaystyle\hat{f}^{y}_{i,j+\frac{1}{2}}$	$\displaystyle:=\frac{1}{6}\left(\frac{f^{y}(q^{n}_{i-\frac{1}{2},j+\frac{1}{2}% })}{6}+\frac{2f^{y}(q^{n}_{i,j+\frac{1}{2}})}{3}+\frac{f^{y}(q^{n}_{i+\frac{1}% {2},j+\frac{1}{2}})}{6}\right)$	(63)
	$\displaystyle\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!+\frac{2}{3}\left(\frac{f^{y}(q^{% n+\frac{1}{2}}_{i-\frac{1}{2},j+\frac{1}{2}})}{6}+\frac{2f^{y}(q^{n+\frac{1}{2% }}_{i,j+\frac{1}{2}})}{3}+\frac{f^{y}(q^{n+\frac{1}{2}}_{i+\frac{1}{2},j+\frac% {1}{2}})}{6}\right)+\frac{1}{6}\left(\frac{f^{y}(q^{n+1}_{i-\frac{1}{2},j+% \frac{1}{2}})}{6}+\frac{2f^{y}(q^{n+1}_{i,j+\frac{1}{2}})}{3}+\frac{f^{y}(q^{n% +1}_{i+\frac{1}{2},j+\frac{1}{2}})}{6}\right)$
	$\displaystyle\simeq\frac{1}{\Delta x}\int_{x_{i-\frac{1}{2}}}^{x_{i+\frac{1}{2% }}}\mathrm{d}x\frac{1}{\Delta t}\int_{t^{n}}^{t^{n+1}}\mathrm{d}t\,f^{y}(q(x,y% _{j+\frac{1}{2}},t))$	(64)

Since this step requires the new point values at $t^{n}+\frac{\Delta t}{2}$ and $t^{n+1}=t^{n}+\Delta t$ , it is performed as the very last.

3.3 Update of the point values

Recall the Euler equations in primitive variables

$\displaystyle\partial_{t}\rho+\mathbf{v}\cdot\nabla\rho+\rho\nabla\cdot\mathbf% {v}$	$\displaystyle=0$	(65)
$\displaystyle\partial_{t}\mathbf{v}+(\mathbf{v}\cdot\nabla)\mathbf{v}+\frac{% \nabla p}{\rho}$	$\displaystyle=0$	(66)
$\displaystyle\partial_{t}p+\mathbf{v}\cdot\nabla p+\rho c^{2}\nabla\cdot% \mathbf{v}$	$\displaystyle=0$	(67)

Contrary to e.g. [ABK25] here the following special choices are made:

1.

The first ingredient of the presented approach is the usage of primitive variables for the point values.
2.

The second is a reconstruction in primitive variables. For point values, the transformation between conserved and primitive variables is immediate. Therefore, here, after transforming the point values to conserved variables, and given the average over the cell (in conserved variables, too), a two-dimensional Simpson rule is inverted to find the value at cell center (in conserved variables). This latter is finally transformed to primitive variables and a biparabolic reconstruction is computed.

An operator splitting is employed, by first solving

$\displaystyle\partial_{t}\rho+\rho\nabla\cdot\mathbf{v}$	$\displaystyle=0$	(68)
$\displaystyle\partial_{t}\mathbf{v}+\frac{\nabla p}{\rho}$	$\displaystyle=0$	(69)
$\displaystyle\partial_{t}p+\rho c^{2}\nabla\cdot\mathbf{v}$	$\displaystyle=0$	(70)

in a locally linearized manner, i.e. actually solving

	$\displaystyle\partial_{t}\mathbf{v}+c_{0}\nabla\left(\frac{p}{\rho_{0}c_{0}}\right)$	$\displaystyle=0$		(71)
	$\displaystyle\partial_{t}\left(\frac{p}{\rho_{0}c_{0}}\right)+c_{0}\nabla\cdot% \mathbf{v}$	$\displaystyle=0$		(72)

with $c_{0}$ , $\rho_{0}$ evaluations of the speed of sound and of the density at the location of the respective point value, and obtaining the density from

\displaystyle\partial_{t}(\rho c_{0}^{2}-p)

\displaystyle=0

(73)

To solve (71)–(72), the exact evolution operator for linear acoustics is used, as described in Section 2. The biparabolic reconstruction (involving values at time $t^{n}$ ) serves as the initial datum.

In a second step,

$\displaystyle\partial_{t}\rho+\mathbf{v}\cdot\nabla\rho$	$\displaystyle=0$	(74)
$\displaystyle\partial_{t}\mathbf{v}+(\mathbf{v}\cdot\nabla)\mathbf{v}$	$\displaystyle=0$	(75)
$\displaystyle\partial_{t}p+\mathbf{v}\cdot\nabla p$	$\displaystyle=0$	(76)

is solved. Local linearization would amount to evolving the data according to

\displaystyle\mathbf{x}\mapsto\mathbf{x}-\mathbf{v}_{0}t

(77)

with, again, $\mathbf{v}_{0}$ the evaluation of the velocity at the location of the respective point value. Here instead, a fully third-order algorithm is used to solve (74)–(76), which is discussed in Section 3.4. The same biparabolic reconstruction (involving values at time $t^{n}$ ) serves as the initial datum.

4.

The two operators are combined in additive operator splitting.

Remark 3.1.

Define solution operators $S_{1},S_{2}\colon\mathbb{R}^{M}\times[0,T]\to\mathbb{R}^{M}$ , such that $S_{k}(q_{0},t)$ , $k=1,2$ is the solution at time $t$ of $\text{IVP}_{k}$ (obtained by splitting some IVP into two operators) with initial data $q_{0}\in\mathbb{R}^{M}$ . Then, recall that multiplicative operator splitting is

\displaystyle S_{2}(S_{1}(q_{0},t),t)

(78)

while

\displaystyle S_{2}(q_{0},t)+S_{1}(q_{0},t)-q_{0}

(79)

is additive operator splitting. Multiplicative operator splitting would require the solution of the first sub-operator to serve as initial data of the second. However, it is not possible to update the average individually with any of the two sub-operators that are considered here, since the presented splitting into acoustics and advection yields non-conservative sub-operators. Additive operator splitting can be applied to the same initial data and allows to circumvent this difficulty. The update of the average then happens using the unsplit equations.

Remark 3.2.

The cell average is evolved in conserved variables by evaluating the flux

\mathbf{f}(q_{\text{prim}})=\left(\begin{array}[]{c}\rho\mathbf{v}\\ \rho\mathbf{v}\otimes\mathbf{v}+p\mathbb{1}\\ \mathbf{v}\left(\frac{\gamma p}{\gamma-1}+\frac{1}{2}\rho|\mathbf{v}|^{2}% \right)\\ \end{array}\right)

using the primitive variables $(\rho,\mathbf{v},p)$ stored at the point values. The method thus is conservative.

Remark 3.3.

The transformation of the point values back and forth between primitive and conserved variables does not entail a big cost. Using primitive variables for the reconstruction is a common approach for Finite Volume methods. More importantly, this extra effort is accepted here for the greater advantage of a significantly simpler usage of the acoustic evolution operator. In principle, however, it seems possible to perform essentially the same algorithm in conserved variables alone.

Due to the usage of the exact evolution operator for acoustics, the linearization of the present scheme amounts to a combination of advection and acoustics, the stability of which has been analyzed in [CHK21]. The shortest separation between two degrees of freedom is $\mathrm{min}(\Delta x,\Delta y)/2$ . The method thus has a maximal CFL number of 0.5.

3.4 Order of accuracy

3.4.1 Overview

It has been shown in [Bar21a] that local linearization (i.e. the evaluation of the Jacobians of the sub-operators (71)–(72), (74)–(76) at the point value that one seeks to update) leads to a second-order accurate method only. In the same work it has been shown that fixpoint iteration on the characteristic equation leads to high-order results for scalar conservation laws. This is because even if they are nonlinear, characteristics remain straight lines. This is no longer true for systems, and [Bar21a] showed that a naive fixpoint iteration does not improve the order beyond second.

For systems, the characteristics are curved since their speed depends on all characteristic variables²²2Riemann invariants are sometimes mistaken for characteristic variables. Riemann invariants have the same value across a simple wave (e.g. contact discontinuity or rarefaction) and they exist for e.g. the Euler equations. Characteristic variables $Q$ are those for which the hyperbolic system can be written in the form $\partial_{t}Q_{k}+\lambda_{k}(Q_{1},\ldots,Q_{N})\partial_{x}Q_{k}=0$ . If they exist, they remain constant on “their” characteristic, which is a stronger statement not restricted to simple waves, but holding true everywhere in the region of smoothness of the solution. However, they exist only if an integrability condition is verified (the matrix $\partial Q/\partial q$ needs to be commuted with the derivative $\partial_{x}$ ). For linear systems, this integrability condition is void, and hence generically across a $k$ -wave all characteristic variables but the $k$ -th are Riemann invariants. This is also the case for nonlinear problems if the integrability condition is true (any $2\times 2$ system, for example, which includes shallow water and isentropic Euler). For the full Euler equations, the integrability condition is false, and characteristic variables do not exist. The explanation therefore is to be understood in the simplified setting of hyperbolic systems with characteristic variables., but only one of them is constant along the characteristic. The speed is thus changing in function of the values brought in by the other characteristic families.

The correct way to achieve a third-order approximate evolution operator in one spatial dimension is thus characterized by a mixing of eigenvalues. For example, for a $2\times 2$ system

	$\displaystyle\partial_{t}q_{1}+\lambda_{1}(q_{1},q_{2})\partial_{x}q_{1}$	$\displaystyle=0$		(80)
	$\displaystyle\partial_{t}q_{2}+\lambda_{2}(q_{1},q_{2})\partial_{x}q_{2}$	$\displaystyle=0$		(81)

with initial data $q_{k}(0,x)=q_{0,k}(x)$ , local linearization produces $x-t\lambda_{1}\big{(}q_{1,0}(x),q_{2,0}(x)\big{)}$ as the estimate of the foot of the 1-characteristic passing through $x$ at time $t$ . The estimate to next order of accuracy is ([Bar21a])

\displaystyle x-t\lambda_{1}\left(q_{1,0}\big{(}x-t\lambda_{1}(x)\big{)},q_{2,% 0}\left(x-t\frac{\lambda_{1}(x)+\lambda_{2}(x)}{2}\right)\right)

(82)

Here, $\lambda_{k}(x)$ is short-hand for $\lambda_{k}(q_{1,0}(x),q_{2,0}(x))$ . Simple characteristic-wise fixpoint iteration will never produce signals traveling at the average speed $\frac{\lambda_{1}(x)+\lambda_{2}(x)}{2}$ . For more complicated hyperbolic systems than (80)-(81) one additionally needs to take into account the changing eigenvectors, see [Bar21a] for the general algorithm in 1-d. Therein, this algorithm has been successfully applied to the Euler equations.

For nonlinear systems in multiple spatial dimensions the question of high-order evolution operators beyond local linearization is open. In practice, local linearization is employed (e.g. in [FR15, CHLM24]), thus formally reducing the order to second. The nonlinear version of advection, however, is amenable to higher-order evolution as discussed in [Mae17]. Here, a different and inexpensive modification for solve (74)–(75) is used, as described below.

Overall, since the reconstruction and all quadratures are third-order, the presented method is third-order accurate in space. To yield an Active Flux method that is also third-order accurate in time, the point updates have to have an error of at most $\mathcal{O}(t^{3})$ ([Bar21a]). While the advective operator is solved with this accuracy (see below), the locally linearized acoustics operator and the operator splitting both have errors $\mathcal{O}(t^{2})$ and thus make the method formally second-order in time. This, at one order higher, finds parallels in the current practice: For high-order methods such as WENO or DG, it is actually customary to save on the effort of an adequately high-order time integration and a method like RK3 is often used. As is shown below, experimental results on the presented method give rise to the conclusion that the global error is not dominated by its second-order component in typical numerical tests. In fact, basing one’s conclusion on certain experimental results one could easily mistake the method for a third-order accurate one.

3.4.2 A high-order evolution operator for nonlinear advection

Consider a hyperbolic system

\displaystyle\partial_{t}q+J^{x}(q)\partial_{x}q+J^{y}(q)\partial_{y}q=0

(83)

with initial data $q(0,\mathbf{x})=q_{0}(\mathbf{x})$ . We denote the components of the initial data as $q_{i,0}(\mathbf{x})$ .

Theorem 3.1.

If $J^{x}=v^{x}\mathbb{1}$ , $J^{y}=v^{y}\mathbb{1}$ and both $v_{x}$ and $v_{y}$ are components of $q$ (as in (74)–(75)), then the solution of (83) (if smooth) is

\displaystyle q(t,\mathbf{x})=q_{\text{loc.lin.}}(t,\mathbf{x})+t^{2}\nabla q_% {i,0}\cdot((\mathbf{v}_{0}\cdot\nabla)\mathbf{v}_{0})+\mathcal{O}(t^{3})

(84)

where $q_{\text{loc.lin.}}(t,\mathbf{x})$ is obtained by considering $J^{x}$ , $J^{y}$ fixed with whatever value they have at $\mathbf{x}$ at $t=0$ , and $\mathbf{v}_{0}=(v^{x}(0,\mathbf{x}),v^{y}(0,\mathbf{x}))$ .

Proof.

We find (summation over repeated indices):

$\displaystyle\partial_{t}q_{i}(t,\mathbf{x})$	$\displaystyle=-J^{x}_{ij}\partial_{x}q_{j}-J^{y}_{ij}\partial_{y}q_{j}$	(85)
$\displaystyle\partial_{t}^{2}q_{i}(t,\mathbf{x})$	$\displaystyle=-J^{x}_{ij}\partial_{x}\partial_{t}q_{j}-J^{y}_{ij}\partial_{y}% \partial_{t}q_{j}-\partial_{t}J^{x}_{ij}\partial_{x}q_{j}-\partial_{t}J^{y}_{% ij}\partial_{y}q_{j}$	(86)
	$\displaystyle=J^{x}_{ij}\partial_{x}(J^{x}_{jk}\partial_{x}q_{k}+J^{y}_{jk}% \partial_{y}q_{k})+J^{y}_{ij}\partial_{y}(J^{x}_{jk}\partial_{x}q_{k}+J^{y}_{% jk}\partial_{y}q_{k})$	(87)
	$\displaystyle\qquad+\frac{\partial J^{x}_{ij}}{\partial q_{k}}(J^{x}_{k\ell}% \partial_{x}q_{\ell}+J^{y}_{k\ell}\partial_{y}q_{\ell})\partial_{x}q_{j}+\frac% {\partial J^{y}_{ij}}{\partial q_{k}}(J^{x}_{k\ell}\partial_{x}q_{\ell}+J^{y}_% {k\ell}\partial_{y}q_{\ell})\partial_{y}q_{j}$
	$\displaystyle=J^{x}_{ij}\left(J^{x}_{jk}\partial_{x}^{2}q_{k}+J^{y}_{jk}% \partial_{x}\partial_{y}q_{k}+\frac{\partial J^{x}_{jk}}{\partial q_{\ell}}% \partial_{x}q_{\ell}\partial_{x}q_{k}+\frac{\partial J^{y}_{jk}}{\partial q_{% \ell}}\partial_{x}q_{\ell}\partial_{y}q_{k}\right)$	(88)
	$\displaystyle\qquad+J^{y}_{ij}\left(J^{x}_{jk}\partial_{x}\partial_{y}q_{k}+J^% {y}_{jk}\partial_{y}^{2}q_{k}+\frac{\partial J^{x}_{jk}}{\partial q_{\ell}}% \partial_{y}q_{\ell}\partial_{x}q_{k}+\frac{\partial J^{y}_{jk}}{\partial q_{% \ell}}\partial_{y}q_{\ell}\partial_{y}q_{k}\right)$
	$\displaystyle\qquad+\frac{\partial J^{x}_{ij}}{\partial q_{k}}(J^{x}_{k\ell}% \partial_{x}q_{\ell}+J^{y}_{k\ell}\partial_{y}q_{\ell})\partial_{x}q_{j}+\frac% {\partial J^{y}_{ij}}{\partial q_{k}}(J^{x}_{k\ell}\partial_{x}q_{\ell}+J^{y}_% {k\ell}\partial_{y}q_{\ell})\partial_{y}q_{j}$

The 4 terms not involving the derivatives of $J^{x}$ , $J^{y}$ are those appearing for a linear problem and hence also those produced by local linearization. The remaining terms are (renaming $j\leftrightarrow k$ in the last parts)

$\displaystyle S_{i}$	$\displaystyle:=J^{x}_{ij}\left(\frac{\partial J^{x}_{jk}}{\partial q_{\ell}}% \partial_{x}q_{\ell}\partial_{x}q_{k}+\frac{\partial J^{y}_{jk}}{\partial q_{% \ell}}\partial_{x}q_{\ell}\partial_{y}q_{k}\right)+J^{y}_{ij}\left(\frac{% \partial J^{x}_{jk}}{\partial q_{\ell}}\partial_{y}q_{\ell}\partial_{x}q_{k}+% \frac{\partial J^{y}_{jk}}{\partial q_{\ell}}\partial_{y}q_{\ell}\partial_{y}q% _{k}\right)$	(89)
	$\displaystyle\qquad+\frac{\partial J^{x}_{ik}}{\partial q_{j}}(J^{x}_{j\ell}% \partial_{x}q_{\ell}+J^{y}_{j\ell}\partial_{y}q_{\ell})\partial_{x}q_{k}+\frac% {\partial J^{y}_{ik}}{\partial q_{j}}(J^{x}_{j\ell}\partial_{x}q_{\ell}+J^{y}_% {j\ell}\partial_{y}q_{\ell})\partial_{y}q_{k}$
	$\displaystyle=\left(J^{x}_{ij}\frac{\partial J^{x}_{jk}}{\partial q_{\ell}}+% \frac{\partial J^{x}_{ik}}{\partial q_{j}}J^{x}_{j\ell}\right)\partial_{x}q_{% \ell}\partial_{x}q_{k}+\left(J^{x}_{ij}\frac{\partial J^{y}_{jk}}{\partial q_{% \ell}}+J^{y}_{ij}\frac{\partial J^{x}_{j\ell}}{\partial q_{k}}+\frac{\partial J% ^{x}_{i\ell}}{\partial q_{j}}J^{y}_{jk}+\frac{\partial J^{y}_{ik}}{\partial q_% {j}}J^{x}_{j\ell}\right)\partial_{x}q_{\ell}\partial_{y}q_{k}$	(90)
	$\displaystyle\qquad+\left(J^{y}_{ij}\frac{\partial J^{y}_{jk}}{\partial q_{% \ell}}+\frac{\partial J^{y}_{ik}}{\partial q_{j}}J^{y}_{j\ell}\right)\partial_% {y}q_{\ell}\partial_{y}q_{k}$

Consider now the special case $(J^{x})_{ij}=v^{x}\delta_{ij}$ , $(J^{y})_{ij}=v^{y}\delta_{ij}$ :

$\displaystyle S_{i}$	$\displaystyle=2v^{x}\frac{\partial v^{x}}{\partial q_{\ell}}\partial_{x}q_{% \ell}\partial_{x}q_{i}+\left(2v^{x}\frac{\partial v^{y}}{\partial q_{\ell}}% \partial_{x}q_{\ell}\partial_{y}q_{i}+2v^{y}\frac{\partial v^{x}}{\partial q_{% k}}\partial_{x}q_{i}\partial_{y}q_{k}\right)+2v^{y}\frac{\partial v^{y}}{% \partial q_{\ell}}\partial_{y}q_{\ell}\partial_{y}q_{i}$	(91)
	$\displaystyle=2(v^{x}\partial_{x}v^{x}+v^{y}\partial_{y}v^{x})\partial_{x}q_{i% }+2(v^{x}\partial_{x}v^{y}+v^{y}\partial_{y}v^{y})\partial_{y}q_{i}$	(92)
	$\displaystyle=2\nabla q_{i}\cdot((\mathbf{v}\cdot\nabla)\mathbf{v})$	(93)

At this point the advantage of primitive variables is clearly visible. The statement is obtained by realizing that the derivatives need to be evaluated at $t=0$ in the Taylor series of the solution. ∎

Theorem 3.2.

If $J^{x}=v^{x}\mathbb{1}$ , $J^{y}=v^{y}\mathbb{1}$ , then the solution of (83) with initial data $q(0,\mathbf{x})=q_{0}(\mathbf{x})$ is

\displaystyle q(t,\mathbf{x})=q_{0}\Big{(}\mathbf{x}-\mathbf{v}_{0}\big{(}x-% \mathbf{v}_{0}(\mathbf{x})t\big{)}t\Big{)}+\mathcal{O}(t^{3})

(94)

where $\mathbf{v}_{0}(\mathbf{x})=(v^{x}(0,\mathbf{x}),v^{y}(0,\mathbf{x}))$ .

This is the same simple iteration as is in use for scalar conservation laws.

Proof.

We need to verify the presence of the term $\Delta t^{2}\nabla q_{i}\Big{|}_{t=0}\cdot((\mathbf{v}_{0}\cdot\nabla)\mathbf{% v}_{0})$ .

$\displaystyle\tilde{q}_{i}(t,\mathbf{x})$	$\displaystyle:=q_{i,0}\Big{(}\mathbf{x}-\mathbf{v}_{0}\big{(}\mathbf{x}-% \mathbf{v}_{0}(\mathbf{x})t\big{)}t\Big{)}$	(95)
$\displaystyle\partial_{t}\tilde{q}_{i}(t,\mathbf{x})$	$\displaystyle=\partial_{k}q_{i,0}\Big{(}\mathbf{x}-\mathbf{v}_{0}\big{(}% \mathbf{x}-\mathbf{v}_{0}(\mathbf{x})t\big{)}t\Big{)}\partial_{t}\big{(}x_{k}-% v^{k}_{0}(\mathbf{x}-\mathbf{v}_{0}(\mathbf{x})t)t\big{)}$	(96)
	$\displaystyle=-\partial_{k}q_{i,0}\Big{(}\mathbf{x}-\mathbf{v}_{0}\big{(}% \mathbf{x}-\mathbf{v}_{0}(\mathbf{x})t\big{)}t\Big{)}\Big{(}-\partial_{\ell}v^% {k}_{0}\big{(}\mathbf{x}-\mathbf{v}_{0}(\mathbf{x})t\big{)}v^{\ell}_{0}(% \mathbf{x})t+v^{k}_{0}\big{(}\mathbf{x}-\mathbf{v}_{0}(\mathbf{x})t\big{)}\Big% {)}$	(97)
$\displaystyle\partial_{t}^{2}\tilde{q}_{i}(t,\mathbf{x})\Big{\|}_{t=0}$	$\displaystyle=\partial_{\ell}\partial_{k}q_{i,0}(\mathbf{x})v^{\ell}_{0}(% \mathbf{x})v^{k}_{0}(\mathbf{x})+2\partial_{k}q_{i,0}(\mathbf{x})v^{\ell}_{0}(% \mathbf{x})\partial_{\ell}v^{k}_{0}(\mathbf{x})$	(98)

Here $\partial_{k}$ denotes $\partial_{0}\equiv\partial_{x}$ , $\partial_{1}\equiv\partial_{y}$ , and analogously for the components $v^{k}$ of $\mathbf{v}$ . The first term is that of local linearization, and the second is precisely the extra term from Theorem 3.1. ∎

In the Active Flux method therefore, instead of solving (74)–(75) via local linearization, two iterations are performed: One obtains first

\displaystyle\tilde{\mathbf{x}}:=\mathbf{x}-\mathbf{v}_{0}\big{(}\mathbf{x}-% \mathbf{v}_{0}(\mathbf{x})t\big{)}t

(99)

and then, the reconstruction of $\rho,\mathbf{v},p$ is evaluated at $\tilde{\mathbf{x}}$ . Recall that this is a very special case, and does not lead to an increase in order in situations of non-commuting Jacobian matrices, i.e. it is not applicable to the acoustic operator. The difference of the present approach to that of [RMF18] is the absence of spatial derivatives. One can draw the following analogy to methods for ODEs: While the approach of [Mae17, RMF18]is similar to the Lax-Wendroff procedure converting higher time derivatives to space derivatives, the approach of [Bar21a] as well as the prescription (99) are more in the spirit of a Runge-Kutta method that generates the right terms in the Taylor series through nested evaluations.

3.5 Bound preservation

3.5.1 Point value update

Since the new point values are required to obtain the fluxes for the average update, bound preservation for the point values is ensured first. Action is deemed necessary, if the new value is NaN (which can happen if the speed of sound has been evaluated at a location with negative density or pressure) or if the new value has a pressure or density that is lower then $\epsilon$ . Here, $\epsilon=10^{-10}$ is chosen. If limiting is needed, a local Lax-Friedrichs method based on the three values

\displaystyle\bar{q}_{ij}

\displaystyle q_{i+\frac{1}{2},j}

\displaystyle\bar{q}_{i+1,j}

(100)

is used (and analogously for the other direction). This amounts to defining control volumes of size $\frac{\Delta x}{2}\times\frac{\Delta y}{2}$ centered at each degree of freedom (see Figure 2), and identifying the average with the point value at cell center. Thus, the bound preserving point value at edge midpoint is

\displaystyle q_{i+\frac{1}{2},j}^{n+1,\text{lim}}=q_{i+\frac{1}{2},j}^{n}-% \frac{\Delta t}{\Delta x/2}\left(\frac{f^{x}(\bar{q}^{n}_{i+1,j})-f^{x}(\bar{q% }^{n}_{ij})}{2}-\frac{\lambda^{\text{HLL}}}{2}(\bar{q}^{n}_{i+1,j}-2q^{n}_{i+% \frac{1}{2},j}+\bar{q}^{n}_{ij})\right)

(101)

An analogous approach is used for the perpendicular edges, while for the node, four neighbours contribute:

	$\displaystyle q_{i+\frac{1}{2},j+\frac{1}{2}}^{n+1,\text{lim}}=q_{i+\frac{1}{2% },j+\frac{1}{2}}^{n}$	$\displaystyle-\frac{\Delta t}{\Delta x/2}\left(\frac{f^{x}(\bar{q}^{n}_{i+1,j+% \frac{1}{2}})-f^{x}(\bar{q}^{n}_{i,j+\frac{1}{2}})}{2}-\frac{\lambda^{\text{% HLL}}}{2}(\bar{q}^{n}_{i+1,j+\frac{1}{2}}-2q^{n}_{i+\frac{1}{2},j+\frac{1}{2}}% +\bar{q}^{n}_{i,j+\frac{1}{2}})\right)$		(102)
		$\displaystyle-\frac{\Delta t}{\Delta y/2}\left(\frac{f^{y}(\bar{q}^{n}_{i+% \frac{1}{2},j+1})-f^{y}(\bar{q}^{n}_{i+\frac{1}{2},j})}{2}-\frac{\lambda^{% \text{HLL}}}{2}(\bar{q}^{n}_{i+\frac{1}{2},j+1}-2q^{n}_{i+\frac{1}{2},j+\frac{% 1}{2}}+\bar{q}^{n}_{i+\frac{1}{2},j})\right)$

Since the point values are in primitive variables, prior to this calculation they need to be converted to conserved ones, and converted back afterwards. Observe that this is an a posteriori limiting, and thus very inexpensive since it occurs rarely.

Finally, for all edges the states

\displaystyle q_{i+\frac{1}{2},j}^{\text{HLL}}=\frac{\bar{q}_{i+1,j}+\bar{q}_{% ij}}{2}-\frac{1}{2\lambda^{\text{HLL}}_{i+\frac{1}{2},j}}\left(f^{x}(\bar{q}_{% i+1,j})-f^{x}(\bar{q}_{ij})\right)

(103)

and the fluxes

\displaystyle\hat{f}_{i+\frac{1}{2},j}^{\text{HLL}}=\frac{f^{x}(\bar{q}_{i+1,j% })+f^{x}(\bar{q}_{ij})}{2}-\frac{\lambda^{\text{HLL}}_{i+\frac{1}{2},j}}{2}% \left(\bar{q}_{i+1,j}-\bar{q}_{ij}\right)

(104)

are computed (analogously for perpendicular edges). These states and fluxes correspond to considering only the grid of the averages, which has spacing $\Delta x$ , and they are needed later for the bound preserving limiting of the averages.

As a practical simplification, the speed $\lambda^{\text{HLL}}_{i+\frac{1}{2},j}$ is computed using the point value at time $t^{n}$ ; a better way is described in [DBK25] and can be applied here as well. It is also taken as the maximum of the wave speed in both directions.

3.5.2 Average update

The procedure adopted here is the one from [DBK25], which goes back to [GP16, Kuz20, Haj21] as well as [PS96, KLT12]: Given a high-order flux $\hat{f}_{i+\frac{1}{2},j}$ (from (61)) and a low-order bound-preserving flux $\hat{f}^{\text{HLL}}_{i+\frac{1}{2},j}$ (from (104)), a convex combination

	$\displaystyle\vartheta\hat{f}_{i+\frac{1}{2},j}+(1-\vartheta)\hat{f}^{\text{% HLL}}_{i+\frac{1}{2},j}$	$\displaystyle=\hat{f}^{\text{HLL}}_{i+\frac{1}{2},j}+\vartheta\Delta\hat{f}_{i% +\frac{1}{2},j}$		(105)
	$\displaystyle\Delta\hat{f}_{i+\frac{1}{2},j}$	$\displaystyle:=\hat{f}_{i+\frac{1}{2},j}-\hat{f}^{\text{HLL}}_{i+\frac{1}{2},j}$		(106)

is found that is optimal, in the sense that it is as close to high-order as possible while maintaining bound-preservation. Since the method is 1-stage, there is no need to rely on an SSP property of the time integrator.

The condition for bound preservation is that the states

\displaystyle q_{i+\frac{1}{2},j}^{\text{HLL}}\pm\vartheta\frac{\Delta\hat{f}_% {i+\frac{1}{2},j}}{\lambda^{\text{HLL}}_{i+\frac{1}{2},j}}

(107)

have positive density and pressure (see [DBK25] for a derivation). For the density, one thus imposes

\displaystyle\vartheta\Delta\hat{f}_{i+\frac{1}{2},j}=\begin{cases}\min(\Delta% \hat{f}_{i+\frac{1}{2},j},(\rho_{i+\frac{1}{2},j}^{\text{HLL}}-\epsilon)% \lambda^{\text{HLL}}_{i+\frac{1}{2},j})&\Delta\hat{f}_{i+\frac{1}{2},j}>0\\ \max(\Delta\hat{f}_{i+\frac{1}{2},j},(\epsilon-\rho_{i+\frac{1}{2},j}^{\text{% HLL}})\lambda^{\text{HLL}}_{i+\frac{1}{2},j})&\text{otherwise}\end{cases}

(108)

and computes the new flux of just the density according to (105).

Next, the positivity of pressure leads to the inequalities

\displaystyle\vartheta^{2}A\pm\vartheta B<C

(109)

which are ensured by enforcing (through $\vartheta\in[0,1]$ )

\displaystyle\vartheta\big{(}\max(0,A)+|B|\big{)}<C

(110)

instead. The values of $A,B,C$ are

$\displaystyle A$	$\displaystyle=\frac{1}{2}\left((\Delta f^{m_{x}})^{2}+(\Delta f^{m_{y}})^{2}% \right)-\Delta f^{e}\Delta f^{\rho}$	(111)
$\displaystyle B$	$\displaystyle=\left(-{\mathbf{m}}^{\text{HLL}}\cdot\Delta f^{\mathbf{m}}+% \Delta f^{e}\rho^{\text{HLL}}+\Delta f^{\rho}e^{\text{HLL}}-\Delta f^{\rho}% \frac{\epsilon}{\gamma-1}\right)\lambda^{\text{HLL}}$	(112)
$\displaystyle C$	$\displaystyle=\left(-\rho^{\text{HLL}}\frac{\epsilon}{\gamma-1}-\frac{1}{2}% \big{(}(m_{x}^{\text{HLL}})^{2}+(m_{y}^{\text{HLL}})^{2}\big{)}+e^{\text{HLL}}% \rho^{\text{HLL}}\right)(\lambda^{\text{HLL}})^{2}$	(113)

where, to save a tree, the indices have been dropped, and $\mathbf{m}$ denotes the momentum. With the $\vartheta$ thus obtained, (105) is used once again to compute the new flux of all the quantities.

Limiting with $\vartheta=0$ is applied straight away if the high-order flux turns out to be NaN.

Deliberately no further approaches to limiting are undertaken, since the aim of this paper is to show how the pure method works, with bound preservation necessary to keep the simulation running. It is clear that it would be possible to endow the method with e.g. shock sensors to limit if necessary in a practical situation.

4 Numerical results

All numerical tests are performed with bound preservation turned on, and a CFL number of 0.45 with respect to $\mathrm{min}(\Delta x,\Delta y)$ . Recall that the smallest separation of two degrees of freedom is half the edge length, and thus $\text{CFL}=0.5$ is the stability limit ([CHK21]). For comparison, many of the test cases are those considered in [ABK25] for the semi-discrete Active Flux method. $\gamma=1.4$ is used throughout.

4.1 Convergence tests

First, convergence tests are performed in order to assess the effect of the combination of a third-order space discretization, and a mixed-order time discretization.

4.1.1 Contact wave

This test is from [CHLM24]: The velocity and the pressure are constant:

$\displaystyle\rho_{0}(x,y)$	$\displaystyle=\frac{5}{2}\exp\left(-40((x-x_{0})^{2}+(y-y_{0})^{2})\right)+% \frac{1}{10}$	(114)
$\displaystyle u_{0}(x,y)$	$\displaystyle=1$	(115)
$\displaystyle v_{0}(x,y)$	$\displaystyle=1$	(116)
$\displaystyle p_{0}(x,y)$	$\displaystyle=1$	(117)

[CHLM24] choose $x_{0}=y_{0}=-0.31$ on a domain $[-1,1]^{2}$ . The exact solution is pure advection. Figure 3 shows the numerical error for various grid sizes at time $t=2$ . One observes third order convergence (due to the high-order advection operator), and for the point values even a machine error for $p$ and $\mathbf{v}$ . This is due to the usage of primitive variables, since both pressure and velocity are constant.

4.1.2 Moving vortex

A standard test case of an isentropic ( $p=\rho^{\gamma}$ ) smooth moving vortex ([HS99], Example 3.3) is studied next. The initial data are

$\displaystyle\left(\begin{array}[]{c}u_{0}(x,y)\\ v_{0}(x,y)\\ \end{array}\right)$	$\displaystyle=\left(\begin{array}[]{c}1\\ 1\\ \end{array}\right)+\frac{\Gamma}{2\pi}\exp\left(\frac{1-r^{2}}{2}\right)\left(% \begin{array}[]{c}-(y-5)\\ x-5\\ \end{array}\right)$	(124)
$\displaystyle T_{0}(x,y):=\frac{p_{0}(x,y)}{\rho_{0}(x,y)}$	$\displaystyle=1-\frac{(\gamma-1)\Gamma^{2}}{8\gamma\pi^{2}}\exp(1-r^{2})$	(125)
$\displaystyle r$	$\displaystyle:=\sqrt{(x-5)^{2}+(y-5)^{2}}$	(126)

with the exact solution being advection with speed $(1,1)$ . However, this time, the pressure is not constant. $\Gamma=5$ is used, as well as a periodic domain of $[0,10]^{2}$ and a final time of $t=10$ , after which the vortex is back where it started from. The errors on different grids are shown in Figure 4. One observes a convergence that is dominated by third order, which starts to deviate slightly towards finer grids.

4.2 Spherical shock tube

Next, the performance of the method is assessed on flows with discontinuities. The first one is the spherical version of the Sod shock tube, i.e.

\displaystyle\rho_{0}(x,y)

\displaystyle=\begin{cases}1&r<0.3\\ 0.125&\text{else}\end{cases}

\displaystyle p_{0}(x,y)

\displaystyle=\begin{cases}1&r<0.3\\ 0.1&\text{else}\end{cases}

\displaystyle\mathbf{v}_{0}(x,y)=0

(127)

Without bound preserving limiting, this test fails; with it, the limiting actually is activated only in the first two time steps (along the initial discontinuity). The results are shown in Fig. 5. One observes some moderate oscillations and some scatter at the contact discontinuity. Part of it might be due to the imprint generated by the limiting at initial time.

4.3 Multi-dimensional Riemann problems

Next, some multi-dimensional Riemann problems from [LL98] are used to assess the ability of the method to cope with more complex interactions between different wave families. The Riemann problems are set up such that outside the interaction region they result in a shock or a simple wave only. In fact, this is well-known to be a difficult situation for numerical methods. As is explained in Section 15.8.4 of [LeV02], whatever the numerical flux used to evolve a discontinuity between two values joined by, say, a Rankine-Hugoniot curve, the very first time step endows the cells in the vicinity with ( $\Delta t$ -dependent) values that are no longer on this curve in general. Thus, the second time step creates waves of all families that propagate away as spurious artifacts. These waves are clearly visible in many simulation results of problems from [LL98] in the literature, and they of course extend into the interaction region. Another reason for their presence is that the initial values given in [LL98] fulfill the jump conditions only to a precision of about $10^{-3}$ .

A selection of multi-d Riemann problems is presented here (Figure 6). The number of grid cells is only half (in each direction) of what [LL98] use. One observes that the features of the interaction region are resolved sharply. Figure 7 shows cuts through the simulation of Configuration 12 for two particular locations, in order to facilitate comparison with [ABK25].

4.4 Low Mach number vortex

Low Mach number flows are of particular interest since numerical methods for compressible flow often require unreasonable grid refinement to resolve subsonic flow. A simple test case is that of a stationary, divergence-free vortex whose Mach number $\epsilon$ can be changed by modifying the background pressure, see [BEK⁺17] for the details of the setup. Results of a first-order Finite Volume method are available e.g. in [Bar21b]. The simulation is evolved until $t=1$ , which is when the quickest part of the flow has completed one revolution, while the sound wave crossing time is about $1/\epsilon$ . Results for various choices of $\epsilon$ on a grid of $50\times 50$ cells are shown in Figure 8. The reader is referred to e.g. [GV99, Del10, Bar21b] for the background on the low Mach number problem. Non low-Mach-compliant methods dissipate the vortex on the sound wave crossing time, while low-Mach-compliant methods do so on an $\epsilon$ -independent time scale in the limit $\epsilon\to 0$ . This is observed in the experimental results: The diffusion (originating from the advection) is asymptotically independent of $\epsilon$ . It is to be noted that this result is by no means trivial and that there is a large literature on “fixing” Finite Volume methods to achieve this. That Active Flux might be low-Mach-compliant has been observed in [BHKR19] as a consequence of the stationarity-preserving property for linear acoustics, and low Mach compliance has been also experimentally observed in [ABK25] for the semi-discrete Active Flux method.

4.5 Kelvin-Helmholtz instability

The Kelvin-Helmholtz instability allows to qualitatively assess the numerical diffusion of a method in subsonic conditions, but in a significantly more complex flow than for an isolated, stationary vortex. Methods that are not low-Mach-compliant are usually artificially stabilizing the setup, preventing the formation of any vortices. The relevant aspect of the next tests is to assess the mere fact of vortex formation. If the initial vorticity is concentrated on a line, then the evolution of the vortices cannot be reasonably studied in the inviscid setting anyway. However, the results can still be compared across different methods in order to qualitatively assess the amount of numerical diffusion. This is the purpose of the test in Section 4.5.2 where the instability is triggered by the passage of a sound wave. The test in Section 4.5.1 has smooth initial vorticity, and therefore until a certain time a well-defined time evolution to which numerical methods converge upon mesh refinement ([LAB⁺24]).

4.5.1 Smooth shear flow

The initial data are as follows ([LAB⁺24]):

	$\displaystyle\rho_{0}(x,y)$	$\displaystyle:=\gamma+R(1-2\eta(y))$	$\displaystyle p_{0}(x,y)$	$\displaystyle:=1$		(128)
	$\displaystyle v_{y,0}(x,y)$	$\displaystyle:=\delta M\sin(2\pi x)$	$\displaystyle v_{x,0}$	$\displaystyle:=M(1-2\eta(y))$		(129)

with

\displaystyle\eta(y)

\displaystyle:=\begin{cases}\frac{1}{2}\left(1+\sin\left(16\pi\left(y+\frac{1}% {4}\right)\right)\right)&-\frac{9}{32}\leq y<-\frac{7}{32}\\ 1&-\frac{7}{32}\leq y<\frac{7}{32}\\ \frac{1}{2}\left(1-\sin\left(16\pi\left(y-\frac{1}{4}\right)\right)\right)&% \frac{7}{32}\leq y<\frac{9}{32}\\ 0&\text{else}\end{cases}

(130)

Here, $R=10^{-3}$ and $\delta=0.1$ . $M$ is a parameter governing the Mach number of the flow, chosen here as $M=0.01$ . The setup is considered on a grid covering $[0,2]\times[-\frac{1}{2},\frac{1}{2}]$ subject to periodic boundaries until a final time of $0.8/M$ . Figure 9 shows the results on several different grids. One observes that the new method is able to resolve well the development of the vortices (see [LAB⁺24] for a comparison with other methods). For low resolution, some spurious features on the interface away from the vortices are visible, which are common among many methods and sometimes even manifest themselves as secondary, spurious vortices. This is not the case here and the features disappear upon mesh refinement.

4.5.2 Instability triggered by the passage of a sound wave

Another type of a Kelvin-Helmholtz instability is a test case proposed in [MRKG03], as a way to assess the ability of the numerical method to cope with acoustic waves in a nearly incompressible flow. A saw-tooth profile of density is hit by an acoustic wave, which induces a shear flow. Since the initial discontinuity in the density translates into a discontinuity in tangential speed, this setup displays ever more features upon grid refinement. Thus, the actual evolution of the vortices cannot be considered as relevant, however, the time scales at which they appear and their minimal size give qualitative information on the numerical diffusion of the method in the subsonic regime. The initial data are given and explained in detail in [BHKR19]. The reader is also referred to this work for comparison of the results with those of the semi-discrete Active Flux, as well as to the original work [MRKG03] (for comparison with Finite Volume methods) and [PP11] (for comparison with a DG method).

Figure 10 shows the results for a rather coarse grid of $200\times 40$ cells. One observes that the vortices of the instability appear much earlier for the present method than for the semi-discrete method from [ABK25], and more homogeneously across the entire interface. This effect appears even stronger if one takes into account the finer resolved grids used in [ABK25]. It thus seems that the present method has significantly less artificial diffusion than its semi-discrete counterpart. Figure 12 shows a simulation on $400\times 80$ cells.

Figure 11 shows the results of a simulation that uses the simple, second-order evolution operator (77) for advection, instead of the higher-order version (99) used otherwise. There are visible differences (for the vortex at $x\simeq 12$ , for example), yet the overall behaviour is similar.

5 Conclusions and outlook

Attempts to develop multi-dimensional Active Flux methods based on evolution operators very quickly reached success for the linear acoustic equations. For the significantly more complicated case of the multi-dimensional Euler equations, several suggestions exist. The main decisions to take in their design, and hence the main differences between existing approaches are whether to use splitting into advection and acoustics (yes), whether to use evolution operators that are exact for the corresponding linearized problems (yes) and whether to use operators that all are sufficiently high-order accurate in time (no) – the terms in brackets indicating the decisions taken in the present work. In practice, the latter question needs to be answered by taking into account the (probably) significant computational effort associated to a fully third-order evolution operator, and the ability to answer it in the affirmative depends on whether it is at all possible to find this operator. Unfortunately, it does not seem to be currently the case, even though a roadmap was traced in [RMF18]. While these considerations apply to the multi-d Euler equations, and also the acoustic suboperator, they do not apply to the advective sub-operator, which has a much simpler nature. As is shown in the present work it can easily be solved to the required order of accuracy. Concerning its accuracy in space, the method is of third order.

The usage of operators that are exact for linearized problems guarantees von Neumann stability of the method. Moreover, since the low Mach compliance for the Euler equations is linked to stationarity preservation for acoustics, a link between the good results of the present method at low Mach number and the corresponding analysis of stationarity preservation in [BHKR19] can be made immediately. The extensive details of an efficient implementation of the exact evolution operator for linear acoustics presented here for the first time hopefully can serve as a reference in future.

Besides the aspects already mentioned, the current work features an example of point values and the reconstruction in primitive variables, as well as a blending of an a priori and a posteriori bound preservation. The usage of primitive variables simplifies the splitting and the corresponding evolution operators greatly, and in certain situations can result in enhanced accuracy.

Numerical results demonstrate that the new method is at least as able to resolve multi-dimensional Riemann problems as the recently developed semi-discrete Active Flux method ([ABK25]), while it might have less diffusion at least in the subsonic regime. It has a higher maximal CFL.

Future work will be devoted to further comparison between the semi-discrete Active Flux and the one based on evolution operators. Other ways to achieve higher order and the possible associated benefits will be investigated.

References

[AB23] Remi Abgrall and Wasilij Barsukow. Extensions of Active Flux to arbitrary order of accuracy. ESAIM: Mathematical Modelling and Numerical Analysis, 57(2):991–1027, 2023.
[Abg22] Rémi Abgrall. A combination of Residual Distribution and the Active Flux formulations or a new class of schemes that can combine several writings of the same hyperbolic problem: application to the 1d Euler equations. Communications on Applied Mathematics and Computation, pages 1–33, 2022.
[ABK25] Rémi Abgrall, Wasilij Barsukow, and Christian Klingenberg. A semi-discrete Active Flux method for the Euler equations on Cartesian grids. Journal of Scientific Computing, 102(2):36, 2025.
[AL24] Rémi Abgrall and Yongle Liu. A new approach for designing well-balanced schemes for the shallow water equations: a combination of conservative and primitive formulations. SIAM Journal on Scientific Computing, 46(6):A3375–A3400, 2024.
[Bar21a] Wasilij Barsukow. The active flux scheme for nonlinear problems. Journal of Scientific Computing, 86(1):1–34, 2021.
[Bar21b] Wasilij Barsukow. Truly multi-dimensional all-speed schemes for the euler equations on cartesian grids. Journal of Computational Physics, 435:110216, 2021.
[BB24] Wasilij Barsukow and Raul Borsche. Implicit active flux methods for linear advection. Journal of Scientific Computing, 98(3):52, 2024.
[BEK⁺17] Wasilij Barsukow, Philipp VF Edelmann, Christian Klingenberg, Fabian Miczek, and Friedrich K Röpke. A numerical scheme for the compressible low-Mach number regime of ideal fluid dynamics. Journal of Scientific Computing, 72(2):623–646, 2017.
[BHKR19] Wasilij Barsukow, Jonathan Hohm, Christian Klingenberg, and Philip L Roe. The active flux scheme on Cartesian grids and its low Mach number limit. Journal of Scientific Computing, 81(1):594–622, 2019.
[BK22] Wasilij Barsukow and Christian Klingenberg. Exact solution and a truly multidimensional Godunov scheme for the acoustic equations. ESAIM: M2AN, 56(1), 2022.
[BKKL24] Wasilij Barsukow, Janina Kern, Christian Klingenberg, and Lisa Lechner. Analysis of the multi-dimensional semi-discrete Active Flux method using the Fourier transform. arXiv preprint arXiv:2412.03477, 2024.
[CHK21] Erik Chudzik, Christiane Helzel, and David Kerkmann. The Cartesian Grid Active Flux method: Linear stability and bound preserving limiting. Applied Mathematics and Computation, 393:125501, 2021.
[CHLM24] Erik Chudzik, Christiane Helzel, and Mária Lukáčová-Medvid’ová. Active Flux methods for hyperbolic systems using the method of bicharacteristics. Journal of Scientific Computing, 99(1):16, 2024.
[DBK25] Junming Duan, Wasilij Barsukow, and Christian Klingenberg. Active Flux methods for hyperbolic conservation laws—Flux Vector splitting and bound-preservation. SIAM Journal on Scientific Computing, 47(2):A811–A837, 2025.
[Del10] Stéphane Dellacherie. Analysis of Godunov type schemes applied to the compressible Euler system at low Mach number. Journal of Computational Physics, 229(4):978–1016, 2010.
[ER11a] Timothy A Eymann and Philip L Roe. Active flux schemes. In 49th AIAA Aerospace Science meeting, 2011.
[ER11b] Timothy A Eymann and Philip L Roe. Active flux schemes for systems. In 20th AIAA computational fluid dynamics conference, 2011.
[ER13] Timothy A Eymann and Philip L Roe. Multidimensional active flux schemes. In 21st AIAA computational fluid dynamics conference, 2013.
[Fan17] Duoming Fan. On the acoustic component of active flux schemes for nonlinear hyperbolic conservation laws. PhD thesis, University of Michigan, Dissertation, 2017.
[FR15] Doreen Fan and Philip L Roe. Investigations of a new scheme for wave propagation. In 22nd AIAA Computational Fluid Dynamics Conference, page 2449, 2015.
[GP16] Jean-Luc Guermond and Bojan Popov. Invariant domains and first-order continuous finite element approximation for hyperbolic systems. SIAM Journal on Numerical Analysis, 54(4):2466–2489, 2016.
[GR63] IS Gradshteyn and IM Ryzhik. Tablitsy integralov, summ, ryadov i proyavleniy.—izd. 4-ye pererab. M.: Fizmatgiz, 1963.
[GV99] Hervé Guillard and Cécile Viozat. On the behaviour of upwind schemes in the low Mach number limit. Computers & fluids, 28(1):63–86, 1999.
[Haj21] Hennes Hajduk. Monolithic convex limiting in discontinuous galerkin discretizations of hyperbolic conservation laws. Computers & Mathematics with Applications, 87:120–138, 2021.
[HKS19] Christiane Helzel, David Kerkmann, and Leonardo Scandurra. A new ADER method inspired by the active flux method. Journal of Scientific Computing, 80(3):1463–1497, 2019.
[HS99] Changqing Hu and Chi-Wang Shu. Weighted essentially non-oscillatory schemes on triangular meshes. Journal of Computational Physics, 150(1):97–127, 1999.
[JS96] Guang-Shan Jiang and Chi-Wang Shu. Efficient implementation of weighted ENO schemes. Journal of computational physics, 126(1):202–228, 1996.
[KLT12] Dmitri Kuzmin, Rainald Löhner, and Stefan Turek, editors. Flux-Corrected Transport: Principles, Algorithms, and Applications. Scientific Computation. Springer Netherlands, Dordrecht, 2012.
[Kuz20] Dmitri Kuzmin. Monolithic convex limiting for continuous finite element discretizations of hyperbolic conservation laws. Computer Methods in Applied Mechanics and Engineering, 361:112804, April 2020.
[LAB⁺24] Giovanni Leidi, R Andrassy, Wasilij Barsukow, J Higl, Philipp VF Edelmann, and Friedrich K Röpke. Performance of high-order godunov-type methods in simulations of astrophysical low mach number flows. Astronomy & Astrophysics, 686:A34, 2024.
[LeV02] Randall J LeVeque. Finite volume methods for hyperbolic problems, volume 31. Cambridge University Press, 2002.
[LL98] Peter D Lax and Xu-Dong Liu. Solution of two-dimensional riemann problems of gas dynamics by positive schemes. SIAM Journal on Scientific Computing, 19(2):319–340, 1998.
[Mae17] Jungyeoul Maeng. On the advective component of active flux schemes for nonlinear hyperbolic conservation laws. PhD thesis, University of Michigan, Dissertation, 2017.
[MRKG03] C-D Munz, Sabine Roller, Rupert Klein, and Karl J Geratz. The extension of incompressible flow solvers to the weakly compressible regime. Computers & Fluids, 32(2):173–196, 2003.
[PP11] Jaime Peraire and Per-Olof Persson. High-order discontinuous Galerkin methods for CFD. In Adaptive high-order methods in computational fluid dynamics, pages 119–152. World Scientific, 2011.
[PS96] Benoit Perthame and Chi-Wang Shu. On positivity preserving finite volume schemes for euler equations. Numerische Mathematik, 73:119–130, 1996.
[RMF18] Philip L Roe, Jungyeoul Maeng, and Doreen Fan. Comparing active flux and Discontinuous Galerkin methods for compressible flow. In 2018 AIAA Aerospace Sciences Meeting, page 0836, 2018.
[Roe17] Philip Roe. Multidimensional upwinding. Handbook of Numerical Analysis, 18:53–80, 2017.
[Roe21] Philip Roe. Designing CFD methods for bandwidth—a physical approach. Computers & Fluids, 214:104774, 2021.
[Vil19] François Vilar. A posteriori correction of high-order discontinuous galerkin scheme through subcell finite volume formulation and flux reconstruction. Journal of Computational Physics, 387:245–279, 2019.
[vL77] Bram van Leer. Towards the ultimate conservative difference scheme. IV. A new approach to numerical convection. Journal of computational physics, 23(3):276–299, 1977.

Appendix A Spherical means

The acoustic evolution operator involves spherical means of the initial data, on which certain differential operators act. Apart from the usual spherical mean of a scalar function, for vector-valued functions there exist additional spherical means of interest that involve the unit normal $\mathbf{n}=\left(\begin{array}[]{c}\cos\varphi\sin\vartheta\\ \sin\varphi\sin\vartheta\\ \cos\vartheta\\ \end{array}\right)$ . For integration over the angular coordinates in three spatial dimensions, $\mathrm{d}\Omega:=\mathrm{d}\varphi\mathrm{d}\vartheta\sin\vartheta$ is used occasionally as an abbreviation.

A.1 General definitions

Definition A.1.

Consider a scalar-valued function $f\colon\mathbb{R}^{3}\to\mathbb{R}$ and a vector valued function $\mathbf{v}\colon\mathbb{R}^{3}\to\mathbb{R}^{3}$ , $\mathbf{v}=(v_{1},\ldots,v_{n})$ .

The spherical mean $M[f]\colon\mathbb{R}^{3}\times\mathbb{R}^{+}_{0}\to\mathbb{R}$ of $f$ is

\displaystyle M[f](\mathbf{x},r)

\displaystyle:=\frac{1}{4\pi}\int_{S^{2}}\mathrm{d}\mathbf{y}\,f(\mathbf{x}+r% \mathbf{y})

(131)

ii.

The spherical mean $M[\mathbf{v}]\colon\mathbb{R}^{3}\times\mathbb{R}^{+}_{0}\to\mathbb{R}^{3}$ of a vector-valued function $\mathbf{v}$ , is defined componentwise:

\displaystyle M[\mathbf{v}](\mathbf{x},r):=\left(M[v_{1}](\mathbf{x},r),\ldots% ,M[v_{n}](\mathbf{x},r)\right)

(132)

In components, this can be expressed as $M[\mathbf{v}]_{i}=M[v_{i}]$ .

iii.

The spherical mean $M[\mathbf{n}f]\colon\mathbb{R}^{3}\times\mathbb{R}^{+}_{0}\to\mathbb{R}^{3}$ of $f$ is

\displaystyle M[\mathbf{n}f](\mathbf{x},r)

\displaystyle:=\frac{1}{4\pi}\int_{S^{2}}\mathrm{d}\mathbf{y}\,\mathbf{y}f(% \mathbf{x}+r\mathbf{y})

(133)

iv.

The spherical mean $M[\mathbf{n}\cdot\mathbf{v}]\colon\mathbb{R}^{3}\times\mathbb{R}^{+}_{0}\to% \mathbb{R}$ of $\mathbf{v}$ is

\displaystyle M[\mathbf{n}\cdot\mathbf{v}](\mathbf{x},r):=\frac{1}{4\pi}\int_{% S^{2}}\mathrm{d}\mathbf{y}\,\mathbf{y}\cdot\mathbf{v}(\mathbf{x}+r\mathbf{y})

(134)

The spherical mean $M[\mathbf{n}(\mathbf{n}\cdot\mathbf{v})]\colon\mathbb{R}^{3}\times\mathbb{R}^{% +}_{0}\to\mathbb{R}^{3}$ of $\mathbf{v}$ is

\displaystyle M[\mathbf{n}(\mathbf{n}\cdot\mathbf{v})](\mathbf{x},r):=\frac{1}% {4\pi}\int_{S^{2}}\mathrm{d}\mathbf{y}\,\mathbf{y}(\mathbf{y}\cdot\mathbf{v}(% \mathbf{x}+r\mathbf{y}))

(135)

Observe that the normalization is such that $M[1]=1$ . For actual computations, the spherical mean is best expressed in spherical coordiantes, e.g.

	$\displaystyle M[f](\mathbf{x},r)$	$\displaystyle=\frac{1}{4\pi}\int_{0}^{2\pi}\mathrm{d}\varphi\int_{0}^{\pi}% \mathrm{d}\vartheta\,\sin\vartheta\,f(x+r\cos\varphi\sin\vartheta,y+r\sin% \varphi\sin\vartheta,z+r\cos\vartheta)$		(136)
	$\displaystyle M[\mathbf{n}f](\mathbf{x},r)$	$\displaystyle=\frac{1}{4\pi}\int_{0}^{2\pi}\mathrm{d}\varphi\int_{0}^{\pi}% \mathrm{d}\vartheta\,\left(\begin{array}[]{c}\cos\varphi\sin\vartheta\\ \sin\varphi\sin\vartheta\\ \cos\vartheta\\ \end{array}\right)\sin\vartheta\,f(x+r\cos\varphi\sin\vartheta,y+r\sin\varphi% \sin\vartheta,z+r\cos\vartheta)$		(140)

etc. The definition for a spherical mean over a smaller angular domain is analogous.

A.2 Spherical means in two spatial dimensions

On two-dimensional domains, integration over subsets of the full domain of $\varphi$ are of interest (because the data might be piecewise defined).

Definition A.2.

The angular domain $W=[\varphi_{\text{min}},\varphi_{\text{max}})\times[0,\pi]$ is called a wedge $W\subset S^{2}$ . Spherical means whose $(\varphi,\vartheta)$ integration domain is $W$ are denoted by $M^{W}$ :

\displaystyle M^{W}[f](\mathbf{x},r)

\displaystyle:=\frac{1}{2(\varphi_{\text{max}}-\varphi_{\text{min}})}\int_{% \varphi_{\text{min}}}^{\varphi_{\text{max}}}\mathrm{d}\varphi\int_{0}^{\pi}% \mathrm{d}\vartheta\,\sin\vartheta\,f(x+r\cos\varphi\sin\vartheta,y+r\sin% \varphi\sin\vartheta)

(141)

with $f:\mathbb{R}^{2}\to\mathbb{R}$ . An analogous definition holds for other types of spherical means such as $M^{W}[\mathbf{n}f]$ .

Remark A.1.

Observe that it only makes sense to define $M^{W}[f]$ for $f$ that do not depend on $z$ .

Remark A.2.

In two spatial dimensions it is natural to only consider the first two components of $\mathbf{n}$ . This will not be made explicit in the notation. This is consistent, since for any $f$ that does not depend on $z$ ,

	$\displaystyle M^{W}[\mathbf{n}f]_{z}$	$\displaystyle=\frac{1}{2(\varphi_{\text{max}}-\varphi_{\text{min}})}\int_{% \varphi_{\text{min}}}^{\varphi_{\text{max}}}\mathrm{d}\varphi\int_{0}^{\pi}% \mathrm{d}\vartheta\,\sin\vartheta\cos\vartheta\,f(x+r\cos\varphi\sin\vartheta% ,y+r\sin\varphi\sin\vartheta)$			(142)
		$\displaystyle\overset{c:=\cos\vartheta}{=}\frac{1}{2(\varphi_{\text{max}}-% \varphi_{\text{min}})}\int_{\varphi_{\text{min}}}^{\varphi_{\text{max}}}% \mathrm{d}\varphi\int_{-1}^{1}\mathrm{d}c\,c\,f(x+r\cos\varphi\sqrt{1-c^{2}},y% +r\sin\varphi\sqrt{1-c^{2}})\$	$\displaystyle=0$		(143)

because an antisymmetric function is integrated over a symmetric interval.

The normalization $M^{W}[1]=1$ is easily confirmed:

\displaystyle\int_{W\subset S^{2}}\mathrm{d}\mathbf{y}=\int_{\varphi_{\text{% min}}}^{\varphi_{\text{max}}}\mathrm{d}\varphi\int_{0}^{\pi}\mathrm{d}% \vartheta\,\sin\vartheta=(\varphi_{\text{max}}-\varphi_{\text{min}})\int_{-1}^% {1}\mathrm{d}(\cos\vartheta)=2(\varphi_{\text{max}}-\varphi_{\text{min}})

(144)

A.3 Spherical means for polynomial data in two space dimensions

In principle, it is possible to evaluate everything analytically. However, in particular for angular domains that are not the entire sphere, the expressions become complicated. Here, some formulae are given that allow to (pre)compute the expressions exactly for polynomial initial data in 2 spatial dimensions.

For polynomial initial data, the key object to compute therefore is the spherical mean

\displaystyle M^{W}[n_{x}^{\alpha}n_{y}^{\beta}x^{a}y^{b}]

(145)

A.3.1 Integration over $\vartheta$

In the following we compute the spherical mean of a polynomial function in two spatial dimensions. For polynomial functions, the dependence on $\vartheta$ will be of the form $\sin^{a}\vartheta$ , which justifies the following definition.

Definition A.3.

Denote by $\displaystyle\eta_{a}:=\int_{0}^{\pi}\mathrm{d}\vartheta\sin^{a}\vartheta$ .

Theorem A.1 ([GR63], 3.621 (3. and 4.)).

One has

\displaystyle\eta_{a}=\sqrt{\pi}\frac{\Gamma\left(\frac{a+1}{2}\right)}{\Gamma% \left(\frac{a+2}{2}\right)}=\frac{(a-1)!!}{a!!}\cdot\begin{cases}\pi&\text{$a>% 0$ even}\\ 2&\text{$a>0$ odd}\end{cases}

(146)

where the double factorial is defined as

\displaystyle A!!:=\begin{cases}A(A-2)(A-4)\cdots 4\cdot 2&\text{$A$ even}\\ A(A-2)(A-4)\cdots 3\cdot 1&\text{$A$ odd}\end{cases}

(147)

This means that $(\eta_{a})_{a=0,1,\ldots}=\pi,2,\frac{\pi}{2},\frac{4}{3},\frac{3\pi}{8},\frac% {16}{15},\frac{5\pi}{16},\frac{32}{35},\ldots$ . Below, we will need the special case

\displaystyle\eta_{i+j+1}=\frac{(i+j)!!}{(i+j+1)!!}\cdot\begin{cases}\pi&\text% {$i+j+1$ even}\\ 2&\text{else}\end{cases}

(148)

A.3.2 Integration over $\varphi$

Definition A.4.

Denote by $\displaystyle\mu^{W}_{ij}:=\int_{\varphi_{\text{min}}}^{\varphi_{\text{max}}}% \mathrm{d}\varphi\cos^{i}\varphi\sin^{j}\varphi$ .

Theorem A.2.

We have

	$\displaystyle\mu^{W}_{ij}$	$\displaystyle=\frac{1}{2^{i+j}}\sum_{m=0}^{i}\sum_{n=0}^{j}{i\choose m}{j% \choose n}(-1)^{n+j}$		(149)
		$\displaystyle\phantom{mmmm}\times\begin{cases}(-1)^{j/2}(\varphi_{\text{max}}-% \varphi_{\text{min}})&i+j=2(m+n),j\text{ even}\\ &\\ \displaystyle\left.\frac{(-1)^{j/2}\sin([i+j-2m-2n]\varphi)}{i+j-2m-2n}\right\|% _{\varphi_{\text{min}}}^{\varphi_{\text{max}}}&i+j\neq 2(m+n),j\text{ even}\\ \displaystyle\left.\frac{(-1)^{\frac{j-1}{2}}\cos([i+j-2m-2n]\varphi)}{i+j-2m-% 2n}\right\|_{\varphi_{\text{min}}}^{\varphi_{\text{max}}}&i+j\neq 2(m+n),j\text% { odd}\\ 0&\text{else}\end{cases}$		(150)

with the notation $f(\varphi)\Big{|}_{\varphi_{\text{min}}}^{\varphi_{\text{max}}}=f(\varphi_{% \text{max}})-f(\varphi_{\text{min}})$ .

Proof.

The integral can be evaluated to

$\displaystyle\mu_{ij}^{W}=\int_{\varphi_{\text{min}}}^{\varphi_{\text{max}}}% \mathrm{d}\varphi\cos^{i}\varphi\sin^{j}\varphi$	$\displaystyle=\int_{\varphi_{\text{min}}}^{\varphi_{\text{max}}}\mathrm{d}% \varphi\left(\frac{\mathrm{e}^{\mathbb{i}\varphi}-\mathrm{e}^{-\mathbb{i}% \varphi}}{2\mathbb{i}}\right)^{j}\left(\frac{\mathrm{e}^{\mathbb{i}\varphi}+% \mathrm{e}^{-\mathbb{i}\varphi}}{2}\right)^{i}$	(151)
	$\displaystyle=\frac{1}{2^{i+j}\mathbb{i}^{j}}\int_{\varphi_{\text{min}}}^{% \varphi_{\text{max}}}\mathrm{d}\varphi\,\,\mathrm{e}^{\mathbb{i}(i+j)\varphi}(% 1-\mathrm{e}^{-2\mathbb{i}\varphi})^{j}(1+\mathrm{e}^{-2\mathbb{i}\varphi})^{i}$	(152)
	$\displaystyle=\frac{1}{2^{i+j}\mathbb{i}^{j}}\sum_{m=0}^{i}\sum_{n=0}^{j}{i% \choose m}{j\choose n}(-1)^{n}\int_{\varphi_{\text{min}}}^{\varphi_{\text{max}% }}\mathrm{d}\varphi\,\,\mathrm{e}^{\mathbb{i}[i+j-2m-2n]\varphi}$	(153)

Now,

\displaystyle\int_{\varphi_{\text{min}}}^{\varphi_{\text{max}}}\mathrm{d}% \varphi\,\,\mathrm{e}^{\mathbb{i}[i+j-2m-2n]\varphi}=\begin{cases}\varphi_{% \text{max}}-\varphi_{\text{min}}&i+j-2m-2n=0\\ &\\ \left.\displaystyle\frac{\mathrm{e}^{\mathbb{i}[i+j-2m-2n]\varphi}}{\mathbb{i}% [i+j-2m-2n]}\right|_{\varphi_{\text{min}}}^{\varphi_{\text{max}}}&\text{else}% \end{cases}

(154)

and thus

\displaystyle\mu_{ij}^{W}

\displaystyle=\frac{1}{2^{i+j}(-1)^{j}}\sum_{m=0}^{i}\sum_{n=0}^{j}{i\choose m% }{j\choose n}(-1)^{n}\begin{cases}\mathbb{i}^{j}(\varphi_{\text{max}}-\varphi_% {\text{min}})&i+j-2m-2n=0\\ &\\ \left.\displaystyle-\frac{\mathbb{i}^{j+1}\mathrm{e}^{\mathbb{i}[i+j-2m-2n]% \varphi}}{i+j-2m-2n}\right|_{\varphi_{\text{min}}}^{\varphi_{\text{max}}}&% \text{else}\end{cases}

(155)

As the integral is real, we only need to include real terms into the sum:

	$\displaystyle\mu_{ij}^{W}$	$\displaystyle=\frac{1}{2^{i+j}}\sum_{m=0}^{i}\sum_{n=0}^{j}{i\choose m}{j% \choose n}(-1)^{n+j}$		(156)
		$\displaystyle\phantom{mmmm}\times\begin{cases}(-1)^{j/2}(\varphi_{\text{max}}-% \varphi_{\text{min}})&i+j=2(m+n),j\text{ even}\\ &\\ \displaystyle\left.\frac{(-1)^{j/2}\sin([i+j-2m-2n]\varphi)}{i+j-2m-2n}\right\|% _{\varphi_{\text{min}}}^{\varphi_{\text{max}}}&i+j\neq 2(m+n),j\text{ even}\\ \displaystyle\left.\frac{(-1)^{\frac{j-1}{2}}\cos([i+j-2m-2n]\varphi)}{i+j-2m-% 2n}\right\|_{\varphi_{\text{min}}}^{\varphi_{\text{max}}}&i+j\neq 2(m+n),j\text% { odd}\\ 0&\text{else}\end{cases}$		(157)

∎

The first few values of $\mu_{ij}^{W}$ are:

$i$	$j$	$\mu_{ij}^{W}\cdot\Delta\varphi$
0	0	$\Delta\varphi$
0	1	$-\Delta\cos(\varphi)$
0	2	$\frac{1}{4}\left(-\Delta\sin(2\varphi)+2\Delta\varphi\right)$
1	0	$\Delta\sin(\varphi)$
1	1	$-\frac{1}{4}\Delta\cos(2\varphi)$
1	2	$\frac{1}{4}\left(-\frac{\Delta\sin(3\varphi)}{3}+\Delta\sin(\varphi)\right)$
2	0	$\frac{1}{4}\left(\Delta\sin(2\varphi)+2\Delta\varphi\right)$
2	1	$-\frac{1}{4}\left(\frac{\Delta\cos(3\varphi)}{3}+\Delta\cos(\varphi)\right)$
2	2	$\frac{1}{16}\left(-\frac{\Delta\sin(4\varphi)}{2}+2\Delta\varphi\right)$

Here, $\Delta\varphi=\varphi_{\text{max}}-\varphi_{\text{min}}$ and $\Delta\sin(n\varphi)=\sin(n\varphi_{\text{max}})-\sin(n\varphi_{\text{min}})$ etc.

A.3.3 Assembling the spherical mean

By linearity and by the assumption that we are concerned with the two-dimensional case only it is sufficient to compute the spherical mean at $\mathbf{x}=(x,y,0)$ of

\displaystyle n_{x}^{\alpha}n_{y}^{\beta}x^{a}y^{b}

(158)

Therefore it is necessary to evaluate the following spherical mean

$\displaystyle M^{W}[n_{x}^{\alpha}n_{y}^{\beta}x^{a}y^{b}](\mathbf{x},r)$	$\displaystyle=\frac{1}{2(\varphi_{\text{max}}-\varphi_{\text{min}})}$	(159)
	$\displaystyle\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!% \!\!\!\!\!\cdot\int_{\varphi_{\text{min}}}^{\varphi_{\text{max}}}\mathrm{d}% \varphi\int_{0}^{\pi}\mathrm{d}\vartheta\,\,\sin\vartheta(\cos\varphi\sin% \vartheta)^{\alpha}(\sin\varphi\sin\vartheta)^{\beta}(x+r\cos\varphi\sin% \vartheta)^{a}(y+r\sin\varphi\sin\vartheta)^{b}$
	$\displaystyle\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!% \!\!\!\!\!=\frac{1}{2(\varphi_{\text{max}}-\varphi_{\text{min}})}\int_{\varphi% _{\text{min}}}^{\varphi_{\text{max}}}\mathrm{d}\varphi\int_{0}^{\pi}\mathrm{d}% \vartheta\,\,\sin^{1+\alpha+\beta}\vartheta\cos^{\alpha}\varphi\sin^{\beta}% \varphi(x+r\cos\varphi\sin\vartheta)^{a}(y+r\sin\varphi\sin\vartheta)^{b}$	(160)

Rewrite now

\displaystyle(x+r\cos\varphi\sin\vartheta)^{a}(y+r\sin\varphi\sin\vartheta)^{b% }=\sum_{k=0}^{a}\sum_{\ell=0}^{b}{a\choose k}{b\choose\ell}x^{k}y^{\ell}r^{a-k% +b-\ell}\sin^{b-\ell}\varphi\cos^{a-k}\varphi\sin^{a-k+b-\ell}\vartheta

(161)

This proves the following

Theorem A.3.

With $\Delta\varphi=\varphi_{\text{max}}-\varphi_{\text{min}}$ , the spherical mean of monomial data reads

\displaystyle M^{W}[n_{x}^{\alpha}n_{y}^{\beta}x^{a}y^{b}](\mathbf{x},r)

\displaystyle=\sum_{k=0}^{a}\sum_{\ell=0}^{b}{a\choose k}{b\choose\ell}x^{k}y^% {\ell}r^{a-k+b-\ell}\frac{\mu^{W}_{a+\alpha-k,b+\beta-\ell}}{2\Delta\varphi}% \eta_{a+\alpha-k+b+\beta-\ell+1}

(162)

We rewrite it as

\displaystyle M^{W}[n_{x}^{\alpha}n_{y}^{\beta}x^{a}y^{b}](\mathbf{x},r)

\displaystyle=\sum_{k=0}^{a}\sum_{\ell=0}^{b}x^{k}y^{\ell}r^{a-k+b-\ell}

(163)

with

\displaystyle F^{(a,b,\alpha,\beta)}_{k,\ell}:={a\choose k}{b\choose\ell}\frac% {\mu^{W}_{a+\alpha-k,b+\beta-\ell}}{2\Delta\varphi}\eta_{a+\alpha-k+b+\beta-% \ell+1}

(164)

In particular, if $x=y=0$ , then the formula simplifies to

\displaystyle M^{W}[n_{x}^{\alpha}n_{y}^{\beta}x^{a}y^{b}](0,r)

\displaystyle=r^{a+b}\cdot\frac{\mu^{W}_{a+\alpha,b+\beta}}{2\Delta\varphi}% \eta_{a+\alpha+b+\beta+1}

(165)

Corollary A.1.

If $\mathbf{v}=(x^{a}y^{b},0,0)$ , then the spherical mean $M^{W}[\mathbf{n}\cdot\mathbf{v}]$ is

	$\displaystyle M^{W}[\mathbf{n}\cdot\mathbf{v}]$	$\displaystyle=M[\mathbf{n}x^{a}y^{b}]_{x}$		(166)
		$\displaystyle=\frac{1}{2\Delta\varphi}\sum_{k=0}^{a}\sum_{\ell=0}^{b}{a\choose k% }{b\choose\ell}x^{k}y^{\ell}r^{a-k+b-\ell}\mu_{a-k+1,b-\ell}(\varphi_{\text{% max}},\varphi_{min})\eta_{a-k+b-\ell+2}$		(167)

Analogously, if $\mathbf{v}=(0,x^{a}y^{b},0)$ , then

	$\displaystyle M^{W}[\mathbf{n}\cdot\mathbf{v}]$	$\displaystyle=M[\mathbf{n}x^{a}y^{b}]_{y}$		(168)
		$\displaystyle=\frac{1}{2\Delta\varphi}\sum_{k=0}^{a}\sum_{\ell=0}^{b}{a\choose k% }{b\choose\ell}x^{k}y^{\ell}r^{a-k+b-\ell}\mu_{a-k,b-\ell+1}(\varphi_{\text{% max}},\varphi_{min})\eta_{a-k+b-\ell+2}$		(169)

Corollary A.2.

If $\mathbf{v}=(x^{a}y^{b},0,0)$ , then the spherical mean $M[\mathbf{n}(\mathbf{n}\cdot\mathbf{v})]$ , restricted to two components, is

\displaystyle M^{W}[\mathbf{n}(\mathbf{n}\cdot\mathbf{v})]=\frac{1}{2\Delta% \varphi}\sum_{k=0}^{a}\sum_{\ell=0}^{b}{a\choose k}{b\choose\ell}x^{k}y^{\ell}% r^{a-k+b-\ell}\left(\begin{array}[]{c}\mu_{a-k+2,b-\ell}(\varphi_{\text{max}},% \varphi_{min})\\ \mu_{a-k+1,b-\ell+1}(\varphi_{\text{max}},\varphi_{min})\\ \end{array}\right)\eta_{a-k+b-\ell+3}

(172)

If $\mathbf{v}=(0,x^{a}y^{b},0)$ , then the spherical mean $M[\mathbf{n}\cdot\mathbf{v}]$ , restricted to two components, is

\displaystyle M^{W}[\mathbf{n}(\mathbf{n}\cdot\mathbf{v})]=\frac{1}{2\Delta% \varphi}\sum_{k=0}^{a}\sum_{\ell=0}^{b}{a\choose k}{b\choose\ell}x^{k}y^{\ell}% r^{a-k+b-\ell}\left(\begin{array}[]{c}\mu_{a-k+1,b-\ell+1}(\varphi_{\text{max}% },\varphi_{min})\\ \mu_{a-k,b-\ell+2}(\varphi_{\text{max}},\varphi_{min})\\ \end{array}\right)\eta_{a-k+b-\ell+3}

(175)

Appendix B Pseudocode implementation of the evolution operator

⬇

1FunctionPolynomial[] timeEvolution(FunctionPolynomial[] initial, double[] position, double /*@$\phi_\text{min}$@*/, double /*@$\phi_\text{max}$@*/){

2/*@ This time evolution yields the result as a polynomial in $r$ that needs to be evaluated at $ct$@*/

4 FunctionPolynomial[] res = new FunctionPolynomial[initial.length];

5 double fraction = Math.abs(/*@$\phi_\text{max}$@*/ - /*@$\phi_\text{min}$@*/) / 2 / Math.PI;

7 FunctionPolynomial R = new FunctionPolynomial(1, 0);

8 R.setCoefficient(1, 0, 1.0); // polynomial

r

10 /*@$M^W[p_0]$@*/ = sphericalMean(/*@$\phi_\text{min}$@*/, /*@$\phi_\text{max}$@*/, initial[PRESS], position);

12 /*@$M^W[n_x p_0]$@*/ = sphericalMean(/*@$\phi_\text{min}$@*/, /*@$\phi_\text{max}$@*/, initial[PRESS], position, 1, 0);

13 /*@$M^W[n_y p_0]$@*/ = sphericalMean(/*@$\phi_\text{min}$@*/, /*@$\phi_\text{max}$@*/, initial[PRESS], position, 0, 1);

15 /*@$M^W[u_0]$@*/ = sphericalMean(/*@$\phi_\text{min}$@*/, /*@$\phi_\text{max}$@*/, initial[XVEL], position);

16 /*@$M^W[v_0]$@*/ = sphericalMean(/*@$\phi_\text{min}$@*/, /*@$\phi_\text{max}$@*/, initial[YVEL], position);

18 /*@$M^W[\vec n \cdot \vec v]$@*/ = sphericalMean(/*@$\phi_\text{min}$@*/, /*@$\phi_\text{max}$@*/, initial[XVEL], position, 1,0).add(

19 sphericalMean(/*@$\phi_\text{min}$@*/, /*@$\phi_\text{max}$@*/, initial[YVEL], position, 0,1) );

21 /*@$M^W[n_x (\vec n \cdot \vec v)]$@*/ = sphericalMean(/*@$\phi_\text{min}$@*/, /*@$\phi_\text{max}$@*/, initial[XVEL], position, 2,0).add(

22 sphericalMean(/*@$\phi_\text{min}$@*/, /*@$\phi_\text{max}$@*/, initial[YVEL], position, 1,1) );

23 /*@$M^W[n_y (\vec n \cdot \vec v)]$@*/ = sphericalMean(/*@$\phi_\text{min}$@*/, /*@$\phi_\text{max}$@*/, initial[XVEL], position, 1,1).add(

24 sphericalMean(/*@$\phi_\text{min}$@*/, /*@$\phi_\text{max}$@*/, initial[YVEL], position, 0,2) );

27 PPterm = /*@$M^W[p_0]$@*/.mult(R).diffX();

28 PUVterm = (/*@$M^W[\vec n \cdot \vec v]$@*/.mult(R).mult(R).diffX()).divideByX();

29 UPterm = ((/*@$M^W[n_x p_0]$@*/.mult(R).mult(R)).diffX()).divideByX();

30 UUVterm = (((/*@$M^W[n_x (\vec n \cdot \vec v)]$@*/.mult(R).mult(R).mult(R).diffX()

31 ).divideByX().sub(/*@$M^W[u_0]$@*/.mult(R))).diffX()).divideByX().integrateX();

32 VPterm = ((/*@$M^W[n_y p_0]$@*/.mult(R).mult(R)).diffX()).divideByX();

33 VUVterm = (((/*@$M^W[n_y (\vec n \cdot \vec v)]$@*/.mult(R).mult(R).mult(R).diffX()

34 ).divideByX().sub(/*@$M^W[v_0]$@*/.mult(R))).diffX()).divideByX().integrateX();

36 res[PRESS] = PPterm.sub(PUVterm).mult(fraction);

37 res[XVEL] = FunctionPolynomial.constantPolynomial(initial[XVEL].eval(position[0], position[1])

38 ).sub(UPterm).add(UUVterm).mult(fraction);

39 res[YVEL] = FunctionPolynomial.constantPolynomial(initial[YVEL].eval(position[0], position[1])

40 ).sub(VPterm).add(VUVterm).mult(fraction);

42 return res;

43}

46FunctionPolynomial sphericalMean(double /*@$\phi_\text{min}$@*/, double /*@$\phi_\text{max}$@*/, FunctionPolynomial /*@$f$@*/, double[] position, int /*@$\alpha$@*/, int /*@$\beta$@*/){

47 /*@ Computes $M^W[n_x^\alpha y^\beta f]$ for polynomial $f$ at @*/position;

48 /*@ Radius $r$ not specified, result of computation is polynomial in $r$ @*/

49 FunctionPolynomial res = new FunctionPolynomial(/*@$f$@*/.degreeX() + /*@$f$@*/.degreeY()); // degree of result

51 for (int i = 0; i <= /*@$f$@*/.degreeX(); i++){

52 for (int j = 0; j <= /*@$f$@*/.degreeY(); j++){

53 res = res.add(sphericalMeanMonomial(/*@$\phi_\text{min}$@*/, /*@$\phi_\text{max}$@*/, i, j, position, /*@$\alpha$@*/, /*@$\beta$@*/).mult(/*@$f$@*/.getCoefficient(i, j)));

54 }

55 }

56 return res;

57}

60FunctionPolynomial sphericalMeanMonomial(double /*@$\phi_\text{min}$@*/, double /*@$\phi_\text{max}$@*/, int /*@a@*/, int /*@b@*/, double[] position, int /*@$\alpha$@*/, int /*@$\beta$@*/){

61 /*@ Computes $M^W[n_x^\alpha y^\beta x^a y^b]$ at @*/position;

62 /*@ Radius $r$ not specified, result of computation is polynomial in $r$ @*/

64 FunctionPolynomial sphericalMean = new FunctionPolynomial(/*@a@*/ + /*@b@*/); // degree of result

66 int phiIntervalIndex = /*@identify the quadrant required (can be passed directly as parameter)@*/

68 for (int /*@$k$@*/ = 0; /*@$k$@*/ <= /*@$a$@*/; /*@$k$@*/++){

69 for (int /*@$\ell$@*/ = 0; /*@$\ell$@*/ <= /*@$b$@*/; /*@$\ell$@*/++){

70 sphericalMean.addToCoefficient(/*@$a+b-\ell-k$@*/, 0,

71 factorInsideEvolution[/*@$a$@*/][/*@$b$@*/][/*@$k$@*/][/*@$\ell$@*/][/*@$\alpha$@*/][/*@$\beta$@*/][phiIntervalIndex]*

72 power(position[0], /*@$k$@*/) *

73 power(position[1], /*@$\ell$@*/));

74 }

75 }

76 return sphericalMean;

77}

80precomputeFactors(double /*@$\phi_\text{min}$@*/, double /*@$\phi_\text{max}$@*/, int phiIntervalIndex){

81 for (int /*@$a$@*/ = 0; /*@$a$@*/ <= /*@$a_\text{max}$@*/; /*@$a$@*/++){

82 for (int /*@$b$@*/ = 0; /*@$b$@*/ <= /*@$b_\text{max}$@*/; /*@$b$@*/++){

83 for (int /*@$k$@*/ = 0; /*@$k$@*/ <= /*@$a$@*/; /*@$k$@*/++){

84 for (int /*@$\ell$@*/ = 0; /*@$\ell$@*/ <= /*@$b$@*/; /*@$\ell$@*/++){

85 for (int /*@$\alpha$@*/ = 0; /*@$\alpha$@*/ <= /*@$\alpha_\text{max}$@*/; /*@$\alpha$@*/++){

86 for (int /*@$\beta$@*/ = 0; /*@$\beta$@*/ <= /*@$\beta_\text{max}$@*/; /*@$\beta$@*/++){

87 factorInsideEvolution[/*@$a$@*/][/*@$b$@*/][/*@$k$@*/][/*@$\ell$@*/][/*@$\alpha$@*/][/*@$\beta$@*/][phiIntervalIndex] =

88 /*@${a \choose k} {b \choose \ell}$@*/ * sphericalMeanIntegrationFactor(/*@$\phi_\text{min}$@*/, /*@$\phi_\text{max}$@*/, /*@$a-k+\alpha$@*/, /*@$b-\ell+\beta$@*/);

89 /*@ This is $F^{(a, b,\alpha,\beta)}_{k,\ell}$. @*/

90 }

91 }

92 }

93 }

94 }

95 }

96}

98double sphericalMeanIntegrationFactor(double /*@$\phi_\text{min}$@*/, double /*@$\phi_\text{max}$@*/, int i, int j){

99 /*@ This function computes $\eta_{i+j+1} \mu^W_{ij} / 2 / \Delta \phi$.@*/

100 double secondFactor = 0.0;

101 double factor = 0.0;

102

103 double firstFactor = /*@$(i+j)!!/(i+j+1)!!$@*/;

104 if ((i+j+1) /*@$\mod 2$@*/ == 0){

105 firstFactor *= /*@$\pi$@*/;

106 } else {

107 firstFactor *= 2.0;

108 }

109 /*@ This is $\eta_{i+j+1}$.@*/

110

111

112 for (int /*@$m$@*/ = 0; /*@$m$@*/ <= i; /*@$m$@*/++){

113 for(int /*@$n$@*/ = 0; /*@$n$@*/ <= j; /*@$n$@*/++){

114 factor = /*@${i \choose m} {j \choose n}$@*/;

115 if (/*@$i+j-2m-2n$@*/ == 0){

116 if(/*@$j\mod 2$@*/==0){

117 factor *= /*@$(-1)^{j/2+n} (\phi_\text{max}-\phi_\text{min})$@*/;

118 } else {

119 factor = 0;

120 }

121 } else {

122 factor *= /*@$(-1)^{n+j} / (i+j-2m-2n)$@*/;

123 if(/*@$j\mod 2$@*/==0){

124 factor *= /*@$(-1)^{j/2}(\sin(\phi_\text{max}(i+j-2m-2n))-\sin(\phi_\text{min}(i+j-2m-2n)))$@*/;

125 } else {

126 factor *= /*@$(-1)^{(j-1)/2} (\cos(\phi_\text{max}(i+j-2m-2n))-\cos(\phi_\text{min}(i+j-2m-2n)))$@*/;

127 }

128 }

129 secondFactor += factor;

130 }

131 }

132 secondFactor /= /*@$2^{i+j}$@*/;

133

134 return firstFactor*secondFactor / 2 / (/*@$\phi_\text{max}$@*/ - /*@$\phi_\text{min}$@*/);

135}

Listing 1: Pseudocode implementation of the evolution operator.

The expressions are simplified tremendously through an implementation of polynomial functions; such an implementation is given in the Appendix C as a reference, where the definitions of functions used in the above listings can be found. It is possible to precompute not just some coefficients, but everything that is known a priori. This is best achieved algorithmically by using the previous implementation to compute necessary coefficients:

⬇

1precomputeEvolutionCoefficients(){

2 double[][] allPositions = /*@where one will want to evaluate the evolution operator, e.g. the 8 points at the boundary of cell@*/;

3/*@ @*/

4 /*@ Positions are assumed to be in $\left[-\frac{\Delta x}{2}, \frac{\Delta x}{2}\right] \times \left[-\frac{\Delta y}{2}, \frac{\Delta y}{2}\right]$@*/

6 evolutionCoefficients = new FunctionPolynomial[3][3][3][allRelativePositions.size()][4][3];

7 /*@ The order is: variable of solution, $a$, $b$, position, quadrant, variable of initial data @*/

8 /*@ A quadrant is the angle between@*/ /*@$\pi$@*//2*angleIndex /*@and@*/ /*@$\pi$@*//2*(angleIndex+1);

10 FunctionPolynomial[] initial;

11 int positionIndex;

13 for (int q = 0; q <= 2; q++){ // variables u,v,p

14 for (int /*@ a @*/ = 0; /*@ a @*/ <= 2; /*@ a @*/++){

15 for (int /*@ b @*/ = 0; /*@ b @*/ <= 2; /*@ b @*/++){

16 initial = new FunctionPolynomial[3];

17 for (int q2 = 0; q2 <= 2; q2++){ initial[q2] = FunctionPolynomial.ZERO(); }

18 initial[q] = new FunctionPolynomial(2, 2);

19 initial[q].setCoefficient(/*@a@*/, /*@b@*/, 1.0);

20 /*@ The q-th variable is now set to $x^a y^b$@*/

22 positionIndex = 0;

23 for (double[] position : allPositions) {

24 for (int angleIndex = 0; angleIndex <= 3; angleIndex++){

25 for (int q2 = 0; q2 <= 2; q2++){ // variables u,v,p

26 evolutionCoefficients[q][/*@a@*/][/*@b@*/][positionIndex][angleIndex][q2] = timeEvolution(initial, position, PI/2*angleIndex, PI/2*(angleIndex+1))[q2];

27 }

28 }

29 positionIndex++;

30 }

31 }

32 }

33 }

34}

36double[] timeEvolutionUsingPrecomputed(FunctionPolynomial[] initial, double[] position, double /*@$t$@*/, double /*@$\phi_\text{min}$@*/, double /*@$\phi_\text{max}$@*/){

37 /*@ positions is assumed to be in $\left[-\frac{\Delta x}{2}, \frac{\Delta x}{2}\right] \times \left[-\frac{\Delta y}{2}, \frac{\Delta y}{2}\right]$@*/

38 double[] res = new double[initial.length];

39 int positionIndex = /*@identify the index corresponding to@*/ position /*@(can be passed directly as parameter)@*/;

40 int[] angleIndices = /*@identify the quadrants required (one or several; can be passed directly as parameter)@*/

42 for (int q = 0; q <= 2; q++){ // variables u,v,p

43 for (int /*@$a$@*/ = 0; /*@$a$@*/ <= 2; /*@$a$@*/++){

44 for (int /*@$b$@*/ = 0; /*@$b$@*/ <= 2; /*@$b$@*/++){

45 for (int angleIndex = 0; angleIndex < angleIndices.length; angleIndex++){

46 for (int q2 = 0; q2 <= 2; q2++){ // variables u,v,p

47 res[q2] += evolutionCoefficients[q][/*@$a$@*/][/*@$b$@*/][positionIndex][angleIndices[angleIndex]][q2].eval(/*@$ct$@*/) * initial[q].getCoefficient(/*@$a$@*/, /*@$b$@*/);

48 }

49 }

50 }

51 }

52 }

53 return res;

54}

Listing 2: Pseudocode implementation of the optimized evolution operator.

⬇

1if (edge.isVertical()){

2 double[] leftValue = timeEvolutionUsingPrecomputed(edge.leftCell().reconstruction(), new double[]{0.5*dx, 0.0}, /*@$\Delta t$@*/, /*@$\pi$@*//2, 3.0*/*@$\pi$@*//2);

3 double[] rightValue = timeEvolutionUsingPrecomputed(edge.rightCell().reconstruction(), new double[]{-0.5*dx, 0.0}, /*@$\Delta t$@*/, -/*@$\pi$@*//2, /*@$\pi$@*//2);

4 return leftValue + rightValue;

5} else if (edge.isHorizontal(){

6 double[] bottomValue = timeEvolutionUsingPrecomputed(edge.bottomCell().reconstruction(), new double[]{0.0, 0.5*dy}, /*@$\Delta t$@*/, -/*@$\pi$@*/, 0);

7 double[] topValue = timeEvolutionUsingPrecomputed(edge.topCell().reconstruction(), new double[]{0.0, -0.5*dy}, /*@$\Delta t$@*/, 0, /*@$\pi$@*/);

8 return bottomValue + topValue;

Listing 3: Example of the usage of the evolution operator.

Appendix C Implementation details

⬇

1public class FunctionPolynomial extends Function{

2 protected double[][] coefficients;

4 public FunctionPolynomial(int degx, int degy) {

5 coefficients = new double[degx+1][degy+1];

6 }

8 protected void resize(int iMax, int jMax){

9 double[][] coefficientsNew = new double[Math.max(iMax + 1, coefficients.length)][Math.max(jMax + 1, coefficients[0].length)];

10 for (int i = 0; i < coefficients.length; i++){

11 for (int j = 0; j < coefficients[0].length; j++){

12 coefficientsNew[i][j] = coefficients[i][j];

13 }

14 }

15 coefficients = coefficientsNew;

16 }

18 public void setCoefficient(int i, int j, double value){

19 if ((i >= coefficients.length) || (j >= coefficients[0].length)){ resize(i, j); }

20 coefficients[i][j] = value;

21 }

23 public void addToCoefficient(int i, int j, double value){

24 if ((i >= coefficients.length) || (j >= coefficients[0].length)){ resize(i, j); }

25 coefficients[i][j] += value;

26 }

28 public FunctionPolynomial integrateX(){

29 FunctionPolynomial res = new FunctionPolynomial(degreeX()+1, 0);

30 res.setCoefficient(0, 0, 0.0);

31 for (int j = 0; j < coefficients[0].length; j++){

32 for (int i = 0; i < coefficients.length; i++){

33 res.setCoefficient(i+1, j, coefficients[i][j]/(i+1));

34 }

35 }

36 return res;

37 }

39 public FunctionPolynomial diffX(){

40 FunctionPolynomial res = new FunctionPolynomial(degreeX(), degreeY());

41 for (int j = 0; j < coefficients[0].length; j++){

42 for (int i = 1; i < coefficients.length; i++){

43 res.setCoefficient(i-1, j, coefficients[i][j]*i);

44 }

45 res.setCoefficient(coefficients.length-1, j, 0);

46 }

47 return res;

48 }

50 public FunctionPolynomial divideByX(){

51 // ATTENTION: the x^0 term is lost!

52 FunctionPolynomial res = new FunctionPolynomial(degreeX(), degreeY());

53 for (int j = 0; j < coefficients[0].length; j++){

54 for (int i = 1; i < coefficients.length; i++){

55 res.setCoefficient(i-1, j, coefficients[i][j]);

56 }

57 res.setCoefficient(coefficients.length-1, j, 0);

58 }

59 return res;

60 }

62 public double eval(double x, double y){

63 double val = 0;

65 double powX, powY;

67 powX = 1.0;

68 for (int i = 0; i < coefficients.length; i++){

69 powY = 1.0;

70 for (int j = 0; j < coefficients[0].length; j++){

71 val += coefficients[i][j] * powX * powY;

72 powY *= y;

73 }

74 powX *= x;

75 }

77 return val;

78 }

80 public Function mult(double val){

81 FunctionPolynomial res = new FunctionPolynomial(degreeX(), degreeY());

82 for (int i = 0; i < coefficients.length; i++){

83 for (int j = 0; j < coefficients[0].length; j++){

84 res.setCoefficient(i, j, val*coefficients[i][j]);

85 }

86 }

87 return res;

88 }

90 public FunctionPolynomial mult(FunctionPolynomial p){

91 FunctionPolynomial res = new FunctionPolynomial(p.degreeX()+degreeX(), p.degreeY()+degreeY());

92 for (int i = 0; i < coefficients.length; i++){

93 for (int j = 0; j < coefficients[0].length; j++){

94 for (int ip = 0; ip <= p.degreeX(); ip++){

95 for (int jp = 0; jp <= p.degreeY(); jp++){

96 res.setCoefficient(ip + i, jp + j, res.getCoefficient(ip+i, jp+j) + getCoefficient(i,j)*p.getCoefficient(ip,jp));

97 }

98 }

99 }

100 }

101 return res;

102 }

103

104 public double getCoefficient(int xPower, int yPower){

105 if ((xPower >= coefficients.length) || (yPower >= coefficients[0].length)){

106 return 0.0;

107 } else {

108 return coefficients[xPower][yPower];

109 }

110 }

111

112 public int[] degree(){

113 return new int[] {degreeX(), degreeY()};

114 }

115

116 public int degreeX(){

117 return coefficients.length-1;

118 }

119

120 public int degreeY(){

121 return coefficients[0].length-1;

122 }

123

124 public FunctionPolynomial add(FunctionPolynomial poly){

125 int[] otherDegree = poly.degree();

126 int[] myDegree = this.degree();

127 FunctionPolynomial res = new FunctionPolynomial(Math.max(otherDegree[0], myDegree[0]), Math.max(otherDegree[1], myDegree[1]));

128 int[] deg = res.degree();

129

130 for (int j = 0; j <= deg[1]; j++){

131 for (int i = 0; i <= deg[0]; i++){

132 res.setCoefficient(i, j, getCoefficient(i, j) + poly.getCoefficient(i,j));

133 }

134 }

135 return res;

136 }

137

138 public FunctionPolynomial sub(FunctionPolynomial poly){

139 int[] otherDegree = poly.degree();

140 int[] myDegree = this.degree();

141 FunctionPolynomial res = new FunctionPolynomial(Math.max(otherDegree[0], myDegree[0]), Math.max(otherDegree[1], myDegree[1]));

142 int[] deg = res.degree();

143

144 for (int j = 0; j <= deg[1]; j++){

145 for (int i = 0; i <= deg[0]; i++){

146 res.setCoefficient(i, j, getCoefficient(i, j) - poly.getCoefficient(i,j));

147 }

148 }

149 return res;

150 }

151

152 public static FunctionPolynomial ZERO(){

153 return constantPolynomial(0);

154 }

155

156 public static FunctionPolynomial constantPolynomial(double val){

157 FunctionPolynomial res = new FunctionPolynomial(0, 0);

158 res.setCoefficient(0, 0, val);

159 return res;

160 }

161

162

163}

Listing 4: An example of an implementation of polynomial functions.

Abstract

1 Introduction

2 Evolution operator for linear acoustics

2.1 Initial data in p𝑝pitalic_p

2.2 Initial data in u𝑢uitalic_u

3 An Active Flux method for the Euler equations

3.1 General remarks

3.2 Update of the average

3.3 Update of the point values

Remark 3.1.

Remark 3.2.

Remark 3.3.

3.4 Order of accuracy

3.4.1 Overview

3.4.2 A high-order evolution operator for nonlinear advection

Theorem 3.1.

Proof.

Theorem 3.2.

Proof.

3.5 Bound preservation

3.5.1 Point value update

3.5.2 Average update

4 Numerical results

4.1 Convergence tests

4.1.1 Contact wave

4.1.2 Moving vortex

4.2 Spherical shock tube

4.3 Multi-dimensional Riemann problems

4.4 Low Mach number vortex

4.5 Kelvin-Helmholtz instability

4.5.1 Smooth shear flow

4.5.2 Instability triggered by the passage of a sound wave

5 Conclusions and outlook

References

Appendix A Spherical means

A.1 General definitions

Definition A.1.

A.2 Spherical means in two spatial dimensions

Definition A.2.

Remark A.1.

Remark A.2.

A.3 Spherical means for polynomial data in two space dimensions

A.3.1 Integration over ϑitalic-ϑ\varthetaitalic_ϑ

Definition A.3.

Theorem A.1 ([GR63], 3.621 (3. and 4.)).

A.3.2 Integration over φ𝜑\varphiitalic_φ

Definition A.4.

Theorem A.2.

Proof.

A.3.3 Assembling the spherical mean

Theorem A.3.

Corollary A.1.

Corollary A.2.

Appendix B Pseudocode implementation of the evolution operator

Appendix C Implementation details

2.1 Initial data in $p$

2.2 Initial data in $u$

A.3.1 Integration over $\vartheta$

A.3.2 Integration over $\varphi$