MercuryDPM  Trunk
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
InsertionBoundary.cc
Go to the documentation of this file.
1 //Copyright (c) 2013-2020, The MercuryDPM Developers Team. All rights reserved.
2 //For the list of developers, see <http://www.MercuryDPM.org/Team>.
3 //
4 //Redistribution and use in source and binary forms, with or without
5 //modification, are permitted provided that the following conditions are met:
6 // * Redistributions of source code must retain the above copyright
7 // notice, this list of conditions and the following disclaimer.
8 // * Redistributions in binary form must reproduce the above copyright
9 // notice, this list of conditions and the following disclaimer in the
10 // documentation and/or other materials provided with the distribution.
11 // * Neither the name MercuryDPM nor the
12 // names of its contributors may be used to endorse or promote products
13 // derived from this software without specific prior written permission.
14 //
15 //THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
16 //ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
17 //WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
18 //DISCLAIMED. IN NO EVENT SHALL THE MERCURYDPM DEVELOPERS TEAM BE LIABLE FOR ANY
19 //DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
20 //(INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
21 //LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
22 //ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
23 //(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
24 //SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
25 
26 #include "InsertionBoundary.h"
27 #include "DPMBase.h"
28 #include "Particles/BaseParticle.h"
29 #include<iostream>
30 
31 #ifdef MERCURY_USE_MPI
32 #include "MpiDataClass.h"
33 #endif
34 
39 {
41  massInserted_ = 0;
42  volumeInserted_ = 0;
43  particleToCopy_ = nullptr;
44  maxFailed_ = 0;
45  isActivated_ = true;
47  initialVolume_ = 0;
50  radMin_ = 0;
51  radMax_ = 0;
52  isManuallyInserting_ = false;
53 }
54 
59  : BaseBoundary(other)
60 {
64  maxFailed_ = other.maxFailed_;
65  isActivated_ = other.isActivated_;
73  radMin_ = other.radMin_;
74  radMax_ = other.radMax_;
76 
77  if (other.particleToCopy_ != nullptr)
78  {
80  }
81  else
82  {
83  particleToCopy_ = nullptr;
84  }
85 }
86 
92 {
93  delete particleToCopy_;
94 }
95 
102 void InsertionBoundary::set(BaseParticle* particleToCopy, unsigned int maxFailed)
103 {
104  if (particleToCopy != nullptr)
105  particleToCopy_ = particleToCopy->copy();
106  maxFailed_ = maxFailed;
107 }
108 
109 
115 {
118  {
119  switch (distribution_)
120  {
121  default:
122  // default is a uniform distribution between radMin_ and radMax_
124  break;
126  // normal distribution
128  static std::mt19937 gen(0);
129  Mdouble particleRadius = 0.5 * (radMax_ + radMin_);
130  Mdouble polydispersity = 0.5 * (radMax_ - radMin_);
131  static std::normal_distribution<> d(particleRadius, polydispersity);
132  static const Mdouble radiusMin = particleRadius - 1.5 * polydispersity;
133  static const Mdouble radiusMax = particleRadius + 1.5 * polydispersity;
134  Mdouble radius = d(gen);
135  while (radius > radiusMax || radius < radiusMin) radius = d(gen);
136  P->setRadius(radius);
137  break;
138  }
139  }
140  // manual insertion routine to insert PSDs as accurate as possible into a given volume
141  else if (isManuallyInserting_)
142  {
143  logger.assert(initialVolume_ > 0.0, "Use setInitialVolume to define the particle insertion volume");
144  Mdouble radius;
145  // getVolumeFlowRate() * time + initialVolume_ - volumeInserted_ lead to more inaccurate results, therfore
146  // -volumeInserted was removed.
148  getHandler()->getDPMBase()->getTime() +
150  P->setRadius(radius);
151  }
152  else
153  {
154  Mdouble radius;
156  P->setRadius(radius);
157  }
158  return P;
159 }
160 
162 {
163  // check if the flow rate limit has been reached
165  {
167  }
168  else
169  {
170  const Mdouble iMax = (Mdouble) variableCumulativeVolumeFlowRate_.size() - 2;
171  const Mdouble i = std::min(time / samplingInterval_, iMax);
172  if (i == iMax)
173  {
174  static unsigned count = 0;
175  if (count == 0) logger(WARN, "Reached end of volume flowrate function");
176  ++count;
177  }
178  const size_t id = i;
179  const Mdouble allowedVolume = variableCumulativeVolumeFlowRate_[id] +
181  variableCumulativeVolumeFlowRate_[id]) * (i - id);
182  return volumeInserted_ < allowedVolume;
183  }
184 }
185 
193 {
194  logger(VERBOSE, "In InsertionBoundary::checkBoundaryBeforeTimeStep\n");
195 
196  if (!isActivated_)
197  return;
198 
199  /* Each timestep, the InsertionBoundary attempts to fill up a region with
200  * particles.
201  *
202  * It first calls generateParticle() to get a randomised particle, subject
203  * to a specified distribution over sizes and species. (The basic class
204  * supports size dispersity only but PolydisperseInsertionBoundary will
205  * support species dispersity.)
206  * Then it repeatedly calls placeParticle(), which gives the particle a
207  * random location (and possibly velocity) in a specified,
208  * geometry-dependent bound. Each time, it checks whether the new particle
209  * would have an interaction with another particle or a wall.
210  *
211  * If it manages to do that within maxFailed_ tries, then:
212  * * the new particle is inserted,
213  * * the failure counter is reset, and
214  * the processes is repeated with a new generateParticle().
215  *
216  * Otherwise, the processes terminates for this timestep.
217  * */
218 
219  // Keep count of how many successive times we have failed to place a new
220  // particle.
221  unsigned int failed = 0;
222  while (failed <= maxFailed_ && insertParticle(md->getNextTime())) // 'generating' loop
223  {
224  /* Generate random *intrinsic* properties for the new particle. */
225  logger(VERBOSE, "about to call generateParticle\n");
226 
227 
228 
229 
231 
232 
233 
234 
235  auto p0 = generateParticle(md->random);
236  // Important for particle generation with a particle size distribution as it generates a particle with zero
237  // radius. If a particle is not allowed to be inserted by the PSD criteria it will generate a particle with
238  // zero diameter. This if statement prevents inserting particles with zero radius, which would else be a problem.
240  if (p0->getRadius() == 0)
241  {
242  logger(VERBOSE, "The PSD for the specified volume is fully set");
243  // free up memory space
244  delete p0;
245  // out of the 'placing' loop
246  failed = maxFailed_ + 1;
247  continue;
248  }
249  logger(VERBOSE, "generated a particle with intrinsics %", p0);
250 
251  while (true) // 'placing' loop
252  {
253  /* Generate extrinsic properties (position and velocity) for this
254  * new particle. */
255 
256 
257 
259 
260 
261 
262 
263 
264  placeParticle(p0, md->random);
265  logger(VERBOSE, "attempting to place particle at %, vel %", p0->getPosition(), p0->getVelocity());
266 
267 #ifdef MERCURY_USE_MPI
268  /* Communicate the new particle's properties by setHandler (note
269  * that this doesn't actually add the particle to the handler). */
270  if (NUMBER_OF_PROCESSORS > 1)
271  {
272  MPIParticle particle;
273  // //Every domain generates a particle (to get the species right etc)
274 
275  //Send particle data from root to other processors to sync the particle properties
276  if (PROCESSOR_ID == 0)
277  {
279  }
280 
282 
283  //Process the received data
284  if (PROCESSOR_ID != 0)
285  {
286  copyDataFromMPIParticleToParticle(&particle, p0, &(md->particleHandler));
287  }
288  }
289 #endif
290  p0->setHandler(&md->particleHandler);
291  /* Check whether the particle has any interactions. */
293  {
294  //Note: in parallel only one of the domains will actually add the particle
295  auto p = md->particleHandler.copyAndAddObject(p0);
296  failed = 0;
297 
299  const double volume = p0->getVolume();
300  volumeInserted_ += volume;
301  massInserted_ += p0->getSpecies()->getDensity() * volume;
302  logger(VERBOSE, "successfully placed a particle %, with position: % after % fails.", p0,
303  p0->getPosition(), failed);
304 
305  /* JMFT: The generateParticle() routine allocates memory, so we should
306  * free it here. (Don't worry, the particle will have been copied to the
307  * particleHandler by this point iff we want it.) */
308  delete p0;
309 
310  break; // out of the 'placing' loop
311  }
312  else
313  {
314  failed++;
315  logger(VERBOSE, "failed to place a particle; have failed % times", failed);
316  }
317 
318  if (failed > maxFailed_)
319  {
320  logger(VERBOSE, "failed too many times; giving up");
321  break; // out of the 'placing' loop (and will leave the 'generating' loop too
322  }
323  }
324  logger(VERBOSE, "failed % times, so breaking out of InsertionBoundary loop for this timestep.", failed);
325  }
326  // logger(INFO, "volumeInserted_ = %", volumeInserted_);
327 }
328 
334 {
336 }
337 
339 {
340  return massInserted_;
341 }
342 
344 {
345  return volumeInserted_;
346 }
347 
352 {
354  massInserted_ = 0;
355  volumeInserted_ = 0;
356 }
357 
359 {
360  isActivated_ = true;
361 }
362 
364 {
365  isActivated_ = false;
366 }
367 
369 {
370  return isActivated_;
371 }
372 
378 void InsertionBoundary::setMaxFailed(unsigned int maxFailed)
379 {
380  maxFailed_ = maxFailed;
381 }
382 
388 {
389  return maxFailed_;
390 }
391 
397 {
398  if (particleToCopy != nullptr)
399  particleToCopy_ = particleToCopy->copy();
400  else
401  logger(ERROR, "Setting particleToCopy to be a null pointer?");
402 }
403 
408 {
410  //if (particleToCopy_==nullptr)
411  // std::cerr << "Error: particleToCopy not set" << std::endl;
412  return particleToCopy_;
413 }
414 
419 void InsertionBoundary::read(std::istream& is)
420 {
421  BaseBoundary::read(is);
422  std::string dummy, type;
423  is >> dummy >> maxFailed_ >> dummy;
424  if (dummy == "volumeFlowRate")
425  is >> volumeFlowRate_ >> dummy;
426  is >> massInserted_;
427  is >> dummy >> volumeInserted_;
428  is >> dummy >> numberOfParticlesInserted_;
429  is >> dummy >> isActivated_;
431  {
433  {
434  is >> dummy >> p.radius;
435  is >> dummy >> p.probability;
436  }
437  }
438  else
439  {
440  is >> dummy >> distribution_;
441  }
443  helpers::readOptionalVariable(is, "checkParticleForInteraction", checkParticleForInteraction_);
444  helpers::readOptionalVariable(is, "initialVolume", initialVolume_);
445  if (helpers::readOptionalVariable(is, "samplingInterval", samplingInterval_))
446  {
447  size_t n;
448  Mdouble flowRate;
449  is >> dummy >> n;
450  //variableCumulativeVolumeFlowRate_.clear();
452  for (size_t i = 0; i < n; ++i)
453  {
454  is >> flowRate;
455  variableCumulativeVolumeFlowRate_.push_back(flowRate);
456  }
457  }
458  is >> dummy;
459  if (dummy != "noParticleToCopy")
460  {
461  delete particleToCopy_;
463  // The .restart file records the index of the particle's species, but
464  // doesn't record the pointer, i.e. the memory address of the species within
465  // the speciesHandler. The latter needs to be reset now.
466  particleToCopy_->setSpecies(getHandler()->getDPMBase()->speciesHandler.getObject(
468  ));
469  }
470 }
471 
476 void InsertionBoundary::write(std::ostream& os) const
477 {
478  logger(VERBOSE, "In InsertionBoundary::write\n");
480  os << " maxFailed " << maxFailed_;
481  if (std::isfinite(volumeFlowRate_))
482  os << " volumeFlowRate " << volumeFlowRate_;
483  os << " massInserted " << massInserted_;
484  os << " volumeInserted " << volumeInserted_;
485  os << " numberOfParticlesInserted " << numberOfParticlesInserted_;
486  os << " isActivated " << isActivated_;
488  {
489  os << " psd " << particleSizeDistribution_.getParticleSizeDistribution().size();
491  {
492  os << " " << p.radius
493  << " " << p.probability;
494  }
495  }
496  else
497  {
498  os << " distribution " << distribution_;
499  }
500  if (checkParticleForInteraction_ == false)
501  {
502  os << " checkParticleForInteraction " << checkParticleForInteraction_;
503  }
504  os << " initialVolume " << initialVolume_;
505  os << " samplingInterval " << samplingInterval_;
506  os << " variableCumulativeVolumeFlowRate " << variableCumulativeVolumeFlowRate_.size();
507  for (const auto flowRate : variableCumulativeVolumeFlowRate_)
508  {
509  os << ' ' << flowRate;
510  }
511  if (particleToCopy_ == nullptr)
512  os << " noParticleToCopy";
513  else
514  {
515  os << " particleToCopy " << *particleToCopy_;
516  }
517 }
518 
520 {
521  return volumeFlowRate_;
522 }
523 
525 {
526  volumeFlowRate_ = volumeFlowRate;
527 }
528 
530 {
531  return initialVolume_;
532 }
533 
535 {
536  initialVolume_ = initialVolume;
537  if (!std::isfinite(volumeFlowRate_)) volumeFlowRate_ = 0;
538 }
539 
540 void InsertionBoundary::setVariableVolumeFlowRate(const std::vector<Mdouble>& variableCumulativeVolumeFlowRate,
541  Mdouble samplingInterval)
542 {
543  logger.assert(samplingInterval > 0, "sampling interval needs to be positive");
544  const Mdouble endTime = variableCumulativeVolumeFlowRate.size() * samplingInterval;
545  logger(INFO, "variable flowrate is defined up to %", endTime);
546  logger.assert_always(getHandler()->getDPMBase()->getTimeMax() < endTime,
547  "variable flowrate is defined up to %, but tMax is set to %", endTime,
548  getHandler()->getDPMBase()->getTimeMax());
549  variableCumulativeVolumeFlowRate_ = variableCumulativeVolumeFlowRate;
550  samplingInterval_ = samplingInterval;
551 }
552 
558 {
560 }
561 
567 {
569 }
570 
575 void InsertionBoundary::setManualInsertion(bool isManuallyInserting)
576 {
577  isManuallyInserting_ = isManuallyInserting;
578 }
579 
581 {
582  distribution_ = distribution;
583 }
584 
589 {
590  return distribution_;
591 }
592 
596 std::ostream& operator<<(std::ostream& os, InsertionBoundary::Distribution type)
597 {
598  os << static_cast<unsigned>(type);
599  return os;
600 }
601 
605 std::istream& operator>>(std::istream& is, InsertionBoundary::Distribution& type)
606 {
607  unsigned uType;
608  is >> uType;
609  type = static_cast<InsertionBoundary::Distribution>(uType);
610  return is;
611 }
std::vector< Mdouble > variableCumulativeVolumeFlowRate_
void copyDataFromParticleToMPIParticle(BaseParticle *p)
bool isManuallyInserting_
A flag to enable a top-down class-by-class manual insertion of a PSD; default is FALSE.
static MPIContainer & Instance()
fetch the instance to be used for communication
Definition: MpiContainer.h:130
The DPMBase header includes quite a few header files, defining all the handlers, which are essential...
Definition: DPMBase.h:72
double getMassOfParticlesInserted() const
Logger< MERCURY_LOGLEVEL > logger("MercuryKernel")
double Mdouble
Definition: GeneralDefine.h:34
virtual void setRadius(Mdouble radius)
Sets the particle's radius_ (and adjusts the mass_ accordingly, based on the particle's species) ...
InsertionBoundary()
Default constructor: set everything to 0/nullptr.
void setPSD(PSD psd)
Sets the range of particle radii that may be generated from a user defined PSD.
Mdouble getVolumeFlowRate() const
const std::complex< Mdouble > i
Definition: ExtendedMath.h:50
unsigned int getMaxFailed() const
Gets the number of times that the boundary may fail to insert a particle.
Mdouble insertManuallyByVolume(Mdouble volume)
Draw sample radius manually per size class and check the volumeAllowed of each size class to insert t...
Definition: PSD.cc:137
bool insertParticle(Mdouble time)
Mdouble radMin_
Minimum and maximum radii of the generated particles.
BaseParticle * particleToCopy_
Particle that will be inserted through the insertion boundary.
Distribution
Defines a custom particle size distribution; distribution_ will always be used, unless particleSizeDi...
void setMaxFailed(unsigned int maxFailed)
Sets the number of times that the wall may fail to insert a particle.
const std::vector< RadiusAndProbability > getParticleSizeDistribution() const
Get the PSD vector.
Definition: PSD.cc:718
void setSpecies(const ParticleSpecies *species)
Mdouble getRandomNumber()
This is a random generating routine can be used for initial positions.
Definition: RNG.cc:143
void read(std::istream &is) override=0
Reads the object's id_ from given istream NB: purely virtual function, overriding the version of Base...
Definition: BaseBoundary.cc:61
virtual BaseParticle * copy() const =0
Particle copy method. It calls to copy constructor of this Particle, useful for polymorfism.
unsigned int maxFailed_
Number of times that the wall may fail to insert a particle.
Distribution distribution_
defines a custom particle size distribution, which by default is uniform
Boundary structure for boundaries used for insertion of particles.
const Mdouble inf
Definition: GeneralDefine.h:44
void setParticleToCopy(BaseParticle *particleToCopy)
Sets the particle that will be inserted through the insertion boundary.
std::istream & operator>>(std::istream &is, InsertionBoundary::Distribution &type)
void copyDataFromMPIParticleToParticle(MPIParticle *bP, BaseParticle *p, ParticleHandler *particleHandler)
Copies data from an MPIParticle class to a BaseParticle and sets the particleHandler and species...
This is a class that generates random numbers i.e. named the Random Number Generator (RNG)...
Definition: RNG.h:52
virtual void placeParticle(BaseParticle *p, RNG &random)=0
Purely virtual function that generates the extrinsic properties (position, velocity) of a particle...
Mdouble drawSample()
Draw a sample radius from a CUMULATIVE_NUMBER_DISTRIBUTION.
Definition: PSD.cc:109
virtual bool checkParticleForInteraction(const BaseParticle &P)
Checks whether a particle P has any interaction with walls or other particles.
Definition: DPMBase.cc:4588
double getVolumeOfParticlesInserted() const
Mdouble volumeInserted_
Total volume of particles inserted.
void checkBoundaryBeforeTimeStep(DPMBase *md) override
Fills the boundary with particles.
unsigned int getNumberOfParticlesInserted() const
Gets the number of particles inserted by the boundary.
ParticleHandler particleHandler
An object of the class ParticleHandler, contains the pointers to all the particles created...
Definition: DPMBase.h:1329
T * getObject(const unsigned int id)
Gets a pointer to the Object at the specified index in the BaseHandler.
Definition: BaseHandler.h:613
std::ostream & operator<<(std::ostream &os, InsertionBoundary::Distribution type)
Data class to send a particle over MPI.
Definition: MpiDataClass.h:81
BaseParticle * getParticleToCopy() const
Gets the particle that will be inserted through the insertion boundary.
void set(BaseParticle *particleToCopy, unsigned int maxFailed)
Sets the particle that will be inserted and the maximum number of times for which insertion may fail...
#define PROCESSOR_ID
Definition: GeneralDefine.h:63
std::enable_if<!std::is_pointer< U >::value, U * >::type copyAndAddObject(const U &object)
Creates a copy of a Object and adds it to the BaseHandler.
Definition: BaseHandler.h:379
void write(std::ostream &os) const override=0
Adds object's id_ to given ostream NB: purely virtual function, overriding the version of BaseObject...
Definition: BaseBoundary.cc:70
BoundaryHandler * getHandler() const
Returns the boundary's BoundaryHandler.
void write(std::ostream &os) const override
Writes the boundary's id_ and maxFailed_.
void setVolumeFlowRate(Mdouble volumeFlowRate_)
bool readOptionalVariable(std::istream &is, const std::string &name, T &variable)
Reads optional variables in the restart file.
Definition: Helpers.h:247
void activate()
Turns on the InsertionBoundary.
void setManualInsertion(bool manualInsertion)
Set the flag for a manual PSD insertion routine.
#define NUMBER_OF_PROCESSORS
For the MPI communication routines this quantity is often required. defining this macro makes the cod...
Definition: GeneralDefine.h:62
~InsertionBoundary() override
Destructor: delete the particle that has to be copied at every insertion.
RNG random
This is a random generator, often used for setting up the initial conditions etc...
Definition: DPMBase.h:1324
std::enable_if< std::is_scalar< T >::value, void >::type broadcast(T &t, int fromProcessor=0)
Broadcasts a scalar from the root to all other processors.
Definition: MpiContainer.h:429
virtual BaseParticle * generateParticle(RNG &random)
Virtual function that generates the intrinsic properties (species, radius) of one particle...
void deactivate()
Turns off the InsertionBoundary.
Contains a vector with radii and probabilities of a user defined particle size distribution (PSD) ...
Definition: PSD.h:62
Mdouble getNextTime() const
Returns the current simulation time.
Definition: DPMBase.cc:805
void setInitialVolume(Mdouble initialVolume)
unsigned int getIndSpecies() const
Returns the index of the species associated with the interactable object.
unsigned int numberOfParticlesInserted_
Number of particles that are already inserted.
Mdouble getInitialVolume() const
DPMBase * getDPMBase()
Gets the problem that is solved using this handler.
Definition: BaseHandler.h:725
Distribution getDistribution()
Gets the range of particle radii that may be generated.
void setDistribution(Distribution distribution)
Sets the range of particle radii that may be generated to custom distributions.
Mdouble massInserted_
Total mass of particles inserted.
void setVariableVolumeFlowRate(const std::vector< Mdouble > &variableCumulativeVolumeFlowRate, Mdouble samplingInterval)
bool isActivated()
Returns whether the InsertionBoundary is activated.
BaseParticle * readAndCreateObject(std::istream &is)
Create a new particle, based on the information provided in a restart file.
void read(std::istream &is) override
Reads the boundary's id_ and maxFailed_.
PSD particleSizeDistribution_
Defines a particle size distribution as an object of the PSD class; if particleSizeDistribution_ is e...
bool isActivated_
The InsertionBoundary is activated by default. If the InsertionBoundary is deactivated, then it introduces no particles (useful for trying to maintain a certain insertion rate).