# Very high order treatment of embedded curved boundaries in compressible flows: ADER discontinuous Galerkin with a space-time Reconstruction for Off-site data

Mirco Ciallella, Stephane Clain, Elena Gaburro, Mario Ricchiuto

► To cite this version:

Mirco Ciallella, Stephane Clain, Elena Gaburro, Mario Ricchiuto. Very high order treatment of embedded curved boundaries in compressible flows: ADER discontinuous Galerkin with a space-time Reconstruction for Off-site data. Computers & Mathematics with Applications, 2024, 175, pp.1-18. 10.1016/j.camwa.2024.08.028. hal-04706736

HAL Id: hal-04706736

https://inria.hal.science/hal-04706736v1

Submitted on 23 Sep 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Very high order treatment of embedded curved boundaries in compressible flows: ADER discontinuous Galerkin with a space-time Reconstruction for Off-site data

Mirco Ciallella<sup>a,*</sup>, Stephane Clain<sup>b</sup>, Elena Gaburro<sup>c</sup>, Mario Ricchiuto<sup>c</sup>

<sup>a</sup> École Nationale Supérieure d'Arts et Métiers, Institut de Mécanique et d'Ingénierie, 33400 Talence, France

<sup>b</sup> CMUC - Mathematical centre, Coimbra University, Coimbra 3001-501, Portugal

<sup>c</sup> Inria, Univ. Bordeaux, CNRS, Bordeaux INP, IMB, UMR 5251, 33405 Talence cedex, France

## Abstract

In this paper we present a novel approach for the design of high order general boundary conditions when approximating solutions of the Euler equations on domains with curved boundaries, using meshes which may not be boundary conformal. When dealing with curved boundaries and/or unfitted discretizations, the consistency of boundary conditions is a well-known challenge, especially in the context of high order schemes. In order to tackle such consistency problems, the so-called Reconstruction for Off-site Data (ROD) method has been recently introduced in the finite volume framework: it is based on performing a boundary polynomial reconstruction that embeds the considered boundary treatment thanks to the implementation of a constrained minimization problem. This work is devoted to the development of the ROD approach in the context of discontinuous finite elements. We use the genuine space-time nature of the local ADER predictors to reformulate the ROD as a single space-time reconstruction procedure. This allows us to avoid a new reconstruction (linear system inversion) at each sub-time node and retrieve a single space-time polynomial that embeds the considered boundary conditions for the entire space-time element. Several numerical experiments are presented proving the consistency of the new approach for all kinds of boundary conditions. Computations involving the interaction of shocks with embedded curved boundaries are made possible through an *a posteriori* limiting technique.

**Keywords:** Curved boundaries, Space-time schemes, Unfitted discretization, Discontinuous Galerkin, compressible flows

## 1. Introduction

The consistency of boundary conditions poses a significant challenge when dealing with high order discretization techniques. High-order methods have demonstrated their potential for providing more precise outcomes with fewer degrees of freedom [1]. However, the management of boundaries remains an open challenging issue. In fact, when curved boundaries are taken into consideration, the precise representation of geometries and the accuracy of simulations become critical concerns. In particular, approximations of boundary conditions must be considered to achieve the desired convergence properties.

Currently, the most traditional method to address these issues employs high order polynomial reconstructions [2] of the same degree used in the mathematical model discretization to estimate

<footnote>*Corresponding author (mirco.ciallella@ensam.eu)</footnote>
the targeted geometry. Additionally, other common approaches include local approximation of the curved geometry or the so-called iso-geometric analysis [3]. Curvilinear mesh generation presents a variety of challenges [4, 5, 6, 7], including nonlinear mappings associated with these elements and the use of complex quadrature formulas. Methodological complexity, along with these features, substantially contributes to increase the computational costs in particular because of the augmented number of required degrees of freedom.

Despite the advancements made in developing numerical methods to handle curvilinear elements [8, 9, 10, 11], obtaining a high-quality curved mesh remains a multifaceted task, particularly in dealing with practical geometric shapes. Alternatively, improving boundary conditions on a straight-faced mesh can be much cheaper, but it would need special care by taking into account the specific characteristics of the geometry in that area. We refer to [12, 13] for the first works on curvature corrected wall boundary conditions and the application to high order schemes. In the previous reference, the boundary normal used when prescribing the slip-wall condition is correctly taken into account, resulting in a recovery of third and, in some cases, fourth order of accuracy on 2D geometries. However, it is well-known that this method can only be formulated for slip-wall conditions, and the research is restricted to 2D geometries.

More recent works on this research topic can be found in [14, 15]. In [14] some of the authors introduced a simple approach to design boundary conditions by reformulating the Shifted Boundary Method (SBM) [16, 17, 18, 19, 20, 21, 22, 23, 24] as a polynomial correction. Although its simplicity to deal with Dirichlet conditions, the development of general Robin boundary conditions has only been addressed recently [25]. With the goal of considering general boundary conditions, some of the authors developed the Reconstruction for Off-site Data (ROD) approach in a high order finite volume framework [15]. The ROD approach replaces the polynomial reconstruction in each boundary cell with a constrained least squares problem. The minimization problem extracts information from neighboring cells and imposes constraints tied to the physical boundary. The work has been extended to other types of boundary conditions due to its generality [26, 27, 28, 29, 30, 31]. In the finite volume framework, the constrained optimization problem is constructed by employing, as constraints, the information provided by the cell averages of the neighboring elements.

In this paper, we investigate the development of the ROD approach within a discontinuous finite element framework [32, 33, 34, 35, 36]. Specifically, by having the local formulation of discontinuous finite element discretizations, a distinct constrained minimization problem can be formulated which only employs the information of the element under consideration. This simplifies certain aspects of the finite volume approach while maintaining its consistency. To achieve high order properties for time-dependent problems, the approach based on spatial polynomial reconstruction requires the solution of this minimization problem at every sub-time step (when coupled with Runge-Kutta schemes) of the considered time integration scheme, resulting in increased computational costs. To further simplify this process, a novel approach is introduced in the space-time ADER framework [37, 38, 39, 40, 41, 42, 43, 44]. In particular, we propose a one-step direct and cheap procedure by employing the space-time basis functions used in the ADER predictor also to write the ROD constrained least squares minimization problem. The obtained system requires only one inversion and provides a single polynomial for the entire space-time element. Due to the modal space-time representation used, the number of modes reconstructed is much reduced compared to the case of several independent reconstructions at each quadrature point in time. This provides a cost effective formulation of ROD for evolutionary problems. Numerical results confirm the accuracy of the proposed approach on unfitted meshes.

The paper is organized as follows. In section 2, we present the studied mathematical model, specifically the hyperbolic system of the Euler equations of gasdynamics. However, it should be
noted that the proposed boundary treatment can be suitably applied to other PDEs. In section 3, we review the high order discretization of conservation laws system using discontinuous finite element [33] with modal polynomials and the one-step ADER method with space-time predictor [37]. A shock limiting approach [45, 46, 47, 48] is also discussed in this section to deal with discontinuous features that may arise in the flow field. Additionally, in section 5, we introduce the boundary treatment based on the Reconstruction for Off-site Data (ROD) technique. We begin by examining the classical approach adapted to the discontinuous Galerkin discretization, which has not yet been explored in the literature. Then, we introduce a new method that utilizes the space-time high order predictor to create a boundary treatment that eliminates the need for polynomial reconstruction at each sub-time node of the ADER correction step. In section 6, we present numerical experiments demonstrating the convergence properties of the methods on conformal linear meshes and unfitted triangulations. The paper ends with section 7, which offers an outlook regarding the future analysis and development of these schemes.

## 2. Mathematical model

The mathematical model considered in this work is the system of Euler equations for compressible gasdynamics in $d = 2$ space dimensions reading:

$$ \frac{\partial \mathcal{U}}{\partial t} + \nabla \cdot \mathcal{F}(\mathcal{U}) = 0, \quad \text{on} \quad \Omega_T = \Omega \times [0, T] \subset \mathbb{R}^d \times \mathbb{R}^+, \qquad (1) $$

with $\mathcal{U} : \mathbb{R}^+ \times \mathbb{R}^d \to \mathbb{R}^{d+2}$ the vector of conserved variables and $\mathcal{F} : \mathbb{R}^{d+2} \to \mathbb{R}^d \times \mathbb{R}^{d+2}$ the nonlinear flux, respectively defined as

$$ \mathcal{U} = \begin{pmatrix} \rho \\ \rho \mathbf{u} \\ \rho E \end{pmatrix}, \quad \mathcal{F}(\mathcal{U}) = \begin{pmatrix} \rho \mathbf{u} \\ \rho \mathbf{u} \otimes \mathbf{u} + p \mathbf{I} \\ \rho H \mathbf{u} \end{pmatrix}, \qquad (2) $$

having denoted by $\rho$ the mass density, by $\mathbf{u}$ the velocity, by $p$ the pressure, and with $E = e + \mathbf{u} \cdot \mathbf{u}/2$ the specific total energy, $e$ being the specific internal energy. Finally, the total specific enthalpy is $H = h + \mathbf{u} \cdot \mathbf{u}/2$ with $h = e + p/\rho$ the specific enthalpy. The relation between the pressure and the internal energy is given by the perfect gas equation of state

$$ p = (\gamma - 1)\rho e, \qquad (3) $$

where $\gamma$ the constant ratio of specific heats ($\gamma = 1.4$ for air).

The validation of the numerical methods are performed by means of the manufactured solution method which needs the discretization of an additional source term in eq. (1) which thus reads as follows

$$ \frac{\partial \mathcal{U}}{\partial t} + \nabla \cdot \mathcal{F}(\mathcal{U}) = \mathcal{S}(\mathcal{U}). \qquad (4) $$

For this reason, we are going to include $\mathcal{S}$ in the space-time discretization presented in section 3.

## 3. High order space-time discretization

### 3.1. Computational domain and data representation

We discretize the domain $\Omega$ as a tessellation $\mathcal{T}$ composed of $\mathcal{N}$ non-overlapping simplicial elements (triangles in 2D). We denote by $K_j$ the generic element, so that the computational
domain is $\Omega_h = \bigcup_{j=1}^N K_j$. Note that in general $\Omega_h \neq \Omega$ and in particular $\partial\Omega_h \neq \partial\Omega$ for most approximations, even conformal, except for very simple geometries (with no curvature) or of isogeometric approaches [3].

The solution $\mathcal{U}$ is approximated by $\mathcal{U}_h$, which belongs to a space of piece-wise polynomials within each triangle $K_j$ and discontinuous across faces, such that in each element we have

$$ \mathcal{U}_h(\mathbf{x}, t)|_{K_j} = \sum_k \psi_k(\mathbf{x}) \hat{\mathcal{U}}_k(t), \qquad (5) $$

where $\psi_k$ is a basis of polynomials of degree $M$. It is well-known that the discontinuous finite element data representation (5) leads to a Finite Volume (FV) scheme if $M = 0$.

## 3.2. Discontinuous Galerkin in space

The elemental semi-discrete discontinuous Galerkin (DG) weak formulation is classically written by projecting each component of eq. (4) on the relevant basis and integrating by parts [33, 49]

$$ \int_{K_j} \psi \frac{d\mathcal{U}_h}{dt} d\mathbf{x} + \int_{\partial K_j} \psi \hat{F}(\mathcal{U}_h^{-}, \mathcal{U}_h^{+}) \cdot \mathbf{n} dS - \int_{K_j} \nabla \psi \cdot \mathcal{F}(\mathcal{U}_h) d\mathbf{x} = \int_{K_j} \psi \mathcal{S}(\mathcal{U}_h) d\mathbf{x}, \quad (6) $$

with $\hat{F}(\mathcal{U}_h^{-}, \mathcal{U}_h^{+})$ a consistent numerical flux which depends on the face values of the internal state $\mathcal{U}_h^{-}$, on the neighboring element state $\mathcal{U}_h^{+}$ and on the face normal $\mathbf{n}$. We recall in particular that a consistent flux is a Lipschitz continuous function whose arguments also verify

$$ \hat{\mathcal{F}}(U, U) \cdot \mathbf{n} = \mathcal{F}(U) \cdot \mathbf{n}. \qquad (7) $$

In this paper, we have used a classical Rusanov-type (local Lax-Friedrichs) flux:

$$ \hat{\mathcal{F}}(\mathcal{U}_h^{-}, \mathcal{U}_h^{+}) \cdot \mathbf{n} = \frac{1}{2} (\mathcal{F}(\mathcal{U}_h^{+}) + \mathcal{F}(\mathcal{U}_h^{-})) \cdot \mathbf{n} - \frac{1}{2} s_{max} (\mathcal{U}_h^{+} - \mathcal{U}_h^{-}), \quad (8) $$

where $s_{max}$ is the largest of the spectral radiuses of the flux Jacobians $\frac{\partial \mathcal{F}}{\partial \mathcal{U}}(\mathcal{U}_h^{+})$ and $\frac{\partial \mathcal{F}}{\partial \mathcal{U}}(\mathcal{U}_h^{-})$.

## 3.3. ADER in time with local space-time predictor

For the time integration, we have used a high order explicit predictor-corrector ADER method [37]. The time domain is discretized in temporal intervals $[t^n, t^{n+1}]$, where $t^{n+1} = t^n + \Delta t$. With the same notation used before, the ADER method can be recast as

$$ \mathcal{U}^{n+1} = \mathcal{U}^n - M^{-1} \int_{t^n}^{t^{n+1}} R(q_h) dt, \qquad (9) $$

where $M$ is the elemental mass matrix, $R$ includes all integrals in eq. (6) but the first, and $q_h = q_h(\mathbf{x}, t)$ is a high order space-time predictor of the solution valid in the time interval $[t^n, t^{n+1}]$. In practice, eq. (9) is replaced by a consistent high order quadrature formula in time (with $\alpha_i$ and $\omega_i$, the quadrature points and weights)

$$ \mathcal{U}^{n+1} = \mathcal{U}^n - M^{-1} \Delta t \sum_i \omega_i R(q_h(t^{n+\alpha_i})). \qquad (10) $$

$q_h(\mathbf{x}, t)$ is defined as a polynomial of degree $M$ in space and time (where time is considered as an additional physical dimension), and $R$ also includes consistent quadrature formulas for the integrals
in space. This polynomial is obtained by means of a genuinely local space-time procedure. In the current implementation this local problem is formulated by means of a modal expansion in the space-time element

$$ q_h(\mathbf{x}, t)|_{K_j \times [t^n, t^{n+1}]} = \sum_{\ell=0}^{Q-1} \theta_{\ell}(\mathbf{x}, t) \hat{q}_{\ell}, \quad Q = \mathcal{L}(M, d+1), \qquad (11) $$

where $\theta_{\ell}(\mathbf{x}, t)$ are $\mathcal{L}(M, d) = \frac{1}{d!} \prod_{m=1}^{d} (M+m)$ modal bases defined as

$$ \theta_{\ell}(\mathbf{x}, t)|_{K_j} = \frac{(x_1 - x_{1,j})^{p_{\ell}} (x_2 - x_{2,j})^{q_{\ell}} (t - t^n)^{r_{\ell}}}{p_{\ell}! h_j^{p_{\ell}} q_{\ell}! h_j^{q_{\ell}} r_{\ell}! h_j^{r_{\ell}}} \quad \ell = 0, \dots, \mathcal{L}(M, d+1), \quad 0 \le p_{\ell} + q_{\ell} + r_{\ell} \le M, \qquad (12) $$

with $(x_{1,j}, x_{2,j})$ the barycenter of element $K_j$, and $h_j$ is the diameter of its incircle.

The values $\hat{q}_{\ell}$ are obtained as the solution of the local space-time weak formulation

$$ \int_{K_j} \theta_k(\mathbf{x}, t^{n+1}) q_h(\mathbf{x}, t^{n+1}) d\mathbf{x} - \int_{K_j} \theta_k(\mathbf{x}, t^n) U_h^n(\mathbf{x}, t^n) dx - \int_{t^n}^{t^{n+1}} \int_{K_j} \frac{\partial \theta_k}{\partial t}(\mathbf{x}, t) q_h(\mathbf{x}, t) d\mathbf{x} dt \\ + \int_{t^n}^{t^{n+1}} \int_{K_j} \theta_k(\mathbf{x}, t) \nabla \cdot \mathcal{F}(q_h^n) d\mathbf{x} dt = \int_{t^n}^{t^{n+1}} \int_{K_j} \theta_k(\mathbf{x}, t) \mathcal{S}(q_h^n) d\mathbf{x} dt, \qquad (13) $$

where $U_h^n$ is the known initial condition at time $t^n$.

Equation (13) is local, meaning that no communication is exchanged between $K_j$ and its neighboring control volumes. The solution to this equation can be obtained independently within each element $K_j$ by means of some iterative procedure. Herein, a simple discrete Picard iteration for each space-time element is performed [41].

### 3.4. A posteriori sub-cell finite volume limiter

In order to properly handle discontinuous solutions we employ the MOOD approach [45, 50, 46], which has already been effectively applied in the ADER framework [47, 40, 39, 48].

The algorithm is based on an *a posteriori* technique. The solution is first evolved from $t^n$ to $t^{n+1}$ using the high order ADER-DG method. Then several admissibility criteria are checked and the solution in all troubled cells (i.e., the cells not satisfying such criteria) is recomputed *a posteriori* using a robust second-order MUSCL-Hancock TVD finite volume scheme, on a sub-triangulation of the initial grid. This allows us to preserve the accuracy of the high order DG scheme also when passing to a lower order but more robust scheme. All aspects of the implementation of this technique are provided in [48] to which we refer for details.

**Remark 1 (CFL constraint).** Concerning the choice of the time step, we have implemented the classical explicit CFL condition

$$ \Delta t < \text{CFL} \left( \frac{|h_{\text{min}}|}{(2M+1)|\lambda_{\text{max}}|} \right), \qquad (14) $$

where $|h_{\text{min}}|$ is the minimum characteristic mesh-size and $|\lambda_{\text{max}}|$ is the maximum eigenvalue of the Jacobian of the flux. For DG on unstructured meshes the CFL stability condition requires a CFL < 1/d to be satisfied (we refer to [51] for more details). Notice that the time step constraint does not need to be modified in presence of troubled cells, because we subdivide each troubled triangle in exactly $(2M+1)^d$ sub-triangles, and then we employ a FV scheme for which (14) holds with $M=0$.
<figure><img src="image_2.png" /><figcaption>Figure 1: The different possibilities available for the discretizations of domains with internal boundaries.</figcaption></figure>

# 4. Boundary conditions

In this section, we firstly discuss the classical approach to enforce boundary conditions in a DG framework. In particular, we focus on Dirichlet-type and slip-wall conditions. However, this technique is consistent as long as the boundary is also discretized with high order polynomials (see fig. 1a). When discretizing it with a linear conformal mesh, as in fig. 1b, this approach is at best second order accurate. Finally, by using fully embedded meshes, as in fig. 1c, the boundary consistency automatically drops to first order if no special treatments are considered. We then consider a more general slip-wall condition to design high order consistent boundary treatments for the cases described in fig. 1b-c.

## 4.1. Classical boundary treatment

When the boundary of element $K_j$ belongs to $\partial\Omega_h$, the normal flux function must account for the appropriate boundary conditions. In this work, the boundary flux is obtained by defining a ghost state $\mathcal{U}^{bc}$, and introducing a numerical flux $\hat{F}(\mathcal{U}_h^{-}, \mathcal{U}^{bc})$ defined by an approximate Riemann solver (8) based on the internal state $\mathcal{U}_h^{-}$ and on $\mathcal{U}^{bc}$. Depending on the condition to be enforced, different definitions of the ghost state are used. In this work we focus on

* Dirichlet-type conditions, where $\mathcal{U}^{bc}$ is fully known and enforced weakly through fluxes;

* reflecting conditions (slip-wall), where we set $\mathbf{u} \cdot \mathbf{n} = 0$, for static walls. In general for conformal meshes, $\mathcal{U}^{bc}$ has the same density, internal energy and tangential velocity of $\mathcal{U}_h^{-}$, and the opposite normal velocity component $\mathbf{u} \cdot \mathbf{n}$.

For conformal meshes, it is then consistent to take simply the face normal $\mathbf{n}$ as the normal of $\partial\Omega_h$ pointing outwards.

## 4.2. General slip-wall boundary treatment

When working on conformal curvilinear meshes, the reflecting boundary condition can be consistently formulated by imposing $\mathbf{u} \cdot \mathbf{n} = 0$ thanks to the high order approximation of the boundary (see fig. 1a). However, when approximating linearly or with general embedded meshes curved boundaries the condition imposed on the normal component of the velocity might not be enough due to the distance between $\partial\Omega$ and $\partial\Omega_h$ (see fig. 1b-c). For this reason, in this work we consider the more general boundary conditions introduced in [52] to have an accurate slip-wall treatment. Herein we recall the basic steps and hypothesis.
Starting from the mass and momentum in eq. (1), under some smoothness assumptions on the conserved variables we can write,

$$ \rho \frac{D\mathbf{u}}{Dt} + \nabla p = 0, \qquad (15) $$

where $\frac{D}{Dt} = \frac{\partial}{\partial t} + (\mathbf{u} \cdot \nabla)$ represents the material derivative.

The first condition to be imposed is on the normal component of the velocity vector for reflecting boundary conditions, meaning

$$ u_n = \mathbf{u} \cdot \mathbf{n} = 0, \qquad (16) $$

where $\mathbf{n}$ is now the normal of the real geometry, which may greatly differ from the normal of the discretized boundary $\partial\Omega_h$.

It should be noticed that in this work the real boundary normal is defined starting from a level set $\phi$ initialized as a signed distance function,

$$ \mathbf{n} := \frac{\nabla\phi}{\|\nabla\phi\|}, $$

and the real tangent vector is defined such that $\mathbf{n} \cdot \boldsymbol{\tau} = 0$, and $\boldsymbol{\tau} \times \mathbf{n} = 1$.

Condition (16) implies that we want to enforce the flow to be tangential to the real boundary $\partial\Omega$, meaning $\mathbf{u} = u_{\tau}\boldsymbol{\tau}$ which allows us to write,

$$ \frac{D\mathbf{u}}{Dt} = \boldsymbol{\tau} \frac{\partial u_{\tau}}{\partial t} + u_{\tau} \frac{\partial \boldsymbol{\tau}}{\partial t} + u_{\tau}(\mathbf{u} \cdot \nabla)\boldsymbol{\tau} = a_{\tau}\boldsymbol{\tau} + u_{\tau} \frac{D\boldsymbol{\tau}}{Dt}, \qquad (17) $$

where $a_{\tau}$ is the tangential acceleration. By using the Frenet-Serret formulae $\frac{D\boldsymbol{\tau}}{Dt} = u_{\tau}k\mathbf{n}$ (where $k$ is the local curvature), and projecting the momentum equation onto the normal direction, we end up with the second condition on the pressure,

$$ \frac{\partial p}{\partial n} = -\rho u_{\tau}^2 k. \qquad (18) $$

The third condition on the density is found by imposing an adiabatic wall, meaning $\partial_n S = 0$ where $S$ denotes the entropy. From $S = p\rho^{-\gamma}$ we can easily obtain $\partial_n S = (\partial_n p)\rho^{-\gamma} - \gamma p\rho^{-\gamma-1}(\partial_n \rho)$, and

$$ \frac{\partial \rho}{\partial n} = \frac{1}{c^2} \frac{\partial p}{\partial n}, \qquad (19) $$

where $c = \sqrt{\gamma p / \rho}$ is the speed of sound.

In order to derive the condition on the tangential component of the velocity, we consider the vorticity $\boldsymbol{\omega} := \nabla \times \mathbf{u}$ and compute it in the local coordinates $(\mathbf{n}, \boldsymbol{\tau})$ considering $\mathbf{u} = u_n\mathbf{n} + u_{\tau}\boldsymbol{\tau}$,

$$ \boldsymbol{\omega} = \frac{\partial u_{\tau}}{\partial n} - \frac{\partial u_n}{\partial \tau} - u_{\tau} k. $$

Finally, by assuming that the flow is irrotational and $u_n = 0$, we can deduce the needed condition on $u_{\tau}$,

$$ \frac{\partial u_{\tau}}{\partial n} = u_{\tau} k. \qquad (20) $$

In [52] it was also shown that the condition on $u_{\tau}$ can be retrieved by imposing $\partial_n H = 0$.
<figure><img src="image_3.png" /><figcaption>Figure 2: Visual representation of a curved boundary discretized with a linear conformal triangulation (a) and a fully embedded (b) triangulation.</figcaption></figure>

**Remark 2 (Equivalent boundary conditions).** From the previous statements, it is straightforward to deduce that the same boundary conditions can be imposed through the more compact system of mixed conditions:

$$u_n = 0, \quad \frac{\partial S}{\partial n} = 0, \quad \frac{\partial H}{\partial n} = 0, \quad \frac{\partial u_{\tau}}{\partial n} - u_{\tau}k = 0. \qquad (21)$$

Notice that by performing a change of variables from the conservative ones to $(u_n, S, H, u_{\tau})$ we can use the mixed Dirichlet-Neumann-Robin conditions (21) to compute the ghost state and impose the slip-wall boundary condition.

# 5. High order ROD boundary reconstruction

In this section, we present the numerical treatments we propose herein to impose general boundary conditions. We begin by presenting a reformulation of the ROD boundary reconstruction in space by employing a set of traditional basis function, which differ from those used for the discretization of the spatial operators. Afterwards, we introduce a novel space-time ROD reconstruction that exploits the space-time basis function used to define the high order local predictor (11). It is then shown that performing the ROD reconstruction by using the same basis functions used for the discretization allows us to greatly simplify the functional of the ROD minimization problem.

## 5.1. Reconstruction for off-site data methods in space

The Reconstruction for Off-site Data (ROD) method consists in discarding the polynomial of the boundary cell $K_j \in \partial\Omega_h$ and replacing it with a new reconstructed one that embeds the considered boundary condition. The ROD polynomial is then used to reconstruct the solution in each quadrature point $\mathbf{x}^*$ of the boundary integral on $\partial\Omega_h$. This reconstruction is performed through a constrained least-squares procedure, which is then solved with Lagrange multipliers.

In particular, for all boundary cells a new polynomial in space is defined as,

$$\varphi(\mathbf{x}, \eta, \mathbf{m}) = P_M(\mathbf{x} - \mathbf{m})\eta = \sum_{\alpha=0}^{M} \sum_{\beta=0}^{M-\alpha} (x_1 - m_1)^{\alpha} (x_2 - m_2)^{\beta} \eta_{(\alpha, \beta)}, \qquad (22)$$

where $\eta$ represents a set of coefficients, and $P_M(\mathbf{x})$ is a vector of polynomial basis of degree $M$. In this section, we consider a set of traditional basis functions that differ from those used for the
<figure><img src="image_4.png" /><figcaption>Figure 3: Visual representation of the ROD stencil: on the left, a constrained spatial reconstruction is performed at each sub-time node of the ADER correction step; on the right, one single space-time constrained polynomial is build by using all the information coming from the space-time predictor. ● – quadrature points in the space-time element; ● – polynomial coefficients of the local space-time predictor; ● – sampling points on the real boundary.</figcaption></figure>

PDE discretization,

$$ 
\begin{align} 
\mathcal{P}_1(x) &= [1 \quad x_1 \quad x_2], \nonumber \\
\mathcal{P}_2(x) &= [1 \quad x_1 \quad x_2 \quad x_1^2 \quad x_1x_2 \quad x_2^2], \tag{23} \\
\mathcal{P}_3(x) &= [1 \quad x_1 \quad x_2 \quad x_1^2 \quad x_1x_2 \quad x_2^2 \quad x_1^3 \quad x_1x_2 \quad x_1x_2^2 \quad x_2^3]. \nonumber
\end{align} 
$$

The point **m** is a reference point used for conditioning purposes, and coincides with the barycenter of the element. It is well-known that a two-dimensional polynomial reconstruction of order *M* requires *N* = (*M* + 1)(*M* + 2)/2 coefficients.

**Remark 3 (Comparison with the original approach).** In the ROD initial formulation for finite volume schemes [15], the size of the stencil was chosen *s* ≈ 1.5 *N*, and the corresponding neighboring elements were found using a searching algorithm.

When working in the DG framework, a local high order polynomial comes from the discretization, hence there is no need of finding the neighboring elements. The reconstructed ROD polynomial can easily embed both the information coming from the high order local polynomial and boundary conditions. It should be noticed that in this case, due to the ROD formulation in space only, the minimization problem must be formulated at each sub-time node of the considered time integration scheme, as shown in fig. 3a.

To build a set of sampling points required for the minimization problem, we exploit the local predictor and evaluate it in the space quadrature/collocation points $\mathbf{x}_q$ of the sub-time node $t^{n+\alpha_i}$ considered (refer again to fig. 3a). This gives us a local stencil of sampling points $\bar{\mathcal{U}}_q$ representing fully the high order behavior of the solution at $t^{n+\alpha_i}$.

Then, the ROD polynomial reconstruction consists in finding the polynomial coefficients $\hat{\eta}$ by minimizing a weighted cost function $\hat{\mathcal{J}}$ subject to certain linear constraints $\hat{\mathcal{G}}$ (the boundary
conditions), that reads

$$ \hat{\mathcal{J}}(\eta) = \frac{1}{2} \sum_q [w_q \mathcal{P}_M(\mathbf{x}_q - \mathbf{m})\eta - w_q \bar{\mathcal{U}}_q]^2, \quad (24) \\ \text{such that } \hat{\eta} = \arg\min_{\eta} \hat{\mathcal{J}}(\eta) \quad \text{subject to} \quad \hat{\mathcal{G}}(\eta) = 0. $$

It should be noticed that the mismatch between the basis functions $\mathcal{P}_M$ used for the ROD reconstruction and those used for the space discretization that define $\bar{\mathcal{U}}$ entails that a least squares problem in terms of quadrature points must be solved. In the unconstrained case, this can be solved in matrix form as follows,

$$ \underbrace{\left[ \sum_q w_q \mathcal{P}_M^T(\mathbf{x}_q - \mathbf{m}) \mathcal{P}_M(\mathbf{x}_q - \mathbf{m}) \right]}_{A_{N \times N}} \eta = \underbrace{\left[ \sum_q w_q \mathcal{P}_M^T(\mathbf{x}_q - \mathbf{m}) \bar{\mathcal{U}}_q \right]}_{b_{N \times 1}} \quad (25) $$

where $w_q := w(d_q, \delta, \gamma)$ is a weight that depends on the position of the sampling point $\mathbf{x}_q$. Usually it is set as the inverse of a distance to a power of $\gamma$:

$$ w(D_q, \delta, \gamma) = \frac{1}{(\delta D_q)^\gamma + 1}, $$

where $\delta$ is a sensibility factor and $D_q = |\mathbf{x}^* - \mathbf{x}_q|$ is the Euclidean distance between point $\mathbf{x}_q$ and the considered quadrature point $\mathbf{x}^*$.

In the constrained case, $\hat{\eta}$ must fulfill $p$ linear constraints, where $0 < p < N$. The ROD polynomial is finally obtained by applying these linear constraints in the matrix form $C\eta = \varrho$, where $C$ represents the general Robin boundary conditions $\alpha\varphi + \beta\partial_n\varphi$ imposed on the polynomials $\mathcal{P}_M$ computed on the real boundary $\Gamma = \partial\Omega$, and $\varrho$ is the set of conditions to enforce. These constraints depend on the nature of the boundary conditions and its geometry through the local normal $n$. Following this formulation, we can easily reformulate the linear constraints as follows,

$$ \underbrace{\left[ \begin{array}{c} \alpha_{\Gamma_1} \mathcal{P}_M(\mathbf{x}_{\Gamma_1} - \mathbf{m}) + \beta_{\Gamma_1} \nabla \mathcal{P}_M(\mathbf{x}_{\Gamma_1} - \mathbf{m}) \cdot \mathbf{n} \\ \vdots \\ \alpha_{\Gamma_p} \mathcal{P}_M(\mathbf{x}_{\Gamma_p} - \mathbf{m}) + \beta_{\Gamma_p} \nabla \mathcal{P}_M(\mathbf{x}_{\Gamma_p} - \mathbf{m}) \cdot \mathbf{n} \end{array} \right]}_{C_{p \times N}} \eta = \underbrace{\left[ \begin{array}{c} g_{\Gamma_1} \\ \vdots \\ g_{\Gamma_p} \end{array} \right]}_{\varrho_{p \times 1}}, \quad (26) $$

where point $\mathbf{x}_{\Gamma_p}$ is defined as the projection of the quadrature point $\mathbf{x}^*$ on the real boundary $\partial\Omega$, along the normal $n$. It is straightforward then to retrieve classical Dirichlet ($\beta_{\Gamma_p} = 0$) and Neumann ($\alpha_{\Gamma_p} = 0$) boundary conditions.

As mentioned above this constrained optimization problem can be solved using Lagrange multipliers. Therefore, we define the functional

$$ \mathcal{L}(\eta, \lambda) = \frac{1}{2} \sum_q [w_q \mathcal{P}_M(\mathbf{x}_q - \mathbf{m})\eta - w_q \bar{\mathcal{U}}_q]^2 + \sum_{k=1}^p \lambda_k [\mathcal{C}\eta - \varrho]. \quad (27) $$

The result of this minimization problem $\hat{\eta}$ is the solution of the linear system

$$ \nabla_{\eta, \lambda} \mathcal{L}(\eta, \lambda) = 0, $$
where the differential operator $\nabla_{\eta,\lambda}$ is built by deriving with respect to each polynomial coefficient and Lagrange multiplier:

$$ \left\{ \begin{array}{l} \nabla_{\eta} \mathcal{L}(\eta, \lambda) = \mathcal{A}\eta - b + \mathcal{C}^T\lambda = 0 \\ \nabla_{\lambda} \mathcal{L}(\eta, \lambda) = \mathcal{C}\eta - \varrho = 0. \end{array} \right. \qquad (28) $$

This problem can be solved in matrix form, as follows

$$ \left[ \frac{\mathcal{A}}{\mathcal{C}} \middle| \frac{\mathcal{C}^T}{0} \right] \left[ \frac{\eta}{\lambda} \right] = \left[ \frac{b}{\varrho} \right] \quad \Rightarrow \quad \left[ \frac{\hat{\eta}}{\hat{\lambda}} \right] = \left[ \frac{\mathcal{A}}{\mathcal{C}} \middle| \frac{\mathcal{C}^T}{0} \right]^{-1} \left[ \frac{b}{\varrho} \right], \qquad (29) $$

where $\hat{\eta}$ represents the polynomial coefficients associated to the Lagrange multipliers $\hat{\lambda}$, which are also a solution of the linear system, that satisfies the boundary conditions.

## 5.2. Reconstruction for off-site data methods for space-time schemes

The original ROD method presented in the previous section is formulated only in space, which entails that the minimization problem must be solved at each sub-time node $t^{n+\alpha_i}$ of the correction step. A similar thing would also happen if dealing with Runge-Kutta methods, meaning that an inversion of the linear system is needed at each sub-time node for every boundary element. In this section, we introduce a novel space-time formulation for the ROD reconstructed polynomial that can easily be integrated in the correction step of the ADER method. Here, we exploit the ADER formulation based on the previous evaluation of the high order space-time predictor (11) to design a new minimization problem that use the same space-time polynomial basis. This idea not only allows to write a much simpler optimization problem but also implies the reconstruction of a single space-time ROD polynomial (meaning a single inversion of the linear system) used for the boundary conditions of the whole space-time element $K_j \times [t^n, t^{n+1}]$.

In this case, the fact that the same polynomial basis of the discretization is also used to perform the ROD boundary reconstruction allows us to write a simplified cost function only with respect to the polynomial coefficients, that no longer depends on the evaluation of the solution within the quadrature points,

$$ \begin{gathered} \hat{\mathcal{J}}(\eta) = \frac{1}{2} [\eta - \hat{q}]^2, \qquad (30) \\ \text{such that} \quad \hat{\eta} = \arg \min_{\eta} \hat{\mathcal{J}}(\eta) \quad \text{subject to} \quad \hat{\mathcal{G}}(\eta) = 0, \end{gathered} $$

where $\hat{q}$ are the predictor coefficients found with eq. (11).

The functional embedding the contraints can now be simplified as

$$ \mathcal{L}(\eta, \lambda) = \frac{1}{2} [\eta - \hat{q}]^2 + \sum_{k=1}^{p} \lambda_k [\mathcal{C}\eta - \varrho], \qquad (31) $$

from which we can deduce the system minimum as,

$$ \left\{ \begin{array}{l} \nabla_{\eta} \mathcal{L}(\eta, \lambda) = \mathcal{1}\eta - \hat{q} + \mathcal{C}^T\lambda = 0 \\ \nabla_{\lambda} \mathcal{L}(\eta, \lambda) = \mathcal{C}\eta - \varrho = 0. \end{array} \right. \qquad (32) $$

The previous system can be solved once again in matrix form, as follows

$$ \left[ \frac{\mathcal{1}}{\mathcal{C}} \middle| \frac{\mathcal{C}^T}{0} \right] \left[ \frac{\eta}{\lambda} \right] = \left[ \frac{\hat{q}}{\varrho} \right] \quad \Rightarrow \quad \left[ \frac{\hat{\eta}}{\hat{\lambda}} \right] = \left[ \frac{\mathcal{1}}{\mathcal{C}} \middle| \frac{\mathcal{C}^T}{0} \right]^{-1} \left[ \frac{\hat{q}}{\varrho} \right]. \qquad (33) $$
By using the space-time polynomial basis we are able to simplify the minimization problem, as now $\mathcal{A} = 1$.

**Remark 4.** It should be noticed that the linear system (33) is slightly larger than the previous one (29), but it is more diagonally-dominant, and has to be solved only once, which makes it more efficient to invert. This is also due to the fact that the total number of modes reconstructed is much lower than in the case of several independent reconstructions. Moreover, by reconstructing a single polynomial for the entire space-time boundary element, the number of required constraints to achieve consistency is halved.

### 5.3. Remarks on the implementation of slip-wall boundary conditions

In this section we address some implementation features needed to develop the slip-wall boundary condition (21) with both the space-only and space-time ROD reconstructions.

#### 5.3.1. Change of variables

As mentioned above, for both the original and space-time ROD reconstructions a change of variables is needed to enforce the mixed Dirichlet-Neumann-Robin boundary conditions, meaning that the ROD polynomial reconstruction must be performed in terms of $\mathcal{V}(\mathcal{U}) = (u_n, S, H, u_\tau)$ so that the conditions can be enforced linearly. When using the original ROD approach, the change of variables is straightforward due to the local evaluation of the sampling points within the cell quadrature points $\mathbf{x}_q$.

However, when working with the space-time formulation only the space-time polynomial coefficients are needed for the minimization problem. In general such coefficients are defined with respect to the vector of conserved variables, meaning that a projection is needed to retrieve them for the new vector of variables $\mathcal{V}$.

An accurate and efficient way to evaluate them is to perform a space-time $L_2$-projection,

$$ \left( \int_{t^n}^{t^{n+1}} \int_{K_j} \theta_i(\mathbf{x}, t) \theta_\ell(\mathbf{x}, t) \, d\mathbf{x}dt \right) \hat{\mathcal{V}}_\ell = \int_{t^n}^{t^{n+1}} \int_{K_j} \theta_i(\mathbf{x}, t) \mathcal{V}(\mathcal{U}_h(\mathbf{x}, t)) \, d\mathbf{x}dt, \quad (34) $$

such that $\mathcal{V}_h$ can be defined as,

$$ \mathcal{V}_h(\mathbf{x}, t)|_{K_j \times [t^n, t^{n+1}]} = \sum_{\ell} \theta_{\ell}(\mathbf{x}, t) \hat{\mathcal{V}}_{\ell}. \quad (35) $$

**Remark 5.** It should be noticed that, in order to solve problem (34), $\mathcal{V}(\mathcal{U}_h(\mathbf{x}, t))$ is evaluated in the space-time element making necessary the knowledge of the boundary normals and tangents also within the computational domain. In our case, since the interface is defined from a level set field, the normal and tangents can be easily computed also far from the boundary.

#### 5.3.2. Shock limiting

For both the ROD formulations, when the a posteriori limiter marks a cell on the wall boundary as troubled, the boundary conditions must be applied differently because the high order polynomials are discarded, and the solution is updated with a second order finite volume scheme.

For simplicity, we adopt a similar flux to the one introduced in [13], which is computed including the information related to the real boundary normal. We take the decomposition $\tilde{\mathbf{n}} = (\tilde{\mathbf{n}} \cdot \mathbf{n})\mathbf{n} + (\tilde{\mathbf{n}} \cdot \boldsymbol{\tau})\boldsymbol{\tau}$, where $\tilde{\mathbf{n}}$ represents the normal of the discretized curved boundary $\partial\Omega_h$ and, $\mathbf{n}$ and $\boldsymbol{\tau}$
represent the real normal and tangent vectors to $\partial\Omega$, respectively. Therefore by writing the flux on the discretized boundary $\mathcal{F}_{\tilde{n}}^{bc} = (\tilde{n} \cdot \mathbf{n})\mathcal{F}_n + (\tilde{n} \cdot \boldsymbol{\tau})\mathcal{F}_{\boldsymbol{\tau}}$, and taking $u_n = 0$, we can recast

$$ \mathcal{F}_{\tilde{n}}^{bc} = \begin{pmatrix} (\tilde{n} \cdot \boldsymbol{\tau}) \rho u_{\tau} \\ (\tilde{n} \cdot \boldsymbol{\tau}) \rho \mathbf{u}_{\tau} + p \tilde{n} \\ (\tilde{n} \cdot \boldsymbol{\tau}) \rho u_{\tau} H \end{pmatrix}. \qquad (36) $$

We remind that eq. (36) is at most third order accurate for general geometries discretized with conformal linear meshes [15], which makes it consistent for troubled cells.

## 6. Numerical experiments

In this section, we present the numerical experiments performed to validate the boundary reconstructions presented in section 5. In the first example we validate the high order convergence properties of the boundary treatments with a smooth manufactured solution characterized by a curved boundary. The adaptation of the space-only ROD reconstruction in the ADER-DG framework will be referred to as ROD-DG; while the novel formulation based on the ADER space-time predictor will be referred to as ROD-ADER-DG. The numerical results obtained without any correction will be simply referred to as DG. For the manufactured solution test we present the numerical solutions obtained by imposing Dirichlet-type boundary conditions on the curved boundary discretized with both conformal linear elements and fully embedded meshes. After that, we present several numerical experiments showing the capabilities of the new approach to deal with both Dirichlet and slip-wall boundary conditions, also when coupled with an a posteriori shock limiter. Once proven for the first example that the two formulations provide very similar discretization errors, for the other test cases we only focus on the results obtained with our new ROD-ADER-DG method.

### 6.1. Manufactured solution: comparisons between the original and space-time reconstruction

The first test case is dedicated to the comparisons between the two ROD formulations described in section 5. In particular we focus on a smooth manufactured solution by considering the two-dimensional inhomogeneous Euler equations (4) with,

$$ S = (0.4 \cos(x + y), 0.6 \cos(x + y), 0.6 \cos(x + y), 1.8 \cos(x + y))^T. $$

This system has the following exact steady state solution,

$$ \rho = 1 + 0.2 \sin(x + y), \quad u = (1, 1), \quad p = 1 + 0.2 \sin(x + y), $$

which is imposed on the real curved boundary as Dirichlet-type boundary conditions. A very coarse mesh is generated and then refined by splitting. The four grids described in table 1 are used to study the convergence properties of the formulations. The boundary conditions are imposed on a complicated curved geometry (see fig. 4), that can be described with the following equation written in polar coordinates:

$$ \begin{pmatrix} x \\ y \end{pmatrix} = r(\alpha, \theta) \begin{pmatrix} \cos \theta \\ \sin \theta \end{pmatrix}, \quad \text{where} \quad r(\alpha, \theta) = r_0 \left( 1 + \frac{1}{10} \sin(\alpha \theta) \right), \qquad (37) $$

where $r_0 = 1$ and $\alpha = 3$.
<figure><img src="image_5.png" /><figcaption>Figure 4: Manufactured solution: setup for the test case presented in section 6.1. We plot the coarsest employed conformal mesh, and the initial density profile.</figcaption></figure>

<table><caption>Table 1: Manufactured solution: characteristics of the conformal linear meshes for the test case presented in section 6.1.</caption><thead><tr><th>Grid level</th><th>Nodes</th><th>Triangles</th><th><i>h</i></th></tr></thead><tbody><tr><td>0</td><td>61</td><td>98</td><td>1.54E-1</td></tr><tr><td>1</td><td>219</td><td>392</td><td>7.72E-2</td></tr><tr><td>2</td><td>829</td><td>1,568</td><td>3.86E-2</td></tr><tr><td>3</td><td>3,225</td><td>6,272</td><td>1.93E-2</td></tr></tbody></table>

<table><caption>Table 2: Manufactured solution: convergence analysis for the test case presented in section 6.1. Numerical results obtained with Dirichlet-type boundary conditions without any correction on the conformal linear meshes described in table 1. One can notice that by imposing naively the boundary conditions, the method is at most second order accurate.</caption><thead><tr><th rowspan="2">Grid level</th><th colspan="2">ρ</th><th colspan="2">ρu</th><th colspan="2">ρv</th><th colspan="2">ρE</th></tr><tr><th>L<sub>2</sub></th><th>ñ</th><th>L<sub>2</sub></th><th>ñ</th><th>L<sub>2</sub></th><th>ñ</th><th>L<sub>2</sub></th><th>ñ</th></tr></thead><tbody><tr><td colspan="9" style="text-align:center; font-style:italic;">DG-<i>P</i><sub>1</sub></td></tr><tr><td>0</td><td>1.7293E-3</td><td>-</td><td>2.3514E-3</td><td>-</td><td>2.2344E-3</td><td>-</td><td>6.6466E-3</td><td>-</td></tr><tr><td>1</td><td>4.4790E-4</td><td>1.95</td><td>6.0925E-4</td><td>1.95</td><td>5.7869E-4</td><td>1.95</td><td>1.7138E-3</td><td>1.96</td></tr><tr><td>2</td><td>1.1276E-4</td><td>1.99</td><td>1.5358E-4</td><td>1.99</td><td>1.4579E-4</td><td>1.99</td><td>4.3195E-4</td><td>1.99</td></tr><tr><td>3</td><td>2.8241E-5</td><td>2.00</td><td>3.8473E-5</td><td>2.00</td><td>3.6512E-5</td><td>2.00</td><td>1.0820E-4</td><td>2.00</td></tr><tr><td colspan="9" style="text-align:center; font-style:italic;">DG-<i>P</i><sub>2</sub></td></tr><tr><td>0</td><td>1.6583E-3</td><td>-</td><td>2.3344E-3</td><td>-</td><td>2.2040E-3</td><td>-</td><td>6.4258E-3</td><td>-</td></tr><tr><td>1</td><td>4.2195E-4</td><td>1.97</td><td>5.9156E-4</td><td>1.98</td><td>5.5908E-4</td><td>1.98</td><td>1.6308E-3</td><td>1.98</td></tr><tr><td>2</td><td>1.0522E-4</td><td>2.00</td><td>1.4777E-4</td><td>2.00</td><td>1.3975E-4</td><td>2.00</td><td>4.0791E-4</td><td>2.00</td></tr><tr><td>3</td><td>2.6187E-5</td><td>2.01</td><td>3.6855E-5</td><td>2.00</td><td>3.4867E-5</td><td>2.00</td><td>1.0170E-4</td><td>2.00</td></tr><tr><td colspan="9" style="text-align:center; font-style:italic;">DG-<i>P</i><sub>3</sub></td></tr><tr><td>0</td><td>1.6817E-3</td><td>-</td><td>2.3809E-3</td><td>-</td><td>2.2562E-3</td><td>-</td><td>6.5431E-3</td><td>-</td></tr><tr><td>1</td><td>4.3095E-4</td><td>1.96</td><td>6.0642E-4</td><td>1.97</td><td>5.7430E-4</td><td>1.97</td><td>1.6640E-3</td><td>1.98</td></tr><tr><td>2</td><td>1.0768E-4</td><td>2.00</td><td>1.5092E-4</td><td>2.01</td><td>1.4280E-4</td><td>2.01</td><td>4.1580E-4</td><td>2.00</td></tr><tr><td>3</td><td>2.6792E-5</td><td>2.01</td><td>3.7467E-5</td><td>2.01</td><td>3.5428E-5</td><td>2.01</td><td>1.0345E-4</td><td>2.01</td></tr></tbody></table>
<table><caption>Table 3: Manufactured solution: characteristics of the fully embedded meshes on a square [-1.5, 1.5] x [-1.5, 1.5] for the test case presented in section 6.1.</caption><thead><tr><td>Grid level</td><td>h</td></tr></thead><tbody><tr><td>5 x 5</td><td>3.79E-1</td></tr><tr><td>10 x 10</td><td>2.02E-1</td></tr><tr><td>15 x 15</td><td>1.33E-1</td></tr><tr><td>20 x 20</td><td>9.93E-2</td></tr><tr><td>35 x 35</td><td>6.02E-2</td></tr></tbody></table>

<figure><img src="image_6.png" /><figcaption>Figure 5: Manufactured solution: setup for the test case presented in section 6.1. We plot the employed embedded linear meshes on a square [-1.5, 1.5] x [-1.5, 1.5], and the initial density profile.</figcaption></figure>

Since the standard Dirichlet boundary condition is enforced onto a curved boundary, discretized with a polygonal mesh, for the DG solution we expect to have at most second order of accuracy no matter the order of the polynomial used. This is clear from table 2 where the $L_2$ error norm is shown for the different experiments.

We then perform the same analysis with both the original and space-time ROD formulations to prove numerically their high order convergence properties. We present the numerical results in tables 4 to 7 obtained with both formulations on conformal linear meshes and fully embedded meshes. In all cases, the expected convergence of $M + 1$ is achieved. Moreover, as mentioned above, the new space-time formulation is able to achieve similar results while being less computationally expensive. In particular, we experience a reduction of ~10% for $P_3$ polynomials, which is remarkable given that the ROD reconstruction only influences the boundary cells. It should be noticed that the convergence trends obtained by the method on fully embedded meshes (described in table 3) are typical of such methods [16], simply because by refining the mesh the discretized boundary gets closer to the real one depending on how it cuts the background mesh (see fig. 5).
<table><caption>Table 4: Manufactured solution: convergence analysis for the test case presented in section 6.1. Numerical results obtained with the ROD-DG reconstruction on the conformal linear meshes described in table 1. One can notice that the high order precision is achieved for all polynomials.</caption><thead><tr><th rowspan="2">Grid level</th><th colspan="2">ρ</th><th colspan="2">ρu</th><th colspan="2">ρv</th><th colspan="2">ρE</th></tr><tr><th>L<sub>2</sub></th><th>ñ</th><th>L<sub>2</sub></th><th>ñ</th><th>L<sub>2</sub></th><th>ñ</th><th>L<sub>2</sub></th><th>ñ</th></tr></thead><tbody><tr><td colspan="9">ROD-DG-P<sub>1</sub></td></tr><tr><td>0</td><td>9.0353E-4</td><td>-</td><td>9.3786E-4</td><td>-</td><td>9.3503E-4</td><td>-</td><td>2.8893E-3</td><td>-</td></tr><tr><td>1</td><td>2.0990E-4</td><td>2.11</td><td>2.0792E-4</td><td>2.17</td><td>2.0892E-4</td><td>2.16</td><td>6.6542E-4</td><td>2.12</td></tr><tr><td>2</td><td>5.0528E-5</td><td>2.05</td><td>4.9278E-5</td><td>2.08</td><td>4.9736E-5</td><td>2.07</td><td>1.6082E-4</td><td>2.05</td></tr><tr><td>3</td><td>1.2438E-5</td><td>2.02</td><td>1.2032E-5</td><td>2.03</td><td>1.2173E-5</td><td>2.03</td><td>3.9724E-5</td><td>2.02</td></tr><tr><td colspan="9">ROD-DG-P<sub>2</sub></td></tr><tr><td>0</td><td>5.8327E-5</td><td>-</td><td>4.0652E-5</td><td>-</td><td>4.5945E-5</td><td>-</td><td>1.7771E-4</td><td>-</td></tr><tr><td>1</td><td>9.5321E-6</td><td>2.61</td><td>5.1266E-6</td><td>2.99</td><td>6.2677E-6</td><td>2.87</td><td>2.8402E-5</td><td>2.65</td></tr><tr><td>2</td><td>1.4695E-6</td><td>2.70</td><td>6.9200E-7</td><td>2.89</td><td>8.6795E-7</td><td>2.85</td><td>4.3499E-6</td><td>2.71</td></tr><tr><td>3</td><td>2.0870E-7</td><td>2.82</td><td>9.7264E-8</td><td>2.83</td><td>1.1487E-7</td><td>2.92</td><td>6.1764E-7</td><td>2.82</td></tr><tr><td colspan="9">ROD-DG-P<sub>3</sub></td></tr><tr><td>0</td><td>1.0030E-6</td><td>-</td><td>1.3555E-6</td><td>-</td><td>1.4008E-6</td><td>-</td><td>3.4352E-6</td><td>-</td></tr><tr><td>1</td><td>4.7755E-8</td><td>4.39</td><td>5.4425E-8</td><td>4.64</td><td>5.5114E-8</td><td>4.67</td><td>1.5203E-7</td><td>4.50</td></tr><tr><td>2</td><td>2.5871E-9</td><td>4.21</td><td>2.7218E-9</td><td>4.32</td><td>2.7433E-9</td><td>4.33</td><td>8.1616E-9</td><td>4.22</td></tr><tr><td>3</td><td>1.5035E-10</td><td>4.10</td><td>1.5323E-10</td><td>4.15</td><td>1.5474E-10</td><td>4.15</td><td>4.7689E-10</td><td>4.10</td></tr></tbody></table>

<table><caption>Table 5: Manufactured solution: convergence analysis for the test case presented in section 6.1. Numerical results obtained with the ROD-DG reconstruction on the fully embedded meshes presented in fig. 5. One can notice that the high order precision is achieved for all polynomials.</caption><thead><tr><th rowspan="2">Grid level</th><th colspan="2">ρ</th><th colspan="2">ρu</th><th colspan="2">ρv</th><th colspan="2">ρE</th></tr><tr><th>L<sub>2</sub></th><th>ñ</th><th>L<sub>2</sub></th><th>ñ</th><th>L<sub>2</sub></th><th>ñ</th><th>L<sub>2</sub></th><th>ñ</th></tr></thead><tbody><tr><td colspan="9">ROD-DG-P<sub>1</sub></td></tr><tr><td>5 × 5</td><td>0.1840E-01</td><td>-</td><td>0.2589E-01</td><td>-</td><td>0.2369E-01</td><td>-</td><td>0.7320E-01</td><td>-</td></tr><tr><td>10 × 10</td><td>0.4015E-02</td><td>2.26</td><td>0.5704E-02</td><td>2.24</td><td>0.5326E-02</td><td>2.21</td><td>0.1566E-01</td><td>2.28</td></tr><tr><td>15 × 15</td><td>0.2141E-02</td><td>1.64</td><td>0.3118E-02</td><td>1.57</td><td>0.3067E-02</td><td>1.44</td><td>0.8588E-02</td><td>1.56</td></tr><tr><td>20 × 20</td><td>0.1732E-02</td><td>0.74</td><td>0.2453E-02</td><td>0.83</td><td>0.2602E-02</td><td>0.57</td><td>0.7137E-02</td><td>0.64</td></tr><tr><td>35 × 35</td><td>0.4966E-03</td><td>2.26</td><td>0.7377E-03</td><td>2.17</td><td>0.7078E-03</td><td>2.36</td><td>0.1996E-02</td><td>2.30</td></tr><tr><td colspan="9">ROD-DG-P<sub>2</sub></td></tr><tr><td>5 × 5</td><td>0.4464E-02</td><td>-</td><td>0.7350E-02</td><td>-</td><td>0.6120E-02</td><td>-</td><td>0.1918E-01</td><td>-</td></tr><tr><td>10 × 10</td><td>0.2611E-03</td><td>4.21</td><td>0.3919E-03</td><td>4.34</td><td>0.3617E-03</td><td>4.19</td><td>0.1029E-02</td><td>4.33</td></tr><tr><td>15 × 15</td><td>0.9730E-04</td><td>2.57</td><td>0.1514E-03</td><td>2.48</td><td>0.1499E-03</td><td>2.29</td><td>0.3833E-03</td><td>2.57</td></tr><tr><td>20 × 20</td><td>0.5464E-04</td><td>2.00</td><td>0.7456E-04</td><td>2.46</td><td>0.8333E-04</td><td>2.04</td><td>0.2213E-03</td><td>1.90</td></tr><tr><td>35 × 35</td><td>0.8168E-05</td><td>3.44</td><td>0.1196E-04</td><td>3.31</td><td>0.1057E-04</td><td>3.73</td><td>0.3169E-04</td><td>3.52</td></tr><tr><td colspan="9">ROD-DG-P<sub>3</sub></td></tr><tr><td>5 × 5</td><td>0.5828E-03</td><td>-</td><td>0.1035E-02</td><td>-</td><td>0.9056E-03</td><td>-</td><td>0.2743E-02</td><td>-</td></tr><tr><td>10 × 10</td><td>0.3038E-04</td><td>4.38</td><td>0.4969E-04</td><td>4.50</td><td>0.4695E-04</td><td>4.38</td><td>0.1252E-03</td><td>4.57</td></tr><tr><td>15 × 15</td><td>0.7848E-05</td><td>3.52</td><td>0.1435E-04</td><td>3.23</td><td>0.1368E-04</td><td>3.21</td><td>0.3360E-04</td><td>3.43</td></tr><tr><td>20 × 20</td><td>0.4749E-05</td><td>1.74</td><td>0.7741E-05</td><td>2.14</td><td>0.7982E-05</td><td>1.87</td><td>0.2019E-04</td><td>1.77</td></tr><tr><td>35 × 35</td><td>0.4305E-06</td><td>4.34</td><td>0.7545E-06</td><td>4.21</td><td>0.7202E-06</td><td>4.35</td><td>0.1822E-05</td><td>4.35</td></tr></tbody></table>
<table><caption>Table 6: Manufactured solution: convergence analysis for the test case presented in section 6.1. Numerical results obtained with the ROD-ADER-DG reconstruction on the conformal linear meshes described in table 1. One can notice that the high order precision is achieved for all polynomials.</caption><thead><tr><th rowspan="2">Grid level</th><th colspan="2">ρ</th><th colspan="2">ρu</th><th colspan="2">ρv</th><th colspan="2">ρE</th></tr><tr><th>L<sub>2</sub></th><th>ñ</th><th>L<sub>2</sub></th><th>ñ</th><th>L<sub>2</sub></th><th>ñ</th><th>L<sub>2</sub></th><th>ñ</th></tr></thead><tbody><tr><td colspan="9">ROD-ADER-DG-P<sub>1</sub></td></tr><tr><td>0</td><td>8.3292E-4</td><td>-</td><td>8.2968E-4</td><td>-</td><td>8.2849E-4</td><td>-</td><td>2.6472E-3</td><td>-</td></tr><tr><td>1</td><td>2.0273E-4</td><td>2.04</td><td>1.9689E-4</td><td>2.07</td><td>1.9822E-4</td><td>2.06</td><td>6.4251E-4</td><td>2.04</td></tr><tr><td>2</td><td>4.9613E-5</td><td>2.03</td><td>4.7793E-5</td><td>2.04</td><td>4.8286E-5</td><td>2.04</td><td>1.5771E-4</td><td>2.03</td></tr><tr><td>3</td><td>1.2295E-5</td><td>2.01</td><td>1.1778E-5</td><td>2.02</td><td>1.1921E-5</td><td>2.02</td><td>3.9159E-5</td><td>2.01</td></tr><tr><td colspan="9">ROD-ADER-DG-P<sub>2</sub></td></tr><tr><td>0</td><td>5.7478E-5</td><td>-</td><td>3.9499E-5</td><td>-</td><td>4.4647E-5</td><td>-</td><td>1.7488E-4</td><td>-</td></tr><tr><td>1</td><td>9.4589E-6</td><td>2.60</td><td>5.0413E-6</td><td>2.97</td><td>6.1664E-6</td><td>2.86</td><td>2.8118E-5</td><td>2.64</td></tr><tr><td>2</td><td>1.4611E-6</td><td>2.70</td><td>6.8355E-7</td><td>2.88</td><td>8.5622E-7</td><td>2.85</td><td>4.3155E-6</td><td>2.70</td></tr><tr><td>3</td><td>2.0771E-7</td><td>2.81</td><td>9.6318E-8</td><td>2.83</td><td>1.1345E-7</td><td>2.92</td><td>6.1345E-7</td><td>2.81</td></tr><tr><td colspan="9">ROD-ADER-DG-P<sub>3</sub></td></tr><tr><td>0</td><td>8.1215E-7</td><td>-</td><td>1.0517E-6</td><td>-</td><td>1.0822E-6</td><td>-</td><td>2.8396E-6</td><td>-</td></tr><tr><td>1</td><td>4.3856E-8</td><td>4.21</td><td>4.9314E-8</td><td>4.41</td><td>5.0113E-8</td><td>4.43</td><td>1.4269E-7</td><td>4.31</td></tr><tr><td>2</td><td>2.4815E-9</td><td>4.14</td><td>2.5941E-9</td><td>4.25</td><td>2.6192E-9</td><td>4.25</td><td>7.9180E-9</td><td>4.17</td></tr><tr><td>3</td><td>1.4723E-10</td><td>4.07</td><td>1.4964E-10</td><td>4.12</td><td>1.5109E-10</td><td>4.11</td><td>4.7004E-10</td><td>4.07</td></tr></tbody></table>

<table><caption>Table 7: Manufactured solution: convergence analysis for the test case presented in section 6.1. Numerical results obtained with the ROD-ADER-DG reconstruction on the fully embedded meshes presented in fig. 5. One can notice that the high order precision is achieved for all polynomials.</caption><thead><tr><th rowspan="2">Grid level</th><th colspan="2">ρ</th><th colspan="2">ρu</th><th colspan="2">ρv</th><th colspan="2">ρE</th></tr><tr><th>L<sub>2</sub></th><th>ñ</th><th>L<sub>2</sub></th><th>ñ</th><th>L<sub>2</sub></th><th>ñ</th><th>L<sub>2</sub></th><th>ñ</th></tr></thead><tbody><tr><td colspan="9">ROD-ADER-DG-P<sub>1</sub></td></tr><tr><td>5 × 5</td><td>2.1244E-2</td><td>-</td><td>2.8921E-2</td><td>-</td><td>2.6676E-2</td><td>-</td><td>8.2175E-2</td><td>-</td></tr><tr><td>10 × 10</td><td>4.7323E-3</td><td>2.39</td><td>6.9117E-3</td><td>2.28</td><td>6.4829E-3</td><td>2.25</td><td>1.9135E-2</td><td>2.32</td></tr><tr><td>15 × 15</td><td>2.5373E-3</td><td>1.51</td><td>3.7614E-3</td><td>1.47</td><td>3.7089E-3</td><td>1.35</td><td>1.0495E-2</td><td>1.45</td></tr><tr><td>20 × 20</td><td>1.9806E-3</td><td>0.83</td><td>2.8678E-3</td><td>0.91</td><td>3.0124E-3</td><td>0.70</td><td>8.3544E-3</td><td>0.76</td></tr><tr><td>35 × 35</td><td>5.7732E-4</td><td>2.46</td><td>8.7147E-4</td><td>2.38</td><td>8.3663E-4</td><td>2.26</td><td>2.3870E-3</td><td>2.50</td></tr><tr><td colspan="9">ROD-ADER-DG-P<sub>2</sub></td></tr><tr><td>5 × 5</td><td>6.4776E-3</td><td>-</td><td>1.0806E-2</td><td>-</td><td>8.9916E-3</td><td>-</td><td>2.8161E-2</td><td>-</td></tr><tr><td>10 × 10</td><td>4.3917E-4</td><td>4.29</td><td>6.7778E-4</td><td>4.41</td><td>6.4805E-4</td><td>4.18</td><td>1.8388E-3</td><td>4.34</td></tr><tr><td>15 × 15</td><td>1.6884E-4</td><td>2.31</td><td>2.6535E-4</td><td>2.26</td><td>2.6897E-4</td><td>2.12</td><td>7.0304E-4</td><td>2.32</td></tr><tr><td>20 × 20</td><td>8.9175E-5</td><td>2.15</td><td>1.2627E-4</td><td>2.50</td><td>1.3957E-4</td><td>2.20</td><td>3.7224E-4</td><td>2.14</td></tr><tr><td>35 × 35</td><td>1.3315E-5</td><td>3.79</td><td>2.0681E-5</td><td>3.61</td><td>1.7738E-5</td><td>4.12</td><td>5.4102E-5</td><td>3.85</td></tr><tr><td colspan="9">ROD-ADER-DG-P<sub>3</sub></td></tr><tr><td>5 × 5</td><td>7.4636E-4</td><td>-</td><td>1.3662E-3</td><td>-</td><td>1.1646E-3</td><td>-</td><td>3.6090E-3</td><td>-</td></tr><tr><td>10 × 10</td><td>7.0291E-5</td><td>3.76</td><td>1.0730E-4</td><td>4.05</td><td>1.0734E-4</td><td>3.80</td><td>2.9700E-4</td><td>3.98</td></tr><tr><td>15 × 15</td><td>1.2500E-5</td><td>4.17</td><td>2.1671E-5</td><td>3.86</td><td>2.2876E-5</td><td>3.73</td><td>5.7446E-5</td><td>3.97</td></tr><tr><td>20 × 20</td><td>7.9376E-6</td><td>1.53</td><td>1.3160E-5</td><td>1.68</td><td>1.2447E-5</td><td>2.05</td><td>3.2858E-5</td><td>1.88</td></tr><tr><td>35 × 35</td><td>8.4487E-7</td><td>4.47</td><td>1.4199E-6</td><td>4.45</td><td>1.3451E-6</td><td>4.44</td><td>3.5280E-6</td><td>4.46</td></tr></tbody></table>
<figure><img src="image_7.png" /><figcaption>Figure 6: Isentropic supersonic vortex: setup for the test case presented in section 6.2. We plot a coarse mesh and the initial density profile.</figcaption></figure>

<table><caption>Table 8: Isentropic supersonic vortex: characteristics of the conformal linear meshes for the test case presented in section 6.2.</caption><thead><tr><th>Grid level</th><th>Nodes</th><th>Triangles</th><th>h</th></tr></thead><tbody><tr><td>0</td><td>71</td><td>92</td><td>1.50E-1</td></tr><tr><td>1</td><td>238</td><td>376</td><td>7.50E-2</td></tr><tr><td>2</td><td>852</td><td>1,504</td><td>3.75E-2</td></tr><tr><td>3</td><td>3,208</td><td>6,016</td><td>1.87E-2</td></tr></tbody></table>

## 6.2. Isentropic supersonic vortex: validation of Dirichlet-type and slip-wall boundary conditions

We now consider an isentropic supersonic flow between two concentric circular arcs of radii $r_i = 1$ and $r_o = 1.384$. The exact density in terms of radius $r$ is given by

$$ \rho = \rho_i \left( 1 + \frac{\gamma + 1}{2} M_i^2 \left( 1 - \left( \frac{r_i}{r} \right)^2 \right) \right)^{\frac{1}{\gamma - 1}}. \quad (38) $$

The velocity and pressure are given by

$$ \|u\| = \frac{c_i M_i}{r}, \quad p = \frac{\rho^\gamma}{\gamma}, \quad (39) $$

where $c_i$ is the speed of sound on the inner circle. The Mach number on the inner circle $M_i$ is set to 2.25 and the density $\rho_i$ to 1.

In this section we present only the results obtained with the novel formulation, which has already shown better performances, for different boundary conditions. Here we present the results obtained on four refined conformal triangulations. In fig. 6 the computational domain and a coarse triangulation are shown.

Since for this test case we have an analytical solution, we present in table 9 the results of the experiment performed by enforcing Dirichlet-type conditions on the internal boundary, and classical slip-wall conditions on the external boundary without any correction. As expected, since
<table><caption>Table 9: Isentropic supersonic vortex: convergence analysis for the test case presented in section 6.2. Numerical results obtained with Dirichlet for the internal boundary and slip-wall for the external boundary without any correction on the conformal linear meshes described in table 8. One can notice that, due to the boundary approximation, the classical slip-wall condition provide convergence trends lower than 2.</caption><thead><tr><th rowspan="2">Grid level</th><th colspan="2">ρ</th><th colspan="2">ρu</th><th colspan="2">ρv</th><th colspan="2">ρE</th></tr><tr><th>L<sub>2</sub></th><th>ñ</th><th>L<sub>2</sub></th><th>ñ</th><th>L<sub>2</sub></th><th>ñ</th><th>L<sub>2</sub></th><th>ñ</th></tr></thead><tbody><tr><td colspan="9" style="text-align:center; font-style:italic;">DG-P<sub>1</sub></td></tr><tr><td>0</td><td>1.0757E-1</td><td>-</td><td>2.1620E-1</td><td>-</td><td>2.1855E-1</td><td>-</td><td>4.9486E-1</td><td>-</td></tr><tr><td>1</td><td>3.9747E-2</td><td>1.44</td><td>8.4226E-2</td><td>1.36</td><td>8.4449E-2</td><td>1.37</td><td>1.8268E-1</td><td>1.44</td></tr><tr><td>2</td><td>1.3584E-2</td><td>1.55</td><td>2.9191E-2</td><td>1.53</td><td>2.9274E-2</td><td>1.53</td><td>6.2762E-2</td><td>1.54</td></tr><tr><td>3</td><td>4.8112E-3</td><td>1.50</td><td>1.0417E-2</td><td>1.49</td><td>1.0448E-2</td><td>1.49</td><td>2.2273E-2</td><td>1.50</td></tr><tr><td colspan="9" style="text-align:center; font-style:italic;">DG-P<sub>2</sub></td></tr><tr><td>0</td><td>1.5596E-1</td><td>-</td><td>3.6487E-1</td><td>-</td><td>3.6579E-1</td><td>-</td><td>7.4800E-1</td><td>-</td></tr><tr><td>1</td><td>6.0207E-2</td><td>1.37</td><td>1.3306E-1</td><td>1.46</td><td>1.3325E-1</td><td>1.46</td><td>2.8177E-1</td><td>1.41</td></tr><tr><td>2</td><td>2.1706E-2</td><td>1.47</td><td>4.7605E-2</td><td>1.48</td><td>4.7647E-2</td><td>1.48</td><td>1.0117E-1</td><td>1.48</td></tr><tr><td>3</td><td>7.8804E-3</td><td>1.46</td><td>1.6924E-2</td><td>1.49</td><td>1.6933E-2</td><td>1.49</td><td>3.6319E-2</td><td>1.48</td></tr><tr><td colspan="9" style="text-align:center; font-style:italic;">DG-P<sub>3</sub></td></tr><tr><td>0</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td></tr><tr><td>1</td><td>1.4987E-1</td><td>-</td><td>1.5568E-1</td><td>-</td><td>1.5817E-1</td><td>-</td><td>5.5100E-1</td><td>-</td></tr><tr><td>2</td><td>6.9222E-2</td><td>1.11</td><td>6.5801E-2</td><td>1.24</td><td>6.6299E-2</td><td>1.25</td><td>2.5165E-1</td><td>1.13</td></tr><tr><td>3</td><td>3.1405E-2</td><td>1.14</td><td>2.8275E-2</td><td>1.22</td><td>2.8553E-2</td><td>1.22</td><td>1.1342E-1</td><td>1.15</td></tr></tbody></table>

both boundaries are discretized with linear meshes the convergence trends cannot be higher than 2. Moreover, it is well-known that the reflecting boundary conditions give rise to severe spurious errors when increasing the polynomial order [49, 13]. For this reason, on the same mesh, increased errors are found when increasing the order of the discretization. Similar trends are also shown in other works in the literature [13, 14]. Notice that due to important spurious disturbances, the simulation performed with $P_3$ on the coarsest mesh does not provide any result, i.e. the simulation accidentally stops before the final time is reached. We validate again the ROD boundary treatment on this experiment, achieving the expected discretization errors and order-of-accuracy (see table 10). Notice that when enforcing slip-wall conditions, the estimated order-of-accuracy is $M$ rather than $M+1$. This is due to the fact that the complex slip-wall treatment needs Neumann and Robin boundary conditions that are imposed using the derivative of the polynomials, hence losing one order of accuracy. This is also typical of other similar approaches [26, 25]. Although the boundary treatment loses one order-of-accuracy, all discretization errors found in our experiments are several orders of magnitude smaller than those obtained with no corrections, which are at most second order. Moreover increasing the polynomial order always allows a reduction in error of at least one order of magnitude, even when the increase of rate of convergence is mild. A trivial way to overcome this issue consists in constructing derivative polynomials of degree $M$ using some informations from neighboring cells. This could be done only for those cells on the Neumann/Robin boundary, as shown in [26]. However, this is out of the scope of this work therefore we will not discuss this case here.

### 6.3. Shock-cylinder interaction

In the last test case, we present the possibility of tackling compressible flows featuring discontinuities with a consistent ROD boundary treatment to improve the accuracy close to internal boundaries. The employed ADER-DG scheme is supplemented with the *a posteriori* sub-cell finite volume limiter discussed in section 3.4. Here we present the numerical results obtained on
<table><caption>Table 10: Isentropic supersonic vortex: convergence analysis for the test case presented in section 6.2. Numerical results obtained with Dirichlet for the internal boundary and slip-wall for the external boundary with the ROD-ADER-DG reconstruction on the conformal linear meshes described in table 8. One can notice that the expected sub-optimal precision given by Neumann conditions is achieved.</caption><thead><tr><th rowspan="2">Grid level</th><th colspan="2">ρ</th><th colspan="2">ρu</th><th colspan="2">ρv</th><th colspan="2">ρE</th></tr><tr><th>L<sub>2</sub></th><th>ñ</th><th>L<sub>2</sub></th><th>ñ</th><th>L<sub>2</sub></th><th>ñ</th><th>L<sub>2</sub></th><th>ñ</th></tr></thead><tbody><tr><td colspan="9">ROD-ADER-DG-P<sub>1</sub></td></tr><tr><td>0</td><td>1.2539E-1</td><td>-</td><td>1.6824E-1</td><td>-</td><td>1.6893E-1</td><td>-</td><td>5.5124E-1</td><td>-</td></tr><tr><td>1</td><td>3.7048E-2</td><td>1.76</td><td>5.4835E-2</td><td>1.62</td><td>5.4921E-2</td><td>1.62</td><td>1.7740E-1</td><td>1.64</td></tr><tr><td>2</td><td>9.5025E-3</td><td>1.96</td><td>1.3583E-2</td><td>2.01</td><td>1.3596E-2</td><td>2.01</td><td>4.5769E-2</td><td>1.95</td></tr><tr><td>3</td><td>2.2951E-3</td><td>2.05</td><td>3.1448E-3</td><td>2.11</td><td>3.1546E-3</td><td>2.10</td><td>1.1120E-2</td><td>2.04</td></tr><tr><td colspan="9">ROD-ADER-DG-P<sub>2</sub></td></tr><tr><td>0</td><td>2.1023E-2</td><td>-</td><td>2.6769E-2</td><td>-</td><td>2.8488E-2</td><td>-</td><td>9.5588E-2</td><td>-</td></tr><tr><td>1</td><td>5.4789E-3</td><td>1.94</td><td>7.1941E-3</td><td>1.90</td><td>7.2213E-3</td><td>1.98</td><td>2.6524E-2</td><td>1.85</td></tr><tr><td>2</td><td>1.2612E-3</td><td>2.12</td><td>1.6497E-3</td><td>2.12</td><td>1.6571E-3</td><td>2.12</td><td>6.1478E-3</td><td>2.11</td></tr><tr><td>3</td><td>2.9615E-4</td><td>2.09</td><td>3.9261E-4</td><td>2.07</td><td>3.9366E-4</td><td>2.07</td><td>1.4588E-3</td><td>2.07</td></tr><tr><td colspan="9">ROD-ADER-DG-P<sub>3</sub></td></tr><tr><td>0</td><td>1.2482E-3</td><td>-</td><td>1.4765E-3</td><td>-</td><td>1.4488E-3</td><td>-</td><td>4.7772E-3</td><td>-</td></tr><tr><td>1</td><td>5.7882E-5</td><td>4.43</td><td>8.1071E-5</td><td>4.18</td><td>8.0220E-5</td><td>4.17</td><td>2.3357E-4</td><td>4.35</td></tr><tr><td>2</td><td>6.6352E-6</td><td>3.12</td><td>6.3173E-6</td><td>3.68</td><td>6.2445E-6</td><td>3.68</td><td>2.3344E-5</td><td>3.32</td></tr><tr><td>3</td><td>1.1183E-6</td><td>2.57</td><td>1.1242E-6</td><td>2.49</td><td>1.1643E-6</td><td>2.42</td><td>4.4317E-6</td><td>2.40</td></tr></tbody></table>

<figure><img src="image_8.png" /><figcaption>Figure 7: Shock-cylinder interaction: characteristics of computational meshes on a computational domain [-2, 6] × [-3, 3] for the test case presented in section 6.3.</figcaption></figure>
a computational domain $[-2, 6] \times [-3, 3]$ discretized with a unstructured conformal triangulation and a fully embedded one (see fig. 7). The cylinder is centered in $(0, 0)$ and has radius 0.5. The initial condition consists in a shock wave traveling at Mach number $M_s = 1.3$ and is then given via the Rankine-Hugoniot conditions. The flow upstream the shock is at rest and is characterized by density and pressure, respectively being $\rho = 1.4$ and $p = 1$.

In fig. 8 we present the numerical results obtained with third order polynomials on the conformal triangulation by imposing naive reflecting boundary conditions (without ROD) and the general slip conditions (with ROD). It is clear that when curved boundaries are considered, a more general boundary treatment that takes into account the real geometry is essential to capture properly the flow dynamics. In the simulation performed without ROD, we notice that spurious disturbances starts from the poor boundary treatment and propagates in the computational domain making the solution less accurate. On the contrary, much smoother iso-contours both close and away from the boundary are observed in the simulation performed with the ROD method.

Finally, in fig. 9 a numerical simulation run on a fully embedded triangulation is presented. Here we show the results obtained by imposing the general slip condition by means of the ROD boundary reconstruction on a surrogate boundary completely detached from the exact geometry as shown in fig. 7. In general the density iso-contours appear globally very similar to the one obtained in the conformal case performed with the ROD method. Some spurious disturbances present in some cells of the discretized boundary are normal due to the use of a less accurate boundary treatment when the *a posteriori* limiter is activated. For this reason, we present side by side also the troubled cells recomputed by means of the low order method to show that disturbances only arise in those regions, while in the non-troubled parts where the ROD treatment is applied the solution appears smooth.

## 7. Conclusions and future investigations

In this work, we have investigated the numerical aspects of the Reconstruction for Off-site Data (ROD) method in a discontinuous Galerkin and space-time ADER framework. The ROD method is formulated as a constraint minimization problem that allows to recover consistent boundary polynomials used for the correct imposition of boundary conditions. Originally developed in a finite volume framework, the ROD idea consisted in managing multiple constraints related to scattered mean values associated with neighboring elements, found with specific algorithms. Thanks to the local discontinuous finite element formulation, we show how the ROD method can be easily adapted by evaluating the boundary cell polynomial within quadrature points, rather than using neighboring elements. By doing so, a classical constrained least squares should be solved with Lagrange multipliers for the boundary integral of the semi-discretization (6). Naturally, it follows that by integrating (6) with classical high order time integration schemes, like one-step Runge-Kutta methods, such minimization problems must be solved at each sub-time node. Although the linear system arising from the constrained least squares is local, an inversion is needed at each boundary cell for each sub-time step, hence impacting the computational costs. For instance, when dealing with $P_2$ polynomials a local matrix of $9 \times 9$ must be inverted each time; for $P_3$ the matrix size becomes $14 \times 14$. In this paper, we also explored a new advantageous perspective for the ROD method when considering the space-time ADER framework. Indeed, we showed that by using the same space-time polynomial coefficients of the ADER predictor (11) for the minimization problem, we obtain a simplified linear system that provides a polynomial consistent with the physical boundary condition for the entire space-time element, hence only one inversion is needed. Several numerical experiments are shown to assess the high order convergence properties
<figure><img src="image_9.png" /><figcaption>Figure 8: Shock-cylinder interaction: density iso-contours at different times for the test case presented in section 6.3. We show the results obtained with classical reflecting boundary conditions (right), and the wall conditions enforced by the ROD boundary polynomial (left) on the conformal mesh presented in fig. 7a.</figcaption></figure>
<figure><img src="image_10.png" /><figcaption>Figure 9: Shock-cylinder interaction: density iso-contours and troubled cells at different times for the test case presented in section 6.3. We show the results obtained with the ROD boundary polynomial on the fully embedded mesh presented in fig. 7b.</figcaption></figure>
of both the spatial and space-time ROD reconstructions when dealing with Dirichlet conditions on a curved boundary using conformal linear meshes and fully embedded triangulations. Application of the novel space-time formulation to general Dirichlet-Neumann-Robin boundary conditions is also presented for the general slip-wall treatment. The latter opens up to many applications when the interaction between a flow and a boundary is crucial. In the last test case, the application of the ROD method to a shock-cylinder interaction is presented, where the discontinuous features of the flow are properly captured using an *a posteriori* MOOD approach.

Several important aspects related to the proposed numerical method need further investigations. In particular, we believe that the introduction of additional constraints within the optimization problem may bring more accurate results and improve the conditioning of the problem and the activation of the shock limiter on boundary cells. For instance, mass conservation, positivity and entropy preservation will be considered in future works.

## Acknowledgments

M.R. is a member of the CARDAMOM team, INRIA and University of Bordeaux research center. E.G. gratefully acknowledges the support and funding received from the European Union with the ERC Starting Grant ALcHyMiA (No. 101114995). Views and opinions expressed are however those of the author(s) only and do not necessarily reflect those of the European Union or the European Research Council Executive Agency. Neither the European Union nor the granting authority can be held responsible for them.

## References

[1] Z. Wang, K. Fidkowski, R. Abgrall, F. Bassi, D. Caraeni, A. Cary, H. Deconinck, R. Hartmann, K. Hillewaert, H. Huynh, N. Kroll, G. May, P.-O. Persson, B. van Leer, M. Visbal, High-order cfd methods: current status and perspective, International Journal for Numerical Methods in Fluids 72 (8) (2013) 811–845.

[2] O. Zienkiewicz, P. Morice, The finite element method in engineering, McGraw-Hill, 1971.

[3] T. Hughes, J. Cottrell, Y. Bazilevs, Isogeometric analysis: Cad, finite elements, nurbs, exact geometry, and mesh refinement, Computer Methods in Applied Mechanics and Engineering 194 (2005) 33–40.

[4] X.-J. Luo, M. S. Shephard, R. M. O’bara, R. Nastasia, M. W. Beall, Automatic p-version mesh generation for curved domains, Engineering with Computers 20 (3) (2004) 273–285.

[5] O. Sahni, X. Luo, K. Jansen, M. Shephard, Curved boundary layer meshing for adaptive viscous flow simulations, Finite Elements in Analysis and Design 46 (1-2) (2010) 132–139.

[6] A. Loseille, L. Rochery, Developments on the $P^2$ cavity operator and Bézier Jacobian correction using the simplex algorithm. arXiv:https://arc.aiaa.org/doi/pdf/10.2514/6.2022-0389, doi:10.2514/6.2022-0389.
URL https://arc.aiaa.org/doi/abs/10.2514/6.2022-0389

[7] M. Couplet, M. Reberol, J.-F. Remacle, Generation of High-Order Coarse Quad Meshes on CAD Models via Integer Linear Programming. arXiv:https://arc.aiaa.org/doi/pdf/10.2514/6.2021-2991, doi:10.2514/6.2021-2991.
URL https://arc.aiaa.org/doi/abs/10.2514/6.2021-2991
[8] S. Dey, R. M. O'Bara, M. S. Shephard, Towards curvilinear meshing in 3d: the case of quadratic simplices, Computer-Aided Design 33 (3) (2001) 199-209.

[9] M. Fortunato, P.-O. Persson, High-order unstructured curved mesh generation using the winslow equations, Journal of Computational Physics 307 (2016) 1-14.

[10] D. Moxey, D. Ekelschot, Ü. Keskin, S. Sherwin, J. Peiró, High-order curvilinear meshing using a thermo-elastic analogy, Computer-Aided Design 72 (2016) 130-139.

[11] A. Veilleux, G. Puigt, H. Deniau, G. Daviller, A stable spectral difference approach for computations with triangular and hybrid grids up to the 6th order of accuracy, Journal of Computational Physics 449 (2022) 110774.

[12] Z. Wang, Y. Sun, Curvature-based wall boundary condition for the euler equations on unstructured grids, AIAA journal 41 (1) (2003).

[13] L. Krivodonova, M. Berger, High-order accurate implementation of solid wall boundary conditions in curved geometries, Journal of computational physics 211 (2) (2006) 492-512.

[14] M. Ciallella, E. Gaburro, M. Lorini, M. Ricchiuto, Shifted boundary polynomial corrections for compressible flows: high order on curved domains using linear meshes, Applied Mathematics and Computation 441 (2023) 127698.

[15] R. Costa, S. Clain, R. Loubère, G. J. Machado, Very high-order accurate finite volume scheme on curved boundaries for the two-dimensional steady-state convection-diffusion equation with dirichlet condition, Applied Mathematical Modelling 54 (2018) 752-767.

[16] A. Main, G. Scovazzi, The shifted boundary method for embedded domain computations. part I: Poisson and Stokes problems, Journal of Computational Physics 372 (2018) 972-995. doi:10.1016/j.jcp.2017.10.026.

[17] A. Main, G. Scovazzi, The shifted boundary method for embedded domain computations. part II: Linear advection-diffusion and incompressible Navier-Stokes equations, Journal of Computational Physics 372 (2018) 996-1026. doi:10.1016/j.jcp.2018.01.023.

[18] T. Song, A. Main, G. Scovazzi, M. Ricchiuto, The shifted boundary method for hyperbolic systems: Embedded domain computations of linear waves and shallow water flows, Journal of Computational Physics 369 (2018) 45-79. doi:10.1016/j.jcp.2018.04.052.

[19] K. Li, N. M. Atallah, G. A. Main, G. Scovazzi, The shifted interface method: A flexible approach to embedded interface computations, International Journal for Numerical Methods in Engineering 121 (3) (2020) 492-518. doi:10.1002/nme.6231.

[20] N. M. Atallah, C. Canuto, G. Scovazzi, The high-order shifted boundary method and its analysis, Computer Methods in Applied Mechanics and Engineering 394 (2022) 114885.

[21] T. Carlier, L. Nouveau, H. Beaugendre, M. Colin, M. Ricchiuto, An enriched shifted boundary method to account for moving fronts, Journal of Computational Physics (2023) 112295.

[22] L. Nouveau, M. Ricchiuto, G. Scovazzi, High-order gradients with the shifted boundary method: An embedded enriched mixed formulation for elliptic pdes, Journal of Computational Physics 398 (2019) 108898. doi:10.1016/j.jcp.2019.108898.
[23] M. Ciallella, M. Ricchiuto, R. Paciorri, A. Bonfiglioli, Extrapolated discontinuity tracking for complex 2d shock interactions, Computer Methods in Applied Mechanics and Engineering 391 (2022) 114543.

[24] A. Assonitis, R. Paciorri, M. Ciallella, M. Ricchiuto, A. Bonfiglioli, Extrapolated shock fitting for two-dimensional flows on structured grids, AIAA Journal 60 (11) (2022) 6301–6312.

[25] J. Visbech, A. P. Engsig-Karup, M. Ricchiuto, A spectral element solution of the poisson equation with shifted boundary polynomial corrections: influence of the surrogate to true boundary mapping and an asymptotically preserving robin formulation, arXiv preprint arXiv:2310.17621 (2023).

[26] R. Costa, J. M. Nóbrega, S. Clain, G. J. Machado, R. Loubère, Very high-order accurate finite volume scheme for the convection-diffusion equation with general boundary conditions on arbitrary curved boundaries, International Journal for Numerical Methods in Engineering 117 (2) (2019) 188–220.

[27] J. Fernández-Fidalgo, S. Clain, L. Ramírez, I. Colominas, X. Nogueira, Very high-order method on immersed curved domains for finite difference schemes with regular cartesian grids, Computer Methods in Applied Mechanics and Engineering 360 (2020) 112782.

[28] S. Clain, D. Lopes, R. M. Pereira, Very high-order cartesian-grid finite difference method on arbitrary geometries, Journal of Computational Physics 434 (2021) 110217.

[29] R. Costa, J. M. Nóbrega, S. Clain, G. J. Machado, Efficient very high-order accurate polyhedral mesh finite volume scheme for 3d conjugate heat transfer problems in curved domains, Journal of Computational Physics 445 (2021) 110604.

[30] R. Costa, S. Clain, G. J. Machado, J. M. Nóbrega, Very high-order accurate finite volume scheme for the steady-state incompressible navier–stokes equations with polygonal meshes on arbitrary curved boundaries, Computer Methods in Applied Mechanics and Engineering 396 (2022) 115064.

[31] R. Costa, S. Clain, G. J. Machado, J. M. Nóbrega, H. B. da Veiga, F. Crispo, Imposing slip conditions on curved boundaries for 3d incompressible flows with a very high-order accurate finite volume scheme on polygonal meshes, Computer Methods in Applied Mechanics and Engineering 415 (2023) 116274.

[32] B. Cockburn, C.-W. Shu, The runge-kutta local projection-discontinuous-galerkin finite element method for scalar conservation laws, ESAIM: Mathematical Modelling and Numerical Analysis 25 (3) (1991) 337–361.

[33] B. Cockburn, C.-W. Shu, The runge-kutta discontinuous Galerkin method for conservation laws v: multidimensional systems, Journal of Computational Physics 141 (2) (1998) 199–224.

[34] B. Cockburn, C.-W. Shu, Runge-kutta discontinuous galerkin methods for convection-dominated problems, Journal of scientific computing 16 (2001) 173–261.

[35] B. Cockburn, G. E. Karniadakis, C.-W. Shu, Discontinuous Galerkin methods: theory, computation and applications, Vol. 11, Springer Science & Business Media, 2012.
[36] A. Mazaheri, C.-W. Shu, V. Perrier, Bounded and compact weighted essentially nonoscillatory limiters for discontinuous Galerkin schemes: Triangular elements, Journal of Computational Physics 395 (2019) 461-488.

[37] M. Dumbser, D. S. Balsara, E. F. Toro, C.-D. Munz, A unified framework for the construction of one-step finite volume and discontinuous Galerkin schemes on unstructured meshes, Journal of Computational Physics 227 (18) (2008) 8209-8253.

[38] W. Boscheri, M. Dumbser, D. Balsara, High-order ader-weno ale schemes on unstructured triangular meshes—application of several node solvers to hydrodynamics and magnetohydrodynamics, International Journal for Numerical Methods in Fluids 76 (10) (2014) 737-778.

[39] W. Boscheri, R. Loubère, High order accurate direct Arbitrary-Lagrangian-Eulerian ADER-MOOD finite volume schemes for non-conservative hyperbolic systems with stiff source terms, Communications in Computational Physics 21 (1) (2017) 271-312.

[40] W. Boscheri, R. Loubere, M. Dumbser, Direct Arbitrary-Lagrangian-Eulerian ADER-MOOD finite volume schemes for multidimensional hyperbolic conservation laws, Journal of Computational Physics 292 (2015) 56-87.

[41] S. Busto, S. Chiocchetti, M. Dumbser, E. Gaburro, I. Peshkov, High order ADER schemes for continuum mechanics, Frontiers in Physics 8 (2020) 32.

[42] F. Fambri, M. Dumbser, O. Zanotti, Space-time adaptive ader-dg schemes for dissipative flows: Compressible navier-stokes and resistive mhd equations, Computer Physics Communications 220 (2017) 297-318.

[43] D. S. Balsara, T. Rumpf, M. Dumbser, C.-D. Munz, Efficient, high accuracy ader-weno schemes for hydrodynamics and divergence-free magnetohydrodynamics, Journal of Computational Physics 228 (7) (2009) 2480-2516.

[44] E. Gaburro, W. Boscheri, S. Chiocchetti, C. Klingenberg, V. Springel, M. Dumbser, High order direct Arbitrary-Lagrangian-Eulerian schemes on moving voronoi meshes with topology changes, Journal of Computational Physics 407 (2020) 109167.

[45] S. Clain, S. Diot, R. Loubère, A high-order finite volume method for systems of conservation laws—multi-dimensional optimal order detection (mood), Journal of computational Physics 230 (10) (2011) 4028-4050.

[46] S. Diot, R. Loubère, S. Clain, The multidimensional optimal order detection method in the three-dimensional case: very high-order finite volume method for hyperbolic systems, International Journal for Numerical Methods in Fluids 73 (4) (2013) 362-392.

[47] R. Loubere, M. Dumbser, S. Diot, A new family of high order unstructured MOOD and ADER finite volume schemes for multidimensional systems of hyperbolic conservation laws, Communications in Computational Physics 16 (3) (2014) 718-763.

[48] E. Gaburro, M. Dumbser, A posteriori subcell finite volume limiter for general PNPM schemes: Applications from gasdynamics to relativistic magnetohydrodynamics, Journal of Scientific Computing 86 (3) (2021) 1-41.
[49] F. Bassi, S. Rebay, High-order accurate discontinuous finite element solution of the 2d euler equations, Journal of computational physics 138 (2) (1997) 251-285.

[50] S. Diot, S. Clain, R. Loubère, Improved detection criteria for the multi-dimensional optimal order detection (mood) on unstructured meshes with very high-order polynomials, Computers & Fluids 64 (2012) 43-63.

[51] N. Chalmers, L. Krivodonova, A robust cfl condition for the discontinuous Galerkin method on triangular meshes, Journal of Computational Physics 403 (2020) 109095. doi:https://doi.org/10.1016/j.jcp.2019.109095.
URL https://www.sciencedirect.com/science/article/pii/S0021999119308009

[52] A. Chertock, A. Coco, A. Kurganov, G. Russo, A second-order finite-difference method for compressible fluids in domains with moving boundaries, Commun Comput Phys 23 (2018) 230-263.