OpenLB 1.7
Loading...
Searching...
No Matches
ParticleDistribution.hh
Go to the documentation of this file.
1/* Lattice Boltzmann sample, written in C++, using the OpenLB
2 * library
3 *
4 * Copyright (C) 2021 Simon Berg
5 * E-mail contact: info@openlb.net
6 * The most recent release of OpenLB can be downloaded at
7 * <http://www.openlb.net/>
8 *
9 * This program is free software; you can redistribute it and/or
10 * modify it under the terms of the GNU General Public License
11 * as published by the Free Software Foundation; either version 2
12 * of the License, or (at your option) any later version.
13 *
14 * This program is distributed in the hope that it will be useful,
15 * but WITHOUT ANY WARRANTY; without even the implied warranty of
16 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
17 * GNU General Public License for more details.
18 *
19 * You should have received a copy of the GNU General Public
20 * License along with this program; if not, write to the Free
21 * Software Foundation, Inc., 51 Franklin Street, Fifth Floor,
22 * Boston, MA 02110-1301, USA.
23 */
24
26#include <utility>
27#ifndef PARTICLE_DISTRIBUTION_H
28#define PARTICLE_DISTRIBUTION_H
29
30namespace olb {
31 template <typename T, typename S>
32 LogNormalDistribution<T, S>::LogNormalDistribution(T standardDeviation, T geometricMean, T scale) :AnalyticalF1D<T, S>(1), _standardDeviation(standardDeviation), _geometricMean(geometricMean), _scale(scale) {
33 this->getName() = "LogNDistribution";
34 };
35
36 template <typename T, typename S>
37 bool LogNormalDistribution<T, S>:: operator()(T output[1], const S x[1]) {
38 output[0] = _scale / (x[0] * _standardDeviation * util::pow((2 * M_PI), 0.5)) * util::exp(-util::pow(util::log(x[0] - _geometricMean), 2.) / (2 * util::pow((_standardDeviation), 2.)));
39
40 return true;
41 };
42
43 template <typename T, typename S>
44 GaussDistribution<T, S>::GaussDistribution(T standardDeviation, T mean, T scale) :AnalyticalF1D<T, S>(1), _standardDeviation(standardDeviation), _mean(mean), _scale(scale) {
45 this->getName() = "GaussDistribution";
46 };
47
48 template <typename T, typename S>
49 bool GaussDistribution<T, S>::operator()(T output[1], const S x[1]) {
50 output[0] = _scale / (_standardDeviation * util::pow((2 * M_PI), 0.5)) * util::exp(-0.5 * util::pow((x[0] - _mean) / (_standardDeviation), 2.));
51 return true;
52 };
53
54 /* When creating a ParticleDistribution a matrix ('_sizeCalculation') is created with evenly spaced
55 * particle categories ranging from 'minRadius' to 'maxRadius' and each concentration is calculated
56 * with the given distribution
57 */
58 template <typename T, typename S>
59 template <typename C>
60 ParticleDistribution<T, S>::ParticleDistribution(FunctorPtr<AnalyticalF1D<T, S>>&& distribution, C numberOfParticles, C numberOfParticleSizeGroups, T minRadius, T maxRadius, T lengthConversion) :
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];
66 }
67 T output[1] = {};
68 T input[1] = {};
69 for (long long i = 0; i < _numberOfParticleSizeGroups; i++) {
70 // calculate size groups
71 _sizeCalculation[i][0] = (_minRadius + delta_d * i) * lengthConversion;
72 input[0] = _sizeCalculation[i][0] / lengthConversion;
73 // calculate concentration
74 _distribution->operator()(output, input);
75 _sizeCalculation[i][1] = output[0];
76 }
78 };
79 /* For the calculation of the particleArray, the total concentration is determined
80 * and for every size group the fraction (relative to all particles) and the number
81 * of particles is calculated.
82 * If the exact particle number has to be met, the function will be called again
83 * with the difference of particles in the particle array to the specified total
84 * number as an offset.
85 */
86 template <typename T, typename S>
87 bool ParticleDistribution<T, S>::calculateParticleArray(bool getExactParticleNumber) {
88 T sum_np = 0;
89 for (long long i = 0; i < _numberOfParticleSizeGroups; i++) {
90 // calculate total particle conc.
91 sum_np += _sizeCalculation[i][1];
92 }
93
94 for (long long i = 0; i < _numberOfParticleSizeGroups; i++) {
95 // calculate fraction
96 _sizeCalculation[i][2] = _sizeCalculation[i][1] / sum_np;
97 // calculate number of whole particles
98 // ofset is added to correct difference in _numberOfParticles to _particleArray..size() from rounding errors
99 _sizeCalculation[i][3] = round(_sizeCalculation[i][2] * (_numberOfParticles + _particleOffset));
100 }
101
102 long long p_count = 0;
103 if (getExactParticleNumber) {
104 for (long long i = 0; i < _numberOfParticleSizeGroups; i++) {
105 p_count += _sizeCalculation[i][3];
106 }
107 }
108 if (p_count == _numberOfParticles || getExactParticleNumber == false) {
109 p_count = 0;
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++) {
113 // write particle radii to array, where every entry represents 1 particle
114 _particleArray.push_back(_sizeCalculation[i][0]);
115 p_count++;
116 }
117 }
118 }
119 return true;
120 }
121 else {
122 _particleArray.clear();
123 // calculate offset and calls function again with offset for convergence, choose
124 // different (larger) divisor if function does not converge
125 _particleOffset += (_numberOfParticles - p_count)/7.;
126
127 return this->calculateParticleArray(true);
128 }
129 };
130 template <typename T, typename S>
131 void ParticleDistribution<T, S>::getParticleArray(std::vector<T>& particleArray) {
132 particleArray = _particleArray;
133 };
134 /* returns the next particle radius in particleArray, loops back to the first entry and/or
135 * shuffles the array before returning first particle if specified.
136 */
137 template <typename T, typename S>
138 T ParticleDistribution<T, S>::nextParticleRadius(bool loopAround, bool shuffleParticles) {
139 if (_currentParticle == 0 && shuffleParticles) {
140 shuffleParticleArray();
141 T Radius = _particleArray[_currentParticle];
142 _currentParticle++;
143 return Radius;
144 }
145 else if (_currentParticle == (_numberOfParticles)) {
146 if (loopAround) {
147 T Radius = _particleArray[_currentParticle];
148 _currentParticle = 0;
149 if (shuffleParticles) {
150 shuffleParticleArray();
151 }
152 return Radius;
153 }
154 else {
155 return 0;
156 }
157 }
158 else {
159 T Radius = _particleArray[_currentParticle];
160 _currentParticle++;
161 return Radius;
162 }
163 };
165 template <typename T, typename S>
167 random_shuffle(&_particleArray[0], &_particleArray[_particleArray.size() - 1]);
168 };
169
170 template <typename T, typename S>
173 OstreamManager clout(std::cout, "printSizeMatrix");
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];
179 }
180 clout << std::endl;
181 }
182 };
183
184 template <typename T, typename S>
187 OstreamManager clout(std::cout, "printParticleArray");
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;
191 }
192 };
193 /* To spread out the particle creation according to a timeDistribution a time array
194 * is calculated by computing the time discrete integral in the interval ['begin', 'end']
195 * and the timeDistribution is evaluated for time-step and added up until a complete particle
196 * is reached. Every time-step a complete particle should be spawned is stored in '_timeArray'.
197 */
198 template <typename T, typename S>
199 template <typename DESCRIPTOR>
201 _timeArray.clear();
202 T input[1] = { };
203 T nSum = 0;
204 T output[1] = { };
205 T integral = 0;
206 // calculate discrete integral
207 for (std::size_t iT = converter.getLatticeTime(begin); iT <= converter.getLatticeTime(end); iT++) {
208 input[0] = converter.getPhysTime(iT);
209 timeDistribution(output, input);
210
211 integral += output[0];
212 }
213 // evaluate timeDistribution for every time-step insde the time intervall
214 for (std::size_t iT = converter.getLatticeTime(begin); iT <= converter.getLatticeTime(end); iT++) {
215 input[0] = converter.getPhysTime(iT);
216 timeDistribution(output, input);
217 output[0] *= ((_numberOfParticles + _particleOffsetTime) / integral);
218 if (output[0] > 0) {
219 // adding up particle parts until a complete particle can be added
220 nSum += output[0];
221 if (nSum >= 1.) {
222 for (int i = 1; i <= nSum; ++i) {
223 // add the value of the time-step for one particle to _timeArray
224 _timeArray.push_back(iT);
225 }
226 nSum = 0;
227 }
228 }
229 }
230 // check if all calculated particles match total number of particles
231 if ((long long)(_timeArray.size()) == _numberOfParticles) {
232 return true;
233 }
234 else {
235 // calculate offset and calls function again with offset for convergence, choose
236 // different (larger) divisor if function does not converge
237 _particleOffsetTime += (_numberOfParticles - _timeArray.size())/7.;
238 calculateTimeArray(timeDistribution, converter, begin, end);
239 return false;
240 }
241 };
242
244 template <typename T, typename S>
246 OstreamManager clout(std::cout, "printTimeArray");
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;
250 }
251 };
252
253 /* this function can be called every time-step to add a new spherical particle when
254 * one should be added to the material number specified by the indicator.
255 * 'loopAround' only affects particle array after reaching last particle
256 * 'shuffleParticles=true' shuffles particle array before adding the first particle and after last
257 */
258 template <typename T, typename S>
259 template <typename C>
260 void ParticleDistribution<T, S>::spawnSphericalParticles(SuperParticleSystem3D<T, Particle3D>& supParticleSystem, IndicatorF3D<S>& indicator, T density, C iT, bool loopAround, bool shuffleParticles) {
261 while (_timeArray[_currentParticle] == iT) {
262 T p_radius = nextParticleRadius(loopAround, shuffleParticles);
263 if (p_radius > 0) {
264 supParticleSystem.addParticle(indicator, 4. / 3. * M_PI * util::pow(p_radius, 3) * density, p_radius, 1);
265 }
266 }
267 };
268
269 /* This function needs to be called at every time-step iT. It checks if according to the pre-calculated time-array a
270 * particle has to be added. It then activates the next particle in line.
271 */
272 template <typename T, typename S>
273 template <typename C>
275 // check if a new particle has to be activated
276 while (_timeArray[_currentParticle] == iT) {
277 // Searches all sub-particle systems
278 for (unsigned int pSystems = 0; pSystems < supParticleSystem.getParticleSystems().size(); pSystems++) {
279 // checks id of particle in sub-particle system
280 for (int particles = 0; particles < supParticleSystem.getParticleSystems()[pSystems]->size(); particles++) {
281 unsigned int currentPiD = supParticleSystem.getParticleSystems()[pSystems]->getParticlesPointer()[particles]->getID();
282 // checks id of particle in sub-particle system
283 if (currentPiD == _currentParticle + 1 + idOfset) {
284 // activates particles
285 supParticleSystem.getParticleSystems()[pSystems]->getParticlesPointer()[particles]->setActive(true);
286 }
287 }
288 }
289 _currentParticle++;
290 // starts at beginning when last particle is reached
291 if(_currentParticle == _numberOfParticles){
292 _currentParticle = 0;
293 }
294 }
295 };
296
298 template <typename T, typename S>
300 // Searches all sub-particle systems
301 for (int pSystems = 0; pSystems < supParticleSystem.numOfPSystems(); pSystems++) {
302 // searches all particles in sub-particle systems
303 for (int particles = 0; particles < supParticleSystem.getParticleSystems()[pSystems]->size(); particles++) {
304 // checks id of particle in sub-particle system
305 int currentPiD = supParticleSystem.getParticleSystems()[pSystems]->getParticlesPointer()[particles]->getID();
306 if (currentPiD <= endID && currentPiD >= beginID) {
307 // deactivates particles
308 supParticleSystem.getParticleSystems()[pSystems]->getParticlesPointer()[particles]->setActive(false);
309 }
310 }
311 }
312 };
313 /*This function adds all particles that are specified by the particle array to the ParticleSystem inside the indicator as spherical particles with
314 * their density. It sets their ID from 'beginID' to 'beginID+_numberOfParticles'. Use different IDs for every particle as they get activated by their
315 * ID. If shuffle=true, the particle array gets shuffled before adding.
316 */
317 template <typename T, typename S>
318 void ParticleDistribution<T, S>::preSpawnSphericalParticles(SuperParticleSystem3D<T, Particle3D>& supParticleSystem, IndicatorF3D<S>& indicator, T density, int beginID, bool shuffle) {
319 // call nextParticleRadius to get particle radius
320 T p_radius = nextParticleRadius(false, shuffle);
321 do {
322 // adds a spherical particle with specified density and '_currentParticle' as its ID
323 supParticleSystem.addTracerParticle(indicator, (T)(_currentParticle), 4. / 3. * M_PI * util::pow(p_radius, 3) * density, p_radius, 1);
324 p_radius = nextParticleRadius(false, shuffle);
325 } while (p_radius > 0);
326 _currentParticle = 0;
327 deactivateParticles(supParticleSystem, beginID, beginID + _numberOfParticles);
328 };
330 template <typename T, typename S>
332 for (long long i = 0; i < _numberOfParticleSizeGroups; i++) {
333 delete[] _sizeCalculation[i];
334 }
335 delete[] _sizeCalculation;
336 };
337};
338#endif
#define M_PI
Smart pointer for managing the various ways of passing functors around.
Definition functorPtr.h:60
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
Definition genericF.hh:51
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)
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)
Definition aDiff.h:475
cpu::simd::Pack< T > pow(cpu::simd::Pack< T > base, cpu::simd::Pack< T > exp)
Definition pack.h:112
ADf< T, DIM > exp(const ADf< T, DIM > &a)
Definition aDiff.h:455
Top level namespace for all of OpenLB.