Embedded Trefftz-DG: Poisson

from ngsolve import *
from ngstrefftz import *
from netgen.occ import *
SetNumThreads(4)

We consider the Poisson equation with Dirichlet 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} \DeclareMathOperator\Ker{ker} \begin{align*} \begin{cases} \Delta u = f &\text{ in } \dom, \\ u=g &\text{ on } \partial \dom. \end{cases} \end{align*}\end{split}\]

We consider again the SIP-DG method

mesh = Mesh(unit_cube.GenerateMesh(maxh=.2))
order = 4
fes = L2(mesh, order=order,  dgjumps=True)

exact = sin(x)*sin(y)*sin(z)
rhs = 3*sin(x)*sin(y)*sin(z)
bndc=exact
alpha = 4
n = specialcf.normal(mesh.dim)
h = specialcf.mesh_size
u = fes.TrialFunction()
v = fes.TestFunction()

jump_u = (u-u.Other())*n
jump_v = (v-v.Other())*n
mean_dudn = 0.5 * (grad(u)+grad(u.Other()))
mean_dvdn = 0.5 * (grad(v)+grad(v.Other()))

a = BilinearForm(fes,symmetric=True)
aterms = grad(u)*grad(v) * dx \
    +alpha*order**2/h*jump_u*jump_v * dx(skeleton=True) \
    +(-mean_dudn*jump_v-mean_dvdn*jump_u) * dx(skeleton=True) \
    +alpha*order**2/h*u*v * ds(skeleton=True) \
    +(-n*grad(u)*v-n*grad(v)*u)* ds(skeleton=True)
a += aterms

f = LinearForm(fes)
fterms = rhs*v * dx \
     +alpha*order**2/h*bndc*v * ds(skeleton=True) \
     +(-n*grad(v)*bndc)* ds(skeleton=True)
f += fterms
with TaskManager():
    a.Assemble()
    f.Assemble()

For \(u_{h,f}\) a particular solution, we are looking for a solution \(u_h\in\IT^p(\Th) + u_{h,f}\) so that

\[\begin{equation} \label{eq:inhom} a_h(u_{h},v_{\IT}) =\ell(v_{\IT}) ~~ \forall~ v_{\IT}\in\IT^p(\Th). \end{equation}\]

After homogenization this means that we are looking for \(u_{\IT}\in\IT^p(\Th)\) that (uniquely) solves

\[\begin{equation} a_h(u_{\IT},v_{\IT}) =\ell(v_{\IT}) - a_h(u_{h,f},v_{\IT}) ~~ \forall~ v_{\IT}\in\IT^p(\Th). \end{equation}\]

This translates to the solution of the linear system

\[\begin{align} \bT^T\bA\bT \bu_{\IT} = \bT^T (\bl-\bA\bu_f). \end{align}\]

To compute a (local) particular solution, on an element \(K\in \mathcal{T}_h\) we assemble \((\bw_K)_{i}=(f,\calG(\be_i))=(f,\phi_i)\) and define \((\bu_{f})_K=\bW^\dagger_K \bw_K.\)

Lap = lambda u : sum(Trace(u.Operator('hesse')))
eps = 10**-7
u,v = fes.TnT()
op = Lap(u)*Lap(v)*dx
lop = -rhs*Lap(v)*dx
with TaskManager():
    emb = TrefftzEmbedding(op,lop,eps=eps)
PP = emb.GetEmbedding()
PPT = PP.CreateTranspose()
uf = emb.GetParticularSolution()
with TaskManager():
    TA = PPT@a.mat@PP
    Tgfu = TA.Inverse(inverse="sparsecholesky")*(PPT*(f.vec-a.mat*uf))
    gfu = GridFunction(fes)
    gfu.vec.data = PP*Tgfu + uf
error = sqrt(Integrate((gfu-exact)**2, mesh))
print("error =",error)

error = 7.234469621652273e-08

from ngsolve.webgui import Draw
Draw(BoundaryFromVolumeCF(gfu),mesh,'gfu')

BaseWebGuiScene

Using EmbeddedTrefftzFES

To use an \(\texttt{EmbeddedTrefftzFES}\) for a problem with a right-hand-side we first assemble the (Bi)linearForm over the Trefftz space, as in the homogeneous case.

etfes = EmbeddedTrefftzFES(emb)
a = BilinearForm(etfes)
a += aterms
f = LinearForm(etfes)
f += fterms
a.Assemble()
f.Assemble()

<ngsolve.comp.LinearForm at 0x7f43dc092bf0>

Computing the residuum

Then we compute the residuum \(\bT^T \bl-\bT^T \bA\bu_f\). We compute \(\bT^T \bA\bu_f\) directly, using a mixed BilinearForm with trial functions in the base space and test functions from the embedded Trefftz space.

res = f.vec.CreateVector()
af = BilinearForm(fes,etfes)
af += aterms
af.Apply(uf,res)
f.vec.data -= res

Now we are ready to solve and then compute the solution from \(u_{\IT} +u_{h,f}\).

tgfu = GridFunction(etfes)
tgfu.vec.data = a.mat.Inverse() * f.vec
gfu.vec.data = emb.Embed(tgfu.vec)
gfu.vec.data += uf

error = sqrt(Integrate((gfu-exact)**2, mesh))
print("error =",error)

error = 7.234469626924987e-08