Embedded Trefftz-DG: Advection

In this section we finally consider an example that is typically not related to Trefftz method: A scalar linear transport problem, the advection equation. It reads as

\[\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} \newcommand{\bb}{\mathbf{b}} \begin{align*} \bb\cdot\nabla u &= f\quad \text{ in } \Omega,\\ u &= u_D \quad \text{ on } \partial \Omega_{\text{in}} :=\{\bx\in\partial \Omega\mid \bb\cdot\nx < 0\}. \end{align*}\end{split}\]

for a given velocity field \(\bb\) which we assume to be divergence-free.

As underlying DG discretization we choose the standard Upwind DG formulation which reads as

\[\begin{split}\begin{align} a_h(u,v) &= \sum_{K\in\Th} \Big\{ \int_K - u ~ \bb\cdot \nabla v\ dV +\int_{\partial K\setminus {\partial \Omega}{\text{in}} } \bb\cdot\nx\hat{u} v~ dS \Big\} \\ \ell(v) &= \sum_{K\in\Th} \int_K f v\ dV - \int_{{\partial \Omega}{\text{in}}} \bb \cdot \nx u_D v~ dS \end{align}\end{split}\]

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

b = CF((-sin(y),cos(x)))
exact = sin(3*x)*sin(3*y)
rhs = b*CF((exact.Diff(x),exact.Diff(y)))

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

bndc=exact
n = specialcf.normal(mesh.dim)
h = specialcf.mesh_size

jump_u = (u-u.Other())*n
jump_v = (v-v.Other())*n

a = BilinearForm(fes,symmetric=False)
a += - u * b * grad(v) * dx
a += b*n*IfPos(b*n,u,u.Other())*v*dx(element_boundary=True)

f = LinearForm(fes)
f += rhs*v * dx \
     - b*n*IfPos(b*n,0,bndc)*v * ds(skeleton=True)
with TaskManager():
    a.Assemble()
    f.Assemble()

eps = 10**-8
op = b*grad(u) * grad(v)[0]*dx
lop = rhs*grad(v)[0]*dx

with TaskManager():
    PP,uf = TrefftzEmbedding(op,fes,lop,eps)
PPT = PP.CreateTranspose()
with TaskManager():
    TA = PPT@a.mat@PP
    Tgfu = TA.Inverse()*(PPT*(f.vec-a.mat*uf))
    gfu = GridFunction(fes)
    gfu.vec.data = PP*Tgfu + uf
error = sqrt(Integrate((gfu-exact)**2, mesh))
print("svdt error : {:13.8e}".format(error))

svdt error : 1.51628610e-07

from ngsolve.webgui import Draw
Draw(gfu)

BaseWebGuiScene