Iteration methods in parallel

103. 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()

103.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.002131411608940556
0.0016785598117133912
0.0014383176229297977
0.0012601136081994502
0.0011163969057570045
0.0009957540133089201
0.000892013039583751
0.0008013860460216729
0.0007213664248458749
0.0006501978632935307
0.0005865868529358019
0.0005295374859014327
0.00047825313447910644
0.0004320763401729044
0.0003904510870202349
0.00035289842915932834
0.00031900021889388745
0.0002883878366131279
0.0002607340736402399
0.00023574705190481803
0.00021316549799011987
0.00019275494792186845
0.00017430461490466548
0.00015762474697784634
0.00014254435985880777
0.00012890926660311873
0.00011658034877248673
0.00010543202872418053
9.53509125215299e-05
8.623457969303072e-05
7.79905007837357e-05
7.053506705294771e-05
6.3792719216495e-05
5.769516408736931e-05
5.218066951299397e-05
4.71934292567164e-05
4.268299050241807e-05
3.860373752689044e-05
3.4914425821440545e-05
3.1577761578216526e-05
2.856002200721243e-05
2.5830712431211527e-05
2.336225652920634e-05
2.112971647279069e-05
1.9110540031668575e-05
1.7284332019465414e-05
1.563264771416978e-05
1.4138806122733465e-05
1.2787721169956143e-05
1.1565749080625242e-05
1.0460550393450774e-05
9.460965197737418e-06
8.556900320854897e-06
7.739227318013283e-06
6.999690227069895e-06
6.330822151353702e-06
5.7258698239108366e-06
5.178725388153858e-06
4.6838647035290175e-06
4.236291551289987e-06
3.831487175498001e-06
3.4653646485897776e-06
3.1342275998439424e-06
2.8347328893352102e-06
2.5638568499710354e-06
2.318864756351698e-06
2.097283211876788e-06
1.896875175059282e-06
1.715617372711863e-06
1.551679871817e-06
1.4034076037199845e-06
1.2693036540219163e-06
1.148014149399825e-06
1.0383145887165753e-06
9.390974803808237e-07
8.493611611120727e-07
7.681996832028874e-07
6.947936681618619e-07
6.284020343837914e-07
5.683545153190265e-07
5.140448925990265e-07
4.649248757945639e-07
4.204985670127125e-07
3.8031745444511146e-07
3.439758843202849e-07
3.111069655442327e-07
2.8137886568143383e-07
2.54491460880751e-07
2.3017330592352966e-07
2.0817889380370092e-07
1.8828617717300159e-07
1.7029432662799943e-07
1.540217032071367e-07
1.3930402462835032e-07
1.259927067541289e-07
1.1395336353993158e-07
1.0306445032160044e-07
9.321603674500005e-08
8.430869694970844e-08
7.625250580230147e-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.046167213571327394     
CG iteration 2, residual = 0.07053899603094559     
CG iteration 3, residual = 0.05070470124514234     
CG iteration 4, residual = 0.0360268233935707     
CG iteration 5, residual = 0.023620592653494303     
CG iteration 6, residual = 0.012187751064729127     
CG iteration 7, residual = 0.0048687652674056035     
CG iteration 8, residual = 0.0036905492669662104     
CG iteration 9, residual = 0.002385951928666893     
CG iteration 10, residual = 0.0014958478667025685     
CG iteration 11, residual = 0.0007978295322313949     
CG iteration 12, residual = 0.00038943081125074753     
CG iteration 13, residual = 0.00020285820875469044     
CG iteration 14, residual = 0.00010870660188901054     
CG iteration 15, residual = 5.770924855637092e-05     
CG iteration 16, residual = 3.4406636395646876e-05     
CG iteration 17, residual = 1.3387850598294938e-05     
CG iteration 18, residual = 5.2616147425749304e-06     
CG iteration 19, residual = 1.7789183584086063e-06     
CG iteration 20, residual = 8.712122195637219e-07     
CG iteration 21, residual = 3.944919015348182e-07     
CG iteration 22, residual = 2.720041089537677e-07     
CG iteration 23, residual = 2.1583250244859396e-07     
CG iteration 24, residual = 9.96610009288972e-08     
CG iteration 25, residual = 5.690711212037661e-08     
CG iteration 26, residual = 3.226654284810232e-08     
CG iteration 27, residual = 1.3834878734993238e-08     
CG iteration 28, residual = 5.759857607292753e-09     
CG iteration 29, residual = 2.4857074071786705e-09     
CG iteration 30, residual = 1.1436471199431326e-09     
CG iteration 31, residual = 4.800144211961727e-10     
CG iteration 32, residual = 2.2081699356439377e-10     
CG iteration 33, residual = 1.0825055599990776e-10     
CG iteration 34, residual = 3.4996245094746733e-11     
CG iteration 35, residual = 1.50355291402918e-11     
CG iteration 36, residual = 4.315483334952038e-12     
CG iteration 37, residual = 1.5661479455200578e-12     
CG iteration 38, residual = 5.460411787704451e-13     
CG iteration 39, residual = 2.256368936213174e-13     
CG iteration 40, residual = 7.046557116957069e-14     
CG iteration 41, residual = 2.471853134469647e-14     

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

c.shutdown(hub=True)
Contents