This is a read only copy of the old FEniCS QA forum. Please visit the new QA forum to ask questions

Setting Function values in parallel

+6 votes

Hi everyone,

I am trying to modify values of function in parallel using the following code:

from dolfin import *                                                            
from mpi4py import MPI as nMPI                                                  

comm = nMPI.COMM_WORLD                                                          

mesh = UnitSquareMesh(10, 10)                                                   
V = FunctionSpace(mesh, "CG", 1)                                                
u = Function(V)                                                                 
v = u.vector()                                                                  
dofmap = V.dofmap()                                                             

rank = comm.Get_rank()                                                          

for cell in cells(mesh):                                                      
  nodes = dofmap.cell_dofs(cell.index())                                      
  for node in nodes:                                                          
    v[node] = rank                                                          

plot(u, interactive=True)  

However, the program deadlocks if run with more than one process. What is the
correct way to this? Thanks for help.

asked Jul 12, 2013 by MiroK FEniCS Expert (80,920 points)

1 Answer

+5 votes
 
Best answer

If you just need to set all the local dofs to (process-wise constant) scalar then v[:] = rank executed by all the processes at the same time should work. But I guess that your snippet is so simple only for demonstration of your issue.

The design of numpy Vector interface in parallel is still open question - check issue 54. Essentially the problem is that Vector.__setitem__() method is collective with the present implementation because of Vector.apply("insert") call hidden therein.

In your case resolution could be to:

  • first collect indices which are to be set
  • prepare array values corresponding to indices
  • finally call v[indices] = values by every process

Problem is that indices and values must be global, i.e. same on all processes to get deterministic results. Alternatively you could prepare local indices and values and call

for i in range(MPI.num_processes()):
    if i == MPI.process_number():
        v[indices] = values
    else:
        v.apply("insert")

I.e. you call __setitem__ by one process while just setting nothing by every other process and repeat this for all the processes.

Preferably you should avoid using numpy interface (until issue 54 is resolved) and use Vector.set_local() method instead if possible. It is collective but accepts array with local values of all the local dofs.

answered Jul 12, 2013 by Jan Blechta FEniCS Expert (51,420 points)
selected Jul 12, 2013 by MiroK

Hi Jan,

thanks for the quick response. I used your suggestion and modified the code as follows:

ranks = []                                                                         
visited = []                                                                    

node_min, node_max = v.local_range()                                            

for cell in cells(mesh):                                                        
  nodes = dofmap.cell_dofs(cell.index())                                        
  for node in nodes:                                                            
    if (node not in visited) and node_min <= node <= node_max:                  
      visited.append(node)                                                      
      ranks.append(float(rank))                                                    

ranks = numpy.array(ranks)                                                            

v.set_local(ranks)

The program now works as expected. Thanks again.

Hi,
at present i have a similar trouble that can be solved using the above solution. But, if i run your proposed code in parallel on fenics 1.6 doesn't work.

Can you point me how to "update" your solution to the actual version of fenics?

Thanks in advance!

Ask a new question with MWE, please.

...