The on-surface radiation conditions (OSRCs), originally developed by Kriegsmann, Taflove and Umashankar Kriegsmann1987, are semi-analytical methods to approximate the solution of wave scattering problems. In general, OSRCs lead to rough approximations, meant to sacrifice accuracy in favor of tremendous computational speed. This property is useful to explore parameter spaces for optimization problems, to produce a good initial guess for iterative methods, and to design inexpensive preconditioners to solve boundary integral equations (BIE) numerically using Krylov subspace methods Antoine2004; Antoine2005; Antoine2007; Darbas2013. The OSRCs have been studied by many researchers including Antoine and collaborators who formulated them in a differential geometric setting in order to handle non-canonical surfaces Antoine1999; Antoine2001; Antoine2001b; Antoine2006; Antoine2008. Other improvements have been accomplished by Jones Jones1988; Jones1990; Jones1992, Ammari Ammari1998; Ammari1998b, Calvo et al. Calvo2004; Calvo2003, Barucq et al. Barucq2003; Barucq2010a; Barucq2012, Chaillat et al. Chaillat2015; Chaillat2017 and Darbas et al. Antoine2004; Antoine2005; Antoine2007; Darbas2013; Darbas2006. See also Atle2007; Acosta2015f; Stupfel1994; Murch1993; Teymur1996; Medvinsky2010; Chniti2016; Alzubaidi2016; Acosta2015f; Acosta2017c.
The underlying assumption in most of the aforementioned works is that the solution radiates outwardly at every point on the surface of the scatterer. This assumption may be violated when the scatterer is not convex since the wave field may exhibit reflections from the scatterer to itself. Hence, one of the main drawbacks of the OSRC method is that its formulation and accuracy are limited to convex scatterers. This issue was addressed by the author in Acosta2015f and by Alzubaidi, Antoine and Chniti in Alzubaidi2016 for multiple scattering problems where the scatterer consists of several disjoint obstacles. The formulations in Acosta2015f; Alzubaidi2016 successfully account for the propagation of waves from one obstacle to another. However, the wave propagation formulas are non-local. Due to the use of boundary integral operators, integration of the wave fields over the surface of the obstacles is required. As a consequence, matrices with dense blocks are obtained upon discretization, which leads to high memory demands and costly matrix inversion, especially in the three-dimensional setting and also for high frequencies.
Our main goal is to modify the formulation of the OSRC for multiple obstacles in such a way that its discretization leads to sparse matrices. This is accomplished by defining a local wave propagation formula to transmit the influence of one obstacle on another using differential rather than integral operators. In Section 2 we provided the mathematical formulation of the multiple scattering problem and decompose it into a system of single-scattering problems. Based on Antoine1999; Acosta2015f; Acosta2017c, we derive a local formula for the propagation of waves which allows us to account for the wave reflections between obstacles. This formula is developed in Section 3. In Section 4 we use the propagation formula to evaluate the far-field pattern using local operations as well. In Section 5, we propose how to numerically implement the OSRC based on a triangulation of the surfaces. A few numerical results and comparisons with exact solutions are shown in Sections 6. Limitations, future work and conclusions are discussed in Section 7.
2 Mathematical formulation
We seek to approximate an outgoing solution to the Helmholtz equation in the exterior to a set of impenetrable convex obstacles. Each obstacle is enclosed by a closed smooth surface that separates the space into a simply connected bounded domain and an exterior domain . We also define
The wave field is assumed to satisfy the following boundary value problem,
where the wavenumber is constant, is the imaginary unit, , and is the boundary of . The proposed approach to approximate the solution of the boundary value problem (2) is based on the following theorem whose objective is to express the solution as the sum of purely-outgoing wave fields each radiating from a single surface for . A proof is found in Balabane2004.
Let the closures of and of be mutually disjoint for , and let be a radiating solution to the Helmholtz equation in . Then can be uniquely decomposed into purely-outgoing wave fields for such that
where radiates only from , that is,
Using Green’s theorem for each of the purely-outgoing wave fields , we obtain a representation of the solution
where is the fundamental solution to the Helmholtz equation. If the exact Dirichlet-to-Neumann (DtN) map on each surface is available, then and we obtain
The representation (6) renders the solution in provided that the Dirichlet data of each purely-outgoing wave field is known. This data can be found by imposing the boundary condition (2b) on each boundary which leads to the following system of equations,
where is the evaluation of on the surface , and where the propagation of the wave field from the surface to the surface for is given by the following operator,
For practical purposes, the evaluation of the propagation operator in (8) has two main challenges. First, the DtN map must be evaluated or at least approximated which is the objective of OSRC in general. Second, the propagation operator involves integration over the boundary of each obstacle. This leads to the appearance of dense blocks once the governing system (7) is discretized. Our objective in the following section is to mitigate this latter problem.
3 Local formula for propagation of wave fields
We seek to derive an analytical formula for the propagation of a wave field away from a generic convex surface . This formula will allow us to account for the interaction between obstacles in order to approximate the system (7) using sparse matrices. This formula is for the propagation of waves on a expansive foliation of parallel surfaces generated by the smooth surface . We follow closely the approach from Antoine1999. The domain exterior to is denoted . If is a convex surface, then it generates a family of parallel surfaces parametrized by defined as follows
is the outward normal vector of the surfaceat the point . Notice . Since is convex, the family of surfaces foliates . For any there exists unique and unique such that . Moreover, the outward normal vector of the surface at a point coincides with the normal vector . See details in (Kress-Book-1999, §6.2), (Gil-Tru-2001, §14.6), (DoCarmo1976, Probl. 11 §3.5) and Antoine1999. The point can be regarded as the projection of y on the surface , and is the distance from x to y (Deu-2001, Ch. 2). The pair satisfies
Now we proceed to write the Helmholtz differential operator in terms of a tangential system of coordinates on the surface
Here is the mean curvature of the surface and is the Laplace–Beltrami operator on the same surface. The differential represents the derivative in the outward normal direction on . See Antoine1999 for a concise review of the differential geometry of surfaces, including the definition of curvatures and the Laplace–Beltrami operator. The mean curvature and Gauss curvature of the surface satisfy the following relations,
where the symbol stands for
Here and are mean and Gauss curvature of the surface , respectively. Also,
is the curvature tensor andis the identity on the tangent plane. The curvature tensor
can be represented as a self-adjoint matrix with eigenvaluesand , called the principal curvatures, such that
The Laplace-Beltrami operator satisfies
where this latter approximation is valid when so that .
Now we seek to decompose the wave field propagating across a surface into incoming and outgoing components using Nirenberg’s factorization theorem Nirenberg1973; Antoine1999. There are two pseudo–differential operators and of order , such that
We refer to and as the outgoing and incoming DtN operators for the radiating boundary value problem defined in the exterior of . The operator admits the following second order approximation derived in Antoine1999,
If is purely-outgoing from the surface , then . Using the approximation (16) for the outgoing DtN operator, we obtain that satisfies the following differential equation in ,
This is a separable differential equation. In order to solve it approximately in closed-form, we need the following expressions
where we have employed the approximation in (14) for the third integral. Hence, we obtain
which leads to
where we have neglected the tangential derivatives of and in order to factor the exponential. Here represents the trace of on the boundary . In order to maintain the locality of the evaluation, we further approximate the exponential of the Laplace–Beltrami operator using an implicit linear approximation valid for large frequency ,
This implicit approximation is chosen to preserve the stability (boundedness) of the operations. Notice that the implementation of (19
) corresponds to solving an elliptic partial differential equation on the surface.
Now we go back to the multi-scattering problem with its notation formulated in Section 2. In order to propagate the wave field from the surface onto the surface , we approximate the propagator (8) as follows,
where . Due to the assumed convexity of , the point and are determined uniquely by satisfying (10). The curvatures and , and the Laplace–Beltrami operator correspond to the surface .
4 Far-Field Pattern
It is commonly of interest to evaluate the radiating field away from the obstacles. The asymptotic behavior is characterized by the far-field pattern of the radiating field which is given by
where and . From the superposition of the purely-outgoing fields, we find that
so that it only remains to find the far-field pattern of each purely-outgoing field . We proceed with the following asymptotic expansions valid as ,
which implies that the far-field pattern corresponding to is given by
5 Discrete Implementation
In this section, we propose a numerical implementation of the on-surface radiation condition defined by the system (7) where the propagation operator is approximated by (20). The approach is based on approximating the Laplace–Beltrami operator, the mean and Gauss curvatures using a triangulation of the surfaces for . Our guiding references for approximating geometrical properties and operators on triangulated surfaces are Xu2006a; Xu2004; Xu2004a; Meyer2003; Magid2007; Bobenko2007.
5.1 Discrete Laplace–Beltrami operator and curvatures
The discrete Laplace–Beltrami operator is described as follows. Let be the collection of vertices of the triangulation of the surface . For a fixed vertex , let be the neighboring vertices of , and let be the edges of the triangulation that connect and . Now for each edge , let and be the angles opposing the edge in the two triangles that share that edge. For a smooth function defined on the surface , we use the following discrete Laplace–Beltrami operator,
where is the area associated with the vertex defined as 13 of the total area of the triangular elements sharing the vertex .
For the discrete mean curvature, we first define the discrete mean curvature vector as
Then the mean curvature is given by
where is the normal vector at the vertex .
The discrete Gauss curvature at each vertex of the triangulation, is defined as follows Xu2006a; Magid2007; Surazhsky2003. Let be the angle between two successive edges sharing the vertex . Then, the approximate Gauss curvature is given by
where is the area associated with the vertex as defined above.
5.2 Discrete wave propagator
Given the triangulation of the surface , we have the necessary definitions in Subsection 5.1 to implement a discrete version of the propagator (20) to pose the system (7) at the discrete level. The discrete operator propagates the waves from the triangulation of the surface to the triangulation of the surface . Let the triangulations and have and vertices, respectively. Then the propagation operator can be represented as a matrix . We proceed to describe the construction of this matrix. Let where is the index in the list of vertices of . Then we can find and as the discrete version of (10), that is,
Then the matrix can be written as
where is an matrix defined by (25) which approximates the Laplace–Beltrami operator, and is an matrix. Each row of this matrix has exactly one non-zero entry. For the row, this non-zero entry is at the column and it is given by
where the total number of degrees of freedom is.
The standard numerical techniques for BIE, such as Nystrom, collocation and boundary element methods, lead to fully populated matrices. Similarly, for the OSRCs developed in Acosta2015f; Alzubaidi2016, matrices with fully populated off-diagonal blocks are obtained. In these cases, high computational costs and large memory demands are needed for matrix-vector multiplication in any iterative solver (e.g. fixed-point or GMRES). By contrast, the proposed OSRC leads to complexity. This follows from the propagator matrix (30) that defines the governing system (32). The implementation of the propagator (30) requires the multiplication by the matrix after the inversion of a system governed by . The former is a sparse matrix, with exactly one non-zero entry per row, leading operations per matrix-vector multiplication. The latter is the solution for the discrete version of a two-dimensional elliptic equation. With proper sorting of the triangulation nodes, this system is both sparse and narrowly banded. Therefore, a direct solver such as LU decomposition requires operations. As a consequence, any iterative solver for the system (32) will require operations per iteration. Figure 1 displays (a) the sparsity pattern of the matrix for the configuration of two spheres from subsection 6.1, (b) the sparsity pattern of the discrete Laplace-Beltrami operator , and (c) CPU time for the application of the propagator (30) as a function of DOF.
For the numerical results presented in the next section, we apply a fixed-point iteration to solve the system (32). Due to the convexity of the surfaces for all , the norm of the propagation operator decreases as the distance between the and obstacles increases. Hence, the system (32) can be solved using the following fixed-point iteration,
for vanishing initial guess at .
6 Numerical Results
In this section we present a few numerical results obtained from the implementation of the discrete formulation defined in Section 5. We consider three convex surfaces for the numerical experiments presented here. These surfaces are the unit sphere, a rounded cube, and a marshmallow–like surface. Their defining equations are as follows,
For comparison, we manufacture an exact solution to the boundary values problem (2), valid for any geometry of the surfaces in . This is accomplished by defining boundary data as the boundary trace of a radiating wave field. We consider the field,
where is the outgoing fundamental solution to the Helmholtz equation with frequency . Hence, represents the superposition of point-sources with respective locations at . If each point is enclosed by the respective surface , for , then the exact solution to the boundary values problem (2) with Dirichlet data is given by .
6.1 Example 1: Two spheres
Now we present some numerical results for the wave radiating problem in the exterior of unit-spheres. These spheres are centered at the follow points:
Table 1 displays the -norm relative error for the far-field pattern for various wavenumbers and mesh refinements. The DOF approximately quadruples with each mesh refinement, which corresponds to halving the number of points per wavelength.
6.2 Example 2: Three different shapes
Now we choose surfaces shown in Figure 2, translated to the respective points
Table 2 displays the -norm relative error for the far-field pattern for various wavenumbers and mesh refinements of these three shapes. The DOF approximately quadruples with each mesh refinement, which corresponds to halving the number of points per wavelength in any tangential direction.
7 Final Remarks
We have formulated an OSRC for multiple scattering problems that involves local operations only. As opposed to Acosta2015f; Alzubaidi2016, integration over the surfaces is avoided which leads to sparse matrices upon discretization and implementation with complexity. The proposed method is a local OSRC that simultaneously accounts for the outgoing behavior of the solution as well as the wave reflections between the multiple obstacles. We expect that this OSRC will provide an inexpensive, physics-based preconditioner to solve multiple scattering problems using boundary integral equations (BIE) via Krylov subspace methods. The derivation of adequate preconditioners for BIE is a mandatory requirement to solve large scale wave scattering problems Antoine2004; Antoine2005; Antoine2007; Darbas2013; Darbas2015; Chaillat2015; Chaillat2017.
The major limitations of the proposed method are similar to those of other OSRCs. The formulation is valid for convex surfaces and the OSRC provides a rough approximation of the solution without a-priori error estimates. Therefore, the OSRC should not be employed if a high degree of accuracy is needed. However, as seen from Table1 and 2, and Figures 6-7 and 10-11, the proposed OSRC can produce qualitatively good approximations to the solution of the boundary value problem (2). Possible future work includes the extension to electromagnetic and elastic waves that govern important engineering applications. It is also possible to increase the order of approximations of the Dirichlet–to–Neumann map to improve the accuracy of the OSRC Acosta2017c. Here we have considered the second order differential approximation (16). The major challenge that still remains to be addressed is the formulation of an OSRC and the associated propagation formula for non-convex obstacles. The same is true for piecewise smooth surfaces in order to handle edges and corners.
This work was partially supported by NSF grant DMS-1712725. The author would like to thank Texas Children’s Hospital for its support and for the research-oriented environment provided by the Predictive Analytics Laboratory.