Hi,
I’m trying to “change” the comm on a submesh from the comm it inherits from the full mesh to a split (or otherwise smaller sub-) comm. The full mesh has been partitioned by region so that certain ranks are guaranteed to have no cells (ghost mode none) and I’d like the submesh comm to only include those ranks with cells. What is the best way of doing this?
I include a somewhat contrived minimal almost working example below, run using the python API with dolfinx 0.8.0 but I can easily upgrade if necessary:
- imports
import gmsh
import dolfinx as df
import numpy as np
from mpi4py import MPI
comm_world = MPI.COMM_WORLD
- create a gmsh mesh
# Create a gmsh mesh of a unit square with various regions marked as physical groups
gmsh.initialize()
gmsh.model.add("square")
lc = 3e-2
gmsh.model.geo.addPoint(0, 0, 0, lc, 1)
gmsh.model.geo.addPoint(1, 0, 0, lc, 2)
gmsh.model.geo.addPoint(1, 0.4, 0, lc, 3)
gmsh.model.geo.addPoint(1, 0.8, 0, lc, 4)
gmsh.model.geo.addPoint(1, 1, 0, lc, 5)
gmsh.model.geo.addPoint(0.5, 1, 0, lc, 6)
gmsh.model.geo.addPoint(0, 1, 0, lc, 7)
gmsh.model.geo.addPoint(0, 0.8, 0, lc, 8)
gmsh.model.geo.addPoint(0, 0.4, 0, lc, 9)
gmsh.model.geo.addPoint(0.5, 0.4, 0, lc, 10)
gmsh.model.geo.addPoint(0.5, 0.8, 0, lc, 11)
gmsh.model.geo.addLine(1, 2, 1)
gmsh.model.geo.addLine(2, 3, 2)
gmsh.model.geo.addLine(3, 4, 3)
gmsh.model.geo.addLine(4, 5, 4)
gmsh.model.geo.addLine(5, 6, 5)
gmsh.model.geo.addLine(6, 7, 6)
gmsh.model.geo.addLine(7, 8, 7)
gmsh.model.geo.addLine(8, 9, 8)
gmsh.model.geo.addLine(9, 1, 9)
gmsh.model.geo.addLine(9, 10, 10)
gmsh.model.geo.addLine(10, 3, 11)
gmsh.model.geo.addLine(10, 11, 12)
gmsh.model.geo.addLine(8, 11, 13)
gmsh.model.geo.addLine(11, 4, 14)
gmsh.model.geo.addLine(11, 6, 15)
gmsh.model.geo.addCurveLoop([1, 2, -11, -10, 9], 1)
gmsh.model.geo.addCurveLoop([11, 3, -14, -12], 2)
gmsh.model.geo.addCurveLoop([14, 4, 5, -15], 3)
gmsh.model.geo.addCurveLoop([13, 15, 6, 7], 4)
gmsh.model.geo.addCurveLoop([10, 12, -13, 8], 5)
gmsh.model.geo.addPlaneSurface([1], 1)
gmsh.model.geo.addPlaneSurface([2], 2)
gmsh.model.geo.addPlaneSurface([3], 3)
gmsh.model.geo.addPlaneSurface([4], 4)
gmsh.model.geo.addPlaneSurface([5], 5)
gmsh.model.geo.synchronize()
gmsh.model.addPhysicalGroup(1, [7, 8, 9], 1)
gmsh.model.addPhysicalGroup(1, [2, 3, 4], 2)
gmsh.model.addPhysicalGroup(1, [1], 3)
gmsh.model.addPhysicalGroup(1, [5, 6], 4)
gmsh.model.addPhysicalGroup(2, [1], 1)
gmsh.model.addPhysicalGroup(2, [2], 2)
gmsh.model.addPhysicalGroup(2, [3], 3)
gmsh.model.addPhysicalGroup(2, [4], 4)
gmsh.model.addPhysicalGroup(2, [5], 5)
# Generate the mesh on rank 0
if comm_world.rank == 0:
gmsh.model.mesh.generate(2)
- get cell region id information
# Get the cell information from the gmsh model
if comm_world.rank == 0:
topologies = df.io.gmshio.extract_topology_and_markers(gmsh.model)
num_cell_types = len(topologies.keys())
cell_information = dict()
cell_dimensions = np.zeros(num_cell_types, dtype=np.int32)
for i, element in enumerate(topologies.keys()):
_, dim, _, num_nodes, _, _ = gmsh.model.mesh.getElementProperties(element)
cell_information[i] = {"id": element, "dim": dim, "num_nodes": num_nodes}
cell_dimensions[i] = dim
# Sort elements by ascending dimension
perm_sort = np.argsort(cell_dimensions)
cell_id = cell_information[perm_sort[-1]]["id"]
num_nodes = cell_information[perm_sort[-1]]["num_nodes"]
num_nodes = comm_world.bcast([num_nodes], root=0)
cells = np.asarray(topologies[cell_id]["topology"], dtype=np.int64)
cell_values = np.asarray(topologies[cell_id]["cell_data"], dtype=np.int32)
else:
num_nodes = comm_world.bcast([None], root=0)[0]
cells = np.empty([0, num_nodes], dtype=np.int32)
cell_values = np.empty((0,), dtype=np.int32)
- setup a crude region partitioner
# convert an array to an adjacency list
def to_adj(cells, dtype):
cflat = []
coff = [0]
for c in cells:
cflat += list(c)
cc = coff[-1] + len(c)
coff += [cc]
adj = df.graph.adjacencylist(np.array(cflat, dtype=dtype),
np.array(coff, dtype=dtype))
return adj
# choose which region ids to group together (if sufficient ranks are available)
rid_groups = [[1], [2, 5], [3, 4]]
# create a partitioner that partitions the mesh depending on the cell region id
def partitioner(comm, commsize, celltype, topo):
if commsize == 1:
# trivial case
al = df.cpp.graph.AdjacencyList_int32(np.zeros(len(topo), dtype=np.int32))
elif commsize == 2:
# special case: not a large enough comm so
# group the last two rid groups together
groups = [rid_groups[0], rid_groups[1] + rid_groups[2]]
procs = np.asarray([i for v in cell_values for i, g in enumerate(groups) if v in g], dtype=np.int32)
al = df.cpp.graph.AdjacencyList_int32(procs)
else:
# arbitrary case
groups = rid_groups
if comm_world.rank == 0:
# figure out which cells belong together
cellindices = [None]*len(topo)
cellgroups = [[] for group in groups]
i = 0
for g, group in enumerate(groups):
for c, v in enumerate(cell_values):
if v in group:
cellgroups[g].append(topo.links(c))
cellindices[c] = i # remember the order
i += 1
groupsizes = [max(1, int(len(cellgroup)*commsize/len(topo))) for cellgroup in cellgroups]
# correct the logic above so that sum(groupsizes) == commsize
while sum(groupsizes) < commsize:
for g in [1, 0, 2]:
groupsizes[g] = groupsizes[g] + 1
if sum(groupsizes) == commsize: break
while sum(groupsizes) > commsize:
for g in [2, 0, 1]:
if groupsizes[g] > 1: groupsizes[g] = groupsizes[g] - 1
if sum(groupsizes) == commsize: break
# communicate the groupsizes
groupsizes = comm.bcast(groupsizes, root=0)
else:
cellgroups = [[] for group in groups]
cellindices = []
groupsizes = comm.bcast([None]*len(groups), root=0)
if sum(groupsizes) != commsize:
raise Exception("sum(groupsizes) != commsize")
values = np.empty((0,), dtype=np.int32)
# for each group, partition the cells using a cell partitioner
for g in range(len(groups)):
lal = df.mesh.create_cell_partitioner(df.mesh.GhostMode.none)(
comm, groupsizes[g], celltype, to_adj(cellgroups[g], np.int64))
values = np.concatenate((values, lal.array + sum(groupsizes[:g])), axis=0)
al = df.graph.adjacencylist(values[cellindices])
return al
- generate the mesh
# convert the gmsh mesh to a dolfinx mesh
mesh, ct, ft = df.io.gmshio.model_to_mesh(gmsh.model, comm_world, 0, gdim=2, partitioner=partitioner)
At this stage the mesh (and on 3 processes the mesh per process) looks like:
- create submeshes:
tdim = mesh.topology.dim
cells_0 = np.concatenate([ct.find(rid) for rid in rid_groups[0]])
submesh_0, _, _, _ = df.mesh.create_submesh(mesh, tdim, cells_0)
cells_1 = np.concatenate([ct.find(rid) for rid in rid_groups[1]])
submesh_1, _, _, _ = df.mesh.create_submesh(mesh, tdim, cells_1)
All the above works, and I’d now like to create submeshes that are on different comms, which include subsets of comm_world. Something like:
7. split comms:
comm_0 = comm_world.Split(submesh_0.topology.index_map(2).size_local>0)
comm_1 = comm_world.Split(submesh_1.topology.index_map(2).size_local>0)
and then create a new (geometrically and topologically identical) submesh_0 with comm comm_0 and similarly for submesh_1 with comm_1.
What would be the best way of doing this?
Many thanks!
Cian