Schwarz preconditioners for high order finite elements

65. Schwarz preconditioners for high order finite elements#

Basis functions for high order finite element methods (aka p-version fem) are associated with vertices, edges, faces and cells of the triangulation. Thus, there is a natural way to define block-Jacobi (or block-Gauss-Seidel) preconditioners by combining dofs sitting on the same geometric entities (called nodes).

from ngsolve import *
from ngsolve.la import EigenValues_Preconditioner
from ngsolve.webgui import Draw
mesh = Mesh(unit_square.GenerateMesh(maxh=0.2))

def DefineProblem(order, printing=True):
    fes = H1(mesh, order=order)
    u,v = fes.TnT()
    a = BilinearForm(grad(u)*grad(v)*dx+10*u*v*dx).Assemble()
    f = LinearForm(x*y*v*dx).Assemble()

    blocks = []
    freedofs = fes.FreeDofs()
    for v in mesh.vertices:
        blocks.append (fes.GetDofNrs(v))
    for ed in mesh.edges:
        blocks.append (fes.GetDofNrs(ed))
    for fa in mesh.faces:
        blocks.append (fes.GetDofNrs(fa))

    if printing:
        print (blocks)
    pre = a.mat.CreateBlockSmoother(blocks, GS=False)

    return a,f,pre

a,f,pre = DefineProblem(2)

[(0,), (1,), (2,), (3,), (4,), (5,), (6,), (7,), (8,), (9,), (10,), (11,), (12,), (13,), (14,), (15,), (16,), (17,), (18,), (19,), (20,), (21,), (22,), (23,), (24,), (25,), (26,), (27,), (28,), (29,), (30,), (31,), (32,), (33,), (34,), (35,), (36,), (37,), (38,), (39,), (40,), (41,), (42,), (43,), (44,), (45,), (46,), (47,), (48,), (49,), (50,), (51,), (52,), (53,), (54,), (55,), (56,), (57,), (58,), (59,), (60,), (61,), (62,), (63,), (64,), (65,), (66,), (67,), (68,), (69,), (70,), (71,), (72,), (73,), (74,), (75,), (76,), (77,), (78,), (79,), (80,), (81,), (82,), (83,), (84,), (85,), (86,), (87,), (88,), (89,), (90,), (91,), (92,), (93,), (94,), (95,), (96,), (97,), (98,), (99,), (100,), (101,), (102,), (103,), (104,), (105,), (106,), (107,), (108,), (109,), (110,), (111,), (112,), (113,), (114,), (115,), (116,), (117,), (118,), (119,), (120,), (121,), (122,), (123,), (124,), (125,), (126,), (127,), (128,), (129,), (130,), (131,), (132,), (), (), (), (), (), (), (), (), (), (), (), (), (), (), (), (), (), (), (), (), (), (), (), (), (), (), (), (), (), (), (), (), (), (), (), (), (), (), (), (), (), (), (), (), (), (), (), (), (), (), (), (), (), (), (), ()]

gfu = GridFunction(a.space)
from ngsolve.krylovspace import CGSolver
inv = CGSolver(mat=a.mat, pre=pre, printrates='\r', maxiter=400)
gfu.vec.data = inv * f.vec
Draw (gfu);

CG iteration 1, residual = 0.03616817538237669     
CG iteration 2, residual = 0.038578886687444205     
CG iteration 3, residual = 0.028786277945169303     
CG iteration 4, residual = 0.02717043861464373     
CG iteration 5, residual = 0.017929160922748995     
CG iteration 6, residual = 0.011843671131637522     
CG iteration 7, residual = 0.008858810330184645     
CG iteration 8, residual = 0.005512321246644729     
CG iteration 9, residual = 0.0035742193069537506     
CG iteration 10, residual = 0.0022787217503585614     
CG iteration 11, residual = 0.0018222848385735526     
CG iteration 12, residual = 0.0011428466557880732     
CG iteration 13, residual = 0.0005434954855783924     
CG iteration 14, residual = 0.0003437652718687495     
CG iteration 15, residual = 0.00017718883102318885     
CG iteration 16, residual = 0.00011375493321065046     
CG iteration 17, residual = 7.490762189941222e-05     
CG iteration 18, residual = 3.84032156624639e-05     
CG iteration 19, residual = 2.3546130827285455e-05     
CG iteration 20, residual = 1.3136911569268165e-05     
CG iteration 21, residual = 7.402546355528151e-06     
CG iteration 22, residual = 4.554507451136723e-06     
CG iteration 23, residual = 2.801171468526849e-06     
CG iteration 24, residual = 2.1678144767741274e-06     
CG iteration 25, residual = 1.5293620323637416e-06     
CG iteration 26, residual = 9.902098365795283e-07     
CG iteration 27, residual = 6.076773468362999e-07     
CG iteration 28, residual = 3.9335995335618345e-07     
CG iteration 29, residual = 2.0838424427178355e-07     
CG iteration 30, residual = 1.1715084968823221e-07     
CG iteration 31, residual = 6.653829414967276e-08     
CG iteration 32, residual = 4.229139870798433e-08     
CG iteration 33, residual = 2.188084777154004e-08     
CG iteration 34, residual = 1.2931264312198394e-08     
CG iteration 35, residual = 8.452500262063084e-09     
CG iteration 36, residual = 5.420344679804673e-09     
CG iteration 37, residual = 2.6351941581130006e-09     
CG iteration 38, residual = 1.2285806360324746e-09     
CG iteration 39, residual = 6.734731997623074e-10     
CG iteration 40, residual = 3.564004473485101e-10     
CG iteration 41, residual = 1.8274090136404232e-10     
CG iteration 42, residual = 8.808562920476457e-11     
CG iteration 43, residual = 4.7976737044634345e-11     
CG iteration 44, residual = 2.367143553464296e-11     
CG iteration 45, residual = 1.2073949988316887e-11     
CG iteration 46, residual = 5.33643118654821e-12     
CG iteration 47, residual = 2.592476359799728e-12     
CG iteration 48, residual = 1.2561165106050925e-12     
CG iteration 49, residual = 5.542308870311771e-13     
CG iteration 50, residual = 2.7131383862528984e-13     
CG iteration 51, residual = 1.1614959668268976e-13     
CG iteration 52, residual = 5.841635350994417e-14     
CG iteration 53, residual = 3.027036794325895e-14     
CG converged in 53 iterations to residual 3.027036794325895e-14

Exercise:

Plot the dependency of the extremal eigenvalues on the polynomial order \(p\).
Compare time to solution. To get comparable results, plot number of dofs per second.

a,f,pre = DefineProblem(3, False)
lam = list(EigenValues_Preconditioner(a.mat, pre))
print ("lammin, lammax=", lam[0], lam[-1])

lammin, lammax= 0.06934594552348519 2.7187250629149218

65.1. Overlapping blocks#

One can improve the condition number by letting blocks overlap. Examples:

One can assign all dofs associated with one element (including the boundary) to a block. The number of blocks is the number of elements.
One can combine all dofs on all elements sharing a vertex to blocks. The number of blocks is now the number of vertics.
One can combine all dofs on all edges sharing a vertex to blocks. The dofs on faces and cells are left as extra blocks.

Exercise:

Implement these definition of blocks.
Find some more possibilities.
Which one gives the best condition number ?
Which one computes the solution in least time ?

Consult NGSolve documentation:

The \(p\)-robust analysis for triangular and tetrahedral meshes can be found in: Additive Schwarz preconditioning for p-version triangular and tetrahedral finite elements

Schwarz preconditioners for high order finite elements

Contents

65. Schwarz preconditioners for high order finite elements#

65.1. Overlapping blocks#