Quasi-Trefftz DG: Wave equation

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

We consider the wave operator

\[\begin{align*} \begin{split} (\square_G f)(\mathbf{x},t):= \Delta f(\mathbf{x},t)-G(\mathbf{x})\partial_t^2 f(\mathbf{x},t). \end{split} \end{align*}\]

with smooth coefficient \(G(\mathbf{x})\). Constructing a basis for a traditional Trefftz space (i.e. a space of functions with \(\square_G f=0\)) is not possible. The crucial idea is that we want to relax the Trefftz porperty to

\[\square_G f=\mathcal{O}(\|(\mathbf{x},t)-(\mathbf{x}_K,t_K)\|^q),\]

with respect to the center of a mesh element \(K\) and up to some \(q\). This leads to the definition of a new quasi-Trefftz space: For an element \(K\) in a space-time mesh let

\[\begin{align*} \begin{split} \mathbb{T}^p(K):=\big\{ f\in\mathbb{P}^p(K) \mid D^{i}\square_G f(\mathbf{x}_K,t_K)=0,\ \forall i\in \mathbb{N}^{n+1}_0, |i|<p-1 \big\}, \qquad p\in \mathbb{N}. \end{split} \end{align*}\]

For this space we are able to construct a basis. We then introduce a space-time DG method with test and trial functions that are locally quasi-Trefftz. The example below shows an acoustic wave propagating through a material with \(G(x,y)=y+1\) and homogeneous Neumann boundary conditions.

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

Test solution

ca=2.5
sig_ex = CoefficientFunction((-ca*(x+y+1)**(ca-1) * exp(-sqrt(2*ca*(ca-1))*z),\
                            -ca*(x+y+1)**(ca-1) * exp(-sqrt(2*ca*(ca-1))*z)) )
v_ex = -sqrt(2*ca*(ca-1))*(x+y+1)**ca * exp(-sqrt(2*ca*(ca-1))*z)
c = CoefficientFunction((x+y+1))
gD=v_ex

Quasi-Trefftz space

Construct a quasi-Trefftz space for the forst order wave equation by setting eq="foqtwave"

fes = trefftzfespace(mesh,order=3,dgjumps=True,eq="foqtwave")
fes.SetCoeff(c)

D = fes.mesh.dim - 1

vsig = fes.TrialFunction()
wtau = fes.TestFunction()
vsigo=vsig.Other()
wtauo=wtau.Other()

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)
dttau = CoefficientFunction(tuple([grad(wtau)[i,D] for i in range(D)]))
gtau = sum([grad(wtau)[i*(D+2)] for i in range(D)])
dtw =grad(wtau)[(D+1)*(D+1)-1]
gw=CoefficientFunction(tuple([grad(wtau)[D,i] for i in range(D)]))
a += (- v*(gtau + pow(c,-2)*dtw) ) * dx
a += ( - sig*(dttau + gw ) ) * dx
a += spacelike * ( pow(c,-2)*IfPos(n_t,v,vo)*jump_wt + IfPos(n_t,sig,sigo)*jump_taut) * dx(skeleton=True)
a += timelike * ( mean_v*jump_taux + mean_sig*jump_wx + alpha*jump_vx*jump_wx + beta*jump_sigx*jump_taux ) * dx(skeleton=True )
a += ( pow(c,-2)*v*w + sig*tau) * IfPos(n_t,1,0)*spacelike * ds(skeleton=True)
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()

gfu = GridFunction(fes, name="uDG")
gfu.vec.data = a.mat.Inverse()*f.vec
Draw(gfu,clipping={"z":-1, "dist":-.33, "function":True})

BaseWebGuiScene

Integrate((gfu-CoefficientFunction((sig_ex,v_ex)))**2,mesh)

1.70851963649639e-07