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

Multigrid for hybrid methods

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

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, smoother="block", blocktype=["vertexpatch"], cycle=1)

with TaskManager(pajetrace=10**8):
    for l in range(6):
        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

f = LinearForm (x*v*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.09459423880925057
CG iteration 2, residual = 0.0006168717231344456
CG iteration 3, residual = 1.0961820396561478e-05
CG iteration 4, residual = 2.9780916114512884e-07
CG iteration 5, residual = 4.311403891010163e-09
CG iteration 6, residual = 6.483214687592656e-11

Draw (gfu.components[0]);

Nearly incompressible materials, H(div)-conforming HDG

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(pajetrace=10**8):
    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.8030966471886978 1.0365605651300314
l = 1 ndof = 8622 0.6859650712535184 1.0490739552090376
l = 2 ndof = 34158 0.6872476152185163 1.0778175315653233
l = 3 ndof = 135918 0.6908366644080608 1.0774921145559064

with TaskManager(pajetrace=10**8):
    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.01974650523880634
CG iteration 2, residual = 0.00020770356871837307
CG iteration 3, residual = 1.3541629720691253e-05
CG iteration 4, residual = 1.2095124817679166e-06
CG iteration 5, residual = 1.1795588285789404e-07
CG iteration 6, residual = 1.2588391516075687e-08
CG iteration 7, residual = 1.489958582040972e-09
CG iteration 8, residual = 1.6707190758288018e-10