Skip to content

Reply To: Particle periodicity running with MPI.

#8891
Rookie
Participant

Dear Jan,

I am very glad that I have fixed the code for the periodic boundary. However, I am not sure if I fixed it correctly. Initially, I thought it was an issue with MPI non-blocking communication, so I searched for information and tried using MPI_Barrier(MPI_COMM_WORLD). I found that the program could run correctly when the number of particles did not exceed the simulation domain, but it would deadlock when particles crossed the boundary. Therefore, I tried commenting out p->setCuboid(newCuboid) because particle allocation was no longer necessary. This confirmed that I correctly implemented the periodic boundary for particles; it prints that particles move from the current rank to the opposite side of the new cuboid. Could you help me confirm if I made the correct modifications? You can quickly complete the simulation by slightly reducing the simulation time and the number of particles.

#ifndef PERIODICBOUNDARY3D_H_
#define PERIODICBOUNDARY3D_H_

#include <math.h>
#include <vector>

namespace olb {

template<typename T, template<typename U> class PARTICLETYPE>
class ParticleSystem3D;

/*
 * Particle boundary based on a cube around the area with material number 1.
 * Only applicable to rectangles since if a particle leaves the area with
 * material number 1 it is moved to the opposing side of the area by
 * newPosition = oldPosition +/- extend(MaterialNumber=1).
 **/

template<typename T, template<typename U> class PARTICLETYPE>
class PeriodicBoundary3D : public Boundary3D<T, PARTICLETYPE> {
public:
  PeriodicBoundary3D(SuperGeometry<T,3>& sg, bool x,
                     bool y, bool z);
  PeriodicBoundary3D(PeriodicBoundary3D<T, PARTICLETYPE>& f);
  virtual ~PeriodicBoundary3D() { };
  virtual void applyBoundary(typename std::deque<PARTICLETYPE<T> >::iterator& p, ParticleSystem3D<T, PARTICLETYPE>& psSys);
  /// Returns number of particles that moved through the periodic boundary
  /// Order: x+, x-, y+, y-, z+, z-
  unsigned int* getJumper();
private:
  //cube extents with origin (0,0,0)
  olb::Vector<T, 3> _minPhys, _maxPhys, _extend;
  bool _x, _y, _z;
  unsigned int _jumper[6];
  CuboidGeometry3D<T>& _cuboidGeometry;
  T _overlap;
};

template<typename T, template<typename U> class PARTICLETYPE>
PeriodicBoundary3D<T, PARTICLETYPE>::PeriodicBoundary3D(
  SuperGeometry<T,3>& sg, bool x, bool y, bool z) : Boundary3D<T, PARTICLETYPE>(),
  _minPhys(sg.getStatistics().getMinPhysR(1)),
  _maxPhys(sg.getStatistics().getMaxPhysR(1)),
  _extend(_maxPhys - _minPhys),
  _x(x),
  _y(y),
  _z(z), _cuboidGeometry(sg.getCuboidGeometry())
{
  _minPhys = sg.getStatistics().getMinPhysR(1);
  _maxPhys = sg.getStatistics().getMaxPhysR(1);
  _extend[0] = _maxPhys[0] - _minPhys[0];
  _extend[1] = _maxPhys[1] - _minPhys[1];
  _extend[2] = _maxPhys[2] - _minPhys[2];
  for (int i=0; i<6; ++i) {
    _jumper[i] = 0;
  }
  _overlap = sg.getOverlap();
}

template<typename T, template<typename U> class PARTICLETYPE>
void PeriodicBoundary3D<T, PARTICLETYPE>::applyBoundary(
  typename std::deque<PARTICLETYPE<T> >::iterator& p,
  ParticleSystem3D<T, PARTICLETYPE>& psSys)
{
  bool crossed_boundary = false;
  
  if (_x) {
    if (p->getPos()[0] > _maxPhys[0]) {
      p->getPos()[0] -= _extend[0];
      ++_jumper[0];
      crossed_boundary = true;
      std::cout << "Particle crossed x+ boundary." << std::endl;
    }
    else if (p->getPos()[0] < _minPhys[0]) {
      p->getPos()[0] += _extend[0];
      ++_jumper[1];
      crossed_boundary = true;
      std::cout << "Particle crossed x- boundary." << std::endl;
    }
  }
  if (_y) {
    if (p->getPos()[1] > _maxPhys[1]) {
      p->getPos()[1] -= _extend[1];
      ++_jumper[2];
      crossed_boundary = true;
      std::cout << "Particle crossed y+ boundary." << std::endl;
    }
    else if (p->getPos()[1] < _minPhys[1]) {
      p->getPos()[1] += _extend[1];
      ++_jumper[3];
      crossed_boundary = true;
      std::cout << "Particle crossed y- boundary." << std::endl;
    }
  }
  if (_z) {
    if (p->getPos()[2] > _maxPhys[2]) {
      p->getPos()[2] -= _extend[2];
      ++_jumper[4];
      crossed_boundary = true;
      std::cout << "Particle crossed z+ boundary." << std::endl;
    }
    else if (p->getPos()[2] < _minPhys[2]) {
      p->getPos()[2] += _extend[2];
      ++_jumper[5];
      crossed_boundary = true;
      std::cout << "Particle crossed z- boundary." << std::endl;
    }
  }

  if (crossed_boundary) {
    int rank;
    MPI_Comm_rank(MPI_COMM_WORLD, &rank);
    int newCuboid = this->_cuboidGeometry.get_iC(p->getPos()[0], p->getPos()[1], p->getPos()[2], _overlap);
    // p->setCuboid(newCuboid);

    std::cout << "Rank " << rank << ": Particle crossed boundary to new cuboid " << newCuboid << std::endl;

    // Ensure synchronization among all processes.
    // MPI_Barrier(MPI_COMM_WORLD);
  }
}

template<typename T, template<typename U> class PARTICLETYPE>
unsigned int* PeriodicBoundary3D<T, PARTICLETYPE>::getJumper()
{
  return _jumper;
}

} // namespace olb

#endif
const T physConvergeTime = 0.1 * charPhysT;                   // time until until statistics sampling in seconds
const T physStatisticsTime = 0.1 * charPhysT;                  // statistics sampling time in seconds
const T particleMaxPhysT = 0.1 * charPhysT;
const T singleMaxPhysT = physConvergeTime + physStatisticsTime;
const T fluidMaxPhysT = physConvergeTime + physStatisticsTime + particleMaxPhysT; // max. simulation time in seconds
const T statisticsSave = 1. / 25.; // time between statistics samples in seconds
const int noOfParticles = 100;
const T checkstatistics = (T)fluidMaxPhysT / 200.;

Best regards,
Rookie

  • This reply was modified 1 week, 3 days ago by Rookie.
  • This reply was modified 1 week, 3 days ago by Rookie.