Add OSRC Magnetic-to-Electric and Electric-to-Magnetic operators

Project.toml

@@ -25,6 +25,7 @@ LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 LinearMaps = "7a12625a-238d-50fd-b39a-03d52299707e"
 NestedUnitRanges = "032820ab-dc03-4b49-91f4-7d58d4da98b3"
 OhMyThreads = "67456a42-1dca-4109-a031-0a68de7e3ad5"
+Polynomials = "f27b6e38-b328-58d1-80ce-0feddd5e7a45"
 ProgressMeter = "92933f4c-e287-5a05-a399-4b506db050ca"
 Requires = "ae029012-a4dd-5104-9daa-d747884805df"
 SauterSchwab3D = "0a13313b-1c00-422e-8263-562364ed9544"
@@ -63,8 +64,9 @@ LiftedMaps = "0.5.1"
 LinearMaps = "3.7 - 3.9, 3.11.2"
 NestedUnitRanges = "0.2.3"
 OhMyThreads = "0.8.3"
-ProgressMeter = "1.11.0"
 PlotlyJS = "0.18.18"
+Polynomials = "4.1.1"
+ProgressMeter = "1.11.0"
 Requires = "1"
 SauterSchwab3D = "0.2"
 SauterSchwabQuadrature = "2.4.0"
@@ -74,5 +76,3 @@ Suppressor = "0.2.8"
 TestItems = "0.1.1, 1"
 WiltonInts84 = "0.2.8"
 julia = "1.10"
-
-

docs/make.jl

@@ -72,6 +72,9 @@ makedocs(;
                 "Low Frequency Stable PMCHWT"=>"projectors/pmchwt_lf.md",
                 "Dense Grid Stable PMCHWT"=>"projectors/pmchwt_theta.md",
             ],
+            "OSRC operators"=>Any[
+                "OSRC operators"=>"operators/OSRC_efie.md",
+            ],
         ],
         "Basis Functions" => Any[
             "Overview"=>"bases/overview.md",

docs/runliterate.jl

@@ -3,4 +3,5 @@ using Literate
 Literate.markdown("./examples/composedoperator.jl","./docs/src/composedoperator"; execute=true)
 Literate.markdown("./examples/pmchwt_theta.jl","./docs/src/projectors"; execute=true)
 Literate.markdown("./examples/pmchwt_lf.jl","./docs/src/projectors"; execute=true)
+Literate.markdown("./examples/OSRC_efie.jl","./docs/src/operators"; execute=true)
 

examples/OSRC_efie.jl

@@ -0,0 +1,128 @@
+# In this example, the preconditioning ability of the On-Surface Radiation Condition (OSRC) Magnetic-to-Electric (MtE) map is shown on the Electric Field Integral Equation (EFIE), on the sphere.
+# # OSRC operators 
+# The MtE map ``\text{V}^{+}`` is a map of the trace of the magnetic field ``n \times h`` to the trace of the electric field ``n \times e`` on a scattering surface.
+# The approximation of the local MtE surface operator as an On-Surface Radiation Condition (OSRC) operator on an arbitrary surface ``\Gamma`` is given by:
+# ```math
+# \begin{equation}
+#     \text{V}^{+} \left( n \times h \right) + n \times e = 0, \quad \text{on} \ \Gamma,
+# \end{equation}
+# ```
+# In the implementation a highly accurate pseudo-differential operator ``\text{V}_{\epsilon}^{+}`` is used to represent the Magnetic-to-Electric map.
+# ```math
+# \begin{aligned}
+#     \text{V}_{\epsilon}^{+}  &:= - \pmb{\Theta}^{-1} \Lambda_{2, \epsilon}^{-1} \Lambda_{1, \epsilon}, \\
+#     \pmb{\Theta} \left( \cdot \right) &:= \left( n \times \cdot \right), \\
+#     \Lambda_{2, \epsilon} &:= I - \textbf{curl}_{\Gamma} \frac{1}{\kappa_{\epsilon}^2} \textup{curl}_{\Gamma}, \\
+#     \Lambda_{1, \epsilon} &:= \left(I + \mathcal{J} \right)^{1/2}, \\
+#     \mathcal{J} &:=\frac{\Delta_{\Gamma}}{\kappa_{\epsilon}^2} = \textbf{Grad}_{\Gamma} \frac{1}{\kappa_{\epsilon}^2} \textup{Div}_{\Gamma} - \textbf{curl}_{\Gamma} \frac{1}{\kappa_{\epsilon}^2} \textup{curl}_{\Gamma}. 
+# \end{aligned}
+# ```
+# More details on this OSRC operator and its implementation are given in
+# [An OSRC Preconditioner for the EFIE (Betcke et al. (2021))](https://arxiv.org/abs/2111.10761)
+# The Electric-to-Magnetic (EtM) OSRC operator is also implemented, as given in 
+# [Approximate local magnetic-to-electric surface operators for time-harmonic Maxwell’s equations (C. Geuzaine et al. (2014))](https://www.researchgate.net/publication/261636307_Approximate_local_magnetic-to-electric_surface_operators_for_time-harmonic_Maxwell's_equations)
+#
+# # OSRC EFIE preconditioning example
+# We begin with loading needed packages
+using BEAST, CompScienceMeshes
+
+using Makeitso
+using LinearAlgebra
+
+using Plots
+using PlotlyDocumenter
+
+# Now, assemble and solve the unpreconditioned EFIE
+
+BLAS.set_num_threads(8)
+
+@target geo (;h) -> begin
+      (; Γ = CompScienceMeshes.meshsphere(radius=1.0, h=h))
+end
+
+@target formulation_EFIE (geo,; κ) -> begin
+    X = raviartthomas(geo.Γ)
+    T = Maxwell3D.singlelayer(wavenumber=κ)
+
+    E = Maxwell3D.planewave(direction=ẑ, polarization=x̂, wavenumber=κ);
+    e = (n × E) × n;
+
+    bx = assemble(e, X)
+    A = assemble(T,X,X); 
+    return (;bilforms=(;A), linforms=(;bx))
+end
+
+@target solution_EFIE (formulation_EFIE,; residual) -> begin
+      (;bilforms, linforms) = formulation_EFIE
+      (;A) = bilforms; (;bx) = linforms;
+      iT = BEAST.GMRESSolver(A; restart=1_500, reltol=residual, maxiter=1_500)
+      u, ch = BEAST.solve(iT, bx)
+      return (;iters=ch.iters, u=u)
+end
+
+# Next, assemble the OSRC MtE operator on the primal grid and use it to precondition the EFIE
+
+@target OSRC_MtE_op (geo,;κ, Np, curvature) -> begin
+    MtE_OSRC_operator = BEAST.MtE_OSRC_op(κ, Np, pi/2, curvature)
+    Nd = BEAST.nedelec(geo.Γ)
+    MtE_map = assemble(MtE_OSRC_operator, Nd, Nd)
+    return (;MtE=MtE_map)
+end
+
+@target solution_OSRC_EFIE (formulation_EFIE, OSRC_MtE_op; residual) -> begin
+      (;bilforms, linforms) = formulation_EFIE
+      (;A) = bilforms; (;bx) = linforms;
+      P_OSRC = OSRC_MtE_op.MtE
+      iT = BEAST.GMRESSolver(A; restart=1_500, reltol=residual, maxiter=1_500, left_preconditioner=P_OSRC)
+      u, ch = BEAST.solve(iT, bx)
+      return (;iters=ch.iters, u=u)
+end
+
+# Finally, assemble and use the Calderón preconditioner, for a comparison with OSRC preconditioning
+
+@target calderon_preconditioner (geo,;κ) -> begin
+    Γ = geo.Γ
+    X = raviartthomas(Γ)
+    Y = BEAST.buffachristiansen(Γ)
+
+    T = Maxwell3D.singlelayer(wavenumber=κ)
+    N = NCross()
+
+    Tyy = assemble(T,Y,Y);
+    Nxy = Matrix(assemble(N,X,Y));
+    iNxy = inv(Nxy);
+    P = iNxy' * Tyy * iNxy
+    return (;P=P)
+end
+
+@target solution_cald_EFIE (formulation_EFIE, calderon_preconditioner; residual) -> begin
+      (;bilforms, linforms) = formulation_EFIE
+      (;A) = bilforms; (;bx) = linforms;
+      P_calderon = calderon_preconditioner.P
+      iT = BEAST.GMRESSolver(A; restart=1_500, reltol=residual, maxiter=1_500, left_preconditioner=P_calderon)
+      u, ch = BEAST.solve(iT, bx)
+      return (;iters=ch.iters, u=u)
+end
+
+# Set the simulation parameters and run for different mesh sizes ``h``
+
+@target benchmark_OSRC (solution_EFIE, solution_OSRC_EFIE, solution_cald_EFIE, ;) -> begin
+      return (;iters_EFIE = solution_EFIE.iters, iters_OSRC_EFIE = solution_OSRC_EFIE.iters, iters_cald_EFIE = solution_cald_EFIE.iters)
+end
+@sweep benchmark_OSRC_sweep (!benchmark_OSRC,; h=[], κ=[], residual=[], Np=[], curvature=[]) -> benchmark_OSRC
+
+h_values = [0.3, 0.2, 0.15]
+κ = 1.0*pi
+residual = 1e-6
+Np = 4
+curvature = 1.0
+
+df = make(benchmark_OSRC_sweep; h=h_values, κ=κ, residual=residual, Np=Np, curvature=curvature)
+
+# Finally, the GMRES iterations for the different methods are plotted
+using Plots 
+plot(df.h, df.iters_EFIE, label="EFIE", l=2)
+plot!(df.h, df.iters_OSRC_EFIE, label="OSRC-EFIE", l=2)
+plot!(df.h, df.iters_cald_EFIE, label="Calderón-EFIE", l=2)
+xlabel!("h")
+ylabel!("GMRES iterations")

src/BEAST.jl

@@ -243,6 +243,7 @@ include("dyadicop.jl")
 include("interpolation.jl")
 include("operators/projectors.jl")
 include("operators/theta.jl")
+include("operators/OSRC.jl")
 
 include("quadrature/rules/momintegrals.jl")
 include("quadrature/doublenumints.jl")