# Hermite-Discontinuous Galerkin Overset Grid Methods for the Scalar Wave Equation

We present high order accurate numerical methods for the wave equation that combines efficient Hermite methods with eometrically flexible discontinuous Galerkin methods by using overset grids. Near boundaries we use thin boundary fitted curvilinear grids and in the volume we use Cartesian grids so that the computational complexity of the solvers approach a structured Cartesian Hermite method. Unlike many other overset methods we do not need to add artificial dissipation but we find that the built in dissipation of the Hermite and discontinuous Galerkin methods is sufficient to maintain stability. By numerical experiments we demonstrate the stability, accuracy, efficiency and applicability of the methods to forward and inverse problems.

## Authors

• 1 publication
• 11 publications
• ### An Energy-Based Discontinuous Galerkin Method with Tame CFL Numbers for the Wave Equation

We extend and analyze the energy-based discontinuous Galerkin method for...
10/14/2021 ∙ by Daniel Appelo, et al. ∙ 0

• ### Sparse Grid Discretizations based on a Discontinuous Galerkin Method

We examine and extend Sparse Grids as a discretization method for partia...
10/25/2017 ∙ by Alexander B. Atanasov, et al. ∙ 0

• ### Fourier Continuation Discontinuous Galerkin Methods for Linear Hyperbolic Problems

Fourier continuation is an approach used to create periodic extensions o...
04/30/2021 ∙ by Daniel Appelo, et al. ∙ 0

• ### Multigrid with Nonstandard Coarsening

We consider the numerical solution of Poisson's equation on structured g...
08/10/2020 ∙ by Kamala Liu, et al. ∙ 0

• ### A Memory-efficient Implementation of Perfectly Matched Layer with Smoothly-varying Coefficients in Discontinuous Galerkin Time-Domain Method

Wrapping a computation domain with a perfectly matched layer (PML) is on...
06/04/2020 ∙ by Liang Chen, et al. ∙ 0

• ### A discontinuous Galerkin method based on a hierarchical orthogonal basis for Lagrangian hydrodynamics on curvilinear grids

We present a new high-order accurate Lagrangian discontinuous Galerkin (...
02/23/2021 ∙ by Xiaodong Liu, et al. ∙ 0

• ### The surrogate matrix methodology: Accelerating isogeometric analysis of waves

The surrogate matrix methodology delivers low-cost approximations of mat...
04/10/2020 ∙ by Daniel Drzisga, et al. ∙ 0

##### This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

## 1 Introduction

Accurate and efficient simulation of waves is important in many areas in science and engineering due to the ability of waves to carry information over large distances. This ability stems from the fact that waves do not change shape in free space. On the other hand when the background medium is changing this induces a change in the wave forms that propagate through the medium and the waves can be used for probing the interior material properties of objects.

In order to preserve the properties of waves from the continuous setting it is preferable to use high order accurate discretizations that are able to control dispersive errors. The development of high order methods for wave propagation problems has been an active area of research for a long time and there are by now many attractive methods. Examples include (but are not limited to) finite difference methods, SBP ; Virta2014 ; Wang2016 ; Petersson2018 ; Hagstrom2012 , embedded boundary finite differences, appelo2012fourth ; FCAD1 ; FCAD2 ; Li2004295 ; Wandzura2004763 , element based methods like discontinuous Galerkin (DG) methods, Wilcox:2010uq ; Upwind2 ; ChouShuXing2014 ; ChungEngquist06 ; ChungEngquist09 ; GSSwave , hybridized discontinuous Galerkin (HDG) methods, Nguyen2011 ; Stanglmeier2016 , cut-cell finite elements STICKO2016364 ; sticko2016higher and Galerkin-difference methods BANKS2016310 .

An advantage of summation-by-parts finite differences and Galerkin type methods is that stability is guaranteed, however this guarantee also comes with some drawbacks. For diagonal norm summation-by-parts finite differences the order of accuracy is reduced to roughly half of that in the interior near boundaries. Further the need for multi-block grids also restricts the geometrical flexibility.

As DG and HDG methods are naturally formulated on unstructured grids they have good geometric flexibility. However, Galerkin based polynomial methods often have the drawback that they require small timesteps (the difference Galerkin and cut-cell finite element methods are less affected by this) when combined with explicit timestepping methods, but on the other hand they preserve high order accuracy all the way up to the boundary and it is easy to implement boundary conditions independent of the order of the method.

The pioneering work by Henshaw and co-authors, see for example chess1990

, describe techniques for generating overset grids as well as how they can be used to solve elliptic and first order time-dependent partial differential equations (PDE) by second order accurate finite differences. In an overset grid method the geometry is discretized by narrow body-fitted curvilinear grids while the volume is discretized on one or more Cartesian grids. The generation of such body-fitted grids is local and typically produces grids of very high quality,

OGEN . The grids overlap (we say that they are overset) so that the solution on an interior (often referred to as non-physical or ghost) boundary can be transferred from the interior of another grid. In chess1990

and in most other overset grid methods the transfer of solutions between grids is done by interpolation. Since the bulk of the domain can be discretized on a Cartesian grid the efficiency asymptotically approaches that of a Cartesian solver but still retains the geometrical flexibility of an unstructured grid method. The same type of efficiency can be expected for embedded boundary and cut-cell finite elements but the errors close to physical boundaries are typically smoother and smaller when body-fitted grids are used.

Here we are concerned with the approximation of the scalar wave equation on overset grids. To our knowledge, high order overset grid methods for wave equations in second order form have been restricted to finite difference discretizations. For example, in henshaw:1730 high order centered finite difference approximations to Maxwell’s equations (written as a system of second order wave equations) was introduced. More recently, in ANGEL2018534 , the upwind discretizations by Banks and Henshaw introduced in BANKS20125854 were generalized to overset grids. A second order accurate overset grid method for elastic waves can be found in smog .

We use the recently introduced dissipative Hermite methods for the scalar wave equation in second order form, secondHermite , for the approximation on Cartesian grids. To handle geometry we use the energy based DG methods of Upwind2 on thin grids that are grown out from physical boundaries. We use projection to transfer the solutions between grids rather than interpolation.

Both the Hermite and DG methods we employ increase the order of accuracy by increasing the number of degrees of freedom on an element or cell. This has practical implications for grid generation as a single grid with minimal overlap can be used independent of order, reducing the complexity of the grid generation step. This can be important for example in problems like optimal shape design, where the boundary changes throughout the optimization. This is different from the finite difference methods where, due to the wider finite difference stencils, the overlap must grow as the order is increased.

The transfer of solutions between overset grids typically causes a perturbation to the discrete operators which, especially for hyperbolic problems, results in instabilities, see smog for example. These instabilities are often weak and can thus be suppressed by a small amount of artificial dissipation. There are two drawbacks of this added dissipation, first it is often not easy to determine the suitable amount needed, i.e. big enough to suppress instabilities but small enough not to reduce the accuracy or timestep too severely. Second, in certain cases the instabilities are strong enough that the dissipation must scale with the discretization parameter (the grid size) in such a way that the order of accuracy of the overall method is reduced by one.

Similar to ANGEL2018534 , we use a dissipative method that has naturally built–in damping that is sufficient to suppress the weak instabilities caused by the overset grids. The order of the hybrid overset grid method is the design order of the Hermite method or DG method, whichever is the smallest.

In the hybrid H–DG overset grid method the Hermite method is used on a Cartesian grid in the interior of the domain, and the discontinuous Galerkin method on another, curvilinear grid at the boundary. The numerical solution is evolved independently on these grids for one timestep of the Hermite method. By using the Hermite method in the interior the strict timestep constraints of the DG method are relaxed by a factor that grows with the order of the method. Asymptotically, as discussed above, the complexity of the hybrid H–DG solver approaches that of the Cartesian Hermite solver secondHermite .

The paper is organized as follows. The Hermite method is described in the next section. We first explain the method in simple one dimensional case and then explain how the method generalizes to two dimensions. The DG method is described in section 3. The details of the overset grids and a hybridization of the DG and the Hermite methods are described in section 4. We illustrate the hybrid H–DG method with numerical simulations in the section 5.

## 2 Dissipative Hermite method for the scalar wave equation

We present the Hermite method in some detail here and refer the reader to the original work secondHermite

for convergence analysis and error estimates.

Consider the one dimensional wave equation in second order form in space and first order in time

 ut = v, (1) vt = c2uxx+f,  x∈Ω,  t∈(0,T). (2)

Here , and for optimal convergence. We refer to as the displacement, and as the velocity. The speed of sound is . We consider boundary conditions of Dirichlet or Neumann type

 u(t,x) = h0(t,x),  x∈∂ΩD, ux(t,x) = h1(t,x),  x∈∂ΩN,

and initial conditions

 u(0,x) = g0(x), v(0,x) = g1(x).

Let the spatial domain be . The domain will be discretized by a primal grid

 xi=a+ih,  h=(b−a)/N,  i=0,…,N,

and a dual grid

 xi=a+ih, i=12,…,N−12.

The use of staggered grids is essential for being able to take large timesteps. In time we discretize using a uniform grid with increments , that is

 tn=nΔt, n=0,1/2,1,….

At each grid point the approximation to the solution is represented by its degrees of freedom (DOF) that approximate the values and spatial derivatives of and . Equivalently, the approximations to and can be represented as polynomials centered at grid points . The Taylor coefficients of these polynomials are scaled versions of the degrees of freedom. To achieve the optimal order of accuracy we require the and first derivatives of and respectively to be stored at each grid point.

At the initial time (which we take to be ) these polynomials are approximations to the initial condition on the primal grid

 u(x,0) v(x,0) ≈m∑l=0^vl(x−xih)l≡qi(x),  i=0,…,N.

The coefficients and are assumed to be accurate approximations to the scaled Taylor coefficients of the initial data. If expressions for the derivatives of the initial data are known we simply set

 ^ul=hll!dlg0dxl∣∣∣x=xi,  ^vl=hll!dlg1dxl∣∣∣x=xi. (3)

Alternatively, if only the functions and are known, we may use a projection or interpolation procedure to find the coefficients in (3).

The numerical algorithm for a single timestep consists of two phases, an interpolation step and an evolution step. First, during the interpolation phase the spatial piecewise polynomials are constructed to approximate the solution at the current time. Then, in the time evolution phase we use the spatial derivatives of the interpolation polynomials to compute time derivatives of the solution using the PDE. We compute new values of the DOF on the next time level by evaluating the obtained Taylor series. We now describe each step separately.

### 2.1 Hermite interpolation

At the beginning of a timestep at time (or at the initial time) we consider a cell and construct the unique local Hermite interpolant of degree for the displacement and degree for the velocity. The interpolating polynomials are centered at the dual grid points and can be written in Taylor form

 pi+12(x) = 2m+3∑l=0^ul,0(x−xi+12h)l,  x∈[xi,xi+1],  i=0,…,N−1, (4) qi+12(x) = 2m+1∑l=0^vl,0(x−xi+12h)l,  x∈[xi,xi+1],  i=0,…,N−1. (5)

The interpolants and are determined by the local interpolation conditions:

 dlpi+12dxl =dlpidxl∣∣∣x=xi,  dlpi+12dxl=dlpi+1dxl∣∣∣x=xi+1,  l=0,…,m+1, dlqi+12dxl =dlqidxl∣∣∣x=xi,  dlqi+12dxl=dlqi+1dxl∣∣∣x=xi+1,  l=0,…,m.

We find the coefficients in (4) and (5) by forming a generalized Newton table as described in hagstrom2015solving .

### 2.2 Time evolution

To evolve the solution in time we further expand the coefficients of and . At each point on the dual grid, we seek temporal Taylor series

 pi+12(x,t) = 2m+3∑l=0κp∑s=0^ul,s(x−xi+12h)l(tΔt)s, (6) qi+12(x,t) = 2m+1∑l=0κq∑s=0^vl,s(x−xi+12h)l(tΔt)s, (7)

where and . The coefficients and are given by the coefficients of (4) and (5). At this time the scaled time derivatives, and , are unknown and must be determined. Once they are determined we may simply evaluate (6) and (7) at to find the solution at the next half timestep.

In Hermite methods the coefficients of temporal Taylor polynomials are determined by collocating the differential equation, good2006 ; secondHermite ; hagstrom2015solving . In particular, by differentiating (1) and (2) in space and time the time derivatives of the solution can be directly expressed in terms of spatial derivatives

 ∂s+1+ru∂ts+1∂xr = ∂s+rv∂ts∂xr, (8) ∂s+1+rv∂ts+1∂xr = c2∂s+r+2u∂ts∂xr+2∂su∂ts∂x2+∂s+rf∂ts∂xr. (9)

Substituting (6) and (7) into (8) and (9) and evaluating at and , we can match the powers of the coefficients to find the recursion relations

 ^ul,s+1 = Δts^vl,s, (10) ^vl,s+1 = c2(l+1)(l+2)h2Δts^ul+2,s+Δts^fl,s. (11)

Here are the coefficients of the Taylor expansion of , or of the polynomial which interpolates in time around . Note that since there are a finite number of coefficients, representing the spatial derivatives at the time , the recursions truncate and only and terms need to be considered.

To complete a half timestep we evaluate the approximation at for the and first derivatives

 ∂lu∂xl(xi+12,tn+12)≈∂lpi+12∂xl(xi+12,tn+12)=l!hlκp∑s=0^ul,s2s,  l=0,…,m+1, (12) ∂lv∂xl(xi+12,tn+12)≈∂lqi+12∂xl(xi+12,tn+12)=l!hlκq∑s=0^vl,s2s,  l=0,…,m. (13)
###### Remark 1

For a piecewise polynomial solution to the wave equation all the terms in the Taylor series are present in the right-hand side of equations (12)-(13). A consequence is that the time evolution is exact whenever the forcing is zero (or a polynomial of degree ) and each cell includes the domain of the dependence of the solution at a dual grid point at time , that is, when

 cΔt2≤h2.

### 2.3 Imposing boundary conditions

Physical boundary conditions are enforced at the half time level, i.e. when the solution on the dual grid is to be advanced back to the primal grid. As there are many degrees of freedom that are located on the boundary the physical boundary condition must be augmented by the differential equation to generate more equations so that the degrees of freedom can be uniquely determined. The basic principle, often referred to as compatibility boundary conditions (see e.g. henshaw:1730 ), is to take tangential derivatives of the PDE and the boundary conditions so that a sufficient number of boundary conditions are obtained.

For example, assume we want to impose the boundary condition

 u(t,0)=g(t). (14)

Then, as is a boundary grid point the Taylor polynomials (6)-(7) centered at should satisfy the boundary condition (14). In addition it should also satisfy the differential equation and derivatives of (14). We thus seek a polynomial outside the domain which together with the polynomial just inside the boundary forms a Hermite interpolant that satisfies the boundary and compatibility conditions.

Precisely, to find the polynomial describing at time we must determine unknowns that specify the polynomial at the boundary. First, this polynomial must interpolate the data describing the current approximation of at , this yields independent linear equations. Second, the first of the remaining independent linear equations can be obtained by requiring that the polynomial coincides with the boundary condition . The next equation is , and so forth.

Once the interpolant is determined on the boundary we evolve it as in the interior (see section 2.2).

###### Remark 2

We note that in the special case of a flat boundary and homogeneous Dirichlet or Neumann boundary conditions then enforcing the boundary conditions reduces to enforcing that the polynomial on the boundary is either odd or even, respectively, in the normal direction. Then the correct odd polynomial can be obtained by constructing the polynomial outside the domain

(often referred as ghost-polynomial) by mirroring the coefficients corresponding to even powers in the normal coordinate variable with a negative sign and the coefficients corresponding to odd powers with the same sign.

Boundary conditions at interior overset grid boundaries are supplied by projection of the known solutions from other grids and will be discussed below.

### 2.4 Higher dimensions

In higher dimensions the approximations to and

take the form of centered tensor product Taylor polynomials. In two dimensions (plus time) the coefficients would be of the form

, with the two first indices representing the powers in the two spatial directions, and the third representing time.

For the scalar wave equation

 ut =v, vt =c2(uxx+uyy),  (x,y)∈Ω,  t>0,

the recursion relations for computing the time derivatives are a straightforward generalization of the one dimensional case

 ck,l,s=Δtsdk,l,s−1, (15) dk,l,s=c2(k+2)(k+1)sΔth2xck+2,l,s−1+c2(l+2)(l+1)sΔth2yck,l+2,s−1. (16)

As noted in secondHermite , using this recursion for all the time derivatives does not produce a method with order independent CFL condition but a method whose time-step size decrease slightly as the order increases. For optimally large timesteps it is necessary to use the special start up procedure

 ck,l,1 =Δtdk,l,0, dk,l,1 =(cΔt)2((k+2)(k+1)h2xcXk+2,l+(l+2)(l+1)h2ycYk,l+2).

Here are the coefficients of the interpolating polynomial of degree in and degree in and are the coefficients of the interpolating polynomial of degree in and degree in . For the remaining coefficients we use (15) and (16) with . Further details of the two dimensional method can be found in secondHermite .

## 3 Energy based discontinuous Galerkin methods for the wave equation

Our spatial discontinuous Galerkin discretization is a direct application of the energy based formulation described for general second order wave equations in Upwind2 ; el_dg_dath ; fluid_solid_DG . Here, our energy based DG method starts from the energy of the scalar wave equation

 H(t)=∫Ωv22+G(x,y,∇u)dΩ,

where is the potential energy density, is the velocity or the time derivative of the displacement, .

Now, the wave equation, written as a second order equation in space and first order in time takes the form

 ut=v,  vt=−δG,

where is the variational derivative of the potential energy

 δG=−∇⋅(c2(x,y)∇u).

For the continuous problem the change in energy is

 dH(t)dt=∫Ωvvt+ut[∇⋅(c2(x,y)∇u]dΩ=[ut(n⋅(c2(x,y)∇u))]∂Ω,

where the last equality follows from integration by parts together with the wave equation.

A variational formulation that mimics the above energy identity can be obtained if the equation is tested with the variational derivative of the potential energy. Let be an element and and be the spaces of tensor product polynomials of degrees and . Then, the variational formulation on that element is:

###### Problem 1

Find , such that for all ,

 ∫Ωj(c2∇ϕ)⋅(∂∇uh∂t−∇vh)dΩ =[(c2∇ϕ)⋅n(v∗−vh)]∂Ωj, (17) ∫Ωjψ∂vh∂t+c2∇ψ⋅∇uhdΩ =[ψ(c2∇u⋅n)∗]∂Ωj. (18)

Let and denote the jump and average of a quantity at the interface between two elements, then, choosing the numerical fluxes as

 v∗ ={vh}−τ1[[c2∇uh⋅n]], (c2∇u⋅n)∗ ={c2∇uh⋅n}−τ2[[vh]],

yields a contribution from each element face to the change of the discrete energy, guaranteeing that

 dHh(t)dt≡ddt∑j∫Ωj(vh)22+G(x,y,∇uh)≤0.

Physical boundary conditions are enforced through the numerical fluxes, see Upwind2 for details.

Note that the above energy estimate follows directly from the formulation (17) - (18) but as the energy is invariant to constants equation (17) must be supplemented by the equation

 ∫Ωj(∂uh∂t−vh)dΩ=0.

Our implementation uses quadrilaterals and approximations by tensor product Chebyshev polynomials of the solution on the reference element . That is, on each quadrilateral we have approximations on the form

 u(x(r,s),y(r,s),tn) ≈qu∑l=0qu∑k=0clkTl(r)Tk(s), v(x(r,s),y(r,s),tn) ≈qv∑l=0qv∑k=0dlkTl(r)Tk(s).

We choose (so called upwind or Sommerfeld fluxes) which result in methods where is observed to be order accurate in space Upwind2 .

### 3.1 Taylor series time-stepping

In order to match the order of accuracy in space and time for the DG method we employ Taylor series time-stepping. Assuming that all the degrees of freedom have been assembled into a vector

we can write the semi-discrete method as with being the matrix representing the spatial discretization. If we know the discrete solution at the time we can advance it to the next time step by the simple formula

 w(tn+Δt) = w(tn)+Δtwt(tn)+(Δt)22!wtt(tn)… = w(tn)+ΔtAw(tn)+(Δt)22!A2w(tn)…

As we use dissipative fluxes this timestepping method is stable as long as the number of stages in the Taylor series is greater than the order of accuracy in space and with the timestep small enough.

## 4 Overset grid methods

In this section we explain how we use the two discretization techniques described above on overset grids to approximate solutions to the scalar wave equation.

The idea behind the overset grid methods is to cover the bulk of the domain with a Cartesian grid, where efficient methods can be employed, and to discretize the geometry with narrow body-fitted grids. In Figure 1 we display two overset grids, a blue Cartesian grid, which we denote , and a red curvilinear grid, which we denote , that are used to discretize a geometry consisting of a circular hole cut out from a square region. Note that the grids overlap, hence the name overset grids. Also, note that the annular grid cuts out a part of the Cartesian grid. This cut of the Cartesian grid creates a internal, non-physical boundary in the blue grid.

Here physical boundary conditions are enforced on the red grid at the black boundary which defines the inner circle and on the outermost boundary on the blue grid.

In order to use the Hermite or DG methods on the grids we will need to supply boundary conditions at the interior boundaries. In the example in Figure 1 this means that we would have to specify the solution on the outer part of the annular grid and on the staircase boundary (marked with filled black circles) that has been cut out from the Cartesian grid.

In most methods that use overset grids, in particular those using finite differences, the communication of the solution on the interior boundaries is done by interpolation, see e.g. chess1990 . For the methods we use here we have found that the stability properties are greatly enhanced if we instead transfer volumetric data (numerical solution) in the elements / gridpoints near the internal boundaries by projection rather than by interpolation. In fact, when we use volume data the resulting methods are stable without adding artificial dissipation, when we use interpolation they are not.

As mentioned above, in a Hermite method, we can think of the degrees of freedom as either being nodal data, consisting of function and derivative values, or as coefficients in a Taylor polynomial. Thus, when transferring data to a grid where a Hermite method is used (like the example in the left subfigure of Figure 1) we must determine a tensor product polynomial centered around a gridpoint local to that grid (the points we would center around are indicated by black points in Figure 1). Below we will explain in detail how we determine this polynomial.

For elements with an internal boundary face (denoted by thick red lines in Figure 1) we could in principle transfer the solution by specifying a numerical flux on that face, however we have found that this approach results in weakly unstable methods. Instead we transfer volumetric data to each element that has an internal boundary face, we give details below. Given the timestep constraints of DG methods we must march the DG solution using much smaller timesteps than those used for the Hermite method. This necessitates the evaluation of the Hermite data not only at the beginning of a Hermite timestep but at many intermediate times.

### 4.1 Determining internal boundary data for the Hermite solver

We first consider the problem of determining internal boundary data required by the Hermite method. An example of how to compute solution data at the gridpoints at the boundary of Cartesian grid (filled black circles) is depicted in Figure 1.

In general, the tensor product polynomial centered around is found by a two step procedure. First we project into a local basis spanned by Legendre polynomials and perform a numerically stable and fast change of basis into the monomial basis. Then we truncate the monomial to the degree required by the Hermite method.

To carry out the projection we introduce a local tensor product Gauss-Legendre-Lobatto (GLL) grid centered around . These points are marked as filled blue circles in the left subfigure of Figure 2. The number of grid points in the local grids are determined by the order of the projection. To maintain the order of the method, the order of the projection should be at least the same as the order of the spatial discretization, thus it is sufficient to have points in each direction. The GLL quadrature nodes are defined on the reference element that maps to a cell defined by the dual gridpoints closest to .

Let be the numerical solution on the red grid. In the first step of the communication we compute the coefficients of a polynomial approximating by projecting on the space of tensor product Legandre polynomials , that is

 ~p(r,s)=2m+3∑l=02m+3∑k=0clkPl(r)Pk(s), clk=(~u,PlPk)∥PlPk∥2. (19)

Here denotes the inner product on and is the norm induced by the inner product. Note that the expression (19) is particularly simple since the Legendre polynomials are orthogonal on the domain of integration. To do this we evaluate at the underlying blue quadrature points in the left subfigure of Figure 2.

Once the polynomial (19) has been found we perform a change of basis into the local monomial used by the Hermite method. Such a change of basis can be done by the fast Vandermonde techniques by Björk and Pereyra, see e.g. bp70 ; Dahlquist:2008fu . At this stage the polynomial is of total degree so the final step is to truncate it to total degree or depending on whether we are considering the displacement or the velocity. With the and degrees of freedom determined everywhere on a Hermite grid we may evolve the solution as described in section 2.

### 4.2 Determining data for DG elements with internal boundary faces

We now consider the problem of determining the data required by the DG method. Here we show how to obtain the data at a single DG element with at least one internal boundary face. As the timesteps of the DG method are significantly smaller than for the Hermite method we must repeat the transfer of data many times. We must also explicitly transfer time derivative data in order to use a Taylor series timestepping approach.

The tensor product polynomials in our implementation of the DG method are composed by the product of Chebyshev polynomials that are expressed on the reference element . Precisely we seek

 p(r,s)=q∑l=0q∑k=0clkTl(r)Tk(s).

To determine such polynomials we perform a projection of the solution , i.e the solution on Cartesian grid,

 clk=(~u,TlTk)C∥TlTk∥2C,

but in this case the weighted inner product is

 (f,g)C=∫1−1∫1−1f(r,s)g(r,s)√1−r2√1−s2drds,

where the Chebyshev polynomials are orthogonal. To carry out this projection we use a local tensor product Chebyshev quadrature nodes, in each dimension, as shown in right subfigure of Figure 2.

Denoting The local time levels used by the DG solver th Hermite timestep are defined to be

 tn,ν=tn,0+νΔtb, ν=0,…NDG,

where and are timesteps taken on grids (Cartesian) and (curvilinear) respectively. For simplicity the starting local time level and the final local time level are equal to consequent timesteps on the Hermite grid, and

 tn,0=tn,  tn,NDG=tn+1.

To transfer the solution values and the time derivatives needed at each of the quadrature points and at each we carry our the following “start up” procedure at . For each of the quadrature points we re-center the Hermite interpolants closest to it and compute the time derivatives precisely by the recursion relations described in section 2. We note that this is an inexpensive computation as the interpolants have already been found as a step in the evolution of the Hermite solution, the only added operation is the re-centering.

## 5 Numerical experiments

The hybrid H–DG method is empirically stable and accurate, and here we demonstrate it with numerical experiments. To test the stability of the method in one dimension we first define the amplification matrix and compute its spectral radius. To test the stability in two dimensions, where the amplification matrix will take too long to compute, we provide the long time simulation and estimate the error growth for multiple refinements. Convergence tests in one and two dimensions are done for the domains where the exact solution is known. In the second half of this section we apply the method to the domain with complex curvilinear boundary in the experiment with wave scattering of the smooth pentagonal object. Finally, in the end of this section we apply the method to the inverse problem of locating the underground cavities as the forward solver.

### 5.1 Numerical stability test

Unlike the Hermite and DG methods, stability of the hybrid H–DG method cannot be shown analytically. As an alternative, the stability can be investigated numerically by looking at the spectrum of the amplification matrix associated with the method. A similar stability analysis was done for finite difference scheme for the wave equation in Strik . To construct the amplification matrix we apply the method to initial data composed of the unit vectors. The vector that is returned after one timestep is then placed as columns in a square matrix. The spectral radius of the amplification has to be smaller then for the method to be stable.

We consider the wave equation (1)-(2) on the unit interval with homogeneous Dirichlet and Neumann boundary conditions at and respectively. We introduce two uniform Cartesian grids which overlap inside a small interval close to one of the boundaries. Precisely, the grids are

 Ωa ={xai=iha,  i=0,…,na}, Ωb ={xbi=1−(nb−i)hb,  i=0,…,nb}.

The Hermite method is used on a grid and the DG method is used on grid . The grids thus overlap inside the interval . Here the ratio of the overlap size and the discretization width is . This ratio is fixed for all values of and . We also fix so that the amount of work done on grid is constant per timestep for all refinements. Fixing the ratio and makes the efficiency of the overall method to asymptotically be the determined by the efficiency the Hermite method.

Let be a vector holding the degrees of freedom of both methods at th timestep, then we may express the complete timestep evolution as where incorporates timestepping and projection. can be expressed as the matrix that can be computed column by column via

 Hk=Hek, (20)

where is the th unit vector. The equation

 wn=Hnw0, (21)

is equivalent to the timesteps of the hybrid H–DG method. Taking -norm of both sides of equation (21) and applying the Cauchy-Schwartz inequality we get

where is the spectral radius of . If we conclude

 ∥wn∥2≤∥∥w0∥∥2,

and the solution will remain bounded.

We consider the case and take the parameters to be

 na=10,20,…,60,  nb=5,  hbha=0.9.

Other parameters are for the DG method and , the number of timesteps done by the DG method during one step of the Hermite method. The parameters and are set so the methods used have the same order of accuracy as the approximation of for the Hermite method

 qu=2m+2, qv=2m+1.

To get an optimal , we take the largest possible timestep for the DG method, so that

 Δtbhb≤0.15/qu,

and

 nDG=ΔtaΔtb,

is an integer. Equivalently, if the Hermite method CFL number is set, we get

 nDG=ΔtaΔtb=⌈CFLqu0.15hahb⌉.

Following the column-by-column construction process (20) described at the start the this section we compute the amplification matrix for the current setting. The spectrum of is shown in Figure 3 for . Displayed results are for the cases and . The CFL numbers set for Hermite method are and . The absolute value of eigenvalues do not exceed . We note that if interpolation is used some eigenvalues of the amplification matrix shift outside of the unit circle. Such unstable modes can possibly be stabilized by numerical dissipation / hyperviscosity but we do not pursue such stabilization here. Instead we observe that when projection is used all eigenvalues are inside the unit circle and the method is stable. Although we only display the results for one problem here the same results were obtained for other grid sizes, various overlap sizes to grid spacing ratios and different CFL numbers set for the Hermite method. We stress that it is possible to make the method unstable if we take the CFL number close to one and if we take to be larger than 3 and thus we only claim that the methods of orders of accuracy up to are stable.

### 5.2 Convergence to an exact solution

Using the same grid setup and boundary conditions as in the example above we test the method for the wave equation (1)-(2), and initial conditions

 u(x,0) =sin(15π2x), (22) v(x,0) =0. (23)

A solution to this problem is the standing wave

 u(x,t)=sin(15π2x)cos(15π2t).

The errors for the solution on the grids are

 εa(x,t)=pi+12(x,t)−u(x,t),x∈xai,xai+1, i=0,…,na, (24)

for the Hermite grid and

 εb(x,t)=uhb(x,t)−u(x,t), (25)

for the DG grid. The maximum error for the total method is

 max(maxx∈[xa0,xana]|εa(x,t)|),maxx∈[xb0,xbnb|εb(x,t)|)).

In Figure 4 we display computed maximum errors as functions of time for the method with (i.e. the order of accuracy is 7). In the left subfigure the CFL number for the Hermite method is set to be and in the right subfigure the CFL number is set to be 0.75. For all Hermite grid sizes, the error growth is linear in time (dashed lines display a least squares fit of a linear function), indicating that the solution is stable for long time computations.

In the left subfigure of Figure 5 the numerical solution and the absolute error are shown for the th order accurate method at time . As can be seen in the lower left subfigure in Figure 5 the error is rather smooth across the overlap indicating that the projection is highly accurate.

To the right in Figure 5 the error at the final time is shown as a function . The dashed lines show the least squares fit with polynomial functions of of order and respectively. The results indicate that the orders of accuracy of the methods are as expected. The parameters (, , , etc.) are the same is in previous example.

### 5.3 Analytical solution in a disk. Rates of convergence

Consider the solution of (1)-(2) with on the unit disk, , with homogeneous Dirichlet boundary conditions. Then the analytical solution can be expressed in polar coordinates as a composition of modes

 uμν(r,θ,t)=Jμ(rκμν)cos(μθ)cosκμνt. (26)

Here is the Bessel function of the first kind of order and is the th zero of . In the following experiment we set , . The initial condition is displayed in the left subfigure of Figure 6.

We setup overset grids as displayed in Figure 7. Grid is a Cartesian grid discretizing a square domain with grid points in each direction and grid spacing . Grid is a curvilinear grid discretizing a thin annulus with radial grid spacing . For all refinements Grid has 7 elements in the radial direction thus the number of elements (or equivalently the number of DOFs of DG method) will grow linearly with the reciprocal of the discretization size . In contrast the number of grid points in the Cartesian grid where the Hermite method will be used grows quadratically with .

To measure the error we evaluate the solution on a finer grid, oversampled with 20 grid points inside each Hermite cell and DG element. The convergence is displayed in the right subfigure of Figure 6. The errors at time as functions of for are displayed as solid lines. The dashed lines show the polynomials in of order . We use in the computations. As can be seen the expected orders of accuracy (3,5 and 7) are observed.

To test the performance of the method we evolve the method over one time period of the solution and measure the CPU time, see the left subfigure of Figure 8. The red curve, displaying the error of the rd, order accurate method only reaches the error in about 1000 seconds while the th and th order accurate methods, using the same compute time, yield errors on the order of and respectively. Clearly the higher order methods are more efficient.

To test the stability of the method we evolve the solution until time which is roughly periods of the solution. We set and test methods with orders of accuracy and . The error growth appears to be linear in time as indicated by dashed lines in the right subfigure of Figure 8.

### 5.4 A wave scattering of a smooth pentagon

In this experiment we study the scattering of a smooth pentagon in free-space. In addition to the use of non-reflecting boundary conditions experiment demonstrates the hybrid Hermite-DG method for the solution which is propagated over many wavelengths. The geometry of the pentagon is defined as the smooth closed parametric curve:

 x(s) = 110(1+110cos(10s))cos(s), (27) y(s) = 110(1+110cos(10s))sin(s),  s∈[0,2π). (28)

The pentagon is placed in a square domain discretized by a Cartesian grid with grid spacing , . The curvilinear grid has 10 elements in the radial direction and the outer boundary is a circle of radius . The overlap width is at most 5 DG elements.

On the boundary of the body we set Dirichlet data

 u(x,y,t)=sin(ωt),  (x,y)∈Γ,  t≥0,  ω=250. (29)

The exterior boundary condition is modeled by truncating the domain using perfectly matched layers governed by the equations, (see appelo2012fourth for derivation)

 utt=∂∂x(ux+σxϕ(1))+∂∂y(uy+σyϕ(2))σ(x)ϕ(3)+σ(y)ϕ(4), (30)

where the auxiliary variables satisfy the equations

 ϕ(1)t+(α+σ(x))ϕ(1)=−ux,ϕ(2)t+(α+σ(y))ϕ(2)=−uy,ϕ(3)t+(α+σ(x))ϕ(3)=−uxx−∂∂x(σ(x)ϕ(1)),ϕ(4)t+(α+σ(y))ϕ(4)=−uyy−∂∂y(σ(y)ϕ(2)). (31)

The damping profiles , are taken as

 σ(z)(z)=σs(tanh(z−z10.7wlay)−tanh(z−z20.7wlay)).