Quasi-Trefftz DG: diffusion-advection-reaction

Consider the following diffusion-advection-reaction BVP

\[\begin{split}\begin{cases} \mathcal{L}u:=\text{div}(-\boldsymbol{K} \nabla u +\boldsymbol{\beta} u) +\sigma u = f &\text{in $\Omega$,}\\ u=g_{D} &\text{on $ \Gamma_{D}$,}\\ - \boldsymbol{K} \nabla u \cdot \mathbf{n} =g_{N} &\text{on $ \Gamma_{N}$.} \end{cases}\end{split}\]

with piecewise-smooth coefficients and right-hand side: $\boldsymbol{K}=\boldsymbol{K}^T$ diffusion matrix, $\boldsymbol{\beta}$ advection vector, $\sigma$ reaction scalar and $ f$ source term.

We want to solve this second-order elliptic boundary value problem using a quasi-Trefftz DG method.

[1]:

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

Constructing a polynomial quasi-Trefftz space

A polynomial quasi-Trefftz space for the diffusion-advection-reaction equation is given by the polynomials that locally belong to

\[\mathbb{Q\!T}^p_f(E):=\big\{ v\in \mathbb{P}^p(E) \mid D^{\boldsymbol{i}} \mathcal{L} v (\boldsymbol{x}^E)=D^{\boldsymbol{i}} f (\boldsymbol{x}^E)\quad \forall \boldsymbol{i}\in \mathbb{N}^d_0,\ |\boldsymbol{i}|\leq p-2\big\} \qquad \boldsymbol{x}^E \in E, \quad p\in \mathbb{N}.\]

We can construct it in NGSolve like so:

[2]:

h = 0.03;
order = 3;
mesh = Mesh(unit_square.GenerateMesh(maxh=h))
fes = trefftzfespace(mesh,order=order,eq="qtelliptic")

Using the eq key word one needs to tell trefftzfespace the operator for which to construct quasi-Trefftz functions.

We will consider the following manufactured example

[3]:

exact = sin(pi*(x+y))
K = CF((1+x+y,0,0,1+x+y),dims=(2,2))
beta = CF((1,0))
sigma = 3/(1+x+y)
rhs = -sum( (K*CF((exact.Diff(x),exact.Diff(y))))[i].Diff(var) for i,var in enumerate([x,y])) + beta*CF((exact.Diff(x),exact.Diff(y))) + sigma*exact

To set the coefficients for the construction of the quasi-Trefftz basis functions of $\mathbb{Q\!T}^p_0(\mathcal{T}_h)$ use

[4]:

with TaskManager():
    fes.SetCoeff(K,beta,sigma)

To set the right-hand side for the construction of the quasi-Trefftz particular approximate solution $u_{h,f}\in \mathbb{Q\!T}^p_f(\mathcal{T}_h)$ use

[5]:

with TaskManager():
    uf = fes.GetEWSolution(rhs)

A possbile DG method is:

\[\text{Find } u_{h} \in \mathbb{Q\!T}^p_f(\mathcal{T}_h) \text{ such that } \mathcal{A}_h^\mathrm{dar}(u_{h},v_h)=L_h(v_h) \quad \forall v_h \in \mathbb{Q\!T}^p_0(\mathcal{T}_h),\]

with the bilinear form

\[\begin{split}\begin{align} \mathcal{A}_h^\mathrm{dar}(u_h,v_h):=&\sum_{E\in \mathcal{T}_h} \int_E \Big(\Big(\boldsymbol{K}\nabla u_h -\boldsymbol{\beta} u_h\Big) \cdot \nabla v_h+\sigma u_h v_h\Big)\\&+ \sum_{F\in\mathcal{F}_h^{\mathrm I}\cup \mathcal{F}_h^{\mathrm D}} \int_F \Big(-\{\!\!\{\boldsymbol{K} \nabla u_h \}\!\!\} \cdot [\![v_h]\!]- [\![u_h]\!] \cdot \{\!\!\{\boldsymbol{K} \nabla v_h \}\!\!\} + \gamma\frac{K_F}{h_F} [\![u_h]\!]\cdot [\![v_h]\!]\Big)\\ &+ \int_{\mathcal{F}_h^{\mathrm I}} \Big( \{\!\!\{\boldsymbol{\beta} u_h\}\!\!\} \cdot [\![v_h]\!]+\frac12|\boldsymbol{\beta}\cdot \mathbf{n}_F|[\![u_h]\!]\cdot[\![v_h]\!]\Big)+ \int_{\mathcal{F}_h^+} u_h\boldsymbol{\beta} \cdot \mathbf{n} v_h, \end{align}\end{split}\]

and the linear form

\[L_h(v_h):=\sum_{E\in\mathcal{T}_h}\int_E f v_h - \int_{\mathcal{F}_h^\mathrm{N}} g_{\mathrm N} v_h + \int_{\mathcal{F}_h^{\mathrm D}} g_{\mathrm D} \Big(-\boldsymbol{K} \nabla v_h \cdot \mathbf{n} + \gamma\frac{K_F}{h_F} v_h%- \bbeta \cdot \mathbf{n} v_h \Big) -\int_{\mathcal{F}_h^-}g_D\boldsymbol{\beta}\cdot\mathbf{n} v_h.\]

[6]:

def dgell(fes,K,beta=CF((0,0)),sigma=0,Dbndc=0,Dbnd=".*",Nbndc=0,Nbnd="",rhs=0,uf=None,alpha=0):
    mesh = fes.mesh
    order = fes.globalorder

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

    alpha = 50*order**2/h if alpha==0 else alpha

    u = fes.TrialFunction()
    v = fes.TestFunction()
    jump = lambda u: (u-u.Other())*n
    mean_d = lambda u: 0.5 * K * (grad(u)+grad(u).Other())
    mean_B = lambda u: 0.5 * beta * (u+u.Other())

    a = BilinearForm(fes)
    a += K*grad(u)*grad(v) * dx \
        +alpha*jump(u)*jump(v) * dx(skeleton=True) \
        +(-mean_d(u)*jump(v)-mean_d(v)*jump(u)) * dx(skeleton=True) \
        +alpha*u*v * ds(skeleton=True,definedon=mesh.Boundaries(Dbnd)) \
        +(-K*grad(u)*n*v-K*grad(v)*n*u)* ds(skeleton=True,definedon=mesh.Boundaries(Dbnd))
    a += (-beta*u*grad(v) + sigma*u*v) * dx \
        + (mean_B(u) * jump(v) + 0.5*sqrt((beta*n)**2)*jump(u)*jump(v)) * dx(skeleton=True) \
        + beta*u*n*v * ds(skeleton=True,definedon=mesh.Boundaries(Nbnd))

    f = LinearForm(fes)
    f += Dbndc * (-K*grad(v)*n + alpha*v - beta*n*v) * ds(skeleton=True,definedon=mesh.Boundaries(Dbnd)) \
         - Nbndc * v * ds(skeleton=True,definedon=mesh.Boundaries(Nbnd)) \
         + rhs*v*dx
    if uf:
        f += -K*grad(uf)*grad(v) * dx \
            -alpha*jump(uf)*jump(v) * dx(skeleton=True) \
            -(-mean_d(uf)*jump(v)-mean_d(v)*jump(uf)) * dx(skeleton=True) \
            -alpha*uf*v * ds(skeleton=True,definedon=mesh.Boundaries(Dbnd)) \
            -(-K*grad(uf)*n*v-K*grad(v)*n*uf)* ds(skeleton=True,definedon=mesh.Boundaries(Dbnd))
        f += (beta*uf*grad(v) - sigma*uf*v) * dx \
            - (mean_B(uf) * jump(v) + 0.5*sqrt((beta*n)**2)*jump(uf)*jump(v)) * dx(skeleton=True) \
            - beta*uf*n*v * ds(skeleton=True,definedon=mesh.Boundaries(Nbnd))

    with TaskManager():
        a.Assemble()
        f.Assemble()

    return a,f

Consider the previous manufactured example and obtain the approximate solution gfu

[7]:

Dbndc = exact
Dbnd = '.*'
Nbndc = 0
Nbnd = ''

a,f = dgell(fes,K,beta,sigma,Dbndc=Dbndc,Dbnd=Dbnd,Nbnd=Nbnd,Nbndc=Nbndc,rhs=rhs,uf=uf)
gfu = GridFunction(fes)

with TaskManager():
    gfu.vec.data = a.mat.Inverse()*f.vec
print("quasi-Trefftz DG error:",sqrt(Integrate((gfu-exact)**2, mesh)))

quasi-Trefftz DG error: 0.00044383471915153015

We can plot the approximate solution like so:

[8]:

from ngsolve.webgui import Draw
Draw(gfu,mesh,"gfu")

[8]:

BaseWebGuiScene

Advection-dominated example

[9]:

nu = 10**(-3)
K = CF((nu,0,0,nu),dims=(2,2))
beta = CF((y+1,-x+2))
sigma = 0
rhs = 0
uf = None
Dbndc = IfPos(IfPos(x-10**-10,0,1) + IfPos(y-10**-10,0,1)*IfPos(x-1/3,0,1) - 0.5, 1, 0)
Nbndc = 0
Dbnd = "bottom|left"
Nbnd = "top|right"

[10]:

with TaskManager():
    fes = trefftzfespace(mesh,order=order,eq="qtelliptic")
    fes.SetCoeff(K,beta,sigma)

alpha = 100*nu/h
a,f = dgell(fes,K,beta,sigma,Dbndc=Dbndc,Dbnd=Dbnd,Nbnd=Nbnd,Nbndc=Nbndc,rhs=rhs,uf=uf,alpha=alpha)
gfu = GridFunction(fes)
with TaskManager():
    gfu.vec.data = a.mat.Inverse() * f.vec
Draw(gfu,mesh,"gfu",deformation=CF((x,y,gfu)))

[10]:

BaseWebGuiScene

Circular obstacle example

[11]:

h = 0.05
order = 3
wp = WorkPlane().RectangleC(2,2) \
    .Circle(0,0,0.25).Reverse()
geo = wp.Face()
mesh = Mesh(OCCGeometry(geo,dim=2).GenerateMesh(maxh=h))
mesh.Curve(order)
mesh.ngmesh.SetBCName(0,"top")
mesh.ngmesh.SetBCName(1,"bottom")
mesh.ngmesh.SetBCName(2,"left")
mesh.ngmesh.SetBCName(3,"right")
mesh.ngmesh.SetBCName(4,"innercircle")

[12]:

nu = 10**(-3)
K = CF((nu,0,0,nu),dims=(2,2))
beta = CF((1+1/16*(y**2-x**2)/((x**2+y**2)**2),-1/8*x*y/((x**2+y**2)**2)))
sigma = 0
Dbndc = IfPos(x+1/2,1,0)
Nbndc = 0
Dbnd="right|innercircle"
Nbnd="left|top|bottom"
rhs = 0
uf = None

with TaskManager():
    fes = trefftzfespace(mesh,order=order,eq="qtelliptic")
    fes.SetCoeff(K,beta,sigma)
    alpha = nu*50*order/h
    a,f = dgell(fes,K,beta,sigma,Dbndc=Dbndc,Dbnd=Dbnd,Nbnd=Nbnd,Nbndc=Nbndc,rhs=rhs,uf=uf,alpha=alpha)
gfu = GridFunction(fes)
with TaskManager():
    gfu.vec.data = a.mat.Inverse() * f.vec

Draw(gfu,mesh,"gfu")

[12]:

BaseWebGuiScene

L-shaped domain example

[13]:

wp = WorkPlane()
wp.MoveTo(0,0.5).LineTo(0.5,0.5).LineTo(0.5,0).LineTo(1,0).LineTo(1,1).LineTo(0,1).Close()
geo = wp.Face()
h = 0.05
mesh = Mesh(OCCGeometry(geo,dim=2).GenerateMesh(maxh=h))

[14]:

nu = 10**(-2)
K = CF((nu,0,0,nu),dims=(2,2))
beta = CF((-y,x))
sigma = 0
Dbndc = IfPos(y,0,1)
Dbnd = ".*"
rhs = 0
uf = None

order = 4
with TaskManager():
    fes = trefftzfespace(mesh,order=order,eq="qtelliptic")
    fes.SetCoeff(K,beta,sigma)
    alpha = nu*50*order/h
    a,f = dgell(fes,K,beta,sigma,Dbndc=Dbndc,Dbnd=Dbnd,Nbnd=Nbnd,Nbndc=Nbndc,rhs=rhs,uf=uf,alpha=alpha)

gfu = GridFunction(fes)
with TaskManager():
    gfu.vec.data = a.mat.Inverse() * f.vec
Draw(gfu,mesh,"gfu",deformation=CF((x,y,gfu)))

[14]:

BaseWebGuiScene