Embedded Trefftz-DG: Stokes problem

\(\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} \DeclareMathOperator{\Div}{\mathrm{div}}\)

The Stokes equations finds a velocity \(u\) and a pressure \(p\) such that

\[\begin{split}\begin{align} -\nu \Delta u + \nabla p & = f && \quad \textrm{in } \Omega, \\ -\Div u &=g && \quad \textrm{in } \Omega,\\ u &= 0 &&\quad \textrm{on } \partial \Omega, \end{align}\end{split}\]

[1]:

from ngsolve import *
from ngstrefftz import *
from netgen.occ import *
from ngsolve.webgui import Draw
SetNumThreads(4)
nu=1

The 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{-\nu \Delta u + \nabla p,w_u}_{0,h} + \inner{\Div u, w_p}_{0,h} \end{align}\]

with test functions

\[\begin{align} (w_u, w_p) \in [\mathbb P ^{k-2}(\mathcal T_h)]^d \times \mathbb P ^{k-1}(\mathcal T_h) \end{align}\]

We are missing a Hesse operator for \(u\) therefore, assuming \(\nu\) constant, we use the equivalent formulation

\[\begin{align} \inner{-\nu \Delta u}_{0,h} = \inner{\nu \nabla u, \nabla w_p}_{0,h} - \sum_{T\in \mathcal T_h} \int_{ \partial T} \nu \partial_n u w_u ds \end{align}\]

where we used the notation \(\partial_n w := \nabla w \cdot n\).

[2]:

def TrefftzStokesEmb(fes,rhs=CF((0,0))):
    mesh = fes.components[0].mesh
    k = fes.components[0].globalorder
    Vs = VectorL2(mesh, order=k - 2)
    Qs = L2(mesh, order=k - 1)
    test_fes = Vs * Qs

    u, p = fes.TrialFunction()
    wu, wp = test_fes.TestFunction()[0:2]

    n = specialcf.normal(mesh.dim)
    op = nu*InnerProduct( Grad(u),Grad(wu) ) * dx \
        - nu*InnerProduct(Grad(u)*n,wu) * dx(element_boundary=True) \
        + InnerProduct(grad(p),wu) * dx + div(u)*wp*dx

    lop = rhs * wu * dx(bonus_intorder=5)
    PP, uf = TrefftzEmbedding(op,fes,lop,eps=10**-6,test_fes=test_fes)
    return PP, uf

We can visualize the Trefftz basis functions via the embedding and drawing the \(\texttt{GridFunction}\) of the corresponding polynomial space

[3]:

def EmbeddedBasisFunctionsAsMultiDim(Emb, fes):
    gfshow = GridFunction(fes, multidim=0)
    gf = GridFunction(fes)
    coefvec = Emb.CreateRowVector()
    for i in range(Emb.width):
        coefvec[:] = 0
        coefvec[i] = 1
        gf.vec.data = Emb * coefvec
        gfshow.AddMultiDimComponent(gf.vec)
    return gfshow

[4]:

order = 3
mesh = Mesh(unit_square.GenerateMesh(maxh=1))

V = VectorL2(mesh, order=order, dgjumps=True)
Q = L2(mesh, order=order - 1, dgjumps=True)
fes = V * Q
PP,uf = TrefftzStokesEmb(fes)

gfshow = EmbeddedBasisFunctionsAsMultiDim(PP,fes)
Draw (gfshow.components[0], mesh, interpolate_multidim=False, animate=False, autoscale=True, vectors = {"grid_size" : 20}, settings={"subdivision":20})

[4]:

BaseWebGuiScene

Stokes-DG formulation

\(\newcommand{\Th}{{\mathcal{T}_h}} \newcommand{\Fh}{{\mathcal{F}_h}} \newcommand{\Fhi}{{\mathcal{F}_h^{\text{i}}}} \newcommand{\Fhb}{{\mathcal{F}_h^{\text{b}}}} \newcommand{\TT}{\mathbb{T}} \newcommand{\X}{\mathbb{X}} \renewcommand{\SS}{\mathbb{S}} \newcommand{\RR}{\mathbb{R}} \newcommand{\PP}{\mathbb{P}} \newcommand{\BDM}{\mathbb{B}\mathbb{D}\mathbb{M}} \newcommand{\RT}{\mathbb{R}\mathbb{T}} \newcommand{\PiN}[1]{\Pi_{\mathcal{N}}^{#1}} \newcommand{\cP}{\mathcal{P}} \DeclareMathOperator{\range}{range} \DeclareMathOperator{\dist}{dist} \DeclareMathOperator{\id}{id} \DeclareMathOperator{\Int}{Int} \DeclareMathOperator{\meas}{meas} \DeclareMathOperator{\Div}{\mathrm{div}} \DeclareMathOperator{\Grad}{\nabla} \DeclareMathOperator\Supp{Supp} \DeclareMathOperator\Span{span} \DeclareMathOperator\Rank{rank} \DeclareMathOperator\Ker{ker} \DeclareMathOperator\Range{range} \DeclareMathOperator\dt{\frac{\partial}{\partial t}} \newcommand\restr[2]{{\left.\kern-\nulldelimiterspace #1 \right|_{#2} }} \newcommand{\avg}[1]{\{\!\!\{ #1\}\!\!\}} \newcommand{\jmp}[1]{[\![#1]\!]} \newcommand{\jump}[1]{[\![#1]\!]}\) As a basis for the following we consider the established DG discretization of the Stokes problem: Find \((u_h, p_h) \in [\PP^k]^d \times \PP^{k-1}\setminus \RR\), s.t.

\[\begin{split}\begin{align} a_h(u_h,v_h) &+ b_h(v_h,p_h) = (f,v_h)_{\Th} && \forall ~ v_h \in [\PP^k]^d,\\ b_h(u_h,q_h) & = (g,q_h)_{\Th} && \forall ~ q_h \in \PP^{k-1}\setminus \RR, \end{align}\end{split}\]

with the bilinear forms

\[\begin{split}\begin{align*} a_h(u_h,v_h) &:= (\nu\! \Grad\! u_h,\Grad\! v_h)_\Th\! - (\avg{ \nu \partial_n u_h },\jmp{v_h})_\Fh\! - (\avg{\nu \partial_n v_h},\jmp{u_h})_\Fh \\ & \qquad \!+ \frac{\alpha \nu }{h} (\jmp{u_h},\jmp{v_h})_\Fh, \\ b_h(v_h,p_h) &:= -(\Div v_h, p_h)_{\Th} + (\jmp{v_h \cdot n}, \avg{p_h})_{\Fh}, \end{align*}\end{split}\]

where the interior penalty parameter \(\alpha = \mathcal{O}(k^2)\) is chosen sufficiently large enough.

[5]:

def StokesDG(fes, nu, rhs):
    k = fes.components[0].globalorder
    u, v = fes.TrialFunction()[0], fes.TestFunction()[0]
    p, q = fes.TrialFunction()[1], fes.TestFunction()[1]

    alpha = 10  # interior penalty param
    n = specialcf.normal(mesh.dim)
    h = specialcf.mesh_size

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

    a = BilinearForm(fes)
    a += nu * InnerProduct(grad(u), grad(v)) * dx
    a += nu * alpha * k**2 / h * jump_u * jump_v * dx(skeleton=True)
    a += nu * (-mean_dudn * jump_v - mean_dvdn * jump_u) * dx(skeleton=True)
    a += nu * alpha * k**2 / h * u * v * ds(skeleton=True)
    a += nu * (-grad(u) * n * v - grad(v) * n * u) * ds(skeleton=True)

    a += (-div(u) * q - div(v) * p) * dx
    a += (mean_p * jump_v + mean_q * jump_u) * dx(skeleton=True)
    a += (p * v * n + q * u * n) * ds(skeleton=True)
    a += 1e-8 * (p * q) * dx
    a.Assemble()

    f = LinearForm(fes)
    f += rhs * v * dx(bonus_intorder=5)
    f.Assemble()
    return a,f

We consider the following exact solution

[6]:

zeta = x**2 * (1-x)**2 * y**2 * (1-y)**2
uexact = CF((zeta.Diff(y), - zeta.Diff(x)))
graduexact = CF((uexact.Diff(x),uexact.Diff(y)),dims=(2,2)).trans
pexact = x**6+y**6 - 2/7

f1 = - nu*uexact[0].Diff(x).Diff(x) - nu*uexact[0].Diff(y).Diff(y) + pexact.Diff(x)
f2 = - nu*uexact[1].Diff(x).Diff(x) - nu*uexact[1].Diff(y).Diff(y) + pexact.Diff(y)
rhs = CF((f1,f2))

[7]:

order = 3
mesh = Mesh(unit_square.GenerateMesh(maxh=.1))

V = VectorL2(mesh, order=order, dgjumps=True)
Q = L2(mesh, order=order - 1, dgjumps=True)
fes = V * Q
a,f = StokesDG(fes, 1, rhs)

PP, uf = TrefftzStokesEmb(fes, rhs)
PPT = PP.CreateTranspose()
TA = PPT @ a.mat @ PP
Tgfu = TA.Inverse(inverse="pardiso") * (PPT * (f.vec - a.mat * uf))
gfu = GridFunction(fes)
gfu.vec.data = PP * Tgfu + uf

[8]:

uh = gfu.components[0]
Draw(uh,mesh,"u")
print(sqrt(Integrate(InnerProduct(uh-uexact,uh-uexact),mesh)))

9.39961407472032e-07