A Sub-Element Adaptive Shock Capturing Approach for Discontinuous Galerkin Methods

In this paper, a new strategy for a sub-element based shock capturing for discontinuous Galerkin (DG) approximations is presented. The idea is to interpret a DG element as a collection of data and construct a hierarchy of low to high order discretizations on this set of data, including a first order finite volume scheme up to the full order DG scheme. The different DG discretizations are then blended according to sub-element troubled cell indicators, resulting in a final discretization that adaptively blends from low to high order within a single DG element. The goal is to retain as much high order accuracy as possible, even in simulations with very strong shocks, as e.g. presented in the Sedov test. The framework retains the locality of the standard DG scheme and is hence well suited for a combination with adaptive mesh refinement and parallel computing. The numerical tests demonstrate the sub-element adaptive behavior of the new shock capturing approach and its high accuracy.


Introduction
For non-linear hyperbolic problems, such as the compressible Euler equations, there are two major sources of instabilities when applying discontinuous Galerkin (DG) methods as a high order spatial discretization.(i) Aliasing is caused by under-resolution of e.g.turbulent vortical structures and can lead to instabilities that may even crash the code.As a cure, de-aliasing mechanisms are introduced in the DG methodology based on e.g.filtering [31,30], polynomial de-aliasing [37,58,38] or analytical integration [32,26], split forms of the non-linear terms that for instance preserve kinetic energy [62,23,25,18], and entropy stability [5,61,48,7,47,45,6,46,44,3,27,22,20].(ii) The second major source of instabilities is the Gibbs phenomenon, i.e. oscillations at discontinuities.It is well known that solutions of non-linear hyperbolic problems may develop discontinuities in finite time, so called shocks.While aliasing issues can be mostly attributed to variational crimes, i.e. a bad design and implementation of the numerical discretization with simplifications that affects aliasing (e.g.collocation), oscillations at discontinuities have their root deeper in the innards of the high order DG approach and are unfortunately an inherent property of the high order polynomial approximation space.Oscillations are fatal when physical constraints and bounds on the solution exist that then break because of over-and undershoots, resulting in nonphysical solutions such as negative density or pressure.

arXiv:2011.03338v1 [math.NA] 6 Nov 2020
A Sub-Element Adaptive Shock Capturing Approach for Discontinuous Galerkin Methods PREPRINT In this work, we focus on a remedy for the second major stability issue, often referred to as shock capturing.There are many shock capturing approaches available in DG literature, for instance based on slope limiters [11,10,41] or (H)WENO limiters [28,51,66], filtering [2], and artificial viscosity [49,9].Slope limiters or (H)WENO limiters are always applied in combination with troubled cell indicators that detect shocks and only flag the DG elements that are affected by oscillations.Only in these elements, the DG polynomial is replaced by an oscillation free reconstruction using data from neighboring elements.Element based slope limiting, e.g., with TVB limiters is effective for low order DG discretizations only, as its accuracy strongly degrades when increasing the polynomial order N , as in some sense the size of the DG elements ∆x gets larger and larger with higher polynomial degree N .The accuracy of the slope limiters is not based on the DG resolution ∼ ∆x N , but only on the element size ∆x.Considering that a finite volume (FV) discretization with high order reconstructions and limiters typically resolve a shock within about 2-3 cells, the shock width for high order DG schemes with large elements would be very wide when relying on element based limited reconstructions.The same behavior is still an issue for high order reconstruction based limiters such as the (H)WENO methodology.While these formally use high order reconstructions, its leading discretization parameter is still ∆x and not ∆x N and hence the shock width still scales with ∆x.Sub-element resolution, i.e. a numerical shock width that scales proportional to ∆x N , can be achieved with e.g.artificial viscosity.The idea is to widen the discontinuity into a sharp, but smooth profile such that high order polynomials can resolve it.Again, some form of troubled cell indicator is introduced to not only flag the element that contains the shock (or that oscillates) but also determine the amount of necessary viscosity, e.g., based on local entropy production [67].It is interesting that less viscosity is needed, i.e. sharper shock profiles are possible, the higher the polynomial degree of the DG discretization.Shocks can be captured in a single DG element if the polynomial order is high enough [49].An issue of artificial dissipation is that a high amount of artificial viscosity is needed for very strong shocks, which makes the overall discretization very stiff with a very small explicit time step restriction.Local time stepping can be used to reduce this issue [24] or specialized many stage Runge-Kutta schemes with optimized coefficients for strongly dissipative operators [33].
An alternative idea to achieve sub-element resolution is to actually replace the elements by sub-elements (lower order DG on finer grid) or maybe even finite volume (FV) subcells.A straight forward, maybe even a natural idea in the DG methodology is to use hp-adaptivity, where the polynomial degree is reduced and simultaneously the grid resolution is increased close to discontinuities.When switching to the lowest order DG, i.e a first order DG, which is nothing but a first order FV method, it is assumed that the inherent artificial dissipation of the first order FV scheme is enough to clear all oscillations at even the strongest shocks.The low accuracy of the first order FV method is compensated by a respective (very strong) increase in local grid resolution through h-adaptivity.An obvious downside to the hpadaptivity approach is the strong dynamic change of the approximation space during the simulation and the associated strong change of the underlying data structures.Thus, as a simpler alternative, subcell based FV discretizations with a first order, second order TVD or high order WENO reconstruction was introduced, e.g., in [60,57,15,63,14].The idea is to completely switch the troubled DG element to a robust FV discretization on a fine subcell grid.This approach keeps the change of the data local to the affected DG element and helps to streamline the implementation.The robust and essential oscillation free FV scheme on the finer subcell grid in the DG element is constructed based on the available DG data.This helps to keep the data structures feasible and reduces the shock capturing again to an element local technique instead of changing the global mesh topology and approximation space.Consequently, this again allows to combine the approach with the idea of a troubled cell mechanism to identify the oscillatory DG elements and replace the high order DG update with the robust FV update.As noted, there are many different variants for the FV subcell scheme with the high order WENO [60,57,15,63,14] variant being especially interesting, as it alleviates the 'stress' on the troubled cell indicator.When accidentally switching the whole high order DG element to the subcell FV scheme in smooth parts of the solution, accuracy may strongly degrade.A first order FV scheme, even on the fine subcell grid, would wash away all fine structure details that the high order DG approach simulated.Thus, from this point of view, a high order subcell WENO FV scheme is highly desired in this context.A downside of a subcell WENO FV scheme is the non-locality of the data dependence due to the high order reconstruction stencils.This is especially cumbersome close to the interface of the high order DG elements, where the reconstruction stencils reach into the neighboring DG elements and thus fundamentally change the data dependency footprint of the resulting implementation, that consequently has a direct impact on the parallelization of the code.
In this paper, the goal is to augment a high order DG method with sub-element shock capturing capabilities to allow the robust simulation of highly supersonic turbulent flows which feature strong shocks, as e.g. in astrophysics, without changing the data dependency footprint.Instead of flagging an element and completely switching from the high order DG scheme to the subcell FV scheme, we aim to smoothly combine different schemes in the flagged element.To achieve this, we firstly interpret the DG element as a collection of available mean value data.As depicted in Fig. 1, local reconstructions allow to define a hierarchy of approximation spaces, from pure piecewise constant approximations PREPRINT (subcell FV) up to a smooth global high order polynomial (DG element), with all piecewise polynomial combinations (sub-element DG) in between.We then associate each possible approximation space with the corresponding FV or DG discretization to define a hierarchy of schemes that is available for the given collection of data.

DG O8
Figure 1: Sketch of the main idea in a one-dimensional setting.Starting with eight mean values (red), we interpret them as a collection of available data.We construct four different approximation spaces (blue): piecewise constant on a subcell grid, piecewise linear and piecewise cubic on a sub-element grid, and the full 7th degree polynomial on the element.
Instead of switching between these different schemes, we aim to smoothly blend them in order to achieve the highest possible accuracy in every DG element.Assuming that the first order subcell FV approximation has enough inherent dissipation to capture shocks in an oscillation free manner, the goal is to give the FV discretization a high enough weight at a shock and then smoothly transition to higher order (sub-element) DG away from the shock.In contrast to the common approaches, where one indicator is used for the whole DG element, we aim to introduce sub-element indicators that are adaptively adjusted inside the DG element.A difficulty that arises when using sub-element indicators and an adaptive blending approach is to retain conservation of the final discretization for arbitrary combinations of blending weights.
The final discretization is a blended scheme, where the weights of the different discretizations may vary inside the DG element.To demonstrate this idea, the resulting behavior of the new approach is depicted in Fig. 2 for the example of a Sedov blast wave in 2D.In this particular setup, we start with 8th order DG as our baseline high order discretization.The grey grid lines indicate the 8th order DG element borders.Re-interpreting the 2D 8th order DG element as 8 2 mean value data, we can construct either a first order FV scheme on 8 2 subcells, a fourth order DG scheme on 2 2 , or a second order DG scheme on 4 2 sub-elements.The blending of the four available discretizations is visualized via a weighted blending factor (introduced in equation ( 57)) in Fig. 2 ranging from 1 (dark red, basically pure first order FV) up to 8 (dark blue, basically pure 8th order DG) and all combinations in between.The gray grid lines represent the 8th order DG elements.The color represents the blending process, where dark red (value 1) corresponds to pure first order FV and dark blue (value 8) corresponds to the pure 8th order DG scheme, with all combination in between.Note, that the plotted value is only for illustration and does not correspond to an actual mixed mathematical order of the final discretization.But the value nicely illustrates how even inside a DG element, the discretization is adaptively adjusted.
Clearly, the new approach adaptively changes on a sub-element level with a very localized low order near the shock front, while a substantial part of the high order DG approximation is recovered away from the shock.We refer to the validation and application sections, Sec. 3 and Sec. 4, for a more detailed discussion.
The remainder of the paper is organized as follows: in Sec. 2 we introduce the numerical scheme with the blending of the different discretizations.Sec. 3 include numerical validations and Sec. 4 shows an application of the presented approach with the conclusion given in the last section.

The Sub-Element Adaptive Blending Approach
The derivations in sections 2.1-2.6 are done in 1D, where the domain Ω is divided into Q non-overlapping elements.Each element q with midpoint x q and size ∆x q is mapped to the reference space as Each element holds a total number of degrees of freedom (DOF) N , where we require N ∈ {2 r , r ∈ N}.For the FV method, the element is divided into N subcells of size ∆xq N leading to a uniform grid with midpoints µ i and interfaces µ i± 1  2 in the reference space at For the sub-element DG scheme of order n = 2 r < N, r ∈ N, the element is divided into N n uniform sub-elements of size n∆xq N , and the mapping of each sub-element s to its (sub-)reference space reads as Fig. 18 in appendix 6.1 provides a visual representation of the reference element and its hierarchical decomposition into sub-elements and subcells. PREPRINT

Finite Volume Method
Consider the general 1D conservation law (4) We approximate the exact solution u(x) in element q with N mean values residing on the uniform subcell grid (2).Then the 1D FV semi-discretization on subcell i reads where f * denotes a consistent numerical two point flux.The adjacent values from neighboring elements q − 1 and q + 1 are given by (ū 0 ) q := (ū N ) q−1 and (ū N +1 ) q := (ū 1 ) q+1 .We define the Finite Volume Difference Matrix and rewrite (6) in compact matrix form with the numerical flux vector f * q ∈ R N +1 .For later use we split ∆ into a volume and a surface operator: Remark: The FV volume operator vanishes when summed along columns, i.e.

Discontinuous Galerkin Method
In this section we briefly outline the Discontinuous Galerkin Spectral Element Method (DGSEM) [1,36].To derive the n-th order DG scheme within sub-element s, (eq.( 3)) residing inside element q, we first rewrite the conservation law (4) into variational form with a test function φ and apply spatial integration-by-parts to the flux divergence We choose the test functions φ as Lagrange functions spanned with the collocation nodes ξ i ∈ − 1 2 , 1 2 , and use the polynomial ansatz where ũ(t) = ũ1 (t), . . ., ũn (t) T and (ξ) = 1 (ξ), . . ., n (ξ) PREPRINT are the vectors of nodal coefficients and Lagrange polynomials.We use collocation of the non-linear flux functions with the same polynomial approximation.Numerical quadrature is collocated as well.The nodes are the Gauss-Legendre quadrature nodes.These quadrature nodes and the interpolation polynomials are also illustrated in Fig. 18 in the appendix 6.1.
Inserting everything into (11) gives the semi-discrete weak-form DG scheme The diagonal mass matrix is with associated numerical quadrature weights ω i , and the differentiation matrix is The DG volume flux fs is computed directly with the nodal values ũs of the polynomial ansatz ( 13) and the surface flux vector is constructed with the numerical flux evaluated at DG sub-element boundaries where b ± are the boundary interpolation operators The boundary evaluation operator can be compactly written as We rewrite (15) and arrive at with Remark: The DG volume operator D vanishes when contracted with the vector of quadrature weights Proof: Here, we use the fact that D 1 = 0 is equivalent to taking the derivative of the constant function u(x) = 1, i.e.

Projection Operator
We consider a given polynomial of degree n − 1 with its nodal data ũ ∈ R n in a (sub-)element, which we want to project to n µ mean values on regular subcells.We define the projection operator P (n→nµ) ∈ R nµ×n component-wise via the mean value of the polynomial in the interval of subcell µ i , The integration of the Lagrange polynomial is done exactly with an appropriate quadrature rule.

Reconstruction Operator
The inverse of the quadratic non-singular projection matrix P (n→n) can be interpreted as a reconstruction operator from n mean values to n nodal values of a polynomial with degree n − 1: However, the naïve reconstruction (28) can lead to invalid data such as negative density or energy.Considering a convex set of permissible states for a given conservation law [65] we want to have a limiting procedure which recovers permissible states while still conserving the element mean value.Provided that the element average is a valid state, we expand the reconstruction process (28) by the following limiting procedure with the dyadic product 1 ⊗ 1 T ∈ R n×n and the "squeezing" parameter β ∈ [0, 1] chosen such that where which shows that the extended reconstruction operator preserves the original element average, i.e.

Remark:
The expanded reconstruction operator fulfills the following relationships: PREPRINT Proof: We begin with the second relation.From (26) we know: The first relation follows as The specific calculation of the parameter β depends on the permissibility constraints of the conservation law at hand.In section 3.1 we solve the compressible Euler equations where the density ρ and pressure p are positive values.Based on [65] we calculate the squeezing parameter β as where ρ , p are the element averages ( 29), ρ min , p min are the minimum values of the unlimited polynomial given in (31) and = min( ρ , p ).

Remark:
We observed that the indispensable positivity limiting can have a potential negative impact on the robustness of the high order DG scheme, especially for high polynomial degrees.It seems that huge artificial jumps and gradients can be generated at element interfaces.To alleviate this problem it is advisable to only allow small squeezing margins between 5% to 10%.If for example the squeezing parameter β falls below 0.95 the reconstructed high order polynomial is decided un-salvageable and the next lower order scheme is considered.

Single-level Blending
We first focus on the case where only two different discretizations per element q are blended: a single-level blending of the low order FV scheme and the high order DG scheme Note that we do not blend the solutions, but directly the discretizations themselves, i.e. the right hand sides which we denote by ∂ t u q .If both schemes operate with different data representations, appropriate transformations ensure compatibility during the blending process.In our case we use the projection operator P (N →N ) introduced in section 2.3 to transfer the nodal output of the DG operator to subcell mean values.The input for the DG scheme on the other hand, i.e. the nodal coefficients ũq ∈ R N , are reconstructed from the given mean values ūq ∈ R N as described in section 2.4.The naïve single-level blended discretization reads as Inserting ( 8) and ( 23) gives The fundamental issue of the naïve single-level blending discretization (35) is that it is not conservative if the blending parameter α q varies from element to element.It turns out that the direct blending of the surface fluxes of both discretizations, namely f , leads to a non-conservative balance of the net fluxes across element interfaces.To remedy this issue the goal is to define unique surface fluxes common to both, the low order and the high order discretization, such that they can be directly blended.We thus compute a blending of the reconstructed interface values to evaluate a unique common surface flux f * q .First, we define two interface vectors of length N + 1 were calculated with (37) where we assumed a blending factor of α = 0.75.
where the outermost interface values are given by A concrete example for the blended boundary values u ± is shown in Fig. 3.
The common surface flux f * q is then evaluated with the interface vectors u + q and u − q as To make the DG surface operator B compatible with f * q , we adapt the notation by inserting an additional column of zeros Corollary: Contracting the DG surface operator B (0) with the quadrature weights ω is equivalent to the contraction of the FV surface operator, i.e.
Proof: We expand the left side of ( 40) and write We replace f * (FV) q and f * (DG) q in (35) with f * q , and additionally, write the DG operators compactly, including the projection operator, as

PREPRINT
The final form of the single-level blended discretization is Remark: It is easy to see that with α q = 0 for all elements q, the pure subcell FV discretization is recovered and with α q = 1 for all elements the blending scheme recovers the pure high order DG method.
Lemma: Given arbitrary blending factors α q ∈ R for each element q the blending scheme ( 42) is conservative.
Proof: We discretely integrate the single-level blended discretization over all elements q with a total number of With ( 9) and ( 41) we get With properties ( 10), ( 24) and ( 26) the volume terms vanish, i.e.
We reformulate the DG surface term with ( 26) and ( 40) as The expansion of the telescopic sum only leaves the outermost surface fluxes , which represent the change due to physical boundary conditions.For instance in case of periodic boundary conditions, these two fluxes would cancel to zero.This concludes the description of the single-level blending scheme.

Multi-level Blending
As mentioned in the introduction, the general idea is the construction of a hierarchy of discretizations, with lower order DG schemes on sub-elements.In principle, it is thus possible to define a multi-level blending discretization, where the low order FV scheme is blended with DG variants of different approximation order.This multi-level extension is driven by the desire to retain as much 'high order DG accuracy' as possible, especially on a sub-element level.
For the discussion, we consider a specific example setup with a DG element with polynomial degree N p = 7.The data of this DG element is collected in form of N = 8 mean values on a regular subcell grid.Our goal is to blend a first order subcell FV scheme, a second and fourth order sub-element DG scheme and the full eight order DG scheme.The construction of the individual schemes follows the discussion in the previous sections 2.1, 2.2, and 2.5.In this section, the extension to a multi-level blending approach is presented.
We get the first data interpretation and its corresponding discretization by directly using the mean values of element q as a first order approximation space with no reconstruction at all ūO(1) q := ūq = (ū 1 , ū2 , . . ., ū8 ) T q .Next, we interpret the eight mean values as four second order DG sub-elements.The reconstruction matrix R O(2) ∈ R 2×2 transforms two adjacent mean values to two nodal values

PREPRINT
The fourth and eight order approximations are constructed analogously as and ũO( 8) 8 with R O(4) ∈ R 4×4 and R O(8) ∈ R 8×8 .Here we used the Kronecker product ⊗ in conjunction with the identity matrix 1 n = diag(1, . . ., 1) ∈ R n×n in order to generate appropriate block-diagonal matrices.Fig. 18 in appendix 6.1 illustrates the hierarchy of the approximation spaces for this example.
We correspond the approximation spaces with the respective discretizations where the DG volume fluxes are computed with the respective polynomial (sub-)element reconstructions of orders n = {2, 4, 8}: Similar to (39) the DG surface operators 1 N/n ⊗ BO(n) had to be slightly adapted in order to be compatible with the common surface flux f * .We list the modified boundary evaluation matrices in the appendix 6.2.Finally, all four candidate discretizations are blended starting from the lowest order up to the highest order discretization where denotes the component-wise multiplication with the vector of blending parameters.Each DG sub-element computes its own blending factor according to the local data, i.e. we get 4 blending factors for the O(2) variant, 2 blending factors for O(4) and 1 for O(8) T .The actual strategy for the blending factors is described in the next section 2.7.
Again, as discussed in the case of a single-blending approach in section 2.5, special care is needed to preserve the conservation of the final discretization.Expanding on the idea in the single-level case, we aim to compute unique surface fluxes for each discretization with a uniquely defined interface value at sub-element or subcell interfaces.We PREPRINT evaluate all high order DG polynomials at the FV subcell interfaces u where the interpolation operator I ±O(n) ∈ R n×n prolongates to the embedded FV subcell interfaces µ i±1/2 of the respective DG sub-element, i.e.
In Fig. 18  We start with the first order interpolation and stack on top the next levels of orders via convex blending To arrive at a complete vector of interface values we append the element interface values from the left and right neighbors q − 1 and q + 1.
These new common interface values allow us to evaluate a common surface flux (38) at each interface, which is again the key to preserve conservation in the final discretization.This concludes the description of the multi-level blending scheme.

Calculation of the Blending Factor α
A good shock indicator for high order methods is supposed to recognize discontinuities, such as weak and strong shocks, in the solution early on and mark the affected elements for proper shock capturing.On the other hand, the indicator should avoid to flag non-shock related fluid features such as shear layers and turbulent flows.There are numerous indicators available for DG schemes [50,15,14,60,2,9,66].
In this work we want to construct an a-priori indicator which relies on the readily available information within each element provided by the multi-level blending framework introduced in section 2.6.The principal idea is to compare a measure of smoothness of the different order reconstructions with each other.Smooth, well-resolved flows are expected to yield rather similar solution profiles compared to data that contain strong variations.
The smoothness measure σ within element q of the n-th order reconstruction ũO(n) q is inferred from the L 1 -norm of its first derivative.For example, if we want to calculate the blending factor α O(8) q for the eighth order we compute and (ξ) dξ PREPRINT utilizing information from the top-level eighth order element and its two fourth order sub-elements.The integrals are evaluated with appropriate quadrature rules.Then the blending factor is calculated by comparing the smoothness of the different orders as where cutoff(x lower , x, x upper ) = min(max(x lower , x), x upper ) and with two parameters τ A , τ S > 0. The design of equation ( 50) ensures that α is always within the unit interval and that for very small σ O(n) q the indicator does not get hypersensitive due to floating point truncation.If not noted otherwise we set the amplification parameter to τ A = 100 and the sensitivity parameter τ S = 0.05.The latter parameter has the biggest influence on the behavior of the indicator and its value demarcates the lower bound where it does not interfere with the convergence tests conducted in section 3.In order to capture all troubling flow features we independently apply the indicator (50) on all primitive state variables and select the smallest of the resulting blending factors.The calculations of the two fourth order blending factors α are done within each sub-element independently together with their respective second order smoothness measures σ , s = 1, . . ., 4.

Remark:
The indicator ( 50) is also used for the single-level blending scheme.

Sketch of the Algorithm
In the last part of this section, a sketch of the algorithm for the blending framework is presented.We present a general outline of the necessary steps for the 1D single-blending scheme evolved with a multi-stage Runge-Kutta time stepping method.We enter at the beginning of a Runge-Kutta cycle and do the following.
(I) Reconstruct the polynomial ũq from the given mean values ūq for each element q as in ( 28).
(II) If the reconstructed polynomial ũq contains non-permissible states, see (31), then calculate the limited version ũ(β) q as in (30) (III) If the squeezing parameter β q is below β L then set α q := 0 else compute the blending factor α q via (50) from the unlimited polynomial ũq .
(IV) Compute the boundary values u ± q via (37) and exchange alongside zone boundaries in case of distributed computing.
(V) Determine the common surface flux f * q via (38).(VI) Calculate the right-hand-side ∂ t ūq : If the blending factor α q above α H then compute ∂ t ūq with the DG-only scheme.
Else if the blending factor α q below α L then compute ∂ t ūq with the FV-only scheme.
Else compute ∂ t ūq with the single-level blending scheme (42).
(VII) Forward in time to the next Runge-Kutta stage and return to step (I).
The switching thresholds are set to α H := 0.99 and α L := 0.01 and the limiter threshold to β L := 0.95.Note that the algorithm only applies the blending procedure where necessary in order to maintain the overall performance of the scheme.
This concludes the presentation of the 1D blending scheme.The description of the blending scheme on 3D Cartesian meshes as well as the algorithm outline for the multi-level blending scheme can be found in appendix 6.3 and 6.4.

PREPRINT 3 Validation
For the computational investigations, the multi-level algorithm with the explicit SSP-RK (5,4) time integrator [59] is implemented in a Fortran-2008 prototype code with a hybrid parallelization strategy based on MPI and OpenMP.Management of the AMR and load balancing is provided by P4EST, a highly efficient Octree library [4].The maximal time step for dimensions d = 1, 2, 3 is estimated by the CFL condition where CF L := 0.8, N := 8 is the number of mean values in each direction of the element q and λq,max is an estimate of the maximum eigenvalue given in (55).For the numerical interface fluxes f * , we use the Harten-Lax-Leer (HLL) approximate Riemann solver [29] with Einfeldt signal speed estimates [16].

Governing Equations
We consider the compressible Euler equations with the vector of conserved quantities u = (ρ, ρ v, E) T , where ρ denotes the density, v the velocity, and E the total energy.We assume a perfect gas equation of state and compute the pressure as If not stated otherwise we choose γ = 1.4.The set of permissible states is given by For the CFL condition (51) the maximum eigenvalue is evaluated on all mean values ūq of element q.Given dimension d = {1, 2, 3} it reads as

Convergence Test
We use the manufactured solution method [52] and validate the 3D multi-level blending scheme on a periodic cube of unit length (L = 1) where the resolution of the mesh is incrementally doubled.We define our manufactured solution in primitive state variables as The final time of the simulation is T = 2 and the center of the domain is refined to introduce non-conforming interfaces in the computational domain.We determine the L ∞ -and L 2 -norms of the errors in the density and total energy.Tables 1 to 4 list the results for the first, second, fourth and eighth order multi-level blending schemes.The initial conditions and the source term of our manufactured solution (56) is in all cases evaluated and applied on the mean values via an appropriate quadrature rule to maintain high order.The results confirm, that the discretizations behave as designed in this assessment.

Conservation Test
The goal in this assessment is to demonstrate that the multi-level-blending discretization is conservative for all choices of blending factors.To do so, we adapt the same setup as in section 3.2 but deactivate the source term.The center of the domain is refined to introduce non-conforming interfaces in the computational domain.Additionally, as a stress test, the blending factors are randomly chosen and changed after each Runge-Kutta stage.As there are multiple blending factors at a given spatial location, we consider the following weighted blending factor to illustrate the distribution of the blending factors in Fig. 4. Note the limiting cases ᾱ = 1 for pure FV and ᾱ = 8 for the eighth order DG scheme.

PREPRINT
The simulation runs to T = 300 performing more than a quarter million timesteps.The result of the test is shown in Fig. 5 where we plot in log-scale the absolute value of the change of bulk ∂ t ūtotal (t), integrated over the whole domain, Q is the total number of elements and |Ω q | = ∆x q ∆y q ∆z q is the volume per element.The results show that the conservation error lies within the range of 64 bit (double precision) floating point truncation and hence confirm that the multi-level blending discretization is fully conservative up to machine precision errors.58)) of each conservative state variable.PREPRINT

1D Shock Tube Problems
We validate the multi-level blending scheme with three well-established 1D problems, namely Sod, Lax and Shu-Osher shock tubes.But first, in order to gain insights into the individual contributions of the different schemes, we introduce a more intuitive reformulation of the blending factors α O(n) .The multi-level blending operation (47) within element q can be interpreted as a linear superposition of four numerical schemes.For that, we collapse the following relations q , into one single term and define blending weights θ O(n) , i.e.
, where θ O(1) θ O( 8) By construction the blending weights have the following property and thus give a proper fraction of each contribution.
The first shock tube problem is the Sod shock tube [56].It is defined on the unit interval Ω = [0, 1] with a diaphragm located at x D = 0.5.The initial condition in primitive state variables reads The resolution is set to 32 elements of N = 8 mean values each.This amounts to 256 total DOF.The result for the density profile (top row) at the final simulation time T = 0.2 is presented in Fig. 6 together with the exact solution.
It shows the correct approximation of the rarefaction wave, contact discontinuity and the forward facing shock front.As designed, only at the shock front the blending scheme gets triggered, which is visible in the bottom row of the plot.The stacked bar chart directly corresponds to the blending weights θ O(n) in ( 59).The vertical dimension of the stacked bars completely fill the unit interval mirroring property (60).2) O( 4) O( 8) Sod Shock Tube Problem: density and blending weights at t = 0.2 The second test case is the Lax shock tube [42] with initial data set to ρ 0 (x), v 0 (x), p 0 (x) = 0.445, 0.689, 3.528 , x < x D , 0.5, 0, 0.571 , x ≥ x D .
All other simulation parameters as well as the resolution are the same as in the Sod test case.The result for the density profile (top row) at final simulation time T = 0.15 is presented in Fig. 7 together with reference solution on a finer grid of 1024 DOF with the third-order piece-wise parabolic method (PPM) [12] implemented in the astrophysics code FLASH (version 4.6, March 2019), see e.g.[21].For comparison, we also included the solution of the PPM with the same grid resolution of 256 DOF.All flow features are resolved correctly and the blending scheme is only triggered in the region around the forward facing shock.Furthermore, this example clearly reveals the adaptive blending on the sub-element level.2) O( 4) O( 8) Lax Shock Tube Problem: density and blending weights at t = 0.15 Here, we compare the results of the multi-level blending DGFV8 to PPM.In order to investigate the importance of the multi-level approach on the accuracy of the DGFV8 result, we additionally perform a simulation where we intentionally deactivate the sub-elements, i.e. α O(2) := 0 and α O(4) := 0. This is identical to the single-level blending approach presented in section 2.5, where only the first order FV scheme is blended with the eighth order DG scheme.The resolution is as before 256 total DOF whereas the reference solution is computed with PPM on a much finer grid of 2048 DOF.The numerical experiments shown in Fig. 8 nicely demonstrate the benefit of using the multilevel approach.Whereas the shocks are about equally resolved, the multi-level blending variant gives the best results in the smooth parts of the solution.Since the resolution is not very high, non-shock related small scale flow features cannot always be resolved by the eighth order DG scheme.This is especially visible in the range x = [1.688,2.25] where the lower order scheme has to completely or at least partially take over.2) O( 4) O (8) Shu-Osher Shock Tube Problem: density and blending weights at t = 1.8  59)) for the single-level blending while the bottom row shows the multi-level blending.The vertical grid lines depict the element boundaries.

2D Riemann Problems
In this section we present a selection of three 2D Riemann problems [43,53,40].The domain for all simulations is set to Ω = [−0.5,0.5] 2 with a uniform grid resolution of 64 2 elements, respectively 512 2 DOF.The setup consists of the four quadrants each initialized with their own constant states.The exact parameters of each Riemann problem are given in [43] where they are addressed with a fixed configuration number.For brevity we omit the setup parameters and refer to [43].In this work we show configurations 3, 4 and 6.We are interested in the change of the numerical solution when we incrementally add one level of blending order from one run to the next.Hence, we present three solutions for each 2D Riemann problem denoted by their respective multi-level schemes: DGFV2, DGFV4 and DGFV8.Note that the schemes are implemented in such a way that they operate on the same element size of N 2 = 8 2 mean values.The results are shown in Fig 9 to 11.The left column shows the density contour and the right column the weighted blending factors as in (57).The general observation is that with increasing order there is more structure visible in the density plots and the blending patterns get more nuanced in tracing the flow structures.[43]) computed with the (from top to bottom) DGFV2, DGFV4 and DGFV8 scheme.The resolution of 512 2 DOF is the same for all runs.Left: density contours in log-scale at final time T = 0.3.Right: Weighted blending factors as defined in (57).Dark blue represents the full eighth order DG and dark red the full first order FV scheme.PREPRINT

2D Sedov Blast
The Sedov blast problem [64,65,13] describes the self-similar evolution of a radially symmetrical blast wave from an initial pressure point (delta distribution) at the center into the surrounding, homogeneous medium.The analytical solution is given in [54,39].In our setup, we approximate the initial pressure point with a smooth Gaussian distribution with the spatial dimension d = 2, the blast energy E = 1 and the width σ such that the initial Gaussian is reasonably resolved.The surrounding medium is initialized with ρ 0 = 1 and p 0 = 10 −14 .Dimensional analysis [54] reveals that the analytical solution of the density right at the shock front is determined by With the adiabatic coefficient γ = 1.4,we investigate how close the numerical results match ρ shock = 6.The Cartesian mesh has a FV equivalent uniform grid resolution of 512 2 DOF, i.e. when computing with the full eighth order DG approximation space, the total number of DG elements is 64 2 .The spatial domain is Ω ∈ [−0.25, 0.25] 2 with the initial blast width σ = 5 × 10 −3 .We compare the accuracy of the results obtained with the single-level and the multi-level (DGFV8) blending discretizations as well as the PPM.Fig. 13 shows the shell-averaged density and pressure profiles at final time T = 0.05.Fig. 12 presents the numerical solution over the whole domain as computed with the multilevel DGFV8.To further illustrate the behavior of the multi-level and single-level blending approach we also show in Fig. 14 the weighted blending factors along the x-axis.The shock front is much sharper for the multi-level blending compared to the single-level blending.It can be observed how the weighted blending factor is dominated by the first order FV scheme (ᾱ ≈ 1) directly at the shock, but transitions quickly to a blended discretization (1 < ᾱ < 8) up to the full eighth order DG (ᾱ ≈ 8) away from the shock, even within a single DG element.This behavior demonstrates the sub-element adaptivity of our novel approach.Again, similarly to the shock tube section 3.4 and the 2D Riemann problem section 3.5 the results for the 2D Sedov blast wave show the benefit of the multi-level approach, with the numerical profiles even slightly closer to ρ shock = 6 compared to the PPM with equal resolution of 512 2 DOF.Sedov Blast with multi-level DGFV8 at t = 0.05 Sedov Blast with single-/multi-level DGFV8 and PPM at t = 0.05   2) O( 4) O (8) Sedov Blast at t = 0.05: blending factor profiles along x-axis  59)) are encoded with stacked bars in their respective colors.The vertical grid lines highlight the eighth order DG element boundaries.For reference the scaled density (ρ/6) and scaled pressure (p/5) are included.PREPRINT 4 Simulation of a Young Supernova Remnant Supernova models have been analyzed and discussed for many decades and since they unite a broad range of features such as strong shocks, instabilities and turbulence, they resemble a good test bed for our novel shock capturing approach in combination with AMR.
The general sequence of events of the presented supernova simulation is like this: We start with a constant distribution of very low density resembling interstellar media (ISM) that typically fills the space between stars.When a star explodes by turning into a supernova it ejects its own mass at very high speeds into the ISM preceded by a strong shock front heating up the ISM.The ejected mass is rapidly decelerated by the swept-up ISM giving rise to a so called reverse shock that travels backwards to the center.The interface, or more precisely the contact discontinuity, between shocked ejecta and shocked ISM is unstable and leads to a layer of slowly growing Rayleigh-Taylor instabilities.This gradually expanding layer, called supernova remnant, is of special interest since this is where astronomical observations reveal a lot of ongoing physics and chemistry, especially driven by mixing and turbulence.
We adapt the setup descriptions in [8,19,17] where we have the initial (internal) blast energy E and the ejecta mass M given in cgs (centimeter-gram-seconds) units.It is beneficial to convert the given units to convenient simulation units reflecting characteristic dimensions of the physical model at hand.Table 5 lists the conversion between cgs and simulation units and table 6 lists the initial parameters used in this simulation.The ambient density ρ a is related to the  mono-atomic particle (hydrogen) density n H via ρ a = m u n H , where m u is 1  12 of the mass of a carbon-12 atom.The ambient pressure p a is calculated from the ideal gas law, i.e. p a = n H k B T with the Boltzmann constant k B .There is some ambiguity regarding the ambient gas temperature T .The literature mentions a warm and neutral interstellar medium which is attributed to temperatures between 6 × 10 3 K and 10 4 K.The heat capacity ratio is fixed to γ = 5/3.The simulation time spans a period from t 0 = 10 yr to T = 500 yr.The expansion of the forward shock (eq.( 64)) is then expected to approximately reach R FS = 5 pc, which determines the size of the computational domain L := 5 pc.Fig. 15 depicts a schematic of the simulation setup.Due to the rotational symmetry of the setup it is sufficient to simulate just one octant of the supernova.The following formulas have been derived in [8] and were adapted to the current setup, i.e. power-law indices of (s, n) = (0, 7).The self-similar solution at initial time t 0 = 10 yr within the power law region, respectively the blast center, is given by We initialize the density as and the total energy with (61) where p 0 = p a , d = 3 and σ = 3 4 r(t 0 ).The initial momentum is (ρ v) 0 = 0. Since we are only interested in the evolution of the instability layer of the supernova we apply the following rules for mesh This allows us to assign an adaptive, high resolution shell of maximal refined elements following the remnant as it expands into the computational domain.The inner and outer radii of the shell are estimated as which have been found adequate via numerical experimentation.Moreover, up to t = 200 yr we enforce R inner = 0, which ensures that the first phase of the explosion is well resolved in any case.The refinement levels range from 2 to 6, which translates to a FV equivalent resolution from 2 2 • 8 = 16 up to 2 6 • 8 = 512 cells in each spatial direction.
We perform three simulations with the 3D multi-level blending schemes DGFV2, DGFV4 and DGFV8, analog to the 2D Riemann problems discussed in section 3.5.Fig. 16 shows the density slice (left column) sketched in Fig. 15 at the final simulation time of T = 500 yr, while Fig. 17 shows the corresponding 3D density contour rendering of the DGFV8 solution.The shock front partially left the domain, which is not considered a problem since the region of interest, namely the instability layer, is still completely covered.The areas of no interest, i.e. the outer region as well as the center, are only coarsely resolved by the AMR scheme.Clearly, an increase in order leads to a much more detailed remnant structure emphasizing the advantage of higher order schemes in resolving small scale turbulence driven by Rayleigh-Taylor instability, see [8].The weighted blending factor is shown in the right column of Fig. 16.The band of highly refined elements is clearly visible following the remnant as intended.Two distinctive lines of blending activity trace the front and reverse shocks.The plots also show a strong qualitative difference of the amount of scales in the instability layer for DGFV2 and DGFV4.Clearly, the DGFV4 result features more scales and finer structures as the DGFV2 simulation with the same FV equivalent grid resolution.The difference in scales between DGFV4 and DGFV8 is less pronounced which is probably related to the extensive blending activity with O4 inside the instability layer of the DGFV8 simulation.In this turbulent part we suspect that the standard collocation eighth order DG scheme might face aliasing instability issues as described in the introduction.There are techniques to reduce aliasing issues available, such as filtering, consistent integration and split-forms.However, this aspect is not the main topic of the paper and it is interesting to see, that the multi-level blending automatically adjusts to cope with these issues as well.

Conclusion
In this work, we introduce an adaptive sub-element based shock capturing approach for DG methods.We interpret an element first as a collection of piecewise constant data.This set of piecewise constant data can be interpreted and reconstructed in a variety of ways.We focus on piecewise polynomial approximation spaces, starting from a pure piecewise constant FV type interpretation (no reconstruction) up to the fully maximum order polynomial reconstruction and every combination in-between, e.g.piecewise linear and piecewise cubic reconstructions.In a second step, we link the data interpretation to a corresponding (high order) DG discretization.Thus, we get a hierarchy of discretizations acting on the same set of data.
The idea is then to adaptively blend this hierarchy of different discretization, where the low order variants are chosen close to discontinuities and the high order variants in smooth or turbulent parts of the simulation.Instead of having an element based troubled cell indicator approach, we use the different data interpretations again to compute sub-element localized indicators, which allows for a sub-element adaptive blending of the discretizations.When the blending can change throughout the element, special care is necessary to preserve exact conservation of the resulting multi-level blended discretization.We achieve exact conservation, by introducing unique, blended reconstruction states at subcell and sub-element interfaces.
In our prototype implementation, we further demonstrate that a natural combination of this sub-element adaptive approach is with adaptive mesh refinement.We base the AMR implementation on the P4EST Octree library, which allows for a straight forward parallelization of the whole computational framework.
We show standard numerical test cases to validate the new approach and a simplified model of a supernova remnant to highlight its high accuracy for challenging test cases with strong shocks and turbulence like structures.

Extension to 3D Cartesian Grids
The extension to higher spatial dimensions with Cartesian conforming grids is relatively straight forward via a tensorproduct strategy.In this section we want to present the most important building blocks for the 3D single-level blending scheme.The steps of the 3D multi-level blending are analogous as in section 2.6 and not detailed any further.Consider the general 3D conservation law We now have an element q with N 3 mean values as available data ūq ∈ R N ×N ×N , and the element sizes ∆x q , ∆y q and ∆z q .With the tensor-product of Lagrange polynomials on Legendre-Gauss nodes ξ i , i = 1, . . ., N , PREPRINT and the reconstruction operator R, which we introduced in section 2.4, we get the nodal coefficients as Here we make use of the n-mode product notation [34].A definition is given in appendix 6.5.In general all vector and matrix operations presented so far can be directly expanded to 3D via the n-mode product.The blending scheme (42) in 3D reads as The DG volume fluxes are calculated from the node values as fq ijk = f (ũ q ) ijk , (g q ) ijk = g (ũ q ) ijk and hq ijk = h (ũ q ) ijk , i, j, k = 1, . . ., N.
, we need to compute common surface fluxes in order to preserve conservation.Therefore, we define two selection operators, s − = (1, 0, . . ., 0, 0) T ∈ R N and s + = (0, 0, . . ., 0, 1) T ∈ R N , mapping the outermost mean values of ūq to the respective interface.Then the 3D analogue of the prolongation procedure (37)  together with the blending factor α q are supposed to be exchanged between processes in case of distributed computing.
Next, we reconstruct again a nodal representation from ū±d we calculate two candidate fluxes at the interface q − 1 2 .In x-direction it reads as f * Next, we determine a common surface flux f * q− 1 2 by blending the two candidate interface fluxes where α q− 1 2 = αq−1+αq 2 . The final step is to calculate the inner interface fluxes analogous to (38) f * q ijk = f * (ū i−1jk , ūijk ), i = 2, . . ., N, j, k = 1, . . ., N and insert the common surface fluxes f * The y-and z-direction are treated analogously.The computation of the blending parameter α for 3D follows the same steps as in section 2.7 where the integrals (49) are rewritten in 3D form.The presented parameters τ a and τ s stay the same.This concludes the 3D blending scheme with conforming interfaces.
We extend the 3D single-blending scheme to non-conforming grids where we assume 4:1 transitions only.See Fig. 20.The steps for the 3D multi-level blending scheme are analogous and not detailed any further.First, we define four matrix operators which allow us to construct refinement and coarsening procedures within the blending framework.We require that N = 2 l , l ∈ N, and define the following expansion and compression operators and PREPRINT coarser level finer level mortar layer with prolongated boundary values Moreover, we need the projection operator P (N →M ) introduced in section 2.3 for mapping N node values to M mean values.We define: mean values to mean values and Z := P (N →2N ) , Q := P (N → N

)
node values to mean values .
We refine the parent element q as rq = (1 The refined block rq is then split into 8 child elements.The reverse operation, namely coarsening, is done by compressing each child element ūr , r = 1, . . ., 8, separately: Afterwards the family of 8 blocks cq is glued together.These blending operations are especially useful for refining or coarsening of highly oscillatory data or at pronounced discontinuities. For the treatment of non-conforming 4:1 interfaces we need a procedure which maps the boundary values from the coarser element to a so-called mortar layer [35] that has the same resolution as the adjacent four smaller elements.See Fig. 20.To compute the mortar layer of the coarser element q, we formulate where i, j = 1, . . ., N and ξ k , η l are the collocation nodes of the tensor-product (67).We translate above equation into matrix notation and get where (L ± ) ij = i 1 2 ξ j ± 1 4 ωj 2 ωi , i, j = 1, . . ., N .Analogous to (75) we determine a common surface flux for the coarse side f *

Figure 2 :
Figure2: Visualization of the sub-element adaptive blending approach for the 2D Sedov blast wave, see Sec. 3 for a detailed discussion.The gray grid lines represent the 8th order DG elements.The color represents the blending process, where dark red (value 1) corresponds to pure first order FV and dark blue (value 8) corresponds to the pure 8th order DG scheme, with all combination in between.Note, that the plotted value is only for illustration and does not correspond to an actual mixed mathematical order of the final discretization.But the value nicely illustrates how even inside a DG element, the discretization is adaptively adjusted.

Figure 3 :
Figure 3: Schematic of four (N = 4) mean values and the reconstructed polynomial on the reference element.The cell boundary values u − q− 1 2 in appendix 6.1 two representatives of µ O(n) i±1/2 are shown which are aligned at the dotted vertical lines, indicating the interfaces of the FV subcells.

Figure 4 :
Figure 4: Conservation test of the 3D multi-level blending scheme: Slice through the computational domain showing weighted blending factor ᾱ according to (57) on a Cartesian non-conforming mesh with refinement levels 3 to 5. The black lines depict the boundaries of the corresponding 8th order DG elements.

Figure 5 :
Figure 5: Conservation test of the 3D multi-level blending scheme: Evolution of the integrated change of bulk log 10 (|∂ t ūtotal |) (eq.(58)) of each conservative state variable.

Figure 6 :
Figure 6: Sod Shock Tube Problem: Numerical solution of the density profile (top row) with the DGFV8 multilevel blending scheme together with the exact solution at final simulation time T = 0.2.The bottom row shows the blending weights θ O(n) (eq.(59)) encoded by stacked bars in different colors.The light grey vertical grid lines depict the element boundaries.

Figure 7 :
Figure 7: Lax Shock Tube Problem: Numerical solution of the density profile (top row) with the DGFV8 (256 DOF) multi-level blending scheme together with the reference solution (PPM, 1024 DOF) and PPM (256 DOF) at final simulation time T = 0.2.The bottom row shows the blending weights θ O(n) (eq.(59)) encoded by stacked bars in different colors.The vertical grid lines depict the element boundaries.

Figure 8 :
Figure 8: Shu-Osher Shock Tube Problem: Numerical solution of the density profile (top row) with the DGFV8 multi-level blending scheme together with the reference solution (PPM on 2048 DOF), PPM and single-level blending scheme all using 256 DOF.The center row shows the blending weights θ O(n) (eq.(59)) for the single-level blending while the bottom row shows the multi-level blending.The vertical grid lines depict the element boundaries.

Figure 9 :
Figure 9: 2D Riemann problem (configuration 3, see[43]) computed with the (from top to bottom) DGFV2, DGFV4 and DGFV8 scheme.The resolution of 512 2 DOF is the same for all runs.Left: density contours at final time T = 0.3.Right: Weighted blending factors as defined in(57).Dark blue represents the full eighth order DG and dark red the full first order FV scheme.

Figure 10 :
Figure 10: 2D Riemann problem (configuration 4, see[43]) computed with the (from top to bottom) DGFV2, DGFV4 and DGFV8 scheme.The resolution of 512 2 DOF is the same for all runs.Left: density contours at final time T = 0.25.Right: Weighted blending factors as defined in(57).Dark blue represents the full eighth order DG and dark red the full first order FV scheme.

Figure 11 :
Figure 11: 2D Riemann problem (configuration 6, see[43]) computed with the (from top to bottom) DGFV2, DGFV4 and DGFV8 scheme.The resolution of 512 2 DOF is the same for all runs.Left: density contours in log-scale at final time T = 0.3.Right: Weighted blending factors as defined in(57).Dark blue represents the full eighth order DG and dark red the full first order FV scheme.

Figure 12 :
Figure 12: 2D Sedov Blast: Numerical solution computed with DGFV8 at final simulation time T = 0.05.Left: density contours.Right: Weighted blending factors as defined in(57).Dark blue represents the full eighth order DG and dark red the full first order FV scheme.

Figure 13 :
Figure 13: 2D Sedov Blast: shell-averaged density and pressure profiles at final simulation time T = 0.05 computed with the single-/multi-level blending scheme (DGFV8) and PPM on a FV equivalent uniform grid resolution of 512 2 DOF.The vertical grid lines depict the element boundaries.

Figure 14 :
Figure 14: 2D Sedov Blast: We compare the blending factor profiles along the x-axis of the single-and multi-level blending schemes at final time T = 0.05.The blending weights θ O(n) (eq.(59)) are encoded with stacked bars in their respective colors.The vertical grid lines highlight the eighth order DG element boundaries.For reference the scaled density (ρ/6) and scaled pressure (p/5) are included.

Figure 15 :
Figure 15: Computational domain (cubic box) covering one octant of the supernova model.The faces at the coordinate axes are set to reflecting walls while the opposite sides are set to outflow.

Figure 16 :Figure 17 :
Figure 16: left column: 2D density slice (see figure 15) showing the instability region respectively remnant of the supernova at T = 500 yr simulated with the (from top to bottom) DGFV2, DGFV4 and DGFV8 scheme.right column: 2D slice of the weighted blending factors of the respective blending schemes at T = 500 yr.The black lines correspond to the element boundaries of the Cartesian non-conforming mesh.

PREPRINT 6 APPENDIX 6 . 1
Four levels of data interpretation of 8 mean values u(χ)

Figure 20 :
Figure 20: Schematic of a non-conforming interface with 4:1 adjacent cells for N = 4.The two kinds of boundary values occurring in (80) are illustrated.Compare with Fig. 19.

ū±d q± 1 2 coarse= ( 1 − 2 coarse
α q ) E ūq × d s ± E T face centered mean values + α q Z nodal boundary values ũq × d b ± ZT face centered mean values ∈ R 2N ×2N .(80)Note the similarity to the prolongation procedure for conforming interfaces (71).The coarse side of the mortar layer ū±d q± 1 is split up to match the faces of the four finer elements r = 1, . . ., 4. The boundary values of the smaller elements are constructed by applying (71) individually.Only the mortar layers together with the associated blending factors α r , r = 1, . . ., 4, are supposed to be exchanged between processes in case of distributed computing.The separate computation of the interface fluxes f * fine,r ∈ R N ×N , r = 1, . . ., 4, follows the formulas (72)-(75).While the resulting four fluxes can be directly copied back to the smaller faces without further treatment, they have to be mapped to the coarser side via L 2 -projection[35] coarse − f * fine,1 φ dξ dη + coarse − f * fine,3 φ dξ dη + coarse − f * fine,4 φ dξ dη = 0.PREPRINTThe integrals are evaluated in the reference space of the coarse face.Applying exact quadrature rules gives

Table 1 :
Experimental order of convergence of the first-order FV variant within the 3D multi-

Table 3 :
Experimental order of convergence of the fourth order DG variant within the 3D multi-level blending framework.

Table 5 :
Conversion from cgs units to simulation units.

Table 6 :
Hydrodynamical parameters in cgs units and simulation units.