Embedded Trefftz-DG: Wave Equation

from ngsolve import *
from ngstrefftz import *
from netgen.occ import *
from ngsolve.TensorProductTools import MakeTensorProductMesh, SegMesh
from ngsolve import *
from ngsolve.webgui import Draw

We aim to solve the first order wave equation

\[\begin{split}\begin{align} \nabla\cdot\sigma + c^{-2}\partial_t v=0 \\ \nabla v + \partial_t \sigma =0\\ v=v_0,\ \sigma=\sigma_0 \text{ at } t=0\\ v=g_D \text{ on }\Gamma_D \end{align}\end{split}\]

and compare to the exact solution given by

from numpy import pi
c=1
sq = sqrt(2.0)
sig_ex = CoefficientFunction((-cos(pi*x)*sin(pi*y)*sin(pi*z*c*sq)/sq,
    -sin(pi*x)*cos(pi*y)*sin(pi*z*c*sq)/sq))
v_ex = sin(pi*x)*sin(pi*y)*cos(pi*z*c*sq)*c
gD = v_ex

Create tensor-product mesh

We consider a space-time mesh with time along the z-Axis.

meshx = Mesh(unit_square.GenerateMesh(maxh=0.2))
mesht = Mesh(SegMesh(5,0,1,periodic=False) )
mesh = Mesh(MakeTensorProductMesh(meshx,mesht))
from ngsolve.webgui import Draw
Draw(mesh)

BaseWebGuiScene

Trefftz-DG formulation

The space-time Trefftz-DG formulation is given by

\[\begin{split}\begin{align} \text{Seek}&\quad (v_{hp},\sigma_{hp})\in V_p (\mathcal{T}_h)\nonumber\\ \text{such that}&\quad \mathcal{A}(v_{hp},\sigma_{hp}; w ,\tau )=\ell( w ,\tau )\quad \forall ( w ,\tau )\in V_p (\mathcal{T}_h), \nonumber\\ \text{where}&\nonumber\\ \mathcal{A}(v_{hp},&\sigma_{hp}; w ,\tau ):= %-\sum_{K\in\mathcal{T}_h} \int_K\bigg(v\Big(\nabla\cdot\tau+c^{-2}\partial_t w \Big) +\sigma_{hp}\cdot\Big(\partial_t \tau +\nabla w \Big)\bigg) %\nonumber\\ +& \int_{F^{space}}\big(c^{-2}v_{hp}^-[[{w}]]_t+\sigma_{hp}^-\cdot[[{\tau}]]_t+v_{hp}^-[[{\tau}]]_N+\sigma_{hp}^-\cdot[[{w}]]_N\big) \nonumber \\ +&\int_{F^{time}}\!\! \big( \{\{{v_{hp}}\}\}[[{\tau }]]_N+\{\{{\sigma_{hp}}\}\}\cdot[[{ w }]]_N +\alpha[[{v_{hp}}]]_N\cdot[[{ w }]]_N+ \beta[[{\sigma_{hp}}]]_N[[{\tau }]]_N \big) \nonumber\\ +&\int_{F^T} (c^{-2}v_{hp} w +\sigma_{hp} \cdot\tau ) +\int_{F^D} \big(\sigma\cdot n_\Omega^x\, w +\alpha v_{hp} w \big) \nonumber\\ \ell( w ,&\tau ):= \int_{F^O} ( c^{-2}v_0 w +\sigma_0\cdot \tau ) +\int_{F^D} g_D\big(\alpha w -\tau\cdot n_\Omega^x\big) \nonumber \end{align}\end{split}\]

Note that for the Trefftz discretization the volume terms vanish. For the embedded Trefftz method we first set up the system on the full polynomial space:

order=3
X = L2(mesh,order=order,dgjumps=True)
fes=FESpace([X for _ in range(mesh.dim)])

D = fes.mesh.dim - 1

vsig = fes.TrialFunction()
wtau = fes.TestFunction()
vsigo = CoefficientFunction(tuple([vsig[i].Other() for i in  range(D+1)]))
wtauo = CoefficientFunction(tuple([wtau[i].Other() for i in  range(D+1)]))

v = vsig[D]
sig = CoefficientFunction(tuple([vsig[i] for i in  range(D)]))
w = wtau[D]
tau = CoefficientFunction(tuple([wtau[i] for i in  range(D)]))

vo = vsigo[D]
sigo = CoefficientFunction(tuple([vsigo[i] for i in  range(D)]))
wo = wtauo[D]
tauo = CoefficientFunction(tuple([wtauo[i] for i in  range(D)]))

h = specialcf.mesh_size
n = specialcf.normal(D+1)
n_t = n[D]/Norm(n)
n_x = CoefficientFunction( tuple([n[i]/Norm(n) for i in  range(D)]) )

mean_v = 0.5*(v+vo)
mean_w = 0.5*(w+wo)
mean_sig = 0.5*(sig+sigo)
mean_tau = 0.5*(tau+tauo)

jump_vx = ( v - vo ) * n_x
jump_wx = ( w - wo ) * n_x
jump_sigx = (( sig - sigo ) * n_x)
jump_taux = (( tau - tauo ) * n_x)

jump_vt = ( v - vo ) * n_t
jump_wt = ( w - wo ) * n_t
jump_sigt = ( sig - sigo ) * n_t
jump_taut = ( tau - tauo ) * n_t

timelike = n_x*n_x # n_t=0
spacelike = n_t**2 # n_x=0

#params
alpha=0.5
beta=0.5

a = BilinearForm(fes)
#space like faces, w/o x jump ASSUME TENSOR MESH
a += spacelike * ( pow(c,-2)*IfPos(n_t,v,vo)*jump_wt + IfPos(n_t,sig,sigo)*jump_taut) * dx(skeleton=True)
#time like faces
a += timelike * ( mean_v*jump_taux + mean_sig*jump_wx + alpha*jump_vx*jump_wx + beta*jump_sigx*jump_taux ) * dx(skeleton=True )        #t=T (or *x)
a += ( pow(c,-2)*v*w + sig*tau) *  IfPos(n_t,1,0)*spacelike * ds(skeleton=True)
#dirichlet boundary 'timelike'
a += ( sig*n_x*w + alpha*v*w ) * timelike *   ds(skeleton=True)

f = LinearForm(fes)
f += ( pow(c,-2)*v_ex*w + sig_ex*tau ) * IfPos(n_t,0,1)*spacelike *  ds(skeleton=True)
f += ( gD * (alpha*w - tau*n_x) ) * timelike *   ds(skeleton=True)

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

To apply the embedded Trefftz framework we take

\[\begin{split} \DeclareMathOperator\dt{\frac{\partial}{\partial t}} \newcommand{\calL}{\mathcal{L}} \begin{align*} \calL= \begin{pmatrix} \nabla\cdot & c^{-2}\dt \\ \nabla & \dt \end{pmatrix} \qquad u= \begin{pmatrix} \mathbf{\sigma}\\ v \end{pmatrix} \end{align*}\end{split}\]

dtsig = CoefficientFunction(tuple([grad(vsig[i])[D] for i in range(D)]))
gsig = sum([grad(vsig[i])[i] for i in range(D)])
dtv =grad(vsig[D])[D]
gv=CoefficientFunction(tuple([grad(vsig[D])[i] for i in range(D)]))
dttau = CoefficientFunction(tuple([grad(wtau[i])[D] for i in range(D)]))
gtau = sum([grad(wtau[i])[i] for i in range(D)])
dtw =grad(wtau[D])[D]
gw=CoefficientFunction(tuple([grad(wtau[D])[i] for i in range(D)]))
op = ( (dtsig+gv)*(dttau+gw)+(dtv+gsig)*(dtw+gtau) )*dx

eps=10**-12
with TaskManager():
    PP = TrefftzEmbedding(op,fes,eps)
PPT = PP.CreateTranspose()
with TaskManager():
    TA = PPT@a.mat@PP
    TU = TA.Inverse()*(PPT*f.vec)
    tgfu = GridFunction(fes)
    tgfu.vec.data = PP*TU

Draw(tgfu.components[0],clipping={"z":-1, "dist":-.11, "function":True})

BaseWebGuiScene

error = Integrate((CF(tgfu.components)-CoefficientFunction((sig_ex,v_ex)))**2,mesh)
print(error)

3.535809015396434e-07