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