Embedded Trefftz-DG: Helmholtz

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

We consider 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{\Vhp}{V^p(\Th)} \newcommand{\bT}{\mathbf{T}} \newcommand{\bW}{\mathbf{W}} \newcommand{\bw}{\mathbf{w}} \newcommand{\bl}{\mathbf{l}} \newcommand{\bM}{\mathbf{M}} \newcommand{\bL}{\mathbf{L}} \newcommand{\bA}{\mathbf{A}} \newcommand{\bU}{\mathbf{U}} \newcommand{\bV}{\mathbf{V}} \newcommand{\calL}{\mathcal{L}} \newcommand{\bu}{\mathbf{u}} \newcommand{\IT}{\mathbb{T}} \newcommand{\calG}{\mathcal{G}} \newcommand{\be}{\mathbf{e}} \newcommand{\bx}{{\mathbf x}} \newcommand{\inner}[1]{\langle #1 \rangle} \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}\]

Standard polynomial Trefftz functions for the Helmholtz equation do not exist, to circumvent this problem, we weaken our condition in the Trefftz space. We introduce a projection \(\Pi\) that is yet to be defined and define the weak Trefftz space and the embedded weak Trefftz DG method:

\[\begin{split}\begin{align} \text{Find }u_{\IT}\in \IT^p(\Th)&,~\text{ s.t. } a_h(u_{\IT},v_{\IT})=\ell(v_{\IT})\qquad \forall v_{\IT}\in \IT^p(\Th)\quad \text{ with } \\ \IT^p(\Th)&=\{v\in V^k(\mathcal T_h),\ \Pi \calL v=0\}. \label{eq:weakTspace} % \\ &\IT^p(\Th)=\{v\in L^2(\dom) \sst \restr{v}{K}\in\IT(K),\forall K\in\Th\} \end{align}\end{split}\]

For the Helmholtz problem we choose \(\Pi:V^{k}(\mathcal T_h)\rightarrow V^{k-2}(\mathcal T_h)\) the \(L^2\) orthogonal projection. This way, we can define the matrix \(\bW\) as

\[\begin{align} \label{def:W3} (\bW)_{ij}&=\inner{\calL\phi_j,\tilde\phi_i}_{0,h}. \end{align}\]

with test functions \(\tilde\phi_i\in V^{k-2}(\mathcal T_h)\) and \(\calL=-\Delta u -\omega^2 u\).

def TrefftzHelmholtzEmb(fes):
    mesh = fes.mesh
    k = fes.globalorder
    u = fes.TrialFunction()

    Lap = lambda u : sum(Trace(u.Operator('hesse')))
    fes2 = L2(mesh, order=order-2, dgjumps=True,complex=True)
    v = fes2.TestFunction()
    op = -u*v*dx - Lap(u)*v*dx
    PP = TrefftzEmbedding(op,fes,test_fes=fes2)
    return PP

def EmbeddedBasisFunctionsAsMultiDim(Emb, fes):
    gfshow = GridFunction(fes, multidim=0)
    gf = GridFunction(fes)
    coefvec = Emb.CreateRowVector()
    for i in range(Emb.width):
        coefvec[:] = 0
        coefvec[i] = 1
        gf.vec.data = Emb * coefvec
        gfshow.AddMultiDimComponent(gf.vec)
    return gfshow

order = 3
mesh = Mesh(unit_square.GenerateMesh(maxh=1))
fes = L2(mesh, order=order, dgjumps=True,complex=True)
PP = TrefftzHelmholtzEmb(fes)
gfshow = EmbeddedBasisFunctionsAsMultiDim(PP,fes)
Draw (gfshow, mesh, interpolate_multidim=False, animate=False, autoscale=True, settings={"subdivision":20})

BaseWebGuiScene

We compare this to the planewave space

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

order = 3
mesh = Mesh(unit_square.GenerateMesh(maxh=1))
fes = trefftzfespace(mesh, order=order, eq="helmholtz",dgjumps=True,complex=True)
gfshow = GridFunction(fes, multidim=0)
gf = GridFunction(fes)
for i in range(fes.ndof):
    gf.vec[:] = 0
    gf.vec[i] = 1
    gfshow.AddMultiDimComponent(gf.vec)
Draw (gfshow, mesh, interpolate_multidim=False, animate=False, autoscale=False, min=0.8,max=1, settings={"subdivision":20})

BaseWebGuiScene

mesh = Mesh(unit_square.GenerateMesh(maxh=.3))
order = 4
fes = L2(mesh, order=order, complex=True, dgjumps=True)
u,v = fes.TnT()

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

jump_u = (u-u.Other())*n
jump_v = (v-v.Other())*n
jump_du = (grad(u)-grad(u.Other()))*n
jump_dv = (grad(v)-grad(v.Other()))*n
mean_u = 0.5 * ((u)+(u.Other()))
mean_du = 0.5 * (grad(u)+grad(u.Other()))
mean_dv = 0.5 * (grad(v)+grad(v.Other()))

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

a += -(jump_u*(mean_dv)+mean_du*(jump_v)) * dx(skeleton=True)
a += -1/(omega*1j)*beta*(jump_du*(jump_dv)) * 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(fes)
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()

PP = TrefftzHelmholtzEmb(fes)
PPT = PP.CreateTranspose()
with TaskManager():
    TA = PPT@a.mat@PP
    TU = TA.Inverse()*(PPT*f.vec)
    tpgfu = GridFunction(fes)
    tpgfu.vec.data = PP*TU
error = sqrt(Integrate((tpgfu-exact)*Conj(tpgfu-exact), mesh).real)
print("error ",error)

error  8.581214563202292e-08