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

GMRES iteration 2, residual = 0.8433834622793179     

GMRES iteration 3, residual = 0.4577402159711756     

GMRES iteration 4, residual = 0.30210716607030697     

GMRES iteration 5, residual = 0.21921235864520172     

GMRES iteration 6, residual = 0.18262653379855692     

GMRES iteration 7, residual = 0.13313336375798615     

GMRES iteration 8, residual = 0.050804977468762766     

GMRES iteration 9, residual = 0.0339219365268719     

GMRES iteration 10, residual = 0.013092407942175605     

GMRES iteration 11, residual = 0.009112513240302019     

GMRES iteration 12, residual = 0.0007753709394000543     

GMRES iteration 13, residual = 0.00010301894074776465     

GMRES iteration 14, residual = 1.9220877064974306e-05     

GMRES iteration 15, residual = 4.005176556115273e-06     

GMRES iteration 16, residual = 1.7426893369405792e-06     

GMRES iteration 17, residual = 4.0853494451285246e-07     

GMRES iteration 18, residual = 1.8137945969065322e-07     

GMRES iteration 19, residual = 8.370894435479248e-08     

GMRES iteration 20, residual = 4.7862205914835166e-08     

GMRES iteration 21, residual = 2.8450592763679095e-08     

GMRES iteration 22, residual = 1.5533527280635323e-08     

GMRES iteration 23, residual = 6.435902470752297e-09     
time= 361.3811936378479

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