A New ADER Method Inspired by the Active Flux Method

We study numerical methods that are inspired by the active flux method of Eymann and Roe and present several new results for one and two-dimensional hyperbolic problems. For one-dimensional linear problems we show that the unlimited active flux method can be interpreted as an ADER method. This interpretation motivates the construction of new third order accurate methods for nonlinear hyperbolic conservation laws. In the two-dimensional case, equivalent methods are only obtained for scalar linear problems. For two-dimensional linear systems the methods are no longer equivalent. For the two-dimensional acoustic equations we compare the accuracy of the two resulting approaches. While commonly used methods for hyperbolic problems are based on discontinuous reconstructions, the active flux method uses a continuous, piecewise quadratic reconstruction. For nonlinear problems we identify a situation in which the continuous reconstruction leads to an unstable approximation. We propose a limiting strategy which overcomes this problem. Our limited version of the active flux method uses the same local stencil as the original method.


Introduction
approximations were also obtained for Burgers' equation. The method has also been used to discretise nonlinear systems of conservation laws, in particular the Euler equations of gas dynamics, by employing a splitting approach [4,14].
Eymann and Roe [4] proposed to reconsider all main components of todays CFD codes. Their goal is to construct methods which are truly multidimensional, more than second order accurate and which have a compact stencil.
The degrees of freedom of the active flux method, which consist of cell average values as well as interface values of the conserved quantities, are chosen in such a way that the method leads to third order accurate approximations using only the most local stencil. Note that interface values are point values of the conserved quantities located at the interfaces. Since the interface values are used by two or more grid cells, this choice of the degrees of freedom reduces memory requirements and computational expense. In the recent paper [14], Roe, Maeng and Fan showed that the active flux method produces accurate results at lower computational costs than other high-order accurate methods. In particular they studied the advection equation, the two-dimensional acoustic equations as well as the pressureless Euler equations.
The choice of the degrees of freedom suggests a continuous reconstruction of the conserved quantities. This continuous reconstruction distinguishes the active flux method from commonly used methods for hyperbolic problems such as Godunov type schemes [8], discontinuous Galerkin schemes [5] or WENO methods [15].
The time evolution of the interface values is performed before the time evolution of the cell average values. Interface value updates should respect the multidimensional physics and the accuracy requirements. A finite volume method is used to evolve average values of the conserved quantities. This guarantees that the method is conservative. The numerical fluxes use the point values at the interfaces at different time levels. The name of the scheme refers to the fact that the interface values are updated independently of the conserved quantities, see [4].
Here we consider methods which use the same degrees of freedom, the same piecewise quadratic reconstruction (whenever possible as explained later) as well as the same finite volume update for the cell average values as suggested in the original active flux method by Eymann and Roe. We show that for linear one-dimensional problems the active flux method can be interpreted as a new version of an ADER scheme, as introduced by Titarev and Toro [17]. We are not aware of an ADER method which uses these same degrees of freedom. In regions where we use the continuous reconstruction, our ADER version does not require the solution of Riemann problems for the conserved quantities but only for certain derivative values of the conserved quantities. The ADER interpretation of the active flux method allows us to construct new third order accurate methods for nonlinear systems of conservation laws, which share many if not all the attractive features of the active flux method. Furthermore, this interpretation provides a strategy for the construction of active flux methods for nonlinear problems which does not rely on a splitting approach. For the Euler equations of gas dynamics, a splitting method might be the best approach for the evolution of the interface values of the active flux method, see Roe [12]. Other authors have also explored splitting methods for the Euler equations, in particular in order to obtain better methods for low Mach number flows, see for example [6,18]. However, in order to construct finite volume methods for general hyperbolic problems, that use the same degrees of freedom as the active flux method, we believe that it is useful to have a general approach for the evolution of the interface values.
Numerical methods for hyperbolic problems use some form of limiting in order to avoid unphysical oscillations near shock waves or steep gradients of the solution profile. In [4], Eymann and Roe suggested a limiting based on neighbouring points of the solution in time as well as space. In [4,12] they point out that the development of appropriate limiting strategies is still ongoing work. In Sect. 5 we present different possibilities for the use of limiters in the active flux method. For linear problems we suggest new continuous reconstructions which use the same stencil as the original active flux method but avoid the occurrence of new minima or maxima in the reconstructed function. These new continuous reconstructions can be used in the original active flux method after some changes, as discussed in Sect. 5.1. They can not be used in the ADER version of the method. For Burgers' equation we identify a situation in which a continuous reconstruction without further modification leads to an unstable approximation. For this situation we propose a discontinuous reconstruction which cures this failure. This limiting approach can be used both in the original active flux method as well as in the ADER version.
It is interesting to note, that 1977 van Leer [20] mentioned a method (scheme V) for the one-dimensional advection equation, which agrees with the active flux method. Although he pointed out advantages in terms of accuracy and efficiency, he wrote: "its value for the ideal compressible flow equations at present seems doubtful". A related but second order accurate finite volume method for advective transport problems, which uses the same compact stencil and distinguishes between the interface update and the update of the cell average values was much more recently proposed by Karabasov and Goloviznin [7]. Furthermore, Zeng [21] proposed related, up to fifth order accurate, hybrid finite difference-finite volume methods.
Our paper is organised as follows: In Sect. 2 we briefly introduce the original form of the active flux method for the advection equation as well as the acoustic equations. In Sect. 3 we show that in the linear case the active flux method can be interpreted as an ADER method. In Sect. 4 we use this interpretation to derive new methods for nonlinear hyperbolic problems and in particular for the Euler equations. In Sect. 5 we present different concepts to introduce limiters for our new one-dimensional schemes. In Sect. 6 we explore the extension of our one-dimensional results to the two-dimensional case. Throughout the paper we show numerous numerical results which illustrate our findings.

Brief Review of the Active Flux Method
We study one-dimensional hyperbolic problems of the form where q : R × R + → R m is a vector of conserved quantities and f : R m → R m is a vector valued flux function. We assume that the system (1) is hyperbolic, i.e. we assume that the flux Jacobian matrix is diagonalisable with real eigenvalues for all q of interest. The degrees of freedom of the active flux method are cell averaged values of the conserved quantities as well as point values of the conserved quantities at the grid cell interfaces. We denote these quantities by Q n i and Q n Since each interface value belongs to two grid cells but is only counted once, the onedimensional active flux method has two degrees of freedom per grid cell. This is the same freedom employed by the discontinuous Galerkin method with linear reconstruction, which can be shown to be the semidiscrete limit of scheme V and active flux respectively, see Roe [11].
In order to describe one time step of the active flux method, we assume that we have third order accurate approximations Q n i+ 1 2 and Q n i at the previous time level t n . The cell averaged values of the conserved quantities at the new time level are obtained via a finite volume approach. On an equidistant grid the method is written in the form where F i+ 1 2 denotes the numerical flux across the grid cell interface located at x i+ 1 2 . The numerical flux F i+ 1 2 is an approximation to the time averaged flux, i.e.
In the active flux method the integral on the right hand side of (3) is approximated by using Simpson's formula We recall that the degree of accuracy of Simpson's formula is p = 3, since p = 3 is the largest integer such that the quadrature formula is exact for x k for all k = 0, . . . , p. This means that we could obtain fourth order accurate approximations of the cell average values if the point values q(x i+ 1 2 , t n + τ ) for τ ∈ {0, t/2, t} would be fourth order accurate. Here the point values are obtained from a piecewise quadratic reconstruction and are thus only third order accurate. This limits the overall accuracy to third order.
Note that the methods considered in this paper only differ in the way how the point values of the conserved quantities are computed. We denote the approximations of the conserved quantities at the grid cell interface by Q n From the cell averaged values Q n i and the interface values Q n i+ 1 2 of the conserved quantities, we can construct a continuous, piecewise quadratic function u(x) for −∞ < x < ∞.
In grid cell i, the reconstructed function has the form )/ x and 0 ≤ ξ ≤ 1. Note that it is the piecewise quadratic reconstruction which limits the method in its current form to third order accuracy. Eymann and Roe [3,4] used characteristic theory in order to compute q(x i+ 1 2 , t n+ 1 2 ) and q(x i+ 1 2 , t n+1 ). For linear hyperbolic problems, they computed exact solutions at (x i+ 1 initial value problem with continuous, piecewise quadratic data u(x) at time t n . We will now illustrate the method for two linear examples.

Advection
For the advection equation with q : R × R + → R and a ∈ R, the conserved quantity at the interface is approximated by For the advection equation with positive advection speed a > 0, the averaged flux at An approximation is only introduced when we replace q(x, t n ) by the piecewise quadratic function u. The spatial integration of the quadratic reconstruction using Simpson's rule provides the exact integral value. Figure 1 illustrates the flux computation.
The active flux method is stable provided that the standard CFL condition a t/ x ≤ 1 is satisfied. In [4], Eymann and Roe explained why this time step restriction is necessary. We are not aware of a formal proof that the method is in fact stable for all according time steps. For the advection equation we can write the method in the form where Q n is a vector which contains all the degrees of freedom at time t n . If we impose periodic boundary conditions, then eigenvalues of the matrix A( t) can be studied in order to understand the stability. While we do not have an analytical formula for the eigenvalues, we plot them for a fixed grid in Fig. 2. All the eigenvalues are inside the unit circle, which experimentally confirms the linear stability of the method. We performed this test for many more CFL numbers in the range 0 < CFL ≤ 1 and show the results of three examples.

Acoustics
In our second example we consider a linear hyperbolic system of the general form q t + Aq x = 0.
As an example for such a linear system we consider the acoustic equations [8], i.e. we set Here p and v denote the pressure and the velocity component in the x-direction, ρ 0 is a constant density, K 0 is called the bulk modulus of compressibility, Z 0 = ρ 0 c 0 is called the impedance and c 0 = √ K 0 /ρ 0 is the speed of sound. Using (5), the piecewise quadratic reconstruction for the acoustic equations is a vector valued function. We use an upper index to denote the two components of u. With the active flux method, the solution of the acoustic equations at (x i+ 1 2 , t n + τ ) is approximated by Formula ( u(x, 0) = 0.
We impose the periodicity condition q(−1, t) = q (1, t). The initial wave package is split into two waves, one moving to the left and one moving to the right with the speed of sound c 0 = 1. Due to the periodicity condition, the exact solution agrees with the initial condition at times t = 2n, n ∈ N. We show numerical results at time t = 7.5, e.g. shortly before the waves overlap in the center of the domain for the fourth time. With the active flux method, the wave package is located at the correct position. This is not the case if we use the Lax-Wendroff method. Furthermore, we note that the numerical results of the active flux method with 200 and 400 grid cells are more accurate than the results of the Lax-Wendroff method with 400 and 800 grid cells, respectively. Recall that the active flux method with 200 grid cells has 400 degrees of freedom. This example clearly shows the advantages of third order accurate methods, in particular if one is interested in the resolution of small scale features. For one-dimensional linear hyperbolic problems, exact evolution formulas can easily be derived. For the nonlinear Euler equations, Eymann and Roe used a splitting approach, which decomposes the Euler flux in a part that describes acoustic wave propagation and a part which describes advective transport. This allowed them to use evolution formulas for the simpler subproblems. We will consider an alternative unsplit approach in Sect. 4.

Alternative Derivation of the Active Flux Method for Linear Problems
In this section we derive a numerical method which shares basic properties with the active flux method. In particular, the method uses the same degrees of freedom and the same piecewise quadratic spatial reconstruction.
While the original active flux method uses exact evolution formulas to approximate q(x i+ 1 2 , t n+ 1 2 ) and q(x i+ 1 2 , t n+1 ), we now instead perform a Taylor series expansion with respect to the time variable and replace time derivatives on the right hand side of (10) by spatial derivatives.

Advection
For the advection equation (6) the time derivatives are replaced by Now we need approximations for q x and q x x at the interface, which we denote by , t n ). Using the reconstruction (5) we compute with Thus, the conserved quantity at the grid cell interface can be approximated by For τ = t/2 and τ = t, we insert these approximations into Simpson's formula and obtain the numerical flux Proof It is enough to show that both approaches lead to the same approximation of q(x i+ 1 2 , t n + τ ). We consider the case a > 0 and assume that aτ/ x ≤ 1. The method from Sect. 2.1 leads to Using (14) and (12), the method from the current section leads to Thus, both approaches are equivalent. The case a < 0 can be considered analogously.

Acoustics
For linear hyperbolic systems of the general form we use the relations to replace the time derivatives in (10) by terms involving spatial derivatives. The vector valued conserved quantities at the grid cell interface are now approximated using the relation For linear systems, the interface values of the derivatives are obtained by solving Riemann problems of the form as well as For the acoustic equations (8), the solution of the Riemann problem with data of the general form consists of three constants states, which are separated by two contact discontinuities moving with speed −c 0 and c 0 . The interface value is given by This formula can also be used to compute the interface values for the Riemann problems (17) and (18). This alternative derivation of the active flux method shares basic properties with both the original active flux method of Eymann and Roe as well as the ADER method of Titarev and Toro [17]. It uses the same degrees of freedom and the same piecewise quadratic reconstruction as well as Simpson's rule for the flux computation as proposed by Eymann and Roe. However, the quantities q(x i+ 1 2 , t n+ 1 2 ) and q(x i+ 1 2 , t n+1 ) are computed in analogy to the ADER approach. It is straightforward to show that both formulations lead to the same numerical results.

Theorem 3.2 For one-dimensional linear hyperbolic systems the ADER version of the active flux method is equivalent to the active flux method that uses the exact evolution formula.
Proof By transforming the system to characteristic variables and using the results for advection, it is straight forward to show that both approaches lead to the same values q(x i+ 1 2 , t n +τ ).
In particular the two formulations of the active flux method for the acoustic equations described in Sects. 2.2 and 3.2 are equivalent. In the next section we will see that the ADER interpretation of the active flux method can be extended to nonlinear problems. The resulting method is no longer equivalent with existing active flux methods.

New Active Flux Methods for Nonlinear Problems
In this section we present new active flux methods for nonlinear problems, which use an ADER approach in order to compute the conserved quantities at the nodes of Simpson's rule.

Burgers' Equation
We start our consideration of nonlinear problems with the Burgers equation, i.e. we consider Our method is again based on the Taylor series expansion of the exact solution (10). For sufficiently smooth solutions, we can replace time derivatives of the conserved quantity by terms involving spatial derivatives, using the relations The conserved quantities at the grid cell interfaces are approximated in the form The derivatives of q at the grid cell interfaces are again obtained by solving generalised Riemann problems. We get with u i (ξ ) and u i (ξ ) as defined in (13).
Evaluations of the right hand side of (20) for τ = t/2 and τ = t provide Q . Now the numerical flux of the resulting finite volume method has the form with f (q) = 1 2 q 2 . Most recently, Roe [12] discussed active flux methods for the one-dimensional Burgers equation. He computes a wave speed at each grid cell interface, using the relation These wave speeds are used to calculate the conserved quantity at the grid cell interface: He compares the approach with a method which instead sets the local wave speed equal to Q n Here we only consider the method with wave speed correction according to (23), since the simpler method without wave speed correction is only second order accurate. We now present numerical results for Burgers' equation with initial data of the form In Fig. 4, we show numerical results at time t = 0.15 for the ADER version of the active flux as well as Roe's version. For this computation we used a grid with 32 grid cells on the interval [0, 1]. The solid line is a highly resolved reference solution obtained by using 2 15 grid cells. The time steps were chosen according to the CFL condition CFL ≤ 0.9, with CFL := max i |Q n In Table 1, we show results of a convergence study using Roe's method with wave speed correction as well as our ADER version of the active flux method. The active flux method with wave-speed correction leads to more accurate results on coarse grids. However, the ADER approach also performs very well and on well resolved grids the accuracy is almost the same.

Euler Equations
In this section we consider approximations of the one-dimensional Euler equations. The Euler equations can be written in the general form of a conservation law (1) with Here ρ, v, p and E denote the density, the velocity component in the x direction, the pressure and the energy, respectively. We use the ideal gas equation of state with adiabatic exponent γ = 1.4.
We explain our approach for the computation of Q , since this is the only part of the method which needs to be adapted to the specific equations. We again want to use the Taylor series expansion (10) in order to compute these interface values of the conserved quantities. For nonlinear hyperbolic systems of the general form (1), we express the time derivatives of the conserved quantities using the relations Thus, the conservative quantities at the grid cell interface can be expressed in the form The exact form of f (q) and f (q) can be found in "Appendix A". The values Q x (x i+ 1 2 , t n ) and Q x x (x i+ 1 2 , t n ) are obtained by solving Riemann problems of the form  as well as To study the accuracy of the resulting method for smooth solutions of the Euler equations, we consider initial values of the form on the interval [0, 1] with periodicity condition. Our time steps satisfy CFL ≤ 0.9. In Fig. 5 we show numerical results for the density at time t = 0.25 on two-different grids together with a highly resolved reference solution. In Table 2 we show the results of our convergence study, which confirm the third order accuracy of the method.

Limiters for the Active Flux Method
Using modified equation analysis, as explained in [19,Chapter 2], it can be shown that methods for the advection equation with an even order of accuracy behave differently at discontinuities or steep gradients than odd-ordered methods. While the leading error term in the modified equation analysis of a second order accurate method creates dispersive waves, the leading error term of a third order accurate method provides damping of spurious oscillations. Thus, in a third order accurate method much less limiting is needed in order to eliminate oscillations. Nevertheless, it is a challenge to find the appropriate amount of limiting. In this section we will discuss different limiting strategies for the active flux method. The piecewise quadratic reconstruction (5), which is used in the unlimited active flux method, can create new minima or maxima. New minima or maxima will occur if u (0)u (1) < 0. Problems need to be expected if there is a relatively large difference between Q n i− 1 2 and Q n i+ 1 2 , while the cell average value Q n i is relatively close to one of these interface values. Such a situation is shown in Fig. 6a.
In order to avoid unphysical oscillations, we want to replace the piecewise quadratic reconstruction, described in Eq. (5), if needed. We will now discuss different reconstructions which are shown in Fig. 6b-d and discuss their use in active flux methods.

Continuous, Piecewise Monotone Reconstructions
In Fig. 6b, c we show reconstructions which are monotone within one grid cell. In Fig.  6b, we show a hyperbola with edge values Q n i− 1 2 and Q n i+ 1 2 and average value Q n i . This reconstruction is motivated by Marquina's work [10], who used a hyperbolic reconstruction in numerical methods for conservation laws. Here, the degrees of freedom of the active flux method are used in the definition of the hyperbola. The reconstruction has the general form The coefficients a, b, c are obtained by solving a nonlinear system of equations, which is obtained from this ansatz by imposing We solved the resulting nonlinear equations using an iterative Newton-type method.
In Fig. 6c we show a continuous, piecewise polynomial reconstruction, which (in each grid cell) consists of a constant part followed by a parabola. This is inspired by a limiter proposed in Roe et al. [13]. The cell average is maintained and the interpolation condition at the two boundary points of the domain is satisfied. This reconstruction is much easier to compute than the hyperbolic reconstruction and can be expressed in the form We now discuss the use of these reconstructions in an active flux method for the advection equation. If the reconstruction does not produce new minima or maxima, the same is true for the exact solution of the advection equation with such data. However, by using the active flux method as described in Sect. 2, the method would in general produce unphysical oscillations. The reason for these oscillations is that Simpson's rule does not integrate the hyperbola or the piecewise polynomial function exactly. The last equality in (7) does not hold if u i is a hyperbola or a piecewise defined function.
For the advection equation we can easily replace Simpson's rule by an exact integration formula, using the relations described in the first two lines of Eq. (7). For the piecewise polynomial reconstruction we can alternatively use a combined Simpson's rule for the two intervals which are divided by ξ * . The new value Q n+1 i+ 1 2 is still computed by using the characteristics but based on the new reconstructions.
These new reconstructions produce very good results as indicated in Fig. 7 (middle) and (right) for a standard advection test problem. In Fig. 7 (left) we show the results of the unlimited active flux method. For this test we used initial data The advection speed has been set to one. We used a grid with 200 grid cells and time steps corresponding to CFL ≤ 0.9. The periodicity condition q(0, t) = q(1, t) is imposed. Thus, at time t = 1 the exact solution agrees with the initial values. Unfortunately, these new reconstructions are not well suited if we use the ADER interpretation of the active flux method. Note that the original form of the active flux method and its ADER interpretation are only equivalent if we use the piecewise parabolic reconstruction (5). For the hyperbolic reconstructions, the derivatives u (k) i (ξ ) do not decrease fast enough for k > 2 and the Taylor series expansion does not converge. Therefore, the ADER approach does not converge. Active flux with exact evolution can still be used. However, the exact flux can no longer be approximated by Simpson's rule, since this would be equivalent to an approach that replaces the hyperbola or the piecewise polynomial reconstruction by a quadratic function which might introduce new minima or maxima.
For the piecewise polynomial reconstruction, the use of the ADER interpretation is only valid up to a time which depends on ξ * . This leads to a reduced time step restriction, with time steps that might become arbitrarily small since ξ * can be located anywhere in the interval (0, 1). The reason for this time step restriction is that the reconstruction is not twice continuously differentiable in each grid cell. Therefore, the Taylor series expansion (10) with (11) is only valid for small time steps.

Discontinuous Quadratic Reconstruction
In Fig. 6d we show a quadratic reconstruction, which has the correct cell average but interpolates only one interface value. We require that the reconstruction uses the interface value which is closer to the cell average value. Furthermore, we require that the first derivative of the reconstruction is equal to zero at that boundary point. This implies that the second boundary point is as close as The discontinuous quadratic reconstruction can be used both for the original version of the active flux method as well as for the ADER version. At interface points where the reconstruction is discontinuous, i.e. at points x i+ 1 2 with u i (1) = u i+1 (0), we need to recompute the interface value Q n i+ 1 2 , since this value is used in Simpson's formula for the flux. This is done by solving the Riemann problem with data The values Q In Fig. 8 we show numerical results for the advection problem from the previous section on three different grids.

Criteria for the Use of a Limited Reconstruction for Linear Problems
For the advection problem we used one of the new reconstructions in grid cell i, if the reconstruction from Eq. (5) produces a new minima or maxima within the grid cell. This is the case, if The only exception in which we continue to use the original piecewise quadratic reconstruction is, if In this situation it is not possible to construct a monotone function with the correct cell average and interface values. Such a situation arises near local extrema of the solution. For a fixed grid, all of our attempts to reduce the number of cells marked for limiting lead to oscillations. In Fig. 9, we show again the results of the advection test with limited,

Limiting Strategies for Burgers' Equation
For the Burgers equation (19), the wave speed at each grid cell interface is determined by the conserved quantity q. This requires a change of the limiting strategy compared to the advection equation. The unlimited active flux method with continuous, piecewise quadratic reconstruction is unstable in situation where q changes sign from positive to negative, as explained in Sect. 5.5.
If u i (0)u i (1) < 0, the piecewise quadratic reconstruction would lead to new minima or maxima. In this situation, we use a discontinuous reconstruction. Near shock waves, i.e. if Q n i− 1 2 > Q n i+ 1 2 , we use a piecewise linear reconstruction of the form with  Fig. 10, we show numerical results using the limited version of the active flux method based on the ADER approach. We computed the solution at different times using time steps corresponding to CFL ≤ 0.8. We compare the numerical solution with results of the 3rd order accurate version of a one-step ADER finite volume method with space-time DG predictor that was kindly provided to us by Michael Dumbser.

Problems of Continuous Reconstructions
We will now describe a situation in which the active flux method with continuous, piecewise quadratic reconstruction fails to approximate the correct solution structure. We consider Burgers' equation and data of the form illustrated in Fig. 11. In this situation, we have a grid cell i with > 0 and Q n update. In the special situation where Q n i = 0 and F i− 1 2 = F i+ 1 2 , the method approximates the correct solution structure. In general, the fluxes F i−1/2 and F i+1/2 will differ and the cell value Q i will increase or decrease indefinitely during the next time steps. Our discontinuous reconstruction offers one possibility to cure this failure. In a personal conversation, Roe pointed out that there might be other possibilities to fix this problem, which still use a continuous reconstruction.

Euler Computations
Finally, we show numerical results for the Euler equations, which require the use of limiters. Again we compare our numerical results with results of the 3rd order accurate version of Dumbser's ADER-DG method.
Our first test problem is Sod's shock tube problem, i.e. we compute solutions of the Euler equations (26) with initial values of the form Our second test problem is the well known Shu-Osher test problem [16], i.e. (26) with initial data In smoothly varying regions near local extrema of the flow, we continue to use the standard reconstruction. Otherwise, we use the discontinuous, piecewise quadratic reconstruction. The discontinuous, quadratic reconstruction will typically be used at the start or the end of a rarefaction wave. The Shu-Osher test confirms that our limited version of the active flux method can capture both the small-scale smooth flow features as well as the shock wave. Furthermore, it compares well with the ADER-DG method of Dumbser. The unlimited version of the active flux method is unstable for the two test problems considered in this section, while the limited version is stable for time steps which satisfy the inequality CFL ≤ 1.

Multidimensional Problems
In this section we discuss possible extensions of our scheme to the two-dimensional case. The equations of interest have the general form where q : R 2 × R + → R m is a vector of conserved quantities and f , g : R m → R m are vector valued flux functions. Again, we assume hyperbolicity of the system. Fig. 16 Left: Configuration in cell (i, j). Right: Area of integration in space-time In this section, we will first describe the setup of our grid, the degrees of freedom of the methods as well as the reconstruction which has been used. We will then see that the equivalence to the active flux method that was derived in the one-dimensional case carries over to the two-dimensional linear advection equation if we solve one-dimensional Riemann problems in a certain way described below. For two-dimensional linear systems we did not find an ADER version which is equivalent with the method proposed by Roe and coauthors. Nevertheless, we will present a third order accurate ADER type evolution procedure which leads to numerical results that compare well with the original method.
Here we will restrict our considerations to two-dimensional Cartesian grids, although other grids such as triangular grids can certainly be used. The original active flux method by Eymann and Roe uses unstructured triangular grids [4,12]. Those grids allow the reconstruction of piecewise quadratic functions with fewer degrees of freedom. Our restriction to Cartesian grids simplifies the comparison of the ADER evolution of interface values with the exact multidimensional evolution that was used by Roe and coauthors. The degrees of freedom that contribute to any given cell will therefore be the four corners and the midpoints of the four edges as well as the cell average. As the edge points are used by two cells and the points in each corner are used by four cells, each grid cell requires the update of four degrees of freedom.
To simplify the calculations, all cells will be mapped to a reference cell [−1, 1] 2 . The cell averages Q n i, j , cell corners Q n with c i j ∈ R, i, j ∈ {0, 1, 2}. This polynomial will more conveniently be expressed in the formû The basis functions N i and coefficients c i , i = 1, . . . , 9 can be found in Table 3. Table 3 Basis functions and coefficients for the two-dimensional reconstruction (37) The resulting piecewise polynomial reconstruction on the whole grid will be denoted by u(x, y). This function is not only continuous in the corner and edge points, but on every point on the edge, as the three points on this edge determine the one-dimensional parabola uniquely. Therefore, we will have continuity of all partial derivatives with respect to x on the horizontal interfaces as well as continuity of all partial derivatives with respect to y on the vertical interfaces. We will use this fact to our advantage later.
In Fig. 16 (right) we show a space time interface for a horizontal grid cell interface. The dots in this figure indicate the nodes of the two-dimensional Simpson formula. The update of the cell average values is described by a finite volume method of the form where F i± 1 2 , j and G i, j± 1 2 are approximations of space-time integrals. At a horizontal grid cell interface, for example, we have These integrals are approximated by a two-dimensional version of Simpson's rule, i.e.
and analogously for F i± 1 2 , j . Note that the same flux approximation is used by the active flux method on triangular shaped grids, see for example [4].
For the advection equation as well as the acoustic equations, exact evolution formulas can be used to update the interface values. Such an approach was used in the original active flux method and will also be used below for comparison.
For general hyperbolic problems exact evolution formulas are not available and we are therefore interested in exploring an ADER type approach. Thus, the update of the cell interface values at the intermediate and new time level is obtained by a Taylor series expansion of the exact solution evaluated at the interface points. Time derivatives are replaced by spatial derivatives using the form of the considered partial differential equation. The spatial derivatives of the conserved quantities at interface points are approximated by solving onedimensional Riemann problems as discussed in more detail below. Especially at the corner points different approximations seem to be possible.

Advection
We start with the advection equation If an exact and easy to evaluate evolution formula is available, then it is of course best to simply use this formula. For general hyperbolic problems such a formula is not available and it is therefore valuable to explore an alternative approach for the approximation of interface values. For the advection equation, Taylor series expansion provides an update formula of the form where the dependencies on the variables in the derivatives were dropped for simplicity. Across a vertical grid cell interface, the reconstructed quantities q, q y and q yy are continuous. Thus, those interface values can directly be used in the Taylor series expansion in order to approximate the interface values at the midpoint (i − 1 2 , j) of this vertical interface. All the other spatial derivatives of the conserved variable are obtained by solving a one-dimensional Riemann problem in the x direction, for example an approximation of q x (x i− 1 2 , y j , t n + 0) is obtained by solving the Riemann problem We obtain Analogously, at a horizontal interface the values q, q x and q x x can be obtained from the continuous reconstruction and all the other spatial derivatives can be computed by solving a one-dimensional Riemann problem. At a corner (i − 1 2 , j + 1 2 ), for each non-mixed derivative, the four initial values of the two-dimensional Riemann problem will collapse into two values due to the same continuity argument. We will therefore approximate this two-dimensional Riemann problem by the resulting one-dimensional Riemann problem as it is done in the edges. We solve two Riemann problem in the x-direction in order to compute q x and q x x . Furthermore, we solve two Riemann problems in the y-direction to compute q y and q yy . The only troublesome terms are the mixed derivatives q xy and q yx . For the advection equation, there is no difference in how these two are treated because of the scalar coefficients. Thus, we only consider q xy . Since we cannot use any continuity of this derivative, we do in fact have four distinct initial values in the general case. We denote these four values by q B L xy := u xy ( Finally, we obtain For the advection equation, the same solution would be obtained by first solving two Riemann problems in the y-direction, followed by a Riemann problem in the x-direction. We recover the correct solution of the two-dimensional Riemann problem. Notice that, as we want to construct a third order method, we typically cut the Taylor series expansion after the second derivative. For the one-dimensional method, all derivatives of the reconstruction past u vanish, i.e. u (n) = 0 , n > 2. This is not the case for our two-dimensional reconstruction. To obtain a slightly more accurate solution, all the non zero spatial derivative terms that arise when we replace q ttt and q tttt by spacial derivatives can be approximated analogously. For our reconstruction on a Cartesian grid, these are the terms q xxy , q xyy and q xxyy . All of these terms can be handled analogously to the terms described above. If all of these terms are included in the Taylor series expansion then the resulting ADER approach is again equivalent to the active flux method with exact evolution formula for the update of the interface values. This is due to the fact that Taylor series expansion recovers the correct characteristic path for our reconstructed values.
In our accuracy test, we solve the advection equation with a = 1 and b = 0.7 on the domain [0, 1] × [0, 1] with periodic boundary conditions. The initial data are given by q(x, y, 0) = sin (4π(x + y)) .
At time t = 1 the numerical solution is compared with the exact solution and the error as a function of the mesh width is shown in Fig. 17. The yellow curve is obtained by the active flux with the exact edge update using time steps that satisfy CFL ≤ 0.9, where the CFL number is given by A larger error, but still 3rd order accuracy, is observed, if we only include those terms of the Taylor series expansion that are shown in Eq. (40). We call the resulting method ADER reduced. The ADER method which uses all the non-zero derivative terms in the Taylor series expansion will be called ADER full. However, it turns out that ADER reduced is less stable. Therefore, we used smaller time steps corresponding to CFL ≤ 0.45. For comparison, we also show the error curve of the active flux method with exact evolution for the same reduced time step size.

Acoustics
The two-dimensional linear acoustic equations have the form where c 0 is the speed of sound, p is the pressure and u, v are the velocity components in the x and y direction, respectively. Let q = ( p, u, v) T denote the vector of the conserved quantities and q 0 = ( p 0 , u 0 , v 0 ) T the initial condition. Then, we can write this linear system in the form Following our approach, the Taylor series expansion for the acoustic equations is given by q(x, y, t n + τ ) = q(x, y, t n ) + τ (−Aq x − Bq y ) The procedure of solving the two-dimensional Riemann problems is almost the same as it is for the advection equation. The only difference is the treatment of the mixed derivative terms. As A and B do not commute, whenever we have to solve a one-dimensional Riemann problem in x-direction (using matrix A) for the term ABq yx , we first multiply the initial states of the Riemann problem with B, so they are now (Bq yx ) L and (Bq yx ) R . The solution to this Riemann problem is then plugged back into (45) and only multiplied with A, and not first with B and then with A, as would have been the case if the Riemann problem was solved for (q yx ) L and (q yx ) R directly. The case B Aq xy is treated in analogy. This concept can be extended to higher derivative terms. For the acoustics problem, the exact two-dimensional solution for data with constant vorticity can be derived by the spherical means formula as was done for example in [4]. This solution formula has the form where R = c 0 t, and This is the formula that was proposed to be used to update the interface values in the active flux method in [4]. Instead of q 0 the reconstruction u at the current time is used and t is replaced by t. The outer integral in (47) is split up into the integrals over the individual cells-two for the edges and four for the corners. This demands that time steps are bounded, such that the circle with radius c 0 t around the midpoint of a grid interface is inside the two neighbouring grid cells. Thus we have a time step restriction of the form CFL ≤ 1, where the CFL number is defined by .
Note that Barsukow et al. [1] have simultaneously to us developed an active flux method for the acoustic equations on Cartesian grids using the spherical mean formula for general initial conditions. In Sect. 6.2, we will compare the active flux method which uses the spherical mean formula with our ADER approach. Upon inspection, the solution (46) is a polynomial in t, as is the Taylor series expansion that we use for our approach. Hence, it is interesting to compare the terms that are obtained by the exact evolution formula with the terms obtained by our approximation. We could not find an approach that was based on the solution of one-dimensional Riemann problems, which matched all the terms that are approximated by the exact evolution formula. Here we do not want to present these long and tedious formulas, which show the differences between the two approaches. Instead we will present some numerical accuracy studies.

Accuracy Study
We consider the acoustic equations (44) with c 0 = 1 and use periodic initial data of the form The solution at time T = 0.2 was computed on the rectangular computational domain 1]. This test problem was suggested by Lukáčová et al. [9], where the exact solution can also be found. We compare the active flux method with exact evolution with our version of the active flux method that uses the ADER approach to update the interface values. All the other components of the method were handled in the same manner.
In Fig. 18 we show the error as a function of the mesh width. Note that for this example the error in u is always equal to the error in v. We again consider two different versions of the ADER approach, one that uses all the non-zero derivative terms in the Taylor series expansion (ADER full) and one version, which uses only those terms of the Taylor series expansion that are needed to get third order accuracy (ADER reduced). In analogy to the advection case, ADER reduced has a reduced stability. Therefore, we used time steps which satisfy CFL ≤ 0.45. The error in the velocity components is almost identical for the three different methods (i.e., exact evolution, ADER full and ADER reduced). The exact evolution produced  Fig. 18 Accuracy study for the two-dimensional acoustic equations. We compare results of our ADER approach with results where the exact evolution formula was used to update the interface values the smallest error, followed by ADER full and ADER reduced. For this test problem, the two ADER versions produced a slightly smaller error in the pressure component than the method which used the exact evolution formula.

Radial Symmetric Test Problem
Now we study a radial symmetric problem, i.e. we solve the acoustic equations with c 0 = 1 and initial values of the form We approximate the problem on the domain [−2, 2] × [−2, 2], set x 0 = y 0 = 0 and use μ = 50. Note that the same problem was studied by Roe and Eymann [4]. There, the active flux method was implemented on a triangular grid and the spherical mean formula was used in order to compute edge values of the conserved quantities during each time step. Again we compare two versions of the active flux which only vary in the update formula of the interface values. We either use the spherical mean formula or the ADER reduced approach. In Figs. 19 and 20 , we show numerical results at time t = 1.25 on grids with 50 × 50 as well as 150 × 150 grid cells. These are scatter plots of the solution, i.e. we plot the interface values over r = (x − x 0 ) 2 + (y − y 0 ) 2 . A scattering of the data indicates grid effects. Both, on the coarse as well as on the fine grid, the two methods produce results that are comparable in accuracy. Note that the accuracy observed on the Cartesian grid compares well with the accuracy documented in [4], where the active flux method was used on triangular grids with comparable resolution.

Conclusion
We have studied the active flux method of Eymann and Roe and proposed several variations of this method for one and two-dimensional hyperbolic problems. The active flux method is very attractive, due to its local stencil as well as its high accuracy, which was documented for several test problems. For one-dimensional linear hyperbolic problems, we showed that the active flux method can be interpreted as a new version of the ADER method. We used this interpretation to propose new active flux methods for nonlinear problems which do not rely on a splitting approach.
While discontinuous reconstructions appear in almost all numerical methods for hyperbolic problems, Roe [12] raised the question as to whether this is really a good idea. Our test computations show that in most situations a continuous reconstruction leads to very accurate results. However, for nonlinear problems we identified a situation where the active flux method in its original version fails to approximate the correct solution structure and becomes unstable. For this situation we propose a discontinuous reconstruction which cures this failure by solving an additional Riemann problem.
Furthermore, we discussed several limiting strategies for the one-dimensional active flux method. For the original active flux method, we proposed two continuous reconstructions,  which led to accurate approximations of linear problems. Unfortunately, these new reconstructions are not appropriate for the ADER implementation of the active flux method. For the ADER version, limiting was introduced via the use of discontinuous, piecewise quadratic reconstructions. Our limiting does not increase the stencil of the active flux method.
Alternatively, we can use a discontinuous, piecewise quadratic reconstruction that is bound preserving as described by Zhang and Shu [22]. Such a limiting approach can be extended to multi-dimensional problems and can be incorporated into the original version of the active flux method as well as into our ADER approach. We are currently investigating this approach for one and two-dimensional nonlinear problems.
In the two-dimensional case, there are many possibilities to define ADER methods which use the same degrees of freedom as the original active flux method. We constructed ADER methods, which use the underlying continuous, piecewise quadratic reconstruction as much as possible in order to reduce the number of Riemann solves. For the advection equation, the resulting ADER method is equivalent to the active flux method that uses the exact evolution formula for the interface update. For the two-dimensional acoustic equations, the exact evolution of interface values could not be obtained by an ADER type approach. However, our accuracy studies show, that we can still obtain very accurate results.
Roe's motivation to introduce the active flux method was the construction of truly multidimensional methods which do not rely on the solution of one-dimensional Riemann problems.
If an exact evolution formula exists, very efficient and accurate active flux methods can be derived. However, for most hyperbolic problems exact evolution formulas are not available. Therefore, we believe that a good understanding of the connection between the active flux method and the ADER method can lead to new third order accurate methods for conservation laws.
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.
In regions with continuous reconstruction, we compute f (Q n i+ 1 2 ), i.e. we evaluate the tensor at the interface value of the conserved quantities, which were computed in the previous time step. At interfaces with a discontinuous reconstruction, we first recompute Q n i+ 1 2 by solving an additional Riemann problem as described in Sect. 5.2. Now the tensor f is evaluated at this recomputed interface value.