How to parallelize code alongside assign()?

Errors
1

I’m trying to manipulate functions using assign() inside of a for-loop, but I want my code to run in parallel, especially when using solve() elsewhere inside of a time loop.

As a simple example, consider the following code code as run using mpirun -np 2

# Define submesh mesh and function spaces
mesh = BoxMesh(Point(0, 0), Point(box_length, box_length, 0.5),
                     64, 64, 32)

P1 = FiniteElement('CG', mesh.ufl_cell(), 1)
W = FunctionSpace(mesh, P1)


#Initialize functions
Aexp = Expression('2', degree=1)
A = interpolate(Aexp, W)

Bexp = Expression('2', degree=1)
B = interpolate(Bexp, W)

Cexp = Expression('0', degree=1)
C = interpolate(C, W)

#time loop
for t in range(0,5):
  C.assign(A + B)
  B.assign(C)  
  print(C.vector().get_local())

This returns the following error message

-------------------------------------------------------------------------
*** Error:   Unable to successfully call PETSc function 'VecCopy'.
*** Reason:  PETSc error code is: 73 (Object is in wrong state).
*** Where:   This error was encountered inside /home/conda/feedstock_root/build_artifacts/fenics-pkgs_1668433332839/work/dolfin/dolfin/la/PETScVector.cpp.
*** Process: 0
*** 
*** DOLFIN version: 2019.1.0
*** Git changeset:  3ea2183fbfe7277de9f16cbe1a9ffaab133ba1fa
*** -------------------------------------------------------------------------

This question seems similar to this question, but that question was never answered and my example is both slightly different in form and returns a different error.

Any help at all would be appreciated, even just an explanation of what is going on with the dofs in parallel, or pointing to the relevant documentation. Clearly, I know that assign itself does not probably need parallelization. Parallelizing other code that might be in the for-loop alongside assign() is my real goal. For instance, when assign() is used in a time loop to update a source field for some variational form to be solved.

2

Surely this is a typo?

3

Even with the fix I mentioned above, I cannot reproduce your issue using the lastest version of DOLFIN (on docker), or the 2019.1.0 dolfin quay image:

from dolfin import *
box_length = 1
# Define submesh mesh and function spaces
mesh = BoxMesh(Point(0, 0), Point(box_length, box_length, 0.5),
                     10,10,10)

P1 = FiniteElement('CG', mesh.ufl_cell(), 1)
W = FunctionSpace(mesh, P1)


#Initialize functions
Aexp = Expression('2', degree=1)
A = interpolate(Aexp, W)

Bexp = Expression('2', degree=1)
B = interpolate(Bexp, W)

Cexp = Expression('0', degree=1)
C = interpolate(Cexp, W)

#time loop
for t in range(0,5):
  C.assign(A + B)
  B.assign(C)  
  print(C.vector().get_local())
1 Like
4

Thanks for your help @dokken! You are correct about the typo.

I’m not sure why you cannot reproduce the error. I’m running the code from the conda-forge version running on a cluster. Could clashes between the conda installation and the cluster’s handling of multithreading or something be going on?

I was able to “fix” things by doing a hacky workaround. From this question I learned that c.assign(a) has different behavior than assign(c, a) when running under mpirun. Thus, the following code has the desired behavior for me:

# Define submesh mesh and function spaces
mesh = BoxMesh(Point(0, 0), Point(1, 1, 0.5),
                     64, 64, 32)

P1 = FiniteElement('CG', mesh.ufl_cell(), 1)
W = FunctionSpace(mesh, P1)

#Define pointwise addition of functions. Returns a numpy array of the pointwise sym of .get_local() values
def Add(functions):
    garrays = []
    for g in functions:
        garrays.append(g.vector().get_local())
    Sum = garrays[0] #initialize the sum variable with the first array in garrays
    
    #successively add each array in garrays together (-1 in loop +1 in index because sum begins with garrays[0])
    for i in range(len(garrays)-1): 
        Sum = Sum+garrays[i + 1]
    result = Function(W)
    result.vector().set_local(Sum)
    return result


Aexp = Expression('2', degree=1)
A = interpolate(Aexp, W)

Bexp = Expression('2', degree=1)
B = interpolate(Bexp, W)


Cexp = Expression('0', degree=1)
C = interpolate(Cexp, W)

for x in range(0,5):
  assign(C, Add([A,B]))
  assign(B, C)  
  print(C.vector().get_local())

While this workaround might “work” for now, I certainly don’t understand what’s going on? Bafflingly, this is the output from the print function using -np 2:

[0. 4. 0. ... 4. 4. 0.]
[0. 6. 0. ... 6. 6. 0.]
[0. 8. 0. ... 8. 8. 0.]
[ 0. 10.  0. ... 10. 10.  0.]
[ 0. 12.  0. ... 12. 12.  0.]
[4. 4. 4. ... 4. 4. 4.]
[6. 6. 6. ... 6. 6. 6.]
[8. 8. 8. ... 8. 8. 8.]
[10. 10. 10. ... 10. 10. 10.]
[12. 12. 12. ... 12. 12. 12.]

One of the processes seems to get it right, while the other doesn’t?

5

You error is in PETSc, which might mean there is an issue with the PETSc installation shipped with conda.

1 Like