A ray-based IPDG method for high-frequency time-domain acoustic wave propagation in inhomogeneous media

TL;DR

This paper introduces a ray-based IPDG method for high-frequency acoustic wave simulation in inhomogeneous media, achieving higher accuracy and stability with fewer unknowns compared to traditional methods.

Contribution

The proposed method combines geometrical optics-inspired basis functions with wavefront tracking to improve high-frequency wave simulation accuracy in inhomogeneous media.

Findings

Significantly higher accuracy than polynomial-based IPDG methods

Moderate growth of relative errors with increasing frequency

Effective dimension reduction and wavefront tracking

Abstract

The numerical approximation of high-frequency wave propagation in inhomogeneous media is a challenging problem. In particular, computing high-frequency solutions by direct simulations requires several points per wavelength for stability and usually requires many points per wavelength for a satisfactory accuracy. In this paper, we propose a new method for the acoustic wave equation in inhomogeneous media in the time domain to achieve superior accuracy and stability without using a large number of unknowns. The method is based on a discontinuous Galerkin discretization together with carefully chosen basis functions. To obtain the basis functions, we use the idea from geometrical optics and construct the basis functions by using the leading order term in the asymptotic expansion. Also, we use a wavefront tracking method and a dimension reduction procedure to obtain dominant rays in each cell. We show numerically that the accuracy of the numerical solutions computed by our method is significantly higher than that computed by the IPDG method using polynomials.Moreover, the relative errors of our method grow only moderately as the frequency increases.

Tables

Table 1. Table 1: In Example 1, the numerical results corresponding to our method.

$\omega$	$1/h$	$\omega h$	Condition number	Number of dof	Relative $L^2$ error (%)
$10\pi$	10	$\pi$	7.59e+07	456	16.856
$20\pi$	20		8.48e+07	1844	14.339
$40\pi$	40		1.08e+08	7376	10.976
$80\pi$	80		1.16e+08	29532	10.128
$160\pi$	160		1.33e+08	118184	11.183

Table 2. Table 2: In Example 1, the numerical results corresponding to the standard IPDG method with bilinear basis.

$\omega$	$1/h$	$\omega h$	Number of dof	Relative $L^2$ error (%)
$10\pi$	100	$\pi /10$	40000	20.11
$20\pi$	200		160000	32.72
$40\pi$	400		640000	57.83
$80\pi$	800		2560000	105.00
$160\pi$	1600		10240000	176.80

Table 3. Table 3: In Example 2, the condition number and the number of degrees of freedom of the original system and the POD systems, and the relative L 2 superscript 𝐿 2 L^{2} error under different settings of ω 𝜔 \omega and meshsize h ℎ h .

$\omega$	$1/h$	$\omega h$	Condition number		Number of dof		Relative
			original	POD	original	POD	$L^2$ error (%)
$10\pi$	10	$\pi$	7.59e+07	4.95e+04	456	442	17.07
	20	$\pi /2$	5.80e+09	6.87e+05	1844	1783	11.56
	40	$\pi /4$	3.52e+12	2.68e+06	7376	7132	5.56
	80	$\pi /8$	6.89e+14	5.47e+06	29532	28540	2.39
	160	$\pi /16$	5.65e+17	5.64e+06	118184	113378	0.85
$20\pi$	20	$\pi$	8.48e+07	5.86e+06	1844	1784	14.45
	40	$\pi /2$	8.20e+09	6.71e+05	7376	7132	7.71
	80	$\pi /4$	1.12e+12	2.70e+06	29532	28549	3.35
	160	$\pi /8$	1.06e+14	5.59e+06	118184	114199	1.18
$40\pi$	40	$\pi$	1.08e+08	5.81e+06	7376	7136	10.99
	80	$\pi /2$	7.17e+09	6.76e+05	29532	28549	5.04
	160	$\pi /4$	6.59e+11	2.92e+06	118184	114238	2.15

Table 4. Table 4: In Example 3, the condition number and the number of degrees of freedom of the original system and the POD systems, and the relative L 2 superscript 𝐿 2 L^{2} error under different settings of ω 𝜔 \omega and meshsize h ℎ h .

$\omega$	$1/h$	$\omega h$	Condition number		Number of dof		Relative
			original	POD	original	POD	$L^2$ error (%)
$10\pi$	10	$\pi$	1.12e+09	2.22e+06	912	884	15.09
$20\pi$	20		3.91e+08	4.23e+06	3688	3566	11.54
$40\pi$	40		4.39e+09	4.43e+06	14752	14263	8.93
$80\pi$	80		3.46e+10	4.56e+06	59064	57090	8.75
$160\pi$	160		8.07e+11	4.54e+06	236368	228443	9.75

Table 5. Table 5: In Example 4, the condition number and the number of degrees of freedom of the original system and the POD systems, and the relative L 2 superscript 𝐿 2 L^{2} error obtained by comparing the solution of the POD systems to the reference solution, under different settings of ω 𝜔 \omega and meshsize h ℎ h .

$\omega$	$1/h$	$\omega h$	Condition number		Number of dof		Relative
			original	POD	original	POD	$L^2$ error (%)
$10\pi$	10	$\pi$	1.37e+10	8.04e+07	1256	940	8.58
$20\pi$	20		6.76e+11	4.99e+07	4880	3654	8.82
$40\pi$	40		1.59e+11	1.46e+08	19792	15200	7.67
$80\pi$	80		5.80e+12	1.85e+08	79244	59535	9.37
$160\pi$	160		8.77e+12	1.82e+08	319592	240167	10.33

Equations

$$u_{tt}(x,t)-c^2(x)\Delta u(x,t)=0,(x,t)\in {\mathbb{R}}^2\times [0,T],$$

$$\displaystyle u(x,t)\approx\text{superposition of }\Big{\{}A_{k}(x,t)e^{i\omega\phi_{k}(x,t)}\Big{\}}_{k=1}^{N},$$

$$({\varphi }_k{)}_t(x,t)\pm \left|\nabla {\varphi }_k(x,t)\right|c(x)=0.$$

$$\Delta v(x)+\frac{{\omega }^2}{c^2(x)}v(x)=0,$$

$$\displaystyle v(x)\approx\text{superposition of }\Big{\{}A_{k}(x)e^{i\omega\phi_{k}(x)}\Big{\}}_{k=1}^{N},$$

$$\begin{array}{rl}{u}_k(x,t)=& \left(A_k(x_0,t)+\nabla A_k(x_0,t)\cdot (x-x_0)\right)e^{i\omega ({\varphi }_k(x_0,t)+\nabla {\varphi }_k(x_0,t)\cdot (x-x_0))}\\ & +\mathcal{O}(h^2+\omega h^2+{\omega }^{-1}),\end{array}$$

$${\Theta }_K:=\left\{\nabla {\varphi }_k(x_K,t):1\le k\le N_K,0\le t\le T\right\}$$

$$V({\Theta }_K)=\text{mbox}span\left\{{\phi }_{j,K}{e}^{i\omega p\cdot (x-x_K)}:p\in {\Theta }_K,j=1,2,3,4\right\}$$

$$V_c={\cup }_{K\in {\mathcal{T}}_h}V({\Theta }_K).$$

$$\left(\frac{{\partial }^2{u}_h}{\partial t^2},v_h\right)+a_h^{\gamma }(u_h,v_h)=0,$$

$$\begin{array}{rl}{a}_h^{\gamma }(u,v):=& {\int }_{\Omega }\nabla u\cdot \overline{\nabla v}dx-\sum _{F\in {\mathcal{F}}_h}{\int }_F\{\nabla u\cdot n\}\overline{\left[v\right]}ds-\sum _{F\in {\mathcal{F}}_h}{\int }_F\left[u\right]\overline{\{\nabla v\cdot n\}}ds\\ & +\frac{\gamma }{h}{\int }_F\left[u\right]\overline{\left[v\right]}ds,\end{array}$$

$$\{u\}:=\frac{1}{2}(u^{+}+u^{-})\text{mbox}and\left[u\right]:=u^{+}{n}^{+}+u^{-}{n}^{-},$$

$$\{\mathbf{\sigma }\}:=\frac{1}{2}({\mathbf{\sigma }}^{+}+{\mathbf{\sigma }}^{-})\text{mbox}and\left[\mathbf{\sigma }\right]:={\mathbf{\sigma }}^{+}\cdot n^{+}+{\mathbf{\sigma }}^{-}\cdot n^{-}.$$

$${\hat{\Theta }}_K:=\left\{p:p=\nabla {\varphi }_k(x_K,t_m),\text{for all }k,t_m\text{ satisfying }1\le k\le N_K,0\le t_m\le T\right\}.$$

$${\tilde{\Theta }}_K:=\left\{p:p=\nabla {\varphi }_k(x_K,t_m),\text{for some }k,t_m\text{ satisfying }1\le k\le N_K,0\le t_m\le T\right\}.$$

$$V({\tilde{\Theta }}_K)=\text{mbox}span\left\{{\phi }_{j,K}{e}^{i\omega p\cdot (x-x_K)}:p\in {\tilde{\Theta }}_K,j=1,2,3,4\right\}$$

$$V={\cup }_{K\in {\mathcal{T}}_h}V({\tilde{\Theta }}_K).$$

$$\left(\frac{u_h^{n+1}-2u_h^n+u_h^{n-1}}{\Delta t^2},v_h\right)+a_h^{\gamma }(u_h^n,v_h)=0,$$

$$\varphi (y(t_0,\alpha ),t_0)\equiv {\varphi }_0\text{mbox}andq(t_0,\alpha )=\nabla \varphi (y(t_0,\alpha ),t_0)$$

$${\varphi }_t(x,t)+c(x)\left|\nabla \varphi (x,t)\right|=0.$$

$$x_t=c(x)\frac{p}{\left|p\right|},p_t=-\left|p\right|\nabla c(x).$$

$$r_j^n:=(y(t_n,{\alpha }_j),q(t_n,{\alpha }_j)).$$

$$I_j:=\max \left\{i:i^{\ast }\le j\right\},$$

$$\mathcal{I}(r_1,r_2,r_3;x_K)={\lambda }_1{p}_1+{\lambda }_2{p}_2+{\lambda }_3{p}_3,$$

$$x_K={\lambda }_1{x}_1+{\lambda }_2{x}_2+{\lambda }_3{x}_3\text{mbox}and{\lambda }_1+{\lambda }_2+{\lambda }_3=1.$$

$$\mathbf{M}^{c}:=\left(\int_{\Omega}c^{-2}\,\psi_{j}(x)\overline{\phi}_{i}(x)\right)\quad\mbox{ and }\quad\mathbf{A}:=\Big{(}a^{\gamma}_{h}(\psi_{j},\psi_{i})\Big{)}.$$

$${\int }_{-1}^1{P}_k(x)e^{i\omega x}dx=i^k{\left(\frac{2\pi }{\omega }\right)}^{1/2}{J}_{k+\frac{1}{2}}(\omega ),$$

$$\rho {|}_K\approx \sum _r\sum _s{\rho }_{rs}{P}_r(x)P_s(y),$$

$$\left|\int g(x)e^{i\omega x}dx-\int \sum a_k{P}_k(x)e^{i\omega x}dx\right|\le {\Vert g-\sum a_k{P}_k\Vert }_1.$$

$${\tilde{\Theta }}_K^{\ast \ast }\leftarrow \text{textsc}PhaseSep({\Theta }_{K,\mathrm{def}},{\tilde{\Theta }}_K^{\ast },\epsilon +2\delta ).$$

Full text

1 Introduction

In this paper, we consider the acoustic wave equation in inhomogeneous media given by

[TABLE]

with initial conditions $u(x,0)=\sum_{l=1}^{L}A_{l}(x)e^{i\omega\phi_{l}(x)}$ and $u_{t}(x,0)=\sum_{l=1}^{L}i\omega B_{l}(x)e^{i\omega\phi_{l}(x)},$ where $c(x)>0$ is the wave speed in the medium, $T>0$ is a given time and $\omega$ is the frequency of the input. Due to the use of the theory of geometrical optics in the construction of basis functions, we assume that the initial condition $u(x,0)$ is a superposition of functions in the form $A(x)e^{i\omega\phi(x)}$ for some smooth functions $A(x)$ and $\phi(x)$ independent of $\omega$, and the initial condition $u_{t}(x,0)$ is a superposition of functions in the form $i\omega B(x)e^{i\omega\phi(x)}$ where the amplitude function $B(x)$ is also a smooth function independent of the frequency $\omega$. We note that, for more general cases, one can use the technique of micro-local analysis to express the initial conditions as a superposition of functions in the form $A(x)e^{i\omega\phi(x)}$ [32, 33]. We are interested in high-frequency solutions propagating in a medium with smooth wave speed $c(x)$, in the sense that the variation within each cell of an underlying domain partition is small. Thus, we will assume that the frequency $\omega\gg 1$ is a very large number.

1.1 Our work

We consider high-frequency solutions of the problem (1) given by a superposition of wave components:

[TABLE]

where $\omega\gg 1$ is the base frequency and $\phi_{k}$ are the phase functions that satisfies

[TABLE]

In our method, we assume that the expansion (2) holds in each cell of an underlying partition of the computational domain, and the number of waves $N$ can vary from cell to cell. One can use evenly-distributed ray directions together with plane wave type basis functions. For good approximations, one needs to use many ray directions but this makes the Galerkin method inefficient and ill-conditioned.

In this paper, we consider a modified wavefront-tracking method[34, 29, 5] to capture possible phases in the solution before the actual simulation. Since the number of phases can be large in general, we assume that there is a small number of phases that contribute most to the solution. Then, we apply a clustering procedure and a dimension reduction procedure to this set of all possible phases and obtain the dominant phases $\phi_{k}(x)$ in the solution. Using these dominant phases, we define our basis functions as $A_{k}(x)e^{i\omega\phi_{k}(x)}$, where $A_{k}$ are polynomials. We remark that these basis functions are obtained for each cell in the partition of the domain. Finally, we use these basis functions together with an interior penalty discontinuous Galerkin (IPDG) method and obtain the approximate solution by solving the resulting linear system. We note that the basis functions can also be used for other DG discretizations such as [8, 9]. We also remark that the wavefront-tracking method can be implemented in a parallel way effortlessly. Moreover, we propose an online/offline scheme to accelerate the simulations for a given medium with many initial conditions, and this situation arises in many practical applications involving inversions. Notice that we will consider smooth media in this paper. This means that the computational mesh is fine enough to capture the variations in the media. For media with more oscillations that cannot be captured by the computational grid, one has to apply some type of multiscale ideas, such as [7, 15, 16, 10]. By using our proposed scheme, we are able to compute the solution of the acoustic wave equation in the high-frequency regime in a very efficient way. Our numerical results show that our scheme is robust on the frequency, i.e. the errors do not grow significantly with the frequency. Furthermore, we compare our scheme with the IPDG scheme that uses the standard polynomial basis. We observe that, with about the same number of unknowns, our scheme performs much better than standard schemes. To the best of our knowledge, our method is the first one that combines ray-based idea and Galerkin method to solve time-dependent wave equations with high frequency solutions.

1.2 Related works

In literature, there are some works that focus on the closely related problem given by the Helmholtz equation

[TABLE]

with frequency $\omega\gg 1$.

On one hand, due to the Shannon’s sampling theorem, the minimum requirement for the number of degrees of freedom of representing the solution is of $\mathcal{O}(\omega^{d})$. On the other hand, the pollution effect haunts the standard Galerkin methods using low order polynomials as the basis [2, 22], i.e. the $L^{2}$-error grows as $\omega$ getting large while keeping $\omega h$ fixed. A common strategy for solving this problem is to incorporate oscillatory functions into the basis for the Galerkin methods. This approach significantly reduces the number of degrees of freedom required for accurately representing the solution.

Mainly, there are two types of methods using this strategy. One type of methods does not assume knowledge of the propagation directions of the wave-field. Instead, they incorporate analytic or approximate solutions with predefined directions into the basis. For example, the generalized finite element method [1, 30] uses basis function given by the product of plane waves with uniformly-spaced directions and finite element basis functions. The Trefftz methods use local solutions of the Helmholtz equation as the basis functions[20, 11]. In particular, a common choice of Trefftz basis for a homogeneous medium is the plane waves, for example, as used in the plane wave discontinuous Galerkin method [17, 19, 21] and the discontinuous enrichment method [14]. For an inhomogeneous medium, some recent papers have proposed methods to construct the basis consisting of approximate local solutions of the form $e^{P(x)}$ with a complex polynomial $P$[24, 25, 23].

Another type of methods relies on an asymptotic approximation, i.e. the WKB ansatz. This ansatz assumes that a high-frequency solution of the Helmholtz equation (3) can be locally approximated by a superposition of wave components,

[TABLE]

where $A_{k}$ and $\phi_{k}$ are non-oscillatory and independent of $\omega$. The function $A_{k}$ and $\phi_{k}$ are called the amplitude and phase, respectively. A comprehensive survey of the techniques arising from this ansatz can be found in [12]. These methods incorporate the phases $\phi_{k}$ in the basis of a Galerkin method and usually called phase-based methods. One fundamental difficulty of these methods is how the phases $\phi_{k}$ should be approximated. For example, a phase-based hybridizable discontinous Galerkin (HDG) method is proposed in [31]. This method uses basis consisting of the product of polynomials and oscillatory functions of the form $e^{i\omega\phi(x)}$, where the phases $\phi$’s are obtained from the eikonal equation or ray-tracing. Later, [28] analyzes the $h$-convergence of a phase-based IPDG method and shows that, under some conditions, these basis functions have the same approximation power as the polynomials in the product. Instead of incorporating the global phases, [4] proposed to use the product of polynomials and plane wave with dominant direction as the basis function for an inhomogeneous medium. Numerical results show that such basis performs significantly better than uniformly-spaced plane waves regarding efficiency and accuracy. Recently, a ray-based finite element method (ray-FEM) is proposed to learn and incorporate dominant directions into the basis in a stable, efficient and systematic way[13]. Numerical tests show that this method achieves asymptotic convergence as $\mathcal{O}(\omega^{-\frac{1}{2}})$ when $\omega\to\infty$.

This paper is organized as follows. In Section 2, we give an overview of our proposed method. In Section 3, we give the detailed algorithms of our method. In addition, we discuss some computational issues, including parallel computation and conditioning of the resulting linear system. In Section 4, we present some numerical experiments to study the behavior of the error as $\omega\to\infty$ and also mesh size $h\to 0$. Finally, we conclude the paper in Section 5.

2 The ray-construction based IPDG method

In this section, we present the key ingredients of our ray-construction based IPDG method for the high frequency acoustic wave propagation problem (1). The algorithmic details are presented in Section 3. The main idea of our approach is to approximate the wave-field via a superposition of plane waves with only dominant ray directions in the solution during the simulation. We will combine the idea of wavefront tracking and a dimensional reduction procedure to construct the dominant ray directions in the solution. We then incorporate these dominant ray directions to form plane wave basis functions for an IPDG method to improve the accuracy and stability for computation of high-frequency solutions.

Consider a cell of an underlying partition of the domain and assume the ansatz (2) in this cell. Note that $N$ is the number of phases and can be large in general. It also varies in time. Recall that the amplitude $A_{k}$ and phase $\phi_{k}$ are independent of the frequency $\omega$. We take a point $x_{0}$ in the center of the cell with size $h$. Using Taylor expansions of $\phi_{k}$ and $A_{k}$ in $x$ around the point $x_{0}$, we have

[TABLE]

for $\left|x-x_{0}\right|<h\ll 1$. When $\omega\to\infty$ and $\omega h=\mathcal{O}(1)$, the asymptotic error decreases as $\mathcal{O}(\omega^{-1})$. This suggest that we can represent each component of the solution locally at $x_{0}$ by the product of a linear function and a plane wave $e^{i\omega\nabla\phi_{k}(x_{0},t)\cdot(x-x_{0})}$. Moreover, we essentially need $\nabla\phi_{k}(x_{0},t)$ to determine the polynomial-modulated plane wave in the local approximation. In this paper, our aim is to find the ray directions $\nabla\phi_{k}(x_{0},t)$ and determine the dominant ray directions within the set of ray directions.

Now, suppose there is a small error $\varepsilon$ in the vector $\nabla\phi_{k}(x_{0},t)$. Then the asymptotic error in (5) becomes $\mathcal{O}(h^{2}+\omega h\varepsilon+\omega h^{2}+\omega^{-1})$. Also, when $\omega\to\infty$ and $\omega h=\mathcal{O}(1)$, the asymptotic error decreases as $\mathcal{O}(\varepsilon+\omega^{-1})$. Our approach is to replace all the $\nabla\phi_{k}(x_{0},t)$ by some approximations, which are easy to compute. In exchange, the number of degrees of freedom could be significantly reduced and the resulting system would be better conditioned. Despite the extra $\mathcal{O}(\varepsilon)$ term in the ansatz, our numerical results in Section 5 confirm that our method still benefits a lot from such approximation.

To simplify our discussion, from now on we restrict our discussion to the domain $\Omega=(0,1)^{2}$, and assume that the solution $u(x,t)$ satisfies the periodic boundary condition. We divide the domain $\Omega$ into $N\times N$ square cells $\mathcal{T}_{h}$, given by $K_{ij}:=[\frac{i-1}{N},\frac{i}{N}]\times[\frac{j-1}{N},\frac{j}{N}],$ for $1\leq i,j\leq N$, where $h=1/N>0$ is the mesh size. We let $\mathcal{F}_{h}$ be the set of faces of the partition. For a given cell $K\in\mathcal{T}_{h}$, we define $x_{K}$ to be the centroid of $K$. We call these points the observation points. For each $x_{K}$, we define the set $\Theta_{K}$ of all possible ray directions by

[TABLE]

where $N_{K}$ is the number of phases, which depends on the cell $K$ and also the time $t$. We remark that we assume the solution contains a finite, and possibly large, number of phases at any time instant.

Next, we use the symmetric interior-penalty type discontinuous Galerkin (IPDG) method [18] to compute the solution of (1). For $j=1,2,3,4$, let $v_{j,K}$ be the four vertices of $K$, and $\varphi_{j,K}$ be the standard Lagrange-type bilinear basis on $K$ such that $\varphi_{j,K}(v_{i,K})=\delta_{ij}$, where $\delta_{ij}$ is the Kronecker delta. We define an enriched local approximation space by

[TABLE]

where functions in $V(\Theta_{K})$ are defined on the cell $K$ only. Then we can define the global approximation space

[TABLE]

Following the derivation of the standard IPDG method, we can write down the following semi-discrete scheme: find $u_{h}\in V_{c}$ such that

[TABLE]

for any $v_{h}\in V_{c}$, where $(u,v)=\int_{\Omega}c^{-2}\,u\overline{v}\,dx$ and the bilinear form $a^{\gamma}_{h}$ is given by

[TABLE]

for any $u,v\in V$, where $\gamma>0$ is the penalty parameter. Here $[\,\cdot\,]$ and $\{\,\cdot\,\}$ are the usual jump and average operators in discontinuous Galerkin methods, which are defined in the following way. Let $K^{\pm}$ be two tiles sharing an edge $F$, $n^{\pm}$ be the outward normal of $K^{\pm}$ on $F$ and $u^{\pm}$ be the two smooth scalar functions on $K^{\pm}$. The average operator $\{\,\cdot\,\}$ and jump operator $\left[\,\cdot\,\right]$ are given by

[TABLE]

respectively. Similarly, for smooth vector fields ${\boldsymbol{\sigma}}^{\pm}$ defined on $K^{\pm}$, respectively, the average operator $\{\,\cdot\,\}$ and jump operator $\left[\,\cdot\,\right]$ are given by

[TABLE]

We note that the set $\Theta_{K}$ contains a continuum of ray directions, and thus cannot be used directly for computations. With this in mind, we will construct a finite set $\widetilde{\Theta}_{K}$, which is a subset of $\Theta_{K}$ and contains all dominant ray directions within the cell $K$. For this purpose, we will consider a partition of $[0,T]$ with time step $\Delta t$, and we denote $t_{m}=m\Delta t$, $m=0,1,\cdots$. Then we consider the set

[TABLE]

We note that the set $\widehat{\Theta}_{K}$ is finite and contains all ray directions resolved by the partition in time. We emphasize that this set has a large dimension, and cannot be used directly in simulations. One key ingredient of our method is to perform a dimension reduction for the set $\widehat{\Theta}_{K}$ and obtain a subset $\widetilde{\Theta}_{K}$ with a much smaller dimension. In particular, the set $\widetilde{\Theta}_{K}$ has the form

[TABLE]

We will give a detail discussion on how to obtain this set in the next section. Next, we define an enriched local approximation space by

[TABLE]

where functions in $V(\Theta_{K})$ are defined on the cell $K$ only. Then we can define the global approximation space

[TABLE]

Using the space $V$, we can write down the fully-discrete scheme: Given $u^{0}_{h}\in V$ and $u^{1}_{h}\in V$, for $n\Delta t<T$, we find $u_{h}^{n+1}\in V$ such that

[TABLE]

for any $v_{h}\in V$. We remark that, by using a standard stability analysis, the time step $\Delta t$ should be chosen such that $\Delta t\|\mathbf{A}\|_{2}<1$, where $\mathbf{A}$ is the stiffness matrix resulting from the bilinear form $a_{\gamma}$ and $\|\cdot\|_{2}$ denotes the $2$-norm.

The above gives a general outline of our scheme. In the next section, we present the detailed implementations.

3 Algorithms

In this section, we present the full algorithm for the ray based IPDG method. We divide it into two conceptual stages:

(ray-construction stage) determining the dominant rays in the solution, i.e. the set $\widetilde{\Theta}_{K}$ at every observation point $x_{K}$. 2. 2.

(time marching stage) solving the IPDG system using the approximate space $V(\widetilde{\Theta}_{K})$.

For the phase-construction stage, we use a modified version of the wavefront tracking method proposed in [34] to determine all ray directions of the solution throughout $0\leq t\leq T$. In particular, we will obtain approximations to the ray directions at the observation point $x_{K}$ for the cell $K$ in the domain partition. We call the set of all these ray directions $\widetilde{\Theta}_{K}^{*}$. Afterward, we apply a clustering procedure on $\widetilde{\Theta}_{K}^{*}$ and obtain our desired set $\widetilde{\Theta}_{K}$ of dominant ray directions.

In the following Sections 3.1–3.4, we discuss the procedures for the phase-construction stage. In Section 3.5, we integrate the ideas from Sections 3.1–3.4 into an algorithm for the ray-construction stage. In Section 3.6, we write down an algorithm for the time marching stage. In Section 3.7, we propose an online/offline scheme based on the ray-construction based IPDG method.

3.1 Wavefront propagation

In this procedure, we compute the wavefront propagation and use it to determine the ray directions. First of all we define wavefronts, which are essentially the level sets of the phase functions, together with the gradients of the phase functions on the level sets. More precisely, we let $\phi$ be one of the phases $\phi_{k}$ in the ansatz (2). A wavefront $W(t_{0},\alpha)$, at a fixed time instant $t_{0}$, corresponding to the phase function $\phi(x,t_{0})$, is essentially the level set curve $\phi(x,t_{0})=\phi_{0}$ for a given constant $\phi_{0}$, which may depend on $t_{0}$. Mathematically, the wavefront $W(t_{0},\alpha)$ is a smooth curve $(y(t_{0},\alpha),q(t_{0},\alpha))\in\mathbb{R}^{2}\times\mathbb{R}^{2}$, parametrized by $\alpha\in J$, such that for any $\alpha$,

[TABLE]

where $J$ denotes an interval of real numbers. That is, the component $y(t_{0},\alpha)\in\mathbb{R}^{2}$ lies on the level set curve $\phi(x,t_{0})=\phi_{0}$, and that the level set curve $\phi(x,t_{0})=\phi_{0}$ is represented by the function $y(t_{0},\alpha)$ and is parametrized by $\alpha$. In addition, the component $q(t_{0},\alpha)\in\mathbb{R}^{2}$ is the gradient vector of the function $\phi(x,t_{0})$ at the point $y(t_{0},\alpha)$.

A related concept is a ray. By inserting $u=A(x,t)e^{i\omega\phi(x,t)}+\mathcal{O}(\omega^{-1})$ into (1) and considering the leading order term, we obtain the following eikonal equation for the phase $\phi(x,t)$:

[TABLE]

A ray $r(t)$ is a bicharacteristic pair $(x(t),p(t))\in\mathbb{R}^{2}\times\mathbb{R}^{2}$ related to the Hamiltonian $H(x,p)=c(x)\left|p\right|$ which is characterized by the following differential equations,

[TABLE]

The functions $x(t)$ and $p(t)$ are called the position and ray direction of the ray $r(t)$ respectively. Rays and wavefronts are related in the following way. If a ray $(x(t),p(t))$ lies on a wavefront $W$ initially, that is, $(x(0),p(0))=W(0,\alpha_{0})\mbox{ for some }\alpha_{0},$ then the ray lies on the same wavefront for any $t>0$, that is, $(x(t),p(t))=W(t,\alpha_{0}).$

We will approximate a wavefront $W$ at the $n$-th time step $t_{n}$ by a finite set of points in $\mathbb{R}^{2}\times\mathbb{R}^{2}$, which is called a discrete wavefront and is denoted by $W^{n}$. We write $W^{n}=\{r^{n}_{j}\}$. We will parametrize $W^{n}$ by a discrete set of parameters $\{\alpha_{j}\}$. For example, using the notations in (7), we will take a discrete set $\{\alpha_{j}\}\subset J$ and define

[TABLE]

The wavefronts at the initial time $t=0$ are given by the initial condition of the problem (1). In particular, we will assume that the phase $\phi(x,0)$ and the ray direction $\nabla\phi(x,0)$ are known at the initial time. We will use a set of level curves of $\phi(x,0)$ to construct the wavefronts at the initial time. The number of these level curves is user-defined. We will assume that the number of these level curves is large enough so that the ray directions of the initial condition are well resolved by the partition of the computational domain. For each of these level curves, we will use (7) to define a wavefront $W^{0}$. For more precise constructions, see Section 4 where numerical examples are presented.

Suppose that all wavefronts are computed at the $n$-th time step. Let $W^{n}=\{r^{n}_{j}\}=\{(y(t_{n},\alpha_{j}),q(t_{n},\alpha_{j}))\}$ be a given wavefront at the time $t_{n}$. To advance the wavefront $W$ to the $(n+1)$-th time step $t_{n+1}$, we treat each $r^{n}_{j}$ as a ray at time $t_{n}$, and then take this as the initial condition of (9) and compute the solution, via a fourth order Runge Kutta scheme, for one time step (that is, from the time $t_{n}$ to the time $t_{n+1}$). We call this solution as $r^{n*}_{j}$, and this will be taken as part of the wavefront at the time $t_{n+1}$. In other words, $r^{n*}_{j}\approx W\left(t_{n+1},\alpha_{j}\right)$. The pseudocode is given in Algorithm 1. Note that, the wavefront at the $(n+1)$-st step will be obtained after the reconstruction procedure presented next.

3.2 Wavefront reconstruction

In this procedure, we insert new rays into a discrete wavefront to maintain the quality of the wavefront. We choose a tolerance function $\mbox{tol}\left(x,p\right)=\alpha_{1}\left|x\right|+\alpha_{2}\left|p\right|$ for some constants $\alpha_{1}$ and $\alpha_{2}$. Whenever $\mbox{tol}(r^{n*}_{j}-r^{n*}_{j+1})\geq 1$ for a pair of neighboring rays $r_{j}^{n*}$ and $r_{j+1}^{n*}$, we insert $\lfloor\mbox{tol}(r^{n*}_{j}-r^{n*}_{j+1})\rfloor$ new rays between them. We take the new rays as equidistant linear interpolation of the corresponding neighboring rays.

By connecting each pair of neighboring rays $r_{j}^{n*}$ and $r_{j+1}^{n*}$ by a line segment, we obtain a linear interpolant in $\mathbb{R}^{4}$. We order all $r^{n*}_{j}$ and the inserted rays along this spline and obtain $W^{n+1}=\{r^{n+1}_{j}\}$. We use $j^{*}$ to denote the new index at the $(n+1)$-th time step corresponding to $r^{n*}_{j}$. In other words, $r^{n*}_{j}$ and $r^{n+1}_{j^{*}}$ are equal. It is clear that $j^{*}\leq\ell\leq(j+1)^{*}$ for any $r^{n+1}_{\ell}$ inserted between $r^{n+1}_{j^{*}}$ and $r^{n+1}_{(j+1)^{*}}$. The pseudocode of this procedure is shown in Algorithm 2. Note that, in this algorithm, we also compute

[TABLE]

which will be used in the next procedure.

3.3 Ray determination

In this procedure, we determine the ray directions of the solution at an observation point when an approximate wavefront passes through the point. Let $x_{K}$ be an observation point. To check whether a wavefront $W$ is passing through $x_{K}$ between the $n$-th and $(n+1)$-th time step, we consider the corresponding discrete wavefronts at the time steps, namely, $W^{n}$ and $W^{n+1}$. We form triangles in $\mathbb{R}^{2}$ using the position of the rays in $W^{n}$ and $W^{n+1}$ and check whether these triangles contains $x_{K}$. More precisely, we write $W^{n}=\{(x^{n}_{j},p^{n}_{j})\}$ and consider the triangles $x^{n}_{j}x^{n+1}_{j^{*}}x^{n+1}_{j^{*}+\ell+1}$ for $\ell\leq(j+1)^{*}-j^{*}-1$, and $x^{n}_{j}x^{n}_{j+1}x^{n+1}_{(j+1)^{*}}$. We call each of these triangles a ray cell.

Whenever $x_{K}$ lies in one of the ray cells, we approximate the ray direction at $x_{K}$ by linear interpolation. Suppose the rays corresponding to the three vertices of the ray cell are given by $r_{1},r_{2},r_{3}$. Then we approximate the ray direction by

[TABLE]

where $\lambda_{\ell}\in\mathbb{R}$ is the barycentric coordinates of $x_{K}$ in the triangle $x_{1}\,x_{2}\,x_{3}$, $x_{\ell}$ and $p_{\ell}$ are the positions and ray directions of the rays $r_{\ell}$ respectively for $\ell=1,2,3$. The barycentric coordinates can be obtained by solving

[TABLE]

We name the set of these approximations to the ray directions at $x_{K}$ by $\widetilde{\Theta}_{K}^{*}$.

Note that the time complexity of Algorithm 3 is $\mathcal{O}(mn)$ for $n$ observation points and $m$ ray cells. However, it is possible to reduce the time complexity by using some heuristics. Supposing that the ray cells are adequately small in size and regular in shape, we may assume that an observation point $x_{K}$ lies in a ray cell only if all the vertices of a ray cell lies in $K$. Then if one of the vertices of a ray cell is $x$, the only possible observation point that lies in the ray cell is given by $\left(\frac{\lfloor Nx\rfloor+\frac{1}{2}}{N},\frac{\lfloor Ny\rfloor+\frac{1}{2}}{N}\right)$, reducing the time complexity to $\mathcal{O}(m)$. Similarly, one can reduce the time complexity of Algorithm 3 for a non-rectangular computational domain to $\mathcal{O}(m)$ whenever there is a region classifier, which maps a points to a region, of time complexity $\mathcal{O}(1)$.

3.4 Ray separation

In this procedure, we find dominant ray directions in each $\widetilde{\Theta}^{*}_{K}$ after all the approximate wavefronts have propagated till time $T$. The aim of this procedure is to obtain a set $\widetilde{\Theta}_{K}$ such that the deviation $d(\widetilde{\Theta}_{K},\widetilde{\Theta}^{*}_{K})$ is not too large while the minimum distance between any two distinct elements in $\widetilde{\Theta}^{*}_{K}$, or namely the separability, is not too small. We will see that using the new set of ray directions $\widetilde{\Theta}_{K}$ will improve the numerical stability in the IPDG formulation, while the approximation power of the new basis do not deteriorate significantly.

We use a parameter $\varepsilon>0$ to control the deviation and the separability. The idea of this procedure is to cover the set $\widetilde{\Theta}^{*}_{K}$ by some balls with radius $\varepsilon$, and make sure the centers of the balls are far enough from each other. We take the dominant ray directions $\widetilde{\Theta}_{K}$ by the centers of these balls. We will also need a predefined set $\Theta_{K,\mathrm{def}}$ to be included in $\widetilde{\Theta}_{K}$ obtained from this procedure. One usage of this predefined set is to include the ray directions corresponding to the initial solutions.

We illustrate this procedure in Figure 2. Let $B_{\varepsilon}(x_{0})$ be the Euclidean ball with radius $\varepsilon$ centered at $x_{0}$. To begin with, we take $\widetilde{\Theta}_{K}$ to be $\Theta_{K,\mathrm{def}}$ and eliminate any ray direction in $\widetilde{\Theta}^{*}_{K}$ that also lies in $\Theta_{K,\mathrm{def}}+B_{\varepsilon}(0)$. Then we repeatedly do the following until $\widetilde{\Theta}^{*}_{K}$ is empty. Firstly, we take a ray direction $p$ from $\widetilde{\Theta}_{K}$. Secondly, we compute the centroid $\bar{q}$ of the set $\widetilde{\Theta}^{*}_{K}\cap B_{\varepsilon}(p)$, and include $\bar{q}$ in $\widetilde{\Theta}_{K}$. Finally, we eliminate any ray direction in $\widetilde{\Theta}^{*}_{K}$ that also lies in $B_{\varepsilon}(\bar{q})$ and repeat. The pseudocode is shown in Algorithm 4.

Next, we discuss some properties of the set $\widetilde{\Theta}_{K}$ obtained from Algorithm 4. We say that a finite set $\Theta\subset\mathbb{R}^{2}$ is $\varepsilon$-separable if the distance between any two distinct elements in $\Theta$ is at least $\varepsilon$. The following theorem gives an estimate for the minimum separation of the set $\widetilde{\Theta}_{K}$.

Theorem 1**.**

If $\Theta_{K,\mathrm{def}}$ is an $(\varepsilon/2)$-separable set, then the set $\widetilde{\Theta}_{K}$ obtained in Algorithm 4 is an $(\varepsilon/{2})$-separable set.

Proof.

Consider $p=(q,0)^{T}$, $q>\varepsilon$, $\left|p_{j}\right|>\varepsilon$ and $\left|p_{j}-p\right|<\varepsilon$ for $j=1,\ldots,m$. Then $p_{j}$ lies in the half-plane $H:=\{(x,y):x>\varepsilon/2\}$. Since $H$ is convex, the convex hull of $\{p_{j}\}_{j=1}^{m}$, which contains $\frac{1}{m}\sum p_{j}$, also lies in $H$. The result follows by translation and rotation of the points in the above argument. ∎

Note that if we replace line 6 in Algorithm 4 by ‘$\bar{q}\leftarrow p$’, this modified algorithm produces an $\varepsilon$-separable set given that $\Theta_{K,\mathrm{def}}$ is $\varepsilon$-separable. Empirically, the set $\widetilde{\Theta}_{K}$ will contain more ray directions in this modified version while the minimum separation of the set of ray directions of the original version is close to $\varepsilon$.

3.5 The ray-construction stage

In this section we combine Algorithms 1–4 to determine the set of dominant ray directions in each cell $K$ throughout the simulation, namely, $\widetilde{\Theta}_{K}$. We give the pseudocode in Algorithm 5. Here we describe the parameters required for this procedure. The parameter $c$ and $T$ are the velocity function and the time of simulation defined in the problem (1), respectively; the parameter $\Delta t$ is the time step in the wavefront propagation; the parameter $\{x_{K}\}$ is the set of observation points in each cell; the parameter $\mbox{tol}(x,p)$ is the tolerance function which controls how many points should be inserted between two neighboring rays in the wavefront reconstruction procedure; the parameter $\{\Theta_{K,\mathrm{def}}\}$ is the collection of the default set of ray directions used in the phase separation procedure; the parameter $\{r^{0}_{j}\}$ is the collection of discrete wavefront corresponding to the initial ray directions; and the parameter $\varepsilon$ controls the deviation and separability in the phase separation procedure. The return sets $\widetilde{\Theta}_{K}$ are our desired dominant ray directions on $K$.

3.6 The time-marching stage

In this stage, we obtain the solution of (1) by solving an IPDG system as described in Section 2. The pseudocode of this stage is given in Algorithm 6. The basis is constructed based on the ray directions obtained in the ray-construction stage. The penalty parameter in the IPDG scheme is denoted by $\gamma$.

We note that it is very time consuming for computing the mass and the stiffness matrices:

[TABLE]

Instead of using numerical quadratures, an alternative method is to use the following formula given by [3]:

[TABLE]

where $J$ is the Bessel function of the first kind and $P_{k}(x)$ is the Legendre polynomial of degree $k$. By writing the integrands into a product of functions in $x$ and in $y$, all the nonzero entries of $\mathbf{A}$ can be computed by equation (10). For the matrix $\mathbf{M}^{c}$, if the mesh size $h$ is small enough such that $\rho(x,y):=1/c^{2}(x,y)$ is approximately given by

[TABLE]

then we can approximate the nonzero entries using (10) with an error estimate of the form

[TABLE]

3.7 Online/offline scheme

Suppose we need to solve the problem (1) with a fixed velocity profile but many initial conditions. From the constructions above, we see that the set of basis functions depends on the pre-defined sets $\{\Theta_{K,\mathrm{def}}\}$, which depend on the initial conditions (c.f. Section 3.4). It is obvious that one does not want to compute the mass and the stiffness matrices for each initial condition. In this section, we propose an offline/online scheme to address this issue by pre-computing the entries of one mass and one stiffness matrix that can be used for all choices of initial conditions. That is, these mass and stiffness matrices are independent of the initial conditions.

First of all, we define a predefined set of ray directions $\Theta_{\mathrm{pre}}$. We assume that this set contains a large set of ray directions that are needed to compute the solution for all initial conditions with good accuracy. Note that this set can be considered as a snapshot space using the terminology in [16]. It is obvious that one can take $\Theta_{\mathrm{pre}}=\mathbb{R}^{2}\backslash\{0\}$, but this choice is not good in practice as it is too large and too expensive to use in the offline stage. Instead, we put a restriction on our solver such that it only handles the case that the ray directions lie in an annulus $R=\{p\in\mathbb{R}^{2}:0<p_{1}<p<p_{2}\}$. Then we choose the set of predefined phases that satisfies $R\subset\Theta_{\mathrm{pre}}+B_{\delta/2}(0)$, for some small $\delta>0$. Next, we present the offline and the online stages. In the offline stage, we compute the matrices $\mathbf{M^{c}}$ and $\mathbf{A}$ in Algorithm 6 using the ray directions in $\Theta_{\mathrm{pre}}$. We note that this computation is performed before the actual simulations and is independent of the initial conditions. We also note that the cost of this offline computation depends on the size of the set $\Theta_{\mathrm{pre}}$.

In the online stage, when the initial condition is given, we will apply a modification of the phase-construction procedure as in Algorithm 5 to obtain the desired ray directions. The modification is discussed as follows. The main idea is that, we will construct the dominant ray directions as before, but we will use the ray directions in the set $\Theta_{\mathrm{pre}}$ that are closest to those dominant ray directions as the basis functions. The advantage is that we can extract sub-matrices of $\mathbf{M^{c}}$ and $\mathbf{A}$ from the offline stage for the simulations. Next, we present in detail how this is performed. Firstly, we apply the original phase separation step with $\varepsilon$ replaced by $\varepsilon+2\delta$.

[TABLE]

Secondly, for each ray direction in $\widetilde{\Theta}^{**}_{K}$, we choose the closest ray direction in $\Theta_{\mathrm{pre}}\cup\Theta_{K,\mathrm{def}}$ to form the set $\widetilde{\Theta}_{K}$. More precisely we take

[TABLE]

Given that $d(\widetilde{\Theta}^{**}_{K},\Theta_{\mathrm{pre}})<\delta/2$, the set $\widetilde{\Theta}_{K}$ is an $(\varepsilon/2)$-separable set such that $d(\widetilde{\Theta}_{K},\widetilde{\Theta}^{*}_{K})<\varepsilon+3\delta$. After we obtain this new $\widetilde{\Theta}_{K}$, we retrieve the corresponding entries from the precomputed $\mathbf{M}^{c}$ and $\mathbf{A}$. Besides, we only need to compute the entries corresponding to $\Theta_{K,\mathrm{def}}\setminus\Theta_{\mathrm{pre}}$ in the online stage. Afterwards, we solve the IPDG scheme as in Algorithm 6.

3.8 Remarks on computation issues

3.8.1 Parallel computation for the ray-construction stage

Suppose we have multiple processing units available for simultaneous calculations. We note that the ray direction and position of a ray do not depend on other rays but only the velocity field. In other words, the computation for the rays can be carried out with duplicated velocity field on each processing unit. Suppose we split a discrete wavefront $W^{n}$ into smaller wavefronts (with overlapping endpoints), solve the ray equations for each smaller wavefront, apply the reconstruction WavefrontRecon, and combine the reconstructed wavefronts into one wavefront, the resulting wavefront is the same as the discrete wavefront $W^{n+1}$ describe in section 3.2. Therefore if we modify Algorithm 5 and perform the procedures WavefrontPropagate, WavefrontRecon and PhaseDet for each smaller wavefronts on different processing units, the resulting set of phases $\widetilde{\Theta}_{K}$ would be the same after performing the PhaseSep procedure. The speedup of this parallel algorithm depends on the maximum number of rays on each wavefront throughout the simulation, which may not be clear before the computation. How to choose the initial splitting of the wavefronts remains a question.

3.8.2 Conditioning of the IPDG system

In this part, we consider the conditioning of the mass matrix. As an illustration, we consider $K=[0,h]^{2}$, $\widetilde{\Theta}_{K}=\{(1,0)^{T},(\cos\theta,\sin\theta)^{T}\}$, and the local approximation space $V(\widetilde{\Theta}_{K})$. In Figure 3, for each mesh size $h$, we compute the condition number of the mass matrix corresponding to the basis $\varphi_{j,K}e^{ip\cdot(x-x_{0})}$, where $\varphi_{j,K}$ is a standard bilinear basis, $p\in\widetilde{\Theta}_{K}$, and $x_{0}$ is the center of $K$. The figure suggests that the resulting system is very ill-conditioned for small $h$ and $\theta\approx 0$.

To form a better conditioned system, in addition to choosing a larger $\varepsilon$, we also consider the proper orthogonal decomposition (POD) basis for the IPDG system [26, 27]. The POD basis is obtained by considering a truncated singular value decomposition of the matrix $\mathbf{M}^{c}=\left(\int\frac{1}{c^{2}}\psi_{i}\psi_{j}\,dx\right)_{ij}$, where the $\psi$’s are the basis in $V(\widetilde{\Theta}_{K})$. Suppose the eigenvalues of $\mathbf{M}^{c}$ are given by $\lambda_{1}\geq\lambda_{2}\ldots\geq\lambda_{N}>0$. For a fixed $\eta$, we choose the minimum $N^{*}$ such that

[TABLE]

Then the POD basis is given by the eigenvectors corresponding to $\lambda_{1},\ldots,\lambda_{N^{*}}$.

4 Numerical experiments

In this section, we present some numerical examples to show the performance of our proposed method. For all the cases, the computational domain is $\Omega:=[0,1]^{2}$, and $\Omega$ is subdivided into $N\times N$ square cells. We define the mesh size $h$ to be $\frac{1}{N}$. We simulate the equation (1) with non-constant wave speed $c(x)$ for the times $0\leq t\leq 1:=T$ using the time step $\Delta t$. The equation is supplemented with the periodic boundary condition and the initial condition will be specified for each of the cases considered below. In all of our examples, we consider two different media with speeds given by

[TABLE]

[TABLE]

We illustrate these wave speed profiles in Figure 4.

To benchmark the performance of our method, we will compare the numerical solutions computed by our method to a reference solution computed by a pseudospectral scheme for a finer mesh.

4.1 Example 1

In this example, we test our method without performing the POD described in Section 3.8.2 and compare it to the standard IPDG method (with polynomial basis) as $\omega\to\infty$ but $\omega h=\mathcal{O}(1)$. We consider the wave speed $c_{1}$. For the initial conditions, we take $L=1$, the phase function $\phi_{1}(x)=x_{1}$ and the amplitude functions $A_{1}(x)=1$ and $B_{1}(x,\omega)=-i\omega$. For the numerical solutions $u_{h}^{0}$ and $u_{h}^{1}$ at the first two time steps, we can compute them by the given initial conditions and a forward Euler scheme in time respectively.

Since $\nabla\phi_{1}(x)=(1,0)^{T}$ for all values of $x$, in the phase-construction stage in Algorithm 5, we can take $\Theta_{K,\mathrm{def}}=\{(1,0)^{T}\}$ for all cells $K$ in the partition of the domain. For the wavefronts, we will use $9$ level set functions defined by $\phi(x)=\beta$ with $\beta=0.1,0.2,\cdots,0.9$. In addition, we take the tolerance function $\mbox{tol}(x,p)=10\left|x\right|+100\left|p\right|$ and the separation parameter $\varepsilon=0.2$. For both solving the equations (9) and the time marching stage in Algorithm 6, we take the time step $\Delta t=h/100$. Finally, we take the penalty parameter $\gamma=10$.

We consider the numerical solution of our scheme using different mesh sizes $h$ and $\omega$ with $\omega h=\pi$. We list the relative errors in Table 1. We also consider the numerical solutions given by standard IPDG scheme with the same time step size $\Delta t=h/100$ and penalty parameter $\gamma=10$. We take more degrees of freedoms and let $\omega h=\pi/10$. The numerical results are shown in Table 2.

In these numerical results, we see that the relative $L^{2}$-errors of the numerical solutions for our method do not increase significantly as the frequency $\omega$ increases, in contrast to the numerical solutions of the standard IPDG method with bilinear basis. We also see that the condition numbers of the mass matrices without using POD are relatively large.

In Figure 5, we show the number of ray directions we captured at each observation point $x_{K}$. The number of rays in different regions of the domain suggests that the phase-construction stage is able to capture the approximate ray directions in a consistent manner. Also, there are maximum 2 phases at each observation point or 8 degrees of freedom in a tile. Therefore, there are $2\times\sqrt{8}\approx 5.66$ points per wavelength for the numerical solutions of our method.

We also demonstrate the quality of the solution by comparing our solution to the reference solutions in Figure 6 and 7. In Figure 8, we compute the absolute difference between our solutions and the reference solutions. In Figure 9, we compare the real parts of our solutions to the reference solutions at $T=1$ and $x_{2}=0.3$.

4.2 Example 2

In this example, we test the $h$-convergence of our method with the use POD basis presented in Section 3.8.2. We consider the same wave speed, initial condition, initial wavefronts and parameters $\Theta_{K,\mathrm{def}},\Delta t,\mbox{tol},\varepsilon,\gamma$ as in Example 1. Instead of solving the ill-conditioned system as $h\to 0$, we consider the POD system with parameter $\eta=10^{-7}$.

In Table 3 we list the relative errors of the approximate solutions. We also compare the condition numbers and degrees of freedom of the POD systems to those of the original systems. We observe that the POD systems has a better conditioning with number of degrees of freedom about the same with the original systems. Our results suggests that the $h$-convergence of the POD systems is first-order. Also, the relative $L^{2}$-error improves as the $\omega$ increases with fixed $\omega h$.

4.3 Example 3

In this example, we consider more phases in the initial condition and the solution. We will observe the behavior of our approximate solution as $\omega\to\infty$ while keeping $\omega h=\mathcal{O}(1)$, using the POD system.

We consider the wave speed $c_{1}$ in Example 1. For the initial condition, we take $L=2$, $\phi_{1}(x)=x_{1}$, $A_{1}=1$, $B_{1}=-i\omega$, $\phi_{2}(x)=x_{2}$, $A_{2}=1$ and $B_{2}=-i\omega$. For the phase-construction stage in Algorithm 5, since $\nabla\phi_{1}=(1,0)^{T}$ and $\nabla\phi_{2}=(0,1)^{T}$, we take $\Theta_{K,\mathrm{def}}=\{(1,0)^{T},(0,1)^{T}\}$ for every cell $K$ in the partition of the domain. For the wavefronts, we will use the level sets $\phi_{1}=\beta$ and $\phi_{2}=\beta$, where $\beta=0.1,0.2,\cdots,0.9$. All other parameters are taken the same as that for Example 2.

In Table 4 we list the relative errors of the approximate solutions, and we see that our method is robust with respect to the frequency $\omega$. We also compare the condition numbers and degrees of freedom of the POD systems to those of the original systems. We observe that the POD systems has a better conditioning with number of degrees of freedom about the same with the original systems. We see that the relative $L^{2}$-error of the numerical solution for the POD system does not increase significantly as the $\omega$ increase.

In Figure 10, we show the number of phases we captured at each observation point $x_{K}$. The number of rays in different region of the domain suggests that the phase-construction stage capture the approximate phases in a consistent manner. Also, there are maximum 4 phases at each observation point or 16 degrees of freedom in a tile. Therefore there are $2\times\sqrt{16}=8$ points per wavelength in the original systems of our method.

We also demonstrate the quality of the solution by comparing our solution to the reference solutions in Figure 11 and 12. In Figure 13, we compute the absolute difference between our solutions and the reference solutions. In Figure 14, we compare the real parts of our solutions to the reference solutions at $T=1$ and $x_{2}=0.3$.

4.4 Example 4

In this example, we consider a periodic velocity and observe the behavior of our approximate solution as $\omega\to\infty$ while keeping $\omega h=\mathcal{O}(1)$, using the POD system. We consider the wave speed $c_{2}$, and take the initial solution so that $L=1$, $\phi_{1}(x)=x_{1}$, $A_{1}(x)=1$ and $B_{1}(x,\omega)=-i\omega$. All other parameters are taken as before.

In Table 5 we list the relative errors of the approximate solutions. We also compare the condition numbers and degrees of freedom of the POD systems to those of the original systems. We observe that the POD systems has a better conditioning with number of degrees of freedom about the same with the original systems. We see that the relative $L^{2}$-error of the numerical solution for the POD system does not increase significantly as the $\omega$ increase.

We also demonstrate the quality of the solution by comparing our solution to the reference solutions in Figure 16 and 17. In Figure 18, we compute the absolute difference between our solutions and the reference solutions. In Figure 19, we compare the real parts of our solutions to the reference solutions at $T=1$ and $x_{2}=0.3$.

5 Conclusion

We proposed a new method for solving the time-domain acoustic wave propagation problem with smoothly varying inhomogeneous media in the high-frequency regime. Our method is based on plane-wave type basis functions, with the phases carefully computed by using ideas from geometrical optics, wavefront tracking, and dimensional reduction. The numerical results show evidence that the accuracy of the solution of our proposed method is dramatically better than standard IPDG method and the proposed method’s accuracy do not worsen significantly as frequency increases, in contrast to standard IPDG method with polynomial basis. In the future, we plan to design numerical schemes for more complicated heterogeneous and multiscale media by combining the approach in this paper and some recent ideas on multiscale schemes [6].

Acknowledgement

Chung is partially supported by Hong Kong RGC General Research Fund (Project: 14301314) and CUHK Direct Grant for Research 2016-17. Qian is partially supported by NSF grants (1522249 and 1614566).