Non-stiff narrow-stencil finite difference approximations of the Laplacian on curvilinear multiblock grids

TL;DR

This paper introduces a new finite difference discretization of the Laplacian on curvilinear multiblock grids that remains stable and non-stiff even with highly skewed grids, enabling larger time steps in simulations.

Contribution

The authors develop stable SATs for the Laplacian that are nearly insensitive to grid skewness, improving computational efficiency for wave propagation modeling.

Findings

The new discretization reduces stiffness on skewed grids.

It enables larger time steps in explicit time integration.

The method applies to variable-coefficient Laplacian problems.

Abstract

The Laplacian appears in several partial differential equations used to model wave propagation. Summation-by-parts--simultaneous approximation term (SBP-SAT) finite difference methods are often used for such equations, as they combine computational efficiency with provable stability on curvilinear multiblock grids. However, the existing SBP-SAT discretization of the Laplacian quickly becomes prohibitively stiff as grid skewness increases. The stiffness stems from the SATs that impose inter-block couplings and Dirichlet boundary conditions. We resolve this issue by deriving stable SATs whose stiffness is almost insensitive to grid skewness. The new discretization thus allows for large time steps in explicit time integrators, even on very skewed grids. It also applies to the variable-coefficient generalization of the Laplacian. We demonstrate the efficacy and versatility of the new method by applying it to acoustic wave propagation problems inspired by marine seismic exploration and infrasound monitoring of volcanoes.

Tables

Table 1. Table 1 : Values of θ H subscript 𝜃 𝐻 \theta_{H} (exact), θ R subscript 𝜃 𝑅 \theta_{R} (approximations, computed in Appendix) and m b subscript 𝑚 𝑏 m_{b} .

	${\theta }_H$	${\theta }_R$	$m_b$
second order	$1/2$	$1.0000$	$2$
fourth order	$17/48$	$0.5776$	$4$
sixth order	$13649/43200$	$0.3697$	$7$

Table 2. Table 2 : CFL numbers used in the classical fourth order Runge–Kutta time-integrator in the application problems, for different orders of spatial accuracy.

	second order	fourth order	sixth order
$CFL$	$0.5$	$0.35$	$0.27$

Table 3. Table 3 : Computed values of θ R subscript 𝜃 𝑅 \theta_{R} for different values of m b subscript 𝑚 𝑏 m_{b} . The chosen ( m b , θ R ) subscript 𝑚 𝑏 subscript 𝜃 𝑅 (m_{b},\theta_{R}) pairs are indicated with boldface font.

	second order	fourth order	sixth order
$m_b=1$	0	0	0
$m_b=2$	1.0000	0	0
$m_b=3$	1.0000	0.1485	0
$m_b=4$	1.0000	0.5776	0
$m_b=5$	1.0000	0.5779	0
$m_b=6$	1.0000	0.5779	0.2318
$m_b=7$	1.0000	0.5779	0.3697
$m_b=8$	1.0000	0.5779	0.3697
$m_b=9$	1.0000	0.5779	0.3697

Equations

$$\begin{array}[]{ll}\ddot{u}=\nabla^{2}u,&\vec{x}\in\Omega^{\prime},\\ u=0,&\vec{x}\in\partial\Omega^{\prime}.\end{array}$$

$$\ddot{\mathbf{u}}=D\mathbf{u}.$$

$$u_{\alpha ,\beta }(r,\theta ,t)=J_{\alpha }(\beta r)\sin \alpha \theta \cos \beta t,$$

$$\Delta t_{\text{new}}/\Delta t_{\text{VM}}=6.69,$$

$$\Delta t_{\text{new}}/\Delta t_{\text{VM}}=6.46,$$

$$\Delta t_{\text{new}}/\Delta t_{\text{VM}}=\mathrm{5.84.}$$

$$u_i{v}_i=\sum _{i=1}^d{u}_i{v}_i,$$

$$\left(u_i,v_i\right)=\sum _{i=1}^d\left(u_i,v_i\right).$$

$$\begin{array}[]{lll}a\ddot{u}=\partial_{x_{i}}b\partial_{x_{i}}u,&\vec{x}\in\Omega^{\prime},&t\in[0,\;T],\\ Lu=0,&\vec{x}\in\Gamma^{\prime},&t\in[0,\;T],\end{array}$$

$${\partial }_{x_i}={\mathcal{K}}_{ij}{\partial }_{{\xi }_j}.$$

$$J{\mathcal{K}}_{ij}{\partial }_{{\xi }_j}={\partial }_{{\xi }_j}J{\mathcal{K}}_{ij}.$$

$${\partial }_{x_i}b{\partial }_{x_i}={\mathcal{K}}_{ij}{\partial }_{{\xi }_j}b{\mathcal{K}}_{ij}{\partial }_{{\xi }_j}=\frac{1}{J}{\partial }_{{\xi }_j}J{\mathcal{K}}_{ij}b{\mathcal{K}}_{ij}{\partial }_{{\xi }_j}.$$

$${\partial }_i={\partial }_{{\xi }_i},{\alpha }_{ijk\ell }={\mathcal{K}}_{ij}J{\mathcal{K}}_{k\ell }.$$

$$a\ddot{u}={\partial }_{x_i}b{\partial }_{x_i}u\iff Ja\ddot{u}={\partial }_j{\alpha }_{ijik}b{\partial }_{k}u.$$

$$\begin{array}[]{lll}Ja\ddot{u}=\partial_{j}\alpha_{ijik}b\partial_{k}u,&\vec{\xi}\in\Omega,&t\in[0,\;T],\\ \widetilde{L}u=0,&\vec{\xi}\in\Gamma,&t\in[0,\;T],\end{array}$$

$${\left(u,v\right)}_{\Omega }=\underset{\Omega }{\int }uv\mathrm{d}\Omega ,{\Vert u\Vert }_{\Omega }^2={\left(u,u\right)}_{\Omega }.$$

$${\left(u,v\right)}_{\Gamma }=\underset{\Gamma }{\int }uv\mathrm{d}\Gamma ,{\Vert u\Vert }_{\Gamma }^2={\left(u,u\right)}_{\Gamma }.$$

$${\left(u,v\right)}_{{\Omega }^{\prime }}={\left(u,Jv\right)}_{\Omega }$$

$${\left(u,v\right)}_{{\Gamma }^{\prime }}={\left(u,\gamma v\right)}_{\Gamma }.$$

$${\left(u,{\partial }_{x_i}v\right)}_{{\Omega }^{\prime }}={\left(u,n_{i}v\right)}_{{\Gamma }^{\prime }}-{\left({\partial }_{x_i}u,v\right)}_{{\Omega }^{\prime }}.$$

$${\left(u,{\partial }_{x_i}v\right)}_{{\Omega }^{\prime }}$$

$$={\left(u,J{\mathcal{K}}_{ij}{\nu }_{j}v\right)}_{\Gamma }-{\left({\partial }_{j}J{\mathcal{K}}_{ij}u,v\right)}_{\Omega }$$

$$={\left(u,J{\mathcal{K}}_{ij}{\nu }_{j}v\right)}_{\Gamma }-{\left(J{\mathcal{K}}_{ij}{\partial }_{j}u,v\right)}_{\Omega }$$

$$={\left(u,\gamma \frac{J}{\gamma }{\mathcal{K}}_{ij}{\nu }_{j}v\right)}_{\Gamma }-{\left({\partial }_{x_i}u,Jv\right)}_{\Omega }$$

$$={\left(u,\frac{J}{\gamma }{\mathcal{K}}_{ij}{\nu }_{j}v\right)}_{{\Gamma }^{\prime }}-{\left({\partial }_{x_i}u,v\right)}_{{\Omega }^{\prime }}.$$

$$n_i=\frac{J}{\gamma }{\mathcal{K}}_{ij}{\nu }_j.$$

$$\frac{\partial u}{\partial \hat{n}}=n_i{\partial }_{x_i}u={\nu }_j\frac{J}{\gamma }{\mathcal{K}}_{ij}{\mathcal{K}}_{ik}{\partial }_{k}u={\nu }_j\frac{{\alpha }_{ijik}}{\gamma }{\partial }_{k}u.$$

$$H_x=h\text{mbox}diag(w_1,\dots ,w_s,1,\dots ,1,w_s,\dots ,w_1),$$

$$H_x{D}_x+D_x^T{H}_x=-e_{\ell }{e}_{\ell }^T+e_r{e}_r^T\text{mbox}(nosumsoverx,\ell ,r),$$

$$e_{\ell }={\left[\begin{array}{c}1,0,\dots ,0\end{array}\right]}^T,e_r={\left[\begin{array}{c}0,\dots ,0,1\end{array}\right]}^T.$$

$$H_x{D}_{xx}(b)$$

$$\text{mbox}(nosumsoverx,\ell ,r),$$

$$\Delta D_x=D_x-{\hat{D}}_x.$$

$$H_x$$

Full text

1 Introduction

Partial differential equations (PDEs) that feature the Laplacian include the second order wave equation. For wave-dominated equations, high-order finite difference methods are computationally efficient [9], assuming that the solution is smooth and the domain not too complex. Finite difference operators with the summation-by-parts (SBP) property [11] lead to energy-stable discretizations when combined with suitable methods for imposing boundary and interface conditions. Weak enforcement via simultaneous approximation terms (SATs) [3] has proven successful in a wide range of applications [6, 23]. A key property of the SBP-SAT method is that it lends itself well to curvilinear multiblock grids, which are necessary in most applications.

For second order differential operators in general, narrow-stencil second-derivative SBP operators [14, 13] provide superior solution accuracy compared to applying a first-derivative operator twice. As a rule of thumb, the global convergence rate is one order higher [14] and the numerical dispersion relation mimics the exact dispersion relation better for marginally resolved modes [10]. Hence, we only consider narrow-stencil operators in this paper. However, the gain in accuracy comes at the cost of more involved stability analysis.

Key properties of the discrete Laplacian that guarantee energy stability for hyperbolic problems are symmetry and negative semidefiniteness with respect to a discrete inner product (see for example [1]). The current state-of-the-art SBP-SAT discretization of the Laplacian was presented by Virta and Mattsson [28], who were the first to derive a symmetric negative semidefinite discretization for Dirichlet boundary conditions on curvilinear multiblock grids. However, as we demonstrate in Section 2, the spectral radius of the discretization by Virta and Mattsson grows rapidly as grid skewness increases. Hence, it forces explicit time integrators to use very small time steps unless the grid is almost Cartesian. We here improve the method by deriving a new discretization whose spectral radius is much less sensitive to grid skewness. The two methods use the same SBP operators and differ only in the SATs for Dirichlet boundary conditions and inter-block couplings. The SATs for Neumann boundary conditions are identical.

Although we use the acoustic wave equation in two dimensions in the numerical experiments, the new discretization is in no way specific to wave equations or two dimensions. Rather, what we present is a symmetric and negative semidefinite discrete Laplacian that may be used to discretize other PDEs, for example the Schrödinger and heat equations. Furthermore, the new discretization is presented for $d$ spatial dimensions and applies to the generalized Laplace operator $\nabla\cdot b(\vec{x})\nabla$.

This paper only considers weak enforcement of boundary and interface conditions using SATs, but it is possible to combine SBP operators with strong enforcement by, for example, injection [4] or ghost points [20, 31]. In fact, for certain classes of boundary and interface conditions, the ghost point technique yields a significantly smaller spectral radius than SATs [31]. We stress that the purpose of this paper is not to find the finite difference discretization of the Laplacian with the smallest possible spectral radius for a given set of boundary conditions. Instead, the purpose is to improve the general-purpose SBP-SAT framework by showing how to impose Dirichlet-type boundary conditions and couple grid blocks with less stiffness than with the existing SBP-SAT method.

The rest of this paper is organized as follows. Section 2 compares the new method with the current state-of-the-art method by Virta and Mattsson [28]. Sections 3 through 8 derive the new method and prove that it is stable. Section 9 demonstrates the efficacy of the new method by solving two application problems inspired by marine seismic exploration and infrasound monitoring of volcanoes.

2 Comparison with Virta and Mattsson, 2014.

This section presents the main result: the new method, derived in subsequent sections of this paper, is significantly more efficient than the method in Virta and Mattsson (VM) [28] when the grid is curved. To compare the two methods, we will consider the second order wave equation with constant coefficients and homogeneous Dirichlet boundary conditions,

[TABLE]

For the spatial discretization we use diagonal-norm SBP operators [13] of interior orders $q=2$, $q=4$ and $q=6$. Near boundaries, the accuracy drops to order $q/2$. As a rule of thumb, one typically observes global convergence rates of order $\min(q,q/2+2)$ for the wave equation [22], but the convergence rate can be both lower and higher in special cases [29, 30]. By the rule of thumb, we expect convergence rates of orders 2, 4, and 5, respectively, in this paper. The boundary conditions are enforced using either the VM SATs or the new SATs.

Section 2.1 shows that the spectral radius of the discrete Laplacian (including SATs) grows much faster with grid skewness with the VM SATs than with the new SATs. This implies that the VM method requires smaller time steps in explicit time steppers, or more iterations for convergence in implicit time steppers combined with iterative solvers.

Section 2.2 compares the two methods in terms of CPU time required to satisfy a given error tolerance. For this comparison we use a mildly curved multiblock grid on the unit disk and an explicit time stepper. Although the VM method is slightly more accurate on a given grid, the new method is significantly more efficient because it allows for much larger time steps. In light of Section 2.1, we anticipate that a more skewed grid would only increase the efficiency gap in favor of the new method.

2.1 Spectral radius as a function of grid skewness

Let $\Omega^{\prime}$ be the parallelogram domain depicted in Figure 1(a). By varying the angle $\varphi$, we obtain a family of parallelograms whose base and height are of unit length and whose bottom and top boundaries are horizontal. The grid lines are chosen to be equidistant and parallel to the domain boundaries. For $\varphi=\pi/2$, $\Omega^{\prime}$ is square and the grid is Cartesian. As $\varphi$ decreases, grid skewness increases. For the experiments in this subsection, we let the grid consist of $31\times 31$ points.

Discretizing (1) in space with the SBP-SAT method yields the semidiscrete equation

[TABLE]

Using the VM SATs and the new SATs derived in later sections of this paper lead to different $D$ matrices, $D_{\text{VM}}$ and $D_{\text{new}}$. Let $\rho(D)$ denote the spectral radius of $D$. Figure 1(b) compares $\rho(D_{\text{VM}})$ and $\rho(D_{\text{new}})$ for the second order accurate discretizations, for $\pi/6\leq\varphi\leq\pi/2$. The spectral radii are normalized by $\rho(D_{\text{new}})$ on the Cartesian grid given by $\varphi=\pi/2$. The new method has a smaller spectral radius for all $\varphi$, but the difference is most apparent for small $\varphi$. As $\varphi$ decreases, the spectral radius of the new method grows much more slowly than that of the VM method.

Figure 2 shows the spectral radii of the fourth and sixth order accurate discretizations, which exhibit similar behavior.

2.2 Computational efficiency on a multiblock grid

We here let $\Omega^{\prime}$ be the unit disk, in which case the model problem (1) has exact solutions

[TABLE]

where $(r,\theta)$ are the polar coordinates, $\alpha\in{\mathbb{Z}}$, $J_{\alpha}$ is the order $\alpha$ Bessel function of the first kind, and $\beta$ is a root of $J_{\alpha}$. We set $\alpha=5$ and $\beta=22.2177998965612$, which is the fifth nonzero root of $J_{5}$. The resulting solution at time $t=0$ is shown in Figure 3(a). We use the multiblock grid configuration in Figure 3(b). We rewrite the semidiscrete system of equations as a first order system in time and use the classical fourth order Runge–Kutta method to integrate until time $t=1$, where we measure the relative $\ell^{2}$ errors. To isolate the accuracy of the spatial discretization, we here use the same time step for the two methods. Let the number of grid points in the center block be $m\times m$ and define the reference grid spacing $h=1/(m-1)$. We here set $\Delta t=0.05h$ for both methods.

Figure 4(a) shows that the two methods converge with the same rates. The VM method consistently yields the smaller error, which is clearly visible in the sixth order case. For second and fourth order, the difference in error is so small that the curves are nearly indistinguishable.

Of more practical relevance than the error as a function of grid spacing is the CPU time required to satisfy a given error tolerance. To compare the CPU times, we let both methods use the largest stable time step (determined by eigenvalue computations). Since the new method has a smaller spectral radius, it ends up using a significantly larger time step. For an example grid, with the configuration in Figure 3(b) and $41\times 41$ grid points in the center block, the time step ratios are

[TABLE]

Figure 4(b) shows the errors as functions of CPU time. For a given order of accuracy, the larger time steps allow the new method to reach the same error level faster, even though the VM method is slightly more accurate on a given grid.

3 Notation conventions

We use the convention to sum over repeated subscript indices so that

[TABLE]

where $d$ denotes the number of spatial dimensions. The summation convention applies to inner products too, i.e.,

[TABLE]

For variable coefficients $a(\vec{x})$, we use the symbol $a$ also in the discrete case, which then is understood to denote a diagonal matrix with the values of $a(\vec{x})$ on the diagonal. The outward unit normal $\hat{n}$ is regarded as a variable coefficient that takes non-zero values only at boundary points. In the discrete setting, the value of $\hat{n}$ at edge and corner points varies with context. When integrating over a face, $\hat{n}$ is understood to denote the unit normal to that face even at edge and corner points.

4 Coordinate transformations

As a model problem that includes the Laplacian, we will consider the second order wave equation,

[TABLE]

where $a>0$ and $b>0$ are variable coefficients, $\Gamma^{\prime}=\partial\Omega^{\prime}$, and the linear operator $L$ represents well-posed boundary conditions. Introduce a one-to-one mapping $x_{i}=x_{i}(\xi_{1},...,\xi_{d})$ from $\Omega=[0,\;1]^{d}$ to $\Omega^{\prime}$. Let $\mathcal{J}$ be the matrix with elements $\mathcal{J}_{ij}=\partial x_{j}/\partial\xi_{i}$, and let $J=\det(\mathcal{J})$. Similarly, let $\mathcal{K}_{ij}=\partial\xi_{j}/\partial x_{i}$. By the chain rule,

[TABLE]

The following metric identities are well known (see [24]):

[TABLE]

Using first (8) and then (9), we have

[TABLE]

For convenience we introduce the notation

[TABLE]

It follows that

[TABLE]

Henceforth, we will consider the transformed wave equation on the unit cube $\Omega$:

[TABLE]

where $\Gamma=\partial\Omega$ and $\widetilde{L}$ is the transformation of $L$.

4.1 Integrals

We will use the $L^{2}$ inner product and norm:

[TABLE]

Similarly, we write

[TABLE]

Note however, that $\left(\cdot,\cdot\right)_{\Gamma}$ is not an inner product but a bilinear form, and $\left\lVert\cdot\right\rVert_{\Gamma}^{2}$ is not a norm but a semi-norm. Further, we define $\gamma$ such that $\mathrm{d}\Gamma^{\prime}=\gamma\,\mathrm{d}\Gamma$ is the surface area element. With this notation, we have

[TABLE]

and

[TABLE]

4.2 Normal derivatives

Let $\hat{n}=[n_{1},\ldots,n_{d}]$ denote the outward unit normal to $\Omega^{\prime}$, and let $\hat{\nu}=[\nu_{1},\ldots,\nu_{d}]$ denote the unit normal to $\Omega$. Using integration by parts, we have

[TABLE]

By instead using the integration-by-parts formula in the reference coordinates, we can also derive

[TABLE]

where we used (9). Since $u$ and $v$ are arbitrary in (18) and (19), it follows that the normal components are related through

[TABLE]

Hence, the normal derivative on $\Omega^{\prime}$ can be expressed as

[TABLE]

5 Summation-by-parts operators

Consider finite-difference approximations of derivatives using an equidistant grid that includes the boundary points of an interval $[x_{\ell},x_{r}]$. In this paper, we consider diagonal-norm SBP operators only. That is, the so-called norm matrix $H_{x}$ has the structure

[TABLE]

where $h$ denotes the grid spacing and the $w_{i}$ are positive dimensionless weights. The first-derivative SBP operators $D_{x}\approx\partial_{x}$ have the property

[TABLE]

where

[TABLE]

The narrow-stencil second derivative operators $D_{xx}(b)\approx\partial_{x}b\partial_{x}$ derived in [13] are based on the same norm matrix and have the property

[TABLE]

where $\hat{D}_{x}$ approximates the first derivative. Note that for the SBP operators considered in this paper, $D_{x}\neq\hat{D}_{x}$. If $D_{x}=\hat{D}_{x}$, then the SBP operators $D_{x}$ and $D_{xx}$ are said to be fully compatible [15]. The difference between the two first-derivative approximations will play an important role in the coming stability analysis and we define

[TABLE]

The SBP operators derived in [13] have $\hat{D}_{x}$ that are accurate of order $q/2+1$, i.e., one order higher than the boundary closure of $D_{x}$. This implies that $\Delta D_{x}{\mathbf{u}}$ is zero to order $q/2$ for all ${\mathbf{u}}$ that are restrictions of smooth functions to the grid.

The matrix $R_{xx}(b)$ is symmetric positive semidefinite and consists of undivided difference approximations in such a way that ${\mathbf{u}}^{T}R_{xx}(b){\mathbf{u}}$ is zero to order $q$ [13].

5.1 Positivity properties

It follows immediately from (22) and (24) that we have

[TABLE]

We define $\theta_{H}=w_{1}$ and write

[TABLE]

or, equivalently,

[TABLE]

Let $x_{\ell,m_{b}}$ and $x_{r,m_{b}}$ denote the $m_{b}$th grid points from the left and right boundaries:

[TABLE]

Let $b_{\ell,min}$ and $b_{r,min}$ denote the minimum of $b$ over the $m_{b}$ leftmost and rightmost grid points:

[TABLE]

The SBP operators considered in this paper also satisfy the following positivity property, which is essential to the new method.

Lemma 1**.**

The second-derivative SBP operators $D_{xx}(b)$ of orders two, four, and six, derived in [13], which can be decomposed as in (25), all satisfy

[TABLE]

for some choice of $m_{b}$ and the dimensionless constant $\theta_{R}>0$.

Proof.

See Appendix. ∎

A property similar but not identical to (32) was used in [5], to prove stability for the elastic wave equation with block interfaces. Because (32) has never before been used in the SBP-SAT literature, we have determined the values of $\theta_{R}$ numerically. We remark that $\theta_{R}$ depends on $m_{b}$. We describe how $m_{b}$ is chosen and how $\theta_{R}$ is computed in Appendix. Table 1 lists the values of all the dimensionless constants for the SBP operators used in this paper.

5.2 Multi-dimensional first-derivative operators

Let operators with subscripts ${\xi_{i}}$ denote one-dimensional operators corresponding to coordinate direction $\xi_{i}$. The multi-dimensional first derivatives $D_{i}\approx\partial_{i}$ and $\hat{D}_{i}\approx\partial_{i}$ are constructed using tensor products:

[TABLE]

[TABLE]

where the $I_{\xi_{i}}$ are one-dimensional identity matrices of appropriate sizes. In analogy with the chain rule (8), we define

[TABLE]

We further define the difference between the two first-derivative approximations:

[TABLE]

The multi-dimensional quadrature is

[TABLE]

Let $\Gamma_{i}^{-}$ and $\Gamma_{i}^{+}$ denote the boundary faces where $\xi_{i}=0$ and $\xi_{i}=1$, respectively. For integration over boundary faces, we define

[TABLE]

Note that $H_{\Gamma_{i}}$ can be used to integrate over $\Gamma_{i}^{+}$ as well as $\Gamma_{i}^{-}$.

For discrete integration over the volume, we write

[TABLE]

Let $e_{S}$ denote a restriction operator that picks out only those solution values that reside on the boundary segment $S$. For discrete integration over the face $\Gamma_{i}^{+}$, for example, we write

[TABLE]

If the surface $S$ is the union of several faces, $S=F_{1}\cup...\cup F_{m}$, we use the convention that

[TABLE]

i.e., the integration is performed over one face at a time. If the integrand contains the unit normal, its values at edges and corners are defined to be the same as on the remainder of that face. As in the continuous case, we define the norm and semi-norm

[TABLE]

With the notation established in this section, we have the discrete integration-by-parts formula

[TABLE]

5.3 Multi-dimensional narrow-stencil second-derivative operators

For any fixed $i$, we construct:

[TABLE]

by using the one-dimensional operator $D_{xx}$ for each grid line. This gives the multi-dimensional SBP property for the second derivative,

[TABLE]

The structure of $R_{ii}(b)$ in (45) is important and requires a closer look. Let $\widetilde{R}_{ii}(b)$ be the operator that performs the action of $R_{xx}(b)$ for each grid line in direction $i$. Further, define

[TABLE]

Then,

[TABLE]

The $R_{ii}$ matrices are symmetric positive semidefinite.

The property (45) can be equivalently written as

[TABLE]

for all ${\mathbf{u}}$ and ${\mathbf{v}}$.

5.4 Combining narrow-stencil derivatives and mixed derivatives

Consider discretizing a term $\partial_{i}b\partial_{j}$. Because we want to use narrow-stencil second derivatives for the non-mixed second derivatives, the following notation will be useful:

[TABLE]

Combining the two integration-by-parts formulas (43) and (48) leads to the integration-by-parts formula

[TABLE]

5.5 Multi-dimensional positivity properties

To suppress unnecessary notation, we assume that the grid spacing in the reference domain is $h$ in each dimension (we stress that the analysis does not rely on this assumption). It follows from (29) that, for any fixed $j$,

[TABLE]

Using (51) we can derive

[TABLE]

which we summarize as

[TABLE]

The multi-dimensional positivity property that follows from (32) is

[TABLE]

where $b_{min}$ denotes the minimum of $b$ over the $m_{b}$ grid points closest to the boundary $\Gamma$ along each grid line. That is, to evaluate $b_{min}$ at a given grid point on $\Gamma$, one needs to evaluate $b$ at the $m_{b}$ first grid points on the same grid line orthogonal to $\Gamma$.

6 Continuous energy balance

We return to the transformed equation (13),

[TABLE]

The energy method, i.e. multiplying (55) by $\dot{u}$ and integrating over $\Omega$, yields

[TABLE]

The volume integrals satisfy

[TABLE]

and

[TABLE]

We define the energy $\mathcal{E}$ as

[TABLE]

which allows us to write (56) as

[TABLE]

7 Semidiscrete energy balance

We discretize the transformed problem (13) in space as

[TABLE]

where $SAT$ imposes the boundary conditions and will be specified later. The discrete energy method amounts to multiplying (61) by $\dot{{\mathbf{u}}}^{T}H$ from the left, which leads to

[TABLE]

We have

[TABLE]

The integration-by-parts formula (50) yields

[TABLE]

We define the discrete normal derivative $D_{\hat{n}}$ as

[TABLE]

and the discrete energy $E$ as

[TABLE]

so that the discrete energy rate reads

[TABLE]

By comparing (65) with (21) and recalling that $\Delta D_{k}{\mathbf{u}}$ is zero to order $q/2$, we conclude that $D_{\hat{n}}$ approximates $\partial/\partial{\hat{n}}$. Further, since the $R_{kk}$ are positive semidefinite and ${\mathbf{u}}^{T}R_{kk}{\mathbf{u}}$ is zero to order $q$, the discrete energy $E$ is nonnegative and approximates the continuous energy $\mathcal{E}$ defined in (59).

7.1 Enforcing Neumann boundary conditions

Homogeneous Neumann boundary conditions, i.e.,

[TABLE]

are energy-conserving for the continuous equation, c.f. (60). By inspecting (67), we see that the discrete energy $E$ will be conserved, and hence the scheme will be energy-stable, if we can find a $SAT$ that is consistent with (68) and satisfies

[TABLE]

Let $F$ denote the set of all faces, i.e. $F=\{\Gamma_{1}^{-},\ldots,\Gamma_{d}^{-},\Gamma_{1}^{+},\ldots,\Gamma_{d}^{+}\}$. We achieve (69) by setting

[TABLE]

7.2 Enforcing Dirichlet boundary conditions

Homogeneous Dirichlet conditions, i.e.,

[TABLE]

are energy-conserving for the continuous problem. As opposed to the case with Neumann conditions, we cannot choose consistent SATs that make the boundary terms in the energy rate vanish. Instead, we will show that a slightly different energy $E_{D}\approx E$ is conserved.

General SATs of the form

[TABLE]

lead to

[TABLE]

We set $B=C-A$, where $A$ is symmetric. Inserting this ansatz in the discrete energy rate (67) yields

[TABLE]

For symmetry we will need

[TABLE]

which implies $C=bD_{\hat{n}}$. We then obtain

[TABLE]

We now obtain the energy balance

[TABLE]

where

[TABLE]

To prove stability, it remains to prove that we can choose $A$ so that $E_{D}$ is a nonnegative quantity. Loosely speaking, we will need $A$ to be sufficiently positive. Note that $E_{D}$ (just like $E$) is a high-order approximation of the continuous energy $\mathcal{E}$ because the boundary integrals are zero to the order of accuracy due to the boundary condition.

To prove that $E_{D}$ is nonnegative, we need to use the positivity of $E$. Since the indefinite term in (78) is a surface integral, we shall bound $E$ from below by relevant surface integrals. The terms in $E$ that stem from the Laplacian are

[TABLE]

Applying the positivity property (53) to Term 1 yields

[TABLE]

Now consider Term 2. To simplify the notation, we define

[TABLE]

Applying the positivity property (54) yields

[TABLE]

Using (80) and (82), we obtain

[TABLE]

Expanding the indefinite term $\left(bD_{\hat{n}}{\mathbf{u}},{\mathbf{u}}\right)_{\Gamma^{\prime},h}$ yields

[TABLE]

We obtain

[TABLE]

where

[TABLE]

[TABLE]

and

[TABLE]

We need to choose $A$ so that $\left(A{\mathbf{u}},{\mathbf{u}}\right)_{\Gamma^{\prime},h}$ compensates for the zero entries in $\widetilde{Y}$ and $\widetilde{Z}_{k}$. We make the ansatz

[TABLE]

where $\tau_{H,R}$ are penalty parameters that we have yet to choose. By noting that

[TABLE]

we can write

[TABLE]

We now obtain

[TABLE]

where

[TABLE]

The matrices $Y$ and $Z_{k}$ are positive semidefinite if

[TABLE]

respectively. We have now proven the following theorem.

Theorem 2**.**

The scheme

[TABLE]

with

[TABLE]

is stable if

[TABLE]

In all simulations in this paper, we set $\tau_{H}=d/(h\theta_{H})$ and $\tau_{R}=1/(h\theta_{R})$, i.e., right on the limit of provable stability. The drawback to using larger values of $\tau_{H,R}$ is that it increases the spectral radius.

We remark that the fundamental difference compared to the approach of Virta and Mattsson [28] is that they did not utilize Lemma 1, and hence not the positivity property (82), which is a consequence of Lemma 1. Hence, they were forced to use a larger $SAT$ to prove that $E_{D}$ is nonnegative, which in turn caused the increase in spectral radius.

If the SBP operators $D_{x}$ and $D_{xx}$ are fully compatible, then $\Delta D_{x}=0$ and hence $\Delta D_{i}=0$, which implies that the scheme in Theorem 2 is stable with $\tau_{R}=0$.

8 Interfaces

Let $\Gamma^{\prime}_{I}$ denote the interface between two domains $\Omega^{\prime}_{u}$ and $\Omega^{\prime}_{v}$. We use superscripts $u$ and $v$ to denote variable coefficients that correspond to the two different sides of the interface. We consider the problem

[TABLE]

augmented with suitable boundary conditions. The system (98) transforms into

[TABLE]

Define the energies

[TABLE]

[TABLE]

Assuming energy-conserving boundary conditions, the energy method yields

[TABLE]

Theorem 3**.**

The scheme

[TABLE]

[TABLE]

where

[TABLE]

is a stable discretization of (99) if

[TABLE]

Proof.

Define the discrete energies

[TABLE]

and

[TABLE]

where

[TABLE]

The discrete energy method yields

[TABLE]

where

[TABLE]

and

[TABLE]

Moving all terms in (110) to the left-hand side yields

[TABLE]

where

[TABLE]

Note that the discrete energy $E_{I}$ approximates the continuous energy $\mathcal{E}_{u}$ + $\mathcal{E}_{v}$ because the interface terms in (114) would be zero if the interface conditions were fulfilled exactly. It remains to show that $E_{I}\geq 0$. We start by expanding the indefinite terms. By (84) we have

[TABLE]

where we have defined $\beta=b\frac{J}{\gamma}$. Similarly, we have

[TABLE]

[TABLE]

[TABLE]

and

[TABLE]

The positivity properties (80) and (82) yield

[TABLE]

and

[TABLE]

Using (90), we can show that, for example,

[TABLE]

and hence

[TABLE]

Define

[TABLE]

After some algebra we obtain

[TABLE]

where

[TABLE]

and

[TABLE]

The matrices $Y$ and $Z_{k}$ are positive semidefinite if $\tau_{H}\geq\frac{d}{4\theta_{H}}$ and $\tau_{R}\geq\frac{1}{4\theta_{R}}$, respectively. ∎

In all experiments with interfaces in this paper, we set $\tau_{H}=d/(4\theta_{H})$ and $\tau_{R}=1/(4\theta_{R})$, i.e., right on the limit of provable stability.

9 Applications

In this section we will use the new method to solve two application problems. To choose the time-step without having to compute eigenvalues of the spatial discretization we use the procedure described below.

Consider first a single-block grid. Let $c=\sqrt{b/a}$ denote the wave speed. Let $h_{i}$ denote the grid spacing in the $\xi_{i}$-direction in the reference domain. We compute approximate physical grid spacings $\widetilde{h}_{j}$ as

[TABLE]

and define

[TABLE]

Note that both $\widetilde{h}$ and $c$ are non-constant on the grid, in general. Let $\mathcal{G}$ denote the set of all grid points. We compute the time step $\Delta t$ as

[TABLE]

where $CFL$ is a dimensionless constant that depends on the spatial order of accuracy and the choice of time-integrator. For a multiblock grid, it only remains to loop over all grid blocks and choose the smallest $\Delta t$ suggested by any block.

For the example applications in this section we use the classical fourth order Runge–Kutta method with the $CFL$ values listed in Table 2. The numbers have been tuned to yield time-steps that are slightly smaller than the largest stable time-step, for the grid configurations that we have tried. While we cannot guarantee that these $CFL$ values yield stability for every possible grid, they should at least provide a reasonable starting point.

9.1 Marine seismic exploration with ocean-bottom nodes

The Marmousi2 velocity model [12] is a fully elastic velocity model. It is an extension of the acoustic Marmousi model [25, 26], which has become a standard benchmark problem for seismic imaging and full waveform inversion. Figure 5 shows the P-wave speed $c_{p}$ in Marmousi2. The red line at $y=y_{I}=-0.45$ km marks the seafloor. Above the red line there is water, and below the red line there is solid Earth.

In this paper, we shall approximate the response of both ocean and solid as acoustic. The acoustic equations can be written on first order form as

[TABLE]

where $p$ is the pressure perturbation, $v_{i}$ is the particle velocity vector, $\rho$ is density, $K$ is bulk modulus, and $F_{i}$ and $G$ are source terms. The Marmousi2 dataset provides density $\rho$ and P-wave speed $c_{p}$. We set the bulk modulus to be $K=\rho c_{p}^{2}$, so that the wave speed in (126) matches the P-wave speed in the elastic model. Both $\rho$ and $K$ are spatially variable within the solid. The system (126) may be written as a second-order equation for $p$ only:

[TABLE]

In seismic exploration, one attempts to solve optimization problems of the form [7, 27]

[TABLE]

where

[TABLE]

and $v_{i}$ and $p$ satisfy (126) and appropriate boundary, interface, and initial conditions. The data $v_{i,data}$ and $p_{data}$ are measured by geophones and hydrophones, respectively, at the locations $\vec{x}_{v}$ and $\vec{x}_{p}$. The weight functions $w_{v,i}$ and $w_{p}$ may be used to window the signal in time. In practice one often regularizes the optimization problem by minimizing $M+S$ rather than $M$, where $S$ is a regularization term. In this paper, we set $S=0$ for simplicity.

We introduce the residuals

[TABLE]

As part of computing the gradient of $M$ with respect to the material parameters, one may solve the adjoint equations [7, 27] for $v_{i}^{\dagger}$ and $p^{\dagger}$,

[TABLE]

where $\tau=T-t$ is reverse time. Notice that the residuals $\widetilde{v}_{i}$ and $\widetilde{p}$ appear as singular source terms at the receiver locations. The second-order equation for $p^{\dagger}$ only is

[TABLE]

Notice that (132) is of the form (7), with added source terms. We use the technique developed in this paper to discretize the operator $\partial_{x_{i}}\rho^{-1}\partial_{x_{i}}$. Pressure and the vertical velocity component are continuous across the fluid-solid interface. The horizontal velocity component, however, may be discontinuous. To avoid differentiating across the discontinuity, we use two grid blocks: one for the water and one for the solid Earth. Expressed in terms of $p^{\dagger}$, the interface conditions are

[TABLE]

where subscripts $w$ and $s$ denote water and solid sides, respectively. The interface conditions are of the same form as those considered in (98) and we may use the interface discretization presented in Theorem 3.

References [16, 21] describe how to discretize $\delta$ and $\partial_{x_{i}}\delta$. When using $q$th order finite difference stencils in the interior, the discretizations of $\delta$ and $\partial_{x_{i}}\delta$ satisfy $q$ and $q+1$ moment conditions, respectively. The distribution $\partial_{x_{i}}\delta$ satisfies

[TABLE]

for all smooth fields $u_{i}$. Let ${\mathbf{u}}_{i}$ denote the restriction of $u_{i}$ to the grid. A discrete approximation $\mathbf{d}_{x_{i},\vec{x}_{v}}^{\prime}\approx\partial_{x_{i}}\delta(\vec{x}-\vec{x}_{v})$ is required to satisfy the moment conditions [21]

[TABLE]

for all $u_{i}$ that are polynomials of degree less than or equal to $q$. To solve (132), we actually need to consider moment conditions for $\partial_{x_{i}}\rho^{-1}\delta(\vec{x}-\vec{x}_{v})$. Integration by parts yields

[TABLE]

The moment conditions on $\mathbf{d}_{x_{i},\vec{x}_{v},\rho^{-1}}^{\prime}\approx\partial_{x_{i}}\rho^{-1}\delta(\vec{x}-\vec{x}_{v})$ are thus

[TABLE]

for all $u_{i}$ that are polynomials of degree less than or equal to $q$. Given $\mathbf{d}_{x_{i},\vec{x}_{v}}^{\prime}$ that satisfies the appropriate moment conditions, we note that we may set

[TABLE]

and satisfy the required moment conditions on $\mathbf{d}_{x_{i},\vec{x}_{v},\rho^{-1}}^{\prime}$. Because we solve the equations on second order form, no smoothness conditions are required.

The far-field boundaries are treated with super-grid-scale absorbing layers [2, 17]. At the water surface, we impose the Dirichlet condition $p=0$. We use the classical fourth order Runge–Kutta method for time integration.

We attempt to set up simulations that are representative of the adjoint equations in marine seismic exploration. Let $W$ denote the Ricker wavelet [18, 19] with peak frequency $f$ centered at time $t_{0}$, i.e.,

[TABLE]

We choose $f=5$ Hz and $t_{0}=0.2$ s, and set

[TABLE]

To illustrate the difference between the two different adjoint sources, we run two separate simulations with only one active source in each. We place the sources at the same location, but on opposite sides of the water-solid interface. That is, we set $\vec{x}_{v}=\begin{bmatrix}10.5\mbox{ km},&y_{I}-\end{bmatrix}$ and $\vec{x}_{p}=\begin{bmatrix}10.5\mbox{ km},&y_{I}+\end{bmatrix}$.

Let $c_{0}=1500$ m/s and $\rho_{0}=1010$ kg/m3 denote reference values for sound speed and density of water. We define non-dimensional pressures and velocities as

[TABLE]

in the case of only hydrophone data, and

[TABLE]

in the case of only geophone data. Figure 6 shows snapshots of $p^{\dagger}_{nd}$, computed with the fourth order ($q=4$) spatial discretization. The left column corresponds to having only hydrophone data ($\widetilde{v}_{i}=0,\,i=1,2$), and the right column corresponds to having only geophone data ($\widetilde{p}=0$). The bottom row shows space-time plots of pressure at the seafloor, on the ocean side ($y=y_{I}+$). Figure 7 shows the horizontal component of $v^{\dagger}_{i,nd}$ in the same way.

Figures 6 and 7 show that, despite placing the singular sources right on the fluid-solid interface, there are no visible artifacts from the source discretizations. Further, Figure 6 shows that the simulated pressure is continuous across the interface, as prescribed by (133). Because the interface conditions are imposed weakly, the discrete pressure values on the two sides of the interface are actually only equal to the order of accuracy, but this small difference is not visible to the eye. Figure 7, on the other hand, shows that a discontinuity in the horizontal velocity component develops at the interface. Despite the discontinuity, there are no signs of numerical artifacts arising from the interface treatment.

9.2 Resonant infrasound tones of Volcán Villarrica

In [8] it was demonstrated how infrasound can be used to track the surface of the lava lake in the crater of Volcán Villarrica (Chile). The position of the lava lake may in turn be used to forecast volcanic eruptions. The nearest infrasound station is located about 4 km northwest of the crater. To obtain a relevant crater topography we use a two-dimensional average of the three-dimensional crater geometry in [8], see $x\in[-173\mbox{ m},\;173\mbox{ m}]$ in Figure 8(b). The conduit, centered at $x=0$ m, is approximated to have perfectly vertical sides. Figure 8(a) shows the full computational domain. Outside of the crater region, we let the topography be described by a degree five polynomial that approximately matches the elevation drop of the Villarrica mountain in the northwest direction.

We assume constant density $\rho=1\,\mbox{kg/m}^{3}$ in the atmosphere and introuduce the velocity potential $\phi$ such that

[TABLE]

We further approximate the sound speed $c=\sqrt{K/\rho}$ as constant and set $c=331$ m/s. The acoustic equations (126) with $F_{i}=0$, $G=0$ are then satisfied if

[TABLE]

The Laplacian in (144) is discretized using the scheme developed in this paper, and we choose fourth order accuracy ($q=4$) in the numerical experiments. We use the classical fourth order Runge–Kutta method for time-integration. The far-field boundaries are again treated with super-grid-scale absorbing layers [2, 17]. The conduit walls and the mountain sides are modeled as rigid walls, resulting in the boundary condition $\vec{v}\cdot\hat{n}=0$. At the conduit bottom we simulate movement of the lava lake by imposing the boundary condition

[TABLE]

with $t_{0}=1$ s and $\sigma=0.13$ s. The same boundary data were used in [8], with the motivation that they include frequencies that are characteristic of lava lake events.

We define a non-dimensional pressure as

[TABLE]

Figure 9 shows snapshots of $p_{nd}$.

Figure 10 shows space-time plots of $p_{nd}$ at position $x=0$, i.e., within and above the conduit. One can observe how the initial pressure pulse, propagating upwards, is partly reflected at the top of the conduit.

Figure 11(a) shows the recorded pressure at the location corresponding to the nearst infrasound station, for different conduit lengths. Figure 11(b) shows the corresponding frequency spectra. The longer the conduit, the lower the crater resonance frequency and the higher the quality factor.

This simulation illustrates the benefits of the curvilinear multiblock finite difference discretization. The method takes topography into account without using staircasing, which might introduce artifacts. Further, it imposes the solid wall boundary condition accurately and stably.

10 Conclusions

We have improved the SBP-SAT method for discretizing the Laplacian on curvilinear multiblock grids in $d$ dimensions. Compared to the previous state-of-the-art method by Virta and Mattsson [28], the new method is significantly less stiff, in particular when the grid is highly skewed. On a simple curvilinear multiblock grid for the unit disk, the largest stable time-step when solving the wave equation with explicit time-stepping is approximately 6 times larger with the new method. Although the method by Virta and Mattsson is typically slightly more accurate on a given grid, our numerical experiments show that the new method is significantly more efficient in terms of CPU time required to achieve a given error tolerance. Further, our results indicate that for more complex grids, which will inevitably contain highly skewed grid cells, the time-step ratio will be even larger in favor of the new method.

The difference compared to the method by Virta and Mattsson lies solely in the SATs for Dirichlet boundary conditions and interface couplings. In essence, the new method is less stiff because the SATs are smaller in magnitude. The reason we were able use smaller SATs and still prove stability is the novel positivity property in Lemma 1, which was not used by Virta and Mattsson.

To illustrate the usefulness of the new method, we have applied it to two problems inspired by marine seismic exploration and infrasound monitoring of volcanoes. The seismic exploration problem uses a Cartesian two-block grid with a discontinuity in material parameters at the ocean floor, where the ocean couples to the seabed. The numerical method handles both the fluid-solid coupling and singular source terms (both $\delta$ functions and derivatives of $\delta$ functions), placed right at the fluid-solid interface, without signs of numerical artifacts. The volcano monitoring problem involves non-trivial topography that requires a curvilinear multiblock grid. Again, the discrete solution shows no sign of numerical artifacts.

MATLAB code that reproduces all figures in this paper is available at

https://bitbucket.org/martinalmquist/laplacian_curvilinear.

Acknowledgements

We thank Leighton Watson for providing the Villarrica crater topography. M. Almquist gratefully acknowledges support from the Knut and Alice Wallenberg Foundation (Dnr. KAW 2016.0498).

Appendix: proof of Lemma 1

By construction, $R_{xx}(a)$ is linear in its argument so that, for all functions $a,b$ and scalars $C$,

[TABLE]

Further, $R_{xx}$ is symmetric positive semidefinite in the sense that, for all $a(x)\geq 0$,

[TABLE]

For detailed information about the structure of $R_{xx}$ we refer to [13]. It follows from linearity and positive semidefiniteness that if $a(x)\geq b(x)$, then

[TABLE]

Given grid endpoints $x_{\ell}$ and $x_{r}$, we define cutoff points at the $m_{b}$th grid point from the boundaries:

[TABLE]

Let $b_{\ell,min}$ and $b_{r,min}$ denote the minimum of $b(x)$ over the $m_{b}$ leftmost and rightmost grid points,

[TABLE]

Let $\chi_{\ell,r}$ denote the indicator functions

[TABLE]

Assuming that there are at least $2m_{b}$ grid points, then $x_{\ell,m_{b}}<x_{r,m_{b}}$ and hence, for all nonnegative functions $b$,

[TABLE]

We thus have

[TABLE]

To determine the constant $\theta_{R}$ in the property

[TABLE]

it is thus sufficient to consider ${\mathbf{u}}^{T}R_{xx}(\chi_{\ell}){\mathbf{u}}$ or, equivalently, ${\mathbf{u}}^{T}R_{xx}(\chi_{r}){\mathbf{u}}$.

We have determined $\theta_{R}$ by computing eigenvalues of

[TABLE]

for different values of the scalar $\zeta$. Since $R_{xx}(\chi_{\ell})$ is only semidefinite, it has one or more zero eigenvalues. We define $\theta_{R}$ as the smallest value such that $N(\theta_{R})$ has negative eigenvalues. In practice, we compute the smallest value such that the most negative eigenvalue of $N(\theta_{R})$ is not within roundoff error from zero. Let $\lambda_{min}(\zeta)$ denote the smallest eigenvalue of $N(\zeta)$. Figure 12 shows $\left|\lambda_{min}(\zeta)\right|+10^{-16}$ on a logarithmic scale, for the fourth order accurate case. The grid spacing is $h=1$ so that roundoff errors in the eigenvalue computations are expected to be of order $10^{-16}$. It is evident that $\lambda_{min}(\zeta)$ is nonzero for $\zeta>0.5776=\theta_{R}$.

We remark that $\theta_{R}$ is independent of $m$, as long as the total number of grid points $m$ is at least $2m_{b}$. However, the value of $\theta_{R}$ increases with $m_{b}$ and approaches an asymptotic value, as shown in Table 3. The drawback to having large $m_{b}$ is that it requires computing the minimum of the coefficient $b$ over many grid points, which can potentially make $b_{\ell,min}$ and $b_{r,min}$ smaller. To minimize the penalty strength required for the stability proof, we would like to maximize the products $b_{\ell,min}\theta_{R}$ and $b_{r,min}\theta_{R}$. We thus choose the smallest value of $m_{b}$ that yields $\theta_{R}$ near the asymptotic value. Based on Table 3 we choose $m_{b}=2,4,7$, for orders two, four and six. The same $m_{b}$ values were chosen in [28] when computing similar constants for the matrix $M_{xx}(b)=D_{x}^{T}HbD_{x}+R_{xx}(b)$.