This page was generated from unit-3.3-scalardg/scalardg.ipynb.

3.3 Discontinuous Galerkin Discretizations for linear transport

We are solving the scalar linear transport problem

Find \(u: [0,T] \to V_D := \{ u \in L^2(\Omega), b \cdot \nabla u \in L^2(\Omega), u|_{\Gamma_{in}} = u_D\}\), s.t. \begin{equation} \int_{\Omega} \partial_t u v + b \cdot \nabla u v = \int_{\Omega} f v \qquad \forall v \in V_0 = \{ u \in L^2(\Omega), b \cdot \nabla u \in L^2(\Omega), u|_{\Gamma_{in}} = 0\} \end{equation}

from ngsolve import *
import netgen.gui

As a first example, we consider the unit square \((0,1)^2\) and the advection velocity \(b = (1,2)\). Accordingly the inflow boundary is \(\Gamma_{in} = \{ x \cdot y = 0\}\).

from netgen.geom2d import SplineGeometry
geo = SplineGeometry()
geo.AddRectangle( (0, 0), (1, 1),
                 bcs = ("bottom", "right", "top", "left"))
mesh = Mesh( geo.GenerateMesh(maxh=0.125))
Draw(mesh)
order = 4

Convection and boundary values:

• To make cases with CoefficientFunctions we use the IfPos-CoefficientFunction. IfPos(a,b,c): a decides on the evaluation. If a is positive b is evaluated, otherwise c.

from math import pi
b = CoefficientFunction((1+sin(4*pi*y),2))
Draw(b,mesh,"wind")
from ngsolve.internal import visoptions; visoptions.scalfunction = "wind:0"
ubnd = IfPos(x-0.125,IfPos(0.625-x,1+cos(8*pi*x),0),0)

We consider an Upwind DG discretization (in space):

Find \(u: [0,T] \to V_h := \bigoplus_{T\in\mathcal{T}_h} \mathcal{P}^k(T)\) so that

\[\sum_{T} \int_T \partial_t u v - b \cdot \nabla v u + \int_{\partial T} b_n \hat{u} v = \int_{\Omega} f v, \quad \forall v \in V_h.\]

Here \(\hat{u}\) is the Upwind flux, i.e. \(\hat{u} = u\) on the outflow boundary part \(\partial T_{out} = \{ x\in \partial T \mid b(x) \cdot n_T(x) \geq 0 \}\) of \(T\) and \(\hat{u} = u^{other}\) else, with \(u^{other}\) the value from the neighboring element.

There is quite a difference in the computational costs (compared to a standard DG formulation) depending on the question if the solution of linear systems is involved or only operator evaluations (explicit method).

• We treat the explicit case here only

• and refer to unit 2.8 for solving linear systems with DG formulations.

Explicit time stepping with a DG formulation

Explicit Euler:

\[\sum_{T} \int_T u^{n+1} v = \sum_{T} \int_T u^{n} v + \Delta t \sum_{T} \left\{ \int_T b \cdot \nabla v u + \int_{\partial T} b_n \hat{u} v \right\} + \Delta t \int_{\Omega} f v, \quad \forall v \in V_h,\]

\[M u^{n+1} = M u^{n} - \Delta t C u^n + \Delta t f\]

In our first example we set \(u_0 = f = 0\).

Computing convection applications \(C u^n\)

• We can define the bilinear form without setting up a matrix with storage. (nonassemble=True)

• A BilinearForm is allowed to be nonlinear in the 1st argument (for operator applications)

VT = L2(mesh,order=order)
u,v = VT.TnT()
c = BilinearForm(VT, nonassemble=True)

• To access the neighbor function we can use u.Other(), cf. unit 2.8.

• To incorporate boundary conditions (\(\hat{u}\) on inflow boundaries), we can use the argument bnd of .Other() (not possible in unit 2.8 )

• If there is no neighbor element (boundary!) the CoefficientFunction bnd is evaluated.

n = specialcf.normal(mesh.dim)
upw_flux = b*n * IfPos(b*n, u, u.Other(bnd=ubnd))
dS = dx(element_boundary=True)
c += - b * grad(v) * u * dx + upw_flux * v * dS

Bilinearform for operator applications

Due to the nonassemble=True in the constructor of the BilinearForm we can use c.mat * x to realize the operator application without assembling a matrix first:

gfu = GridFunction(VT);
Draw(gfu,mesh,"u",sd=4,min=-0.1,max=2.1,autoscale=False)
res = gfu.vec.CreateVector()

#Operator application (equivalent to assemble and mult but faster)
res.data = c.mat * gfu.vec

# this does not work:
c2 = BilinearForm(VT)
res.data = c2.mat * gfu.vec

---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
<ipython-input-1-05fb163505bf> in <module>
      1 # this does not work:
      2 c2 = BilinearForm(VT)
----> 3 res.data = c2.mat * gfu.vec

TypeError: matrix not ready - assemble bilinearform first

Solving mass matrix problems (see also unit-2.11)

• Need to invert the mass matrix

• For DG methods the mass matrix is block diagonal, often even diagonal.

• FESpace offers (if available):

  • Mass(rho): mass matrix as an operator (not a sparse matrix)

  • Mass(rho).Inverse(): inverse mass matrix as an operator (not a sparse matrix)

invm = VT.Mass(1).Inverse()
print(invm)

Print base-matrix

The (simple) time loop

%%time
#better: %%timeit
t = 0; gfu.vec[:] = 0
dt = 2e-4
tend = 0.6

while t < tend-0.5*dt:
    res.data = invm @ c.mat * gfu.vec
    gfu.vec.data -= dt * res
    t += dt
    Redraw(blocking=True)

CPU times: user 5.72 s, sys: 7.33 ms, total: 5.73 s
Wall time: 5.73 s

Skeleton formulation (cf. unit 2.8)

So far we considered the DG formulation with integrals on the boundary of each element. Instead one could formulate the problem in terms of facet integrals where every facet appears only once.

Then, the corresponding formulation is:

Find \(u: [0,T] \to V_h := \bigoplus_{T\in\mathcal{T}_h} \mathcal{P}^k(T)\) so that

\[\begin{split}\sum_{T} \int_T \partial_t u v - b \cdot \nabla v u + \sum_{F\in\mathcal{F}^{inner}} \int_{F} b_n \hat{u} (v - v^{neighbor}) + \\ \sum_{F\in\mathcal{F}^{bound}} \int_{F} b_n \hat{u} v = \int_{\Omega} f v , \quad \forall v \in V_h.\end{split}\]

Here \(\mathcal{F}^{inner}\) is the set of interior facets and \(\mathcal{F}^{bound}\) is the set of boundary facets.

integral types:

The facet integrals are divided into inner facets and boundary facets. To obtain these integrals we combine the DifferentialSymbol dx or ds with skeleton=True.

dskel_inner  = dx(skeleton=True)
dskel_bound  = ds(skeleton=True)
# or:
dskel_inflow = ds(skeleton=True, definedon=mesh.Boundaries("left|bottom"))
dskel_outflow = ds(skeleton=True, definedon=mesh.Boundaries("right|top"))

Skeleton formulation of \(C\)

c = BilinearForm(VT,nonassemble=True)
c += -b * grad(v) * u * dx
c += upw_flux * (v-v.Other()) * dskel_inner

c += b*n * IfPos(b*n,u,ubnd) * v * dskel_bound
#alternatively (if you know in/out beforehand):
#c += b*n * ubnd * v * dskel_inflow
#c += b*n * u * v * dskel_outflow

The (simple) time loop

%%time
t = 0; gfu.vec[:]=0
while t < tend-0.5*dt:
    gfu.vec.data -= dt * invm @ c.mat * gfu.vec
    t += dt
    Redraw(blocking=True)

CPU times: user 4.42 s, sys: 45.2 ms, total: 4.46 s
Wall time: 4.4 s

TraceOperator (Purpose: avoid neighbor access Other())

The trace operator maps from an FESpace to a compatible FacetFESpace by evaluating the traces.

The stored value on a facet is the sum of the traces of the aligned elements.

No need to access the neighbor (.Other()) later on.

VF = FacetFESpace(mesh, order=order)
V = FESpace( [VT, VF] )

trace = VT.TraceOperator(VF, False)
print("The space VT has {} dofs.".format(VT.ndof))
print("The trace space VF has {} dofs.".format(VF.ndof))
print("trace represents an {} x {} operator".format(trace.height, trace.width))

The space VT has 2160 dofs.
The trace space VF has 1160 dofs.
trace represents an 1160 x 2160 operator

gfuF = GridFunction(VF)
gfuF.vec.data = trace * gfu.vec

To make use of the trace operator, we split the previous (element_boundary)formulation into three parts:

\[\sum_{T} \int_T - b \cdot \nabla v u + \sum_{T} \int_{\partial T} b_n \hat{u}^* v + \sum_{F\in\mathcal{F}^{inflow}} \int_{F} b_n u^{inflow} v \quad \text{ with } \quad \hat{u}^* = f(u^{left}+u^{right},u)\]

• part 1: element local volume integrals (no interaction with neighbor elements) (on VT)

• part 2: couplings on the facets (on VT \(\times\) V)

• part 3: (rhs + ) boundary terms (on VT)

uT, vT = VT.TnT()
c1 = BilinearForm(space=VT, nonassemble=True)                     # part 1
c1 += -uT * b * grad(vT) * dx
uT,uF = V.TrialFunction()
c2 = BilinearForm(trialspace=V, testspace=VT, nonassemble=True)   # part 2
# here uf-u = u_me + u_other - u_me is the neighbor value (0 at bnd)
c2 += b*n * IfPos(b*n, uT, uF-uT) * vT * dx(element_boundary=True)
rhs = LinearForm(VT)                                              # part 3
rhs += b*n * IfPos(b*n, 0, ubnd) * vT * dskel_bound # dskel_inflow
rhs.Assemble()

<ngsolve.comp.LinearForm at 0x7f5854d77a70>

Note that c1 and c2 act on different spaces. To make them compatible, we combine: * TraceOperator: VT \(\to\) VF * Embeddings: VT \(\to\) VT \(\times\) VF and VF \(\to\) VT \(\times\) VF (cf. unit 2.10)

The Embeddings are (with \(n_F = \operatorname{dim}(V_F)\), \(n_T = \operatorname{dim}(V_T)\), \(n_V=n_F+n_T\))

\[\begin{split}\texttt{embV}: V_T \to V_T \times V_F, \qquad \texttt{embV} = \underbrace{\left( \begin{array}{c} I \\ 0 \end{array} \right)}_{\texttt{embT} \in \mathbb{R}^{n_V \times n_T}} + \underbrace{\left( \begin{array}{c} 0 \\ I \end{array} \right)}_{\texttt{embF} \in \mathbb{R}^{n_V \times n_F}} \underbrace{\operatorname{tr}}_{\texttt{trace}\in \mathbb{R}^{n_F \times n_T}}\end{split}\]

embT = Embedding(V.ndof, V.Range(0))
embF = Embedding(V.ndof, V.Range(1))
embV = embT + embF @ trace

The operators can now be combined easily. Recall: * c1: VT \(\to\) VT’ \(\quad\) \(\bullet\) c2: VT \(\times\) VF \(\to\) VT’ \(\quad\) \(\bullet\) \(\Rightarrow\) c: VT \(\to\) VT’

c = c1.mat + c2.mat @ embV

The time loop:

%%time
t = 0; gfu.vec[:] = 0
invm = VT.Mass(1).Inverse()
r = gfu.vec.CreateVector(); r.data = invm * rhs.vec # <- the boundary data
while t < tend-0.5*dt:
    gfu.vec.data -= dt * (invm @ c * gfu.vec + r)
    t += dt
    Redraw(blocking=True)

CPU times: user 4.64 s, sys: 62.3 ms, total: 4.71 s
Wall time: 4.64 s

Careful with one-liners!

The following will not work correctly: (gfu appears left and right: order of ops is important!)

%%time
t = 0; gfu.vec[:] = 0
while t < tend-0.5*dt:
    gfu.vec.data -= dt * (r + invm @ c * gfu.vec)
    t += dt;
    Redraw(blocking=True)

CPU times: user 4.51 s, sys: 62.4 ms, total: 4.58 s
Wall time: 4.51 s

In doubt? Create a temporary vector:

%%time
tmp = gfu.vec.CreateVector()
t = 0; gfu.vec[:] = 0
while t < tend-0.5*dt:
    tmp.data = r + invm @ c * gfu.vec
    gfu.vec.data -= dt * tmp
    t += dt;
    Redraw(blocking=True)

CPU times: user 4.55 s, sys: 0 ns, total: 4.55 s
Wall time: 4.55 s

Speed up with geom_free integrals

With \(u(x) = \hat{u}(\hat{x})\), \(v(x) = \hat{v}(\hat{x})\) (directly mapped) and \(\mathbf{w}(x) = \frac{1}{|\operatorname{det}(F)|} F \cdot \hat{\mathbf{w}}(\hat{x})\) (Piola mapped) there holds

\[\int_T \mathbf{w} \cdot \nabla u \ v = \int_{\hat{T}} \hat{\mathbf{w}} \cdot \hat{\nabla} \hat{u} \ \hat{v}\]

and similarly for the facet integrals.

Integrals are independent of the “physical element” (up to ordering) similar to integrals in unit-2.11.

\(\leadsto\) allows to prepare intermediate computations for a large class of elements at once.

Piola wind:

gfwind = GridFunction(HDiv(mesh, order=order))
gfwind.Set(b)
b=gfwind; # <- overwrite old name "b"

Same as before, but with geom_free=True flag:

# part 1
uT, vT = VT.TnT()
c1 = BilinearForm(space=VT, nonassemble=True, geom_free=True)
c1 += -uT * b * grad(vT) * dx
# part 2
uT,uF = V.TrialFunction()
c2 = BilinearForm(trialspace=V, testspace=VT, nonassemble=True, geom_free=True)
# here uf-u = u_me + u_other - u_me is the neighbor value
c2 += b*n * IfPos(b*n, uT, uF-uT) * vT * dx(element_boundary=True)

c = c1.mat + c2.mat @ embV

%%time
t = 0; gfu.vec[:] = 0
while t < tend-0.5*dt:
    gfu.vec.data -= dt * (invm @ c * gfu.vec + r)
    t += dt
    Redraw(blocking=True)

CPU times: user 2.47 s, sys: 28.7 ms, total: 2.5 s
Wall time: 2.48 s