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.9737767045911394 0.9999650900396271
3 11001 0.9731893341903823 0.9999569893696165
4 43601 0.9738508899391483 0.9999593243765181
5 173601 0.9740047338843161 0.9999581308399157

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.5006659907705164
CG iteration 2, residual = 0.024383987764233152
CG iteration 3, residual = 0.0001879979181935111
CG iteration 4, residual = 1.1484794277924998e-06
CG iteration 5, residual = 9.252209915639308e-09
CG iteration 6, residual = 6.018344897359637e-11
CG iteration 7, residual = 4.580227880382381e-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.8873455539760282 0.9998975910633994
3 11001 0.8957401274225065 0.9998894014559885
4 43601 0.8945786078982361 0.999886328500714
5 173601 0.8943861692924864 0.9998960435342534
6 692801 0.8953852278217929 0.9998858093863481

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.4886087977865685
CG iteration 2, residual = 0.1234544421908208
CG iteration 3, residual = 0.003639088091883737
CG iteration 4, residual = 9.080367025015582e-05
CG iteration 5, residual = 2.783128212652333e-06
CG iteration 6, residual = 7.696308282558781e-08
CG iteration 7, residual = 2.172768027790504e-09
CG iteration 8, residual = 6.173415107276088e-11
CG iteration 9, residual = 1.888718288842808e-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.3132568089255071 0.9999007841868173
3 14146 0.1975018103233654 0.9999100509629262
4 55938 0.19549161267820164 0.9999047978307516
5 222466 0.19549539298198004 0.9998999587479189
6 887298 0.19556232072463403 0.9998955800582647

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.5451163631461933 1.2567450360086476
3 3618 0.3588082330609082 1.4293139465314442
4 14146 0.3055401019239069 1.5320915176938796
5 55938 0.2798589825176764 1.6946834108509634

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

CG iteration 1, residual = 0.020553570901506772
CG iteration 2, residual = 0.0024505017135056824
CG iteration 3, residual = 0.0006507984597103664
CG iteration 4, residual = 0.00020038866210370223
CG iteration 5, residual = 7.105863986308642e-05
CG iteration 6, residual = 3.0287581801163978e-05
CG iteration 7, residual = 1.3957421251276822e-05
CG iteration 8, residual = 6.026536683803503e-06
CG iteration 9, residual = 2.7439210373692124e-06
CG iteration 10, residual = 1.218630532465995e-06
CG iteration 11, residual = 4.964351575231132e-07
CG iteration 12, residual = 1.9398388083486507e-07
CG iteration 13, residual = 7.68667870228005e-08
CG iteration 14, residual = 3.304683227230027e-08
CG iteration 15, residual = 1.4352549305528452e-08
CG iteration 16, residual = 6.0211650878863604e-09
CG iteration 17, residual = 2.5562823884925753e-09
CG iteration 18, residual = 1.0915843646826387e-09
CG iteration 19, residual = 4.519581656232971e-10
CG iteration 20, residual = 1.7588626620529382e-10
CG iteration 21, residual = 6.957906924198661e-11
CG iteration 22, residual = 2.850210754902431e-11
CG iteration 23, residual = 1.1345979806610844e-11
CG iteration 24, residual = 4.16916211936322e-12
CG iteration 25, residual = 1.5090611127106195e-12
CG iteration 26, residual = 5.394782139902078e-13
CG iteration 27, residual = 1.987123416729369e-13
CG iteration 28, residual = 7.249987931852906e-14
CG iteration 29, residual = 2.7920178284432005e-14
CG iteration 30, residual = 1.1077828823726836e-14

Draw (gfu);

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

from Dissertation J. Schöberl, page 116.

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

CG iteration 1, residual = 156.0216977186489
CG iteration 2, residual = 4.696993518461411
CG iteration 3, residual = 0.441269791601204
CG iteration 4, residual = 0.06908617512251757
CG iteration 5, residual = 0.029433444799943795
CG iteration 6, residual = 0.014480167237856648
CG iteration 7, residual = 0.008734149812429095
CG iteration 8, residual = 0.005471780853158728
CG iteration 9, residual = 0.002040640466466324
CG iteration 10, residual = 0.000867522773558759
CG iteration 11, residual = 0.0003033426394281592
CG iteration 12, residual = 0.00014955848925266905
CG iteration 13, residual = 5.951615094306691e-05
CG iteration 14, residual = 2.711796490136618e-05
CG iteration 15, residual = 1.04423570214825e-05
CG iteration 16, residual = 5.09553102170522e-06
CG iteration 17, residual = 2.3425876326679256e-06
CG iteration 18, residual = 9.935176562467443e-07
CG iteration 19, residual = 3.8165395459728584e-07
CG iteration 20, residual = 1.5921102516466537e-07
CG iteration 21, residual = 7.00246203208351e-08
CG iteration 22, residual = 3.024142833369374e-08
CG iteration 23, residual = 1.4336012132499179e-08
CG iteration 24, residual = 6.296140098886312e-09
CG iteration 25, residual = 2.5090395548228084e-09
CG iteration 26, residual = 1.0887999980541186e-09
CG iteration 27, residual = 4.765297122863209e-10
CG iteration 28, residual = 2.1085634936610977e-10
CG iteration 29, residual = 1.1199062981922105e-10

Draw (gfu);