## 1 Introduction

Solutions to the initial value problem for hyperbolic systems of conservation laws

(1) | ||||||

(2) |

generically develop discontinuities in finite time, even if the initial data are smooth. The natural setting therefore are weak solutions. For convergence to the weak solution of (1), a numerical method needs to be conservative (Lax-Wendroff theorem).

A popular way to derive numerical methods (due to Godunov) is to introduce discontinuities at every cell interface (reconstruction step), and to evolve piecewise constant data over a short period of time (Riemann solver). Even if a more complicated reconstruction is used inside the cell, or if an approximate evolution is used instead of the exact one, the main idea of Godunov’s approach remains to introduce discontinuities at every cell interface. This approach has also been adapted to Finite Element methods, which has led to the development of Discontinuous Galerkin (DG) methods. These methods are conservative.

However, convergence to the weak solution is not enough in practice. In view of the big computational effort associated with grid refinement (particularly in multi-d), there is ongoing interest in guaranteeing properties of numerical solutions for coarse grids already. Some of them are essential for the simulation to continue running, e.g. non-negativity of the density. Other shortcomings might not cause a simulation to crash, but still are a big problem in practice, in particular if they require the computational grid to be much finer than reasonable. For example it might turn out to be necessary to resolve a vortex with a number of grid cells that depends on how long one intends to run the simulation for. This can occur in simulations of fluid flow phenomena that are smooth (low Mach number/incompressible limit, vortices, …). They are not well approximated by standard Godunov methods on coarse grids because the effects of numerical diffusion become more and more pronounced as the simulation runs. Grid refinement increases the time scales on which diffusion becomes overwhelming, but – besides slowing down the simulation – does not solve the problem in a fundamental way.

The difficulties of the Godunov method in regions of smooth flow are not surprising, as it introduces shocks everywhere in the computational domain. For example, the failure of the Godunov method to resolve the low Mach number limit on coarse grids has been traced back directly to the acoustic waves introduced at every time step at every cell interface (in [GM04]). For low Mach number flow, modifications of Godunov methods have been suggested which restore the correct behaviour, but they are ad hoc and cannot guarantee adequate treatment of other aspects of smooth multi-dimensional flow (for example vortices). Note at this point that in the quest for a new numerical method we expect it to be able to work well in presence of shocks as well as when they are absent.

The argumentation of [GM04] assumes that the multi-dimensional numerical method is constructed using a one-dimensional Riemann solver. The fluxes through each interface are computed while ignoring the influence of the corners of the cell. Such an approach shall be called directionally-split. Truly multi-dimensional phenomena, where contributions from different direction balance each other, tend to be a challenge for dimensionally split methods. The low Mach number/incompressible limit is one example of such a truly multi-dimensional phenomenon, because the divergence-free condition is trivial in one spatial dimension. One might therefore conjecture that the difficulties do not originate in the Godunov method per se, but in the fact that quasi-one-dimensional solvers are applied in a directionally split way, and thus multi-dimensionality is not taken into account appropriately. One might think that perhaps it would suffice to take into account all the multi-dimensional interactions in multi-dimensional Riemann problems to obtain a Godunov-type method that resolves smooth multi-dimensional problems accurately on coarse grids.

Unfortunately, the situation is not that simple. In [BK22] this line of thought has been developed completely for the linear acoustic equations, which have interesting multi-dimensional features that are not captured by standard dimensionally-split Godunov-type methods on coarse grids (e.g. an analogue of the low Mach number limit, and the involution of vorticity). The complete solution of the four-quadrant Riemann problem for linear acoustics has been obtained and a two-dimensional Godunov method on Cartesian grids constructed in [BK22], whose evolution step was thus exact and fully multi-dimensional. None of the desirable properties (low Mach number compliance, vorticity preservation) was found to hold true for this method. Taking into account all multi-dimensional interactions from multi-dimensional Riemann problems is not enough (and, of course, quite an effort).

Similar observations described in [Roe17] have sparked the development of the Active Flux method by Roe and collaborators ([ER11], based on [vL77]), a numerical method with a continuous reconstruction, whose evolution step thus requires a short-time solution of the IVP for (1) with *continuous* data.

This does not mean that this method is only applicable in the subsonic regime. Continuous data can self-steepen over the course of a time step, and in such cases, the update needs to be performed carefully, before the discontinuity is projected onto a continuous reconstruction again. In a typical simulation, however, shocks are located on sets of codimension 1 (e.g. in 1-d on countably many points), and thus in this sense almost everywhere in the computational domain the assumption of continuity is correct. Compare this to the approach of a Godunov method solving a rarefaction: instead of maintaining the continuous solution, the rarefaction is cut apart in a staircase shape at every time step, and over the course of the evolution step every jump of the staircase is evolving again into a small rarefaction.

Leaving aside conceptual questions, from a practical point of view, the answers to the following three questions matter most: Is this method stable, possibly even under an explicit time integration? Does the numerical solution obtained with this method converge to the weak solution of (1), i.e. is there an analogue of the Lax-Wendroff-Theorem? Does this method have favorable properties on smooth, multi-dimensional solutions? The answer to all three questions is yes, which makes it an interesting alternative to Godunov-type methods.

The degrees of freedom of the traditional Active Flux ([vL77, ER11]) are cell averages and, additionally, point values located at cell interfaces and shared by adjacent cells. This allows for a continuous (but not necessarily differentiable) reconstruction, which passes through the given point values, and whose average agrees with the given one. The evolution of the point value requires the solution of the IVP for (1) over a short time step for continuous initial data. This building block replaces the Riemann solver in Godunov-type methods. The evolution of the average is obtained as usual by integrating the conservation law over the cell and over a time step, and by applying Gauss’ theorem. The fluxes at cell interfaces can be immediately evaluated using the available point values. The evolution of the averages therefore is conservative, and an analogue of a Lax-Wendroff theorem for Active Flux has been shown in [Abg20]. Active Flux decouples the problem of average evolution from the computation of interface values, which are declared independent degrees of freedom, hence the name Active Flux, although it would be more accurate to say “Active Point Value at Cell Interface”. This allows for more flexibility. The update of point values is not subject to the constraint of conservation, because conservation applies only to averages. Therefore it is possible to even switch variables and update the point values in, say, primitive variables. This has been pointed out in [Abg20]. A more thorough description to the traditional Active Flux method is given in Section 2, as well as in [ER11, ER13, BHKR19, HKS19, Bar21].

The development of a high order method requires chiseling out the distinctive features of the low order method. The aim of the present paper therefore is not only to propose practical high order extensions of Active Flux. These extensions are tied to possible interpretations of Active Flux that allow to place it amidst other families of numerical methods. The three extensions that we propose allow to see Active Flux either as a coupled Finite Volume/Finite Difference method (Section 3), or as an enriched Finite Volume method (Section 4), or as a coupled Finite Element/Finite Difference method (Section 5). Which of the three interpretations, or maybe another, will turn out to be the most fruitful one, only future can tell. Therefore at this stage we content ourselves with developing them to a level that allows to solve the Euler equations, and with stating their respective advantages and disadvantages, whereby we cannot single out any of the three extensions as overperforming the others in all respects.

This paper focuses on the one-dimensional case and defers the case of multiple dimensions to a forthcoming publication in order to keep the size of the paper reasonable. Additionally, in multiple spatial dimensions there are further aspects that need to be taken into account, and the numerical discretization also is different for Cartesian grids and unstructured ones.

Here, the grid is assumed equidistant with cells where , . As customary, indices referring to time are denoted as superscripts, i.e. is the -th time level. denotes the class of times continuously differentiable functions, and the set of polynomials of degree less or equal to .

## 2 Review: Third order Active Flux

In this section, an overview of the existing Active Flux method(s) in one spatial dimension shall be given. At the same time we would like to define Active Flux in the largest possible generality.

Active Flux has two distinctive features. First, the mixed type of degrees of freedom: a cell average, reminiscent of Finite Volume methods, and point values, of which at least some are located at cell boundaries. The numerical solution is thus given by the two sets

(3) |

with the interpretations

(4) |

The same formulae are used to initialize the degrees of freedom given initial data . Traditional degrees of freedom in Active Flux are one cell average per cell, and one (shared) point value located at every cell interface.

Second, the point values are shared between adjacent cells, contrary to the approach of e.g. DG methods. The point values at cell interfaces allow to immediately write down a conservative update of the cell average. Indeed, integrating the conservation law (1) over a computational cell and applying Gauss’ law yields

(5) |

The evolution of the point values can happen in many ways, such that for the moment we content ourselves with the following general definition:

###### Definition 2.1 (Traditional Active Flux).

Inside cell , consider a reconstruction

(6) |

with the properties

(7) | ||||

(8) |

Define the global reconstruction

(9) |

(10) |

The *traditional Active Flux method* is the following semi-discretization

(11) |

of (1) with the interpretations

(12) |

In the following, we often drop the dependence of on its parameters other than . When the time step is important, we denote by the reconstruction obtained by using values of time step .

Usually, but limiting might require other choices, discussed in Section 2.2. An Active Flux method with this choice of reconstruction is at most third order accurate.

In [Abg20], a closely related definition was given, which, however, does not rely on a reconstruction:

###### Definition 2.2 (Semidiscrete Active Flux from [Abg20]).

Its order of accuracy depends on the approximation order of (and eventually on the order of accuracy of the time integration).

Later, (Sections 4 and 5) we discuss extensions of Active Flux to higher order that include other types of degrees of freedom. In these cases, the definition is modified accordingly in order to include their respective updates. This is discussed in the corresponding Sections.

### 2.1 Time integration

Only third-order methods have been suggested so far in the literature, with two different time discretizations for the two methods (11) and (13). In [Abg20], system (13) is integrated with a Runge-Kutta method. The update of the point value is obtained by introducing a point value at the center of the cell and expressing the averages as . This allows to choose for a standard finite difference formula from [Ise82] on a grid of half the spacing:

(14) | ||||

(15) |

For linear advection , the equation for the point value at cell interface then reads

(16) |

To the same order of accuracy, it can be expressed in the form (13) with , as

(17) |

The stability bound for this method is (see Section A).

This method has been successfully applied to the Euler equations in [Abg20] by replacing the prefactor in (16) by the Jacobian and using two finite-difference formulae (biased in different directions) in order to include upwinding.

The earliest version of Active Flux ([vL77, ER11]) employed a leap-frog-type time integrator for the traditional Active Flux (11). First, approximations

(18) |

are computed by solving the IVP (1) as described in Definition 2.1 over half the time step and over the full time step (from the same initial data ), and then these values are used to compute fluxes in the update of the average:

(19) |

Here, are the quadrature weights of Simpson’s rule, which is the adequate quadrature for a third order method. This is a leap-frog-type of integration, because the new point values are already used for computing the new value of the average.

The original publication [vL77] considered Active Flux only for linear advection , in which case one has for positive

(20) | ||||

(21) |

Given the unique parabolic reconstruction which fulfills

(22) |

Equation (20) reads

(23) | ||||

(24) |

with .

This method is stable up to

. It has been extended to nonlinear scalar conservation laws and to hyperbolic systems of conservation laws by means of an estimate of the characteristic directions in

[Bar21]. The main difficulty is to obtain a sufficiently accurate approximate solution of the IVP. In multiple spatial dimensions, the concept of characteristics (along which some quantity is constant) is generally to be replaced by that of characteristic cones, which do not allow for simple transformations to characteristic variables. For linear acoustics, an exact evolution operator has been obtained in [ER13, BK22]. It has been then applied to Active Flux on Cartesian two-dimensional grids in [BHKR19], and structure preservation properties of the resulting method have been demonstrated: The Active Flux method has been found to be vorticity preserving / low Mach compliant. Despite promising efforts towards approximate evolution operators for the multi-dimensional Euler equations in e.g. [Fan17], further research is required for hyperbolic systems in multiple spatial dimensions that is beyond the scope of the present paper.### 2.2 Limiting

Limiting is necessary for Active Flux because it is a high order method, and because it is linear when applied to linear problems. The usual philosophy of limiting in the context of finite volume methods is to modify the reconstruction such that its value at the cell interface satisfies a maximum principle. Examples of such approaches are slope limiters, or the approach in [ZS10] where the reconstruction polynomial is scaled around its average. For Active Flux, it is not possible to combine such ways of limiting with continuity of the reconstruction because the values at cell interfaces are prescribed. Limiting strategies for Active Flux that give up continuity can be found in e.g. [HKS19, CHK21].

Here, we prefer to maintain continuity. The traditional Active Flux method uses a reconstruction to evolve the point values, and it is a natural approach to limiting to replace an oscillative reconstruction by a monotone, or less oscillative function. Such limiting strategies based on piecewise defined functions (e.g. a constant and a parabola) have been proposed in [RLM15, BB20]. A simple limiting strategy has been introduced in [Bar21] where the parabolic reconstruction is replaced by a power law. This shall be taken as inspiration here.

Recall the following result from [Bar21]:

###### Theorem 2.1.

*From left to right*: The cell average is larger than the two point values, there exists is no continuous monotone reconstruction, and the parabola is accepted. The cell average is between the two point values, but the parabola is not monotone, we choose one of the power laws. The average is inside the region given by the dashed lines, we choose the (monotone) parabola. The average is below the region, we choose the other power law. The average is below both point values, we choose the parabola again.

Thus, parabolic interpolants can violate the maximum principle if the average is close to one of the point values . The power law reconstructions, proposed in [Bar21]

(27) | ||||

(28) | ||||

are interpolating the same data while being monotone whenever the data are, and are thus a good replacement for the parabolic reconstruction whenever the latter fails to be monotone. Observe that the expression appearing in (26) is just the power in (27). Thus, when

(29) |

the power law (27) reduces to the parabolic interpolant. The same is true for (28) for

(30) |

If , one can switch from the parabola to the power-law (28), and at to the the power law (27):

###### Theorem 2.2.

For data that are not monotone (i.e. violating (25)), a parabolic reconstruction is used. See Figure 1 for an overview.

A rational interpolation (suggested in [HKS19]) is similar in spirit to the power law reconstruction, but its coefficients cannot be computed analytically because of the condition on its average, which requires integration of a rational function.

Note that in all suggestions of limiting available in the literature, only the update of the point values is limited. This provides satisfactory results, but does not exclude the appearance of oscillations, as the finite volume step

(31) |

remains unlimited. A limiting of the finite volume step is subject of future work.

### 2.3 Transonic upwinding

In [HKS19, Bar21], the following problematic situation has been observed with early versions of Active Flux. Consider initial data

(32) |

for Burgers’ equation, such that the discrete data initially are (see also Figure 2)

(33) | ||||||

(34) | ||||||

(35) |

The value of is irrelevant for the discussion.

The exact solution of (32) would be a shock moving at speed . In the update of , if its value 2 is used to estimate the upwind direction, information would be taken from the left, where the reconstruction is identically equal to 2. would thus remain stationary. The same happens to , whose initial value is negative, and thus information is taken from the right. In fact, the reconstruction thus never contributes. The numerical method keeps all point values stationary, while the average in cell 0 will grow, because the flux difference is .

The error is the wrong estimate of the upwind direction at . An exact time evolution of the continuous reconstruction would see it form a shock, and this shock would first move left, and then right, reaching before the end of the time step. At this moment, the direction of information flow at would be reversed, and information would no longer be taken from the right (see Figure 3). Thus, the deficiency of the early versions of Active Flux is due to an erroneous estimate of the upwind direction in the point value update. This estimate needs to take into account the fact that a continuous reconstruction might self-steepen before the end of the time step.

This error is most visible in the transonic case, because the shock has time to form during the time step. Depending on the CFL condition, setups close to transonic might be affected as well. In the worst case, the consequence is a shock that is not moving (and thus violating the Rankine-Hugoniot conditions), with the cell average in its vicinity growing without bound. In setups that are not strictly transonic, an artificially large average at the location of the shock is visible in the form of a spike (see examples in [Bar21]).

In [Bar21] it has been suggested to estimate the upwind direction more carefully. In particular, when solving for the foot point of the characteristic (which, generally, requires to solve a nonlinear algebraic equation), it has been suggested to initialize the iteration with not only the value at the respective cell interface, but also with values from its neighbours. If a shock forms over the time step, one would thus obtain several candidate characteristics, while in a smooth setup the iteration would converge to the unique characteristic. A simple strategy that works well in practice is then to select from the candidate characteristics the one whose speed is largest in modulus. In [Bar21, BB20] it has been demonstrated that this strategy works well even for systems. Some authors (e.g. [CHK21]) suggest to average to values obtained from the candidate characteristics.

## 3 High order via larger stencils

### 3.1 Discretization in space

The first strategy for high order accuracy that we propose involves the traditional degrees of freedom (cell averages and shared point values at cell interfaces), and uses larger stencils. The starting point therefore is the semidiscrete Active Flux of Definition 2.2, with a Runge-Kutta time integration. Observe again that the update of the averages

(36) |

is exact, and that the order of accuracy of the method is given by the order of accuracy of the right-hand side of the point value update in (13) and of the order of accuracy of the time integration.

Consider linear finite difference -th order approximations to the derivative of the form

(37) |

We shall adopt the following tableau notation for such a finite difference formula that involves point values at cell interfaces and cell averages

###### Notation 3.1.

(38) | |||

The vertical lines indicate cell interfaces (such that coefficients of cell averages are written “inside” the cell and coefficients of point values are associated to cell interfaces). The double vertical line indicates the cell interface at which the finite difference is supposed to provide an approximation to the derivative (of the corresponding order of accuracy). Coefficients not marked explicitly in the tableau are assumed zero. The notation will be slightly expanded in Section 5, where the superscript on the coefficients of the cell averages will become clear.

Here are some examples ( are free parameters):

(39) | ||||

(40) | ||||

(41) | ||||

(42) |

Note once more that the double vertical line indicates the cell interface at which the finite difference is providing the desired approximation of the derivative. E.g. for FD1,

(43) | ||||

while for FD3a, | ||||

(44) |

The finite difference formulas are involve one degree of freedome more than it would be necessary to obtain the desired order of accuracy; the remaining free parameter is used later to optimize stability. These finite difference formulas are different from standard finite differences because they involve point values and averages. Table 1 summarizes some of the possible choices. Given a finite difference formula

(45) |

define the flipped formula

(46) |

It is an approximation of the derivative at the same location and of the same order of accuracy as , and arises because the spatial derivative changes sign upon reflection .

We propose to use these finite differences in the update of the point values by writing

(47) |

where is any of the finite difference formulae (39)–(42). In the scalar case the positive and negative parts are given by

(48) |

In the case of systems,

is the Jacobian and they are defined via its eigenvalues:

(49) |

For an Active Flux method of order , it is sufficient to use a derivative approximation of order . This offset is due to the offset in the definition of “order”. For example, a 3rd order Active Flux method is exact on parabolic data, while already a 2nd order derivative approximation is exact for parabolae.

It remains to define . For scalar conservation laws (), if does not switch sign in the computational domain, then the natural choice works. Further details are given in the next section.

The use of finite differences in the point value update is inspired by the approach in [Abg20], but different in the following crucial aspect. In [Abg20], a point value at the cell center is first estimated through Simpson’s rule:

(50) |

One thus obtains a grid of half spacing, on which standard finite differences can be used. (This is described in Section 2.1.) While for third order, Simpson’s rule is accurate enough, for higher orders the natural expression of the cell average in terms of point values would require solving a linear system, i.e. it would become nonlocal. Indeed, a formula such as

(51) |

cannot be as easily solved for the midpoint values as (50). One could use a quadrature for the average that involves precisely one midpoint value and correspondingly more of the point values at cell interfaces:

(52) |

It is clear, however, that the form (37) includes this case, and is more general.

### 3.2 Time integration and stability conditions

The algorithm consists of the ODE system (36) and (47). We propose to solve it with a Runge-Kutta time integrator, e.g. a strong stability preserving method of order 3 (SSP-RK3). For linear problems, of course, higher order Runge-Kutta methods are easily available. (Compare Figures 11 and 12, where stability of FD5b is shown for RK3 and RK5, respectively, with the maximum CFL number being only slightly higher for RK5.). For nonlinear problems, often a lower order time integrator is used with a time step small enough, such that the global error is dominated by the spatial accuracy. We thus content ourselves with results on RK3. The finite difference formulas considered above are typically one cell larger than it would be necessary to obtain the desired order of accuracy, such that the free parameter can be used to optimize stability. The results of the linear stability analysis (see Section A) can thus be depicted in a diagram whose axes show the CFL number and the value of the free parameter. Some of these diagrams are shown in Figures 6–15. Table 1, shows the maximum CFL numbers for parameter values optimizing stability.

name | FD1 | FD3a | FD3b | FD3c | FD4a | FD4b |

1 (2) | 3 | 3 | 3 | 4 | 4 | |