Parallel usage of vertex_to_dof_map

1

Hello,

I have an issue which I can’t figure out why it might be happening. The problem is as follows:

In an MPI parallelised code (For simplicity, lets says 2 processes P_0,P_1), suppose I have an array on each process A_0, A_1 which will always have a length matching the number of dofs and vertices on that process. (I will only be using CG1 elements so this should always hold true).

I also have two corresponding arrays which store the local vertex at which each value in A_k should lie on. Lets call these v_k.

Essentially, now I want to map the values of A_k onto a function u, but using v_k to ensure that they are mapped to the correct vertices.

If I am not mistaken, I can use vertex_to_dof_map to do this. However, each time I have tried the vertex_to_dof_map it gives me dof indices which are higher than that of the number of dofs owned by a process k.

Consider the MWI below:

from fenics import *
import numpy as np

# MPI communicator
comm = MPI.comm_world
rank = comm.Get_rank()

# Generate a mesh and functionspace
grid = RectangleMesh(Point(0,0),Point(10,10),10,10,"crossed")
V = FunctionSpace(grid,'CG',1)
u = Function(V)

# Connectivity
dofs = V.dofmap().dofs()
v2d = vertex_to_dof_map(V)

# Generate a random array on each process
A = np.random.uniform(0,1,len(dofs))

# Vertex mapping array
v = np.linspace(0,len(dofs)-1,len(dofs)).astype(int)
np.random.shuffle(v)

# Assign to the function vector
u.vector()[:] = A[v2d[v]]

The above code returns the errors (running with 2 processes):

IndexError: index 113 is out of bounds for axis 0 with size 109
IndexError: index 114 is out of bounds for axis 0 with size 112

If i investigate the ranges of the dofmap using:

print('Process:',rank,max(v2d),len(dofs),V.dofmap().ownership_range())

I get:

Process: 0 117 109 (0,109)
Process: 1 118 112 (109,221)

Clearly, some dof indices returned from vertex_to_dof_map exceed that of the number of dofs on any process k - Am i missing something regarding the mapping used here? As far I as know, they should all be process (local) indices here and should work fine.

If run with 1 process, it all works out fine.

Any help would be greatly appreciated,
Thanks

2

I would check the documentation to see if you’re handling the ghost entries correctly.

3

HI nate, thanks for your reply.

I had a look, and if I understand correctly - In a parallel setting, dofs which lie on the boundary between processes are marked as shared and are owned by both processes? Thus they both have a local index essentially mapping to the same global dof?

I cannot, however, find where the indexing scheme comes from - it appears to me that shared dofs have a different numbering system when they show up on a local process?

How would one go about assigning (in parallel) a value to a shared dof?

Please correct me if my understanding is incorrect,
Thanks

4

Since dolfin stores parts of the mesh on each process, dofs on process boundaries are owned by one process, and shared with the other(s). If you want to assign a value to only the local dofs, they are the first n dofs in V.dofmap().dofs().
The following code shows you how to get the number of dofs owned by each process, and how to get the global number of each degree of freedom not owned by the process that are in v2d:

from dolfin import *
import numpy as np

grid = RectangleMesh(Point(0,0),Point(10,10),1,1,"crossed")
V = FunctionSpace(grid,'CG',1)
u = Function(V)

# Connectivity
dofs = V.dofmap().dofs()
ownership_range = V.dofmap().ownership_range()
num_dofs_local = ownership_range[1] - ownership_range[0]
global_unowned = V.dofmap().local_to_global_unowned()
v2d = vertex_to_dof_map(V)

print(MPI.comm_world.rank, f"Num owned dofs {num_dofs_local}, Num dofs on process {len(v2d)}  Unowned dofs (global) {global_unowned}")

which gives

1 Num owned dofs 4, Num dofs on process 4  Unowned dofs (global) []
0 Num owned dofs 1, Num dofs on process 4  Unowned dofs (global) [4 3 1]

I would recommend only assigning data to dofs owned by the process

2 Likes
5

Hi @dokken,

Thanks so much for this insight - everything makes sense now, and I got my application working. Cheers!

6

@dokken If I’m looking at dof_to_vertex_map(V) how do I know which indices correspond to which dofs in V.dofmap().dofs()? Because there are fewer local dofs than there are indices in the map, it’s not clear to me how the map corresponds to the dofs. I want to initialize the values of a function with data that I have, but I have to map the correct data to the proper vertices. I’ve been using dof_to_vertex_map(V) to do this but I’m struggling to get it to work in parallel:

V = fenics.FunctionSpace(mesh, 'CG', 1)
loc0, loc1 = V.dofmap().ownership_range()
dofs = V.dofmap().dofs()
d_to_vert = dof_to_vertex_map(V)
# T_0 is an array of raster data. 
T_local = T_0.flatten()[loc0:loc1][d_to_vert]

d_to_vert is too large and has indices that are too large to be used as a mask here. How do I adapt the mask to only map the local dofs to their vertices?

Thanks!

7

dof_to_vertex_map includes ghosted dofs, as easily illustrated by:

from dolfin import *

n = 10
mesh = UnitSquareMesh(n, n)
V = FunctionSpace(mesh, "Lagrange", 1)
u = Function(V)


loc0, loc1 = V.dofmap().ownership_range()
d_to_vert = dof_to_vertex_map(V)
num_ghosts = len(V.dofmap().local_to_global_unowned())
print(
    f"Number of local dofs {loc1-loc0}, Number of dofs (including ghosts) {len(d_to_vert)}, Number of ghosts {num_ghosts}")

returning the following when executed on three processes

Number of local dofs 40, Number of dofs (including ghosts) 48, Number of ghosts 8
Number of local dofs 39, Number of dofs (including ghosts) 45, Number of ghosts 6
Number of local dofs 42, Number of dofs (including ghosts) 47, Number of ghosts 5

Thus you should only access the loc1-loc1 first entities in dof_to_vertex_map(V) if you only want to work with local entries.

2 Likes
8

@dokken How does that work exactly? If I access the first loc1-loc0 entries in dof_to_vertex_map(V) then I’ll still have indices that are out of bounds for something of the size of dofs. I guess I’m interpreting what you’re saying to be something like verts = dofs[d_to_vert[:loc1-loc0]], but that obviously doesn’t work with the indexing issue. I feel like I’m still missing something basic here.

9

You do not need the last part of your function,
i.e.

should be verts=d_to_vert[:loc1-loc0]
This gives you the local vertex indices on the current process that are associated with the dofs owned by this process.

1 Like
10

Thanks a ton! However, I still don’t understand how that vertex indexing works. If I want to initialize a vector with data, I would think I need to do something like this:

local_data = data[loc0:loc1]
local_vert_idxs = [d_to_vert][:loc1-loc0]
T_n.vector().set_local(local_data[local_vert_idxs])

But local_vert_idxs here has indices that are out of range for local_data.

Thanks for being patient with me; I’m really trying to get this.

11

The vertices you get out is using local numbers.
See the following code on how to map global data to the local dofs:

from dolfin import *
from mpi4py import MPI as _MPI
import numpy as np

n = 10
mesh = UnitSquareMesh(n, n)
V = FunctionSpace(mesh, "Lagrange", 1)
u = Function(V)

num_dofs = V.dim()
num_vertices = mesh.num_entities_global(0)
assert num_dofs == num_vertices
global_data = np.arange(num_dofs, dtype=np.int32)

loc0, loc1 = V.dofmap().ownership_range()
d_to_vert = dof_to_vertex_map(V)

global_vertex_numbers = mesh.topology().global_indices(0)
global_vertices = global_vertex_numbers[d_to_vert[:loc1-loc0]]
print(
    f"Process {MPI.comm_world.rank} Num vertices with owned dofs: {len(global_vertices)}")

print(f"Process {MPI.comm_world.rank} Total number of dofs {MPI.comm_world.allreduce(len(global_vertices), op=_MPI.SUM)}")

local_data = global_data[global_vertices]
assigned_data = MPI.comm_world.gather(local_data, root=0)

if MPI.comm_world.rank == 0:
    all_data = np.hstack(assigned_data)
    unique_data = np.unique(all_data)
    print(f"Unique vertices: {len(unique_data)}")
1 Like
12

@dokken Thanks! I’ve been trying to work with this, but I’m not sure how to recover the data with a vert_to_d mask or something. I’m initializing data as below:

loc0, loc1 = V.dofmap().ownership_range()
global_vertex_numbers = mesh.topology().global_indices(0)
global_vertices = global_vertex_numbers[d_to_vert[:loc1-loc0]]

T_0 = 2d_np_array_spatial_data
T_local = T_0.flatten()[global_vertices]
T_n.vector().set_local(T_local)
T_n.vector().apply('insert')
as_backend_type(T_n.vector()).vec().ghostUpdate()

Then in a for loop I step through time solving for new T values and at certain times I want to save the numpy array of T values. How can I convert back to the original spatial coordinates to save the new values? I’ve been trying to play around with the global_vertices you showed before to be able to something like

result = T_n.vector().gather_on_zero()
if mpi_rank == 0:
    vert_to_d = ? 
    perimeters[j] = result[vert_to_d]

but everything I’ve tried has failed and I don’t understand how I should go about this. Any help appreciated!

I would also add, I had a reason for using dolfin instead of dolfinx, but it no longer applies, so if this is easier to do/advisable I’ll happily upgrade to dolfinx and begin learning.

13

Here is an example on how to map from a function to the mesh geometry:

from dolfin import *
from mpi4py import MPI as _MPI
import numpy as np

n = 3
mesh = UnitSquareMesh(n, n)
V = FunctionSpace(mesh, "Lagrange", 1)
u = Function(V)
u.interpolate(Expression("2*x[0]", degree=1))

num_dofs = V.dim()
num_vertices = mesh.num_entities_global(0)
assert num_dofs == num_vertices
global_data = np.arange(num_dofs, dtype=np.int32)

loc0, loc1 = V.dofmap().ownership_range()
d_to_vert = dof_to_vertex_map(V)

global_vertex_numbers = mesh.topology().global_indices(0)


global_data = np.zeros(num_vertices, dtype=np.float64)
global_data[global_vertex_numbers[d_to_vert[:loc1-loc0]]
            ] = u.vector().get_local()[:loc1-loc0]

# print(V.tabulate_dof_coordinates()[:loc1-loc0])
global_data = MPI.comm_world.allreduce(global_data, op=_MPI.SUM)
global_vertices = np.zeros((num_vertices, mesh.geometry().dim()))
vertices = MPI.comm_world.gather(mesh.coordinates(), root=0)
vertex_order = MPI.comm_world.gather(mesh.topology().global_indices(0), root=0)

if MPI.comm_world.rank == 0:
    for vertex, order in zip(vertices, vertex_order):
        global_vertices[order] = vertex
    for coordinate, value in zip(global_vertices, global_data):
        assert(2*coordinate[0] == value)
2 Likes
14

That’s perfect! Thanks so much for all your help!