This page was generated from unit-3.6-opsplit/opsplit.ipynb.

3.6 DG/HDG splitting

When solving unsteady problems with an operator-splitting method it might be benefitial to consider different space discretizations for different operators.

For a problem of the form

$${\partial }_{t}u+Au+Cu=0$$

We consider the operator splitting:

• 1.Step: Given $u^0$, solve $t^n\to t^{n+1}$: ${\partial }_{t}u+Cu=0\Rightarrow u^{\frac{1}{2}}$

• 2.Step: Given $u^{\frac{1}{2}}$, solve $t^n\to t^{n+1}$: ${\partial }_{t}u+Au=0\Rightarrow u^1$

This splitting is only first order accurate but allows to take different time discretizations for the substeps 1 and 2.

In this example we consider the Navier-Stokes problem again (cf. unit 3.2) and want to combine

• an $H(div)$ -conforming Hybrid DG method (which is a very good discretization for Stokes-type problems) and

• a standard upwind DG method for the discretization of the convection

The weak form is: Find $(\mathbf{u},p):[0,T]\to (H_{0,D}^1{)}^d\times L^2$, s.t. $$\begin{array}{rlrl}{\int }_{\mathrm{\Omega }}{\partial }_t\mathbf{u}\cdot v+{\int }_{\mathrm{\Omega }}\nu \mathrm{\nabla }\mathbf{u}\mathrm{\nabla }\mathbf{v}+\mathbf{u}\cdot \mathrm{\nabla }\mathbf{u}\mathbf{v}-{\int }_{\mathrm{\Omega }}\mathrm{div}(\mathbf{v})p& =\int f\mathbf{v}& & \mathrm{\forall }\mathbf{v}\in (H_{0,D}^1{)}^d,\\ -{\int }_{\mathrm{\Omega }}\mathrm{div}(\mathbf{u})q& =0& & \mathrm{\forall }q\in L^2,\\ \mathbf{u}(t=0)& ={\mathbf{u}}_0\end{array}$$

Again, we consider the benchmark setup from http://www.featflow.de/en/benchmarks/cfdbenchmarking/flow/dfg_benchmark2_re100.html . The geometry:

The viscosity is set to $\nu ={10}^{-3}$.

from ngsolve import *
import netgen.gui
from netgen.geom2d import SplineGeometry
geo = SplineGeometry()
geo.AddRectangle((0, 0), (2, 0.41), bcs = ("wall", "outlet", "wall", "inlet"))
geo.AddCircle((0.2, 0.2), r=0.05, leftdomain = 0, rightdomain = 1, bc = "cyl")
mesh = Mesh( geo.GenerateMesh(maxh = 0.08) ); Draw(mesh)
order = 3
mesh.Curve(order)

For the HDG formulation we use the product space with

• $BD{M}_k$: $H(div)$ conforming FE space (degree k)

• Vector facet space: facet functions of degree k (vector valued and only in tangential direction)

• piecewise polynomials up to degree $k-1$ for the pressure

HDG spaces

V1 = HDiv ( mesh, order = order, dirichlet = "wall|cyl|inlet" )
V2 = FESpace ( "vectorfacet", mesh, order = order,
              dirichlet = "wall|cyl|inlet" )
Q = L2( mesh, order = order-1)
V = FESpace ([V1,V2,Q])

u, uhat, p = V.TrialFunction()
v, vhat, q = V.TestFunction()

Parameter and directions (normal / tangential projection):

nu = 0.001 # viscosity
n = specialcf.normal(mesh.dim)
h = specialcf.mesh_size
def tang(vec):
    return vec - (vec*n)*n

Stokes discretization

The bilinear form to the HDG (see unit-2.8-DG for scalar HDG) discretized Stokes problem is: $$\begin{array}{rl}{K}_h(({\mathbf{u}}_T,{\mathbf{u}}_F,p),({\mathbf{v}}_T,{\mathbf{v}}_F,q):=& {\displaystyle \sum _{T\in \mathcal{T}}\{{\int }_T\nu \mathrm{\nabla }{\mathbf{u}}_T:\mathrm{\nabla }{\mathbf{v}}_{T}dx-{\int }_{\mathrm{\Omega }}\mathrm{div}(\mathbf{u})qdx-{\int }_{\mathrm{\Omega }}\mathrm{div}(\mathbf{v})pdx}\\ & -{\int }_{\partial T}\nu \frac{\partial {\mathbf{u}}_T}{\partial n}[\mathbf{v}{]}_{t}ds-{\int }_{\partial T}\nu \frac{\partial {\mathbf{v}}_T}{\partial n}[\mathbf{u}{]}_{t}ds+{\int }_{\partial T}\nu \frac{\alpha k^2}{h}[\mathbf{u}{]}_t[\mathbf{v}{]}_{t}ds\}\end{array}$$

where $[\cdot {]}_t$ is the tangential projection of the jump $(\cdot {)}_T-(\cdot {)}_F$.

alpha = 4  # SIP parameter
dS = dx(element_boundary=True)
a = BilinearForm ( V, symmetric=True)
a += nu * InnerProduct ( Grad(u), Grad(v) ) * dx
a += nu * InnerProduct ( Grad(u)*n, tang(vhat-v) ) * dS
a += nu * InnerProduct ( Grad(v)*n, tang(uhat-u) ) * dS
a += nu * alpha*order*order/h * InnerProduct(tang(vhat-v), tang(uhat-u)) * dS
a += (-div(u)*q - div(v)*p) * dx
a.Assemble()

<ngsolve.comp.BilinearForm at 0x7ff9680060b0>

The mass matrix is simply

$$\begin{array}{r}{M}_h(({\mathbf{u}}_T,{\mathbf{u}}_F,p),({\mathbf{v}}_T,{\mathbf{v}}_F,q):={\int }_{\mathrm{\Omega }}{\mathbf{u}}_T\cdot {\mathbf{v}}_{T}dx\end{array}$$

m = BilinearForm(V , symmetric=True)
m += u * v * dx
m.Assemble()

<ngsolve.comp.BilinearForm at 0x7ff970081670>

Initial conditions

As initial condition we solve the Stokes problem:

gfu = GridFunction(V)

U0 = 1.5
uin = CoefficientFunction( (U0*4*y*(0.41-y)/(0.41*0.41),0) )
gfu.components[0].Set(uin, definedon=mesh.Boundaries("inlet"))

invstokes = a.mat.Inverse(V.FreeDofs(), inverse="umfpack")
gfu.vec.data += invstokes @ -a.mat * gfu.vec

Draw( gfu.components[0], mesh, "velocity" )
Draw( gfu.components[2], mesh, "pressure" )
Draw( Norm(gfu.components[0]), mesh, "absvel(hdiv)")

Application of convection

In the operator splitted approach we want to apply only operator applications for the convection part. Further, we want to do this in a usual DG setting. As a model problem we use the following procedure:

• Given $(\mathbf{u},p)$ in HDG space: project into $\hat{\mathbf{u}}$ in usual DG space

• Solve ${\partial }_t\hat{\mathbf{u}}=C\hat{\mathbf{u}}$ by explicit scheme (involves convection evaluations and mass matrix operations only)

• Solve Unsteady Stokes step with r.h.s. from convection sub-problem. To this end evaluate mixed mass matrix $\int \hat{\mathbf{u}}\cdot \mathbf{v}$ to obtain a functional on the HDG space

For the projection steps we use mixed mass matrices:

• $M_m$ : $HDG\times DG\to \mathbb{R}$

• $M_m^T$ : $DG\times HDG\to \mathbb{R}$

• $M_{DG}$ : $DG\times DG\to \mathbb{R}$ (block diagonal)

To set up mixed mass matrices we use a bilinear form with two different FESpaces.

We do not assemble the matrices as we will only need the matrix-vector applications of $M_m$, $M_m^T$ and $M_{DG}^{-1}$.

There is a more elegant version (ConvertL2Operator, similar to TraceOperator) that we discuss below.

mixed mass matrices

VL2 = VectorL2(mesh, order=order, piola=True)
uL2, vL2 = VL2.TnT()
bfmixed = BilinearForm ( V, VL2, nonassemble=True )
bfmixed += vL2*u*dx
bfmixedT = BilinearForm ( VL2, V, nonassemble=True)
bfmixedT += uL2*v*dx

convection operator

• convection operation with standard Upwinding

• No set up of the matrix (only interested in operator applications)

• for the advection velocity we use the $H(div)$-conforming velocity (which is pointwise divergence free).

vel = gfu.components[0]
convL2 = BilinearForm(VL2, nonassemble=True )
convL2 += (-InnerProduct(Grad(vL2) * vel, uL2)) * dx
un = InnerProduct(vel,n)
upwindL2 = IfPos(un, un*uL2, un*uL2.Other(bnd=uin))

dskel_inner  = dx(skeleton=True)
dskel_bound  = ds(skeleton=True)

convL2 += InnerProduct (upwindL2, vL2-vL2.Other()) * dskel_inner
convL2 += InnerProduct (upwindL2, vL2) * dskel_bound

Solution of convection steps as an operator (BaseMatrix)

We now define the solution of the convection problem for an initial data $u$ in the HDG space. The return value (“res”) is $M_m\hat{u}$ where $\hat{u}$ is the solution of several explicit Euler steps of the convection problem

$${\partial }_t\hat{u}=-M_{DG}^{-1}C\hat{u}$$

class PropagateConvection(BaseMatrix):
    def __init__(self,tau,steps):
        super(PropagateConvection, self).__init__()
        self.tau = tau; self.steps = steps
        self.h = V.ndof; self.w = V.ndof # operator domain and range
        self.mL2 = VL2.Mass(Id(mesh.dim)); self.invmL2 = self.mL2.Inverse()
        self.vecL2 = bfmixed.mat.CreateColVector() # temp vector
    def Mult(self, x, y):
        self.vecL2.data = self.invmL2 @ bfmixed.mat * x # <- project from Hdiv to L2
        for i in range(self.steps):
            self.vecL2.data -= self.tau/self.steps * self.invmL2 @ convL2.mat * self.vecL2
        y.data = bfmixedT.mat * self.vecL2
    def CreateColVector(self):
        return CreateVVector(self.h)

The operator splitted method (Yanenko splitting)

Setup of $M^{\ast }$:

t = 0; tend = 0
tau = 0.01; substeps = 10

mstar = m.mat.CreateMatrix()
mstar.AsVector().data = m.mat.AsVector() + tau * a.mat.AsVector()
inv = mstar.Inverse(V.FreeDofs(), inverse="umfpack")

The time loop:

%%time
tend += 1
res = gfu.vec.CreateVector()
convpropagator = PropagateConvection(tau,substeps)
while t < tend:
    gfu.vec.data += inv @ (convpropagator - mstar) * gfu.vec
    t += tau
    print ("\r  t =", t, end="")
    Redraw(blocking=True)

  t = 1.0000000000000007CPU times: user 9.53 s, sys: 19.2 ms, total: 9.55 s
Wall time: 9.48 s

ConvertL2Operator

Instead of doing the conversion between the spaces “by hand”, we can use the convenience operator ConvertL2Operator that realizes $ M_{DG}^{-1} :nbsphinx-math:`cdot `M_m $:

HDG_to_L2 = V1.ConvertL2Operator(VL2) @ Embedding(V.ndof, V.Range(0)).T

The transpose of $ M_{DG}^{-1} \cdot `M_m $is$ M_m^T :nbsphinx-math:cdot M_{DG}^{-1}$ * ``V1.ConvertL2Operator(VL2)`: V1 $\to$ VL2 * V1.ConvertL2Operator(VL2).T: VL2${}^{\prime }$ $\to$ V1${}^{\prime }$

Overwrite the application and use the ConvertL2Operator instead ot the mixed mass matrices:

def NewMult(self, x, y):
    self.vecL2.data = HDG_to_L2 * x # <- project from Hdiv to L2
    for i in range(self.steps):
        self.vecL2.data -= tau/self.steps*self.invmL2@convL2.mat*self.vecL2
    y.data = HDG_to_L2.T @ self.mL2 * self.vecL2

PropagateConvection.Mult = NewMult

The time loop:

%%time
tend += 1
convpropagator = PropagateConvection(tau,substeps)
while t < tend:
    gfu.vec.data += inv @ (convpropagator - mstar) * gfu.vec
    t += tau
    print ("\r  t =", t, end="")
    Redraw(blocking=True)

  t = 2.0000000000000013CPU times: user 9.36 s, sys: 28.6 ms, total: 9.39 s
Wall time: 9.34 s