Trefftz + Tent pitching

Tent pitching techniques generate a space-time mesh, which complies with the causality properties of the hyperbolic PDE. The resulting mesh consists of tent shaped objects, the height of rach tent is complient with the local wavespeed. The PDE can be solved explicitly in each tent and independant tents can be solves in parallel.

[1]:
from ngstrefftz import *
from ngsolve.TensorProductTools import *
from netgen.occ import *
from ngsolve import *
from ngsolve.webgui import *

Constructing a tent-pitched mesh

We use ngstents to construct tent-pitched meshes.

[2]:
def Lmesh(maxh):
    mu = 0.4
    mp = MeshingParameters (maxh = maxh)
    minh = math.pow(maxh,1/mu)
    refpoints = int(1/minh)
    for i in range(0, refpoints+1):
        for j in range(0, refpoints+1):
            xk = (i/refpoints-0.5)*2
            yk = (j/refpoints-0.5)*2
            r = sqrt(xk*xk+yk*yk)
            mp.RestrictH (x=xk, y=yk, z=0, h=max(minh, maxh*math.pow(r,1-mu)) )

    wp = WorkPlane()
    wp.MoveTo(-1,-1).LineTo(0,-1).LineTo(0,0).LineTo(1,0).LineTo(1,1).LineTo(-1,1).Close()
    geo = wp.Face()
    geo = OCCGeometry(geo,dim=2)

    return Mesh(geo.GenerateMesh(mp=mp))

maxh = 0.13
mesh = Lmesh(maxh)
Draw(mesh)
[2]:
BaseWebGuiScene
[3]:
dt = 0.05   # tent slab height (in time)
wavespeed = 1
ts = TentSlab(mesh, method="edge")
ts.SetMaxWavespeed(wavespeed)
ts.PitchTents(dt=dt, local_ct=True)
print("max slope", ts.MaxSlope())
print("n tents", ts.GetNTents())
max slope 1.0
n tents 780

Trefftz-DG for the acoustic wave equation + Tents

[4]:
for i in range(0,len(mesh.GetBoundaries())):
    mesh.ngmesh.SetBCName(i,"neumann")

order=4
delta = 500
ee = exp(-(x+1)**2*delta )
initc=((ee,-2*(x+1)*delta * ee, 0, 0))
bdd = CoefficientFunction(( 0,0,0,0))

TT=TWave(order,ts,CoefficientFunction(wavespeed))
TT.SetInitial(initc)
TT.SetBoundaryCF(bdd)

V = L2(mesh,order=order)
u = GridFunction(V,"u")
gfshow = GridFunction(V, multidim=0)
t=0
SetNumThreads(4)
with TaskManager():
    while t < 1.5:
        TT.Propagate()
        t += dt
        TT.GetWave(u)
        gfshow.AddMultiDimComponent(u.vec)
Draw (gfshow,mesh,interpolate_multidim=True,deformation=True,animate=True,autoscale=False,min=0,max=0.3)
[4]:
BaseWebGuiScene

Quasi-Trefftz + Tents

if a non-constant coefficient function as wavespeed is passed to \(\texttt{TWave}\) the quasi-Trefftz method will be used, see also the notebook on the quasi-Treffty method.

[5]:
order = 4
delta=500
dt=0.02
ee = exp(-(x)**2*delta )
initc=((ee,-2*(x)*delta * ee, 0, 0))
bdd = CoefficientFunction(( 0,0,0,0))
mesh = Mesh(unit_square.GenerateMesh(maxh = 0.05))
for i in range(0,len(mesh.GetBoundaries())):
    mesh.ngmesh.SetBCName(i,"neumann")

ts = TentSlab(mesh, method="edge")
ts.SetMaxWavespeed(y+1)
ts.PitchTents(dt=dt, local_ct=True)
TT=TWave(order,ts,CoefficientFunction(y+1))
TT.SetInitial(initc)
TT.SetBoundaryCF(bdd)

V = L2(mesh,order=order)
u = GridFunction(V,"u")
gfshow = GridFunction(V, multidim=0)
t=0
SetNumThreads(4)
with TaskManager():
    while t < 0.8:
        TT.Propagate()
        t += dt
        TT.GetWave(u)
        gfshow.AddMultiDimComponent(u.vec)
Draw (gfshow,mesh,interpolate_multidim=True,deformation=True,animate=True,autoscale=False,min=0,max=0.3)
[5]:
BaseWebGuiScene