This page was generated from unit-2.1.7-facetspace-mg/facetspace-mg.ipynb.

2.1.7 Multigrid for hybrid methods

Mixed methods for second order problems can often be reduced to the mesh facet, so called hybrid mixed methods. Simiar, hybrid DG methods introduce new variables on the facets, such that the bulk of element variables can be condensed out.

We show how to setup a multigrid preconditioner for hybrid methods. Interesting applications are nearly incompressible materials, or Stokes, discretized by \(H(\operatorname{div})\)-conforming HDG or hybrid mixed methods.

from ngsolve import *
from ngsolve.webgui import Draw
from ngsolve.la import EigenValues_Preconditioner

The hybrid DG method:

On \(V_h = V_T \times V_F = P^k({\mathcal T}) \times P^k ({\mathcal F})\) we define the bilinear form

\[\DeclareMathOperator{\Div}{div} \sum_T \int_T \nabla u \nabla v - \sum_T \int_{\partial T} \tfrac{\partial u}{\partial n} (v-\widehat v) - \sum_T \int_{\partial T} \tfrac{\partial v}{\partial n} (u-\widehat u) + \frac{\alpha p^2}{h} \sum_F \int_F (u-\widehat u)(v-\widehat v)\]

Element variables can be condensed out, which leads to a system reduced to the Skeleton.

When splitting a large triangle \(T_H\) into small trianles, the functions on \(\partial T_H\) have a canonical representation on the facets of the fine triangles. However, facet variables on internal edges of \(T_H\) are not defined by embedding. The HarmonicProlongation provides the energy optimal extension to the internal edges. To define energy optimal we need the energy defined by a bilinear form.

ngmesh = unit_square.GenerateMesh(maxh=2)
mesh = Mesh(ngmesh)

order = 3
fes = L2(mesh, order=order) * FacetFESpace(mesh, order=order, hoprolongation=True, dirichlet=".*")
(u,uhat), (v,vhat) = fes.TnT()

n = specialcf.normal(2)
h = specialcf.mesh_size
dS = dx(element_vb=BND)

HDGform = u*v*dx+ grad(u)*grad(v)*dx - n*grad(u)*(v-vhat)*dS - n*grad(v)*(u-uhat)*dS + 5*(order+1)**2/h*(u-uhat)*(v-vhat)*dS
bfa = BilinearForm(HDGform, condense=True).Assemble()
fes.SetHarmonicProlongation(bfa, inverse="sparsecholesky")
pre = preconditioners.MultiGrid(bfa, blocktype=["vertexpatch"], cycle=1)

with TaskManager():
    for l in range(7):
        mesh.Refine()
        bfa.Assemble()
        # pre.Update()
        lam = EigenValues_Preconditioner(bfa.mat, pre)
        print ("l =", l, "ndof =", fes.ndof, "lam_min/lam_max = ", lam[0], lam[-1])

l = 0 ndof = 144 lam_min/lam_max =  0.9999999999999993 0.9999999999999993
l = 1 ndof = 544 lam_min/lam_max =  0.9573848616699543 0.9999918601826863
l = 2 ndof = 2112 lam_min/lam_max =  0.9261189386281288 0.9999692236038735
l = 3 ndof = 8320 lam_min/lam_max =  0.9325349652513784 0.9998939800629592
l = 4 ndof = 33024 lam_min/lam_max =  0.92945510895588 0.9997230751035944
l = 5 ndof = 131584 lam_min/lam_max =  0.9282844459834527 0.9996349409051806
l = 6 ndof = 525312 lam_min/lam_max =  0.9280114490358107 0.9995334765161776

f = LinearForm (x*v*dx).Assemble()
gfu = GridFunction(fes)
gfu.vec[:]=0
with TaskManager():
    Solve (bfa*gfu==f, pre, lin_solver=solvers.CGSolver, printrates=True)

CG iteration 1, residual = 0.09459461645928992
CG iteration 2, residual = 0.0005944516866628238
CG iteration 3, residual = 7.695100604718358e-06
CG iteration 4, residual = 2.0699972197723346e-07
CG iteration 5, residual = 3.8945463326402356e-09
CG iteration 6, residual = 4.9504671345225695e-11
CG iteration 7, residual = 1.0763303055548899e-12
CG iteration 8, residual = 2.096450968451275e-14

Draw (gfu.components[0]);

Hybrid-mixed methods:

Find \(\sigma, u, \widehat u \in \Sigma_h \times V_h \times F_h\):

\[\begin{split}\DeclareMathOperator{\Div}{div} \begin{array}{ccccccll} \int a \sigma \tau & + & \sum_T \int_T \Div \tau \, u & + & \sum_F \int_F [\tau_n] \widehat u & = & 0 & \forall \, \tau \in \Sigma \\ \int \Div \sigma \, v &&&&& = & \int f v & \forall \, v \in V_h \\ \int [ \sigma_n ] \, \widehat v &&&&& = & \int_{\Gamma_n} g \widehat v & \forall \, \widehat v \in F_h \end{array}\end{split}\]

where \(\Sigma_h\) is an discontinuous \(H(div)\) finite element space, \(V_h\) a sub-space of \(L_2\), and \(F_h\) consists of polynomials on every facet.

ngmesh = unit_square.GenerateMesh(maxh=0.2)
mesh = Mesh(ngmesh)

order = 2

fesSigma = PrivateSpace(Discontinuous(HDiv(mesh, order=order, RT=True)))
fesL2 = L2(mesh, order=order)
fesFacet = FacetFESpace(mesh, order=order, hoprolongation=True, dirichlet=".*")
fes = fesSigma*fesL2*fesFacet

(sigma, u,uhat), (tau, v,vhat) = fes.TnT()
n = specialcf.normal(2)
dS = dx(element_vb=BND)
mixedform = -sigma*tau*dx
mixedform += div(sigma)*v*dx - sigma*n*vhat*dS
mixedform += div(tau)*u*dx - tau*n*uhat*dS

bfa = BilinearForm(mixedform, condense=True).Assemble()
fes.SetHarmonicProlongation(bfa, inverse="sparsecholesky")
pre = preconditioners.MultiGrid(bfa, blocktype=["vertexpatch"], cycle=1)

for l in range(5):
    mesh.Refine()
    bfa.Assemble()
    lam = EigenValues_Preconditioner(bfa.mat, pre)
    print ("l =", l, "ndof =", fes.ndof, "lam_min/lam_max = ", lam[0], lam[-1])

l = 0 ndof = 2328 lam_min/lam_max =  0.8830690135570511 1.000082660916307
l = 1 ndof = 9192 lam_min/lam_max =  0.8695401255366382 0.9995622254215344
l = 2 ndof = 36528 lam_min/lam_max =  0.8480639402311638 0.9985740487773203
l = 3 ndof = 145632 lam_min/lam_max =  0.830706421034568 0.9979298221317969
l = 4 ndof = 581568 lam_min/lam_max =  0.8389819990459044 0.9976337192728126

f = LinearForm(x*v*dx).Assemble()
gfu = GridFunction(fes)
Solve(bfa*gfu==f, pre, solvers.CGSolver, printrates=True)
Draw (gfu.components[1]);

CG iteration 1, residual = 0.09771383816444208
CG iteration 2, residual = 0.0011964806411029606
CG iteration 3, residual = 3.13755506370263e-05
CG iteration 4, residual = 9.5216727952247e-07
CG iteration 5, residual = 3.248449810413196e-08
CG iteration 6, residual = 1.649970071284961e-09
CG iteration 7, residual = 8.248735706462719e-11
CG iteration 8, residual = 3.983130055402826e-12
CG iteration 9, residual = 2.0772868157657282e-13
CG iteration 10, residual = 1.2019646043967861e-14

Nearly incompressible materials, H(div)-conforming HDG

[Lehrenfeld+Schöberl, 2016]

ngmesh = unit_square.GenerateMesh(maxh=0.3)
mesh = Mesh(ngmesh)

order = 3

fesT = HDiv(mesh, order=order, hoprolongation=True, dirichlet=".*")
fesF = TangentialFacetFESpace(mesh, order=order, hoprolongation=True, highest_order_dc=True, dirichlet=".*")
fes = fesT*fesF

(u,uhat), (v,vhat) = fes.TnT()
n = specialcf.normal(2)
def tang(v): return v-(v*n)*n
h = specialcf.mesh_size
dS = dx(element_vb=BND)

HDGform = InnerProduct(Grad(u),Grad(v))*dx - (Grad(u)*n)*tang(v-vhat)*dS - (Grad(v)*n)*tang(u-uhat)*dS \
    + 1*(order+1)**2/h*tang(u-uhat)*tang(v-vhat)*dS

bfa = BilinearForm(HDGform + 1e3*div(u)*div(v)*dx, condense=True).Assemble()
fes.SetHarmonicProlongation(bfa)
pre = preconditioners.MultiGrid(bfa, smoother="block", smoothingsteps=1, blocktype=["vertexpatch"], cycle=1)

with TaskManager():
    for l in range(4):
        mesh.Refine()
        bfa.Assemble()
        lam = EigenValues_Preconditioner(bfa.mat, pre)
        print ("l =", l, "ndof =", fes.ndof, lam[0], lam[-1])

l = 0 ndof = 2190 0.8031555871568468 1.0389867675724136
l = 1 ndof = 8622 0.6881351203081549 1.0480292176047623
l = 2 ndof = 34158 0.6819544900463768 1.0702830065295164
l = 3 ndof = 135918 0.6926495255892722 1.0752680621326927

with TaskManager():
    f = LinearForm ((0.5-y)*v[0]*dx).Assemble()
gfu = GridFunction(fes)
gfu.vec[:]=0
with TaskManager(pajetrace=10**8):
    solvers.BVP(bfa, f, gfu, pre, print=True)

CG iteration 1, residual = 0.01974650523880605
CG iteration 2, residual = 0.00020770356871082298
CG iteration 3, residual = 1.3541629719963928e-05
CG iteration 4, residual = 1.2095124816793476e-06
CG iteration 5, residual = 1.179558828526886e-07
CG iteration 6, residual = 1.258839151557343e-08
CG iteration 7, residual = 1.4899585819634964e-09
CG iteration 8, residual = 1.6707190757074536e-10

Draw (gfu.components[0])

BaseWebGuiScene

Nearly incompressible materials / Stokes in 3D

ngmesh = unit_cube.GenerateMesh(maxh=2)
mesh = Mesh(ngmesh)

order = 2

fesT = HDiv(mesh, order=order, hoprolongation=True, dirichlet=".*")
fesF = TangentialFacetFESpace(mesh, order=order, hoprolongation=True, highest_order_dc=True, dirichlet=".*")
fes = fesT*fesF

(u,uhat), (v,vhat) = fes.TnT()
n = specialcf.normal(3)
def tang(v): return v-(v*n)*n
h = specialcf.mesh_size
dS = dx(element_vb=BND)

HDGform = 0.001*u*v*dx+InnerProduct(Grad(u),Grad(v))*dx - (Grad(u)*n)*tang(v-vhat)*dS - (Grad(v)*n)*tang(u-uhat)*dS \
    + 5*(order+1)**2/h*tang(u-uhat)*tang(v-vhat)*dS

bfa = BilinearForm(HDGform + 1e3*div(u)*div(v)*dx, condense=True).Assemble()
fes.SetHarmonicProlongation(bfa)
pre = preconditioners.MultiGrid(bfa, smoother="block", smoothingsteps=3, blocktype=["edgepatch"], cycle=1)

with TaskManager():
    for l in range(3):
        mesh.Refine()
        bfa.Assemble()
        # pre.Update()
        lam = EigenValues_Preconditioner(bfa.mat, pre)
        print ("l =", l, "ndof =", fes.ndof, "lam min/max = ", lam[0], lam[-1])

l = 0 ndof = 5574 lam min/max =  0.5060508739232161 1.2112630153266313
l = 1 ndof = 43470 lam min/max =  0.2580696754603477 1.2905858121466474
l = 2 ndof = 342894 lam min/max =  0.2376744824447379 1.4487474122388502

with TaskManager():
    f = LinearForm ((0.5-y)*v[0]*dx).Assemble()
gfu = GridFunction(fes)
gfu.vec[:]=0

with TaskManager():
    Solve(bfa*gfu==f, pre, solvers.CGSolver, printrates=True)

CG iteration 1, residual = 0.016109906113163878
CG iteration 2, residual = 0.0036631197303151736
CG iteration 3, residual = 0.0013763841180133103
CG iteration 4, residual = 0.0005505700530464946
CG iteration 5, residual = 0.00023667946218137526
CG iteration 6, residual = 0.0001001439542047559
CG iteration 7, residual = 4.349404812113597e-05
CG iteration 8, residual = 1.6590186209442128e-05
CG iteration 9, residual = 6.187485725236571e-06
CG iteration 10, residual = 2.3330442207180076e-06
CG iteration 11, residual = 9.509298878743054e-07
CG iteration 12, residual = 3.7770775958317634e-07
CG iteration 13, residual = 1.6068815303140225e-07
CG iteration 14, residual = 6.952194561493456e-08
CG iteration 15, residual = 2.782312459468385e-08
CG iteration 16, residual = 1.0924337753518158e-08
CG iteration 17, residual = 4.487871076211575e-09
CG iteration 18, residual = 1.8627230253709616e-09
CG iteration 19, residual = 7.402872804341187e-10
CG iteration 20, residual = 2.9942322000888977e-10
CG iteration 21, residual = 1.2726933117910823e-10
CG iteration 22, residual = 5.4302541953643354e-11
CG iteration 23, residual = 2.2027345474574765e-11
CG iteration 24, residual = 8.934086296736379e-12
CG iteration 25, residual = 3.678442279210992e-12
CG iteration 26, residual = 1.457910563469523e-12
CG iteration 27, residual = 5.88632994935434e-13
CG iteration 28, residual = 2.568065072287564e-13
CG iteration 29, residual = 1.0908752441007383e-13
CG iteration 30, residual = 4.422174497439918e-14
CG iteration 31, residual = 1.80457472127406e-14
CG iteration 32, residual = 7.616932077299142e-15

clipping = { "function" : True,  "pnt" : (0.5,0.5,0.51), "vec" : (0,0,-1) }
Draw (gfu.components[0], order=2, clipping=clipping);

Flow channel in 3D

from ngsolve import *
from netgen.occ import *
from ngsolve.webgui import Draw
from ngsolve.krylovspace import CGSolver

box = Box((0,0,0), (2.5,0.41,0.41))
box.faces.name="wall"
box.faces.Min(X).name="inlet"
box.faces.Max(X).name="outlet"
cyl = Cylinder((0.5,0.2,0), Z, h=0.41,r=0.05)
cyl.faces.name="cyl"
shape = box-cyl

mesh = shape.GenerateMesh(maxh=0.2).Curve(3)
Draw (mesh);

order = 2

fesT = HDiv(mesh, order=order, hoprolongation=True, dirichlet="wall|inlet|cyl")
fesF = TangentialFacetFESpace(mesh, order=order, hoprolongation=True, highest_order_dc=True, dirichlet=".*")
fes = fesT*fesF

(u,uhat), (v,vhat) = fes.TnT()
n = specialcf.normal(3)
def tang(v): return v-(v*n)*n
h = specialcf.mesh_size
dS = dx(element_vb=BND)

HDGform = 0.001*u*v*dx+InnerProduct(Grad(u),Grad(v))*dx - (Grad(u)*n)*tang(v-vhat)*dS - (Grad(v)*n)*tang(u-uhat)*dS \
    + 5*(order+1)**2/h*tang(u-uhat)*tang(v-vhat)*dS

with TaskManager():
    bfa = BilinearForm(HDGform + 1e3*div(u)*div(v)*dx, condense=True).Assemble()
    fes.SetHarmonicProlongation(bfa)
    pre = preconditioners.MultiGrid(bfa, inverse="sparsecholesky", smoother="block", smoothingsteps=1, blocktype=["edgepatch"], cycle=1)

with TaskManager():
    for l in range(1):
        mesh.Refine()
        bfa.Assemble()
        lam = EigenValues_Preconditioner(bfa.mat, pre)
        print ("l =", l, "ndof =", fes.ndof, "lam min/max = ", lam[0], lam[-1])

l = 0 ndof = 1779310 lam min/max =  0.07395292018165522 3.336861388881027

gfu = GridFunction(fes)

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

inv = CGSolver(bfa.mat, pre.mat, printrates=True, tol=1e-5)

with TaskManager():
    gfu.vec.data -= inv@bfa.mat * gfu.vec
gfu.vec.data += bfa.harmonic_extension * gfu.vec

CG iteration 1, residual = 14.705439018174086
CG iteration 2, residual = 0.6535781335123068
CG iteration 3, residual = 0.2949122809986667
CG iteration 4, residual = 0.1501249968250976
CG iteration 5, residual = 0.11534735471253654
CG iteration 6, residual = 0.0720507592698503
CG iteration 7, residual = 0.051982089875520214
CG iteration 8, residual = 0.038709221148398185
CG iteration 9, residual = 0.022859269501356475
CG iteration 10, residual = 0.016697123231952635
CG iteration 11, residual = 0.01089695408606514
CG iteration 12, residual = 0.007552851135448333
CG iteration 13, residual = 0.005539446505053609
CG iteration 14, residual = 0.0041399961601813725
CG iteration 15, residual = 0.0030604930093198024
CG iteration 16, residual = 0.0021714417553159214
CG iteration 17, residual = 0.001557073093463609
CG iteration 18, residual = 0.0011741926334623994
CG iteration 19, residual = 0.0008572928809040887
CG iteration 20, residual = 0.0006182722955589243
CG iteration 21, residual = 0.0005100286526489225
CG iteration 22, residual = 0.00036400544225980043
CG iteration 23, residual = 0.00027695011270463613
CG iteration 24, residual = 0.00020615710374360176
CG iteration 25, residual = 0.0001450445075884344

clipping = { "function" : True,  "pnt" : (1,0.2,0.2 ), "vec" : (0,0,-1.0) }
Draw (gfu.components[0], mesh, order=2, clipping=clipping);

Pressure:

Draw (div(gfu.components[0]), mesh, order=2, clipping=clipping, draw_surf=False);