This page was generated from unit-2.1.6-highorder-mg/highorder-multigrid.ipynb.

2.1.6. Multigrid for high order finite element spaces

Officially released with NGSolve-2504, most high order spaces provide high order polynomial preserving prolongation operators. Current restrictions:

• require simplicial meshes

• refinement by bisection

• require uniform polynomial order

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

Many FESpaces provide now a high-order accurate prolongation operator. It has to be enabled by the flag hoprolongation=True, maybe becoming default in future.

mesh = Mesh(unit_square.GenerateMesh(maxh=0.2))
gfu = GridFunction(H1(mesh, order=10, hoprolongation=True))
gfu.Set (sin(50*x*y))
Draw (gfu, order=10);

for l in range(2):
    mesh.Refine()
Draw (gfu, order=5);

Multigrid preconditioners for high order spaces

If the high order prolongation is enabled, the multigrid preconditioner uses the high order discretization on the mesh hierarchy. If not, the coarse grid spaces use the lowest order spaces in the background.

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

fes = H1(mesh, order=5, hoprolongation=True, dirichlet=".*")
u,v = fes.TnT()

a = BilinearForm(grad(u)*grad(v)*dx).Assemble()
pre = preconditioners.MultiGrid(a, blocktype=["vertexpatch"])
# pre = preconditioners.MultiGrid(a, smoother="block", blocktype=["vertexedge","face"])
# pre = preconditioners.MultiGrid(a, smoother="block", blocktype=["vertpatch"])

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

2 2801 0.9737767048916046 0.9999650900459265
3 11001 0.9731893341714941 0.9999569893694966
4 43601 0.9738508899616769 0.9999593243766862
5 173601 0.9740047337504241 0.9999581308391128

f = LinearForm(x*v*dx).Assemble()
gfu = GridFunction(fes)

Solve(a * gfu == f, dirichlet=x*(1-x), lin_solver=solvers.CGSolver, pre=pre, printrates=True)

CG iteration 1, residual = 2.500665990770517
CG iteration 2, residual = 0.024383987764233152
CG iteration 3, residual = 0.00018799791819350655
CG iteration 4, residual = 1.1484794277901268e-06
CG iteration 5, residual = 9.252209912710351e-09
CG iteration 6, residual = 6.018344974750805e-11
CG iteration 7, residual = 4.580233475866812e-13

ea = { "euler_angles" : (-70, 0,-55) }
Draw (gfu, deformation=True, **ea);

High order with static condensation

For high order methods, static condensation may save a lot of computation. However, canonical prolongation for the skeleton variables only does not preserve high order polynomials. Here the HarmonicProlongation comes into play: It prolongates functions on the boundaries of the coarse grid elements, and then solves local Dirichlet problems for the dofs on the fine-grid skeleton inside the coarse grid elements.

The Dirichlet problem is solved for the problem-specific bilinear form, which has to be provided when enabling the HarmonicProlongation:

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

fes = H1(mesh, order=5, hoprolongation=True, dirichlet=".*")
u,v = fes.TnT()

a = BilinearForm(grad(u)*grad(v)*dx, condense=True).Assemble()
fes.SetHarmonicProlongation(a)
pre = preconditioners.MultiGrid(a, blocktype="vertexpatch")

for l in range(5):
    mesh.Refine()
    a.Assemble()
    lam = EigenValues_Preconditioner(a.mat, pre)
    print (mesh.levels, fes.ndof, lam[0], lam[-1])

2 2801 0.8873455539760278 0.9998975910633985
3 11001 0.8957401274225056 0.9998894014559885
4 43601 0.8945786078982365 0.9998863285007136
5 173601 0.894386169292487 0.9998960435342532
6 692801 0.8953852278217878 0.9998858093863445

f = LinearForm(x*v*dx).Assemble()
gfu = GridFunction(fes)

Solve(a * gfu == f, dirichlet=x*(1-x), lin_solver=solvers.CGSolver, pre=pre, printrates=True)

CG iteration 1, residual = 3.4886087977865676
CG iteration 2, residual = 0.12345444219082062
CG iteration 3, residual = 0.0036390880918837146
CG iteration 4, residual = 9.08036702501544e-05
CG iteration 5, residual = 2.7831282126523136e-06
CG iteration 6, residual = 7.696308282558869e-08
CG iteration 7, residual = 2.172768027790556e-09
CG iteration 8, residual = 6.173415107275903e-11
CG iteration 9, residual = 1.8887182888427588e-12

Draw (gfu, deformation=True, **ea);

This example shows the result of a HarmonicProlongation:

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

fes = H1(mesh, order=8, hoprolongation=True)
u,v = fes.TnT()

a = BilinearForm(grad(u)*grad(v)*dx, condense=True).Assemble()
fes.SetHarmonicProlongation(a)
pre = preconditioners.MultiGrid(a, smoother="block", blocktype="vertexpatch")

gfu = GridFunction(fes)
gfu.Set(sin(10*x))
Draw(gfu, order=5)
mesh.Refine()
Draw(gfu, order=5)
gfu.vec.data += a.harmonic_extension * gfu.vec
Draw (gfu, order=5);

Nearly incompressible elasticity and Stokes

The Scott-Vogelius element:

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

fes = VectorH1(mesh, order=4, hoprolongation=True, dirichlet=".*")
u,v = fes.TnT()

eq = InnerProduct(Grad(u),Grad(v)) * dx + 1e4 * div(u)*div(v)* dx
a = BilinearForm(eq).Assemble()

# a = BilinearForm( (Grad(u)|Grad(v)) * dx).Assemble()

pre = preconditioners.MultiGrid(a, blocktype=["vertexpatch"])

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

2 3618 0.3132568089258626 0.999900784186667
3 14146 0.19750181032339126 0.9999100509628502
4 55938 0.19549161267829734 0.9999047978306761
5 222466 0.19549539298189234 0.9998999587479052
6 887298 0.1955623207246145 0.9998955800582556

A robust method with reduced integration

We define the bilinear form as

\[A(u,v) = \int_\Omega \nabla u : \nabla v \, dx + \frac{1}{\varepsilon} \int_\Omega P_{L_2}^Q \operatorname{div} u \operatorname{div} v \, dx,\]

where \(P_{L_2}^Q\) is an element-wise \(L_2\)-projector into a lower order space.

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

fesp = L2(mesh, order=0)
fes = VectorH1(mesh, order=2, hoprolongation=True, dirichlet=".*")
u,v = fes.TnT()

eq = InnerProduct(Grad(u),Grad(v)) * dx + 1e4 * Interpolate(div(u),fesp)*div(v)* dx
a = BilinearForm(eq).Assemble()

# a = BilinearForm( (Grad(u)|Grad(v)) * dx).Assemble()
fes.SetHarmonicProlongation(a)
pre = preconditioners.MultiGrid(a, blocktype=["vertexpatch"])

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

2 946 0.5451163631463438 1.2567450360089776
3 3618 0.35880823306083937 1.4293139465317155
4 14146 0.30554010192394687 1.5320915176928378
5 55938 0.27985898251764896 1.6946834108546738

f = LinearForm((x-0.5)*v[1]*dx).Assemble()
gfu = GridFunction(fes)
Solve (a*gfu==f, pre, CGSolver, printrates=True)

GMRes iteration 1, residual = 0.23493908809833197
GMRes iteration 2, residual = 0.007800992971133687
GMRes iteration 3, residual = 0.0017934210623298725
GMRes iteration 4, residual = 0.0005261276687800984
GMRes iteration 5, residual = 0.00015073318275287846
GMRes iteration 6, residual = 5.7899132958387704e-05
GMRes iteration 7, residual = 2.6799606698463287e-05
GMRes iteration 8, residual = 1.2709143556341768e-05
GMRes iteration 9, residual = 6.0059690801987565e-06
GMRes iteration 10, residual = 2.7884138826742696e-06
GMRes iteration 11, residual = 1.1107640406255571e-06
GMRes iteration 12, residual = 4.4088565769485225e-07
GMRes iteration 13, residual = 1.79013146736778e-07
GMRes iteration 14, residual = 7.690860207089606e-08
GMRes iteration 15, residual = 3.307757519032511e-08
GMRes iteration 16, residual = 1.3612173714185154e-08
GMRes iteration 17, residual = 5.618571264560557e-09
GMRes iteration 18, residual = 2.408103474479244e-09
GMRes iteration 19, residual = 1.0124601177825438e-09
GMRes iteration 20, residual = 4.05305008489e-10
GMRes iteration 21, residual = 1.5830654940261236e-10
GMRes iteration 22, residual = 6.353606152351412e-11
GMRes iteration 23, residual = 2.5324646225872506e-11
GMRes iteration 24, residual = 9.580036115876123e-12
GMRes iteration 25, residual = 3.471055740981902e-12
GMRes iteration 26, residual = 1.2635465365592798e-12
GMRes iteration 27, residual = 4.708786735524909e-13
GMRes iteration 28, residual = 1.709423721114709e-13

Draw (gfu)

BaseWebGuiScene

from netgen.occ import *
shape = MoveTo(0,0).LineTo(1,0,"in").LineTo(1,1).LineTo(2,1).LineTo(3,0).LineTo(4,1).LineTo(5,1) \
    .LineTo(5,2,"out").LineTo(4,2).LineTo(3,1).LineTo(2,2).LineTo(1,2).Rotate(180).Arc(1,90).Close().Face()
mesh = shape.GenerateMesh(dim=2, maxh=0.25).Curve(3)
Draw (mesh)
print (mesh.GetBoundaries())

fesp = L2(mesh, order=0)
fes = VectorH1(mesh, order=2, hoprolongation=True, dirichlet="in|default")
u,v = fes.TnT()

eq = InnerProduct(Grad(u),Grad(v)) * dx + 1e4 * Interpolate(div(u),fesp)*div(v)* dx
a = BilinearForm(eq).Assemble()

# a = BilinearForm( (Grad(u)|Grad(v)) * dx).Assemble()
fes.SetHarmonicProlongation(a)
pre = preconditioners.MultiGrid(a, blocktype=["vertexpatch"])

for l in range(4):
    mesh.Refine()
    a.Assemble()
    print (mesh.levels, fes.ndof)

('in', 'default', 'default', 'default', 'default', 'default', 'out', 'default', 'default', 'default', 'default', 'default', 'default')
2 3482
3 13426
4 52706
5 208834

see this Dissertation, page 116.

f = LinearForm(fes)
gfu = GridFunction(fes)
Solve (a*gfu==f, pre, CGSolver, dirichlet=CF((0,x*(1-x))) | mesh.Boundaries("in"), printrates=True)

GMRes iteration 1, residual = 31.925481011977798
GMRes iteration 2, residual = 2.988227598636033
GMRes iteration 3, residual = 0.7390615672896171
GMRes iteration 4, residual = 0.17058254115885987
GMRes iteration 5, residual = 0.05935867997739947
GMRes iteration 6, residual = 0.02499133540663386
GMRes iteration 7, residual = 0.01406142437421026
GMRes iteration 8, residual = 0.00901487970381001
GMRes iteration 9, residual = 0.004613270577852951
GMRes iteration 10, residual = 0.0018629688936926646
GMRes iteration 11, residual = 0.0006782704227368849
GMRes iteration 12, residual = 0.0003013641480589808
GMRes iteration 13, residual = 0.00013405661307782702
GMRes iteration 14, residual = 5.9566650574728934e-05
GMRes iteration 15, residual = 2.412102486912539e-05
GMRes iteration 16, residual = 1.0668702568317427e-05
GMRes iteration 17, residual = 4.916939332289962e-06
GMRes iteration 18, residual = 2.0799757960403246e-06
GMRes iteration 19, residual = 8.381498341359801e-07
GMRes iteration 20, residual = 3.5130547815437323e-07
GMRes iteration 21, residual = 1.4994936269763095e-07
GMRes iteration 22, residual = 6.554696017925733e-08
GMRes iteration 23, residual = 3.1287151779121895e-08
GMRes iteration 24, residual = 1.3882484126552414e-08
GMRes iteration 25, residual = 5.562545781132151e-09
GMRes iteration 26, residual = 2.3376106097852802e-09
GMRes iteration 27, residual = 9.927167451062273e-10
GMRes iteration 28, residual = 4.0206267854772215e-10
GMRes iteration 29, residual = 1.66105785580514e-10
GMRes iteration 30, residual = 7.081503078766824e-11
GMRes iteration 31, residual = 2.970810703311859e-11

Draw (gfu);