Maxwell SL with preconditioner

from ngsolve import *
from netgen.occ import *
from ngsolve.webgui import Draw
from ngsolve.bem import *

face = Glue(Sphere((0,0,0),1).faces)
mesh = Mesh(OCCGeometry(face).GenerateMesh(maxh=0.1)).Curve(4)
Draw (mesh);

visplane = WorkPlane(Axes((0,0,0), Z, X)).RectangleC(5,5).Face()
vismesh = Mesh(OCCGeometry(visplane).GenerateMesh(maxh=0.2))
Draw (vismesh);

kappa = 3*pi
d = CF((0,-kappa,0))
E_inc = exp(1j*d*CF((x,y,z))) * CF((1,0,0))

Draw (E_inc[0], vismesh, animate_complex=True, order=3);

fes = HDivSurface(mesh, order=3, complex=True)
u,v = fes.TnT()

fespot = H1(mesh, order=4, complex=True)
upot, vpot = fespot.TnT()

with TaskManager():
    V1 = HelmholtzSL(u.Trace()*ds(bonus_intorder=4), kappa)*v.Trace()*ds(bonus_intorder=4)
    V2 = HelmholtzSL(div(u.Trace())*ds(bonus_intorder=4), kappa)*div(v.Trace())*ds(bonus_intorder=4)

lhs = (1j*kappa) * V1.mat + 1/(1j*kappa) * V2.mat

rhs = LinearForm(E_inc*v.Trace()*ds(bonus_intorder=3)).Assemble()

preconditioner is $\( (\cdot \times n) @ V @ (\cdot \times n) \)$

invMHd = BilinearForm(u.Trace()*v.Trace()*ds).Assemble().mat.Inverse(inverse="sparsecholesky")
n = specialcf.normal(3)
with TaskManager():
    Vrot = HelmholtzSL(Cross(u.Trace(),n)*ds(bonus_intorder=2), kappa)*Cross(v.Trace(),n)*ds(bonus_intorder=2)
    Vpot = HelmholtzSL(upot*ds(bonus_intorder=2), kappa)*vpot*ds(bonus_intorder=2)

surfcurl = BilinearForm(Cross(grad(upot).Trace(), n)*v.Trace()*ds).Assemble().mat
invMH1 = BilinearForm(upot*vpot*ds).Assemble().mat.Inverse(inverse="sparsecholesky")
                       
pre = invMHd @ ( kappa*Vrot.mat - 1/kappa*surfcurl@invMH1@Vpot.mat@invMH1@surfcurl.T)  @ invMHd

gfj = GridFunction(fes)
from time import time
ts = time()
with TaskManager():
    gfj.vec[:] = solvers.GMRes(A=lhs, b=rhs.vec, pre=pre, maxsteps=500, tol=1e-8)
te = time()
print ("time=", te-ts)

GMRes iteration 1, residual = 2.0872552459581626     

GMRes iteration 2, residual = 0.8433834616981585     

GMRes iteration 3, residual = 0.4577402158638526     

GMRes iteration 4, residual = 0.30210716687713585     

GMRes iteration 5, residual = 0.21921235909534736     

GMRes iteration 6, residual = 0.18262653379842578     

GMRes iteration 7, residual = 0.13313336368706075     

GMRes iteration 8, residual = 0.05080497747827218     

GMRes iteration 9, residual = 0.03392193638853567     

GMRes iteration 10, residual = 0.013092407954775282     

GMRes iteration 11, residual = 0.009112512981696528     

GMRes iteration 12, residual = 0.000775370575912981     

GMRes iteration 13, residual = 0.00010301881897470279     

GMRes iteration 14, residual = 1.921792249505858e-05     

GMRes iteration 15, residual = 3.987345458297536e-06     

GMRes iteration 16, residual = 1.6856198507385223e-06     

GMRes iteration 17, residual = 4.029576030425146e-07     

GMRes iteration 18, residual = 1.698942364374051e-07     

GMRes iteration 19, residual = 6.544785768660022e-08     

GMRes iteration 20, residual = 2.8702373462610824e-08     

GMRes iteration 21, residual = 1.6334233160269858e-08     

GMRes iteration 22, residual = 1.1209235550105296e-08     

GMRes iteration 23, residual = 6.368012713722374e-09     
time= 835.7865998744965

Draw( gfj[0], mesh, draw_vol=False, order=3, min=-2, max=2, animate_complex=True);

fes_screen = VectorH1(vismesh, order=3, complex=True)
E_screen = GridFunction(fes_screen)

E_screen.Set(1j*kappa*HelmholtzSL(u.Trace()*ds(bonus_intorder=4), kappa)(gfj) \
             -1/(1j*kappa)*grad(HelmholtzSL(div(u.Trace())*ds(bonus_intorder=4), kappa)(gfj)) \
             , definedon=vismesh.Boundaries(".*"))

Draw ( (E_inc-E_screen)[0], vismesh, animate_complex=True, min=-1, max=1)
Draw ( (E_inc)[0].real, vismesh, animate_complex=True, min=-1, max=1)
Draw ( (E_screen)[0].real, vismesh, animate_complex=True, min=-1, max=1)

BaseWebGuiScene