Schwarz preconditioners for high order finite elements

81. 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,), (), (), (), (), (), (), (), (), (), (), (), (), (), (), (), (), (), (), (), (), (), (), (), (), (), (), (), (), (), (), (), (), (), (), (), (), (), (), (), (), (), (), (), (), (), (), (), (), (), (), (), (), (), ()]

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.03600107735785368     
CG iteration 2, residual = 0.038372133430393085     
CG iteration 3, residual = 0.029110059484005368     
CG iteration 4, residual = 0.027141112208440944     
CG iteration 5, residual = 0.017569446437950133     
CG iteration 6, residual = 0.011526175793908976     
CG iteration 7, residual = 0.008223670026623633     
CG iteration 8, residual = 0.005175756512492302     
CG iteration 9, residual = 0.0033326056989442495     
CG iteration 10, residual = 0.0021695790281934047     
CG iteration 11, residual = 0.0017192307844422971     
CG iteration 12, residual = 0.0011194899718497886     
CG iteration 13, residual = 0.0005551689787361585     
CG iteration 14, residual = 0.00033185029452326164     
CG iteration 15, residual = 0.00020548463318327418     
CG iteration 16, residual = 0.00012501023681255603     
CG iteration 17, residual = 8.249278113708203e-05     
CG iteration 18, residual = 4.194808314485239e-05     
CG iteration 19, residual = 2.774519246951208e-05     
CG iteration 20, residual = 1.629950611838762e-05     
CG iteration 21, residual = 1.0683041626248798e-05     
CG iteration 22, residual = 7.829498488416673e-06     
CG iteration 23, residual = 4.998589544127889e-06     
CG iteration 24, residual = 3.1419610329681682e-06     
CG iteration 25, residual = 1.774162749131244e-06     
CG iteration 26, residual = 1.0621529093398554e-06     
CG iteration 27, residual = 5.923827837653877e-07     
CG iteration 28, residual = 3.0868547517303654e-07     
CG iteration 29, residual = 1.6372588267151164e-07     
CG iteration 30, residual = 8.453675310914573e-08     
CG iteration 31, residual = 4.827118044289648e-08     
CG iteration 32, residual = 2.4503550791698903e-08     
CG iteration 33, residual = 1.336017892188534e-08     
CG iteration 34, residual = 6.922662690227991e-09     
CG iteration 35, residual = 3.974036650498149e-09     
CG iteration 36, residual = 2.436704957055982e-09     
CG iteration 37, residual = 1.3020788349876099e-09     
CG iteration 38, residual = 7.428345203403012e-10     
CG iteration 39, residual = 4.4932781620273573e-10     
CG iteration 40, residual = 2.1902309893115748e-10     
CG iteration 41, residual = 1.432282758621592e-10     
CG iteration 42, residual = 6.154601025499428e-11     
CG iteration 43, residual = 3.055221663315438e-11     
CG iteration 44, residual = 1.723177834905711e-11     
CG iteration 45, residual = 9.365286235132443e-12     
CG iteration 46, residual = 5.019742972592409e-12     
CG iteration 47, residual = 2.3167393994567743e-12     
CG iteration 48, residual = 1.1284855487475755e-12     
CG iteration 49, residual = 5.355145750599808e-13     
CG iteration 50, residual = 2.8434096010959235e-13     
CG iteration 51, residual = 1.3755390110821064e-13     
CG iteration 52, residual = 6.29889613302593e-14     
CG iteration 53, residual = 2.867016439654294e-14     
CG converged in 53 iterations to residual 2.867016439654294e-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.07035098554201802 2.7338478939072934

81.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

81. Schwarz preconditioners for high order finite elements#

81.1. Overlapping blocks#