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
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. 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.
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}\).
[1]:
from netgen import gui
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) ); Draw(mesh)
order = 3
mesh.Curve(order)
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¶
[2]:
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):
[3]:
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:
where \([\cdot]_t\) is the tangential projection of the jump \((\cdot)_T - (\cdot)_F\).
[4]:
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()
The mass matrix is simply
[5]:
m = BilinearForm(V , symmetric=True)
m += SymbolicBFI( u * v )
m.Assemble()
Initial conditions¶
As initial condition we solve the Stokes problem:
[6]:
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¶
[7]:
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).
[8]:
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
[9]:
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^*\):
[10]:
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:
[11]:
%%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)
---------------------------------------------------------------------------
AttributeError Traceback (most recent call last)
<ipython-input-11-ac18c6825f32> in <module>
----> 1 get_ipython().run_cell_magic('time', '', 'tend += 1\nres = gfu.vec.CreateVector()\nconvpropagator = PropagateConvection(tau,substeps)\nwhile t < tend:\n gfu.vec.data += inv @ (convpropagator - mstar) * gfu.vec\n t += tau\n print ("\\r t =", t, end="")\n Redraw(blocking=True)\n')
/usr/local/lib/python3.6/dist-packages/IPython/core/interactiveshell.py in run_cell_magic(self, magic_name, line, cell)
2356 with self.builtin_trap:
2357 args = (magic_arg_s, cell)
-> 2358 result = fn(*args, **kwargs)
2359 return result
2360
</usr/local/lib/python3.6/dist-packages/decorator.py:decorator-gen-61> in time(self, line, cell, local_ns)
/usr/local/lib/python3.6/dist-packages/IPython/core/magic.py in <lambda>(f, *a, **k)
185 # but it's overkill for just that one bit of state.
186 def magic_deco(arg):
--> 187 call = lambda f, *a, **k: f(*a, **k)
188
189 if callable(arg):
/usr/local/lib/python3.6/dist-packages/IPython/core/magics/execution.py in time(self, line, cell, local_ns)
1284 if len(expr_ast.body) > 1 :
1285 expr_val=expr_ast.body[-1]
-> 1286 code_val = self.shell.compile(ast.Expression(expr_val.value)
1287 , '<timed eval>'
1288 , 'eval')
AttributeError: 'While' object has no attribute 'value'
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 $:
[12]:
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\('\) \(\to\) V1\('\)
Overwrite the application and use the ConvertL2Operator instead ot the mixed mass matrices:
[13]:
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:
[14]:
%%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)
---------------------------------------------------------------------------
AttributeError Traceback (most recent call last)
<ipython-input-14-f5c073bdb4fa> in <module>
----> 1 get_ipython().run_cell_magic('time', '', 'tend += 1\nconvpropagator = PropagateConvection(tau,substeps)\nwhile t < tend:\n gfu.vec.data += inv @ (convpropagator - mstar) * gfu.vec\n t += tau\n print ("\\r t =", t, end="")\n Redraw(blocking=True)\n')
/usr/local/lib/python3.6/dist-packages/IPython/core/interactiveshell.py in run_cell_magic(self, magic_name, line, cell)
2356 with self.builtin_trap:
2357 args = (magic_arg_s, cell)
-> 2358 result = fn(*args, **kwargs)
2359 return result
2360
</usr/local/lib/python3.6/dist-packages/decorator.py:decorator-gen-61> in time(self, line, cell, local_ns)
/usr/local/lib/python3.6/dist-packages/IPython/core/magic.py in <lambda>(f, *a, **k)
185 # but it's overkill for just that one bit of state.
186 def magic_deco(arg):
--> 187 call = lambda f, *a, **k: f(*a, **k)
188
189 if callable(arg):
/usr/local/lib/python3.6/dist-packages/IPython/core/magics/execution.py in time(self, line, cell, local_ns)
1284 if len(expr_ast.body) > 1 :
1285 expr_val=expr_ast.body[-1]
-> 1286 code_val = self.shell.compile(ast.Expression(expr_val.value)
1287 , '<timed eval>'
1288 , 'eval')
AttributeError: 'While' object has no attribute 'value'