Parallel Mesh Generation Using mshr and Legacy FEniCS

General
, ,
1

Hi everyone,

I’m using FEniCS 2019.1.0 with mshr. From what I understand, mshr is serial only - however, is there a way to modify meshes generated with mshr so that I can use MPI parallelization?

Consider as an example the rather trivial program below, where we solve the Poisson equation on a patterned surface:

import fenics as fe
import numpy as np
import mshr
 
comm = fe.MPI.comm_world
rank = fe.MPI.rank(comm)

fe.parameters["form_compiler"]["optimize"] = True
fe.parameters["form_compiler"]["cpp_optimize"] = True

L_x = 1
L_y = 1

surface_amplitude = 0.
surface_freq = L_x/(8*np.pi)
def surfaceExpr(x):
    
    return surface_amplitude*np.cos(surface_freq*(x - L_x/2) )

def surfExprDeriv(x):
    return - surface_amplitude*surface_freq*np.sin(surface_freq*(x-L_x/2))

domain_n_points = 60
domain_points = []
for n in range(domain_n_points + 1):
    x = n*L_x/domain_n_points
    domain_points.append(fe.Point(x, surfaceExpr(x) ))
domain_points.append(fe.Point(L_x,  surface_amplitude))
domain_points.append(fe.Point(L_x, L_y))
domain_points.append(fe.Point(0., L_y))
domain_points.append(fe.Point(0., surface_amplitude))
domain1 = mshr.Polygon(domain_points)
if rank == 0:
    mesh = mshr.generate_mesh(domain1, domain_n_points)

#mesh = mshr.generate_mesh(domain1, domain_n_points)

V = fe.FunctionSpace(mesh, "Lagrange", 1)


# Define boundary condition
u_D = fe.Expression('1 + x[0]*x[0] + 2*x[1]*x[1]', degree=2)

def boundary(x, on_boundary):
    return on_boundary
bc = fe.DirichletBC(V, u_D, boundary)
# Define variational problem
trialFn = fe.TrialFunction(V)
testFn = fe.TestFunction(V)
f = fe.Constant(-6.0)
a = fe.dot(fe.grad(trialFn), fe.grad(testFn))*fe.dx
L = f*testFn*fe.dx
# Compute solution
u = fe.Function(V)
fe.solve(a == L, u, bc)
2

You can save the mesh to disk on rank 0 with COMM_SELF using a parallel capable format (e.g., XDMF), and then reload it for parallel computation with all ranks using COMM_WORLD.

3

This doesn’t seem to work. Consider the block of code pasted below -

import mshr

comm = MPI.comm_world
rank = MPI.rank(comm)

L_x = 1
L_y = 1
print("Entering rank-0 block", flush=True)
if rank == 0:
    Domain = mshr.Rectangle(Point(0.0, 0.0), Point(L_x, L_y))
    print("Domain has been defined", flush=True)
    
    mesh_serial = mshr.generate_mesh(Domain, 1)
    print("Mesh has been generated", flush=True)

    with XDMFFile(MPI.comm_self, meshFile) as xdmf:
        xdmf.write(mesh_serial)
    print("Mesh has been written on rank 0", flush=True)

MPI.barrier(comm)


mesh = Mesh(comm)
with XDMFFile(comm, meshFile) as xdmf:
    xdmf.read(mesh)
    
print(rank, mesh.num_cells())

Upon running this with mpirun -np N python test.py, N > 2, the program hangs before mesh_serial = mshr.Rectangle(Point(0.0, 0.0), Point(L_x, L_y)) is executed. The following is what is written to console -

Entering rank-0 block
Entering rank-0 block
Domain has been defined
4

Urgh. I have a horrible feeling using mshr may be locking you into using a default global MPI_COMM_WORLD.

If this is the case, then my proposal won’t work.

You can try to figure this out by looking at the legacy code: GitHub - martinal/mshr · GitHub

You may have to resort to generating a mesh in a separate script. That script can be run independently via the user or spawning a new process for its execution.

You could also try modifying mshr’s legacy code to require an MPI_COMM, but the GPL license may be unfriendly for you.

Ultimately I’d recommend avoiding mshr if you can since it’s not been maintained for a very very long time, and look into something like gmsh.