Structure aware Runge–Kutta time stepping for spacetime tents

We introduce a new class of Runge–Kutta type methods suitable for time stepping to propagate hyperbolic solutions within tent-shaped spacetime regions. Unlike standard Runge–Kutta methods, the new methods yield expected convergence properties when standard high order spatial (discontinuous Galerkin) discretizations are used. After presenting a derivation of nonstandard order conditions for these methods, we show numerical examples of nonlinear hyperbolic systems to demonstrate the optimal convergence rates. We also report on the discrete stability properties of these methods applied to linear hyperbolic equations.


Introduction
For simulating wave phenomena, the state of the art relies heavily on efficient and accurate numerical solution techniques for hyperbolic systems. This paper is concerned with those solution techniques that proceed by subdividing the spacetime into tentshaped subregions satisfying a causality condition. Just as light cones are often used to delineate what is causally possible and impossible in the spacetime, tent-shaped spacetime regions are natural to impose causality when numerically solving hyperbolic equations. By constraining the height of the tent pole, erected vertically in an increasing time direction, one can ensure that the tent encloses the domain of dependence of all its points. This constraint on the tent pole height is a causality condition that a numerical scheme using such tents should satisfy. The spacetime subdivision into tents may be unstructured, thus allowing such schemes to advance in time by different amounts at different spatial locations, i.e., local time stepping can be naturally built in while subdividing the spacetime into tents.
The main contribution of this paper is a new explicit Runge-Kutta type time stepping scheme for solving hyperbolic systems within a spacetime tent. Standard time stepping methods cannot be directly applied on tents, since tents are generally not a tensor product of a spatial domain with a time interval. A non-tensor product spacetime tent can be mapped to a tensor product spacetime cylinder using a degenerate Duffy-like transformation. This is the basis of the Mapped Tent Pitching (MTP) schemes that we introduced previously in Ref. [5]. As shown there, a spacetime Piola map can be used to pull back the hyperbolic system from the tent to the spacetime cylinder. Being a tensor product domain, the spacetime cylinder, admits the use of standard explicit time stepping schemes, like the classical RK4 scheme. However, as we shall show here, expected convergence rates are not observed when such standard explicit schemes are used. The cause of this problem can be traced back to the degeneracy of the map. After illustrating this problem, we shall introduce a new Structure Aware Runge-Kutta (SARK) scheme, which overcomes this problem.
We first reported the above-mentioned order reduction in [6], where a fix was proposed for linear hyperbolic systems, called the Structure Aware Taylor (SAT) scheme. In contrast, the new SARK schemes of this paper are applicable to both linear and nonlinear hyperbolic systems.
Prior work on tent-based methods spans both the computational engineering literature [2,11] and numerical analysis literature [4,12]. These works have clearly articulated the promise of tent-based schemes, including local time stepping, even with higher order spatio-temporal discretizations, and the opportunities to utilize concurrency. Recent advances include tent-based Trefftz methods [13] and the use of asynchronous SDG (spacetime discontinuous Galerkin) methods to new engineering applications [1]. One can find explicit methods for conservation laws in the literature, like the finite volume approach in [10], which incorporate local time stepping without the usage of a tent mesh. In contrast to the presented MTP schemes, the extension to high order is yet to be done for those methods.
To place the present contribution in the perspective of these existing works, a few words regarding our focus on explicit time stepping are in order. The ratio of memory movements to flops is very low for explicit schemes, making them highly suitable for the newly emerging many-core processors. However, before the introduction of MTP schemes in Ref. [5], it was not clear that such advantages of explicit time stepping could be brought to any tent-based method. Now that we have an algorithmic avenue to perform explicit time stepping within tent-based schemes, we turn to the study of accuracy and convergence orders. Having encountered the unexpected roadblock of the above-mentioned convergence order reduction, we have been focusing on developing time stepping techniques to overcome it. This paper is an outgrowth of these studies.
In the next section, we quickly review the construction of MTP schemes, showing how the main system of ordinary differential equations (ODEs) that is the subject of this paper arises. In Sect. 3, we show why one should not use standard Runge-Kutta schemes for solving this ODE system. Then, in Sect. 4, we propose our new SARK schemes for solving the ODE system. Sect. 5 derives order conditions for these schemes. In Sect. 6, we study the discrete stability of SARK schemes. Section 7 reports on the good performance of the new schemes when applied to some standard nonlinear hyperbolic systems.

Construction of mapped tent pitching schemes
In this section we give a brief overview of MTP schemes. Let T be a simplicial conforming spatial mesh of a bounded spatial domain X 0 & R N : Spacetime tents are built atop this spatial mesh by vertically erecting tent poles at its vertices in a sequence of steps. At the ith step, a tent K i is added. It takes the form where x V is the vertex patch made up of all elements in T connected to a mesh vertex V. In (1), the function s i ðxÞ is a continuous function of the spatial coordinate x that is piecewise linear with respect to T . The graph of s i represents the advancing spacetime front at the ith step, so the tent K i may be thought of as the spacetime domain between these advancing fronts-see Fig. 1. Within K i , the distance the central vertex V can advance in time is restricted by the causality constraint where c max is an upper bound for the local wavespeed on K i (clarified further below) of a hyperbolic system under consideration. Such an upper bound is readily computable for linear systems. It is also possible to compute such bounds for some nonlinear systems, including the Euler system (without computing the full solution), as described in Ref. [8].
Once we have c max , there are multiple algorithmic options for generating the advancing fronts s i satisfying (2). One such algorithm is written out in detail in our prior work [5].
The hyperbolic systems we have in mind are general systems with L unknowns in N spatial dimensions, posed on the spacetime cylinder X :¼ X 0 Â ð0; t max Þ for some final time t max . To recall the standard settings in the literature, we follow [3], and assume we are given some g : X Â R L ! R L and f : X Â R L ! R LÂN . The standard hyperbolic problem is to find u : X ! R L such that o t gðx; t; uÞ þ div x f ðx; t; uÞ ¼ 0 where o t ¼ o=ot denotes the time derivative and div x denotes the row-wise divergence operator. We assume that the system (3) is hyperbolic in t-direction, as defined in Ref. [3]. This implies, in particular, that the first order derivatives of components of g and f exist. These derivatives, together with a direction vector, are then used to form a generalized eigenproblem, whose eigenvalues are real whenever (3)  SN Partial Differential Equations and Applications of these eigenvalues, which depend on the direction vector and x, t, u, give wavespeeds. Let c(x, t, u) denote the maximum of these wavespeeds over all directions. We choose the quantity c max in (2) to be any upper bound for these c(x, t, u) for all ðx; tÞ 2 K i . MTP schemes proceed by mapping each of the tents arising above to a spacetime cylinder. To define the mapping, we consider a general tent K over any given vertex patch x V , defined by The functions u b and u t are continuous functions that are piecewise linear on the vertex patch and may be identified as the bottom and top advancing fronts restricted to the vertex patch x V . To map the tent K to the spacetime cylinderK :¼ x V Â ð0; 1Þ, we define the transformation U :K ! K by Uðx;tÞ :¼ ðx; uðx;tÞÞ; where uðx;tÞ :¼ ð1 ÀtÞu b ðxÞ þtu t ðxÞ: ð5Þ Fðx; t; uÞ :¼ ½f ðx; t; uÞ; gðx; t; uÞ 2 R LÂðNþ1Þ ; we may write (3) as div x;t Fðx; t; uÞ ¼ 0: The spacetime divergence div x;t is a row-wise operator which applies the spatial derivative to the first N components and the temporal derivative to the last component. The wellknown Piola transformation of F, defined byF ¼ det U 0 ð Þ F U ð Þ U 0 ð Þ ÀT can be simplified after calculating the derivatives of U tô where dðxÞ ¼ u t ðxÞ À u b ðxÞ and I 2 R NÂN is the identity matrix. By the properties of the Piola map, we then immediately have div x;tF ðx;t;ûÞ ¼ 0; withû ¼ u U and the spacetime divergence div x;t on Hence we useK. Finally, as described in Ref. [5], writingF in terms of f and g, we find that (6) is equivalent to ot gðx;t;ûÞ À f ðx;t;ûÞruðtÞ ð Þ þ div x dðxÞf ðx;t;ûÞ ð Þ¼0; ð7Þ which is again a conservation law. For readability, we omit the spatial variable x and pseudo-timet from the arguments of functions in (7) and simply write which describes the evolution ofû along pseudo-time fromt ¼ 0 tot ¼ 1. Since uðx;tÞ ¼ ð1 ÀtÞu b ðxÞ þtu t ðxÞ ¼ u b ðxÞ þtdðxÞ; we may split gðûÞ À f ðûÞru into parts with and without explicit dependence on pseudotime, allowing us to rewrite (8) as This equation is the starting point for our spatial discretization. We use a discontinuous Galerkin method based on V h ¼ fv : vj K is a polynomial of degree p on all spatial elements T 2 T g: for all v h 2 V h and allt 2 ½0; 1. Here and throughout, every facet F is assigned a unit normal, simply denoted by n, whose direction is arbitrarily fixed, except when F & oX, in which case it points outward. The tracesû þ andû À ofû from either side are defined bŷ In (9), we also used a numerical flux f n on each facet F (that takes values in R L depending on valuesû þ ;û À from either side) and the jump sv h t :¼v þ h Àv À h . In these definitions, wheneverû þ falls outside X, it is prescribed using some given boundary conditions.
Let m ¼ dim V h ðx V Þ and let w i , i ¼ 1; . . .; m denote any standard local basis for V h ðx V Þ. Introducing U : ½0; 1 ! R m , consider the basis expansion Equation (9) leads to an ODE system for UðtÞ as follows. Define (possibly nonlinear) for all v h 2 V h ðx V Þ, whereû in all terms on the right hand sides above is to be expanded in terms of U using (10).
With these notations, problem (9) becomes the problem of finding U : ½0; 1 ! R m satisfying d dt Mðt; UðtÞÞ ¼ AðUðtÞÞ; given some Uð0Þ 3 Difficulty with standard time stepping In this section, we describe the problem we must overcome, thus setting the stage for the new schemes proposed in Sect. 4. The problem is that standard Runge-Kutta methods when applied to the tent system (12)-after a standard reformulation-do not give expected orders of convergence.
A standard approach to numerically solve (12) proceeds by introducing a new variable YðtÞ :¼ Mðt; UðtÞÞ: Then, using the inverse mapping M À1 ðt; ÁÞ, the primary variable is UðtÞ ¼ M À1 ðt; YðtÞÞ. Substituting this expression for U on the right hand side of (12), we can bring (12) to the standard form Standard ODE solvers, such as the classical explicit RK4 method, may now be directly applied to (14). Unfortunately, this leads to reduced convergence order, as we shall now see.
Consider the example of the one-dimensional Burgers' equation with initial values set by an inflow boundary condition an x ¼ 0, and an outflow boundary condition at x ¼ 1. Let u h ðxÞ be the numerical solution of (15) at t ¼ t max obtained using the DG spatial discretization with p ¼ 2 in (11) and solving the resulting semidiscretized ODE system using one of the standard RK methods as mentioned above. The final time t max ¼ 0:1 is chosen such that the exact solution is still a smooth function (before the onset of shock). Therefore no regularization or limiting is expected to be essential to witness high order convergence.
The exact solution at t max shown in Fig. 2a is obtained by the method of characteristics together with a Newton method. Thus one would expect the error to go to zero at a rate of Oðh 3 Þ. However our observations in Fig. 2b run counter to that expectation. Figure 2b reports the rates we observed when two standard time stepping schemes were used to solve (14), namely the two-stage (RK2) and the four-stage (RK4) explicit Runge-Kutta time stepping schemes. Although we see third order convergence for the first few refinement steps, the rate eventually drops to first order for both methods. We shall return to this example in Sect. 7.1 after developing a method without this convergence reduction.

Structure aware Runge-Kutta type methods
In this section, we develop specialized Runge-Kutta type schemes that do not show the above mentioned order loss of classical Runge-Kutta schemes.
We motivate the definition of the new scheme by reformulating (12) in terms of two variables ZðtÞ and YðtÞ, defined by YðtÞ ¼ Mðt; UðtÞÞ ¼ ZðtÞ ÀtM 1 ðUðtÞÞ: together with the initial conditions Yð0Þ ¼ Zð0Þ ¼ M 0 ðU 0 Þ: Here and throughout we use primes ( 0 ) to abbreviate d=dt. The key idea is to avoid the inversion of the time-dependent M at allt, limiting the inversion to just that of M 0 . Assuming we can compute the timeindependent inverse M À1 0 , we definẽ Then, (16) yields the following ODE system for Y and Z on 0\t\1: Convergence rates of the spatial error at t max for Ralston's method (RK2) and the classical RK4 method.  (15) where Integrating the equations of (17) from 0 to s, we obtain The new scheme, defined below, may be thought of as motivated by quadrature approximations to the integrals above. Note that we are only interested in such quadratures that result in explicit schemes. Moreover, we must also approximate sM 1 ðZðsÞÞ by an extrapolation formula that uses prior values of Z, in order to keep the scheme explicit.
Definition 1 Given an initial condition Y 0 , an s-stage SARK method for (17) computes This explicit method is determined by the coefficient matrices Hence we use instead of the standard Butcher tableau to express our scheme. Here we restrict ourselves to schemes where c 2 R s is set by the consistency condition In the next section, we shall develop a theory to choose appropriate values of a ij ; d ij ; and b i . There, Sect. 5.4 contains some specific examples of SARK scheme tableaus.

Order conditions for the scheme
Appropriate values of a ij ; d ij ; and b i can be found by order conditions obtained by matching terms in the Taylor expansions of the exact solution YðsÞ and the discrete solution Y s . To derive these order conditions we follow the general methodology laid out in Ref. [9]. For this, we need to first compute the derivatives of the exact flow (in Sect. 5.1), then the derivatives of the discrete flow (in Sect. 5.2), followed by the formulation of resulting order conditions (in Sect. 5.3).

Derivatives of the exact solution
Continuing to use primes ( 0 ) for total derivatives with respect to a single variable like d=ds, to ease the tedious calculations below, we shall also employ the nth order Frechet derivative of a function g : D & R m ! V, for some vector space V. It is denoted by g ðnÞ ðzÞ : R m Â Á Á Á Â R m ! V and defined by the symmetric multilinear form . .½v n i n for any v 1 ; . . .; v n 2 R m : Whenever g and z : ð0; 1Þ ! R m are sufficiently smooth for the derivatives below to exist continuously, we have the following formulae.

Derivatives of the discrete flow
The next task is to compute the coefficients of the Taylor expansion of the function Y s defined in (19b). The arguments Z i in (19b) are also functions of s, as given by (19a).
Therefore, in what follows, we first differentiate Z i Z i ðsÞ and then Y s . Obviously, Z i ð0Þ and Z(0) coincide, so we will focus on the first and higher derivatives of Z i at s ¼ 0. To this end, we differentiate (19a) k times to get Using (21) for k ¼ 1; 2; 3, then (20), and evaluating at s ¼ 0 we obtain ÁÀ ðd jk l ð1Þ þ a jk a ð1Þ Þðd kl l þ a lk aÞ Á : Next, we focus on Y s . By (19b), Using (20), and evaluating the resulting terms at s ¼ 0 by means of (24), we obtain ij l ð1Þ þ a ij a ð1Þ Þðd jk l þ a jk aÞ Á :

Formulation of order conditions
To obtain a specific method, we find values for a ij ; d ij and b i by matching the coefficients in the Taylor expansions of YðsÞ and Y s . Note that Y s ð0Þ ¼ Y 0 ¼ Yð0Þ, so the 0th order coefficients match. The next terms in the Taylor expansions will match if Y 0 ð0Þ ¼ Y 0 s ð0Þ. For this it is sufficient that because of (23a) and (25a). To match the third terms in the Taylor expansions, equating (23b) and (25b), a ð1Þ ðaÞ þ a ð1Þ ðlÞ ¼ X s i¼1 X j\i 2b i d ij a ð1Þ ðlÞ þ 2b i a ij a ð1Þ ðaÞ: Equating the coefficients of a ð1Þ ðaÞ and a ð1Þ ðlÞ, we conclude that Y 00 ð0Þ ¼ Y 00 If one desires to further match the next higher order terms, Y 000 s ð0Þ ¼ Y 000 ð0Þ, then the expressions in (23c) and (25c) must be equated, i.e., Thus, we have proved the following result, which summarizes our discussions on order conditions.

Examples of methods up to third order
Observe that the standard order conditions of Runge-Kutta methods are a subset of the order conditions derived in Sect. 5.3. Thus we base our SARK methods on existing Runge-Kutta methods. Below, we shall refer to an s-stage SARK method based on an existing Runge-Kutta method called ''RKname'' as ''SARK(s, RKname)''. A second order two-stage SARK method can be derived from a second order Runge-Kutta method once we find d ij satisfying the additional condition which was introduced in (27). For example, one may start with the standard explicit midpoint rule and select d 21 ¼ 1=2 to satisfy (29), thus arriving at the ''SARK(2, midpoint)'' method, listed first in Table 1. The table continues on to display further such methods obtained from other well-known second order Runge-Kutta schemes. The third order SARK methods in Table 2 are based on known third order Runge-Kutta methods with three stages. The additional coefficients d ij are chosen, such that (27)-(28) are satisfied.

Application of multiple steps within a tent
Recall that the ODE system we need to solve within one mapped tent is (17) for 0\t\1: Since thet interval is not small, we subdivide it into r subintervals and use the previously described s-stage SARK scheme within each subinterval, as described next.
We subdivide the unit interval [0, 1] into r subintervals ½t k ;t kþ1 ; k ¼ 0; 1; . . .; r À 1; wheret k ¼ k r ; and apply (1) within each subinterval as described next.  First observe that the above splitting of the unitt-interval corresponds to subdividing the original tent K, as given by (4), into r ''subtents'' (see Fig. 3) of the form where u ½k ¼ uðt k Þ. Clearly u ½0 ¼ u b and u ½r ¼ u t . We then apply (1) to each of these subtents. Accordingly, let M 0;½k be defined by (11a) after replacing u b by u ½k . Keeping the same definition of A and M 1 , letM 0;½k , and s ½k ¼t kþ1 Àt k . Then the application of (1) on each interval ½t k ;t kþ1 results in the following algorithm.  Fig. 3 Illustration of the subtent K k (shaded) defined in (30). It is the image under U of the tensor product a ijÃ ½k ðZ ½k j Þ: 3. Set Y r;s 1 ¼ Y ½r : Output this as the approximation to Y(1) at the tent top.
We conclude this section by defining the propagation operators of the above algorithm, which we shall use later. At step k, we define the (generally nonlinear) partial propagation operator T ½kþ1 : V h ðx V Þ ! V h ðx V Þ, using the intermediate quantities in the algorithm: Let the total propagation operator on the tent T : Clearly, the input and output of the algorithm are related to T by

Investigation of discrete stability
This section is devoted to remarks on the stability of the new SARK schemes. While it is common to study stability of ODE solvers by applying them to a simple scalar ODE, keeping our application of spatially varying hyperbolic solutions in mind, we consider changes in an energy-like measure on the solution UðtÞ. Recall that UðtÞ 2 R m is the coefficient vector of the basis expansion of the mapped finite element solution uðx;tÞ 2 V h ðx V Þ, as defined by (10). We limit ourselves to the case where the energy-like quantity is a norm and (the generally nonlinear operators) M; M 0 and M 1 defined in (13), (11a) and (11b), respectively, are linear, so that we may rewrite Mðt; UÞ ¼ MðtÞU using the linear operator MðtÞ :¼ M 0 ÀtM 1 : R m ! V h ðx V Þ. For many standard linear hyperbolic systems, the causality condition-see (2)-can be used to easily show that MðtÞ is identifiable with a symmetric positive definite matrix so that (33) indeed defines a norm. In the special case of gðvÞ ¼ v, we note that on flat advancing fronts, where uðx;tÞ is independent of x for some fixedt, (33) reduces to so kUðtÞk MðtÞ becomes the familiar spatial L 2 norm ofûðÁ;tÞ.

Our procedure to study linear stability
Stability of the scheme within a tent can be understood by studying the discrete analogue of the ratio kUð1Þk Mð1Þ =kUð0Þk Mð0Þ for all possible initial data U(0). This amounts to studying the norm of the discrete propagation operator for U, which we proceed to formulate. First, recall the connection between U and Y, namely YðtÞ ¼ MðtÞUðtÞ. Algorithm 1, takes as input an approximation U 0 to U(0) at the tent bottom and outputs Y r;s 1 , an approximation to Y(1) at the tent top. Hence the associated approximation to U(1) is Next, recall the discrete propagation operator defined by (31). It is now a linear operator that maps Y 0 ¼ Mð0ÞU 0 to Y r;s 1 according to (32). Define the tent propagation matrix S : R m ! R m by Clearly, (32) implies that The discrete analogue of kUð1Þk Mð1Þ =kUð0Þk Mð0Þ is kU r;s 1 k Mð1Þ =kU 0 k Mð0Þ which can be bounded using the following norm of S: It is immediate from (35) that kU r;s 1 k Mð1Þ kSk LðMð0Þ;Mð1ÞÞ kU 0 k Mð0Þ . Thus the study of stability of SARK schemes is reduced to computing estimates for the norm of S.
We now describe how we computed the norm of S for some examples below. Writinĝ vðx;tÞ ¼ P m i¼1 V i ðtÞw i ðxÞ andŵðx;tÞ ¼ P m i¼1 W i ðtÞw i ðxÞ, in analogy with the basis expansion ofû in (10), let Mt be the m Â m symmetric positive definite matrix satisfying Thus, to investigate the stability of a scheme, we computed T ½k from the scheme's Butcher-like tableau, then T by (31b), followed by S per (34), and finally, the square root of the spectral radius of M À1 0 ðS > M 1 SÞ, which equals kSk LðMð0Þ;Mð1ÞÞ as shown above. We expand on the first of these steps in the next few subsections by displaying T ½k for some SARK schemes and end this section by reporting our numerical estimates for kSk LðMð0Þ;Mð1ÞÞ for an example.

Propagation operator of two-stage SARK methods
For an arbitrary two-stage SARK method the only non-zero coefficients are b 1 ; b 2 ; a 21 ; d 21 . For a given Y ½k ¼ M 0;½k U ½k we obtain with the identity matrix I 2 R mÂm . The propagation fromt k tot kþ1 reads where we used the order conditions (26) and (27) for second order methods. This results in the propagation matrix

Propagation operator of three-stage SARK methods
A similar calculation for three-stage SARK methods, using the order conditions (26)-(28), leads to the propagation matrix with X 0 ¼ ½0; 1 2 ; t max ¼ 0:05, the flux field b ¼ ð1; 1Þ > and periodic boundary conditions. The time slab X ¼ X 0 Â ð0; t max Þ is filled with tents. Within each such tent K i , let C i denote the norm kSk LðMð0Þ;Mð1ÞÞ computed with S, M(0), and M(1) specific to that tent. We expect C i to be close to one for a stable method. Let C :¼ max where the maximum is taken over all tents in the time slab. To gain an understanding of practical stability, we examine the values of C as a function of the number of SARK stages (s), polynomial degree (p), and more importantly, the number of substeps per tent (r).
In all our numerical experiments, we observed that on each tent, for a fixed s, the norm kSk LðMð0Þ;Mð1ÞÞ tends to 1 with increasing number of substeps r, and moreover, we discovered a dependence of the following form kSk LðMð0Þ;Mð1ÞÞ ¼ 1 þ O r Às ð Þ on each tent K i . Therefore, we organize our report on numerical stability observations into plots of values of C as a function of r. We do so for two SARK methods, one with s ¼ 2 and another with s ¼ 3. The results are displayed in Fig. 4. After a prominent preasymptotic region, we observe that C, as a function of r, exhibits the rate O r Às ð Þ in all cases, except one.
The exceptional case is the case p ¼ 2 in Fig. 4b, where the stability measure approaches the ideal value of 1 much faster. We do not have an explanation for this observation.
Note that all the plotted curves in Fig. 4 shift to the top and right as p increases, i.e., the number of substeps r required to keep the same stability measure C increases with p. This behavior is akin to the p-dependence of the CFL-conditions of standard time stepping schemes.

Numerical results
In this section, we collect our observations on the performance of the new SARK schemes, on the one-dimensional Burgers' equation (in Sect. 7.1) and the two-dimensional Euler system (in Sects. 7.2-7.3). While Sect. 7.2 focuses on the study of convergence rates for a smooth Euler solution, Sect. 7.3 presents the application of SARK scheme on the computationally challenging problem of simulating a Mach 3 wind tunnel with a forwardfacing step.

Convergence rates for Burgers' equation
Let us begin by returning to the one-dimensional model problem of Sect. 3 to show that the SARK methods do not suffer from the previously described convergence order reduction. For this discussion, the equation and error e h are as in (15). We apply Algorithm 1 with SARK schemes of s ¼ 2 and s ¼ 3 stages, collect values of e h for various h and plot them in Fig. 5. The data shown in Fig. 5 was generated with the polynomial order p ¼ 2 in space and h ¼ 2 Ài =10 for i ¼ 0. . .12. The tents were built so that (2) is satisfied with c max ¼ 2. Algorithm 1 is applied with r ¼ 4 and r ¼ 10 substeps within each tent for s ¼ 2 and s ¼ 3 respectively. As h decreases, in Fig. 5a we eventually see quadratic convergence for the two-stage SARK method (although the convergence rate seems to be slightly higher in a preasymptotic regime), while the rate of the underlying standard Runge-Kutta method drops to first order. The three-stage SARK method in Fig. 5b shows cubic convergence while the rate of the underlying standard Runge-Kutta method drops to first order again. These plots clearly show the benefit of using SARK scheme over the corresponding standard Runge-Kutta scheme.  Table 1b) and the standard Ralston method.  Table 2b) and the standard Heun scheme.

Convergence rates for a 2D Euler system
Now we apply SARK methods to the Euler system. Similar to the Burgers' example, which we discussed in the previous section, we choose smooth initial data and fix a final time before the onset of shock so that no limiting is needed.
Recall that the Euler system fits into (3) with Here the functions q : X 0 ! R, m : X 0 ! R 2 and E : X 0 ! R denote the density, momentum, and total energy of a perfect gas in the spatial domain X 0 ¼ ½0; 1 2 . Furthermore, we use P ¼ 1 2 qT for the pressure, T ¼ 4 for the temperature and d ¼ 5 denotes the degrees of freedom of the gas particles. The initial values are set by q 0 ¼ 1 þ e À100ððxÀ0:5Þ 2 þðyÀ0:5Þ 2 Þ ; ð37bÞ P 0 ¼ 1 þ e À100ððxÀ0:5Þ 2 þðyÀ0:5Þ 2 Þ ; ð37dÞ and the final time t max ¼ 0:1. The data shown in Fig. 6 was generated with polynomial degree p ¼ 2 in space and mesh sizes h ¼ 0:1 Â 2 Ài , for i ¼ 0. . .6. For the tent generation c max in (2) was set to 8 and the number of substeps r ¼ 4. Since we do not have an exact solution in closed form, we compare the numerical solution computed using c max with a ''reference solution''  Table 1b) and the standard Ralston method.  Table 2b) and the standard Heun scheme. Fig. 6 Error e h as defined in (37e) over spatial degrees of freedom (dof) for SARK and standard RK methods applied to the Euler equation on tents as described in (37) computed with the higher characteristic speed 2 Á c max . The latter requires many more tents to reach the final time. Let the former and latter approximations to uðÁ; t max Þ be denoted by u h and u ref h , respectively. We define the error by This is the quantity that is plotted in Fig. 6. The errors of the two-stage SARK method and the underlying RK method is seen to diverge already for the first refinement level in Fig. 6a. While the SARK method shows the expected second convergence order, the rate of the RK method drops to first order. For the three-stage methods in Fig. 5b, we see cubic convergence for both method for the first few refinements. The convergence rate of the RK method eventually drops to first order while the SARK converges at third order.

Mach 3 wind tunnel
We conclude with the well-known benchmark example [14] of the wind tunnel with a forward-facing step onto which gas flows at Mach 3. The situation is modeled by the already described Euler system (37a), but now with the initial values q 0 ¼ 1:4; m 0 ¼ q 0 ð3; 0Þ > ; P 0 ¼ 1 ð38Þ on a spatial domain X 0 with a re-entrant corner at the edge of the forward-facing stepthe domain and the boundary conditions are exactly as illustrated in numerous previous works, see e.g., [5, Fig. 4(a)]. Our numerical experience with this problem shows that it is beneficial to use high order local time stepping. As in our prior study [5], we use a spatially refined mesh near the re-entrant corner and let the tents adapt, providing automatic local time stepping. In contrast to the standard time stepping used in Ref. [5], we now use one of the newly proposed SARK schemes. We shall apply the SARK(3, Heun) method. Unlike the study in Sect. 7.2, now we must handle multiple shocks that develop over time, so it is necessary to add some stabilization to the system. This is done by adding artificial viscosity based on the entropy residual as suggested by [7]-details of this stabilization on tents are exactly as already described in Ref. [5], so we omit them here.
One of the components of the computed solution is shown in Fig. 7. This was generated with polynomial order p ¼ 4 in space, maximal characteristic speed c max ¼ 10 and r ¼ 16 substeps within each tent. Figure 8 shows the spatial mesh with the locally refined corner. The zoom in illustrates the local refinement of the tents which comes in naturally through the causality constraint while pitching the tents. The solution component (logarithmic density) shown in Fig. 7a is comparable with the solution we previously obtained using standard methods in Ref. [5], but now due to the higher accuracy of the new SARK time integration, we obtained a similar quality solution faster (with the overall simulation time on the same processor reduced by a factor of 10). We also observed that the entropy residuals calculated off the computed solution with SARK schemes led to a significantly reduced addition of artificial viscosity. The artificial viscosity coefficients generated by the entropy residual are shown in Fig. 7b, which is about half the size of what is shown in the corresponding plot in our earlier work [5, Fig. 5].