Parallel domain decomposition for particle methods: Part 3
Introduction
In Part 2 of this series we showed the direct way of communicating ghost particles between patches. That approach requires 26 communication steps per patch in three-dimensions. In this article we discuss the approach suggested by Steve Plimpton (“Fast parallel algorithms for short-range molecular dynamics”, Sandia Report SAND91-1144.UC-405, 1993).
Plimpton’s paper has been cited almost 15,000 times since its publication. Among other things, the paper explains how the number of communication steps can be reduced to six in three dimensions.
Plimpton’s scheme for exchanging particles
If you run the animation below, you will notice that the left-right ghost particles are exchanged first. The ghost regions in the up-down directions are then extended to include the left-right ghost regions. The particles in these enlarged regions are then transferred in the up-down directions. Therefore, there are only four communication steps in two-dimensions. The same process can be used in three-dimensions, leading to only six communication steps.
MPI implementation
In Part 2 we defined
a PatchNeighborComm
struct and a Patch
struct. We can keep the PatchNeighborComm
struct in the same form, with the possible addition of a method of two. However, the Patch
struct becomes considerably simplified as show below.
Patch struct
The Patch
struct now needs only six PatchNeighborComm
objects but three waitToFinish
and combineReceivedParticles
methods.
struct Patch {
int d_rank;
double d_ghostWidth, d_tolerance;
IntVec d_patchMPICoords;
Vec d_lower, d_upper;
PatchNeighborComm d_xMinus, d_yMinus, d_zMinus, d_xPlus, d_yPlus, d_zPlus;
Patch(MPI_Comm& cartComm,
int rank, const IntVec& mpiCoords, const Vec& lower, const Vec& upper,
double ghostWidth, double tolerance);
void setXMinus(MPI_Comm& cartComm); // A patch boundary is at the domain boundary
void setXPlus(MPI_Comm& cartComm);
void setYMinus(MPI_Comm& cartComm);
// .....
void setZPlus(MPI_Comm& cartComm);
// Step 1: Send and receive data from x+ and x-
void sendRecvGhostXMinus(boost::mpi::communicator& boostWorld,
const ParticlePArray& patchParticles);
void sendRecvGhostXPlus(boost::mpi::communicator& boostWorld,
const ParticlePArray& patchParticles);
void waitToFinishX();
void combineReceivedParticlesX(ParticlePArray& patchParticles);
// Step 2: Send and receive data from y+ and y-
void sendRecvGhostYMinus(boost::mpi::communicator& boostWorld,
const ParticlePArray& patchParticles);
void sendRecvGhostYPlus(boost::mpi::communicator& boostWorld,
const ParticlePArray& patchParticles);
void waitToFinishY();
void combineReceivedParticlesY(ParticlePArray& patchParticles);
// Step 3: Send and receive data from z+ and z-
void sendRecvGhostZMinus(boost::mpi::communicator& boostWorld,
const ParticlePArray& patchParticles);
void sendRecvGhostZPlus(boost::mpi::communicator& boostWorld,
const ParticlePArray& patchParticles);
void waitToFinishZ();
void combineReceivedParticlesZ(ParticlePArray& patchParticles);
};
The new implementation of the Patch
struct is shown below. The initialization of
the struct is the same as before; as are the setXMinus
etc. methods. Also, the
sendRecvGhostXMinus
and sendRecvGhostXPlus
methods are the same as before.
Patch::Patch(//...) {
//..... Same as for direct communication
}
void Patch::setXMinus(MPI_Comm& cartComm) {
//..... Same as for direct communication
}
//......
Next we do the communication of the x+ and x- neighbors and combine the received particles with the particles in the current patch.
// Step 1: Do the x+ and x- communication
void Patch::sendRecvGhostXMinus(//....) {
//..... Same as for direct communication
}
void Patch::sendRecvGhostXPlus(//....) {
//..... Same as for direct communication
}
void Patch::waitToFinishX() {
if (d_xMinus.d_boundary == PatchBoundary::inside) {
d_xMinus.waitToFinish();
}
if (d_xPlus.d_boundary == PatchBoundary::inside) {
d_xPlus.waitToFinish();
}
}
void Patch::combineReceivedParticlesX(ParticlePArray& patchParticles) {
received.clear();
d_xMinus.combineReceivedParticles(patchParticles);
d_xPlus.combineReceivedParticles(patchParticles);
}
Now that the patchParticles
vector has been updated, we can repeat the
process for the y+ and y- directions. Note that the size of the ghost
region has been expanded in the negative and positive x-direction.
// Step 2: Do the y+ and y- communication
void Patch::sendRecvGhostYMinus(boost::mpi::communicator& boostWorld, const ParticlePArray& patchParticles)
{
if (d_yMinus.d_boundary == PatchBoundary::inside) {
Vec ghostLower = d_lower - Vec(d_ghostWidth, 0.0, 0.0);
Vec ghostUpper = d_upper + Vec(d_ghostWidth, 0.0, 0.0);
ghostUpper.setY(d_lower.y() + d_ghostWidth);
Box ghostBox(ghostLower, ghostUpper);
d_yMinus.asyncSendRecv(boostWorld, ghostBox, d_tolerance, patchParticles);
}
}
void Patch::sendRecvGhostYPlus(boost::mpi::communicator& boostWorld, const ParticlePArray& patchParticles)
{
if (d_yPlus.d_boundary == PatchBoundary::inside) {
Vec ghostLower = d_lower - Vec(d_ghostWidth, 0.0, 0.0);
ghostLower.setY(d_upper.y() - d_ghostWidth);
Vec ghostUpper = d_upper + Vec(d_ghostWidth, 0.0, 0.0);
Box ghostBox(ghostLower, ghostUpper);
d_yPlus.asyncSendRecv(boostWorld, ghostBox, d_tolerance, patchParticles);
}
}
void Patch::waitToFinishY() {
if (d_yMinus.d_boundary == PatchBoundary::inside) {
d_yMinus.waitToFinish();
}
if (d_yPlus.d_boundary == PatchBoundary::inside) {
d_yPlus.waitToFinish();
}
}
void Patch::combineReceivedParticlesY(ParticlePArray& patchParticles) {
received.clear();
d_yMinus.combineReceivedParticles(patchParticles);
d_yPlus.combineReceivedParticles(patchParticles);
}
Finally, we do the third stage of communication in the z-direction. The ghost-regions have now been expanded to contain for the x-, x+ and y-, y+ extensions.
// Step 2: Do the z+ and z- communication
void Patch::sendRecvGhostZMinus(boost::mpi::communicator& boostWorld, const ParticlePArray& patchParticles)
{
if (d_zMinus.d_boundary == PatchBoundary::inside) {
Vec ghostLower = d_lower - Vec(d_ghostWidth, d_ghostWidth, 0.0);
Vec ghostUpper = d_upper + Vec(d_ghostWidth, d_ghostWidth, 0.0);
ghostUpper.setZ(d_lower.z() + d_ghostWidth);
Box ghostBox(ghostLower, ghostUpper);
d_zMinus.asyncSendRecv(boostWorld, ghostBox, d_tolerance, patchParticles);
}
}
void Patch::sendRecvGhostZPlus(boost::mpi::communicator& boostWorld, const ParticlePArray& patchParticles)
{
if (d_zPlus.d_boundary == PatchBoundary::inside) {
Vec ghostLower = d_lower - Vec(d_ghostWidth, d_ghostWidth, 0.0);
ghostLower.setZ(d_upper.z() - d_ghostWidth);
Vec ghostUpper = d_upper + Vec(d_ghostWidth, d_ghostWidth, 0.0);
Box ghostBox(ghostLower, ghostUpper);
d_zPlus.asyncSendRecv(boostWorld, ghostBox, d_tolerance, patchParticles);
}
}
void Patch::waitToFinish() {
if (d_zMinus.d_boundary == PatchBoundary::inside) {
d_zMinus.waitToFinish();
}
if (d_zPlus.d_boundary == PatchBoundary::inside) {
d_zPlus.waitToFinish();
}
}
void Patch::combineReceivedParticlesY(ParticlePArray& patchParticles) {
received.clear();
d_zMinus.combineReceivedParticles(patchParticles);
d_zPlus.combineReceivedParticles(patchParticles);
}
The particle exchange function
The Plimpton particle exchange function the main simulation can then be simplified to the following.
void ParticleCode::exchangeGhostParticles() {
// Initialize list of all particles in patch + ghost
mergeParticleVec.clear();
mergeParticleVec = particleVec;
// First the x+/x- directions
d_patchP->sendRecvGhostXMinus(boostWorld, mergeParticleVec);
d_patchP->sendRecvGhostXPlus(boostWorld, mergeParticleVec);
d_patchP->waitToFinishX();
d_patchP->combineReceivedParticlesX(mergeParticleVec);
// Next the y+/y- directions
d_patchP->sendRecvGhostYMinus(boostWorld, mergeParticleVec);
d_patchP->sendRecvGhostYPlus(boostWorld, mergeParticleVec);
d_patchP->waitToFinishX();
d_patchP->combineReceivedParticlesX(mergeParticleVec);
// Next the z+/z- directions
d_patchP->sendRecvGhostZMinus(boostWorld, mergeParticleVec);
d_patchP->sendRecvGhostZPlus(boostWorld, mergeParticleVec);
d_patchP->waitToFinishX();
d_patchP->combineReceivedParticlesX(mergeParticleVec);
}
In this case the number of communication steps is much smaller. However, there is a wait period at the end of each communication step that may reduce the benefits of the Plimpton approach in some situations where the particles are unevenly distributed.
Remarks
Plimpton’s scheme is attractive for its simplicity in communicating ghost particle positions. However, there are two more important communication steps that need to be considered - the computation of interparticle forces and the migration of particles between patches. In the next part of this series, we will discuss the migration of particles when we use the Plimpton scheme.