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=10
order=4

screen = WorkPlane(Axes( (0,0,0), Z, X)).RectangleC(15,15).Face()

sp = Fuse(Sphere( (0,0,0), pi).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 =  17010 ndof_screen = 31110

with TaskManager():
    # V = HelmholtzSingleLayerPotentialOperator(fes_sphere, fes_sphere, kappa=kappa, intorder=10)
    # K = HelmholtzDoubleLayerPotentialOperator(fes_sphere, fes_sphere, kappa=kappa, intorder=10)
    # C = HelmholtzCombinedFieldOperator(fes_sphere, fes_sphere, kappa=kappa, intorder=10)
    C = HelmholtzCF(u*ds("sphere"), kappa)*v*ds
    u,v  = fes_sphere.TnT()
    Id = BilinearForm(u*v*ds).Assemble()

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 = 60.52440579389443
GMRes iteration 2, residual = 24.055821835237605
GMRes iteration 3, residual = 10.866381612939415
GMRes iteration 4, residual = 5.207724837026282
GMRes iteration 5, residual = 2.5606451105593298
GMRes iteration 6, residual = 1.2928262518532612
GMRes iteration 7, residual = 0.6567675233068863
GMRes iteration 8, residual = 0.35921392193432006
GMRes iteration 9, residual = 0.20001557577759385
GMRes iteration 10, residual = 0.10890790559448278
GMRes iteration 11, residual = 0.05889051992112595
GMRes iteration 12, residual = 0.03274406127301673
GMRes iteration 13, residual = 0.013561377443551775
GMRes iteration 14, residual = 0.00728664076375631
GMRes iteration 15, residual = 0.0036319265850675315
GMRes iteration 16, residual = 0.0020364922932962254
GMRes iteration 17, residual = 0.001213947151580113
GMRes iteration 18, residual = 0.0007228533895894088
GMRes iteration 19, residual = 0.0004301332817236075
GMRes iteration 20, residual = 0.0002583410832707025
GMRes iteration 21, residual = 0.00014417783739586982
GMRes iteration 22, residual = 7.770652618062729e-05
GMRes iteration 23, residual = 4.468014818822589e-05
GMRes iteration 24, residual = 2.211722906069809e-05
GMRes iteration 25, residual = 1.1247609798544885e-05
GMRes iteration 26, residual = 5.988321048133798e-06
GMRes iteration 27, residual = 3.2611582435840055e-06
GMRes iteration 28, residual = 1.7721145247161923e-06
GMRes iteration 29, residual = 1.0683563619743509e-06
GMRes iteration 30, residual = 5.955571812727601e-07
GMRes iteration 31, residual = 3.4514286396785986e-07
GMRes iteration 32, residual = 1.997649360837229e-07
GMRes iteration 33, residual = 1.1199302079006898e-07
GMRes iteration 34, residual = 6.802286055501477e-08
GMRes iteration 35, residual = 4.016335375213108e-08
GMRes iteration 36, residual = 1.8963293310895724e-08
GMRes iteration 37, residual = 1.0014341697466003e-08
GMRes iteration 38, residual = 5.59642497612889e-09

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

prostprocessing on screen

uscat = GridFunction(fes_screen)
with TaskManager(pajetrace=10**8):
    # uscat.Set(1j*kappa*V.GetPotential(gfu)-K.GetPotential(gfu) , definedon=mesh.Boundaries("screen"))
    # uscat.Set(HelmholtzCF(u*ds("sphere"), kappa)(gfu)  , definedon=mesh.Boundaries("screen"))
    uscat.Set(1j*kappa*HelmholtzSL(u*ds("sphere"),kappa)(gfu) - HelmholtzDL(u*ds("sphere"), kappa)(gfu), \
              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