Reply To: Particle periodicity running with MPI.

[Rookie]

Dear Jan,

I have successfully run it in the latest version as well. When you directly copy this code and compile it, you may encounter an error related to the particle period. In that case, replace the period code in this file with the one I sent earlier.

/olb-1.7r0/src/particles/subgrid3DLegacyFramework/boundaries/periodicBoundary3D.h

#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)
  // std::vector<T> _minPhys, _maxPhys, _extend;
  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)
{
  if (_x) {
    if (p->getPos()[0] > _maxPhys[0]) {
      p->getPos()[0] -= _extend[0];
      ++_jumper[0];
      int C = this->_cuboidGeometry.get_iC(p->getPos()[0], p->getPos()[1], p->getPos()[2], _overlap);
      p->setCuboid(C);
    }
    else if (p->getPos()[0] < _minPhys[0]) {
      p->getPos()[0] += _extend[0];
      ++_jumper[1];
      int C = this->_cuboidGeometry.get_iC(p->getPos()[0], p->getPos()[1], p->getPos()[2], _overlap);
      p->setCuboid(C);
    }
  }
  if (_y) {
    if (p->getPos()[1] > _maxPhys[1]) {
      p->getPos()[1] -= _extend[1];
      ++_jumper[2];
      int C = this->_cuboidGeometry.get_iC(p->getPos()[0], p->getPos()[1], p->getPos()[2], _overlap);
      p->setCuboid(C);
    }
    else if (p->getPos()[1] < _minPhys[1]) {
      p->getPos()[1] += _extend[1];
      ++_jumper[3];
      int C = this->_cuboidGeometry.get_iC(p->getPos()[0], p->getPos()[1], p->getPos()[2], _overlap);
      p->setCuboid(C);
    }
  }
  if (_z) {
    if (p->getPos()[2] > _maxPhys[2]) {
      p->getPos()[2] -= _extend[2];
      ++_jumper[4];
      int C = this->_cuboidGeometry.get_iC(p->getPos()[0], p->getPos()[1], p->getPos()[2], _overlap);
      p->setCuboid(C);
    }
    else if (p->getPos()[2] < _minPhys[2]) {
      p->getPos()[2] += _extend[2];
      ++_jumper[5];
      int C = this->_cuboidGeometry.get_iC(p->getPos()[0], p->getPos()[1], p->getPos()[2], _overlap);
      p->setCuboid(C);
    }
  }
}

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

}

#endif

Best regards,
Rookie