diff options
Diffstat (limited to 'src/particles/superParticleSystem3D.hh')
-rw-r--r-- | src/particles/superParticleSystem3D.hh | 2257 |
1 files changed, 2257 insertions, 0 deletions
diff --git a/src/particles/superParticleSystem3D.hh b/src/particles/superParticleSystem3D.hh new file mode 100644 index 0000000..f6450f7 --- /dev/null +++ b/src/particles/superParticleSystem3D.hh @@ -0,0 +1,2257 @@ +/* This file is part of the OpenLB library + * + * Copyright (C) 2016 Thomas Henn + * E-mail contact: info@openlb.net + * The most recent release of OpenLB can be downloaded at + * <http://www.openlb.net/> + * + * 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 SUPERPARTICLESYSTEM_3D_HH +#define SUPERPARTICLESYSTEM_3D_HH + +#define shadows + +#include <utility> +#include <string> +#include <array> +#include <iomanip> +#include <ios> +#include <iostream> +#include <random> +#include "io/fileName.h" +#include "superParticleSystem3D.h" + +#ifndef M_PI +#define M_PI 3.14159265358979323846 +#endif + +namespace olb { + +template<typename T, template<typename U> class PARTICLETYPE> +SuperParticleSystem3D<T, PARTICLETYPE>::SuperParticleSystem3D( + CuboidGeometry3D<T>& cuboidGeometry, LoadBalancer<T>& loadBalancer, + SuperGeometry3D<T>& superGeometry) + : SuperStructure3D<T>(cuboidGeometry, loadBalancer, + superGeometry.getOverlap()), + clout(std::cout, "SuperParticleSystem3d"), + _superGeometry(superGeometry), + _overlap(0) +{ + init(); +} + +template<typename T, template<typename U> class PARTICLETYPE> +SuperParticleSystem3D<T, PARTICLETYPE>::SuperParticleSystem3D( + SuperGeometry3D<T>& superGeometry) + : SuperStructure3D<T>(superGeometry.getCuboidGeometry(), + superGeometry.getLoadBalancer(), + superGeometry.getOverlap()), + clout(std::cout, "SuperParticleSystem3d"), + _superGeometry(superGeometry), + _overlap(0) +{ + init(); +} + +template<typename T, template<typename U> class PARTICLETYPE> +void SuperParticleSystem3D<T, PARTICLETYPE>::init() +{ + int rank = 0; + int size = 0; + + _stopSorting = 0; + +#ifdef PARALLEL_MODE_MPI + rank = singleton::mpi().getRank(); + size = singleton::mpi().getSize(); +#endif + for (int i = 0; i < this->_cuboidGeometry.getNc(); ++i) { + //clout << i << std::endl; + if (this->_loadBalancer.rank(i) == rank) { + auto dummy = new ParticleSystem3D<T, PARTICLETYPE>(i, _superGeometry); + this->_cuboidGeometry.get(i).getOrigin().toStdVector()[0]; + std::vector<T> physPos = this->_cuboidGeometry.get(i).getOrigin().toStdVector(); + std::vector<T> physExtend(3, 0); + T physR = this->_cuboidGeometry.get(i).getDeltaR(); + for (int j = 0; j < 3; j++) { + physPos[j] -= .5 * physR; + physExtend[j] = (this->_cuboidGeometry.get(i).getExtend()[j] + 1) + * physR; + } + dummy->setPosExt(physPos, physExtend); + _pSystems.push_back(dummy); + } + } + //singleton::exit(0); + +#ifdef PARALLEL_MODE_MPI + for (int i = 0; i < this->_cuboidGeometry.getNc(); ++i) { + if (this->_loadBalancer.rank(i) == rank) { + std::vector<int> dummy; + this->getCuboidGeometry().getNeighbourhood(i, dummy, 3); + _rankNeighbours.insert(_rankNeighbours.end(), dummy.begin(), dummy.end()); + _cuboidNeighbours.push_back(dummy); + } + } + for (auto& N : _rankNeighbours) { + N = this->_loadBalancer.rank(N); + } +#endif + + /* Ein jeder ist sein eigener Nachbar*/ + if (rank == 0) { + _rankNeighbours.push_back(size - 1); + } + if (rank == size - 1) { + _rankNeighbours.push_back(0); + } + _rankNeighbours.push_back(rank); + _rankNeighbours.sort(); + _rankNeighbours.unique(); +} + +template<typename T, template<typename U> class PARTICLETYPE> +SuperParticleSystem3D<T, PARTICLETYPE>::SuperParticleSystem3D( + SuperParticleSystem3D<T, PARTICLETYPE>& spSys) + : SuperStructure3D<T>(spSys._cuboidGeometry, spSys._loadBalancer, + int(spSys._overlap)), + clout(std::cout, "SuperParticleSystem3d"), + _pSystems(spSys._pSystems), + _superGeometry(spSys._superGeometry), + _rankNeighbours(spSys._rankNeighbours), + _overlap(spSys._overlap) +{ +} + +template<typename T, template<typename U> class PARTICLETYPE> +SuperParticleSystem3D<T, PARTICLETYPE>::SuperParticleSystem3D( + SuperParticleSystem3D<T, PARTICLETYPE> const& spSys) + : SuperStructure3D<T>(spSys._cuboidGeometry, spSys._loadBalancer, + int(spSys._overlap)), + clout(std::cout, "SuperParticleSystem3d"), + _pSystems(spSys._pSystems), + _superGeometry(spSys._superGeometry), + _rankNeighbours(spSys._rankNeighbours), + _overlap(spSys._overlap) +{ +} + +template<typename T, template<typename U> class PARTICLETYPE> +SuperParticleSystem3D<T, PARTICLETYPE>::SuperParticleSystem3D( + SuperParticleSystem3D<T, PARTICLETYPE> && spSys) + : SuperStructure3D<T>(spSys._cuboidGeometry, spSys._loadBalancer, + int(spSys._overlap)), + clout(std::cout, "SuperParticleSystem3d"), + _superGeometry(spSys._superGeometry), + _rankNeighbours(spSys._rankNeighbours), + _overlap(spSys._overlap) +{ + _pSystems = std::move(spSys._pSystems); +} + +template<typename T, template<typename U> class PARTICLETYPE> +void SuperParticleSystem3D<T, PARTICLETYPE>::print() +{ + int no = globalNumOfParticles(); + int active = globalNumOfActiveParticles(); + clout << "activeParticles= " << active << " (" << no << ") " << std::endl; + // cout << "[SuperParticleSystem3D] " << _pSystems.size() + // << " pSystems on rank " << singleton::mpi().getRank() << "\n"; + // for (auto pS : _pSystems) { + // cout << pS->_particles.size() << " "; + // } +} + +template<typename T, template<typename U> class PARTICLETYPE> +void SuperParticleSystem3D<T, PARTICLETYPE>::printDeep(std::string message) +{ + clout << "=========================================================================================================================" << std::endl; + clout << "printDeep diagnostic tool" << message << std::endl; + clout << "=========================================================================================================================" << std::endl; + for (auto pS : _pSystems) { + pS->printDeep ( std::to_string(_pSystems.size()) + " pSystems on rank " + std::to_string(singleton::mpi().getRank()) + ": " ); + } +} + +template<typename T, template<typename U> class PARTICLETYPE> +void SuperParticleSystem3D<T, PARTICLETYPE>::print(std::list<int> mat) +{ + std::list<int>::iterator _matIter; + int no = globalNumOfParticles(); + clout << "globalNumOfParticles=" << no; + int active = globalNumOfActiveParticles(); + clout << "; activeParticles=" << active; + for (_matIter = mat.begin(); _matIter != mat.end(); _matIter++) { + clout << "; material" << *_matIter << "=" << countMaterial(*_matIter); + } + clout << std::endl; +} + +template<typename T, template<typename U> class PARTICLETYPE> +void SuperParticleSystem3D<T, PARTICLETYPE>::captureEscapeRate( + std::list<int> mat) +{ + std::list<int>::iterator _matIter; + T sum = T(); + // capture rate + for (_matIter = mat.begin(); _matIter != mat.end(); _matIter++) { + sum += (T) countMaterial(*_matIter); + } + clout << "captureRate=" << 1. - sum / globalNumOfParticles() + << "; escapeRate=" << sum / globalNumOfParticles() << std::endl; +} + +template<typename T, template<typename U> class PARTICLETYPE> +void SuperParticleSystem3D<T, PARTICLETYPE>::diffEscapeRate( + std::list<int> mat, int& globalPSum, int& pSumOutlet, T& diffEscapeRate, T& maxDiffEscapeRate, + int iT, int iTConsole, T genPartPerTimeStep) +{ + std::list<int>::iterator _matIter; + T pSumOutletNew = T(); + T maxDiffEscapeRateNew = maxDiffEscapeRate; + T diffERtmp = T(); // temporal differencial escape rate + + // count particles at outlet + for (_matIter = mat.begin(); _matIter != mat.end(); _matIter++) { + pSumOutletNew += (T) countMaterial(*_matIter); + } + + // calculate diff. escape rate + if (globalNumOfParticles() > globalPSum) { + T diffERtmp = (T) (pSumOutletNew - pSumOutlet) / (globalNumOfParticles() - globalPSum); + diffEscapeRate += diffERtmp; + } else { + if (genPartPerTimeStep != 0.) { + diffERtmp = (T) (pSumOutletNew - pSumOutlet) / (genPartPerTimeStep); + diffEscapeRate += diffERtmp; + } + } + // calculate max. diff. escape rate + if (diffERtmp > maxDiffEscapeRateNew) { + maxDiffEscapeRateNew = diffERtmp; + } + // console output + if (iT % iTConsole == 0) { + diffEscapeRate /= iTConsole; + if (globalNumOfParticles() > globalPSum) { + clout << "diffEscapeRate = " << diffEscapeRate << std::endl; + } else { + if (genPartPerTimeStep != 0.) { + clout << "no particles in feedstream, continue calculation of diffEscapeRate with theoretical " + << genPartPerTimeStep << " generated particles per phys. time step" + << " diffEscapeRate = " << diffEscapeRate << std::endl; + } else { + clout << "no particles in feedstream, calculation of diffEscapeRate not possible" << std::endl; + } + } + if (maxDiffEscapeRateNew > maxDiffEscapeRate) { + clout << "maxDiffEscapeRate = " << maxDiffEscapeRateNew << std::endl; + } + diffEscapeRate = 0. ; + } + maxDiffEscapeRate = maxDiffEscapeRateNew; + pSumOutlet = pSumOutletNew; + globalPSum = globalNumOfParticles(); +} + +template<typename T, template<typename U> class PARTICLETYPE> +void SuperParticleSystem3D<T, PARTICLETYPE>::diffEscapeRate( + std::list<int> mat, int& globalPSum, int& pSumOutlet, T& diffEscapeRate, T& maxDiffEscapeRate, + int iT, int iTConsole, T genPartPerTimeStep, + T& avDiffEscapeRate, T latticeTimeStart, T latticeTimeEnd) +{ + std::list<int>::iterator _matIter; + T pSumOutletNew = T(); + T maxDiffEscapeRateNew = maxDiffEscapeRate; + T avDiffEscapeRateNew = T(); + + // count particle at outlet + for (_matIter = mat.begin(); _matIter != mat.end(); _matIter++) { + pSumOutletNew += (T) countMaterial(*_matIter); + } + // calculate diff. escape rate + if (globalNumOfParticles() > globalPSum) { + avDiffEscapeRateNew = (T) (pSumOutletNew - pSumOutlet) / (globalNumOfParticles() - globalPSum); + diffEscapeRate += avDiffEscapeRateNew; + } else { + if (genPartPerTimeStep != 0.) { + avDiffEscapeRateNew = (T) (pSumOutletNew - pSumOutlet) / (genPartPerTimeStep); + diffEscapeRate += avDiffEscapeRateNew; + } + } + // calculate max. diff. escape rate + if (avDiffEscapeRateNew > maxDiffEscapeRateNew) { + maxDiffEscapeRateNew = avDiffEscapeRateNew; + } + // calculate average diff. escape rate between tStart and tEnd + if ((iT >= latticeTimeStart) && (iT <= latticeTimeEnd)) { + avDiffEscapeRate += avDiffEscapeRateNew; + } + if (iT == latticeTimeEnd) { + avDiffEscapeRate /= (latticeTimeEnd - latticeTimeStart); + clout << "average diffEscapeRate between t = " << latticeTimeStart << "s and t = " << latticeTimeEnd << "s : " + << avDiffEscapeRate << std::endl; + } + // console output + if (iT % iTConsole == 0) { + diffEscapeRate /= iTConsole; + if (globalNumOfParticles() > globalPSum) { + clout << "diffEscapeRate = " << diffEscapeRate << std::endl; + } else { + if (genPartPerTimeStep != 0.) { + clout << "no particles in feedstream, continue calculation of diffEscapeRate with theoretical " + << genPartPerTimeStep << " generated particles per phys. time step" << std::endl; + clout << "diffEscapeRate = " << diffEscapeRate << std::endl; + } else { + clout << "no particles in feedstream, calculation of diffEscapeRate not possible" << std::endl; + } + } + if (maxDiffEscapeRateNew > maxDiffEscapeRate) { + clout << "maxDiffEscapeRate = " << maxDiffEscapeRateNew << std::endl; + } + diffEscapeRate = 0. ; + } + maxDiffEscapeRate = maxDiffEscapeRateNew; + pSumOutlet = pSumOutletNew; + globalPSum = globalNumOfParticles(); +} + +// Output_txt.file for tracer particles (particles that have an id != 0 +/* _filename is filename of output txt-file + * _file is contact to output file + * */ +template<typename T, template<typename U> class PARTICLETYPE> +void SuperParticleSystem3D<T, PARTICLETYPE>::getOutput(std::string _filename, + int it, T conversionFactorTime, unsigned short _properties ) +{ + + enum particleProperties + : unsigned short {position = 1, + velocity = 2, + radius = 4, + mass = 8, + force = 16, + storeForce = 32, + }; + + /// initialize Parameters for the desired ParticleOutput: + /// ID, Radius, Position, Velocity, Force, storeForces + std::vector < T > id, rad; + std::vector<std::array<T, 3>> pos, vel, forces, storeForces; + /// initialize iteration Parameters for the loops + int k, j, numOfTracerParticles = globalNumOfTracerParticles(); + int m = 0; + /// setPrecision for the decimal places in the txtOutputFile + std::setprecision(9); + /// set size of list in correlation to the number of TracerParticles + for (k = 0; k < numOfTracerParticles; k++) { + id.push_back(T()); + rad.push_back(T()); + pos.push_back( { T(), T(), T() }); + vel.push_back( { T(), T(), T() }); + forces.push_back( { T(), T(), T() }); + storeForces.push_back( { T(), T(), T() }); + } + + k = 0; + /// get Information about each Particle in each particleSystem + std::vector<ParticleSystem3D<T, PARTICLETYPE>*> _pSystems = getPSystems(); + for (unsigned int pS = 0; pS < _pSystems.size(); ++pS) { + std::deque<PARTICLETYPE<T>*> particles = + _pSystems[pS]->getParticlesPointer(); + for (auto p : particles) { + if (p->getID() != 0) { + /// check whether it is a tracerParticle(!=0) or not (==0) + id[k] = p->getID(); + rad[k] = p->getRad(); + for (j = 0; j < 3; j++) { /// loop for each coordinate [x,y,z] + pos[k][j] = p->getPos()[j]; + vel[k][j] = p->getVel()[j]; + forces[k][j] = p->getForce()[j]; + storeForces[k][j] = p->getStoreForce()[j]; + } + // p->resetStoreForce(); + k++; + p++; + } + } + } + +#ifdef PARALLEL_MODE_MPI + /// prepare data vectors for parallel computing (mpi) + for (m = 0; m < numOfTracerParticles; m++) { + singleton::mpi().reduceAndBcast(id[m], MPI_SUM); + singleton::mpi().reduceAndBcast(rad[m], MPI_SUM); + for (j = 0; j < 3; j++) { /// loop for each coordinate [x,y,z] + singleton::mpi().reduceAndBcast(pos[m][j], MPI_SUM); + singleton::mpi().reduceAndBcast(vel[m][j], MPI_SUM); + singleton::mpi().reduceAndBcast(forces[m][j], MPI_SUM); + singleton::mpi().reduceAndBcast(storeForces[m][j], MPI_SUM); + } + } +#endif + + /// Particle Output to txtFile just on main process ID (==0) + if (!singleton::mpi().getRank()) { + std::ofstream _file; + + int i = 0; + if (it == 0) { + /// write headers of each column at timeStep zero + /// after header name, there is the number of the column + _file.open(_filename, std::ios::out | std::ios::trunc); + _file << "Timestep" << i + 1 << " " << "physTime" << i + 2; + i = i + 2; + + for (m = 0; m < numOfTracerParticles; m++) { + _file << " id" << i + 1; + i = i + 1; + _file << " rad" << i + 1; + i = i + 1; + if (_properties & particleProperties::position) { + _file << " pos0_" << i + 1 << " pos1_" << " pos2_" << i + 3; + i = i + 3; + } + if (_properties & particleProperties::force) { + _file << " forc0_" << i + 1 << " forc1_" << i + 2 << " forc2_" + << i + 3; + i = i + 3; + } + if (_properties & particleProperties::velocity) { + _file << " vel0_" << i + 1 << " vel1_" << i + 2 << " vel2_" << i + 3; + i = i + 3; + } + if (_properties & particleProperties::storeForce) { + _file << " hforc0_" << i + 1 << " hforc1_" << i + 2 << " hforc2_" + << i + 3; + i = i + 3; + } + } + _file << "\n"; + _file.close(); + } + + if (it >= 0) { + /// write the results in regart to the above defined headers + _file.open(_filename, std::ios::out | std::ios::app); + _file << it << " " << conversionFactorTime*it << " "; + + for (m = 0; m < numOfTracerParticles; m++) { + _file << id[m] << " "; + _file << rad[m] << " "; + if (_properties & particleProperties::position) { + _file << pos[m][0] << " " << pos[m][1] << " " << pos[m][2] << " "; + } + if (_properties & particleProperties::force) { + _file << forces[m][0] << " " << forces[m][1] << " " << forces[m][2] + << " "; + } + if (_properties & particleProperties::velocity) { + _file << vel[m][0] << " " << vel[m][1] << " " << vel[m][2] << " "; + } + if (_properties & particleProperties::storeForce) { + _file << storeForces[m][0] << " " << storeForces[m][1] << " " + << storeForces[m][2] << " "; + } + } + } + _file << "\n"; + _file.close(); + } +} + + +template<typename T, template<typename U> class PARTICLETYPE> +std::vector<ParticleSystem3D<T, PARTICLETYPE>*>& SuperParticleSystem3D<T, + PARTICLETYPE>::getPSystems() //+* +{ + return _pSystems; +} + +template<typename T, template<typename U> class PARTICLETYPE> +void SuperParticleSystem3D<T, PARTICLETYPE>::simulate(T dT, bool scale) +{ + + // for (auto pS : _pSystems) { + // pS->velocityVerlet1(dT); + // if (_overlap > 0) { + // pS->makeTree(); + // cout << "makeTree" << endl; + // } + // } + // updateParticleDistribution(); + // for (auto pS : _pSystems) { + // pS->velocityVerlet2(dT); + // } + + // updateParticleDistribution(); + // for (auto pS : _pSystems) { + // pS->computeForce(); + // pS->predictorCorrector1(dT); + // } + // + // updateParticleDistribution(); + // for (auto pS : _pSystems) { + // pS->computeForce(); + // pS->predictorCorrector2(dT); + // } + + // for (auto pS : _pSystems) { + // pS->rungeKutta4_1(dT); + // } + // updateParticleDistribution(); + // for (auto pS : _pSystems) { + // pS->rungeKutta4_2(dT); + // } + // updateParticleDistribution(); + // for (auto pS : _pSystems) { + // pS->rungeKutta4_3(dT); + // } + // updateParticleDistribution(); + // for (auto pS : _pSystems) { + // pS->rungeKutta4_4(dT); + // } + // updateParticleDistribution(); + for (auto pS : _pSystems) { + time_t delta = clock(); + pS->_contactDetection->sort(); + _stopSorting += clock() - delta; + pS->simulate(dT, scale); + pS->computeBoundary(); + + } + updateParticleDistribution(); +} + +template<typename T, template<typename U> class PARTICLETYPE> +void SuperParticleSystem3D<T, PARTICLETYPE>::simulateWithTwoWayCoupling_Mathias ( T dT, + ForwardCouplingModel<T,PARTICLETYPE>& forwardCoupling, + BackCouplingModel<T,PARTICLETYPE>& backCoupling, + int material, int subSteps, bool resetExternalField, bool scale ) +{ + // reset external field + if (resetExternalField) + backCoupling.resetExternalField(material); + + for (auto pS : _pSystems) { + time_t delta = clock(); + pS->_contactDetection->sort(); + _stopSorting += clock() - delta; + pS->simulateWithTwoWayCoupling_Mathias(dT, forwardCoupling, backCoupling, material, subSteps, scale); + pS->computeBoundary(); + + } + updateParticleDistribution(); +} + +template<typename T, template<typename U> class PARTICLETYPE> +void SuperParticleSystem3D<T, PARTICLETYPE>::simulateWithTwoWayCoupling_Davide ( T dT, + ForwardCouplingModel<T,PARTICLETYPE>& forwardCoupling, + BackCouplingModel<T,PARTICLETYPE>& backCoupling, + int material, int subSteps, bool resetExternalField, bool scale ) +{ + // reset external field + if (resetExternalField) + backCoupling.resetExternalField(material); + + for (auto pS : _pSystems) { + time_t delta = clock(); + pS->_contactDetection->sort(); + _stopSorting += clock() - delta; + pS->simulateWithTwoWayCoupling_Davide(dT, forwardCoupling, backCoupling, material, subSteps, scale); + pS->computeBoundary(); + } + updateParticleDistribution(); +} + +// multiple collision models +template<> +void SuperParticleSystem3D<double, MagneticParticle3D>::simulate(double dT, std::set<int> sActivityOfParticle, bool scale) +{ + for (auto pS : _pSystems) { + time_t delta = clock(); + if (pS->getIGeometry() == singleton::mpi().getRank()) { + pS->_contactDetection->sort(); + } + _stopSorting += clock() - delta; + if (pS->getIGeometry() == singleton::mpi().getRank()) { + pS->simulate(dT, sActivityOfParticle, scale); + pS->computeBoundary(); + } + } + updateParticleDistribution(); +} + +template<> +bool SuperParticleSystem3D<double, MagneticParticle3D>::particleSActivityTest(int sActivity) +{ + for (auto pS : _pSystems) { + for (auto p : pS->_particles) { + if (p.getSActivity() == sActivity) { + return false; + } + } + } + return true; +} + +template<typename T, template<typename U> class PARTICLETYPE> +void SuperParticleSystem3D<T, PARTICLETYPE>::setOverlap(T overlap) +{ + if ( int(overlap) + 1 > _superGeometry.getOverlap() ) { + clout << "Overlap of SuperParticleSystem3D should be < overlap " + "of SuperStructure3D" << std::endl; + exit(-1); + } + _overlap = overlap; +} + +template<typename T, template<typename U> class PARTICLETYPE> +T SuperParticleSystem3D<T, PARTICLETYPE>::getOverlap() +{ + return _overlap; +} + +template<typename T, template<typename U> class PARTICLETYPE> +int SuperParticleSystem3D<T, PARTICLETYPE>::numOfPSystems() +{ + return _pSystems.size(); +} + +template<typename T, template<typename U> class PARTICLETYPE> +int SuperParticleSystem3D<T, PARTICLETYPE>::globalNumOfParticles() +{ +#ifdef PARALLEL_MODE_MPI + // cout << "return1" << endl; + int buffer = rankNumOfParticles(); + // cout << "return2" << endl; + singleton::mpi().reduceAndBcast(buffer, MPI_SUM); + // cout << "return3" << endl; + return buffer; +#else + // cout << "return4" << endl; + return rankNumOfParticles(); +#endif +} + +template<typename T, template<typename U> class PARTICLETYPE> +int SuperParticleSystem3D<T, PARTICLETYPE>::globalNumOfShadowParticles() +{ +#ifdef PARALLEL_MODE_MPI + int buffer = rankNumOfShadowParticles(); + singleton::mpi().reduceAndBcast(buffer, MPI_SUM); + return buffer; +#else + return rankNumOfShadowParticles(); +#endif +} + +template<typename T, template<typename U> class PARTICLETYPE> +int SuperParticleSystem3D<T, PARTICLETYPE>::rankNumOfParticles() +{ + int num = 0; + for (auto pS : _pSystems) { + num += pS->size(); + } + return num; +} + +template<typename T, template<typename U> class PARTICLETYPE> +int SuperParticleSystem3D<T, PARTICLETYPE>::rankNumOfShadowParticles() +{ + int num = 0; + for (auto pS : _pSystems) { + num += pS->_shadowParticles.size(); + } + return num; +} + +template<typename T, template<typename U> class PARTICLETYPE> +int SuperParticleSystem3D<T, PARTICLETYPE>::rankNumOfActiveParticles() +{ + int num = 0; + for (auto pS : _pSystems) { + num += pS->numOfActiveParticles(); + } + return num; +} + +template<typename T, template<typename U> class PARTICLETYPE> +int SuperParticleSystem3D<T, PARTICLETYPE>::countLocMaterial(int mat) +{ + int num = 0; + for (auto pS : _pSystems) { + num += pS->countMaterial(mat); + } + return num; +} + +template<typename T, template<typename U> class PARTICLETYPE> +int SuperParticleSystem3D<T, PARTICLETYPE>::countMaterial(int mat) +{ +#ifdef PARALLEL_MODE_MPI + int buffer = countLocMaterial(mat); + singleton::mpi().reduceAndBcast(buffer, MPI_SUM); + return buffer; +#else + return countLocMaterial(mat); +#endif +} + +template<typename T, template<typename U> class PARTICLETYPE> +int SuperParticleSystem3D<T, PARTICLETYPE>::globalNumOfActiveParticles() +{ +#ifdef PARALLEL_MODE_MPI + int buffer = rankNumOfActiveParticles(); + singleton::mpi().reduceAndBcast(buffer, MPI_SUM); + return buffer; +#else + return rankNumOfActiveParticles(); +#endif +} + +template<typename T, template<typename U> class PARTICLETYPE> +int SuperParticleSystem3D<T, PARTICLETYPE>::globalNumOfTracerParticles() +{ +#if PARALLEL_MODE_MPI + int buffer = rankNumOfTracerParticles(); + singleton::mpi().reduceAndBcast(buffer, MPI_SUM); + return buffer; +#else + return rankNumOfTracerParticles(); +#endif +} + +template<typename T, template<typename U> class PARTICLETYPE> +int SuperParticleSystem3D<T, PARTICLETYPE>::rankNumOfTracerParticles() +{ + int num = 0; + for (auto pS : _pSystems) { + std::deque<PARTICLETYPE<T>*> particles = pS->getParticlesPointer(); + int pSNum = 0; + for (auto p : particles) { + if (p->getID() != 0) { + pSNum++; + } + } + num += pSNum; + } + return num; +} + +// TODO class olb::ParticleSystem3D<double, olb::Particle3D>’ has no member named ‘radiusDistribution +/* + template<typename T, template<typename U> class PARTICLETYPE> + std::map<T, int> SuperParticleSystem3D<T, PARTICLETYPE>::rankRadiusDistribution() + { + std::map<T, int> distrAll; + typename std::map<T, int>::iterator ita; + for (auto pS : _pSystems) { + std::map<T, int> distrPs; + distrPs = pS->radiusDistribution(); + for (typename std::map<T, int>::iterator itp = distrPs.begin(); + itp != distrPs.end(); ++itp) { + T radKeyP = itp->first; + int countP = itp->second; + if (distrAll.count(radKeyP) > 0) { + ita = distrAll.find(radKeyP); + ita->second = ita->second + countP; + } else { + distrAll.insert(std::pair<T, int>(radKeyP, countP)); + } + } + } + return distrAll; + } + */ + +template<typename T, template<typename U> class PARTICLETYPE> +std::vector<int> SuperParticleSystem3D<T, PARTICLETYPE>::numOfForces() +{ + std::vector<int> dummy; + for (auto pS : _pSystems) { + dummy.push_back(pS->numOfForces()); + } + return dummy; +} + +template<typename T, template<typename U> class PARTICLETYPE> +template<typename DESCRIPTOR> +void SuperParticleSystem3D<T, PARTICLETYPE>::setVelToFluidVel( + SuperLatticeInterpPhysVelocity3D<T, DESCRIPTOR>& fVel) +{ + for (auto pS : _pSystems) { + pS->setVelToFluidVel(fVel); + } +} + +template<typename T, template<typename U> class PARTICLETYPE> +void SuperParticleSystem3D<T, PARTICLETYPE>::setVelToAnalyticalVel( + AnalyticalConst3D<T, T>& aVel) +{ + for (auto pS : _pSystems) { + pS->setVelToAnalyticalVel(aVel); + } +} + +template<typename T, template<typename U> class PARTICLETYPE> +bool SuperParticleSystem3D<T, PARTICLETYPE>::findCuboid(PARTICLETYPE<T> &p) +{ + return findCuboid(p, int(_overlap)); +} + +template<typename T, template<typename U> class PARTICLETYPE> +bool SuperParticleSystem3D<T, PARTICLETYPE>::findCuboid(PARTICLETYPE<T> &p, + int overlap) +{ + int C = this->_cuboidGeometry.get_iC(p.getPos()[0], p.getPos()[1], + p.getPos()[2], overlap); + if (C != this->_cuboidGeometry.getNc()) { + p.setCuboid(C); + return 1; + } else { + clout << "Lost Particle! Pos: " << p.getPos()[0] << " " << p.getPos()[1] + << " " << p.getPos()[2] << " Vel: " << p.getVel()[0] << " " + << p.getVel()[1] << " " << p.getVel()[2] << " Cuboid: " << C << std::endl; + p.setActive(false); + return 0; + } +} + +template<typename T, template<typename U> class PARTICLETYPE> +bool SuperParticleSystem3D<T, PARTICLETYPE>::checkCuboid(PARTICLETYPE<T>& p, + T overlap) +{ + return checkCuboid(p, overlap, p.getCuboid()); +} + +template<typename T, template<typename U> class PARTICLETYPE> +bool SuperParticleSystem3D<T, PARTICLETYPE>::checkCuboid(PARTICLETYPE<T>& p, + T overlap, int iC) +{ + // Checks whether particle is contained in the cuboid + // extended with an layer of size overlap*delta + return this->_cuboidGeometry.get(iC).physCheckPoint(p.getPos()[0], + p.getPos()[1], + p.getPos()[2], overlap); +} + +template<typename T, template<typename U> class PARTICLETYPE> +void SuperParticleSystem3D<T, PARTICLETYPE>::updateParticleDistribution() +{ + // Find particles on wrong cuboid, store in relocate and delete + // maps particles to their new rank + // std::multimap<int, PARTICLETYPE<T> > relocate; + _relocate.clear(); +#ifdef shadows + _relocateShadow.clear(); +#endif + for (unsigned int pS = 0; pS < _pSystems.size(); ++pS) { + + _pSystems[pS]->_shadowParticles.clear(); + auto par = _pSystems[pS]->_particles.begin(); + + while (par != _pSystems[pS]->_particles.end()) { + + //Check if particle is still in its cuboid + if (checkCuboid(*par, 0)) { +#ifdef shadows + // Check if its inside boundary layer of its cuboid + if (!(checkCuboid(*par, -_overlap))) { + + std::set<int> sendTo; + // Run through all cuboids to search the one that + // shares its boundary with the cuboid containing the particle + for (int iC = 0; iC < this->_cuboidGeometry.getNc(); iC++) { + // Check if iC is neighbor cuboid in which overlap the particle is + // located + if (par->getCuboid() != iC && checkCuboid(*par, _overlap, iC)) { + // rank is the processing thread of the global cuboid number iC + int rank = this->_loadBalancer.rank(iC); + // insert rank into sendTo, if sendTo does not already contain + // rank for the actual particle + if (!sendTo.count(rank)) { + _relocateShadow.insert(std::make_pair(rank, (*par))); + sendTo.insert(rank); + } + } + } + } +#endif + //If not --> find new cuboid + } else { + // check to which cuboid particle is located + // and give new cuboid to par + findCuboid(*par, 0); + // gives particle the new cuboid number where it is now located to + _relocate.insert( + std::make_pair(this->_loadBalancer.rank(par->getCuboid()), (*par))); +#ifdef shadows + // If the new particle position is in boundary layer + // If yes --> check if its inside boundary layer + if (!(checkCuboid(*par, -_overlap))) { + // yes, particle is in boundary layer + std::set<int> sendTo; + for (int iC = 0; iC < this->_cuboidGeometry.getNc(); iC++) { + // checks if iC is neighbor cuboid in whose overlap the particle + // is located + if (par->getCuboid() != iC && checkCuboid(*par, _overlap, iC)) { + int rank = this->_loadBalancer.rank(iC); + if (!sendTo.count(rank)) { + // shadow particle (copy of original) is inserted for iC + // but just with information of rank, not of cuboid! + _relocateShadow.insert(std::make_pair(rank, (*par))); + sendTo.insert(rank); + } + } + } + } +#endif + par = _pSystems[pS]->_particles.erase(par); + par--; + } + par++; + } + } + /* Communicate number of Particles per cuboid*/ +#ifdef PARALLEL_MODE_MPI + singleton::MpiNonBlockingHelper mpiNbHelper; + //mpiNbHelper.allocate(_rankNeighbours.size()); + + int k = 0; + + /* Serialize particles */ + // it is: std::map<int, std::vector<double> > _send_buffer + _send_buffer.clear(); + + // fill std::map<int, std::vector<double> > _send_buffer with + // particle properties of particles in _relocate + + // buffer size equals number of particle properties + T buffer[PARTICLETYPE<T>::serialPartSize]; + // insert all regular particles of _relocate into _send_buffer + for (auto rN : _relocate) { + // second element in rN is (*par) + // write particle properties into buffer + rN.second.serialize(buffer); + // first element in rN is rank + // it is: std::map<int, std::vector<double> > _send_buffer; + // enlarge the _send_buffer element (rank), which is of vector type, with + // new elements of buffer (begin: buffer, end: buffer + serialPartSize) + _send_buffer[rN.first].insert(_send_buffer[rN.first].end(), + buffer, buffer + PARTICLETYPE<T>::serialPartSize); + } + + /*Send Particles */ + k = 0; + // noSends contains number of neighboring cuboids + int noSends = 0; + // it is: std::list<int> _rankNeighbours, rank of neighboring cuboids + for (auto rN : _rankNeighbours) { + // if there are particles contained in _send_buffer, increase noSends + if (_send_buffer[rN].size() > 0) { + ++noSends; + } + } + // int noSends = _send_buffer.size(); + if (noSends > 0) { + mpiNbHelper.allocate(noSends); + // cout << mpiNbHelper.get_size() << std::endl; + for (auto rN : _rankNeighbours) { + if (_send_buffer[rN].size() > 0) { + singleton::mpi().iSend<double>(&_send_buffer[rN][0], + _relocate.count(rN)*PARTICLETYPE<T>::serialPartSize, + rN, &mpiNbHelper.get_mpiRequest()[k++], 1); + } + } + } + + /*Receive and add particles*/ + singleton::mpi().barrier(); + k = 0; + int flag = 0; + MPI_Iprobe(MPI_ANY_SOURCE, 1, MPI_COMM_WORLD, &flag, MPI_STATUS_IGNORE); + if (flag) { + for (auto rN : _rankNeighbours) { + MPI_Status status; + int tmpFlag = 0; + MPI_Iprobe(rN, 1, MPI_COMM_WORLD, &tmpFlag, &status); + if (tmpFlag) { + int number_amount = 0; + MPI_Get_count(&status, MPI_DOUBLE, &number_amount); + T recv_buffer[number_amount]; + singleton::mpi().receive(recv_buffer, number_amount, rN, 1); + + for (int iPar = 0; iPar < number_amount; iPar += PARTICLETYPE<T>::serialPartSize) { + PARTICLETYPE<T> p; + p.unserialize(&recv_buffer[iPar]); + if (singleton::mpi().getRank() == this->_loadBalancer.rank(p.getCuboid())) { + // particle added to pSystem with local cuboid number on actual thread + _pSystems[this->_loadBalancer.loc(p.getCuboid())]->addParticle(p); + |