Stationarity Preservation Properties of the Active Flux Scheme on Cartesian Grids

Hyperbolic systems of conservation laws in multiple spatial dimensions display features absent in the one-dimensional case, such as involutions and non-trivial stationary states. These features need to be captured by numerical methods without excessive grid refinement. The active flux method is an extension of the finite volume scheme with additional point values distributed along the cell boundary. For the equations of linear acoustics, an exact evolution operator can be used for the update of these point values. It incorporates all multi-dimensional information. The active flux method is stationarity preserving, i.e


Introduction
Hyperbolic conservation laws in several space dimensions show a number of phenomena which are absent in one-dimensional situations. E.g., in case of the Euler equations, these additional phenomena include vortices, nontrivial stationary states (themselves governed by PDEs), involutions, and the incompressible (low Mach number) limit. Design principles for numerical methods in one spatial dimension have been subject of successful investigation over several decades (correct approximation of weak solutions, stability, entropy, ⋯ ). Numerical methods able to capture truly multi-dimensional phenomena on coarse grids still lack a comprehensive set of design principles.

3
There exist examples of numerical methods that can cope with particular multi-dimensional phenomena without requiring excessive refinement of the computational grid. Numerical methods able to reproduce the low Mach number limit are referred to as low Mach number compliant, numerical methods that keep stationary a discretization of an involution-involution-preserving, numerical methods with stationary states being discretizations of all the stationary states of the PDE-stationarity preserving.
A standard approach for constructing numerical methods in multiple spatial dimensions is to use a one-dimensional method in a direction-by-direction fashion, i.e., solving a sequence of one-dimensional problems in different directions. Such methods generally are unable to resolve multi-dimensional features on coarse grids and only achieve this by means of a "fix" (e.g., [2,3,6]). Such a fix generally affects other properties, e.g., the stability.
Several ways of incorporating truly multi-dimensional information in a numerical method have been suggested. They differ in the kind of information included and in the abilities of the resulting method. For instance including the solution of multi-dimensional Riemann problems, exactly or approximately, does improve the resolution of multi-dimensional shock interactions, but has been shown in [5] not to yield a stationarity preserving or involution preserving method in general. Retaining the enlarged stencil, but modifying the way how cells contribute "over the corner" has, however, in certain cases proven itself a successful way to obtain involution preserving schemes, e.g., in [2,10,12,13].
An entirely different way to incorporate multi-dimensional information into a numerical method has been suggested in [7] and given the name active flux. It is conceptually based on the so-called "Scheme V" from [14], which is a one-dimensional finite volume scheme with additional point values located at cell boundaries. These additional pointwise degrees of freedom are evolved using an exact evolution operator, in the case of [14] solving the initial value problem for one-dimensional linear advection. The initial data for this update are provided by a conservative reconstruction of the data at the previous time step which interpolates the point values at cell boundaries.
The scheme from [14] has been extended to multiple spatial dimensions on triangular grids and to the equations of linear acoustics in [7]. The exact solution operator for linear acoustics is used for the update of the point values, and thus, the entire (multidimensional) information is retained. In [4], this method has been extended to Cartesian grids and it has been shown that it then yields a stationarity preserving and vorticity preserving scheme. Since then, the method has been applied to other equations ( [1,8,9,11]), in particular to those for which the IVP cannot easily be solved exactly. In these cases, an approximate evolution operator has to be used for the update of the pointwise degrees of freedom. To still guarantee structure preservation properties for an active flux method endowed with an approximate evolution operator, it is important to understand the way active flux achieves structure preservation with an exact evolution operator. This then shall allow to derive conditions on the approximate evolution operators.
This paper contributes to this aim by providing an experimental study of the discrete stationary states of active flux for linear acoustics. The paper is organized as follows. Section 2 introduces the equations of linear acoustics, Sect. 3 describes the active flux scheme in detail, and in Sect. 4, the discrete stationary states of the active flux scheme for linear acoustics are presented. The behavior of setups that correspond to discrete stationary states, and of those that do not, is analyzed in Sect. 5 using carefully chosen numerical test cases.

3 2 The Initial Value Problem for Linear Acoustics
Consider a hyperbolic system of m conservation laws in d spatial dimensions: and the corresponding initial value problem with q(0, ) = q 0 ( ) . One example of a hyperbolic system (2.1) is the Euler equations of ideal hydrodynamics 1 The system is closed with the equation of state of an ideal gas: Upon linearization around the state of constant density ̄= const , constant pressure p = const , and vanishing velocity ̄ = 0 , the following system arises: with c = √p ∕̄ . The last two equations form a closed system. Replacing p∕̄↦ p yields the acoustic equations: They form a hyperbolic system with eigenvalues ±c and 0. They also possess the involution ∇ × , which always remains stationary: Stationary states of (2.9)-(2.10) are given by p = const and ∇ ⋅ = 0 . Thus, both the involution and the stationary states involve differential operators. A numerical involution should involve some discretization of the curl, and numerical stationary states should be characterized by some discretization of the divergence. This mimics the case of more complicated equations, such as the Euler equations, or the equations of magnetohydrodynamics.
In particular, the latter are famous for possessing an involution ( ∇ ⋅ = 0 ) which is hard to capture numerically.
The acoustic system (2.9)-(2.10) can be solved exactly. The exact solution operator is derived in [5], and the result is reproduced below. Definition 1 (Spherical mean) The spherical mean M f ( , r) of an integrable function f that depends on = (x, y, z) ∈ ℝ 3 is given by with the outward normal vector given by Theorem 1 (Solution operator) The solution to the initial value problem for (2.9)-(2.10) is given by The acoustic system thus offers the possibility to study aspects of multi-dimensional hyperbolic systems while having access to an exact evolution operator. It can also be used in the active flux method.

The Active Flux Scheme on Cartesian Grids
The following focuses on the two-dimensional case d = 2 . Consider a Cartesian grid covering the computational domain with cells

Degrees of Freedom
Upon integration of (2.1) over a time step [t n , t n+1 ] and over a computational cell C ij , the Gauss theorem allows to write where the latter is the flux average through an edge e ⊂ C ij .
Active flux additionally evolves independent point values at the cell boundaries. In this paper, they are chosen to be distributed at the following locations: The pointwise degrees of freedom are correspondingly denoted by q N ij , q EH ij , q EV ij (see Fig. 1). The pointwise degrees of freedom are shared between the cells. Having four nodes and four edges, every cell thus has access to 9 degrees of freedom (including the cell average). The continuity across cell boundaries leads to a reconstruction procedure that differs from the usual finite volume approach.

Reconstruction
To update the point values in time, already in [14], it has been suggested to solve an initial value problem at their locations, with a reconstruction of the data at time t n serving as the initial datum. However, [14] focused on the one-dimensional linear advection only.
In multiple spatial dimensions, the 9 degrees of freedom accessible to every cell make a biparabolic reconstruction inside every cell natural. The reconstruction interpolates the pointwise degrees of freedom and its average matches the average in the cell. The explicit formula can be found in, e.g., [4].
Along every edge of the cell, this reconstruction reduces to a parabola. Every edge involves three pointwise degrees of freedom which are interpolated, and they suffice to uniquely define this parabola. Therefore, the reconstruction is continuous across cell boundaries, but, in general, has discontinuous derivatives.

Update of Point Values
For certain, in particular linear equations, an exact evolution operator is available. It then is used to find the exact solution of the initial value problem at the location of the pointwise degree of freedom. Nonlinear problems generally require approximate evolution operators (see, e.g., [1,9]). In this paper, the exact evolution operator of linear acoustics, as discussed in Sect. 2, is used.

Update of the Averages
After the update of the point values, they are available at times t n , t n+ 1 2 , and t n+1 . These values serve as quadrature points in an approximation of the flux according to (3.3): Simpson's rule) and having written = (f x , f y ). This allows to update the averages according to (3.1). In the context of the active flux method, the usual problem of finding a numerical flux encountered for finite volume schemes is thus replaced by the problem of evolving the additional pointwise degrees of freedom in time.

Summary
Given both the averages and the point values at cell interfaces at time t n , the algorithm of the active flux thus consists of the following steps.
(i) Reconstruct the data in every cell, such that the reconstruction is conservative and interpolates the pointwise degrees of freedom. (ii) Use an evolution operator for the update of the pointwise degrees of freedom to times t n+ 1 2 and t n+1 (for acoustics, (2.13)-(2.14) is used). (iii) Compute intercell flux via quadrature (equation (3.6)).
(iv) Update the average as in a standard finite volume scheme (equation (3.1)).

Properties of the Active Flux Method
Certain properties of the active flux are entirely different from usual finite volume schemes. They are briefly addressed in the following.
-The usage of independently evolved point values makes the computation of the flux through the cell boundary very easy (hence, the name active flux). -Point values shared by adjacent cells imply a continuous reconstruction and discontinuities are approximated by steep gradients. For nonlinear problems, the continuous data might self-steepen into a shock-here, care must be taken in the design of the (approximate) evolution operator (see [1] for examples). However, this is not an issue for the linear problems discussed here. -The active flux as described above is a time-explicit method. It is stable, because the evolution of the point values respects the characteristic directions and thus naturally provides upwinding. Additional upwinding in the computation of the intercell flux is not necessary. -The active flux scheme as described above is linear and has a compact stencil. This means that the update of the degrees of freedom of each cell (point values and the cell average) is a linear function of the degrees of freedom in a neighborhood of the cell. The corresponding coefficients (which follow from the reconstruction and the exact evolution operator) need to be computed only once. The update thus can be implemented very efficiently. -The active flux scheme is third-order accurate in both space and time (see also [4]), because the reconstruction is parabolic, and a half-step in time is used.
Nonlinear problems make approximate evolution operators necessary. They introduce additional error into the scheme, and the properties of the scheme will generally depend on the approximation chosen for the evolution operator. This paper focuses on structure preservation properties. A first step of the analysis should address the question whether structure preserving properties are attainable at all, i.e., even upon the usage of the exact evolution operator. If so, then this result requires any approximate evolution operator not to spoil these properties. Therefore, equations which admit exact evolution operators are valuable to separate the approximation introduced by the active flux in general from the possible influence of an approximate evolution operator. This paper focuses on the equations of linear acoustics.

Stationarity Preservation
Numerical methods achieve stability by adding the diffusion which manifests itself in a decay of Fourier modes. The von Neumann stability criterion, therefore, enforces every Fourier mode to be nonincreasing in time. Fourier modes characterizing numerical stationary states (of which, e.g., the uniformly constant state is a simple example) neither decay nor increase in time. In multiple spatial dimensions, stationary states can be very diverse (e.g., all divergence-free vector fields for linear acoustics). The analysis of stationary Fourier modes has been recognized in [2] to be an efficient way of finding all numerical stationary states of a numerical method. Whenever structures (e.g., stationary states) are themselves given as differential operators, the question arises whether a numerical method can be at all capable of "preserving" them. A numerical method only has access to a finite amount of information and surely cannot achieve this preservation exactly. A comparison of the numerical stationary states to the stationary states of the PDE has two possible outcomes. The numerical stationary states should be a discretization of all the stationary states of the PDE. However, often, they only discretize a subset of all the stationary states of the PDE. For example, certain numerical methods keep only spatially constant states stationary (see [2] for a precise definition of trivial stationary states). They thus entirely fail to capture all the others. The method in this case adds too much diffusion and this diffusion acts even on states which should remain stationary. Numerical schemes with numerical stationary states discretizing the full set of stationary states of the PDE are called stationarity preserving.
A linear PDE with nontrivial stationary states possesses an involution, i.e., a functional of the solution which remains stationary even if the solution does not. Stationarity-preserving schemes mimic this feature: such a numerical scheme always possesses a discrete involution. In the case of linear acoustics (2.9)-(2.10), the involution is the vorticity ∇ × and, thus, stationarity preserving schemes are also vorticity preserving. For more details, see [2].
Thus, stationarity preserving numerical methods keep those discrete data stationary that fulfill a particular discrete counterpart to div = 0 (in the acoustic case, for example). The question whether any other discretization is also kept stationary must generally be answered in the negative. This leads to the question what happens when the discrete data are not fulfilling this particular discretization. Note that this is the general case, as no preprocessing or projection of the initial data is ever assumed. This is discussed in Sect. 5, while the remainder of this section presents the discrete stationary states of the active flux.
A prominent feature of the active flux method in this paper is the usage of the exact solution operator (2.14)-(2.13) in the point value update. Obviously, the exact evolution operator keeps divergence-free velocities and a constant pressure stationary. Therefore, the point values remain stationary if the reconstruction of the pressure is constant, and the reconstruction of the velocity is divergence-free. This implies conditions on the values of the discrete degrees of freedom. In [4], the usage of the exact evolution operator in the active flux method has been recognized as the crucial ingredient for a drastic simplification of the study of the numerical stationary states. Now, only the reconstruction needs to be dealt with, and the evolution operator is merely involved through the fact that its stationary states are known.
It is easiest to study numerical stationary states using the discrete Fourier transform. Define t x ∶= exp( k x Δx) , t y ∶= exp( k y Δy) with the imaginary unit and = (k x , k y ) ∈ ℝ 2 the wave vector characterizing the spatial frequency of the Fourier mode. Then, apply the Fourier transform to the four sets of degrees of freedom: Recall that (2.9)-(2.10) in two spatial dimensions form a linear system of three equations, and therefore, the total set of Fourier modes is

Theorem 2
If Q is such that then p recon = 0 and div recon = 0 uniformly.
Proof See [4]. These equations are discretizations of div = 0 , and (by inverting the Fourier transform) can be written as the following finite difference formulae:

Corollary 1 From
Here, the following notation has been used: The notation is combined for different directions, e.g., After the conditions are clarified under which the point values remain stationary, the remaining question concerns the stationarity of the cell averages.

Theorem 3
If the discrete degrees of freedom on the grid fulfill the conditions of Theorem 2, then the update (3.1) of the cell averages leaves the averages stationary, as well.
The proof can be found in [4]. Thus, the active flux scheme as described above on Cartesian grids is stationarity preserving.

Numerical Study of Stationary States
Having identified the numerical stationary states does not yet provide an answer to the question how initial data behave that are sampled from a stationary state of the PDE, but fail to satisfy the discrete relations (4.9)-(4.12). As pointed out in [2] for finite difference methods, such initial data will undergo a transition towards one of the discrete stationary states. For stationarity preserving schemes, the smaller Δx , the closer do the discrete stationary states approximate the stationary states of the PDE. Therefore, the transition can be expected to be exponentially quick and the discrete stationary state onto which the simulation settles to be an equally recognizable representation of the stationary state of the PDE. An experimental demonstration of this transition for the active flux scheme using carefully chosen setups is subject of Sect. 5.2.
When the initial data satisfy (4.13)-(4.14), the active flux scheme should remain stationary up to the machine precision. This statement can be observed in simulations and Sect. 5.1 serves to illustrate this behavior.

Discrete Stationary State
Consider a Fourier mode Q exp( k x iΔx + k y jΔy) . Theorem 2 states a choice of Q (which via t x , t y depends on ) that remains stationary. This mode is complex valued, but adding up two such modes with opposite signs of yields a real-valued grid function:

3
To see how this mode arises, consider, e.g., û N , ignore Δx and Δy for the moment, and observe that k ↦ −k maps t x ↦ 1∕t x . Then, Here, k x = 2 , k y = 10 ⋅ 2 is set; the pressure is chosen to vanish identically. This initial datum is now evolved in a 100 × 100 grid covering [0, 1] 2 with periodic boundaries. Figure 2 shows the initial datum. Figure 3 shows that it, indeed, remains stationary up to the machine precision, and that the discrete divergences (4.13)-(4.14) remain at the level of machine precision.

General Stationary States
Numerical initial data originating from a stationary state of the PDE generally do not fulfill (4.13)-(4.14), and, therefore, generally do not remain exactly stationary. A stable scheme applies diffusion to all non-stationary Fourier modes (von Neumann stability). Therefore, one can expect the numerical evolution to stationarize on data fulfilling (4.13)-(4.14). In practical applications, numerical data originating from a stationary state of the PDE evolve to a numerical stationary state which is virtually indistinguishable, although, of course, the transition can be noticed if measured carefully. The numerical stationary state is not diffused further and persists forever.
Consider, therefore, the following initial data of a divergence-free vortex: with r = √ (x − 0.5) 2 + (y − 0.5) 2 and w = 0.2 . This setup is solved on a 10 × 10 computational grid covering [0, 1] 2 with zero-gradient boundaries. Figure 4 shows the decay of the discrete divergences (4.13)-(4.14) and the stationarization of the setup. One observes that the discrete divergences decay exponentially to machine level. To demonstrate how special the decay to machine zero of the discrete divergences (4.13)-(4.14) is, consider now an alternative discretization of the divergence: The time evolution of this discrete divergence is shown in Fig. 5. One observes that it becomes stationary, as the entire setup, but it does not decay.   Fig. 4. Measuring another discrete divergence (5.6) instead of (4.13). One observes that the discrete divergence becomes stationary (as is the case for the entire setup), but that it does not reach (machine) zero

Conclusion and Outlook
The active flux method is an extension of the finite volume scheme which incorporates point values distributed along the cell boundary. This paper deals with the active flux method on Cartesian grids for the equations of linear acoustics. They admit an exact evolution operator which is used to update the pointwise degrees of freedom and, thus, incorporates all multi-dimensional information. This is presumably the reason why the active flux method is found to resolve multi-dimensional features much better than usual finite volume schemes. In particular, it is stationarity preserving, i.e., it discretizes all the stationary states of the PDE. Note that, in multiple spatial dimensions, the stationary states of a hyperbolic system of PDEs themselves are given by PDEs. Numerical schemes in general add diffusion, which, in many cases, affects even states which should remain stationary. The fact that a scheme is stationarity preserving can, thus, be reinterpreted as a careful choice of the numerical diffusion. This paper demonstrates experimental evidence for the discrete stationary states of the active flux method. In particular, the numerical evolution of such a state is demonstrated to be only due to the machine error. In general, the initial data that originate from a stationary state of the PDE do not fulfill the particular discretization of a numerical stationary state. However, by stability, the instationary part is diffused away and the setup transitions exponentially quickly to a numerical stationary state which then persists forever. For stationarity preserving schemes, this new stationary state is an equally recognizable discrete representation of the stationary state of the PDE. This is demonstrated experimentally in this paper.
The active flux scheme thus seems to allow to construct truly multi-dimensional structure preserving schemes in a natural way. Nonlinear systems of conservation laws, however, generally do not admit exact solutions in the closed form. Future work shall be devoted to an understanding of the conditions that lead to structure preservation for the active flux scheme endowed with an approximate evolution operator.
Funding Open access funding provided by the University of Zurich.

Compliance with ethical standards
Conflict of interest The authors declare that they have no conflict of interest.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creat iveco mmons .org/licen ses/by/4.0/.