Geometry-based approximation of waves in complex domains

Abstract

We consider wave propagation problems over 2-dimensional domains with piecewise-linear boundaries, possibly including scatterers. We assume that the wave speed is constant, and that the initial conditions and forcing terms are radially symmetric and compactly supported. We propose an approximation of the propagating wave as the sum of some special space-time functions. Each term in this sum identifies a particular field component, modeling the result of a single reflection or diffraction effect. We describe an algorithm for identifying such components automatically, based on the domain geometry. To showcase our proposed method, we present several numerical examples, such as waves scattering off wedges and waves propagating through a room in presence of obstacles. Software implementing our numerical algorithm is made available as open-source code.

Keywords: wave propagation, surrogate modeling, scattering, geometrical theory of diffraction.

AMS subject classifications: 35L05, 35Q60, 65M25, 78A45, 78M34.

1 Introduction

The discretization of numerical models for the simulation of complex phenomena results in high-dimensional systems to be solved, usually at an extremely high cost in terms of computational time and storage memory. Among these models, wave propagation problems represent an extremely interesting topic: relevant applications can be found, e.g., in the field of array imaging, where acoustic, electromagnetic, and elastic waves in scattering media are modeled by the reflectivity coefficient, which is often unknown. Some examples in this direction can be found in [1, 2, 3, 4], where inverse scattering problems are used to infer the reflectivity of one or more scatterers embedded either in a known and smooth medium, or in a randomly inhomogeneous medium. Another example of application of wave propagation problems is numerical acoustics, where the goal is to simulate the propagation of sound in a room, in presence of obstacles and walls with different absorbing and/or reflecting properties. See, e.g., [5, 6] for some frequency-domain examples of this.

Our focus here are problems in the time domain, whose numerical simulation is expensive, mainly because one needs to use both a fine spatial mesh and a carefully chosen time step in order to satisfy the CFL condition [7, 8]. In the interest of making these simulations feasible, model order reduction (MOR) [9, 10, 11] represents a promising framework, whose goal is to reduce the computational cost of solving the problem of interest.

In this context, it is well known [12] that wave propagation problems are characterized by a slowly decaying Kolmogorov $n$-width. Because of this, classical linear-subspace MOR methods are not able to reproduce the behavior of the wave propagation without relying on a very high-dimensional linear manifold. This makes linear surrogate models unappealing, since they do not yield significant speed-ups. In recent years, methods that rely on nonlinear and/or hybrid space-time approaches have been proposed. See, e.g., [13, 14, 15] and references therein.

In this work, we focus on wave propagation over 2-dimensional spatial domains, possibly including obstacles. We aim at proposing a new approximation framework, which satisfies the main goal of MOR techniques, that is, making the computational cost of the numerical simulations more feasible. We limit our investigation to domains with piecewise-linear boundaries and a constant wave speed. The initial conditions and forcing terms are assumed to be compactly supported and radially symmetric around a “source point”. Under these assumptions, and following ideas similar to those in [5, 6], we propose to construct a surrogate model for the solution of the target problem as a superposition of some special space-time terms, which we call “field components”.

Each field component models a reflection or diffraction effect, and is characterized by:

• a space-time propagation term, related to the free-space radially symmetric solution of the wave equation;

• a spatial indicator function, determining the spatial support of each component;

• a nonlinear spatial term encoding the angular modulation of the component, which is crucial when modeling diffraction effects.

To compute the “angular modulation” term, we leverage concepts from the frequency-domain geometric theory of diffraction [16, 17]. However, in order to obtain a time-domain diffraction model that is suitable for our purposes, we first need to adapt the available results, effectively developing a novel time-domain diffraction model in the process.

The number of field components appearing in the surrogate is determined by the number of reflection and diffraction effects that are required to faithfully approximate the target wave, which ultimately depends on the geometry of the physical domain.

Among the advantages of the proposed approach, we mention the fact that each field component is separable into time-radial and angular components (in the “polar coordinates” sense). As we will see, we can leverage this to drastically reduce the computational cost and storage memory requirements of our approximation strategy.

The rest of the paper is structured as follows. In Section 1.1 we present the problem of interest. In Section 2 we introduce the main ingredients of our method, and we present our surrogate modeling algorithm. In Sections 3 and 4 we detail how we model reflection and diffraction effects, respectively. In Section 5 we present some numerical results to showcase our method. Both simple benchmarks (wedges) and more complicated tests (2D room model with scatterers) are considered. Some final considerations follow in Section 6.

1.1 Target problem

We are interested in the numerical approximation of the solution of the wave equation in complex domains. In this work, we consider 2-dimensional domains only. However, most of our discussion generalizes to 3D. We defer a discussion on this till Section 6.

We denote by $\Omega\subset\mathbb{R}^{2}$ the physical domain in which the wave equation is considered. We assume that $\Omega$ is either a closed polygon or a set-subtraction of polygons (to allow for multiply connected domains). We denote by $n_{e}$ and $n_{v}$ the number of edges and vertices of $\partial\Omega$, respectively. We study the propagation of waves in $\Omega$ over a given time interval of interest $[0,T]$. The model problem is the wave equation with constant (unit) wave speed:

$$\begin{cases}\partial_{tt}u(\bm{x},t)=\Delta u(\bm{x},t)+f(\bm{x},t)\quad&\text{for }(\bm{x},t)\in\Omega\times(0,T),\\ u(\bm{x},0)=u_{0}(\bm{x}),\partial_{t}u(\bm{x},0)=u_{1}(\bm{x})\quad&\text{for }\bm{x}\in\Omega,\\ \partial_{\nu}u(\bm{x},t)=0\quad&\text{for }(\bm{x},t)\in\partial\Omega\times(0,T],\end{cases}$$ (1)

with $\Delta$ the Laplacian operator, defined, in 2 dimensions, as $\Delta=\sum_{j=1}^{2}\partial_{x_{j}x_{j}}$. The homogeneous Neumann condition (i.e., the last equation above) models the whole boundary $\partial\Omega$ as sound-hard [7]. More generally, all or parts of $\partial\Omega$ may be modeled as sound-soft via a Dirichlet-type condition: $u(\bm{x},t)=0$. Depending on the target application, impedance boundary conditions may also be appropriate: $Z\partial_{\nu}u(\bm{x},t)+\partial_{t}u(\bm{x},t)=0$, with $Z\geq 0$.

We assume that the initial conditions $u_{0}$ and $u_{1}$, as well as the forcing term $f$, have radial symmetry around a given point. Without loss of generality, we will take such point to be the origin of $\mathbb{R}^{2}$:

$$u_{0}(\bm{x})=\eta_{0}(\left\|\bm{x}\right\|),\;u_{1}(\bm{x})=\eta_{1}(\left\|\bm{x}\right\|),\;f(\bm{x},t)=\eta_{2}(\left\|\bm{x}\right\|,t),$$ (2)

for all $(\bm{x},t)\in\Omega\times(0,T)$, with $\left\|\bm{x}\right\|^{2}=\sum_{j=1}^{2}x_{j}^{2}$. We further assume that the functions $\eta_{j}$ have compact support, namely, that there exist $R>0$ such that $\eta_{j}(\rho)=0$ for all $\rho>R$ and $j=0,1,2$. Moreover, to avoid incompatibilities with the boundary conditions, for simplicity we will only consider the situation where the supports of the functions $\eta_{j}$ are fully contained in $\Omega$.

2 Approximation framework

Before we can model boundary effects (reflection and diffraction), we need to understand how the solution $u$ would behave if no boundary were present. To this aim, we consider the wave equation in free space

$$\begin{cases}\partial_{tt}U(\bm{x},t)=\Delta U(\bm{x},t)+f(\bm{x},t)\quad&\text{for }(\bm{x},t)\in\mathbb{R}^{2}\times(0,\infty),\\ U(\bm{x},0)=u_{0}(\bm{x})\quad&\text{for }\bm{x}\in\mathbb{R}^{2},\\ \partial_{t}U(\bm{x},0)=u_{1}(\bm{x})\quad&\text{for }\bm{x}\in\mathbb{R}^{2},\end{cases}$$ (3)

which we have obtained from Eq. 1 by replacing $\Omega$ with the whole plane.

Due to radial symmetry (of the initial conditions and of the forcing term), we can recast the problem in polar coordinates. This allows us to define the free-space solution $\Psi$ in the radial coordinate, as the solution of

$$\begin{cases}\partial_{tt}\Psi(\rho,t)=\widehat{\Delta}\Psi(\rho,t)+\eta_{2}(\rho,t)\quad&\text{for }(\rho,t)\in(0,\infty)\times(0,\infty),\\ \Psi(\rho,0)=\eta_{0}(\rho)\quad&\text{for }\rho\in[0,\infty),\\ \partial_{t}\Psi(\rho,0)=\eta_{1}(\rho)\quad&\text{for }\rho\in[0,\infty),\end{cases}$$ (4)

where $\widehat{\Delta}$ is the Laplace operator in polar coordinates (under radial symmetry), i.e., $\widehat{\Delta}=\partial_{\rho\rho}+\frac{1}{\rho}\partial_{\rho}$, and $U(\bm{x},t)=\Psi(\left\|\bm{x}\right\|,t)$ for all $\bm{x}\in\mathbb{R}^{2}$. For closure, it may also be necessary to provide a condition at $\rho=0$, e.g., $\partial_{\rho}\Psi(0,\cdot)=0$ if the data $\eta_{j}$ are sufficiently smooth [18]. Note that, by the compact support of the initial conditions and of the forcing term, and by the finite (unit) speed of propagation of the wave equation, we have $\Psi(\rho,t)=0$ whenever $\rho>t+R$.

Remark 2.1.

Generally, the free-space solution $\Psi$ is not available analytically, except for very simple choices of initial conditions and forcing term. Accordingly, in most applications, the function $\Psi$ will need to be replaced with a suitable approximation. To this effect, one can discretize Eq. 4, e.g., with a finite element approximation (in space) and some timestepping scheme (in time). See Section 5 for more details on this.

Our goal is to approximate, for all $(\bm{x},t)\in\Omega\times[0,T]$, the solution $u(\bm{x},t)$ of the wave equation problem in Eq. 1 with a surrogate, defined as a sum of special functions $\widetilde{u}_{n}$ (field components):

$$u(\bm{x},t)\approx\widetilde{u}(\bm{x},t)=\sum_{n=1}^{N}\underbrace{\Psi(\left\|\bm{x}-\bm{\xi}_{n}\right\|+r_{n},t)\mathbbm{1}_{\Omega_{n}}(\bm{x})\zeta_{n}(\bm{x}-\bm{\xi}_{n})}_{\widetilde{u}_{n}(\bm{x},t)}.$$ (5)

Therein, $\Psi$ is the above-mentioned free-space radially symmetric solution of Eq. 4, and $\mathbbm{1}_{A}$ denotes the indicator function with support $A$, i.e.,

$$\mathbbm{1}_{A}(y)=\begin{cases}1\quad&\text{if }y\in A,\\ 0\quad&\text{if }y\notin A.\end{cases}$$ (6)

Moreover, in Eq. 5, we have introduced the following quantities:

• $N$ is the number of field components used in the approximation.

• $\bm{\xi}_{n}$ is the location of the point source from which $\widetilde{u}_{n}$ originates.

• $r_{n}\geq 0$ is a spatial delay, which is used for the synchronization of diffraction effects.

• $\Omega_{n}\subset\Omega$ is the spatial support of $\widetilde{u}_{n}$.

• $\zeta_{n}:\mathbb{R}^{2}\setminus\{\bm{0}\}\to\mathbb{R}$ is a weight function describing the angular modulation. We assume $\zeta_{n}$ to be a positive-homogeneous function, i.e., $\zeta_{n}(\bm{y})=\zeta_{n}(\bm{y}/\left\|\bm{y}\right\|)$ for all $\bm{y}\in\mathbb{R}^{2}\setminus\{\bm{0}\}$.

Remark 2.2.

We refer to $\zeta_{n}$ as “angular modulation” since, in 2D, the positive-homogeneity assumption is equivalent to requiring $\zeta_{n}(\bm{y})$ to depend only on the angle between the vector $\bm{y}$ and some reference direction, e.g., the positive $x_{1}$-axis.

Note that, due to the finite speed of propagation of the free-space solution $\Psi$, we have that a generic term $\widetilde{u}_{n}(\bm{x},t)$ is zero whenever $t<\left\|\bm{x}-\bm{\xi}_{n}\right\|+r_{n}-R$, i.e., for $t$ small enough, depending on $\bm{x}$.

The number of field components $N$ in Eq. 5 will be determined based on how many boundary effects (reflections and diffractions) need to be included in $\widetilde{u}$ in order to have a good approximation of the target wave $u$. We describe a strategy for automatically identifying a good $N$ in the next section.

2.1 Building the low-rank skeleton

Recalling that $u$ solves the wave equation in Eq. 1 in the domain $\Omega$, we use the first term in Eq. 5, namely, $\widetilde{u}_{1}$, to approximate the outgoing component of $u$, ignoring any effect due to the boundary $\partial\Omega$, except for shadows. Then, given such $\widetilde{u}_{1}$, we use the other terms $\widetilde{u}_{2},\ldots,\widetilde{u}_{N}$ to correct this first approximation. Each extra term models a single effect due to a certain portion of the boundary, specifically, an edge (reflection off that edge) or a vertex (diffraction about that vertex).

Going back to the first field component $\widetilde{u}_{1}$, let us define it, by providing its “ingredients” $\bm{\xi}_{1}$, $r_{1}$, $\Omega_{1}$, and $\zeta_{1}$, cf. Eq. 5. We set $\bm{\xi}_{1}=\bm{0}$, the center of the initial condition, as well as $r_{1}=0$, since no delay is necessary for this first term. Then, leveraging symmetry, we set $\zeta_{1}\equiv 1$, which corresponds to the (physical) assumption that the propagation of $\widetilde{u}_{1}$ is purely radial. Finally, we set $\Omega_{1}$ (the spatial support of the first field component around $\bm{0}$) as the set of points that can be reached from $\bm{0}$ via a straight line without going outside $\partial\Omega$, i.e.,

$$\Omega_{1}=\left\{\bm{x}\in\Omega\ :\ \tau\bm{x}\in\Omega\ \ \forall 0\leq\tau\leq 1\right\}.$$ (7)

In summary, the first term of $\widetilde{u}$ is

$$\widetilde{u}_{1}(\bm{x},t)=\Psi(\left\|\bm{x}\right\|,t)\mathbbm{1}_{\Omega_{1}}(\bm{x}).$$ (8)

Remark 2.3.

From a practical viewpoint, the implementation of spatial supports like $\Omega_{1}$ is not trivial. In our code, we verify that an arbitrary point $\bm{x}\in\Omega$ belongs to $\Omega_{1}$ by checking if there are no intersections between the line segment from $\bm{x}$ to $\bm{\xi}_{1}$ and any of the edges composing $\partial\Omega$. In some sense, this may be considered akin to “ray-tracing” [6]. A similar approach, with minor modifications, works for the other spatial supports $\Omega_{n}$ related to reflected or diffracted waves, cf. Eqs. 14 and 15 below.

Then we can move to the subsequent terms $\widetilde{u}_{n}$, $n\geq 2$. Their expressions depend on our choice of reflection and diffraction modeling, and will be provided in the upcoming sections. Instead, in the rest of the present section we focus on understanding how large $N$ should be, in order for $\widetilde{u}$ to provide a faithful approximation of $u$. Equivalently, we want to count the number of times the wave gets reflected or diffracted at the boundary $\partial\Omega$. This is done incrementally, starting from the initial value $N=1$ (no boundary effects) and then updating this guess as more and more boundary effects get “discovered”.

Figure 1: Computation of some timetable entries. The boundary $\partial\Omega$ has 11 sides, so that, e.g., $(\bm{a}_{1})_{14}$ is related to $\bm{y}_{3}$ and $(\bm{a}_{1})_{17}$ is related to $\bm{y}_{6}$. The shadowed area $\Omega\setminus\Omega_{1}$ is in darker grey.

To help us in this endeavor, we employ what we call a timetable, which, in this work, is simply a list of vectors, each with size $n_{e}+n_{v}$. (Note that similar ideas can can also be found in [19, 5, 6] and references therein.) The timetable is built incrementally starting from an empty list, appending one new vector every time a new term is added in the sum in Eq. 5, starting from $\widetilde{u}_{1}$. The entries of the $n$-th timetable vector are the waiting times before $\widetilde{u}_{n}$ comes in contact with an edge or a vertex of $\partial\Omega$. If it is impossible for $\widetilde{u}_{n}$ to “cast light” (along a straight path) onto a certain edge or vertex (e.g., $\gamma_{3}$ and $\bm{y}_{4}$ in Fig. 1), then the corresponding entry in the timetable is set to $\infty$. After this, it suffices to look for the smallest not-yet-explored entry of the timetable to identify what the next term of the approximation $\widetilde{u}$ should be. Once the entry in the timetable has been explored, its value is set to $\infty$.

We start by describing how the first vector $\bm{a}_{1}\in\mathbb{R}^{n_{e}+n_{v}}$ of the timetable (corresponding to $\widetilde{u}_{1}$) is computed, and how $\bm{a}_{1}$ allows us to identify the (geometric) features of $\widetilde{u}_{2}$. The vector $\bm{a}_{1}$ can be partitioned into edges-related part (the first $n_{e}$ entries) and vertices-related part (the last $n_{v}$ entries).

• Edge times. Given a generic edge $\gamma_{j}\subset\partial\Omega$ ($j=1,\ldots,n_{e}$) belonging to the domain boundary, we define the corresponding entry of $\bm{a}_{1}$ as

  $$\left(\bm{a}_{1}\right)_{j}=\begin{cases}r_{1}+\inf\left\{\left\|\bm{x}-\bm{\xi}_{1}\right\|:\bm{x}\in\gamma_{j}\cap\overline{\Omega}_{1}\right\}&\text{if the set is non-empty},\\ \infty&\text{otherwise}.\end{cases}$$ (9)

  Note that we have taken the shortest path from $\bm{\xi}_{1}$ to $\gamma_{j}$, and that we have denoted the closure of $\Omega_{1}$ as $\overline{\Omega}_{1}$.

• Vertex times. Given a generic vertex $\bm{y}_{j}\subset\partial\Omega$ ($j=1,\ldots,n_{v}$) of the domain boundary, we set

  $$\left(\bm{a}_{1}\right)_{n_{e}+j}=\begin{cases}r_{1}+\left\|\bm{y}_{j}-\bm{\xi}_{1}\right\|\quad&\text{if }\bm{y}_{j}\in\overline{\Omega}_{1},\\ \infty\quad&\text{otherwise}.\end{cases}$$ (10)

Note that we have included the delay $r_{1}$ (which is actually zero here) as a way to streamline Eqs. 9 and 10 for the upcoming section. See Fig. 1 for a diagram showcasing these formulas.

The smallest entry of $\bm{a}_{1}$ is the time at which the first “boundary event” (reflection or diffraction) can happen. The index of the smallest entry tells us whether the event is a reflection (index $1\leq j\leq n_{e}$) or a diffraction (index $n_{e}+1\leq j\leq n_{e}+n_{v}$), and also what edge/vertex causes the event. From here, we use the models described in Sections 3 and 4 to build $\widetilde{u}_{2}$, by computing $\bm{\xi}_{2}$, $r_{2}$, $\Omega_{2}$, and $\zeta_{2}$.

Then, the second timetable vector $\bm{a}_{2}$ can be computed by replacing all subscripts “1” by “2” in Eqs. 9 and 10. This is followed by the construction of $\widetilde{u}_{3}$, and so on. The process continues until all not-yet-explored entries of the timetable are larger than $T+R$, where, as mentioned above, $R$ is the half-width of the support of forcing term and initial conditions. Indeed, starting from this time instant, the would-be next terms of $\widetilde{u}$ do not affect the approximation anymore, since, due to the finite speed of wave propagation, they only act (on $\Omega$) after the end of the time horizon, i.e., for $t>T$. The total number of field components $N$ is simply the number of vectors in the timetable.

We summarize the overall procedure for the construction of the terms $\widetilde{u}_{n}$ in Algorithm 1. For ease of presentation, once an entry of the timetable has been explored, it is set to $\infty$ as a way for the algorithm to ignore it from that point forward.

Algorithm 1 Step-by-step construction of the surrogate model

  Set $N\leftarrow 1$, find $\Omega_{1}$ as in Eq. 7, and define $\widetilde{u}_{1}$ as in Eq. 8

  Define $\bm{a}_{1}\in\mathbb{R}^{n_{e}+n_{v}}$ using Eqs. 9 and 10

  Set $i\leftarrow 1$ and $j\leftarrow\operatorname*{arg\,min}_{j=1,\ldots,n_{e}+n_{v}}\left(\bm{a}_{1}\right)_{j}$

  while $\left(\bm{a}_{i}\right)_{j}\leq T+R$ do

     Set $\left(\bm{a}_{i}\right)_{j}\leftarrow\infty$ and $N\leftarrow N+1$

     if $j\leq n_{e}$ then

        Find $\bm{\xi}_{N}$, $r_{N}$, $\Omega_{N}$, and $\zeta_{N}$ as in Section 3 $\leftarrow$ Reflection from edge $j$

     else

        Find vertex index $j^{\prime}\leftarrow j-n_{e}$

        Find $\bm{\xi}_{N}$, $r_{N}$, $\Omega_{N}$, and $\zeta_{N}$ as in Section 4 $\leftarrow$ Diffraction from vertex $j^{\prime}$

     end if

     Define $\widetilde{u}_{N}$ from $\bm{\xi}_{N}$, $r_{N}$, $\Omega_{N}$, and $\zeta_{N}$, as in Eq. 5

     Define $\bm{a}_{N}\in\mathbb{R}^{n_{e}+n_{v}}$ using Eqs. 9 and 10, with “$N$” replacing “$1$” in subscripts

     Set $(i,j)\leftarrow\operatorname*{arg\,min}_{i=1,\ldots,N,j=1,\ldots,n_{e}+n_{v}}\left(\bm{a}_{i}\right)_{j}$

  end while

In trapping domains, see, e.g., Section 5.2, the number of terms $N$ might be rather large due to waves repeatedly “bouncing back and forth” between two or more edges/vertices. A large $N$, although necessary for a good approximation of all wavefronts, is undesirable since it increases the computational cost of both the construction of the surrogate $\widetilde{u}$ and its evaluation.

In [19], it is suggested to set a fixed upper bound (e.g., 1 or 2) on the maximum number of successive diffraction events that are taken into account. Here, we propose and employ an alternative method that allows for a finer control on whether a certain field component is included or not, based on its magnitude. Specifically, we remove all terms $\widetilde{u}_{n}$ that are smaller (in magnitude) than a certain tolerance tol, uniformly over $\bm{x}$ and $t$. This can be done as a post-processing step (thus speeding up the evaluation of $\widetilde{u}$ but not its construction) or even while building the surrogate itself. This can be achieved with a simple modification of Algorithm 1, by introducing a test on the magnitude of each soon-to-be-added wave contribution $\widetilde{u}_{n}$, discarding terms that are too small.

3 Modeling reflection

We now present the strategy for modeling reflection due to an edge $\gamma$ of the domain boundary $\partial\Omega$. We rely on the well-known “geometrical optics” model, which describes wave propagation in terms of rays, not accounting for any diffraction [20]. We assume that we are adding a new term $\widetilde{u}_{n}$ to the surrogate model in Eq. 5, due to a reflection phenomenon caused by the field component $\widetilde{u}_{i}$, $i<n$. Specifically, we assume that a wave originating at $\bm{\xi}_{i}$ hits the edge $\gamma\subset\partial\Omega$, which, in particular, requires $\gamma\cap\overline{\Omega_{i}}\neq\emptyset$. We need to prescribe several ingredients.

Spatial correction $r_{n}$

We just transfer $r_{n}$ over from the incoming wave: $r_{n}=r_{i}$. Indeed, as we will see in Section 4, we require the term $r_{n}$ only when modeling diffraction.

Source point $\bm{\xi}_{n}$

We use the method of images, which gives the position of $\bm{\xi}_{n}$ as the reflection of $\bm{\xi}_{i}$ with respect to the edge $\gamma$:

$$\bm{\xi}_{n}=2\operatorname*{arg\,min}_{\bm{z}\in\widetilde{\gamma}}\left\|\bm{z}-\bm{\xi}_{i}\right\|-\bm{\xi}_{i},$$ (11)

where $\widetilde{\gamma}\subset\mathbb{R}^{2}$ is the straight line on which edge $\gamma$ lies. See Fig. 2 (left).

Figure 2: Graphical representation of a reflection off edge $\gamma$. On the left, the law of reflection prescribes $\theta_{r}=\theta_{i}$. We show the straight line $\widetilde{\gamma}$ supporting $\gamma$ with a dotted stroke. For a given observation point $\bm{x}$, $\bm{y}(\bm{x})$ denotes the point of incidence of the reflected field component. On the right, computation of the spatial support $\Omega_{n}$ (light-grey area) and its complementary shadow zone $\Omega\setminus\Omega_{n}$ (dark-grey area) for the reflected field component, in the presence of a rectangular obstacle. The dashed portion of edge $\gamma$ denotes the shadow $\gamma\setminus\gamma^{(i)}$. The shadow zone consists of two connected components.

Weight function $\zeta_{n}$

Let $\bm{x}-\bm{\xi}_{n}$ be a generic point where we wish to evaluate the weight function $\zeta_{n}$. We define the incidence point $\bm{y}(\bm{x})$ as the intersection (if any) between edge $\gamma$ and the segment from $\bm{\xi}_{n}$ to $\bm{x}$. See Fig. 2 (left). According to the method of images, the amplitude of the reflected wave is equal (up to sign) to the amplitude of the incoming wave:

$$\zeta_{n}(\bm{x}-\bm{\xi}_{n})=\pm\zeta_{i}(\bm{y}(\bm{x})-\bm{\xi}_{i}).$$ (12)

Above and throughout this section, the sign $\pm$ depends on the type of boundary conditions on the edge $\gamma$ (“$+$” for Neumann, “$-$” for Dirichlet).

Now, recall that we are assuming all weight functions to be positive-homogeneous. According to Remark 2.2, this means that $\zeta_{i}(\bm{x}-\bm{\xi}_{i})$ is only a function of the angle $\phi_{i}(\bm{x})$ between $\bm{x}-\bm{\xi}_{i}$ and the positive $x_{1}$-axis. See Fig. 2 (left) for a graphical depiction. Specifically, with an abuse of notation, let $\zeta_{i}(\bm{x}-\bm{\xi}_{i})=\zeta_{i}(\phi_{i}(\bm{x}))$ and $\zeta_{n}(\bm{x}-\bm{\xi}_{n})=\zeta_{n}(\phi_{n}(\bm{x}))$, where the “new” angle-dependent functions $\zeta_{i}$ and $\zeta_{n}$ are $2\pi$-periodic. By Eq. 12, we deduce the property

$$\zeta_{n}(\phi_{n}(\bm{x}))=\pm\zeta_{i}(\phi_{i}(\bm{y}(\bm{x})))= \pm\zeta_{i}(2\beta_{\gamma}-\phi_{n}(\bm{y}(\bm{x})))$$

$$= \pm\zeta_{i}(2\beta_{\gamma}-\phi_{n}(\bm{x})),$$ (13)

where $\beta_{\gamma}$ is the angle between edge $\gamma$ and the positive $x_{1}$-axis. This uniquely identifies $\zeta_{n}$, given $\zeta_{i}$ and $\beta_{\gamma}$.

Spatial support $\Omega_{n}$

We first identify what portion of $\gamma$ is actually “lit” by $\widetilde{u}_{i}$: $\gamma^{(i)}=\gamma\cap\overline{\Omega_{i}}$. Note that we may have $\gamma\neq\gamma^{(i)}$, for instance when obstacles are present between $\bm{\xi}_{i}$ and $\gamma$. See Fig. 2 (right) for an illustration. Then, roughly speaking, we define the new support $\Omega_{n}$ as the union of all line segments from $\bm{\xi}_{n}$ that pass through $\gamma^{(i)}$. To be more precise, given $\bm{x}\in\Omega$, let $\bm{y}(\bm{x})$ be the intersection (if any) between $\gamma$ and the line segment from $\bm{\xi}_{n}$ to $\bm{x}$. Also, if $\bm{y}(\bm{x})$ exists, we define $\tau_{0}(\bm{x})=\left\|\bm{y}(\bm{x})-\bm{\xi}_{n}\right\|/\left\|\bm{x}-\bm{\xi}_{n}\right\|\in(0,1)$, which satisfies $\bm{y}(\bm{x})=\bm{\xi}_{n}+\tau_{0}(\bm{x})(\bm{x}-\bm{\xi}_{n})$. The new support is defined as

$$\Omega_{n}=\left\{\bm{x}\in\Omega\ :\ \bm{y}(\bm{x})\in\gamma^{(i)}\text{ and }\bm{\xi}_{n}+\tau(\bm{x}-\bm{\xi}_{n})\in\Omega\ \ \forall\tau_{0}(\bm{x})<\tau\leq 1\right\}.$$ (14)

Figure 3: Example of reflection off an edge in the presence of an obstacle, from Fig. 2. Neumann conditions are imposed on all edges. Source wave (left), reflected wave (middle), and superimposition of the two (right). Note how the obstacle creates a shadow zone for source and reflected waves. For simplicity, in this plot we are not showing any reflection or diffraction effects due to the rectangular obstacle, since they would be modeled at different stages of the algorithm.

Figure 3 represents a possible output of the numerical algorithm. In this case, we simulate only the reflections, thus discarding, for the time being, any effect due to diffraction. It is clear that, by modeling reflection effects only, we may obtain a discontinuous approximation of the solution of our target problem, with discontinuities arising at the boundaries of the spatial supports identified so far. As we will see in the next section, introducing diffraction in our approximation will allow us to obtain a continuous approximation $\widetilde{u}$.

4 Modeling diffraction

Here, we describe a strategy for modeling waves diffracted by a vertex of the domain boundary $\partial\Omega$. This is required in building a new field component $\widetilde{u}_{n}$ whenever the smallest unexplored entry of the timetable is related to a vertex, i.e., $j>n_{e}$ in Algorithm 1, cf. Section 2.

In our modeling, we take inspiration from the uniform theory of diffraction (UTD) [17, 20, 21]. UTD provides an effective description of diffraction in the frequency domain, which we need to specialize for our time-domain modeling. This task is difficult for several reasons. On the one hand, standard approaches for frequency-to-time-domain conversion involve the inversion of the Fourier transform, which is rather costly in a numerical setting, where integration/convolution must be replaced by quadrature. (This issue actually affects also time-domain diffraction modeling based on impulse responses, see, e.g., [22, 23].) On the other hand, UTD is an asymptotic method, whose accuracy relies on assuming the frequency to be sufficiently large. As such, an inverse Fourier transform is not guaranteed to provide sound results, whenever the signal bandwidth includes low frequencies.

Both these reasons behoove us to develop a novel time-domain version of UTD. Specifically, we wish our diffraction modeling to fit in the framework of Eq. 5, so that:

• Diffraction is modeled as a wave outgoing from a point source at some $\bm{\xi}_{n}$. This is actually in line with standard approaches, e.g., [24, 6], where diffraction is modeled through virtual sources on the vertex (in 2D) or edge (in 3D) that causes the diffraction event. This motivates the choice of the center $\bm{\xi}_{n}:=\bm{y}_{j^{\prime}}=\bm{y}_{j-n_{e}}$, the diffraction vertex (we are employing the notation of Algorithm 1).

• Diffraction has a spatial support $\Omega_{n}$, which (again in accordance to diffraction theory) we define as the set of all points that are visible (along straight-line paths) from $\bm{\xi}_{n}$, i.e.,

  $$\Omega_{n}=\left\{\bm{x}\in\Omega\ :\ \bm{\xi}_{n}+\tau(\bm{x}-\bm{\xi}_{n})\in\Omega\ \ \forall 0<\tau\leq 1\right\}.$$ (15)

• The space-time dependence of the diffraction amplitude should be separable into an angular component $\zeta_{n}$, and a radial-temporal component, which, in fact, is assumed to coincide with the free-space wave propagation profile $\Psi$.

The last step is the most critical one, as it involves a simplifying modeling choice, specifically made to fit the approximate diffraction wave within our framework. Still, there are strong theoretical foundations for this, as we proceed to explain.

4.1 Modeling the angular modulation

From a phenomenological point of view, the angular modulation $\zeta_{n}$ can be related to the diffraction coefficient appearing in geometrical diffraction theory [16, 20]. Indeed, the diffraction coefficient $D$ has exactly the desired role of scaling factor for the diffraction amplitude, which depends on the angle that the target measurement point $\bm{x}$ forms with the edges adjacent the diffraction vertex. As such, we set $\zeta_{n}$ equal to the UTD diffraction coefficient $D$.

We define $D$ in full detail in Appendix A. Here, for our discussion, we only need the following basic facts: in 2-dimensional UTD, the diffraction coefficient $D$ depends on the above-mentioned angle between $\bm{x}$ and $\bm{\xi}_{n}$, but also on (i) the domain geometry locally around the diffraction point, (ii) the location of the source that causes the diffraction, (iii) the frequency of the incident signal, and (iv) the distance $\left\|\bm{x}-\bm{\xi}_{n}\right\|$.

The last two above-mentioned items are problematic in our framework: since we work in the time domain, we do not have a single incident frequency $k$; moreover, due to our separability assumption, we require the angular modulation to be independent of $\left\|\bm{x}-\bm{\xi}_{n}\right\|$. However, frequency and radial distance always appear together, as a product, in the UTD diffraction coefficient. Specifically, such product $\mu:=k\left\|\bm{x}-\bm{\xi}_{n}\right\|$ can be interpreted as the distance between $\bm{x}$ and the diffraction point $\bm{\xi}_{n}$, measured in wavelength units.

In practice, the value of $\mu$ determines how quickly the magnitude of $D$ decays to $0$, as one moves away from a so-called “shadow boundary”, i.e., the boundary of the spatial support of either the incident or the reflected wave(s). Equivalently, $\mu$ determines the width of the transition region, which encompasses any shadow boundaries. See Appendices A and 15. Also, note that UTD requires $\mu$ to be large enough ($\geq 1$), and that $\mu\to\infty$ yields the geometrical-optics setting, i.e., a diffraction-free model, with $\zeta_{n}\equiv 0$.

Since the dependence of $D$ on the two troublesome terms happens only through $\mu$, we propose the following strategy for defining a $k$- and $\left\|\bm{x}-\bm{\xi}_{n}\right\|$-independent angular modulation $\zeta_{n}$: we set $\zeta_{n}$ as the diffraction coefficient $D$ obtained for some value of $\mu=\overline{\mu}$ that is fixed a priori. Such value $\overline{\mu}$ should be sufficiently large not to cause issues in UTD, but also sufficiently small to avoid the geometrical-optics pitfall.

Obviously, the specific choice of $\overline{\mu}$ should be based on the ultimate approximation target: minimizing the discrepancy between $u$ and $\widetilde{u}$. However, it turns out that $D$ depends relatively mildly on the value of $\mu$. See Appendix A for a theoretical justification, and Section 5.1 for some numerical tests. In our experiments, we settle for the value $\overline{\mu}=10$. A more careful quantitative investigation of the role of $\mu$ on the diffraction approximation accuracy is envisioned as a future research direction.

Instead of using frequency-domain UTD, the time-domain diffraction coefficient could be computed by convolution of the incoming wave with a “diffraction impulse response”, as done in [23, 25]. Such a convolution would need to be approximated (by quadrature) over and over, at each diffraction event. Notably, since the diffraction impulse response is not a Dirac delta, a convolution with it naturally introduces a radial modulation too. This entails that the diffraction wave does not behave like an angular modulation of the free-space wave $\Psi$, as assumed in our model.

Our modeling of diffraction through a diffraction coefficient that does not depend on time nor radius, although biased in the above-described way, is very efficient, as it avoids the expensive computation of convolutions with the diffraction impulse response. As we will show in our numerical experiments, this increased efficiency comes at the cost of a relatively small modeling error. As such, our approach remains justified.

4.2 Estimating the modeling error

We can assess the quality of our diffraction model a priori, by comparing its resulting approximation of the diffraction wave ($\widetilde{u}_{n}$) with that obtained via the above-mentioned convolution approach [23], which we denote by $\mathring{u}_{n}$ here.

As before, assume that field component $\widetilde{u}_{i}$ impinges on vertex $\bm{\xi}_{n}$, causing a diffraction event emanating from there. Given an arbitrary space-time position $(\bm{x},t)$, the local discrepancy between the diffraction waves obtained with the two methods can be expressed by inverse Fourier transform, as

$$\widetilde{u}_{n}(\bm{x},t)-\mathring{u}_{n}(\bm{x},t)= \widetilde{D}_{\overline{\mu}}\int_{\mathbb{R}}\mathcal{F}[\widetilde{u}_{i}](\bm{\xi}_{n},k)e^{ikt}\textrm{d}k-\int_{\mathbb{R}}D_{k}\mathcal{F}[\widetilde{u}_{i}](\bm{\xi}_{n},k)e^{ikt}\textrm{d}k$$

$$= \int_{\mathbb{R}}\left(\widetilde{D}_{\overline{\mu}}-D_{k}\right)\mathcal{F}[\widetilde{u}_{i}](\bm{\xi}_{n},k)e^{ikt}\textrm{d}k.$$

Above, $\mathcal{F}$ denotes the Fourier transform operator, whereas $D_{k}=D_{k}(\bm{x}-\bm{\xi}_{n})$ is the UTD diffraction coefficient, which depends on both the position $\bm{x}$ and the wavenumber $k$, cf. Appendix A. On the other hand, $\widetilde{D}_{\overline{\mu}}$ is the diffraction coefficient used in our modeling, which is $k$-independent and depends on $\bm{x}$ only through the direction of $\bm{x}-\bm{\xi}_{n}$. (Note that, in geometrical optics, one would set $\widetilde{D}_{\overline{\mu}}=0$.)

It follows that, to achieve a small error, one should pick a value of $\widetilde{D}_{\overline{\mu}}$ close to that of $D_{k}$, specifically at wavenumbers $k$ that lie closest to the Fourier spectrum of $\widetilde{u}_{i}(\bm{\xi}_{n},\cdot)$. In this way, we can expect that our approach, despite being based on a heuristic choice of $\overline{\mu}$, independently of $k$, should do significantly better than geometrical optics.

This being said, our method has limitations, which are better seen by looking at the following global-in-wavenumber bound on the error magnitude:

$$\left|\widetilde{u}_{n}(\bm{x},t)-\mathring{u}_{n}(\bm{x},t)\right|\leq\sup_{k}\left|\widetilde{D}_{\overline{\mu}}-D_{k}\right|\left\|\mathcal{F}[\widetilde{u}_{i}](\bm{\xi}_{n},\cdot)\right\|_{L^{1}(\mathbb{R})}.$$ (16)

Our method may lead to gross errors if $D_{k}$ has large variations in $k$. In this respect, it is worth nothing that, thanks to the availability of explicit formulas for $D_{k}$, cf. Appendix A, it is possible to obtain an upper-bound for $|D_{k}-\widetilde{D}_{\overline{\mu}}|$ by direct inspection. Still, bound Eq. 16 is pessimistic, as we will show in our numerical experiments.

As a practical alternative, we propose a heuristic a posteriori error indicator, with the objective of obtaining better (but still approximate) information on the magnitude of the diffraction modeling error. To this aim, we apply the following multi-fidelity argument. Let $\widetilde{u}$ be the approximation obtained through our method, including also diffraction terms, modeled as in Section 4.1. On the other hand, let $\widetilde{u}_{\text{GO}}$ be the approximation obtained with our approach, disregarding diffraction completely (“GO” stands for “geometrical optics”). As indicator of the approximation error $\left|\widetilde{u}-u\right|$, we simply use the cross-fidelity discrepancy $|\widetilde{u}_{\text{GO}}-\widetilde{u}|$. The justification behind this follows.

First, we note that, since reflection effects are modeled exactly in our approach, cf. Section 3, $u-\widetilde{u}_{\text{GO}}$ coincides with the superposition of all diffraction-caused effects. Second, we make the heuristic assumption that diffraction effects in $\widetilde{u}$ are approximated with at most a $50\%$ relative error. (Of course, we have no way of verifying this in practice, as we would need access to the exact solution.) Accordingly, we have that $2|\widetilde{u}-u|\leq|\widetilde{u}_{\text{GO}}-u|$ and, by the triangular inequality,

$$\left|\widetilde{u}-u\right|\leq\left|\widetilde{u}_{\text{GO}}-u\right|-\left|\widetilde{u}-u\right|\leq\left|\widetilde{u}_{\text{GO}}-\widetilde{u}\right|.$$

We assess the reliability of this estimate in our numerical experiments in the next section.

5 Numerical results

In our experiments, we require a “reference” solution of Eq. 1 to validate our results. To this effect, we use the solution $u_{\textup{FE}}$ obtained by discretizing Eq. 1 with:

• the P1-finite element (FE) method with mass-lumping, over a regular triangulation (mesh) of the physical domain $\Omega$;

• explicit leapfrog timestepping with a uniform time step that satisfies the CFL condition on the chosen mesh.

See [7, 8] for more details on this discretization strategy.

If the domain $\Omega$ is unbounded, we first need to truncate it in such a way that reflections from the non-physical truncation boundary do not affect the solution in the region of interest for $t<T$. Recalling that the problem data are supported in a ball of radius $R$ and center $\bm{0}$, this can be done, e.g., by truncating $\Omega$ at the sphere with radius $R+T$ (we recall that we are assuming a unit wave speed) and center $\bm{0}$. In our tests, we rely on FEniCS [26] to carry out the FE discretization on 2-dimensional domains $\Omega$. Note that, instead, unbounded geometries are allowed in our proposed approach, making domain truncations unnecessary.

All our tests are performed in Python 3.8 on a machine with an 8-core 4.70 GHz Intel^® processor. For reproducibility, our code is made available at https://github.com/pradovera/ray-wave-2d.

5.1 Some simple wedges

As a way to assess our proposed method in simple settings, we consider four different “wedge” domains. We define $\Omega$ to be one of the portions of the plane $\mathbb{R}^{2}$ delimited by two straight lines intersecting at a point, with $\alpha$ being the outer angle at such point. The specific choices of wedge angles $\alpha$ are reported in Table 1 for the four cases.

Figure 4: Initial conditions for the wedge examples, indexed #1 through #4 from left to right. The (dashed) distance between the center of the Gaussian and the boundary vertex is 4 units in all cases.

Example	exterior	incidence	number $N$ of	$\sup_{t\in[0,T]}\textup{err}(t)$, cf. Eq. 17
index	angle $\alpha$	angle $\theta$	components $\widetilde{u}_{n}$
#1	$1.5\pi$	$0.313\pi$	5	$0.146\%$
#2	$1.62\pi$	$0.192\pi$	8	$0.828\%$
#3	$0.879\pi$	$0.434\pi$	4	$3.14\%$
#4	$0.879\pi$	$1.04\pi$	3	$4.68\%$

Table 1: Setup for the four wedge examples. The angle $\theta$ is as in Fig. 4. The last two columns refer to the results of our algorithm.

We set up a wave-propagation problem like Eq. 1, with $u_{0}$ an isotropic Gaussian with standard deviation $0.2$. The center of $u_{0}$ is at a point located at a 4-unit distance from the wedge vertex, in the direction determined by the “incidence angle” $\theta$. See Fig. 4 for a representation of the initial conditions in the four cases. We set $u_{1}=f=0$, we enforce Neumann boundary conditions on the whole $\partial\Omega$, and we seek the solution at the final time $T=5$, i.e., $1$ time unit after the wave crest has reached the wedge vertex.

Figure 5: Free-space solution $\Psi$. The dashed line denotes the upper bound of the “causality cone” of $\Psi$, i.e., $\rho=t+R$, with $R=1$.

Figure 6: Results for the four wedge examples. Each row pertains to a different example. In each row, from left to right: surrogate solution, FE solution, and error. The color scales for the first two columns are the same. All results are shown at the final time $t=T$.

To this aim, we employ our proposed approach, see Section 2. First, we compute an approximation of the free-space solution $\Psi$, which solves Eq. 4, by employing the P1-FE method with explicit leapfrog timestepping. Note that, since Eq. 4 is cast in polar coordinates, we only need to discretize a 1D interval with P1-FE. Since the initial condition $u_{0}$ is supported within the unit disk, we have $R=1$, and it suffices to approximate $\Psi(\rho,t)$ for $(\rho,t)\in[0,T+R]\times[0,T]$. Since this space-time domain is only 2-dimensional, we can afford a very fine discretization. In our experiments, we employ a $1001\times 2001$ uniform Cartesian space-time grid, i.e., the mesh size is $\delta x=\frac{T+R}{1000}$ and the time step is $\delta t=\frac{T}{2000}$. This satisfies the CFL condition. We show the resulting $\Psi$ (which, in fact, we should denote by $\Psi_{\textup{FE}}$) in Fig. 5. Note that, since we are solving the wave equation in 2D, the magnitude of $\Psi$ decays like $\mathcal{O}(t^{-1/2})$ as $t\to\infty$, which corresponds to a spatial decay of $\mathcal{O}(\rho^{-1/2})$. In 3D, the decay would be quicker, namely, $\mathcal{O}(t^{-1})$ in time and $\mathcal{O}(\rho^{-1})$ in space.

After this preliminary step, we use the timetable-based strategy from Section 2 to identify reflection and scattering effects, which are then added up to give the final approximation $\widetilde{u}$. We show the resulting $\widetilde{u}(\cdot,T)$ in Fig. 6. In this figure, we also display a reference solution $u_{\textup{FE}}(\cdot,T)$, which we obtain by direct discretization of Eq. 1 by P1-FE and leapfrog timestepping, as described at the beginning of Section 5.

In all four examples, we see that $\widetilde{u}$ and $u_{\textup{FE}}$ seem qualitatively close. Notably, we can observe a good representation of the most prominent wavefronts, which are due to propagation of either the main “free-space” wave or to its reflections. Indeed, those wave contributions are reconstructed exactly: the only errors are the ones due to FE approximation and timestepping, which affect both $u_{\textup{FE}}$ and $\widetilde{u}$ (the latter through the approximation of $\Psi$). Instead, some differences are present when comparing diffraction effects, which arise as circular waves about the wedge vertex. We can quantitatively observe this in the last column of both Tables 1 and 6.

In example #1, we observe a very small error, which, in fact, is simply the (FE and timestepping) discretization error. This is related to the fact that the wedge index $\nu=\pi/(2\pi-\alpha)=2$ is an integer, thus making diffraction unnecessary in approximating the wave $u$ [17].

In the other examples, diffraction effects are necessary to correctly identify $u$. While a good qualitative behavior can be observed in Fig. 6, we can see in Table 1 that a modest error is present. Specifically, as error measure, we use the supremum over $t\in[0,T]$ of the relative $L^{2}(\Omega)$-approximation error, defined as

$$\textup{err}(t):=\left(\int_{\Omega}\left(\widetilde{u}(\bm{x},t)-u_{\textup{FE}}(\bm{x},t)\right)^{2}\textup{d}\bm{x}\leavevmode\nobreak\ \bigg{/}\int_{\Omega}u_{\textup{FE}}(\bm{x},t)^{2}\textup{d}\bm{x}\right)^{1/2}.$$ (17)

We see the largest error in example #4, where the relative $L^{2}(\Omega)$-approximation error amounts to about $5\%$. This example corresponds to a case of “almost grazing” incidence, with the source point being located rather close to one of the wedge’s edges. In turn, this leads to two shadow boundaries located very close to each other, with overlapping transition regions, cf. Appendix A. We interpret this result as evidence of the fact that our diffraction model from Section 4 is less accurate with grazing than with non-grazing incidence. Despite this, the error remains small in all the above tests, and our diffraction model may be considered satisfactory overall.

5.1.1 Building a cavity out of wedges

As a slightly more complicated example, we now combine the four wedges from the previous section to obtain the open cavity represented in Fig. 7. In this case, more reflection and diffraction effects arise, due to the trapping nature of the domain. Our initial conditions and forcing term are the same as before, but now all edges are sound-soft. Accordingly, we model them using Dirichlet boundary conditions. The time horizon is $T=9$.

Figure 7: Graph view of the field components for the cavity example. Markers are used to denote source points $\bm{\xi}_{n}$. Dotted line segments denote causality effects: a black star is linked to a blue square (resp., red circle) if the field component emanating from the former causes a reflection (resp., diffraction) emanating from the latter. To avoid clutter, the field components have been separated into three plots according to causality: in the left plot we show all $7$ events (reflections or diffractions) directly caused by $\widetilde{u}_{1}$, in the middle plot we show all $20$ events directly caused by the $7$ field components in the left plot, etc.

Using our strategy from Section 2, we build the approximation $\widetilde{u}$, which contains $45$ wave terms ($1$ source wave, $32$ reflected waves, and $12$ diffraction waves). As visual reference, we display the source points of all field components $\widetilde{u}_{n}$ in Fig. 7. Each dotted line segment is a “causality link”, showing how one field component causes the next, by either reflection or diffraction.

Figure 8: Results for the cavity domain. Each row corresponds to a different time instant $t\in\{3,6,9\}$, from top to bottom. In each row, from left to right: surrogate solution $\widetilde{u}(\cdot,t)$, FE solution $u_{\textup{FE}}(\cdot,t)$, and error $\widetilde{u}(\cdot,t)-u_{\textup{FE}}(\cdot,t)$. The color scales for the first two columns are the same.

In Fig. 8, we compare the approximation $\widetilde{u}$ with the reference solution $u_{\textup{FE}}$, obtained as described at the beginning of Section 5, on a spatial mesh with around $10^{6}$ degrees of freedom. Once more, we see a good qualitative agreement between $\widetilde{u}$ and $u_{\textup{FE}}$, with the most important features of $u$ being identified well. We see “error waves” of small amplitude propagating from the 3 vertices of $\Omega$ that generate diffraction effects. These correspond to errors in diffraction modeling. More quantitatively, the relative $L^{2}(\Omega)$-approximation errors are roughly $0.05\%$ (at $t=3$), $7.21\%$ (at $t=6$), and $12\%$ (at $t=9$).

The approximation $\widetilde{u}$ is built in approximately $500$ms, out of which $420$ms are used to compute the free-space solution $\Psi$. This is a surprisingly small time, compared to the computation of the reference $u_{\textup{FE}}$, which takes about a minute, with the 2D meshing alone taking around $20$s. Note that we are including also the meshing time in order to have a fair comparison of the costs incurred by the two considered methods, with the target being obtaining a reasonably accurate approximate solution of the wave equation.

As an additional experiment to validate our diffraction modeling, we study how the approximation error depends on the choice of $\overline{\mu}$, which affects the diffraction coefficient by being (roughly) inversely proportional to the width of the transition regions around shadow boundaries, cf. Section 4. We compute the $L^{2}(\Omega)$-approximation error at times $t=6$ and $t=9$, for values of $\overline{\mu}$ between $1$ and $100$. Note that every value of $\overline{\mu}$ involves the computation of a new diffraction coefficient $D$ for each diffraction event, but otherwise does not require retraining the surrogate model $\widetilde{u}$: only the angular weights $\zeta_{n}$ vary, while all other terms in Eq. 5 remain the same.

Figure 9: Relative $L^{2}(\Omega)$-approximation error at $t=6$ and $t=9$ for different values of $\overline{\mu}$.

In Fig. 9, we show how the errors have a minimum around $\overline{\mu}=7$ for $t=6$, and around $\overline{\mu}=11$ for $t=9$. This shift in the minimum point seems to suggest that larger values of $\overline{\mu}$ should be used for approximations at longer times, an observation that may be especially relevant for simulations over long time horizons. That being said, the error appears to vary only mildly around the locations of such minima: for instance, at $t=6$, the relative error obtained for $\overline{\mu}=10$ is $7.21\%$, only slightly more than the minimum error, $7.06\%$. (Note also the narrowness of the vertical scale of the plot.) This provides an empirical justification for choosing a single value of $\overline{\mu}$ (around $10$) for all diffraction events, and for all times $t$.

5.2 A tall room

We now move to a simplified sound propagation problem in a room. For simplicity, we consider a 2-dimensional problem, thus assuming an infinitely tall room, and modeling line sources in the $z$-direction (e.g., an array of loudspeakers) as point sources.

The complicated domain $\Omega\subset\mathbb{R}^{2}$ is depicted in Fig. 10. It is composed of two communicating “rooms” with sound-hard walls, as well as of a third large room (above), which is modeled as infinitely large. In the main room, three sound-soft triangular obstacles are also present.

Setting once more $u_{1}=f=0$, we are interested in modeling the propagation of an initial condition $u_{0}$ modeled as a Ricker wavelet centered at $\bm{0}$, see Fig. 10 (top left), over the time horizon $t\in[0,T]$, with $T=20$. To this aim, we employ our proposed method from Section 2.

As in the previous example, we start by computing an approximation of the free-space solution $\Psi=\Psi(\rho,t)$ for $(\rho,t)\in[0,T+R]\times[0,T]$, see Eq. 4, with $R$ being the radius of the support of the initial condition $u_{0}$. Again, we use P1-FE with leapfrog timestepping for this.

Since many reflective surfaces face each other, the domain $\Omega$ is trapping. Accordingly, we expect the number $N$ of waves in the approximation $\widetilde{u}$ to be rather large. In the interest of reducing the number of such terms, we employ the on-the-fly tolerance-based strategy described in Section 2.1, removing all wave terms $\widetilde{u}_{n}$ whose magnitude is smaller than $\textup{tol}=2.5\cdot 10^{-2}$. After this, $N=798$ terms are left. Although this value of $N$ may seem large, the evaluation of the corresponding surrogate $\widetilde{u}$ is rather quick, due to the explicit nature of each wave contribution (and to the fact that their supports are smaller than the whole $\Omega$).

Figure 10: 2-dimensional domain $\Omega$ modeling a room. Top left plot: initial condition $\widetilde{u}(\cdot,0)=u(\cdot,0)=u_{0}$, a Ricker wavelet; we also show the point $\bm{x}_{\textup{trace}}$ as a cross. Top right plot: intermediate solution $\widetilde{u}(\cdot,7.5)$. Bottom left plot: intermediate solution $\widetilde{u}(\cdot,15)$. Bottom right plot: final solution $\widetilde{u}(\cdot,20)$.

Figure 11: Value of solution at point $\bm{x}_{\textup{trace}}=(-1,-2)$.

We show the resulting $\widetilde{u}(\cdot,t)$ for the four times $t\in\{0,7.5,15,20\}$ in Fig. 10. There, we can see why so many terms are necessary for the approximation of $u$: we must model many reflection and diffraction effects. Since energy escapes the system only through the top “door”, the wave will persist for quite a long time. Accordingly, a larger $T$ will make a larger $N$ necessary.

In order to better inspect this effect, we show the trace of the solution at the arbitrarily chosen point $\bm{x}_{\textup{trace}}=(-1,-2)$ in Fig. 11. We notice that oscillations persist for $t>10$. We use this last plot also to validate our results. To this aim, we compare:

• The surrogate $\widetilde{u}$ obtained as described above, with $\textup{tol}=2.5\cdot 10^{-2}$.

• The surrogate $\widetilde{u}$ obtained with our strategy, but with $\textup{tol}=10^{-3}$. This leads to an increased number of rays $N\approx 9.6\cdot 10^{3}$.

• The reference solution $u_{\textup{FE}}$ obtained by the P1-FE with leapfrog timestepping, as described at the beginning of Section 5. The mesh size must be chosen small enough to well resolve both the initial condition and the domain $\Omega$, as well as propagation of the solution itself. In our case, we found that a mesh with approximately $1.4\cdot 10^{6}$ elements leads to sufficient accuracy for a comparison with both above surrogate models $\widetilde{u}$. To satisfy the CFL condition on this mesh, we choose a time step $\Delta t\approx 7\cdot 10^{-3}$.

We can observe that the two surrogates obtained with our approach give very similar results. Indeed, the cutoff tolerance tol affects the results only for sufficiently large $t$, due to the accumulation of “small” waves that are excluded from the coarser surrogate but included in the finer one.

Moreover, taking the FE solution as reference, we see that most of the peaks of the surrogates are aligned with the FE ones (i.e., the “phase” of the wave is well approximated), but there are some noticeable discrepancies in their amplitudes. This is due to the fact that, in our approach, reflection is modeled exactly, whereas the magnitudes of the diffraction waves are only approximated. For this reason, we should not expect the amplitude error to get smaller if we reduce tol. The only effective way of improving the approximation would be using a more accurate diffraction model.

As a final result, we also report:

• The “construction” time, i.e., the time required to compute the numerical solution. For $\widetilde{u}$, this means executing Algorithm 1. For $u_{\textup{FE}}$, this means building the mesh, assembling the FE stiffness and (lumped) mass matrices, and carrying out the timestepping.

• The “evaluation” time, i.e., the time required to evaluate the numerical solution ($\widetilde{u}$ or $u_{\textup{FE}}$) at a single $(\bm{x},t)$-point.

They can be found in Table 2.

	$\widetilde{u}$ ($\textup{tol}=2.5\cdot 10^{-2}$)	$\widetilde{u}$ ($\textup{tol}=10^{-3}$)	$u_{\textup{FE}}$
Construction	$1.03\cdot 10^{1}$	$1.39\cdot 10^{2}$	$9.21\cdot 10^{1}$
Evaluation	$3.38\cdot 10^{-4}$	$4.71\cdot 10^{-3}$	$3.89\cdot 10^{-5}$

Table 2: Computational times (in seconds) for the room test case. To obtain more statistically significant results, each displayed time is the average over $3$ (resp. $10^{3}$) runs of the construction (resp. evaluation) phase with identical parameters.

We can observe the increased construction and evaluation times that result from decreasing tol. Moreover, we see that, in this example, our proposed approach is more competitive in the construction phase, but less so in the evaluation phase. The larger evaluation time of our surrogate is ultimately due to the nonlinearity of the functions that appear in $\widetilde{u}$. In this context, it may be surprising to note that, in our implementation, the most expensive step in evaluating $\widetilde{u}$ (taking about half of the online time) is determining whether an evaluation point is in the spatial supports $\Omega_{n}$ or not. On the other hand, evaluating the FE solution at a space-time point is an extremely cheap operation, essentially corresponding to a vector dot product.

This being said, we can notice a clear advantage of our approach in this context. Once the approximation $\widetilde{u}$ has been built, it can easily be evaluated at arbitrary locations in space and time, namely, not only at the final time $t=T$. This is because all the terms defining the approximate expansion Eq. 5 ($\Psi$, $\bm{\xi}_{n}$, etc.) are cheap to store in memory. One cannot usually do the same with the FE solution. Indeed, once the time-stepping has been carried out, the evaluation of $u_{\textup{FE}}$ at arbitrary instants before $t=T$ (by space-time interpolation) is possible only if the whole sequence of solution at all time-steps has been stored, a practice that is commonly avoided due to the (often) unfeasible memory requirements.

5.2.1 A time-harmonic source

One of the advantages of our approach is that it allows changing the source terms of the problem in a seamless way. Notably, under minor technical constraints (e.g., the spatial support of the new source term should not be larger than the old one), this kind of change does not require training a new surrogate.

To showcase this, we approximate the wave propagating from a time-harmonic point source at $\bm{x}=\bm{0}$ with angular frequency $\omega>0$. In our tests, we pick $\omega\in\{2\pi,10\pi\}$. To this aim, we define $u$ as the solution of the following ($\omega$-dependent) problem:

$$\begin{cases}\partial_{tt}u(\bm{x},t)=\Delta u(\bm{x},t)-\omega^{2}\sin(\omega t)g(\bm{x})\quad&\text{for }(\bm{x},t)\in\Omega\times(0,T),\\ u(\bm{x},0)=0\quad&\text{for }\bm{x}\in\Omega,\\ \partial_{t}u(\bm{x},0)=0\quad&\text{for }\bm{x}\in\Omega,\\ \partial_{\nu}u(\bm{x},t)=0\quad&\text{for }(\bm{x},t)\in\partial\Omega\times(0,T],\end{cases}$$ (18)

where $g$ denotes a narrow 2-dimensional Gaussian, centered at $\bm{0}$ and with standard deviation $0.05$.

As usual, we define $\Psi=\Psi(\rho,t)$ as the ($\omega$-dependent) solution of the free-space version of Eq. 18 in radial-temporal coordinates. To obtain an approximation for the wave $u$ generated by the time-harmonic source for an arbitrary $\omega$, it suffices to plug the corresponding $\Psi$ in each term of the surrogate $\widetilde{u}$ from the previous section! We show the results of our approximation in Fig. 12, where we can observe a very good agreement between approximation and reference FE solution.

Figure 12: Value of solution at point $\bm{x}_{\textup{trace}}=(-1,-2)$ for excitation frequencies $\omega=2\pi$ (top) and $\omega=10\pi$ (bottom). The reference FE solution is also included in both plots as a dashed line.

Figure 13: Approximation error at $\bm{x}_{\textup{trace}}=(-1,-2)$ for excitation frequencies $\omega=2\pi$ (top) and $\omega=10\pi$ (bottom). We also include the a priori error bound and the heuristic a posteriori error indicator using dotted and dashed lines, respectively. Note that, to avoid clutter, we show a “filtered” version of the error, using a window of width $0.1$: $\max_{|\delta|\leq 0.05}|\widetilde{u}(\bm{x}_{\textup{trace}},t+\delta)-u(\bm{x}_{\textup{trace}},t+\delta)|$.

We display the approximation error in Fig. 13. There, we also display an a priori error bound based on Eq. 17, as well as the heuristic a posteriori error indicator $|\widetilde{u}_{\text{GO}}-\widetilde{u}|$, with $\widetilde{u}_{\text{GO}}$ being the approximation obtained with our approach under the assumption of geometrical optics. The a priori error bound is piecewise-constant since time-independent terms are added to the error whenever new diffraction waves reach $\bm{x}_{\textup{trace}}$.

For $t<6$, no diffraction is present, the error indicators are both zero, since the only errors are due to discretization of the PDE. For $t>6$, the a priori estimator provides a reliable but pessimistic upper bound for the error. On the other hand, the a posteriori indicator manages to stay much closer to the actual error. Although not substantiated by theory, this error indicator seems quite effective at giving quantitative information on the error magnitude, without the need to compute a reference solution.

Note that using FE to approximate the wave $u$ requires carrying out a new simulation from scratch for every frequency to be studied. To this end, one needs to choose a mesh with $\omega$-dependent resolution: the mesh size should be small enough for the well-known pollution effect (see, e.g., [27]) to be absent.

A constraint on the mesh resolution is present also in our proposed approach. However, it only applies to the problem defining the free-space solution $\Psi$, which is 1-dimensional in space. Hence, having to refine the mesh represents a much smaller obstacle to efficiency. In particular, as the frequency $\omega$ increases, the computation of $\widetilde{u}$ becomes comparatively more and more efficient, with respect to the computation of $u_{\textup{FE}}$.

6 Conclusions

We have presented a surrogate modeling strategy for approximating waves propagating through complex 2-dimensional domains with polygonal boundaries. Our method relies on the automatic identification of reflection and diffraction effects caused by the domain geometry. Each effect is modeled through a relatively simple nonlinear expression. Reflection-related components are built using geometrical optics, whereas diffraction-related components are modeled by a novel ad hoc modification of the geometrical theory of diffraction. We have also provided estimators for the quality of the approximation achieved by our method.

In our numerical tests, we have observed a good accuracy, with the main features of the target wave being well identified. Notably, our diffraction model has proven to be fairly effective. Still, it relies on the parameter $\mu$, which, in some sense, determines the strength of surrogate diffraction effects. Although we have presented some heuristics for choosing $\mu$, more refined strategies for selecting $\mu$ and, more generally, a thorough validation of our diffraction model remain open issues.

In terms of complexity, our method requires the solution of a simplified 1D-in-space problem, much simpler than the original 2D-in-space one. Another favorable aspect of our algorithm is its potential to be run on parallel architectures, since the computation of different rays can be carried out independently. This is not the case for standard timestepping-based discretizations, due to their intrinsically sequential nature.

6.1 Open questions and possible extensions

We now present some open questions concerning our method, which give cues for further research on the subject.

Additional physics

A first question is whether our algorithm may be applied to wave equations incorporating further physical effects, like dissipation. In such setting, information can propagate faster than the wave speed due to diffusion. In our opinion, it should still be possible to apply our method to some effect, as long as dissipation is kept low, so as to maintain at least a partial notion of causality.

3D

Concerning other extensions, we have already mentioned that our geometry-based approach to wave propagation generalizes from 2D to 3D. In 3D, reflection and diffraction happen at facets and edges, as opposed to edges and vertices, respectively [20]. Moving from polygons to polyhedra leads to an increase in the “optical entities” to be considered (e.g., triangles have 3 edges and 3 vertices, whereas tetrahedra have 4 facets and 6 edges), especially if complex 2D surfaces need to be meshed. This leads to an increase in the number of field components to be considered. If the number of successive reflections/diffractions is large, this effect is further amplified. For instance, if at most 3 (resp., 5) successive reflection or diffraction events happen involving tetrahedral scatterers, we may expect the cost of our method to increase by a factor of 5 (resp., 13) with respect to a similar simulation with triangular scatterers in 2D.

Moreover, some additional implementation issues arise in 3D due to the higher spatial dimension. For instance, when computing the spatial supports $\Omega_{n}$, our current strategy, cf. Remark 2.3, could be rather expensive if implemented naively, especially if the number of facets of $\partial\Omega$ is large. More sophisticated strategies to represent the sets $\Omega_{n}$, e.g., involving beam-tracing [6] as opposed to ray-tracing, may be required to handle 3 dimensions more effectively. That being said, besides the geometrical side, our method still only requires solving a 1D FEM problem in spherical symmetry, even though the physical dimension of the underlying problem has increased.

Coupling with other methods

Our method is not too well suited to approximate wave propagation in trapping media, since the number $N$ of terms in the approximation may increase dramatically. To alleviate this issue, one might consider splitting the domain into non-trapping and trapping parts. Our approach could be then applied to approximate wave propagation in non-trapping subdomains, while standard MOR approaches (based, e.g., on projection onto a basis of standing waves) may deliver an efficient approximation over the trapping region. As in other MOR methods based on domain decomposition [28, 29], the most complex step is the (dynamic) coupling between subdomains. This remains an open issue for our method. We also note that, more generally, developing such an approach could also enable applying our method within multi-physics problems, e.g., when modeling the interaction between propagating wave and structures.

Parametrized problems

Finally, we recall that, in many applications, the ultimate target is understanding how the wave $u$ solving Eq. 1 depends on underlying parameters $\mathbf{p}$, e.g., the forcing term $f$, the shape of the domain $\Omega$, etc. In this setting, MOR tries to construct a surrogate model of the form $\widetilde{u}=\widetilde{u}(\bm{x},t;\mathbf{p})$, providing a good approximation of $u$ over a whole range of parameter values. Even though our technique was presented here in the non-parametric setting, we believe that it potentially allows incorporating the parameter dependence in a natural and efficient way. In our opinion, this might be achievable by leveraging the simple and interpretable structure of the field components (free-space solution, spatial support, and angular modulation). As a simple preliminary example, we showcased this in Section 5.2.1 for a parametric source term, with the parameter being the frequency. We are currently investigating how to extend our method to more complicated parametric problems.

Appendix A UTD diffraction coefficient

Consider the setup shown in Fig. 14: a wedge with exterior angle $\alpha$ and vertex $\bm{\xi}_{n}$ is hit by a wave coming from $\bm{\xi}$, located at an angular position $\phi=\theta$. The angle $\phi$ is measured from one of the two sides of the wedge. For simplicity, we assume that each of the wedge’s sides is either sound-soft or sound-hard. For a discussion of the UTD coefficients in the case of impedance boundary conditions, we refer to [30].

Let $\bm{x}$ be a point located at distance $s=\left\|\bm{x}-\bm{\xi}_{n}\right\|$ from $\bm{\xi}_{n}$, at an angle $\phi=\phi(\bm{x})$. The diffraction coefficient $D=D(\bm{x})$ represents the magnitude (and sign) of the diffraction wave at $\bm{x}$, assuming that the incident wave, with origin at $\bm{\xi}$, is time-harmonic with unit amplitude and frequency $k$. UTD [17] predicts a diffraction coefficient

$$D=D_{1}+D_{2}\pm(D_{3}+D_{4}),$$

with the sign $\pm$ depending on the type of boundary conditions (“$+$” for Neumann, “$-$” for Dirichlet). Given the wedge index $\nu:=\pi/(2\pi-\alpha)$, all contributions $D_{j}$ have the form

$$D_{j}=-\frac{\nu}{2\sqrt{2\pi ks}}\cot(\widehat{\phi}_{j})F\left(2ks\cos(\widetilde{\phi}_{j})^{2}\right),$$

with

$$\widehat{\phi}_{1}=\frac{\pi+\phi-\theta}{2}\nu,\quad\widetilde{\phi}_{1}= N_{1}(2\pi-\alpha)-\frac{\phi-\theta}{2},$$

$$\widehat{\phi}_{2}=\frac{\pi-\phi+\theta}{2}\nu,\quad\widetilde{\phi}_{2}= N_{2}(2\pi-\alpha)-\frac{\phi-\theta}{2},$$

$$\widehat{\phi}_{3}=\frac{\pi+\phi+\theta}{2}\nu,\quad\widetilde{\phi}_{3}= N_{3}(2\pi-\alpha)-\frac{\phi+\theta}{2},$$

$$\widehat{\phi}_{4}=\frac{\pi-\phi-\theta}{2}\nu,\quad\widetilde{\phi}_{4}= N_{4}(2\pi-\alpha)-\frac{\phi+\theta}{2},$$

and

$$F(x)=2\sqrt{x}\sqrt{\left(\int_{\sqrt{x}}^{\infty}\cos(y^{2})\textrm{d}y\right)^{2}+\left(\int_{\sqrt{x}}^{\infty}\sin(y^{2})\textrm{d}y\right)^{2}}.$$

The integers $N_{j}$ are defined as

$$N_{1}=\left[\frac{\nu}{2}\right],\quad N_{2}=\left[-\frac{\nu}{2}\right],\quad N_{3}=\left[\frac{1+\nu}{2}\right],\quad\text{and}\quad N_{4}=\left[\frac{1-\nu}{2}\right],$$

with square brackets denoting the rounding operator, which yields the integer that is closest to its argument.

Figure 14: Diagrams for the two cases of scattering for concave corners ($0<\alpha<\pi$): without (left plot) and with shadow zone (right plot). The dashed lines are reflection shadow boundaries. The dash-dotted line is an incident shadow boundary. Shadow regions are absent if and only if $\pi-\alpha\leq\theta\leq\pi$. The angular coordinate $0<\phi<2\pi-\alpha$ is measured starting from one of the two adjacent edges of $\partial\Omega$.

The diffraction coefficient has unit jumps (with suitable signs) at the shadow boundaries, i.e., at the angular positions corresponding to boundaries of the spatial support of the source wave (in which case, we have an “incident shadow boundary”, ISB) or of one of its reflections (in which case, we have a “reflection shadow boundary”, RSB). The locations of the shadow boundaries can be identified geometrically. For instance, for convex wedges ($0<\alpha<\pi$), shadow boundaries are at $\phi=\left|\pi-\theta\right|$ (which is an ISB if $\theta>\pi$, and an RSB otherwise) and at $\phi=2\pi-\alpha-\left|\pi-\alpha-\theta\right|$ (which is an ISB if $\theta<\pi-\alpha$, and an RSB otherwise).

The diffraction coefficient for $\alpha=\frac{\pi}{3}$, $\theta=\frac{\pi}{5}$, Neumann conditions, and $ks\in\{1,10,100\}$ are shown in Fig. 15 (left). Two jumps happen, at an RSB (at $\phi=\frac{4}{5}\pi$) and at an ISB (at $\phi=\frac{6}{5}\pi$). Similar results are obtained for Dirichlet conditions, cf. Fig. 15 (right).

We can observe that the value of $ks$ determines how quickly $D$ decays to 0 around each shadow boundaries: larger values of $ks$ yield a narrower “transition region”, which, loosely speaking, is the support of $D$ around each discontinuity. Indeed, asymptotically in $ks$, the width of each transition region is proportional to $1/\sqrt{ks}$ [20].

Figure 15: UTD diffraction coefficients for Neumann (left) and Dirichlet (right) wedges.