Embedded Trefftz-DG: Navier-Stokes

For the stationary incompressible Navier-Stokes equations,

\[-\nu\Delta u + (u\cdot\nabla)u + \nabla p = f, \qquad \nabla\cdot u = 0,\]

the nonlinear part is the convective velocity in \((u\cdot\nabla)u\). The embedded Trefftz construction enters after linearization: freeze one velocity field and solve an Oseen problem in a Trefftz space adapted to that frozen operator.

We use the following DG formulation

\[\begin{split}\newcommand{\jump}[1]{\left[\!\left[#1\right]\!\right]} \newcommand{\avg}[1]{\left\{\!\left\{#1\right\}\!\right\}} \begin{aligned} a_h^w((u,p),(v,q)) &= \nu\int_{\Omega}\nabla u:\nabla v\,\mathrm{d}x - \int_{\Omega} \left( q\,\nabla\cdot u + p\,\nabla\cdot v \right) \,\mathrm{d}x \\ &\quad + \int_{\Omega} \left( (w\cdot \nabla u)\cdot v + \frac12(\nabla\cdot w)\,u\cdot v \right) \,\mathrm{d}x \\ &\quad - \int_{\mathcal F_h^{\mathrm{int}}} \left( \nu\avg{\nabla u}n\cdot\jump{v} + \nu\avg{\nabla v}n\cdot\jump{u} - \nu\frac{\alpha k^2}{h} \jump{u}\cdot\jump{v} \right) \,\mathrm{d}s \\ &\quad - \int_{\mathcal F_h^{\mathrm{int}}} \left( \avg{p}\,n\cdot\jump{v} - \avg{q}\,n\cdot\jump{u} \right) \,\mathrm{d}s \\ &\quad + \int_{\mathcal F_h^{\mathrm{int}}} \left( (\avg{w}\cdot n) \jump{u}\cdot\avg{v} + \frac12(\jump{w}\cdot n) \avg{u\cdot v} \right) \,\mathrm{d}s \\ &\quad +\text{ suitable boundary terms}. \end{aligned}\end{split}\]

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

def VLap(u, dim):
    hesse = u.Operator("hesse")
    return CF(tuple(sum(hesse[i, j * (dim + 1)] for j in range(dim))for i in range(dim)))

def SolveOseen(mesh,k,nu,rhs=None,wn=None,trefftz=True,ubnd=None,bndname=".*",inletname=None,alpha=50):
    if hasattr(wn, "Operators") and "div" in wn.Operators():
        divwn = div(wn)
    elif wn is not None:
        if mesh.dim == 2:
            divwn = wn[0].Diff(x) + wn[1].Diff(y)
        elif mesh.dim == 3:
            divwn = wn[0].Diff(x) + wn[1].Diff(y) + wn[2].Diff(z)

    bint = 16
    if inletname is None:
        inletname = bndname
    stab = 1e-7

    V = VectorL2(mesh, order=k, dgjumps=True)
    Q = L2(mesh, order=k - 1, dgjumps=True)
    basefes = V * Q

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

    if trefftz:
        Vs = VectorL2(mesh, order=k - 2)
        Qs = L2(mesh, order=k - 1)
        test_fes = Vs * Qs
        u, p = basefes.TrialFunction()[0:2]
        wu, wp = test_fes.TestFunction()[0:2]
        oseen_residual = -nu * VLap(u, mesh.dim) + grad(p)
        if wn is not None:
            oseen_residual += grad(u) * wn
        op = oseen_residual * wu * dx(bonus_intorder=bint)
        op += div(u) * wp * dx(bonus_intorder=bint)

        emb = TrefftzEmbedding(op)
        fes = EmbeddedTrefftzFES(emb)
        if rhs is not None:
            upf = GridFunction(basefes)
            lop = rhs * wu * dx(bonus_intorder=bint)
            upf.vec.data = emb.GetParticularSolution(lop)
    else:
        fes = basefes

    u, v = fes.TrialFunction()[0], fes.TestFunction()[0]
    p, q = fes.TrialFunction()[1], fes.TestFunction()[1]

    def jump(v):
        return v - v.Other()

    def mean(v):
        return 0.5 * (v + v.Other())

    def mean2(u, v):
        return 0.5 * (u * v + u.Other() * v.Other())

    a = BilinearForm(fes)
    def dgterms(u,p,v,q):
        a = nu * InnerProduct(grad(u), grad(v)) * dx
        a += nu * alpha * k**2 / h * jump(u) * jump(v) * dx(skeleton=True)
        a += nu * (-mean(grad(u)) * n * jump(v) - mean(grad(v)) * n * jump(u)) * dx(skeleton=True)
        a += nu * alpha * k**2 / h * u * v * ds(skeleton=True, definedon=mesh.Boundaries(bndname))
        a += nu * (-grad(u) * n * v - grad(v) * n * u) * ds(skeleton=True, definedon=mesh.Boundaries(bndname))
        a += (mean(p) * n * jump(v) + mean(q) * n * jump(u)) * dx(skeleton=True)
        a += (p * v * n + q * u * n) * ds(skeleton=True, definedon=mesh.Boundaries(bndname))
        a += (-div(u) * q - div(v) * p) * dx
        a += -stab * p * q * dx
        if wn is not None:
            a += (grad(u) * wn) * v * dx(bonus_intorder=bint)
            a += -mean(wn) * n * jump(u) * mean(v) * dx(skeleton=True, bonus_intorder=bint)
            a += 0.5 * divwn * u * v * dx(bonus_intorder=bint)
            a += -0.5 * jump(wn) * n * mean2(u, v) * dx(skeleton=True, bonus_intorder=bint)
            a += -0.5 * wn * n * u * v * ds(skeleton=True,bonus_intorder=bint,definedon=mesh.Boundaries(bndname),)
        return a
    a += dgterms(u,p,v,q)

    f = LinearForm(fes)
    if rhs is not None:
        f += rhs * v * dx(bonus_intorder=bint)
    if ubnd is not None:
        f += nu * alpha * k**2 / h * ubnd * v * ds(skeleton=True,definedon=mesh.Boundaries(inletname),bonus_intorder=bint,)
        f += nu * (-grad(v) * n * ubnd) * ds(skeleton=True,definedon=mesh.Boundaries(inletname),bonus_intorder=bint,)
        f += q * ubnd * n * ds(skeleton=True,definedon=mesh.Boundaries(inletname),bonus_intorder=bint,)
        if wn is not None:
            f += -0.5 * wn * n * ubnd * v * ds(skeleton=True,definedon=mesh.Boundaries(inletname),bonus_intorder=bint,)
    if trefftz and rhs is not None:
        mupf = GridFunction(basefes)
        mupf.vec.data = -1*upf.vec
        muf = mupf.components[0]
        mpf = mupf.components[1]
        f += dgterms(muf, mpf, v, q)


    a.Assemble()
    f.Assemble()

    gfu = GridFunction(fes)
    gfu.vec.data = a.mat.Inverse() * f.vec

    if trefftz:
        polygfu = GridFunction(basefes)
        polygfu.vec.data = emb.Embed(gfu.vec)
        gfu = GridFunction(basefes)
        if rhs is not None:
            gfu.vec.data = polygfu.vec + upf.vec
        else:
            gfu = polygfu

    uh, ph = gfu.components[0:2]

    return uh, ph, fes.ndof

Picard Linearization

Given a velocity \(w\), the Oseen operator is

\[\mathcal O_w(u,p) = -\nu\Delta u + (w\cdot\nabla)u + \nabla p, \qquad \nabla\cdot u=0.\]

The Trefftz embedding uses as local operator.

Picard iteration freezes the convection field at the previous velocity:

\[\begin{split}\begin{aligned} \mathcal O_{u^{m}}(u^{m+1},p^{m+1}) &= f,\\ \nabla\cdot u^{m+1} &= 0. \end{aligned}\end{split}\]

In each loop we solve the discrete problem with

\[(u\cdot\nabla)u \quad\leadsto\quad (u^m\cdot\nabla)u^{m+1}.\]

The embedded Trefftz space is updated in each iteration step.

wn = None
for step in range(maxiter):
    uh, ph = solve_oseen(mesh, nu, wn=wn, trefftz=True)
    if wn is not None and norm(uh - wn) < tol:
        break
    wn = uh

The next Oseen solve, and therefore the next embedded Trefftz space, is rebuilt for this new convection field.

Oseen Trefftz Space

On one element \(K\), start with the discontinuous polynomial trial space

\[V_h^k(K) = [\mathbb P^k(K)]^d \times \mathbb P^{k-1}(K)\]

for velocity and pressure. The residual is tested in

\[Q_h^k(K) = [\mathbb P^{k-2}(K)]^d \times \mathbb P^{k-1}(K).\]

For fixed convection \(w\), the embedded Oseen Trefftz space is

\[\begin{split}\mathbb T_{w,h}^k(K) = \left\{ (u_h,p_h)\in V_h^k(K): \begin{array}{l} \left(-\nu\Delta u + (w\cdot\nabla)u + \nabla p, v\right)_K=0\\ (\nabla\cdot u,q)_K =0 \end{array} \quad\forall (v_h,q_h)\in Q_h^k(K) \right\},\end{split}\]

mesh = Mesh(unit_square.GenerateMesh(maxh=0.45))
k = 3
nu = 1.0
wn = CF((1, 0))

V = VectorL2(mesh, order=k, dgjumps=True)
Q = L2(mesh, order=k - 1, dgjumps=True)
base = V * Q

R = VectorL2(mesh, order=k - 2) * L2(mesh, order=k - 1)
u, p = base.TrialFunction()[0:2]
v, q = R.TestFunction()[0:2]

op = (-nu * VLap(u, mesh.dim) + grad(p) + grad(u) * wn) * v * dx(bonus_intorder=16)
op += div(u) * q * dx(bonus_intorder=16)

emb = TrefftzEmbedding(op)
tfes = EmbeddedTrefftzFES(emb)

print(f"full Oseen polynomial dofs: {base.ndof}")
print(f"embedded Oseen Trefftz dofs: {tfes.ndof}")

full Oseen polynomial dofs: 156
embedded Oseen Trefftz dofs: 84

Schäfer-Turek Benchmark

We now run the DFG 2D-1 Schäfer-Turek flow in the channel with a circular obstacle,

\[\Omega=(0,2.2)\times(0,0.41)\setminus B_{0.05}(0.2,0.2).\]

The inflow profile is parabolic,

\[u_{\mathrm{in}}(y) = \left( \frac{4U_{\max}y(0.41-y)}{0.41^2}, 0 \right), \qquad U_{\max}=0.3, \qquad \nu=10^{-3}.\]

def solve_ns(mesh,k=4,nu=0.001,ubnd=None,bndname="inlet|wall|cyl",inletname="inlet",maxiter=100,tol=1e-8,alpha=50,trefftz=True):
    wn = CF((0,0))
    uh_old = wn

    for step in range(1, maxiter + 1):
        uh, ph, ndof = SolveOseen(mesh,k=k,nu=nu,wn=wn,trefftz=trefftz,ubnd=ubnd,bndname=bndname,inletname=inletname,alpha=alpha,)
        update = sqrt(Integrate(InnerProduct(uh - uh_old, uh - uh_old), mesh))
        print(f"Step: {step}, residual: {update}")
        if update < tol:
            return uh, ph, ndof, step

        uh_old = uh
        wn = uh

    return uh, ph, ndof, maxiter

ubnd = CF((4 * .3 * y * (0.41 - y) / (0.41 * 0.41), 0))
shape = Rectangle(2.2, 0.41).Circle(0.2, 0.2, 0.05).Reverse().Face()
shape.edges.name = "cyl"
shape.edges.Min(X).name = "inlet"
shape.edges.Max(X).name = "outlet"
shape.edges.Min(Y).name = shape.edges.Max(Y).name = "wall"
mesh = Mesh(OCCGeometry(shape, dim=2).GenerateMesh(maxh=0.05)).Curve(4)

with TaskManager():
    uh, ph, ndof, steps = solve_ns(mesh,k=4,ubnd=ubnd,trefftz=True)

Step: 1, residual: 0.2100373796626424
Step: 2, residual: 0.031504370815592556
Step: 3, residual: 0.013030499750240468
Step: 4, residual: 0.004684564927230454
Step: 5, residual: 0.0017049626177180696
Step: 6, residual: 0.0005757952588127792
Step: 7, residual: 0.00017496011231973453
Step: 8, residual: 5.029517045522804e-05
Step: 9, residual: 1.4265290644866134e-05
Step: 10, residual: 4.105605887995219e-06
Step: 11, residual: 1.2290775430418004e-06
Step: 12, residual: 3.8070306855857596e-07
Step: 13, residual: 1.1937128902954822e-07
Step: 14, residual: 3.753694679165403e-08
Step: 15, residual: 1.1790913454175637e-08
Step: 16, residual: 3.6753806689757763e-09

Draw(Norm(uh), mesh, "speed");
Draw(ph - Integrate(ph, mesh) / Integrate(1, mesh), mesh, "pressure");