Accuracy analysis for explicit-implicit finite volume schemes on cut cell meshes

The solution of time-dependent hyperbolic conservation laws on cut cell meshes causes the small cell problem: standard schemes are not stable on the arbitrarily small cut cells if an explicit time stepping scheme is used and the time step size is chosen based on the size of the background cells. In [J. Sci. Comput. 71, 919-943 (2017)], the mixed explicit implicit approach in general and MUSCL-Trap in particular have been introduced to solve this problem by using implicit time stepping on the cut cells. Theoretical and numerical results have indicated that this might lead to a loss in accuracy when switching between the explicit and implicit time stepping. In this contribution we examine this in more detail and will prove in one dimension that the specific combination MUSCL-Trap of an explicit second-order and an implicit second-order scheme results in a fully second-order mixed scheme. As this result is unlikely to hold in two dimensions, we also introduce two new versions of mixed explicit implicit schemes based on exchanging the explicit scheme. We present numerical tests in two dimensions where we compare the new versions with the original MUSCL-Trap scheme.


Introduction
Cartesian embedded boundary meshes for computing flow problems involving complex geometries have become very popular as mesh generation is fully automatic and fairly cheap: the geometry is simply cut out of Cartesian background cells.This results in cut cells where the object intersects the background mesh.Cut cells can have various shapes and can in particular be arbitrarily small.This causes various issues when standard methods are used to solve partial differential equations (PDEs) on these meshes.The specific issue depends on the kind of PDE to be solved.
In the context of time-dependent hyperbolic conservation laws the main issue is what is referred to as the small cell problem: Typically, explicit time stepping schemes are used.By the CFL condition the time step size is coupled to the size of the cells.One wants to choose the time step based on the size of the larger background cells and use the same time step on the arbitrarily small cut cells.For standard schemes, this approach causes the values on the small cut cells to explode.
One approach to solve this issue is to use cell merging or cell agglomeration, see, e.g., [26,38,39].In that approach cut cells that are too small are merged / combined with neighboring cells, which eliminates the problem.The downside is that the complexity is moved back into mesh generation process.The alternative is to develop algorithmic solution approaches to overcome the small cell problem.Two very well-established approaches are the flux redistribution method [11,13] and the h-box method [8,23,7].More recent approaches include the dimensionally split approach [25,21], the extension of the active flux method to cut cells [24], and the state redistribution scheme (SRD) [6].All of these schemes are based on finite volume approaches (if one counts the active flux method as a finite volume scheme).A few years ago, the development of algorithmic solution approaches in the context of discontinuous Galerkin (DG) methods has started.Existing solution approaches here are the DoD stabilization [15,34], the extension of the ghost penalty stabilization [9] to time-dependent first-order hyperbolic problems [19,18], as well as the extension of the SRD scheme to a DG setting [20].
In this contribution we will consider the mixed explicit implicit scheme in more detail, which was introduced by May and Berger [29,31,33] in the context of a finite volume setting to overcome the small cell problem.The idea is quite simple: cut cells are treated implicitly for stability but cells away from the cut cells use a standard explicit time stepping scheme to keep the cost low.The switch happens by means of flux bounding.The authors combined a second-order explicit scheme with a second-order implicit scheme.Numerical experiments [33] have shown that the resulting mixed scheme converges with second order in the L 1 norm.In the L ∞ norm, one numerically sees full second order in one dimension but in two and three dimensions orders between 1 and 2 were observed.
In this contribution we want to examine the error behavior of the mixed explicit implicit scheme in more detail and also test possible cures.In one dimension, we will first analyze the one step error and then show that the mixed scheme indeed converges with second order.As numerical results indicate that we will not be able to prove the same result in higher dimensions, we will then introduce two new versions of the mixed explicit implicit scheme with an improved transition error.These will be designed so that on a Cartesian mesh the transition between explicit and implicit time stepping has a third-order one step error.We will include results from a test involving cut cells to see whether this increases the overall accuracy of the scheme.Some of the material presented here has been part of the master thesis of Fabian Laakmann [27].
The idea of combining explicit and implicit time stepping to gain stability has been used by other authors as well and has in particular become more popular in recent years, see, e.g., [14,36,37,17] and the references cited therein.Here, we will focus on the approach as introduced in [29,33] in the context of cut cells.We note though that this approach can easily be extended to other problem settings as well.Recently, it has been applied in the context of two phase flow problems: if one uses a sharp interface method, then the creation of new phases (nucleation or cavitation) results in the creation of new tiny cells.In [35], the newly created tiny cells have successfully been treated with the first-order version of the mixed explicit implicit scheme, extended to the isothermal Euler equations.
This contribution is structured as follows: in section 2, we focus on the situation in one dimension.We introduce the mixed explicit implicit scheme MUSCL-Trap as used in [33] and examine its error.We then present two different variants of the mixed explicit implicit scheme, which have a reduced transition error compared to the original scheme.In section 3, we extend the considerations to two dimensions.We first formulate the new variants in two dimensions and then compare the three schemes numerically.We conclude with an outlook in section 4. x Figure 1: 1d model problem: Equidistant grid with one small cell of length αh labeled as cell 0.
2 Mixed explicit implicit schemes in 1d We consider the linear advection equation with initial data s(t 0 , •) = s 0 .The standard model mesh for developing cut cell schemes in 1d is shown in figure 1: an equidistant mesh with mesh width h contains one cell of length αh, labeled as cell 0, in the middle.Here, α ∈ (0, 1] denotes the volume fraction -the ratio of the small cell to full cell volume.

The MUSCL-Trap scheme
The idea behind mixed explicit implicit schemes for cut cell meshes as used in [29,31,33] is to employ an implicit scheme on the cut cell 0 for stability.Away from the cut cell an explicit scheme is used to keep the cost low.In [33], May and Berger developed what we will refer to as MUSCL-Trap, the combination of the explicit MUSCL scheme and the implicit Trapezoidal scheme by means of flux bounding.MUSCL (Monotonic Upwind Scheme for Conservation Laws) [40,12] is a common explicit finite volume scheme.The scheme is second-order accurate in space and time.In one dimension, the MUSCL scheme for the linear advection equation (1) on an equidistant grid is given by with S n i,x ≈ ∂ x s(t n , x i ), computed using central difference gradients and a standard slope limiter [28], with CFL number λ = u∆t h .The scheme is stable for 0 < λ ≤ 1 and would therefore become unstable on the cell 0 for α < 1.
The implicit Trapezoidal rule in time in combination with a suitable slope reconstruction in space is given by for a non-equidistant mesh with cells of length h i .
In [29,33], May and Berger examined various ways for combining an explicit scheme with an implicit scheme and found flux bounding to be the best way.The idea is illustrated in Step 1: Update of explicit cells.
Step 2: Update of cut cell and transition cells Figure 2: Flux bounding: First the cells away from the cut cell are updated using an explicit scheme (edges marked in green).Then the neighborhood of the cut cell is updated.The Cartesian neighbors of the cut cell are transition cells.The update on these cells employs both explicit fluxes (edges marked in green) and implicit fluxes (edges marked in light blue), see also [33], [30,Fig. 6].
figure 2: One first updates all Cartesian cells that are not direct neighbors of the cut cell using the explicit scheme, see figure 2(a).In a second step, one updates the value on cells −1, 0, 1, compare figure 2(b).Here, cell 0 is treated fully implicitly to guarantee stability on the small cut cell, indicated by the light blue edges.The cells −1 and 1 are transition cells: When updating cell −2 in Step 1 with the fully explicit scheme, one assumed that a certain amount of mass, represented by , has left cell −2.For the full scheme to be conservative, one needs to ensure that the mass arrives in cell −1.This is achieved by re-using the explicit flux for the edge −3/2 in Step 2. Note that transition cells employ both explicit and implicit fluxes.Cut cells are treated fully implicitly, and Cartesian cells that are not neighbors of cut cells are treated fully explicitly.Flux bounding is set up in a symmetric manner, i.e., it does not distinguish between u > 0 and u < 0, allowing the extension to more general problems, such as Euler equations [35].Flux bounding has the following properties: 1.It is conservative (by construction).
2. It satisfies a TVD property in a suitable setting.May and Berger [33,Thm 1] showed TVD stability for a mixed scheme combining MUSCL (with minmod limiter [28]) as explicit scheme with implicit Euler and piecewise constant data as implicit scheme for 0 < λ ≤ 1, independent of the size of α ∈ (0, 1]. 3. It has a natural extension to 2d and 3d (see section 3 as well as [33]).
This results in the following formulae for MUSCL-Trap in 1d: For the slope reconstruction on cells i ≤ −2 and i ≥ 2 unlimited central differences are used.
On the cut cell 0 and the transition cells −1 and 1 a least squares approach is employed [5].This extends in a straight-forward way to higher dimensions [4,32].Note that the least squares slope S 1,x enters the computation of the MUSCL flux . Except for this special case, only central difference slopes are used for the computation of MUSCL fluxes.Throughout this paper we will assume all slopes to be unlimited.
Remark 1.This contribution focuses on the accuracy of mixed schemes in general and the accuracy of MUSCL-Trap in particular.We will not discuss limiting here.For limiting on cut cells, see, e.g., [32].

The one step error of MUSCL-Trap
The one step error of scheme (4) has roughly been analyzed in [33,30].We will do this here more thoroughly.Let the MUSCL-Trap scheme be given by , depending on i.For the computation of the one step error, we replace the input arguments of Φ i with the true cell averages at time t n and t n+1 , given by sn i and sn+1 i , resulting in Remark 2. The one step error L(s n , sn+1 ) i measures the error made in one time step t n → t n+1 on cell i.Typically, i.e., on uniform meshes, the order of convergence of the one step error is one order higher than the overall order of the scheme.
The numerically observed convergence order of the one step error for a smooth test problem is summarized in figure 3. We observe a third-order one step error on cells that are sufficiently far away from the cut cell.These cells use the full MUSCL scheme on a uniform mesh.In the neighborhood of the cut cell, we observe a second-order one step error.This is due to the following four error sources.
Error sources 2.1.(1) The switch in the time stepping scheme between MUSCL and Trapezoidal rule results in an error term of size O(h 2 ).The examination of its propagation is one key aspect of this research.
(2) The irregularity of the cut cell causes several errors if α < 1: (a) the slope reconstruction, which uses a least squares approach, is only of first order, resulting in second-order error terms that we summarize in a i h 2 for cell i; (this is the only reason for the error of size O(h 2 ) on cell 2;) (b) in the computation of the fluxes, we use the variable S n i (which is supposed to be an approximation to the cell average sn i ) as approximation to the point value s(t n , x i ); this results in second-order error terms that we summarize in b i h 2 ; (c) for a third-order one step error, we would need to use a quadratic reconstruction in space instead of a linear one; we summarize the resulting error terms on cell i in c i h 2 .
We note that on an equidistant mesh, i.e., α = 1, all error terms listed under (2) reduce to third-order errors, partly due to cancellation effects with neighboring cells.Using this notation, we observe on cells −1, 0, and 1 the following errors: Note that the error that we make when switching from MUSCL to Trapezoidal on cell −1 is to leading order the negative of the error that we make when switching back from Trapezoidal to MUSCL on cell 1.

Numerical results for MUSCL-Trap
In [33], the authors presented numerical results for the model problem shown in figure 1, which we will start with as well to keep the presentation mostly self-contained., α = 10 −4 , and final time is T = (1 + αh) so that the test function is back to its original position.The L 1 error has been normalized to account for the changing domain length.We solve s t + s x = 0 with λ = 0.8, independent of α.
The results for the one step error and the error at time T are shown in Table 1.
For the one step error in the L ∞ norm we observe second-order convergence, as expected.The third-order convergence in the L 1 norm can be explained by a simple counting argument: let M denote the number of Cartesian cells.Then, there are M -3 cells (out of M +1 total cells) with an error of size O(h 3 ) and only 4 cells with an error of size O(h 2 ).This results in the third-order convergence in the L 1 norm.The results for the error at time T are more surprising: we observe full second-order convergence in the L 1 and L ∞ norm despite the one step error only converging with second order in L ∞ .However, there is a reasonable explanation for this behavior.Let us consider error accumulation from a Lagrangian view point: let us fix one cell, e.g., the cell that contains the peak of the sine curve, put a tracer there, and follow the tracer during the simulation.If we assume that our scheme is stable (and that error propagates with norm smaller/equal than 1), the error at time T essentially corresponds to the sum of all the one step errors of all the cells that our tracer has passed.Given that most cells have a one step error of third order and that there are only 4 isolated cells with a one step error of second order, this consideration actually implies to expect a second-order error on all cells at time T .This is consistent with the numerical results.
This consideration motivates the following new test, Test 2, which we expect to be more challenging: Instead of using only one cut cell, we double the number of cut cells as we halve the grid size h.Then, the number of cells with one step error O(h 2 ) increase with O( 1 h ) and should affect the convergence order at time T .
Test 2: The test setup is shown in figure 4. We use K = 40 and L = 4 on the coarsest mesh.The overall domain length is given by [0, 1 + 4αh 0 ] with h 0 denoting the coarsest mesh (in our case h 0 = 1 160 ).Different to Test 1, the domain length stays constant.The test function is sin 2π(x+0.36)1+4αh 0 , α = 10 −4 , and T = (1 + 4αh 0 ) so that the test function is back to its original position.The results for the one step error and the error at time T are shown in Table 2.
Despite the one step error being only of second order in both the L 1 and L ∞ norm, the error after one period also converges with second order both in L 1 and L ∞ .This behavior ) and one cut cell of length αH in the middle.When the mesh is refined by a factor of 2, we keep K (and α) fixed but double L. Analogous to figure 1, we refer to cut cells as (cells of type) 0, to the left Cartesian neighbors of cut cells as (cells of type) -1, and to the right Cartesian neighbor as (cells of type) 1. of the error accumulation behaving differently on non-uniform meshes has been observed and analyzed before [41].We will follow the proof from [41] to show second order accuracy for MUSCL-Trap in theorem 1 below.Similar ideas have also been used to show the second order accuracy of the h-box method [8].
Different to Test 1, we observe for the one step error measured in the L 1 norm, a convergence order of 2 (compared to 3 before).This can be explained as follows: let e i = S n i − sn i denote the error on cell i, let K+1 be the number of cells per block and let there be L blocks.As before, M denotes the total number of Cartesian cells.This implies M + L = L • (K + 1).Then, the L 1 error can be computed as (with c denoting a generic constant) Therefore, for keeping K (the number of cells per block) fixed, we get an L 1 error of O(h 2 ).This theoretical consideration coincides with our numerical observation.Note that in Test 1, K increased with O( 1 h ), explaining the better convergence rate of 3 in the L 1 norm for the one step error.

Proof for MUSCL-Trap being second-order accurate
For the proof of second-order convergence, we require stability of the method.Examining this stability mathematically is very challenging due to using a non-uniform mesh and a mixed scheme.Applying, for example, von Neumann stability analysis is not feasible.We therefore formulate this as assumption 2.1 in the proof below.In our numerical tests, we found that this was true for the L 1 and L ∞ norm.
Theorem 1.Let Ass.2.1 below hold true.Consider the MUSCL-Trap scheme given by (4), but assume that slopes are computed by means of an unlimited forward difference quotient.Let 0 < λ ≤ 1 independent of α.Then the scheme is second-order accurate with respect to the norm used in Ass.2.1 for the linear advection equation for model problem 4 for sufficiently smooth initial data s 0 .
Proof.For the proof, we use an idea that goes back to Wendroff and White [41], and which has also been used to show second-order accuracy of the h-box method [8]: Assume that we are able to construct a grid function with cell averages wn i such that (i) the new grid function wn i is sufficiently close to sn i and (ii) the one step error is of third order for all cells.To be precise, wn i is supposed to satisfy Assumption 2.1.We assume that the MUSCL-Trap scheme is stable with respect to the norm ∥•∥ in the following sense: there holds for the error propagation.
For the error grid function E n = S n − wn there holds by linearity of Φ From Ass. 2.1 and assumption (8b), we can conclude with N = T ∆t and ∆t = O(h).In other words: we have 'normal' error accumulation with respect to the new solution wn .Together with property (8a) we get by means of the triangle inequality which implies global, second-order convergence with respect to the true solution s.
It remains to find such a suitable grid function wn .We note that this is the essence of this proof.
Case α = 1: We first consider the case of a uniform mesh.The main error source that we like to examine is the switch in time stepping.Note though that different to (7), we now assume that forward differences are used.Therefore, the actual error terms look slightly different and are given by We define A direct computation then shows L(γ n , γ n+1 ) 0 = O(h 3 ) as well as This implies by linearity L( wn , wn+1 Therefore, assumptions (8a) and (8b) are satisfied.
Case α < 1: In this case we need to address all four error sources.First, a direct computation shows that for forward differences there holds with Note in particular that for this setup, the one step error on the small cell is surprisingly of third order as errors cancel nicely.We define with This results in L(γ n , γ n+1 ) 0 = O(h 3 ) as well as in which implies the claim.
Remark 3. The advantage of forward differences compared to central differences lies in the reduced coupling of cells.For central differences, we have not been able to find suitable coefficients γ n when using periodic boundary conditions.But the numerical results in section 2.1.2imply that the same order of convergence holds true when one employs central differences instead of forward ones.
Remark 4. We note that we addressed two different error sources in this proof: the error caused by the switch in the time stepping scheme and the error caused by the irregular length of the cut cell.The reason that the former error does not accumulate is simply that we make an error of size − λ 3 4 h 2 s xx (t n , x −1 ) on cell −1 and an error of size λ 3 4 h 2 s xx (t n , x 1 ) on cell 1.In other words: the error that we made on cell −1 cancels to leading order with the error that we make two cells later.
Finally, we verify our results numerically.We repeat Test 2 using forward difference quotients but this time we compare our discrete solution S n i to our new grid function w given by (11) and (12).The result is given in Table 3 and shows the expected convergence rates.

New variants of mixed explicit implicit schemes
Our theoretical and numerical considerations above imply that the mixed scheme MUSCL-Trap is second-order accurate in 1d, despite only having a second-order one step error.The numerical results in 2d though for MUSCL-Trap presented in [33] showed convergence rates between 1.3 and 1.6 measured in the L ∞ norm.This implies that in 2d the error does accumulate to some extent.Therefore, we cannot expect to be able to find a new grid function wn with properties (i) and (ii) in 2d.
In order to find a scheme with a third-order one step error, we would need to address all four error sources.As the error sources 2(a)-(c) are very difficult to examine in 2d (due to varying sizes and shapes of cut cells) and are generally shared by most cut cell schemes, we focus here on the error caused by switching the time stepping scheme.Further, we expect the error for switching the time stepping scheme to have a very different effect in 1d and 2d: in 1d, characteristics are typically set up in a way that an error caused by switching from explicit to implicit is followed very briefly afterwards by an error caused by switching from implicit back to explicit.In 2d however, see the setup of the ramp test below in figure 6, this effect of switching back and forth might not be there.Instead there can be characteristics that mainly go through transition cells.

Developing mixed time stepping schemes with better transition error
We will discuss the issue of reducing the error by switching the time stepping scheme in a sightly more general setting.For finite volume schemes, one typically uses the conservative update formula with F i±1/2 being an appropriate approximation to the flux at the edge x i±1/2 during the time step.In the case of the linear advection equation (1) (with WLOG u = 1), we try to approximate To only focus on the temporal error, we will reinterpret time stepping schemes as quadrature formulae for approximating the integral in (14), assuming that we know all values at the edges x = x i±1/2 that we need.To give an example: using explicit Euler in time would correspond to approximating F true i+1/2 by the point evaluation s(t n , x i+1/2 ).For MUSCL and Trapezoidal time stepping we have the following results: • MUSCL: The corresponding time stepping scheme is the explicit midpoint rule, which corresponds to the following approximation of F true i+1/2 : This results in ).
• Trapezoidal: The Trapezoidal time stepping scheme approximates F true i+1/2 by This results in the error Therefore, assuming that the reconstruction in space is done sufficiently accurately, when using Trapezoidal rule time stepping, we approximate F true i+1/2 with an error of O(∆t 2 ).However, when we use Trapezoidal on both edges, we make to leading order the same (systematic) error for the approximation of F true i−1/2 .When applying the update formula ( 13), the leading order error term cancels and there holds Therefore, for h = O(∆t), we get an error of order O(∆t 3 ) (for the time approximation).However, if we use the explicit midpoint rule ( = MUSCL) for F i−1/2 and Trapezoidal time stepping for F i+1/2 , we get This is the transition error that we saw before.
Note that these considerations have another important implication: the seemingly easiest way to improve the transition error that we observed for MUSCL-Trap would be to just keep MUSCL and use a third-order implicit scheme.Or to keep the Trapezoidal time stepping scheme and upgrade the explicit time stepping scheme to third order.However, both versions would not show the result one hoped for.The easiest way to see this is by combining, for example, the MUSCL flux F (M ) i−1/2 with the true flux F true i+1/2 .Then, when applying the conservative update formula (13), we would be left with an error of O(∆t 2 ) despite the scheme being apparently more accurate than when applying MUSCL fluxes on both edges.
Therefore, the two options to get a scheme with a third-order one step error (with respect to the time error) are to (a) either use both an explicit and an implicit time stepping scheme that approximate F true i±1/2 with third order or (b) to find a combination of a second-order explicit and a second-order implicit time stepping scheme such that the leading order error terms match.
Besides approach (a) being more expensive, third-order time stepping schemes also typically involve several stages, which makes the coupling very complicated.We therefore follow approach (b): we fix the implicit Trapezoidal rule and look for a new explicit time stepping scheme.We will consider two different modifications of the MUSCL scheme.

The MUSCLmod-Trap scheme
The first approach is to change the MUSCL scheme given in (2) to use The idea is to account for the error term ∆t 4 s tt (t n , x 1+1/2 ) that we saw above when switching between Trapezoidal and explicit midpoint rule but the details slightly vary as we generally do not have the exact values on the edges.The precise formulation simply comes out of the error analysis for combining MUSCL with Trapezoidal.The truncation error analysis shows that • using MUSCLmod on an equidistant mesh as fully explicit scheme, there holds L(s n+1 i ) = O(h 3 ); • using MUSCLmod in a mixed MUSCLmod-Trap scheme (with α = 1), there holds A von Neumann stability analysis for the new MUSCL version shows a CFL condition of 0 < λ ≤ 1.
In table 4 we show the result for the new mixed MUSCLmod-Trap scheme for Test 2, but using α = 1, i.e., we switch the time stepping between explicit and implicit at the corresponding cells but do that on an equidistant mesh to avoid the error sources 2(a)-(c).As expected we observe third order for the one step error and second error for the error at time T .We use a standard second-order accurate difference quotient to evaluate the second derivatives.

The MPRKC-Trap scheme
For the second variant, consider the explicit Trapezoidal rule, also known as standard secondorder SSP RK (strong stability preserving Runge-Kutta) scheme [22], and given by (for the ODE y t = g(y(t))): ) .
If we compare that with the implicit Trapezoidal rule, we see that there is a O(∆t 2 ) difference due to using a second-order predictor step y (1) instead of y n+1 .Therefore, to get an explicit scheme that matches the leading order error term of i+1/2 up to third order, we need to replace the predictor step (computation of y (1) ), which currently uses explicit Euler, by a more accurate predictor.We choose the explicit midpoint rule for this purpose.Together with a suitable space discretization (the new predictor step will correspond to the MUSCL scheme), this results in a new explicit finite volume scheme, which we will present now.Note that by construction this new scheme will show a third-order one step transition error when coupled to (implicit) Trapezoidal time stepping with slope reconstruction.
The new explicit scheme, which we call MPRKC (for MUSCL Predictor RK Corrector) scheme, is given by with Remark 5.In terms of cost, the new scheme is roughly twice as expensive as MUSCL for taking one time step.Compared to using the standard second-order SSP RK scheme with slope reconstruction in space, the cost is roughly the same.We just exploit the information from the slope reconstruction in stage 1 more carefully.
Proof.The claim follows by means of a von Neumann stability analysis.We analytically deduce a very lengthy and complicated expression for the amplification factor |G| (not given here).
We combine the explicit MPRKC scheme with the implicit Trapezoidal scheme (3) using flux bounding, resulting in MPRKC-Trap.This is sketched in figure 5. Due to MPRKC being a two-stage scheme, the implicit region has become bigger compared to MUSCL-Trap.The results for the mixed scheme on an equidistant mesh with α = 1 but with switching the scheme as indicated in Test 2 are shown in table 5.As expected, the one step error converges with third order and the error at time T with second order.
Figure 5: The new MPRKC-Trap scheme: We compute S (1) i using the MUSCL scheme for cells i ≤ −2 and i ≥ 2 as sketched in figure 2(a).Then we compute the explicit flux F n+1/2,ET given by ( 16) for all edges marked in green and update the fully explicitly treated cells i ≤ −4 and i ≥ 4. We combine this using flux bounding with an implicit treatment (edges marked in light blue) of cells −3, . . ., 3. Note that now cells −3 and 3 correspond to the transition cells.3 Mixed explicit implicit schemes in 2d In the following we will discuss mixed explicit implicit schemes in 2d: we will first briefly introduce the 2d version of the MUSCL-Trap scheme as presented in [33] and examine its error.Then, we will introduce variants of MUSCLmod-Trap and MPRKC-Trap in 2d.We will conclude with a comparison of the resulting schemes for both a fully Cartesian mesh and a cut cell mesh.

The MUSCL-Trap scheme
In 2d, we solve the linear advection equation Here, s(t, x, y) denotes a scalar field.For simplicity, we will only consider a constant velocity field, i.e., u and v are assumed to be constant.This is sufficient to show our main findings.
As explicit scheme, an unsplit, two-dimensional MUSCL scheme [3] is used.For brevity reasons, we will not present the details here.We note though that this MUSCL scheme uses corner-coupling to capture proper dependencies in 2d and is therefore stable under the CFL condition 0 < ν ≤ 1 in combination with the time step computation Flux bounding is used to couple the explicit MUSCL scheme with an implicit scheme.The idea is sketched in the left graphic of figure 6: Full Cartesian cells that are edge neighbors of cut cells are marked as transition cells.Cut cells are treated fully implicit.Therefore, fluxes between cut cells and transition cells use an implicit scheme.All remaining fluxes, in particular fluxes between two transition cells, use an explicit scheme.
This results in the following two-step algorithm, which is analogous to the 1d situation: Given S n ij , (i) compute all explicit fluxes using the MUSCL scheme and update all fully explicitly treated cells to S n+1 ij ; (ii) compute all implicit fluxes and update cut cells and transition cells to S n+1 ij .
As implicit scheme, Trapezoidal time stepping in combination with a least squares formulation for slope reconstruction in space [4,32] is used.On Cartesian cells that are not transition cells standard central differences are used for slope reconstruction.No limiting is applied.
We note that on cut cells, one reconstructs to the midpoint of the cut cell edges, not to the midpoint of the Cartesian edges.To be more precise, the update on a cut cell is given by where α ij ∈ (0, 1) is the volume fraction and Here (x i+1/2,j , y i+1/2,j ) denotes the location of the edge midpoint of face (i + 1/2, j), (x ij , y ij ) denotes the centroid of cut cell (i, j), and S ij,x and S ij,y refer to the reconstructed unlimited x-and y-slopes respectively in cell (i, j).Further, β i+1/2,j ∈ [0, 1] represents the area fraction of cut cell edge (i + 1/2, j) compared to a full Cartesian edge.The fluxes G n+1/2,T i,j+1/2 are defined analogously.This concludes the brief description of the MUSCL-Trap scheme.More details can be found in [33].
Remark 6 (Cost of MUSCL-Trap).In each time step, one needs to solve an implicit system that couples the cut cells and transitions cells.Since cut cells occur only at the boundary of the embedded object, the size of this system is one dimension lower than the overall number of cells.In 2d on a Cartesian grid with N cells in each direction, one expects O(N ) unknowns in the implicit system.Therefore, using this mixed explicit implicit approach is significantly cheaper than using an implicit scheme everywhere.
In [33], May and Berger presented numerical results for the (unlimited) MUSCL-Trap scheme in 2d.The L 1 error converged with second order.The L ∞ error however showed convergence rates between 1.3 and 1.6.This is not in line with the second-order accuracy that we saw in 1d.Newer results with DG codes on cut cell meshes for piecewise linear polynomials however have also showed reduced convergence rates in the L ∞ norm of 1.5 to 1.6 [15,20].This raises the general question of whether we can expect to see full second-order convergence in the L ∞ norm on 2d cut cell meshes at all but this is the goal.
In the following we want to examine the accuracy of the scheme in 2d more closely than done in [33].Here, we will focus on the error caused by switching the time stepping scheme.Due to the complexity of cut cells in 2d, our examination will be mostly based on numerical experiments.
Remark 7.For our numerical tests in 2d, we use the code setup from [33] as starting point: our implementation is based on BoxLib [16], a library for massively parallel AMR applications.For the generation of the cut cells, we use patchCubes, a variant of the cubes mesh generator that is part of the Cart3D package [10,2].The solution of the resulting implicit system is done using umfpack [1].

The transition error for MUSCL-Trap (in absence of cut cells)
We would like to examine the time stepping transition error in the absence of error sources (2)(a)-(c) caused by the irregularity of the cut cells.Therefore, we consider the test setup shown in figure 7: we create a setup where we switch between explicit and implicit time stepping but all cells are full Cartesian cells.As on a Cartesian mesh error sources (2)(a)-(c) drop out, the remaining error should be dominated by the switch from MUSCL to Trapezoidal.
For the theoretical considerations, we focus on the case of 45 • and choose u = v = 1 and ∆x = ∆y.The update on a transition cell is given by A lengthy computation, which involves the usage of the equation s t + s x + s y = 0 to re-express derivatives, leads to the following formula for the one step error on transition cells: To do so, we intersect a Cartesian grid with a line with angle β.All cells cut by this line (marked in pink) are treated as 'fake' (as Cartesian) cut cells.All 'fake' cut cell edges as well as all edges of cells below the 'fake' cut cells use the implicit scheme (marked in red).All remaining edges (marked in blue), i.e., all edges of cells above the 'fake' cut cells use the explicit scheme.The switch takes place on the cells marked with a 'T'.
Therefore, we expect in general a one step error of second order.For the special case s xx = s yy , we expect a third-order one step error.(We confirmed this in numerical experiments not shown here.)In numerical tests, see e.g.Test 3 below, we often observed (for a variety of angles) that on a full Cartesian mesh, we have a second-order one step error in L ∞ due to this transitionbut that this transition error accumulates 'only' to a scheme of order O(h 1.5 ) for the L ∞ error at time T .

New mixed explicit implicit schemes
We now present the extensions of MUSCLmod-Trap and MPRKC-Trap to two dimensions.

The MUSCLmod-Trap scheme
To derive the new explicit MUSCLmod scheme, we essentially do a truncation error analysis of Trapezoidal rule and MUSCL scheme on a mesh as shown in figure 7 and design MUSCLmod to contain the second-order transition error terms.One needs to be a bit careful in this computation as the original MUSCL scheme is based on an unsplit version, involving corner coupling and a more evolved evaluation of transverse derivatives.This then results in the following change of the reconstructed value S M,n+1/2 i+1/2,j (the approximation to the solution at the midpoint of edge (i + 1/2, j) at time t n+1/2 used in the MUSCL scheme) for the case u, v > 0 The second derivatives are computed using standard second-order difference quotients on Cartesian cells away from cut cells.On transition cells we fit a quadratic polynomial to compute them.
Numerical tests show for MUSCLmod as fully explicit scheme on a Cartesian mesh a thirdorder one step error and standard second-order accuracy at time T .The numerical stability limit also seems to be very similar to the original unsplit MUSCL.

The MPRKC-Trap scheme
We now extend the explicit MPRKC scheme, given by ( 15)-( 16) in 1d, to 2d.We use the unsplit MUSCL scheme described by Almgren et al. [3], which we also use for the MUSCL-Trap scheme, to compute the predictor S (1) .Afterwards, we use the two-dimensional analogue of ( 16) for the computation of the fluxes F n+1/2,ET i±1/2,j and G n+1/2,ET i,j±1/2 .The classic two-stage second-order SSP RK scheme in combination with central differences for ∆x = ∆y = h is stable with ν ≤ 1.0 under the CFL condition [7] due to using a split approach.The new MPRKC scheme uses the unsplit MUSCL as predictor (instead of the split upwind scheme) but the split explicit Trapezoidal scheme as corrector.
We therefore expect MPRKC to be at least as stable as the classic two-stage second-order SSP RK scheme and to be in particular stable under the CFL condition (20).We confirmed this in numerical tests.Next, we briefly sketch the extension of the mixed MPRKC-Trap to 2d.If we compare the 1d sketches of MUSCL-Trap, see figure 2(b), and MPRKC-Trap, see figure 5, we observe that the implicit zone has been extended by two cells.This holds also true in 2d: Given S n ij , (i) compute explicit fluxes using the MUSCL scheme and update all cells that have been treated fully explicitly by MUSCL-Trap to S (ii) use the predicted values S (1) ij for taking the second step of the explicit MPRKC scheme; due to slope reconstruction, this will 'cost' two layers of cells; therefore, the position of the transition cells has been shifted by two cell layers to the interior of the flow domain, compare figure 6; (iii) compute all implicit fluxes for the extended implicit zone; update cut cells, fully implicitly treated implicit cells in the extended implicit zone, and transition cells to S n+1 ij .
Note that the MPRKC-Trap scheme has been constructed to have a third-order one step error on Cartesian meshes.

Numerical results in 2d
We will consider two different tests in the following: We will first compare the various mixed schemes on a mesh that contains only Cartesian cells.This is to confirm our analytical considerations above.Afterwards we will compare the mixed schemes on a mesh that contains cut cells.)) so that the peak of the Gaussian aligns with where we switch the scheme.We choose the velocity field parallel to the angle β, i.e., γ = β, and set (u, v) T = (2, 2 tan(β)).We use T = 0.15.The initial data and the solution at the final time is shown in figure 8.We use ν = 0.8.The time step ∆t for the mixed schemes involving MUSCL and MUSCLmod is computed using ( 18), for the mixed scheme involving MPRKC we use (20).We compare MUSCL-Trap, MUSCLmod-Trap, and MPRKC-Trap.We also include results using the explicit versions of the scheme everywhere (i.e., in that case we do not switch to an implicit scheme).
In figure 9 we show the error after taking one time step, both in the L 1 norm and in the L ∞ norm.As expected, all fully explicit versions show a third-order one step error.For the mixed schemes, both new versions (MUSCLmod-Trap and MPRKC-Trap) show third order in L ∞ whereas MUSCL-Trap only converges with second order.
In figure 10 we show the error at time T .We observe second-order convergence for the L 1 norm for all 6 schemes.The error for MUSCL is smallest, followed by MUSCL-Trap.For the error in the L ∞ norm, all schemes except for MUSCL-Trap show second order.MUSCL-Trap converges on the finest meshes with order 1.5.So there is some error accumulation from the second-order one step error but we only seem to lose half an order.We have observed an order of 1.5 for MUSCL-Trap on Cartesian meshes in several tests.In terms of absolute error sizes, we observe the smallest errors for MUSCL, MUSCLmod, and MUSCLmod-Trap, with them being essentially identical on the finest mesh.

Test on cut cell mesh
Our discussion of the schemes in 2d so far has focused on the transition error between the explicit and implicit schemes in the absence of cut cells.The goal was to examine and eliminate error source (1) separately from error sources (2)(a)-(c).We now add cut cells again and conclude this contribution with a numerical comparison involving the following three schemes: MUSCL-Trap, MUSCLmod-Trap, MPRKC-Trap.The goal is to examine whether the reduction of the transition error leads to more accurate results and better convergence orders for the cut cell situation as well.
Test 4: We consider the ramp setup shown in figure 6.We again choose γ = β.We will consider 4 different angles for this test: 10 • , 20 • , 30 • , and 40 • .The setup is similar to Test 3 with the main difference being that we now have cut cells along the ramp.The Cartesian mesh is chosen to cover [0, 1] 2 with ∆x = ∆y.We use as test function again s 0 = 1 + exp(−120 * ((x − x 0 ) 2 + (y − y 0 ) 2 )).The starting point of the ramp and the centering of the test function (x 0 , y 0 ) varies for the different angles.They are chosen in such a way that at the initial time, the cut cells are located around the peak, see, e.g., figure 11 for the solution at the final time for ramp angle 30 • .
Besides the accuracy of the 3 different schemes, we also want to examine the effect of the inaccurate slope reconstruction on cut cells and transition cells.Typically, we use a least squares fit to construct gradients on cut cells and transition cells, see, e.g., [32], which is only first-order accurate.As an easy way of obtaining slopes with higher accuracy, we will also include results that use analytic slopes on cut cells and transition cells (together with standard  second-order slope reconstruction on Cartesian cells that are updated fully explicitly).These will be marked with 'ana' in the legend whereas 'LS' implies that the least squares slope reconstruction has been used.The tests will focus on comparing the accuracy of the different schemes.In terms of cost, for all schemes we need to solve implicit systems of size O(N ).The costs for MUSCL-Trap and MUSCLmod-Trap are pretty comparable.MPRKC-Trap is somewhat more expensive due to applying the explicit scheme twice, the potentially reduced time step length caused by the different CFL conditions, and the extended implicit zone, which is also reflected in somewhat longer running times.
In figure 11 we show the results for the error in the L 1 norm for the 30 • ramp.The results for the other angles are very comparable.We observe second order convergence for all schemes.Generally, using the original MUSCL scheme as explicit scheme leads to slightly smaller L 1 errors for this test.We have also included results for versions of MUSCL-Trap and MUSCLmod-Trap with an extended implicit zone, where the implicit zone has been extended by 2 cells.These are marked as 'ext2'.As a result the set of implicitly treated cells is now the same as for MPRKC-Trap.
In figure 12 we present the results for the error in the L ∞ norm for all 4 angles.Recall that the test function was chosen such that we expect the bigger errors on cut cells, which is actually the case.It is very common for errors on cut cells in the L ∞ norm to show a zig-zag behavior.In table 6 we therefore also present slopes that we get from fitting straight lines through the data shown in figure 12 by means of a least squares approach.While the details vary a bit for the different ramp angles, we can generally observe that • When using the LS slopes, the 3 different schemes (MUSCL-Trap, MUSCLmod-Trap, and MPRKC-Trap) show similar absolute errors and also similar convergence orders (with the exception of the 40 • ramp).Having improved the transition error does not result in smaller errors overall for this test.
• Using analytic slopes reduces the error (in terms of its size) significantly, with factors varying between 4 and 10.
• For MPRKC-Trap, using analytic slopes instead of LS slopes improves the convergence orders, within a range of 0.15 to 0.42.
• For MUSCL-Trap and MUSCLmod-Trap, we do not see significant improvement in convergence orders when using analytic slopes instead of LS slopes.For extending the

Conclusions and Future Plans
In this contribution, we analyzed the accuracy of the mixed explicit implicit scheme consisting of MUSCL as explicit scheme and Trapezoidal rule with slope reconstruction as implicit scheme.
For the one step error, we identified several second-order error sources linked to the irregular size of the cut cells as well as a second-order transition error when switching from explicit to implicit schemes.For linear advection in 1d, we can show that these errors do not accumulate in the usual way and that the resulting scheme is second-order accurate.This is not the case in 2d.We therefore introduced two new mixed schemes, MUSCLmod-Trap and MPRKC-Trap, with improved transition errors that are based on exchanging the explicit scheme.When using the mixed scheme on a fully Cartesian mesh, this led to improved convergence orders in 2d.When using these new schemes on a test involving cut cells however there was no significant difference to using the original MUSCL-Trap.Using analytic slopes led to a slight improvement of convergence orders and a significant improvement in the actual error size.Note though that the improved convergence orders are not that different from the newer results with DG codes on cut cell meshes for piecewise linear polynomials, which show reduced convergence orders in the L ∞ norm of 1.5 to 1.6 [15,20] as well.So maybe it is to much to aim for full second order in L ∞ for the ramp test.
More intensive numerical tests are needed to make a definite statement.It currently seems that it might only pay off to reduce the transition error if also the other errors sources (caused by the irregularity of the cut cells) are taken care of.The first step would be to upgrade the slope reconstruction on cut cells and neighbors of cut cells to second-order.This can be achieved by fitting quadratic polynomials.Our initial attempts of implementing this have led to irritating results.It is currently not clear whether there is a bug in the implementation or whether there is some kind of weird interaction going on with the solves necessary for the implicit time stepping.We currently use an approximate Newton scheme for that with the approximate Jacobian being based on a first-order discretization.

Test 1 :
Consider the model problem shown in figure 1.One cut cell is located in the interval [0.5, 0.5 + αh].The overall domain length is given by [0, 1 + αh].Periodic boundary conditions are used.The test function is sin 2π(x+0.36)1+αh

Figure 4 :
Figure 4: New 1d model problem: The grid consists of L blocks with K + 1 cells that each contain K cells of length H (hereH = h, h 2 , h4) and one cut cell of length αH in the middle.When the mesh is refined by a factor of 2, we keep K (and α) fixed but double L. Analogous to figure 1, we refer to cut cells as (cells of type) 0, to the left Cartesian neighbors of cut cells as (cells of type) -1, and to the right Cartesian neighbor as (cells of type) 1.

Figure 6 :
Figure 6: Switching between schemes in 2d: Cut cells are marked by pink color.Implicit fluxes are marked by a bold, red line, explicit fluxes by blue.Left: MUSCL-Trap: 'Transition' cells (marked with 'T') are full Cartesian cells that share an edge with a cut cell.Right: MPRKC-Trap: The shaded cells corresponds to the layer of cells closest to the embedded object that can be updated using the MUSCL predictor step.Transition cells (marked with 'T') have been shifted by 2 cell layers compared to MUSCL-Trap.

Figure 7 :
Figure 7: Setup for testing transition error in 2d: All cells are Cartesian cells but we simulatethe switch in time stepping.To do so, we intersect a Cartesian grid with a line with angle β.All cells cut by this line (marked in pink) are treated as 'fake' (as Cartesian) cut cells.All 'fake' cut cell edges as well as all edges of cells below the 'fake' cut cells use the implicit scheme (marked in red).All remaining edges (marked in blue), i.e., all edges of cells above the 'fake' cut cells use the explicit scheme.The switch takes place on the cells marked with a 'T'.

Figure 11 :
Figure 11: Results for Test 4 for 30 • ramp: Left: Zoom on the solution at time T using MUSCL-Trap for N = 512.Right: L 1 error for various schemes.

Table 1 :
Result of Test 1: MUSCL-Trap for the model problem in Fig. 1.

Table 2 :
Result of Test 2: Error for MUSCL-Trap for the model problem shown in Fig. 4.

Table 3 :
Result of Test 2: Error for MUSCL-Trap using forward difference quotients for the model problem shown in figure4.The new grid function w given by (11) is used for initialization and error computation.

Table 4 :
Result of Test 2: Error for MUSCLmod-Trap for model problem 4 with α = 1 (i.e., equidistant mesh but including change in time stepping scheme).

Table 5 :
Result of Test 2: Error for MPRKC-Trap for model problem 4 with α = 1 (i.e., equidistant mesh but including change in time stepping scheme).