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.0872552459587936     

GMRes iteration 2, residual = 0.8433834616996533     

GMRes iteration 3, residual = 0.4577402158645292     

GMRes iteration 4, residual = 0.3021071668771748     

GMRes iteration 5, residual = 0.2192123590957193     

GMRes iteration 6, residual = 0.18262653379883698     

GMRes iteration 7, residual = 0.13313336368773615     

GMRes iteration 8, residual = 0.050804977477992466     

GMRes iteration 9, residual = 0.03392193638851947     

GMRes iteration 10, residual = 0.013092407954857908     

GMRes iteration 11, residual = 0.00911251298161594     

GMRes iteration 12, residual = 0.0007753705757960248     

GMRes iteration 13, residual = 0.00010301881926497622     

GMRes iteration 14, residual = 1.9217923338645667e-05     

GMRes iteration 15, residual = 3.9873503025119514e-06     

GMRes iteration 16, residual = 1.6856364719722588e-06     

GMRes iteration 17, residual = 4.0295819605350466e-07     

GMRes iteration 18, residual = 1.6989529441468275e-07     

GMRes iteration 19, residual = 6.54489473000369e-08     

GMRes iteration 20, residual = 2.8704023125516513e-08     

GMRes iteration 21, residual = 1.633546726668888e-08     

GMRes iteration 22, residual = 1.120982242004986e-08     

GMRes iteration 23, residual = 6.368004702641874e-09     
time= 245.01112818717957

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)

WebGLScene