This page was generated from unit-5a.1-mpi/poisson_mpi.ipynb.

5.1 Poisson Equation in Parallel

NGSolve can be executed on a cluster using the MPI message passing interface. You can download poisson_mpi.py and run it as

mpirun -np 4 python3 poisson_mpi.py

This computes and saves a solution which we can visualize it with drawsolution.py

netgen drawsolution.py

For proper parallel execution, Netgen/NGSolve must be configured with -DUSE_MPI=ON. Recent binaries for Linux and Mac are built with MPI support (?). If you are unsure, when your Netgen/NGSolve supports MPI, look for outputs like "Including MPI version 3.1" during Netgen startup.

MPI-parallel execution using ipyparallel

For the MPI jupyter-tutorials we use ipyparallel module. Please consult https://ipyparallel.readthedocs.io for installation.

We start an MPI server:

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

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

We use mpi4py https://mpi4py.readthedocs.io/ for issuing MPI calls from Python. The %%px syntax magic causes parallel execution of that cell:

%%px
from mpi4py.MPI import COMM_WORLD as comm
print (comm.rank, comm.size)

[stdout:1] 1 4

[stdout:3] 3 4

[stdout:0] 0 4

[stdout:2] 2 4

The master process (rank==0) generates the mesh, and distributes it within the group of processes defined by the communicator. All other ranks receive a part of the mesh. The function mesh.GetNE(VOL) returns the local number of elements:

%%px
from ngsolve import *
from netgen.occ import unit_square

if comm.rank == 0:
    mesh = Mesh(unit_square.GenerateMesh(maxh=0.1).Distribute(comm))
else:
    mesh = Mesh(netgen.meshing.Mesh.Receive(comm))
print (mesh.GetNE(VOL))

[stdout:3] 57

[stdout:1] 55

[stdout:0] 56

[stdout:2] 62

We can define spaces, bilinear / linear forms, and gridfunctions in the same way as in sequential mode. But now, the degrees of freedom are distributed on the cluster following the distribution of the mesh. The finite element spaces define how the dofs match together.

%%px
fes = H1(mesh, order=3, dirichlet=".*")
u,v = fes.TnT()

a = BilinearForm(grad(u)*grad(v)*dx)
pre = Preconditioner(a, "local")
a.Assemble()

f = LinearForm(1*v*dx).Assemble()
gfu = GridFunction(fes)

from ngsolve.krylovspace import CGSolver
inv = CGSolver(a.mat, pre.mat, printrates=comm.rank==0, maxiter=200, tol=1e-8)
gfu.vec.data = inv*f.vec

[stdout:0] CG iteration 1, residual = 0.05313636107212468
CG iteration 2, residual = 0.07324483851673458
CG iteration 3, residual = 0.06077337326925861
CG iteration 4, residual = 0.049656759272644624
CG iteration 5, residual = 0.04628399294099896
CG iteration 6, residual = 0.03051563331540459
CG iteration 7, residual = 0.025898317569357752
CG iteration 8, residual = 0.016994550808865836
CG iteration 9, residual = 0.013010287656153195
CG iteration 10, residual = 0.01143208876099897
CG iteration 11, residual = 0.007816924495134901
CG iteration 12, residual = 0.0035098066609772533
CG iteration 13, residual = 0.0017387130849062845
CG iteration 14, residual = 0.001150230939682908
CG iteration 15, residual = 0.0008184957095722701
CG iteration 16, residual = 0.000541963076272078
CG iteration 17, residual = 0.00038670700524591587
CG iteration 18, residual = 0.00028367421577579023
CG iteration 19, residual = 0.0001777343614941448
CG iteration 20, residual = 0.00013341157334192987
CG iteration 21, residual = 9.379683067760162e-05
CG iteration 22, residual = 6.550835827505457e-05
CG iteration 23, residual = 5.123579298412112e-05
CG iteration 24, residual = 3.5362251474140825e-05
CG iteration 25, residual = 2.2782716580339947e-05
CG iteration 26, residual = 1.567745927863883e-05
CG iteration 27, residual = 1.1117149587134473e-05
CG iteration 28, residual = 7.709451159725371e-06
CG iteration 29, residual = 5.474421176425092e-06
CG iteration 30, residual = 3.68767966574082e-06
CG iteration 31, residual = 2.0833695750954497e-06
CG iteration 32, residual = 1.6564423874703208e-06
CG iteration 33, residual = 1.3736999959996538e-06
CG iteration 34, residual = 8.539631243750991e-07
CG iteration 35, residual = 4.6174975273231643e-07
CG iteration 36, residual = 2.9507200081461955e-07
CG iteration 37, residual = 1.9900657697864932e-07
CG iteration 38, residual = 1.4268161595535362e-07
CG iteration 39, residual = 9.868137345841459e-08
CG iteration 40, residual = 7.020839769088627e-08
CG iteration 41, residual = 4.5770714733282975e-08
CG iteration 42, residual = 2.962279920594035e-08
CG iteration 43, residual = 1.9267023819983864e-08
CG iteration 44, residual = 1.368498086053997e-08
CG iteration 45, residual = 9.117459309938521e-09
CG iteration 46, residual = 5.918193053148567e-09
CG iteration 47, residual = 3.739344480179684e-09
CG iteration 48, residual = 2.4133694034436254e-09
CG iteration 49, residual = 1.7270530802503105e-09
CG iteration 50, residual = 1.1180335858497373e-09
CG iteration 51, residual = 7.628606214875802e-10
CG iteration 52, residual = 5.294659179776453e-10

The solution is distributed over the processes. Gather() is a collective operation that assembles the whole mesh and solution on the process with rank=0, while pickling (what c[:][...] does) transfers the local parts:

%%px
gfu_global = gfu.Gather()   # collective: whole solution on rank 0, None elsewhere

We pull the gathered solution from rank 0 to the client and draw it:

from ngsolve.webgui import Draw
Draw (c[0]["gfu_global"]);

Pulling gfu from all engines gives the local parts; drawing gfu[n] shows the part that the process with rank=n possesses.

gfu = c[:]["gfu"]   # the local parts
Draw (gfu[3]);

We can also visualize the sub-domains obtained by the automatic partitioning, without using any computed solution, as follows.

%%px
fesL2 = L2(mesh, order=0)
gfL2 = GridFunction(fesL2)
gfL2.vec.local_vec[:] = comm.rank

%%px
gfL2_global = gfL2.Gather()

Draw (c[0]["gfL2_global"]);