WENO schemes on unstructured meshes using a relaxed a posteriori MOOD limiting approach

In this paper a relaxed formulation of the a posteriori Multi-dimensional Optimal Order Detection (MOOD) limiting approach is introduced for weighted essentially non-oscillatory (WENO) finite volume schemes on unstructured meshes. The main goal is to minimise the computational footprint of the MOOD limiting approach by employing WENO schemes—by virtue of requiring a smaller number of cells to reduce their order of accuracy compared to an unlimited scheme. The key characteristic of the present relaxed MOOD formulation is that the Numerical Admissible Detector (NAD) is not uniquely defined for all orders of spatial accuracy, and it is relaxed when reaching a 2nd-order of accuracy. The augmented numerical schemes are applied to the 2D unsteady Euler equations for a multitude of test problems including the 2D vortex evolution, cylindrical explosion, double-Mach reflection, and an implosion. It is observed that in many events, the implemented MOOD paradigm manages to preserve symmetry of the forming structures in simulations, an improvement comparing to the non-MOOD limited counterparts which cannot be easily obtained due to the multi-dimensional reconstruction nature of the schemes. c ⃝ 2020 The Author(s). Published by Elsevier B.V. This is an open access article under the CC BY license (http://creativecommons.org/licenses/by/4.0/).


Introduction
Coupled with increasingly larger computational power, modern schemes are mature enough to attack intensive gasdynamic and aerodynamic problems, as well as simulations of relativistic fluids with magnetic fields [1]. Nevertheless, the solution of partial differential equations in the presence of strong shocks and stiff gravitational source terms constitutes a very challenging computational problem. The advances on many numerical techniques strive for robust, yet high-order accurate solutions, which generally entail unfortunate increases in terms of computational cost. In the last decade an overwhelming number of numerical methods have been developed for the compressible Navier-Stokes equations on unstructured meshes in the context of large-eddy simulations (LES). The common goal of all the numerical methods in this context is to achieve high-order of accuracy in a computationally efficient manner, while obtaining physically meaningful results free of spurious oscillations. Finding the ideal balance between these three facets (accuracy, efficiency, non-oscillatory behaviour) is far from obvious. Therefore the numerical methods employ some form of adaptivity in order to adapt their spatial accuracy at the presence of sharp gradients that have the potential to contaminate the solution. The majority of these methods are developed in the finite-volume (FV) , the discontinuous Galerkin (DG) [16,[23][24][25][26][27][28][29][30][31], the spectral finite-volume (SFV) [32][33][34][35][36][37][38][39], and the Flux Reconstruction (FR) [23,[40][41][42][43] mathematical frameworks that have made novel efforts in trying to address these challenges for unstructured meshes.
One of the interesting and promising techniques that have been used in various numerical frameworks is the MOOD introduced by Diot et al. [44]. Since then, it has been finely tuned by Diot, Clain, Loubére, Dumbser, Nogueira and others [44][45][46][47][48][49][50][51] and has been successfully applied to a wide range of challenging flow problems. The key benefit of the MOOD approach is that it increases the robustness of the numerical framework employed since the solutions at any cell that do not conform to certain, pre-selected criteria are recomputed with spatially lower-order numerical schemes until said criteria are universally satisfied. This solution-checking strategy is performed in an a posteriori fashion once the spatial discretisation and fluxes have been computed and the solution is advanced in time. Essentially this a posteriori strategy is more capable of detecting solutions that violate certain criteria, compared to an a priori limiting approach. However, once a candidate solution is deemed invalid for a cell (based on the selected criteria) and prompts the reduction in spatial order of accuracy, the neighbours of this cell also need to have their solutions recomputed with a lower-order approximation from the considered cell. In other words the fluxes computation have already contaminated the neighbours of the considered cell with a solution which is not valid and we need to recompute the solution at the considered cell and its neighbours.
The first distinctive drawback of this technique is that revisiting (i.e. reconstruction, fluxes, solution update) the cells without a valid solution and their corresponding neighbours has a significant computational footprint. In the parallel nature of any numerical framework, this corresponds to an additional layer of communication across processes necessary for capturing neighbouring cells that might have been contaminated by invalid solutions from neighbouring cells belonging to another process.
The second drawback is that the definition of solution validity criteria are either non-adaptive for various types of schemes, or that they are following the same principle as a priori limiters where the candidate solution is compared against the minimum and maximum variables of the solution of the neighbouring elements. The latter being linked with a fine balance between monotonicity and high-order of accuracy, especially when taking into account that for high-order FV schemes it is not just the direct-side neighbours that contribute to the reconstruction process but all the elements in the reconstruction's stencil. Omitting the values from all these elements from the definition of the bounding limits can have a profound effect on the sensitivity and accuracy of the scheme [52].
The scope of the present study is to improve the computational performance of the MOOD approach for finitevolume methods using unstructured meshes. Unstructured meshes coupled with the MOOD techniques present an excellent opportunity to enhance the robustness of numerical frameworks that are regularly employed for industrialscale engineering applications. More importantly we are also presented with a window of opportunity to study the numerical characteristics of finite volume and unstructured meshes for problems not traditionally attempted with such as instability-driven flows.
Our first goal in reducing the computational footprint of the MOOD approach is directly linked to reducing the number of cells with an invalid solution. To pursue this instead of employing an unlimited reconstruction process, a WENO scheme will be used since apart from the solution validity criteria there is the potential of obtaining fewer cells with solution violations compared to an unlimited scheme. The classical WENO scheme offers an a priori limiting, facilitated by nonlinear weights compounding seven (in 2D) or nine (in 3D) polynomial reconstructions per cell per variable to yield a desirable reconstruction which is essentially-non-oscillatory [53]. This makes the WENO paradigm computationally expensive, and often a very-high order WENO reconstruction may still suffer from oscillations. Unchecked these oscillations can propagate and contaminate the information carried by the neighbouring cells. A WENO scheme, while non-oscillatory in principle, may be plagued by minute oscillations which can contaminate the solution in one way or another. This is especially true for very high-order schemes, and/or superfine meshes, not to mention that an unstructured solver groups more elements per stencil than a structured one, making propagation of even the smallest oscillation undesirable, justifying the use of a cascade of auxiliary schemes. It is of prime interest to reliably detect whence the solution is affected, and recover some of its original accuracy wherever possible. The degeneration of the solution within the computational domain is distinguished by threshold values declaring the candidate solutions as valid (smooth and physical) or not.
The second goal is the investigation of the performance of WENO schemes for instability driven flows such as Richtmyer-Meshkov (RM) flows. The existence of perturbations along fluid interfaces like those developing from the Richtmyer-Meshkov instability form explicitly by the interaction of shock-waves with dints along material interfaces [54]. The driven flows undergo linear and non-linear evolution, and transiently transcend to turbulence. The influence of different initial conditions, material properties, and geometric configurations on the motion of the interface are investigated. Of particular interest is the link between the aggravated interface perturbations caused by the unstructured grid (and possibly the multidimensional reconstruction) itself and the breakdown of symmetry of structures at a much coarser resolution compared to structured grid results, and how the MOOD technique can influence that. The mesh itself, the numerical dissipation of the scheme, and the information propagation between the discretised terms are factors usually left unaccounted for. The UCNS3D [55,56] solver employed in this study by   [56] provides a robust OpenMP-MPI parallel framework for the multi-dimensional unsplit FV numerical computation. This permits the simulation of the late-time mixing of, and the strongly-nonlinear evolution of single-and multi-mode RMIs, and Kelvin-Helmholtz instabilities (KHIs), in all spatial dimensions. A certain limitation of the solver is that the adiabatic index is assumed constant for both fluids by the solver, making it a single component or "single-gamma" which can quantitatively affect the solution. The mixing process may not be severely affected by this, but the implementation may mispredict quantities of the flow fluid like the transmitted and reflected shock velocities, the time of reshock, or the velocity of the initialised perturbed interface. Addressing these will be the subject of future work.
To the best of our knowledge this is the first published work which investigates the effects of an unstructured solver for such flow problems and the first that while doing so employs a MOOD technique aimed at improving the robustness of the numerical frameworks. The outline of this paper is as follows: In Section 2 we describe the numerical framework in detail, in particular we outline briefly the reconstruction, fluxes and time advancement operators. The revised MOOD operator is detailed in the same section. In Section 3 we present several typical unsteady two-and three-dimensional test problems of computational fluid dynamics (CFD). Our computational results indicate that WENO-MOOD coupling could improve the robustness of the employed framework while decreasing the computational footprint of the MOOD approach, and it may be potentially applied for industrial-scale LES on unstructured meshes in the near future. The conclusions drawn from the present study are summarised in Section 4.

General computational framework
The conservative form of the three-dimensional Euler equations considered reads: ⎛ Here, the velocity components are u, v, and w, and the total energy E is per unit volume. The total energy in three dimensions is similar to its 2-D counterpart, the energy per unit area; in both cases, it is the sum of the internal energy e and the kinetic energy E = ρe + 0.5ρ(u 2 + v 2 + w 2 ), where e = p(γ − 1) −1 ρ −1 . Pressure p and the conservative quantities are related by an equation of state (EOS). For simplicity, in this work the sound speed c s = √ γ p ρ is used. For the two-dimensional formulation, the z component of the velocity and the associated terms are eliminated, resulting in a simpler expression for both the momentum and the energy equations.

Numerical framework
The discretisation in a domain Ω is achieved by combining conforming arbitrary shaped elements of volume |V i | which can be quadrilateral, triangles, tetrahedra, hexahedra, prisms or pyramids. Integrating Eq. (1) over a mesh element using the finite volume formulation the ordinary differential equation is obtained as follows: where U i is the volume averaged conserved variable vector, F n,l is the flux vector in the direction normal to the face l, N f is the number of faces per element, N q p is the number of quadrature points used for approximating the surface integrals. |A l | is the surface area of the corresponding face, and α corresponds to different Gaussian integration points x α and weights ω α over the face. The interface fluxes are computed based on the boundary extrapolated reconstructed values, which are obtained by a polynomial reconstruction from element-averaged data.

Spatial discretisation
The present reconstruction algorithm is based upon the approaches of [5,6,10], that has been successful applied various challenging flow problems [9,[57][58][59][60][61][62][63][64][65]. The reconstruction is performed in a reference space (ξ, η, ζ ) in order to remove the scaling effects in stencils consisting of various elements sizes, and details regarding the decomposition and transformation can be found in the work of Tsoutsanis et al. [5]. The reconstruction goal is to build a high-order polynomial p i (x, y, z) of arbitrary order r , for each considered element V i that has the same average as a general quantity U i this yields: Employing the inverse mapping the considered element V i can be transformed to the element V ′ i in the reference coordinate system, and the spatial average of the conserved variable U i does not change during coordinate transformation.
A reconstruction stencil as shown in Fig. 1, is built by recursively adding neighbouring elements, consisting of M + 1 elements, in order to perform the reconstruction on the target element where the index m refers to the local numbering of the elements in the stencil, with the element with index 0 being the considered element i. The r th order reconstruction polynomial at the transformed element V ′ 0 is an expansion over local polynomial basis functions φ k (ξ, η, ζ ) as given by: where U 0 is a general quantity at the considered element i. a k are the degrees of freedom and the upper index in the summation of expansion K corresponds to the number of the degrees of freedom and is related to the order of the polynomial r by K = 1 6 (r + 1)(r + 2)(r + 3) − 1 for 3D and K = 1 2 (r + 1)(r + 2) − 1. For computing the degrees of freedom a k , a minimum of K element are required within the stencil, in addition to the target element. For improving the robustness of the algorithm we use M = 2 · K as described in [5,6,66].
The unknown degrees of freedom a k for each element m of the stencil are determined by the cell average of the reconstruction polynomial p(ξ, η, ζ ) being equal to the element average of the solution U m : Where V ′ m is the volume of the element m in the stencil, in the reference coordinate system. The basis functions φ k need to satisfy the constraint of Eq. (3). The basis functions φ k for all the elements in the stencil are defined as follows: The basis functions employed are based on the Legendre type polynomials. Denoting the integrals of the basis function k over the considered element m in the stencil, and the vector of right-hand side by A mk and b respectively as given by: and by rewriting the equations for degrees of freedom a k in a matrix form as The matrix A mk can be precomputed since it is based on the geometry of the stencil's elements, which remains unchanged during the simulations. The resulting linear system is solved by a QR decomposition by making use of a pseudo-inverse matrix as detailed in [56].

The WENO scheme
The WENO scheme employs a non-linear combination of various reconstruction polynomials from the central stencil and the directional stencils; which are shown in Fig. 1, where the weight of each polynomial is dependent to the smoothness of its solution, and it is based on the approaches of [6,18]. The polynomials are given as where s t is the total number of stencils. Substituting back to Eq. (6) for p s (ξ, η, ζ ), we obtain the following expression: Using the condition that the sum of all weights is unity, yields whereã k are the reconstructed degrees of freedom; and the non-linear weight ω m is defined as The smoothness indicator SI m is given by: where β is a multi-index, r is the polynomial's order, and λ m is the linear weight. The central stencil is assigned a large linear weight of λ 1 = 1000 and the value to prevent division by zero of ε = 10 −6 is used in this study, and D being the derivative operator. The reader is referred to [5,10] for the definition of geometrical sectors and a detailed explanation of the different set of geometrical conditions. The smoothness indicator is a quadratic function of the degrees of freedom (a s k ) and Eq. (15) can be rewritten as: where the oscillation indication matrix OI kq is given by: and can be easily precomputed and stored at the initialisation stage of the simulation.

The MUSCL scheme
In the MUSCL framework employed in this work the high-order variation of the solution within every cell is approximated by the polynomial arising from the central stencil only. The scheme can be written as: where U i j,α is the extrapolated reconstructed solution at face j, and at quadrature point α, U i is the value for the conserved variable of element i, and (ξ a , η a ) are the coordinates of the quadrature point at the i j edge. All polynomials, for all the faces and for all the quadrature points are then limited by the limiter φ i which is valid for cell i to prevent any spurious oscillations from contaminating the solution. In the present work the slope limiter of Barth and Jespersen [67] is used and the reader is refered to [67] for more details.

Fluxes approximation & temporal discretisation
For the evaluation of the intercell numerical flux function F n l the approximate HLLC (Harten-Lax-van Leer-Contact) Riemann solver of Toro [68] is employed. The solution is advanced in time by the explicit Total Variation Diminishing (TVD), Runge-Kutta 3rd-order method. The influence of the Riemann solver and time-stepping methods within the MOOD paradigm is a topic of further study in the future. A C F L of 0.8 is used for all the test-cases in the present study, unless otherwise noted. For all the test problems involved the initial condition is approximated with a 7th-order Gaussian quadrature rule, to ensure that the initial condition is of high-accuracy, and during the simulation the volume/surface/line integrals are approximated by Gaussian quadrature rule suitable for the order of polynomial employed.

The MOOD paradigm
The key factor here, is the use of a succession of auxiliary, robust, schemes that are increasingly more dissipative than the principal numerical scheme employed by the solver. The underlying schemes contrast with the original MOOD (MOOD • ), where a single scheme ranging from a single, arbitrary high-order to first order is used. In principle, decreasing the order of accuracy of the solution for one cell can start at any order, and may occur multiple times within a single RK step; nevertheless, this cascade concludes after a finite number of iterations, and has to be short, otherwise the overhead might eliminate any potential benefits offered by the a posteriori approach. The ideal cascade provides stability with minimal footprint, whilst the principal scheme which is substantially more expensive, offers superior accuracy, especially for the smooth parts of the solution. The end result is that the solution for a cell is a posteriori detected (at time level n + 1) as either valid for a polynomial degree greater than 0, or the zeroth degree is reached. The present the moderated MOOD framework uses three different schemes. For any spatial order higher than 2 a WENO scheme is employed, that can be cascaded down in one step to 2nd-order MUSCL scheme, and if the solution is still not valid cascade to a first-order upwind scheme. For instance if the user selects a 6th-order WENO scheme, if the solution for a cell is not valid it will cascade to 2nd-order MUSCL, and if the solution is still not valid it will drop to first-order upwind scheme. The key difference between the original MOOD • approach and the relaxed revised MOOD R approach is that the latter has cascade adaptive criteria for admissibility of solutions as it will be further detailed later.
Offering an a posteriori stratagem to the nuisance of limiting, the allure of this paradigm is to yield a candidate solution from the principal numerical scheme applied by the solver, and assess the validity of said candidate solution using a set of criteria for physical and numerical admissibility. Each criterion has a predefined sensitivity. The solution of the computational domain is done by a WENO scheme, as shown in the flowchart of Fig. 2. Computing a candidate solution U * the beginning of each Runge-Kutta step, a cluster of detectors evaluate the candidate solution against the most pertinent properties (positivity, monotonicity, etc.). Those cells that pass the tests are deemed valid, and have their solutions admitted to the next Runge-Kutta stage. MOOD loops the next RK step, and so on. Unless the margins of the detectors are broken through (in which case the solution fails), the solution advances to the next time step. Consequently a stricken cell along with its central stencil around it (its von Neumann neighbourhood) are marked as 'problematic' for the current RK stage, and have their solutions recomputed locally with a lower degree polynomial reconstruction. The utility of the direct-side neighbours is that the candidate value of a flagged cell remains within the confines of the local minimum and local maximum on the previous time step, and ensures that no local extrema form, preserving monotonicity.
The "MOOD R -loop" is layered with one buffer step before 0th order, and recalculates any high-order solution with the built in second-order MUSCL scheme, alongside an appropriate (preferably strict) limiter. Making use of Fig. 3. Flow-chart of the MOOD-augmented UCNS3D code. In its current configuration, the implemented WENO scheme of UCNS3D "piggybacks" the a posteriori MOOD algorithm. The switch of accuracy criteria of the PAD/NAD sensors for the parachuting and the bulletproof scheme is also shown as steps 14, 15a, and 15b. The steps with dark orange outlines refer to the transition from a higher order of accuracy to the 2nd order MUSCL scheme, the first auxiliary, while those outlined with red refer to parachuting the solution to the bulletproof first order accurate Godunov's flux. A solution decremented to first order is exempted of the PAD/NAD due to its conservativity property.
only one central stencil 1 the process reconstructs a new second-order accurate candidate solution, which is tested anew for its physical and numerical validity. As a hedge, the numerical detector (19) is run on the MUSCL scheme with tempered numerical criteria (21), to tamp the aggressiveness of the previous criteria (20). The premise of this relaxation is that the 2nd order MUSCL scheme is more robust and dissipative that a higher order WENO reconstruction especially when a strict limiter is applied, thus less prone to oscillations. If the new candidate solution violates the bounds of a criterion again, the order is further demoted. This constitutes parachuting, the worst case scenario and end point of the cascade, which calls up the so called bulletproof scheme; the offending cell is updated with the robust, and stable, first-order accurate Godunov finite volume scheme and produces in principle valid (monotonicity-and positivity-wise) solutions.
Concordantly a candidate solution can be evaluated, and can detect the cluster of failed cells from the preceding tests. In principle, blending the MOOD loop within the framework of the UCNS3D is doubly interesting, from a CPU time and memory management standpoint; whilst the WENO scheme maintains the global high-order accuracy of the solution, the MOOD R contingent engineers a more efficient and robust framework, alleviating the Weighted Essentially Non-Oscillatory scheme from inheriting non-physical or numerically degenerate information, and thus halting in situ their propagation throughout the domain. Additionally, pitfalls of the a priori criteria by the WENO and MUSCL schemes are effectively canceled-out. The integration of the MOOD loop within the UCNS3D entails the development of robust detectors, capable of discerning the non-smooth solutions and shock regions; see Fig. 3 for more information on its incorporation into the UCNS3D, and its interactions therein. The concept presented in this work may not be restricted to the concert of numerical schemes found here, or their accuracy, or compactness thereof, since the applied principle can emerge as combination of any (f.e.) finite volume scheme.

MOOD admissibility criteria
An in depth analysis of the usually employed detectors can be found in [69], here we will only describe the ones we employ in our formulation, outlined in Fig. 3. The collection of points V i represents the set of first neighbours of the point in consideration. The cells that do not fulfil any of the employed detectors, are recalculated along with a stencil around them, as explained in the previous section, and illustrated in Fig. 3. In its present form, the NAD criterion assesses the conservative variables' vector as suggested by [3]: 1. Physical Admissible Detector (PAD): This detector checks that the solution is physical, that is, all points must have positive density and positive pressure at all times. In practice, this detector also identifies points with NaN values. 2. Numerical Admissible Detector (NAD) [48]: relaxed version of the Discrete Maximum Principle (DMP) [20].
It checks that the solution is monotonic and new extrema are not created. It compares the candidate solution with the solution obtained in the previous Runge-Kutta step. We remark that the superscript ( n ) indicates here the previous RK step, not the previous time step: The original DMP-relaxed margins, employed as introduced by [50], read: while the tempered relaxed criteria, used to dampen the transition from second to first order of accuracy were ad hoc selected as:

Applications
We present the numerical simulations employed to assess the performance of the WENO schemes with the relaxed MOOD approach in terms of robustness, accuracy and computational efficiency for the solution of the Euler equations in 2D setups. Several benchmark test problems have been performed and are listed below: • 2D Vortex Evolution. This test problem provides an assessment of the accuracy and computational footprint of the method.
• 2D Cylindrical Explosion. This test problem provides an assessment of the non-oscillatory properties of the method.
• 2D Double Mach Reflection. This test problem provides an assessment of the accuracy and robustness of the method in resolving all the flow features presented in this problem without any artefacts.
• 2D Implosion. This test problem is employed to understand the evolution of the nonlinearities for different Atwood numbers, and the onset of activation of the MOOD technique.

2D vortex evolution
The 2D vortex evolution test problem of Balsara and Shu [70] is used to assess the impact of the MOOD paradigm in a smooth flow problem. An isentropic vortex propagates at supersonic Mach number at 45 • across the domain. The computational domain is given by [0, 10] × [0, 10] with periodic boundary conditions applied on all sides. The unperturbed domain has an initial condition (ρ, u, v, p) = (1, 1, 1, 1), where temperature and density are defined as T = p/ρ, and S = p/r γ and the vortex perturbations are given by: The vortex strength ϵ = 5 and the adiabatic gas constant γ = 1.4. The e L 2 and the e L ∞ errors are computed as follows:   Mesh points where U c ( x, t f ) and U e ( x, t f ) are the computed and exact solutions at the end of the simulation t f . The exact solution U e ( x, t f ) is given by the initial condition itself at t 0 . Hybrid unstructured meshes consisting of quadrilateral and triangular elements are used for this test problem of 16, 32, 64, and 128 nodes per side resolution, and the simulation is run for a time of t f = 10 s with 3rd-order, 4th-order and 5th-order WENO schemes with the original and the relaxed MOOD algorithm.
From the results shown in Table 1, it is evident that the WENO schemes without the MOOD algorithm achieve convergence rates close to their theoretical ones. When employing the MOOD algorithm in its original formulation the computational cost increases and the convergence rates are far from the theoretical ones. The accuracy gets even worse as the order of the scheme increases since there are more cells that their solution is invalid and their spatial order is reduced to first-order. However the relaxed MOOD algorithm has a significantly smaller computational footprint compared to the original MOOD approach, and at the same time it recovers the theoretical order of accuracy of the WENO schemes, since there are fewer cells that will drop their solution to 1st-order. Additionally from Fig. 4 it can be noticed that the relaxed scheme has fewer flagged cells, and that their lowest spatial order is 2nd-order compared to the original MOOD scheme where the majority of the flagged cells are 1storder. Therefore it is encouraging that this relaxed MOOD algorithm for smooth flow problems has a significantly reduced computational footprint (2.5%-12% increase in computational time) compared to the original MOOD approach (7%-34% increase in computational time), and more importantly that it allows the WENO schemes to achieve their theoretical order of accuracy.

Cylindrical explosion
This test is Sod's shock tube problem, a form of cylindrical explosion. The 2-D Euler equations are solved in a circular domain 2.0 × 2.0 in the x-y plane, using an Advancing Front unstructured mesh of 223,455 triangularquadrilateral elements and a solution time of t = 0.2 s. Two fluid regions are separated at t = 0 s by a circular discontinuity [1] and the initial conditions read: (1, 0, 0, 1), r ∈ (0, 0.5), (0.125, 0, 0, 0.14), r ∈ (0.5, 1), Fig. 4. Snapshot of flagged cells coloured by the permitted spatial order of accuracy using the original (left) and relaxed (right) MOOD criteria, in an unstructured mesh for the 2D vortex evolution problem at t = 10 s using a WENO 5th-order scheme. It can be noticed that the relaxed scheme has fewer flagged cells, and that their lowest spatial order is 2nd-order compared to the original MOOD scheme where the majority of the flagged cells are 1st-order.
Two families of schemes are considered in this test problem, the WENO and the unlimited one. The unlimited one refers to the scheme using only one central stencil for reconstruction without WENO or MUSCL at any stage during the cascade of the MOOD. The main purpose for testing this setting was to assess if a much simpler and computationally cheaper scheme compared to WENO could be employed for these types of flow problems.
In agreement with published data by Toro 2009 [1], the results exhibit a circular shock wave that recedes from the origin (0, 0). A companion circular contact radiates from the origin in the form of a continuous surface. The contact discontinuity shows some local oscillations, as seen in Fig. 5). Finally, a circular rarefaction can be discerned, which appears to retreat towards the origin. In Fig. 5 the results for WENO4 and WENO6, with the MOOD • and MOOD R versions of the same algorithm enabled are presented. Lest the reader forgets, the MOOD • has the original cascade criteria applied to both the high-order to MUSCL transition, and the MUSCL-to-Godunov switch, which leads a lot more elements to Godunov's flux; its relaxed counterpart retains the MUSCL scheme. Since the aforementioned minor oscillations persist even if they are smoother, they provide an indirect assertion that the MOOD algorithm restricts itself with a highly localised footprint. It can be seen in Figs. 7 and 6, that the WENO scheme is also prone to some very low-amplitude oscillations which is to be expected for this test problem and resolution employed (which they can be further reduced by tuning the ε parameter for preventing division by zero), and this is true for the original and relaxed MOOD criteria. The unlimited scheme with the relaxed MOOD technique on the other hand has significantly more oscillations compared to the WENO schemes as highlighted in Fig. 7, and therefore is not considered suitable for these types of problems.

Double Mach reflection
The deflection of a shock wave from a wedge can give rise to a series of interesting and complex flow phenomena. The most important part of the simulation, that is also useful for the evaluation of different methods, is the generation of a dense jet along the lower wall. According to Woodward and Collela (1984) [71], and Liska and Wendroff (2003) [72] this structure, akin to a shaped charge jet is very sensitive to the numerical diffusion of a scheme. By definition, the wedge located upstream the shock wave is inclined at an angle θ . While such a "ramp" might be of special interest in specialised applications like ramjets, or easier to set up in a real world experiment, it requires plenty of computational work for coding the appropriate boundary conditions [1]. In this simplified condition used, the shock wave is given an angle θ . The initial conditions involve a moving shock wave with Mach number M = 10 at an inclination of θ = 60 • . The conditions ahead of the shock is at rest with uniform density and pressure ρ = 1.4 Fig. 5. Close up of the shock front region of Sod's 2D shock-tube problem. The higher-order WENO6 scheme presents a similar pattern of low mode oscillations as WENO4, albeit with a slightly larger amplitude. MOOD R marks significantly fewer elements than MOOD • ; an even smaller number of cells was reconstructed to first order using Godunov's scheme.  The WENO 3rd-order scheme is used for the present test problem, with and without the relaxed MOOD algorithm.
In the domain two Mach stems form together with two contact discontinuities. The first contact appears well resolved, and all triple points are shown. At the same time, the second Mach reflection is somewhat underresolved, and gradually diminishes upon reaching the contact discontinuity generated by the first Mach shock. Wrinkly islands in the domain are permitted oscillations in the solution forming as the scheme reacts to machine noise, indicating that the schemes have a tendency to oscillate. In all schemes, the Mach stem appears diffracted close to the tip of the vortex "cocoon". For the WENO3, the MOOD R has a remedial effect to carbuncle, as it can be seen in Fig. 8. At the same time, the fact that the MOOD R is not fully suppressing carbuncles demonstrates that coupling the algorithm with the well-established WENO and MUSCL frameworks is not making additional and potentially unwarranted physical assumptions, since the phenomenon is not a mere failing of numerics (n.b. mass-flux corrections may help) may be a special class of entropy solutions, and valid vanishing viscosity limits to the Euler equations [73], since it has an intrinsic correspondence with an instability mechanism of the Euler equations [74,75].
The cascade capture the vortex generated in the primary slip line area, in agreement to [71,76] and the jet is well represented and indicates the pressure gradient is building up in the region where the contact discontinuity nears the wall. As expected the results are in agreement with the higher-resolution results obtained with the WENO7 schemes on a finer mesh with UCNS3D [66], although the shear-waves are not captured with the same resolution of the WENO7 scheme as seen in Fig. 9.
An additional, yet minor, source of error is also detected at the very first time step of the solution. The schemes branch off the initial shock, which immediately creates a structure like a shadow trailing the real shock. This "starting-error" is generated as the numerical representation of the shock has a latency to catch up with the steady shock profile, resulting in an "overheated region" [71]. The only observable trail of this error in the domain can be seen atop the peak of the reflected shock.
Contact discontinuities, which have special importance in this test case, can be readily surveyed in the density contour plots of Fig. 8. The more interesting and elaborate structures are confined in the restricted region of the double Mach reflection [71,76].   9. Zoomed on the wave interaction zone for the double Mach reflection problem, where the key differences can be noticed between the WENO3 with a posteriori MOOD R -limiting, and the results obtained with a WENO 7th-order scheme on a finer mesh as reproduced from [66]. The structures related to the shear waves appear more pronounced in the finer-resolution results with WENO7 as expected, and the figure on the right highlights the locations of the cells for which the MOOD R -limiting approach has reduced their spatial accuracy.

2D implosion at different Atwood numbers
Another class of numerical setups for the RMI study the instability in converging geometries. Relatively recent examples of this approach is the work of [77] and [78]. Our study focuses on Case 2 of the former work which is also tested in more detail in the latter publication. The setups tested in this work are almost identical to those of [78], but with M = 1.22 as seen in Fig. 10. Four different Atwood numbers were tested: A = −0.33, −0.67, −0.82, −0.90: (1.3764, −0.4021, 1.5698) r < 1.5, (1, 0, 1) , r ∈ [1 + α 0 · sin(nθ ), 1.5], (0.2, 0, 1) , r < 1 + α 0 · sin(nθ ); An unstructured mesh was created out of two mesh domains; an interior domain, extending to 1 cm radius, bounded by two semi-circular connectors with 901 nodes apiece, and a exterior mesh, radiating to r = 3 cm, enclosed by two boundaries with 301 grid points each. The inner, refined region was meshed with the Advancing Front algorithm for quadrilateral cells, and the exterior topology utilised the Advancing Front Orthogonal algorithm with both triangular, and quadrilateral elements, and a boundary decay of 0.99, with a reflective wall boundary condition applied at the outer boundary. The WENO 5th-order scheme was used with the relaxed MOOD algorithm.
Nonlinearity affects the converging implosions more than their planar counterparts. The impinging shock bifurcates into a transmitted shock which propagates within the lighter fluid and an outbound rarefaction wave. Phase inversion, as expected, occurs immediately after the interaction [77]. Fig. 11 shows that the briefly symmetric spikes and bubble structures rapidly lose their symmetry for higher Atwood numbers. The nonlinear evolution is in agreement with the results of [78], and the reader is referred to the two animations included for this test-problem Implosion1, and Implosion2 which depict the evolution of density and density gradient magnitude respectively for all the Atwood numbers. The grid resolution offers a detailed vista for implosion and reshock, and the time-integrated histories of the kinetic energy shown in Fig. 12 demonstrate the time of peak compression (where KE is minimum).
The comparison between the four different Atwood numbers in Fig. 11 shows that there is indeed a correspondence between the Atwood number and the penetration depth of the lighter fluid by the spikes. The larger the density ratio, the more elongated and finer the spike. The MOOD R -coupled solver again proves efficient, with an average MOOD R activation rate of 0.03%. The peaks between t = 2 s and t = 3 s in the time-series of Fig. 13 correspond to KE ramp-ups of Fig. 12, highlighting that the MOOD algorithm is an essential component for increasing the robustness of the framework at the required instances. From Fig. 12 where the evolution of the spike position is plotted, the trends obtained agree with the ones noticed in [78]. As expected the use of higher-grid resolution and less dissipative methods employed in the present study compared to [78] results in increased spike growth after the linear regime (up to t = 0.6 s), and they are driven further into the quiescent fluid, which agrees with previous studies [79].  11. Left: Density contour plots at time t = 3 s for all four Atwood numbers, showing the spike penetration with increasing density ratio. The ratio increases from top left (denser inner gas) to bottom right (least dense inner gas). The surrounding density is kept the same, and the motion of the interface asymptotically converges to the solid-vacuum case. Right: Density gradient magnitude for the same time and Atwood numbers. The formation of the KHI is readily apparent, and the finer the spikes, the easier to break up.

Conclusions
In this work we have developed a new relaxed MOOD algorithm, whose solution admissibility criteria are adaptive through the cascade of spatial order reduction. This algorithm has been employed with the WENO and MUSCL FV schemes on unstructured meshes. The key ingredient of the algorithm is that the solution validity criteria are relaxed once a 2nd-order MUSCL reconstruction is reached based on the Barth & Jespersen limiter, and as a consequence the number of cells that reduce their spatial order to 1st-order of accuracy is greatly reduced. Therefore the computational footprint of the method is significantly smaller compared to the original MOOD algorithm. This framework has been implemented within the 2D/3D MPI+OpenMP parallel UCNS3D [55,56] code tailored for unstructured meshes. We have shown on a series of test problems using mixed element unstructured meshes that the present WENO+MOOD R framework: • attains and maintains the theoretical order of accuracy of the WENO schemes for smooth flow problems, while reducing the computational cost compared to the original MOOD algorithm • improves the non-oscillatory properties of the employed WENO FV method targeting only the affected regions • reduces the artefacts in the solutions (such as carbuncle) that would otherwise lead to non-physical results • improves the robustness of the numerical framework at various locations of the test problems examined such as waves emitted by shock-interface interactions It should be noted that a very small cell number has to drop their order of accuracy for realistic test-problems (generally 0.03%). This is crucial for reducing the computational footprint and overhead associated with the MOOD algorithm, especially in a parallel CFD code. In the future we plan to extend the MOOD R algorithm, to viscous Navier-Stokes equations and in particular for iLES turbulent simulations of compressible flows. In such cases an even more affordable version of the WENO schemes such as the compact WENO (CWENO) might be needed, to further reduce the computational cost of the framework.