OpenLB 1.7
Loading...
Searching...
No Matches
superParticleSystem3D.hh
Go to the documentation of this file.
1/* This file is part of the OpenLB library
2 *
3 * Copyright (C) 2016 Thomas Henn, Davide Dapelo
4 * E-mail contact: info@openlb.net
5 * The most recent release of OpenLB can be downloaded at
6 * <http://www.openlb.net/>
7 *
8 * This program is free software; you can redistribute it and/or
9 * modify it under the terms of the GNU General Public License
10 * as published by the Free Software Foundation; either version 2
11 * of the License, or (at your option) any later version.
12 *
13 * This program is distributed in the hope that it will be useful,
14 * but WITHOUT ANY WARRANTY; without even the implied warranty of
15 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
16 * GNU General Public License for more details.
17 *
18 * You should have received a copy of the GNU General Public
19 * License along with this program; if not, write to the Free
20 * Software Foundation, Inc., 51 Franklin Street, Fifth Floor,
21 * Boston, MA 02110-1301, USA.
22 */
23
24#ifndef SUPERPARTICLESYSTEM_3D_HH
25#define SUPERPARTICLESYSTEM_3D_HH
26
27#include <cassert>
28#define shadows
29
30#include <utility>
31#include <string>
32#include <array>
33#include <iomanip>
34#include <ios>
35#include <iostream>
36#include <random>
37#include "io/fileName.h"
38#include "superParticleSystem3D.h"
39
40#ifndef M_PI
41#define M_PI 3.14159265358979323846
42#endif
43
44namespace olb {
45
46 template<typename T, template<typename U> class PARTICLETYPE>
47 SuperParticleSystem3D<T, PARTICLETYPE>::SuperParticleSystem3D(
48 CuboidGeometry3D<T>& cuboidGeometry, LoadBalancer<T>& loadBalancer,
49 SuperGeometry<T,3>& superGeometry)
50 : SuperStructure<T,3>(cuboidGeometry, loadBalancer,
51 superGeometry.getOverlap()),
52 clout(std::cout, "SuperParticleSystem3d"),
53 _superGeometry(superGeometry),
54 _overlap(0)
55 {
56 init();
57 }
58
59 template<typename T, template<typename U> class PARTICLETYPE>
60 SuperParticleSystem3D<T, PARTICLETYPE>::SuperParticleSystem3D(
61 SuperGeometry<T,3>& superGeometry)
62 : SuperStructure<T,3>(superGeometry.getCuboidGeometry(),
63 superGeometry.getLoadBalancer(),
64 superGeometry.getOverlap()),
65 clout(std::cout, "SuperParticleSystem3d"),
66 _superGeometry(superGeometry),
67 _overlap(0)
68 {
69 init();
70 }
71
72 template<typename T, template<typename U> class PARTICLETYPE>
73 void SuperParticleSystem3D<T, PARTICLETYPE>::init()
74 {
75 int rank = 0;
76 int size = 0;
77
78 _stopSorting = 0;
79
80#ifdef PARALLEL_MODE_MPI
81 rank = singleton::mpi().getRank();
82 size = singleton::mpi().getSize();
83#endif
84 for (int i = 0; i < this->_cuboidGeometry.getNc(); ++i) {
85 if (this->_loadBalancer.rank(i) == rank) {
86 auto dummy = new ParticleSystem3D<T, PARTICLETYPE>(i, _superGeometry);
87 std::vector<T> physPos = toStdVector(this->_cuboidGeometry.get(i).getOrigin());
88 std::vector<T> physExtend(3, 0);
89 T physR = this->_cuboidGeometry.get(i).getDeltaR();
90 for (int j = 0; j < 3; j++) {
91 physPos[j] -= .5 * physR;
92 physExtend[j] = (this->_cuboidGeometry.get(i).getExtent()[j] + 1)
93 * physR;
94 }
95 dummy->setPosExt(physPos, physExtend);
96 _pSystems.push_back(dummy);
97 }
98 }
99
100#ifdef PARALLEL_MODE_MPI
101 for (int i = 0; i < this->_cuboidGeometry.getNc(); ++i) {
102 if (this->_loadBalancer.rank(i) == rank) {
103 std::vector<int> dummy;
104 this->getCuboidGeometry().getNeighbourhood(i, dummy, 3);
105 _rankNeighbours.insert(_rankNeighbours.end(), dummy.begin(), dummy.end());
106 _cuboidNeighbours.push_back(dummy);
107 }
108 }
109 for (auto& N : _rankNeighbours) {
110 N = this->_loadBalancer.rank(N);
111 }
112#endif
113
114 /* Ein jeder ist sein eigener Nachbar*/
115 if (rank == 0) {
116 _rankNeighbours.push_back(size - 1);
117 }
118 if (rank == size - 1) {
119 _rankNeighbours.push_back(0);
120 }
121 _rankNeighbours.push_back(rank);
122 _rankNeighbours.sort();
123 _rankNeighbours.unique();
124 }
125
126 template<typename T, template<typename U> class PARTICLETYPE>
127 SuperParticleSystem3D<T, PARTICLETYPE>::SuperParticleSystem3D(
128 SuperParticleSystem3D<T, PARTICLETYPE>& spSys)
129 : SuperStructure<T,3>(spSys._cuboidGeometry, spSys._loadBalancer,
130 int(spSys._overlap)),
131 clout(std::cout, "SuperParticleSystem3d"),
132 _pSystems(spSys._pSystems),
133 _superGeometry(spSys._superGeometry),
134 _rankNeighbours(spSys._rankNeighbours),
135 _overlap(spSys._overlap)
136 {
137 }
138
139 template<typename T, template<typename U> class PARTICLETYPE>
140 SuperParticleSystem3D<T, PARTICLETYPE>::SuperParticleSystem3D(
141 SuperParticleSystem3D<T, PARTICLETYPE> const& spSys)
142 : SuperStructure<T,3>(spSys._cuboidGeometry, spSys._loadBalancer,
143 int(spSys._overlap)),
144 clout(std::cout, "SuperParticleSystem3d"),
145 _pSystems(spSys._pSystems),
146 _superGeometry(spSys._superGeometry),
147 _rankNeighbours(spSys._rankNeighbours),
148 _overlap(spSys._overlap)
149 {
150 }
151
152 template<typename T, template<typename U> class PARTICLETYPE>
153 SuperParticleSystem3D<T, PARTICLETYPE>::SuperParticleSystem3D(
154 SuperParticleSystem3D<T, PARTICLETYPE> && spSys)
155 : SuperStructure<T,3>(spSys._cuboidGeometry, spSys._loadBalancer,
156 int(spSys._overlap)),
157 clout(std::cout, "SuperParticleSystem3d"),
158 _superGeometry(spSys._superGeometry),
159 _rankNeighbours(spSys._rankNeighbours),
160 _overlap(spSys._overlap)
161 {
162 _pSystems = std::move(spSys._pSystems);
163 }
164
165 template<typename T, template<typename U> class PARTICLETYPE>
166 void SuperParticleSystem3D<T, PARTICLETYPE>::print()
167 {
168 int no = globalNumOfParticles();
169 int active = globalNumOfActiveParticles();
170 clout << "activeParticles= " << active << " (" << no << ") " << std::endl;
171 // cout << "[SuperParticleSystem3D] " << _pSystems.size()
172 // << " pSystems on rank " << singleton::mpi().getRank() << "\n";
173 // for (auto pS : _pSystems) {
174 // cout << pS->_particles.size() << " ";
175 // }
176 }
177
178 template<typename T, template<typename U> class PARTICLETYPE>
179 void SuperParticleSystem3D<T, PARTICLETYPE>::printDeep(std::string message)
180 {
181 for (auto pS : _pSystems) {
182 pS->printDeep ( std::to_string(_pSystems.size()) + " pSystems on rank " + std::to_string(singleton::mpi().getRank()) + ": " );
183 }
184 }
185
186 template<typename T, template<typename U> class PARTICLETYPE>
187 void SuperParticleSystem3D<T, PARTICLETYPE>::print(std::list<int> mat)
188 {
189 std::list<int>::iterator _matIter;
190 int no = globalNumOfParticles();
191 clout << "globalNumOfParticles=" << no;
192 int active = globalNumOfActiveParticles();
193 clout << "; activeParticles=" << active;
194 for (_matIter = mat.begin(); _matIter != mat.end(); _matIter++) {
195 clout << "; material" << *_matIter << "=" << countMaterial(*_matIter);
196 }
197 clout << std::endl;
198 }
199
200 template<typename T, template<typename U> class PARTICLETYPE>
201 void SuperParticleSystem3D<T, PARTICLETYPE>::captureEscapeRate(
202 std::list<int> mat)
203 {
204 std::list<int>::iterator _matIter;
205 T sum = T();
206 // capture rate
207 for (_matIter = mat.begin(); _matIter != mat.end(); _matIter++) {
208 sum += (T) countMaterial(*_matIter);
209 }
210 clout << "captureRate=" << 1. - sum / globalNumOfParticles()
211 << "; escapeRate=" << sum / globalNumOfParticles() << std::endl;
212 }
213
214 template<typename T, template<typename U> class PARTICLETYPE>
215 void SuperParticleSystem3D<T, PARTICLETYPE>::diffEscapeRate(
216 std::list<int> mat, int& globalPSum, int& pSumOutlet, T& diffEscapeRate, T& maxDiffEscapeRate,
217 int iT, int iTConsole, T genPartPerTimeStep)
218 {
219 std::list<int>::iterator _matIter;
220 T pSumOutletNew = T();
221 T maxDiffEscapeRateNew = maxDiffEscapeRate;
222 T diffERtmp = T(); // temporal differencial escape rate
223
224 // count particles at outlet
225 for (_matIter = mat.begin(); _matIter != mat.end(); _matIter++) {
226 pSumOutletNew += (T) countMaterial(*_matIter);
227 }
228
229 // calculate diff. escape rate
230 if (globalNumOfParticles() > globalPSum) {
231 T diffERtmp = (T) (pSumOutletNew - pSumOutlet) / (globalNumOfParticles() - globalPSum);
232 diffEscapeRate += diffERtmp;
233 }
234 else {
235 if (genPartPerTimeStep != 0.) {
236 diffERtmp = (T) (pSumOutletNew - pSumOutlet) / (genPartPerTimeStep);
237 diffEscapeRate += diffERtmp;
238 }
239 }
240 // calculate max. diff. escape rate
241 if (diffERtmp > maxDiffEscapeRateNew) {
242 maxDiffEscapeRateNew = diffERtmp;
243 }
244 // console output
245 if (iT % iTConsole == 0) {
246 diffEscapeRate /= iTConsole;
247 if (globalNumOfParticles() > globalPSum) {
248 clout << "diffEscapeRate = " << diffEscapeRate << std::endl;
249 }
250 else {
251 if (genPartPerTimeStep != 0.) {
252 clout << "no particles in feedstream, continue calculation of diffEscapeRate with theoretical "
253 << genPartPerTimeStep << " generated particles per phys. time step"
254 << " diffEscapeRate = " << diffEscapeRate << std::endl;
255 }
256 else {
257 clout << "no particles in feedstream, calculation of diffEscapeRate not possible" << std::endl;
258 }
259 }
260 if (maxDiffEscapeRateNew > maxDiffEscapeRate) {
261 clout << "maxDiffEscapeRate = " << maxDiffEscapeRateNew << std::endl;
262 }
263 diffEscapeRate = 0. ;
264 }
265 maxDiffEscapeRate = maxDiffEscapeRateNew;
266 pSumOutlet = pSumOutletNew;
267 globalPSum = globalNumOfParticles();
268 }
269
270 template<typename T, template<typename U> class PARTICLETYPE>
271 void SuperParticleSystem3D<T, PARTICLETYPE>::diffEscapeRate(
272 std::list<int> mat, int& globalPSum, int& pSumOutlet, T& diffEscapeRate, T& maxDiffEscapeRate,
273 int iT, int iTConsole, T genPartPerTimeStep,
274 T& avDiffEscapeRate, T latticeTimeStart, T latticeTimeEnd)
275 {
276 std::list<int>::iterator _matIter;
277 T pSumOutletNew = T();
278 T maxDiffEscapeRateNew = maxDiffEscapeRate;
279 T avDiffEscapeRateNew = T();
280
281 // count particle at outlet
282 for (_matIter = mat.begin(); _matIter != mat.end(); _matIter++) {
283 pSumOutletNew += (T) countMaterial(*_matIter);
284 }
285 // calculate diff. escape rate
286 if (globalNumOfParticles() > globalPSum) {
287 avDiffEscapeRateNew = (T) (pSumOutletNew - pSumOutlet) / (globalNumOfParticles() - globalPSum);
288 diffEscapeRate += avDiffEscapeRateNew;
289 }
290 else {
291 if (genPartPerTimeStep != 0.) {
292 avDiffEscapeRateNew = (T) (pSumOutletNew - pSumOutlet) / (genPartPerTimeStep);
293 diffEscapeRate += avDiffEscapeRateNew;
294 }
295 }
296 // calculate max. diff. escape rate
297 if (avDiffEscapeRateNew > maxDiffEscapeRateNew) {
298 maxDiffEscapeRateNew = avDiffEscapeRateNew;
299 }
300 // calculate average diff. escape rate between tStart and tEnd
301 if ((iT >= latticeTimeStart) && (iT <= latticeTimeEnd)) {
302 avDiffEscapeRate += avDiffEscapeRateNew;
303 }
304 if (iT == latticeTimeEnd) {
305 avDiffEscapeRate /= (latticeTimeEnd - latticeTimeStart);
306 clout << "average diffEscapeRate between t = " << latticeTimeStart << "s and t = " << latticeTimeEnd << "s : "
307 << avDiffEscapeRate << std::endl;
308 }
309 // console output
310 if (iT % iTConsole == 0) {
311 diffEscapeRate /= iTConsole;
312 if (globalNumOfParticles() > globalPSum) {
313 clout << "diffEscapeRate = " << diffEscapeRate << std::endl;
314 }
315 else {
316 if (genPartPerTimeStep != 0.) {
317 clout << "no particles in feedstream, continue calculation of diffEscapeRate with theoretical "
318 << genPartPerTimeStep << " generated particles per phys. time step" << std::endl;
319 clout << "diffEscapeRate = " << diffEscapeRate << std::endl;
320 }
321 else {
322 clout << "no particles in feedstream, calculation of diffEscapeRate not possible" << std::endl;
323 }
324 }
325 if (maxDiffEscapeRateNew > maxDiffEscapeRate) {
326 clout << "maxDiffEscapeRate = " << maxDiffEscapeRateNew << std::endl;
327 }
328 diffEscapeRate = 0. ;
329 }
330 maxDiffEscapeRate = maxDiffEscapeRateNew;
331 pSumOutlet = pSumOutletNew;
332 globalPSum = globalNumOfParticles();
333 }
334
335 // Output_txt.file for tracer particles (particles that have an id != 0
336 /* _filename is filename of output txt-file
337 * _file is contact to output file
338 * */
339 template<typename T, template<typename U> class PARTICLETYPE>
340 void SuperParticleSystem3D<T, PARTICLETYPE>::getOutput(std::string _filename,
341 int it, T conversionFactorTime, unsigned short _properties )
342 {
343
344 enum particleProperties
345 : unsigned short {position = 1,
346 velocity = 2,
347 radius = 4,
348 mass = 8,
349 force = 16,
350 storeForce = 32,
351 };
352
355 std::vector < T > id, rad;
356 std::vector<std::array<T, 3>> pos, vel, forces, storeForces;
358 int k, j, numOfTracerParticles = globalNumOfTracerParticles();
359 int m = 0;
361 std::setprecision(9);
363 for (k = 0; k < numOfTracerParticles; k++) {
364 id.push_back(T());
365 rad.push_back(T());
366 pos.push_back( { T(), T(), T() });
367 vel.push_back( { T(), T(), T() });
368 forces.push_back( { T(), T(), T() });
369 storeForces.push_back( { T(), T(), T() });
370 }
371
372 k = 0;
374 std::vector<ParticleSystem3D<T, PARTICLETYPE>*> _pSystems = getPSystems();
375 for (unsigned int pS = 0; pS < _pSystems.size(); ++pS) {
376 std::deque<PARTICLETYPE<T>*> particles =
377 _pSystems[pS]->getParticlesPointer();
378 for (auto p : particles) {
379 if (p->getID() != 0) {
381 id[k] = p->getID();
382 rad[k] = p->getRad();
383 for (j = 0; j < 3; j++) {
384 pos[k][j] = p->getPos()[j];
385 vel[k][j] = p->getVel()[j];
386 forces[k][j] = p->getForce()[j];
387 storeForces[k][j] = p->getStoreForce()[j];
388 }
389 // p->resetStoreForce();
390 k++;
391 p++;
392 }
393 }
394 }
395
396#ifdef PARALLEL_MODE_MPI
398 for (m = 0; m < numOfTracerParticles; m++) {
399 singleton::mpi().reduceAndBcast(id[m], MPI_SUM);
400 singleton::mpi().reduceAndBcast(rad[m], MPI_SUM);
401 for (j = 0; j < 3; j++) {
402 singleton::mpi().reduceAndBcast(pos[m][j], MPI_SUM);
403 singleton::mpi().reduceAndBcast(vel[m][j], MPI_SUM);
404 singleton::mpi().reduceAndBcast(forces[m][j], MPI_SUM);
405 singleton::mpi().reduceAndBcast(storeForces[m][j], MPI_SUM);
406 }
407 }
408#endif
409
411 if (!singleton::mpi().getRank()) {
412 std::ofstream _file;
413
414 int i = 0;
415 if (it == 0) {
418 _file.open(_filename, std::ios::out | std::ios::trunc);
419 _file << "Timestep" << i + 1 << " " << "physTime" << i + 2;
420 i = i + 2;
421
422 for (m = 0; m < numOfTracerParticles; m++) {
423 _file << " id" << i + 1;
424 i = i + 1;
425 _file << " rad" << i + 1;
426 i = i + 1;
427 if (_properties & particleProperties::position) {
428 _file << " pos0_" << i + 1 << " pos1_" << " pos2_" << i + 3;
429 i = i + 3;
430 }
431 if (_properties & particleProperties::force) {
432 _file << " forc0_" << i + 1 << " forc1_" << i + 2 << " forc2_"
433 << i + 3;
434 i = i + 3;
435 }
436 if (_properties & particleProperties::velocity) {
437 _file << " vel0_" << i + 1 << " vel1_" << i + 2 << " vel2_" << i + 3;
438 i = i + 3;
439 }
440 if (_properties & particleProperties::storeForce) {
441 _file << " hforc0_" << i + 1 << " hforc1_" << i + 2 << " hforc2_"
442 << i + 3;
443 i = i + 3;
444 }
445 }
446 _file << "\n";
447 _file.close();
448 }
449
450 if (it >= 0) {
452 _file.open(_filename, std::ios::out | std::ios::app);
453 _file << it << " " << conversionFactorTime*it << " ";
454
455 for (m = 0; m < numOfTracerParticles; m++) {
456 _file << id[m] << " ";
457 _file << rad[m] << " ";
458 if (_properties & particleProperties::position) {
459 _file << pos[m][0] << " " << pos[m][1] << " " << pos[m][2] << " ";
460 }
461 if (_properties & particleProperties::force) {
462 _file << forces[m][0] << " " << forces[m][1] << " " << forces[m][2]
463 << " ";
464 }
465 if (_properties & particleProperties::velocity) {
466 _file << vel[m][0] << " " << vel[m][1] << " " << vel[m][2] << " ";
467 }
468 if (_properties & particleProperties::storeForce) {
469 _file << storeForces[m][0] << " " << storeForces[m][1] << " "
470 << storeForces[m][2] << " ";
471 }
472 }
473 }
474 _file << "\n";
475 _file.close();
476 }
477 }
478
479
480 template<typename T, template<typename U> class PARTICLETYPE>
481 std::vector<ParticleSystem3D<T, PARTICLETYPE>*>& SuperParticleSystem3D<T,
482 PARTICLETYPE>::getPSystems() //+*
483 {
484 return _pSystems;
485 }
486
487 template<typename T, template<typename U> class PARTICLETYPE>
488 void SuperParticleSystem3D<T, PARTICLETYPE>::simulate(T dT, bool scale)
489 {
490
491 // for (auto pS : _pSystems) {
492 // pS->velocityVerlet1(dT);
493 // if (_overlap > 0) {
494 // pS->makeTree();
495 // cout << "makeTree" << std::endl;
496 // }
497 // }
498 // updateParticleDistribution();
499 // for (auto pS : _pSystems) {
500 // pS->velocityVerlet2(dT);
501 // }
502
503 // updateParticleDistribution();
504 // for (auto pS : _pSystems) {
505 // pS->computeForce();
506 // pS->predictorCorrector1(dT);
507 // }
508 //
509 // updateParticleDistribution();
510 // for (auto pS : _pSystems) {
511 // pS->computeForce();
512 // pS->predictorCorrector2(dT);
513 // }
514
515 // for (auto pS : _pSystems) {
516 // pS->rungeKutta4_1(dT);
517 // }
518 // updateParticleDistribution();
519 // for (auto pS : _pSystems) {
520 // pS->rungeKutta4_2(dT);
521 // }
522 // updateParticleDistribution();
523 // for (auto pS : _pSystems) {
524 // pS->rungeKutta4_3(dT);
525 // }
526 // updateParticleDistribution();
527 // for (auto pS : _pSystems) {
528 // pS->rungeKutta4_4(dT);
529 // }
530 // updateParticleDistribution();
531 for (auto pS : _pSystems) {
532 time_t delta = clock();
533 pS->_contactDetection->sort();
534 _stopSorting += clock() - delta;
535 pS->simulate(dT, scale);
536 pS->computeBoundary();
537
538 }
539 updateParticleDistribution();
540 }
541
542 template<typename T, template<typename U> class PARTICLETYPE>
543 void SuperParticleSystem3D<T, PARTICLETYPE>::simulateWithTwoWayCoupling_Mathias ( T dT,
544 ForwardCouplingModel<T,PARTICLETYPE>& forwardCoupling,
545 BackCouplingModel<T,PARTICLETYPE>& backCoupling,
546 int material, int subSteps, bool resetExternalField, bool scale )
547 {
548 // reset external field
549 if (resetExternalField) {
550 backCoupling.resetExternalField(material);
551 }
552
553 for (auto pS : _pSystems) {
554 time_t delta = clock();
555 pS->_contactDetection->sort();
556 _stopSorting += clock() - delta;
557 pS->simulateWithTwoWayCoupling_Mathias(dT, forwardCoupling, backCoupling, material, subSteps, scale);
558 pS->computeBoundary();
559
560 }
561 updateParticleDistribution();
562 }
563
564 template<typename T, template<typename U> class PARTICLETYPE>
565 void SuperParticleSystem3D<T, PARTICLETYPE>::simulateWithTwoWayCoupling_Davide ( T dT,
566 ForwardCouplingModel<T,PARTICLETYPE>& forwardCoupling,
567 BackCouplingModel<T,PARTICLETYPE>& backCoupling,
568 int material, int subSteps, bool resetExternalField, bool scale )
569 {
570 for (auto pS : _pSystems) {
571 time_t delta = clock();
572 pS->_contactDetection->sort();
573 _stopSorting += clock() - delta;
574 pS->simulateWithTwoWayCoupling_Davide(dT, forwardCoupling, backCoupling, material, subSteps, scale);
575 pS->computeBoundary();
576 pS->eraseInactiveParticles();
577 }
578 updateParticleDistribution();
579 backCoupling.communicate();
580 }
581
582 // multiple collision models
583 template<>
584 void SuperParticleSystem3D<double, MagneticParticle3D>::simulate(double dT, std::set<int> sActivityOfParticle, bool scale)
585 {
586 for (auto pS : _pSystems) {
587 time_t delta = clock();
588 if (pS->getIGeometry() == singleton::mpi().getRank()) {
589 pS->_contactDetection->sort();
590 }
591 _stopSorting += clock() - delta;
592 if (pS->getIGeometry() == singleton::mpi().getRank()) {
593 pS->simulate(dT, sActivityOfParticle, scale);
594 pS->computeBoundary();
595 }
596 }
597 updateParticleDistribution();
598 }
599
600 template<>
601 bool SuperParticleSystem3D<double, MagneticParticle3D>::particleSActivityTest(int sActivity)
602 {
603 for (auto pS : _pSystems) {
604 for (auto p : pS->_particles) {
605 if (p.getSActivity() == sActivity) {
606 return false;
607 }
608 }
609 }
610 return true;
611 }
612
613 template<typename T, template<typename U> class PARTICLETYPE>
614 void SuperParticleSystem3D<T, PARTICLETYPE>::setOverlap(T overlap)
615 {
616 assert( int(overlap) + 1 <= _superGeometry.getOverlap() );
617 _overlap = overlap;
618 }
619
620 template<typename T, template<typename U> class PARTICLETYPE>
621 T SuperParticleSystem3D<T, PARTICLETYPE>::getOverlap()
622 {
623 return _overlap;
624 }
625
626 template<typename T, template<typename U> class PARTICLETYPE>
627 int SuperParticleSystem3D<T, PARTICLETYPE>::numOfPSystems()
628 {
629 return _pSystems.size();
630 }
631
632 template<typename T, template<typename U> class PARTICLETYPE>
633 int SuperParticleSystem3D<T, PARTICLETYPE>::globalNumOfParticles()
634 {
635#ifdef PARALLEL_MODE_MPI
636 // cout << "return1" << std::endl;
637 int buffer = rankNumOfParticles();
638 // cout << "return2" << std::endl;
639 singleton::mpi().reduceAndBcast(buffer, MPI_SUM);
640 // cout << "return3" << std::endl;
641 return buffer;
642#else
643 // cout << "return4" << std::endl;
644 return rankNumOfParticles();
645#endif
646 }
647
648 template<typename T, template<typename U> class PARTICLETYPE>
649 int SuperParticleSystem3D<T, PARTICLETYPE>::globalNumOfShadowParticles()
650 {
651#ifdef PARALLEL_MODE_MPI
652 int buffer = rankNumOfShadowParticles();
653 singleton::mpi().reduceAndBcast(buffer, MPI_SUM);
654 return buffer;
655#else
656 return rankNumOfShadowParticles();
657#endif
658 }
659
660 template<typename T, template<typename U> class PARTICLETYPE>
661 int SuperParticleSystem3D<T, PARTICLETYPE>::rankNumOfParticles()
662 {
663 int num = 0;
664 for (auto pS : _pSystems) {
665 num += pS->size();
666 }
667 return num;
668 }
669
670 template<typename T, template<typename U> class PARTICLETYPE>
671 int SuperParticleSystem3D<T, PARTICLETYPE>::rankNumOfShadowParticles()
672 {
673 int num = 0;
674 for (auto pS : _pSystems) {
675 num += pS->_shadowParticles.size();
676 }
677 return num;
678 }
679
680 template<typename T, template<typename U> class PARTICLETYPE>
681 int SuperParticleSystem3D<T, PARTICLETYPE>::rankNumOfActiveParticles()
682 {
683 int num = 0;
684 for (auto pS : _pSystems) {
685 num += pS->numOfActiveParticles();
686 }
687 return num;
688 }
689
690 template<typename T, template<typename U> class PARTICLETYPE>
691 int SuperParticleSystem3D<T, PARTICLETYPE>::countLocMaterial(int mat)
692 {
693 int num = 0;
694 for (auto pS : _pSystems) {
695 num += pS->countMaterial(mat);
696 }
697 return num;
698 }
699
700 template<typename T, template<typename U> class PARTICLETYPE>
701 int SuperParticleSystem3D<T, PARTICLETYPE>::countMaterial(int mat)
702 {
703#ifdef PARALLEL_MODE_MPI
704 int buffer = countLocMaterial(mat);
705 singleton::mpi().reduceAndBcast(buffer, MPI_SUM);
706 return buffer;
707#else
708 return countLocMaterial(mat);
709#endif
710 }
711
712 template<typename T, template<typename U> class PARTICLETYPE>
713 int SuperParticleSystem3D<T, PARTICLETYPE>::globalNumOfActiveParticles()
714 {
715#ifdef PARALLEL_MODE_MPI
716 int buffer = rankNumOfActiveParticles();
717 singleton::mpi().reduceAndBcast(buffer, MPI_SUM);
718 return buffer;
719#else
720 return rankNumOfActiveParticles();
721#endif
722 }
723
724 template<typename T, template<typename U> class PARTICLETYPE>
725 int SuperParticleSystem3D<T, PARTICLETYPE>::globalNumOfTracerParticles()
726 {
727#if PARALLEL_MODE_MPI
728 int buffer = rankNumOfTracerParticles();
729 singleton::mpi().reduceAndBcast(buffer, MPI_SUM);
730 return buffer;
731#else
732 return rankNumOfTracerParticles();
733#endif
734 }
735
736 template<typename T, template<typename U> class PARTICLETYPE>
737 int SuperParticleSystem3D<T, PARTICLETYPE>::rankNumOfTracerParticles()
738 {
739 int num = 0;
740 for (auto pS : _pSystems) {
741 std::deque<PARTICLETYPE<T>*> particles = pS->getParticlesPointer();
742 int pSNum = 0;
743 for (auto p : particles) {
744 if (p->getID() != 0) {
745 pSNum++;
746 }
747 }
748 num += pSNum;
749 }
750 return num;
751 }
752
753 // TODO class olb::ParticleSystem3D<double, olb::Particle3D>’ has no member named ‘radiusDistribution
754 /*
755 template<typename T, template<typename U> class PARTICLETYPE>
756 std::map<T, int> SuperParticleSystem3D<T, PARTICLETYPE>::rankRadiusDistribution()
757 {
758 std::map<T, int> distrAll;
759 typename std::map<T, int>::iterator ita;
760 for (auto pS : _pSystems) {
761 std::map<T, int> distrPs;
762 distrPs = pS->radiusDistribution();
763 for (typename std::map<T, int>::iterator itp = distrPs.begin();
764 itp != distrPs.end(); ++itp) {
765 T radKeyP = itp->first;
766 int countP = itp->second;
767 if (distrAll.count(radKeyP) > 0) {
768 ita = distrAll.find(radKeyP);
769 ita->second = ita->second + countP;
770 } else {
771 distrAll.insert(std::pair<T, int>(radKeyP, countP));
772 }
773 }
774 }
775 return distrAll;
776 }
777 */
778
779 template<typename T, template<typename U> class PARTICLETYPE>
780 std::vector<int> SuperParticleSystem3D<T, PARTICLETYPE>::numOfForces()
781 {
782 std::vector<int> dummy;
783 for (auto pS : _pSystems) {
784 dummy.push_back(pS->numOfForces());
785 }
786 return dummy;
787 }
788
789 template<typename T, template<typename U> class PARTICLETYPE>
790 template<typename DESCRIPTOR>
791 void SuperParticleSystem3D<T, PARTICLETYPE>::setVelToFluidVel(
792 SuperLatticeInterpPhysVelocity3D<T, DESCRIPTOR>& fVel)
793 {
794 for (auto pS : _pSystems) {
795 pS->setVelToFluidVel(fVel);
796 }
797 }
798
799 template<typename T, template<typename U> class PARTICLETYPE>
800 void SuperParticleSystem3D<T, PARTICLETYPE>::setVelToAnalyticalVel(
801 AnalyticalConst3D<T, T>& aVel)
802 {
803 for (auto pS : _pSystems) {
804 pS->setVelToAnalyticalVel(aVel);
805 }
806 }
807
808 template<typename T, template<typename U> class PARTICLETYPE>
809 bool SuperParticleSystem3D<T, PARTICLETYPE>::findCuboid(PARTICLETYPE<T> &p)
810 {
811 return findCuboid(p, int(_overlap));
812 }
813
814 template<typename T, template<typename U> class PARTICLETYPE>
815 bool SuperParticleSystem3D<T, PARTICLETYPE>::findCuboid(PARTICLETYPE<T> &p,
816 int overlap)
817 {
818 int C = this->_cuboidGeometry.get_iC(p.getPos()[0], p.getPos()[1],
819 p.getPos()[2], overlap);
820 if (C != this->_cuboidGeometry.getNc()) {
821 p.setCuboid(C);
822 return 1;
823 }
824 else {
825 clout << "Lost Particle! Pos: " << p.getPos()[0] << " " << p.getPos()[1]
826 << " " << p.getPos()[2] << " Vel: " << p.getVel()[0] << " "
827 << p.getVel()[1] << " " << p.getVel()[2] << " Cuboid: " << C << std::endl;
828 p.setActive(false);
829 return 0;
830 }
831 }
832
833 template<typename T, template<typename U> class PARTICLETYPE>
834 bool SuperParticleSystem3D<T, PARTICLETYPE>::checkCuboid(PARTICLETYPE<T>& p,
835 T overlap)
836 {
837 return checkCuboid(p, overlap, p.getCuboid());
838 }
839
840 template<typename T, template<typename U> class PARTICLETYPE>
841 bool SuperParticleSystem3D<T, PARTICLETYPE>::checkCuboid(PARTICLETYPE<T>& p,
842 T overlap, int iC)
843 {
844 // Checks whether particle is contained in the cuboid
845 // extended with an layer of size overlap*delta
846 return this->_cuboidGeometry.get(iC).physCheckPoint(p.getPos()[0],
847 p.getPos()[1],
848 p.getPos()[2], overlap);
849 }
850
851 template<typename T, template<typename U> class PARTICLETYPE>
852 void SuperParticleSystem3D<T, PARTICLETYPE>::updateParticleDistribution()
853 {
854 // Find particles on wrong cuboid, store in relocate and delete
855 // maps particles to their new rank
856 // std::multimap<int, PARTICLETYPE<T> > relocate;
857 _relocate.clear();
858#ifdef shadows
859 _relocateShadow.clear();
860#endif
861 for (unsigned int pS = 0; pS < _pSystems.size(); ++pS) {
862
863 _pSystems[pS]->_shadowParticles.clear();
864 auto par = _pSystems[pS]->_particles.begin();
865
866 while (par != _pSystems[pS]->_particles.end()) {
867
868 //Check if particle is still in its cuboid
869 if (checkCuboid(*par, 0)) {
870#ifdef shadows
871 // Check if its inside boundary layer of its cuboid
872 if (!(checkCuboid(*par, -_overlap))) {
873
874 std::set<int> sendTo;
875 // Run through all cuboids to search the one that
876 // shares its boundary with the cuboid containing the particle
877 for (int iC = 0; iC < this->_cuboidGeometry.getNc(); iC++) {
878 // Check if iC is neighbor cuboid in which overlap the particle is
879 // located
880 if (par->getCuboid() != iC && checkCuboid(*par, _overlap, iC)) {
881 // rank is the processing thread of the global cuboid number iC
882 int rank = this->_loadBalancer.rank(iC);
883 // insert rank into sendTo, if sendTo does not already contain
884 // rank for the actual particle
885 if (!sendTo.count(rank)) {
886 _relocateShadow.insert(std::make_pair(rank, (*par)));
887 sendTo.insert(rank);
888 }
889 }
890 }
891 }
892#endif
893 //If not --> find new cuboid
894 }
895 else {
896 // check to which cuboid particle is located
897 // and give new cuboid to par
898 findCuboid(*par, 0);
899 // gives particle the new cuboid number where it is now located to
900 _relocate.insert(
901 std::make_pair(this->_loadBalancer.rank(par->getCuboid()), (*par)));
902#ifdef shadows
903 // If the new particle position is in boundary layer
904 // If yes --> check if its inside boundary layer
905 if (!(checkCuboid(*par, -_overlap))) {
906 // yes, particle is in boundary layer
907 std::set<int> sendTo;
908 for (int iC = 0; iC < this->_cuboidGeometry.getNc(); iC++) {
909 // checks if iC is neighbor cuboid in whose overlap the particle
910 // is located
911 if (par->getCuboid() != iC && checkCuboid(*par, _overlap, iC)) {
912 int rank = this->_loadBalancer.rank(iC);
913 if (!sendTo.count(rank)) {
914 // shadow particle (copy of original) is inserted for iC
915 // but just with information of rank, not of cuboid!
916 _relocateShadow.insert(std::make_pair(rank, (*par)));
917 sendTo.insert(rank);
918 }
919 }
920 }
921 }
922#endif
923 par = _pSystems[pS]->_particles.erase(par);
924 par--;
925 }
926 par++;
927 }
928 }
929 /* Communicate number of Particles per cuboid*/
930#ifdef PARALLEL_MODE_MPI
931 singleton::MpiNonBlockingHelper mpiNbHelper;
932 //mpiNbHelper.allocate(_rankNeighbours.size());
933
934 int k = 0;
935
936 /* Serialize particles */
937 // it is: std::map<int, std::vector<double> > _send_buffer
938 _send_buffer.clear();
939
940 // fill std::map<int, std::vector<double> > _send_buffer with
941 // particle properties of particles in _relocate
942
943 // buffer size equals number of particle properties
944 T buffer[PARTICLETYPE<T>::serialPartSize];
945 // insert all regular particles of _relocate into _send_buffer
946 for (auto rN : _relocate) {
947 // second element in rN is (*par)
948 // write particle properties into buffer
949 rN.second.serialize(buffer);
950 // first element in rN is rank
951 // it is: std::map<int, std::vector<double> > _send_buffer;
952 // enlarge the _send_buffer element (rank), which is of vector type, with
953 // new elements of buffer (begin: buffer, end: buffer + serialPartSize)
954 _send_buffer[rN.first].insert(_send_buffer[rN.first].end(),
955 buffer, buffer + PARTICLETYPE<T>::serialPartSize);
956 }
957
958 /*Send Particles */
959 k = 0;
960 // noSends contains number of neighboring cuboids
961 int noSends = 0;
962 // it is: std::list<int> _rankNeighbours, rank of neighboring cuboids
963 for (auto rN : _rankNeighbours) {
964 // if there are particles contained in _send_buffer, increase noSends
965 if (_send_buffer[rN].size() > 0) {
966 ++noSends;
967 }
968 }
969 // int noSends = _send_buffer.size();
970 if (noSends > 0) {
971 mpiNbHelper.allocate(noSends);
972 // cout << mpiNbHelper.get_size() << std::endl;
973 for (auto rN : _rankNeighbours) {
974 if (_send_buffer[rN].size() > 0) {
975 singleton::mpi().iSend<double>(&_send_buffer[rN][0],
976 _relocate.count(rN)*PARTICLETYPE<T>::serialPartSize,
977 rN, &mpiNbHelper.get_mpiRequest()[k++], 1);
978 }
979 }
980 }
981
982 /*Receive and add particles*/
983 singleton::mpi().barrier();
984 k = 0;
985 int flag = 0;
986 MPI_Iprobe(MPI_ANY_SOURCE, 1, MPI_COMM_WORLD, &flag, MPI_STATUS_IGNORE);
987 if (flag) {
988 for (auto rN : _rankNeighbours) {
989 MPI_Status status;
990 int tmpFlag = 0;
991 MPI_Iprobe(rN, 1, MPI_COMM_WORLD, &tmpFlag, &status);
992 if (tmpFlag) {
993 int number_amount = 0;
994 MPI_Get_count(&status, MPI_DOUBLE, &number_amount);
995 //T recv_buffer[number_amount];
996 //singleton::mpi().receive(recv_buffer, number_amount, rN, 1);
997 std::vector<T> recv_buffer;
998 recv_buffer.reserve(number_amount);
999 singleton::mpi().receive(recv_buffer.data(), number_amount, rN, 1);
1000
1001 for (int iPar = 0; iPar < number_amount; iPar += PARTICLETYPE<T>::serialPartSize) {
1002 PARTICLETYPE<T> p;
1003 p.unserialize(&recv_buffer[iPar]);
1004 if (singleton::mpi().getRank() == this->_loadBalancer.rank(p.getCuboid())) {
1005 // particle added to pSystem with local cuboid number on actual thread
1006 _pSystems[this->_loadBalancer.loc(p.getCuboid())]->addParticle(p);
1007 }
1008 }
1009 }
1010 }
1011 }
1012 if (noSends > 0) {
1013 singleton::mpi().waitAll(mpiNbHelper);
1014 }
1015#ifdef shadows
1016 /**************************************************************************************************************/
1017 //Same Again for shadowParticles
1018 mpiNbHelper.allocate(_rankNeighbours.size());
1019 k = 0;
1020 /* Serialize particles */
1021 _send_buffer.clear();
1022 // T buffer[PARTICLETYPE<T>::serialPartSize];
1023 for (auto rN : _relocateShadow) {
1024 // std::vector<T> buffer = rN.second.serialize();
1025 rN.second.serialize(buffer);
1026 _send_buffer[rN.first].insert(_send_buffer[rN.first].end(),
1027 buffer, buffer + PARTICLETYPE<T>::serialPartSize);
1028 }
1029
1030 /*Send Particles */
1031 k = 0;
1032 noSends = _send_buffer.size();
1033 if (noSends > 0) {
1034 mpiNbHelper.allocate(noSends);
1035 for (auto rN : _rankNeighbours) {
1036 if (_send_buffer[rN].size() > 0) {
1037 singleton::mpi().iSend<double>(&_send_buffer[rN][0],
1038 _relocateShadow.count(rN)*PARTICLETYPE<T>::serialPartSize,
1039 rN, &mpiNbHelper.get_mpiRequest()[k++], 4);
1040 }
1041 }
1042 }
1043 /*Receive and add particles*/
1044 singleton::mpi().barrier();
1045 k = 0;
1046 flag = 0;
1047 MPI_Iprobe(MPI_ANY_SOURCE, 4, MPI_COMM_WORLD, &flag, MPI_STATUS_IGNORE);
1048 if (flag) {
1049 for (auto rN : _rankNeighbours) {
1050 MPI_Status status;
1051 int tmpFlag = 0;
1052 MPI_Iprobe(rN, 4, MPI_COMM_WORLD, &tmpFlag, &status);
1053 // cout << "Message from " << rN << " found on " << singleton::mpi().getRank() << std::endl;
1054 if (tmpFlag) {
1055 int number_amount = 0;
1056 MPI_Get_count(&status, MPI_DOUBLE, &number_amount);
1057 // cout << "Contains " << number_amount << " infos" << std::endl;
1058 std::vector<T> recv_buffer;
1059 recv_buffer.reserve(number_amount);
1060 singleton::mpi().receive(recv_buffer.data(), number_amount, rN, 4);
1061 for (int iPar = 0; iPar < number_amount; iPar += PARTICLETYPE<T>::serialPartSize) {
1062 // std::cout << "Particle unserialized" << std::endl;
1063 // par contains information of shadow particle
1064 // par is already shadow particle !!
1065 PARTICLETYPE<T> p;
1066 p.unserialize(&recv_buffer[iPar]);
1067 addShadowParticle(p);
1068 }
1069 }
1070 }
1071 }
1072 /* WARNING:
1073 * Performance warning: mpi barrier added - Asher 04/01/2017
1074 * MPI hanging on application /apps/asher/particle/simpleParticleExample
1075 * apparent communication error, tested with mpirun -np 6 when particle is
1076 * in overlap regions. This could be due to tag value of 4 coinciding with
1077 * tags used in fluid parallelization and conditional blocking receives.
1078 * Unaffected threads could continue on to collision step where these buffers
1079 * are altered/overwritten. When tested with tag of 4000007 no such error arises,
1080 * implying such a conflict.
1081 */
1082 singleton::mpi().barrier();
1083 if (noSends > 0) {
1084 singleton::mpi().waitAll(mpiNbHelper);
1085 }
1086#endif
1087 mpiNbHelper.free();
1088#else
1089 for (auto& par : _relocate) {
1090 // _relocate contains make_pair(rank, (*par))
1091 // therefore par.second is pointer to particle par
1092 addParticle(par.second);
1093 }
1094 for (auto& par : _relocateShadow) {
1095 // _relocateShadow contains make_pair(rank, (*par))
1096 // therefore par.second is pointer to particle par
1097 addShadowParticle(par.second);
1098 }
1099#endif
1100 }
1101
1102 template<typename T, template<typename U> class PARTICLETYPE>
1103 void SuperParticleSystem3D<T, PARTICLETYPE>::addParticle(PARTICLETYPE<T>& p)
1104 {
1105 if (findCuboid(p)) {
1106#ifdef PARALLEL_MODE_MPI
1107 if (singleton::mpi().getRank() == this->_loadBalancer.rank(p.getCuboid())) {
1108 _pSystems[this->_loadBalancer.loc(p.getCuboid())]->addParticle(p);
1109 }
1110 else {
1111 // clout << "Particle not found on Cuboid: " << p.getCuboid() << std::endl;
1112 // clout << "Ppos: " << p.getPos()[0] << " " << p.getPos()[1] << " " << p.getPos()[2]<< std::endl;
1113 }
1114#else
1115 _pSystems[this->_loadBalancer.loc(p.getCuboid())]->addParticle(p);
1116#endif
1117 }
1118 }
1119
1120 template<typename T, template<typename U> class PARTICLETYPE>
1121 void SuperParticleSystem3D<T, PARTICLETYPE>::addParticle(
1122 IndicatorF3D<T>& ind, T mas, T rad, int no, std::vector<T> vel)
1123 {
1124 srand(clock());
1125 // srand(rand());
1126 std::vector<T> pos(3, 0.);
1127 bool indic[1] = { false };
1128
1129 no += globalNumOfParticles();
1130 while (globalNumOfParticles() < no) {
1131 pos[0] = ind.getMin()[0]
1132 + (T) (rand() % 100000) / 100000. * (ind.getMax()[0] - ind.getMin()[0]);
1133 pos[1] = ind.getMin()[1]
1134 + (T) (rand() % 100000) / 100000. * (ind.getMax()[1] - ind.getMin()[1]);
1135 pos[2] = ind.getMin()[2]
1136 + (T) (rand() % 100000) / 100000. * (ind.getMax()[2] - ind.getMin()[2]);
1137
1138#ifdef PARALLEL_MODE_MPI
1139 singleton::mpi().bCast(&*pos.begin(), 3);
1140#endif
1141
1142 int x0, y0, z0, C;
1143 std::vector<int> locLat(4, 0);
1144 if (this->_cuboidGeometry.getFloorLatticeR(pos, locLat)) {
1145 C = locLat[0];
1146 if (this->_loadBalancer.rank(C) == singleton::mpi().getRank()) {
1147 x0 = locLat[1];
1148 y0 = locLat[2];
1149 z0 = locLat[3];
1150 if (_superGeometry.get(C, x0, y0, z0) == 1
1151 && _superGeometry.get(C, x0, y0 + 1, z0) == 1
1152 && _superGeometry.get(C, x0, y0, z0 + 1) == 1
1153 && _superGeometry.get(C, x0, y0 + 1, z0 + 1) == 1
1154 && _superGeometry.get(C, x0 + 1, y0, z0) == 1
1155 && _superGeometry.get(C, x0 + 1, y0 + 1, z0) == 1
1156 && _superGeometry.get(C, x0 + 1, y0, z0 + 1) == 1
1157 && _superGeometry.get(C, x0 + 1, y0 + 1, z0 + 1) == 1
1158 && ind(indic, &pos[0])) {
1159 PARTICLETYPE<T> p(pos, vel, mas, rad);
1160 addParticle(p);
1161 }
1162 }
1163 }
1164 }
1165 }
1166
1167 template<typename T, template<typename U> class PARTICLETYPE>
1168 void SuperParticleSystem3D<T, PARTICLETYPE>::addParticle(std::set<int>
1169 material, int no, T mas,
1170 T rad, std::vector<T> vel)
1171 {
1172 srand(time(nullptr));
1173 std::vector<T> pos(3, 0.), min(3, std::numeric_limits<T>::max()),
1174 max(3, std::numeric_limits<T>::min());
1175 olb::Vector<T, 3> tmpMin(0.), tmpMax(0.);
1176 std::set<int>::iterator it = material.begin();
1177 for (; it != material.end(); ++it) {
1178 tmpMin = _superGeometry.getStatistics().getMinPhysR(*it);
1179 tmpMax = _superGeometry.getStatistics().getMaxPhysR(*it);
1180 max[0] = util::max(max[0], tmpMax[0]);
1181 max[1] = util::max(max[1], tmpMax[1]);
1182 max[2] = util::max(max[2], tmpMax[2]);
1183 min[0] = util::min(min[0], tmpMin[0]);
1184 min[1] = util::min(min[1], tmpMin[1]);
1185 min[2] = util::min(min[2], tmpMin[2]);
1186 }
1187 // cout << "Min: " << min[0] << " " << min[1] << " " << min[2] << std::endl;
1188 // cout << "Max: " << max[0] << " " << max[1] << " " << max[2] << std::endl;
1189 for (int i = 0; i < 3; i++) {
1190 min[i] -= this->_cuboidGeometry.get(0).getDeltaR();
1191 max[i] += 2 * this->_cuboidGeometry.get(0).getDeltaR();
1192 }
1193 no += globalNumOfParticles();
1194 while (globalNumOfParticles() < no) {
1195 pos[0] = min[0] + (T) (rand() % 100000) / 100000. * (max[0] - min[0]);
1196 pos[1] = min[1] + (T) (rand() % 100000) / 100000. * (max[1] - min[1]);
1197 pos[2] = min[2] + (T) (rand() % 100000) / 100000. * (max[2] - min[2]);
1198
1199#ifdef PARALLEL_MODE_MPI
1200 singleton::mpi().bCast(&*pos.begin(), 3);
1201#endif
1202
1203 std::vector<int> locLat(4, 0);
1204 if (this->_cuboidGeometry.getFloorLatticeR(pos, locLat)) {
1205 if (this->_loadBalancer.rank(locLat[0]) == singleton::mpi().getRank()) {
1206 if (material.find(_superGeometry.get(locLat[0], locLat[1], locLat[2], locLat[3]))
1207 != material.end()
1208 && material.find(_superGeometry.get(locLat[0], locLat[1], locLat[2] + 1,
1209 locLat[3])) != material.end()
1210 && material.find(_superGeometry.get(locLat[0], locLat[1], locLat[2],
1211 locLat[3] + 1)) != material.end()
1212 && material.find(_superGeometry.get(locLat[0], locLat[1], locLat[2] + 1,
1213 locLat[3] + 1)) != material.end()
1214 && material.find(_superGeometry.get(locLat[0], locLat[1] + 1, locLat[2],
1215 locLat[3])) != material.end()
1216 && material.find(_superGeometry.get(locLat[0], locLat[1] + 1, locLat[2] + 1,
1217 locLat[3])) != material.end()
1218 && material.find(_superGeometry.get(locLat[0], locLat[1] + 1, locLat[2],
1219 locLat[3] + 1)) != material.end()
1220 && material.find(_superGeometry.get(locLat[0], locLat[1] + 1, locLat[2] + 1,
1221 locLat[3] + 1)) != material.end()) {
1222 PARTICLETYPE<T> p(pos, vel, mas, rad);
1223 addParticle(p);
1224 }
1225 }
1226 }
1227 }
1228 singleton::mpi().barrier();
1229 }
1230
1231 template<typename T, template<typename U> class PARTICLETYPE>
1232 void SuperParticleSystem3D<T, PARTICLETYPE>::addParticle(
1233 IndicatorF3D<T>& ind, std::set<int> material, T mas, T rad, int no, std::vector<T> vel)
1234 {
1235 srand(clock());
1236 std::vector<T> pos(3, 0.);
1237 bool indic[1] = { false };
1238
1239 no += globalNumOfParticles();
1240 while (globalNumOfParticles() < no) {
1241 pos[0] = ind.getMin()[0]
1242 + (T) (rand() % 100000) / 100000. * (ind.getMax()[0] - ind.getMin()[0]);
1243 pos[1] = ind.getMin()[1]
1244 + (T) (rand() % 100000) / 100000. * (ind.getMax()[1] - ind.getMin()[1]);
1245 pos[2] = ind.getMin()[2]
1246 + (T) (rand() % 100000) / 100000. * (ind.getMax()[2] - ind.getMin()[2]);
1247
1248#ifdef PARALLEL_MODE_MPI
1249 singleton::mpi().bCast(&*pos.begin(), 3);
1250#endif
1251
1252 int x0, y0, z0;
1253 std::vector<int> locLat(4, 0);
1254 if (this->_cuboidGeometry.getFloorLatticeR(pos, locLat)) {
1255 if (this->_loadBalancer.rank(locLat[0]) == singleton::mpi().getRank()) {
1256 x0 = locLat[1];
1257 y0 = locLat[2];
1258 z0 = locLat[3];
1259 if (_superGeometry.get(locLat[0], x0, y0, z0) == 1
1260 && _superGeometry.get(locLat[0], x0, y0 + 1, z0) == 1
1261 && _superGeometry.get(locLat[0], x0, y0, z0 + 1) == 1
1262 && _superGeometry.get(locLat[0], x0, y0 + 1, z0 + 1) == 1
1263 && _superGeometry.get(locLat[0], x0 + 1, y0, z0) == 1
1264 && _superGeometry.get(locLat[0], x0 + 1, y0 + 1, z0) == 1
1265 && _superGeometry.get(locLat[0], x0 + 1, y0, z0 + 1) == 1
1266 && _superGeometry.get(locLat[0], x0 + 1, y0 + 1, z0 + 1) == 1
1267 && ind(indic, &pos[0])) {
1268 if (material.find(
1269 _superGeometry.get(locLat[0], locLat[1], locLat[2], locLat[3]))
1270 != material.end()
1271 && material.find(
1272 _superGeometry.get(locLat[0], locLat[1], locLat[2] + 1,
1273 locLat[3])) != material.end()
1274 && material.find(
1275 _superGeometry.get(locLat[0], locLat[1], locLat[2],
1276 locLat[3] + 1)) != material.end()
1277 && material.find(
1278 _superGeometry.get(locLat[0], locLat[1], locLat[2] + 1,
1279 locLat[3] + 1)) != material.end()
1280 && material.find(
1281 _superGeometry.get(locLat[0], locLat[1] + 1, locLat[2],
1282 locLat[3])) != material.end()
1283 && material.find(
1284 _superGeometry.get(locLat[0], locLat[1] + 1, locLat[2] + 1,
1285 locLat[3])) != material.end()
1286 && material.find(
1287 _superGeometry.get(locLat[0], locLat[1] + 1, locLat[2],
1288 locLat[3] + 1)) != material.end()
1289 && material.find(
1290 _superGeometry.get(locLat[0], locLat[1] + 1, locLat[2] + 1,
1291 locLat[3] + 1)) != material.end()) {
1292 PARTICLETYPE<T> p(pos, vel, mas, rad);
1293 addParticle(p);
1294 }
1295 }
1296 }
1297 }
1298 }
1299 }
1300
1301 template<typename T, template<typename U> class PARTICLETYPE>
1302 void SuperParticleSystem3D<T, PARTICLETYPE>::addHL3DParticle(
1303 IndicatorF3D<T>& ind, std::set<int> material, T mas, T rad, int no, std::vector<T> vel, T surface, T volume)
1304 {
1305 std::cout << "addHL3DParticle begin" <<std::endl;
1306 srand(clock());
1307 std::vector<T> pos(3, 0.);
1308 bool indic[1] = { false };
1309
1310 no += globalNumOfParticles();
1311 while (globalNumOfParticles() < no) {
1312 pos[0] = ind.getMin()[0]
1313 + (T) (rand() % 100000) / 100000. * (ind.getMax()[0] - ind.getMin()[0]);
1314 pos[1] = ind.getMin()[1]
1315 + (T) (rand() % 100000) / 100000. * (ind.getMax()[1] - ind.getMin()[1]);
1316 pos[2] = ind.getMin()[2]
1317 + (T) (rand() % 100000) / 100000. * (ind.getMax()[2] - ind.getMin()[2]);
1318
1319#ifdef PARALLEL_MODE_MPI
1320 singleton::mpi().bCast(&*pos.begin(), 3);
1321#endif
1322
1323 int x0, y0, z0;
1324 std::vector<int> locLat(4, 0);
1325 if (this->_cuboidGeometry.getFloorLatticeR(pos, locLat)) {
1326 if (this->_loadBalancer.rank(locLat[0]) == singleton::mpi().getRank()) {
1327 x0 = locLat[1];
1328 y0 = locLat[2];
1329 z0 = locLat[3];
1330 if (_superGeometry.get(locLat[0], x0, y0, z0) == 1
1331 && _superGeometry.get(locLat[0], x0, y0 + 1, z0) == 1
1332 && _superGeometry.get(locLat[0], x0, y0, z0 + 1) == 1
1333 && _superGeometry.get(locLat[0], x0, y0 + 1, z0 + 1) == 1
1334 && _superGeometry.get(locLat[0], x0 + 1, y0, z0) == 1
1335 && _superGeometry.get(locLat[0], x0 + 1, y0 + 1, z0) == 1
1336 && _superGeometry.get(locLat[0], x0 + 1, y0, z0 + 1) == 1
1337 && _superGeometry.get(locLat[0], x0 + 1, y0 + 1, z0 + 1) == 1
1338 && ind(indic, &pos[0])) {
1339 if (material.find(
1340 _superGeometry.get(locLat[0], locLat[1], locLat[2], locLat[3]))
1341 != material.end()
1342 && material.find(
1343 _superGeometry.get(locLat[0], locLat[1], locLat[2] + 1,
1344 locLat[3])) != material.end()
1345 && material.find(
1346 _superGeometry.get(locLat[0], locLat[1], locLat[2],
1347 locLat[3] + 1)) != material.end()
1348 && material.find(
1349 _superGeometry.get(locLat[0], locLat[1], locLat[2] + 1,
1350 locLat[3] + 1)) != material.end()
1351 && material.find(
1352 _superGeometry.get(locLat[0], locLat[1] + 1, locLat[2],
1353 locLat[3])) != material.end()
1354 && material.find(
1355 _superGeometry.get(locLat[0], locLat[1] + 1, locLat[2] + 1,
1356 locLat[3])) != material.end()
1357 && material.find(
1358 _superGeometry.get(locLat[0], locLat[1] + 1, locLat[2],
1359 locLat[3] + 1)) != material.end()
1360 && material.find(
1361 _superGeometry.get(locLat[0], locLat[1] + 1, locLat[2] + 1,
1362 locLat[3] + 1)) != material.end()) {
1363 //PARTICLETYPE<T> p(pos, vel, mas, rad);
1364 PARTICLETYPE<T> p( pos, mas, rad, volume, surface);
1365 addParticle(p);
1366
1367 }
1368 }
1369 }
1370 }
1371 }
1372 }
1373
1374 template<>
1375 void SuperParticleSystem3D<double, MagneticParticle3D>::addParticle(
1376 IndicatorF3D<double>& ind, double mas, double rad, int no, int id,
1377 std::vector<double> vel, std::vector<double> dMoment, std::vector<double> aVel, std::vector<double> torque, double magnetisation,
1378 int sActivity)
1379 {
1380 std::vector<double> pos(3, 0.);
1381 bool indic[1] = { false };
1382
1383 no += globalNumOfParticles();
1384 while (globalNumOfParticles() < no) {
1385 pos[0] = ind.getMin()[0]
1386 + (double) (rand() % 100000) / 100000. * (ind.getMax()[0] - ind.getMin()[0]);
1387 pos[1] = ind.getMin()[1]
1388 + (double) (rand() % 100000) / 100000. * (ind.getMax()[1] - ind.getMin()[1]);
1389 pos[2] = ind.getMin()[2]
1390 + (double) (rand() % 100000) / 100000. * (ind.getMax()[2] - ind.getMin()[2]);
1391
1392#ifdef PARALLEL_MODE_MPI
1393 singleton::mpi().bCast(&*pos.begin(), 3);
1394#endif
1395
1396 int x0, y0, z0, C;
1397 std::vector<int> locLat(4, 0);
1398 if (this->_cuboidGeometry.getFloorLatticeR(pos, locLat)) {
1399 C = locLat[0];
1400 if (this->_loadBalancer.rank(C) == singleton::mpi().getRank()) {
1401 x0 = locLat[1];
1402 y0 = locLat[2];
1403 z0 = locLat[3];
1404 if (_superGeometry.get(C, x0, y0, z0) == 1
1405 && _superGeometry.get(C, x0, y0 + 1, z0) == 1
1406 && _superGeometry.get(C, x0, y0, z0 + 1) == 1
1407 && _superGeometry.get(C, x0, y0 + 1, z0 + 1) == 1
1408 && _superGeometry.get(C, x0 + 1, y0, z0) == 1
1409 && _superGeometry.get(C, x0 + 1, y0 + 1, z0) == 1
1410 && _superGeometry.get(C, x0 + 1, y0, z0 + 1) == 1
1411 && _superGeometry.get(C, x0 + 1, y0 + 1, z0 + 1) == 1
1412 && ind(indic, &pos[0])) {
1413 MagneticParticle3D<double> p(pos, vel, mas, rad, id, dMoment, aVel, torque, magnetisation, sActivity);
1414 id++;
1415 addParticle(p);
1416 }
1417 }
1418 }
1419 }
1420 }
1421
1422 template<>
1423 void SuperParticleSystem3D<double, MagneticParticle3D>::addParticle(IndicatorF3D<double>& ind, std::set<int> material, double mas, double rad, int no, int id,
1424 std::vector<double> vel, std::vector<double> dMoment, std::vector<double> aVel, std::vector<double> torque, double magnetisation,
1425 int sActivity)
1426
1427 {
1428 std::vector<double> pos(3, 0.);
1429 bool indic[1] = { false };
1430
1431 no += globalNumOfParticles();
1432 while (globalNumOfParticles() < no) {
1433 pos[0] = ind.getMin()[0]
1434 + (double) (rand() % 100000) / 100000. * (ind.getMax()[0] - ind.getMin()[0]);
1435 pos[1] = ind.getMin()[1]
1436 + (double) (rand() % 100000) / 100000. * (ind.getMax()[1] - ind.getMin()[1]);
1437 pos[2] = ind.getMin()[2]
1438 + (double) (rand() % 100000) / 100000. * (ind.getMax()[2] - ind.getMin()[2]);
1439
1440#ifdef PARALLEL_MODE_MPI
1441 singleton::mpi().bCast(&*pos.begin(), 3);
1442#endif
1443
1444 int x0, y0, z0;
1445 std::vector<int> locLat(4, 0);
1446 if (this->_cuboidGeometry.getFloorLatticeR(pos, locLat)) {
1447 if (this->_loadBalancer.rank(locLat[0]) == singleton::mpi().getRank()) {
1448 x0 = locLat[1];
1449 y0 = locLat[2];
1450 z0 = locLat[3];
1451 if (_superGeometry.get(locLat[0], x0, y0, z0) == 1
1452 && _superGeometry.get(locLat[0], x0, y0 + 1, z0) == 1
1453 && _superGeometry.get(locLat[0], x0, y0, z0 + 1) == 1
1454 && _superGeometry.get(locLat[0], x0, y0 + 1, z0 + 1) == 1
1455 && _superGeometry.get(locLat[0], x0 + 1, y0, z0) == 1
1456 && _superGeometry.get(locLat[0], x0 + 1, y0 + 1, z0) == 1
1457 && _superGeometry.get(locLat[0], x0 + 1, y0, z0 + 1) == 1
1458 && _superGeometry.get(locLat[0], x0 + 1, y0 + 1, z0 + 1) == 1
1459 && ind(indic, &pos[0])) {
1460 if (material.find(
1461 _superGeometry.get(locLat[0], locLat[1], locLat[2], locLat[3]))
1462 != material.end()
1463 && material.find(
1464 _superGeometry.get(locLat[0], locLat[1], locLat[2] + 1,
1465 locLat[3])) != material.end()
1466 && material.find(
1467 _superGeometry.get(locLat[0], locLat[1], locLat[2],
1468 locLat[3] + 1)) != material.end()
1469 && material.find(
1470 _superGeometry.get(locLat[0], locLat[1], locLat[2] + 1,
1471 locLat[3] + 1)) != material.end()
1472 && material.find(
1473 _superGeometry.get(locLat[0], locLat[1] + 1, locLat[2],
1474 locLat[3])) != material.end()
1475 && material.find(
1476 _superGeometry.get(locLat[0], locLat[1] + 1, locLat[2] + 1,
1477 locLat[3])) != material.end()
1478 && material.find(
1479 _superGeometry.get(locLat[0], locLat[1] + 1, locLat[2],
1480 locLat[3] + 1)) != material.end()
1481 && material.find(
1482 _superGeometry.get(locLat[0], locLat[1] + 1, locLat[2] + 1,
1483 locLat[3] + 1)) != material.end()) {
1484
1485 MagneticParticle3D<double> p(pos, vel, mas, rad, id, dMoment, aVel, torque, magnetisation, sActivity);
1486 id++;
1487 addParticle(p);
1488 }
1489 }
1490 }
1491 }
1492 }
1493 }
1494
1495 template<typename T, template<typename U> class PARTICLETYPE>
1496 template<typename DESCRIPTOR>
1497 void SuperParticleSystem3D<T, PARTICLETYPE>::generateParticlesCircleInletMassConcentration(
1498 IndicatorCircle3D<T>& indicatorCircle, T particleMassConcentration, T charPhysVelocity,
1499 T conversionFactorTime, SuperLatticeInterpPhysVelocity3D<T, DESCRIPTOR>& getVel,
1500 PARTICLETYPE<T>& p, std::set<int> material, int iT, T& particlesPerPhyTimeStep,
1501 std::vector<T>& inletVec, std::deque<std::vector<T>>& posDeq, int deqSize)
1502 {
1503
1504 std::vector<T> pos(3, 0.);
1505 std::vector<T> vel(3, 0.);
1506
1507 PARTICLETYPE<T> pCopy(p);
1508
1509 // r can be initialized with arbitrary values
1510 // r is calculated to lie in the same plane as indicatorCircle
1511 Vector<T, 3> r(inletVec);
1512 // s is calculated to lie in the same plane as indicatorCircle
1513 // s is orthogonal to r
1514 Vector<T, 3> s(0., 0., 0.);
1515
1516 T fluidVel[3] = {0., 0., 0.};
1517
1518 // calculation of r in the first step of iteration
1519 if (iT == 0) {
1520 T fluxDensity = charPhysVelocity * M_PI * util::pow(indicatorCircle.getRadius(), 2.) ;
1521 T particleVolume = 4. / 3.* M_PI * util::pow(p.getRad(), 3);
1522 T particleDensity = p.getMass() / particleVolume;
1523 T particleVolumeConcentration = particleMassConcentration / particleDensity;
1524 T particleVolumeFlux = fluxDensity * particleVolumeConcentration;
1525 T particleFlux = particleVolumeFlux / particleVolume;
1526 particlesPerPhyTimeStep = particleFlux * conversionFactorTime;
1527
1528 bool b = false;
1529 for (int i = 0; i < 3; i++) {
1530 if (indicatorCircle.getNormal()[i] == 0.) {
1531 b = true;
1532 }
1533 }
1534 if (b == true) {
1535 for (int i = 0; i < 3; i++) {
1536 if (indicatorCircle.getNormal()[i] == 0.) {
1537 r[i] = 1.;
1538 }
1539 else {
1540 r[i] = 0.;
1541 }
1542 }
1543 }
1544 else {
1545 r[0] = -(indicatorCircle.getNormal()[1] + indicatorCircle.getNormal()[2]) / indicatorCircle.getNormal()[0] ;
1546 r[1] = 1.;
1547 r[2] = 1.;
1548 }
1549 normalize(r) ;
1550 for (int i = 0; i <= 2; i++) {
1551 inletVec[i] = r[i];
1552 }
1553 }
1554
1555 // norm of r
1556 T r_max = indicatorCircle.getRadius();
1557 T r_min = -1. * r_max;
1558
1559 s = crossProduct3D(r, indicatorCircle.getNormal());
1560
1561 // Non-deterministic random number generator
1562 std::random_device rd;
1563 // Pseudo-random number engine: Mersenne Twister 19937 generator
1564 std::mt19937 engine(rd());
1565 int id = this->globalNumOfParticles();
1566
1567 while (this->globalNumOfParticles() < (iT + 1) * particlesPerPhyTimeStep) {
1568
1569 gt_mark:
1570 normalize(r);
1571 normalize(s);
1572
1573 // Random number distribution that produces floating-point values according to a uniform distribution
1574 std::uniform_real_distribution<T> distR(r_min, r_max);
1575
1576 // r_norm is between r_min and r_max
1577 T r_norm = distR(engine);
1578 r *= r_norm ;
1579
1580 T s_max = util::sqrt(util::pow(indicatorCircle.getRadius(), 2.) - util::pow(r_norm, 2.)) ;
1581 T s_min = -1. * s_max ;
1582 std::uniform_real_distribution<T> distS(s_min, s_max);
1583
1584 // s_norm is between s_min and s_max
1585 T s_norm = distS(engine);
1586 s *= s_norm ;
1587
1588 std::vector<T> posVecTmp(3, 0.);
1589 posVecTmp[0] = indicatorCircle.getCenter()[0] + r[0] + s[0];
1590 posVecTmp[1] = indicatorCircle.getCenter()[1] + r[1] + s[1];
1591 posVecTmp[2] = indicatorCircle.getCenter()[2] + r[2] + s[2];
1592
1593 for (auto a : posDeq) {
1594 T dist = util::sqrt(util::pow(a[0] - posVecTmp[0], 2.)
1595 + util::pow(a[1] - posVecTmp[1], 2.)
1596 + util::pow(a[2] - posVecTmp[2], 2.));
1597 if (dist <= 3 * p.getRad()) {
1598 goto gt_mark;
1599 }
1600 }
1601
1602 pos[0] = posVecTmp[0];
1603 pos[1] = posVecTmp[1];
1604 pos[2] = posVecTmp[2];
1605
1606 std::vector<int> latticeRoundedPos(4, 0);
1607
1608 if (this->_cuboidGeometry.getFloorLatticeR(pos, latticeRoundedPos)) {
1609 int globCuboid = latticeRoundedPos[0]; // is global cuboid number
1610 if (this->_loadBalancer.rank(globCuboid) == singleton::mpi().getRank()) {
1611
1612 if (material.find(_superGeometry.get(latticeRoundedPos[0], latticeRoundedPos[1], latticeRoundedPos[2], latticeRoundedPos[3]))
1613 != material.end()
1614 && material.find(_superGeometry.get(latticeRoundedPos[0], latticeRoundedPos[1], latticeRoundedPos[2] + 1,
1615 latticeRoundedPos[3])) != material.end()
1616 && material.find(_superGeometry.get(latticeRoundedPos[0], latticeRoundedPos[1], latticeRoundedPos[2],
1617 latticeRoundedPos[3] + 1)) != material.end()
1618 && material.find(_superGeometry.get(latticeRoundedPos[0], latticeRoundedPos[1], latticeRoundedPos[2] + 1,
1619 latticeRoundedPos[3] + 1)) != material.end()
1620 && material.find(_superGeometry.get(latticeRoundedPos[0], latticeRoundedPos[1] + 1, latticeRoundedPos[2],
1621 latticeRoundedPos[3])) != material.end()
1622 && material.find(_superGeometry.get(latticeRoundedPos[0], latticeRoundedPos[1] + 1, latticeRoundedPos[2] + 1,
1623 latticeRoundedPos[3])) != material.end()
1624 && material.find(_superGeometry.get(latticeRoundedPos[0], latticeRoundedPos[1] + 1, latticeRoundedPos[2],
1625 latticeRoundedPos[3] + 1)) != material.end()
1626 && material.find(_superGeometry.get(latticeRoundedPos[0], latticeRoundedPos[1] + 1, latticeRoundedPos[2] + 1,
1627 latticeRoundedPos[3] + 1)) != material.end()) {
1628 getVel(fluidVel, &pos[0], globCuboid);
1629 vel[0] = fluidVel[0];
1630 vel[1] = fluidVel[1];
1631 vel[2] = fluidVel[2];
1632
1633 pCopy.setPos(pos);
1634 pCopy.setVel(vel);
1635 pCopy.setID(id);
1636 this->addParticle(pCopy);
1637 id++;
1638
1639 if (posDeq.size() <= (unsigned)deqSize) {
1640 posDeq.push_front(posVecTmp);
1641 }
1642 else {
1643 posDeq.push_front(posVecTmp);
1644 posDeq.pop_back();
1645 }
1646
1647 }
1648 }
1649 }
1650 }
1651 normalize(r);
1652 }
1653
1654
1655 template<>
1656 void SuperParticleSystem3D<double, MagneticParticle3D>::setMagneticParticlesdMomRandom()
1657 {
1658
1659 for (auto pS : _pSystems) {
1660 std::deque<MagneticParticle3D<double>*> particles = pS->getParticlesPointer();
1661
1662 for (auto p : particles) {
1663 std::vector<double> dMoment = { 0., 0., 0. };
1664 for (int i = 0; i < 3; i++) {
1665 dMoment[i] = rand() % (9 - (-9) + 1) + (-9);
1666 }
1667
1668 double dMoment_norm = util::sqrt(util::pow(dMoment[0], 2.) + util::pow(dMoment[1], 2.) + util::pow(dMoment[2], 2.)) ;
1669
1670 for (int i = 0; i < 3; i++) {
1671 dMoment[i] /= dMoment_norm ;
1672 }
1673
1674 p->setMoment(dMoment);
1675 }
1676 }
1677 }
1678
1679 template<typename T, template<typename U> class PARTICLETYPE>
1680 void SuperParticleSystem3D<T, PARTICLETYPE>::setParticlesVelRandom(T velFactor)
1681 {
1682
1683 for (auto pS : _pSystems) {
1684 std::deque<PARTICLETYPE<T>*> particles = pS->getParticlesPointer();
1685
1686 for (auto p : particles) {
1687 std::vector<T> vel = { 0., 0., 0. };
1688 for (int i = 0; i < 3; i++) {
1689 vel[i] = rand() % (9 - (-9) + 1) + (-9);
1690 }
1691
1692 T vel_norm = util::sqrt(util::pow(vel[0], 2.) + util::pow(vel[1], 2.) + util::pow(vel[2], 2.)) ;
1693
1694 for (int i = 0; i < 3; i++) {
1695 vel[i] /= vel_norm ;
1696 vel[i] *= velFactor ;
1697 }
1698
1699 p->setVel(vel);
1700 }
1701 }
1702 }
1703
1704
1705 template<typename T, template<typename U> class PARTICLETYPE>
1706 void SuperParticleSystem3D<T, PARTICLETYPE>::setParticlesPosRandom(T posFactor)
1707 {
1708
1709 for (auto pS : _pSystems) {
1710 std::deque<PARTICLETYPE<T>*> particles = pS->getParticlesPointer();
1711
1712 for (auto p : particles) {
1713 std::vector<T> pos = { 0., 0., 0. };
1714 for (int i = 0; i < 3; i++) {
1715 pos[i] = rand() % (9 - (-9) + 1) + (-9);
1716 }
1717
1718 T pos_norm = util::sqrt(util::pow(pos[0], 2.) + util::pow(pos[1], 2.) + util::pow(pos[2], 2.)) ;
1719
1720 for (int i = 0; i < 3; i++) {
1721 pos[i] /= pos_norm ;
1722 pos[i] *= posFactor ;
1723 }
1724
1725 for (int i = 0; i < 3; i++) {
1726 p->getPos()[0] += pos[0] ;
1727 p->getPos()[1] += pos[1] ;
1728 p->getPos()[2] += pos[2] ;
1729 }
1730 }
1731 }
1732 }
1733
1734 template<typename T, template<typename U> class PARTICLETYPE>
1735 void SuperParticleSystem3D<T, PARTICLETYPE>::setParticlesPosRandom(T posFactorX, T posFactorY, T posFactorZ)
1736 {
1737
1738 for (auto pS : _pSystems) {
1739 std::deque<PARTICLETYPE<T>*> particles = pS->getParticlesPointer();
1740
1741 for (auto p : particles) {
1742 std::vector<T> pos = { 0., 0., 0. };
1743 for (int i = 0; i < 3; i++) {
1744 pos[i] = rand() % (9 - (-9) + 1) + (-9);
1745 }
1746
1747 T pos_norm = util::sqrt(util::pow(pos[0], 2.) + util::pow(pos[1], 2.) + util::pow(pos[2], 2.)) ;
1748
1749 for (int i = 0; i < 3; i++) {
1750 pos[i] /= pos_norm ;
1751 }
1752
1753 pos[0] *= posFactorX ;
1754 pos[1] *= posFactorY ;
1755 pos[2] *= posFactorZ ;
1756
1757 for (int i = 0; i < 3; i++) {
1758 p->getPos()[0] += pos[0] ;
1759 p->getPos()[1] += pos[1] ;
1760 p->getPos()[2] += pos[2] ;
1761 }
1762 }
1763 }
1764 }
1765
1766 template<>
1767 void SuperParticleSystem3D<double, MagneticParticle3D>::setMagneticParticles(std::vector<double> dMoment,
1768 std::vector<double> vel, std::vector<double> aVel, std::vector<double> torque, double magnetisation)
1769 {
1770 int i = 0;
1771 for (auto pS : _pSystems) {
1772 std::deque<MagneticParticle3D<double>*> particles = pS->getParticlesPointer();
1773
1774 for (auto p : particles) {
1775
1776 p->setMoment(dMoment);
1777 p->setVel(vel);
1778 p->setAVel(aVel);
1779 p->setTorque(torque);
1780 p->setMagnetisation(magnetisation);
1781 p->setID(i);
1782 i++;
1783
1784 }
1785 }
1786 }
1787
1788 template<>
1789 void SuperParticleSystem3D<double, MagneticParticle3D>::setMagneticParticles(std::vector<double> dMoment,
1790 std::vector<double> vel, std::vector<double> aVel, std::vector<double> torque, double magnetisation, int sActivity)
1791 {
1792 int i = 0;
1793 for (auto pS : _pSystems) {
1794 std::deque<MagneticParticle3D<double>*> particles = pS->getParticlesPointer();
1795
1796 for (auto p : particles) {
1797
1798 p->setMoment(dMoment);
1799 p->setVel(vel);
1800 p->setAVel(aVel);
1801 p->setTorque(torque);
1802 p->setMagnetisation(magnetisation);
1803 p->setID(i);
1804 i++;
1805 p->setSActivity(sActivity);
1806 }
1807 }
1808 }
1809
1810 template<>
1811 void SuperParticleSystem3D<double, MagneticParticle3D>::prepareAgglomerates()
1812 {
1813 for (auto pS : _pSystems) {
1814 std::list<MagneticParticle3D<double>*> particlesList;
1815 pS->_Agglomerates.push_back(particlesList) ;
1816 }
1817 }
1818
1819 template<>
1820 void SuperParticleSystem3D<double, MagneticParticle3D>::initAggloParticles()
1821 {
1822 for (auto pS : _pSystems) {
1823 pS->initAggloParticles() ;
1824 }
1825 }
1826
1827 template<>
1828 void SuperParticleSystem3D<double, MagneticParticle3D>::findAgglomerates(int iT, int itVtkOutputMagParticles)
1829 {
1830 int pSi = 0;
1831 for (auto pS : _pSystems) {
1832 pS->findAgglomerates() ;
1833
1834 if (iT % itVtkOutputMagParticles == 0) {
1835
1836 clout << "Particlesystem number: " << pSi << std::endl;
1837 clout << "Number of non agglomerated particles" << ": " << pS->_Agglomerates[0].size() << std::endl;
1838 clout << "Number of agglomerated particles" << ": " << pS->size() - pS->_Agglomerates[0].size() << std::endl;
1839 clout << "Proportion of agglomeratet particles" << ": "
1840 << double(pS->size() - pS->_Agglomerates[0].size()) / double(pS->size()) * 100. << "%" << std::endl;
1841 clout << "Number of agglomerates" << ": " << pS->_Agglomerates.size() - 1 << std::endl;
1842 }
1843 pSi++;
1844 }
1845 }
1846
1847 template<typename T, template<typename U> class PARTICLETYPE>
1848 void SuperParticleSystem3D<T, PARTICLETYPE>::addParticleEquallyDistributed(
1849 IndicatorCuboid3D<T>& cuboid, T pMass, T pRad,/* number of particles on x, y, z axis*/
1850 int nox, int noy, int noz, std::vector<T> vel)
1851 {
1852 std::vector < T > pos(3, 0.);
1853 Vector<T, 3> minPos(cuboid.getMin());
1854 bool indic[1] = { false };
1855
1856 clout << "Number of particles to create: nox*noy*noz = "
1857 << nox * noy * noz << std::endl;
1858
1859 T xlength = cuboid.getMax()[0] - cuboid.getMin()[0];
1860 T ylength = cuboid.getMax()[1] - cuboid.getMin()[1];
1861 T zlength = cuboid.getMax()[2] - cuboid.getMin()[2];
1862 int modNox = nox - 1, modNoy = noy - 1, modNoz = noz - 1;
1863 if (nox == 1) {
1864 modNox = 1;
1865 }
1866 if (noy == 1) {
1867 modNoy = 1;
1868 }
1869 if (noz == 1) {
1870 modNoz = 1;
1871 }
1872
1873 for (int i = 0; i < nox; ++i) {
1874 pos[0] = minPos[0] + (T) (i) * xlength / modNox;
1875 for (int j = 0; j < noy; ++j) {
1876 pos[1] = minPos[1] + (T) (j) * ylength / modNoy;
1877 for (int k = 0; k < noz; ++k) {
1878 pos[2] = minPos[2] + (T) (k) * zlength / modNoz;
1879
1880 if (cuboid(indic, &pos[0])) {
1881 PARTICLETYPE<T> p(pos, vel, pMass, pRad);
1882 addParticle(p);
1883 }
1884 }
1885 }
1886 }
1887 clout << "Number of created particles = "
1888 << globalNumOfParticles() << std::endl;
1889 }
1890
1891 template<typename T, template<typename U> class PARTICLETYPE>
1892 void SuperParticleSystem3D<T, PARTICLETYPE>::addParticleEquallyDistributed(
1893 IndicatorCuboid3D<T>& cuboid, int nox, int noy, int noz, PARTICLETYPE<T>& p)
1894 {
1895 std::vector < T > pos(3, 0.);
1896 int id = 0;
1897 Vector<T, 3> minPos(cuboid.getMin());
1898 PARTICLETYPE<T> pCopy(p);
1899 bool indic[1] = { false };
1900
1901 clout << "Number of particles to create: nox*noy*noz = "
1902 << nox * noy * noz << std::endl;
1903
1904 T xlength = cuboid.getMax()[0] - cuboid.getMin()[0];
1905 T ylength = cuboid.getMax()[1] - cuboid.getMin()[1];
1906 T zlength = cuboid.getMax()[2] - cuboid.getMin()[2];
1907 int modNox = nox - 1, modNoy = noy - 1, modNoz = noz - 1;
1908 if (nox == 1) {
1909 modNox = 1;
1910 }
1911 if (noy == 1) {
1912 modNoy = 1;
1913 }
1914 if (noz == 1) {
1915 modNoz = 1;
1916 }
1917
1918 for (int i = 0; i < nox; ++i) {
1919 pos[0] = minPos[0] + (T) (i) * xlength / modNox;
1920 for (int j = 0; j < noy; ++j) {
1921 pos[1] = minPos[1] + (T) (j) * ylength / modNoy;
1922 for (int k = 0; k < noz; ++k) {
1923 pos[2] = minPos[2] + (T) (k) * zlength / modNoz;
1924
1925 if (cuboid(indic, &pos[0])) {
1926 pCopy.setPos(pos);
1927 pCopy.setID(id);
1928 addParticle(pCopy);
1929 id++;
1930 }
1931 }
1932 }
1933 }
1934 clout << "Number of created particles = "
1935 << globalNumOfParticles() << std::endl;
1936 }
1937
1938 template<typename T, template<typename U> class PARTICLETYPE>
1939 void SuperParticleSystem3D<T, PARTICLETYPE>::addTracerParticle(
1940 IndicatorF3D<T>& ind, int idTP, T mas, T rad, int noTP, std::vector<T> vel)
1941 {
1942 srand(clock());
1943 std::vector < T > pos(3, 0.);
1944 bool indic[1] = { false };
1945
1946 noTP += globalNumOfTracerParticles();
1947 while (globalNumOfTracerParticles() < noTP) {
1948 pos[0] = ind.getMin()[0]
1949 + (T) (rand() % 100000) / 100000. * (ind.getMax()[0] - ind.getMin()[0]);
1950 pos[1] = ind.getMin()[1]
1951 + (T) (rand() % 100000) / 100000. * (ind.getMax()[1] - ind.getMin()[1]);
1952 pos[2] = ind.getMin()[2]
1953 + (T) (rand() % 100000) / 100000. * (ind.getMax()[2] - ind.getMin()[2]);
1954
1955#ifdef PARALLEL_MODE_MPI
1956 singleton::mpi().bCast(&*pos.begin(), 3);
1957#endif
1958
1959 int x0, y0, z0, C;
1960 std::vector<int> locLat(4, 0);
1961 if (this->_cuboidGeometry.getFloorLatticeR(pos, locLat)) {
1962 C = locLat[0];
1963 if (this->_loadBalancer.rank(C) == singleton::mpi().getRank()) {
1964 x0 = locLat[1];
1965 y0 = locLat[2];
1966 z0 = locLat[3];
1967 if (_superGeometry.get(C, x0, y0, z0) == 1
1968 && _superGeometry.get(C, x0, y0 + 1, z0) == 1
1969 && _superGeometry.get(C, x0, y0, z0 + 1) == 1
1970 && _superGeometry.get(C, x0, y0 + 1, z0 + 1) == 1
1971 && _superGeometry.get(C, x0 + 1, y0, z0) == 1
1972 && _superGeometry.get(C, x0 + 1, y0 + 1, z0) == 1
1973 && _superGeometry.get(C, x0 + 1, y0, z0 + 1) == 1
1974 && _superGeometry.get(C, x0 + 1, y0 + 1, z0 + 1) == 1
1975 && ind(indic, &pos[0])
1976 && !util::nearZero(idTP)) {
1977 PARTICLETYPE<T> p(pos, vel, mas, rad, idTP, mas);
1978 addParticle(p);
1979 }
1980 }
1981 }
1982 }
1983 }
1984
1985
1986 template<typename T, template<typename U> class PARTICLETYPE>
1987 void SuperParticleSystem3D<T, PARTICLETYPE>::addParticleBoxMuller(
1988 IndicatorF3D<T>& ind, T partRho, T mu, T sigma, int no, T appProb, std::vector<T> vel)
1989 {
1990
1991 srand(clock());
1992
1993 //Randomization for Appearance-Likelyhood
1994 T rndmApp[1] = {(T) (rand() % 100000) / 100000.};
1995
1996#ifdef PARALLEL_MODE_MPI
1997 singleton::mpi().bCast(rndmApp, 1);
1998#endif
1999
2000 if (rndmApp[0] <= appProb ) {
2001
2002 std::vector<T> pos(3, 0.);
2003 T rad;
2004 T mas;
2005
2006 bool indic[1] = { false };
2007
2008 no += globalNumOfParticles();
2009 while (globalNumOfParticles() < no) {
2010 pos[0] = ind.getMin()[0]
2011 + (T) (rand() % 100000) / 100000. * (ind.getMax()[0] - ind.getMin()[0]);
2012 pos[1] = ind.getMin()[1]
2013 + (T) (rand() % 100000) / 100000. * (ind.getMax()[1] - ind.getMin()[1]);
2014 pos[2] = ind.getMin()[2]
2015 + (T) (rand() % 100000) / 100000. * (ind.getMax()[2] - ind.getMin()[2]);
2016
2017 //Normally Distributed Particle Radius (Box-Muller Method)
2018 T u1 = (T) (rand() % 100000) / 100000.;
2019 T u2 = (T) (rand() % 100000) / 100000.;
2020
2021 T x = util::cos(2 * M_PI * u1) * util::sqrt(-2 * util::log(u2));
2022 rad = mu + x * sigma;
2023 mas = 4. / 3. * M_PI * util::pow( rad, 3 ) * partRho;
2024
2025#ifdef PARALLEL_MODE_MPI
2026 singleton::mpi().bCast(&*pos.begin(), 3);
2027#endif
2028
2029 int x0, y0, z0, C;
2030 std::vector<int> locLat(4, 0);
2031 if (this->_cuboidGeometry.getFloorLatticeR(pos, locLat)) {
2032 C = locLat[0];
2033 if (this->_loadBalancer.rank(C) == singleton::mpi().getRank()) {
2034 x0 = locLat[1];
2035 y0 = locLat[2];
2036 z0 = locLat[3];
2037 if (_superGeometry.get(C, x0, y0, z0) == 1
2038 && _superGeometry.get(C, x0, y0 + 1, z0) == 1
2039 && _superGeometry.get(C, x0, y0, z0 + 1) == 1
2040 && _superGeometry.get(C, x0, y0 + 1, z0 + 1) == 1
2041 && _superGeometry.get(C, x0 + 1, y0, z0) == 1
2042 && _superGeometry.get(C, x0 + 1, y0 + 1, z0) == 1
2043 && _superGeometry.get(C, x0 + 1, y0, z0 + 1) == 1
2044 && _superGeometry.get(C, x0 + 1, y0 + 1, z0 + 1) == 1
2045 && ind(indic, &pos[0])) {
2046 PARTICLETYPE<T> p(pos, vel, mas, rad);
2047 addParticle(p);
2048 }
2049 }
2050 }
2051 }
2052 }
2053 }
2054
2055
2056 template<typename T, template<typename U> class PARTICLETYPE>
2057 void SuperParticleSystem3D<T, PARTICLETYPE>::addParticleWithDistance(
2058 IndicatorCuboid3D<T>& ind, T pMass, T pRad, std::vector<T> vel,
2059 T conc, T minDist, bool checkDist)
2060 {
2061
2062 srand(clock());
2063 bool indic[1] = { false };
2064 std::vector < T > pos(3, 0.);
2065 int size = T();
2066 T dist = T();
2067 Vector<T, 3> diff;
2068 int C;
2069
2070 T indicatorVol = (ind.getMax()[0] - ind.getMin()[0])
2071 * (ind.getMax()[1] - ind.getMin()[1])
2072 * (ind.getMax()[2] - ind.getMin()[2]);
2073
2074 int noParticles = (int) (conc * indicatorVol / (4. / 3 * M_PI * util::pow(pRad, 3)));
2075
2076 if (checkDist == true && (noParticles * util::pow(minDist, 3) * 4 / 3. * M_PI > indicatorVol) ) {
2077 std::cout << "Error: minDist too large" << std::endl;
2078 exit(-1);
2079 }
2080
2081 std::cout << " noparticles " << noParticles << std::endl;
2082
2083 noParticles += globalNumOfParticles();
2084
2085 while (globalNumOfParticles() < noParticles) {
2086
2087 pos[0] = ind.getMin()[0]
2088 + (T) (rand() % 100000) / 100000. * (ind.getMax()[0] - ind.getMin()[0]);
2089 pos[1] = ind.getMin()[1]
2090 + (T) (rand() % 100000) / 100000. * (ind.getMax()[1] - ind.getMin()[1]);
2091 pos[2] = ind.getMin()[2]
2092 + (T) (rand() % 100000) / 100000. * (ind.getMax()[2] - ind.getMin()[2]);
2093
2094#ifdef PARALLEL_MODE_MPI
2095 singleton::mpi().bCast(&*pos.begin(), 3);
2096#endif
2097
2098 std::vector<int> locLat(4, 0);
2099 if (this->_cuboidGeometry.getFloorLatticeR(pos, locLat)) {
2100 C = locLat[0];
2101 //related cuboidID of util::floor lattice position
2102 if (this->_loadBalancer.rank(C) == singleton::mpi().getRank()) {
2103
2104 if (ind(indic, &pos[0])) {
2105
2106 int psno = 0;
2107 int globIC = 0;
2108 for (unsigned int pS = 0; pS < _pSystems.size(); ++pS) {
2109 globIC = this->_loadBalancer.glob(pS);
2110 if (globIC == C) {
2111 size = _pSystems[pS]->sizeInclShadow();
2112 psno = pS;
2113 }
2114 }
2115 if (size == 0) {
2116 PARTICLETYPE<T> p(pos, vel, pMass, pRad);
2117 addParticle(p);
2118 }
2119 else {
2120 for (int j = 0; j < size;) {
2121 diff[0] = _pSystems[psno]->operator[](j).getPos()[0]
2122 - pos[0];
2123 diff[1] = _pSystems[psno]->operator[](j).getPos()[1]
2124 - pos[1];
2125 diff[2] = _pSystems[psno]->operator[](j).getPos()[2]
2126 - pos[2];
2127 dist = norm(diff);
2128
2129 if (dist < minDist) {
2130 goto marke;
2131 }
2132 if (j == (size - 1)) {
2133 PARTICLETYPE<T> p(pos, vel, pMass, pRad);
2134 addParticle(p);
2135 j += 1;
2136 }
2137 else {
2138 j += 1;
2139 }
2140 }
2141 }
2142 marke:
2143 ;
2144 }
2145 }
2146 }
2147 }
2148 }
2149
2150 template<typename T, template<typename U> class PARTICLETYPE>
2151 void SuperParticleSystem3D<T, PARTICLETYPE>::addShadowParticle(
2152 PARTICLETYPE<T>& p)
2153 {
2154 for (unsigned int pS = 0; pS < _pSystems.size(); ++pS) {
2155 int globIC = this->_loadBalancer.glob(pS);
2156 if (globIC != p.getCuboid() && checkCuboid(p, _overlap, globIC)
2157 && !checkCuboid(p, 0, globIC)) {
2158 _pSystems[pS]->addShadowParticle(p);
2159 }
2160 }
2161 }
2162
2163 template<typename T, template<typename U> class PARTICLETYPE>
2164 std::vector<ParticleSystem3D<T, PARTICLETYPE>*> SuperParticleSystem3D<T,
2165 PARTICLETYPE>::getParticleSystems()
2166 {
2167 return _pSystems;
2168 }
2169
2170 template<typename T, template<typename U> class PARTICLETYPE>
2171 ParticleSystem3D<T, PARTICLETYPE>& SuperParticleSystem3D<T, PARTICLETYPE>::operator[](
2172 int i)
2173 {
2174 return *(_pSystems[i]);
2175 }
2176
2177 template<typename T, template<typename U> class PARTICLETYPE>
2178 void SuperParticleSystem3D<T, PARTICLETYPE>::addForce(
2179 std::shared_ptr<Force3D<T, PARTICLETYPE> > f)
2180 {
2181 for (auto pS : _pSystems) {
2182 pS->addForce(f);
2183 }
2184 }
2185
2186 template<typename T, template<typename U> class PARTICLETYPE>
2187 void SuperParticleSystem3D<T, PARTICLETYPE>::addBoundary(
2188 std::shared_ptr<Boundary3D<T, PARTICLETYPE> > b)
2189 {
2190 for (auto pS : _pSystems) {
2191 pS->addBoundary(b);
2192 }
2193 }
2194
2195 template<typename T, template<typename U> class PARTICLETYPE>
2196 void SuperParticleSystem3D<T, PARTICLETYPE>::addParticleOperation(
2197 std::shared_ptr<ParticleOperation3D<T, PARTICLETYPE> > o)
2198 {
2199 for (auto pS : _pSystems) {
2200 pS->addParticleOperation(o);
2201 }
2202 }
2203
2204 template<typename T, template<typename U> class PARTICLETYPE>
2205 void SuperParticleSystem3D<T, PARTICLETYPE>::saveToFile(std::string name)
2206 {
2207 std::string fullName = createFileName(name) + ".particles";
2208
2209 int rank = 0;
2210
2211#ifdef PARALLEL_MODE_MPI
2212 int size = 1;
2213 size = singleton::mpi().getSize();
2214 rank = singleton::mpi().getRank();
2215#endif
2216
2217 if (rank == 0) {
2218 std::ofstream fout(fullName.c_str(), std::ios::trunc);
2219 if (!fout) {
2220 clout << "Error: could not open " << fullName << std::endl;
2221 }
2222 fout.close();
2223 }
2224#ifdef PARALLEL_MODE_MPI
2225 if (rank > 0) {
2226 int prev = rank - 1;
2227 int buffer = 0;
2228 MPI_Status status;
2229 MPI_Recv(&buffer, 1, MPI_INT, prev, 0, MPI_COMM_WORLD, &status);
2230 }
2231#endif
2232 for (auto pS : _pSystems) {
2233 pS->saveToFile(fullName);
2234 }
2235#ifdef PARALLEL_MODE_MPI
2236 if (rank < size - 1) {
2237 int next = rank + 1;
2238 int buffer = 0;
2239 MPI_Send(&buffer, 1, MPI_INT, next, 0, MPI_COMM_WORLD);
2240 }
2241#endif
2242 }
2243
2244 template<typename T, template<typename U> class PARTICLETYPE>
2245 void SuperParticleSystem3D<T, PARTICLETYPE>::addParticlesFromFile(
2246 std::string name, T mass, T radius)
2247 {
2248 std::string fullName = createFileName(name) + ".particles";
2249 std::ifstream fin(fullName.c_str());
2250
2251 std::string line;
2252 while (std::getline(fin, line)) {
2253 std::istringstream iss(line);
2254 T buffer[PARTICLETYPE<T>::serialPartSize];
2255 for (int i = 0; i < PARTICLETYPE<T>::serialPartSize; i++) {
2256 iss >> buffer[i];
2257 }
2258 PARTICLETYPE<T> p;
2259 p.unserialize(buffer);
2260 if ( !util::nearZero(radius) ) {
2261 p.setRad(radius);
2262 }
2263 if ( !util::nearZero(mass) ) {
2264 p.setMass(mass);
2265 }
2266 addParticle(p);
2267 }
2268 }
2269
2270 template<typename T, template<typename U> class PARTICLETYPE>
2271 void SuperParticleSystem3D<T, PARTICLETYPE>::addParticlesFromParaviewFile(
2272 std::string name)
2273 {
2274 std::string fullName = createFileName(name) + ".csv";
2275 std::ifstream fin(fullName.c_str());
2276
2277 std::string line;
2278 std::getline(fin, line);//first line is header
2279 //std::cout << line << std::endl;
2280 int counter = 0;
2281 T line_size = 13;
2282 while (std::getline(fin, line)) {
2283 //std::cout << counter << std::endl;
2284 std::istringstream iss(line);
2285
2286 T para_buffer[13];
2287 std::string A;
2288 iss >> A;
2289 T buffer[PARTICLETYPE<T>::serialPartSize];
2290 for (unsigned int i = 0; i < PARTICLETYPE<T>::serialPartSize; i++) {
2291
2292 buffer[i]=1.;
2293 }
2294 for (unsigned int i = 1; i < line_size; i++) {
2295 iss >> para_buffer[i];
2296 }
2297 for (unsigned int i=0;i<3;i++)
2298 {
2299 buffer[i] = para_buffer[i+5];
2300 buffer[i+3] = para_buffer[i+10];
2301 }
2302 buffer[9] = para_buffer[4];
2303 buffer[12] = para_buffer[9];
2304 buffer[14] = para_buffer[3];
2305 buffer[15] = counter;
2306 PARTICLETYPE<T> p;
2307 p.unserialize(buffer);
2308
2309 /*
2310 if ( !util::nearZero(radius) ) {
2311 p.setRad(radius);
2312 }
2313 if ( !util::nearZero(mass) ) {
2314 p.setMass(mass);
2315 }*/
2316 addParticle(p);
2317 if ( !util::nearZero(radius) ) {
2318 p.setRad(0.000001);
2319 }
2320 if ( !util::nearZero(mass) ) {
2321 p.setMass(1.);
2322 }
2323 counter++;
2324 }
2325 }
2326
2327
2328
2329 template<typename T, template<typename U> class PARTICLETYPE>
2330 void SuperParticleSystem3D<T, PARTICLETYPE>::clearParticles()
2331 {
2332 for (auto pS : _pSystems) {
2333 pS->clearParticles();
2334 }
2335 }
2336
2337 template<typename T, template<typename U> class PARTICLETYPE>
2338 void SuperParticleSystem3D<T, PARTICLETYPE>::setContactDetection(
2339 ContactDetection<T, PARTICLETYPE>& contactDetection)
2340 {
2341 clout << "Setting ContactDetectionAlgorithm " << contactDetection.getName() << std::endl;
2342
2343 for (auto pS : _pSystems) {
2344 pS->setContactDetection(contactDetection);
2345 }
2346 }
2347
2348 template<typename T, template<typename U> class PARTICLETYPE>
2349 void SuperParticleSystem3D<T, PARTICLETYPE>::setContactDetectionForPSys(
2350 ContactDetection<T, PARTICLETYPE>& contactDetection, int pSysNr)
2351 {
2352 clout << "Setting ContactDetectionAlgorithm for pSys: " << pSysNr
2353 << " = " << contactDetection.getName() << std::endl;
2354
2355 _pSystems[pSysNr]->setContactDetection(contactDetection);
2356 // this->getParticleSystems()[pSysNr]->setContactDetection(contactDetection);
2357 }
2358
2359
2360 //template<typename T, template<typename U> class PARTICLETYPE>
2361 //template<typename DESCRIPTOR>
2362 //void SuperParticleSystem3D<T, PARTICLETYPE>::particleOnFluid(
2363 // SuperLattice<T, DESCRIPTOR>& sLattice, T eps,
2364 // SuperGeometry<T,3>& sGeometry) {
2365 // for (unsigned int i = 0; i < _pSystems.size(); ++i) {
2366 // _pSystems[i]->particleOnFluid(
2367 // // sLattice.getBlock(i),
2368 // sLattice.getBlock(i),
2369 // sLattice.getCuboidGeometry().get(this->_loadBalancer.glob(i)),
2370 // sLattice.getOverlap(), eps, sGeometry.getBlockGeometry(i));
2371 // }
2372 //}
2373 //
2374 //template<typename T, template<typename U> class PARTICLETYPE>
2375 //template<typename DESCRIPTOR>
2376 //void SuperParticleSystem3D<T, PARTICLETYPE>::resetFluid(
2377 // SuperLattice<T, DESCRIPTOR>& sLattice) {
2378 // for (unsigned int i = 0; i < _pSystems.size(); ++i) {
2379 // _pSystems[i]->resetFluid(
2380 // // sLattice.getBlock(i),
2381 // sLattice.getBlock(i),
2382 // sLattice.getCuboidGeometry().get(this->_loadBalancer.glob(i)),
2383 // sLattice.getOverlap());
2384 // }
2385 //}
2386
2387 //template<typename T>
2388 //class SuperRotParticleSystem3D : public SuperParticleSystem3D<T, RotatingParticle3D> {
2389 // public:
2390 // SuperRotParticleSystem3D(SuperGeometry<T,3>& sg, UnitConverter<T,DESCRIPTOR>& conv): SuperParticleSystem3D<T, RotatingParticle3D>(sg, conv)
2391 // {};
2392 // void simulate(T dT) {
2393 // for (auto pS : this->_pSystems) {
2394 // pS->_contactDetection->sort();
2395 // pS->simulate(dT);
2396 // pS->computeTorque();
2397 // pS->computeBoundary();
2398 // }
2399 // this->updateParticleDistribution();
2400 // }
2401 //};
2402
2403}//namespace olb
2404
2405#endif /* SUPERPARTICLESYSTEM_3D_HH */
M_PI
#define M_PI
Definition analyticalF.hh:43
olb::AnalyticalConst
AnalyticalConst: DD -> XD, where XD is defined by value.size()
Definition analyticalF.h:78
olb::BackCouplingModel
Abstact base class for BaseBackCouplingModel.
Definition backCouplingModels.h:40
olb::BackCouplingModel::resetExternalField
virtual void resetExternalField(int material)=0
Resets external field.
olb::BackCouplingModel::communicate
virtual void communicate()=0
Communicates POPULATION and FORCE fields if the model is non-local.
olb::Boundary3D
Prototype for all particle boundaries.
Definition boundary3D.h:43
olb::CuboidGeometry3D
A cuboid geometry represents a voxel mesh.
Definition cuboidGeometry3D.h:72
olb::Force3D
Prototype for all particle forces.
Definition force3D.h:43
olb::ForwardCouplingModel
Abstact base class for all the forward-coupling models Its raison d'etre consists of not being temple...
Definition forwardCouplingModels3D.h:41
olb::IndicatorCircle3D
indicator function for a 3D circle
Definition indicatorF3D.h:77
olb::IndicatorCircle3D::getCenter
Vector< S, 3 > const & getCenter() const
Definition indicatorF3D.hh:111
olb::IndicatorCircle3D::getNormal
Vector< S, 3 > const & getNormal() const
Definition indicatorF3D.hh:117
olb::IndicatorCuboid3D
indicator function for a 3d-cuboid, parallel to the planes x=0, y=0, z=0.
Definition indicatorF3D.h:238
olb::IndicatorF3D
IndicatorF3D is an application from .
Definition indicatorBaseF3D.h:42
olb::IndicatorF3D::getMax
virtual Vector< S, 3 > & getMax()
Definition indicatorBaseF3D.hh:51
olb::IndicatorF3D::getMin
virtual Vector< S, 3 > & getMin()
Definition indicatorBaseF3D.hh:45
olb::LoadBalancer
Base class for all LoadBalancer.
Definition loadBalancer.h:59
olb::ParticleSystem3D::addForce
void addForce(std::shared_ptr< Force3D< T, PARTICLETYPE > > pF)
Add a force to ParticleSystem.
Definition particleSystem3D.hh:239
olb::SuperGeometry
Representation of a statistic for a parallel 2D geometry.
Definition superGeometry.h:70
olb::SuperParticleSystem3D
The class superParticleSystem is the basis for particulate flows within OpenLB.
Definition superParticleSystem3D.h:79
olb::SuperParticleSystem3D::addParticleOperation
void addParticleOperation(std::shared_ptr< ParticleOperation3D< T, PARTICLETYPE > > o)
Add an operation to system.
Definition superParticleSystem3D.hh:2196
olb::SuperParticleSystem3D::saveToFile
void saveToFile(std::string name)
Save Particles to file. Add using addParticlesFromFile(std::string name, T mass, T radius);.
Definition superParticleSystem3D.hh:2205
olb::SuperParticleSystem3D::setParticlesPosRandom
void setParticlesPosRandom(T posFactor)
Changes particle positions randomly.
Definition superParticleSystem3D.hh:1706
olb::SuperParticleSystem3D::setMagneticParticlesdMomRandom
void setMagneticParticlesdMomRandom()
Gives random dipolemoment orientation to all MagneticParticle3D.
Definition superParticleSystem3D.h:199
olb::SuperParticleSystem3D::addBoundary
void addBoundary(std::shared_ptr< Boundary3D< T, PARTICLETYPE > > b)
Add a boundary to system.
Definition superParticleSystem3D.hh:2187
olb::SuperParticleSystem3D::setOverlap
void setOverlap(T)
Set overlap of ParticleSystems, overlap has to be in lattice units particle system _overlap+1 <= _sup...
Definition superParticleSystem3D.hh:614
olb::SuperParticleSystem3D::clearParticles
void clearParticles()
Removes all particles from System.
Definition superParticleSystem3D.hh:2330
olb::SuperParticleSystem3D::updateParticleDistribution
void updateParticleDistribution()
Redistribute particles on compute nodes.
Definition superParticleSystem3D.hh:852
olb::SuperParticleSystem3D::globalNumOfShadowParticles
int globalNumOfShadowParticles()
Get global number of shadow particles (particles hold in overlap)
Definition superParticleSystem3D.hh:649
olb::SuperParticleSystem3D::setParticlesVelRandom
void setParticlesVelRandom(T velFactor)
Gives random velocity to all particles.
Definition superParticleSystem3D.hh:1680
olb::SuperParticleSystem3D::addParticleEquallyDistributed
void addParticleEquallyDistributed(IndicatorCuboid3D< T > &cuboid, T pMass, T pRad, int nox, int noy, int noz, std::vector< T > vel={0., 0., 0.})
Add a number of identical particles equally distributed in a given Material Number.
Definition superParticleSystem3D.hh:1848
olb::SuperParticleSystem3D::getOverlap
T getOverlap()
Get overlap of ParticleSystems.
Definition superParticleSystem3D.hh:621
olb::SuperParticleSystem3D::printDeep
void printDeep(std::string message="")
Definition superParticleSystem3D.hh:179
olb::SuperParticleSystem3D::simulateWithTwoWayCoupling_Davide
void simulateWithTwoWayCoupling_Davide(T dT, ForwardCouplingModel< T, PARTICLETYPE > &forwardCoupling, BackCouplingModel< T, PARTICLETYPE > &backCoupling, int material, int subSteps=1, bool resetExternalField=true, bool scale=false)
Definition superParticleSystem3D.hh:565
olb::SuperParticleSystem3D::addParticle
void addParticle(PARTICLETYPE< T > &p)
Add a Particle to SuperParticleSystem.
Definition superParticleSystem3D.hh:1103
olb::SuperParticleSystem3D::init
void init()
Init the SuperParticleSystem.
Definition superParticleSystem3D.hh:73
olb::SuperParticleSystem3D::addForce
void addForce(std::shared_ptr< Force3D< T, PARTICLETYPE > > f)
Add a force to system.
Definition superParticleSystem3D.hh:2178
olb::SuperParticleSystem3D::prepareAgglomerates
void prepareAgglomerates()
Agglomerate detection functions: Todo: enable for parallel mode Initializes an empty agglomerate list...
Definition superParticleSystem3D.h:208
olb::SuperParticleSystem3D::rankNumOfShadowParticles
int rankNumOfShadowParticles()
Get number of shadow particles computed on this node.
Definition superParticleSystem3D.hh:671
olb::SuperParticleSystem3D::addParticlesFromParaviewFile
void addParticlesFromParaviewFile(std::string name)
Add particles form a File. Save using saveToFile(std::string name)
Definition superParticleSystem3D.hh:2271
olb::SuperParticleSystem3D::addParticleWithDistance
void addParticleWithDistance(IndicatorCuboid3D< T > &ind, T pMass, T pRad, std::vector< T > vel, T conc, T minDist, bool checkDist)
Generates particles with specific volume concentration conc equally and randomly distributed in given...
Definition superParticleSystem3D.hh:2057
olb::SuperParticleSystem3D::captureEscapeRate
void captureEscapeRate(std::list< int > mat)
console output of escape (E), capture (C) rate for material numbers mat
Definition superParticleSystem3D.hh:201
olb::SuperParticleSystem3D::findAgglomerates
void findAgglomerates(int iT, int itVtkOutputMagParticles)
Detects and manages particle agglomerates.
Definition superParticleSystem3D.h:212
olb::SuperParticleSystem3D::checkCuboid
bool checkCuboid(PARTICLETYPE< T > &p, T overlap)
Check if particle is still on cuboid.
Definition superParticleSystem3D.hh:834
olb::SuperParticleSystem3D::setVelToAnalyticalVel
void setVelToAnalyticalVel(AnalyticalConst3D< T, T > &)
Set particle velocity to analytical velocity (e.g. as initial condition.
Definition superParticleSystem3D.hh:800
olb::SuperParticleSystem3D::findCuboid
bool findCuboid(PARTICLETYPE< T > &, int overlap)
Find the cuboid the particle is on.
Definition superParticleSystem3D.hh:815
olb::SuperParticleSystem3D::addParticleBoxMuller
void addParticleBoxMuller(IndicatorF3D< T > &ind, T partRho, T mu, T sigma, int no=1, T appProb=1., std::vector< T > vel={0., 0., 0.})
Add a number of nonidentical particles with normally distributed radius (Box-Muller Method) in a give...
Definition superParticleSystem3D.hh:1987
olb::SuperParticleSystem3D::getOutput
void getOutput(std::string filename, int iT, T conversionFactorTime, unsigned short particleProperties)
Get Output of particleMovement Write the data of the particle movement into an txtFile.
Definition superParticleSystem3D.hh:340
olb::SuperParticleSystem3D::globalNumOfActiveParticles
int globalNumOfActiveParticles()
Get global number of active particles.
Definition superParticleSystem3D.hh:713
olb::SuperParticleSystem3D::diffEscapeRate
void diffEscapeRate(std::list< int > mat, int &globalPSum, int &pSumOutlet, T &diffEscapeRate, T &maxDiffEscapeRate, int iT, int iTConsole, T genPartPerTimeStep=0)
Console output of differential escape rate for material numbers mat (e.g.
Definition superParticleSystem3D.hh:215
olb::SuperParticleSystem3D::addShadowParticle
void addShadowParticle(PARTICLETYPE< T > &p)
Add a shadow particle to system.
Definition superParticleSystem3D.hh:2151
olb::SuperParticleSystem3D::countMaterial
int countMaterial(int mat)
Get number of particles in the vicinity of material number mat.
Definition superParticleSystem3D.hh:701
olb::SuperParticleSystem3D::simulateWithTwoWayCoupling_Mathias
void simulateWithTwoWayCoupling_Mathias(T dT, ForwardCouplingModel< T, PARTICLETYPE > &forwardCoupling, BackCouplingModel< T, PARTICLETYPE > &backCoupling, int material, int subSteps=1, bool resetExternalField=true, bool scale=false)
Integrate on Timestep dT with two-way coupling, scale = true keeps the particle velocity in stable ra...
Definition superParticleSystem3D.hh:543
olb::SuperParticleSystem3D::rankNumOfParticles
int rankNumOfParticles()
Get number of particles computed on this node.
Definition superParticleSystem3D.hh:661
olb::SuperParticleSystem3D::addTracerParticle
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...
Definition superParticleSystem3D.hh:1939
olb::SuperParticleSystem3D::setVelToFluidVel
void setVelToFluidVel(SuperLatticeInterpPhysVelocity3D< T, DESCRIPTOR > &)
Set particle velocity to fluid velocity (e.g. as initial condition.
Definition superParticleSystem3D.hh:791
olb::SuperParticleSystem3D::simulate
void simulate(T dT, bool scale=false)
Integrate on Timestep dT, scale = true keeps the particle velocity in stable range.
Definition superParticleSystem3D.hh:488
olb::SuperParticleSystem3D::rankNumOfTracerParticles
int rankNumOfTracerParticles()
Get number of TracerParticles computed on this node.
Definition superParticleSystem3D.hh:737
olb::SuperParticleSystem3D::setMagneticParticles
void setMagneticParticles(std::vector< T > dMoment, std::vector< T > vel, std::vector< T > aVel, std::vector< T > torque, T magnetisation)
Gives specific attributes to all MagneticParticle3D.
Definition superParticleSystem3D.h:201
olb::SuperParticleSystem3D::particleSActivityTest
bool particleSActivityTest(int sActivity)
Tests if particles with specific sActivity exist.
Definition superParticleSystem3D.h:215
olb::SuperParticleSystem3D::globalNumOfTracerParticles
int globalNumOfTracerParticles()
Get number of TracerParticles computed on this node.
Definition superParticleSystem3D.hh:725
olb::SuperParticleSystem3D::numOfForces
std::vector< int > numOfForces()
Get number of linked Forces.
Definition superParticleSystem3D.hh:780
olb::SuperParticleSystem3D::SuperParticleSystem3D
SuperParticleSystem3D(CuboidGeometry3D< T > &cuboidGeometry, LoadBalancer< T > &loadBalancer, SuperGeometry< T, 3 > &)
Constructor for SuperParticleSystem.
Definition superParticleSystem3D.hh:47
olb::SuperParticleSystem3D::rankNumOfActiveParticles
int rankNumOfActiveParticles()
Get number of active particles computed on this node.
Definition superParticleSystem3D.hh:681
olb::SuperParticleSystem3D::operator[]
ParticleSystem3D< T, PARTICLETYPE > & operator[](int i)
Get a ParticleSystem.
Definition superParticleSystem3D.hh:2171
olb::SuperParticleSystem3D::generateParticlesCircleInletMassConcentration
void generateParticlesCircleInletMassConcentration(IndicatorCircle3D< T > &indicatorCircle, T particleMassConcentration, T charPhysVelocity, T conversionFactorTime, SuperLatticeInterpPhysVelocity3D< T, DESCRIPTOR > &getVel, PARTICLETYPE< T > &p, std::set< int > material, int iT, T &particlesPerPhyTimeStep, std::vector< T > &inletVec, std::deque< std::vector< T > > &posDeq, int deqSize)
Generates particle at a circle shaped inlet, amount given by mass concentration in feedstream.
Definition superParticleSystem3D.hh:1497
olb::SuperParticleSystem3D::setContactDetection
void setContactDetection(ContactDetection< T, PARTICLETYPE > &contactDetection)
Set contact detection algorithm for particle-particle contact. Not yet implemented.
Definition superParticleSystem3D.hh:2338
olb::SuperParticleSystem3D::globalNumOfParticles
int globalNumOfParticles()
Get global number of particles.
Definition superParticleSystem3D.hh:633
olb::SuperParticleSystem3D::addParticlesFromFile
void addParticlesFromFile(std::string name, T mass, T radius)
Add particles form a File. Save using saveToFile(std::string name)
Definition superParticleSystem3D.hh:2245
olb::SuperParticleSystem3D::addHL3DParticle
void addHL3DParticle(IndicatorF3D< T > &ind, std::set< int > material, T mas, T rad, int no=1, std::vector< T > vel={0., 0., 0.}, T surface=1., T volume=1.)
Definition superParticleSystem3D.hh:1302
olb::SuperParticleSystem3D::setContactDetectionForPSys
void setContactDetectionForPSys(ContactDetection< T, PARTICLETYPE > &contactDetection, int pSysNr)
Set contact detection algorithm for particle-particle contact. Not yet implemented.
Definition superParticleSystem3D.hh:2349
olb::SuperParticleSystem3D::initAggloParticles
void initAggloParticles()
Adds new generated particles to the list of non agglomerated Particles.
Definition superParticleSystem3D.h:210
olb::SuperParticleSystem3D::numOfPSystems
int numOfPSystems()
Get number of ParticleSystems.
Definition superParticleSystem3D.hh:627
olb::Vector
Plain old scalar vector.
Definition vector.h:47
olb::singleton::MpiManager::bCast
void bCast(T *sendBuf, int sendCount, int root=0, MPI_Comm comm=MPI_COMM_WORLD)
Broadcast data from one processor to multiple processors.
olb::singleton::MpiManager::getSize
int getSize() const
Returns the number of processes.
Definition mpiManager.cpp:65
olb::singleton::MpiManager::iSend
void iSend(T *buf, int count, int dest, MPI_Request *request, int tag=0, MPI_Comm comm=MPI_COMM_WORLD)
Sends data at *buf, non blocking.
olb::singleton::MpiManager::reduceAndBcast
void reduceAndBcast(T &reductVal, MPI_Op op, int root=0, MPI_Comm comm=MPI_COMM_WORLD)
Reduction operation, followed by a broadcast.
olb::singleton::MpiManager::getRank
int getRank() const
Returns the process ID.
Definition mpiManager.cpp:70
olb::singleton::MpiManager::waitAll
void waitAll(MpiNonBlockingHelper &mpiNbHelper)
Complete a series of non-blocking MPI operations.
Definition mpiManager.cpp:1303
olb::singleton::MpiManager::barrier
void barrier(MPI_Comm comm=MPI_COMM_WORLD)
Synchronizes the processes.
Definition mpiManager.cpp:93
olb::singleton::MpiManager::receive
void receive(T *buf, int count, int source, int tag=0, MPI_Comm comm=MPI_COMM_WORLD)
Receives data at *buf, blocking.
olb::singleton::MpiNonBlockingHelper
Helper class for non blocking MPI communication.
Definition mpiManager.h:51
olb::singleton::MpiNonBlockingHelper::allocate
void allocate(unsigned i)
Allocates memory.
Definition mpiManager.cpp:1323
olb::singleton::MpiNonBlockingHelper::get_mpiRequest
MPI_Request * get_mpiRequest(int i=0) const
Get the specified request object.
Definition mpiManager.cpp:1341
fileName.h
These functions help you to create file names.
olb::descriptors::m
constexpr T m(unsigned iPop, unsigned jPop, tag::MRT)
Definition mrtLatticeDescriptors.h:394
olb::singleton::mpi
MpiManager & mpi()
Definition mpiManager.cpp:29
olb::util::sqrt
cpu::simd::Pack< T > sqrt(cpu::simd::Pack< T > value)
Definition pack.h:100
olb::util::min
cpu::simd::Pack< T > min(cpu::simd::Pack< T > rhs, cpu::simd::Pack< T > lhs)
Definition pack.h:124
olb::util::log
ADf< T, DIM > log(const ADf< T, DIM > &a)
Definition aDiff.h:475
olb::util::max
cpu::simd::Pack< T > max(cpu::simd::Pack< T > rhs, cpu::simd::Pack< T > lhs)
Definition pack.h:130
olb::util::pow
cpu::simd::Pack< T > pow(cpu::simd::Pack< T > base, cpu::simd::Pack< T > exp)
Definition pack.h:112
olb::util::nearZero
bool nearZero(const ADf< T, DIM > &a)
Definition aDiff.h:1087
olb::util::cos
ADf< T, DIM > cos(const ADf< T, DIM > &a)
Definition aDiff.h:578
olb
Top level namespace for all of OpenLB.
Definition boundaryPostProcessors2D.h:34
olb::createFileName
std::string createFileName(std::string name)
for .pvd masterFile
Definition fileName.hh:34
olb::toStdVector
constexpr std::vector< T > toStdVector(const ScalarVector< T, D, IMPL > &a)
Copies data into a standard vector.
Definition scalarVector.h:88
olb::norm
constexpr T norm(const ScalarVector< T, D, IMPL > &a)
Euclidean vector norm.
Definition scalarVector.h:66
olb::crossProduct3D
constexpr Vector< T, 3 > crossProduct3D(const ScalarVector< T, 3, IMPL > &a, const ScalarVector< T, 3, IMPL_ > &b) any_platform
Definition vector.h:224
olb::normalize
constexpr Vector< T, D > normalize(const ScalarVector< T, D, IMPL > &a, T scale=T{1})
Definition vector.h:245
Generated on Tue Mar 5 2024 09:57:52 for OpenLB by doxygen 1.10.0