Admissibility Preserving Subcell Limiter for Lax–Wendroff Flux Reconstruction

Abstract

Lax-Wendroff Flux Reconstruction (LWFR) is a single-stage, high order, quadrature free method for solving hyperbolic conservation laws. We develop a subcell based limiter by blending LWFR with a lower order scheme, either first order finite volume or MUSCL-Hancock scheme. While the blending with a lower order scheme helps to control spurious oscillations, it may not guarantee admissibility of discrete solution, e.g., positivity property of quantities like density and pressure. By exploiting the subcell structure and admissibility of lower order schemes, we devise a strategy to ensure that the blended scheme is admissibility preserving for the mean values and then use a scaling limiter to obtain admissibility of the polynomial solution. For MUSCL-Hancock scheme on non-cell-centered subcells, we develop a slope limiter, time step restrictions and suitable blending of higher order fluxes, that ensures admissibility of lower order updates and hence that of the cell averages. By using the MUSCL-Hancock scheme on subcells and Gauss-Legendre points in flux reconstruction, we improve small-scale resolution compared to the subcell-based RKDG blending scheme with first order finite volume method and Gauss-Legendre-Lobatto points. We demonstrate the performance of our scheme on compressible Euler’s equations, showcasing its ability to handle shocks and preserve small-scale structures.

Appendices

Appendix A: Admissibility of MUSCL-Hancock Scheme for General Grids

For the conservation law (1), define $\sigma \left( u_1, u_2 \right)$ as

$$\begin{aligned} \sigma \left( u_1, u_2 \right) = \max \{ \rho (f'( u_\lambda )): u_\lambda = \lambda u_1 + (1-\lambda ) u_2, \quad 0 \le \lambda \le 1\} \end{aligned}$$

where $\rho (A)$ denotes the spectral radius of matrix A. For the 2-D hyperbolic conservation law

$$\begin{aligned} u_t + f_x + g_y = 0 \end{aligned}$$

(34)

where (f, g) are Cartesian components of the flux vector; the wave speed estimates in x, y directions are defined as follows

$$\begin{aligned} \sigma _x \left( u_1, u_2 \right) = \max \{ \rho (f'( u_\lambda )) : u_\lambda = \lambda u_1 + (1-\lambda ) u_2, \quad 0 \le \lambda \le 1\} \\ \sigma _y \left( u_1, u_2 \right) = \max \{ \rho (g'( u_\lambda )) : u_\lambda = \lambda u_1 + (1-\lambda ) u_2, \quad 0 \le \lambda \le 1\} \end{aligned}$$

We assume that the admissibility set $\mathcal {U}_{\text {ad}}$ of the conservation law is a convex subset of $\mathbb {R}^d$ which can be written as (2). The following assumption is made concerning the admissibility of first order finite volume scheme.

Admissibility of first order finite volume scheme. Under the time step restriction

$$\begin{aligned} {\max _{j} \frac{\Delta t}{\Delta x_j} \sigma (u_j, u_{j+1})} \le 1 \end{aligned}$$

(35)

the first order finite volume method

$$\begin{aligned} u_j^{n+1} = u_j^n -\frac{\Delta t}{\Delta x_j}\left( f(u_j^n,u_{j+1}^n) - f(u_{j-1}^n, u_j^n) \right) \end{aligned}$$

is admissibility preserving, i.e., $u_j^n \in \mathcal {U}_{\text {ad}}$ for all j implies that $u_j^{n+1} \in \mathcal {U}_{\text {ad}}$ for all j.

1.1 A.1 Review of MUSCL-Hancock Scheme

Here we review the MUSCL-Hancock scheme for general uniform grids that need not be cell-centered (Fig. 22) in the sense that

$$\begin{aligned} x_{j+\frac{1}{2}}- x_j \ne x_j - x_{j-\frac{1}{2}}, \end{aligned}$$

(36)

for some j where $x_j$ is the solution point in finite volume element $(x_{j-\frac{1}{2}}, x_{j+\frac{1}{2}})$. The grid used in the blending limiter (Fig. 2) is a special case of (36).

Fig. 22

Non-uniform, non-cell-centered finite volume grid

For the $j^{\text {th}}$ finite volume element $(x_{j - \frac{1}{2}}, x_{j + \frac{1}{2}})$, the constant state is denoted $u_j^n$ and the linear approximation will be denoted $r^n_j (x)$. For conservative reconstruction, the linear reconstruction is given by

$$\begin{aligned} r^n (x) = u_j^n + (x - x_j) \delta _j, \qquad x \in \left( x_{j - \frac{1}{2}}, x_{j + \frac{1}{2}} \right) \end{aligned}$$

The values on left and right faces will be computed as

$$\begin{aligned} u^{n, -}_j = u_j^n + (x_{{j-\frac{1}{2}}} - x_j) \delta _j, \qquad u_j^{n, +} = u_j + (x_{j + \frac{1}{2}} - x_j) \delta _j \end{aligned}$$

(37)

We use Taylor’s expansion to evolve the solution to $t_n + \frac{1}{2}\Delta t$

$$\begin{aligned} \begin{aligned} u_j^{{n+\frac{1}{2}},-}&= u_j^{n,-}-\frac{\Delta t}{2 \Delta x_j}(f(u_j^{n, +})-f(u_j^{n, -})) \\ u_j^{{n+\frac{1}{2}},+}&= u_j^{n,+}-\frac{\Delta t}{2 \Delta x_j}(f(u_j^{n, +})-f(u_j^{n, -})) \end{aligned} \end{aligned}$$

(38)

where $\Delta x_j = x_{{j+\frac{1}{2}}} - x_{{j-\frac{1}{2}}}$. The final update is performed by using an approximate Riemann solver on the evolved quantities

$$\begin{aligned} u_j^{n + 1} = u_j^n - \frac{\Delta t}{\Delta x_j} \left( f_{j + \frac{1}{2}}^{n + \frac{1}{2}} - f^{n + \frac{1}{2}}_{j - \frac{1}{2}} \right) \end{aligned}$$

(39)

where

$$\begin{aligned} f_{j + \frac{1}{2}}^{n + \frac{1}{2}} = f\left( u_j^{n + \frac{1}{2}, +}, u_{j + 1}^{n + \frac{1}{2}, -} \right) \end{aligned}$$

is some numerical flux function. The key idea of the proof is to write the evolution $u_{j}^{{n+\frac{1}{2}}, \pm }$ from (38) as a convex combination of exact solution of some Riemann problem and the final update $u_j^{n+1}$ from (39) as a convex combination of first order finite volume updates on appropriately chosen subcells.

1.2 A.2 Primary Generalization for Proof

For the uniform, cell-centered case, Berthon [7] defined $u_j^{*,\pm }$ to satisfy

$$\begin{aligned} \frac{1}{2} u_j^{n, -} + u_j^{*, \pm } + \frac{1}{2} u_j^{n, +} = 2 u_j^{n, \pm } \end{aligned}$$

We generalize it for non-cell centred grids (36)

$$\begin{aligned} \mu _-u_j^{n, -} + u_j^{*, \pm } + \mu _+u_j^{n, +} = 2 u_j^{n, \pm } \end{aligned}$$

where

$$\begin{aligned} \mu _-= \frac{x_{{j+\frac{1}{2}}} - x_j}{x_{{j+\frac{1}{2}}} - x_{{j-\frac{1}{2}}}}, \qquad \mu _+= \frac{x_j - x_{{j-\frac{1}{2}}}}{x_{{j+\frac{1}{2}}} - x_{{j-\frac{1}{2}}}} \end{aligned}$$

(40)

This choice was made to keep the natural extension of $u_j^{*,\pm }$ in the conservative reconstruction case:

$$\begin{aligned} u_j^{*,\pm } = u_j^n + 2(x_{j\pm \frac{1}{2}} - x_j)\delta _j \end{aligned}$$

noting that $u_j^{n, \pm }$ are given by (37).

1.3 A.3 Proving Admissibility

The following lemma about conservation laws will be crucial in the proof.

Lemma 1

Consider the 1-D Riemann problem

$$\begin{aligned} u_t + f(u)_x&= 0\\ u(x, 0)&= \left\{ \begin{array}{lll} u_l, &{} \quad &{} x < 0\\ u_r, &{} &{} x > 0 \end{array}\right. \end{aligned}$$

in $[-h,h] \times [0,\Delta t]$ where

$$\begin{aligned} \frac{\Delta t}{h} \sigma (u_l, u_r) \le 1 \end{aligned}$$

(41)

Then, for all $t \le \Delta t$,

$$\begin{aligned} \int _{- h}^h u(x,t) \text {d}x = h (u_l + u_r) - t (f (u_r) - f (u_l)) \end{aligned}$$

(42)

Proof

Integrate the conservation law over $(- h, 0) \times (0, t)$

$$\begin{aligned} 0&= \int _{- h}^0 u(x, t) \text {d}x - h u_l + \int _0^t (f(u(0^-,t)) - f(u(-h, t))) \text {d}t\\&= \int _{- h}^0 u(x,t) \text {d}x - h u_l + t (f(\tilde{u} (0^-)) - f(u_l)) \end{aligned}$$

where, by self-similarity of solution of Riemann problem, $\tilde{u}$ is defined so that $u(x,t) = \tilde{u}(x/t)$ and $f(u(-h,t)) = f(u_l)$ is obtained as characteristics from [0, h] do not reach $x=-h$ due to the time restriction (41). Rewriting gives

$$\begin{aligned} \int _{- h}^0 u(x,t) \text {d}x = h u_l - t (f(\tilde{u} (0^-)) - f(u_l)) \end{aligned}$$

Similarly,

$$\begin{aligned} \int _0^h u(x,t) \text {d}x = h u_r - t (f(u_r) - f(\tilde{u}(0^+))) \end{aligned}$$

If $\tilde{u}$ is discontinuous at $x=0$, by Rankine-Hugoniot conditions, we will have a stationary jump at $x/t=0$ and obtain $f(\tilde{u}(0^+)) = f(\tilde{u}(0^-))$. The same trivially holds if $\tilde{u}$ is continuous at $x/t=0$. Thus, we can sum the previous two identities to get (42). $\square$

We will now give a criterion under which we can prove $u_j^{{n+\frac{1}{2}},\pm } \in \mathcal {U}_{\text {ad}}$, i.e., the evolution step (38) preserves $\mathcal {U}_{\text {ad}}$.

Lemma 2

Define $\mu _\pm$ by (40) and pick $u_j^{*, \pm }$ to satisfy

$$\begin{aligned} \frac{\mu _-}{2} u_j^{n, -} + \frac{1}{2} u_j^{*, \pm } + \frac{\mu _+}{2} u_j^{n, +} = u_j^{n, \pm } \end{aligned}$$

(43)

Assume $u_j^{n, \pm }, u_j^{*, \pm } \in \mathcal {U}_{\text {ad}}$ and the CFL restrictions

$$\begin{aligned} \max _{j} \frac{\Delta t}{\mu _-\Delta x_j} \sigma \left( u_j^{n, -}, u_j^{*, \pm } \right) \le 1, \qquad \max _{j} \frac{\Delta t}{\mu _+\Delta x_j} \sigma \left( u_j^{*, \pm }, u_j^{n, +} \right) \le 1 \end{aligned}$$

(44)

are satisfied. Then, $u_j^{n + \frac{1}{2}, \pm }$ given by the first step (38) of the MUSCL-Hancock scheme is in $\mathcal {U}_{\text {ad}}$.

Proof

We will prove that $u_j^{n + \frac{1}{2}, +} \in \mathcal {U}_{\text {ad}}$, and the proof for $u_j^{n + \frac{1}{2}, -}$ shall follow similarly. The key idea is to write $u_j^{n + \frac{1}{2}, {+}}$ as the exact solution of some Riemann problems. Define $u^h (x, t): (x_{j-\frac{1}{2}}, x_{j+\frac{1}{2}}) \times (0, \Delta t/2) \rightarrow \mathcal {U}_{\text {ad}}$ to be the weak solution of the Cauchy problem with initial data

$$\begin{aligned} u^h (x, 0) = {\left\{ \begin{array}{ll} u_j^{n, -}, \quad &{} {\text {if}} x \in (x_{{j-\frac{1}{2}}}, x_{j - 1 / 4})\\ u_j^{*, +}, &{} {\text {if}} x \in (x_{j - 1 / 4}, x_{j + 1 / 4})\\ u_j^{n, +}, &{} {\text {if}} x \in (x_{j + 1 / 4}, x_{{j+\frac{1}{2}}}) \end{array}\right. } \end{aligned}$$

where

$$\begin{aligned} x_{j - \frac{1}{4}} = \frac{1}{2} (x_{{j-\frac{1}{2}}} + \tilde{x}_j), \qquad x_{j + \frac{1}{4}} = \frac{1}{2} (\tilde{x}_j + x_{{j+\frac{1}{2}}}), \qquad \tilde{x}_j = x_{j - \frac{1}{2}} + \mu _-\Delta x_j \end{aligned}$$

Under our time step restrictions (44), the solution $u^h$ at time $\frac{\Delta t}{2}$ is made up of non-interacting Riemann problems centered at $x_{j \pm \frac{1}{4}}$, see Fig. 23. We take the projection of $u^h (x, \Delta t / 2)$ on piecewise-constant functions

$$\begin{aligned} \tilde{u}_j^{n + \frac{1}{2}, +}:= \frac{1}{\Delta x_j} \int _{x_{j-\frac{1}{2}}} ^{x_{j+\frac{1}{2}}} u^h \left( x, \frac{\Delta t}{2} \right) \text {d}x \end{aligned}$$

Since we assumed that the conservation law preserves $\mathcal {U}_{\text {ad}}$, we get $\tilde{u}_j^{n + \frac{1}{2}, +} \in \mathcal {U}_{\text {ad}}$. If we prove $\tilde{u}_j^{n + \frac{1}{2}, +} = u_j^{n + \frac{1}{2}, +}$, we will have our claim. Applying Lemma 42 to the two non-interacting Riemann problems, we get

$$\begin{aligned} \widetilde{u}_j^{n + \frac{1}{2}, +}&=\frac{1}{\Delta x_j} \left( \int _{x_{{j-\frac{1}{2}}}}^{\tilde{x}_j} u^h \left( x, \frac{\Delta t}{2} \right) \text {d}x + \int _{\tilde{x}_j}^{x_{{j+\frac{1}{2}}}} u^h \left( x, \frac{\Delta t}{2} \right) \text {d}x \right) \\&= \frac{1}{\Delta x_j} \left[ \frac{\tilde{x}_j - x_{{j-\frac{1}{2}}}}{2} u_j^{n, -} \!+ \!\frac{\Delta x_j}{2} u_j^{*, +} \!+\! \frac{x_{{j+\frac{1}{2}}} - \tilde{x}_j}{2} u_j^{n, +}\!-\! \frac{\Delta t}{2} \left( f \left( u_j^{n, +} \right) \!-\! f \left( u_j^{n, -} \right) \right) \right] \\&= \frac{1}{2} \left( \mu _-u_j^{n, -} + u_j^{*, +} + \mu _+u_j^{n, +} \right) - \frac{\Delta t / 2}{\Delta x_j} \left( f \left( u_j^{n, +} \right) - f \left( u_j^{n, -} \right) \right) \\&= u_j^{n, +} - \frac{\Delta t / 2}{\Delta x_j} \left( f \left( u_j^{n, +} \right) - f \left( u_j^{n, -} \right) \right) , \quad \text {using (43)}\\&= u_j^{n + \frac{1}{2}, +}, \quad \text {by (38)} \end{aligned}$$

Fig. 23

Two non-interacting Riemann problems

This proves our claim. $\square$

Now, we introduce a new variable $u_j^{n + \frac{1}{2}, *}$ defined as follows:

$$\begin{aligned} \mu _-u_j^{n + \frac{1}{2}, -} + u_j^{n + \frac{1}{2}, *} + \mu _+u_j^{n + \frac{1}{2}, +} = 2 u_j^n \end{aligned}$$

(45)

Fig. 24

Finite volume evolution

As illustrated in Fig. 24, we evolve each state according to the associated first order scheme to define the following

$$\begin{aligned} \begin{aligned} u_j^{n + 1, -}&= u_j^{n + \frac{1}{2}, -} - \frac{\Delta t}{\mu _-\Delta x_j / 2} \left( f \left( u_j^{n + \frac{1}{2}, -}, u_j^{n + \frac{1}{2}, *} \right) - f \left( u_{j - 1}^{n + \frac{1}{2}, +}, u_j^{n + \frac{1}{2}, -} \right) \right) \\ u_j^{n + 1, *}&= u_j^{n + \frac{1}{2}, *} - \frac{\Delta t}{\Delta x_j / 2} \left( f \left( u_j^{n + \frac{1}{2}, *}, u_j^{n + \frac{1}{2}, +} \right) - f \left( u_j^{n + \frac{1}{2}, -}, u_j^{n + \frac{1}{2}, *} \right) \right) \\ u_j^{n + 1, +}&= u_j^{n + \frac{1}{2}, +} - \frac{\Delta t}{\mu _+\Delta x_j / 2} \left( f \left( u_j^{n + \frac{1}{2}, +}, u_{j + 1}^{n + \frac{1}{2}, -} \right) - f \left( u_j^{n + \frac{1}{2}, *}, u_j^{n + \frac{1}{2}, +} \right) \right) \end{aligned} \end{aligned}$$

(46)

Recall that (39) is

$$\begin{aligned} u_j^{n + 1} = u_j^n - \frac{\Delta t}{\Delta x_j} \left( f\left( u_j^{n + \frac{1}{2}, +}, u_{j + 1}^{n + \frac{1}{2}, -} \right) - f\left( u_{j - 1}^{n + \frac{1}{2}, +}, u_j^{n + \frac{1}{2}, -} \right) \right) \end{aligned}$$

Using (45) and (46), we get

$$\begin{aligned} \frac{\mu _-}{2} u_j^{n + 1, -} + \frac{1}{2} u_j^{n + 1, *} + \frac{\mu _+}{2} u_j^{n + 1, +} = u_j^{n + 1} \end{aligned}$$

Thus, assuming $u_j^{n + \frac{1}{2}\pm }, u_j^{n + \frac{1}{2}, *} \in \mathcal {U}_{\text {ad}}$ for all j, and since $\frac{1}{2}\mu _-+ \frac{1}{2}\mu _+= 1$, we get $u_j^{n+1} \in \mathcal {U}_{\text {ad}}$ under the following time step restrictions arising from the assumed time step requirement (35) for admissibility of the first order finite volume method

$$\begin{aligned} \begin{aligned}&\max _j \frac{\Delta t}{\mu _-\Delta x_j / 2} \sigma \left( u_j^{n + \frac{1}{2}, -}, u_j^{n + \frac{1}{2}, *} \right) \le 1,\\&\max _j \frac{\Delta t}{\mu _-\Delta x_j / 2} \sigma \left( u_{j - 1}^{n + \frac{1}{2}, +}, u_j^{n + \frac{1}{2}, -} \right) \le 1,\\&\max _j \frac{\Delta t}{\Delta x_j / 2} \sigma \left( u_j^{n + \frac{1}{2}, *}, u_j^{n + \frac{1}{2}, +}\right) \le 1, \end{aligned} \qquad \begin{aligned}&\max _j \frac{\Delta t}{\Delta x_j / 2} \sigma \left( u_j^{n + \frac{1}{2}, -}, u_j^{n + \frac{1}{2}, *}\right) \le 1 \\&\max _j \frac{\Delta t}{\mu _+\Delta x_j / 2} \sigma \left( u_j^{n + \frac{1}{2}, +}, u_{j + 1}^{n + \frac{1}{2}, -} \right) \le 1\\&\max _j \frac{\Delta t}{\mu _+\Delta x_j / 2} \sigma \left( u_j^{n + \frac{1}{2}, *}, u_j^{n + \frac{1}{2}, +} \right) \le 1 \end{aligned}\nonumber \\ \end{aligned}$$

(47)

This can be summarised in the following Lemma.

Lemma 3

Assume that the states $\left\{ u_j^{n + \frac{1}{2}, \pm } \right\} _j$, $\left\{ u_j^{n + \frac{1}{2}, *} \right\} _j$ belong to $\mathcal {U}_{\text {ad}}$, where $u_j^{n + \frac{1}{2}, *}$ is defined as in (45). Then, the updated solution $u_j^{n+1}$ of MUSCL-Hancock scheme (37–39) is in $\mathcal {U}_{\text {ad}}$ under the CFL conditions (47).

Since Lemma 2 states that $u_j^{n + \frac{1}{2}, \pm } \in \mathcal {U}_{\text {ad}}$ if $u_j^{*, \pm } \in \mathcal {U}_{\text {ad}}$, the only new condition pertains to $u_j^{n + \frac{1}{2}, *}$. Our goal now is to understand this condition, and ultimately prove that it follows from the requirement that $u_j^{*,\pm } \in \mathcal {U}_{\text {ad}}$ in case of conservative reconstruction.

Recall that $u_j^{n + \frac{1}{2}, *}$ was defined by (45); expanding the definition of $u_j^{n + \frac{1}{2}, \pm }$ given by (38) yields

$$\begin{aligned} u_j^{n + \frac{1}{2}, *} = 2 u_j^n - \left( \mu _-u_j^{n, -} + \mu _+u_j^{n, +} \right) - \frac{\Delta t}{2 \Delta x_j} (f (u_j^{n, -}) - f (u_j^{n, +})) \end{aligned}$$

(48)

This identity (48) will be seen as an evolution update similar to (38) with $u_j^{n, +}$ and $u_j^{n, -}$ being swapped and $u_j^n$ replaced with $2 u_j^n - \left( \mu _-u_j^{n, -} + \mu _+u_j^{n, +} \right)$. The admissibility of $u_j^{n + \frac{1}{2}, *}$ will be studied by adapting the proof of admissibility for (38), accounting for the differences in the case of (48). Define $u_j^{*,*}$ so that

$$\begin{aligned} \frac{\mu _-}{2} u_j^{n,-} + \frac{1}{2}u_j^{*,*} + \frac{\mu _+}{2} u_j^{n,+} = 2u_j^n-(\mu _-u_j^{n,-}+\mu _+u_j^{n,+}) \end{aligned}$$

(49)

i.e.,

$$\begin{aligned} u_j^{*,*} = 4 u_j^n - 3(\mu _-u_j^{n,-}+\mu _+u_j^{n,+}) \end{aligned}$$

(50)

The following Lemma extends the proof of Lemma 2 to obtain conditions for $u_j^{n + \frac{1}{2}, *} \in \mathcal {U}_{\text {ad}}$.

Lemma 4

Assume that $u_j^n \in \mathcal {U}_{\text {ad}}$ for all j. Consider the reconstructions $u_j^{n,\pm }$ and the $u_j^{*, *}$ defined in (49). Assume $u_j^{n, \pm }, u_j^{*, *} \in \mathcal {U}_{\text {ad}}$ and the time step restrictions

$$\begin{aligned} \max _{j} \frac{\Delta t}{\mu _-\Delta x_j} \sigma \left( u_j^{*, *}, u_j^{n, -} \right) \le 1, \qquad \max _{j} \frac{\Delta t}{\mu _+\Delta x_j} \sigma \left( u_j^{n, +}, u_j^{*, *} \right) \le 1 \end{aligned}$$

(51)

Then $u_j^{n + \frac{1}{2}, *} \in \mathcal {U}_{\text {ad}}$.

Proof

We will use the identity which follows from (48, 49)

$$\begin{aligned} u_j^{n + \frac{1}{2}, *} = \frac{\mu _-u_j^{n,-} + u_j^{*,*} + \mu _+u_j^{n,+}}{2} - \frac{\Delta t}{2 \Delta x_j} (f (u_j^{n, -}) - f (u_j^{n, +})) \end{aligned}$$

(52)

to fall back to previous case of Lemma 2.

Define $u^h (x, t): (x_{j-\frac{1}{2}}, x_{j+\frac{1}{2}}) \times (0, \Delta t/2) \rightarrow \mathcal {U}_{\text {ad}}$ to be the weak solution of the Cauchy problem with initial data

$$\begin{aligned} u^h (x, 0) = {\left\{ \begin{array}{ll} u_j^{n, +}, \quad &{} {\text {if}} x \in (x_{j-\frac{1}{2}}, x_{j - 1 / 4})\\ u_j^{*, *}, &{} {\text {if}} x \in (x_{j - \frac{1}{4}}, x_{j + 1 / 4})\\ u_j^{n, -}, &{} {\text {if}} x \in (x_{j + \frac{1}{4}}, x_{j+\frac{1}{2}}) \end{array}\right. } \end{aligned}$$

where

$$\begin{aligned} x_{j - \frac{1}{4}} = \frac{1}{2} (x_{j-\frac{1}{2}}+ x_j), \qquad x_{j + \frac{1}{4}} = \frac{1}{2} (x_j + x_{j+\frac{1}{2}}) \end{aligned}$$

Note that we have already accounted for the swapped $u_j^{n,-}$ and $u_j^{n,+}$ while defining this initial condition, see Fig. 25.

Under the assumed CFL conditions (51), the solution $u^h$ at time $\frac{\Delta t}{2}$ is made up of non-interacting Riemann problems centered at $x_{j \pm \frac{1}{4}}$. Take the projection of $u^h (x, t / 2)$ on piecewise-constant functions

$$\begin{aligned} \widetilde{u}_j^{n + \frac{1}{2}, *}:= \frac{1}{\Delta x_j} \int _{x_{j - \frac{1}{2}}}^{x_{{j+\frac{1}{2}}}} u^h \left( x, \frac{\Delta t}{2} \right) \text {d}x \in \mathcal {U}_{\text {ad}}\end{aligned}$$

Fig. 25

Two non-interacting Riemann problems

As in Lemma 2, we will show $u_j^{n + \frac{1}{2}, *} \in \mathcal {U}_{\text {ad}}$ by showing $u_j^{n + \frac{1}{2}, *} = \tilde{u}_j^{n + \frac{1}{2}, *}$. Applying Lemma 42 to the two non-interacting Riemann problems, we get

$$\begin{aligned} \widetilde{u}_j^{n + \frac{1}{2}, *}&= \frac{1}{\Delta x_j} \left( \int _{x_{{j-\frac{1}{2}}}}^{x_j} u^h \left( x, \frac{\Delta t}{2} \right) \text {d}x + \int _{x_j}^{x_{{j+\frac{1}{2}}}} u^h \left( x, \frac{\Delta t}{2} \right) \text {d}x \right) \\&= \frac{1}{\Delta x_j} \left( \frac{x_j - x_{{j-\frac{1}{2}}}}{2} u_j^{n, +} + \frac{\Delta x_j}{2} u_j^{*, *} + \frac{x_{{j+\frac{1}{2}}} - x_j}{2} u_j^{n, -} - \frac{\Delta t}{2} \left( f \left( u_j^{n, -} \right) - f \left( u_j^{n, +} \right) \right) \right) \\&= \frac{1}{2} \left( \mu _+u_j^{n, +} + u_j^{*, *} + \mu _-u_j^{n, -} \right) - \frac{\Delta t / 2}{\Delta x_j} \left( f \left( u_j^{n, -} \right) - f \left( u_j^{n, +} \right) \right) \\&=u_j^{n + \frac{1}{2}, *},\quad \text {by (52)} \end{aligned}$$

This proves our claim. $\square$

For conservative reconstruction,

$$\begin{aligned} \mu _-u_j^{n, -} + \mu _+u_j^{n, +} = u_j^n \end{aligned}$$

and thus by (50), $u_j^{*, *} = u_j^n$. The previous lemma can thus be specialized as follows.

Lemma 5

Assume that $u_j^n \in \mathcal {U}_{\text {ad}}$ and $u_j^{n, \pm } \in \mathcal {U}_{\text {ad}}$ for all j with conservative reconstruction. Also assume the CFL restrictions

$$\begin{aligned} \max _{j} \frac{\Delta t}{\mu _-\Delta x_j} \sigma \left( u_j^n, u_j^{n, -} \right) \le 1, \qquad \max _{j} \frac{\Delta t}{\mu _+\Delta x_j} \sigma \left( u_j^{n, +}, u_j^n \right) \le 1 \end{aligned}$$

(53)

where $\mu ^{\pm }$ are defined in (40). Then, $u_j^{n + \frac{1}{2}, *}$ defined in (45) is in $\mathcal {U}_{\text {ad}}$.

Combining Lemmas 2, 3, 5, we obtain the final criterion for admissibility preservation of MUSCL-Hancock with conservative reconstruction in the following Theorem 3.

Theorem 3

Let $u_j^n \in \mathcal {U}_{\text {ad}}$ for all j and $u_j^{n, \pm }$ be the conservative reconstructions defined as

$$\begin{aligned} u_j^{n, +} = {u_j^n} + (x_{{j+\frac{1}{2}}} - x_j) \delta _j, \qquad u_j^{n, -} = {u_j^n} + (x_{{j-\frac{1}{2}}} - x_j) \delta _j \end{aligned}$$

so that $u_j^{*, \pm }$ defined in (43) is also given by

$$\begin{aligned} u_j^{*, \pm } = u_j^n + 2 ( x_{j \pm \frac{1}{2}} - x_j ) \delta _j \end{aligned}$$

(54)

Assume the slope $\delta _j$ is chosen such that $u_j^{*, \pm } \in \mathcal {U}_{\text {ad}}$ and the CFL restrictions (44, 47, 53) hold. Then, the updated solution $u_j^{n + 1}$, defined by MUSCL-Hancock scheme (39) is in $\mathcal {U}_{\text {ad}}$.

Proof

Once we obtain $u_j^{n, \pm } \in \mathcal {U}_{\text {ad}}$, the claim follows from Lemmas 2–5. To prove that $u_j^{n, \pm }$ is indeed in $\mathcal {U}_{\text {ad}}$, we make the straight forward observation that

$$\begin{aligned} u_j^{n, \pm } = \frac{1}{2} u_j^{*, \pm } + \frac{1}{2} u_j^n \end{aligned}$$

Since $u_j^{*, \pm }$ and $u_j^n$ are in $\mathcal {U}_{\text {ad}}$, the proof is completed by the convex property of $\mathcal {U}_{\text {ad}}$. $\square$

Remark 2

The strictest time step restriction for admissibility of the MUSCL-Hancock scheme is imposed by (47). Thus, we can find the CFL coefficient for grid used by subcell-based blending scheme (10) by minimizing the denominator in (47) which is given by

$$\begin{aligned} \frac{1}{2} \min _{j=0,\dots ,N}\left( \xi _j - \sum _{k=0}^{j-1} w_k\right) w_j = \frac{1}{2} \xi _0w_0 \end{aligned}$$

where $\xi _0,w_0$ are the first Gauss-Legendre quadrature point (3) and weight in [0, 1]. This coefficient is less than half of the optimal CFL coefficient that arises from Fourier stability analysis of the LWFR scheme with D2 dissipation, see Table 1 of [5].

1.4 A.4 Non-conservative Reconstruction

To maintain the simple admissibility criterion (Theorem 3), we have restricted ourselves to conservative reconstruction in this work. In this section, we explain the complexities that will arise in enforcing admissibility if we perform reconstruction with non-conservative variables $v$ defined by the change of variables formula

$$\begin{aligned} v= \kappa \left( u\right) \end{aligned}$$

The linear approximation is given by

$$\begin{aligned} r^n (x) = v^n_j + (x - x_j)\delta _j, \quad x \in [x_{{j-\frac{1}{2}}}, x_{{j+\frac{1}{2}}}] \end{aligned}$$

and thus the trace values are

$$\begin{aligned} v_j^{n, \pm } = v_j^n + (x_{{j\pm \frac{1}{2}}} - x_j)\delta _j \end{aligned}$$

Since the arguments of proof of admissibility depend on constraints on the conservative variables, we have to take the inverse map on our reconstructions. For example, conservative variables at the face are obtained as

$$\begin{aligned} u_j^{n, \pm } = \kappa ^{- 1} (v_j^{n, \pm }) \end{aligned}$$

(55)

Due to the non-linearity of the map $\kappa$, unlike the conservative case, we have

$$\begin{aligned} \mu _-u_j^{n, -} + \mu _+u_j^{n, +} \ne u_j^n \end{aligned}$$

which is why several reductions of admissibility constraints will fail. The admissibility criteria for non-conservative reconstruction is stated in Theorem 4.

Theorem 4

Assume that $u_j^n \in \mathcal {U}_{\text {ad}}$ for all j. Consider $u_j^{n,\pm }$ defined in (55), $u_j^{*,\pm }$ defined in (43) and $u_j^{*,*}$ defined so that

$$\begin{aligned} \frac{\mu _-}{2}u_j^{n,-} + \frac{1}{2} u_j^{*,*} + \frac{\mu _+}{2} u_j^{n,+} = 2u_j^n-(\mu _-u_j^{n,-} + \mu _+u_j^{n,+}) \end{aligned}$$

Assume that the slope $\delta _j$ is chosen so that $u_j^{n,\pm }, u_j^{*,\pm }, u_j^{*,*} \in \mathcal {U}_{\text {ad}}$ and that the CFL restrictions (44, 47, 51) are satisfied. Then the updated solution $u_j^{n+1}$ of MUSCL-Hancock scheme (39) is in $\mathcal {U}_{\text {ad}}$.

1.5 A.5 MUSCL-Hancock Scheme in 2-D

Consider the 2-D hyperbolic conservation law (34) with fluxes $f, g$. For simplicity, assume that the reconstruction is performed on conservative variables. Thus, the linear reconstructions are given by

$$\begin{aligned} r^n_{ij} (x, y) = u_{ij}^n + (x - x_i) \delta _i^x + (y - y_j) \delta _j^y, \end{aligned}$$

and the approximations at the face $u^{n, +x}, u^{n, -x}, u^{n, +y}, u^{n, -y}$ are

$$\begin{aligned} \begin{aligned} u^{n, \pm x}_{ij}&= r^n_{ij} (x_{{i\pm \frac{1}{2}}}, y_j) = {u}_{ij}^n + (x_{{i\pm \frac{1}{2}}} - x_i)\delta _i^x\\ u^{n, \pm y}_{ij}&= r^n_{ij} (x_i, y_{{j\pm \frac{1}{2}}}) = {u}_{ij}^n + (y_{{j\pm \frac{1}{2}}} - y_j) \delta _j^y \end{aligned} \end{aligned}$$

(56)

and the derivative approximations are given by

$$\begin{aligned}&\partial _x f_{i j}:= \frac{1}{\Delta x_i} \left( f\left( u^{n, +x}_{{ij}} \right) - f\left( u^{n, -x}_{ij} \right) \right) , \partial _y g_{i j}:= \frac{1}{\Delta y_j} \left( g\left( u^{n, +y}_{{ij}} \right) - g\left( u^{n, -y}_{ij} \right) \right) \\&\partial _t u_{i j}^n:= - \partial _x f_{i, j} - \partial _y g_{i, j} \end{aligned}$$

The evolutions to time level ${n+\frac{1}{2}}$ are given by

$$\begin{aligned} u^{{n+\frac{1}{2}}, \pm x}_{i j} = u^{n, \pm x}_{i j} + \frac{\Delta t}{2} \partial _t u_{i j}^n,\qquad u^{{n+\frac{1}{2}}, \pm y}_{i j} = u^{n, \pm y}_{i j} + \frac{\Delta t}{2} \partial _t u_{i j}^n \end{aligned}$$

(57)

and then the final update is performed as

$$\begin{aligned} u_{i j}^{n + 1} = u_{i j}^n - \frac{\Delta t}{\Delta x_i} ( f_{{i+\frac{1}{2}}, j}^{n+\frac{1}{2}}- f_{{i-\frac{1}{2}}, j}^{{n+\frac{1}{2}}}) - \frac{\Delta t}{\Delta y_j} ( g_{i, {j+\frac{1}{2}}}^{n+\frac{1}{2}}- g_{i,{j-\frac{1}{2}}}^{{n+\frac{1}{2}}}) \end{aligned}$$

(58)

where the numerical fluxes are computed as

$$\begin{aligned} f_{{i+\frac{1}{2}}, j}^{n+\frac{1}{2}}= f\left( u^{{n+\frac{1}{2}}, +x}_{i j}, u^{{n+\frac{1}{2}}, -x}_{i+1, j} \right) , \qquad g_{i, {j+\frac{1}{2}}}^{n+\frac{1}{2}}= g\left( u^{{n+\frac{1}{2}}, +y}_{i j}, u^{{n+\frac{1}{2}}, -y}_{i, j+1} \right) \end{aligned}$$

1.5.1 A.5.1 First Evolution Step

As in 1-D, define $u^{*, \pm x}_{i j}, u^{*, \pm y}_{ij}$ so that

$$\begin{aligned} \begin{aligned} \mu _{+x}u^{n, +x}_{ij} + u^{*, \pm x}+ \mu _{-x}u^{n, -x}_{ij}&= 2 u^{n, \pm x}_{ij} \\ \mu _{+y}u^{n, +y}_{ij} + u^{*, \pm y}+ \mu _{-y}u^{n, -y}_{ij}&= 2 u^{n, \pm y}_{ij} \end{aligned} \end{aligned}$$

(59)

where

$$\begin{aligned} \begin{aligned} \mu _{+x}&= \frac{x_i - x_{{i-\frac{1}{2}}}}{x_{{i+\frac{1}{2}}} - x_{{i-\frac{1}{2}}}}, \qquad \mu _{-x}= \frac{x_{{i+\frac{1}{2}}} - x_i}{x_{{i+\frac{1}{2}}} - x_{i - 1 /2}}\\ \mu _{+y}&= \frac{y_j - y_{{j-\frac{1}{2}}}}{y_{{j+\frac{1}{2}}} - y_{{j-\frac{1}{2}}}}, \qquad \mu _{-y}= \frac{y_{{j+\frac{1}{2}}} - y_j}{y_{{j+\frac{1}{2}}} - y_{{j-\frac{1}{2}}}} \end{aligned} \end{aligned}$$

(60)

Since we assume conservative reconstruction

$$\begin{aligned} \mu _{+x}u^{n, +x}_{i j} + \mu _{-x}u^{n, -x}_{i j} = \mu _{+y}u^{n, +y}_{i j} + \mu _{-y}u^{n, -y}_{i j} = u_{i j}^n \end{aligned}$$

Thus, we have

$$\begin{aligned} u^{*, \pm x}_{i j} = u_{i j} + 2 (x_{i\pm \frac{1}{2}}- x_i) \delta _i^x, \qquad u^{*, \pm y}_{i j} = u_{i j} + 2 (y_{{j\pm \frac{1}{2}}} - x_j) \delta _j^y \end{aligned}$$

We will particularly discuss admissibility of the updates

$$\begin{aligned} u^{{n+\frac{1}{2}}, +x}_{i j} = u^{n, +x}_{i j} - \frac{\Delta t / 2}{\Delta x_i} \left( f\left( u^{n, +x}_{i j} \right) - f\left( u^{n, -x}_{ij} \right) \right) - \frac{\Delta t / 2}{\Delta y_j} \left( g\left( u^{n, +y}_{i j} \right) - g\left( u^{n, -y}_{i j} \right) \right) \end{aligned}$$

(61)

Admissibility of the other three updates $u^{{n+\frac{1}{2}}, -x}_{i j}, u^{{n+\frac{1}{2}}, \pm y}_{i j}$ will follow similarly. For some $k_x, k_y$ chosen such that $k_x + k_y = 1,$ we write (61) as

$$\begin{aligned} u^{{n+\frac{1}{2}}, +x}_{i j} = k_x \mathcal {\theta }_{ij}^{+x}+ k_y \mathcal {\theta }_{ij}^{+y}\end{aligned}$$

where

$$\begin{aligned} \mathcal {\theta }_{ij}^{+x}:=u^{n, +x}_{i j} - \frac{\Delta t / 2}{k_x \Delta x_i} \left( f\left( u^{n, +x}_{i j} \right) - f\left( u^{n, -x}_{i j} \right) \right) \end{aligned}$$

(62)

and

$$\begin{aligned} \mathcal {\theta }_{ij}^{+y}:= u^{n, +x}_{i j} - \frac{\Delta t / 2}{k_y \Delta y_j} \left( g\left( u^{n, +y}_{i j} \right) - g\left( u^{n, -y}_{i j} \right) \right) \end{aligned}$$

(63)

We will choose the slopes $\delta ^x_i, \delta ^y_j$ and time step $\Delta t$ so that $\mathcal {\theta }_{ij}^{+x}, \mathcal {\theta }_{ij}^{+y}\in \mathcal {U}_{\text {ad}}$. Then, we can take convex combinations of the two terms to obtain admissibility of $u^{{n+\frac{1}{2}}, +x}_{i j}$. The choice of $k_x, k_y$ will not influence the slope restriction, but only the time step restriction required to obtain admissibility. In this work, we only use the Fourier CFL restriction imposed by the Lax–Wendroff scheme (32) and observe admissibility preservation in all our numerical experiments and thus do not study the choice of $k_x,k_y$. However, in a similar context, [78] proposed the choice of

$$\begin{aligned} k_x = \frac{a_x / \Delta x_i}{a_x / \Delta x_i + a_y / \Delta y_j}, \qquad k_y = \frac{a_y / \Delta y_j }{a_x / \Delta x_i + a_y / \Delta y_j} \end{aligned}$$

(64)

where

$$\begin{aligned} a_x = \sigma _x(u^{n, -x}_{i j}, u^{n, +x}_{i j}), \qquad a_y = \sigma _y(u^{n, -y}_{i j}, u^{n, +y}_{i,j}) \end{aligned}$$

(65)

In [23], it was shown that the time step restriction imposed by the above decomposition is suboptimal and optimal decompositions were proposed. After choosing $k_x, k_y$, following the 1-D procedures from Sect. 1, the slopes $\delta ^x_i, \delta ^y_j$ will be limited to enforce admissibility of $\mathcal {\theta }_{ij}^{+x}, \mathcal {\theta }_{ij}^{+y}$ (62, 63). The admissibility preservation of $\mathcal {\theta }_{ij}^{+x}$ (62) follows directly from the arguments used in Lemma 2, enforcing slope restriction so that $u^{n, \pm x}_{i j}$ and $u^{*, +x}_{i j}$ are admissible, and appropriate time step restrictions. For admissibility of $\mathcal {\theta }_{ij}^{+y}$ (63), we define $u^{*, +xy}_{ij}$ so that

$$\begin{aligned} \mu _{+y}u^{n, +y}_{i j} + u^{*, +xy}_{ij} + \mu _{-y}u^{n, -y}_{i j} = 2 u^{n, +x}_{i j} \end{aligned}$$

Thus, the proof of Lemma 2 shall apply as in 1-D under the assumption of admissibility of $u^{n, \pm y}_{i j}, u^{*, +xy}_{ij}$ and some CFL conditions. Thus, we will have admissibility of $\mathcal {\theta }_{ij}^{+y}\in \mathcal {U}_{\text {ad}}$. We obtain further simplifications because of conservative reconstructions

$$\begin{aligned} u^{*, +xy}_{ij} = u^{*, +x}_{ij} \end{aligned}$$

(66)

and thus the slope limiting for enforcing admissibility of $u^{*, +x}_{ij}$ will suffice. We note the precise slope and CFL restrictions are in Lemma 6.

Lemma 6

For $\mu _{\pm x}, \mu _{\pm y}$ defined in (60), $u^{n, \pm x}_{ij}, u^{n, \pm y}_{ij}$ reconstructed in (56), $u^{*, \pm x}_{ij}, u^{*, \pm y}_{ij}$ picked as in (59), assume

$$\begin{aligned} u^{n, \pm x}_{ij}, u^{n, \pm y}_{ij}, u^{*, \pm x}_{ij}, u^{*, \pm y}_{ij} \in \mathcal {U}_{\text {ad}}\end{aligned}$$

and the CFL restrictions

$$\begin{aligned} \begin{aligned} \max _{i,j} \frac{\lambda _{x_i}}{\mu _{-x}} \sigma _x \left( u^{n, -x}_{ij}, u^{*, \pm x}_{ij} \right) \le 1, \qquad \max _{i,j} \frac{\lambda _{x_i}}{\mu _{+x}} \sigma _x\left( u^{*, \pm x}_{ij}, u^{n, +x}_{ij} \right) \le 1 \\ \max _{i,j} \frac{\lambda _{y_j}}{\mu _{-y}} \sigma _y \left( u^{n, -y}_{ij}, u^{*, \pm x}_{ij} \right) \le 1, \qquad \max _{i,j} \frac{\lambda _{y_j}}{\mu _{+y}} \sigma _y \left( u^{*, \pm x}_{ij}, u^{n, +y}_{ij} \right) \le 1 \\ \max _{i,j} \frac{\lambda _{y_j}}{\mu _{-y}} \sigma _y \left( u^{n, -y}_{ij}, u^{*, \pm y}_{ij} \right) \le 1, \qquad \max _{i,j} \frac{\lambda _{y_j}}{\mu _{+y}} \sigma _y \left( u^{*, \pm y}_{ij}, u^{n, +y}_{ij} \right) \le 1 \\ \max _{i,j} \frac{\lambda _{x_i}}{\mu _{-x}} \sigma _x \left( u^{n, -x}_{ij}, u^{*, \pm y}_{ij} \right) \le 1, \qquad \max _{i,j} \frac{\lambda _{x_i}}{\mu _{+x}} \sigma _x \left( u^{*, \pm y}_{ij}, u^{n, +x}_{ij} \right) \le 1 \end{aligned} \end{aligned}$$

(67)

where $\lambda _{x_i} = \frac{\Delta t}{k_x \Delta x_i}, \lambda _{y_j} = \frac{\Delta t}{k_y \Delta y_j}$ for all i, j and $k_x + k_y = 1$. Then, the updates $u^{{n+\frac{1}{2}}, \pm x}_{i j}, u^{{n+\frac{1}{2}}, \pm y}_{i j}$ (61) of the first step of 2-D MUSCL-Hancock scheme are admissible.

1.5.2 A.5.2 Finite volume step

The final update is given by

$$\begin{aligned} u_{i j}^{n + 1} = u_{i j}^n - \frac{\Delta t}{\Delta x_i} ( f_{{i+\frac{1}{2}}, j}^{n + \frac{1}{2}} - f_{{i-\frac{1}{2}}, j}^{n + \frac{1}{2}}) - \frac{\Delta t}{\Delta y_j} ( g_{i, {j+\frac{1}{2}}}^{n + \frac{1}{2}} - g_{i, {j-\frac{1}{2}}}^{n + \frac{1}{2}}) \end{aligned}$$

(68)

where the numerical fluxes are computed as

$$\begin{aligned} f_{{i+\frac{1}{2}}, j}^{n + \frac{1}{2}} = f\left( u^{{n+\frac{1}{2}}, +x}_{ij}, u^{{n+\frac{1}{2}}, -x}_{i+1,j} \right) , \qquad g_{i, {j+\frac{1}{2}}}^{n + \frac{1}{2}} = g\left( u^{{n+\frac{1}{2}}, +y}_{i,j},u^{{n+\frac{1}{2}}, -y}_{i,j+1} \right) \end{aligned}$$

As in the previous step, the expression (68) is split into a convex combination

$$\begin{aligned} u_{i j}^{n + 1} = k_x \zeta _{ij}^{x} + k_y \zeta _{ij}^{y} \end{aligned}$$

where

$$\begin{aligned} \zeta _{ij}^{x}:= u_{ij}^n - \frac{\Delta t}{k_x \Delta x_i} ( f_{{i+\frac{1}{2}}, j}^{n + \frac{1}{2}} - f_{{i-\frac{1}{2}}, j}^{n + \frac{1}{2}}), \qquad \zeta _{ij}^{y}:= u_{ij}^n - \frac{\Delta t}{k_y \Delta y_j} ( g_{i, {j+\frac{1}{2}}}^{n + \frac{1}{2}} - g_{i, {j-\frac{1}{2}}}^{n + \frac{1}{2}}) \end{aligned}$$

for some $k_x,k_y \ge 0$ with $k_x+k_y=1.$ The admissibility of $\zeta _{ij}^{x}$ and $\zeta _{ij}^{y}$ will imply the admissibility of $u_{i j}^{n + 1}.$ The admissibility of $\zeta _{ij}^{x}, \zeta _{ij}^{y}$ will follow exactly as from the procedure in 1-D (Lemma 3) with appropriate time step restrictions and assumption of admissibility of terms $u^{{n+\frac{1}{2}}, \pm x}_{ij}, u^{{n+\frac{1}{2}}, \pm y}_{ij}, u^{{n+\frac{1}{2}}, *x}_{ij}, u^{{n+\frac{1}{2}}, *y}_{ij}$ for $u^{{n+\frac{1}{2}}, *x}_{ij}, u^{{n+\frac{1}{2}}, *y}_{ij}$ defined as

$$\begin{aligned} \begin{aligned} \mu _{-x}u^{{n+\frac{1}{2}}, -x}_{ij} + u_{ij}^{n+\frac{1}{2},*x} + \mu _{+x}u^{{n+\frac{1}{2}}, +x}_{ij}&= 2 u_{ij}^n \\ \mu _{-y}u^{{n+\frac{1}{2}}, -y}_{ij} + u_{ij}^{n + \frac{1}{2}, *y} + \mu _{+y}u^{{n+\frac{1}{2}}, +y}_{ij}&= 2 u_{ij}^n \end{aligned} \end{aligned}$$

(69)

The precise CFL restrictions and admissibility constraints are in the following Lemma 7.

Lemma 7

Assume that the states $\left\{ u^{{n+\frac{1}{2}}, \pm x}_{ij}, u^{{n+\frac{1}{2}}, \pm y}_{ij}, u^{{n+\frac{1}{2}}, *x}_{ij}, u^{{n+\frac{1}{2}}, *y}_{ij} \right\} _{i,j}$ belong to $\mathcal {U}_{\text {ad}}$, where $u^{{n+\frac{1}{2}}, *x}_{ij}, u^{{n+\frac{1}{2}}, *y}_{ij}$ are defined as in (69). Then, the updated solution $u_{ij}^{n+1}$ of MUSCL-Hancock scheme is in $\mathcal {U}_{\text {ad}}$ under the CFL conditions

$$\begin{aligned} \begin{aligned}&\frac{2\lambda _{x_i}}{\mu _{-x}} \sigma _x \left( u^{{n+\frac{1}{2}}, -x}_{ij}, u^{{n+\frac{1}{2}}, *x}_{ij} \right) \le 1, \quad 2\lambda _{x_i} \sigma _x \left( u^{{n+\frac{1}{2}}, *x}_{ij}, u^{{n+\frac{1}{2}}, +x}_{ij}\right) \\&\quad \le 1, \quad \frac{2\lambda _{x_i}}{\mu _{+x}} \sigma _x \left( u^{{n+\frac{1}{2}}, +x}_{ij}, u^{{n+\frac{1}{2}}, -x}_{i+1,j} \right) \le 1 \\&\frac{2\lambda _{x_i}}{\mu _{-x}} \sigma _x \left( u^{{n+\frac{1}{2}}, +x}_{i-1,j}, u^{{n+\frac{1}{2}}, -x}_{ij} \right) \le 1, \quad 2\lambda _{x_i} \sigma _x \left( u^{{n+\frac{1}{2}}, -x}_{ij}, u^{{n+\frac{1}{2}}, *x}_{ij}\right) \\&\quad \le 1, \quad \frac{2\lambda _{x_i}}{\mu _{+x}} \sigma _x \left( u^{{n+\frac{1}{2}}, *x}_{ij}, u^{{n+\frac{1}{2}}, +x}_{ij} \right) \le 1 \end{aligned} \nonumber \\ \begin{aligned}&\frac{2 \lambda _{y_j}}{\mu _{-y}} \sigma _y \left( u^{{n+\frac{1}{2}}, -y}_{ij}, u^{{n+\frac{1}{2}}, *y}_{ij} \right) \le 1, \quad 2 \lambda _{y_j} \sigma _y \left( u^{{n+\frac{1}{2}}, *y}_{ij}, u^{{n+\frac{1}{2}}, +y}_{ij}\right) \\&\quad \le 1, \quad \frac{2 \lambda _{y_j}}{\mu _{+y}} \sigma _y \left( u^{{n+\frac{1}{2}}, +y}_{ij}, u^{{n+\frac{1}{2}}, -y}_{i,j+1} \right) \le 1 \\&\frac{2 \lambda _{y_j}}{\mu _{-y}} \sigma _y \left( u^{{n+\frac{1}{2}}, +y}_{i,j-1}, u^{{n+\frac{1}{2}}, -y}_{ij} \right) \le 1, \quad 2 \lambda _{y_j} \sigma _y \left( u^{{n+\frac{1}{2}}, -y}_{ij}, u^{{n+\frac{1}{2}}, *y}_{ij}\right) \\&\quad \le 1, \quad \frac{2 \lambda _{y_j}}{\mu _{+y}} \sigma _y \left( u^{{n+\frac{1}{2}}, *y}_{ij}, u^{{n+\frac{1}{2}}, +y}_{ij} \right) \le 1 \end{aligned} \end{aligned}$$

(70)

where $\lambda _{x_i} = \frac{\Delta t}{k_x \Delta x_i}, \lambda _{y_j} = \frac{\Delta t}{k_y \Delta y_j}$ for all i, j.

As in 1-D, we now show that admissibility of $u_{i j}^{n + \frac{1}{2}, *x}, u_{i j}^{n + \frac{1}{2}, *y}$ can also be reduced to admissibility of $u^{*, \pm x}_{i j}, u^{*, \pm y}_{ij}$, similar to Lemma 5. Expanding the definition of $u^{{n+\frac{1}{2}}, *y}_{ij}$ gives us

$$\begin{aligned}{} & {} u_{ij}^{{n+\frac{1}{2}}, *y} = 2 u_{ij}^n -(\mu _{-y}u^{n, -y}_{i j} + \mu _{+y}u^{n, +y}_{ij}) - \frac{\Delta t}{\Delta x_i} \left( f(u^{n, -x}_{i j}) - f(u^{n, +x}_{ij}) \right) \nonumber \\{} & {} \quad - \frac{\Delta t}{\Delta y_j} \left( g(u^{n, -y}_{i j}) - g(u^{n, +y}_{ij}) \right) \end{aligned}$$

(71)

If we obtain the admissibility of

$$\begin{aligned} \mathcal \eta _{ij}^{* {yx}}:= 2 u_{i j}^n - (\mu _{-y}u^{n, -y}_{i j} + \mu _{+y}u^{n, +y}_{i j}) - \frac{\Delta t}{k_x \Delta x_i} \left( f\left( u^{n, -x}_{i j} \right) - f\left( u^{n, +x}_{i j} \right) \right) \end{aligned}$$

(72)

and

$$\begin{aligned} \mathcal \eta _{ij}^{* {yy}}:= 2 u_{i j}^n - (\mu _{-y}u^{n, -y}_{i j} + \mu _{+y}u^{n, +y}_{i j}) - \frac{\Delta t}{k_y \Delta y_j} \left( g\left( u^{n, -y}_{i j} \right) - g\left( u^{n, +y}_{i j} \right) \right) \end{aligned}$$

(73)

for some $k_x,k_y \in [0,1]$ with $k_x+k_y=1,$ then the admissibility of $u_{i j}^{{n+\frac{1}{2}}, *y}$ follows as we can write it as a convex combination

$$\begin{aligned} u_{i j}^{{n+\frac{1}{2}}, *y} = k_x \eta _{ij}^{*{yx}}+ k_y\eta _{ij}^{*{yx}} \end{aligned}$$

and obtain the admissibility of $u_{ij}^{{n+\frac{1}{2}}, *y}$. Thus, we need to limit the slope so that (72, 73) are admissibile. To that end, define $u_{i j}^{**{y x}}, u_{ij}^{**{y y}}$ to satisfy

$$\begin{aligned} \mu _{-x}u^{n, -x}_{i j} + u_{i j}^{**{y x}} + \mu _{+x}u^{n, +x}_{i j} =&2 \left( 2 u_{i j}^n - (\mu _{-y}u^{n, -y}_{i j} + \mu _{+y}u^{n, +y}_{i j}) \right) \\ \mu _{-y}u^{n, -y}_{i j} + u_{i j}^{**{y y}} + \mu _{+y}u^{n, +y}_{i j} =&2 \left( 2 u_{i j}^n - (\mu _{-y}u^{n, -y}_{i j} + \mu _{+y}u^{n, +y}_{i j}) \right) \end{aligned}$$

respectively. Consequently,

$$\begin{aligned} \eta _{ij}^{* {yx}}&= \frac{1}{2} ( \mu _{-x}u^{n, -x}_{i j} + u_{i j}^{**{y x}} + \mu _{+x}u^{n, +x}_{i j})-\frac{\Delta t}{ k_x\Delta x_i}\left( f\left( u^{n, -x}_{i j} \right) - f\left( u^{n, +x}_{i j} \right) \right) \\ \eta _{ij}^{* {yy}}&= \frac{1}{2}( \mu _{-y}u^{n, -y}_{i j} + u_{i j}^{**{y y}} + \mu _{+y}u^{n, +y}_{i j})-\frac{\Delta t}{k_y \Delta y_j} \left( g\left( u^{n, -y}_{i j} \right) - g\left( u^{n, +y}_{i j} \right) \right) \end{aligned}$$

Then, assuming the admissibility of $u_{ij}^{**{y x}}, u_{i j}^{**{y y}}$ and proceeding as in the proof of Lemma 4, we can ensure that $\eta _{ij}^{* {yx}}, \eta _{ij}^{* {yy}} \in \mathcal {U}_{\text {ad}}$ and thus $u_{i j}^{{n+\frac{1}{2}}, *y} \in \mathcal {U}_{\text {ad}}$. Furthermore, since the reconstruction is conservative

$$\begin{aligned} \mu _{-y}u^{n, -y}_{i j} + \mu _{+y}u^{n, +y}_{ij} = \mu _{-x}u^{n, -x}_{i j} + \mu _{+x}u^{n, +x}_{i j} = u_{i j}^n \end{aligned}$$

Thus, admissibility of $u_{ij}^{**{yx}}, u_{ij}^{**{yy}}$ is obtained as

$$\begin{aligned} u_{i j}^{**{yx}} = u_{ij}^{**{yy}} = u_{i j}^n \end{aligned}$$

The arguments for admissibility of $u_{i j}^{{n+\frac{1}{2}}, *x}$ are similar. The admissibility criteria of $u_{i j}^{{n+\frac{1}{2}}, *x}, u_{i j}^{{n+\frac{1}{2}}, *y}$ are summarised in the following lemma.

Lemma 8

Assume that ${u_{ij}^n} \in \mathcal {U}_{\text {ad}}$ and $u^{n, \pm x}_{ij}, u^{n, \pm y}_{ij} \in \mathcal {U}_{\text {ad}}$ for all i, j with conservative reconstruction. Also assume the CFL restrictions

$$\begin{aligned} \begin{aligned} \max _{i,j} \frac{\lambda _{x_i}}{\mu _{-x}} \sigma _x \left( u_{ij}^n, u^{n, -x}_{ij} \right) \le 1, \qquad \max _{i,j} \frac{\lambda _{x_i}}{\mu _{+x}} \sigma _x \left( u^{n, +x}_{ij}, u_{ij}^n \right) \le 1 \\ \max _{i,j} \frac{\lambda _{y_j}}{\mu _{-y}} \sigma _y \left( u_{ij}^n, u^{n, -y}_{ij} \right) \le 1, \qquad \max _{i,j} \frac{\lambda _{y_j}}{\mu _{+y}} \sigma _y \left( u^{n, +y}_{ij}, u_{ij}^n \right) \le 1 \end{aligned} \end{aligned}$$

(74)

where $\lambda _{x_i} = \frac{\Delta t}{k_x \Delta x_i}, \lambda _{y_j} = \frac{\Delta t}{k_y \Delta y_j}$ and $\mu _{\pm x}, \mu _{\pm y}$ are defined in (60). Then, $u_{i j}^{{n+\frac{1}{2}}, *x}, u_{i j}^{{n+\frac{1}{2}}, *y}$ defined in (69) are in $\mathcal {U}_{\text {ad}}$.

Combining Lemmas 6, 7, 8, we will have the 2-D result corresponding to Theorem 3 with the same proof.

Theorem 5

Let $u_{ij}^n \in \mathcal {U}_{\text {ad}}$ for all i, j and $u^{n, \pm x}_{ij}, u^{n, \pm y}_{ij}$ be the conservative reconstructions defined as

$$\begin{aligned} u^{n, \pm x}_{ij} = u_{ij}^n + (x_{{i\pm \frac{1}{2}}} - x_i) \delta ^x_i, \qquad u^{n, \pm y}_{ij} = u_{ij}^n + (y_{j\pm \frac{1}{2}}- y_j) \delta ^y_j \end{aligned}$$

so that $u^{*, \pm x}_{ij}, u^{*, \pm y}_{ij}$ (59) are given by

$$\begin{aligned} u^{*, \pm x}_{ij} = u_{ij}^n + 2(x_{{i\pm \frac{1}{2}}} - x_i) \delta ^x_i,\qquad u^{*, \pm y}_{ij} = u_{ij}^n + 2(y_{j\pm \frac{1}{2}}- y_j) \delta ^y_j \end{aligned}$$

Assume that the slopes $\delta ^x_i, \delta ^y_j$ are chosen to satisfy $u^{*, \pm x}_{ij}, u^{*, \pm y}_{ij} \in \mathcal {U}_{\text {ad}}$ for all i, j and that the CFL restrictions (67, 70, 74) are satisfied. Then the updated solution $u_{ij}^{n+1}$ of MUSCL-Hancock procedure is in $\mathcal {U}_{\text {ad}}$.

B Limiting Numerical Flux in 2-D

Consider the 2-D hyperbolic conservation law (34). The Lax–Wendroff update is

$$\begin{aligned} (u_{ij}^e)^{n+1} = (u_{ij}^e)^n - \Delta t \left[ \frac{1}{\Delta x_e} \frac{\partial F_h^\textbf{e}}{\partial \xi } (\xi _i, \xi _j) + \frac{1}{\Delta y_e} \frac{\partial G_h^\textbf{e}}{\partial \eta }(\xi _i, \xi _j) \right] , \qquad 0 \le i,j \le N \end{aligned}$$

where $F^\textbf{e}_h, G^\textbf{e}_h$ are continuous time-averaged fluxes (5) in the x, y directions for the grid element $\textbf{e}= (e_x, e_y)$. The reader is referred to Section 10 of [5] for more details of the 2-D Lax–Wendroff Flux Reconstruction (LWFR) scheme. Since the 2-D scheme is formed by taking a tensor product of the 1-D scheme, Theorem 2 applies, i.e., the scheme will be admissibility preserving in means (Definition 2) if we choose the blended numerical flux such that the lower order updates are admissible at solution points adjacent to the interfaces. Thus, we now explain the process of constructing the numerical flux where, to minimize storage requirements and memory reads, we will perform the correction within the interface loop where only one of x or y flux will be available in one iteration. Thus theoretical justification for the algorithm comes from breaking 2-D lower order updates into 1-D convex combinations. The general structure of the LWFR Algorithm 1 will remain the same. Here, we justify Algorithm 2 for construction of blended x flux with knowledge of only the x flux. The algorithm for blended y fluxes will be analogous.

We consider the calculation of the blended numerical flux for a corner solution point of the element, as this situation differs from 1-D, due to the fact that a corner solution point is adjacent to both x and y interfaces. In particular, we consider the bottom-left corner point $\textbf{0}= (0,0)$ and show that the procedure in Algorithm 2 ensures admissibility at such points. The same justification applies to other corner and non-corner points. For the element $\textbf{e}= (e_x, e_y)$, denoting interfaces along x, y directions as ${(e_x \pm \frac{1}{2}, e_y)}$, ${(e_x, e_y\pm \frac{1}{2})}$, we consider the update at the bottom left corner $\textbf{0}= (0,0)$, suppressing the local solution point index $i=0$ or $j=0$ when considering the FR interface fluxes. The lower order update is given by

$$\begin{aligned} \tilde{u}_{\textbf{0}}^{n + 1}&= (u_\textbf{0}^{\textbf{e}})^n - \frac{\Delta t}{\Delta x_e w_0} (f_{(\frac{1}{2}, 0)}^\textbf{e}- \tilde{F}_{(e_x- \frac{1}{2}, e_y)}) - \frac{\Delta t}{\Delta y_e w_0} (g_{(0, \frac{1}{2})}^\textbf{e}- \tilde{G}_{(e_x, {e_y-\frac{1}{2}})}) \end{aligned}$$

where $\tilde{F}_{(e_x- \frac{1}{2}, e_y)}, \tilde{G}_{(e_x, {e_y-\frac{1}{2}})}$ are heuristically guessed candidates for the blended numerical flux (12). Pick $k_x, k_y > 0$ such that $k_x + k_y = 1$ and

$$\begin{aligned} \begin{aligned} \tilde{u}_x^{\text {low}, n + 1}&= (u_\textbf{0}^{\textbf{e}})^n - \frac{\Delta t}{k_x \Delta x_e w_0} \left( f^\textbf{e}_{\left( \frac{1}{2}, 0\right) } - f_{\left( e_x- \frac{1}{2}, e_y\right) } \right) \\ \tilde{u}_y^{\text {low}, n + 1}&= (u_\textbf{0}^{\textbf{e}})^n - \frac{\Delta t}{k_y \Delta y_e w_0} \left( g^\textbf{e}_{\left( 0, \frac{1}{2}\right) } - g_{\left( e_x, {e_y-\frac{1}{2}}\right) }\right) \end{aligned} \end{aligned}$$

(75)

satisfy

$$\begin{aligned} \tilde{u}_x^{\text {low}, n + 1}, \tilde{u}_y^{\text {low}, n + 1} \in \mathcal {U}_{\text {ad}}\end{aligned}$$

(76)

Such $k_x, k_y$ exist because the lower order scheme with lower order flux at element interfaces is admissibility preserving. The choice of $k_x, k_y$ should be made so that (76) is satisfied with the least time step restriction, but we have found the Fourier stability restriction imposed by (32) to be sufficient even with the most trivial choice of $k_x = k_y = \frac{1}{2}$. The discussion of literature for the optimal choice of $k_x,k_y$ is the same as the one made for the 2-D MUSCL Hancock scheme (64) and is not repeated here. After the choice of $k_x, k_y$ is made, if we repeat the same procedure as in the 1-D case, we can perform slope limiting to find $F_{e_x - \frac{1}{2}, e_y}$, $F_{e_x, e_y - \frac{1}{2}}$ such that

$$\begin{aligned} \tilde{u}_x^{n + 1}&= (u_\textbf{0}^{\textbf{e}})^n - \frac{\Delta t}{k_x \Delta x_e w_0} \left( f_{\left( \frac{1}{2},0\right) }^\textbf{e} - F_{\left( e_x - \frac{1}{2}, e_y\right) }\right) \end{aligned}$$

(77)

$$\begin{aligned} \tilde{u}_y^{n + 1}&= (u_\textbf{0}^{\textbf{e}})^n - \frac{\Delta t}{k_y \Delta y_e w_0} \left( g^\textbf{e}_{\left( 0, \frac{1}{2}\right) } - G_{\left( e_x, e_y - \frac{1}{2}\right) }\right) \end{aligned}$$

(78)

are also in the admissible region. Then, we will get

$$\begin{aligned} k_x \tilde{u}_x^{n+1} + k_y \tilde{u}_y^{n+1} = \tilde{u}_{\textbf{0}}^{n+1} \end{aligned}$$

(79)

We now justify Algorithm 2 as follows. Algorithm 2 corrects the numerical fluxes during the loop over x interfaces to enforce admissibility of $\tilde{u}_x^{n + 1}$ (77) at all solution points neighbouring x interfaces including the corner solution points, and the analogous algorithm for y interfaces will ensure admissibility of $\tilde{u}_y^{n + 1}$ (78) at all solution points neighbouring y interfaces including the corner points. At the end of the loop over interfaces, (79) will ensure that lower order updates at all solutions points neighbouring interfaces are admissible and Algorithm 1 will be an admissibility preserving Lax–Wendroff scheme for 2-D if we compute the blended numerical fluxes $F_{(e_x+\frac{1}{2}, e_y)}, F_{(e_x, e_y+\frac{1}{2})}$ using Algorithm 2 and its counterpart in the y direction.

Algorithm 2

Computation of blended flux $F_{e_x+\frac{1}{2}, e_y, j}$ where $(e_x+\frac{1}{2}, e_y)$ are the interface indices and $j \in \{0,\dots , N\}$ is the solution point index on the interface