From 3d3f09550ec20de17a3e2e0cce833268c04b5ed8 Mon Sep 17 00:00:00 2001
From: Kristof Cools <krcools@gmail.com>
Date: Wed, 23 Oct 2024 13:09:07 +0200
Subject: [PATCH] efie convergence examples

---
 examples/efie_convergence.jl        | 55 +++++++++++++++++++++++++
 examples/efie_convergence_cuboid.jl | 63 +++++++++++++++++++++++++++++
 2 files changed, 118 insertions(+)
 create mode 100644 examples/efie_convergence.jl
 create mode 100644 examples/efie_convergence_cuboid.jl

diff --git a/examples/efie_convergence.jl b/examples/efie_convergence.jl
new file mode 100644
index 00000000..ff22ece0
--- /dev/null
+++ b/examples/efie_convergence.jl
@@ -0,0 +1,55 @@
+using CompScienceMeshes
+using BEAST
+
+using BakerStreet
+using DataFrames
+
+function payload(;h)
+    # Γ = readmesh(joinpath(dirname(pathof(BEAST)),"../examples/sphere2.in"))
+    @show h
+    Γ = meshsphere(radius=1.0, h=h)
+    # @show length(Γ)
+    X = raviartthomas(Γ)
+
+    κ, η = 1.0, 1.0
+    t = Maxwell3D.singlelayer(wavenumber=κ)
+    s = -Maxwell3D.singlelayer(gamma=κ)
+    E = Maxwell3D.planewave(direction=ẑ, polarization=x̂, wavenumber=κ)
+    e = (n × E) × n
+
+    @hilbertspace j
+    @hilbertspace k
+
+    A = assemble(t[k,j], X, X)
+    S = assemble(s[k,j], X, X)
+    b = assemble(e[k], X)
+
+    Ai = BEAST.GMRESSolver(A; restart=1500, abstol=1e-5, maxiter=10_000)
+    u = Ai * b
+    u_Snorm = real(sqrt(dot(u, S*u)))
+    return (;u, X, u_Snorm)
+end
+
+α = 0.8
+h = collect(0.4 * α.^(0:10))  
+# h = collect(logrange(0.4, 0.05, length=8))
+@runsims payload h
+
+df = @collect_results
+
+using LinearAlgebra
+using Plots
+plot(df.h, real.(df.u_Snorm), marker=:.)
+
+# α = df.h[1] / df.h[2]
+W1 = reverse(df.u_Snorm)
+W2 = Iterators.drop(W1,1)
+W3 = Iterators.drop(W1,2)
+p = [log((w2-w1)/(w3-w2)) / log(1/α) for (w1,w2,w3) in zip(W1,W2,W3)]
+
+# efie = @discretise t[k,j]==e[k]  j∈X k∈X
+# u, ch = BEAST.gmres_ch(efie; restart=1500)
+
+# include("utils/postproc.jl")
+# include("utils/plotresults.jl")
+
diff --git a/examples/efie_convergence_cuboid.jl b/examples/efie_convergence_cuboid.jl
new file mode 100644
index 00000000..8f4fa5b7
--- /dev/null
+++ b/examples/efie_convergence_cuboid.jl
@@ -0,0 +1,63 @@
+using CompScienceMeshes
+using BEAST
+
+using BakerStreet
+using DataFrames
+
+function payload(;h)
+    # Γ = readmesh(joinpath(dirname(pathof(BEAST)),"../examples/sphere2.in"))
+    # @show h
+    d = 1.0
+    Γ = meshcuboid(d, d, d/2, h)
+    d = 1.2
+    Γ′ = meshcuboid(d, d, d/2, h)
+    # @show length(Γ)
+    X = raviartthomas(Γ)
+    X'
+
+    κ, η = 1.0, 1.0
+    t = Maxwell3D.singlelayer(wavenumber=κ)
+    s = -Maxwell3D.singlelayer(gamma=κ)
+    E = Maxwell3D.planewave(direction=ẑ, polarization=x̂, wavenumber=κ)
+    e = (n × E) × n
+
+    @hilbertspace j
+    @hilbertspace k
+
+    A = assemble(t[k,j], X, X)
+    S = assemble(s[k,j], X, X)
+    b = assemble(e[k], X)
+
+    Ai = BEAST.GMRESSolver(A; restart=1500, abstol=1e-5, maxiter=10_000)
+    u = Ai * b
+    u_Snorm = real(sqrt(dot(u, S*u)))
+    return (;u, X, u_Snorm)
+end
+
+α = 0.8
+h = collect(0.4 * α.^(0:13))  
+# h = collect(logrange(0.4, 0.05, length=8))
+@runsims payload h
+
+df = @collect_results
+
+using LinearAlgebra
+using Plots
+plot(df.h, real.(df.u_Snorm), marker=:.)
+
+# α = df.h[1] / df.h[2]
+W1 = reverse(df.u_Snorm)
+# W1 = reverse((df.h).^2.231)
+W2 = Iterators.drop(W1,1)
+W3 = Iterators.drop(W1,2)
+p = [log((w2-w1)/(w3-w2)) / log(1/α) for (w1,w2,w3) in zip(W1,W2,W3)]
+plot(collect(Iterators.drop(reverse(df.h),2)), p, marker=:.)
+
+refsol = df.u_Snorm[1]
+plot(log.(df.h), log.(abs.(df.u_Snorm .- refsol)), marker=:.)
+plot!(log.(df.h), 1.75 * log.(df.h))
+
+fcr, geo = facecurrents(df.u[1], df.X[1])
+import Plotly
+Plotly.plot(patch(geo, log10.(norm.(fcr))))
+