1 /* ----------------------------------------------------------------------
2 This is the
3
4 ██╗ ██╗ ██████╗ ██████╗ ██████╗ ██╗ ██╗████████╗███████╗
5 ██║ ██║██╔════╝ ██╔════╝ ██╔════╝ ██║ ██║╚══██╔══╝██╔════╝
6 ██║ ██║██║ ███╗██║ ███╗██║ ███╗███████║ ██║ ███████╗
7 ██║ ██║██║ ██║██║ ██║██║ ██║██╔══██║ ██║ ╚════██║
8 ███████╗██║╚██████╔╝╚██████╔╝╚██████╔╝██║ ██║ ██║ ███████║
9 ╚══════╝╚═╝ ╚═════╝ ╚═════╝ ╚═════╝ ╚═╝ ╚═╝ ╚═╝ ╚══════╝®
10
11 DEM simulation engine, released by
12 DCS Computing Gmbh, Linz, Austria
13 http://www.dcs-computing.com, office@dcs-computing.com
14
15 LIGGGHTS® is part of CFDEM®project:
16 http://www.liggghts.com | http://www.cfdem.com
17
18 Core developer and main author:
19 Christoph Kloss, christoph.kloss@dcs-computing.com
20
21 LIGGGHTS® is open-source, distributed under the terms of the GNU Public
22 License, version 2 or later. It is distributed in the hope that it will
23 be useful, but WITHOUT ANY WARRANTY; without even the implied warranty
24 of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. You should have
25 received a copy of the GNU General Public License along with LIGGGHTS®.
26 If not, see http://www.gnu.org/licenses . See also top-level README
27 and LICENSE files.
28
29 LIGGGHTS® and CFDEM® are registered trade marks of DCS Computing GmbH,
30 the producer of the LIGGGHTS® software and the CFDEM®coupling software
31 See http://www.cfdem.com/terms-trademark-policy for details.
32
33 -------------------------------------------------------------------------
34 Contributing author and copyright for this file:
35 Christoph Kloss (JKU Linz, DCS Computing GmbH, Linz)
36 Richard Berger (JKU Linz)
37
38 Copyright 2012- DCS Computing GmbH, Linz
39 Copyright 2009-2015 JKU Linz
40 ------------------------------------------------------------------------- */
41
42 #include <cmath>
43 #include <algorithm>
44 #include <stdlib.h>
45 #include <string.h>
46 #include "atom.h"
47 #include "atom_vec.h"
48 #include "force.h"
49 #include "update.h"
50 #include "comm.h"
51 #include "modify.h"
52 #include "domain.h"
53 #include "random_park.h"
54 #include "memory.h"
55 #include "error.h"
56 #include "fix_multisphere.h"
57 #include "fix_particledistribution_discrete.h"
58 #include "fix_template_sphere.h"
59 #include "fix_property_atom.h"
60 #include "irregular.h"
61 #include "fix_insert.h"
62 #include "math_extra_liggghts.h"
63 #include "mpi_liggghts.h"
64 #include "vector_liggghts.h"
65
66 #include "probability_distribution.h"
67 #include "region_neighbor_list.h"
68
69 using namespace LAMMPS_NS;
70 using namespace FixConst;
71
72 #define EPSILON 0.001
73
74 #define LMP_DEBUGMODE_FIXINSERT false //(667 == update->ntimestep)// true
75 #define LMP_DEBUG_OUT_FIXINSERT screen
76
77 /* ---------------------------------------------------------------------- */
78
FixInsert(LAMMPS * lmp,int narg,char ** arg)79 FixInsert::FixInsert(LAMMPS *lmp, int narg, char **arg) :
80 Fix(lmp, narg, arg),
81 neighList(*new RegionNeighborList<interpolate_no>(lmp))
82 {
83 if (narg < 7) error->fix_error(FLERR,this,"not enough arguments");
84
85 restart_global = 1;
86
87 setup_flag = false;
88
89 fix_distribution = NULL;
90 fix_multisphere = NULL;
91 multisphere = NULL;
92
93 compress_flag = false ;
94
95 // required args
96 iarg = 3;
97
98 if(strcmp(arg[iarg++],"seed")) error->fix_error(FLERR,this,"expecting keyword 'seed'");
99 // random number generator, seed depends on proc
100 random = new RanPark(lmp, arg[iarg++], true);
101 seed = random->getSeed();
102 if (seed <= 0) error->fix_error(FLERR,this,"illegal seed");
103
104 // set defaults
105 init_defaults();
106
107 // parse args
108
109 #ifdef SUPERQUADRIC_ACTIVE_FLAG
110 check_obb_flag = 1;
111 #endif
112
113 bool hasargs = true;
114 while(iarg < narg && hasargs)
115 {
116 hasargs = false;
117
118 if(strcmp(arg[iarg],"distributiontemplate") == 0) {
119 if (iarg+2 > narg) error->fix_error(FLERR,this,"");
120 int ifix = modify->find_fix(arg[iarg+1]);
121 if(ifix < 0 || strncmp(modify->fix[ifix]->style,"particledistribution/discrete",29))
122 error->fix_error(FLERR,this,"Fix insert requires you to define a valid ID for a fix of type particledistribution/discrete");
123 fix_distribution = static_cast<FixParticledistributionDiscrete*>(modify->fix[ifix]);
124 iarg += 2;
125 hasargs = true;
126 } else if (strcmp(arg[iarg],"maxattempt") == 0) {
127 if (iarg+2 > narg) error->fix_error(FLERR,this,"");
128 maxattempt = atoi(arg[iarg+1]);
129 iarg += 2;
130 hasargs = true;
131 } else if (strcmp(arg[iarg],"nparticles") == 0) {
132 if (iarg+2 > narg) error->fix_error(FLERR,this,"");
133 if(strcmp(arg[iarg+1],"INF") == 0)
134 ninsert_exists = 0;
135 else ninsert = atof(arg[iarg+1]);
136 iarg += 2;
137 hasargs = true;
138 } else if (strcmp(arg[iarg],"mass") == 0) {
139 if (iarg+2 > narg) error->fix_error(FLERR,this,"");
140 if(strcmp(arg[iarg+1],"INF") == 0)
141 ninsert_exists = 0;
142 else massinsert = atof(arg[iarg+1]);
143 iarg += 2;
144 hasargs = true;
145 } else if (strcmp(arg[iarg],"massrate") == 0) {
146 if (iarg+2 > narg) error->fix_error(FLERR,this,"");
147 massflowrate = atof(arg[iarg+1]);
148 iarg += 2;
149 hasargs = true;
150 } else if (strcmp(arg[iarg],"particlerate") == 0) {
151 if (iarg+2 > narg) error->fix_error(FLERR,this,"");
152 nflowrate = atof(arg[iarg+1]);
153 iarg += 2;
154 hasargs = true;
155 } else if (strcmp(arg[iarg],"insert_every_time") == 0 || strcmp(arg[iarg],"insert_every") == 0 || strcmp(arg[iarg],"every") == 0) {
156 if (iarg+2 > narg) error->fix_error(FLERR,this,"");
157 if(strcmp(arg[iarg+1],"once") == 0) insert_every = 0;
158 else if(strcmp(arg[iarg],"insert_every_time") == 0)
159 {
160 if(!update->timestep_set)
161 error->fix_error(FLERR,this,"need so set 'timestep' before");
162 insert_every = static_cast<int>(atof(arg[iarg+1])/update->dt);
163 }
164 else
165 insert_every = atoi(arg[iarg+1]);
166 if(insert_every < 0) error->fix_error(FLERR,this,"insert_every must be >= 0");
167 iarg += 2;
168 hasargs = true;
169 } else if (strcmp(arg[iarg],"start") == 0) {
170 if (iarg+2 > narg) error->fix_error(FLERR,this,"");
171 first_ins_step = atoi(arg[iarg+1]);
172 if(first_ins_step < update->ntimestep + 1 && !modify->fix_restart_in_progress())
173 error->fix_error(FLERR,this,"'start' step can not be before current step");
174 iarg += 2;
175 hasargs = true;
176 } else if (strcmp(arg[iarg],"overlapcheck") == 0) {
177 if (iarg+2 > narg) error->fix_error(FLERR,this,"");
178 if(strcmp(arg[iarg+1],"yes")==0) check_ol_flag = 1;
179 else if(strcmp(arg[iarg+1],"no")==0) check_ol_flag = 0;
180 else error->fix_error(FLERR,this,"");
181 iarg += 2;
182 hasargs = true;
183 } else if (strcmp(arg[iarg],"all_in") == 0) {
184 if (iarg+2 > narg) error->fix_error(FLERR,this,"");
185 if(strcmp(arg[iarg+1],"yes")==0) all_in_flag = 1;
186 else if(strcmp(arg[iarg+1],"no")==0) all_in_flag = 0;
187 else error->fix_error(FLERR,this,"");
188 iarg += 2;
189 hasargs = true;
190 } else if (strcmp(arg[iarg],"set_property") == 0) {
191 if (iarg+3 > narg) error->fix_error(FLERR,this,"");
192 int n = strlen(arg[iarg+1]) + 1;
193 property_name = new char[n];
194 strcpy(property_name,arg[iarg+1]);
195 fix_property_value = force->numeric(FLERR,arg[iarg+2]);
196 iarg += 3;
197 hasargs = true;
198 } else if (strcmp(arg[iarg],"random_distribute") == 0) {
199 if (iarg+2 > narg) error->fix_error(FLERR,this,"");
200 if(strcmp(arg[iarg+1],"uncorrelated")==0) exact_number = 0;
201 else if(strcmp(arg[iarg+1],"exact")==0) exact_number = 1;
202 else error->fix_error(FLERR,this,"");
203 iarg += 2;
204 hasargs = true;
205 } else if (strcmp(arg[iarg],"verbose") == 0) {
206 if (iarg+2 > narg) error->fix_error(FLERR,this,"");
207 if(strcmp(arg[iarg+1],"no")==0) print_stats_during_flag = 0;
208 else if(strcmp(arg[iarg+1],"yes")==0) print_stats_during_flag = 1;
209 else error->fix_error(FLERR,this,"");
210 iarg += 2;
211 hasargs = true;
212 } else if (strcmp(arg[iarg],"compress_tags") == 0) {
213 if (iarg+2 > narg) error->fix_error(FLERR,this,"not enough arguments for compress_tags");
214 if(0 == strcmp(arg[iarg+1],"yes"))
215 compress_flag = true;
216 else if(0 == strcmp(arg[iarg+1],"no"))
217 compress_flag = false;
218 else
219 error->fix_error(FLERR,this,"expecting 'yes' or 'no' after 'compress_tags'");
220 iarg += 2;
221 hasargs = true;
222 } else if (strcmp(arg[iarg],"vel") == 0) {
223 if (iarg+5 > narg) error->fix_error(FLERR,this,"not enough keyword for 'vel'");
224 if (strcmp(arg[iarg+1],"constant") == 0) {
225 v_insert[0] = atof(arg[iarg+2]);
226 v_insert[1] = atof(arg[iarg+3]);
227 v_insert[2] = atof(arg[iarg+4]);
228 iarg += 5;
229 } else if (strcmp(arg[iarg+1],"uniform") == 0) {
230 if (iarg+8 > narg) error->fix_error(FLERR,this,"not enough keyword for 'uniform'");
231 v_randomSetting = RANDOM_UNIFORM;
232 v_insert[0] = atof(arg[iarg+2]);
233 v_insert[1] = atof(arg[iarg+3]);
234 v_insert[2] = atof(arg[iarg+4]);
235 v_insertFluct[0] = atof(arg[iarg+5]);
236 v_insertFluct[1] = atof(arg[iarg+6]);
237 v_insertFluct[2] = atof(arg[iarg+7]);
238 iarg += 8;
239 } else if (strcmp(arg[iarg+1],"gaussian") == 0) {
240 if (iarg+8 > narg) error->fix_error(FLERR,this,"not enough keyword for 'gaussian'");
241 v_randomSetting = RANDOM_GAUSSIAN;
242 v_insert[0] = atof(arg[iarg+2]);
243 v_insert[1] = atof(arg[iarg+3]);
244 v_insert[2] = atof(arg[iarg+4]);
245 v_insertFluct[0] = atof(arg[iarg+5]);
246 v_insertFluct[1] = atof(arg[iarg+6]);
247 v_insertFluct[2] = atof(arg[iarg+7]);
248 iarg += 8;
249 } else
250 error->fix_error(FLERR,this,"expecting keyword 'constant' or 'uniform' or 'gaussian' after keyword 'vel'");
251 hasargs = true;
252 } else if (strcmp(arg[iarg],"omega") == 0) {
253 if (iarg+5 > narg) error->fix_error(FLERR,this,"");
254 if (strcmp(arg[iarg+1],"constant") == 0)
255 {
256 omega_insert[0] = atof(arg[iarg+2]);
257 omega_insert[1] = atof(arg[iarg+3]);
258 omega_insert[2] = atof(arg[iarg+4]);
259 } else error->fix_error(FLERR,this,"expecting keyword 'constant' after keyword 'omega'");
260 iarg += 5;
261 hasargs = true;
262 } else if (strcmp(arg[iarg],"orientation") == 0) {
263 if (iarg+2 > narg)
264 error->fix_error(FLERR,this,"not enough arguments for 'orientation'");
265 iarg++;
266 if(strcmp(arg[iarg],"random") == 0)
267 {
268 quat_random_ = true;
269 iarg++;
270 }
271 else if(strcmp(arg[iarg],"template") == 0)
272 {
273 quat_random_ = false;
274 iarg++;
275 }
276 else if (strcmp(arg[iarg],"constant") == 0)
277 {
278 iarg++;
279 if (iarg+4 > narg) error->fix_error(FLERR,this,"");
280 quat_insert[0] = atof(arg[iarg++]);
281 quat_insert[1] = atof(arg[iarg++]);
282 quat_insert[2] = atof(arg[iarg++]);
283 quat_insert[3] = atof(arg[iarg++]);
284 } else error->fix_error(FLERR,this,"expecting 'random', template' or 'constant' after keyword 'quat'");
285 hasargs = true;
286 }
287
288 #ifdef SUPERQUADRIC_ACTIVE_FLAG
289 else if (strcmp(arg[iarg],"check_obb") == 0) {
290 if (iarg+2 > narg) error->fix_error(FLERR,this,"");
291 if(strcmp(arg[iarg+1],"yes")==0) check_obb_flag = 1;
292 else if(strcmp(arg[iarg+1],"no")==0) check_obb_flag = 0;
293 else error->fix_error(FLERR,this,"");
294 if(check_ol_flag==0) check_obb_flag = 0;
295 iarg += 2;
296 hasargs = true;
297 }
298 #endif
299 else if(strcmp(style,"insert") == 0) error->fix_error(FLERR,this,"unknown keyword");
300 }
301
302 // memory not allocated initially
303 ninsert_this_max_local = 0;
304
305 // check for missing or contradictory settings
306 sanity_check();
307
308 //min/max type to be inserted, need that to check if material properties defined for all materials
309 type_max = fix_distribution->max_type();
310 type_min = fix_distribution->min_type();
311
312 // allgather arrays
313 MPI_Comm_rank(world,&me);
314 MPI_Comm_size(world,&nprocs);
315 recvcounts = new int[nprocs];
316 displs = new int[nprocs];
317
318 // set next reneighbor
319 force_reneighbor = 1;
320 next_reneighbor = first_ins_step;
321 most_recent_ins_step = -1;
322
323 vector_flag = 1;
324 size_vector = 2;
325 global_freq = 1;
326
327 print_stats_start_flag = 1;
328
329 irregular = new Irregular(lmp);
330
331 // calc max insertion radius
332 int ntypes = atom->ntypes;
333 maxrad = 0.;
334 minrad = 1000.;
335 for(int i = 1; i <= ntypes; i++)
336 {
337 maxrad = std::max(maxrad,max_rad(i));
338 minrad = std::min(minrad,min_rad(i));
339 }
340 }
341
342 /* ---------------------------------------------------------------------- */
343
~FixInsert()344 FixInsert::~FixInsert()
345 {
346 delete random;
347 delete [] recvcounts;
348 delete [] displs;
349 delete &neighList;
350 if(property_name) delete []property_name;
351
352 if(irregular) delete irregular;
353 irregular = 0;
354 }
355
356 /* ---------------------------------------------------------------------- */
357
setup(int vflag)358 void FixInsert::setup(int vflag)
359 {
360
361 // do this only once
362 if(setup_flag) return;
363 else setup_flag = true;
364
365 // calculate ninsert, insert_every, ninsert_per
366 calc_insertion_properties();
367
368 // calc last step of insertion
369 if(ninsert_exists)
370 {
371 if(ninsert <= ninsert_per)
372 final_ins_step = first_ins_step;
373 else
374 final_ins_step = first_ins_step +
375 static_cast<int>(static_cast<double>(ninsert)/ninsert_per) * static_cast<double>(insert_every);
376
377 if(final_ins_step < 0)
378 error->fix_error(FLERR,this,"Particle insertion: Overflow - need too long for particle insertion. "
379 "Please decrease # particles to insert or increase insertion rate");
380 if(ninsert < 0)
381 error->fix_error(FLERR,this,"Particle insertion: Overflow - too many particles for particle insertion. "
382 "Please decrease # particles to insert.");
383 }
384 else
385 final_ins_step = -1;
386
387 // print statistics
388 print_stats_start();
389
390 }
391
392 /* ---------------------------------------------------------------------- */
393
init_defaults()394 void FixInsert::init_defaults()
395 {
396 // default is that total # of particles to insert by this command is known
397 ninsert_exists = 1;
398
399 ninsert = ninserted = 0;
400 massinsert = massinserted = 0.;
401 nflowrate = massflowrate = 0.;
402
403 insert_every = -1;
404 ninsert_per = 0.;
405
406 // 1st insertion on next timestep is default
407 first_ins_step = update->ntimestep + 1;
408
409 maxattempt = 50;
410
411 check_ol_flag = 1;
412 all_in_flag = 0;
413
414 exact_number = 1;
415
416 v_randomSetting = RANDOM_CONSTANT;
417 vectorZeroize3D(v_insert);
418 vectorZeroize3D(v_insertFluct);
419 vectorZeroize3D(omega_insert);
420
421 quatIdentity4D(quat_insert);
422 quat_random_ = false;
423
424 print_stats_during_flag = 1;
425 warn_boxentent = true;
426
427 property_name = 0;
428 fix_property = 0;
429 fix_property_value = 0.;
430 }
431
432 /* ---------------------------------------------------------------------- */
433
sanity_check()434 void FixInsert::sanity_check()
435 {
436 if(fix_distribution == NULL)
437 error->fix_error(FLERR,this,"have to define a 'distributiontemplate'");
438
439 if(MathExtraLiggghts::abs(vectorMag4DSquared(quat_insert)-1.) > 1e-10)
440 error->fix_error(FLERR,this,"quaternion not valid");
441
442 if(ninsert > 0 && massinsert > 0.)
443 error->fix_error(FLERR,this,"must not define both 'nparticles' and 'mass'");
444 if(nflowrate > 0. && massflowrate > 0.)
445 error->fix_error(FLERR,this,"must not define both 'particlerate' and 'massrate'");
446
447 if(insert_every == 0 && (massflowrate > 0. || nflowrate > 0.))
448 error->fix_error(FLERR,this,"must not define 'particlerate' or 'massrate' for 'insert_every' = 0");
449
450 if(0 == comm->me)
451 {
452
453 std::vector<int> seeds;
454 seeds.push_back(random->state());
455 seeds.push_back(fix_distribution->random_state());
456 for(int itemplate = 0; itemplate < fix_distribution->n_particletemplates(); itemplate++)
457 {
458
459 seeds.push_back(fix_distribution->particletemplates()[itemplate]->random_insertion_state());
460 }
461
462 std::sort(seeds.begin(),seeds.end());
463
464 if(std::unique(seeds.begin(),seeds.end()) !=seeds.end() )
465 {
466 char errstr[1024];
467 sprintf(errstr,"Fix %s, ID %s: Random number generation: It is required that all the random seeds of this fix insert/*, \n"
468 " the random seed of particle distribution fix (id %s) template and all random seeds of the \n"
469 " fix particletemplate/* commands used by particle distribution fix (id %s) are different\n"
470 " Hint: possible valid (different) seeds would be the following numbers:\n"
471 " 15485863, 15485867, 32452843, 32452867, 49979687, 49979693, 67867967, 67867979, 86028121, 86028157",
472 style,id,fix_distribution->id,fix_distribution->id);
473
474 if(input->seed_check_throw_error())
475 error->one(FLERR,errstr);
476 else
477 error->warning(FLERR,errstr);
478 }
479 }
480 }
481
482 /* ---------------------------------------------------------------------- */
483
print_stats_start()484 void FixInsert::print_stats_start()
485 {
486 if (me == 0 && print_stats_start_flag) {
487
488 if(ninsert_exists)
489 {
490 if (screen)
491 fprintf(screen ,"INFO: Particle insertion %s: %f particles every %d steps - particle rate %f (mass rate %e)\n"
492 " %d particles (mass %e) within %d steps\n",
493 id,ninsert_per,insert_every,nflowrate,massflowrate,ninsert,massinsert,final_ins_step-first_ins_step);
494
495 if (logfile)
496 fprintf(logfile,"INFO: Particle insertion %s: %f particles every %d steps - particle rate %f, (mass rate %e)\n"
497 " %d particles (mass %e) within %d steps\n",
498 id,ninsert_per,insert_every,nflowrate,massflowrate,ninsert,massinsert,final_ins_step-first_ins_step);
499 }
500 else if(massflowrate > 0.)
501 {
502 if (screen)
503 fprintf(screen ,"INFO: Particle insertion %s: %f particles every %d steps - particle rate %f (mass rate %e)\n",
504 id,ninsert_per,insert_every,nflowrate,massflowrate);
505
506 if (logfile)
507 fprintf(logfile,"INFO: Particle insertion %s: %f particles every %d steps - particle rate %f, (mass rate %e)\n",
508 id,ninsert_per,insert_every,nflowrate,massflowrate);
509 }
510 else
511 {
512 if (screen)
513 fprintf(screen ,"INFO: Particle insertion %s: inserting every %d steps\n",id,insert_every);
514
515 if (logfile)
516 fprintf(logfile ,"INFO: Particle insertion %s: inserting every %d steps\n",id,insert_every);
517 }
518 }
519 }
520
521 /* ---------------------------------------------------------------------- */
522
print_stats_during(int ninsert_this,double mass_inserted_this)523 void FixInsert::print_stats_during(int ninsert_this, double mass_inserted_this)
524 {
525 bigint step = update->ntimestep;
526
527 if (me == 0 && print_stats_during_flag)
528 {
529 if (screen)
530 fprintf(screen ,"INFO: Particle insertion %s: inserted %d particle templates (mass %e) at step " BIGINT_FORMAT "\n"
531 " - a total of %d particle templates (mass %e) inserted so far.\n",
532 id,ninsert_this,mass_inserted_this,step,ninserted,massinserted);
533
534 if (logfile)
535 fprintf(logfile,"INFO: Particle insertion %s: inserted %d particle templates (mass %e) at step " BIGINT_FORMAT "\n"
536 " - a total of %d particle templates (mass %e) inserted so far.\n",
537 id,ninsert_this,mass_inserted_this,step,ninserted,massinserted);
538 }
539 }
540
541 /* ---------------------------------------------------------------------- */
542
setmask()543 int FixInsert::setmask()
544 {
545 int mask = 0;
546 mask |= PRE_EXCHANGE;
547 return mask;
548 }
549
550 /* ---------------------------------------------------------------------- */
551
init()552 void FixInsert::init()
553 {
554 int ntimestep = update->ntimestep;
555
556 if (!atom->radius_flag || !atom->rmass_flag)
557 error->fix_error(FLERR,this,"Fix insert requires atom attributes radius, rmass");
558 if (domain->triclinic)
559 error->fix_error(FLERR,this,"Cannot use with triclinic box");
560 if (domain->dimension != 3)
561 error->fix_error(FLERR,this,"Can use fix insert for 3d simulations only");
562
563 fix_multisphere = static_cast<FixMultisphere*>(modify->find_fix_style("multisphere", 0));
564 if(!fix_multisphere) multisphere = NULL;
565 else multisphere = &fix_multisphere->data();
566
567 // in case of new fix insert in a restarted simulation, have to add current time-step
568 if(next_reneighbor > 0 && next_reneighbor < ntimestep)
569 {
570
571 error->fix_error(FLERR,this,"'start' step can not be before current step");
572 }
573
574 if(property_name)
575 {
576 fix_property = static_cast<FixPropertyAtom*>(modify->find_fix_property(property_name,"property/atom","scalar",1,1,this->style,true));
577 }
578 }
579
580 /* ---------------------------------------------------------------------- */
581
reset_timestep(bigint newstep,bigint oldstep)582 void FixInsert::reset_timestep(bigint newstep,bigint oldstep)
583 {
584
585 next_reneighbor += (newstep-oldstep);
586
587 }
588
589 /* ---------------------------------------------------------------------- */
590
min_type()591 int FixInsert::min_type()
592 {
593 return type_min;
594 }
595
596 /* ---------------------------------------------------------------------- */
597
max_type()598 int FixInsert::max_type()
599 {
600 return type_max;
601 }
602
603 /* ---------------------------------------------------------------------- */
604
max_rad(int type)605 double FixInsert::max_rad(int type)
606 {
607 return fix_distribution->max_rad(type);
608 }
609
610 /* ---------------------------------------------------------------------- */
611
min_rad(int type)612 double FixInsert::min_rad(int type)
613 {
614 return fix_distribution->min_rad(type);
615 }
616
617 /* ---------------------------------------------------------------------- */
618
max_r_bound()619 double FixInsert::max_r_bound()
620 {
621 return fix_distribution->max_r_bound();
622 }
623
624 /* ---------------------------------------------------------------------- */
625
extend_cut_ghost()626 double FixInsert::extend_cut_ghost()
627 {
628
629 if(!fix_multisphere)
630 return 0.;
631
632 return 2.*fix_distribution->max_r_bound();
633 }
634
635 /* ---------------------------------------------------------------------- */
636
calc_ninsert_this()637 int FixInsert::calc_ninsert_this()
638 {
639 if(ninsert_per == 0.) error->fix_error(FLERR,this,"ninsert_per == 0.");
640
641 // number of bodies to insert this timestep
642 int ninsert_this = static_cast<int>(ninsert_per + random->uniform());
643 if (ninsert_exists && ninserted + ninsert_this > ninsert) ninsert_this = ninsert - ninserted;
644
645 return ninsert_this;
646 }
647
648 /* ----------------------------------------------------------------------
649 perform particle insertion
650 ------------------------------------------------------------------------- */
651
pre_exchange()652 void FixInsert::pre_exchange()
653 {
654
655 int ninsert_this, ninsert_this_local; // global and local # bodies to insert this time-step
656
657 // just return if should not be called on this timestep
658
659 if (next_reneighbor != update->ntimestep || most_recent_ins_step == update->ntimestep) return;
660 most_recent_ins_step = update->ntimestep;
661
662 // things to be done before inserting new particles
663
664 if(!pre_insert())
665 return;
666
667 // number of particles to insert this timestep
668 ninsert_this = calc_ninsert_this();
669
670 // limit to max number of particles that shall be inserted
671 // to avoid that max # may be slightly exceeded by random processes
672 // in fix_distribution->randomize_list, set exact_number to 1
673 if (ninsert_exists && ninserted + ninsert_this >= ninsert)
674 {
675 ninsert_this = ninsert - ninserted;
676 if(ninsert_this < 0)
677 ninsert_this = 0;
678 exact_number = 1;
679 }
680
681 // distribute ninsert_this across processors
682 ninsert_this_local = distribute_ninsert_this(ninsert_this);
683
684 // re-allocate list if necessary
685
686 if(ninsert_this_local > ninsert_this_max_local)
687 {
688 init_list(ninsert_this_local);
689 ninsert_this_max_local = ninsert_this_local;
690 }
691
692 // generate list of insertions
693 // number of inserted particles can change if exact_number = 0
694
695 ninsert_this_local = generate_list(ninsert_this_local,groupbit,exact_number);
696
697 MPI_Sum_Scalar(ninsert_this_local,ninsert_this,world);
698
699 if(ninsert_this == 0)
700 {
701 // warn if flowrate should be fulfilled
702 if((nflowrate > 0. || massflowrate > 0.) && comm->me == 0)
703 error->warning(FLERR,"Particle insertion: Inserting no particle - check particle insertion settings");
704
705 // schedule next insertion
706 if (insert_every && (!ninsert_exists || ninserted < ninsert))
707 next_reneighbor += insert_every;
708
709 else if(0 == insert_every)
710 next_reneighbor = -1;
711
712 return;
713 }
714 else if(ninsert_this < 0)
715 {
716
717 error->fix_error(FLERR,this,"Particle insertion: Internal error");
718 }
719
720 double min_subbox_extent;
721 int min_dim;
722 domain->min_subbox_extent(min_subbox_extent,min_dim);
723
724 if(warn_boxentent && min_subbox_extent < 2.2 *max_r_bound())
725 {
726 char msg[256];
727 sprintf(msg,"Particle insertion on proc %d: sub-domain is smaller than the bounding radius of insert particles to insert: \nMax. bounding "
728 "sphere diameter is %f, sub-domain extent in %s direction is only %f ",
729 comm->me,2.*max_r_bound(),0==min_dim?"x":(1==min_dim?"y":"z"),min_subbox_extent);
730 error->warning(FLERR,msg);
731 }
732
733 // warn if max # insertions exceeded by random processes
734 if (ninsert_exists && ninserted + ninsert_this > ninsert)
735 {
736 error->warning(FLERR,"INFO: Particle insertion: Number of particles to insert was slightly exceeded by random process");
737 }
738
739 // fill xnear array with particles to check overlap against
740
741 // add particles in insertion volume to xnear list
742 neighList.reset();
743
744 if(check_ol_flag)
745 load_xnear(ninsert_this_local);
746
747 // insertion counters in this step
748 int ninserted_this = 0, ninserted_spheres_this = 0;
749 int ninserted_this_local = 0, ninserted_spheres_this_local = 0;
750 double mass_inserted_this = 0.;
751 double mass_inserted_this_local = 0.;
752
753 // randomize insertion positions and set v, omega
754 // also performs overlap check via xnear if requested
755 // returns # bodies and # spheres that could actually be inserted
756 x_v_omega(ninsert_this_local,ninserted_this_local,ninserted_spheres_this_local,mass_inserted_this_local);
757
758 // actual particle insertion
759
760 fix_distribution->pre_insert(ninserted_this_local,fix_property,fix_property_value);
761
762 ninserted_spheres_this_local = fix_distribution->insert(ninserted_this_local);
763
764 // warn if max # insertions exceeded by random processes
765 if (ninsert_exists && ninserted + ninsert_this > ninsert)
766 {
767 error->warning(FLERR,"INFO: Particle insertion: Number of particles to insert was slightly exceeded by random process");
768 }
769
770 // set tag # of new particles beyond all previous atoms, reset global natoms
771 // if global map exists, reset it now instead of waiting for comm
772 // since deleting atoms messes up ghosts
773 int step = update->ntimestep;
774
775 if (atom->tag_enable)
776 {
777
778 //force all tags to be reset by setting them to zero
779 if(compress_flag)
780 {
781 if(comm->me == 0)
782 printf("FixInsertStream: resetting tags @ step %d. \n", step);
783 int *tag = atom->tag;
784 for (int i = 0; i < atom->nlocal; i++)
785 tag[i] = 0;
786 }
787
788 atom->tag_extend();
789 atom->natoms += static_cast<double>(ninserted_spheres_this);
790 if (atom->map_style)
791 {
792 atom->nghost = 0;
793 atom->map_init();
794 atom->map_set();
795 }
796 }
797
798 // give particle distributions the chance to do some wrap-up
799
800 fix_distribution->finalize_insertion();
801
802 // give derived classes the chance to do some wrap-up
803
804 finalize_insertion(ninserted_spheres_this_local);
805
806 // tally stats
807 MPI_Sum_Scalar(ninserted_this_local,ninserted_this,world);
808 ninserted += ninserted_this;
809 MPI_Sum_Scalar(mass_inserted_this_local,mass_inserted_this,world);
810 massinserted += mass_inserted_this;
811 print_stats_during(ninserted_this,mass_inserted_this);
812
813 if(ninserted_this < ninsert_this && comm->me == 0)
814 error->warning(FLERR,"Particle insertion: Less insertions than requested");
815
816 if (irregular->migrate_check())
817 irregular->migrate_atoms();
818
819 // next timestep to insert
820 if (insert_every && (!ninsert_exists || ninserted < ninsert)) next_reneighbor += insert_every;
821 else next_reneighbor = 0;
822
823 }
824
825 /* ----------------------------------------------------------------------
826 distribute insertions across processors
827 ------------------------------------------------------------------------- */
828
distribute_ninsert_this(int ninsert_this)829 int FixInsert::distribute_ninsert_this(int ninsert_this)
830 {
831 int me, nprocs, ngap, ninsert_this_local, *ninsert_this_local_all;
832 double fraction_local, fraction_local_all_sum, *fraction_local_all, *remainder, r, rsum;
833
834 me = comm->me;
835 nprocs = comm->nprocs;
836
837 fraction_local = insertion_fraction();
838
839 if(!exact_number)
840 return static_cast<int>(fraction_local*static_cast<double>(ninsert_this) + random->uniform());
841
842 // for exact_number==1, have to allgather to exactly match ninsert_this
843
844 fraction_local_all = new double[nprocs];
845 remainder = new double[nprocs];
846 ninsert_this_local_all = new int[nprocs];
847
848 // allgather local fractions
849 MPI_Allgather(&fraction_local,1,MPI_DOUBLE,fraction_local_all,1,MPI_DOUBLE,world);
850
851 // proc0 calculates ninsert_this_local for all processes
852
853 if(me == 0)
854 {
855 // remove fractions < 2% / nprocs
856 // have to normalize so not all portions get cancelled away for higher proc counts
857 // normalize fraction_local_all so sum across processors is 1
858
859 double lower_thresh = 0.02 / static_cast<double>(nprocs);
860
861 fraction_local_all_sum = 0.;
862 for(int iproc = 0; iproc < nprocs; iproc++)
863 {
864 if(fraction_local_all[iproc] < lower_thresh)
865 fraction_local_all[iproc] = 0.;
866 fraction_local_all_sum += fraction_local_all[iproc];
867 }
868
869 if(fraction_local_all_sum == 0.)
870 error->one(FLERR,"Internal error distributing particles to processes");
871
872 for(int iproc = 0; iproc < nprocs; iproc++)
873 fraction_local_all[iproc] /= fraction_local_all_sum;
874
875 rsum = 0.;
876 for(int iproc = 0; iproc < nprocs; iproc++)
877 {
878 ninsert_this_local_all[iproc] = static_cast<int>(fraction_local_all[iproc]*static_cast<double>(ninsert_this));
879 remainder[iproc] = fraction_local_all[iproc]*static_cast<double>(ninsert_this) - ninsert_this_local_all[iproc];
880 rsum += remainder[iproc];
881
882 }
883
884 ngap = round(rsum);
885
886 for(int i = 0; i < ngap; i++)
887 {
888 r = random->uniform() * static_cast<double>(ngap);
889 int iproc = 0;
890 rsum = remainder[iproc];
891
892 while(iproc < (nprocs-1) && rsum < r)
893 {
894 iproc++;
895 rsum += remainder[iproc];
896 }
897 ninsert_this_local_all[iproc]++;
898 }
899 }
900
901 // Bcast the result
902 MPI_Bcast(ninsert_this_local_all,nprocs, MPI_INT,0,world);
903 ninsert_this_local = ninsert_this_local_all[me];
904
905 delete []fraction_local_all;
906 delete []remainder;
907 delete []ninsert_this_local_all;
908
909 return ninsert_this_local;
910 }
911
912 /* ----------------------------------------------------------------------
913 fill neighbor list with nearby particles
914 ------------------------------------------------------------------------- */
915
load_xnear(int ninsert_this_local)916 int FixInsert::load_xnear(int ninsert_this_local)
917 {
918 // load up neighbor list with local and ghosts
919
920 double **x = atom->x;
921 double *radius = atom->radius;
922 const int nall = atom->nlocal + atom->nghost;
923
924 BoundingBox bb = getBoundingBox();
925 neighList.reset();
926 #ifdef SUPERQUADRIC_ACTIVE_FLAG
927 neighList.set_obb_flag(check_obb_flag);
928 #endif
929
930 if(neighList.setBoundingBox(bb, maxrad,true,true))
931 {
932 for (int i = 0; i < nall; ++i)
933 {
934
935 if (is_nearby(i) && neighList.isInBoundingBox(x[i]) )
936 {
937 #ifdef SUPERQUADRIC_ACTIVE_FLAG
938 if(atom->superquadric_flag and check_obb_flag)
939 neighList.insert_superquadric(x[i], radius[i], atom->quaternion[i], atom->shape[i], atom->blockiness[i]);
940 else
941 neighList.insert(x[i], radius[i]);
942 #else
943 neighList.insert(x[i], radius[i]);
944 #endif
945 }
946 }
947 }
948
949 return neighList.count();
950 }
951
952 /* ----------------------------------------------------------------------
953 generate random velocity based on random setting
954 ------------------------------------------------------------------------- */
955
generate_random_velocity(double * velocity)956 void FixInsert::generate_random_velocity(double * velocity) {
957 switch(v_randomSetting) {
958 case RANDOM_UNIFORM:
959 velocity[0] = v_insert[0] + v_insertFluct[0] * 2.0 * (random->uniform()-0.50);
960 velocity[1] = v_insert[1] + v_insertFluct[1] * 2.0 * (random->uniform()-0.50);
961 velocity[2] = v_insert[2] + v_insertFluct[2] * 2.0 * (random->uniform()-0.50);
962 break;
963
964 case RANDOM_GAUSSIAN:
965 velocity[0] = v_insert[0] + v_insertFluct[0] * random->gaussian();
966 velocity[1] = v_insert[1] + v_insertFluct[1] * random->gaussian();
967 velocity[2] = v_insert[2] + v_insertFluct[2] * random->gaussian();
968 }
969 }
970
971 /* ----------------------------------------------------------------------
972 pack entire state of Fix into one write
973 ------------------------------------------------------------------------- */
974
write_restart(FILE * fp)975 void FixInsert::write_restart(FILE *fp)
976 {
977 int n = 0;
978 double list[5];
979 list[n++] = static_cast<double>(random->state());
980 list[n++] = static_cast<double>(ninserted);
981 list[n++] = static_cast<double>(first_ins_step);
982 list[n++] = static_cast<double>(next_reneighbor);
983 list[n++] = massinserted;
984
985 if (comm->me == 0) {
986 int size = n * sizeof(double);
987 fwrite(&size,sizeof(int),1,fp);
988 fwrite(list,sizeof(double),n,fp);
989 }
990 }
991
992 /* ----------------------------------------------------------------------
993 use state info from restart file to restart the Fix
994 ------------------------------------------------------------------------- */
995
restart(char * buf)996 void FixInsert::restart(char *buf)
997 {
998 int n = 0;
999 double *list = (double *) buf;
1000 bigint next_reneighbor_re;
1001
1002 seed = static_cast<int> (list[n++]) + comm->me;
1003 ninserted = static_cast<int> (list[n++]);
1004 first_ins_step = static_cast<int> (list[n++]);
1005 next_reneighbor_re = static_cast<bigint> (list[n++]);
1006 massinserted = list[n++];
1007
1008 random->reset(seed);
1009
1010 // in order to be able to continue pouring with increased number of particles
1011 // if insert was already finished in run to be restarted
1012 if(next_reneighbor_re != 0) next_reneighbor = next_reneighbor_re;
1013
1014 }
1015
1016 /* ----------------------------------------------------------------------
1017 output
1018 ------------------------------------------------------------------------- */
1019
compute_vector(int index)1020 double FixInsert::compute_vector(int index)
1021 {
1022 if(index == 0) return static_cast<double>(ninserted);
1023 if(index == 1) return massinserted;
1024 return 0.0;
1025 }
1026