Pass vector value to function when using MPI

General
1

Dear all! I try to pass value of a vector to a function. However, when using MPI, there will be errors.

Here is my code:

from mpi4py import MPI
import ufl
from dolfinx import mesh, fem
import numpy as np
domain = mesh.create_unit_square(MPI.COMM_WORLD, 10, 10, mesh.CellType.triangle)
V = fem.functionspace(domain, ("Lagrange", 1))
a = np.ones(11*11)
l = fem.Function(V)
l.x.array[:] = a

When I use the command python test.py, everything is fine. However, when I use mpirun -n 2 python test.py , there will be errors like:

l.x.array[:] = a
ValueError: could not broadcast input array from shape (121,) into shape (75,)

Sometimes I will need to restore the vector a in a numpy array and copy them to another function, so I hope that someone can help me with this error.

2

When running with MPI, your mesh is distributed. This means that all processes do not contain all cells, which also means that all processes do not have all degrees of freedom.

As it seems like you are assuming a 1-1 relationship between the input vertices and your degrees of freedom, you should consider: How to apply a boundary condition to a range of DOFs? - #9 by dokken

3

Thanks for your code. But I am still confused about how to pass values to a functions. With the code I only get the local index, but I do not know how to pass values. Suppose I define a function space V and I want to pass values a to a function l in V, how to make it use MPI without errors? I just do not know how to change the function input_indices_to_dof such that it will return the local nodes and the corresponding global nodes.

4

input_indices_to_dofs takes a range of your original input indices, i.e. say you have data marked from 0 to 10 000, and maps it to the local indices on the current process.

I.e. if you take the aforementioned data into a system where the mesh has been partitioned in 4 pieces (0-2500, 2501-5000, 5001-7500, 7501-1000), it will give you the local indices for each of these on the given process.

As you get the input_global_indices of the current process through the mesh geometry, you can use these to get the values a[i] that you want to add to the process.

Furthermore, as your example above is simply passing 1 everywhere, this can be done with l.x.array[:] = 1.

5

Thanks for your reply. I just give a simple example for a, where sometimes a will be much more complicated. I try to use you code to pass values. It will work when I do not use mpirun. But there is still some errors when I use mpirun.

Here is a simple code:

from petsc4py import PETSc
from mpi4py import MPI
import dolfinx
import numpy as np


def input_indices_to_dof(V, global_dofs):
    mesh = V.mesh
    mesh.topology.create_connectivity(0, mesh.topology.dim)
    # Create a mapping from global to local node indices
    num_nodes_global = mesh.geometry.index_map().size_global
    local_to_input_global = mesh.geometry.input_global_indices
    global_to_local = np.full(num_nodes_global, -1, dtype=np.int32)
    global_to_local[local_to_input_global] = np.arange(len(local_to_input_global), dtype=np.int32)
    

    geom_layout =  mesh.geometry.cmap.create_dof_layout()
    V_layout = V.dofmap.dof_layout
    for i in range(mesh.topology.dim):
        num_entities_per_cell = dolfinx.cpp.mesh.cell_num_entities(
            mesh.topology.cell_type, i
    )
        for j in range(num_entities_per_cell):
            gl = geom_layout.entity_dofs(i, j)
            vl = V_layout.entity_dofs(i, j)
            assert vl == gl, "Mismatch in dof layout"
    geom_to_dof_layout = np.empty(V.dofmap.index_map.size_local+ V.dofmap.index_map.num_ghosts, dtype=np.int32)
    geom_to_dof_layout[mesh.geometry.dofmap.reshape(-1)] = V.dofmap.list.reshape(-1)

    # Remove all nodes that are not on process
    local_nodes = global_to_local[global_dofs]
    local_nodes = local_nodes[local_nodes >= 0]
    return geom_to_dof_layout[local_nodes]

if __name__ == "__main__":
    mesh = dolfinx.mesh.create_unit_square(MPI.COMM_WORLD, 10, 10)
    V = dolfinx.fem.functionspace(mesh, ("Lagrange", 1))
    u = dolfinx.fem.Function(V)
    num_nodes_global = mesh.geometry.index_map().size_global
    np.random.seed(0)
    global_indices = np.random.randint(0, num_nodes_global, 11*11)  
    values = np.random.rand(11*11)  
    u.x.array[:]=values[input_indices_to_dof(V, global_indices)]
6

Thanks for your advice but I still do not know how to do it. The problem is that I will restore the matrix a from a function u and I want to pass values a to another function l. There are two obstacles. One is that I do not know how to obtain the local index on the certain process and their corresponding global index. The other one is that once I know the local index and the corresponding global index, how to pass values? Do I only simply use the code l.x.array[ local index]=a[global index]? I really hope that you can help me with that.