underworldcode / underworld2

underworld2: A parallel, particle-in-cell, finite element code for Geodynamics.
http://www.underworldcode.org/
Other
161 stars 58 forks source link

Gathering All Mesh Data in Parallel #660

Open gp37 opened 1 year ago

gp37 commented 1 year ago

Hi, I am new to parallel computing and I am running into a problem with gathering data on the root processor. I need to gather all the data together as I need to find a certain value in each column in my model. I dont get any errors with the following code but it does not give the right answer. For example, in parallel when I run the code below I get np.sum(recvbufT) = 118502.35 however I know the correct answer is 103315.66. Any help would be appreciated!

sendbufy = np.array(mesh.data[:,1]) sendbufx = np.array(mesh.data[:,0]) sendbufT = np.array(temperatureField.data)

size of data

sendcountsy = np.array(MPI.COMM_WORLD.gather(len(temperatureField.data), 0)) sendcountsx = np.array(MPI.COMM_WORLD.gather(len(temperatureField.data), 0)) sendcountsT = np.array(MPI.COMM_WORLD.gather(len(temperatureField.data), 0))

if rank == 0:

prepare to receive all this data

recvbufy = np.empty(sum(sendcountsy), dtype=float) recvbufx = np.empty(sum(sendcountsx), dtype=float) recvbufT = np.empty(sum(sendcountsT), dtype=float)

else: recvbuf = None recvbufy = None recvbufx = None recvbufT = None

Gather up all the data and put it in recvbuf

MPI.COMM_WORLD.Gatherv(sendbuf=sendbufy, recvbuf=(recvbufy, sendcountsy), root=0) MPI.COMM_WORLD.Gatherv(sendbuf=sendbufx, recvbuf=(recvbufx, sendcountsx), root=0) MPI.COMM_WORLD.Gatherv(sendbuf=sendbufT, recvbuf=(recvbufT, sendcountsT), root=0)

if rank == 0: print(np.sum(recvbufT))

bknight1 commented 1 year ago

Hi @gp37,

UW has a built in function to reduce values to root. You can try the following:

import underworld.function as fn
### reduce x, y and temperatureField to  root proc
x = fn.input()[0].evaluate_global(mesh.data)
y = fn.input()[1].evaluate_global(mesh.data)
T = temperatureField.evaluate_global(mesh.data)

### print the sum of the temperatureField
if rank == 0:
    print(np.sum(T))
gp37 commented 1 year ago

Thanks Ben,

I tried to implement the code above and I don't get any errors, however, the code seems to take a long time to run or is getting stuck in a loop or something. For example, I ran the above code for 5 minutes and still, I was not able to finish evaluating for x. This is true on multiple processors and just running in serial. Is this normal for the evaluate_global function?

gp37 commented 1 year ago

Giving an update, I actually get this error when running x = fn.input()[0].evaluate_global(mesh.data)

File "/opt/venv/lib/python3.10/site-packages/underworld/function/_function.py", line 744, in evaluate_global total_output[incoming_positions] = incoming_data IndexError: index 6633 is out of bounds for axis 0 with size 6633 Abort(1) on node 0 (rank 0 in comm 496): application called MPI_Abort(comm=0x84000002, 1) - process 0

and I am not sure how to fix that.

bknight1 commented 1 year ago

Hi @gp37, can you share the entire block of code? It looks like you're trying to update an array in certain locations, but the incoming data should contain the entire array (if using the code above).

The evaluate_global does take a while unfortunately, we could use MPI routines but the evaluate_global is much more straight forward to use

gp37 commented 1 year ago

Hi, Attached is a script with a simple mantle convection model in which I cannot get evaluate_global to work.

MantleConvectionParallel.ipynb.zip

bknight1 commented 1 year ago

Hi @gp37, I've had a look and come up with the following using MPI:

### on local nodes
x = np.ascontiguousarray(mesh.data[:,0])
y = np.ascontiguousarray(mesh.data[:,1])
T = np.ascontiguousarray(temperatureField.data[:,0])

### Collect local array sizes using the high-level mpi4py gather
sendcounts = np.array(comm.gather(len(x), root=0))

if rank == 0:
### creates dummy data on all nodes to store the data that is gathered
    x_data = np.zeros((sum(sendcounts)), dtype='float64')
    y_data = np.zeros((sum(sendcounts)), dtype='float64')
    T_data = np.zeros((sum(sendcounts)), dtype='float64')
else:
    x_data = None
    y_data = None
    T_data = None

comm.barrier()

### gather data on root proc (0)
comm.Gatherv(sendbuf=x, recvbuf=(x_data, sendcounts), root=0)
comm.Gatherv(sendbuf=y, recvbuf=(y_data, sendcounts), root=0)
comm.Gatherv(sendbuf=T, recvbuf=(T_data, sendcounts), root=0)

if rank == 0:
    '''do some stuff'''
    print(x_data.shape)

Let me know how it goes