Code Verification

This package is well suited to verify the correctness of more involved numerical techniques to determine the scattering behavior of real-world objects, such as boundary element methods. Two examples are given below, both employing the packages BEAST and CompScienceMeshes: a PEC sphere and a sound-soft and a sound-hard sphere, each excited by a plane wave.

Electromagnetic: PEC Sphere

The scattering from PEC objects can be determined by solving the electric field integral equation (EFIE) by the method of moments (MoM):

using BEAST, CompScienceMeshes

# --- parameters
f = 1e8             # frequency
c = 2.99792458e8    # speed of light
μ = 4π * 1e-7       # permeability
κ = 2π * f / c      # wavenumber

# --- obtain a triangulation of the sphere
spRadius = 1.0 # radius of sphere
Γ = meshsphere(spRadius, 0.25)

# --- define basis functions on the triangulation
RT = raviartthomas(Γ)

# --- excitation by plane wave
𝐸 = Maxwell3D.planewave(; direction=ẑ, polarization=x̂, wavenumber=κ)
𝑒 = n × 𝐸 × n

# --- integral operator
𝑇 = Maxwell3D.singlelayer(; wavenumber=κ)

# --- assemble matrix and RHS of the LSE
e = -assemble(𝑒, RT)
T = assemble(𝑇, RT, RT)

# --- solve the LSE
u = T \ e

942-element Vector{ComplexF64}:
   -0.0410797762571886 - 0.27608593748490046im
   0.06928485262595385 + 0.2770897702960765im
   0.09785418567547545 + 0.21879174276068888im
 -0.002678944507175422 - 0.01512788714999675im
   0.03530860897605319 - 0.43197120769481995im
  -0.02575857884534293 + 0.4305875884769166im
 0.0018068627371974743 + 0.0324134118695019im
   0.08414041525797701 - 0.41141589712303483im
  0.007442032775828352 + 0.046321363033134874im
    0.0648896313437135 + 0.20307900953006489im
                       ⋮
  0.014169422118546453 + 0.034253543437254265im
   0.03698742671590207 - 0.1600123353247938im
   0.00950059219655661 - 0.10529002384030085im
  0.049467495796120343 + 0.07935593345991158im
   0.13395756098447714 + 0.08279348256040814im
  -0.05905764310030015 - 0.03248677609805959im
    0.0802761586240356 - 0.171956652941191im
  -0.06203225036121507 + 0.0668622731330171im
  0.031799655743877234 + 0.07451277275211148im

Now that we know the expansion coefficients for the basis functions we can compute the scattered fields:

# --- define points where the fields are computed
using SphericalScattering # use this package to get points on a spherical grid
points_cartFF, points_sphFF = sphericalGridPoints()
points_cartNF, points_sphNF = sphericalGridPoints(r=5.0)

# --- compute the fields
EF_MoM = potential(MWSingleLayerField3D(𝑇), points_cartNF, u, RT)
HF_MoM = potential(BEAST.MWDoubleLayerField3D(wavenumber=κ), points_cartNF, u, RT) / (c * μ)
FF_MoM = -im * f / (2 * c) * potential(MWFarField3D(𝑇), points_cartFF, u, RT)

The fields for the same scenario can be obtained (more accurately) by this package:

using SphericalScattering

sp = PECSphere(radius=spRadius)
ex = planeWave(frequency=f)

EF = scatteredfield(sp, ex, ElectricField(points_cartNF))
HF = scatteredfield(sp, ex, MagneticField(points_cartNF))
FF = scatteredfield(sp, ex, FarField(points_cartFF))

The agreement between the two solutions can be determined as a worst case relative error of all evaluated points:

using LinearAlgebra

# --- relative worst case errors in percent
diff_EF = round(maximum(norm.(EF - EF_MoM) ./ maximum(norm.(EF))) * 100, digits=2)
diff_HF = round(maximum(norm.(HF - HF_MoM) ./ maximum(norm.(HF))) * 100, digits=2)
diff_FF = round(maximum(norm.(FF - FF_MoM) ./ maximum(norm.(FF))) * 100, digits=2)

print("E-field error: $diff_EF %\n")
print("H-field error: $diff_HF %\n")
print("Far-field error: $diff_FF %\n")

E-field error: 1.29 %
H-field error: 1.28 %
Far-field error: 1.37 %

Acoustic: Sound-Soft and Sound-Hard Sphere

The same approach applies to the acoustic scatterers. Both packages employ the time convention $\mathrm{e}^{\mathrm{j}\omega t}$, so that the incident plane wave $u^\mathrm{i} = \mathrm{e}^{-\mathrm{j} \kappa \hat{\bm e}_z \cdot \bm r}$ and the Green's function $\mathrm{e}^{-\mathrm{j} \kappa R} / (4 \pi R)$ agree. As in the electromagnetic example, the sphere is embedded in the default medium Medium(ε0, μ0), so that the wavenumber follows from the frequency as $\kappa = 2 \pi f / c$ with $c = 1 / \sqrt{\varepsilon_0 \mu_0}$:

using BEAST, CompScienceMeshes
using SphericalScattering

# --- parameters
f = 7e7                 # frequency
c = 1 / sqrt(ε0 * μ0)   # wave speed in the default medium Medium(ε0, μ0)
κ = 2π * f / c          # wavenumber

# --- obtain a triangulation of the sphere
spRadius = 1.0 # radius of sphere
Γ = meshsphere(spRadius, 0.15)

# --- excitation by a plane wave travelling in z-direction
uⁱ = Helmholtz3D.planewave(; direction=ẑ, wavenumber=κ) # BEAST
ex = Acoustic.planeWave(; frequency=f)                  # this package

# --- points where the scattered fields are compared
points_cartNF, points_sphNF = sphericalGridPoints(r=5.0)

On a sound-soft sphere the total pressure $u = u^\mathrm{i} + u^\mathrm{s}$ vanishes. Its Neumann trace $\sigma = \partial_n u$ solves the first-kind equation $\mathcal{S} \sigma = u^\mathrm{i}$ with the single-layer operator $\mathcal{S}$, and it radiates the scattered field $u^\mathrm{s} = -\mathcal{S}[\sigma]$:

# --- sound-soft: piecewise constant Neumann trace
X = lagrangecxd0(Γ)
S = assemble(Helmholtz3D.singlelayer(; wavenumber=κ), X, X)
σ_MoM = S \ assemble(strace(uⁱ, Γ), X)
us_soft_MoM = -potential(HH3DSingleLayerNear(; wavenumber=κ), points_cartNF, σ_MoM, X; type=ComplexF64)

# --- the same quantities by this package, the trace at the centroids of the triangles
sp = SoftSphere(radius=spRadius)
σ = field(sp, ex, PressureNormalGradient(X.pos))
us_soft = scatteredfield(sp, ex, Pressure(points_cartNF))

On a sound-hard sphere the normal derivative of the total pressure vanishes. Its Dirichlet trace $u$ solves the first-kind equation $\mathcal{W} u = \partial_n u^\mathrm{i}$ with the hypersingular operator $\mathcal{W}$, and it radiates the scattered field $u^\mathrm{s} = \mathcal{D}[u]$ by the double-layer potential:

# --- sound-hard: piecewise linear Dirichlet trace
Y = lagrangec0d1(Γ)
W = assemble(Helmholtz3D.hypersingular(; wavenumber=κ), Y, Y)
u_MoM = W \ assemble(∂n(uⁱ), Y)
us_hard_MoM = potential(HH3DDoubleLayerNear(; wavenumber=κ), points_cartNF, u_MoM, Y; type=ComplexF64)

# --- the same quantities by this package, the trace at the vertices
sp = HardSphere(radius=spRadius)
u = field(sp, ex, PressureTrace(Y.pos))
us_hard = scatteredfield(sp, ex, Pressure(points_cartNF))

The centroids and the vertices of the mesh do not lie exactly on the sphere; as the traces depend only on the direction of a location, they are evaluated on the sphere nevertheless, see Quantities. The agreement is again determined as a worst case relative error of all evaluated points:

# --- relative worst case errors in percent
relerr(a, b) = round(maximum(abs.(a - b)) / maximum(abs.(b)) * 100, digits=2)

print("sound-soft, Neumann trace error:   $(relerr(σ_MoM, σ)) %\n")
print("sound-soft, scattered field error: $(relerr(us_soft_MoM, us_soft)) %\n")
print("sound-hard, Dirichlet trace error: $(relerr(u_MoM, u)) %\n")
print("sound-hard, scattered field error: $(relerr(us_hard_MoM, us_hard)) %\n")

sound-soft, Neumann trace error:   3.54 %
sound-soft, scattered field error: 0.33 %
sound-hard, Dirichlet trace error: 0.47 %
sound-hard, scattered field error: 0.62 %