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= f &\text{ in } \dom, \\ \frac{\partial u}{\partial \nx} - i\omega 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,omega=1):
    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 = -omega**2*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

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}} - i\frac{\beta}{\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}} + i\frac{\delta}{\omega}{\nabla u\cdot\nx}{\overline{\nabla v}\cdot\nx} \right) dS \\ \ell(v) &= \int_\Omega fv\ dV + \int_{\Fh^\text{bnd}}\left( (1-\delta)g\overline{v} - i\frac{\delta}{\omega}g\overline{\nx\cdot\nabla v}\right) dS \end{align}\end{split}\]

def dghelmholtz(fes,omega,rhs=0,bndc=0,bnd=".*"):
    u,v = fes.TnT()
    h = specialcf.mesh_size
    order = fes.globalorder
    alpha = order**2/h
    beta = h/order
    delta = 0.1*omega*h/order
    IP = lambda u,v: InnerProduct(u,v)

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

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

    f = LinearForm(fes)
    f += rhs * v * dx
    f += (1-delta)*bndc*v * ds(bnd,skeleton=True)
    f += -1j*delta/omega*bndc*grad(v)*n * ds(bnd,skeleton=True)

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

mesh = Mesh(unit_square.GenerateMesh(maxh=.3))
order = 4
fes = L2(mesh, order=order, complex=True, dgjumps=True)
a,f = dghelmholtz(fes,omega,bndc=bndc)
PP = TrefftzHelmholtzEmb(fes)
PPT = PP.CreateTranspose()
with TaskManager():
    TA = PPT@a.mat@PP
    TU = TA.Inverse()*(PPT*f.vec)
    gfu = GridFunction(fes)
    gfu.vec.data = PP*TU
error = sqrt(Integrate((gfu-exact)*Conj(gfu-exact), mesh).real)
print("error ",error)

error  6.622323484588101e-08

air = Circle((0.5, 0.5), 0.8).Face()
air.edges.name = 'outer'
scatterer = MoveTo(0.7, 0.3).Rectangle(0.05, 0.4).Face()
scatterer.edges.name = 'scat'
geo = OCCGeometry(air - scatterer, dim=2)
mesh = Mesh(geo.GenerateMesh(maxh=0.05))
omega = 25
rhs = 3e2*exp(-(40**2)*((x-0.5)*(x-0.5) + (y-0.5)*(y-0.5)))

fes = L2(mesh, order=order, complex=True, dgjumps=True)
a,f = dghelmholtz(fes,omega,rhs=rhs,bnd="outer")
PP = TrefftzHelmholtzEmb(fes,omega)
PPT = PP.CreateTranspose()
with TaskManager():
    TA = PPT@a.mat@PP
    TU = TA.Inverse()*(PPT*f.vec)
    gfu = GridFunction(fes)
    gfu.vec.data = PP*TU

Draw(gfu, mesh, animate_complex=True,deformation=True);