A comparison of interpolation techniques for non-conformal high-order discontinuous Galerkin methods

The capability to incorporate moving geometric features within models for complex simulations is a common requirement in many fields. The fluid mechanics within aeronautical applications, for example, routinely feature rotating (e.g. turbines, wheels and fan blades) or sliding components (e.g. in compressor or turbine cascade simulations). With an increasing trend towards the high-fidelity modelling of these cases, in particular combined with the use of high-order discontinuous Galerkin methods, there is therefore a requirement to understand how different numerical treatments of the interfaces between the static mesh and the sliding/rotating part impact on overall solution quality. In this article, we compare two different approaches to handle this non-conformal interface. The first is the so-called mortar approach, where flux integrals along edges are split according to the positioning of the non-conformal grid. The second is a lesser-documented point-to-point interpolation method, where the interior and exterior quantities for flux evaluations are interpolated from elements lying on the opposing side of the interface. Although the mortar approach has advantages in terms of its numerical properties, in that it preserves the local conservation properties of DG methods, in the context of complex 3D meshes it poses significant implementation difficulties which the point-to-point method handles more readily. In this article we examine the numerical properties of each method, focusing not only on observing convergence orders for smooth solutions, but also how each method performs in under-resolved simulations of linear and nonlinear hyperbolic problems, to inform the use of these methods in implicit large-eddy simulations.


Introduction
Problems containing features that move or deform are found in many research areas, but are particularly prevalent in the study of various fluid dynamics phenomena. In particular, aeronautical applications commonly feature rotating or sliding geometries, with typical examples in this area including turbomachinery [1,2], unmanned aerial vehicles [3,4], insect and avian flight aerodynamics [5,6], and HVAC (heating, ventilating, and air conditioning) [7,8]. Being able to accurately model these moving geometries and their subsequent impact on the underlying flow physics is highly important: for example, predicting how the profiles of turbine or compressor blades impact on propulsion efficiency, or how wing profiles affect the performance of wind turbines. Moreover, the cost and difficulty of performing full-scale experimental testing of such geometries can be challenging from the perspective of both instrumentation and expense. For these reasons, computational fluid dynamics (CFD) is now commonplace in the design and modelling process. If a CFD method is to be regarded as universally useful in these application areas, then it must be capable of accurately modelling moving geometry and ideally provide high-fidelity results beyond the scope of physical field tests alone.
Most leading software for CFD is based around lower-order finite volume or finite element methods, typically leveraging the computationally-cheap Reynolds Averaged Navier-Stokes (RANS) equations in combination with a turbulence closure model. However this approach has natural limitations in studying the aforementioned problems at very high levels of fidelity [9]. With the large increases in computational power in recent years, a more recent trend is to instead consider transient simulations that leverage implicit large-eddy simulation (ILES) or under-resolved direct numerical simulation (uDNS) [10,11]. This approach is more computationally expensive than RANS, but also provides greater accuracy and enables high-fidelity simulations of the complex geometries that lie in this regime [12]. The combination of LES with less common high-order methods, either based on continuous or discontinuous Galerkin (CG/DG) methods, has seen significant interest in recent years, particularly in aeronautics applications [13]. From a numerical perspective, high-order methods possess far lower levels of numerical diffusion and dispersion, making them ideally suited to resolving features across long time-and length-scales. This can overcome a significant bottleneck when considering these simulations at lower orders, since very fine grid resolutions are required to overcome the effects of numerical errors [14]. Additionally, from a computational perspective, the larger number of floating-point operations that are required per degree-of-freedom as the polynomial order is increased means that high-order methods can be used to overcome the memory bandwidth bottlenecks that are common in modern computational hardware [15]. The combination of these effects means that high-order methods can achieve higher accuracy per degree-of-freedom at equivalent or lower computational cost to lower-order methods.
However, these methods are somewhat less well-explored in the simulation of problems involving rotating or sliding geometries, which require the treatment  of non-conformal interfaces between elements. In this article, we explore two common approaches to the handling of non-conformal interfaces and compare their numerical performance in a range of linear and nonlinear problems.

Requirements for moving geometry simulations
One approach to tackling the problem of moving geometry is the sliding mesh method, where the mesh is separated into two or more separate regions, and during the simulation the regions will slide relative to one another. This provides a way to prescribe simple mesh motion via rotation or translation. The most basic problem case is to employ a stationary outer region, with a rotating circular region within it, which is found in many applications, for example modelling flow in a stirred tank [16]. An exaggerated example of a the sliding mesh method with a non-conforming interface is shown in Fig. 1. It is clear that this process results in a non-conformal mesh: i.e. a mesh where elements do not connect to precisely one other element through one of their edges or faces. Most CFD simulations make use of conformal meshes, where each edge (in 2D) or face (in 3D) of an element has precisely one neighbouring element. As such, techniques need to be developed in order to accurately preserve solution quality across the non-conformal interface.
In the 'classical' spectral element method, where C 0 continuity is imposed between elements in the continuous Galerkin formulation, three main techniques have been evaluated for use in non-conformal meshes. Possibly the most wellknown of these arises when performing h-adaptation in an octree-like manner, so that 2-to-1 element subdivisions are obtained in the resulting mesh. In this case, hanging nodes are generated and their values can be constrained through the analytic definition of the basis functions lying along an edge or face, together with the assembly mapping that is used to construct mass and stiffness matrices [17,18]. However, in the sliding mesh case, elements may overlap at arbitrary positions along their edges and faces, making this approach infeasible. Possibly the most widely-adopted approach to implementing generic non-conformal interfaces in the CG setting is the mortar technique [19]. In this setting, one augments the traditional C 0 function spaces for each conformal domain with functions defined on mortar elements at the interface between two domains. The weak form of the problem is then augmented to incorporate a penalty for the jump across the interface in an appropriate manner, so that the convergence order of the scheme is retained. This approach is visualised in Fig. 2, where we note that the mortar elements are constructed at the common intersection points of each element.
An alternative approach to imposing non-conformal conditions, and which is perhaps less commonly-used, is to instead adopt a point-to-point interpolation across the elemental interface. In this setting, no attempt to construct mortar elements is made and the function space is defined in the usual manner for each conformal domain. However, when values within elements are desired at the left-hand side of the interface, they are obtained by performing a polynomial interpolation from the values on the right-hand side, and vice versa. This approach was first implemented and tested for geophysical problems in [20] in the CG setting, where it was shown to demonstrate convergence-order preserving properties. A sample visualisation of this approach is shown in Fig. 3, where dotted arrows denote the evaluation of the high-order polynomial defined by points on the edge of element Ω A to obtain their values within the boundaries of elements in ∂Ω B and ∂Ω C . This interpolation process can be built into the assembly operation that is used to construct mass and stiffness matrices.

Non-conformal techniques for the discontinuous Galerkin method
At present, there is a significant interest in the development of high-order fluid dynamics solvers for iLES/uDNS based around the DG method due to its favourable stability properties in this regime [21,22,23]. In the context of sliding mesh simulations, DG also offers an easier route to the accurate treatment of non-conformal interfaces between elements across the sliding interface, since elements are naturally disconnected as part of the formulation of the method. Additionally, approaches to impose non-conformal interfaces in DG are perhaps less well-explored than in CG. In the DG formulation, connectivity between elements is imposed through a flux term, which may be either an upwind-type solver for linear problems or a more complex Riemann problem for more general nonlinear hyperbolic systems. These fluxes are computed on integrals across each edge of an element and take the form where Γ e is an edge of element Ω e , u + and u − is a vector of conserved variables on the exterior and interior of the element respectively,f is the numerical flux function and n is an outwards-facing normal. The question then is how one computes these integrals, given that the exterior values u + may now lie across more than one element on the other side of the interface.
The mortaring approach has been investigated in a number of works from the DG perspective. First, we note that unlike the CG setting which requires modifications to the function space and weak form of the problem, in DG by 'mortaring' we only refer to the act of constructing mortar elements on which to compute the flux integral. That is, the integral above is split into multiple integrals, one for each mortar element. This approach was first investigated by Kopriva et al. in the study of both fluid dynamics [24] and electromagnetics problems [25]. The same approach has been used fairly extensively for problems involving sliding meshes; for example by [26] in a hybrid DG-Fourier pseudospectral solver for the incompressible Navier-Stokes equations, in the construction of a spectral difference solver for the compressible Navier-Stokes equations that incorporates sliding grids by Zhang and Liang [27] and more recently in the hyperbolic solver FLEXI [28]. This approach has the significant advantage that it preserves the local conservation property of the DG method, which is important from the perspective of obtaining accurate results that conserve mass (in the case of CFD). However, although mortaring is straightforward in two dimensions, a significant challenge in the use of the mortar approach for general three-dimensional problems is the generation of the mortar elements themselves. In general this could involve the re-meshing of the non-conformal interface between domains at each timestep in order to generate an appropriate mortar space, as adopted by Aguerre et al. [29] based on the supermesh construction of Farrell et al. [30]. We note that these works consider linear problems -this approach would therefore be even more complex in the context of high-order methods where the interface is curved.
The alternative approach is therefore to consider the point-to-point interpolation method in the DG context, since implementation is relatively straightforward by comparison as it does not require the construction of mortar elements. However, neither the implementation, performance or robustness of this approach for DG has been thoroughly investigated in the literature to date. In particular, potential issues may arise from the discontinuity of fluxes between elements: for smooth solutions and at high polynomial orders, the interpolation between neighbouring non-conformal elements will likely introduce very little error into the resulting solution. However, in the presence of under-resolved simulations, which are more prone to admitting discontinuities in flow solution between elements, the discontinuity may introduce additional numerical error that warrants further study. A prototypical example which demonstrates this in an illustrative manner is shown in Fig, 4. On the left side of the interface, the two discontinuous solutions from elements Ω A and Ω B must be sampled at integration points on the skeleton of Ω C . If the two functions are sufficiently discontinuous, the interpolation procedure could result in spurious noise introduced into the interior of Ω C .

Aim of this work
To date, the point-to-point method has not been well-studied in the literature. A study by Kopera and Giraldo [31] is one of the very few references, to the best of the authors' knowledge, that consider the point-to-point interpolation approach in DG, where CG and DG implementations of the interpolation technique are examined and their mass conservation properties are reported. However we note that in this case, only hanging-node type vs. more generic non-conformal interfaces are considered. Additionally, this work was performed in well-resolved cases which may not be the case for more general iLES/uDNStype problems. In this article, we therefore aim to address this gap in the understanding of the performance of these approaches by performing a comparative study of the mortar and point-to-point techniques. We consider several aspects, including a validation of convergence order for both approaches, the performance of each method in terms of numerical diffusion for a linear transport equation, at varying degrees of underresolution, and the behaviour of each method when considering the nonlinear problem of the compressible Euler equations across long time periods.
The remainder of the paper is structured as follows. In section 2 we outline the theoretical framework of the two formulations and outline our implementation strategy within the spectral/hp element framework Nektar++ [32,33]. Section 3 presents the results of our studies from a linear transport equation and the nonlinear compressible Euler equations. Finally, in section 4, we draw some brief conclusions and discuss the key performance characteristics of each method.

The DG formulation of the spectral/hp element method
In this section, we briefly introduce the discontinuous Galerkin (DG) discretisation of the spectral/hp element method. A more thorough overview can be found in several other works, e.g. [33,34]. The starting point for the DG formulation is the same as any other typical finite element problem, in that we consider a domain Ω comprised of non-overlapping elements Ω e such that Ω = e Ω e . Given a general hyperbolic conservation law for conserved variables u taking the form ∂u ∂t we follow the standard Galerkin approach and, on a single element, construct the weak form via multiplication by a test function v and integrating by parts where (u, v) Ω e = Ω e uv dx and u, v ∂Ω e = ∂Ω e uv ds denote inner products on the volume and surface, respectively. Moreover,f defines a numericallycalculated flux term which, as explained in the previous section, may take the form of a general Riemann problem, and which depends on the element-exterior and interior velocities u + and u − , respectively. Within each element, we represent u using an expansion of high-order polynomials, so that In this expression, we note that the approximation is defined with the use of a standard (reference) element Ω st , with φ n denoting an appropriate set of basis functions. An isoparametric mapping χ e : Ω st → Ω e defines a possibly curvilinear element Ω e , so that x = χ e (ξ) for ξ ∈ Ω st . We additional equip the standard element with a distribution of quadrature points ξ q and weight w q , so that upon selecting test functions v = φ n we then evaluate the terms in eq. (2) as finite summations, i.e.
In this study, we consider only two-dimensional elements and select tensor products of Gauss-Lobatto-Legendre points to evaluate quadrature. As basis functions we adopt the hierarchical modified basis of Karniadakis & Sherwin [34].
Similarly to the classical Lagrange basis, these basis functions have the beneficial property of boundary-interior decomposition, which makes the addition of flux terms into the overall elemental degrees of freedom a straightforward addition operation. In particular we note that the flux integral terms can be considered along each edge i of Ω e , which we denote by Γ e i , as the integral where now ψ n denotes a basis function with support along edge i. In particular, we note that the solution variables along y ∈ Γ e i can be written as a polynomial expansion u(y) = nû n ψ n (y) As alluded to in the introduction, the central focus of this work is to understand how different evaluations of the flux term in the presence of a non-conformal mesh influence the overall properties and stability of the DG method. In the following sections, we outline the formulation of both the point-to-point interpolation method and the mortar method.

The point-to-point interpolation method
In the point-to-point interpolation method, the interface is handled using a direct interpolation from one side to the other. That is, when we require the values of exterior conserved variables u + at a spatial position x ∈ Γ e i , we adopt the following procedure: • determine a corresponding element Ω f that contains the point x along an edge Γ f j ; • perform a polynomial interpolation at that position using eq.(3) in order to determine u + .
This interpolation is performed for every integration point along Γ e i , as shown diagrammatically in Fig. 3. Once the trace space (i.e. the collection of all edges in the interface) has been fully populated by interpolation, the DG solver can continue as usual with a Riemann solver to calculate the correct numerical flux to then be added into elemental coefficient spaces.
In order to determine a corresponding element that contains the point y, we require the ability to determine the distance of a desired point from any given edge Γ f j . For edges that are straight-sided, this translates into a simple geometric problem which may be solved analytically. However, for curvilinear elements, we must instead utilise the parametric mapping x = ξ e (ξ). In particular, for each edge in the non-conformal interface, we minimise an objective function d(ξ; y) = x − y 2 2 = χ e (ξ) − y 2 2 , i.e. the square of the Euclidean norm · 2 between a point ξ within the edge and the target point y. This then allows us to determine the corresponding reference space point ξ that aligns with the target point x which has minimum distance to y. In our implementation, this is solved via a gradient-descent method utilising a quasi-Newton search direction and backtracking line search, but other Newton-type methods will provide similar convergence properties. Since this is additionally an expensive operation to be performed for every edge within the interface, we make use of an r-tree structure to reduce the initial search space. The octants that are used to construct the r-tree are defined as the bounding box for each curvilinear edge. In this manner, the r-tree can first be interrogated to determine a subset of possible edges under which to then perform the nonlinear optimisation of distance, which further reduces computational cost.
Finally, we require the evaluation of each polynomial expansion (3) at an arbitrary point within the reference element. Although this can be computed directly from eq.(3), this would require the evaluation of each basis function itself at an arbitrary point in the reference element. Instead, therefore, we leverage the use of the representation of the function at quadrature points and rewrite (3) as a summation in terms of Lagrange interpolants, so that  Classically one would then generate an interpolation matrix I as outlined in [34], and perform a dot product to evaluate the desired quantity. However, our timings demonstrate that the use of fast summation based on barycentric interpolation techniques described in [35] yield far better performance for this operation.

The mortar method
The second approach we will consider in this paper is the mortar method which maintains the local conservation properties of DG by constructing mortar elements as visualised in fig. 2. This method applied to the spectral element method was originally developed by Maday et al. [36], and has been used for both incompressible flow [19] and compressible flow problems [24]. We note again that 'mortar' in this sense refers to the act of construction of mortar elements so that flux integrals may be expressed as where M is the number of mortars on edge i, and Ξ m denotes each mortar element. We then construct a polynomial expansion on each mortar element of the same polynomial order. The mortar method is realised by projecting variables from across the interface onto its corresponding mortar element, solving the Riemann problem on the mortars, and then performing an L 2 projection in order to consolidate the contributions from each mortar element. The number of mortars connected to a single interface edge and their relative size is arbitrary, allowing for a wide range of varying mesh circumstances. To give a more concrete definition of the method, we utilise the notation prevalent in Zhang and Liang [27] and Kopriva et al. [25], labelling the two contributing interface segments 'L' and 'R' as shown in Fig. 5.
First we recall that each edge in the interface Γ e i may be represented on a standard segment −1 ≤ ξ ≤ 1 and then mapped using the isoparametric mapping χ e . Similarly, each mortar element has a similar mapping −1 ≤ z ≤ 1 and, in particular, we may write the relationship between the two as where o is the offset of the centre of the mortar relative to the centre of the interface edge, and s is the relative scale factor. The solution on an interface edge can be represented by eq. (3), so that where we consider now only a single scalar quantity u for clarity. We can similarly define the solution on the mortar element, Ξ, as To project the solutions from the element onto the mortar we minimise the norm in the L 2 sense, i.e.
When evaluated at all quadrature points, this can be expressed in matrix form asû where M is the standard elemental mass matrix, and S Ω→Ξ are constructed as To apply the mortar method to the DG formulation, we therefore adopt the following approach: • Construct both the left and right solutions u + and u − onto the mortar using the projection matrices P L→Ξ and P R→Ξ as shown in Fig. 5a.
• Once the solutions are on the mortar, the Riemann solver can be used to compute the numerical fluxf .
• Projecting from the mortars back onto the interface element requires minimising the trace quantities norm in the L 2 sense. For M mortars to the interface element Ω this is as follows In matrix form the solution to this iŝ where S Ξ→Ω is the transpose of S Ω→Ξ taking care to include the respective scale factors.
It is also worth noting that for where the geometry of the interface element is identical to the mortar, for example between Ω A and its corresponding mortar in Fig. 2, the projection matrix is merely the identity matrix soû Ξ =û Ω and f Ω =f Ξ . This can be used to reduce computational costs in this case.

Results
In this section, we report on the results of a number of tests taken using both linear and nonlinear problems to evaluate the efficacy of both the mortar and point-to-point interpolation technique. At each stage, we use conformal grids of similar resolutions to provide a benchmark against which to compare. Each method has been implemented within the Nektar++ spectral/hp element framework [33,32]. In all cases, we consider only explicit timestepping methods with the use of a standard 4th-order Runge-Kutta time integration scheme. The timestep used for each case is reported separately.

Convergence order
In this first case, we test the correctness of our implementation by performing a standard hp-convergence study. For this, we select a standard linear transport equation within a domain Ω = [−5, 5] 2 , so that in eq. (1), F (u) = vu for a constant velocity v = (1, 0). We select an initial condition that is non-polynomial, so that u(x, 0) = sin(πx) sin(πy) together with periodic boundary conditions on all edges, so that the solution propagates indefinitely. Regular grids are constructed using between 81 and 13, 452 quadrilateral elements in the conformal case. The non-conformal case incorporates two interfaces to ensure that the periodic boundaries are conformal to one another for ease of implementation. This results in three domains, of which we shift the central one vertically by half a cell length to create the non-conformal grid similar to Fig. 8b, meaning that the non-conformal mesh will have slightly higher element numbers than its equivalent conformal counterpart. Polynomial orders of P = 3 through P = 11 are considered for each grid, and we select Q = P + 2 quadrature points in each coordinate direction so as to exactly integrate the mass matrix and remove any spurious aliasing error due to the use of numerical integration. We select a timestep size of ∆t = 0.001 and measure the error after one tenth of a cycle (i.e. t = 1) so that error due to timestepping is reduced. Fig. 6 highlights the convergence properties in the L 2 sense of the two nonconformal methods, together with the conformal interface. To ensure clarity the results have been trimmed to remove subsequent points from each polynomial   order after the minimum L 2 error of ∼ 10 −7 has been reached. The results of this study show that, for smooth solutions, all three methods exhibit the expected convergence orders. We note that although there is a slight shifting in global degrees of freedom for the two non-conformal interface handling methods, this is solely due to an increase in elements present in the non-conformal cases (owing to the increased resolution used to offset the grids) and this does not affect the overall trend in convergence properties. These results validate both that the solvers are implemented correctly and, moreover, that in the presence of smooth solutions both methods yield the same convergence properties.

Decay properties
In order to more robustly validate each method, we now consider a more challenging problem at varying degrees of resolution. In order to evaluate the numerical diffusion that is introduced by the presence of an interface, we consider the rotation of a Gaussian in a circular manner using the transport equation. More precisely we utilise the same transport equation as the previous setting but now consider the velocity v(x, y) = (−x, y), so that the initial scalar Gaussian field u(x), 0 = exp(− x − x 0 2 /σ 2 ) is rotated around the origin. We consider a domain Ω = [−2, 2], using an initial starting point x 0 = (−0.625, −0.625) with σ = 0.1. The mesh used in this test consists of 16 × 16 uniform quadrilateral cells for the conformal case, while in the non-conformal cases the right-hand half of the grid has been displaced by half a cell vertically along the central interface in relation to the left-hand side, as shown in Fig. 7. Constructing the mesh in this way aims to keep a consistent cell density by ensuring the half cell height sections are on the extreme ends of the interface, distant from where the peak crosses the interface. We also note that the selection of x 0 is designed to place the peak in the centre of a cell to the left of the interface, as well as to ensure minimal interaction with the domain boundaries which all have a homogeneous Dirichlet condition imposed on them.
The peak starts at t = 0 with a maximum value of 1. Unlike the exact solution, which precisely preserves this peak indefinitely, the non-polynomial nature of the solution field means that we can expect the peak to decrease every rotational cycle due to numerical diffusion introduced by each method. We then measure the L ∞ error precisely through the solution of a minimisation problem -i.e. we do not solely sample the error at quadrature points, as at lower orders very few quadrature points are used within each element, and this may lead to a significant difference in the observed error. We select a timestep size of ∆t = 0.001, and for each combination of polynomial order and interface handling method, we measure the L ∞ after 100 cycles of the Gaussian (i.e. t = 200π). We note that at lower polynomial orders, the solution will be underresolved by design -our aim in this series of simulation is to examine how this affects numerical stability across the methods and/or if there are significant differences in performance of the methods. The results of these experiments are shown in Table 1.  Two trends are evident from the presented results. At the lowest polynomial order of P = 4, we see a reasonable level of difference in the point-to-point method vs. the mortar and conformal grids. The oscillations of numerical error at this order are clearly evident in Fig. 7b which shows the 4 th order interpolation method after 100 cycles. Curiously, the values observed at P = 4 through P = 7 are higher for the point-to-point method than both the mortar/conformal cases, indicating that the point-to-point method is somewhat better able to resolve the peak of the Gaussian. It is also clear that as the polynomial order increases this difference in the maximum value decreases, so that at P ≥ 8 the results are essentially identical. Broadly speaking, however, the performance of the methods is reasonably comparable across the range of polynomial orders.

Long-time advection of an isentropic vortex
In order to examine the non-conformal methods in more realistic problems, whilst still considering their long-term stability and diffusion properties, we now move on to consider a nonlinear problem. In particular, we consider the compressible Euler equations in two dimensions. In this case, the conserved variables are given as u = [ρ, ρu, ρv, E] with ρ being the density, (u, v) the fluid velocity, E the specific total energy, and where p is the pressure. To close the system we need to specify an equation of state; in this case we use the ideal gas law p = ρRT where T is the temperature and R is the gas constant.
To consider long-term stability, we opt to study an isentropic vortex that is advected at constant velocity through periodic boundaries. This is a commonly used benchmark when testing numerical discretisation of the compressible Euler equations, particularly for higher-order codes, as it is one of the few problems that admits an exact solution calculable at all times whilst also being relatively simple to implement [37,38,39,40].
For our purposes, we consider a domain Ω = [−5, 5] 2 . At any given time t, the solution for the isentropic vortex is given by the equations We select an initial vortex position (x 0 , y 0 ) = (0, 0) with strength β = 5 and γ = 1.4, and advect the velocity in the x-direction with velocity (u 0 , v 0 ) = (1, 0). The initial vortex size and location can be seen in Fig. 8a. The concept behind this series of simulations is much the same as the rotating Gaussian peak; i.e. we wish to cycle the vortex through the domain a number of times, and compare the error as a function of time for each interface method. This will be undertaken for polynomial orders, P = 3 through P = 7 on a grid of fixed size, where the lower orders are expected to be under-resolved and the higher orders somewhat more resolved. To impose this, a single pair of periodic conditions at the constant x boundaries were used so that u(−5, y) = u(5, y), while the constant y boundaries were set to free-stream conditions. Although it may seem more natural to impose periodic boundaries on all of the edges of the domain, in a similar fashion to [39] we found that this leads to a gradual accumulation of numerical error, which left unchecked eventually causes simulations to diverge. The solution, proposed in [39] and [40], is to impose farfield conditions at constant y boundaries, which allows recirculated  waves of accumulated numerical error to escape the domain and avoid premature divergence. This is particularly important in this case, as in order to further reduce sources of artificial dissipation, we elect to use the exact Riemann solver of Toro [41] to calculate the numerical fluxf (u + , u − ). We note that although this is computationally expensive, cheaper solvers such as the Roe solver may introduce additional numerical diffusion [42].
One additional consideration that needs to be taken in this nonlinear regime is the order of integration used to evaluate integrals in the weak form of eq. (2). Aliasing errors are a well known phenomenon in this regime, due to the cubic nonlinearity that arises in the definition of the Euler equations, as well as the non-polynomial flux term which calculated between elements [37]. Typically it is necessary to use a higher order of quadrature than is used for linear problems, in order to remove sources of aliasing error due to under-integration of these terms. For this reason, we consider two different numbers of quadrature points with Q = P + 2 and Q = 2P + 2, respectively.
The conformal case consists of a singular domain made up of a 21 × 21 regular quadrilateral mesh as shown in Fig. 8a. The resulting non-conformal mesh consists of the domains 7 × 21, 7 × 22 and 7 × 21 elements, as shown in Fig. 8b.
The periodic boundary also allows us to conveniently express the time in cycles, where one cycle is the length of time taken for the vortex to propagate through the domain and return to its initial position. In our case, with a propagation speed of u 0 = 1 and domain length L = 10, this leads to the exact solution every t = 10. We can calculate the exact solution at any time, t, by moving the vortex centre by u 0 t in the x-direction and making sure to account for the periodic condition. We select a timestep size of ∆t = 0.001 and use the same explicit 4th-order Runge-Kutta timestepping scheme as in previous results.
In Fig. 9, we visualise the L 2 error of the density field ρ, denoted by L 2 ρ , for a simulation spanning 100 cycles of the vortex through the domain. This figure yields a number of interesting features that warrant further discussion. Firstly, as validation of our results, we note that the broad characteristics of the conformal error broadly agree with those seen in other work and, in particular, those of [39]. More generally, we observe that the conformal method and mortar method yield extremely close results for all polynomial and quadrature orders under observation, which we would expect given the similar levels of resolution and the local conservation properties of the mortar method. However, when considering the point-to-point interpolation method, there are indeed clear differences in comparison to the mortar and conformal methods. Perhaps the most obvious peculiarity is the possible relationship between odd numbers of quadrature points and the long term stability of the interpolation method; P3Q5 and P5Q7 in Fig. 9 show significant divergence at low cycle counts. The errors in this case appear to be related to aliasing error: as the integration order is increased to Q = 8 and Q = 12 respectively, the results remain consistent with those found by the mortar method and the reference conformal case.
To investigate the effect of aliasing and integration order further, additional point-to-point interpolation simulations were run at P = 5 with quadrature orders ranging between Q = 7 and Q = 13. Fig. 10 depicts the L 2 ρ error for these cases. The pronounced abnormality at Q = 7 is clearly visible, and indeed at Q = 8, there is a sudden increase in error after ∼ 70 cycles which is indicative of further long-time increases in error. However, for Q ≥ 9, we observe much more consistent trends and better agreement with the mortar and conformal cases.
More generally then, we can state that so long as appropriate levels of aliasing are used so that Q = 2P + 2, the interpolation method closely follows the same trend as the mortar method and the benchmark conformal case, with the same reduction in error as polynomial order increases. These cases are all visualised in a single Fig. 11 to highlight this more clearly.

Conclusions
In this paper, we have compared the numerical performance of the pointto-point interpolation and mortar techniques, together with equivalent conformal cases, for a number of linear and non-linear hyperbolic conservation law problems. For problems that admit smooth solutions (i.e. which are adequately resolved in space), it is clear that either method is capable of performing equally well, both in terms of preserving the high-order convergence properties of the DG method, and also when considering the advection of structures across very long time periods. Likewise, when considering problems that are marginallyor under-resolved, it is equally clear that the mortar technique yields the most consistently accurate results when compared to the point-to-point interpolation approach.
Although there were relatively minor differences between the point-to-point and mortar methods for the linear Gaussian hump case in the presence of under-    resolution, the isentropic vortex case clearly highlights the care that must be taken when using the point-to-point method in such a regime. From the results we observe here, aliasing can have a significant impact on the ability of this method to accurately resolve flow features across long time periods. However, from this perspective, it would be unusual for higher-order fluid dynamics simulations being performed in an implicit LES or under-resolved DNS regime without a significant level of dealiasing. As demonstrated in [43], running either compressible Euler or Navier-Stokes simulations without a comparable level of dealiasing to that we present here can yield inaccurate results and potentially lead to instability. Indeed, it is worth emphasising that in realistic fluid dynamics simulations, most problems consist of inflow-outflow setups in which structures would be naturally removed from the domain within a far shorter time period than we consider here. In this context, the point-to-point approach would likely be sufficient for most application purposes. On balance then, the authors believe that with appropriate precautions, either method is sufficiently capable of handling the non-conformal grids that appear as a result of sliding or moving meshes. This has important implications in terms of implementation, however. For three-dimensional geometries, the implementation challenge of constructing mortar elements across an arbitrary interface at high-order presents a significant challenge, which makes the point-to-point method an attractive alternative, particularly in the context of highly parallel simulations. Future work in this area could therefore consider the application of the point-to-point method in the context of a challenging three-dimensional simulation involving sliding meshes, and comparing this performance against normal conformal simulations as appropriate.