1 /* ----------------------------------------------------------------------
2    SPARTA - Stochastic PArallel Rarefied-gas Time-accurate Analyzer
3    http://sparta.sandia.gov
4    Steve Plimpton, sjplimp@sandia.gov, Michael Gallis, magalli@sandia.gov
5    Sandia National Laboratories
6 
7    Copyright (2014) Sandia Corporation.  Under the terms of Contract
8    DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
9    certain rights in this software.  This software is distributed under
10    the GNU General Public License.
11 
12    See the README file in the top-level SPARTA directory.
13 ------------------------------------------------------------------------- */
14 
15 #include "math.h"
16 #include "stdlib.h"
17 #include "string.h"
18 #include "create_particles.h"
19 #include "update.h"
20 #include "particle.h"
21 #include "mixture.h"
22 #include "grid.h"
23 #include "modify.h"
24 #include "comm.h"
25 #include "domain.h"
26 #include "region.h"
27 #include "input.h"
28 #include "variable.h"
29 #include "random_mars.h"
30 #include "random_knuth.h"
31 #include "math_const.h"
32 #include "memory.h"
33 #include "error.h"
34 
35 using namespace SPARTA_NS;
36 using namespace MathConst;
37 
38 enum{UNKNOWN,OUTSIDE,INSIDE,OVERLAP};   // same as Grid
39 
40 #define EPSZERO 1.0e-14
41 
42 /* ---------------------------------------------------------------------- */
43 
CreateParticles(SPARTA * sparta)44 CreateParticles::CreateParticles(SPARTA *sparta) : Pointers(sparta) {}
45 
46 /* ---------------------------------------------------------------------- */
47 
command(int narg,char ** arg)48 void CreateParticles::command(int narg, char **arg)
49 {
50   if (!grid->exist)
51     error->all(FLERR,"Cannot create particles before grid is defined");
52 
53   particle->exist = 1;
54 
55   if (narg < 1) error->all(FLERR,"Illegal create_particles command");
56 
57   imix = particle->find_mixture(arg[0]);
58   if (imix < 0) error->all(FLERR,"Create_particles mixture ID does not exist");
59   particle->mixture[imix]->init();
60 
61   // style arg
62 
63   bigint np = 0;
64   single = 0;
65 
66   int iarg = 1;
67   if (strcmp(arg[iarg],"n") == 0) {
68     if (iarg+2 > narg) error->all(FLERR,"Illegal create_particles command");
69     np = ATOBIGINT(arg[iarg+1]);
70       if (np < 0) error->all(FLERR,"Illegal create_particles command");
71       iarg += 2;
72   } else if (strcmp(arg[iarg],"single") == 0) {
73     if (iarg+8 > narg) error->all(FLERR,"Illegal create_particles command");
74     single = 1;
75     mspecies = particle->find_species(arg[iarg+1]);
76     if (mspecies < 0)
77       error->all(FLERR,"Create_particles species ID does not exist");
78     xp = atof(arg[iarg+2]);
79     yp = atof(arg[iarg+3]);
80     zp = atof(arg[iarg+4]);
81     vx = atof(arg[iarg+5]);
82     vy = atof(arg[iarg+6]);
83     vz = atof(arg[iarg+7]);
84     iarg += 8;
85   } else error->all(FLERR,"Illegal create_particles command");
86 
87   // optional args
88 
89   int globalflag = 0;
90   twopass = 0;
91   region = NULL;
92   speciesflag = densflag = velflag = tempflag = 0;
93   sstr = sxstr = systr = szstr = NULL;
94   dstr = dxstr = dystr = dzstr = NULL;
95   tstr = txstr = tystr = tzstr = NULL;
96   vxstr = vystr = vzstr = vstrx = vstry = vstrz = NULL;
97 
98   while (iarg < narg) {
99     if (strcmp(arg[iarg],"global") == 0) {
100       if (iarg+2 > narg) error->all(FLERR,"Illegal create_particles command");
101       if (strcmp(arg[iarg+1],"no") == 0) globalflag = 0;
102       else if (strcmp(arg[iarg+1],"yes") == 0) globalflag = 1;
103       else error->all(FLERR,"Illegal create_particles command");
104       iarg += 2;
105     } else if (strcmp(arg[iarg],"region") == 0) {
106       if (iarg+2 > narg) error->all(FLERR,"Illegal create_particles command");
107       int iregion = domain->find_region(arg[iarg+1]);
108       if (iregion < 0)
109         error->all(FLERR,"Create_particles region does not exist");
110       region = domain->regions[iregion];
111       iarg += 2;
112     } else if (strcmp(arg[iarg],"species") == 0) {
113       if (iarg+5 > narg) error->all(FLERR,"Illegal create_particles command");
114       speciesflag = 1;
115       sstr = arg[iarg+1];
116       if (strcmp(arg[iarg+2],"NULL") == 0) sxstr = NULL;
117       else sxstr = arg[iarg+2];
118       if (strcmp(arg[iarg+3],"NULL") == 0) systr = NULL;
119       else systr = arg[iarg+3];
120       if (strcmp(arg[iarg+4],"NULL") == 0) szstr = NULL;
121       else szstr = arg[iarg+4];
122       iarg += 5;
123     } else if (strcmp(arg[iarg],"density") == 0) {
124       if (iarg+5 > narg) error->all(FLERR,"Illegal create_particles command");
125       densflag = 1;
126       dstr = arg[iarg+1];
127       if (strcmp(arg[iarg+2],"NULL") == 0) dxstr = NULL;
128       else dxstr = arg[iarg+2];
129       if (strcmp(arg[iarg+3],"NULL") == 0) dystr = NULL;
130       else dystr = arg[iarg+3];
131       if (strcmp(arg[iarg+4],"NULL") == 0) dzstr = NULL;
132       else dzstr = arg[iarg+4];
133       iarg += 5;
134     } else if (strcmp(arg[iarg],"temperature") == 0) {
135       if (iarg+5 > narg) error->all(FLERR,"Illegal create_particles command");
136       tempflag = 1;
137       tstr = arg[iarg+1];
138       if (strcmp(arg[iarg+2],"NULL") == 0) txstr = NULL;
139       else txstr = arg[iarg+2];
140       if (strcmp(arg[iarg+3],"NULL") == 0) tystr = NULL;
141       else tystr = arg[iarg+3];
142       if (strcmp(arg[iarg+4],"NULL") == 0) tzstr = NULL;
143       else tzstr = arg[iarg+4];
144       iarg += 5;
145     } else if (strcmp(arg[iarg],"velocity") == 0) {
146       if (iarg+7 > narg) error->all(FLERR,"Illegal create_particles command");
147       velflag = 1;
148       if (strcmp(arg[iarg+1],"NULL") == 0) vxstr = NULL;
149       else vxstr = arg[iarg+1];
150       if (strcmp(arg[iarg+2],"NULL") == 0) vystr = NULL;
151       else vystr = arg[iarg+2];
152       if (strcmp(arg[iarg+3],"NULL") == 0) vzstr = NULL;
153       else vzstr = arg[iarg+3];
154       if (strcmp(arg[iarg+4],"NULL") == 0) vstrx = NULL;
155       else vstrx = arg[iarg+4];
156       if (strcmp(arg[iarg+5],"NULL") == 0) vstry = NULL;
157       else vstry = arg[iarg+5];
158       if (strcmp(arg[iarg+6],"NULL") == 0) vstrz = NULL;
159       else vstrz = arg[iarg+6];
160       iarg += 7;
161     } else if (strcmp(arg[iarg],"twopass") == 0) {
162       if (iarg+1 > narg) error->all(FLERR,"Illegal create_particles command");
163       twopass = 1;
164       iarg += 1;
165     } else error->all(FLERR,"Illegal create_particles command");
166   }
167 
168   if (globalflag)
169     error->all(FLERR,"Create_particles global option not yet implemented");
170 
171   // error checks and further setup for variables
172 
173   if (speciesflag) {
174     svar = input->variable->find(sstr);
175     if (svar < 0)
176       error->all(FLERR,"Variable name for create_particles does not exist");
177     if (!input->variable->equal_style(svar))
178       error->all(FLERR,"Variable for create_particles is invalid style");
179     if (sxstr) {
180       sxvar = input->variable->find(sxstr);
181       if (sxvar < 0)
182         error->all(FLERR,"Variable name for create_particles does not exist");
183       if (!input->variable->internal_style(sxvar))
184         error->all(FLERR,"Variable for create_particles is invalid style");
185     }
186     if (systr) {
187       syvar = input->variable->find(systr);
188       if (syvar < 0)
189         error->all(FLERR,"Variable name for create_particles does not exist");
190       if (!input->variable->internal_style(syvar))
191         error->all(FLERR,"Variable for create_particles is invalid style");
192     }
193     if (szstr) {
194       szvar = input->variable->find(szstr);
195       if (szvar < 0)
196         error->all(FLERR,"Variable name for create_particles does not exist");
197       if (!input->variable->internal_style(szvar))
198         error->all(FLERR,"Variable for create_particles is invalid style");
199     }
200   }
201 
202   if (densflag) {
203     dvar = input->variable->find(dstr);
204     if (dvar < 0)
205       error->all(FLERR,"Variable name for create_particles does not exist");
206     if (!input->variable->equal_style(dvar))
207       error->all(FLERR,"Variable for create_particles is invalid style");
208     if (dxstr) {
209       dxvar = input->variable->find(dxstr);
210       if (dxvar < 0)
211         error->all(FLERR,"Variable name for create_particles does not exist");
212       if (!input->variable->internal_style(dxvar))
213         error->all(FLERR,"Variable for create_particles is invalid style");
214     }
215     if (dystr) {
216       dyvar = input->variable->find(dystr);
217       if (dyvar < 0)
218         error->all(FLERR,"Variable name for create_particles does not exist");
219       if (!input->variable->internal_style(dyvar))
220         error->all(FLERR,"Variable for create_particles is invalid style");
221     }
222     if (dzstr) {
223       dzvar = input->variable->find(dzstr);
224       if (dzvar < 0)
225         error->all(FLERR,"Variable name for create_particles does not exist");
226       if (!input->variable->internal_style(dzvar))
227         error->all(FLERR,"Variable for create_particles is invalid style");
228     }
229   }
230 
231   if (tempflag) {
232     tvar = input->variable->find(tstr);
233     if (tvar < 0)
234       error->all(FLERR,"Variable name for create_particles does not exist");
235     if (!input->variable->equal_style(tvar))
236       error->all(FLERR,"Variable for create_particles is invalid style");
237     if (txstr) {
238       txvar = input->variable->find(txstr);
239       if (txvar < 0)
240         error->all(FLERR,"Variable name for create_particles does not exist");
241       if (!input->variable->internal_style(txvar))
242         error->all(FLERR,"Variable for create_particles is invalid style");
243     }
244     if (tystr) {
245       tyvar = input->variable->find(tystr);
246       if (tyvar < 0)
247         error->all(FLERR,"Variable name for create_particles does not exist");
248       if (!input->variable->internal_style(tyvar))
249         error->all(FLERR,"Variable for create_particles is invalid style");
250     }
251     if (tzstr) {
252       tzvar = input->variable->find(tzstr);
253       if (tzvar < 0)
254         error->all(FLERR,"Variable name for create_particles does not exist");
255       if (!input->variable->internal_style(tzvar))
256         error->all(FLERR,"Variable for create_particles is invalid style");
257     }
258   }
259 
260   if (velflag) {
261     if (vxstr) {
262       vxvar = input->variable->find(vxstr);
263       if (vxvar < 0)
264         error->all(FLERR,"Variable name for create_particles does not exist");
265       if (!input->variable->equal_style(vxvar))
266         error->all(FLERR,"Variable for create_particles is invalid style");
267     }
268     if (vystr) {
269       vyvar = input->variable->find(vystr);
270       if (vyvar < 0)
271         error->all(FLERR,"Variable name for create_particles does not exist");
272       if (!input->variable->equal_style(vyvar))
273         error->all(FLERR,"Variable for create_particles is invalid style");
274     }
275     if (vzstr) {
276       vzvar = input->variable->find(vzstr);
277       if (vzvar < 0)
278         error->all(FLERR,"Variable name for create_particles does not exist");
279       if (!input->variable->equal_style(vzvar))
280         error->all(FLERR,"Variable for create_particles is invalid style");
281     }
282     if (vstrx) {
283       vvarx = input->variable->find(vstrx);
284       if (vvarx < 0)
285         error->all(FLERR,"Variable name for create_particles does not exist");
286       if (!input->variable->internal_style(vvarx))
287         error->all(FLERR,"Variable for create_particles is invalid style");
288     }
289     if (vstry) {
290       vvary = input->variable->find(vstry);
291       if (vvary < 0)
292         error->all(FLERR,"Variable name for create_particles does not exist");
293       if (!input->variable->internal_style(vvary))
294         error->all(FLERR,"Variable for create_particles is invalid style");
295     }
296     if (vstrz) {
297       vvarz = input->variable->find(vstrz);
298       if (vvarz < 0)
299         error->all(FLERR,"Variable name for create_particles does not exist");
300       if (!input->variable->internal_style(vvarz))
301         error->all(FLERR,"Variable for create_particles is invalid style");
302     }
303   }
304 
305   // calculate Np if not set explicitly
306 
307   if (single) np = 1;
308   else if (np == 0) {
309     Grid::ChildCell *cells = grid->cells;
310     Grid::ChildInfo *cinfo = grid->cinfo;
311     int nglocal = grid->nlocal;
312 
313     double flowvolme = 0.0;
314     for (int icell = 0; icell < nglocal; icell++) {
315       if (cells[icell].nsplit > 1) continue;
316       if (cinfo[icell].type != INSIDE)
317         flowvolme += cinfo[icell].volume / cinfo[icell].weight;
318     }
319     double flowvol;
320     MPI_Allreduce(&flowvolme,&flowvol,1,MPI_DOUBLE,MPI_SUM,world);
321     np = particle->mixture[imix]->nrho * flowvol / update->fnum;
322   }
323 
324   // generate particles
325   // NOTE: invoke local or global option here
326 
327   if (comm->me == 0)
328     if (screen) fprintf(screen,"Creating particles ...\n");
329 
330   MPI_Barrier(world);
331   double time1 = MPI_Wtime();
332 
333   bigint nprevious = particle->nglobal;
334   if (single) create_single();
335   else if (!globalflag) {
336     if (twopass) create_local_twopass(np);
337     else create_local(np);
338   }
339   //else create_global(np);
340 
341   MPI_Barrier(world);
342   double time2 = MPI_Wtime();
343 
344   // error check
345   // only if no region and no variable species/density specified
346 
347   bigint nglobal;
348   bigint nme = particle->nlocal;
349   MPI_Allreduce(&nme,&nglobal,1,MPI_SPARTA_BIGINT,MPI_SUM,world);
350   if (!region && !speciesflag && !densflag && !tempflag &&
351       nglobal-nprevious != np) {
352     char str[128];
353     sprintf(str,"Created incorrect # of particles: "
354 	    BIGINT_FORMAT " versus " BIGINT_FORMAT,
355 	    nglobal-nprevious,np);
356     error->all(FLERR,str);
357   }
358   bigint ncreated = nglobal-nprevious;
359   particle->nglobal = nglobal;
360 
361   // print stats
362 
363   if (comm->me == 0) {
364     if (screen) {
365       fprintf(screen,"Created " BIGINT_FORMAT " particles\n",ncreated);
366       fprintf(screen,"  CPU time = %g secs\n",time2-time1);
367     }
368     if (logfile) {
369       fprintf(logfile,"Created " BIGINT_FORMAT " particles\n",ncreated);
370       fprintf(logfile,"  CPU time = %g secs\n",time2-time1);
371     }
372   }
373 }
374 
375 /* ----------------------------------------------------------------------
376    create a single particle
377    find cell it is in, and store on appropriate processor
378 ------------------------------------------------------------------------- */
379 
create_single()380 void CreateParticles::create_single()
381 {
382   double x[3],v[3],vstream[3];
383   double *lo,*hi;
384 
385   x[0] = xp;  x[1] = yp;  x[2] = zp;
386   v[0] = vx;  v[1] = vy;  v[2] = vz;
387   double temp_thermal = particle->mixture[imix]->temp_thermal;
388   double temp_rot = particle->mixture[imix]->temp_rot;
389   double temp_vib = particle->mixture[imix]->temp_vib;
390   vstream[0] = vstream[1] = vstream[2] = 0.0;
391 
392   if (domain->dimension == 2 && x[2] != 0.0)
393     error->all(FLERR,"Create_particles single requires z = 0 "
394 	       "for 2d simulation");
395 
396   Grid::ChildCell *cells = grid->cells;
397   Grid::ChildInfo *cinfo = grid->cinfo;
398   int nglocal = grid->nlocal;
399 
400   int iwhich = -1;
401   for (int icell = 0; icell < nglocal; icell++) {
402     if (cinfo[icell].type != OUTSIDE) continue;
403     lo = cells[icell].lo;
404     hi = cells[icell].hi;
405     if (x[0] >= lo[0] && x[0] < hi[0] &&
406 	x[1] >= lo[1] && x[1] < hi[1] &&
407 	x[2] >= lo[2] && x[2] < hi[2]) iwhich = icell;
408   }
409 
410   // insure that exactly one proc found cell to insert particle into
411 
412   int flag,flagall;
413   if (iwhich < 0) flag = 0;
414   else flag = 1;
415   MPI_Allreduce(&flag,&flagall,1,MPI_INT,MPI_SUM,world);
416   if (flagall != 1)
417     error->all(FLERR,"Could not create a single particle");
418 
419   // nfix_update_custom = # of fixes with update_custom() method
420 
421   particle->error_custom();
422   modify->list_init_fixes();
423   int nfix_update_custom = modify->n_update_custom;
424 
425   // add the particle
426 
427   RanKnuth *random = new RanKnuth(update->ranmaster->uniform());
428 
429   if (iwhich >= 0) {
430     int id = MAXSMALLINT*random->uniform();
431     double erot = particle->erot(mspecies,temp_rot,random);
432     double evib = particle->evib(mspecies,temp_vib,random);
433     particle->add_particle(id,mspecies,iwhich,x,v,erot,evib);
434     if (nfix_update_custom)
435       modify->update_custom(particle->nlocal-1,temp_thermal,
436                            temp_rot,temp_vib,vstream);
437   }
438 
439   delete random;
440 }
441 
442 /* ----------------------------------------------------------------------
443    create Np particles in parallel
444    every proc creates fraction of Np for cells it owns
445    only insert in cells uncut by surfs
446    account for cell weighting
447    attributes of created particle depend on number of procs
448 ------------------------------------------------------------------------- */
449 
create_local(bigint np)450 void CreateParticles::create_local(bigint np)
451 {
452   int dimension = domain->dimension;
453 
454   int me = comm->me;
455   RanKnuth *random = new RanKnuth(update->ranmaster->uniform());
456   double seed = update->ranmaster->uniform();
457   random->reset(seed,me,100);
458 
459   Grid::ChildCell *cells = grid->cells;
460   Grid::ChildInfo *cinfo = grid->cinfo;
461   int nglocal = grid->nlocal;
462 
463   // volme = volume of grid cells I own that are OUTSIDE surfs
464   // skip cells entirely outside region
465   // Nme = # of particles I will create
466   // MPI_Scan() logic insures sum of nme = Np
467 
468   double *lo,*hi;
469   double volone;
470 
471   double volme = 0.0;
472   for (int i = 0; i < nglocal; i++) {
473     if (cinfo[i].type != OUTSIDE) continue;
474     lo = cells[i].lo;
475     hi = cells[i].hi;
476     if (region && region->bboxflag && outside_region(dimension,lo,hi))
477       continue;
478 
479     if (dimension == 3) volone = (hi[0]-lo[0]) * (hi[1]-lo[1]) * (hi[2]-lo[2]);
480     else if (domain->axisymmetric)
481       volone = (hi[0]-lo[0]) * (hi[1]*hi[1]-lo[1]*lo[1])*MY_PI;
482     else volone = (hi[0]-lo[0]) * (hi[1]-lo[1]);
483     volme += volone / cinfo[i].weight;
484   }
485 
486   double volupto;
487   MPI_Scan(&volme,&volupto,1,MPI_DOUBLE,MPI_SUM,world);
488 
489   double *vols;
490   int nprocs = comm->nprocs;
491   memory->create(vols,nprocs,"create_particles:vols");
492   MPI_Allgather(&volupto,1,MPI_DOUBLE,vols,1,MPI_DOUBLE,world);
493 
494   // nme = # of particles for me to create
495   // gathered Scan results not guaranteed to be monotonically increasing
496   // can cause epsilon mis-counts for huge particle counts
497   // enforce that by brute force
498 
499   for (int i = 1; i < nprocs; i++)
500     if (vols[i] != vols[i-1] &&
501         fabs(vols[i]-vols[i-1])/vols[nprocs-1] < EPSZERO)
502       vols[i] = vols[i-1];
503 
504   bigint nstart,nstop;
505   if (me > 0) nstart = static_cast<bigint> (np * (vols[me-1]/vols[nprocs-1]));
506   else nstart = 0;
507   nstop = static_cast<bigint> (np * (vols[me]/vols[nprocs-1]));
508   bigint nme = nstop-nstart;
509 
510   memory->destroy(vols);
511 
512   // nfix_update_custom = # of fixes with update_custom() method
513 
514   particle->error_custom();
515   modify->list_init_fixes();
516   int nfix_update_custom = modify->n_update_custom;
517 
518   // loop over cells I own
519   // only add particles to OUTSIDE cells
520   // skip cells entirely outside region
521   // ntarget = floating point # of particles to create in one cell
522   // npercell = integer # of particles to create in one cell
523   // basing ntarget on accumulated volume and nprev insures Nme total creations
524   // particle species = random value based on mixture fractions
525   // particle velocity = stream velocity + thermal velocity
526 
527   int *species = particle->mixture[imix]->species;
528   double *cummulative = particle->mixture[imix]->cummulative;
529   double *vstream = particle->mixture[imix]->vstream;
530   double *vscale = particle->mixture[imix]->vscale;
531   int nspecies = particle->mixture[imix]->nspecies;
532   double temp_thermal = particle->mixture[imix]->temp_thermal;
533   double temp_rot = particle->mixture[imix]->temp_rot;
534   double temp_vib = particle->mixture[imix]->temp_vib;
535 
536   int npercell,ncreate,isp,ispecies,id;
537   double x[3],v[3],vstream_variable[3];
538   double ntarget,scale,rn,vn,vr,theta1,theta2,erot,evib;
539 
540   double tempscale = 1.0;
541   double sqrttempscale = 1.0;
542 
543   double volsum = 0.0;
544   bigint nprev = 0;
545 
546   for (int i = 0; i < nglocal; i++) {
547     if (cinfo[i].type != OUTSIDE) continue;
548     lo = cells[i].lo;
549     hi = cells[i].hi;
550     if (region && region->bboxflag && outside_region(dimension,lo,hi))
551       continue;
552 
553     if (dimension == 3) volone = (hi[0]-lo[0]) * (hi[1]-lo[1]) * (hi[2]-lo[2]);
554     else if (domain->axisymmetric)
555       volone = (hi[0]-lo[0]) * (hi[1]*hi[1]-lo[1]*lo[1])*MY_PI;
556     else volone = (hi[0]-lo[0]) * (hi[1]-lo[1]);
557     volsum += volone / cinfo[i].weight;
558 
559     ntarget = nme * volsum/volme - nprev;
560     npercell = static_cast<int> (ntarget);
561 
562     if (random->uniform() < ntarget-npercell) npercell++;
563     ncreate = npercell;
564 
565     if (densflag) {
566       scale = density_variable(lo,hi);
567       ntarget *= scale;
568       ncreate = static_cast<int> (ntarget);
569       if (random->uniform() < ntarget-ncreate) ncreate++;
570     }
571 
572     for (int m = 0; m < ncreate; m++) {
573       rn = random->uniform();
574 
575       isp = 0;
576       while (cummulative[isp] < rn) isp++;
577       ispecies = species[isp];
578 
579       x[0] = lo[0] + random->uniform() * (hi[0]-lo[0]);
580       x[1] = lo[1] + random->uniform() * (hi[1]-lo[1]);
581       x[2] = lo[2] + random->uniform() * (hi[2]-lo[2]);
582       if (dimension == 2) x[2] = 0.0;
583 
584       if (region && !region->match(x)) continue;
585       if (speciesflag) {
586         isp = species_variable(x) - 1;
587         if (isp < 0 || isp >= nspecies) continue;
588         ispecies = species[isp];
589       }
590 
591       if (tempflag) {
592         tempscale = temperature_variable(x);
593         sqrttempscale = sqrt(tempscale);
594       }
595 
596       vn = vscale[isp] * sqrttempscale * sqrt(-log(random->uniform()));
597       vr = vscale[isp] * sqrttempscale * sqrt(-log(random->uniform()));
598       theta1 = MY_2PI * random->uniform();
599       theta2 = MY_2PI * random->uniform();
600 
601       if (velflag) {
602         velocity_variable(x,vstream,vstream_variable);
603         v[0] = vstream_variable[0] + vn*cos(theta1);
604         v[1] = vstream_variable[1] + vr*cos(theta2);
605         v[2] = vstream_variable[2] + vr*sin(theta2);
606       } else {
607         v[0] = vstream[0] + vn*cos(theta1);
608         v[1] = vstream[1] + vr*cos(theta2);
609         v[2] = vstream[2] + vr*sin(theta2);
610       }
611 
612       erot = particle->erot(ispecies,temp_rot*tempscale,random);
613       evib = particle->evib(ispecies,temp_vib*tempscale,random);
614 
615       id = MAXSMALLINT*random->uniform();
616 
617       particle->add_particle(id,ispecies,i,x,v,erot,evib);
618 
619       if (nfix_update_custom)
620         modify->update_custom(particle->nlocal-1,temp_thermal,
621                              temp_rot,temp_vib,vstream);
622     }
623 
624     // increment count without effect of density variation
625     // so that target insertion count is undisturbed
626 
627     nprev += npercell;
628   }
629 
630   delete random;
631 }
632 
633 /* ----------------------------------------------------------------------
634    create Np particles in parallel in two passes
635    every proc creates fraction of Np for cells it owns
636    only insert in cells uncut by surfs
637    account for cell weighting
638    attributes of created particle depend on number of procs
639 ------------------------------------------------------------------------- */
640 
create_local_twopass(bigint np)641 void CreateParticles::create_local_twopass(bigint np)
642 {
643   int dimension = domain->dimension;
644 
645   int me = comm->me;
646   RanKnuth *random = new RanKnuth(update->ranmaster->uniform());
647   double seed = update->ranmaster->uniform();
648   random->reset(seed,me,100);
649 
650   Grid::ChildCell *cells = grid->cells;
651   Grid::ChildInfo *cinfo = grid->cinfo;
652   int nglocal = grid->nlocal;
653 
654   // volme = volume of grid cells I own that are OUTSIDE surfs
655   // skip cells entirely outside region
656   // Nme = # of particles I will create
657   // MPI_Scan() logic insures sum of nme = Np
658 
659   double *lo,*hi;
660   double volone;
661 
662   double volme = 0.0;
663   for (int i = 0; i < nglocal; i++) {
664     if (cinfo[i].type != OUTSIDE) continue;
665     lo = cells[i].lo;
666     hi = cells[i].hi;
667     if (region && region->bboxflag && outside_region(dimension,lo,hi))
668       continue;
669 
670     if (dimension == 3) volone = (hi[0]-lo[0]) * (hi[1]-lo[1]) * (hi[2]-lo[2]);
671     else if (domain->axisymmetric)
672       volone = (hi[0]-lo[0]) * (hi[1]*hi[1]-lo[1]*lo[1])*MY_PI;
673     else volone = (hi[0]-lo[0]) * (hi[1]-lo[1]);
674     volme += volone / cinfo[i].weight;
675   }
676 
677   double volupto;
678   MPI_Scan(&volme,&volupto,1,MPI_DOUBLE,MPI_SUM,world);
679 
680   double *vols;
681   int nprocs = comm->nprocs;
682   memory->create(vols,nprocs,"create_particles:vols");
683   MPI_Allgather(&volupto,1,MPI_DOUBLE,vols,1,MPI_DOUBLE,world);
684 
685   // nme = # of particles for me to create
686   // gathered Scan results not guaranteed to be monotonically increasing
687   // can cause epsilon mis-counts for huge particle counts
688   // enforce that by brute force
689 
690   for (int i = 1; i < nprocs; i++)
691     if (vols[i] != vols[i-1] &&
692         fabs(vols[i]-vols[i-1])/vols[nprocs-1] < EPSZERO)
693       vols[i] = vols[i-1];
694 
695   bigint nstart,nstop;
696   if (me > 0) nstart = static_cast<bigint> (np * (vols[me-1]/vols[nprocs-1]));
697   else nstart = 0;
698   nstop = static_cast<bigint> (np * (vols[me]/vols[nprocs-1]));
699   bigint nme = nstop-nstart;
700 
701   memory->destroy(vols);
702 
703   // nfix_update_custom = # of fixes with update_custom() method
704 
705   modify->list_init_fixes();
706   int nfix_update_custom = modify->n_update_custom;
707 
708   // loop over cells I own
709   // only add particles to OUTSIDE cells
710   // skip cells entirely outside region
711   // ntarget = floating point # of particles to create in one cell
712   // npercell = integer # of particles to create in one cell
713   // basing ntarget on accumulated volume and nprev insures Nme total creations
714   // particle species = random value based on mixture fractions
715   // particle velocity = stream velocity + thermal velocity
716 
717   int *species = particle->mixture[imix]->species;
718   double *cummulative = particle->mixture[imix]->cummulative;
719   double *vstream = particle->mixture[imix]->vstream;
720   double *vscale = particle->mixture[imix]->vscale;
721   int nspecies = particle->mixture[imix]->nspecies;
722   double temp_thermal = particle->mixture[imix]->temp_thermal;
723   double temp_rot = particle->mixture[imix]->temp_rot;
724   double temp_vib = particle->mixture[imix]->temp_vib;
725 
726   int npercell,ncreate,isp,ispecies,id;
727   double x[3],v[3],vstream_variable[3];
728   double ntarget,scale,rn,vn,vr,theta1,theta2,erot,evib;
729 
730   double tempscale = 1.0;
731   double sqrttempscale = 1.0;
732 
733   double volsum = 0.0;
734   bigint nprev = 0;
735 
736   int* ncreate_values;
737   memory->create(ncreate_values, nglocal, "create_particles:ncreate");
738 
739   for (int i = 0; i < nglocal; i++) {
740     if (cinfo[i].type != OUTSIDE) continue;
741     lo = cells[i].lo;
742     hi = cells[i].hi;
743     if (region && region->bboxflag && outside_region(dimension,lo,hi))
744       continue;
745 
746     if (dimension == 3) volone = (hi[0]-lo[0]) * (hi[1]-lo[1]) * (hi[2]-lo[2]);
747     else if (domain->axisymmetric)
748       volone = (hi[0]-lo[0]) * (hi[1]*hi[1]-lo[1]*lo[1])*MY_PI;
749     else volone = (hi[0]-lo[0]) * (hi[1]-lo[1]);
750     volsum += volone / cinfo[i].weight;
751 
752     ntarget = nme * volsum/volme - nprev;
753     npercell = static_cast<int> (ntarget);
754     if (random->uniform() < ntarget-npercell) npercell++;
755     ncreate = npercell;
756 
757     if (densflag) {
758       scale = density_variable(lo,hi);
759       ntarget *= scale;
760       ncreate = static_cast<int> (ntarget);
761       if (random->uniform() < ntarget-ncreate) ncreate++;
762     }
763 
764     ncreate_values[i] = ncreate;
765 
766     // increment count without effect of density variation
767     // so that target insertion count is undisturbed
768 
769     nprev += npercell;
770   }
771 
772   for (int i = 0; i < nglocal; i++) {
773     if (cinfo[i].type != OUTSIDE) continue;
774     lo = cells[i].lo;
775     hi = cells[i].hi;
776     if (region && region->bboxflag && outside_region(dimension,lo,hi))
777       continue;
778     ncreate = ncreate_values[i];
779 
780     for (int m = 0; m < ncreate; m++) {
781       rn = random->uniform();
782       isp = 0;
783       while (cummulative[isp] < rn) isp++;
784       ispecies = species[isp];
785 
786       x[0] = lo[0] + random->uniform() * (hi[0]-lo[0]);
787       x[1] = lo[1] + random->uniform() * (hi[1]-lo[1]);
788       x[2] = lo[2] + random->uniform() * (hi[2]-lo[2]);
789       if (dimension == 2) x[2] = 0.0;
790 
791       if (region && !region->match(x)) continue;
792       if (speciesflag) {
793         isp = species_variable(x) - 1;
794         if (isp < 0 || isp >= nspecies) continue;
795         ispecies = species[isp];
796       }
797 
798       if (tempflag) {
799         tempscale = temperature_variable(x);
800         sqrttempscale = sqrt(tempscale);
801       }
802 
803       vn = vscale[isp] * sqrttempscale * sqrt(-log(random->uniform()));
804       vr = vscale[isp] * sqrttempscale * sqrt(-log(random->uniform()));
805       theta1 = MY_2PI * random->uniform();
806       theta2 = MY_2PI * random->uniform();
807 
808       if (velflag) {
809         velocity_variable(x,vstream,vstream_variable);
810         v[0] = vstream_variable[0] + vn*cos(theta1);
811         v[1] = vstream_variable[1] + vr*cos(theta2);
812         v[2] = vstream_variable[2] + vr*sin(theta2);
813       } else {
814         v[0] = vstream[0] + vn*cos(theta1);
815         v[1] = vstream[1] + vr*cos(theta2);
816         v[2] = vstream[2] + vr*sin(theta2);
817       }
818 
819       erot = particle->erot(ispecies,temp_rot*tempscale,random);
820       evib = particle->evib(ispecies,temp_vib*tempscale,random);
821 
822       id = MAXSMALLINT*random->uniform();
823 
824       particle->add_particle(id,ispecies,i,x,v,erot,evib);
825       if (nfix_update_custom)
826         modify->update_custom(particle->nlocal-1,temp_thermal,
827                              temp_rot,temp_vib,vstream);
828     }
829   }
830 
831   memory->destroy(ncreate_values);
832 
833   delete random;
834 }
835 
836 /* ----------------------------------------------------------------------
837    return 1 if grid cell with lo/hi is entirely outside region bounding box
838    else return 0
839 ------------------------------------------------------------------------- */
840 
outside_region(int dim,double * lo,double * hi)841 int CreateParticles::outside_region(int dim, double *lo, double *hi)
842 {
843   int flag = 1;
844   if (hi[0] > region->extent_xlo &&
845       lo[0] < region->extent_xhi) flag = 0;
846   if (hi[1] > region->extent_ylo &&
847       lo[1] < region->extent_yhi) flag = 0;
848   if (dim == 3) {
849     if (hi[2] > region->extent_zlo &&
850         lo[2] < region->extent_zhi) flag = 0;
851   }
852   return flag;
853 }
854 
855 /* ----------------------------------------------------------------------
856    use particle position in svar variable to generate particle species
857    first plug in particle x,y,z values into sxvar,syvar,szvar
858 ------------------------------------------------------------------------- */
859 
species_variable(double * x)860 int CreateParticles::species_variable(double *x)
861 {
862   if (sxstr) input->variable->internal_set(sxvar,x[0]);
863   if (systr) input->variable->internal_set(syvar,x[1]);
864   if (szstr) input->variable->internal_set(szvar,x[2]);
865 
866   double value = input->variable->compute_equal(svar);
867   int isp = static_cast<int> (value);
868   return isp;
869 }
870 
871 /* ----------------------------------------------------------------------
872    use grid cell center in dvar variable to generate density scale factor
873    first plug in grid x,y,z values into dxvar,dyvar,dzvar
874 ------------------------------------------------------------------------- */
875 
density_variable(double * lo,double * hi)876 double CreateParticles::density_variable(double *lo, double *hi)
877 {
878   double center[3];
879   center[0] = 0.5 * (lo[0]+hi[0]);
880   center[1] = 0.5 * (lo[1]+hi[1]);
881   center[2] = 0.5 * (lo[2]+hi[2]);
882 
883   if (dxstr) input->variable->internal_set(dxvar,center[0]);
884   if (dystr) input->variable->internal_set(dyvar,center[1]);
885   if (dzstr) input->variable->internal_set(dzvar,center[2]);
886 
887   double scale = input->variable->compute_equal(dvar);
888   return scale;
889 }
890 
891 /* ----------------------------------------------------------------------
892    use grid cell center in tvar variable to generate temperature scale factor
893    first plug in grid x,y,z values into txvar,tyvar,tzvar
894 ------------------------------------------------------------------------- */
895 
temperature_variable(double * x)896 double CreateParticles::temperature_variable(double *x)
897 {
898 
899   if (txstr) input->variable->internal_set(txvar,x[0]);
900   if (tystr) input->variable->internal_set(tyvar,x[1]);
901   if (tzstr) input->variable->internal_set(tzvar,x[2]);
902 
903   double scale = input->variable->compute_equal(tvar);
904   return scale;
905 }
906 /* ----------------------------------------------------------------------
907    use particle position in vxvar,vyvar,vzvar variables to generate vel stream
908    first plug in particle x,y,z values into vvarx,vvary,vvarz
909 ------------------------------------------------------------------------- */
910 
velocity_variable(double * x,double * vstream,double * vstream_variable)911 void CreateParticles::velocity_variable(double *x, double *vstream,
912                                         double *vstream_variable)
913 {
914   if (vstrx) input->variable->internal_set(vvarx,x[0]);
915   if (vstry) input->variable->internal_set(vvary,x[1]);
916   if (vstrz) input->variable->internal_set(vvarz,x[2]);
917 
918   if (vxstr) vstream_variable[0] = input->variable->compute_equal(vxvar);
919   else vstream_variable[0] = vstream[0];
920   if (vystr) vstream_variable[1] = input->variable->compute_equal(vyvar);
921   else vstream_variable[1] = vstream[1];
922   if (vzstr) vstream_variable[2] = input->variable->compute_equal(vzvar);
923   else vstream_variable[2] = vstream[2];
924 }
925 
926 /* ----------------------------------------------------------------------
927    create Np particles in serial
928    every proc generates all Np coords, only keeps those in cells it owns
929    created particle attributes should be independent of number of procs
930 ------------------------------------------------------------------------- */
931 
932 /*
933 void CreateParticles::create_all(bigint n)
934 {
935   int dimension = domain->dimension;
936   double xlo = domain->boxlo[0];
937   double ylo = domain->boxlo[1];
938   double zlo = domain->boxlo[2];
939   double xprd = domain->xprd;
940   double yprd = domain->yprd;
941   double zprd = domain->zprd;
942 
943   int me = comm->me;
944   RanKnuth *random = new RandomPark(update->ranmaster->uniform());
945 
946   int icell,id;
947   double x,y,z;
948 
949   // loop over all N particles
950 
951   for (bigint m = 0; m < n; m++) {
952     x = xlo + random->uniform()*xprd;
953     y = ylo + random->uniform()*yprd;
954     z = zlo + random->uniform()*zprd;
955     if (dimension == 2) z = 0.0;
956 
957     // which_cell() returns global grid cell index the particle is in
958     // if I own that grid cell, add particle
959 
960     icell = grid->which_cell(x,y,z);
961     id = MAXSMALLINT*random->uniform();
962 
963     if (grid->cells[icell].proc == me) {
964       particle->add_particle(id,1,icell,x,y,z);
965     }
966   }
967 
968   delete random;
969 }
970 */
971