This page was generated from unit-11.3-bem-Helmholtz/BrakhageWerner.ipynb.

11.3.1 Helmholtz solver using Brakhage-Werner formulation

Combined field integral equations combine single and double layer integral operators, one simple option is the Brakhage-Werner formulation.

The solution is represented as

\[u = (i \kappa S - D) \phi,\]

where \(\phi\) solve the boundary integral equation

\[\big( \tfrac{1}{2} + K + i \kappa V \big) \phi = u_{in} \qquad \text{on} \, \Gamma\]

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

kappa=20
order=4

screen = WorkPlane(Axes( (0,0,0), Z, X)).RectangleC(15,15).Face()
sphere = Sphere( (0,0,0), pi)
screen = screen - sphere
sp = Fuse(sphere.faces)
screen.faces.name="screen"
sp.faces.name="sphere"
shape = Compound([screen,sp])

mesh = shape.GenerateMesh(maxh=5/kappa).Curve(order)
Draw (mesh);

fes_sphere = Compress(SurfaceL2(mesh, order=order, complex=True, definedon=mesh.Boundaries("sphere")))
u,v = fes_sphere.TnT()
fes_screen = Compress(SurfaceL2(mesh, order=order, dual_mapping=True, complex=True, definedon=mesh.Boundaries("screen")))
print ("ndof_sphere = ", fes_sphere.ndof, "ndof_screen =", fes_screen.ndof)

ndof_sphere =  64530 ndof_screen = 98325

with TaskManager():
    C = HelmholtzCF(u*ds("sphere"), kappa)*v*ds
    Id = BilinearForm(u*v*ds).Assemble()

with TaskManager():
    lhs = 0.5 * Id.mat + C.mat
    source = exp(1j * kappa * x)
    rhs = LinearForm(-source*v*ds).Assemble()

gfu = GridFunction(fes_sphere)
pre = BilinearForm(u*v*ds, diagonal=True).Assemble().mat.Inverse()
with TaskManager():
    gfu.vec[:] = solvers.GMRes(A=lhs, b=rhs.vec, pre=pre, maxsteps=40, tol=1e-8)

GMRES iteration 1, residual = 119.67056861548505
GMRES iteration 2, residual = 54.52890573982241
GMRES iteration 3, residual = 27.43583083346146
GMRES iteration 4, residual = 14.630009765901004
GMRES iteration 5, residual = 8.056881046815842
GMRES iteration 6, residual = 4.509288098389363
GMRES iteration 7, residual = 2.564403593372173
GMRES iteration 8, residual = 1.4899517634512522
GMRES iteration 9, residual = 0.902481528936717
GMRES iteration 10, residual = 0.575820167038915
GMRES iteration 11, residual = 0.35606560095017303
GMRES iteration 12, residual = 0.23167701960538362
GMRES iteration 13, residual = 0.14055733273567553
GMRES iteration 14, residual = 0.08643232124400872
GMRES iteration 15, residual = 0.047855814477777474
GMRES iteration 16, residual = 0.025506372247077018
GMRES iteration 17, residual = 0.014746624601715663
GMRES iteration 18, residual = 0.008510791472470533
GMRES iteration 19, residual = 0.00551036915384327
GMRES iteration 20, residual = 0.0036187074818088173
GMRES iteration 21, residual = 0.002430613803273291
GMRES iteration 22, residual = 0.0015919307009793363
GMRES iteration 23, residual = 0.0010083373431220167
GMRES iteration 24, residual = 0.0006321666858694492
GMRES iteration 25, residual = 0.00040090381921028417
GMRES iteration 26, residual = 0.00026043699290784607
GMRES iteration 27, residual = 0.00017182267969723236
GMRES iteration 28, residual = 0.00010868824103024936
GMRES iteration 29, residual = 6.329089442291172e-05
GMRES iteration 30, residual = 3.863213309648279e-05
GMRES iteration 31, residual = 2.1392678447595102e-05
GMRES iteration 32, residual = 1.229005257860321e-05
GMRES iteration 33, residual = 7.293316710365911e-06
GMRES iteration 34, residual = 4.7962049336238685e-06
GMRES iteration 35, residual = 3.114162424918259e-06
GMRES iteration 36, residual = 2.00842578952865e-06
GMRES iteration 37, residual = 1.2865040351225297e-06
GMRES iteration 38, residual = 7.999753526082447e-07
GMRES iteration 39, residual = 5.195825990318129e-07
GMRES iteration 40, residual = 3.36270833792524e-07
WARNING: GMRES did not converge to TOL

Draw (gfu, order=5, min=-1, max=1);

prostprocessing on screen

uscat = GridFunction(fes_screen)
with TaskManager(pajetrace=10**8):
    intop = 1j*kappa*HelmholtzSL(u*ds("sphere"),kappa) + (-1)*HelmholtzDL(u*ds("sphere"), kappa)
    uscat.Set(intop(gfu, mesh.Boundaries("screen")), definedon=mesh.Boundaries("screen"))

print ("Scattered field")
Draw (uscat, mesh, min=-1,max=1, animate_complex=True, order=4);

Scattered field

uin = mesh.BoundaryCF( {"screen": source }, default=0)
print ("Total field")
Draw (uin-uscat, mesh, min=-1,max=1, animate_complex=True, order=4);

Total field

Scattering from sphere with \(D = 50 \lambda\). About 5 min on Macbook Apple M4 Pro