Asymptotic properties of a class of linearly implicit schemes for weakly compressible Euler equations

In this paper we derive and analyse a class of linearly implicit schemes which includes the one of Feistauer and Kučera (J Comput Phys 224:208–221, 2007) as well as the class of RS-IMEX schemes (Schütz and Noelle in J Sci Comp 64:522–540, 2015; Kaiser et al. in J Sci Comput 70:1390–1407, 2017; Bispen et al. in Commun Comput Phys 16:307–347, 2014; Zakerzadeh in ESAIM Math Model Numer Anal 53:893–924, 2019). The implicit part is based on a Jacobian matrix which is evaluated at a reference state. This state can be either the solution at the old time level as in Feistauer and Kučera (2007), or a numerical approximation of the incompressible limit equations as in Zeifang et al. (Commun Comput Phys 27:292–320, 2020), or possibly another state. Subsequently, it is shown that this class of methods is asymptotically preserving under the assumption of a discrete Hilbert expansion. For a one-dimensional setting with some limitations on the reference state, the existence of a discrete Hilbert expansion is shown.


Introduction
We consider multi-dimensional systems of hyperbolic conservation laws that depend on a parameter ε ∈ (0, ε 0 ], ε 0 > 0 fixed, ∂ t w(x, t, ε) + ∇ • f (w(x, t, ε), ε) = 0, (1) which are stiff as ε tends to 0.Here (x, t) ∈ Ω × R + ⊂ R d × R + are the space-time variables, and is the solution vector, consisting of the conserved quantities.Here N ⊂ R m is a suitable image space depending on the problem at hand, e.g., taking into account positivity of density and the like.The function is the flux matrix.We assume that for any unit vector n ∈ R d and any w ∈ N , the Jacobian matrix f (w, ε) • n is real diagonalizable with eigenvalues λ 1 (w, ε, n), . . ., λ m (w, ε, n), and that for fixed w and n, min j=1,...m as ε → 0. A classical example is low Mach number Euler equations of gas dynamics, which is also the system that we will consider in the sequel.A key issue is the choice of time discretization.For explicit schemes, the CFL condition imposes a small time step of order O(ε∆x).This might be feasible for a very fast, highly parallel solver such as [13] for some given ε, but there exists a threshold on ε such that for any value smaller than this threshold, the restriction on ∆t becomes too demanding.Fully implicit schemes, on the other hand, necessitate solving large systems of nonlinear equations, whose condition number deteriorates as the parameter ε tends to zero.Our focus here is on IMEX (implicit-explicit) schemes [1,2,6,22], which attempt to split the system into a fast part (treated implicitly) and a slow part (treated explicitly).
Besides the questions of accuracy and efficiency, there is also a qualitative issue of change of type of the system of conservation laws as ε tends to zero.For instance, weakly compressible solutions become incompressible in this limit.An important question is whether this property holds also for the numerical approximation.
The literature on numerical methods for singularly perturbed hyperbolic conservation laws is huge.The interest of this paper is on IMEX schemes for the Euler equations; those schemes necessitate a splitting of the function f (w(x, t, ε), ε) into stiff and non-stiff parts.Possible splittings have been introduced in, e.g., [4,7,10,12,18,21,28], see also the references in the cited papers.
The linearly implicit scheme presented in [9], building heavily on the work of [8], is not of the IMEX type; it is presented for the dimensional Euler equations, so there is no (explicit) ε−dependency.Our interest here is to compare asymptotic properties of the scheme [9], which we call Dolejší-Feistauer-Kučera in this work, with the RS-IMEX scheme presented in [28].Also the latter scheme is a linearly implicit scheme, see [4,10,28].
This research has been motivated through the following observation: Although different in type, numerically, both schemes perform very well in the ε → 0 limit.For the RS-IMEX scheme, a formal asymptotic consistency analysis has been given in [16]; no such analysis has been presented for the Dolejší-Feistauer-Kučera scheme.Even more, the Dolejší-Feistauer-Kučera scheme is not designed to work with the nondimensionalized equations.Nevertheless, consider the convergence results shown in Fig. 1.These results show error of a travelling vortex computation for the isentropic Euler equations, see [3, page 122] for details on the vortex.The equations are ε−dependent, and so is the vortex.For ε → 0, the equations converge towards the incompressible isentropic Euler equations.The precise definition of solver parameters are not of importance here, we refer to [9,27] for details.It should be mentioned that they are certainly not comparable (different linear solvers, different triangulations, different numerical fluxes, different representative mesh sizes and so Representative mesh parameter Figure 1: Numerical results for the RS-IMEX scheme (left) and the Dolejší-Feistauer-Kučera scheme (right); errors are in density and momentum.It can be seen that the errors are apparently independent of ε, which is typically a good indicator for a scheme being asymptotically consistent.Errors and mesh sizes have been scaled (independently on ε), so that they begin at (1, 1).To account for the dimensions in the Dolejší-Feistauer-Kučera scheme, error in density is scaled by ε −2 .Please note that quantities on the left and on the right cannot be compared right away.
on).The key observation is that both schemes, and not only the RS-IMEX, perform very well for ε → 0, which is typically a clear indicator for a scheme being AP.
The contribution of the present paper comes in three parts: • First, we present a unified framework of the RS-IMEX (RS for reference solution) and a class of linearly implicit schemes.
After having unified the schemes, we solely focus on the full Euler equations of gas dynamics given in form (1), for the ease of presentation formulated in two dimensions, with Throughout the paper, ε denotes a reference Mach number.Here u has been defined as the velocity vector u := (u, v); the equations come along with the dimensionless equation of state: It is known that for ε → 0, the solution w converges towards the solution of the incompressible equations if initial and boundary data are so-called well-prepared, see Def. 3.8, see [23]; see also [24] for a generalization and review of the existing results and [20] for a discussion in the case of more generalized initial conditions.
• Assuming the existence of an asymptotic expansion of the discretization, we show that the semidiscrete-in-time algorithm converges for ε → 0 to a consistent discretization of the incompressible Euler equations.This property has been named in [14] asymptotic preserving (AP), we refer a reader to [19] where this property has been firstly studied.See also [11] for the so-called unified preserving schemes.
• Subsequently, we show under some restrictions that there exists an asymptotic expansion of the semidiscrete-in-time discretization.
In this paper, we work with the strong form of the equations.It is hence very important to state the following assumption: Assumption 1.1.Throughout the paper, we consider initial data and final times in such a way that w remains sufficiently smooth.
The paper is organized as follows: In Sec. 2 we introduce the so-called RS-IMEX schemes.We write them as a class of linearly implicit schemes and show that Dolejší-Feistauer-Kučera is a particular, and canonical, member of this class.Sec. 3 shows that this class of schemes is asymptotically preserving assuming the existence of a discrete Hilbert expansion.Under some restrictions on the reference state, we show in Sec. 4 that this discrete Hilbert expansion exists in one spatial dimension.Sec. 5 offers conclusion and outlook.

Linearly implicit schemes based on a reference state
In this section, we formulate a unified framework containing both the Dolejší-Feistauer-Kučera and the RS-IMEX scheme.For simplicity of exposition, we suppress the dependence on ε and rewrite (1) as be the stiff and non-stiff fluxes.Note that for fixed w R and ε, the stiff flux f (w; w R ) is linear in w.
The underlying idea is that the Jacobian matrix f contains all singular eigenvalues (of order ε −1 ), and f will hence be discretized implicitly.The Jacobian f contains eigenvalues of order ε 0 , and f will hence be discretized explicitly.We call f the stiff and f the non-stiff flux.
In the following, we introduce the RS-IMEX scheme, which is based on a reference state that is a function depending on time and space: Then the RS-IMEX scheme is given by Next we introduce a variant of the RS-IMEX scheme which is not based on a reference solution w R (t), but on a reference state w n R which is constant in the time interval [t n , t n+1 ) (but possibly variable in space): Definition 2.3 (IMEX time-discretization based on a reference state).Here we suppose that w R (t) ≡ w n R is constant in time on the interval [t n , t n+1 ).We call w n R : Ω → N the reference state.Then the RS-IMEX scheme based on a reference state is given by The following lemma considerably simplifies the form of the scheme (12).It also provides a convenient basis for a DG space discretization: Lemma 2.4 (Linearly implicit scheme based on a reference state).The scheme (12) is equivalent to the linearly implicit scheme Proof.From ( 9) and (10), Remark 2.5.Taking the reference state to be the discretization at time level n, i.e., w n R = w n , then (13) reduces to the classical linear implicit scheme If, in addition, the flux is homogeneous of degree one, i.e. f (w) = f (w)w, then This is at the basis of the Dolejší-Feistauer-Kučera scheme, proposed in [8,9] for the Euler equations of gas dynamics.
Remark 2.6.In his dissertation [15], Kaiser observed that for the full Euler equations in multiple space dimensions, the Jacobian of the non-stiff flux, f (w n R ), may have complex eigenvalues if the tangential velocities are large enough compared with the normal velocities.This was remedied in [28] by removing terms of order ε 2 from the linearized equation of state.
Remark 2.7.To simplify the notation, we will usually omit the bar at w R and write simply whenever this does not lead to confusion.From now on, it is assumed that w R is of the form given in Def.

2.3.
In Section 3, we study the asymptotic consistency of the RS-IMEX scheme given in Definition 2.3 for the two-dimensional Euler equations of gas dynamics.In Section 4, we specialize to the one-dimensional case and a constant reference solution w R and prove the existence of an asymptotic expansion for our class of linearly implicit schemes.

AP analysis
Considering the Euler fluxes (6) and defining w := (w 1 , w 2 , w 3 , w 4 ) T , one can write the two Euler fluxes in terms of w as Using this notation (1) reads Jacobi matrices of f 1 and f 2 with respect to w (written in terms of the physical variables density ρ, momentum ρu and energy E) are given by We fix the boundary conditions as follows: Assumption 3.1.In the following we assume either periodic boundary conditions or slip (wall) boundary conditions for the velocity: u• n = 0 on ∂Ω, where n is the unit outer normal to Ω.

Formal expansion of the scheme
We make the following formal assumption on the existence of a Hilbert expansion.For the validity of this assumption, we refer the reader to Sec. 4.
Assumption 3.2.We assume that the physical quantities ρ, u, E and p on each time level have a formal Hilbert expansion of the form (written e.g. for ρ n ) similarly, this is assumed for the reference state w R .
Remark 3.3.It is trivial that w R used in the RS-IMEX [28] has a Hilbert expansion, because it does not depend on ε.For the Dolejší-Feistauer-Kučera scheme [9], however, this is not clear, as w R is the solution from the previous time iterate.
Substituting the Hilbert expansions into the expressions (19) and (20) gives the expansion for s = 1, 2, where and f s,(−1) (w) = 0 for s = 1, 2. Finally, since due to the Taylor expansion, we have Taking all the expansions ( 21) -( 26) and substituting into the linearized problem (13), we gather terms according to the powers of ε.For ε −2 and ε −1 we get the following lemma.
Collecting the ε 0 terms of the mass equation from (13) gives Similarly, from the momentum equation we get and Finally from the energy equation we get We note that if we assume periodic or slip boundary conditions e.g. for u n , then the same boundary conditions hold for the individual terms in its Hilbert expansion.This can be seen (e.g. in the case of slip boundary conditions) by taking the limit ε → 0 in the boundary condition • n = 0, etc. Lemma 3.5.Assuming either slip boundary conditions for u R and u n for all n or periodic boundary conditions, the functions E n (0) and p n (0) are constant in space and independent of n.
Proof.We integrate (33) over Ω and apply Green's theorem.Since E n (0) and E n+1 (0) are constant by Lemma 3.4, we get where E corresponds to the terms under the divergence symbol in (33).Since each of these terms contains either u R,(0) , u n (0) or u n+1 (0) , all of which have zero normal component on ∂Ω, the whole boundary integral in (34) vanishes.This is the case of slip-boundary conditions, for periodic boundary conditions, the boundary integral vanishes due to spatial periodicity of all the terms.Altogether, (34) then implies E n+1 (0) = E n (0) and (29) implies p n+1 (0) = p n (0) .

Asymptotic preserving property
In this section we prove that the zero order variables from the Hilbert expansion satisfy the incompressible Euler equations.First, we start with the incompressibility.
Lemma 3.6.Assume either slip boundary conditions for u R and u n for all n or periodic boundary conditions.Let ρ n (0) and ρ R,(0) be constant in space and let is also constant in space, and ∇• u n+1 (0) = 0.
4. There holds: The same is true for the inverse of θ.
The proof of the lemma is rather straightforward, which is why we omit it here.
Using the operator θ it is possible to express the momentum at time level n + 1 as a function of p L .What we are doing here is very similar to the work of [3], in the discrete case, it could be interpreted as a Gaussian elimination procedure.Lemma 4.8.There holds: with δ(ρu) * * being a quantity that possesses a Hilbert expansion.
Proof.There holds Plugging this into the momentum equation yields (note that, again, we have omitted the time level n + 1 on the right-hand side) By δ(ρu) + we denote terms that are known to have a Hilbert expansion in ε.Rearranging terms yields Exploiting the properties of θ formulated in Lemma 4.7 yields the claim.
Based on this lemma, we can find that p L fulfills a third-order differential equation: Lemma 4.9.Let p L be given as in (62).Then p L satisfies at time level n + 1 the equation with the constants ω i being defined by and p * L ∈ H ∞ being a function that possesses a Hilbert expansion.
Proof.The proof consists of lengthy and tedious, but rather straightforward computations.The important steps are the following: • First, write δE n+1 explicitly based on (61).Use (64) and (63) to express all quantities δρ and δ(ρu) in terms of p L .Substitute E n+1 on the right-hand side by using the definition of p L in (62).Then, apply θ to the equation, which results in As above, δE * * is a smooth term having a Hilbert expansion.The constants ω l i are given by • Second, write δE n+1 explicitly, this time based on the definition of p L in (62), substitute δρ and δ(ρu) accordingly.Applying θ on both sides then yields Again, δE * * * is a smooth term with a Hilbert expansion.The constants ω r i are given by , ω r 3 =0.
Proof.Assume that ω 2 = 0 and ω 3 = 0. Then there holds Hence, The only roots of this equation are γ = −3 and γ = 0, they are hence outside the range of γ.Proof.Note that p L fulfills the equation Plugging this into (65) yields the algebraic equation for p L (k) where p * L (k) denotes the Fourier coefficients of the right-hand side.Because we know that the right-hand side has the Hilbert expansion, we also know that there exists a Hilbert expansion for p * L (k).In particular, with respect to ε, we have p * L (k) = O(1).The Fourier coefficients of p L are hence given by while for k = 0, there holds (note that ω 2 and ω 3 are not zero simultaneously!) = O(1) and the properties of θ −1 , see Lemma 4.7, also δ(ρu) n+1 can be written in terms of a Hilbert expansion.Due to (64) this property carries over to δρ n+1 .Now, as p L , δρ and δ(ρu) have the Hilbert expansions, it is clear that also δE has the Hilbert expansion, too, due to (62).This proves the claim.

Conclusion and Outlook
In this work we have introduced and analysed a class of linearly implicit methods for the discretization of the full Euler equation that unifies several already existing schemes, in particular the Dolejší-Feistauer-Kučera and the RS-IMEX scheme.We have shown that this class of methods is asymptotically consistent and exhibits a phenomenon that we call superconsistency, i.e., the consistency of the flux approximation is higher than expected.Furthermore, for a prototype example, we have shown that this unified class of methods possesses the Hilbert expansion in the case of the full Euler equations which is, to the best of our knowledge, a novel contribution.
Ongoing work focuses on the extension of the analysis, in particular the existence of the Hilbert expansion, to more general situations in multiple dimensions.It is unclear whether the Fourier analysis is then still a suitable framework, as the straightforward extension of the approach we presented here is severely more complicated and it is restricted to the periodic boundary conditions.Finally, it remains to investigate numerically the efficiency and accuracy of the proposed splittings in general experiments.

Theorem 4 . 11 .
Let γ ≥ 1. Furthermore (as in this whole section), assume that Assumptions 4.1 and 4.2 hold.Then p L fulfilling the equation (65) has a Hilbert expansion, in particular it holds p L = const +O(ε 2 ).
65); with p * L having a Hilbert expansion.Due to Lemma 4.10 ω 2 and ω 3 cannot be zero simultaneously.Because we are operating under periodic boundary conditions, we apply the Fourier expansion for p L p L (x) := k∈Z p L (k)e ikx .