The numerical simulation of time-harmonic waves scattered from perfect electric conductors (PECs) is of fundamental importance across the spectrum of electromagnetic applications.
Denote by an incident field. We are looking for the solution of the exterior scattering problem, satisfying
Here, denotes the wavenumber of the problem, with denoting the frequency and and the electric permeability and magnetic permittivity in vacuum. The PEC object is denoted by and it is enclosed by a smooth boundary , also denotes the propagation medium. Frequently, the incident field is a plane wave given by , where
is a non-zero vector representing the polarisation of the wave andis a unit vector perpendicular to that gives the direction of the plane wave.
The simulation of time-harmonic waves scattered from perfect electric conductors (PECs) is a key application area for boundary element methods as they naturally deal with the infinite domain problem by reducing the computation of the numerical solution of an integral equation on a finite dimensional surface. In the low-frequency regime a direct solution of the electric field integral equation (EFIE) is well established (see [Buffa2003]). As the mesh width decreases, however, the formulation suffers from ill-conditioning, hence a regulariser, is needed:
Here, a direct formulation of the EFIE has been presented, with as the electric field integral operator and as the magnetic field integral operator. In the most simple form is just an inverse diagonal or an algebraic preconditioner, but more sophisticated implementations are based on the idea of choosing as the discretisation of an operator that regularises the electric field operator in continuous spaces [Bagci2009].
The most common example of is the Calderón Multiplicative Preconditioner [andriulli2008multiplicative] and difficulty here is the evaluation of operator products. Assume a function in some Hilbert space, and operators and in compatible Hilbert spaces. In order to evaluate the product through Galerkin discretisations of the operators and we need to compute a finite dimensional matrix product of the form where we have used rhe subscript to denote finite dimensional quantities after discretisation [Betcke20]. The matrix is a mass matrix, which contains the inner product of the test space of and the domain space of . The difficulty is that this mass matrix is numerically singular for the standard choice of RWG spaces in electromagnetic scattering [buffa2007dual], a problem that does not occure for simple P1 spaces used in many acoustic problems. In order to overcome this problem one can use so-called Buffa-Christiansen (BC) spaces as test spaces [andriulli2008multiplicative]. However, their use is expensive due to the requirement of barycentric refinements. One aim of this research is to build a preconditioner that avoids the use of BC spaces and that can simply be defined using RWG spaces in the primal grid.
Aiming to avoid the use of BC spaces, the optimal operator to use as a regulariser has been considered to be the Magnetic-to-Electric (MtE) operator. It turns out that with this choice the left-hand side of (2) becomes a second kind integral operator, thus in theory, making it much more amenable to iterative solvers. However, a direct evaluation of the MtE operator has a similar complexity to solving the original scattering problem. It is therefore not practical. In the early 2000s Antoine, Darbas, and their co-authors investigated in the acoustic scattering case a class of approximations to the Dirichlet-to-Neumann (DtN) map based on Padé expansions of high-frequency approximations of this operator. It turns out that this OSRC method (On-Surface-Radiation Condition) allows an accurate approximation of the DtN map and its inverse simply by solving sparse surface linear systems of equations on the boundary of the scatterer (see [antoine2006improved, antoine2007generalized, darbas13], and [kriegsmann1987new] for earlier work by Kiregsmann, Taflove and Umashankar).
In recent work by Bouajaji, Antoine, and Geuzaine this OSRC concept was extended to the MtE operator in three dimensional electromagnetic scattering, allowing it to be well approximated by the solution of sparse linear systems arising from surface differential operators [el2014approximate].
In this paper we follow this approach for approximating the MtE operator and investigate the suitability of this approximation as regulariser for the EFIE. We also demonstrate details of the implementation and the resulting preconditioning performance with respect to time and matrix-vector products.
The paper is organised as follows: in section 2 an overview of physical problem and function spaces for Maxwell boundary integral operators is presented. Section 3 shows the MtE preconditioner and its continuous implementation and approximations, while section 4 shows the discrete implementation of MtE as a preconditioner. Section 5 presents a validation of this new preconditioner and some tests on its performance, with conclusion sgiven in section 6.
2 Problem setting, Tangential Sobolev Spaces and Surface Operators
2.1 Function Spaces
To solve the problem in (19), consider the definition of generic Sobolev spaces on :
Here is a pseudodifferential operator, the sub-fix stands for the assumption of integrability within all compact subsets in (which is not necessary when is bounded. In that case, is replaced by . However, in order to formulate boundary value problems for Maxwell’s equations, it is necessary to define spaces on the boundary, hence the tangential, magnetic and Neuman traces are given respectively by:
for , , and .
As stated in [Buffa2003], Green’s formula implies the definition of the tangential trace on :
with , . Since is surjective, then [Buffa2003] defines the following space for smooth surfaces:
and from (3) it is natural to define as the dual space of . Hence , is well-defined and continuous. Also, the scalar surface divergence of , for is defined as . Therefore, on a smooth domain we have: . Since this is not obvious for non-smooth domains [Buffa2003], in order to build a the Galerkin formulation, it is assumed that the grid over consists of a finite union of smooth faces . Hence, according to [Buffa2003, definition 1] is replaced by , which preserves the definition in (4): (here defines the tangential space of square integrable functions).
Notice that for , considering the dual pairing in (3), it is possible to define as its self dual with respect to as a pivot space. For the same reason becomes the dual of respect to the standard dual product. To complete the description of the properties of the tangential and magnetic traces, it is necessary to remark that and are continuous and subjective ( [Buffa2003, Lemma 3] ) on finite domains.
To solve (1), the Stratton-Chu representation formulas [monk2003finite] must be considered for any :
where are defined as:
With , .
By applying magnetic and tangential traces to the Electric and Magnetic field potential operators the electric and magnetic BIOs can be obtained:
Then, applying electric and magnetic traces to the representation formula (5), the following can be derived:
This is the so called Calderon projector, which allows to find the electric and magnetic fields in the space, from a boundary condition, or . This Projector has some important properties, for instance: from which the following can be deduced:
This is the basis for Calderón Preconditioning. It basically states that equals the identity plus a compact operator (a second kind Fredholm operator [kress1989linear]). Since the spectrum of compact operators clusters at 0, the addition of the identity and this compact operator results in a third operator that in terms of spectral properties will remain stable when discretised. The difficulties of building the discrete version operator arise when trying to build gram matrices to implement the discrete product . These matrices are unstable due to the nature of the discrete spaces used in their construction, so to overcome this problem a CMP has been proposed in [andriulli2008multiplicative] and has become standard when building a CP. However, as it was pointed out in the introduction, this technique requires the use of BC basis functions, which translates into an increase of computational complexity of the algorithm, hence the intention to build alternative preconditioners.
3 Construction of an OSRC Preconditioner
In this section we show that the exact Magnetic to Electric (MtE) operator and its inverse (EtM) constitute perfect preconditioners and with this, establish the basic idea of the behavior of an approximate OSRC operator as a preconditioner.
3.1 OSRC operator as a preconditioner for the EFIE
Starting with the exact EtM (Electric to Magnetic) operator, this can be derived from the first row of the Calderon Projector (6):
hence the EtM and its inverse (MtE, Magnetic to Electric) operator are given by:
And from the second row :
a second version of the EtM exact operator and its inverse (MtE) is obtained:
An important fact to bear in mind when deriving these expressions is that it is assumed that either or are invertible operators when required.
It is necessary to remark that the discrete versions of the pairs (8), (11) and (9), (10) are not necessarily the same, since they might be defined on different discrete spaces. For example, if the electric data is represented in an RWG space and the magnetic data in a BC space [scroggs2017software] possible domain and range pairing of (8), (11) are (RWG, BC) and (BC, RWG) respectively and for (9), (10) are (BC, RWG) and (RWG, BC), also respectively. 111Here RWG stands for Rao Wilton Glisson basis functions. To discriminate these versions of the operator we have used the subscripts and and we state the following:
Assuming that is smooth and that is not an interior eigenvalue, the following holds for the regularised EFIE:
is not an interior eigenvalue, the following holds for the regularised EFIE:
where is the inverse MtE operator. In other words, is a preconditioner.
The aplication of the inverse MtE operator to the EFIE operator by the right leads directly to an operator with better spectral properties:
On the other hand, the inverse MtE applied from the left gives:
From (7) we have that . Then:
As a conclusion, and are equivalent to second kind Fredholm operators.
From now on, we drop the subscripts and we refer to just as .
From the last theorem, the MtE operator would be an ideal alternative as a preconditioner, :
However, the construction of is as expensive as solving the EFIE, so finding a good approximation is essential.
3.2 Approximation of the MtE operator ()
[el2014approximate] shows that an approximation for the EtM on smooth surfaces is given by
where is given by
Hence . However, in order to implement the MtE operator as a preconditioner for the CFIE as in (12), it is required to modify the mapping properites of the MtE operator in such a way such that it maps from to . Therefore, it is necessary to define two auxiliary maps,
can be used as a preconditioner. Such maps can be defined conveniently, depending on the discrete spaces used as it will be shown in next sections.
From the definition of (14), the inverse of this operator can also be found: .
3.3 Padé approximation of the MtE operator ()
Regarding the implementation of , is a local symmetrical operator that can be easily approximated. However, is a non-local pseudodifferential operator whose calculation is not given straightforward. To overcome this problem, a local approximation of can be obtained from a rotating branch cut Padé Approximation:
Renaming as , the Padé approximation of can be applied on by solving the following ():
Now, by introducing , the original weak formulation for reads:
Which in conjunction with (16), it can be rewritten as a system of equations:
The calculation of , follows from composing the weak form of with the weak form of :
As an example of how similar is to as , figure 1 shows how the spectrum of converges to as grows.
4 Implementation of the discrete MtE preconditioner
To build the discrete approximation of , consider a polyhedral approximation of with a triangulation . The following discrete spaces are introduced: as the usual Nédeléc’s space of linear edge finite element defined in [monk2003finite], as the usual Rao Wilton Glisson’s space of linear edge finite element defined in [rao1982electromagnetic] and (space of scalar linear basis functions).
Resorting to the same notation used in [el2014approximate], some discrete operators are needed:
where and belong to , belongs to and , belong to .
From (17) the discrete formulation of () can be written in matricial form as:
Therefore, a block matrix inversion can be performed to get:
Which gives a first insight for the implementation of .
Finally, the discrete preconditioned formulation is given by:
We must be careful with the operator products in the discretisation. Naively, we would introduce Gram matrices to evaluate the product between and . However, due to the nature of RWG spaces, some of these are singular, but luckily, it turns that they are not needed in practice. Notice that when defining one must pay attention to the fact that maps from the Nédélec space, , to SNC. Since the test space of is also SNC, then becomes just a standard identity, .
on the other hand is a map from the SNC space to the RWG space. Having in mind that then is just .
Finally, we need to map back to the trial space of space of (also the test space of ) so the discrete formulation remains coherent. Thus, we perform the product and the discrete implementation reduces to:
There are 2 ways to solve the discrete problem:
Proceed as in [el2014approximate] and build a block sparse system. Unfortunately, the system generated by this method also had poor spectral properties, so it was discarded as an option.
Implement a GMRES iteration using the following algorithm as the matrix-vector product:
(the presence of mass matrices has been omitted here). Notice that in the previous algorithm, all the matrices to be inverted are sparse, except for , in the third step. Recalling (19):
There are different approaches to complete the construction of the GMRES iteration:
The first method consists of inverting using either expression (1) or (2). Notice that solving (1) needs the inversion of a sparse matrix of dimension ( stands for number of edges and for number of vertices on the grid), but (2) requires the inversion of a dense matrix of dimension , which is as expensive as solving the non regularised EFIE. Hence in 2 a LU decomposition of (1) is passed as an input and a forward-backward substitution is performed in step 3.
The second method consists of implementing (2), based in the following heuristic, inspired in figure2 which illustrated the general behavior of Padé coefficients for any .
Notice that when is small, then the value of becomes big and in this case . On the other hand, when is big, becomes small and then can be neglected.
Hence, we can approximate
where is a subset of indexes in such that . Therefore, the calculation of (2) reduces to a partial sum of Padé coefficients and the calculation of .
Finally, the following summarizes the notation used for each kind of preconditioner:
: MtE approximate operator.
: Padé MtE approximate operator of first type, of order .
: Padé MtE approximate operator of second type, of order .
5 Numerical Experiments
This section demonstrates the performance of the proposed OSRC preconditioner and some comparisons with other existing preconditioners. All the tests were performed using BEM++ software [scroggs2017software] on a spherical grid of radius .
The OSRC operator was validated using the same method as in [el2014approximate], where bistatic RCS calculated using:
A direct formulation of the EFIE ().
An approximation of calculated by computing the square roots of the eigenvalues of the matrix that generates .
By applying .
In order to validate the construction of the calculations of numerical solutions using this operator were compared to analytical solutions on the unit sphere calculated using spherical harmonics.
The first test (figure 3) replicates the results obtained in [el2014approximate] and shows the bistatic RCS obtained from the scattering problem of an incident electromagnetic plane wave by a PEC unit sphere, for and . In both cases, the solution calculated using is a good approximation of the analytic solution.
5.2 Performance Comparison
In order to compare the performance of the OSRC preconditioner to existing preconditioners, the following attributes were benchmarked:
GMRES number of iterations.
Relative errors between numerical and analytical solutions
Assembly and solving times for:
Pure direct formulation of the EFIE, denoted by (always calculated in the primal grid).
EFIE preconditioned using the Calderon Multiplicative Preconditioner (CMP), denoted by .
EFIE preconditioned with , with .
EFIE preconditioned with , with . In this case, for we have just one Padé coefficient. When we also only keep just one Padé coefficient because the second coefficient in the expansion is sufficiently small in size to be discarded. We note that this is not identical to the case , since the definition of Padé coefficients in [el2014approximate] shows that these differ when the degree of the expansion () changes.
These tests were performed using -matrix compressions, with
Figure 4 shows performs better than as the number of iterations that takes to solve the system grows at a lower rate. As expected, shows (slightly) better results than , because the first is a better approximation than the latter.
Tables 1 shows that is more efficient in terms of assembly times than the CMP, since only requires to assemble s