Helmholtz

The Helmholtz equation with Robin boundary conditions

\[\begin{split}\newcommand{\Th}{{\mathcal{T}_h}} \newcommand{\Fh}{\mathcal{F}_h} \newcommand{\dom}{\Omega} \newcommand{\jump}[1]{[\![ #1 ]\!]} \newcommand{\tjump}[1]{[\![{#1} ]\!]_\tau} \newcommand{\avg}[1]{\{\!\!\{#1\}\!\!\}} \newcommand{\nx}{n_\mathbf{x}} \newcommand{\IT}{\mathbb{T}} \newcommand{\bx}{\mathbf{x}} \newcommand{\sst}{\;\text{s.t.}\;} \begin{align*} \begin{cases} -\Delta u - \omega^2 u= 0 &\text{ in } \dom, \\ \frac{\partial u}{\partial \nx} + i u = g &\text{ on } \partial \dom. \end{cases} \end{align*}\end{split}\]

from ngsolve import *
from ngstrefftz import *
from netgen.occ import *

A Trefftz space for the Helmholtz equation in two dimensions is given by the (non-polynomial) space of plane wave functions

\[\begin{align*} \IT^p=\{e^{-i\omega(d_j\cdot \bx)} \sst j=-p,\dots,p\}. \end{align*}\]

mesh = Mesh(unit_square.GenerateMesh(maxh=.3))
fes = trefftzfespace(mesh,order=3,eq="helmholtz",complex=True,dgjumps=True)
fes2 = trefftzfespace(mesh,order=3,eq="helmholtzconj",complex=True,dgjumps=True)

omega=1
n = specialcf.normal(2)
exact = exp(1j*sqrt(0.5)*(x+y))
gradexact = CoefficientFunction((sqrt(0.5)*1j*exact, sqrt(0.5)*1j*exact))
bndc = gradexact*n + 1j*omega*exact
eps = 10**-7

We consider the DG-scheme given by

\[\begin{split}\begin{align} a_h(u,v) &= \sum_{K\in\Th}\int_K \nabla u\nabla v-\omega^2 uv\ dV -\int_{\Fh^\text{int}}\left(\avg{\nabla u}\jump{v}+\jump{u} \avg{\overline{\nabla v}} \right) dS \nonumber \\ &\qquad+\int_{\Fh^\text{int}} \left( i\alpha \omega\jump{u}\jump{\overline{v}} - \frac{\beta}{i\omega}\jump{\nabla u}\jump{\overline{\nabla v}} \right) dS -\int_{\Fh^\text{bnd}}\delta\left(\nx\cdot\nabla u \overline{v}+u \overline{\nx\cdot\nabla v}\right) dS\\ \nonumber &\qquad+\int_{\Fh^\text{bnd}} \left( i(1-\delta)\omega{u}{\overline{v}} - \frac{\delta}{i\omega}{\nabla u}{\overline{\nabla v}} \right) dS \\ \ell(v) &= \int_{\Fh^\text{bnd}}\left( (1-\delta)g\overline{v} - \frac{\delta}{i\omega}g\overline{\nx\cdot\nabla v}\right) dS \end{align}\end{split}\]

h = specialcf.mesh_size
alpha = 1/(omega*h)
beta = omega*h
delta = omega*h

u = fes.TrialFunction()
v = fes.TestFunction()
if fes2 is not None:
    v = fes2.TestFunction()
jump = lambda u: (u-u.Other())*n
mean = lambda u: 0.5 * ((u)+(u.Other()))

a = BilinearForm(fes,fes2)
a += grad(u)*(grad(v))*dx - omega**2*u*(v)*dx

a += -(jump(u)*(mean(grad(v)))+mean(grad(u))*jump(v)) * dx(skeleton=True)
a += -1/(omega*1j)*beta*jump(grad(u))*(jump(grad(v))) * dx(skeleton=True)
a += omega*1j*alpha*jump(u)*jump(v) * dx(skeleton=True)

a += -delta*(u*(grad(v))*n+grad(u)*n*(v)) * ds(skeleton=True)
a += -1/(omega*1j)*delta*(grad(u)*n)*((grad(v))*n) * ds(skeleton=True)
a += omega*1j*(1-delta)*u*(v) * ds(skeleton=True)

f = LinearForm(fes2)
f += -1/(omega*1j)*delta*bndc*(grad(v))*n*ds(skeleton=True)
f += (1-delta)*bndc*(v)*ds(skeleton=True)

with TaskManager():
    a.Assemble()
    f.Assemble()

gfu = GridFunction(fes)
with TaskManager():
    gfu.vec.data = a.mat.Inverse() * f.vec

error = sqrt(Integrate((gfu-exact)*Conj(gfu-exact), mesh).real)
print("error ",error)

error  1.4124364281169749e-06

from ngsolve.webgui import Draw
Draw(gfu)

BaseWebGuiScene