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.

[1]:

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.

[2]:

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)

[3]:

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

[4]:

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

[5]:

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.

[6]:

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.0000826609163065
l = 1 ndof = 9192 lam_min/lam_max =  0.8695401255366377 0.9995622254215335
l = 2 ndof = 36528 lam_min/lam_max =  0.848063940231164 0.9985740487773187
l = 3 ndof = 145632 lam_min/lam_max =  0.830706421034594 0.9979298221317767
l = 4 ndof = 581568 lam_min/lam_max =  0.8389819990458989 0.9976337192727811

[7]:

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

CG iteration 1, residual = 0.1856429027154921
CG iteration 2, residual = 0.0023298507638458642
CG iteration 3, residual = 5.4897196584525635e-05
CG iteration 4, residual = 1.8801116118647993e-06
CG iteration 5, residual = 6.651127368304594e-08
CG iteration 6, residual = 3.242660560201135e-09
CG iteration 7, residual = 1.7540201305378962e-10
CG iteration 8, residual = 8.506846723262754e-12
CG iteration 9, residual = 4.50287845997025e-13
CG iteration 10, residual = 2.581624946882175e-14

Nearly incompressible materials, H(div)-conforming HDG¶

[8]:

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)

[9]:

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.8031555871568414 1.0389867675724214
l = 1 ndof = 8622 0.6881351203081374 1.0480292176047754
l = 2 ndof = 34158 0.6819544900463874 1.0702830065295148
l = 3 ndof = 135918 0.6926495255892826 1.0752680621327033

[10]:

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.01974650523880636
CG iteration 2, residual = 0.00020770356871556786
CG iteration 3, residual = 1.3541629720175803e-05
CG iteration 4, residual = 1.2095124816987853e-06
CG iteration 5, residual = 1.1795588285255381e-07
CG iteration 6, residual = 1.258839151560123e-08
CG iteration 7, residual = 1.4899585819894852e-09
CG iteration 8, residual = 1.670719075746596e-10

Nearly incompressible materials / Stokes in 3D¶

[11]:

ngmesh = unit_cube.GenerateMesh(maxh=0.3)
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)

[12]:

with TaskManager():
    for l in range(2):
        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 = 93831 lam min/max =  0.09744034470926972 1.3584288477036652
l = 1 ndof = 738607 lam min/max =  0.11159161337425566 1.5639580133845095

[13]:

with TaskManager(pajetrace=10**8):
    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.01649864156546645
CG iteration 2, residual = 0.002880221500837048
CG iteration 3, residual = 0.0012973281127649434
CG iteration 4, residual = 0.0006291767071514194
CG iteration 5, residual = 0.00034212404496442473
CG iteration 6, residual = 0.00018986161398905436
CG iteration 7, residual = 0.00010495480172690592
CG iteration 8, residual = 5.6792676001633144e-05
CG iteration 9, residual = 3.223619624399929e-05
CG iteration 10, residual = 1.7889721196786516e-05
CG iteration 11, residual = 9.932461473408207e-06
CG iteration 12, residual = 5.702882845951065e-06
CG iteration 13, residual = 3.4772537242985423e-06
CG iteration 14, residual = 2.2030120819785107e-06
CG iteration 15, residual = 1.3845181990633264e-06
CG iteration 16, residual = 8.452023058454922e-07
CG iteration 17, residual = 4.865032668541236e-07
CG iteration 18, residual = 2.7214663873327504e-07
CG iteration 19, residual = 1.4658773546049505e-07
CG iteration 20, residual = 7.661756638640799e-08
CG iteration 21, residual = 4.1094123304675296e-08
CG iteration 22, residual = 2.2676336133507237e-08
CG iteration 23, residual = 1.2750052703273745e-08
CG iteration 24, residual = 7.156068104141158e-09
CG iteration 25, residual = 4.064413378474485e-09
CG iteration 26, residual = 2.320562805494629e-09
CG iteration 27, residual = 1.3027530001561051e-09
CG iteration 28, residual = 7.299511520183448e-10
CG iteration 29, residual = 4.12605356640944e-10
CG iteration 30, residual = 2.3526123467024506e-10
CG iteration 31, residual = 1.3444504125613385e-10
CG iteration 32, residual = 7.610937986484883e-11
CG iteration 33, residual = 4.3867287103046e-11
CG iteration 34, residual = 2.5753023412050017e-11
CG iteration 35, residual = 1.4921637032046386e-11
CG iteration 36, residual = 8.346406504216869e-12
CG iteration 37, residual = 4.499350043747646e-12
CG iteration 38, residual = 2.3647568485386997e-12
CG iteration 39, residual = 1.2285664041302538e-12
CG iteration 40, residual = 6.437836701077966e-13
CG iteration 41, residual = 3.352332863064103e-13
CG iteration 42, residual = 1.7399475039986366e-13
CG iteration 43, residual = 9.162203999145559e-14
CG iteration 44, residual = 4.734363660329852e-14
CG iteration 45, residual = 2.4659671922482638e-14
CG iteration 46, residual = 1.2791468177817011e-14

[14]:

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

[ ]: