The on-surface radiation conditions (OSRCs), originally developed by Kriegsmann, Taflove and Umashankar [1], are semi-analytical formulations to roughly approximate the solution of wave scattering problems. The approximate nature of an OSRC prevents it from being used as a satisfactory solution procedure. However, this semi-analytical construct yields tremendous computational speed which is useful to explore feasibility regions 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 fixed-point or Krylov subspace methods [2][3][4][5]. 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 [6][7][8][9][10]. Other improvements have been accomplished by Jones [11][12][13], Ammari [14,15], Calvo et al. [16,17], , Chaillat et al. [21,22] and Darbas et al. [2][3][4][5]23]. See also [24,25,[25][26][27][28][29][30][31][32].
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 [25] and by Alzubaidi, Antoine and Chniti in [31] for multiple scattering problems where the scatterer consists of several disjoint obstacles. The formulations in [25,31] 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 [6,25,32], 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 as well. In Section 6, we implement the OSRC as a preconditioner to reduce the spectral radius of a block-Jacobi iteration matrix to solve a multiple scattering problem using a boundary element method. Limitations, future work and conclusions are discussed in Section 7.
We seek to approximate an outgoing solution to the Helmholtz equation in the exterior to a set of J impenetrable convex obstacles. Each obstacle is enclosed by a closed smooth surface ∂Ω j that separates the space into a simply connected bounded domain Ω - j and an exterior domain Ω + j . We also define
The wave field is assumed to satisfy the following boundary value problem,
where the wavenumber k > 0 is constant, ι is the imaginary unit, r = |x|, and ∂Ω is the boundary of Ω + . The proposed approach to roughly approximate the solution u of the boundary value problem ( 2) is based on the following theorem whose objective is to express the solution u as the sum of purely-outgoing wave fields u j each radiating from a single surface ∂Ω j for j = 1, 2, . . . , J . A proof is found in [33].
Theorem 1. Let the closures of Ω - j and of Ω - i be mutually disjoint for i ̸ = j, and let u be a radiating solution to the Helmholtz equation in Ω + . Then u can be uniquely decomposed into purely-outgoing wave fields u j for j = 1, 2, . . . , J such that
where u j radiates only from ∂Ω j , that is,
Using Green's theorem for each of the purely-outgoing wave fields u j , we obtain a representation of the solution
where Φ is the fundamental solution to the Helmholtz equation. If the exact outgoing Dirichlet-to-Neumann (DtN) map Λ j on each surface ∂Ω j is available, then ∂ ν u j = Λ j u j and we obtain
The representation (6) renders the solution in Ω + provided that the Dirichlet data of each purely-outgoing wave field u j is known. This data can be found by imposing the boundary condition (2b) on each boundary ∂Ω j which leads to the following system of equations,
where f i is the evaluation of f on the surface ∂Ω i , and where the propagation of the wave field u j from the surface ∂Ω j to the surface ∂Ω i for j ̸ = i 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 Λ j must be evaluated or at least roughly 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.
We seek to derive an analytical formula to roughly approximate the propagation of a wave field away from a generic convex surface Γ . It is conceptually similar to the beam propagation method developed in [34] for twodimensional scattering problems. 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 an expansive foliation of parallel surfaces generated by the smooth surface Γ . We follow closely the approach from [6]. The domain exterior to Γ is denoted Ω . This family of surfaces Γ s parametrized by s ≥ 0 is defined as follows
where n(x) is the outward normal vector of the surface Γ at the point x ∈ Γ . Notice Γ = Γ 0 . Since Γ is convex, the family of surfaces Γ s foliates Ω . For any y ∈ Ω there exists unique x ∈ Γ and unique s > 0 such that y = x+sn(x). Moreover, the outward normal vector n(y) of the surface Γ s at a point y = x + sn(x) coincides with the normal vector n(x).
Now we proceed to write the Helmholtz differential operator in terms of a tangential system of coordinates on the surface Γ s
Here H s is the mean curvature of the surface Γ s and ∆ Γ s is the Laplace-Beltrami operator on the same surface. The differential ∂ s represents the derivative in the outward normal direction on Γ s . See [6] for a concise review of the differential geometry of surfaces, including the definition of curvatures and the Laplace-Beltrami operator. The mean curvature H s and Gauss curvature K s of the surface Γ s satisfy the following relations,
where the symbol µ s stands for
Here H and K are mean and Gauss curvature of the surface Γ , respectively. Also, R is the curvature tensor and I is the identity on the tangent plane. The curvature tensor R can be represented as a self-adjoint matrix with eigenvalues κ 1 and κ 2 , called the principal curvatures, such that
The Laplace-Beltrami operator satisfies
where this latter approximation is valid when κ 1 ≈ κ 2 so that µ s (I + sR) -2 ≈ I. Now we seek to decompose the wave field propagating across a surface Γ s into incoming and outgoing components using Nirenberg's factorization theorem [6,39]. There are two pseudo-differential operators Λ -(s) and Λ + (s) of order +1, such that
We refer to Λ + (s) and Λ -(s) as the outgoing and incoming DtN operators for the radiating boundary value problem defined in the exterior of Γ s . The operator Λ + s admits the following second order approximation derived in [6],
If w is purely-outgoing from the surface Γ , then ∂ s w = Λ + (s)w. Using the approximation ( 16) for the outgoing DtN operator, we obtain that w satisfies the following differential equation in s > 0,
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 ln
which leads to
where we have neglected the tangential derivatives of H and K in order to factor the exponential. Here w o represents the trace of w 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 k,
This implicit approximation is chosen to preserve the stability (boundedness) of the operations. In fact, the spectrum of the operator in (19) lies in the unit circle because the Laplace-Beltrami ∆ Γ operator has real eigenvalues. Now we go back to the multi-scattering problem with its notation formulated in Section 2. In order to propagate the wave field u j from the surface ∂Ω j onto the surface ∂Ω i , we approximate the propagator (8) as follows,
where y = x+sn(x). Due to the assumed convexity of ∂Ω j , the point x ∈ ∂Ω j and s > 0 are determined uniquely by y ∈ ∂Ω i satisfying (10). The curvatures H and K, and the Laplace-Beltrami operator ∆ ∂Ω j correspond to the surface ∂Ω j . Regarding the boundedness of (20), notice that the two exponentials have purely imaginary arguments, the operator ( 19) is non-expansive, and the factor µ s (x) > 1 for s > 0 and for convex ∂Ω j (H > 0, K > 0). Therefore each propagator P i j can be shown to be a contraction mapping. This renders the OSRC system (7) uniquely solvable by the Neumann series theorem in an L 2 framework.
It is commonly of interest to evaluate the radiating field u away from the obstacles. The asymptotic behavior is characterized by the far-field pattern u ∞ of the radiating field u which is given by
where r = |y| and ŷ = y/|y|. From the superposition of the purely-outgoing fields, we find that
so that it only remains to find the far-field pattern u ∞ j of each purely-outgoing field u j . We proceed with the following asymptotic expansions valid as r → ∞,
Using ( 18) and ( 23), we find that the purely-outgoing wave field u j has the following asymptotic behavior
] which implies that the far-field pattern corresponding to u j is given by
where x ∈ ∂Ω j is uniquely determined by ŷ such that n(x) = ŷ. Hence, once the Dirichlet data of the field u j is found by solving the system (7), then plugging ( 24) into ( 22) renders an approximation for the far-field pattern u ∞ .
In this section, we propose a numerical implementation of the on-surface radiation condition defined by the system (7) where the propagation operator P i, j is roughly approximated by (20). The approach is based on discretizing the Laplace-Beltrami operator, the mean and Gauss curvatures using a triangulation of the surfaces ∂Ω j for j = 1, 2, . . . , J . Our guiding references for approximating geometrical properties and operators on triangulated surfaces are [40][41][42][43][44][45].
The discrete Laplace-Beltrami operator is described as follows. Let {x j } J j be the collection of vertices of the triangulation T of the surface Γ . For a fixed vertex x j , let {x j(i) } i=1,2,... be the neighboring vertices of x j , and let {e j(i) } i=1,2,... be the edges of the triangulation that connect x j and x j(i) . Now for each edge e j(i) , let α j(i) and β j(i) be the angles opposing the edge e j(i) in the two triangles that share that edge. For a smooth function u defined on the surface Γ , we use the following discrete Laplace-Beltrami operator,
where A j is the area associated with the vertex x j defined as 1 /3 of the total area of the triangular elements sharing the vertex x j .
For the discrete mean curvature, we first define the discrete mean curvature vector as
Then the mean curvature is given by
where n(x j ) is the normal vector at the vertex x j .
The discrete Gauss curvature at each vertex x j of the triangulation, is defined as follows [40,44,46]. Let θ i be the angle between two successive edges sharing the vertex x j . Then, the discrete Gauss curvature is given by
where A j is the area associated with the vertex x j as defined above.
Given the triangulation T j of the surface ∂Ω j , we have the necessary definitions in Section 5.1 to implement a discrete version of the propagator (20) to pose the system (7) at the discrete level. The discrete operator P i, j propagates the waves from the triangulation T j of the surface ∂Ω j to the triangulation T i of the surface ∂Ω i . Let the triangulations T j and T i have N j and N i vertices, respectively. Then the propagation operator can be represented as a matrix P i, j : R N j → R N i . We proceed to describe the construction of this matrix. Let y n ∈ T i where n is the index in the list of vertices of T i . Then we can find x m ∈ T j and s m > 0 as the discrete version of (10), that is,
Then the matrix P i, j can be written as
where ∆ T j is an N j × N j matrix defined by (25) which approximates the Laplace-Beltrami operator, and A i, j is an N i × N j matrix. Each row of this matrix has exactly one non-zero entry. For the nth row, this non-zero entry is at the mth column and it is given by
where m depends on n by satisfying (29)
where the total number of degrees of freedom is
The standard numerical techniques for BIE, such as Nystrom, collocation and Galerkin boundary element methods, lead to fully populated matrices. Similarly, for the OSRCs developed in [25,31], matrices with fully populated off-diagonal blocks are obtained. In these cases, high computational costs O(N 2 ) and large memory demands O(N 2 ) are needed for matrix-vector multiplication in any iterative solver (e.g. fixed-point or GMRES). By contrast, the proposed OSRC leads to O(N ) 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 A i, j after the inversion of a system governed by I + ∆ T j /(2ιkH). The former is a sparse matrix, with exactly one non-zero entry per row, leading O(N i ) 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 O(N j ) operations. As a consequence, any iterative solver for the system (32) will require O(N ) operations per iteration. 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 ∂Ω j for all j = 1, 2, . . . , J , the norm of the propagation operator P i, j decreases as the distance between the ith and jth obstacles increases. Hence, the system (32) can be solved using the following fixed-point iteration,
for vanishing initial guess at m = 0.
In this section we present a few numerical results obtained from the implementation of the discrete formulation defined in the previous subsections. To obtain these numerical results, we implemented the fixed-point method (33) with M = 10 iterations. 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, Sphere:
Rounded Cube:
Marshmallow:
Coarse triangulations of these three surfaces are shown in Fig. 2. The discrete mean and Gauss curvatures of the rounded cube and of the marshmallow, defined by ( 26)-( 27), are shown in Fig. 3.
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 f as the boundary trace of a radiating wave field. We consider the field,
where Φ(z) = e ιk|z| /(4π |z|) is the outgoing fundamental solution to the Helmholtz equation with frequency k > 0. Hence, F represents the superposition of J point-sources with respective locations at c j . If each point c j is enclosed by the respective surface ∂Ω j , for j = 1, 2, . . . , J , then the exact solution to the boundary values problem (2) with Dirichlet data f = F| ∂Ω is given by F| Ω + .
Table 1 L 2 -norm relative error for the far-field pattern for various wavenumbers and mesh refinements. Multiple-scattering from two spheres centered at the points (36). The DOF approximately quadruples with each mesh refinement.
Two spheres Now we present some numerical results for the wave radiating problem in the exterior of J = 2 unit-spheres. These spheres are centered at the follow points:
Table 1 displays the L 2 -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.
Figs. 4 and 5 show the exact and numerical far-field patterns and the error between them, for wavenumbers k = π and k = 2π , respectively.
For visual comparison, cross-sections of the far-field pattern are shown Figs. 6 and7 for wavenumbers k = π and k = 2π , respectively.
Now we choose J = 3 surfaces shown in Fig. 2, translated to the respective points c 1 = 4 (1, 0, 0), c 2 = 4 (cos 2π /3, sin 2π /3, 0), and c 3 = 4 (cos 4π /3, sin 4π /3, 0).
Table 2 displays the L 2 -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. Figs. 8 and9 show the exact and numerical far-field patterns and the error between them, for wavenumbers k = π and k = 2π, respectively.
For this example, we also display some cross-sections of the far-field pattern as shown in Figs. 10 and 11 for wavenumbers k = π and k = 2π , respectively. L 2 -norm relative error for the far-field pattern for various wavenumbers and mesh refinements. Multiple scattering from three surfaces (34) centered at the points (37) The DOF approximately quadruples with each mesh refinement.
We demonstrate the utility of the proposed local OSRC formulation to precondition the Galerkin boundary element method (BEM) for a multiple scattering problem. This is one of the main applications for OSRCs [2][3][4][5]21,22,47,48]. The implementation of the Galerkin BEM and the OSRC was done using GYPSILAB, an open source MATLAB library developed by Alouges and Aussal [49]. We choose a multiple scattering problem in threedimensional space with two sound-soft unit spheres centered at (x, y, z) = (2, 1, 0) and (x, y, z) = -(2, 1, 0). The incident field u inc is a plane wave of wavenumber k = 4π propagating along the x-axis in the positive direction. We first pose the multiple scattering problem on the single-layer direct integral operator framework [50][51][52] governed
where ρ j are the unknown boundary densities on ∂Ω j , f j = -u inc | ∂Ω j is the Dirichlet data on each sphere ∂Ω j , and the single-layer operators S i j are defined as At the continuous level, the system (38) can be solved using a block-Jacobi fixed-point iterative method [47,50,51] of this form
where v j = S j j ρ j . We precondition this iterative method using the proposed OSRC formulation (32) in order to reduced the spectral radius (contraction mapping constant) of the iteration matrix for the fixed-point algorithm. The preconditioned algorithm is written as follows
Notice that this iterative method requires the inversion of the system (32) which is accomplished running the iterative algorithm (33) nested within each iteration of the algorithm (41).
For the BEM numerical implementation, we first discretize each sphere using triangular meshes. We ran experiment using two meshes, one with N = 2500 vertices (coarse mesh) and the other with N = 5000 vertices (fine mesh). The operators in (40) and (41) were discretized using the Galerkin projection on continuous piecewise linear functions. The boundary element implementation of the OSRC is easily obtained with the GYPSILAB toolbox. The discretization in the weak form of the propagator (30) corresponds to invert a two-dimensional elliptic system governed by the Laplace-Beltrami operator on each surface, and the function A i j in (31) is pre-computed offline on the Gauss quadrature points of the meshes. The numerical approximations for the total acoustic field on the x y cross-sectional plane using the BEM and OSRC formulation are displayed in Fig. 12. We note that the result from the OSRC formulation is a good approximation to the BEM solution, but not entirely accurate especially in the shadow regions of the spheres. Hence, instead of using the OSRC formulation to obtain a final solution, we implemented it as a right-preconditioner as shown in the algorithm (41). The convergence pattern and computational cost for both (40) and (41) are displayed in Fig. 13. This procedure reduced the spectral radius of the iteration matrix Fig. 12. Evaluation on the x y-plane of the numerical solutions using the BEM and OSRC formulation for the scattering the plane wave from two sound-soft spheres. Fig. 13. Performance and computational cost for the block-Jacobi fixed point iterations (Jacobi) and OSRC-based preconditioned block-Jacobi fixed point iterations (OSRC+Jacobi) to solve the boundary element system. The reported CPU time is normalized to the CPU time needed for one iteration of the preconditioned block-Jacobi method. by a factor of 300 approximately, which leads to a significant improvement on the rate of convergence due to the OSRC preconditioner.
For the coarse and fine meshes, the OSRC preconditioner adds about 27% and 7% more computational work, respectively, to each iteration of the block-Jacobi method. This reduction in the relative added computational burden as the mesh is refined is consistent with the linear complexity O(N ) of the OSRC formulation (described in Section 5.3). This extra marginal cost pays off in terms of faster convergence. See the left panel in Fig. 13. In fact, to reach a desired level of error, the preconditioned method requires less than half of the total computational work needed for the unpreconditioned block-Jacobi method to reach the same error level. See the right panel in Fig. 13.
We have formulated an OSRC for multiple scattering problems that involves local operations only. Integration over surfaces is avoided which leads to sparse matrices upon discretization and implementation with O(N ) complexity. We expect that this OSRC will provide an inexpensive initial guess or preconditioner to solve multiple scattering problems using boundary integral equations (BIE) or a fast approximate method to explore parameter spaces for optimization problems. The derivation of adequate preconditioners for BIE is a mandatory requirement to solve large scale wave scattering problems because the discretization of BIE leads to notoriously ill-conditioned, non-hermitian dense matrices [2][3][4][5]21,22,47,48]. In contrast to some algebraic preconditioners for dense systems (such as splitting, incomplete factorization, and sparsification), the OSRC formulation is based on the physics of the problem, including the wave reflections between the multiple obstacles. This renders the OSRC as an effective preconditioner for this specific problem due to its capacity to account for the global interdependence patterns among entries of the system.
The major limitations of the proposed formulation are similar to those of other OSRCs. This OSRC provides a rough approximation of the solution without a-priori error estimates and lacking a methodical procedure for convergence. Moreover, to obtain closed-form formulas in Section 3, we adopted several assumptions: (i) convexity of the surfaces bounding the obstacles, (ii) that the principal curvatures of these surfaces are similar, (iii) that their tangential derivatives are small, and (iv) that the frequency is sufficiently high. Therefore, we expect this OSRC to lose accuracy for surfaces with skewed or rapidly changing curvatures, or for low-frequency problems. As a consequence, the OSRC should not be employed as a viable solution procedure, but rather as a fast pre-processing step for preconditioning or exploration.
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 approximation of the Dirichlet-to-Neumann map to improve the accuracy of the OSRC [32]. 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 surfaces. The same is true for piecewise smooth surfaces in order to handle edges and corners. However, recent progress has been reported in this area by Modave, Geuzaine and Antoine [53].