Energy conservative SBP discretizations of the acoustic wave equation in covariant form on staggered curvilinear grids

TL;DR

This paper introduces an energy conserving numerical method for the acoustic wave equation on curvilinear grids, utilizing covariant basis decomposition and high-order SBP operators to enhance accuracy and stability.

Contribution

The paper presents a novel covariant basis decomposition approach combined with SBP operators for energy conserving discretization on curvilinear grids, improving rotational invariance and stability.

Findings

Outperforms Cartesian basis on rotated grids

Provides a conditionally stable discretization method

Offers bounds for stability evaluation

Abstract

We develop a numerical method for solving the acoustic wave equation in covariant form on staggered curvilinear grids in an energy conserving manner. The use of a covariant basis decomposition leads to a rotationally invariant scheme that outperforms a Cartesian basis decomposition on rotated grids. The discretization is based on high order Summation-By-Parts (SBP) operators and preserves both symmetry and positive definiteness of the contravariant metric tensor. To improve accuracy and decrease computational cost, we also derive a modified discretization of the metric tensor that leads to a conditionally stable discretization. Bounds are derived that yield a point-wise condition that can be evaluated to check for stability of the modified discretization. This condition shows that the interpolation operators should be constructed such that their norm is close to one.

Tables

Table 1. Table 1: Errors and convergence rates for two discretizations of the contravariant metric tensor. Results for the metric tensor discretization defined in ( 26 ) are labeled ( G ) 𝐺 (G) and the modified discretization defined in ( 28 ) are labeled ( G ~ ) ~ 𝐺 (\widetilde{G}) . Here, q h subscript 𝑞 ℎ q_{h} and q ∞ subscript 𝑞 q_{\infty} denote the convergence rates in the l 2 subscript 𝑙 2 l_{2} and l ∞ subscript 𝑙 l_{\infty} norms, respectively.

N	${\Vert {\text{Error}}^{\left(G\right)}\Vert }_h$	${\Vert {\text{Error}}^{\left(G\right)}\Vert }_{\infty }$	$q_h^{\left(G\right)}$	$q_{\infty }^{\left(G\right)}$	${\Vert {\text{Error}}^{\left(\tilde{G}\right)}\Vert }_h$	${\Vert {\text{Error}}^{\left(\tilde{G}\right)}\Vert }_{\infty }$	$q_h^{\left(\tilde{G}\right)}$	$q_{\infty }^{\left(\tilde{G}\right)}$
16	2.49e-2	1.05e-1	-	-	2.44e-2	1.06e-1	-	-
32	2.88e-3	1.49e-2	3.11	2.81	2.00e-3	1.56e-2	3.61	2.76
64	4.03e-4	2.89e-3	2.84	2.37	2.31e-4	2.39e-3	3.11	2.70
128	6.32e-5	6.89e-4	2.67	2.07	3.59e-5	4.16e-4	2.69	2.52
256	1.06e-5	1.77e-4	2.58	1.96	6.24e-6	9.84e-5	2.53	2.08

Table 2. Table 2: Errors and convergence rates for boundary characteristic variables (c) and non-boundary characteristic variables (nc) are reported. Errors are measured for the grid points in the SBP boundary-modified region of the bottom boundary.

h	${\Vert {\text{Err}}^{\left(c\right)}\Vert }_h$	${\Vert {\text{Err}}^{\left(c\right)}\Vert }_{\infty }$	$q_h^{\left(c\right)}$	$q_{\infty }^{\left(c\right)}$	${\Vert {\text{Err}}^{\left(nc\right)}\Vert }_h$	${\Vert {\text{Err}}^{\left(nc\right)}\Vert }_{\infty }$	$q_h^{\left(nc\right)}$	$q_{\infty }^{\left(nc\right)}$
6.25e-2	3.81e-3	9.82e-3	-	-	9.73e-3	2.43e-2	-	-
3.12e-2	6.12e-4	1.27e-3	2.64	2.95	1.10e-3	4.27e-3	3.14	2.51
1.56e-2	9.49e-5	3.95e-4	2.69	1.69	1.15e-4	7.00e-4	3.26	2.61
7.81e-3	1.29e-5	1.14e-4	2.87	1.79	1.41e-5	1.17e-4	3.02	2.58
3.91e-3	2.31e-6	3.05e-5	2.48	1.90	2.04e-6	2.04e-5	2.80	2.52

Table 3. Table 3: Errors and convergence rates obtained for the disc geometry with a cavity. Results for the covariant (co) and Cartesian (ca) formulations are reported separately.

N	${\Vert {\text{Err}}^{\left(co\right)}\Vert }_h$	${\Vert {\text{Err}}^{\left(co\right)}\Vert }_{\infty }$	$q_2^{\left(co\right)}$	$q_{\infty }^{\left(co\right)}$	${\Vert {\text{Err}}^{\left(ca\right)}\Vert }_h$	${\Vert {\text{Err}}^{\left(ca\right)}\Vert }_{\infty }$	$q_h^{\left(ca\right)}$	$q_{\infty }^{\left(ca\right)}$
16	1.69e-2	4.78e-2	-	-	6.65e-2	2.71e-1	-	-
32	6.74e-4	2.43e-3	4.18	3.87	2.72e-3	1.28e-2	4.15	3.96
64	4.68e-5	3.29e-4	3.67	2.75	2.04e-4	1.41e-3	3.56	3.04
128	5.59e-6	6.67e-5	2.99	2.25	2.68e-5	3.20e-4	2.86	2.09
256	8.62e-7	1.55e-5	2.67	2.08	4.41e-6	7.87e-5	2.57	2.00

Table 4. Table 4: Relative errors and convergence rates obtained for a point source placed on top of a curved boundary. Errors are computed using a numerical reference solution with N λ = 160 subscript 𝑁 𝜆 160 N_{\lambda}=160 grid points per minimum wavelength.

$N_{\lambda }$	${\Vert \Delta p\Vert }_{𝐫}$	$q_h^{\left(p\right)}$	${\Vert \Delta v^1\Vert }_{𝐫}$	$q_h^{\left(v^1\right)}$	${\Vert \Delta v^2\Vert }_{𝐫}$	$q_h^{\left(v^2\right)}$
5	1.20e-01	2.75	9.93e-02	3.14	1.20e-01	2.77
10	6.32e-03	4.25	4.77e-03	4.38	6.40e-03	4.22
20	3.16e-04	4.32	2.52e-04	4.24	3.31e-04	4.28
40	1.64e-05	4.26	1.45e-05	4.12	1.77e-05	4.22

Equations

$$\frac{\partial p}{\partial t}+\nabla \cdot \mathbf{v}=0,$$

$$\frac{\partial \mathbf{v}}{\partial t}+\nabla p=0,$$

$$0\le r^1\le 1\text{mbox}and0\le r^2\le 1,$$

$${\mathbf{a}}_1=\left[\begin{array}{c}\frac{\partial X}{\partial r^1}\\ \frac{\partial Y}{\partial r^1}\end{array}\right]\text{mbox}and{\mathbf{a}}_2=\left[\begin{array}{c}\frac{\partial X}{\partial r^2}\\ \frac{\partial Y}{\partial r^2}\end{array}\right].$$

$${\mathbf{a}}^i\cdot {\mathbf{a}}_j={\delta }_j^i=\{\begin{array}{ll}1,& i=j\\ 0,& i\ne j\end{array}.$$

$$\mathbf{v}=v^1{\mathbf{a}}_1+v^2{\mathbf{a}}_2,$$

$$\left[\begin{array}{c}{v}^1\\ v^2\end{array}\right]=\left[\begin{array}{cc}{g}^{11}& g^{12}\\ g^{12}& g^{22}\end{array}\right]\left[\begin{array}{c}{v}_1\\ v_2\end{array}\right].$$

$$\left[\begin{array}{c}{v}_1\\ v_2\end{array}\right]=\left[\begin{array}{cc}{g}_{11}& g_{12}\\ g_{12}& g_{22}\end{array}\right]\left[\begin{array}{c}{v}^1\\ v^2\end{array}\right].$$

$$\nabla \cdot \mathbf{v}=\frac{1}{J}\frac{\partial J{v}^1}{\partial r^1}+\frac{1}{J}\frac{\partial J{v}^2}{\partial r^2}.$$

$$\nabla p={\mathbf{a}}^1\frac{\partial p}{\partial r^1}+{\mathbf{a}}^2\frac{\partial p}{\partial r^2}.$$

$$J\frac{\partial p}{\partial t}+\frac{\partial }{\partial r^1}\left(J{v}^1\right)+\frac{\partial }{\partial r^2}\left(J{v}^2\right)$$

$$\frac{\partial }{\partial t}\left[\begin{array}{c}{v}^1\\ v^2\end{array}\right]+\left[\begin{array}{cc}{g}^{11}& g^{12}\\ g^{12}& g^{22}\end{array}\right]\left[\begin{array}{c}\frac{\partial p}{\partial r^1}\\ \frac{\partial p}{\partial r^2}\end{array}\right]$$

$$e(t)=\frac{1}{2}\int p^{2}Jd{r}^{1}d{r}^2+\frac{1}{2}\int (v_1{v}^1+v_2{v}^2)Jd{r}^{1}d{r}^2.$$

$$e(t)=\frac{1}{2}\int p^{2}Jd{r}^{1}d{r}^2+\frac{1}{2}\sum _{i,j}\int g_{ij}{v}^i{v}^{j}Jd{r}^{1}d{r}^2.$$

$$\frac{de}{dt}$$

$$=-\int \left(\frac{\partial (pJ{v}^1)}{\partial r^1}+\frac{\partial (pJ{v}^2)}{\partial r^2}\right)d{r}^{1}d{r}^2$$

$$=-\int {\left[pJ{v}^1\right]}_{r^1=0}^{1}d{r}^2-\int {\left[pJ{v}^2\right]}_{r^2=0}^{1}d{r}^1.$$

$$\frac{de}{dt}=0.$$

$$\mathbf{x}=[x_0{x}_1\dots x_N{]}^T\text{mbox}and\hat{\mathbf{x}}=[{\hat{x}}_0{\hat{x}}_1\dots {\hat{x}}_{N+1}{]}^T.$$

$$x_0$$

$${\hat{x}}_0$$

$$\hat{D}\mathbf{\xi }=2\hat{\mathbf{\xi }},D\hat{\mathbf{\xi }}=2\mathbf{\xi },\mathbf{\xi }=[x_0^2{x}_1^2\dots x_N^2{]}^T,\hat{\mathbf{\xi }}=[{\hat{x}}_0^2{\hat{x}}_1^2\dots {\hat{x}}_{N+1}^2{]}^T.$$

$${\int }_{x=0}^{1}u\frac{\partial v}{\partial x}dx$$

$${\mathbf{u}}^{T}MD\hat{\mathbf{v}}$$

$$P\hat{\mathbf{x}}=\mathbf{x},\hat{P}\mathbf{x}=\hat{\mathbf{x}}.$$

$${\mathbf{u}}^{T}MP\hat{\mathbf{v}}=(\hat{P}\mathbf{u}{)}^T\hat{M}\hat{\mathbf{v}}\Rightarrow MP={\hat{P}}^T\hat{M},$$

$$\mathbf{p}$$

$${\mathbf{v}}^1$$

$${\mathbf{v}}^2$$

$$\hat{N}=(N+2)(N+2),N_1=(N+1)(N+2),N_2=(N+2)(N+1),$$

$$D_1$$

$${\hat{D}}_1$$

$$H_1$$

$$\begin{array}{r}{P}_{1c}=(P\otimes \hat{I})\in {\mathbb{R}}^{N_1\times \hat{N}},P_{c1}=(\hat{P}\otimes \hat{I})\in {\mathbb{R}}^{\hat{N}\times N_1},\\ P_{2c}=(\hat{I}\otimes P)\in {\mathbb{R}}^{N_2\times \hat{N}},P_{c2}=(\hat{I}\otimes \hat{P})\in {\mathbb{R}}^{\hat{N}\times N_2}.\end{array}$$

$$H_1{P}_{1c}=P_{c1}^T\hat{H}\text{mbox}and{H}_2{P}_{2c}=P_{c2}^T\hat{H}.$$

Full text

1 Introduction

Many applications featuring wave propagation require numerical methods capable of treating complex geometry. A popular approach for finite difference methods is to solve the governing equations on curvilinear grids. Typically, the governing equations are first formulated in terms of Cartesian vector components followed by transforming the spatial derivatives with respect to the curvilinear coordinates. An alternative approach utilizes the invariance of the governing equation with respect to the choice of coordinate basis. In this paper we develop a staggered grid discretization of the acoustic wave equation where the components of the velocity vector are expressed with respect to the covariant basis, which varies throughout space.

Discretization of wave equations on general staggered curvilinear grids is challenging due to the emergence of off-diagonal metric terms. On a Cartesian grid, all fields and their components can be positioned in the grid so that all of the terms can be discretized by staggered difference approximations. However, a curvilinear coordinate transformation introduces additional terms into the governing equations that cannot be discretized without interpolation. One can directly discretize these additional terms by collocating some of the components, but this collocation increases the computational cost and memory storage. Despite this increase in computational cost, the approach has been used in computational seismology for solving the elastic wave equation [2, 34, 23]. An alternative solution is to solve the covariant form of the governing equations on orthogonal grids. In this case, the metric tensor becomes diagonal. This approach has been used in computational electrodynamics for solving Maxwell’s equations in covariant form [37, 33]. Unfortunately, the generation of orthogonal grids is a non-trivial task, in particular for many practical applications that feature complex 3D geometries. A third solution is to use interpolation in the off-diagonal components. So far this approach has only been developed with one-dimensional low order interpolation [9, 10, 11], resulting in a scheme that may be susceptible to long-term instabilities, requiring ad-hoc stabilization [8]. In the present work, we develop a provable stable numerical method that solves the acoustic wave equation in covariant form on a general staggered curvilinear grids. The proposed scheme is high-order-accurate, energy conservative and preserves positive definiteness of the discretized metric tensor.

To obtain a provably stable method, we discretize the acoustic wave equation by applying the Summation-By-Parts (SBP) principle. The SBP principle provides a way to analyze and derive provably stable schemes by constructing energy estimates for semi-discrete approximations. SBP methods were originally designed to obtain provably stable high order finite difference approximations for first order hyperbolic systems on collocated grids [13, 32, 20, 21]. Since that time, there have been numerous developments that, for example, extend the SBP methodology to non-conforming grids [16, 25, 3] and to equations with second and higher order derivatives [31, 26, 15]. SBP methods have also been extend beyond finite differences, e.g., to finite volume, discontinuous Galerkin, and flux reconstruction schemes [18, 5, 29]. Methods have also been devised for coupling different SBP discretizations [19, 14, 12, 4].

Our proposed curvilinear staggered SBP discretization builds upon our previous work on staggered grid methods for wave equations in Cartesian geometries [22, 28, 17]. Our main theoretical result is that stability can be guaranteed on a curvilinear grid under two conditions. First, the interpolation operators must be compatible with the SBP norm weights such that the scalar product between a staggered and a node-centered grid function gives the same result, regardless of which function is being interpolated. Secondly, the discretization must preserve the positive definiteness of the metric tensor. The latter condition is essential because the velocity field is decomposed with respect to the covariant basis. As a result, the symmetric and positive definite covariant metric tensor appears in the kinetic energy term. The discretization of this tensor involves interpolation, which can destroy positive definiteness. We show how to preserve positive definiteness of this tensor in the discretized equations. While it is possible to preserve positive definiteness for any non singular curvilinear grid, such a discretization has to involve interpolation in the diagonal components of the contravariant metric tensor, which reduces accuracy and increases computational cost. By modifying the metric tensor discretization in its diagonal components, it is possible to overcome this reduction in accuracy and increase in computational cost. However, positive definiteness can no longer be guaranteed for all non singular grids. We derive a condition that can be checked in a point-wise manner to ensure positive definiteness of the metric tensor discretization. This condition shows that the interpolation operators should be constructed such that their norm is close to one.

The remainder of the paper is organized as follows. The covariant form of the acoustic wave equation and curvilinear mappings are reviewed in Section 2. In Section 3 we develop SBP finite difference for staggered grids, first in one dimension and then in two dimensions for the acoustic wave equation. In Section 4, we show that the discretization is energy conservative. Two alternative ways of discretizing the metric tensor such that the positive definiteness is preserved are presented in Section 5. In Section 6, we conduct numerical experiments. We first investigate the accuracy of the scheme for the two different metric tensor discretizations. We then demonstrate that the covariant formulation can outperform a Cartesian formulation due to the loss of rotational invariance in the latter. The section is concluded by demonstrating that the solution is free from numerical artifacts when the acoutic wave propagation is driven by a discretized point force applied to the boundary. Finally, in Section 7, the study is summarized and conclusions are drawn.

2 Problem formulation

Consider the 2-D acoustic wave equation in dimensionless form (all quantities are scaled, including space and time)

[TABLE]

posed on a 2-D bounded domain $(x,y)\in\Omega$ for $t\geq 0$. Here, $p$ is the pressure and $\boldsymbol{v}$ is the velocity vector. We consider cases where the domain $\Omega$ can be defined through a curvilinear mapping from the unit square in parameter space [35]. We define the curvilinear coordinates

[TABLE]

and the continuously differentiable mapping $x=X(r^{1},r^{2})$ and $y=Y(r^{1},r^{2})$ between the curvilinear and Cartesian coordinates. We assume that the mapping is non singular. By differentiating the mapping with respect to the curvilinear coordinates, one obtains the covariant basis vectors,

[TABLE]

The corresponding covariant metric tensor is defined by $g_{ij}={\mathbf{a}}_{i}\cdot{\mathbf{a}}_{j}$. It is symmetric and positive definite.

The contravariant basis vectors ${\mathbf{a}}^{1}$ and ${\mathbf{a}}^{2}$ can be defined by the orthogonality condition

[TABLE]

The contravariant metric tensor satisfies $g^{ij}={\mathbf{a}}^{i}\cdot{\mathbf{a}}^{j}$. It is symmetric and positive definite, and the Jacobian of the curvilinear mapping, $J=1/\sqrt{|g^{ij}|}$, is assumed to be bounded. Here, $|g^{ij}|$ denotes the determinant of the tensor $g^{ij}$.

The velocity field is decomposed with respect to the covariant basis,

[TABLE]

where $v^{1}$ and $v^{2}$ are the contravariant velocity components. The contravariant metric tensor can be used to transform between covariant and contravariant velocity components,

[TABLE]

The inverse relationship is

[TABLE]

Contravariant quantities are always denoted by superscripts, whereas covariant quantities are always denoted by subscripts.

Next, we derive the covariant formulation of the acoustic wave equation. The divergence of the velocity vector field satisfies [35, 6]

[TABLE]

The pressure gradient satisfies

[TABLE]

By inserting (5) into (1) and (6) into (2), followed by taking the dot product of (2) with the contravariant basis, results in

[TABLE]

The total energy in the system is the sum of the acoustic and the kinetic energies

[TABLE]

The second integral expresses the kinetic energy as an invariant formed by contracting the covariant velocity components with the contravariant components. Alternately, the covariant velocity components can be transformed into the contravariant velocity components using (4), yielding

[TABLE]

In this form, we see that (10) is positive since the covariant metric tensor $g_{ij}$ is symmetric and positive definite. By differentiating the energy (10) with respect to time and inserting the governing equations (7)-(8), we obtain the energy rate

[TABLE]

The terms in brackets on the right hand side must be bounded by imposing one boundary condition per side of the unit square in parameter space. With the homogeneous pressure boundary condition $p=0$ on all sides,

[TABLE]

As expected, the acoustic wave equation in covariant formulation conserves the total energy in the system.

3 SBP operators

Before we can present the numerical scheme, we need to introduce some definitions. First, we review SBP staggered grid finite difference methods in one spatial dimension, and then in two spatial dimensions [22]. The SBP operators form building blocks for constructing our provably stable staggered grid scheme. To make it easier to understand how the construction works, we share some of our codes in the repository github.com/ooreilly/sbp. In particular, this repository contains the SBP operators derived in this work and a prototype implementation of the scheme.

We introduce the grid vectors $\mathbf{x}$ (node-centered) and $\hat{\mathbf{x}}$ (cell-centered) that discretize the unit interval using $N+1$ and $N+2$ grid points, respectively,

[TABLE]

The grid points are

[TABLE]

and $h=1/N$ is the grid spacing, see Figure 1. On the $\mathbf{x}$-grid, the function $u(x)$ is approximated by grid function $\mathbf{u}=[u_{0}\ u_{1}\ \ldots\ u_{N}]^{T}$, where $u_{i}\approx u(x_{i})$. Similarly, on the $\hat{\mathbf{x}}$-grid, the function ${v}(x)$ is approximated by the grid function $\hat{\mathbf{v}}=[\hat{v}_{0}\ \hat{v}_{1}\ \ldots\ \hat{v}_{N+1}]^{T}$, where $\hat{v}_{i}\approx v(\hat{x}_{i})$.

We introduce the SBP staggered grid difference operators $D$ and $\hat{D}$ that approximate the first derivative, see Figure 1. These operators are accurate to order $2s$ in the interior and to order $s$ for a few points near each boundary. The operator $D$ approximates the derivative of a grid function on the $\hat{\mathbf{x}}$-grid at the $\mathbf{x}$-grid. In contrast, the operator $\hat{D}$ approximates the derivative of a grid function on the $\mathbf{x}$-grid at the $\hat{\mathbf{x}}$-grid. For example, for all difference operators that are at least second order accurate on the boundary ($s\geq 2$), quadratic monomials $\mathbf{\xi}$ and $\hat{{\boldsymbol{\xi}}}$ are differentiated exactly,

[TABLE]

Staggered SBP difference operators satisfy a summation by parts property, corresponding to integration by parts,

[TABLE]

The property (13) is key to proving energy stability of the semi-discrete approximation. The matrices $M$ and $\hat{M}$ are positive definite diagonal matrices that define $2s-1$ order accurate quadrature rules on the grids $\bf x$ and $\hat{\mathbf{x}}$, respectively. Since the grid points of the cell-centered and nodal grids overlap, the SBP property (13) gives boundary terms involving unknowns on the boundary. Alternate approaches that avoid overlapping boundary points either interpolates or extrapolates the numerical solution to the boundary [3], or strongly impose the boundary condition [28].

In addition to difference operators, we also need to define SBP staggered interpolation operators. These operators are $2s$ order accurate in the interior and $s$ order accurate on the boundary. The interpolation operator $P$ interpolates a grid function on the $\hat{\mathbf{x}}$-grid to the $\mathbf{x}$-grid, and the interpolation operator $\hat{P}$ interpolates in the opposite direction. For example, if the interpolation operator is at least second order accurate on the boundary ($s\geq 2$), then

[TABLE]

The interpolation operators satisfy the SBP property

[TABLE]

for all grid functions ${\mathbf{u}}$ and ${\mathbf{v}}$. The property (14) states that when evaluating the integral $\int uvdx$ using the SBP quadrature, the result must be the same whether ${\mathbf{u}}$ is interpolated to the $\hat{\mathbf{x}}$-grid, or $\hat{{\mathbf{v}}}$ is interpolated to the ${\mathbf{x}}$-grid.

3.1 Two-dimensional SBP Operators

In two spatial dimensions, we introduce a staggered grid that discretizes the pressure field $p$ and the contravariant velocity components ($v^{1}$, $v^{2}$). This grid is parameterized by the curvilinear coordinates $r^{1}$ and $r^{2}$, see Figure 2.

On the staggered grid, pressure $p$ is located at cell-centers and the contravariant velocity components $v^{1}$ and $v^{2}$ are located on the edges with $r^{1}=\mbox{const.}$ or $r^{2}=\mbox{const.}$, respectively.

Without loss of generality, we use the same number of cells $N$ (and grid spacing $h$) in both coordinate directions. We collect the spatial discretizations of the pressure and velocity component fields in the vectors $\mathbf{p}$, $\mathbf{v}^{1}$ and $\mathbf{v}^{2}$. These column vectors are ordered as follows:

[TABLE]

where $p_{ij}\approx p(\hat{r}^{1}_{i},\hat{r}^{2}_{j},t)$, $v^{1}_{ij}\approx v^{1}(r^{1}_{i},\hat{r}^{2}_{j},t)$, and $v^{2}_{ij}\approx v^{2}(\hat{r}^{1}_{i},r^{2}_{j},t)$. Furthermore, the length of each grid vector is, respectively,

[TABLE]

i.e., $N_{1}=N_{2}$ in this particular case.

We extend the differentiation and interpolation operators to two spatial dimensions by applying them in a line-by-line manner using Kronecker products. The difference operators become

[TABLE]

where $\hat{I}$ is the $\hat{N}\times\hat{N}$ identity matrix. In this notation, $D_{1}$ is a difference operator that acts on a cell-centered quantity and approximates the first derivative $\partial/\partial r^{1}$ on the edge-1 grid, see Figure 3. We also introduce the following norm weight matrices,

[TABLE]

and interpolation operators

[TABLE]

The interpolation operator $P_{1c}$ interpolates from the cell-centered grid to the edge-1 grid, and so forth. By applying the SBP interpolation property (14) to (15), we find

[TABLE]

The extension of the SBP property (13) to 2D becomes,

[TABLE]

where ${\mathbf{p}}_{R},\ {\mathbf{p}}_{L},\ {\mathbf{p}}_{T},\ {\mathbf{p}}_{B}\in\mathbb{R}^{N+2}$, contain the grid values of ${\mathbf{p}}$ along the right, left, top, and bottom boundaries, respectively (and similarly for ${\mathbf{v}}^{1}$ and ${\mathbf{v}}^{2}$).

3.2 The semi-discrete approximation

After discretizing the spatial derivatives on the staggered grid, the semi-discretization of the pressure equation (7) becomes

[TABLE]

In (18), $\hat{J}$, $J_{1}$, and $J_{2}$ are diagonal matrices that hold the Jacobian evaluated at the cell-centers and the edge grid points, respectively.

The semi-discrete approximation of the equations for the contravariant velocity components (8) become

[TABLE]

In (19), $G$ is the discrete approximation of the contravariant metric tensor $g^{ij}$. The explicit form of $G$ will be presented in the next section. Recall that the contravariant metric tensor gives the transformation from covariant to contravariant components. The discrete counterpart of (3) is

[TABLE]

Since the contravariant components ${\mathbf{v}}^{1}$ and ${\mathbf{v}}^{2}$ are associated with different locations of the staggered grid, $G$ must be constructed using interpolation, unless the curvilinear grid is orthogonal.

4 The energy method

In this section, we apply the energy method to show that the scheme (18)-(19) is stable. To simplify the presentation, we consider periodic boundary conditions and refer to Appendix A for a stability proof with weakly imposed boundary conditions using the SAT penalty method.

In the following, the matrix $G^{-1}$ is obtained by inverting $G$ in (19). The matrix $G^{-1}$ approximates the covariant metric tensor $g_{ij}$. In practice, we never have to explicitly compute the inverse. Consider the following discrete approximation of the energy (10), $\mathcal{E}=\mathcal{E}_{a}+\mathcal{E}_{k}$,

[TABLE]

In (21), we have defined

[TABLE]

Recall that the first term on the right-hand side of (10) is the acoustic energy and the two remaining terms represent the kinetic energy. The discrete approximation of the acoustic energy, $\mathcal{E}_{a}$, defines a norm of ${\mathbf{p}}$ since the matrices $\hat{H}$ and $\hat{J}$ are diagonal with positive entries. In order for the kinetic energy, $\mathcal{E}_{k}$, to be a norm, the matrix $HJG^{-1}$ must be positive definite. Since $H_{1}$, $H_{2}$, $J_{1}$, and $J_{2}$ are positive diagonal matrices, $HJ$ is a positive diagonal matrix. The construction of $G$ such that $HJG^{-1}$ is symmetric and positive definite is discussed in the next section. Given such a $G$, we obtain an energy stable scheme, as is shown in the following lemma.

Lemma 1**.**

Assume that the grid metric and solution is periodic in each direction. Then the energy (21) of the semi-discrete approximation (18)-(19) is conserved.

Proof.

Differentiating the acoustic energy in (21) with respect to time and inserting (18) into the result yields

[TABLE]

Similarly, differentiating the kinetic energy in (21) with respect to time and inserting (19), and assuming $HJG^{-1}=(HJG^{-1})^{T}$, yields

[TABLE]

The energy rate is the sum of (23) and (24),

[TABLE]

On a two dimensional periodic grid, the boundary terms in (17) cancel out and therefore the SBP property simplifies to $\widehat{H}\widehat{D}_{1}=-D_{1}^{T}H_{1}$ and $\widehat{H}\widehat{D}_{2}=-D_{2}^{T}H_{2}$. Thus, the right-hand side of (25) is identically zero and the energy of the semi-discrete approximation is conserved. ∎

5 Discretization of the contravariant metric tensor

To define a norm for $({\mathbf{v}}^{1},{\mathbf{v}}^{2})$, note that the second matrix in (21) can be written as

[TABLE]

Because $HJ$ is diagonal and positive definite, it is sufficient to show that $HJG$ is symmetric and positive definite.

There are many ways to construct the matrix $G$ to form the approximation $\mathcal{E}_{k}$ of the kinetic energy in (21), but care is needed to preserve energy positivity. What complicates matters is the fact that ${\mathbf{v}}^{1}$ and ${\mathbf{v}}^{2}$ are defined on different grids. Therefore, the construction of $G$ must involve interpolation operators at least in its off-diagonal blocks. However, as we will see in Section 5.2, positive definiteness may be lost when only the off-diagonal blocks are interpolated.

5.1 Unconditionally energy positivity preserving construction

To obtain a matrix $G$ that gives a positive definite $HJG$, we evaluate the contravariant metric tensor and Jacobian at the cell-centered grid points, and then interpolate back and forth to each velocity grid. This approach results in

[TABLE]

The reason why $J$ appears in the construction of $G$ will become clear next. To show the following result, we first need to make an assumption.

Assumption 1**.**

The null-space of the SBP interpolation operators in (15) is empty.

We have verified this assumption through numerical experiments. We are now ready to state the following result.

Lemma 2**.**

Let $HJ$ be given by (21) and let $G$ be defined by (26). The matrix $HJG$ is positive definite if Assumption 1 holds, $\hat{J}>0$ and $\hat{g}^{ij}$ is positive definite.

Proof.

Multiplying (26) by $HJ$ results in

[TABLE]

Using the SBP interpolation property (16),

[TABLE]

By factoring out the interpolation and norm operators as well as the Jacobians, we find

[TABLE]

This factorization relies on the fact that the matrices $\hat{J}$, $\hat{H}$, $\hat{g}^{11}$, $\hat{g}^{12}$, and $\hat{g}^{22}$ are diagonal with positive entries, and hence commute. Next, we show that the matrix $V$ is positive definite. We can permute $V$ so that it consists of positive definite $2\times 2$ blocks, one per grid point. Let $R$ be a permutation matrix that is obtained by permuting the rows of the identity matrix. Then $\tilde{V}=R^{T}VR$ is the permutation of $V$ that consists of $2\times 2$ positive definite blocks. Each block holds the contravariant metric tensor corresponding to one grid point, which is non singular by assumption. Since permutation matrices satisfy $RR^{T}=I$, we have that $V=RR^{T}VRR^{T}=R\tilde{V}R^{T}$ is positive definite. If Assumption 1 holds, then $HJG$ is symmetric and positive definite. ∎

5.2 Conditionally energy positivity preserving

construction

While the particular choice of $G$ in (26) yields a positive kinetic energy for all non singular mappings, it may be undesirable to use this choice in practice. The reason is that blocks along the diagonal of $G$ contain interpolation operators that interpolate back and forth between two different grids. This interpolation adds extra computation and may also reduce the accuracy in the numerical solution. Next, we show that it is possible to avoid these extra computations and improve the accuracy. However, care has to be taken to ensure that the contravariant metric tensor is discretized such that it remains positive definite.

Consider the following discretization,

[TABLE]

This choice preserves the symmetry property $HJ\widetilde{G}=(HJ\widetilde{G})^{T}$, but may cause a loss of positive definiteness. To avoid having to perform expensive eigenvalue computations to assert whether $HJ\widetilde{G}$ is positive definite, we derive and analyze an approximate condition that can be evaluated in a point-wise manner.

Introduce the coefficients $\alpha>0$ and $\beta>0$. We use the SBP interpolation relations (16) to decompose the matrix (28) according to

[TABLE]

where

[TABLE]

and

[TABLE]

We can guarantee that $HJ\widetilde{G}$ is positive definite if the matrix $A(\alpha,\beta)$ is positive definite and $B(\alpha,\beta)$ is positive semi-definite. The matrix $A(\alpha,\beta)$ becomes indefinite if $\alpha<=0$ or $\beta<=0$. We must therefore assume that the coefficients $\alpha$ and $\beta$ are positive. On the other hand, the matrix $B(\alpha,\beta)$ becomes indefinite if $\alpha$ or $\beta$ are too large. Our strategy is therefore to first derive upper bounds on $\alpha$ and $\beta$ that guarantee that $B(\alpha,\beta)$ is positive semi-definite. The definiteness of $A(\alpha,\beta)$ can then be verified by solving a $2\times 2$ eigenvalue problem at each cell center.

The matrix $B(\alpha,\beta)$ in (32) is clearly positive semi-definite if $B_{11}(\alpha)$ and $B_{22}(\beta)$ are positive semi-definite, which can be guaranteed by bounding $\alpha$ and $\beta$ as is shown in the following lemma.

Lemma 3**.**

The matrices $B_{11}(\alpha)$ and $B_{22}(\beta)$ are positive definite if the coefficients $\alpha$ and $\beta$ satisfy the bounds

[TABLE]

Here, $\max\lambda^{(1)}$ and $\max\lambda^{(2)}$ are the largest eigenvalues of the respective generalized eigenvalue problems,

[TABLE]

Proof.

Both bounds can be derived in the same fashion and we focus on the bound for $\alpha$. From the definitions above,

[TABLE]

The matrices ${\mathbb{X}}^{(1)}$ and ${\mathbb{Y}}^{(1)}$ are real, symmetric and positive definite. All eigenvalues $\lambda^{(1)}_{j}$ of the generalized eigenvalue problem (34) are therefore real and positive. Furthermore, the corresponding eigenvectors $\boldsymbol{y}_{j}$ form a complete set and can be normalized such that

[TABLE]

Any vector $\boldsymbol{z}$ can therefore be expanded according to

[TABLE]

Thus,

[TABLE]

The matrix $B_{11}(\alpha)$ is therefore positive semi-definite if

[TABLE]

The inequality must be satisfied for all eigenvalues $\lambda^{(1)}_{j}$. We have $1/\max(\lambda^{(1)})\leq 1/\lambda^{(1)}_{j}\leq 1/\min(\lambda^{(1)})$, which shows that $B_{11}(\alpha)$ is positive semi-definite as long as

[TABLE]

∎

In Appendix C we show how the above 2-D generalized eigenvalue problems decouples into a number of regular eigenvalue problems along each 1-D grid line.

When using the modified metric tensor discretization, the norm $\|P\hat{P}\|_{2}$ should be close to one to avoid instability. During the construction of the SBP interpolation operators, there can be undetermined coefficients remaining after satisfying the SBP property (14) and accuracy conditions. These coefficients can be determined by choosing an objective function that minimizes $\|P\hat{P}\|_{2}$. In this work we optimize for accuracy and manually check that $\|P\hat{P}\|_{2}\approx 1$. Please see [22, 28] for further construction details.

To understand why $\|P\hat{P}\|_{2}\approx 1$ is desirable, assume the mapping is linear. Then the metrics $J$ and $g^{ij}$ are constant, and the generalized eigenvalue problem (34) simplifies to

[TABLE]

Thereafter, by applying the SBP property (16), we get

[TABLE]

By observing that $P_{1c}=(P\otimes\hat{I})$ and $P_{c1}=(\hat{P}\otimes\hat{I})$, we obtain an eigenvalue problem along a single grid line in the $r^{1}$-direction, $P\hat{P}\mathbf{x}=\lambda^{(1)}\mathbf{x}$. The maximum eigenvalue $\max(\lambda^{(1)})$ satisfies $1\leq\max(\lambda^{(1)})\leq\|P\hat{P}\|_{2}$. The lower bound comes from observing that the constant vector ${\mathbf{1}}=[1\ 1\ \ldots\ 1]^{T}$ is interpolated exactly, and is therefore an eigenvector of $P\hat{P}$.

To investigate how the size of $\|P\hat{P}\|_{2}$ influences the stability of the scheme, we construct one pair of interpolation operators with $\|P\hat{P}\|_{2}=1.04$ and another pair with $\|P\hat{P}\|_{2}=8.53$. For the test, we consider an almost square domain where three sides are straight and the fourth side is a Gaussian hill,

[TABLE]

We increase the amplitude $\gamma$ of the Gaussian hill and monitor the minimum eigenvalue of the symmetric matrix $HJ\widetilde{G}$ in the kinetic energy. This matrix must be positive definite to guarantee stability. We expect that for a sufficiently large amplitude of the Gaussian hill the minimum eigenvalue of $HJ\widetilde{G}$ becomes negative. In the test, the grids are coarse with $17\times 17$ grid points. Figure 4 show the grids at which $HJ\widetilde{G}$ become positive semi-definite for each set of interpolation operators. As can been seen, if $\|P\hat{P}\|_{2}\approx 1$, the grid can undergo significant deformation before $HJ\widetilde{G}$ becomes positive semi-definite (Figure 4a). On the other hand, if $\|P\hat{P}\|_{2}\gg 1$, the grid can only undergo small deformations before $HJ\widetilde{G}$ becomes positive semi-definite (Figure 4b).

Since it is computationally expensive to compute the minimum eigenvalue $\min(\lambda)$ of $HJ\tilde{G}$, we demonstrate how to estimate $\min(\lambda)$ using Lemma 3. Note that the constants $\alpha$ and $\beta$ can be bounded using Lemma 3 without much computational work because these bounds only involve solving 1-D eigenvalue problems. Once $\alpha$ and $\beta$ are bounded, the matrix $C(\alpha,\beta)$ in (29) can be evaluated in a point-wise manner. This approach provides an efficient way of estimating whether $HJ\widetilde{G}$ is positive definite. If $C(\alpha,\beta)$ is positive definite at all grid points, then $HJ\widetilde{G}$ is positive definite. However, if $C(\alpha,\beta)$ has a negative eigenvalue at any grid point, then the test for definiteness of $HJ\widetilde{G}$ is inconclusive.

Figure 5 shows how accurate the approximate method is at estimating if $HJ\widetilde{G}$ is positive definite by computing the minimum eigenvalue of $C(\alpha,\beta)$ for all grid points. In Figure 5a, the approximate method estimates that $HJ\widetilde{G}$ is positive definite for $\gamma<0.5$. (the actual amplitude at which $HJ\widetilde{G}$ becomes positive semi-definite is $\gamma\approx 2.5$). In Figure 5b, the approximate method estimates that $HJ\widetilde{G}$ is positive definite for amplitudes $\gamma<0.02$ (the actual amplitude at which $HJ\widetilde{G}$ becomes positive semi-definite is $\gamma\approx 0.07$).

In conclusion, large values of $\|P\hat{P}\|_{2}$ can severely limit the amount of grid skewness that the scheme can tolerate without becoming unstable. This restriction may be so severe that stable computations on realistic grids become impossible. Therefore, if one wants to use the approximation (28) of the contravariant metric tensor, it is important to construct interpolation operators such that $\|P\hat{P}\|_{2}\approx 1$.

6 Numerical experiments

We use the method of manufactured solutions to conduct convergence rate studies [30]. The following manufactured solution is designed to satisfy the acoustic wave equation (1)-(2) in Cartesian velocity components ($\boldsymbol{v}=\tilde{v}_{x}\boldsymbol{e}_{x}+\tilde{v}_{y}\boldsymbol{e}_{y}$),

[TABLE]

To apply the manufactured solution, we impose the inhomogeneous pressure boundary condition, $p(x,y,t)=\tilde{p}(x,y,t)$ for all points $(x,y)$ on the boundary. We use the classical fourth order Runge-Kutta scheme to advance the solution in time. Since a naive implementation of time dependent boundary data reduces the convergence rate to second order for Runge-Kutta schemes, we apply the boundary correction procedure described in [1].

The error is the difference between the numerical and manufactured solution, and it is measured using either the discrete $l_{2}$ and $l_{\infty}$ (max) norm for each field. The $l_{2}$ errors are

[TABLE]

where we used grid function notation $p_{ij}\approx p(\hat{x}_{i},\hat{y}_{j},t)$ etc. The sum of errors for all fields is denoted by

[TABLE]

where $l=h,\infty$ (depending on the choice of norm). Let $q_{h}$ and $q_{\infty}$ denote the convergence rates in the $l_{2}$ and $l_{\infty}$ norms, respectively. In (37), the contravariant components of the manufactured solution are obtained by the formulas

[TABLE]

where ${\mathbf{e}}_{x}=(1,0)$ and ${\mathbf{e}}_{y}=(0,1)$ are the Cartesian basis vectors. In all numerical experiments, we use staggered SBP interpolation and differentiation operators having fourth order interior accuracy and second order boundary accuracy. To construct the metrics in the scheme, we use the same staggered SBP interpolation and differentiation operators for discretizing the covariant basis vectors ${\mathbf{a}}_{1}$ and ${\mathbf{a}}_{2}$, Jacobian $J$. The discretized contravariant basis vectors ${\mathbf{a}}^{1}$ and ${\mathbf{a}}^{2}$ are obtained from (2).

6.1 Comparison of metric tensors

As discussed in section 5, the discretization of the contravariant metric tensor $G$ in (26) results in a provably stable scheme. However, as we proceed to show, the modified discretization $\widetilde{G}$ in (28) gives better accuracy.

For this test, we use linear transfinite interpolation [35] to define a mapping by specifying the following parameterization of the four boundary segments

[TABLE]

where $a=0.05$, $k=2\pi$, $x_{0}=y_{0}=0.2$.

Table 1 lists the errors and convergence rates for the energy positivity preserving discretization of $G$ in (26) and the modified $\widetilde{G}$ in (28). In most cases, the error is smaller for the scheme using $\widetilde{G}$ compared to $G$. Note that the convergence rates for the velocity components on the finer grids seem to drop to order $2.5$ in $l_{2}$.

6.2 Characteristic boundary error

In this section, we continue analyzing the error from the experiment in Section 6.1 with the modified metric tensor discretization (28). In particular, we observed that the convergence rate dropped to order $2.5$ in the $l_{2}$ norm on the finer grids. This convergence rate is consistent with the truncation error analysis using the energy method, see [36], and generalizes for higher order operators as $p+1/2$ where $p$ is the order of the truncation error of the SBP operator near the boundary. According to the general theory for SBP finite difference approximations of first order hyperbolic problems [7], the convergence rate is in many cases one order higher than the order of the truncation error at the boundary, measured in the $l_{2}$ norm. However, as is pointed out in [27], the general theory assumes a non-characteristic boundary. Since the acoustic wave equation contains zero-speed characteristic variables, the general theory does not apply.

In order to examine the error in more detail, we symmetrize and diagonalize the acoustic wave equation with respect to the direction normal to the boundary, see Appendix B. For the bottom boundary, we obtain the characteristic variables

[TABLE]

where $\eta=\sqrt{g_{11}g_{22}-g_{12}^{2}}$. These variables correspond to the eigenvalues

[TABLE]

The characteristic variable $w^{+}$ corresponds to the positive eigenvalue and is ”incoming”, meaning that it propagates into the domain, whereas $w^{-}$ corresponds to the negative eigenvalue and is ”outgoing”, implying that it propagates out of the domain. However, the eigenvalue corresponding to the variable $w^{0}$ is zero, making it boundary characteristic, meaning that it does not propagate in the direction normal to the boundary. It is the presence of the boundary characteristic variable that violates the assumptions of the general theory [7].

To examine the error in the characteristic variables, we focus on a cross-section that starts in the middle of the bottom boundary ($r^{1}=0.5$) and extends $10\%$ into the domain, i.e., $0\leq r^{2}\leq 0.1$. We use the SBP interpolation operators to compute the incoming and outgoing characteristic variables at the cell-centers. Note that in (42), only the boundary characteristic variable $w^{0}$ depends on $v^{1}$. For simplicity, we therefore define

[TABLE]

as the norm of the error in the characteristic and non-characteristic boundary variables, respectively.

Figure 6 shows the absolute value of the solution error at these grid points for the two finest grids in Table 1. Note that the solution error at the first four grid points is much larger than the error at the interior grid points. These four grid points correspond to the region in which the SBP operators have a second order accurate truncation error. Table 2 lists the errors in the $l_{2}$ norm and $l_{\infty}$ norm measured at these grid points for each grid. The convergence rates of the non-characteristic variables show close to third order convergence in $l_{2}$, which is consistent with the standard theory [7]. The convergence rate of the boundary characteristic variable is reduced compared to the non-characteristic boundary variables. This drop in the convergence rate can also be seen in Figure 6.

6.3 Rotational invariance

A desirable property of the numerical scheme is that its accuracy and stability should not depend on the orientation of the coordinate system. Since we have chosen to decompose the velocity field with respect to the covariant basis, the discretization is local with respect to the orientation of the curved grid lines, and therefore rotationally invariant. If we instead the decompose the velocity field with respect to the Cartesian basis, the discretization does not preserve rotational invariance. In this section we compare the accuracy of these two approaches.

When the velocity is decomposed with respect to the Cartesian basis, the metric tensor that transforms from Cartesian components to contravariant components is defined by

[TABLE]

This transformation is explicitly given by

[TABLE]

When a grid is rotated with respect to the Cartesian basis vectors such that the vector pairs $(\mathbf{{\mathbf{e}}}_{x},\ {{\mathbf{a}}}^{1})$ and $(\mathbf{{\mathbf{e}}}_{y},\ {\mathbf{a}}^{2})$ become orthogonal, the diagonal entries in $A$ vanish and only off-diagonal entries remain. These off-diagonal entries must be treated using interpolation. It is therefore to be expected that the Cartesian formulation will suffer from accuracy degradation for such grids.

Before proceeding to the numerical experiment, we give some more details. When transformed to the Cartesian basis, the governing equations become

[TABLE]

The semi-discrete approximation (excluding SAT penalty terms) becomes

[TABLE]

When discretizing the block matrix $\tilde{A}$ in (46), it is necessary to use interpolation operators for the off-diagonal entries,

[TABLE]

where the diagonal matrix $\tilde{A}_{ij,k}$ is obtained by discretizing entry $A_{ij}$ in $A$ on the edge-$k$ grid. Here, the interpolation operators are $P_{12}=(P\otimes\hat{P})$ and $P_{21}=(\hat{P}\otimes P)$.

When the velocity field is decomposed with respect to the Cartesian basis, the mixed terms that arise in the acoustic wave equation depend on the angle between the contravariant basis and Cartesian basis vectors. Since the mixed terms must be discretized by a combination of interpolation and differentiation operators, the resulting operators become wide. These operators may be less accurate than the staggered difference operators used for the diagonal terms.

To test how the accuracy depends on these angles, we use the manufactured solution (37) for the following disc geometry with a cavity,

[TABLE]

where the radial coordinate is stretched using linear interpolation from the inner radius $R_{0}$ to the outer radius $R_{1}$ by the function

[TABLE]

The stretching parameter $a$ is taken to be $a=4\pi/n_{2}$ ($n_{2}$ is the number of cells in the angular, $r^{2}$, direction). The inner and outer radii are set to $R_{0}=0.3$ and $R_{1}=1$. Periodic boundary conditions are used in the angular direction and an offset $\phi=0.2\pi$ is used to avoid aligning the periodic boundary with any of the Cartesian axes.

On the coarsest grid we use $n_{1}=16$ grid cells in the radial direction and $n_{2}=48$ grid cells in the angular direction. Finer grids are obtained by a factor of two grid refinement at each level. The solution is advanced in time until $T=0.5$ using the time step $\Delta t=0.015625$. The results in Table 3 show that the scheme using the covariant representation of the velocity field is always more accurate compared to the Cartesian basis decomposition.

6.4 Point source on the boundary

In computational acoustics, point sources are frequently used to model explosions. Spurious oscillation can arise when source terms are introduced in schemes discretized by central finite differences. The triggering of spurious oscillations can be avoided by using Cartesian staggered grids, or by smoothing out the source discretization by imposing smoothness conditions, see [24]. However, there is no guarantee that spurious oscillations are avoided in our curvilinear staggered scheme.

The point source we consider is placed on the boundary $r^{2}=1$ by enforcing the pressure boundary condition

[TABLE]

In the above, $\delta(r^{1}-r_{*})$ is the Dirac distribution, centered at the source location $r_{*}$, and $s(t)$ is a given source time function. For a detailed discussion on how to discretize the Dirac distribution for hyperbolic problems, see [24]. To place the source on the boundary, we treat (48) as a boundary condition and impose it using a SAT term, see [22] for details.

In this test, only the top boundary is curved, and its shape is given by the Gaussian hill,

[TABLE]

where $\sigma=0.105$. A 1-D coordinate stretching is defined to generate the mapping. The source is placed on the Gaussian hill at the location ${\mathbf{r}}_{s}=(0.45,1)$ in parameter space, or $\mathbf{x}_{s}=(4.5,5.79711)$ in physical space, see Figure 7. The source time function is a Ricker wavelet, i.e.,

[TABLE]

where the time delay is $t_{0}=1.7$, to prevent an abrupt startup. We define the minimum wavelength that must be resolved as $\lambda_{min}=c_{min}/f_{max}$, where $c_{min}=1$ is the minimum wave speed (here set to 1 due to solving the acoustic wave equation in dimensionless form) and $f_{max}=2.5$. This estimate corresponds to the frequency where the Ricker wavelet reaches $5\%$ of the peak amplitude in the frequency domain.

Since we do not have an analytic solution for this test, we compare the numerical solutions against a reference solution that has been computed on a very fine mesh. When estimating convergence rates, we take into account that the reference solution itself contains an error. The resolution of the reference solution is $4096\times 2048$ grid points, which corresponds to $N_{\lambda}=\lambda_{min}/h\approx 160$ grid points per shortest wavelength. The error is defined as the difference between the reference solution and a numerical solution measured in time at the interior receiver location ${\mathbf{r}}_{r}=(0.5,0.5)$, corresponding to $\mathbf{x}_{r}=(5,3)$. Since this point does not align with a grid point for each field, we use cubic interpolation of the neighboring values to estimate the solution.

We advance the solution in time using the time step $\Delta t=0.03125$ (for the coarsest grid) until the final time $T=7.8125$. Figure 7 shows that there are no visible artifacts from placing a point source on the boundary when discretizing using the modified metric tensor (28). On the other hand, the provably stable metric tensor discretization (26) causes spurious oscillations. For the same simulation, we also show the pressure as function of time at the receiver location and the estimated error (Figure 8). The initial oscillations correspond to the direct arrival of the pressure wave from the source; oscillations at later times are due to reflections against the non-planar boundary.

To investigate the error in more detail, let $\|\Delta p\|_{{\mathbf{r}}}$ denote the error in pressure measured in the $l_{2}$ norm over time at the fixed point ${\mathbf{r}}$, and let the two velocity components be defined in a similar manner. Table 4 lists the relative error and convergence rates for both pressure and velocity as well as an estimate of the number of grid points per shortest wavelength to obtain the desired error. We note that fourth order convergence rates are obtained for all components. In this case, we see no indication that the error due to the characteristic boundary variable has polluted the solution at this interior location.

6.5 Computational Efficiency

After having established that the staggered scheme with the modified metric tensor is capable of handling point sources on the boundary, we now investigate its computational performance and compare it to a collocated finite difference method. This test case is quite similar to the previous test, and can be reproduced by using the code found in the repository github.com/ooreilly/sbp/.

In the collocated scheme the governing equations are discretized using the classical fourth order SBP operators derived by B. Strand in [32]. One can convert the staggered scheme into a collocated scheme by replacing the interpolation operators with identity matrices and the derivatives with Strand’s SBP operators.

The primary difference in terms of computational work between the collocated and staggered scheme stems from the use of interpolation. In the interior of the staggered scheme, the use of interpolation operators in the velocity equations introduces two stencils with $4\times 7$ points, whereas the remaining terms involve only four stencils with $4$ points each. In contrast, the entire collocated scheme consists of only six stencils with 4 points each. Despite these large differences in computational work, the staggered scheme can outperform the collocated scheme when running on modern hardware.

The performance of each scheme can be highly dependent on the implementation. For prototyping purposes, it can be quite convenient to implement the schemes as sparse matrix and dense vector products. In this way, the implementation looks similar to the discretizations presented in Section 3. Unfortunately, this approach results in excessive computational overhead due explicitly storing the large sparse matrix that holds the spatial discretization. Instead, significant performance improvements can be achieved if the implementation is matrix-free, or stencil-based.

We measure the performance of the collocated and staggered scheme implemented in a stencil-based manner. Both these compute kernels are executed on a Nvidia Geforce RTX 2080 Ti card in single precision. In our timings, we only investigate the performance of the interior computation (in most practical cases, the boundary only consitutes a small fraction of the total cost). No particular effort has been invested into optimizing the kernels performance. Both kernels are memory bound, achieving $75\%$ of peak memory bandwidth (staggered) and $85\%$ of peak memory bandwidth (collocated). Compute utilization of the GPU streaming multiprocessors (SM)s reaches $45\%$ (staggered) and $25\%$ (collocated). This performance analysis was conducted using a grid size of 1024 x 512 grid points.

The solution is advanced in time using the same settings as before. Before the final time is reached, acoustic waves have reflected against the exterior boundaries and these reflections have reached the receiver location. As can be seen in Figure 9, the staggered scheme outperforms the collocated scheme on all grids tested, and its performance improves as the error levels decrease. On the finest grids, the difference in relative error approaches an order of magnitude.

7 Conclusions

We have shown how to construct a provably stable high order accurate finite difference approximation of the acoustic wave equation in covariant form on general curvilinear staggered grids. To obtain an energy conserving scheme, the discretization uses SBP interpolation and differentiation operators and imposes the boundary conditions weakly. After transforming the acoustic wave equation to curvilinear coordinates, the covariant metric tensor emerges in the kinetic energy. To make the scheme stable, the discretized kinetic energy must form a discrete norm of the dependent variables. The weights in the norm are represented by a symmetric matrix that depends on the discretization of the metric tensor. Stability follows if this matrix is positive definite. We offer two alternative ways of making this matrix positive definite. The first option uses interpolation of both the diagonal and off-diagonal terms in the metric tensor. While this approach leads to a stable scheme for all non singular meshes, it results in a slight degradation in accuracy and is more computationally expensive compared to our second alternative. In that approach, a modified discretization of the metric tensor is derived where only the off-diagonal elements are interpolated. Here, positive definiteness of the norm matrix depends on the mesh quality. Estimates are theoretically derived to efficiently determine if the discretization is stable for a given mesh. Numerical experiments illustrate the benefits of using interpolation operators with norms close to one and demonstrate the accuracy and stability of the proposed methodology.

While this work has focused on the acoustic wave equation in covariant form, future extensions to Maxwell’s equations or the elastic wave equation should also be possible. However, these extensions might be non-trivial and warrant further investigations. For Maxwell’s equations one must contend with the discretization of the curl operator in a stable and consistent manner. Some additional difficulties with solving the elastic wave equation in covariant form are the emergence of Christoffel symbols, and transformation of the stress tensor from covariant to contravariant form.

Acknowledgment

We would like to thank Dr. Martin Almquist for many fruitful discussions during the initial developments of this work, as well as his insightful comments and feedback on an early version of this manuscript.

This research was supported by the Southern California Earthquake Center (Contribution No. 9223), with contributions from NSF Cooperative Agreements EAR-1600087, USGS Cooperative Agreement G17AC00047, NSF-ACI Award 1450451, and Pacific Gas and Electric.

This work performed under the auspices of the U.S. Department of Energy by Lawrence Livermore National Laboratory under Contract DE-AC52-07NA27344. This is contribution LLNL-JRNL-780017.

Appendix A Weakly imposed pressure boundary condition

Here, we demonstrate how to weakly impose the homogeneous pressure boundary condition method in the covariant formulation of our scheme. To simplify the analysis, we will only consider the left boundary $r^{1}=0$ (the other boundaries are treated in a analogous manner).

When the grid is not periodic, the SBP property (17) applied to the right-hand side of the energy rate (25) leads to

[TABLE]

where $\hat{J}_{L}$ contains the values of the Jacobian along the edge-1 grid on left boundary, and the contributions from the other boundaries have been neglected. Unless a boundary condition is specified, the energy rate is indefinite. For simplicity, we impose the homogeneous boundary condition $p=0$ using the Simultaneous Approximation Term (SAT) penalty method.

Consider the velocity equations (19), augmented with a SAT penalty term on the right-hand side,

[TABLE]

The penalty vector $\mathbf{S}$ restricts the pressure field to the left boundary, and can be written as

[TABLE]

where the matrix $L$ is defined such that $L^{T}{\mathbf{v}}^{1}={\mathbf{v}}^{1}_{L}$, and ${\mathbf{f}}_{L}(t)$ is the boundary data. Since we focus on homogeneous boundary conditions, we set ${\mathbf{f}}_{L}=0$. When computing the energy rate including the penalty vector $\mathbf{S}$ in (49), we find

[TABLE]

By repeating the same procedure for the remaining boundaries, the energy of the semi-discrete approximation is conserved.

Appendix B Symmetrization and diagonalization of the governing equations

We show that the 2D acoustic wave equation in covariant form (7)-(8) can be written as a symmetric hyperbolic system. Recall the governing equations,

[TABLE]

Multiplying (51) by the covariant metric tensor gives the system

[TABLE]

This system can be written in matrix form as

[TABLE]

where

[TABLE]

The matrices $A$ and $B$ are symmetric and $W=W^{T}$ is positive definite. Therefore, (52) is a symmetric hyperbolic system. The energy in the system is given by

[TABLE]

and is the same as the energy (10).

Next, we compute the characteristic variables for (52) with respect to the direction normal to the boundary. First, we apply a transform to remove $W$ in front of the time derivative term in (52). Note that since $W$ is positive definite, it can be factorized into $W=FF^{T}$ using the Cholesky factorization. The matrix $F$ is the unique lower triangular matrix given by

[TABLE]

where $\eta=\sqrt{g_{11}g_{22}-g_{12}^{2}}>0$ ($\eta$ is positive since the metric tensor $g_{ij}$ is positive definite). By inserting $W=FF^{T}$, $FF^{-1}=I$ and $F^{-T}F^{T}=I$, equation (52) is transformed to

[TABLE]

We introduce the definitions,

[TABLE]

where

[TABLE]

Then (55) becomes,

[TABLE]

By applying the product rule, we find

[TABLE]

To investigate boundary errors near the bottom boundary $r^{2}=0$, where the normal is $(n_{1},n_{2})=(0,-1)$, we consider $A^{\prime}n_{1}+B^{\prime}n_{2}^{\prime}=-B^{\prime}$. Since $B^{\prime}$ is symmetric, it has the eigendecomposition $B^{\prime}=X_{B}\Lambda_{B}X_{B}^{T}$, where $X_{B}$ is an orthonormal eigenvector basis, and $\Lambda_{B}$ is a diagonal matrix. This eigendecomposition is given by

[TABLE]

Finally, in the direction $(n_{1},n_{2})=(0,-1)$, the characteristic properties of (56) can be studied by the model problem,

[TABLE]

Because $X_{B}$ is constant, we can change variables to ${\mathbf{w}}=X_{B}^{T}{\mathbf{q}}^{\prime}$, which leads to the diagonal system

[TABLE]

The characteristic variables are

[TABLE]

Appendix C Solving the generalized eigenvalue problems for $B(\alpha,\beta)$

The conditions for the matrix $HJ\widetilde{G}$ to be positive definite relies on the solution of two generalized eigenvalue problems (34) and (35), see Lemma 3. In this section we use Kronecker product identities to show how the eigenvalues of the generalized eigenvalue problems can be calculated by instead solving a number of decoupled 1-D eigenvalue problems.

Due to the ordering of the 2-D dependent variables and SBP operators, we can directly use Kronecker product identities to analyze the generalized eigenvalue problem (35). The analysis of (34) follows in a corresponding way after the dependent variables and SBP operators have been permuted.

For all matrices $A$, $B$, $C$ and $D$ where the sizes are such that the products $AC$ and $BD$ are well-defined, it is well-known that $(A\otimes B)(C\otimes D)=(AC)\otimes(BD)$. Furthermore, for all matrices $A$ and $B$, we have $(A\otimes B)^{T}=A^{T}\otimes B^{T}$.

We proceed by analyzing the generalized eigenvalue problem (35). To simplify the notation in this section we drop the superscripts on the matrices ${\mathbb{X}}$ and ${\mathbb{Y}}$. Because the matrix ${\mathbb{X}}$ is positive definite, the generalized eigenvalue problem

[TABLE]

has the same eigenvalues as

[TABLE]

and the eigenvectors are related by $\widetilde{\boldsymbol{y}}={\mathbb{X}}^{1/2}\boldsymbol{y}$.

From the definitions of the 2-D SBP operators in Section 3.1,

[TABLE]

Denote the components of the norm weight matrix by $\hat{M}=\mbox{diag}(\hat{m}_{0},\hat{m}_{1},\ldots,\hat{m}_{N+1})$. Furthermore, let the diagonal matrices $J_{2}g_{2}^{22}$ and $\hat{J}_{2}\hat{g}_{2}^{22}$ have the blocked matrix representations,

[TABLE]

In this case, $K^{(2)}_{k}$ is a diagonal $(N+1)\times(N+1)$ matrix and $\hat{K}^{(2)}_{k}$ is a diagonal $(N+2)\times(N+2)$ matrix. They hold the values of the corresponding metric coefficients along a nodal or a cell-centered grid line in the $r^{2}$-direction, with $r^{1}=\hat{x}_{k}$ (cf. (12) for a definition of $\hat{x}_{k}$).

From the definition of a Kronecker product,

[TABLE]

Therefore,

[TABLE]

In a similar way,

[TABLE]

and the above blocked representation of $J_{2}g_{2}^{22}$ gives

[TABLE]

Note that $\mathbb{X}$ is a positive definite diagonal matrix.

Because of the diagonal blocked structure of the matrices ${\mathbb{X}}$ and ${\mathbb{Y}}$, the eigenvalue problem (57) reduces to $N+2$ decoupled eigenvalue problems

[TABLE]

where,

[TABLE]

Note that $(K^{(2)}_{k})^{-1/2}$ and $M^{-1/2}$ commute because they are both diagonal. To summarize, we have shown that the eigenvalues of the 2-D generalized eigenvalue problem (35) can be calculated by solving $N+2$ decoupled 1-D eigenvalue problems (58).

In a similar fashion, calculating the eigenvalues of the generalized eigenvalue problem (35) reduces to $N+2$ decoupled eigenvalue problems

[TABLE]

where,

[TABLE]

In this case, the diagonal matrices $K^{(1)}_{j}$ and $\hat{K}^{(1)}_{j}$ are the blocks of the diagonal matrices

[TABLE]

Here, $K^{(1)}_{j}$ and $\hat{K}^{(1)}_{j}$ hold the values of the corresponding metric coefficients along a nodal or cell-centered grid line in the $r^{1}$-direction with $r^{2}=\hat{x}_{j}$.