Iteration methods in parallel

105. Iteration methods in parallel#

from ipyparallel import Cluster
c = await Cluster(engines="mpi").start_and_connect(n=4, activate=True)
c.ids

Starting 4 engines with <class 'ipyparallel.cluster.launcher.MPIEngineSetLauncher'>

[0, 1, 2, 3]

%%px
from mpi4py import MPI
comm = MPI.COMM_WORLD

from ngsolve import *

comm = MPI.COMM_WORLD
mesh = Mesh(unit_square.GenerateMesh(maxh=0.1, comm=comm))

setup a standard problem, in parallel. We use a Jacobi preconditioner: it extracts local diagonals, and cumulates them over identified dofs.

%%px
fes = H1(mesh, dirichlet='.*')
u,v = fes.TnT()
a = BilinearForm(grad(u)*grad(v)*dx).Assemble()
f = LinearForm(1*v*dx).Assemble()
gfu = GridFunction(fes)

pre = Preconditioner(a, "local") # Jacobi preconditioner
pre.Update()

105.1. Richardson iteration#

The right hand side vector f.vec is distributed, the vector gfu.vec of the gridfunction is consistent. Create help vectors of the same type:

  • residual calculation r = f - A u is purely local

  • the preconditioning step w = pre * r includes type conversion from a distributed input to a consistent output

  • the inner product for the error acts on two vectors of opposite type

  • solution vector update can be done purely local

%%px
r = f.vec.CreateVector()    # a distributed vector
w = gfu.vec.CreateVector()  # a consistent vector

for it in range(100):
    r.data = f.vec - a.mat * gfu.vec
    w.data = pre * r
    err = InnerProduct(w, r)
    if comm.rank==0: print (err)
    gfu.vec.data += w

[stdout:0] 0.0021064513366757624
0.0016649942356400777
0.0014293622352480324
0.0012540554438783874
0.001112404612257642
0.0009933390479475918
0.000890831942976833
0.0008011690168185088
0.0007218924281157399
0.0006512848690148516
0.0005880862634490914
0.0005313295098339835
0.00048024242401788295
0.0004341879116486988
0.0003926266502992656
0.0003550932023919776
0.00032118024145896124
0.00029052774906968394
0.0002628153108169191
0.00023775638776913217
0.00021509388242576453
0.0001945965815081428
0.00017605621543874845
0.00015928496924243752
0.00014411333724077553
0.0001303882493063521
0.00011797141848740605
0.00010673787378974326
9.657465096663952e-05
8.737962020820057e-05
7.906043377785088e-05
7.153357960521118e-05
6.472352903316176e-05
5.856196859246741e-05
5.298710700251122e-05
4.794304967465525e-05
4.33792338922704e-05
3.924991860282375e-05
3.551372341262978e-05
3.2133211945081825e-05
2.9074515223568345e-05
2.630699118224849e-05
2.380291680011571e-05
2.15372097032181e-05
1.9487176392098893e-05
1.7632284531086348e-05
1.595395698683949e-05
1.4435385528847577e-05
1.3061362307207032e-05
1.1818127405416772e-05
1.0693230930278576e-05
9.675408249165827e-06
8.754467118573381e-06
7.921185568486566e-06
7.1672195160009044e-06
6.485019179969889e-06
5.867753457298018e-06
5.309241501770909e-06
4.803890818858559e-06
4.346641255507124e-06
3.9329143232140194e-06
3.5585673462821717e-06
3.2198519756196523e-06
2.913376652277429e-06
2.636072644562713e-06
2.38516331841581e-06
2.1581363331689332e-06
1.9527184841408617e-06
1.7668529400544409e-06
1.5986786472714488e-06
1.4465116945531416e-06
1.3088284517023017e-06
1.1842503132148085e-06
1.071529894150113e-06
9.695385399768453e-07
8.77255025312514e-07
7.937553283851727e-07
7.182033788183627e-07
6.498426860896932e-07
5.879887648335622e-07
5.320222811376854e-07
4.813828512053036e-07
4.355634302857997e-07
3.9410523568941653e-07
3.565931530484706e-07
3.226515798276342e-07
2.9194066446357063e-07
2.641529034764591e-07
2.390100624797615e-07
2.1626039025863877e-07
1.9567609802136595e-07
1.7705107858376246e-07
1.6019884264875789e-07
1.4495065151739598e-07
1.3115382753416282e-07
1.1867022534926646e-07
1.073748486908884e-07
9.715459879719703e-08
8.79071419763558e-08
7.953988495554161e-08

from ngsolve.webgui import Draw

gfu = c[:]['gfu']
Draw(gfu[0]);

Very similar for other iteraton methods, such as Conjugate Gradients. The matrix operation goes from consistent to distributed without communication, the preconditioner does the cumulation. The inner products are between different vector types:

%%px
from ngsolve.krylovspace import CGSolver
inv = CGSolver(a.mat, pre, printrates=comm.rank==0)

gfu.vec.data = inv * f.vec

[stdout:0] CG iteration 1, residual = 0.0458960928258143     
CG iteration 2, residual = 0.07023551604308398     
CG iteration 3, residual = 0.050613080017471855     
CG iteration 4, residual = 0.036110304721396196     
CG iteration 5, residual = 0.024371225619145134     
CG iteration 6, residual = 0.012037144772437378     
CG iteration 7, residual = 0.003957451163626997     
CG iteration 8, residual = 0.002808522019112092     
CG iteration 9, residual = 0.0019576676071666097     
CG iteration 10, residual = 0.0012105518357746058     
CG iteration 11, residual = 0.0006678458492479842     
CG iteration 12, residual = 0.00034994343416230693     
CG iteration 13, residual = 0.00015505713813410548     
CG iteration 14, residual = 9.053904926414744e-05     
CG iteration 15, residual = 4.1301684511661955e-05     
CG iteration 16, residual = 2.681164276691173e-05     
CG iteration 17, residual = 1.3058747054232963e-05     
CG iteration 18, residual = 4.920983843351047e-06     
CG iteration 19, residual = 1.8602999450351314e-06     
CG iteration 20, residual = 1.0070126786667637e-06     
CG iteration 21, residual = 5.936469705558383e-07     
CG iteration 22, residual = 4.061853379897207e-07     
CG iteration 23, residual = 3.073398845261929e-07     
CG iteration 24, residual = 1.2811327744940765e-07     
CG iteration 25, residual = 7.645607659519623e-08     
CG iteration 26, residual = 3.576003073963308e-08     
CG iteration 27, residual = 1.3675779003025412e-08     
CG iteration 28, residual = 5.429352347662427e-09     
CG iteration 29, residual = 2.7388776138098953e-09     
CG iteration 30, residual = 1.492972650904689e-09     
CG iteration 31, residual = 6.559774367270432e-10     
CG iteration 32, residual = 2.2925027188366396e-10     
CG iteration 33, residual = 9.769065823884045e-11     
CG iteration 34, residual = 3.873546483208607e-11     
CG iteration 35, residual = 1.087987525169979e-11     
CG iteration 36, residual = 3.638827062697718e-12     
CG iteration 37, residual = 1.0472900892991234e-12     
CG iteration 38, residual = 2.905348173962424e-13     
CG iteration 39, residual = 7.585404739032898e-14     
CG iteration 40, residual = 2.59351193019184e-14     

gfu = c[:]['gfu']
Draw(gfu[0]);

c.shutdown(hub=True)
Contents