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

Multigrid for high order finite element spaces

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

• require simplicial meshes

• require uniform polynomial order

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

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

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

0.9737767048916046 0.9999650900459265
0.9731893341714941 0.9999569893694966
0.9738508899616769 0.9999593243766862
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

Draw (gfu, deformation=True);

High order with static condensation

mesh = Mesh(unit_square.GenerateMesh(maxh=0.2))
fes = H1(mesh, order=5, hoprolongation=True, dirichlet=".*")

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

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

1 2801 0.8873455539760278 0.9998975910633985
2 11001 0.8957401274225056 0.9998894014559885
3 43601 0.8945786078982365 0.9998863285007136
4 173601 0.894386169292487 0.9998960435342532
5 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);

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)

BaseWebGuiScene

gfu.vec.data += a.harmonic_extension * gfu.vec

Draw (gfu, order=6)

BaseWebGuiScene