On using MPI

Question

(I don't have much experience in parallel programming, so apologies if the following does not make a lot of sense or is trivial.)

I wrote a object-oriented Python program, where each PDE that I encounter is a class on its own. In a very simplified way it looks like this:

program.py

import subprocess
from dolfin import *

class Main(object):
    def __init__(self):
        self.geometry = Geometry()
        self.poissonpde = Poisson(self.geometry)

    def do_stuff(self):
        fct_list = [Expression("x[0]*x[0] + x[1]"), Expression("x[0] + 7 * x[1]")]
        for fct in fct_list:
            solution = self.poissonpde.solve(fct)
            # plot(solution, interactive=True)
            print solution.vector().array()[:10]


class Geometry(object):
    def __init__(self):
        subprocess.call(["gmsh", "-2", "square.geo"])
        subprocess.call(["dolfin-convert", "square.msh", "square.xml"])
        self.mesh = Mesh("square.xml")
        self.regions = MeshFunction('size_t', self.mesh, "square_physical_region.xml")
        self.boundaries = MeshFunction('size_t', self.mesh, "square_facet_region.xml")
        self.dx = Measure('dx')[self.regions]
        self.ds = Measure('ds')[self.boundaries]


class Poisson(object):
    def __init__(self, geometry):
        self.geometry = geometry
        self.V = FunctionSpace(self.geometry.mesh, 'CG', 1)
        self.u = TrialFunction(self.V)
        self.v = TestFunction(self.V)
        self.sol = Function(self.V)
        self.bcs = [DirichletBC(self.V, Constant(0.0), self.geometry.boundaries, 2),
                    DirichletBC(self.V, Constant(0.0), self.geometry.boundaries, 3)]
        self.solver = LUSolver()

    def solve(self, sourcefct):
        a = inner(grad(self.u), grad(self.v)) * self.geometry.dx
        L = sourcefct * self.v * self.geometry.dx
        A, L = assemble_system(a, L, self.bcs)
        self.solver.set_operator(A)
        self.solver.solve(self.sol.vector(), L)
        return self.sol


main = Main()
main.do_stuff()

square.geo

cl = 1.0;
Point(1) = {0, 0, 0, cl};
Point(2) = {0, 20, 0, cl};
Point(3) = {20, 20, 0, cl};
Point(4) = {20, 0, 0, cl};
Line(1) = {1, 2};
Line(2) = {2, 3};
Line(3) = {3, 4};
Line(4) = {4, 1};
Line(5) = {1, 3};
Line Loop(6) = {1, 2, -5};
Plane Surface(7) = {6};
Line Loop(8) = {5, 3, 4};
Plane Surface(9) = {8};
Physical Line(1) = {1, 3};
Physical Line(2) = {2};
Physical Line(3) = {4};
Physical Surface(0) = {7};
Physical Surface(1) = {9};

... where in practice, this code consists of multiple modules. Now I would like to speed things up using MPI, but the examples all assume one python script.

Hence my question what is the recommended way to solve the PDE's using MPI, with as few modifications as possible?

---

EDIT Based on the discussion below, splitting program.py into the following two files works (running serial.py first, then parallel.py):

serial.py

import subprocess
from dolfin import *

subprocess.call(["gmsh", "-2", "square.geo"])
subprocess.call(["dolfin-convert", "square.msh", "square.xml"])
mesh = Mesh("square.xml")
regions = MeshFunction('size_t', mesh, "square_physical_region.xml")
boundaries = MeshFunction('size_t', mesh, "square_facet_region.xml")
hdf = HDF5File(mesh.mpi_comm(), "file.h5", "w")
hdf.write(mesh, "/mesh")
hdf.write(regions, "/regions")
hdf.write(boundaries, "/boundaries")

parallel.py

from dolfin import *

class Main(object):
    def __init__(self):
        self.geometry = Geometry()
        self.poissonpde = Poisson(self.geometry)

    def do_stuff(self):
        fct_list = [Expression("x[0]*x[0] + x[1]"), Expression("x[0] + 7 * x[1]")]
        for fct in fct_list:
            solution = self.poissonpde.solve(fct)
            # plot(solution, interactive=True)
            print solution.vector().array()[:10]


class Geometry(object):
    def __init__(self):
        self.mesh = Mesh()
        hdf = HDF5File(self.mesh.mpi_comm(), "file.h5", "r")
        hdf.read(self.mesh, "/mesh", False)
        self.regions = CellFunction("size_t", self.mesh)
        hdf.read(self.regions, "/regions")
        self.boundaries = FacetFunction("size_t", self.mesh)
        hdf.read(self.boundaries, "/boundaries")
        self.dx = Measure('dx')[self.regions]
        self.ds = Measure('ds')[self.boundaries]

class Poisson(object):
    def __init__(self, geometry):
        self.geometry = geometry
        self.V = FunctionSpace(self.geometry.mesh, 'CG', 1)
        self.u = TrialFunction(self.V)
        self.v = TestFunction(self.V)
        self.sol = Function(self.V)
        self.bcs = [DirichletBC(self.V, Constant(0.0), self.geometry.boundaries, 2),
                    DirichletBC(self.V, Constant(0.0), self.geometry.boundaries, 3)]
        self.solver = LUSolver()

    def solve(self, sourcefct):
        a = inner(grad(self.u), grad(self.v)) * self.geometry.dx
        L = sourcefct * self.v * self.geometry.dx
        A, L = assemble_system(a, L, self.bcs)
        self.solver.set_operator(A)
        self.solver.solve(self.sol.vector(), L)
        return self.sol


main = Main()
main.do_stuff()

Answer 1

If your real code is close enough to the simplified version you showed (mesh import, PDE solution, solution export) mpi should work without any modification, just running mpirun -n 2 python script.py.

Comments

In that case it will just run everything twice, no? I simplified the code here quite a bit; a lot of code should only be executed once (e.g. creating and writing to files). Ideally only solving the PDE's happens in parallel.

---

I changed the example above to a minimal working example: just put both files in the same directory and it will work (but not in parallel).

---

Your code works fine in parallel if you omit the system calls (you may need to switch loading your mesh as dolfin XML to HDF in newer versions of FEniCS though to read the regions/boundaries properly in parallel though).

You might be confused by the fact that the output is different from running on a single core. This is because parallelization is performed by distributing nodes across separate cores, and so each process will only have access to a piece of the solution vector (so when you print solution.vector().array()[:10], you will get several different results). Try plotting the solution using plot(solution), it should make this behaviour more clear.

If you want to analyze the solution as a whole, you can output to a file (such as a .pvd) and visualize/postprocess in Paraview, or use the various MPI-related functions in FEniCS to shuttle/interact with data across multiple cores (although this is not as trivial as the first option).

---

Thanks, I updated the code above and it works now. Some questions remain, though:

• In the part that is not executed in parallel, what it the role of the communicator in hdf = HDF5File(mesh.mpi_comm(), "file.h5", "w") ?

• To make things complicated, the boundary conditions in my problem are solutions of an eigenvalue problem on (part of) the boundarymesh. It is probably easier if I can keep this part as it is now, i.e. let every node solve this eigenvalue problem again on the entire mesh. Can I ask dolfin to work with a distributed mesh and a complete one at the same time, for different PDE's?

• I need access to the entire solution for post-processing in my code (e.g. taking the integral of the solution over a domain). What methods are there to do so?

• Given all the problems that I anticipate in doing so, is there an easier way that runs the code in serie and spawns a new process every time when PETSc is called?

---

For the first bullet, in serial.py mpi_comm() will be only one processor.
For the third bullet, I think integration of the entire domain should work.
For the second bullet, can you open a new question for this with a very minimal example? Thanks

---

Regarding the second bullet point: it was quite a challenge to use functions defined on a boundary (i.e. on a submesh of a boundarymesh) in the boundary conditions of the problem. It involved manually mapping the degrees of freedom between meshes. I doubt a minimal example is possible.

To make it a bit more clear what I am after, assume there would be another class in the example above:

class EigenvalueProblem(object):
    def __init__(self, geometry):
        self.geometry = geometry
        self.V = FunctionSpace(self.geometry.serial_mesh, 'CG', 1)
        self.u = TrialFunction(self.V)
        self.v = TestFunction(self.V)
        self.bcs = [DirichletBC(self.V, Constant(0.0), self.geometry.boundaries, 2),
                    DirichletBC(self.V, Constant(0.0), self.geometry.boundaries, 3)]
        self.A = df.PETScMatrix()
        self.B = df.PETScMatrix()
        self.eigensolver = df.SLEPcEigenSolver(self.A, self.B)

    def solve(self, sourcefct):
        a = inner(grad(self.u), grad(self.v)) * self.geometry.dx
        B = sourcefct * self.u * self.v * self.geometry.dx
        A = assemble(a, tensor=self.A)
        B = assemble(a, tensor=self.B)
        for bc in bcs:
            bc.apply(A)
        self.eigensolver.solve(1)
        val, ival, vect, ivect = self.eigensolver.get_eigenpair(0)
        return val, vect

Note that here I use self.geometry.serial_mesh to define a functionspace; in Poisson I would use self.geometry.distributed_mesh instead.
I am not sure if this is possible, but the aim is that each node computes the same eigenvalue problem again entirely, whereas the Poisson problem will be divided among nodes.