MPI and Subdomains

mesh
1

I am trying do fluid structure interaction simulations, for which I use the set_subdomain function when creating my mesh in order to ensure that the boundaries of the mesh coincide with the boundaries between the fluid and solid domains. However, if I try to run this using (for example) mpirun -np 2 python3 script.py I get the following error:

*** -------------------------------------------------------------------------

*** Error: Unable to extract local mesh data.

*** Reason: Marked subdomains are not yet supported.

*** Where: This error was encountered inside LocalMeshData.cpp.

*** Process: 0

***

*** DOLFIN version: 2018.1.0

*** Git changeset: 948dc42cc4e06ed9227d0201ad50f94ac94cbf9f

*** -------------------------------------------------------------------------

A minimal example that produces this error, but works if it is run with python3 script.py is below:

from dolfin import *
from mshr import *

L, H = 2.5, 0.41
l, h = 0.35, 0.02
C, r = (0.2,0.2), 0.05
SEGMENTS = 25
RESOLUTION = 50
domain = Rectangle(Point(0.0,0.0),Point(L,H)) - \
        Circle(Point(C[0],C[1]),r,SEGMENTS)
domain.set_subdomain(1,Rectangle(Point(C[0],C[1]-h/2), \
        Point(C[0]+r+l,C[1]+h/2)))
mesh = generate_mesh(domain, RESOLUTION)

I am running the latest stable version of FEnICS available on docker, using a mac. Is there any plan to implement subdomains for MPI? Or does anyone know of a workaround I could use?

2

If it is okay to generate the mesh in serial, then solve in parallel with a different script, you could try something like the following. First, run this in serial:

from dolfin import *
from mshr import *

L, H = 2.5, 0.41
l, h = 0.35, 0.02
C, r = (0.2,0.2), 0.05
SEGMENTS = 25
RESOLUTION = 50
domain = Rectangle(Point(0.0,0.0),Point(L,H)) - \
        Circle(Point(C[0],C[1]),r,SEGMENTS)
domain.set_subdomain(1,Rectangle(Point(C[0],C[1]-h/2), \
        Point(C[0]+r+l,C[1]+h/2)))
mesh = generate_mesh(domain, RESOLUTION)

# Convert the markers into a `MeshFunction` that can be used in parallel:
d = mesh.topology().dim()
mf = MeshFunction("size_t",mesh,d,0)
for cell in cells(mesh):
    mf[cell] = mesh.domains().markers(d)[cell.index()]

# Write out both the `Mesh` and the `MeshFunction`:
File("mesh.xml") << mesh
File("mf.xml") << mf

# (It is better to use, say, HDF5 instead of xml, but I don't remember the
# API off the top of my head.)

Then, run this in parallel:

from dolfin import *

# Read the `Mesh` and `MeshFunction` back in, possibly in parallel:
mesh = Mesh("mesh.xml")
mf = MeshFunction("size_t",mesh,"mf.xml")

# Set the `dx` volume `Measure` to use the `MeshFunction` as subdomain data:
dx = dx(subdomain_data=mf)

# Test:  Results are independent of number of processors.
area0 = assemble(Constant(1.0)*dx(0,domain=mesh))
area1 = assemble(Constant(1.0)*dx(1,domain=mesh))
area = assemble(Constant(1.0)*dx(domain=mesh)) 
if(MPI.rank(MPI.comm_world)==0):
    print("Area of region 0: "+str(area0))
    print("Area of region 1: "+str(area1))
    print("Total area: "+str(area1))
    print("Difference of total and marked: "+str(area0 + area1 - area))
1 Like