Particle periodicity running with MPI.

[Rookie]

Dear OpenLB developers,

/* This file is part of the OpenLB library
*
* Copyright (C) 2014 Thomas Henn, Mathias J. Krause
* E-mail contact: info@openlb.net
* The most recent release of OpenLB can be downloaded at
* <http://www.openlb.net/&gt;
*
* This program is free software; you can redistribute it and/or
* modify it under the terms of the GNU General Public License
* as published by the Free Software Foundation; either version 2
* of the License, or (at your option) any later version.
*
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public
* License along with this program; if not, write to the Free
* Software Foundation, Inc., 51 Franklin Street, Fifth Floor,
* Boston, MA 02110-1301, USA.
*/

#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

This is the code for particle periodic boundaries that I want to use. This code can only be used in serial mode, and when run in parallel, it encounters issues due to particles crossing block boundaries defined by MPI core count. To illustrate the problem, I have set up two particles, A and B, starting at positions A1 and B1 respectively. When they move to positions A2 and B2, the simulation terminates because particle A cannot enter the next block. I understand that this is because the code lacks relevant MPI settings. I have also found your documentation on the new settings for particle periodic boundaries (olb-1.6r0/src/particles/communication/relocation.h), which explains particle serialization and MPI non-blocking communication settings. The crux of the issue is that my code lacks MPI communication-related settings. I hope you can provide me with some guidance on how to modify it. I have tried adding MPI communication settings myself, but without success. I am very grateful for your assistance.

https://postimg.cc/F7vmCHDM

Best regards,
Rookie

[jan]

Dear Rookie,

are you using the subgrid3DLegacyFramework in your setup? Unfortunately, I cannot provide support in that context, so please correct me if I’m wrong.

Here a few thoughts:

I don’t see why the explained error should be caused by missing communications (I don’t think you have to add any communication, without the periodic boundary the particles move from one to the other block without problem, don’t they?). Because, if I understand it correctly, then the particles are duplicated when they pass through a block boundary in the middle of the domain. The periodic boundary shouldn’t be called there, should it?

I could imagine that the hardcoded //cube extents with origin (0,0,0) doesn’t fit in your setup. Apparently, this part of the legacy code assumes that the simulation domain’s origin is at (0,0,0). Perhaps you should try to change your simulation setup to fit that assumption, if this doesn’t align with your setup.

In general, consider to check these values (inside the applyBoundary method) for correctness, that is, are the values would they should be?

_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];

However, I’m aware that this doesn’t explain why it’s working in sequential mode, but not in parallel. Are you sure that it worked in sequential? Do the above values differ between sequential and parallel mode? I believe that here the positions are merely updated and that the communication should take place somewhere else.

You could also try to change the places where the particle position is updated (in case some change to CuboidGeometry broke the logic, for example:

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);

Consider following movePositionToEnd and movePositionToStart of the current code: https://www.openlb.net/DoxyGen/html/d9/d46/namespaceolb_1_1particles_1_1communication.html#a1647f5b5b28d226caa8ae64e54f0320b
(min and max area already read from the SuperGeometry at the top and position refers to the current position)

Best regards,
Jan

• This reply was modified 4 days, 9 hours ago by jan.

[jan]

~You could also try to change the places where the particle position is updated~
You could also try to change the code where the particle position is updated

[Rookie]

Dear Jan,

Thank you for answering my question. I am indeed using the subgrid3DLegacyFramework in my code. Although these codes are not easy to use, they do contain the functionality I need. Why aren’t you using this part of the code anymore? If this makes you uncomfortable, I won’t bother you with this issue again. I want to thank you for telling me not to add additional communication. From the results of inserting two particles, it seems that particles cannot move from one block to another. The problem may be that when particles pass through the boundary of the domain, they are not copied. Particle periodic boundaries should indeed not be called there. The question about whether the origin (0,0,0) for fluid and particle periodicity is set consistently. I set the origin of the fluid simulation domain to (0,0,0), but because of the fluid periodicity setting, an extra layer is added outside the geometric region, so it becomes (-Δx, -Δy, -Δz). I output the extend and minPhys of the particle periodicity, and these values ​​do not change in both serial and parallel modes. I only set periodicity in the x and y directions, so the extend values ​​for these two directions are correct. I want to emphasize that periodicity can run successfully in serial mode, so could this be the reason for termination in parallel mode?

[SuperGeometryStatistics3D] materialNumber=1; count=39520; minPhysR=(0,0,0.00266667); maxPhysR=(0.250667,0.0826667,0.0346667)
[SuperGeometryStatistics3D] materialNumber=2; count=6080; minPhysR=(0,0,0); maxPhysR=(0.250667,0.0826667,0.0373333)
[prepareGeometry] Dimension of the channel in meters: x = 0.250667 ; y = 0.0826667 ; z = 0.0373333
[prepareGeometry] Prepare Geometry … OK
[main] noOfCuboid=8
[prepareLattice] Prepare Lattice …
[setInitialConditions] Set initial conditions …
[setInitialConditions] Set initial conditions … OK
[prepareLattice] Prepare Lattice … OK
[main] Prepare Particles …
_extend[0]=0.250667
_extend[1]=0.0826667
_extend[2]=0.032
(_minPhys[0],_minPhys[1],_minPhys[2])=(0,0,0.00266667)

Best regards,
Rookie

[mathias]

You should change to the new code. It is much better. It should not be a big deal to transfer the functionallites you need from the old one to the new code. Once, you have done that, please consider contributing it to the OpenLB community by adding it at https://gitlab.com/openlb/release .

[Author]

Posts