MercuryDPM  Trunk
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
ChuteBottom.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 
27 #include <algorithm>
28 #include "ChuteBottom.h"
32 #include "Walls/InfiniteWall.h"
33 #include "Logger.h"
34 
39 {
40  constructor();
41 }
42 
51  : DPMBase(other), Chute(other)
52 {
53  constructor();
54 }
55 
64  : DPMBase(other), Chute(other)
65 {
66  constructor();
67 }
68 
77  : DPMBase(other), Chute(other)
78 {
79  constructor();
80 }
81 
90  : DPMBase(other), Chute(other)
91 {
92  constructor();
93 }
94 
101  : DPMBase(other), Chute(other),
102  thickness_(other.thickness_), isBottomPeriodic_(other.isBottomPeriodic_)
103 {
104 
105 }
106 
111 {
112  setName("roughbottom");
113  fStatFile.setFileType(FileType::NO_FILE); //set to 0 for no data creation
118  setThickness(2.4);
119  setIsBottomPeriodic(true);
120 }
121 
140 {
141  // set up mini-simulation with particle inflow on the left end and horizontal chute
142  // set all parameters that should be different from the original chute
143  setChuteAngle(0.0);
145  //~ setInflowHeight(45.*getInflowParticleRadius());
146  // note: Changing the Inflow height was an attempt to make the bottom density homogeneous, but it did not have the desired effect
149 
150  auto species = dynamic_cast<LinearViscoelasticNormalSpecies*>(speciesHandler.getObject(0));
151  if (species != nullptr)
152  {
155  //this number is chosen to be consistent with the old implementation (with contained a bug and was replaced); originally, a factor of 10.0 (and a restitution of 0.2) was chosen.
156  const Mdouble collisionTime = 13.2863 * species->getCollisionTime(mass);
157  species->setCollisionTimeAndRestitutionCoefficient(collisionTime, 0.343008, mass);
158  setTimeStep(0.02 * 9.41823 * collisionTime);
159  logger(INFO, "Time step: %", getTimeStep());
160  //logger(INFO, "Species: %", *species);
161  }
162  else
163  {
165  logger(WARN, "[ChuteBottom::makeRoughBottom()] species type does not allow setting the parameters.");
166  }
167 
168  auto species2 = dynamic_cast<SlidingFrictionSpecies*>(speciesHandler.getObject(0));
169  if (species2 != nullptr)
170  species2->setSlidingFrictionCoefficient(0);
171 
172  // set the simulation to run for 2000 time steps
173  setTimeMax(getTimeStep() * 2000);
174  setSaveCount(100);
175 
176  // run the simulation
177  solve();
178 
179  //Find the Z-position of the highest particle in the system
180  Mdouble height = 0;
181  for (BaseParticle* const p : particleHandler)
182  {
183  height = std::max(height, p->getPosition().Z);
184  }
185 
186  // Next, all particles are removed from the system, except those with a Z-position s.t.:
187  // hmax = height - maxInflowParticleRadius_;
188  // hmax - thickness_ * maxInflowParticleRadius_ < Z-position < hmax
189  // note that hmax = height - rMax
190  logger(INFO, "[ChuteBottom::makeRoughBottom()] Thickness: %", thickness_);
191 
192  for (std::vector<BaseParticle*>::iterator it = particleHandler.begin(); it != particleHandler.end(); ++it)
193  {
194  if ((*it)->getPosition().Z < height - (1.0 + thickness_) * getMaxInflowParticleRadius() ||
195  (*it)->getPosition().Z > height - getMaxInflowParticleRadius())
196  {
197  //delete particles outside the given range
198  particleHandler.removeObject((*it)->getIndex());
199  --it;
200  }
201  }
202  //The remaining particles must be moved downwards and fixed.
203  // They also must have the species as specified by the chute.
204  for (BaseParticle* p : particleHandler)
205  {
206  p->move(Vec3D(0.0, 0.0, getMaxInflowParticleRadius() - height));
207  p->fixParticle();
208  p->setSpecies(chute.speciesHandler.getObject(0));
209  }
210 
211  //copy the rough bottom over
212  logger(INFO, "[ChuteBottom::makeRoughBottom()] Chute bottom finished, consisting of % particles",
213  particleHandler.getNumberOfObjects());
215 
216 }
217 
223 {
224 
225  particleHandler.setStorageCapacity(static_cast<unsigned int>(std::min(
226  3.0 * getXMax() * getYMax() * getZMax() / mathsFunc::cubic(2.0 * getInflowParticleRadius()), 1e6)));
227 
228  createBottom();
229 
235  wallHandler.clear();
237  if (isBottomPeriodic_)
238  {
239  InfiniteWall w0;
241  w0.set(Vec3D(0.0, 0.0, -1.0), Vec3D(0, 0, getZMin() - getInflowParticleRadius()));
243  PeriodicBoundary b0;
244  b0.set(Vec3D(1.0, 0.0, 0.0), getXMin(), getXMax());
246  b0.set(Vec3D(0.0, 1.0, 0.0), getYMin(), getYMax());
248  }
249  else
250  {
251  InfiniteWall w0;
253  w0.set(Vec3D(0.0, 0.0, -1.0), Vec3D(0, 0, getZMin() - getInflowParticleRadius()));
255  w0.set(Vec3D(-1.0, 0.0, 0.0), Vec3D(getXMin(), 0, 0));
257  w0.set(Vec3D(1.0, 0.0, 0.0), Vec3D(getXMax(), 0, 0));
259  w0.set(Vec3D(0.0, -1.0, 0.0), Vec3D(0, getYMin(), 0));
261  w0.set(Vec3D(0.0, 1.0, 0.0), Vec3D(0, getYMax(), 0));
263  }
264 
265  // add particles
272  unsigned int failed = 0;
273  const unsigned int max_failed = 500;
274 
275  SphericalParticle inflowParticle_;
276  inflowParticle_.setSpecies(speciesHandler.getObject(0));
277  inflowParticle_.setHandler(&particleHandler);
278  inflowParticle_.setOrientation({1, 0, 0, 0});
279  inflowParticle_.setAngularVelocity(Vec3D(0.0, 0.0, 0.0));
280 
281  // try max_failed times to find new insertable particle
282  // every time an insertion succeeds the counter failed is set to 0 again.
283  while (failed <= max_failed)
284  {
285  inflowParticle_.setRadius(getFixedParticleRadius());
286  //inflowParticle_.computeMass();
287 
288  // The position components are first stored in a Vec3D, because if you pass
289  // them directly into setPosition the compiler is allowed to change the order
290  // in which the numbers are generated
291  Vec3D position;
292  position.X = random.getRandomNumber(inflowParticle_.getRadius(), getXMax() - inflowParticle_.getRadius());
293  position.Y = random.getRandomNumber(inflowParticle_.getRadius(), getYMax() - inflowParticle_.getRadius());
294  position.Z = random.getRandomNumber(2 * inflowParticle_.getRadius(), getZMax() - inflowParticle_.getRadius());
295  inflowParticle_.setPosition(position);
296  inflowParticle_.setVelocity(Vec3D(0.0, 0.0, 0.0));
297 
298  //if the volume we want to insert the particle is still free, insert the particle
299  if (checkParticleForInteraction(inflowParticle_))
300  {
301  particleHandler.copyAndAddObject(inflowParticle_);
302  failed = 0;
303  }
304  else
305  {
306  failed++;
307  }
308  }
309  //set_Nmax(particleHandler.getNumberOfObjects());
310  logger(INFO, "[ChuteBottom::setupInitialConditions()] Number of particles created: %",
312 
313  //fix hgrid (there is still an issue when particles are polydispersed)
314  //assume 1-2 levels are optimal (which is the case for mono and bidispersed) and set the cell size to min and max
315  // !this is not optimal for polydispersed
316  Mdouble minCell = 2. * std::min(getFixedParticleRadius(), getMinInflowParticleRadius());
317  Mdouble maxCell = 2. * std::max(getFixedParticleRadius(), getMaxInflowParticleRadius());
318  if ((minCell == maxCell) | (minCell == 0.))
320  else
322  //set_HGRID_cell_to_cell_ratio (1.0000000001*maxCell/minCell);
323  //optimize number of buckets
324  //setHGridNumberOfBucketsToPower(particleHandler.getNumberOfObjects() * 1.5);
325  //end: fix hgrid
326 
327  //~ write(std::cout,false);
328 }
329 
335 {
336 }
337 
344 {
345  return thickness_;
346 }
347 
354 {
355  logger.assert_always(new_ > 0.0, "[ChuteBottom::setThickness()] thickness % not positive.", new_);
356  thickness_ = new_;
357 }
358 
364 {
365  return isBottomPeriodic_;
366 }
367 
372 void ChuteBottom::setIsBottomPeriodic(bool isBottomPeriodic)
373 {
374  isBottomPeriodic_ = isBottomPeriodic;
375 }
void setInflowHeight(Mdouble inflowHeight)
Sets maximum inflow height (Z-direction)
Definition: Chute.cc:892
void set(Vec3D normal, Mdouble distanceLeft, Mdouble distanceRight)
Defines a PeriodicBoundary by its normal and positions.
void solve()
The work horse of the code.
Definition: DPMBase.cc:3968
The DPMBase header includes quite a few header files, defining all the handlers, which are essential...
Definition: DPMBase.h:72
Mdouble X
the vector components
Definition: Vector.h:65
ChuteBottom()
This is the default constructor. All it does is set sensible defaults.
Definition: ChuteBottom.cc:38
A basic particle.
void setVelocity(const Vec3D &velocity)
set the velocity of the BaseInteractable.
void setTimeMax(Mdouble newTMax)
Sets a new value for the maximum simulation duration.
Definition: DPMBase.cc:840
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) ...
unsigned int getNumberOfObjects() const override
Returns the number of objects in the container. In parallel code this practice is forbidden to avoid ...
void setupInitialConditions() override
Sets up initial conditions before running a chute simulation.
Definition: ChuteBottom.cc:222
Mdouble getZMax() const
If the length of the problem domain in z-direction is ZMax - ZMin, then getZMax() returns ZMax...
Definition: DPMBase.h:617
Mdouble getXMin() const
If the length of the problem domain in x-direction is XMax - XMin, then getXMin() returns XMin...
Definition: DPMBase.h:586
void setStorageCapacity(const unsigned int N)
Sets the storage capacity of this BaseHandler.
Definition: BaseHandler.h:669
Used by Chute::createBottom to create an unordered particle layer.
Definition: ChuteBottom.h:39
void setHandler(ParticleHandler *handler)
Sets the pointer to the particle's ParticleHandler.
Mdouble getYMin() const
If the length of the problem domain in y-direction is YMax - YMin, then getYMin() returns YMin...
Definition: DPMBase.h:599
File interactionFile
File class to handle in- and output into .interactions file. This file hold information about interac...
Definition: DPMBase.h:1396
void setChuteAngle(Mdouble chuteAngle)
Sets gravity vector according to chute angle (in degrees)
Definition: Chute.cc:703
void setSpecies(const ParticleSpecies *species)
Mdouble getRandomNumber()
This is a random generating routine can be used for initial positions.
Definition: RNG.cc:143
Mdouble getMassFromRadius(Mdouble radius) const
T cubic(const T val)
calculates the cube of a number
Definition: ExtendedMath.h:113
void setRoughBottomType(RoughBottomType roughBottomType)
Sets the type of rough bottom of the chute.
Definition: Chute.cc:649
Defines a pair of periodic walls. Inherits from BaseBoundary.
Mdouble getXMax() const
If the length of the problem domain in x-direction is XMax - XMin, then getXMax() returns XMax...
Definition: DPMBase.h:593
Mdouble getInflowParticleRadius() const
Returns the average radius of inflow particles.
Definition: Chute.cc:864
void constructor()
This is the actual constructor METHOD; it is called by all constructors above (except the default cop...
Definition: ChuteBottom.cc:110
Mdouble thickness_
Thickness of the multilayer chute rough bottom. See also documentation of ChuteBottom::makeRoughBotto...
Definition: ChuteBottom.h:123
bool checkParticleForInteraction(const BaseParticle &P) final
Checks if given BaseParticle has an interaction with a BaseWall or other BaseParticle.
Definition: MercuryBase.cc:588
Creates chutes with different bottoms. Inherits from Mercury3D (-> MercuryBase -> DPMBase)...
Definition: Chute.h:64
This is the base class for both Mercury2D and Mercury3D. Note the actually abstract grid is defined i...
Definition: MercuryBase.h:125
Mdouble getMaxInflowParticleRadius() const
Returns the maximum radius of inflow particles.
Definition: Chute.cc:882
void makeRoughBottom(Chute &chute)
Makes a multilayered rough bottom with thickness thickness_.
Definition: ChuteBottom.cc:139
File dataFile
An instance of class File to handle in- and output into a .data file.
Definition: DPMBase.h:1370
file will not be created/read
LinearViscoelasticNormalSpecies contains the parameters used to describe a linear elastic-dissipative...
File fStatFile
An instance of class File to handle in- and output into a .fstat file.
Definition: DPMBase.h:1375
void setIsBottomPeriodic(bool isBottomPeriodic)
Sets whether the bottom should be periodic in Y.
Definition: ChuteBottom.cc:372
BoundaryHandler boundaryHandler
An object of the class BoundaryHandler which concerns insertion and deletion of particles into or fro...
Definition: DPMBase.h:1344
Mdouble getFixedParticleRadius() const
Returns the particle radius of the fixed particles which constitute the (rough) chute bottom...
Definition: Chute.cc:634
Mdouble getThickness() const
Returns the thickness of the multilayer rough bottom.
Definition: ChuteBottom.cc:343
This adds on the hierarchical grid code for 3D problems.
Definition: Mercury3D.h:36
ParticleHandler particleHandler
An object of the class ParticleHandler, contains the pointers to all the particles created...
Definition: DPMBase.h:1329
all data will be written into/ read from a single file called name_
T * getObject(const unsigned int id)
Gets a pointer to the Object at the specified index in the BaseHandler.
Definition: BaseHandler.h:613
void setSaveCount(unsigned int saveCount)
Sets File::saveCount_ for all files (ene, data, fstat, restart, stat)
Definition: DPMBase.cc:398
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 setSlidingFrictionCoefficient(Mdouble new_mu)
Allows the (dynamic) Coulomb friction coefficient to be changed; also sets mu_s by default...
Mdouble getRadius() const
Returns the particle's radius.
Definition: BaseParticle.h:345
void setHGridMaxLevels(unsigned int HGridMaxLevels)
Sets the maximum number of levels of the HGrid in this MercuryBase.
Definition: MercuryBase.cc:470
SpeciesHandler speciesHandler
A handler to that stores the species type i.e. LinearViscoelasticSpecies, etc.
Definition: DPMBase.h:1319
void hGridActionsBeforeTimeLoop() override
This sets up the broad phase information, has to be done at this stage because it requires the partic...
Definition: MercuryBase.cc:94
Mdouble getYMax() const
If the length of the problem domain in y-direction is YMax - YMin, then getYMax() returns XMax...
Definition: DPMBase.h:605
Mdouble Y
Definition: Vector.h:65
SlidingFrictionSpecies contains the parameters used to describe sliding friction. ...
void setFileType(FileType fileType)
Sets the type of file needed to write into or read from. File::fileType_.
Definition: File.cc:216
WallHandler wallHandler
An object of the class WallHandler. Contains pointers to all the walls created.
Definition: DPMBase.h:1339
void setFixedParticleRadius(Mdouble fixedParticleRadius)
Sets the particle radius of the fixed particles which constitute the (rough) chute bottom...
Definition: Chute.cc:616
RNG random
This is a random generator, often used for setting up the initial conditions etc...
Definition: DPMBase.h:1324
void setPosition(const Vec3D &position)
Sets the position of this BaseInteractable.
virtual void createBottom()
Creates the chute bottom, which can be either flat or one of three flavours of rough.
Definition: Chute.cc:331
void setName(const std::string &name)
Allows to set the name of all the files (ene, data, fstat, restart, stat)
Definition: DPMBase.cc:412
Mdouble getMinInflowParticleRadius() const
returns the minimum radius of inflow particles
Definition: Chute.cc:873
void set(Vec3D normal, Vec3D point)
Defines a standard wall, given an outward normal vector s.t. normal*x=normal*point for all x of the w...
Mdouble getZMin() const
If the length of the problem domain in z-direction is ZMax - ZMin, then getZMin() returns ZMin...
Definition: DPMBase.h:611
void setTimeStep(Mdouble newDt)
Sets a new value for the simulation time step.
Definition: DPMBase.cc:1195
File eneFile
An instance of class File to handle in- and output into a .ene file.
Definition: DPMBase.h:1380
This is a class defining walls.
Definition: InfiniteWall.h:47
File restartFile
An instance of class File to handle in- and output into a .restart file.
Definition: DPMBase.h:1385
void setOrientation(const Quaternion &orientation)
Sets the orientation of this BaseInteractable.
void setThickness(Mdouble thickness)
Sets the thickness of the multilayer rough bottom.
Definition: ChuteBottom.cc:353
Definition: Vector.h:49
Mdouble getTimeStep() const
Returns the simulation time step.
Definition: DPMBase.cc:1211
Mdouble Z
Definition: Vector.h:65
void setAngularVelocity(const Vec3D &angularVelocity)
set the angular velocity of the BaseInteractble.
bool isBottomPeriodic_
TRUE if the bottom is periodic in Y.
Definition: ChuteBottom.h:128
virtual void clear()
Empties the whole BaseHandler by removing all Objects and setting all other variables to 0...
Definition: BaseHandler.h:528
void actionsBeforeTimeStep() override
Performs all necessary actions before the start of a time step (none in this case) ...
Definition: ChuteBottom.cc:334
void setSpecies(const ParticleSpecies *species)
Defines the species of the current wall.
Definition: BaseWall.cc:171
bool getIsBottomPeriodic() const
Returns TRUE if the bottom is periodic in Y.
Definition: ChuteBottom.cc:363