The Time Domain Linear Sampling Method for Determining the Shape of Multiple Scatterers Using Electromagnetic Waves

A The Discontinuous Galerkin Forward Solver

Recalling our assumption that c 0 = 1 , Maxwell’s equations for a regularized electric dipole at 𝐱 i with polarization 𝒑 can be written in the conservation form as (see [21, Section 10.5])

Here (A.1), Q is a block matrix given by

and

where 𝐞 ℓ is the ℓ -th Cartesian unit vector. In addition, the right-hand side of (A.1) is written

The coefficient β = β ⁢ ( 𝐱 ) is used to add additional damping in the zone next to the absorbing condition and δ 𝐱 i is the Dirac delta function at 𝐱 i . We assume that the materials are isotropic and piecewise homogeneous, so the coefficients corresponding to the relative electric permittivity ϵ r and magnetic permeability μ r are identified piecewise with real scalars.

The spatial derivatives in (A.1) are discretized using the nodal discontinuous Galerkin method [21], while the time integration is done by the low-storage explicit Runge–Kutta method [7]. In the discretized version, we assume that the computational domain Ω ~ ⊂ ℝ 3 is divided into N K tetrahedral elements, Ω ~ = ⋃ k = 1 N K D k . The boundary of element D k is denoted by Γ k . We assume that the elements are aligned with material discontinuities. Furthermore, for any element D k the superscript “ - ” refers to interior information while “ + ” refers to exterior information.

We multiply (A.1) by a local test function ϕ k and integrate by parts twice to obtain an elementwise variational formulation

𝐪 k is the restriction of 𝐪 to the element D k and ℱ * is the numerical flux across neighboring element interfaces. For the numerical flux ℱ * along the normal 𝝂 , we use the upwind [20]

where

In this work, we apply impedance and perfect electric conductor (PEC) boundary conditions. On the exterior boundary, the PEC condition is recovered from (A.2) and (A.3) by setting ϵ r + = ϵ r - , μ r + = μ r - ,

The impedance boundary condition is obtained from (A.2) and (A.3) by setting ϵ r + = ϵ r - = ϵ r bc , μ r + = μ r - = μ r bc , and

The impedance boundary condition reduces to the Silver–Müller absorbing (SMA) by setting parameters ϵ r bc = ϵ r - and μ r bc = μ r - (i.e. the physical values of the interior element). Unfortunately, the SMA condition is not perfect and some unwanted reflections will happen at the outflow boundaries if the incoming wave is not parallel with the boundary. In this paper, we couple the absorbing boundary condition with a sponge layer (SL) that damps the wave. To do so, the variable β introduced in (A.1) is non-zero for the regions next to the SMA condition. Moreover, the SL is coupled with the grid stretching.

To illustrate the functioning of the SL, let us consider an example in which the layer is applied on one coordinate axis only. Now, for a node coordinate x 1 ( ℓ ) ∈ [ x 1 ( 0 ) , x 1 ( 0 ) + L ] , where x 1 ( 0 ) denotes a starting location of the SL and L its thickness, the grid stretched coordinate x ^ 1 ( ℓ ) is defined as

where g max denotes the maximum value given for the grid stretching. Similarly, the β value at x ^ 1 ( ℓ ) in the SL is

where β max is the maximum value given for the damping coefficient.

The current version of the wave solver is written in the C/C++ programming language and is integrated with the Open Concurrent Compute Abstraction (OCCA) [28] library and message passing interface to enable parallel computations both on CPU and GPU clusters. Currently, the solver uses only constant order basis functions and the computational load between different elements is balanced by the parMetis software [22].

Due to the assumption of using a magnetic dipole as a source, and the fact that the current DG-based wave solver assumes an electric dipole, we use the magnetic field ℋ as data for the inverse solver. Because of the constant coefficients in Maxwell’s equations, this corresponds to the electric field due to a magnetic dipole.