This page was generated from unit-3.5.1-dgapply/dgapply-scalar.ipynb.

3.5.1 Operator applications for Discontinuous Galerkin discretizations - linear transport

In the previous units we already saw DG discretizations. However, we only solved stationary problems in these cases so that DG discretizations have only been treated in terms of solving linear systems (where HDG has some advantages over classical DG discretizations).

Here, however, we treat DG discretizations in the case where operator applications are more relevant than the solution of linear systems, as e.g. in explicit time integration schemes. Note that for continuous Galerkin methods, in unit 3.2 we already saw simple ways on how to use operator applications within an IMEX scheme. We will develop this further and take a look at efficient ways to do it for DG methods.

In this unit we start with a model PDE problem: The unsteady 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}

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

# import libraries, define geometry and generate mesh
from ngsolve import *
from ngsolve.webgui import *
from netgen.occ import *
shape = Rectangle(1,1).Face()
shape.edges.Min(X).name="left"
shape.edges.Max(X).name="right"
shape.edges.Min(Y).name="bottom"
shape.edges.Max(Y).name="top"
geom = OCCGeometry(shape, dim=2)
mesh = Mesh(geom.GenerateMesh(maxh=0.1))

We consider higher order discretizations and set \(k=5\) here:

order = 4

Warning

• We will make some timing measurements.

• To avoid the impact of initialization overheads you may want to increase \(k\) and decrease the mesh size.

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.

b = CF((1+sin(4*pi*y),2))
ubnd = 0.1 * IfPos(x-0.125,IfPos(0.625-x,1+cos(8*pi*x),0),0)

Draw(b,mesh,vectors = {"grid_size" : 32}, min = 2, max = 3, autoscale=False,
     height = "30vh")

BaseWebGuiScene

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.

In the following we only treat the explicit time integration case here (take a look at unit 2.8 for solving linear systems with DG formulations).

Explicit time stepping with a DG formulation

Upwind-in-space combined with an explicit Euler scheme in time yields:

\[\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\).

# roughly estimated time step:
dt = 1e-3 / (order+1)

Computing convection applications \(C u^n\)

We recall:

• 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, settings = { "subdivision" : 10 },
#     min=-0.1,max=2.1,autoscale=False)

#Operator application (equivalent to assemble and mult but faster)
res = gfu.vec.CreateVector()
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)
Cell In[10], line 3
      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)

We also have to deal with the mass matrix:

• 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() # not a sparse matrix
print(invm)

Print base-matrix

The (simple) time loop

def SolveAndVisualize(conv_op, t=0, tend=0.6, dt=dt, dt_sample=0.06, rhs_vec = None):
    gfu.vec[:] = 0; cnt = 0;
    sample_rate = int(floor(dt_sample*tend/dt))
    gfu_t = GridFunction(gfu.space,multidim=0)
    gfu_t.AddMultiDimComponent(gfu.vec)
    if rhs_vec:
        inv_rhs = invm * rhs_vec
    while t < tend-0.5*dt:
        if rhs_vec:
            res = invm @ conv_op * gfu.vec + inv_rhs
        else:
            res = invm @ conv_op * gfu.vec
        gfu.vec.data -= dt * res
        t += dt; cnt += 1; print("\r t=",t,end="")
        if (cnt+1) % sample_rate == 0:
            gfu_t.AddMultiDimComponent(gfu.vec)
    return gfu_t

%%time
gfu_t = SolveAndVisualize(c.mat)

 t= 0.5999999999999666CPU times: user 9.38 s, sys: 33.7 ms, total: 9.41 s
Wall time: 9.38 s

Draw(gfu_t,mesh,interpolate_multidim=True,animate=True,
     min=0, max=0.2, autoscale=False, deformation=True, height = "60vh")

BaseWebGuiScene

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

Next run

%%time
gfu_t = SolveAndVisualize(c.mat)
#Draw(gfu_t,mesh,interpolate_multidim=True,animate=True,
#     min=0, max=0.2, autoscale=False, deformation=True)

 t= 0.5999999999999666CPU times: user 7.52 s, sys: 49.4 ms, total: 7.57 s
Wall time: 7.53 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 = VT*VF

trace = VT.TraceOperator(VF, False)
#info output about the trace operator
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 3480 dofs.
The trace space VF has 1840 dofs.
trace represents an 1840 x 3480 operator

example application of the trace 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)

Part 1:

uT, vT = VT.TnT()
c1 = BilinearForm(space=VT, nonassemble=True)                     # part 1
c1 += -uT * b * grad(vT) * dx

Part 2:

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
# on bnd: uf = u_me and hence uf-u = 0
c2 += b*n * IfPos(b*n, uT, uF-uT) * vT * dx(element_boundary=True)

Right hand side (Part 3):

rhs = LinearForm(VT)                                              # part 3
rhs += b*n * IfPos(b*n, 0, ubnd) * vT * dskel_bound # dskel_inflow
rhs.Assemble()

<ngsolve.comp.LinearForm at 0x7fe8206544f0>

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
gfu_t = SolveAndVisualize(c, rhs_vec=rhs.vec)

 t= 0.5999999999999666CPU times: user 8.13 s, sys: 52.6 ms, total: 8.19 s
Wall time: 8.15 s

#Draw(gfu_t,mesh,interpolate_multidim=True,animate=True,
#     min=0, max=0.2, autoscale=False, deformation=True)

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
gfu_t = SolveAndVisualize(c, rhs_vec=rhs.vec)

 t= 0.5999999999999666CPU times: user 4.54 s, sys: 35.4 ms, total: 4.58 s
Wall time: 4.54 s

#Draw(gfu_t,mesh,interpolate_multidim=True,animate=True,
#     min=0, max=0.2, autoscale=False, deformation=True)

And now with the TaskManager:

%%time
with TaskManager():
    gfu_t = SolveAndVisualize(c, rhs_vec=rhs.vec)

 t= 0.5999999999999666CPU times: user 7.38 s, sys: 956 ms, total: 8.34 s
Wall time: 2.08 s