27#ifndef PARTICLE_DISTRIBUTION_H
28#define PARTICLE_DISTRIBUTION_H
31 template <
typename T,
typename S>
33 this->
getName() =
"LogNDistribution";
36 template <
typename T,
typename S>
43 template <
typename T,
typename S>
45 this->
getName() =
"GaussDistribution";
48 template <
typename T,
typename S>
58 template <
typename T,
typename S>
61 _distribution(std::move(distribution)), _numberOfParticles(numberOfParticles), _numberOfParticleSizeGroups(numberOfParticleSizeGroups), _minRadius(minRadius), _maxRadius(maxRadius) {
62 _sizeCalculation =
new T * [_numberOfParticleSizeGroups];
63 const T delta_d = (_maxRadius - _minRadius) / _numberOfParticleSizeGroups;
64 for (
long long i = 0; i < _numberOfParticleSizeGroups; i++) {
65 _sizeCalculation[i] =
new T[4];
69 for (
long long i = 0; i < _numberOfParticleSizeGroups; i++) {
71 _sizeCalculation[i][0] = (_minRadius + delta_d * i) * lengthConversion;
72 input[0] = _sizeCalculation[i][0] / lengthConversion;
74 _distribution->operator()(output, input);
75 _sizeCalculation[i][1] = output[0];
86 template <
typename T,
typename S>
89 for (
long long i = 0; i < _numberOfParticleSizeGroups; i++) {
91 sum_np += _sizeCalculation[i][1];
94 for (
long long i = 0; i < _numberOfParticleSizeGroups; i++) {
96 _sizeCalculation[i][2] = _sizeCalculation[i][1] / sum_np;
99 _sizeCalculation[i][3] = round(_sizeCalculation[i][2] * (_numberOfParticles + _particleOffset));
102 long long p_count = 0;
103 if (getExactParticleNumber) {
104 for (
long long i = 0; i < _numberOfParticleSizeGroups; i++) {
105 p_count += _sizeCalculation[i][3];
108 if (p_count == _numberOfParticles || getExactParticleNumber ==
false) {
110 for (
long long i = 0; i < _numberOfParticleSizeGroups; i++) {
111 if (_sizeCalculation[i][3] > 0) {
112 for (
long long j = 1; j <= _sizeCalculation[i][3]; j++) {
114 _particleArray.push_back(_sizeCalculation[i][0]);
122 _particleArray.clear();
125 _particleOffset += (_numberOfParticles - p_count)/7.;
127 return this->calculateParticleArray(
true);
130 template <
typename T,
typename S>
132 particleArray = _particleArray;
137 template <
typename T,
typename S>
139 if (_currentParticle == 0 && shuffleParticles) {
140 shuffleParticleArray();
141 T Radius = _particleArray[_currentParticle];
145 else if (_currentParticle == (_numberOfParticles)) {
147 T Radius = _particleArray[_currentParticle];
148 _currentParticle = 0;
149 if (shuffleParticles) {
150 shuffleParticleArray();
159 T Radius = _particleArray[_currentParticle];
165 template <
typename T,
typename S>
167 random_shuffle(&_particleArray[0], &_particleArray[_particleArray.size() - 1]);
170 template <
typename T,
typename S>
174 clout << std::setw(7) <<
"GroupNo" << std::setw(12) <<
"Radius" << std::setw(12) <<
"concentr." << std::setw(12) <<
"partiton" << std::setw(12) <<
"p_number" << std::endl;
175 for (
long long i = 0; i < _numberOfParticleSizeGroups; i++) {
176 clout << std::setw(7) << i + 1;
177 for (
long long j = 0; j <= 3; j++) {
178 clout << std::setw(12) << _sizeCalculation[i][j];
184 template <
typename T,
typename S>
188 clout << std::setw(7) <<
"p_number" << std::setw(12) <<
"radius" << std::endl;
189 for (
long long i = 0; i < _particleArray.size(); i++) {
190 clout << std::setw(7) << i << std::setw(12) << _particleArray[i] << std::endl;
198 template <
typename T,
typename S>
199 template <
typename DESCRIPTOR>
209 timeDistribution(output, input);
211 integral += output[0];
216 timeDistribution(output, input);
217 output[0] *= ((_numberOfParticles + _particleOffsetTime) / integral);
222 for (
int i = 1; i <= nSum; ++i) {
224 _timeArray.push_back(iT);
231 if ((
long long)(_timeArray.size()) == _numberOfParticles) {
237 _particleOffsetTime += (_numberOfParticles - _timeArray.size())/7.;
238 calculateTimeArray(timeDistribution, converter, begin, end);
244 template <
typename T,
typename S>
247 clout << std::setw(7) <<
"p_number" << std::setw(12) <<
"time" << std::endl;
248 for (
long long i = 0; i < _timeArray.size(); i++) {
249 clout << std::setw(7) << i << std::setw(12) << _timeArray[i] << std::endl;
258 template <
typename T,
typename S>
259 template <
typename C>
261 while (_timeArray[_currentParticle] == iT) {
262 T p_radius = nextParticleRadius(loopAround, shuffleParticles);
272 template <
typename T,
typename S>
273 template <
typename C>
276 while (_timeArray[_currentParticle] == iT) {
278 for (
unsigned int pSystems = 0; pSystems < supParticleSystem.
getParticleSystems().size(); pSystems++) {
280 for (
int particles = 0; particles < supParticleSystem.
getParticleSystems()[pSystems]->size(); particles++) {
281 unsigned int currentPiD = supParticleSystem.
getParticleSystems()[pSystems]->getParticlesPointer()[particles]->getID();
283 if (currentPiD == _currentParticle + 1 + idOfset) {
285 supParticleSystem.
getParticleSystems()[pSystems]->getParticlesPointer()[particles]->setActive(
true);
291 if(_currentParticle == _numberOfParticles){
292 _currentParticle = 0;
298 template <
typename T,
typename S>
301 for (
int pSystems = 0; pSystems < supParticleSystem.
numOfPSystems(); pSystems++) {
303 for (
int particles = 0; particles < supParticleSystem.
getParticleSystems()[pSystems]->size(); particles++) {
305 int currentPiD = supParticleSystem.
getParticleSystems()[pSystems]->getParticlesPointer()[particles]->getID();
306 if (currentPiD <= endID && currentPiD >= beginID) {
308 supParticleSystem.
getParticleSystems()[pSystems]->getParticlesPointer()[particles]->setActive(
false);
317 template <
typename T,
typename S>
320 T p_radius = nextParticleRadius(
false, shuffle);
324 p_radius = nextParticleRadius(
false, shuffle);
325 }
while (p_radius > 0);
326 _currentParticle = 0;
327 deactivateParticles(supParticleSystem, beginID, beginID + _numberOfParticles);
330 template <
typename T,
typename S>
332 for (
long long i = 0; i < _numberOfParticleSizeGroups; i++) {
333 delete[] _sizeCalculation[i];
335 delete[] _sizeCalculation;
Smart pointer for managing the various ways of passing functors around.
GaussDistribution(T standardDeviation, T mean, T scale=1)
bool operator()(T output[1], const S x[1])
std::string & getName()
read and write access to name
IndicatorF3D is an application from .
LogNormalDistribution(T standardDeviation, T geometricMean, T scale=1)
bool operator()(T output[1], const S x[1])
class for marking output with some text
void spawnSphericalParticles(SuperParticleSystem3D< T, Particle3D > &supParticleSystem, IndicatorF3D< S > &indicator, T density, C iT, bool loopAround=false, bool reshuffleAfterLoop=false)
T nextParticleRadius(bool loopAround=false, bool reshuffleAfterLoop=false)
~ParticleDistribution()
destructor
void deactivateParticles(SuperParticleSystem3D< T, Particle3D > &supParticleSystem, int beginID, int endID)
Sets the particle state 'active' to 'false' for particles with ids from 'beginID' to 'endID'.
void printParticleArray()
prints the particle Radii to the console
ParticleDistribution(FunctorPtr< AnalyticalF1D< T, S > > &&distribution, C numberOfParticles, C numberOfParticleSizeGroups, T minRadius, T maxRadius, T lengthConversion=1)
void printSizeMatrix()
prints the calculation matrix to the console
void shuffleParticleArray()
shuffles particle array using random_shuffle
bool calculateTimeArray(AnalyticalF1D< T, S > &timeDistribution, UnitConverter< T, DESCRIPTOR > &converter, T begin, T end)
bool calculateParticleArray(bool getExactParticleNumber=false)
void timeActivateParticles(SuperParticleSystem3D< T, Particle3D > &supParticleSystem, C iT, C idOfset=0)
void printTimeArray()
prints out all time-steps a new particle should appear
void preSpawnSphericalParticles(SuperParticleSystem3D< T, Particle3D > &supParticleSystem, IndicatorF3D< S > &indicator, T density, int beginID=1, bool shuffle=false)
void getParticleArray(std::vector< T > &particleArray)
The class superParticleSystem is the basis for particulate flows within OpenLB.
void addParticle(PARTICLETYPE< T > &p)
Add a Particle to SuperParticleSystem.
void addTracerParticle(IndicatorF3D< T > &ind, int idTP, T mas, T rad, int noTP=1, std::vector< T > vel={0., 0., 0.})
Add a number of particles with a certain ID (TracerParticle) equally distributed in a given Indicator...
std::vector< ParticleSystem3D< T, PARTICLETYPE > * > getParticleSystems()
Get ParticleSystems.
int numOfPSystems()
Get number of ParticleSystems.
Conversion between physical and lattice units, as well as discretization.
constexpr size_t getLatticeTime(T physTime) const
conversion from physical to lattice time
constexpr T getPhysTime(size_t latticeTime) const
conversion from lattice to physical time
ADf< T, DIM > log(const ADf< T, DIM > &a)
cpu::simd::Pack< T > pow(cpu::simd::Pack< T > base, cpu::simd::Pack< T > exp)
ADf< T, DIM > exp(const ADf< T, DIM > &a)
Top level namespace for all of OpenLB.