AENO: a Novel Reconstruction Method in Conjunction with ADER Schemes for Hyperbolic Equations

In this paper, we present a novel spatial reconstruction scheme, called AENO, that results from a special averaging of the ENO polynomial and its closest neighbour, while retaining the stencil direction decided by the ENO choice. A variant of the scheme, called m-AENO, results from averaging the modified ENO (m-ENO) polynomial and its closest neighbour. The concept is thoroughly assessed for the one-dimensional linear advection equation and for a one-dimensional non-linear hyperbolic system, in conjunction with the fully discrete, high-order ADER approach implemented up to fifth order of accuracy in both space and time. The results, as compared to the conventional ENO, m-ENO and WENO schemes, are very encouraging. Surprisingly, our results show that the L 1 -errors of the novel AENO approach are the smallest for most cases considered. Crucially, for a chosen error size, AENO turns out to be the most efficient method of all five methods tested.


Introduction
Circumventing Godunov's theorem [8] is a 60-year old challenge to numerical mathematicians concerned with the development of high-order methods for solving hyperbolic equations. The invention of total variation diminishing (TVD) methods in the 1970s and 1980s 1 3 [11,16,23,29,43] represented an important step forward in facing the challenge. However such methods, even if very effective in practice, are limited to, at most, second-order accuracy for smooth solutions. Currently, the alternative approaches relying on non-linear spatial reconstruction pioneered by Harten and others [12,13], such as ENO, m-ENO and WENO [14,15,[26][27][28] have become the dominant approach. Unlike TVD methods, however, ENO, m-ENO and WENO based methods cannot guarantee the absence of spurious oscillations for scalar problems. However, these methods guarantee that such oscillations vanish in the limit of mesh refinement. This is not the case for linear methods, for which, no matter how fine the mesh is, the amplitude of the oscillations will remain unchanged. The initially popular ENO method has, in recent years, been largerly discarded by most researchers and practitioners. The reconstruction scenery is currently dominated by the WENO approach, in a variety of forms. See for example [5,31] and [6].
In this paper, we propose a novel spatial reconstruction method that is akin to both ENO and WENO. The method, called AENO, results from averaging two polynomials, the classical ENO polynomial and its closest neighbour, while the search for the stencil remains commanded by ENO. A variant of the scheme, called m-AENO, results from averaging the modified ENO (m-ENO) polynomial of Shu [26] and its closest neighbour; m-AENO is equivalent to AENO for 2nd and 3rd order of accuracy. Here, the reconstruction schemes are applied in conjunction with the fully discrete high-order ADER approach [4,7,30,34,40]. See [36] for an up-to-date review of ADER and references to the many contributions to its development. Schemes of up to 5th order of accuracy in space and time are implemented and tested, first for the linear advection equation and then for a non-linear hyperbolic system, namely the blood flow equations. For both problems, we first carry out experiments to compare numerical solutions with exact solutions for short and long evolution times. Results for three reconstruction methods are compared, namely ENO, m-ENO, WENO and the novel AENO and m-AENO schemes. Then we carry out a convergence rate study for both types of problems, the linear advection equation and the blood flow equations. Overall, the results are encouraging. For most problems, all five reconstruction methods give comparable results. Surprisingly, our results show that the L 1 -errors of the novel AENO/m-AENO schemes are the smallest, for most cases considered.
The rest of this paper is structured as follows. In Sect. 2 we review the fully discrete ADER methodology to solve hyperbolic equations, first for the linear scalar case and then for non-linear systems. In Sect. 3 we briefly introduce necessary background on the ENO reconstruction method. In Sect. 4 we describe the new spatial reconstruction technique AENO and its variant m-AENO; in Sect. 5 we show results for the linear advection equation. In Sect. 6 we apply the methods to the blood flow equations and in Sect. 7 we draw conclusions.

Review of the ADER Methodology for Hyperbolic Equations
Here, we briefly review the ADER approach to construct one-step, fully discrete high-order Godunov methods for solving hyperbolic balance laws, whose building block is the generalized Riemann problem, a piece-wise smooth data Cauchy problem including stiff or non-stiff source terms, rather than the conventional piece-wise constant data homogeneous Riemann problem. In this section we review the ADER methodology as applied to onedimensional systems of hyperbolic balance laws. We first deal with the scalar linear case.

ADER for the Scalar Linear Case
Before proceeding, we recall the finite volume framework for a general scalar balance law.

The Finite Volume Framework
Consider the one-dimensional balance law solved by a finite volume method, with space-time volumes Exact integration of (1) on V in (2) gives with definitions and

ADER for Linear Advection-Reaction
The ADER method [38] departs from the finite volume formula (3), understood as resulting from approximations to the integrals in (4). The one-step updating formula (3) updates the cell average q n i to q n+1 i , for which definitions for the numerical flux f i+ 1 2 and the numerical source s i are required. The rest of the presentation of the ADER approach will be carried out with (1) as the linear advection-reaction equation, in which f (q) = q and s(q) = q , with a constant wave propagation speed and constant reaction rate satisfying ⩽ 0.

The ADER Flux and the Generalized Riemann Problem
To determine the numerical flux f i+ 1 2 in (3) we compute a high-order accurate approximation to the flux integral in (4), which in turn requires a high-order approximation to the integrand q(x i+ 1 2 , t) . This is achieved by solving, approximately, the following generalized Riemann problem (1) t q + x f (q) = s(q), q(x, t n )dx, , t))dt, ).
in which the polynomials p i (x) and p i+1 (x) are reconstruction polynomials of arbitrary degree K. To solve (6) we adopt here the Toro-Titarev [41] method, which starts from the LeFloch-Raviart expansion [18] for the approximate solution The leading term of the expansion. The leading term q(0, 0 + ) in (7) results from solving the conventional piece-wise-constant data Riemann problem in which the initial condition consists of the boundary extrapolated values p i (0) from the left and p i+1 (0) from the right. The position of the interface x i+ 1 2 corresponds to 0, in local coordinates.
The higher-order terms of the expansion. To compute the higher-order terms in (7) one requires the determination of the time derivatives (k) t q(0, 0 + ) . In the first version of ADER [38] it was proposed to apply the Cauchy-Kovalevskaya procedure to the governing equation in (6). Just for a moment let us assume the linear homogenous case in (6) ( = 0 ). Then, the Cauchy-Kovalevskaya procedure gives In fact, it is easily shown that for an arbitrary positive integer k, Similarly, for the inhomogeneous case with non-zero source term in (6), it can also be shown that the time derivatives can still be expressed as functions of spacial derivatives, as follows: Now the pending problem is that of determining (l) x q(0, 0 + ) in (7), for which a useful step is the construction of evolution equations for the spatial derivatives. For the homogenous case, it can be easily shown that (l) x q(x, t) obeys the same evolution equation as q(x, t), that is (6) PDE: t q + x q = q, PDE: t q + x q = 0, (9) t q + x q = 0, t q = − x q, (2) t q = 2 (2) x q.

3
The next problem is to explore the possibility of formulating initial value problems for (12). In fact this is possible if we find suitable initial conditions. If for each computational cell we have a spatial polynomial representation p i (x) for q(x, t n ) , via a spatial reconstruction technique for example, then it is possible to pose classical Riemann problems for spatial derivatives as follows: The initial condition for these classical Riemann problems (piece-wise constant data) consists of intercell-boundary extrapolated spatial derivatives of the reconstruction polynomials, p i (x) from the left and p i+1 (x) from the right. Solutions of classical Riemann problems for the spatial derivatives, as in (13), will give the spatial derivatives (l) x q(0, 0 + ) in (7), of any order l. This in turn determines the time partial derivatives in (11) and hence in (7). We have therefore found the complete solution of the homogeneous version of the generalized Riemann problem (6), which is given as The solution for the inhomogeneous case is similar. The numerical flux results from an approximation to the flux integral in (4), which is now obtained as

Exact integration gives
Generally, the numerical flux in (15) will require numerical integration of the appropriate order. Inserting the numerical flux (16) into the conservative finite volume formula (3) will give a numerical approximation of order K + 1 in both space and time. In the presence of source terms in the governing equations, one must include the source terms when performing the Cauchy-Kovalevskaya procedure, as these enter in the calculation of the numerical flux. In addition, one must compute the numerical source to the appropriate order of accuracy. This requires an approximation of the source integral in (4) in the entire control volume V in (2), for which a high-order approximation to q(x, t) within V is needed. Such approximation is found in a manner similar to that of the generalized Riemann problem.

ADER for Non-linear Systems
Consider a general one-dimensional system of hyperbolic balance laws The vectors , ( ) and ( ) define conserved variables, fluxes and sources, respectively. Assuming the space/time domain to be discretised by finite volumes V in (2), exact integration of (17)

ADER and the Godunov Method
The classical Godunov method [8] results from (18)- (19) by assuming ( ) = and (x∕t) to be the similarity solution of the piece-wise constant data, homogeneous, Riemann problem Figure 1, top frame, depicts the initial condition (left) and the structure of the solution (right) of (20). In this simplified case the numerical flux i+ 1 2 in (18) results from the exact integration in (19), namely God For exact and approximate Riemann solvers for the first-order Godunov method see [35]. The ADER extension of Godunov's method [38] includes three steps.
i) Replacing piece-wise constant data by piece-wise smooth data, as resulting from reconstruction procedures, for example; see Fig. 1, bottom frame. This leads naturally to the generalized Riemann problem ii) Computing LR ( ) , the time-dependent solution of (22) at the interface, followed by evaluation of the numerical flux iii) Defining a function i (x, t) within each volume V, followed by evaluation of the numerical source (22) PDEs: Early communications on the ADER methodology include [38] and [25]. For an introduction to ADER methods see Chaps. 19 and 20 of [35]. The next section deals with solvers for the generalised Riemann problem (22).

Solvers for the Generalized Riemann Problem
Several methods for solving the generalised Riemann problem are currently available. Here we limit ourselves to brief descriptions for two of them and give appropriate references for other solvers.
The Toro-Titarev solver for the generalized Riemann problem This solver [41] seeks an approximate solution LR ( ) of (22) expressed as a truncated series expansion following [19]. The leading term is defined as and is computed by solving a conventional Riemann problem, see (8). The challenge is to determine the coefficients (k) t (0, 0 + ) in (25) for the higher-order terms. A precursor to the Toro-Titarev solver [41] results from a modification of the Ben-Artzi/ Falcovitz second-order GRP solver [1] that departs from LR ( ) = (0, 0 + ) + t (0, 0 + ) , with (0, 0 + ) defined as in (26). See Chap. 14 of the 1997 edition of [35]. The Cauchy-Kovalevskaya procedure is then used to determine t (0, 0 + ) , that is and the solution then reads The pending problem in (28) is to compute the spatial derivative x (0, 0 + ) . First it is noted that for the linear homogeneous case the following evolution equations for all spatial derivatives are valid: One can then define the classical Riemann problem for (30), with k = 1 , to find x (0, 0 + ) in (28). The numerical flux (23) then follows.
The Toro-Titarev solver for the GRP k [35] departs from the LeFloch-Raviart expansion (25) and the Cauchy-Kovalevskaya procedure, leading to The functionals (k) are specific to the particular system (17); their arguments are spatial derivatives, yet to be found. To determine these, evolution equations are invoked Then, simplified classical Riemann problems for spatial derivatives are posed Strong simplifications have been made to the evolution equations in order to arrive at (33); the source terms (k) have been neglected and the advection term has been linearised around the leading term of the expansion (25). The coefficient matrix (0) LR results from evaluating the Jacobian at the leading term of the expansion. The complete solution (25) has then been determined and the numerical flux follows from (23). A similar but simpler procedure is used to determine the source term in (24), see [41].

The Montecinos-Toro implicit GRP solver
This solver [21,40] is implicit and can deal with stiff source terms. The key step is the following lemma.

Lemma 1 Let (x, ) be an analytical function of . Then
This expression is the implicit counterpart of the explicit Raviart-LeFloch expansion (25). Now the implicit Cauchy-Kovalevskaya procedure leads to As for the Toro-Titarev solver, a key step is to determine the arguments of functionals (k) in (35). For details see [40] for the general case. Example 1 Second-order scheme. In this special case, see [21] for details, we have The implicit Cauchy-Kovalevskaya procedure gives Then, applying the implicit expansion (34) we have and Therefore Finally, the solution (0, ) = * ( ) satisfies the non-linear algebraic system x (x, t)) = , ICs:      x (0, 0 + ). and the numerical flux follows. This gives a second-order, locally implicit method suitable for balance laws with stiff source terms. Other k solvers We have described two schemes to solve the generalized Riemann problem. The first solver is explicit, while the second one is implicit, leading to a locally implicit ADER method. The implicit method can deal with stiff source terms and the corresponding ADER methods are then able to reconcile stiffness and high accuracy, for smooth solutions. These two methods seek the solution right at the interface and extend the (explicit) Ben-Artzi/Falcovitz second-order method. Another class of solvers emerges from a re-interpretation of the high-order method of Harten et al. [13] that extends the MUSCL-Hancock second-order method. In this approach, the initial data is evolved in time and interactions at the interface are sought at selected time-integration points via classical Riemann problems [2]. The already well-established method of Dumbser et al. [4] is a numerical version of the Harten method, see also [2], whereby the data is evolved numerically via an implicit space-time discontinuous Galerkin approach. The resulting ADER scheme is suitable for stiff source terms, with the additional benefit of avoiding the Cauchy-Kovalevskaya procedure. Other recent solvers for the generalized Riemann problem are due to Götz and Iske [10], Götz and Dumbser [9], and Dematté et al. [3].
We have reviewed the ADER approach to extend Godunov's method to high order of accuracy. Very high-order methods are essential for the following reasons: (i) there is a growing trend to use PDEs to understand the physics they embody; (ii) only very accurate solutions of PDEs will achieve this and also reveal limitations of mathematical models (PDEs themselves) and uncertainty on parameters of the problem; (iii) efficiency: given an error deemed acceptable, then high-order methods attain that error, if small, much more efficiently than low-order methods on fine meshes, by orders of magnitude. For very long time evolution simulations high-order methods are essential.
In the presentation of ADER we assumed that the solution at a time t n in any computing cell was represented by a polynomial resulting from some spatial reconstruction procedure. In principle, it can be any kind of reconstruction procedure, centred, biased, linear, non-linear. In what follows we review a non-linear reconstruction procedure.

Review of ENO Polynomial Reconstruction Methodology
To circumvent Godunov's theorem when constructing high-order methods, it is necessary to resort non-linear schemes and hence avoid or reduce spurious oscillations in the vicinity of large gradients of the solution. A popular method to do so is the essentially non-oscillatory (ENO) reconstruction method of Harten and collaborators [12], which we briefly review here. The presentation is succinct and we follow [17] to construct the polynomial with the required properties. More popular than ENO is the WENO method, which we also use in this paper for comparison, but a review of WENO is omitted as the reader can consult the classical references [15,27,28]; see also [31]. There is a variant of ENO due to Shu [26], called here m-ENO. We omit all details of m-ENO here but encourage the reader to consult the original reference.

ENO Polynomial of Arbitrary Degree
We assume a set of cell integral averages of the unknown q(x, t) denoted as {q j } i j=0 and adopt the primitive-function approach [17,45] to construct an ENO polynomial with samples q(x i+1∕2 ) at the intercell boundaries x i+1∕2 . Recall that a function q(x) is the primitive function of f(x), or antiderivative, if where x −1∕2 is the left most interface.Then, the samples are defined as That is, the partial sums of cell averages are the samples of the primitive function q(x). Now the task is to construct the interpolating polynomial P(x) of degree N + 1 passing through (x i+1∕2 , q(x i+1∕2 )) , of one degree higher than the sought reconstruction polynomial p(x), with For convenience, the final sought polynomial p(x) of degree N is written in Taylor form as follows: The coefficients a 0 , a 1 , ⋯ , a N depend on the stencil of the polynomial. Next, we find the Newton divided differences. Recalling from (43) point added, to left or to the right, depending on the relative size of divided differences in absolute value. The general algorithm is as follows. Start with the two points l 0 (i) = i + 1∕2 at the cell interface x i+ 1 2 and l 1 (i) = i − 1 2 at the cell interface x i− 1 2 . The most left point is l 1 (i) = i − 1 2 . Then the next point in the sequence is l 2 (i) . In general the new point to be added is chosen recursively as follows.
For m = 1, ⋯ , N, The coefficients a j of the polynomial (45) are as follows: with Therefore the sought polynomial (45) has been determined. For details see for example [17].
In what follows we give details for ENO polynomials of second to fourth degree.

Piecewise-Linear ENO Reconstruction
The reconstruction polynomial of degree 1 on the cell ] is or on the stencils

and the coefficients are
As stated in Sect. 2, the ADER method requires the value of the polynomial, as well as its spatial derivatives, at the cell interfaces. The sought values of the polynomial at interface x i+1∕2 are

Piecewise-Quadratic ENO Reconstruction
The reconstruction on the cell

Piecewise-Cubic ENO Reconstruction
The reconstruction on the cell The stencils are The ENO choice of indexes is as follows: The coefficients read The polynomial values at the interface x i+ 1 2 are (62)

Piecewise-Quartic ENO Reconstruction
The reconstruction on the cell The stencils are The ENO choice of indexes is To determine the coefficients we first set (67) Then the coefficients are The polynomial values at interface x i+ 1 2 are (71) As already pointed out, the evaluation of the polynomial and its derivatives at the cell interface x i+ 1 2 , from the left and the right, is needed for solving the generalized Riemann problem to determine the numerical flux. To determine the numerical source we must evaluate the volume integral (24), for which one makes use of the polynomial p(x) and its derivatives in cell i, at nodes x l of the integration formula used, e.g., a Gauss-Legendre quadrature.
We have reviewed the ENO method as implemented in the fully discrete high-order ADER methodology. The ENO reconstruction method has been criticised for a number of reasons. One of them is the abrupt change of stencil simply due to small changes in the divided differences that determine the coefficients of the ENO polynomial. In fact, even round-off errors may decide the choice of stencil. Shu [26] addressed this issue by proposing a modified ENO method, called m-ENO here. In the next section, we propose a simple averaged ENO-type method, in which the classical ENO polynomial is averaged with its closest neighbour. A variant m-AENO results from substituting ENO by m-ENO, for schemes of 4th and greater order of accuracy. The results are encouraging.

Novel Polynomial Reconstruction: AENO
In this section, we present a non-linear polynomial reconstruction procedure that is akin to both the existing ENO and WENO procedures [12-15, 27, 28]. We first present an example that motivates this section.

A Motivating Example: Second-Order ADER Method
We begin by considering the simplest, first-degree, polynomial reconstruction ) is the cell centre and i is the slope of p i (x) , still to be determined. Figure 2 depicts two polynomials of the form (74) and the results are respectively determined by choosing the slope thus (73) Assume that the slope i is chosen as the weighted average It is easy to show that the second-order ADER method with the linear reconstruction (74) reproduces the classical second-order schemes of Lax-Wendroff, Warming-Beam and Fromm if in (76) is chosen as follows: The classical ENO [12] second-order method chooses In other words, the classical ENO second-order method switches non-linearly between the Lax-Wendroff and the Warming-Beam methods, depending on the relative sizes in the absolute value of the slopes i− 1 2 and i+ 1 2 . A simple third-order method. It is interesting to note that, still in the frame of the second-order ADER method, if the weight parameter is chosen as then, a third-order accurate method results, in both space and time. Here c is the Courant number. This scheme was first presented in [33], see Eq. (13.39), and also in [37], ] . These result from choosing the slope as in (75) in the setting of the MUSCL-Hancock method [44]. However, it has been shown that for the homogeneous case, the MUSCL-Hancock method is entirely equivalent to the ADER method in second-order mode (ADER2), see [20]. We have again verified this in the setting of ADER methods. The details of the analysis are omitted. Figure 3 displays computed results for the linear advection equation, for which a Gaussian profile is used as initial condition, where comparison is made with ADER2 in conjunction with ENO reconstruction (74), (75) and (78). Table 1 shows empirical convergence rates for the linear advection equation in which it is verified that the ADER2 (formally second-order) scheme with reconstruction (74) with (76) and given by (79) is third-order accurate in both space and time.
The AENO reconstruction method presented in this paper, in fact emerges from averaging of two polynomials in the frame of the ENO scheme to choose the reconstruction stencil. The construction of ADER schemes of up to fifth order in space and time, based on averages of the type (76) were constructed by MSc student Andrea Santacá in his Master thesis [24]. Additionally, in the present paper we include a variant of AENO, called m-AENO, to be explained later.  Table 1 Empirical convergence rates for ADER2 with the particular choice of slope

A Joining Function for Averaging Two Polynomials
Motivated by the example of the previous Sect. 4.1 we seek an averaging procedure that can be applied to two polynomials of arbitrary degree. First we define the function and ∈ R , a constant, with | | < 1 . Consider now the averaging function J( (x)) joining smoothly two constant, real states q L and q R , as follows: Note the analogy between the weighted average (76) and (81). While in (76) the weight is constant, in (81) (x) is obviously variable. Figure 4 shows the weight function (x) for 2 = 0.2 . For x < 0 the function (x) approaches 1 rapidly, while for x > 0 , (x) approaches −1 . Note the following properties of (x): Note also the associated properties of J( (x)):  Figure 5 shows the distribution of J( (x)) for q L = 3 , q R = −2 and 2 = 0.2 . It is seen that That is the constant states q L and q R are approached asymptotically by J(x). However we note that these states are approached very rapidly. Consider a positive integer N, then For large N, That is, at a distance of 10 from the origin, J(x) attains the constant states with an error of less than 1% . Finally, we note that the use of 2 in (87) (N ) ≈ −sign( ) = −1.  the definition of (x) is convenient, as it permits a precise evaluation of J(x) at a distance from the origin. In particular applications, a translation might be necessary. In this paper we use J(x) for averaging two polynomials, through the averaging of their coefficients.
When applying the above averaging function (81) to two states q L and q R we proceed as follows. First we define where TOL is a small positive quantity chosen to avoid division by zero, e.g., TOL = 10 −6 . Then we redefine It is easy to see that for the joining function (81) we have In other words, in the joining function (81) acting on two states q L and q R , the smaller argument in absolute value receives the larger weight and the larger argument in absolute value receives the smaller weight. This averaging procedure will be applied to two polynomials, the ENO polynomial and that closest to it, by applying it to each pair of respective coefficients. The coefficients of the ENO polynomial will receive the larger weight, as we shall see in the next section.

AENO: a Two-Polynomial Average Method
In Sect. 3 we reviewed the ENO spatial reconstruction method. In Sect. 4.1 we showed that averaging of the ENO polynomial with its closest neighbour is capable of producing encouraging results. In particular in Sect. 4.1, we showed that for the formally secondorder ADER method there is an averaging for which such formally second-order ADER method is actually third-order accurate in both space and time.
In this section, we propose a reconstruction method based on a two-polynomial averaging procedure. Two schemes are proposed. Scheme 1, called AENO, results from a weighted average of the classical ENO polynomial and its closest neighbour. Scheme 2, called m-AENO, follows a similar approach and results from a weighted average of the classical m-ENO polynomial of Shu [26] and its closest neighbour. We first describe AENO.
Here we introduce a general averaging procedure as applied to two neighbouring polynomials of arbitrary degree. Assume an Nth degree polynomial that has been obtained by differentiating a polynomial P(x) of degree N + 1 . By redefining the coefficients as we write For the sake of simplicity, in what follows the hat is omitted.
Recall that with a reconstruction polynomial of degree N it is possible to construct an ADER numerical method of order (N + 1) in both space and time. Then in the context of the ENO method, we need N + 1 stencils and following the ENO approach we find a single, best polynomial, the ENO polynomial. In the AENO approach, instead of choosing just the single ENO polynomial we consider another polynomial, the closest to the ENO polynomial and take an average of these two. Moreover, the search of the stencils for the two polynomials is commanded by the conventional ENO approach and the coefficients of the ENO polynomial take the largest weight.
Example 2 In the case of a third-order method we just need second-degree polynomials, for which we have three possible candidates, namely In the ENO context we choose one of these polynomials according to the ENO algorithm presented earlier. In this example, we need two stages. In the AENO case we perform just one step of the algorithm; basically, we choose the direction based on the first-order slope, which is either left or right. If the direction is left then we take an average of the coefficients of p L (x) and p M (x) ; otherwise we take an average of the coefficients of p M (x) and p R (x) . Suppose the direction is the left one, then the sought polynomial is with and (91) The coefficient ã 0 is chosen so as to comply with the conservation property, that is The parameter 2 is a positive constant, while TOL is a small positive tolerance to avoid division by zero. In practice we take TOL = 10 −6 .
Generalization: The AENO procedure is then generalized as follows. The m-AENO scheme follows a similar approach to the AENO scheme presented above. The difference is that, instead of relying on the classical ENO polynomial we rely on the classical m-ENO polynomial of Shu [26]. Then the m-AENO scheme results from a weighted averaged of m-ENO polynomial and its closest neighbour. We omit the details. For background of the m-ENO polynomial the reader is referred to the original paper of Shu [26]. (97) In the next section, we carry out a thorough assessment of the newly proposed spatial reconstruction procedure, as applied to the linear advection equation.

Results for the Linear Advection Equation
Here we assess the performance of ADER with the newly proposed AENO reconstruction procedure. For the assessment we compare numerical solutions against the exact solutions and against numerical solutions with ENO, m-ENO and WENO, also in conjunction with ADER. To this end we consider four test problems, namely IVP1, IVP2, IVP3 and IVP4 given below The exact solution in each IVP is given as

Solution Profiles for Long-Time Evolution
The purpose of this section is to illustrate the performance of the new schemes for linear IVPs for long-time evolution, which is known to be particularly challenging for low accuracy schemes. We compare results to exact solutions and to numerical solutions obtained with established reconstruction methods, namely ENO, m-ENO and WENO. To this end, we use IVP1 with smooth initial conditions and IVP4 with a combination of smooth and discontinuous parts. The second problem IVP4 is a more realistic representation of (104) practical problems for hyperbolic equations that involve both smooth parts as well as discontinuities. Figures 6,7,8,9 and 10 show computed results (symbols) for IVP1 (104) compared to the exact solution (line) for the ADER high-order scheme with various reconstruction methods, from 2nd to 5th orders of accuracy in both space and time. For the new methods AENO and m-AENO we used 2 = 1 2 . Obviously, for all methods, increasing the order of accuracy progressively improves the agreement between the numerical solution and the exact solution. The reconstruction method used has a visible effect on the computed solution; this is patently obvious for 2nd and 3rd orders of accuracy but it is still visible for higher orders of accuracy. Long-time evolution also contributes to differentiate between methods in terms of accuracy. Judging from the figures, AENO is clearly more accurate than both ENO and m-ENO; this is especially clear for 2nd, 3rd and 4th orders of accuracy. In fact, the new AENO schemes of this paper compare well with the sophisticated WENO; it is even reasonable to state that the new AENO schemes have a small advantage over WENO; compare carefully the profiles for 3rd and 4th orders of accuracy. For second order, WENO is superior to all other schemes. Figures 11,12,13,14 and 15 show computed results for the multiple wave test IVP4 for a long evolution time T out = 2 000 units and a coarse mesh of just M = 100 cells. This is an exceedingly demanding test problem, as it contains the conflicting requirements of high accuracy for smooth parts and high resolution without spurious oscillations in the vicinity of discontinuities. Figure 11 shows results for ENO in conjunction with ADER for schemes of 2nd to 5th order in space and time. Clearly, for the long chosen output time and the coarse mesh used, even the 5th order scheme gives large errors; even higher order of accuracy would be required to obtain acceptable results. The challenges of this test are representative of practical computational problems involving long time evolution, such as in acoustics, seimic waves and tsunamic waves. Figure 12 shows results for m-ENO in conjunction with ADER for schemes of 2nd to 5th order in space and time. Results are comparable to those of ENO in Fig. 11. Figure 13 shows results for the new AENO scheme. Compared to the classical ENO and m-ENO shown in Figs. 11 and 12, the new AENO scheme is clearly superior; this is more evident for 3rd and 4th orders of accuracy. Surprisingly, this observation is also true for WENO, comparing Figs. 13 and 15. Only for the 5th order case, AENO and WENO are, roughly, comparable; WENO displays some more visible undershoots and a spurious oscillation on the right-hand side. Results for the second version of our averaged ENO, m-AENO, are shown in Fig. 14. For 2nd and 3rd order accuracy, these results are identical to those of AENO. In fact, the schemes ENO and m-AENO are identical in these two cases. For 4th and 5th order cases the performance of m-AENO is visibly inferior to that of AENO. In fact, m-AENO is also inferior to the classical ENO, m-ENO and WENO schemes for this test problem and for the 4th and 5th order cases. Table 2 shows the errors in the computed solutions for the multi-wave test problem.
For the multi-wave test IVP4, it is the discontinuous components of the profile that pose the most severe challenge to the numerical methods. This is a recurrent feature of high-order methods for hyperbolic equations. To highlight this feature, we carried out computations by keeping the square wave only, in the initial condition of IVP4. Results are shown in Figs. 16 and 17. Figure 16 shows four frames corresponding to ADER schemes of orders 2, 3, 4 and 5. For each order we plot results for all reconstruction  schemes, noting that for orders 4 and 5 we have five distinct reconstruction schemes. For second-order ADER (ADER-2), WENO gives the best results, whereas for all higher-order ADER schemes, AENO gives the best results, by an appreciable margin, most evident for 3rd and 4th orders of accuracy. Note the undershoots in the WENO results for the 5th order scheme. Figure 17 shows four frames corresponding to four reconstruction schemes, including the two new ones of this paper, namely AENO and m-AENO. For each reconstruction schemes we see the effect of increasing the order of accuracy, from 2 to 5. Note the peculiar behaviour of m-AENO for the 4th order scheme. At this stage we note, however, that when it comes to smooth solutions, even with large derivatives, the modified AENO scheme, namely m-AENO based on m-ENO, performs very well indeed, as we shall see in the convergence rates results in Sect. 5.2.

Convergence Rate Study for the Linear Advection Equation
In this section, we carry out a convergence rate study for the linear advection equation through IVP2 (105) and IVP3 (106), for the schemes of 2nd to 5th order of accuracy in space and time. We compare the newly presented schemes AENO and m-AENO, with the established ENO, m-ENO and WENO reconstruction methods, all of them used with the fully discrete ADER approach. IVP2 (105) could be described as a smooth test of the kind commonly used in the literature to assess convergence rates. IVP3 (106) also serves the same purpose but it is recognised as an exceedingly severe test problem, for which many commonly used high-order schemes fail to give the expected rates. Tables 3, 4 Tables 3, 4, 5, 6 and 7. All schemes attain the theoretically expected convergence rates, for all orders, in the L 1 norm. WENO, with the exception of the second-order scheme, attains the theoretically expected convergence rates also in the L 2 and L ∞ norms. The remaining schemes show poor performance in these norms, especially in the L ∞ norm. Generally, ENO and m-ENO show comparable performance, perhaps with small advantage to ENO. The new AENO scheme does not perform satisfactory in the L 2 and L ∞ norms; as a matter of fact in these norms AENO is inferior to ENO and m-ENO. The new m-AENO scheme improves with respect to AENO for the 5th order scheme. Figure 18 shows computed L 1 -errors for IVP2 (105), as the mesh is refined, for the ADER scheme with all reconstruction schemes of this paper: ENO, m-ENO, WENO, AENO and m-AENO. With the exception of the second-order case, for which WENO is the   most accurate reconstruction method, the AENO schemes of this paper outperform all the other reconstruction methods. Results for IVP3. As stated previously IVP3 is an exceedingly severe test problem, for which many of the commonly used high-order schemes fail to give the expected convergence rates. Results for the schemes of this paper are shown in Tables 8, 9, 10, 11 and 12, where errors and corresponding orders of accuracy are shown for the three norms L 1 , L 2 and L ∞ .

Results for IVP2. Results are shown in
For the ENO scheme, results are shown in Table 8. The scheme attains the expected convergence rates sub-optimally for the second and third order cases in the L 1 norm, while failing in the L 2 and L ∞ norms. The 4th order ENO scheme fails to attain the expected rates in all three norms, while for the 5th order case computed rates are close to the expected ones, but are sub-optimal. In summary, ENO fails the convergence rate test for IVP3. For the modified ENO (m-ENO) scheme of Shu [26] we show the computed results in Table 9, Table 3 ENO results for IVP2 (105): convergence-rate study for the ADER scheme with the ENO reconstruction for the linear advection with = 1 , T out = 0.5 and C cfl = 0.9 Mesh where significant improvements with respect to ENO are seen; however, computed rates are sub-optimal and failure is seen in some cases for the L 2 and L ∞ norms. For the new AENO scheme of this paper, results are shown in Table 10. For the lowerorder cases, 2nd and 3rd order, results are satisfactory, but the scheme fails to attain the expected convergence rates for the higher-order schemes, 4th and 5th order. For the new variant of AENO, called m-AENO, results are shown in Table 11. Note that m-AENO is identical to AENO for the second-and third-order schemes, see Table 10 and corresponding comments. The performance of m-AENO for the 4th and 5th order schemes is very satisfactory, attaining the expected rates in the L 1 and L 2 norms, even if sub-optimally for the L 2 norm. In the 4th order case the scheme does not attain the expected rate in the L ∞ norm. In summary, m-AENO constitutes a significant improvement over AENO, for the higher-order range. Table 4 m-ENO results for IVP2 (105): convergence-rate study for the ADER scheme with the m-ENO reconstruction for the linear advection with = 1 , T out = 0.5 and C cfl = 0.9 Mesh Results for the classical WENO scheme in conjunction with the fully discrete ADER approach are shown in Table 10. Generally, the WENO results are satisfactory for this severe test problem. For the second-order case convergence in the L 1 norm is attained, being sub-optimal for the L 2 norm and failing for the L ∞ norm. For the higher-order schemes, convergence is attained in all norms, except for the 4th order scheme, where convergence is sub-optimal in the L 1 norm and fails in the L 2 and L ∞ norms. In summary, out of the five reconstruction schemes, WENO gives the best convergence rate performance for this demanding IVP3 test problem. Figure 19 shows computed L 1 -errors for IVP3 (106), as the mesh is refined, for the ADER scheme with all five reconstructions methods: ENO, m-ENO, WENO, AENO and m-AENO. For the second-order scheme, WENO has the smallest errors, followed by AENO and then by ENO. For the third-order schemes, AENO has the smallest error, followed by WENO, m-ENO and ENO. For the 4th and 5th order schemes, m-AENO has Table 5 AENO results for IVP2 (105): convergence-rate study for the ADER schemes with the AENO reconstruction with 2 = 0.5 for the linear advection for = 1 , T out = 0.5 , C cfl = 0.9 Mesh the smallest error, while AENO's performance deteriorates visibly as the order of accuracy increases. This observation is consistent with the convergence rate study for IVP3 discussed above.
In the next section, we assess the new AENO reconstruction method as applied to a nonlinear hyperbolic system.

Results for the Blood Flow Equations for Arteries
In this section, we assess the methods as applied to a non-linear hyperbolic system. The ADER method is implemented with the Toro-Titarev solver [41] for the generalized Riemann problem, as described in Sect. 2.2.2. All three reconstruction methods are implemented in terms of characteristic variables, for all orders of accuracy. Table 6 m-AENO results for IVP2 (105): convergence-rate study for the ADER schemes with the m-AENO reconstruction with 2 = 0.5 for the linear advection for = 1 , T out = 0.5 , C cfl = 0.9 Mesh 1 3

The Governing Equations
Now we extend and test the new spatial reconstruction AENO method to a non-linear system of hyperbolic equations, namely the blood flow equations for arteries. The system, in general conservation-law form reads The first two equations are statements of conservation of mass and momentum, while the third equation represents the advection of a tracer (x, t) transported with the blood velocity u(x, t). Here A(x, t) is the blood vessel cross-sectional area. In the source term R is resistance, which is neglected here. The system is closed with a tube law for arteries, given as Here A 0 is the equilibrium cross-sectional area, E is the Young modulus of the vessel wall and h 0 is the vessel-wall thickness. These three quantities are parameters of the problem and in general they vary along the length of the vessel. In (110) we have assumed that all parameters of the problem are constant and where is the constant density of the blood. We remark that for blood flow in veins, a different tube law applies [39]. We assess the method in terms of two problems. First, we consider a Riemann problem with exact solution containing discontinuities. The purpose is to assess the performance of the non-linear reconstruction method at discontinuities, expecting absent or much reduced Table 9 m-ENO results for IVP3 (106): convergence-rate study for the ADER scheme with the m-ENO reconstruction for the linear advection with = 1 , T out = 4 and C cfl = 0.9 Mesh spurious oscillations. The second problem has a smooth exact solution and is used to evaluate the accuracy, expecting to attain the desired convergence rates.

Riemann Problem Solution Profiles
Here we assess the methods as applied to a Riemann problem test with the exact solution.
For this purpose, we solve the homogeneous equations (109)-(110) with initial conditions as given in Table 13. The solution of this problem is discontinuous and therefore the objective of the exercise is to test the expected ENO character of the non-linear reconstruction methods. More particularly, to test the newly presented AENO method. We expect good resolution of smooth parts of the flow, sharp discontinuities and absence, or much reduced, spurious oscillations behind the shock. Table 10 AENO results for IVP3 (106): convergence-rate study for the ADER schemes with the AENO reconstruction with 2 = 0.5 for the linear advection for = 1 , T out = 4 , C cfl = 0.9 Mesh Results are shown in Fig. 20 for ENO, Fig. 21 for m-ENO, Fig. 22 for AENO, Fig. 23 for m-AENO and Fig. 24 for WENO. The results show that the numerical solution approximates well the exact one for every order of accuracy and the discontinuity is handled well. For all orders, no visible spurious oscillations are present, except for WENO, and the approximation of the waves improves as the order of accuracy increases. To plotting accuracy, all methods give virtually identical results. The new method performs as well as the well-established ENO and WENO methods. As already noted, WENO differs slightly from ENO, mENO, AENO and m-AENO, in that small spurious oscillations are seen behind the shock for the 4th and 5th order methods. Table 11 m-AENO results for IVP3 (106): convergence-rate study for the ADER schemes with the m-AENO reconstruction with 2 = 0.5 for the linear advection for = 1 , T out = 4 , C cfl = 0.9 Mesh  Table 12 WENO results for IVP3 (106): convergence-rate study for ADER schemes with the WENO reconstruction for the linear advection with = 1 , T out = 4 and C cfl = 0.9 Mesh

Variables
Initial value

Efficiency Assessment on a Riemann Problem Test
A criterion of fundamental importance in assessing the performance of numerical methods is the cost/benefit relation. In other words, accuracy against cost, or equivalently error against cost. This is conventionally achieved by computing the solution to a test problem with the exact solution on a sequence of successively refined meshes M k , and computing the corresponding errors E k on a given norm. Associated to each pair (M k , E k ) there is a computing cost CPU k . Then, one usually fits a least-square line to the set of points (CPU k , E k ) on log-log axes. The procedure is performed for each method to be assessed, which produces a graph of several straight lines, one for each method, as shown in Fig. 25.
For the efficiency exercise here, we choose the Riemann problem, with the exact solution, given in Table 14.
Results are shown in Fig. 25 for third-order schemes for all reconstruction methods. The top frame shows results for the error in cross-sectional area A, while the bottom frame shows results for the error in velocity u. We have used the sequence of meshes M = 25, 50, 100, 200, 400, 800, 1 600 ; see symbols. Then, as is customary in the literature, we have fitted least square straight lines for each scheme. From the results, it is clearly seen that as the mesh is refined, the CPU cost increases and the error decreases. The straight lines corresponding to the schemes diverge as the error decreases. It is seen that for a given cost, the smallest error is obtained from the newly proposed AENO scheme, while the largest error is given by the WENO scheme. The error lines from the remaining schemes lie in between AENO and WENO.
To assess the error/cost relation we fix an error and draw the corresponding horizontal line. For example, in Fig. 25 we first choose a constant error denoted by E AENO 1 600 , this is the error of the AENO scheme for the finest mesh used with M = 1 600 cells. Then the intercept of this horizontal line with the least square line corresponding to AENO gives the CPU cost. For the cross-sectional area we obtain 12 827 s for AENO and 54 648 s for WENO. This means that AENO is 4.3 times more efficient than WENO.
Obviously, for even smaller errors, this difference in efficiency will increase. By extrapolating to smaller errors, assuming straight lines, we gain a clearer idea of the difference in efficiency of all methods assessed. Extrapolating the straight lines to the chosen small errors of E = 10 −8 for area A and of E = 10 −4 for velocity u we obtain estimates for the cost of the schemes; see corresponding arrows in Fig. 25. From the given numbers in Fig. 25 it is seen that the efficiency gain of AENO over WENO is 9.7 for area, while for velocity, bottom plot, this gain is 9.13. In other words, for small errors, the AENO reconstruction method is estimated to be one order of magnitude more efficient than WENO.

Convergence-Rate Study
To carry out an empirical convergence rate study we consider the method of manufactured solutions. First we prescribe a smooth vector-value function ̃ (x, t) [22] given as follows: The function ̃ (x, t) in (113) is the exact solution of the new balance law (115). We remark that for this test we considered the equations without the passive scalar.
Convergence rates are calculated with respect to the cross-sectional area A(x, t) and flow q(x, t) = A(x, t) u(x, t) . Results are shown in Tables 15,16,17,18,19,20,21,22,23 and 24. The overall results differ somehow from those for the linear advection equation. Now, essentially, the expected convergence rates are reached for all orders in all three norms, by all three reconstruction methods. In the linear advection case, the AENO method reached the expected rate only in the L 1 norm, being sub-optimal in the L 2 norm; the expected rate in the L ∞ norm was not reached.
Regarding errors, Figs. 26 and 27 show a comparison of L 1 errors for all three methods for all orders and all meshes. Figure 26 shows results for the second (top) and third (bottom) order methods. For the second-order case all three reconstruction schemes give essentially equivalent results. For the third-order schemes, surprisingly, ENO exhibits smaller errors than AENO. Figure 27 shows results for the fourth (top) and fifth   Table 15 Convergence-rate study for the blood flow equations. ADER schemes from 2nd to 5th order of accuracy with the ENO reconstruction. Rates are calculated for the cross-sectional area A. Computational parameters are: T out = 0.1 and C cfl = 0.9 Mesh (bottom) order methods. For the fourth-order schemes, all results are virtually coincident. For the fifth-order case, the AENO schemes give the smaller errors.

Shock/Turbulence Interaction Problem
Here we assess the performance of the methods for a blood-flow analogue of the so-called shock/turbulence interaction problem for gas dynamics proposed in [14]. A variation of this test that is more suitable to assess methods for long evolution times was proposed in [42] and [32]. Our test problem here is inspired on the modified version of the problem. We solve the blood flow equations with the initial condition (Figs. 28 [25,50], showing on top the shock wave and below the domain portion [28,40] displaying the physical oscillations upstream of the shock. As the order of the numerical method increases, the resolution of the physical oscillations of the solution is improved. For second order methods and all the reconstructions techniques, only two oscillations just behind the shock are captured. The WENO reconstruction method shows a small advantage over the other reconstruction schemes, for the second-order case. However, for the fourth and fifth order cases, the WENO numerical solution presents spurious oscillations and over/under shoots. For the third-order case, AENO shows a small advantage over the other methods. For the fourth-order case, AENO is comparable to ENO and m-ENO; both AENO and m-AENO show a small advantage over the remaining methods. For the fifth-order case, all methods are comparable, with the exception of WENO, which shows the worse performance, due to the spurious oscillations and over/under shoots. Table 19 Convergence-rate study for the blood flow equations. ADER schemes from 2nd to 5th order of accuracy with the AENO reconstruction with 2 = 0.5 . Rates are calculated for the cross-sectional area A. Computational parameters are: T out = 0.1 and C cfl = 0.9 Mesh  Table 20 Convergence-rate study for the blood flow equations. ADER schemes from 2nd to 5th order of accuracy with the AENO reconstruction with 2 = 0.5 . Rates are calculated for the flow Au. Computational parameters are: T out = 0.1 and C cfl = 0.9 Mesh  Table 22 Convergence-rate study for the blood flow equations. ADER schemes from 2nd to 5th order of accuracy with the m-AENO reconstruction with 2 = 0.5 . Rates are calculated for the flow Au. Computational parameters are: T out = 0.1 and C cfl = 0.9 Mesh

Conclusions
In this paper, we have presented a novel, non-linear spatial reconstruction scheme called AENO, with two variants. The method is akin to the established methods ENO, m-ENO and WENO. All reconstruction schemes have been implemented in conjunction with the fully discrete ADER scheme, from second to fifth order of accuracy, in both space and time, along with the Toro-Titarev solver for the generalized Riemann problem. The methods have been tested thoroughly for the linear advection equation and for a non-linear hyperbolic system, namely the blood flow equations. The assessment of the methods has been through comparison of profiles with exact solutions, convergence rate studies and efficiency studies. Overall, the results of the new AENO reconstruction methods are comparable to the ENO, m-ENO and WENO. However, AENO shows a distinctive advantage over ENO for long-time evolution problems; this is more obvious for second-and third-order methods, but will also be apparent for high-order methods on coarse meshes. Crucially, AENO turns out to be the most efficient of all schemes tested. Estimates reveal AENO to be up to one order of magnitude more efficient than WENO, for sought small errors. Desirable future developments are implementation and assessment of the AENO schemes in conjunction with semi-discrete methods in one and multiple space dimensions.
Funding Open access funding provided by University of Helsinki including Helsinki University Central Hospital.

Compliance with Ethical Standards
Conflicts of Interest The authors declare that they have no conflict of interest.  Fig. 34 Shock interaction test problem: ADER schemes from 2nd to 5th order with m-AENO reconstruction