The Active Flux scheme for nonlinear problems

The Active Flux scheme is a finite volume scheme with additional point values distributed along the cell boundary. It is third order accurate and does not require a Riemann solver. Instead, given a reconstruction, the initial value problem at the location of the point value is solved. The intercell flux is then obtained from the evolved values along the cell boundary by quadrature. Whereas for linear problems an exact evolution operator is available, for nonlinear problems one needs to resort to approximate evolution operators. This paper presents such approximate operators for nonlinear hyperbolic systems in one dimension and nonlinear scalar equations in multiple spatial dimensions. They are obtained by estimating the wave speeds to sufficient order of accuracy. Additionally, an entropy fix is introduced and a new limiting strategy is proposed. The abilities of the scheme are assessed on a variety of smooth and discontinuous setups.


Introduction
Hyperbolic m × m systems of conservation laws in d spatial dimensions have the form The function f is called the flux. Exact solutions of these equations are unavailable in general, and one needs to resort to numerical methods. Cell based methods consider the computational domain to be partitioned into cells. A certain number of discrete degrees of freedom are associated with every cell: e.g. finite volume methods store the average of the dependent variable and spectral/Galerkin methods store coefficients of a decomposition in some basis.
In order to evolve the cell average, finite volume methods require the knowledge of the flux through the intercell boundary (see section 2.1 for a derivation). This flux cannot generally be approximated by a symmetric average of fluxes associated to the values in the two adjacent cells, because this results in an unstable method. Instead, the fact that hyperbolic PDEs have certain preferred directions of information propagation needs to be reflected in the numerical method. The choice of the numerical flux as an asymmetric average of the neighbouring values is referred to as upwinding. It has been suggested in [God59] to use an exact short-time solution as a building block in order to find a numerical flux that leads to a stable scheme. First, a piecewise polynomial function is found, such that it is continuous in every cell and its average agrees with the given average. The discontinuities at cell interfaces present so-called Riemann Problems, which then are solved over a time interval that does not allow them to interact. The exact flux at the location of the cell interface then is used as a numerical flux in the finite volume method. To save computation time, an approximate solution of the Riemann Problem can be used, see e.g. [Roe81,HLL83,JX95] as well as [LeV02,Tor09] for more details. Higher order of accuracy is achieved by widening the stencil ([VL77, CW84,TT02]).
Galerkin methods represent the numerical solution in e.g. a polynomial basis. Every basis coefficient is then evolved using the weak formulation of (1). Again, for hyperbolic equations this requires modification in order to achieve a stable method, one of which is the discontinuous Galerkin method ( [CS98], but see also [BH82]). The basis functions are piecewise polynomial, and in order to deal with the jumps across cell interfaces a Riemann solver is invoked. Higher order of accuracy is achieved by retaining more coefficients of the basis decomposition.
Even in one spatial dimension, conservation laws (1) therefore pose a number of challenges to numerical methods. This does not only include the necessity of upwinding. It is also known that continuous solutions do not generally exist for all times, and thus numerical methods need to be designed in such a way that they can capture discontinuities (weak solutions). Weak solutions in one spatial dimension only become unique upon additional conditions (entropy conditions), and numerical methods need to fulfill a discrete counterpart of these conditions (entropy stability, see e.g. [Tad03]). These aspects have been subject of numerous investigations, see e.g. [LeV02] for an introduction.
In multiple spatial dimensions, systems of conservation laws have a rich phenomenology which is absent in the one-dimensional case. In the context of the Euler equations these are vortices (e.g. created by Kelvin-Helmholtz instabilities), multi-dimensional shock interactions, the low Mach number/incompressible limit and many more. The easiest way of extending a one-dimensional numerical method to multiple dimensions is directional splitting, i.e. the problem is replaced by a number of one-dimensional problems. This, however, has been demonstrated to require excessive grid refinement in order to capture truly multi-dimensional features even for systems much simpler than the Euler equations (e.g. [MR01,GM04,Bar19,BK20]). It has been found that numerical methods should reflect essential properties of the solution at discrete level in order to avoid expensive grid refinement. Such methods are called structure preserving. So far, modifications of existing schemes have been suggested, but it is largely unexplored how such schemes can be derived from first principles.
The Active Flux scheme is a new scheme ( [ER13], an extension of [VL77]) that combines a finite volume scheme with additional, independently evolved degrees of freedom which are interpreted as point values. These point values are located at cell boundaries and Active Flux thus uses a continuous reconstruction. This is a major difference to finite volume schemes. The Active Flux scheme is nevertheless able to resolve shocks (which are approximated by steep gradients), as can be seen below. The reconstruction is parabolic and therefore Active Flux is third order accurate. It has been shown in [BHKR19] for the equations of linear acoustics that the scheme is vorticity preserving without any fix which makes it a good candidate as a structure preserving method for more complicated multi-dimensional problems.
So far, the Active Flux scheme has been studied in great detail for linear equations ([VL77, ER13, BHKR19]). As explained in section 2, the essential ingredient is an approximate solution operator for the initial value problem which is used to update the pointwise degrees of freedom. For linear equations, the point values can be updated using an exact evolution operator. For nonlinear problems, an approximate evolution operator is required. By exploiting special properties the Active Flux scheme has been applied to Burgers' equation ( [ER11b,ER11a,Roe17]) and Euler equations ( [ER11b,Fan17,Mae17,HKS19]). Some of these extensions lose the order of convergence when applied to nonlinear equations. The approximate evolution operator has to be of sufficiently high order, e.g. local linearization as used in [ER11b] is not sufficient to yield an overall third order scheme. In [HKS19], a solution operator based on the Cauchy-Kovalevskaya/Lax-Wendroff procedure has been suggested which can be applied to general nonlinear hyperbolic systems in one spatial dimension and has shown the correct order of convergence when applied to one-dimensional Euler equations in practice. However, both the procedure of [HKS19] itself and the evaluation of the higher order spatial derivatives can be rather complicated. In particular, the derivatives are required at locations where the reconstruction is not differentiable.
The aim of this paper is to provide a simpler solution operator that allows to apply the Active Flux scheme to a large class of hyperbolic conservation laws. The general idea is to keep the structure of a characteristic-based (or in multi-d characteristic-cone-based, see [BHKR19]) evolution operator but to estimate carefully the wave speeds -which are not constant in the nonlinear case. This also includes an estimate on whether a shock has occurred by self-steepening. In section 3 such approximate evolution operators are provided for scalar conservation laws in one and several spatial dimensions, and in section 4 -for hyperbolic systems of conservation laws in one spatial dimension. This leads to algorithms significantly different from the scalar case. This paper is the first part of a sequence of papers devoted to the application of the Active Flux scheme to nonlinear problems. The case of multi-dimensional systems shall form the content of a forthcoming work. Although the examples presented here are all computed on Cartesian grids, the approximate evolution operators can be immediately applied to unstructured grids. A detailed experimental study concerning unstructured grids, however, is subject of future work.
Section 5 is describing a limiting procedure. As the Active Flux scheme is of higher order, spurious oscillations can appear. The continuous reconstructions employed in the Active Flux scheme do not allow to make immediate use of the same limiting strategies as in the case of usual finite volume schemes. Several limiters for Active Flux have been suggested in the literature: in [RLM15,HKS19] the parabolic reconstruction inside a cell is replaced by several parabolae joined in a continuous and monotone manner. This, however, might add implementational and computational complexity. The same is true for the hyperbolic reconstruction considered in [HKS19] as its parameters cannot be computed analytically. Additionally, discontinuous reconstructions have been considered in [ER11b, Eym13, HKS19] as limiting strategies. However, favorable properties have been deduced from the continuous reconstruction in [BHKR19, Bar19], and these discontinuous limiting strategies violate the principle of Active Flux and move it again closer to usual finite volume schemes. This shows the need for a simple limiter that keeps the reconstruction continuous. Here, such a limiter is presented -it is optimally monotone (see Theorem 5.2 for precise statement) and is at the same time computationally efficient.
The paper thus is organized as follows: Approximate evolution operators are presented in section 3 for scalar conservation laws and in section 4 for systems. Limiting is discussed in section 5 and numerical examples for problems in one and two spatial dimensions are shown in section 6.
2 The Active Flux scheme 2.1 Finite volume scheme Consider the computational domain to be divided into (polyhedral) computational cells C ⊂ R d and discretize the time into points t n , n ∈ N 0 separated by (not necessarily equal) time steps ∆t. Recall that in order to solve (1), finite volume schemes use cell averages 2 as discrete degrees of freedom. The cell average is updated in time using fluxes f e through cell boundaries (edges e):q Here,q n C denotes the value of the average at time t n . Applying Gauss' law to (1), f e can be given the interpretation of approximating with n e the outward normal to edge e. In the one-dimensional (d = 1) case the cells are indexed by a finite subset of the integers. Thenq i denotes the averages in cell ], which is centered around x i . The size |C i | of a cell is for simplicity denoted by ∆x. Most of the results remain valid when the size of the cells varies smoothly.

Pointwise degrees of freedom
The Active Flux scheme is an extension of the finite volume scheme (3) for (1). Active flux uses point values q x located at points x ∈ R d along the cell boundary as additional discrete degrees of freedom (recall that indices never denote derivatives in this paper). They approximate the value q(t, x). So far, the following choices have been considered in the literature (see also Figure 1  Active flux is not a staggered-grid finite volume method. Staggered grids consider offset grids of averages, each for a different variable; Active Flux stores additional point values of all the variables inside every cell. This approach is closer to Lagrange-basis spectral/Galerkin methods, with the difference that the cell average is retained as one of the degrees of freedom. A particularity of the Active Flux method is the exclusive distribution of the point values along the cell boundary. For the evolution of the point values at cell boundaries the Active Flux scheme considers an initial value problem: a reconstruction q recon (x) plays the role of the initial data, and the evolution in time can be either exact or approximate. This is explained in more detail in the next sections.
Once the time evolution of the point values is known, the numerical flux is obtained using a quadrature of (4) in time and along the edge. In order to obtain a third order scheme it is necessary to compute the point values also at half the time step (see [BHKR19] for further implementation details). Only the update of the average needs to be conservative, there is no notion of a conservative update for a point value.

Reconstruction
The reconstruction is interpolating the point values and the average: The reconstruction thus is conservative. The difference to reconstructions in the context of finite volume schemes is the fact that the reconstruction is continuous at the locations of the pointwise degrees of freedom. Additionally, the choices used for the reconstruction so far in the literature were such that the reconstruction is continuous everywhere in the computational domain. These particular choices are briefly reviewed next: • In one spatial dimension, the reconstruction is chosen piecewise parabolic in [VL77]. This is a natural choice, as (5)-(6) amount to three conditions in each cell. It reads The reconstruction is continuous everywhere.
• In two spatial dimensions, in [ER13,BHKR19] the pointwise degrees of freedom are placed at endpoints and at the midpoints of every edge. The reconstruction is chosen always to reduce to a parabola along any edge and, as a parabola is uniquely defined by three points, the reconstruction is thus continuous across any edge, and thus everywhere.

Evolution of pointwise degrees of freedom
Active flux is a time-explicit method and thus subject to a CFL condition ∆t < L min λ max (9) In the following, the time step is chosen based on the maximum value λ max of the characteristic speed at the location of pointwise degrees of freedom. The shortest length L min in the one-dimensional case is the size of the cell, and in the two-dimensional case half the edge length (as there is a pointwise degree of freedom located at its midpoint). The update procedure of the Active Flux scheme for the point value is the (exact or approximate) solution of the initial value problem at its location. The initial data are given by the reconstruction. When the Active Flux scheme is applied to linear equations (as in [VL77,ER13,BHKR19]) an exact evolution is easily available. For nonlinear equations it is necessary to devise approximate evolution operators. This is the topic of section 3 (for scalar nonlinear conservation laws) and section 4 (for systems of conservation laws). Here, only a general statement shall be given that concerns the necessary accuracy of an approximate evolution operator. First, the following result is needed: Lemma 2.1. For f : R × R → R and g : R → R, both analytic and n 1 , n 2 ∈ N, assume Proof. Expand Recall that f ∈ Θ(g) means that asymptotically c 1 |g| ≤ |f | ≤ c 2 |g| for some c 1 , c 2 > 0.
Theorem 2.1. Assume a hyperbolic CFL condition ∆x ∈ Θ(∆t) as ∆t → 0. If the approximate evolutionq(t, x) for fixed x ∈ R approximates the exact solution q(t, x) at least asq and the quadrature rules used to approximate (4) yield the exact value up to an error of O(∆t α ∆x β ), α + β ≥ 3 then Active Flux formally achieves third order accuracy.
Proof. Denote by T t [q 0 ] the exact evolution operator applied to initial data q 0 and evolving them to a time t, and byT t [q 0 ] its corresponding approximation.
Assume point values of q(t n , x) to be used in the reconstruction. Then, because the reconstruction is an interpolation, and taking x to be the location of one of the point values (where the interpolation is exact) This statement can also be understood as follows: the interpolation matches the Taylor series of q(t n , x + δx) in δx to sufficiently high powers of δx. At the same time, the derivatives if q that appear as coefficients in this Taylor series are approximated by finite differences, which carry error terms O(∆x β ), β > 1. Because they all use the same point values in the approximation, lower order derivatives are approximated better. The approximate evolution operator uses initial data from the neighbouring cells at a distance O(∆t) from some fixed x. Thereforẽ For the average update, the numerical flux is obtained using a quadrature of (4), such that the numerical flux differs from the exact one by the same error. Thus, using the assumption of a hyperbolic CFL constraint, lemma 2.1 implies that the leading errors cancel when the fluxes at x + ∆x and x are subtracted. One is left with This on total gives a numerical method of third order.

Overview of the algorithm
The overall algorithm of Active Flux is as follows: 1. Given cell averages and point values, compute a reconstruction according to section 2.3.
2. Use the reconstruction as initial data in the update of the point values (section 2.4). Approximate evolution operators for scalar nonlinear problems are discussed in section 3 and for nonlinear systems in one spatial dimensions in section 4 below.
3. Given the updated point values along the cell interfaces, compute the intercell fluxes via quadrature of (4). Here, a space-time Simpson rule is used.

Scalar nonlinear equations
Consider the initial value problem for the following scalar (i.e. m = 1) conservation law Assume the flux function to be smooth and convex.

Fix-point iteration
In the absence of shocks (19) can be rewritten as with a(q) = ∂ q f (q). The characteristics ξ : R + 0 → R d are straight lines on which the solution is constant. They fulfill The exact solution is found by evaluating the initial data at the footpointξ = ξ(0) of the characteristic as q remains constant along it. This also allows to write x =ξ + a(q 0 (ξ))t This equation can be solved forξ efficiently using a fixpoint iteration: , given recursively bŷ for t ≥ 0, formally approximatesξ to n-th order, i.e.ξ which proves the assertion.
This iteration seems related to the Picard iteration, but it is exact for linear problems for any initial data after one step. Therefore the above iteration is even more powerful than a standard Picard iteration.
In view of Theorem 2.1, an evolution operator for the discrete degree of freedom q x located at x is and instead of q 0 the reconstruction q recon based on values at time t n would be used in the fixpoint iteration.

Comparison to previous results
Before turning to questions regarding the possible presence of shocks in the solution, compare this evolution to similar approaches available in the literature. Note that (32) estimates the speed of the characteristic as Local linearization would correspond to taking the evolution operator q recon (ξ (1) ), and thus estimate the characteristic speed simply by a(q 0 (x)). For the special case of Burgers' equation, in [ER11b] it is suggested to estimate the characteristic speed in one spatial dimension by However, this approach does not lead to an increase in the order of convergence (as can be shown by direct computation) and thus is not fundamentally superior to local linearization. In [Roe17] the exact speed of the characteristic for linear data is used as an estimate.
Usage of this formula as an evolution operator for the pointwise degrees of freedom requires the evaluation of the derivative at a location where the data are not differentiable. Also, equations with more complicated wave speeds lead to a lot more complicated formulae.

Modification of the fixpoint iteration in order to account for shocks
It is well-known that nonlinear hyperbolic equations develop shocks even when the initial data are smooth. Therefore, when studying the time evolution of the reconstruction, the assumption that no shocks appear cannot always be true. However, the reconstruction is continuous. An initial value problem with continuous data does not develop a shock immediately. The shock can only appear only after a time t s > 0. Whenever the time step happens to be small enough (∆t < t s ), the reconstruction did not have time to develop a shock and (32) is a good estimate. This gives an explanation why in certain cases even Riemann problems can be successfully computed with the Active Flux scheme endowed with (32). Fig. 2 shows such a successful computation of a Riemann problem between values q high = 11 and q low = 1 for Burgers' equation. Recall that the initial data in the cell containing the discontinuity are still reconstructed continuously. One can estimate its self-steepening time in this situation as . The CFL condition involves the maximum speed q high in this case. Thus, for a Riemann problem with uniformly positive values the time step is always smaller than the estimate of the self-steepening time.
Riemann problems involving both positive and negative values do show artefacts. In [HKS19] it has been shown, that on such Riemann problems for Burgers' equation evolution operators like the one from [Roe17] fail. (32) suffers from very similar problems. In [HKS19] it is suggested to revert to a discontinuous reconstruction in this case. However, the failure can be explained by an insufficiently accurate evolution operator rather than tracing it back to continuity of the reconstruction. In order to do this, consider an even simpler Riemann problem for Burgers' equation: The exact solution is a shock moving at speed 1 2 . However, the evolution operator using (32) leaves these data stationary! Indeed, for x > 0 the fixpoint iteration is initialized with zero speed. For x < 0 the fixpoint iteration converges after one iteration toξ (n) =ξ (1) = x − t and q 0 (x − t) = 1 ∀t.
Additionally, transonic rarefactions show non-entropic artefacts similar to the ones observed for finite volume schemes with Riemann solvers. Examples of such are shown in Figure 10.
All these problems are removed by modifying the initialization of the fixpoint iteration (32) as follows:ξ On two-dimensional Cartesian grids, and q n+1 The reasoning behind this algorithm is the following: Shock formation occurs because of crossing characteristics. Iteration (25)-(26) converges to both footpointsξ + andξ − if its initial estimatesξ (0) ± are chosen appropriately. The above algorithm (37)-(41) initializes the iteration with the two locations x±∆x, placed symmetrically around x. In order to find the solution, one thus needs to estimate which of the characteristics will have survived until time t (and not gone into the shock). The choice of (41) is to use the value transported by the quicker characteristic. This choice is inspired by the above example (36) of a Riemann problem, where it makes information flow into the right direction. Of course in general it remains an approximation.
Note that the order of the approximation is not modified, as the modification affects only the initial step of the iteration and is of the order O(∆x). On smooth solutions, characteristics do not cross, and the two initializations are expected to converge to the same final result.
Despite its simplicity, in experiments the modification has proven itself able to reliably cure both the artificially stationary shocks and the non-entropic features at transonic rarefactions. One thus may speak of the modification as an entropy fix. To actually prove a statement on the discrete entropy is subject of future work. Instead here a number of different test cases are shown: self-steepening and Riemann problems resulting in strong and weak shocks and (transonic) rarefactions (see sections 6.1-6.2).

Nonlinear systems
Consider now an m × m nonlinear hyperbolic system of conservation laws in one spatial dimension: The Jacobian matrix is denoted by J(q) := ∇ q f . Hyperbolicity guarantees that J has real eigenvalues.
In certain cases a variable change from conservative to characteristic variables q → Q can be found, such that in the absence of shocks (42) can be rewritten as with λ 1 , . . . , λ m the eigenvalues of J. Denote the initial data as Q i,0 (x) = Q i (0, x), i = 1, . . . , m.
In the linear case this can be solved by solving an advection equation in every component (see e.g. [ER11b]). In the nonlinear case λ i is, in general, a function of all the components of Q: Therefore, in general the characteristics are curved. This is also a fundamental difference to the nonlinear scalar case, where the characteristics remain straight. This is why applying the fixpoint iteration (32) to every component of (44) does not lead to sufficient order of accuracy (as can be checked by direct computation). A different approximate evolution operator is necessary in the case of systems, which takes into account the curvature of characteristics. Sections 4.1 and 4.2 describe two such approaches. They yield comparable results, but the strategies and resulting algorithms are fundamentally different. In particular, the algorithm in section 4.1 does not assume a transformation to characteristic variables. In view of future extensions, e.g. to multiple spatial dimensions, so far it is not clear which of them would be most suitable. They also differ in the nature of necessary computations. Therefore both are presented here.
Even in case the approximate evolution operator is formulated in characteristic variables, the reconstruction still uses conservative variables (as it requires a cell average). At the locations where the initial data need to be evaluated, a transformation to characteristic variables is performed. After obtaining the result of the approximate evolution operator in characteristic variables, they are transformed back to conservative variables.

Estimating curved characteristics
It can be shown by explicit calculation that a straightforward extension of iteration (25)-(26) to (44) does not allow to prove a statement analogous to Theorem 3.1 -the higher order terms are not correct. This is due to the fact that characteristics are now curved. However, other ways of obtaining a high order estimate can be found. The first is presented in this section, the second -in section 4.2 Consider (42) and diagonalize RJR −1 = Λ = diag(λ 1 , . . . , λ m ). Note that, in general, R and all λ i depend on q. To make this explicit, write R(q) and Λ(q). The equation becomes Although for some systems it is possible to find Q(q) such that ∂ t Q = R(q)∂ t q (as mentioned in the introduction to this section), this is not assumed in the following algorithm.
Denote by F (k) = R −1 diag(0, . . . , 0, k 1, 0 . . . , 0)R the projector associated with the k-th eigenvalue. Then obviously In the following, matrix indices are frequently made explicit, e.g. q i denotes the i-th component of the vector q, and R ij the element of R found in its i-th row and j-th column. It is also sometimes useful to write λ i (q 1 , . . . , q m ) instead of λ i (q).
where λ i (x) is shorthand for λ i (q 1,0 (x), . . . , q m,0 (x)) ∀i and analogously for F (k) (x). Then, with the approximate solution operator approximates the exact solution of (42) with J = R −1 diag(λ 1 , . . . , λ m )R as Note: The inconspicuously looking equation R * ij := R ij (q (i) ) is non-trivial. It states that the rows of R * are evaluated independently, each on a different predictor value.
Proof. Wherever the summation is from 1 to m, it is omitted for the sake of readability. Recall and note that (The prime denotes differentiation with respect to the unique argument.) In order to compare the leading order terms in the Taylor series, differentiate the approximate evolution operator with respect to time: On the one hand then Here λ * i | t=0 = λ i , R * | t=0 = R was used. On the other hand The term in brackets can now be expanded using the definitions of R * and λ * : and with (56) and h,s was used.
On the other hand, by performing the Cauchy-Kovalevskaya/Lax-Wendroff procedure on the PDE, which proves the assertion.
It makes sense to express the matrices R in variables which make the computation simple. In the numerical examples for the full Euler equations, (ρ, v, p) are used: The transformation matrix R (from (ρ, v, p) to the eigenspace (+, 0, −)) reads This gives Proof. Assume p = const, v = const uniformly. Then , and the index is taken from the symbolic set i ∈ {+, −, 0}.
This completes the proof.

Runge-Kutta scheme
Recall the second order Runge Kutta method for an ordinary differential equatioṅ x * = x(0) + αtλ(0, x(0)) (79) For α = 1 2 this can be simplified to the midpoint method For simplicity of presentation, consider m = 2. The Runge-Kutta integration can be applied to the characteristic relations that govern the time evolution of the characteristic curves ξ i : R + 0 → R, i = 1, 2.

Limiting in one spatial dimension
The Active Flux scheme uses a conservative reconstruction in order to evolve the pointwise degrees of freedom. It has to fulfill several conditions (i.e. (5)-(6)). In one spatial dimension conservation and interpolation of the two point values at cell boundaries amount to three conditions. The natural choice therefore is a parabola (e.g. in [VL77]). However, polynomials in general do not fulfill a maximum principle: The maximum of the reconstruction q recon (x) in cell i can exceed max(q i , q i+ 1 2 , q i− 1 2 ). To correct this (whenever possible) is the objective of the limiting procedure introduced below.
The starting point is an analysis of the failure of the parabolic reconstruction to be monotone. Assume in the following that q i− 1 2 < q i+ 1 2 ; for the opposite situation analogous statements are true. ii) The reconstruction (8) has a maximum inside [− ∆x 2 , ∆x 2 ] if the average is too close to the point values. The maximum is located at . Equating this to ± ∆x 2 yields the bounds. Thus, there are three possible cases, which are shown in Figure 3.
: no continuous monotone reconstruction exists: use the parabolic reconstruction.
Proof. Monotonicity and (5) are obvious. For (6) compute 1 ∆x Thus, using the limiter amounts to replacing the parabolic reconstruction by (99) in the region B: The effect is illustrated in Figure 4. . Outside this range no monotone function can be found, and parabolic reconstruction is still used. Right: A version of the limiting symmetric with respect to N → 1 N according to (103).

Notes:
i) Ifq i > q i− 1 2 + r then N < 2. Therefore inside the region C (q i− 1 2 + r <q i < q i+ 1 2 − r) where the parabolic reconstruction (8) is monotone, the power law would have a formal approximation order less than the parabola.
ii) Usage of the parabolic reconstruction whenever an overshoot (undershoot) is unavoidable (region A) is making the limiter affect maxima (minima) as little as possible, and thus to avoid clipping. In practice, limiting is discarded if it would imply max(N, 1 N ) > 50.
iii) As can be seen from Figure 4 (center), the curves for N and 1/N are not symmetric. One thus can consider as a limiting strategy instead of (102). The result is shown in Fig. 4 (right).
iv) For q i− 1 2 > q i+ 1 2 the parabola is monotone ifq i fulfills but the formula (99) remains unchanged.
In Figure 5 the effect of limiting is shown for linear advection. The initial data are compared to the numerical solution on a periodic grid after one revolution. The approximate evolution operators in sections 3 and 4 approximate the solution at time t by evaluating the initial data at some particular location. The above limiting procedure thus guarantees that the point value update satisfies the maximum principle in all cases when a monotone reconstruction is at all possible. This does not mean, however, that the full numerical solution will be free of spurious oscillations. The finite volume update of the average might still lead to the appearance of oscillations. In practice, for linear advection one does observe oscillations for large CFL numbers, but they are much smaller than without limiting. A limiting of the finite volume step for Active Flux has not yet been considered in the literature and is subject of future work.
For systems, limiting is applied to conservative quantities individually. Numerical examples of limiting for nonlinear problems are shown in section 6.

Numerical examples
The Active Flux scheme endowed with the approximate solution operators of the above section is now applied to several problems in order to assess its abilities experimentally. First, scalar equations and systems in one spatial dimension are considered (sections 6.1-6.4); in section 6.6 the scheme is applied to multi-dimensional scalar conservation laws.
The CFL condition is applied using the maximum absolute value of the characteristic speed (eigenvalue of the Jacobian in case of systems) evaluated at the point values (edge midpoint values in the multi-dimensional case).
Third order accuracy in time requires the computation of the point values at both the full and the half time step. Simpson's rule in time then is used to compute the fluxes necessary for the cell averages update (3).

Burgers' equation
Here, Burgers' equation is solved with the Active Flux scheme. In all cases the approximate evolution operator (41) is used with a(q) = q. Additionally, section 6.1.3 shows the artefacts appearing upon usage of the simple, unmodified fixpoint iteration (32), and their absence when using (41).

Convergence study
In order to assess the experimental convergence rate, Gaussian initial data are evolved on grids of different refinement. The simulation is stopped before a shock occurs. Figure 6 (right) shows the setup and the solution at final time on a grid with ∆x = 1/100.

Self-steepening
Figures 7 and 8 show continuous initial data which self-steepen into a shock. In both cases the arising shock connects a positive and a negative value -a situation in which the simple fixpoint iteration (32) is known to fail. Figure 7 (left) shows such a failure: the shock is stationary, which violates the Rankine-Hugoniot condition (this problem has been first reported in [HKS19]). Figure 7 (right) demonstrates that the modified iteration (41) yields the correct shock speed 0.75. Moreover, in figure 9 it is shown that the simple fixpoint iteration (32) produces a spurious oscillation when a shock appears. This is not the usual oscillation due to the high order of the method, as it is removed upon usage of the modified fixpoint iteration (41), even if limiting is not applied. The reason for the appearance of the oscillation is related to the stationary-shock artefact of figure 7 (left). As is emphasized in [HKS19], stationary point values imply fluxes which do not change in time. At the location of the shock, the two fluxes of the cell are different and the average keeps growing. This pile-up is observed as an oscillation in figure 9, although the non-zero value on the left side of the shock is enough to make it travel at the right speed.

Riemann problems
Finally, a number of Riemann problems are used in order to assess the properties of the suggested modification of the fixpoint iteration. Fig. 10 shows an initial setup chosen to include many interesting cases at once, including strong and weak shocks and rarefactions, and shocks and rarefactions where the discontinuity crosses 0. These latter are known to be the most problematic. On the right also some of the constant states are 0. The boundaries are periodic. In the same Figure, a solution at time t = 0.1 is shown which is using the simple, unmodified fixpoint iteration (32). One observes that only uniformly positive or uniformly For comparison, Fig. 11 shows the evolution of the same setup using the modified fixpoint iteration (41). Now all the shocks have the correct speed, and the rarefactions do not contain additional non-entropic shocks. The strong rarefaction with values symmetric around 0 seems particularly difficult to capture. For large CFL numbers, little artefacts have been observed, which, however, vanish upon grid refinement.
To the author's knowledge the performance of Active Flux for nonlinear problems has never been studied on transonic rarefactions or strong shocks in the literature.

Other scalar equations
Consider the following scalar conservation law The speed of a shock joining q L and q R is E.g. a shock joining 1 and −5 is moving at speed −26.
The self-similar rarefaction is given by 3 x t . Fig. 12 shows again a selection of shocks and rarefaction together with the numerical and the exact solution. Usage of the modified fixpoint iteration (41) allows to resolve all the shocks and rarefactions correctly.

p-system
The p-system in case of smooth solutions can be rewritten in the form (44) as with c = p (ρ). In the following, γ = 1.4 is used.
As the eigenvalues cannot switch sign, and never vanish, Active Flux solving the psystem seems less prone to artefacts.

Isentropic Euler equations
Consider the isentropic Euler equations On smooth solutions this system is equivalent to with c = γp(ρ)/ρ. In the following, K = 1, γ = 1.4 are used.
6.4.1 Convergence study Figure 15 shows the setup of a Gaussian initial density and no velocity. Its evolution (before a shock forms) is used to study the convergence of the method. The reference solution is the setup solved with the iteration from section 4.2 on a grid of 16384 cells. Third order is confirmed experimentally.

Riemann problems
The scheme is now applied to Riemann problems. Here it becomes important to use the modification (95)-(98). Figure 16 shows a setup with a transonic rarefaction, accurately resolved, and Figure 17 shows a strong shock and a double rarefaction. The two tests cover the cases of the two eigenvalues of the system having different and same sign. The double shock case shows little artefacts which vanish upon refining the grid.

Full Euler equations
The full Euler equations with the ideal equation of state e = p γ − 1 + 1 2 ρv 2 γ > 1 (118) form a hyperbolic system of conservation laws, but do not admit characteristic variables. Thus, the most general solution algorithm (52) is required. In the following, γ = 1.4 is used.

Convergence study
To demonstrate the convergence of the method, the test from [HKS19] is run: The numerical results are compared at t = 0.25 to a reference solution obtained on a grid of 2048 points. In Figure 18, one observes third order convergence, in agreement with the theoretical expectation.

Riemann problem
To assess the performance of the numerical method on discontinuous problems, two Riemann problems are computed: the Sod shock tube ( [Sod78], Figure 19, left) and the Lax shock tube ( [Lax54], Figure 19, right). One observes the good perfromance of Active Flux even on these discontinuous setups.

Interaction between a shock and a sound wave
Finally, to demonstrate the performance of the algorithm on a more challenging setup, the Shu-Osher test ( [SO89]) is shown in Figure 20. One observes that due to its high order, Active Flux is able to capture the details of the interaction on poorly resolved grids without difficulty (compare e.g. to [CLS89]).

Multi-dimensional scalar equations
As a last test case, consider the multi-dimensional Burgers' equation and a 4-quadrant Riemann problem setup as follows: Figure 21 shows the solution at t = 0.3 using Active Flux (with fixpoint iteration (41)) along with the exact solution taken from [GPP11] (p. 4258). No limiting is used.

Conclusion and outlook
The Active Flux scheme is a finite volume scheme with additional pointwise degrees of freedom located at the cell boundary. It involves a continuous reconstruction and thus does not make use of Riemann solvers. Instead, an evolution operator for the pointwise degrees of freedom is required. Once their evolution is obtained, the update of the cell average follows the usual finite volume/Godunov scheme idea; the intercell flux is obtained by evaluating the flux function on the point values along the boundary (and using quadrature). The Active Flux scheme has initially ([VL77, ER13]) been employing exact evolution operators, because the problems onto which the scheme was applied admitted closed form solution operators.
In particular it has been shown in [BHKR19] that for linear acoustics the Active Flux method is low Mach number compliant and stationarity preserving without the need for any fix. This makes Active Flux an interesting candidate for a class of methods, which are structure preserving by construction. Inspired by the finding for linear acoustics, this paper serves as a stepping stone towards deriving structure preserving Active Flux methods for multi-dimensional nonlinear systems.
Upon an extension of the Active Flux method to nonlinear problems usage of exact evolution operators cannot be maintained. Approximate evolution operators suggested so far in the literature either were reducing the order of the scheme or involved complicated expressions such as the Lax-Wendroff expansion and subsequent solution of Riemann problems (ADER). This paper shows how approximate evolution operators can be found which are not costly and allow to maintain third order of accuracy. The cases considered here are nonlinear scalar conservation laws in one and two spatial dimensions as well as nonlinear hyperbolic systems of conservation laws in one spatial dimension.
As Active Flux is very different from standard finite volume or Galerkin methods, many aspects need to be reconsidered, and many questions well-studied for other methods were still open. Continuous reconstruction might raise doubts about the applicability of Active Flux to problems involving shock formation. It is found that in certain setups too simple an evolution operator fails to correctly recognize the self-steepening. This then leads to artefacts that resemble entropy glitches encountered with certain finite volume schemes. As a cure, in this paper a simple strategy is presented which allows to take into account the fact that characteristics may cross. This "entropy fix" is found to lead to accurate numerical evolutions without artefacts. By means of numerical examples it is shown that Active Flux, for example, is able to accurately solve Riemann problems for one-dimensional systems of nonlinear conservation laws, such as the Euler equations.
As a numerical method of higher order is prone to overshoots around discontinuities, a limiting procedure needs to be in place. Here, a simple limiting is suggested which modifies the reconstruction whenever an avoidable overshoot/undershoot is recognized. Whereas the Active Flux limiters available in the literature either require joining several polynomials inside the cell or re-introduce discontinuities, the suggested monotone reconstruction is simple to compute and retains a continuous reconstruction.
The correct approximation of the entropy solution and limiting in one spatial dimension may not outperform currently available methods of third and higher order. However, all these are necessary ingredients for an extension to multiple spatial dimensions that so far were open, or at least insufficiently studied for the Active Flux method.
Future work shall be devoted to multi-dimensional hyperbolic systems. The approximate evolution operators presented here shall be extended to multiple spatial dimensions and thus combined with the favorable properties of Active Flux in multiple dimensions. This hopefully will pave the way towards a powerful structure-preserving method for multidimensional systems of conservation laws.