Parallel domain decomposition for particle methods: Part 4

Introduction

The Plimpton scheme of communicating ghost information between patches was described in Part 3 of this series. Let us now see how a similar approach can be used to migrate particles that have moved across patches.

In the animation below we just move the particles within each patch randomly. To make the identity of the particles clear, we have used solid circles for the patch particles and three-quarter circle for the particles in the ghost regions. As you can see, some of the particles have moved outside the patches and need either to be deleted (if they have left the computational domain - assuming that the domain size remains unchanged) or they need to be moved to adjacent patches.

Plimpton’s scheme for migrating particles

If we run Plimpton’s scheme in reverse order, we can move the particles to the appropriate patches with only four communication steps in 2D and six in 3D. Notice in the animation below that we start with a search region in the x-direction that contains the top and bottom patches along with the right (or left) patch. We relocate particles in this region first and then need to move particles only in the top and bottom patches. Note also that the ghost particles have been moved back to their original locations, indicating that we can ignore these during the migration process. Depending on the requirements of the problem, we may either delete particles that have left the domain, introduce them back in a periodic manner, or extend the domain itself.

MPI implementation

The implementation of the migration process is similar to that for ghost exchange. A typical migrateParticles function can have the following form:

using ParticleIDHashMap = std::unordered_set<ParticleID>;
void migrateParticle(...., const Vec& patchWidth, ParticlePArray& patchParticles)
{
  ParticleIDHashMap sentParticles;  // sent particles per process
  ParticlePArray    recvParticles;  // received particles per process
  // First migrate in the x-direction
  d_patchP->sendRecvMigrateXMinus(boostWorld, patchWidth, patchParticles);
  d_patchP->sendRecvMigrateXPlus(boostWorld, patchWidth, patchParticles);
  d_patchP->waitToFinishX();
  d_patchP->combineSentParticlesX(sentParticles);
  d_patchP->combineReceivedParticlesX(recvParticles);
  d_patchP->deleteSentParticles(sentParticles, patchParticles);
  d_patchP->addReceivedParticles(recvParticles, patchParticles);
  // Next migrate in the y-direction
  sentParticles.clear();
  recvParticles.clear();
  d_patchP->sendRecvMigrateYMinus(boostWorld, patchWidth, patchParticles);
  d_patchP->sendRecvMigrateYPlus(boostWorld, patchWidth, patchParticles);
  d_patchP->waitToFinishY();
  d_patchP->combineSentParticlesY(sentParticles);
  d_patchP->combineReceivedParticlesY(recvParticles);
  d_patchP->deleteSentParticles(sentParticles, patchParticles);
  d_patchP->addReceivedParticles(recvParticles, patchParticles);
  // Next migrate in the z-direction
  sentParticles.clear();
  recvParticles.clear();
  d_patchP->sendRecvMigrateZMinus(boostWorld, patchWidth, patchParticles);
  d_patchP->sendRecvMigrateZPlus(boostWorld, patchWidth, patchParticles);
  d_patchP->waitToFinishZ();
  d_patchP->combineSentParticlesZ(sentParticles);
  d_patchP->combineReceivedParticlesZ(recvParticles);
  d_patchP->deleteSentParticles(sentParticles, patchParticles);
  d_patchP->addReceivedParticles(recvParticles, patchParticles);
  // delete outgoing particles (if needed by the problem)
  d_patchP->removeParticlesOutsidePatch(patchParticles);
}

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 described in Part 3 now has a few more methods. Let us see how some of these new functions may be implemented.

The first new function is sendRecvMigrateXMinus which is the equivalent of sendRecvGhostXMinus for th emigration process. Note that the only difference between these two function is the definition of the search box.

void 
Patch::sendRecvMigrateXMinus(boost::mpi::communicator& boostWorld, 
                             const Vec& neighborWidth,
                             const ParticlePArray& particles) 
{
  if (d_xMinus.d_boundary == PatchBoundary::inside) {
    Vec neighborLower = d_lower - neighborWidth;
    Vec neighborUpper = d_upper + neighborWidth;
    neighborUpper.setX(d_lower.x());
    Box neighborBox(neighborLower, neighborUpper);
    d_xMinus.asyncSendRecv(boostWorld, 
                           neighborBox, d_tolerance,
                           particles);
  }
}

The next new method is combineSentParticlesX which is defined as:

void 
Patch::combineSentParticlesX(ParticleIDHashMap& sent) 
{
  d_xMinus.combineSentParticles(sent);
  d_xPlus.combineSentParticles(sent);
}

The combineSentParticles method in PatchNeighborComm is defined as

void 
PatchNeighborComm::combineSentParticles(ParticleIDHashMap& sent) 
{
  if (!d_sentParticles.empty()) {
    for (const auto& particle : d_sentParticles) {
      sent.insert(particle->getId());
    }
  }
}

One also needs to delete the sent particles from the patch, using the method deleteSentParticles; this is where the use of a hashmap becomes handy.

void
Patch::deleteSentParticles(const ParticleIDHashMap& sent,
                           ParticlePArray& particles)
{
  if (sent.size() > 0) {
    particles.erase(
      std::remove_if(
        particles.begin(), particles.end(),
        [&sent](const ParticleP& particle) {
          return (sent.find(particle->getId()) != std::end(sent));
        }),
      std::end(particles));
  }
}

Finally, we add the received particles to the list of particles in the patch using addReceivedParticles:

void
Patch::addReceivedParticles(const ParticlePArray& received,
                            ParticlePArray& particles)
{
  particles.insert(particles.end(), received.begin(), received.end());
}

In some special cases, we will also need to remove particles outside the domain. One approach is to use removeParticlesOutsidePatch, but this step is typically not recommended in general as it is costly and often not necessary.

void 
Patch::removeParticlesOutsidePatch(ParticlePArray& particles)
{
  Box box(d_lower, d_upper);
  double epsilon = d_tolerance;
  particles.erase(
    std::remove_if(
      particles.begin(), particles.end(),
      [&box, &epsilon](ParticleP particle) {
        if (box.inside(particle->currentPos(), epsilon)) {
          return false;
        }
        return true;
      }),
    particles.end());
}

Remarks

In the next part of this series we will explore how information about forces can be communicated across patches.