PnumpyLaplacian

Use case: applying the Laplacian operator to a field

As an example of use case relying on nearest neighbor communication, we apply a finite difference version of the Laplace operator to a 2D field. We will use one ghost.

# create and set the values of the distributed array field
zda = pnumpy.ghZeros( zz.shape, zz.dtype, numGhosts=1 )
zda[:] = numpy.sin(numpy.pi*xx) * numpy.cos(2*numpy.pi*yy)

Next, we apply the Laplacian operator inside the domain, where no communication is required

# compute the star Laplacian in the interior, this does not require
# any communication

laplaceZ = 4 * zda[:]

# local neighbor contributions, no communication
laplaceZ[1:  , :] -= zda[0:-1,:]
laplaceZ[0:-1, :] -= zda[1:  ,:]
laplaceZ[:, 1:  ] -= zda[:,0:-1]
laplaceZ[:, 0:-1] -= zda[:,1:  ]

Let's know tackle the edges of the sub-domains. We need to know what our neighbors are to the north, south, east, and west:

noProc = dc.getNeighborProc(rk, ( 1,  0), periodic = (False, True))
soProc = dc.getNeighborProc(rk, (-1,  0), periodic = (False, True))
eaProc = dc.getNeighborProc(rk, ( 0,  1), periodic = (False, True))
weProc = dc.getNeighborProc(rk, ( 0, -1), periodic = (False, True))

Note that we assume that the east and west wrap around, periodicity is set to True. We can now access the ghost data to complete our operator:

weZData = zda.get(weProc, winID=(0, +1))
eaZData = zda.get(eaProc, winID=(0, -1))
soZData = zda.get(soProc, winID=(+1, 0))
noZData = zda.get(noProc, winID=(-1, 0))

weSlc = zda.getEllipsis(winID=(0, -1))
eaSlc = zda.getEllipsis(winID=(0, +1))
soSlc = zda.getEllipsis(winID=(-1, 0))
noSlc = zda.getEllipsis(winID=(+1, 0))

if weProc is not None:
    laplaceZ[weSlc] -= weZData
if eaProc is not None:
    laplaceZ[eaSlc] -= eaZData
if soProc is not None:
    laplaceZ[soSlc] -= soZData
if noProc is not None:
    laplaceZ[noSlc] -= noZData

Finally, we need to clean up the windows

zda.free()

Some comments are in order. First we need to account for the possiblity that there are no neighbors to the north and south of the sub-array by checking for soProc is None and similarly for noProc. Second, it is good practice to write the slicing operators in such a way that they don't depend on the resolution of the local grid, this is why we're using the getEllipsis. Third, communication takes place only in the get method. The user is required to explicitly call this method, making it easier to get good perfromance. Fourth, communication is global, across all the processes. Passing the None value for the processor rank will bypass the communication call.