High Order ADER Schemes for Continuum Mechanics

In this paper we first review the development of high order ADER finite volume and ADER discontinuous Galerkin schemes on fixed and moving meshes, since their introduction in 1999 by Toro et al. We show the modern variant of ADER based on a space-time predictor-corrector formulation in the context of ADER discontinuous Galerkin schemes with a posteriori subcell finite volume limiter on fixed and moving grids, as well as on space-time adaptive Cartesian AMR meshes. We then present and discuss the unified symmetric hyperbolic and thermodynamically compatible (SHTC) formulation of continuum mechanics developed by Godunov, Peshkov, and Romenski (GPR model), which allows to describe fluid and solid mechanics in one single and unified first order hyperbolic system. In order to deal with free surface and moving boundary problems, a simple diffuse interface approach is employed, which is compatible with Eulerian schemes on fixed grids as well as direct Arbitrary-Lagrangian-Eulerian methods on moving meshes. We show some examples of moving boundary problems in fluid and solid mechanics.

In this paper we first review the development of high order ADER finite volume and ADER discontinuous Galerkin schemes on fixed and moving meshes, since their introduction in 1999 by Toro et al. We show the modern variant of ADER based on a space-time predictor-corrector formulation in the context of ADER discontinuous Galerkin schemes with a posteriori subcell finite volume limiter on fixed and moving grids, as well as on space-time adaptive Cartesian AMR meshes. We then present and discuss the unified symmetric hyperbolic and thermodynamically compatible (SHTC) formulation of continuum mechanics developed by Godunov, Peshkov, and Romenski (GPR model), which allows to describe fluid and solid mechanics in one single and unified first order hyperbolic system. In order to deal with free surface and moving boundary problems, a simple diffuse interface approach is employed, which is compatible with Eulerian schemes on fixed grids as well as direct Arbitrary-Lagrangian-Eulerian methods on moving meshes. We show some examples of moving boundary problems in fluid and solid mechanics.

INTRODUCTION AND REVIEW OF THE ADER APPROACH
The development of high order numerical schemes for hyperbolic conservation laws has been one of the major challenges of numerical analysis for the last decades. Godunov [1] proved that for the linear advection equation no monotone linear schemes of second or higher order of accuracy can be constructed. Therefore, even if physical viscosity is considered, a linear high order scheme will present spurious oscillations near discontinuities, as it can be seen, for instance for the Lax-Wendroff scheme, Lax and Wendroff [2]. A first idea to circumvent this theorem has been proposed in Kolgan [3], where limited slopes are employed to produce a non-linear scheme of second order of accuracy in space. Since then, many high order numerical methods have been developed like the Total Variation Disminishing methods (TVD) and Flux limiter methods (see, for instance, [4][5][6][7][8][9]). Despite these methodologies being already well-established at the end of the last century, their major drawback was that they just provided global second order of accuracy and reduced locally to first order in the vicinity of smooth extrema.
More advanced non-linear methods for advection dominated problems involve the family of ENO and WENO schemes, see Harten and Osher [10], Harten et al. [11], and Shu [12]. In particular, the method of Harten et al. [11] is a fully discrete high order scheme that can be re-interpreted in terms of the solution of a generalized Riemann problem (GRP), see Castro and Toro [13]. Moreover, it can be seen as a generalization of the MUSCL-Hancock method of van Leer, see van Leer [8], Toro [9], and Berthon [14].
Following the idea of solving a generalized Riemann problem (GRP), see also Ben-Artzi and Falcovitz [15], LeFloch and Tatsien [16], Ben-Artzi et al. [17], and Han et al. [18], the ADER approach (Arbitrary high order DErivative Riemann problem) has been first put forward for the linear advection equation with constant coefficients by Millington et al. [19] and Toro et al. [20]. The first step of the methodology involves piecewise polynomial data reconstruction, where a non-linear ENO reconstruction is applied in order to avoid spurious oscillations of the numerical solution. Then, a GRP is defined at each cell interface. Classically, the initial condition for the GRP was given as piece-wise linear polynomials and second order schemes could be obtained by constructing a space-time integral of the solution in an appropriate control volume [21,22], or following a MUSCL approach, van Leer [23] and Colella [24]. An alternative methodology proposed in Ben-Artzi and Falcovitz [25] consists in expressing the solution of the GRP as a Taylor series expansion in time. The ADER approach obtains the high order time derivatives of the GRP solution at the cell interface via the Cauchy-Kovalevskaya procedure, which replaces time derivatives by spatial derivatives using repeated differentiation of the differential form of the PDE. The spatial derivatives, which may also jump at the interface, are defined via the solution of linearized Riemann problems for the derivatives, where linearization is carried out about the Godunov state obtained from the classical Riemann problem between the boundary extrapolated values at the interface. In Figure 1, the classical piece-wise constant polynomials are plotted against a high order reconstruction and the similarity solutions for both cases are sketched. Finally, these similarity solutions are used to construct the numerical flux. The resulting schemes are arbitrary high order accurate in both space and time, in the sense that they have no theoretical accuracy barrier.
Since their introduction in Toro et al. [20] and Millington et al. [19], many extensions of the ADER methodology have been proposed. Regarding 2D linear PDEs, one may refer to Schwartzkopff et al. [26] and their simplification for the particular case of structured grids in Schwartzkopff et al. [27]. Moreover, non-linear systems have been initially addressed in Toro and Titarev [28] and Titarev and Toro [29]. Further applications of ADER on non-Cartesian meshes have been presented in Käser [30], Käser and Iske [31], Dumbser et al. [32], and Castro and Toro [13]. One should also mention the development of ADER schemes in the framework of discontinuous Galerkin (DG) finite element methods, see Qiu et al. [33], Dumbser and Munz [34] and Gassner et al. [35]. One of the main advantages of using DG is that the reconstruction step of classical ADER finite volume (ADER-FV) schemes can be skipped, since the discrete solution is already given by high order piecewise polynomials that can be directly evolved during each time step. Furthermore, ADER-DG schemes avoid the use of classical Runge-Kutta time stepping and thus provide efficient communication-avoiding schemes for parallel computing, see Fambri et al. [36] and allow for simple and natural time-accurate local time stepping (LTS), see Dumbser et al. [37].
An important step forward in the development of more general ADER schemes was achieved in Dumbser et al. [38], where a new class of ADER-FV methods has been introduced. The main contribution of this paper consists in the introduction of a new element-local space-time DG predictor, which allows at the same time the treatment of stiff source terms, as well as the replacement of the cumbersome Cauchy-Kovalevskaya procedure. First, a high order WENO method is employed to compute a polynomial reconstruction of the data inside each spatial element; then, an element-local weak formulation of the conservation law is considered in space-time and the predictor is applied to construct the time evolution of the WENO polynomials within each cell. Note that, in this step, the integration by parts is performed only in time, which differs from global space-time DG schemes [39,40], which are globally implicit. Finally, the cell averages are updated with an explicit fully discrete one-step scheme, considering the integral form of the equations. As a result, the proposed methodology maintains arbitrary high order of accuracy, while avoiding the issues related to the use of a Taylor series expansion in time. As already mentioned above, it naturally provides an approach for the treatment of stiff source terms [for further details on this topic, see [41] and references therein].
The above methodology can also be applied in the discontiuous Galerkin framework as presented in Dumbser et al. [42], where, a unified P N P M framework for arbitrary high order one-step finite volume and DG schemes has been introduced. For other reconstruction-based DG schemes, see e.g., Luo et al. [43,44]. Afterwards, the methodology has been extended to solve a wide variety of different PDE systems, such as the resistive relativistic MHD equations, Dumbser and Zanotti [45]; non-conservative hyperbolic systems found in geophysical flows, Dumbser et al. [46] in which a wellbalanced and path-conservative version of the scheme has been developed; compressible multi-phase flows Dumbser et al. [47], the compressible Navier-Stokes equations, Dumbser [48]; the compressible Euler equations and divergence-free schemes for MHD, Balsara et al. [49], and Balsara and Dumbser [50], where ADER schemes were used in combination with genuinely multidimensional Riemann solvers. The last extensions concern the special and general relativistic MHD equations, see Zanotti et al. [51], and Fambri et al. [36], as well as the Einstein field equations of general relativity [52,53].
Later, ADER schemes have been extended to adaptive mesh refinement on Cartesian grids (AMR), in combination with time accurate local time stepping (LTS). This technique has initially been introduced in Dumbser et al. [54,55] for conservative and non-conservative hyperbolic systems, respectively. Moreover, the schemes of the ADER family were the first high order methods to be applied for the numerical solution of the unified first order hyperbolic formulation of continuum mechanics by Godunov, Peshkov and Romenski [56][57][58], see Dumbser et al. [59][60][61]. In the rest of this paper, we will refer to the Godunov-Peshkov-Romenski model of continuum mechanics as GPR model.
The ADER approach has also been extended to the direct Arbitrary-Lagrangian-Eulerian framework (ALE), where the mesh moves with an arbitrary velocity, taken as close as possible to the local fluid velocity. Initially developed for one space dimension, it has been soon extended to the case of the two and three dimensional Euler equations on unstructured meshes, Boscheri and Dumbser [62,63], including the discretization of non-conservative products. Further works in this area involve the use of local timestepping techniques, [64,65]; coupling with multidimensional HLL Riemann solvers, Boscheri et al. [66]; solution of magnetohydrodynamics problems (MHD), [67,68]; development of a quadrature-free approach to increase the computational efficiency of the overall method, Boscheri and Dumbser [69]; use of curvilinear unstructured meshes, Boscheri and Dumbser [70]; or extension to solve the GPR model, Boscheri et al. [71] and Peshkov et al. [72]. Furthermore, in Gaburro et al. [73] a novel algorithm to deal with moving non-conforming polygonal grids has been presented. The methodology reduces the typical mesh distortion arising in shear flows and provides high quality elements even for long-time simulations. An exactly well-balanced path-conservative version of this approach for the Euler equations with gravity can be found in Gaburro et al. [74]. Still in the ALE framework, within this article, we will present new results for the family of ADER-FV and ADER-DG schemes on moving unstructured Voronoi meshes [75], as recently introduced in Gaburro [76] and Gaburro et al. [77].
It is well-known that when dealing with high order schemes special care must be paid to the limiting methodology employed. In most of the previous referenced papers classical a priori limiters have been used, such as WENO reconstruction. Nevertheless, some alternative contributions to this topic can be found in the series of papers [51,[77][78][79][80][81][82][83][84][85], where a novel a posteriori sub-cell FV limiter of high order DG schemes, based on the MOOD paradigm of Clain et al. [86] and Diot et al. [87,88], has been employed.
Besides the references given above, which focus on the development of the ADER methodology with a local space-time Galerkin predictor, many recent papers have been devoted to the development of other families of ADER schemes, like the classical ADER finite volume methods. Without pretending to be exhaustive, we may refer to Castro et al. [ [101], and Dematté et al. [102] and references therein.
In this paper, as a promising application of the family of ADER schemes, we solve a diffuse interface formulation of the GPR model of continuum mechanics. In comparison with existing continuum mechanics models, the novel feature of the GPR model is in that it incorporates the two main branches of continuum mechanics, fluid and solid mechanics, in one single unified PDE system. Recall that traditionally fluid and solid mechanics are described by PDE systems of different types, i.e., parabolic (viscous fluids) and hyperbolic (linear elasticity and hyperelasticity), which imposes many theoretical and technical difficulties if one wishes to model natural and industrial processes involving co-existence of the fluid and solid states such as in fluid-structure interaction (FSI) problems, modeling of general solid-fluid transition such as in melting and solidification processes, e.g., additive manufacturing, see for example Francois et al. [103], flows of granular media [104], viscoplastic flows, e.g., debris flows, avalanches, mantle convection, flows of many industrial Bingham-type fluids, see Balmforth et al. [105]. Due to the unified treatment of fluids and solids, the GPR model thus has a great potential for simplifying the modeling process and code development for solving the aforementioned problems. Yet, before to be applied to practical problems, the GPR model may require a coupling with an interface tracking/capturing technique for the modeling of moving material boundaries such as in free surface flows or solid body motion. In particular, in this paper, we couple the GPR model with a simple diffuse interface approach, see Tavelli et al. [85], Dumbser [106], Gaburro et al. [107], Kemm et al. [108]. For example, very interesting computational results with similar diffuse interface approaches and level set techniques for compressible multi-material flows have been obtained for example in Gavrilyuk et al. [109], Favrie et al. [110], Favrie and Gavrilyuk [111], Ndanou et al. [112], de Brauer et al. [113], Michael and Nikiforakis [114], Jackson and Nikiforakis [115], and Barton [116]. Finally, we demonstrate that the ADER family of schemes is capable to resolve the GPR model in both solid and fluid regimes.
The paper is organized as follows. In section 2 we present the family of ADER finite volume and ADER discontinuous Galerkin finite element schemes on fixed Cartesian and moving polygonal meshes in two space dimensions. Next, in section 3 we introduce the diffuse interface formulation of the GPR model. In section 4 we show some computational results obtained with different kinds of ADER schemes (ADER-FV and ADER-DG) on different mesh topologies, including moving unstructured Voronoi meshes, as well as fixed and adaptive Cartesian grids. The paper is rounded off by some concluding remarks and an outlook to future work in section 5.
Here, we briefly describe the key features of our numerical scheme, keeping the notation as general as possible, and referring to the literature for further details. We start by introducing the general form of our governing PDE system and a moving unstructured discretization of two-dimensional domains (sections 2.1 and 2.2); next, in section 2.3 we describe the data representation of the discrete solution. Then, we explain how to obtain high order of accuracy in space: this is available by construction in the DG case, and obtained via some variants of the well-known WENO procedure [32,[121][122][123][124][125] for the FV approach. Finally, we focus on the predictor-corrector version of the ADER scheme that allows to achieve arbitrary high order of accuracy in space and time. Since it is out of the scope of this paper to recall all the details, a general overview is given in sections 2.5 and 2.7, and an inedited proof of the convergence of the predictor for a non-linear conservation law is presented in section 2.6.
We would like to emphasize that, besides this novel convergence proof, other progress has been introduced within this work. Indeed, up to our knowledge, it is the first time that: (i) the ADER approach is used to solve a diffuse interface formulation of the GPR model that addresses the free surface problem in both solid and fluid mechanics context (previously, a similar formulation was used only in the solid dynamics context [112,126,127]); (ii) non-conservative products are taken into account in the high order direct ALE scheme of Gaburro et al. [77], where they have to be integrated also on degenerate spacetime control volumes (see section 2.5.2).

Governing PDE System
In this paper we consider high order fully-discrete schemes for non-linear systems of hyperbolic PDE with non-conservative products and algebraic source terms of the form where Q = Q(x, t) ∈ Q ⊂ R m is the state vector, t ∈ R + 0 is the time, x ∈ ⊂ R d is the spatial coordinate, d is the number of space dimensions, Q is the so-called state space or phase space, F(Q) is the non-linear flux tensor, B(Q)·∇Q is a non-conservative product and S(Q) is a purely algebraic source term. Introducing the system matrix A(Q) = ∂F/∂Q + B(Q) the above system can also be written in quasi-linear form as The system is said to be hyperbolic if for all n = 0 and for all Q ∈ Q the matrix A(Q) · n has m real eigenvalues and a full set of m linearly independent right eigenvectors. The system (1) needs to be provided with an initial condition Q(x, 0) = Q 0 (x) and appropriate boundary conditions on ∂ .
In this paper we focus on a particular, but very general, example of a first-order system (1) describing elastic and viscoplastic heat-conducting media; it will be discussed in section 3.

Domain Discretization
In the general ALE case, we consider a moving two-dimensional (d = 2) domain (t) and we cover it using an unstructured mesh made of N P non-overlapping polygons P i , i = 1, . . . N P . The mesh is first built at time t = 0 and then it is rearranged at each time step t n : elements and nodes are moved following the local fluid velocity and when necessary, in order to prevent mesh distortion, also the mesh topology (i.e., the shape of the elements and their connectivities) is changed.
Given a polygon P n i we denote by V(P n i ) = {v n i 1 , . . . , v n i j , . . . , v n i N n V i } the set of its N n V i Voronoi neighbors (the neighbors that share with P n i at least a vertex), and by E(P n i ) = {e n i 1 , . . . , e n i j ,. . . ,e n i N n V i } the set of its N n V i edges, and by D(P n i ) = {d n i 1 , . . . , d n i j ,. . . ,d n i N n V i } the set of its N n V i vertexes, consistently ordered counterclockwise. Finally, the barycenter of P n i is noted as x n b i = (x n b i , y n b i ). When necessary, by connecting x n b i with each vertex of D(P i ) we can subdivide a polygon P n i in N n V i subtriangles denoted as T (P n i ) = {T n i 1 , . . . , T n i j , . . . , T n The coordinates of each node at time t n are denoted by x n k , and V n k represents the velocity at which it is supposed to move, Frontiers in Physics | www.frontiersin.org so that its new coordinates at time t n+1 are given from the following relation More details on how to obtain V can be found in Boscheri et al. [68], Boscheri and Dumbser [63,119] for what concerns classical direct ALE schemes on conforming unstructured grids, in Gaburro et al. [73,74] for non-conforming unstructured grids, in Boscheri and Dumbser [70] for curvilinear meshes, and we refer in particular to section 2.4 and 2.5 of Gaburro et al. [77] for what concerns moving unstructured polygonal grids allowing for topology changes, which indeed is the ALE case considered in this paper (see case B below). Moreover, working in the ALE framework, we are allowed to take V = 0, i.e., we can also work in a fixed Eulerian system where the initial mesh is never modified.
In particular, in this paper we will consider the following two situations for our domain discretization: A. A fixed Cartesian mesh made of N P quadrilaterals elements, which is not moved during the simulation, but which can be successively refined, with a general space-tree-type data structure that allows element-by-element refinement with a general refinement factor r ≥ 2, in order to increase the resolution in the areas of interest, as can be seen in Figure 2 (for the details on the refinement procedure we refer to Dumbser et al. [54] and Fambri et al. [36]). To ease the description of the numerical method, we will associate to each quadrilateral element P n i , a set of indices that refer to its Cartesian coordinates, j, k , such that P n jk : B. A moving polygonal grid as the one described in Gaburro et al. [77] that (i) moves with the fluid flow in order to reduce the numerical dissipation associated with transport terms and (ii) also allows for topology changes at any time step in order to maintain always a high quality of the moving mesh; in this case we remark that our method is also able to deal with degenerate space time control volumes at arbitrary high order of accuracy.

Space-Time Connectivity
To better understand the context of moving meshes we refer the reader to Figure 3: note that the tessellation at time t n has been evolved resulting in a slightly different tessellation at time t n+1 ; for each element P n i the new vertex coordinates x n+1 k , k = 1, . . . , N n V i , are connected to the old coordinates x n+1 k via straight line segments, yielding the multidimensional space-time control volume C n i , that involves N n,st V i +2 space-time sub-surfaces. Specifically, the space-time volume C n i is bounded on the bottom and on the top by the element configuration at the current time level P n i and at the new time level P n+1 i , respectively, while it is closed with a total number of N n,st V i lateral space-time surfaces ∂C n i j , j = 1, . . . , N n,st V i that are given by the evolution of each edge e n i j of element P n i within the time step t = t n+1 − t n . A priori, ∂C n i j are not parallel to the time direction: thus to be treated numerically they can be mapped to a reference square by using a FIGURE 2 | Sketch of the mesh refinement structure of three AMR levels with refinement factor r = 3. Solid lines indicate active cells, whereas the dashed ones are the virtual cells allowing interpolation between the coarse and the refined mesh, needed in the case of high order WENO reconstruction.
set of of bilinear basis functions (see Boscheri and Dumbser [62]). To resume, the space-time volume C n i is bounded by its surface ∂C n i which is given by Note that in the fixed Cartesian case, C n i reduces to a right parallelepiped with four lateral space-time surfaces ∂C n i j parallel to the time-direction, so many simplifications are possible.
We close this part by emphasizing that the family of direct ALE schemes proposed in this work, based on the ADER predictor-corrector approach, is based on the integration of the governing Equation (1) in space and in time directly over these space-time control volumes, see section 2.7. Note that this procedure, which is more evident when C n i is an oblique prism, is also hidden when C n i is just a right parallelepiped.

Data Representation
The conserved variables Q in (1) are discretized in each polygon P n i at the current time t n via piecewise polynomials of arbitrary high order N, denoted by u n h (x, t n ) and defined as where in the last equality we have employed the classical tensor index notation based on the Einstein summation convention, which implies summation over two equal indices. The functions ϕ ℓ (x, t n ) can be either: i. Nodal spatial basis functions given by a set of Lagrange interpolation polynomials of maximum degree N with the property where {x m GL } are the set of the Gauss-Legendre (GL) quadrature points on P n i (see Stroud [128] for the multidimensional case). In particular, when employing these basis functions on a Cartesian grid, each quadrilateral P n i is easily mapped to a reference square, we only need the tensor product of the GL quadrature points in the unit interval [0, 1], and the ϕ ℓ are simply generated by multiplying one-dimensional nodal basis functions, i.e., with ϕ ℓ i satisfying (6) with d = 1, and x = x j− 1 2 + ξ x j , y = y k− 1 2 +η y k being the set of reference coordinates related to P n i . In this case, the total number of GL quadrature points per polygon, as well as the total number of basis functions {ϕ ℓ } and expansion coefficientsû n ℓ,i , the so-called degrees of freedom (DOF), is N = (N + 1) d . These basis functions are used on Cartesian grids, i.e., for Case A. ii. Modal spatial basis functions written through a Taylor series of degree N in the variables x = (x, y) directly defined on the physical element P n i , expanded about its current barycenter x n b i and normalized by its current characteristic length h i h i being the radius of the circumcircle of P n i . In this case the (N + m). We employ this kind of basis functions in the moving unstructured polygonal Case B.
The discontinuous finite element data representation (5) leads naturally to discontinuous Galerkin (DG) schemes if N > 0, but also to finite volume (FV) schemes in the case N = 0. This indeed means that for N = 0 we have ϕ ℓ (x) = 1, with ℓ = 0 and (5) reduces to the classical piecewise constant data that are typical of finite volume methods. In the case N > 0 (DG) the form given by (5) already provides a spatially high order accurate data representation with accuracy N + 1, where instead for the case N = 0 (FV), if we are interested in increasing the spatial order of accuracy, up to M + 1 for examle, we need to perform a spatial reconstruction. With this notation, our method falls within the more general class of P N P M schemes introduced in Dumbser et al. [42] for fixed unstructured meshes.

Data Reconstruction
In this section we focus on the reconstruction procedure needed in the finite volume context (N = 0, M > 0) in order to obtain order of accuracy M + 1 in space starting from the piecewise constant values of u n h (x, t n ) in P n i and its neighbors, i.e., in order to obtain a high order polynomial of degree M representing our solution in each P n where the ψ ℓ functions simply coincide with the ϕ ℓ basis functions of (5). Our reconstruction procedures are based on the WENO algorithm in its polynomial formulation as presented in Dumbser et al. [38], Dumbser and Käser [32,123], Titarev et al. [129], Tsoutsanis et al. [130], Levy et al. [131], Dumbser et al. [132], and Semplice et al. [133], and not based on the original version of WENO proposed in Jiang and Shu [121], Balsara and Shu [122], Hu and Shu [134], and Zhang and Shu [124] which provides only point values. For each P n i , the basic idea consists in (i) selecting a central stencil of elements S 0 i with a total number of elements, containing the cell P n i itself, its first layer of Voronoi neighbors V(P n i ) and filled by recursively adding neighbors of those elements that have been already included in the stencil, and in (ii) using the cell-average values of the elements of S 0 i to reconstruct a polynomial of degree M by imposing the integral conservation criterion, i.e., by requiring that its average on each cell match the known cell average. If f > 1 (which occurs in the unstructured case, where we take f = 1.5), this of course leads to an overdetermined linear system, which is solved using a constrained least-squares technique (CLSQ) [123], i.e., the reconstructed polynomial has exactly the cell averageû n 0,i on the polygon P n i and matches all the other cell averages of the remaining stencil elements in the least-square sense.
However, as well-known thanks to the Godunov theorem [1], the use of only one central stencil (which is indeed a linear procedure) would introduce oscillations in the presence of shock waves or other discontinuities. So, in order to make the reconstruction procedure non-linear, we will compute the final reconstruction polynomial as a non-linear combination or more than only one reconstruction polynomial, each one defined on a different reconstruction stencil S s i . We refer to the cited literature for further details, and here we just highlight the main characteristics of the two reconstruction procedures adopted in this work.

Case A: Cartesian Mesh
In Case A, of a fixed Cartesian mesh, we employ the polynomial WENO procedure given in Dumbser et al. [54], which is implemented in a dimension by dimension fashion. For each cell, we define its related sets of one-dimensional reconstruction stencils as where L = {M, s} and R = {M, s} denote the order and stencil dependent spatial extension of the stencil to the left and to the right. For odd order schemes we consider three stencils, one central, one fully left-sided, and one fully right-sided stencil in each space dimension (see Figure 4 for a graphical interpretation for M = 2), while for even order schemes we have four stencils, two of which are central, while the remaining two are again given by the fully left-sided and fully right-sided in each space dimension. In both cases the total amount of elements in each stencil is always n e = M + 1, the order of the scheme. Focusing on the reconstruction procedure on the x direction, given a element P n i , we start by expressing the first coordinate of the reconstruction polynomial at each stencil in terms of one dimensional basis functions, Then, we integrate on the stencil elements obtaining an algebraic system on the polynomial coefficients: withū n mk the average value obtained by integrating the solution at the previous time step on the cell P mk . Once the coefficients, and thus the polynomials, related to all the stencils are obtained, we compute a reconstruction polynomial in the x direction as the data-dependent non-linear combination of these, where n s is the number of stencils, n s = 3 if M =2 and n s = 4 otherwise; and ω s denote the non-linear weights (see Dumbser et al. [54] for further details).
To complete the reconstruction polynomial, we now repeat the above procedure in the y direction for each degree of freedom w n jk,ℓ 1 . First, we write the reconstruction polynomial in terms of the basis functions, Then, we solve the algebraic system 1 y m Finally, we get the WENO reconstruction polynomial In order to enforce bounds on the WENO reconstruction polynomial, such as the condition 0 ≤ α ≤ 1 on the volume fraction function of for example (56a), we rescale the reconstruction coefficientsŵ n jk,ℓ 1 ℓ 2 around the cell average as follows:ŵ * where the scaling factor ϕ jk is computed via the Barth and Jespersen limiter (see Barth and Jespersen [135]) applied to the volume fraction function α in all Gauss-Legendre and Gauss-Lobatto quadrature nodes, i.e., ϕ jk = min(ϕ jk,p ) is the global minimum in each element, with the nodal limiter values given by Here α max = 1 − ε ≤ 1 is the upper bound of the volume fraction function and α min = ε ≥ 0 is its lower bound;ᾱ denotes the cell average of α and α p denotes the node value of α in the quadrature point x p under consideration. As already mentioned above, this strategy is inspired from the Barth and Jespersen limiter [135], but also from the new bound-preserving polynomial approximation introduced in Després [136] and Campos-Pinto et al. [137]. Since the physical solution of α must satisfy 0 ≤ α ≤ 1, the above bound preserving limiter does not reduce the formal order of accuracy of the reconstruction, as proven in Després [136].

Case B: Moving Polygonal Mesh
In Case B of our moving and topology changing polygonal mesh we adopt a CWENO reconstruction algorithm, first introduced in Levy et al. [138][139][140] and Semplice et al. [133], and which can be cast in the general framework described in Cravero et al. [141]. We closely follow the work outlined in Dumbser et al. [132] and Boscheri et al. [142] for unstructured triangular and tetrahedral meshes, and extended it to moving polygonal grids in Gaburro et al. [77]. We emphasize that the main advantages of such a procedure is that only one stencil (the central one) is required to contain the total amount of elements stated in (10) and only this one is used to construct a polynomial of degree M; the other ones are used to compute polynomials of lower degree. In particular, we consider N n V i stencils S s i , each of them containing exactlŷ n e = (d + 1) cells, i.e., the central cell P n i and two consecutive neighbors belonging to V(P n i ). Refer to Figure 5 for a graphical description of the stencils. For each stencil S s i we compute a linear polynomial by solving a simple reconstruction system which is not overdetermined. According to the above mentioned literature, the reconstructed polynomial obtained via a nonlinear combination of the polynomial of degree M, computed over S s 0 , and of the N n V i linear polynomials, computed over S s i , maintains the order of convergence of the method and avoids unwanted spurious oscillations. In particular, in the case of moving meshes with topology changes, where the set of neighbors may change at any time step, the use of smaller so-called sectorial stencils significantly speeds up computations.
For the sake of uniform notation, in the DG case, i.e., when N > 0 and M = N, we trivially impose that the reconstruction polynomial is given by the DG polynomial, i.e., w n h (x, t n ) = u n h (x, t n ), which automatically implies that in the case N = M the reconstruction operator is simply the identity.

Space-Time Predictor Step
In this section we focus on the key feature, the element-local space-time predictor step, of our ADER FV-DG schemes: this part of the algorithm (the predictor) produces a high order approximation in both space and time of Q in all P n i . This allows to obtain a fully discrete one-step scheme that is uniformly high order accurate in both space and time.
The predictor step consists in a completely local procedure which solves the governing PDE (1) in the small, see Harten et al. [11], inside each space-time element C n i , and it only considers the geometry of volume C n i , the initial data w n h on P n i and the governing Equations (1), without taking into account any interaction between C n i and its neighbors. Because of this absence of communications, we refer to it as local. The procedure finally provides, for each C n i , a space-time polynomial data representation q n h , which serves as a predictor solution, only valid inside C n i , to be used for evaluating the numerical fluxes, the non-conservative products and the algebraic source terms when integrating the PDE in the final corrector step (see section 2.7) of the ADER scheme.
The predictor q n h is a polynomial of degree M, which takes the following form where θ ℓ (x, t) can be either i. For fixed and adaptive Cartesian grids (Case A), nodal space-time basis functions of degree M given by the product of one-dimensional nodal basis functions verifying (6) (with d = 1 ), two of them mapped to the unit interval [0, 1] as in (7) and with the time coordinate mapped to the reference time τ ∈ [0, 1] via t = t n + τ t. In this case, the total number of GL quadrature points per cell, as well as the total number of DOF is Q = (M + 1) d+1 , see also  In the other panels we report two of the N n Vi = 5 sectorial stencils containing the element itself and two consecutive neighbors belonging to V(P n i ). dimensions plus time) are used, which read with the total number of DOF Q = 1 Since we are only interested in an element local predictor solution, i.e., we do not need to consider the interactions with the neighbors, we do not yet take into account the jumps of q n h across the space-time lateral surfaces, because this will be done in the final corrector step (section 2.7). Instead, we insert the known discrete solution w n h (x, t n ) at time t n in order to introduce a weak initial condition for solving our PDE; note that w n h (x, t n ) uses information coming from the past only (following an upwinding approach) in such a way that the causality principle is correctly respected. To this purpose, the first term is integrated by parts in time. This leads to Equation (25) results in an element-local non-linear system for the unknown degrees of freedomq n ℓ of the space-time polynomials q n h . The solution of (25) can be found via a simple and fast converging fixed point iteration (a discrete Picard iteration) as detailed e.g., in Dumbser et al. [42] and Hidalgo and Dumbser [41]. For linear homogeneous systems, the discrete Picard iteration converges in a finite number of at most N + 1 steps, since the involved iteration matrix is nilpotent, see Jackson [143]. Moreover a proof of the convergence of this procedure in the case of a non-linear homogeneous conservation law in 1D is given in next section 2.6.

Simplification in the Case of a Fixed Cartesian Mesh
The space-time predictor step formerly presented can be simplified in the case of a Cartesian mesh with nodal basis functions resulting in a more efficient algorithm. Under these assumptions the governing PDE (1), can be rewritten as with Now, by substituting the discrete space-time predictor solution q n h with its expansion on the nodal basis and after integrating by parts in time, we obtain 1 0 1 0 1 0 θ k (ξ , η, 1) θ ℓ (ξ , η, 1)q n ℓ dξ dηdτ To recover the value of the unknown degrees of freedomq n ℓ , it is sufficient to solve the above equation locally for each element. One important advantage of using the nodal Gauss-Legendre basis is that the terms in (29) can be evaluated in a dimensionby-dimension fashion.

Space-Time Predictor for Sliver Space-Time Elements
When a topology change occurs, some space-time sliver elements, as those shown on the right side of Figure 8, are originated (see Gaburro et al. [77]), and the predictor procedure over them needs particular care. The problem connected with sliver elements is the fact that their bottom face, which consists only in a line segment, is degenerate, hence the spatial integral over P n i vanishes, i.e., there is no possibility to introduce an initial condition for the local Cauchy problem at time t n into their predictor. Thus, in order to couple however (24) with some known data from the past, we will end up with a formula different from (25). We underline that we first carry out the space-time predictor for all standard elements using, which can be computed independently of each other, and only subsequently we process the remaining space-time sliver elements. Then, when FIGURE 8 | Space time connectivity with topology changes and sliver element. Left: at time t n the polygons P n 2 and P n 3 are neighbors and share the highlighted edge, instead at time t n+1 they do not touch each other; the opposite situation occurs for polygons P n 1 and P n 4 . This change of topology causes the appearance of degenerate elements of different types (refer to Gaburro et al. [77] for all the details). In particular, so-called space-time sliver elements (right) need to be taken into account when considering the space-time framework, so the predictor and the corrector step have to be a adapted to their special features. Sliver elements (right) are indeed completely new control volumes which do neither exist at time t n , nor at time t n+1 , since they coincide with an edge of the tessellation and, as such, have zero areas in space. However, they have a non-negligible volume in space-time. The difficulties associated to this kind of element are due to the fact that w h is not clearly defined for it at time t n (thus the predictor has to be modified) and that contributions across it should not be lost at time t n+1 in order to guarantee conservation (thus the corrector has to be modified).
considering a sliver, we use the upwinding in time approach on the entire space-time surface ∂C n i that closes a sliver control volume, and again respecting the causality principle, we take the information to feed the predictor only from the past, i.e., only from those space-time neighbors C n j whose common surface ∂C n ij exhibit a negative time component of the outward pointing space-time normal vector (ñ t < 0). In this way, we can introduce information from the past into the space-time sliver elements.
As a consequence, the predictor solution q n h is again obtained by means of (24), but by treating the entire ∂C n i with the upwind in time approach, i.e., by considering also the jump terms between the still unknown predictor of the slivers (call it q n,− h ) and the already known predictors of its neighbors (call them q n,+ h ), where ∂C − i = ∂C n i withñ t < 0 is the part of the space-time boundary that has a negative time component of the space-time normal vector. Note that here we have taken into account also the jump of the non-conservative terms, and that these contributions have been added entirely [i.e., not only half of them, as in (49)]. Indeed, in (49) half of the jump contribution goes to one element, while the other half goes to the neighboring element; here instead, since the interaction between neighbors is only computed from the side of the sliver element, the entire jump contributes to the predictor in the sliver element.

Convergence Proof of the Predictor
Step for a Non-linear Conservation Law In this section, the convergence proof of the predictor for a non-linear conservation law is given. The proof is provided, for simplicity, in the case of a fixed mesh in one space dimension, following the nomenclature already employed in section 2.5.1, but it still holds in higher dimensions. Let us consider a general hyperbolic system of conservation laws of the form Then, the corresponding space-time DG predictor used in the ADER-DG framework reads For convenience, all derivatives and integrals in (32) have been transformed to the reference space-time element [0, 1] 2 . Moreover, the discrete solution is given by q h = θ l (ξ , τ )q ℓ , and the flux is expanded in the same basis as f h = θ ℓ (ξ , τ )f ℓ . When using a nodal basis, we can compute the degrees of freedom for the flux interpolant f h simply asf ℓ = f q ℓ . We also recall that the initial condition given by the DG scheme at time t n reads w h = ϕ ℓ (ξ )ŵ ℓ . Then, integration of the first term in (32) by parts in time yields and insertion of the definitions of the discrete solution leads to The iterative scheme employed to find the solution for the spacetime degrees of freedomq, at any Picard iteration r, can therefore be rewritten in compact matrix-vector notation as with where we have dropped the indices to ease the notation. After inverting K 1 (this matrix is built using the linearly independent basis functions so that it is invertible), we obtain the explicit iteration formulâ To prove that the former iterative formula will converge, we introduce the operator and the induced matrix norm Furthermore, we assume the flux to be Lipschitz continuous with Lipschitz constant L > 0 so that We now need to show that the operator ϕ is a contraction: The operator is therefore a contraction under the CFL-type condition on the time step t which connects the Lipschitz constant L with the mesh spacing x and the matrix norm of K −1 1 K ξ . Since the operator is contractive under the above assumptions, the Banach fixed point theorem, Banach [144], guarantees convergence of the iterative method.
In the previous reasoning, we have assumed that the inequality in the right hand side of (43) be strict. Thus, to conclude the proof, let us assume that the equality holds, this is true if and only if K −1 1 K ξ = 0. By taking into account the definition of the induced matrix norm (40), it implies K −1 1 K ξ x = 0 for any x in the metric space. Thus, K −1 1 K ξ = 0. Direct substitution in (38) gives so that no iterative procedure is done. Note: The matrix K −1 1 K ξ has been proven to be nilpotent and thus all its eigenvalues are zero, see Jackson [143], which guarantees convergence to the exact solution in a finite number of steps for linear homogeneous PDE.

Corrector Step
The corrector step is the last step of our path-conservative ADER FV-DG scheme, where the update of the solution from time t n up to time t n+1 can take place in a single step procedure thanks to the use of the predictor q n h . The update formula is recovered starting from the space-time divergence form of the PDẼ ∇ ·F(Q) +B(Q) ·∇Q = S(Q),F = (F, Q), which is multiplied by a set of space-time test functionsφ k and integrated over each space-time control volume C n Frontiers in Physics | www.frontiersin.org Note that the employed test functionsφ k coincide with the θ k of (22) for the Cartesian Case A. Instead, for the moving polygonal Case B, they need to be tied to the motion of the barycenter x b i (t) and must be moved together with P i (t) in such a way that at time t = t n they refer to the current barycenter x n b i and at time t = t n+1 they refer to the new barycenter x n+1 b i , thus they are defined as follows These moving modal basis functions are essential to the moving approach presented in Gaburro et al. [77] and used in this paper. They naturally allow for topology changes, without the need of any remapping steps, which we want to avoid in a direct ALE formulation. Now, (46) by applying the Gauss theorem to the fluxdivergence term and by splitting the non-conservative products into their volume and surface contribution, becomes where Q on P n+1 i is represented by the unknown u n+1 h , on P n i is taken to be the current representation of the conserved variables u n h , in the interior of C n i is given by the predictor q n h and on the space-time lateral surfaces ∂C n ij is given by q n,− h and q n,+ h which are the so-called boundary-extrapolated data, i.e., the values assumed respectively by the predictors of the two neighbor elements C n i and C n j on the shared space-time lateral surface ∂C n ij . Furthermore, we have employed a two-point path-conservative numerical flux function of Rusanov-type where s max is the maximum eigenvalue of the ALE Jacobian matrices A V n (q n,+ h ) and A V n (q n,− h ) being and the path = (q − h , q + h , s) is a straight-line segment path connecting q n,− h and q n,+ h which allow to treat the jump of the non-conservative products following the theory introduced in Dal Maso et al. [145], Parés [146], and Castro et al. [147], and extended to ADER FV-DG schemes of arbitrary high order in Dumbser et al. [46] and Dumbser and Toro [148]. Despite in this paper we only consider the Rusanov flux, the above methodology can be extended to different flux functions, adapting to the new flux splitting techniques like the ones presented in Toro and Vázquez-Cendón [149]. Finally, the time step size t is given by (52) where h min is the minimum characteristic mesh-size, ℓ i j is the length of the edge j of P n i and |λ max | is the spectral radius of the Jacobian of the flux F. Stability on unstructured meshes is guaranteed by the satisfaction of the inequality CFL < 1 d , see Dumbser et al. [42].
We close this section by remarking that the integration of the governing PDE over closed space-time volumes C n i automatically satisfies the geometric conservation law (GCL) for all test functionsφ k . This simply follows from the Gauss theorem and we refer to Boscheri and Dumbser [63] for a complete proof.

A Posteriori Subcell Finite Volume Limiter
Up to now, we have presented a family of FV and DG type schemes which achieves arbitrary high order of accuracy in space and time; the main difference between the FV and the DG approach lies in the fact that FV schemes, thanks to the WENOtype non-linear reconstruction procedure, are robust in the presence of shocks and discontinuities, while the DG formulation as presented so far, being linear in the sense of Godunov, is subject to the appearance of spurious oscillations. Thus, in order to employ a DG scheme in the context of solving hyperbolic partial differential equations, where usually discontinuities are developed, a technique that is able to limit spurious oscillations (called limiter) should be introduced. Several attempts in that direction can be found in the literature. For example, we could recall the artificial viscosity technique used in Hartmann and Houston [150],Persson and Peraire [151], and Cesenek et al. [152] which consists in adding a small parabolic term in the equation in order to smooth out the discontinuities.
Here, instead, we follow a different approach based on exploiting the respective strengths of FV and DG schemes, i.e., the resolution of DG in smooth regions and the robustness of FV across discontinuities. Thus, we first evolve the solution everywhere by using our DG scheme; then, we check a posteriori, at the end of each time step, if the obtained DG solution in each cell respects or not some criteria [as density and pressure positivity, a relaxed discrete maximum principle, specific physical bounds, or more elaborate choices as those of Guermond et al. [153]], and we mark as troubled those cells where the obtained DG solution is marked as not acceptable. Only for these troubled cells we repeat the time step using, instead of the DG scheme, a second order TVD FV method, which always assures a robust solution.
This idea is founded on works as those of Cockburn and Shu [154], Qiu and Shu [155,156], Balsara et al. [157], Luo et al. [158], Krivodonova [159], Zhu et al. [160], Zhu and Qiu [161], Clain et al. [86], Diot et al. [87,88], Loubére et al. [79], Boscheri et al. [162],and Boscheri and Loubére [83]; but in particular, here, we adopt a so-called subcell approach aimed at not losing the resolution of the DG scheme when switching to the FV method, as forwarded in Sonntag and Munz [163], Dumbser et al. [78], Zanotti et al. [80], Dumbser and Loubére [81], Boscheriand Dumbser [119], Fambri et al. [84], Rannabauer et al. [164], de la Rosa and Munz [165], and Boscheri et al. [142]. Indeed, at the beginning of the time step we project the DG solution u n h of a troubled cell P n i on a subdivision of it in sub-cells s n i,α obtaining a value for the cell averages on s n i,α at time t n v n i,α (x, t n ) = We evolve the cell averages up to time t n+1 using a classical TVD FV scheme, obtaining v n+1 i,α (x, t n+1 ). Finally, we recover a DG polynomial representation of the solution at time t n+1 over P n+1 i using the values on the sub-grid level v n+1 i,α and by applying a reconstruction operator as where the reconstruction is imposed to be conservative on the main cell P n+1 i yielding the additional linear constraint Thus, the limited solution on a troubled cell is robust thanks to the use of a TVD scheme and accurate thanks to the subcell resolution. For all the details of the a posteriori subcell FV limiter used in this work, we refer to Dumbser et al. [78] and Fambri et al. [36] for the fixed Cartesian Case A and to Gaburro et al. [77] for the moving polygonal Case B.

Governing PDE System
A simplified diffuse interface formulation of the unified continuum fluid and solid mechanics model [57,59,60,166], which can be used for modeling moving boundary problems of fluids and solids of arbitrary geometry, is given by the following PDE system (throughout this paper we make use of the Einstein summation convention over repeated indices) Here, (56a) is the evolution equation for the color function α that is needed in the diffuse interface approach as introduced in Tavelli et al. [85] for the description of linear elastic solids of arbitrary geometry and as used in Dumbser [106] and Gaburro et al. [107] for a simple diffuse interface method for the simulation of non-hydrostatic free surface flows. We assume that the color function α equals to 1 in the regions of the computational domain occupied by the material and 0 outside these regions. In the computational code, α = 1 − ε inside of the material and α = ε outside the material. Here, ε is a small parameter ε ≪ 1, see section 4. Then, inside of the diffuse interface, α may take any values between 0 and 1 (between ε and 1 − ε in the computational code). Equation (56b) is the mass conservation law and ρ is the material density; (56c) is the momentum conservation law, where v i is the velocity field and g i is the gravity vector; (56d) is the evolution equation for distortion field A ik (non-holonomic basis triad, see Peshkov et al. [167]); (56e) is the evolution equation for the specific thermal impulse J k constituting the heat conduction in the matter via a hyperbolic (non-Fourier-type) model. Finally, (56f) is the entropy balance equation and (56g) is the energy conservation law. Other thermodynamic parameters are defined via the total energy potential E = E(α, ρ, S, v, A, J): ik = pδ ik − σ ik is the total stress tensor (δ ik is the Kronecker delta); p = ρ 2 E ρ is the thermodynamic pressure; σ ik = −ρA jk E A ji is the non-isotropic part of the stress tensor, T = E S is the temperature, and the notations such as E ρ , E A ik , etc. stand for the partial derivatives of the energy potential, e.g., , etc. The dissipation in the medium includes two relaxation processes: the shear stress relaxation characterized by the scalar function θ 1 (τ 1 ) > 0 depending on the relaxation time τ 1 and thermal impulse relaxation characterized by θ 2 (τ 2 ) > 0 depending on the relaxation time τ 2 . Both these relaxation processes then contribute to the entropy production term [the source on the right hand-side of (56f)] which is positive because it is quadratic in E A ik and E J k .
From the mathematical standpoint, the unification of the model (56) consists in the use of only first-order hyperbolic equations for both dissipative and non-dissipative processes in contrast to the classical continuum mechanics relying on the mixed hyperbolic-parabolic formulations such as the famous Navier-Stokes-Fourier equations, for example. From the physical standpoint, the unification of Equations (56) consists in treating solid and fluid states of matter from the solid-dynamics viewpoint. Indeed, as discussed in Peshkov and Romenski [57] and Dumbser et al. [59,166], similarly to standard continuum solid-dynamics, the distortion field introduces additional degrees of freedom (in comparison to the classical continuum fluid mechanics) which characterizes deformation and rotational degrees of freedom of the continuum particles, represented not as scaleless mathematical points but characterized by a finite length scale, or equivalently, time scale τ 1 , e.g., see Dumbser et al. [166]. In such a formulation, solid-type behavior corresponds to relaxation times τ 1 such that T problem ≪ τ 1 , while the fluidtype behavior corresponds to τ 1 ≪ T problem , where T problem is the characteristic time scale of the problem under consideration.
In order to close system (56), that is, in order to define pressure p = ρ 2 E ρ , stresses σ ik = −ρA jk E A ji , temperature T = E S , and the dissipative source terms, one needs to provide the energy potential E. In this paper, we rely on a rather simple choice of E, which is, however, enough to deal with Newtonian fluids and simple hyperelastic solids. Thus, we assume that the specific total energy can be written as a sum of three contributions as (57) with the specific internal energy given by the ideal gas equation of state in the case of gases, and given by either the so-called stiffened gas equation of state or the well-known Mie-Grüneisen equation of state in the case of solids and liquids. Here, c v is the specific heat capacity at constant volume, γ is the ratio of the specific heats, p 0 is the reference (atmospheric) pressure, ρ 0 is the reference material density, and Ŵ 0 , and s are some material parameters. The specific energy stored in material deformations and in the thermal impulse is where • Gij= G ij − 1 3 G kk δ ij is the trace-free part of the metric tensor G ij = A ki A kj , which is induced by the mapping from Eulerian coordinates to the current stress-free reference configuration. The coefficientsc s (α) andc h (α) in (61) are the characteristic velocities for propagation of shear and thermal perturbations accordingly. In the present diffuse interface model, we choose the following simple linear mixture rule for the computation of the shear sound speed and of the heat wave propagation as a function of the volume fraction ᾱ where c s and c h are the material parameters inside the continuum and c g h ≪ 1 and c g s ≪ 1 are free parameters that can be chosen for the region outside the continuum. The specific kinetic energy is contained in the third contribution to the total energy and reads E 3 (v k ) = 1 2 v i v i . With the equation of state chosen above, we get the following expressions for the stress tensor, the heat flux and the dissipative sources E A ik and E J k present in the relaxation source terms: The functions θ 1 and θ 2 are chosen in such a way that a constant viscosity and heat conduction coefficient are obtained in the stiff relaxation limit, see Dumbser et al. [59] for a formal asymptotic analysis, Thus, following the procedure detailed in Dumbser et al. [59], one can show via formal asymptotic expansion that in the stiff relaxation limit τ 1 → 0, τ 2 → 0, the stress tensor and the heat flux reduce to and that is the effective shear viscosity and effective heat conductivity of model (56) are with ρ 0 and T 0 are reference density and temperature, see Dumbser et al. [59], where also an explanation has been provided of how the relaxation times τ could be obtained experimentally via ultrasound measurements.

Symmetric Godunov Form of the Model
It is important to note an interesting structural feature of Equations (56) that may affect future developments of the ADER schemes in an attempt to respect such structural properties at the discrete level that may help to improve physical consistency of the numerical solution. Thus, as many PDE systems studied in some other of our papers [59,60,168,169], system (56) belongs to the class of so-called Symmetric Hyperbolic Thermodynamically Compatible (SHTC) PDE systems originally studied by Godunov [170,171] and later by Godunov and Romenski [172], Godunov et al. [173], Romenski [168] and Romensky [174]. Indeed, by simply rescaling the quantitiesρ = αρ,p = αp =ρ 2 Eρ, andσ ik = ασ ik = −ρA jk E A ji and replacing the non-conservative Equation (56a) by an equivalent (on smooth solutions) conservative form (69a), system (56) can be written as where we have omitted the energy equation. Now, this system looks exactly as the system studied in Dumbser et al. [59], apart from the additional Equation (69a) which has the same structure as (69b) and does not change the essence. Then, after denoting E =ρE and introducing new variables P = (̺ 1 , ̺ 2 , v i , α ik , i , σ ) which are thermodynamically conjugate to the conservative variables Q = (αρ,ρ,ρv i , A ik ,ρJ i ,ρS), and a new thermodynamic potential L(P) = Q · E Q − E = Q · P − E, system (69) can be written in a symmetric form In this PDE system, the first two terms in each equation form the canonical Godunov form introduced in Godunov [170] which can be immediately written as a quasilinear symmetric form, e.g., see Peshkov et al. [169], Romenski [168], and Romensky [174]. The other (non-conservative) terms obviously form a symmetric matrix. Therefore, the entire system (71) can be written in a symmetric quasi-linear form and hence, it is a symmetric hyperbolic system if the thermodynamic potential L is convex. We note that the understanding of the structural properties of the continuous equations might be beneficial for developing of so-called structure-preserving numerical integrators (e.g., symplectic integrators). Thus, the energy conservation law (56g) is in fact a consequence of the other Equations (56) or (71), e.g., see Dumbser et al. [59] and Peshkov et al. [169], and can be viewed as a constraint of the system (71). Its non-violation at the discrete level cannot be guaranteed by the general purpose ADER family of schemes studied in this paper and hence, usually, as well as in our implementation, it is included into the set of discretized PDEs instead of the entropy equation. In principle, a structure-preserving scheme which satisfies all SHTC properties [169] of the continuous equations at the discrete level should guarantee the automatic satisfaction of the energy conservation law, without its explicit discretization. We hope to cover this topic in future work.

NUMERICAL RESULTS
In this section, we present some numerical results in order to illustrate the capabilities and potential applicability of the proposed numerical approach in non-linear continuum mechanics. The first three test problems are carried out without making explicit use of the diffuse interface approach, i.e., setting α = 1 everywhere in the entire computational domain. The last three test problems illustrate the full potential of the diffuse interface extension of the GPR model in the context of moving free boundary problems. Gravity effects are neglected in all test cases, apart from the dambreak problem shown in subsection 4.6. Whenever values for ν = µ/ρ 0 and c s are provided, the corresponding relaxation time τ 1 is computed according to (68).

Numerical Convergence Studies in the Stiff Relaxation Limit
In order to verify the high order property of our ADER schemes in both space and time in the stiff relaxation limit, we first represent the numerical convergence study that was already carried out in Dumbser et al. [59] on a smooth unsteady flow, for which an exact analytical solution is known for the compressible Euler equations, i.e., in the stiff relaxation limit τ 1 → 0 and τ 2 → 0 of the GPR model. The problem setup is the one of the classical isentropic vortex, see Hu and Shu [175]. The initial condition consists in a stationary isentropic vortex, whose exact solution can easily be found by solving the compressible Euler equations in cylindrical coordinates. Due to the Galilean invariance of the Euler equations and of the GPR model, one can then simply superimpose a constant velocity field to this stationary vortex solution in order to get an unsteady version of the test problem. The vortex strength is chosen as ε = 5 and the perturbation of entropy S = p ρ γ is assumed to be zero. For details of the setup, see Hu and Shu [175] and Dumbser et al. [59]. In this test we set the distortion field initially to  Table 1, together with the chosen values for the effective viscosity µ and the effective heat conductivity coefficient κ. From Table 1 one can observe that high order of convergence of the numerical method is achieved also in the stiff limit of the governing PDE system.

Circular Explosion Problem in a Solid
In this Section, we simulate a circular explosion problem in an ideal elastic solid. We compare the results obtained with a third order ADER-WENO finite volume scheme on moving unstructured Voronoi meshes with possible topology changes, Gaburro et al. [77], with those obtained with a fourth order ADER discontinuous Galerkin finite element scheme on a very fine uniform Cartesian mesh composed of 512 × 512 elements, which will be taken as the reference solution for this benchmark.  1 | Experimental errors and order of accuracy at time t = 1 for the density ρ for ADER-DG schemes applied to the GPR model (c s = 0.5, α = 1) in the stiff relaxation limit (µ ≪ 1, κ ≪ 1). mesh with 82 919 control volumes. The computational results obtained with the unstructured ADER-WENO ALE scheme and those obtained with the high order Eulerian ADER-DG scheme are presented and compared with each other in Figure 9.
We can note a very good agreement between the two results. The high quality of the ADER-WENO finite volume scheme on coarse grids is mainly due to the natural mesh refinement around the shock, which is typical for Lagrangian schemes. Furthermore, Lagrangian schemes are well-known to capture material interfaces and contact discontinuities very well, since the mesh is moving with the fluid and thus numerical dissipation at linear degenerate fields moving with the fluid velocity is significantly lower than with classical Eulerian schemes.

Rotor Test Problem
A second solid mechanics benchmark consists in the simulation of a plate on which a rotational impulse is initially impressed, in a circular region centered with respect to the computational domain. This rotor will initially move according to the rotational impulse, while emitting elastic waves which ultimately determine the formation of a set of concentric rings with alternating direction of rotation. The test is analogous to the rotor problem shown in Peshkov et al. [72], but with a weakened material in order to show stronger motion of the Voronoi grid.
The results of the third order ADER-WENO finite volume method on a moving Voronoi grid with variable connectivity, composed of 150 561 cells, are compared against a reference  solution obtained with a fourth order ADER discontinuos Galerkin scheme on a very fine uniform Cartesian mesh counting 512 × 512 elements, for a total of over four million spatial degrees of freedom.
The computational domain is the square = [−1, 1] × [−1, 1] and the final simulation time is set to t = 0.5. With exception made for the velocity field, all variables are initially constant throughout the domain. Specifically we set α = 1, ρ = 1, p = 1, A = I, J = 0, while the velocity field is v = [−y/R, x/R, 0] if r = x 2 + y 2 ≤ R, and v = 0 otherwise, that is, outside of the circle of radius R = 0.2; this way, the initial tangential velocity at r = R is one. The solid is taken to be elastic (τ 1 → ∞), heat wave propagation is neglected (c h = 0), and the characteristic speed of shear waves is c s = 0.25. The constitutive law is chosen to be the stiffened-gas EOS with γ = 1.4 and p 0 = 0. We can see in Figure 10 that, although some of the finer features are lost (specifically the small central counterclockwise-rotating ring) due to the lower resolution of the finite volume method on a coarser grid, the shear waves travel outwards with the correct velocity and the moving Voronoi finite volume simulation can be said to be in agreement with the high resolution discontinuous Galerkin results. Also in Figure 10, it is shown that the central region of the computational grid has undergone significant motion but thanks to the absence of constraints on the connectivity between elements, the Voronoi control volumes have not been stretched excessively as would instead happen for a similar moving unstructured grid, but with fixed connectivity.

Elastic Vibrations of a Beryllium Plate
The first benchmark for our new diffuse interface version of the GPR model consists in the purely elastic vibrations of a beryllium plate, subject to an initial velocity distribution, see for example Sambasivan et al. [176], Maire et al. [177], Burton et al. [178], Boscheri et al. [71], and Peshkov et al. [72] for a setup of the same test problem in the framework of Lagrangian and ALE schemes.
Unlike in Lagrangian schemes, no boundary conditions need to be imposed on the surface of the bar. We simply use transmissive boundaries on ∂ . The entire computational domain is initialized with the reference density for beryllium as ρ(x, 0) = ρ 0 , while the pressure is set to p(x, 0) = 0. The distortion field is initialized with A = I. According to Burton et al. [178], the final time is set to t f = 53.25 so that it corresponds approximately to two complete flexural periods. The simulations are carried out with a third order ADER-WENO scheme on two uniform Cartesian meshes composed of 256 × 128 and 512 × 256 elements, respectively.
For the fine grid simulation in Figure 11, we present the temporal evolution of the color contour map of the volume fraction function α, which represents the moving geometry of the bar. Here, dark gray color is used to indicate the regions with α > 0.5 and white color is used for the regions of α < 0.5. In the same figure, we also depict the pressure field in the region α > 0.5 at times t = 5, t = 14, t = 23, and t = 28. These time instants cover approximately one flexural period. The time evolution of the vertical velocity component v(0, 0, t) in the origin is depicted in Figure 12. For comparison, in the same figure we also show the results obtained on the coarse mesh for the same test problem with a fourth order ADER-DG scheme with second order TVD subcell finite volume limiter (red line).
Our computational results compare visually well against available reference solutions in the literature, see Sambasivan et al. [176], Maire et al. [177], Burton et al. [178], Boscheri et al. [71], and Peshkov et al. [72], which were all carried out with pure Lagrangian or Arbitrary-Lagrangian-Eulerian schemes on moving meshes, while here we use a diffuse interface approach on a fixed Cartesian grid.
The Taylor bar impact problem is a classical benchmark for an elasto-plastic aluminium projectile that hits a rigid solid wall, see Sambasivan et al. [176], Maire et al. [177], Dobrev et al. [181], and Boscheri et al. [71]. The projectile is initially moving with velocity v = (−0.015, 0) toward a wall located at x = 0. This velocity field is imposed within the subregion b , while in the rest of the domain we set v = 0. The remaining initial conditions are chosen as ρ = ρ 0 , p = p 0 , A = I, J = 0 and with the parameters τ 0 = 1 and m = 20 for the computation of the relaxation time (73). Unlike in Lagrangian schemes, we do not need to set any boundary conditions on the free surface of the moving bar. We only apply reflective slip wall boundary conditions on the wall in x = 0. According to Maire et al. [177], Dobrev et al. [181], and Boscheri et al. [71] the final time of the simulation is t = 5, 000. The computational domain is discretized on a regular Cartesian grid composed of 512 × 256 elements using a third order ADER-WENO finite volume scheme. As in Boscheri et al. [71] we FIGURE 13 | Geometry of the Taylor bar at time t = 1, 000 (top) and at the final time t = 5, 000 (bottom) obtained with a third order ADER-WENO finite volume scheme applied to the diffuse interface GPR model. We plot the contour colors of the volume fraction function α, where black regions denote α > 0.5 and white regions α < 0.5. employ a classical source splitting for the treatment of the stiff sources that arise in the regions of plastic deformations, i.e., when σ ≫ σ 0 . In Figure 13, we show the computational results at t = 1000 and at the final time t = 5, 000. The obtained solution is in agreement with the results presented in Maire et al. [177], Boscheri et al. [71], and Peshkov et al. [72]. At time t = 5, 000, we measure a final length of the projectile of L f = 456, which fits the results achieved in Maire et al. [177] and Boscheri et al. [71] up to 2%.

Dambreak Problem
In this last section on numerical test problems, we solve a twodimensional dambreak problem with different relaxation times in order to show the entire range of potential applications of the GPR model. For this purpose, we also activate the gravity source term, setting the gravity vector to g = (0, −g) with g = 9.81. The computational domain is chosen as = [0, 4] × [0, 2] and is discretized with a fourth order ADER discontinuous Galerkin finite element scheme with polynomial approximation degree N = 3 and a posteriori subcell TVD finite volume limiter. Computations are run on a uniform Cartesian mesh composed of 128 × 64 elements until the final time t = 0.5. The initial condition is chosen as follows: we set ρ = ρ 0 , v = 0, A = I and J = 0 in the entire computational domain. We impose the slip boundary condition on the bottom. In the subdomain d = [0, 2] × [0, 1], we set α = 1 − ε, and p = ρ 0 g(y − 1), while in the rest of the domain we set α = ε and p = 0. In this test problem we set ε = 10 −2 and use a stiffened gas FIGURE 15 | Dambreak problem at t = 0.4, simulated with a fourth order ADER-DG scheme using a space-time adaptive Cartesian AMR mesh applied to the GPR model with with ν = 10 −3 (Top), and reference solution, computed with a third order ADER-WENO finite volume scheme on a very fine uniform Cartesian grid, solving the inviscid and barotropic reduced Baer-Nunziato approach presented in [106,107] (Bottom). equation of state with parameters ρ 0 = 1, 000, p 0 = 5 × 10 4 , γ = 2, c h = 0 and a shear sound speed c s = 6. Simulations are run in three different regimes, only characterized by a different choice of the strain relaxation time τ 1 . In the first simulation, we set τ 1 so that a kinematic viscosity ν = µ/ρ 0 = 10 −3 is reached in the stiff relaxation limit, i.e., the GPR model in this case describes an almost inviscid fluid. In the second simulation we choose τ 1 so that ν = 0.1, i.e., a high viscosity Newtonian fluid behavior is reached. In the last simulation we set τ 1 → ∞, i.e., the strain relaxation term is switched off so that an ideal elastic solid with low shear resistance is described, similar to a jelly-type medium. In all cases, we apply solid slip wall boundary conditions on the left and on the right of the computational domain, while on the right and upper boundary, transmissive boundary conditions are set. The temporal evolution of the volume fraction function α, together with the coarse mesh used in this simulation, are depicted in Figure 14. The results for the almost inviscid fluid agree qualitatively well with those shown in Ferrari et al. [182], Dumbser et al. [106], and Gaburro et al. [107] for nonhydrostatic dambreak problems. In order to corroborate this statement quantitatively, we now repeat the simulation with ν = 10 −3 using a fourth order ADER-DG scheme on a coarse AMR grid composed of only 32 × 16 elements on the level zero grid. We then apply two levels of AMR refinement with refinement factor r = 3, i.e., we employ a general space-tree, rather than a simple quad-tree. We note that the simulations on the AMR grid are run in combination with time-accurate local time stepping (LTS), which is trivial to implement in high order ADER-DG and ADER-FV schemes, due to their fully-discrete one-step nature. For details on LTS, see Dumbser et al. [37,54],Dumbser [64] and Gaburro et al. [65]. As a reference solution of this almost inviscid flow problem, we solve the reduced barotropic and inviscid Baer-Nunziato model introduced in Dumbser [106] and Gaburro et al. [107], using a third order ADER-WENO finite volume scheme on a very fine uniform Cartesian grid composed of 1024 × 512 elements. The direct comparison of the two simulations at time t = 0.4 is shown in Figure 15. Overall we can indeed note an excellent agreement between the behavior of the diffuse interface GPR model in the stiff relaxation limit and the weakly compressible inviscid non-hydrostatic free surface flow model of Dumbser [106] and Gaburro et al. [107].

CONCLUSIONS AND OUTLOOK
In the first part of this paper we have provided a review of the ADER approach, whose development started about 20 years ago with the seminal works of Toro et al. [20] Millington et al. [19], Titarev and Toro [29], and Toro and Titarev [28] in the context of approximate solvers for the generalized Riemann problem (GPR). The ADER method provides fully discrete explicit onestep schemes that are in principle arbitrary high order accurate in both space and time. The most recent developments include ADER schemes for stiff source terms, as well as ADER finite volume and discontinuous Galerkin finite element schemes on fixed and moving meshes, which are all based on a spacetime predictor-corrector approach. The fact that ADER schemes are fully discrete makes the implementation of time accurate local time stepping (LTS) particularly simple, both on adaptive Cartesian AMR meshes [54], as well as in the context of Lagrangian schemes on moving grids [64,65]. The fully discrete space-time formulation also allows the treatment of topology changes during one time step in a very natural way [77]. In the second part of the paper we have then shown several applications of high order ADER finite volume and discontinuous Galerkin finite element schemes to the novel unified hyperbolic model of continuum mechanics (GPR model) proposed by Godunov, Peshkov and Romenski [56,57,59]. The presented test problems cover the entire range of continuum mechanics, from ideal elastic solids over plastic solids to viscous fluids. The use of a diffuse interface approach allows also to simulate moving boundary problems on fixed Cartesian meshes. Future developments will concern the extension of the mathematical model to non-Newtonian fluids [183] and to free surface flows with surface tension, see Schmidmayer et al. [184] and Chiocchetti et al. [185], as well as to the conservative multi-phase model of Romenski et al. [186,187]. In future work we will also consider the use of novel all speed schemes [188] and semi-implicit space-time discontinuous Galerkin finite element schemes [189][190][191] for the diffuse interface version of the GPR model used in this paper.

DATA AVAILABILITY STATEMENT
The datasets generated for this study are available on request to the corresponding author.

AUTHOR CONTRIBUTIONS
The governing PDE system was developed by IP. The numerical method and the computer codes were developed by MD, EG, and SC. The test problems were computed by MD, EG, and SC. The analysis of the method was performed by SB. All authors discussed the results and contributed to the final manuscript.