Reply To: Particle periodicity running with MPI.
OpenLB – Open Source Lattice Boltzmann Code › Forums › on OpenLB › General Topics › Particle periodicity running with MPI. › Reply To: Particle periodicity running with MPI.
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