Embedded Trefftz method

\(\newcommand{\Th}{\mathcal{T}_h} \newcommand{\Vhp}{V^p(\Th)} \newcommand{\bT}{\mathbf{T}} \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}\)

Let \(\{\phi_i\}\) be the set of basis functions of \(\Vhp\), \(N =\operatorname{dim}(\Vhp)\) and \(\calG: \mathbb{R}^n \to \Vhp\) be the Galerkin isomorphism, \(\calG(\bx) = \sum_{i=1}^N \bx_i \phi_i\). With \(\be_i,~i=1,..,N\) the canonical unit vectors in \(\mathbb{R}^n\), so that \(\calG(\be_i) = \phi_i\), we define the following matrices and vector for \(i,j = 1,\dots, N\)

\[\begin{split}\begin{align} (\bA)_{ij}&=a_h(\calG(\be_j),\calG(\be_i))=a_h(\phi_j,\phi_i), \qquad (\bl)_i = \ell(\calG(\be_i)) = \ell(\phi_i) \\ (\bW)_{ij}&=\inner{\calL\phi_j,\calL\phi_i}_{0,h}, \label{def:W} \end{align}\end{split}\]

where \(\inner{\cdot,\cdot}_{0,h}=\sum_{K\in\Th}\inner{\cdot,\cdot}_K\) is the element-wise \(L^2\)-inner product. We are interested in the kernel of \(\calL\) (in an element-wise and pointwise sense) in \(\Vhp\) as this is the part where the Trefftz DG method operates. We note that

\[\begin{equation} \Ker(\calL) = \calG(\Ker(\bW)) \end{equation}\]

and hence are looking for a basis of \(\Ker(\bW)\) which . We can determine the kernel of \(\bW\) element-wise, for \(\bW_K\) collect a set of orthogonal basis vectors in a matrix \(\bT_K\in\mathbb{R}^{N_{\!K}\!\times\! M_{\!K}}\) so that \(\ker(\bW_K) = \bT_K \cdot \mathbb{R}^{M_K}\).

We can then implicitly solve using Trefftz test and trial functions by solving:

Find \(\bu_\IT\) so that

\[\begin{equation} \label{eq:trefftzlinearsys} \bT^T\bA\bT ~ \bu_\IT = \bT^T \bl. \end{equation}\]

The solution in the full polynomial space is then given by \(\bu=\bT\bu_\IT\)

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

Laplace problem

We are looking to solve

\[\begin{split}\begin{align*} \begin{cases} -\Delta u = 0 &\text{ in } \Omega, \\ u=g &\text{ on } \partial \Omega, \end{cases} \end{align*}\end{split}\]

with the embedded Trefftz method, now $ \calL `= :nbsphinx-math:Delta`$.

exact = exp(x+y)*sin(sqrt(2)*z)
order = 4
mesh = Mesh(unit_cube.GenerateMesh(maxh=.2))
fes = L2(mesh, order=order,  dgjumps=True)

A suiteable DG method is given by

\[\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}} \begin{align} a_h(u,v) &= \int_\dom \nabla u\nabla v\ dV -\int_{\Fh^\text{int}}\left(\avg{\nabla u}\jump{v}+\avg{\nabla v}\jump{u} - \frac{\alpha p^2}{h}\jump{u}\jump{v} \right) dS \\ &\qquad -\int_{\Fh^\text{bnd}}\left(\nx\cdot\nabla u v+\nx\cdot\nabla v u-\frac{\alpha p^2}{h} u v \right) dS\\ \ell(v) &= \int_{\Fh^\text{bnd}}\left(\frac{\alpha p^2}{h} gv -\nx\cdot\nabla vg\right) dS. \end{align}\end{split}\]

def lapdg(fes,bndc):
    alpha = 4
    n = specialcf.normal(mesh.dim)
    h = specialcf.mesh_size
    u = fes.TrialFunction()
    v = fes.TestFunction()

    jump = lambda u: u-u.Other()
    mean_dn = lambda u: 0.5*n * (grad(u)+grad(u).Other())

    a = BilinearForm(fes,symmetric=True)
    a += grad(u)*grad(v) * dx \
        +alpha*order**2/h*jump(u)*jump(v) * dx(skeleton=True) \
        +(-mean_dn(u)*jump(v)-mean_dn(v)*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.Assemble()

    f = LinearForm(fes)
    f += alpha*order**2/h*bndc*v * ds(skeleton=True) \
         +(-n*grad(v)*bndc)* ds(skeleton=True)
    f.Assemble()
    return a,f

a,f = lapdg(fes,exact)
with TaskManager():
    a.Assemble()
    f.Assemble()

Trefftz embedding

We use the function \(\texttt{TrefftzEmbedding}\) and pass the operator

\[\newcommand{\Vhp}{V^p(\Th)} \newcommand{\bT}{\mathbf{T}} \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} \inner{\calL \cdot,\calL \cdot}_{0,h} \end{align}\]

When numerically computing the kernel of a matrix, we use a SVD. The diagonal matrix of the SVD would have \(M_K\) zeros assuming exact arithmetics. Due to inexact computer arithmetics this will not be the case exactly and hence we use a truncation parameter \(\varepsilon>0\) to determine which values are considered as (numerical) zeros.

eps=10**-9
Lap = lambda u : sum(Trace(u.Operator('hesse')))
u,v = fes.TnT()
op = Lap(u)*Lap(v)*dx
with TaskManager():
    PP = TrefftzEmbedding(op,fes,eps)
PPT = PP.CreateTranspose()
with TaskManager():
    TA = PPT@a.mat@PP
    TU = TA.Inverse(inverse='sparsecholesky')*(PPT*f.vec)
    tgfu = GridFunction(fes)
    tgfu.vec.data = PP*TU

error = sqrt(Integrate((tgfu-exact)**2, mesh))
print("embedded Trefftz error ",error)

embedded Trefftz error  6.901461102254937e-07

Trefftz embedding type 2

We can also construct the embedding by finding the kernel of

\[\begin{align} (\bW)_{ij}&=\inner{\calL\phi_j,\tilde\phi_i}_{0,h},\qquad \forall \tilde\phi_i\in V^{p-2}(\Th) \end{align}\]

eps=10**-9
Lap = lambda u : sum(Trace(u.Operator('hesse')))
u = fes.TrialFunction()
fes2 = L2(mesh, order=order-2,  dgjumps=True)
v = fes2.TestFunction()
op = Lap(u)*v*dx
with TaskManager():
    PP = TrefftzEmbedding(op,fes,test_fes=fes2)
PPT = PP.CreateTranspose()
with TaskManager():
    TA = PPT@a.mat@PP
    TU = TA.Inverse(inverse='sparsecholesky')*(PPT*f.vec)
    tgfu = GridFunction(fes)
    tgfu.vec.data = PP*TU

error = sqrt(Integrate((tgfu-exact)**2, mesh))
print("embedded Trefftz error ",error)

embedded Trefftz error  6.901461109735476e-07

Trefftz embedding via FESpace

The following improves the runtime of the embedded Trefftz method by avoiding the full assembly of the larger system matrix. Instead we directly assemble the matrix \(\bT^T\bA\bT\) using the local embedding matrices.

We use EmbeddedTrefftzFES to construct an embedded space from an underlying L2 space.

with TaskManager():
    fes = L2(mesh, order=order, dgjumps=True)
    u,v = fes.TnT()
    op = Lap(u)*Lap(v)*dx

    etfes = EmbeddedTrefftzFES(fes)
    etfes.SetOp(op,eps=eps)

    a,f = lapdg(etfes,exact)
    a.Assemble()
    f.Assemble()

    inv = a.mat.Inverse(inverse="sparsecholesky")
    gfu = GridFunction(etfes)
    gfu.vec.data = inv * f.vec

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

embedded Trefftz error  nan