DG/HDG splitting¶

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

For a problem of the form $\partial_t u + A u + C u = 0$ We consider the operator splitting:

1.Step: Given $u^0$ , solve $t^n \to t^{n+1}$ : $\quad \partial_t u + C u = 0 \Rightarrow u^{\frac12}$
2.Step: Given $u^{\frac12}$ , solve $t^n \to t^{n+1}$ : $\quad \partial_t u + A u = 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. 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 $(u,p):[0,T] \to (H_{0,D}^1)^d \times L^2$ , s.t.

$\begin{align} \int_{\Omega} \partial_t u \cdot v + \int_{\Omega} \nu \nabla u \nabla v + u \cdot \nabla u v - \int_{\Omega} \operatorname{div}(v) p &= \int f v && \forall v \in (H_{0,D}^1)^d, \\ - \int_{\Omega} \operatorname{div}(u) q &= 0 && \forall q \in L^2, \\ \quad u(t=0) & = u_0 \end{align}$

Again, we consider the benchmark setup from http://www.featflow.de/en/benchmarks/cfdbenchmarking/flow/dfg_benchmark2_re100.html . The geometry is a 2D channel with a circular obstacle which is positioned (only slightly) off the center of the channel. The geometry: SchaeferTurekBenchmarkSetup

The left part (red) is the inflow boundary, The right part (green) is the outflow boundary. The viscosity is set to $\nu = 10^{-3}$ .

import netgen.gui
%gui tk
import tkinter
from math import pi
from ngsolve import *
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) )
order = 3
mesh.Curve(order)

bcs =  ('wall', 'outlet', 'wall', 'inlet')

For the HDG formulation we use the product space with

$BDM_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()

Stokes discretization / initial conditions¶

The bilinear form to the HDG discretized Stokes problem is:

$\begin{align*} K_h((u_T,u_F,p),(v_T,v_F,q) := & \displaystyle \sum_{T\in\mathcal{T}} \left\{ \int_{T} \nu {\nabla} {u_T} \! : \! {\nabla} {v_T} \ d {x} \right. \\ & \qquad \left. - \int_{\partial T} \nu \frac{\partial {u_T}}{\partial {n} } [ v ]_t \ d {s} - \int_{\partial T} \nu \frac{\partial {v_T}}{\partial {n} } [ u ]_t \ d {s} + \int_{\partial T} \nu \frac{\lambda k^2}{h} [u]_t [v]_t \ d {s} \right\}\\ & - \int_{\Omega} \operatorname{div}(u) q \ d {x} - \int_{\Omega} \operatorname{div}(v) p \ d {x} \end{align*}$

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

The mass matrix is simply

$\begin{align*} M_h((u_T,u_F,p),(v_T,v_F,q) := \int_{\Omega} u_T \cdot v_T d {x} \end{align*}$

nu = 0.001
alpha = 4

gradu = CoefficientFunction ( (grad(u),), dims=(2,2) )
gradv = CoefficientFunction ( (grad(v),), dims=(2,2) )

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

def tang(vec):
    return vec - (vec*n)*n

a = BilinearForm ( V, symmetric=True)
a += SymbolicBFI ( nu*InnerProduct ( gradu, gradv) )
a += SymbolicBFI ( nu*InnerProduct ( gradu.trans * n,  tang(vhat-v) ), element_boundary=True )
a += SymbolicBFI ( nu*InnerProduct ( gradv.trans * n,  tang(uhat-u) ), element_boundary=True )
a += SymbolicBFI ( nu*alpha*order*order/h * InnerProduct ( tang(vhat-v),  tang(uhat-u) ), element_boundary=True )
a += SymbolicBFI ( -div(u)*q -div(v)*p )
a.Assemble()

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

f = LinearForm ( V )

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")

res = f.vec.CreateVector()
res.data = f.vec - a.mat*gfu.vec
gfu.vec.data += invstokes * res

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

Application of convection¶

In the operator split 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 $(u,p)$ in HDG space: project into $\hat{u}$ in usual DG space
Solve $\partial_t \hat{u} = C \hat{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{u} 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)

mixed mass matrices¶

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}$ .

VL2 = L2(mesh, dim=mesh.dim, order=order, flags = { "dgjumps" : True } )
uL2 = VL2.TrialFunction()   # 1 x dim matrix
vL2 = VL2.TestFunction()    

gfuL2 = GridFunction(VL2)

bfmixed = BilinearForm ( V, VL2, flags = { "nonassemble" : True } )
bfmixed += SymbolicBFI ( vL2*u )

bfmixedT = BilinearForm ( VL2, V, flags = { "nonassemble" : True } )
bfmixedT += SymbolicBFI ( uL2*v )

convection operator¶

For the convection operation we use a standard Upwind scheme. Again, we do not set up the matrix as we are only interested in operator applications. For the advection velocity we use the $H(div)$ -conforming velocity (which is point-wise divergence free).

vel = gfu.components[0]
convL2 = BilinearForm(VL2, flags = { "nonassemble" : True } )
convL2 += SymbolicBFI( (-InnerProduct(grad(vL2).trans * vel, uL2.trans)) )
un = InnerProduct(vel,n)
upwindL2 = IfPos(un, un*uL2, un*uL2.Other(bnd=uin))
convL2 += SymbolicBFI( InnerProduct (upwindL2, vL2-vL2.Other()), VOL, skeleton=True )
convL2 += SymbolicBFI( InnerProduct (upwindL2, vL2), BND, skeleton=True )

solution of convection steps¶

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}$

def SolveConvectionSteps(gfuvec, res, tau, steps):
    bfmixed.Apply (gfuvec, gfuL2.vec) 
    VL2.SolveM(CoefficientFunction(1), gfuL2.vec)
    conv_applied = gfuL2.vec.CreateVector()
    for i in range(steps):
        convL2.Apply(gfuL2.vec,conv_applied)
        VL2.SolveM(CoefficientFunction(1), conv_applied)
        gfuL2.vec.data -= tau/steps * conv_applied
        #Redraw()    
    bfmixedT.Apply (gfuL2.vec, res)
    
#SolveConvectionSteps(gfu,res,0.01,1)    
#Draw(gfuL2, mesh, "velocity(L2)")
#Draw(Norm(gfuL2), mesh, "absvel(L2)")

Operator splitting method¶

# initial values again:
res.data = f.vec - a.mat*gfu.vec
gfu.vec.data += invstokes * res
t = 0

tend = 1
substeps = 10
tau = 0.01

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

while t < tend:
    SolveConvectionSteps(gfu.vec, res, tau, substeps)
    res.data -= mstar * gfu.vec
    gfu.vec.data += inv * res
    t += tau
    print ("t=", t)
    Redraw(blocking=True)

t= 0.01
t= 0.02
t= 0.03
t= 0.04
t= 0.05
t= 0.060000000000000005
t= 0.07
t= 0.08
t= 0.09
t= 0.09999999999999999
t= 0.10999999999999999
t= 0.11999999999999998
t= 0.12999999999999998
t= 0.13999999999999999
t= 0.15
t= 0.16
t= 0.17
t= 0.18000000000000002
t= 0.19000000000000003
t= 0.20000000000000004
t= 0.21000000000000005
t= 0.22000000000000006
t= 0.23000000000000007
t= 0.24000000000000007
t= 0.25000000000000006
t= 0.26000000000000006
t= 0.2700000000000001
t= 0.2800000000000001
t= 0.2900000000000001
t= 0.3000000000000001
t= 0.3100000000000001
t= 0.3200000000000001
t= 0.3300000000000001
t= 0.34000000000000014
t= 0.35000000000000014
t= 0.36000000000000015
t= 0.37000000000000016
t= 0.38000000000000017
t= 0.3900000000000002
t= 0.4000000000000002
t= 0.4100000000000002
t= 0.4200000000000002
t= 0.4300000000000002
t= 0.4400000000000002
t= 0.45000000000000023
t= 0.46000000000000024
t= 0.47000000000000025
t= 0.48000000000000026
t= 0.49000000000000027
t= 0.5000000000000002
t= 0.5100000000000002
t= 0.5200000000000002
t= 0.5300000000000002
t= 0.5400000000000003
t= 0.5500000000000003
t= 0.5600000000000003
t= 0.5700000000000003
t= 0.5800000000000003
t= 0.5900000000000003
t= 0.6000000000000003
t= 0.6100000000000003
t= 0.6200000000000003
t= 0.6300000000000003
t= 0.6400000000000003
t= 0.6500000000000004
t= 0.6600000000000004
t= 0.6700000000000004
t= 0.6800000000000004
t= 0.6900000000000004
t= 0.7000000000000004
t= 0.7100000000000004
t= 0.7200000000000004
t= 0.7300000000000004
t= 0.7400000000000004
t= 0.7500000000000004
t= 0.7600000000000005
t= 0.7700000000000005
t= 0.7800000000000005
t= 0.7900000000000005
t= 0.8000000000000005
t= 0.8100000000000005
t= 0.8200000000000005
t= 0.8300000000000005
t= 0.8400000000000005
t= 0.8500000000000005
t= 0.8600000000000005
t= 0.8700000000000006
t= 0.8800000000000006
t= 0.8900000000000006
t= 0.9000000000000006
t= 0.9100000000000006
t= 0.9200000000000006
t= 0.9300000000000006
t= 0.9400000000000006
t= 0.9500000000000006
t= 0.9600000000000006
t= 0.9700000000000006
t= 0.9800000000000006
t= 0.9900000000000007
t= 1.0000000000000007