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

GMRES iteration 2, residual = 0.8433834622783686     

GMRES iteration 3, residual = 0.4577402159706648     

GMRES iteration 4, residual = 0.30210716607024823     

GMRES iteration 5, residual = 0.2192123586448595     

GMRES iteration 6, residual = 0.18262653379843952     

GMRES iteration 7, residual = 0.13313336375759927     

GMRES iteration 8, residual = 0.05080497746867199     

GMRES iteration 9, residual = 0.03392193652672886     

GMRES iteration 10, residual = 0.013092407942040654     

GMRES iteration 11, residual = 0.009112513240253705     

GMRES iteration 12, residual = 0.0007753709392715611     

GMRES iteration 13, residual = 0.00010301894044587265     

GMRES iteration 14, residual = 1.9220876415732964e-05     

GMRES iteration 15, residual = 4.005172772596476e-06     

GMRES iteration 16, residual = 1.7426778712070648e-06     

GMRES iteration 17, residual = 4.0853479856449564e-07     

GMRES iteration 18, residual = 1.8137917151343438e-07     

GMRES iteration 19, residual = 8.370880836771866e-08     

GMRES iteration 20, residual = 4.7861878871679306e-08     

GMRES iteration 21, residual = 2.845053334601645e-08     

GMRES iteration 22, residual = 1.553356429777655e-08     

GMRES iteration 23, residual = 6.4359503487042576e-09     
time= 298.69758677482605

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