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 "mpi.h"
16 #include "math.h"
17 #include "string.h"
18 #include "stdlib.h"
19 #include "ctype.h"
20 #include "particle.h"
21 #include "grid.h"
22 #include "update.h"
23 #include "comm.h"
24 #include "mixture.h"
25 #include "collide.h"
26 #include "random_mars.h"
27 #include "random_knuth.h"
28 #include "memory.h"
29 #include "error.h"
30 #include "fix_vibmode.h"
31
32 using namespace SPARTA_NS;
33
34 enum{PKEEP,PINSERT,PDONE,PDISCARD,PENTRY,PEXIT,PSURF}; // several files
35 enum{NONE,DISCRETE,SMOOTH}; // several files
36 enum{INT,DOUBLE}; // several files
37
38 #define DELTA 16384
39 #define DELTASPECIES 16
40 #define DELTAMIXTURE 8
41 #define MAXLINE 1024
42
43 // customize by adding an abbreviation string
44 // also add a check for the keyword in 2 places in add_species()
45
46 #define AIR "N O NO"
47
48 /* ---------------------------------------------------------------------- */
49
Particle(SPARTA * sparta)50 Particle::Particle(SPARTA *sparta) : Pointers(sparta)
51 {
52 MPI_Comm_rank(world,&me);
53
54 exist = sorted = 0;
55 nglobal = 0;
56 nlocal = maxlocal = 0;
57 particles = NULL;
58
59 nspecies = maxspecies = 0;
60 species = NULL;
61 maxvibmode = 0;
62
63 //maxgrid = 0;
64 //cellcount = NULL;
65 //first = NULL;
66 maxsort = 0;
67 next = NULL;
68
69 // create two default mixtures
70
71 nmixture = maxmixture = 0;
72 mixture = NULL;
73
74 char **newarg = new char*[1];
75 newarg[0] = (char *) "all";
76 add_mixture(1,newarg);
77 newarg[0] = (char *) "species";
78 add_mixture(1,newarg);
79 delete [] newarg;
80
81 // custom per-particle vectors/arrays
82
83 ncustom = 0;
84 ename = NULL;
85 etype = esize = ewhich = NULL;
86
87 ncustom_ivec = ncustom_iarray = 0;
88 icustom_ivec = icustom_iarray = NULL;
89 eivec = NULL;
90 eiarray = NULL;
91 eicol = NULL;
92
93 ncustom_dvec = ncustom_darray = 0;
94 icustom_dvec = icustom_darray = NULL;
95 edvec = NULL;
96 edarray = NULL;
97 edcol = NULL;
98
99 custom_restart_flag = NULL;
100
101 // RNG for particle weighting
102
103 wrandom = NULL;
104
105 copy = copymode = 0;
106 }
107
108 /* ---------------------------------------------------------------------- */
109
~Particle()110 Particle::~Particle()
111 {
112 if (copy || copymode) return;
113
114 memory->sfree(species);
115 for (int i = 0; i < nmixture; i++) delete mixture[i];
116 memory->sfree(mixture);
117
118 memory->sfree(particles);
119 //memory->destroy(cellcount);
120 //memory->destroy(first);
121 memory->destroy(next);
122
123 for (int i = 0; i < ncustom; i++) delete [] ename[i];
124 memory->sfree(ename);
125 memory->destroy(etype);
126 memory->destroy(esize);
127 memory->destroy(ewhich);
128
129 for (int i = 0; i < ncustom_ivec; i++) memory->destroy(eivec[i]);
130 for (int i = 0; i < ncustom_iarray; i++) memory->destroy(eiarray[i]);
131 for (int i = 0; i < ncustom_dvec; i++) memory->destroy(edvec[i]);
132 for (int i = 0; i < ncustom_darray; i++) memory->destroy(edarray[i]);
133
134 memory->destroy(icustom_ivec);
135 memory->destroy(icustom_iarray);
136 memory->sfree(eivec);
137 memory->sfree(eiarray);
138 memory->destroy(eicol);
139 memory->destroy(icustom_dvec);
140 memory->destroy(icustom_darray);
141 memory->sfree(edvec);
142 memory->sfree(edarray);
143 memory->destroy(edcol);
144
145 delete wrandom;
146 }
147
148 /* ---------------------------------------------------------------------- */
149
init()150 void Particle::init()
151 {
152 // check for errors in custom particle vectors/arrays
153
154 error_custom();
155
156 // initialize mixtures
157
158 for (int i = 0; i < nmixture; i++) mixture[i]->init();
159
160 // RNG for particle weighting
161
162 if (!wrandom) {
163 wrandom = new RanKnuth(update->ranmaster->uniform());
164 double seed = update->ranmaster->uniform();
165 wrandom->reset(seed,me,100);
166 }
167
168 // if first run after reading a restart file,
169 // delete any custom particle attributes that have not been re-defined
170 // use nactive since remove_custom() may alter ncustom
171
172 if (custom_restart_flag) {
173 int nactive = ncustom;
174 for (int i = 0; i < nactive; i++)
175 if (custom_restart_flag[i] == 0) remove_custom(i);
176 delete [] custom_restart_flag;
177 custom_restart_flag = NULL;
178 }
179
180 // if vibstyle = DISCRETE,
181 // all species with vibdof = 4,6,8 must have info read from a species.vib file
182
183 if (collide && collide->vibstyle == DISCRETE) {
184 for (int isp = 0; isp < nspecies; isp++) {
185 if (species[isp].vibdof <= 2) continue;
186 if (species[isp].vibdiscrete_read == 0) {
187 char str[128];
188 sprintf(str,"Discrete vibrational info for species %s not read in",
189 species[isp].id);
190 error->all(FLERR,str);
191 }
192 }
193 }
194
195 // reallocate cellcount and first lists as needed
196 // NOTE: when grid becomes dynamic, will need to do this in sort()
197
198 //if (maxgrid < grid->nlocal) {
199 // maxgrid = grid->nlocal;
200 // memory->destroy(cellcount);
201 //memory->destroy(first);
202 //memory->create(first,maxgrid,"particle:first");
203 //memory->create(cellcount,maxgrid,"particle:cellcount");
204 // }
205 }
206
207 /* ----------------------------------------------------------------------
208 compress particle list to remove particles with indices in mlist
209 mlist indices MUST be in ascending order
210 overwrite deleted particle with particle from end of nlocal list
211 inner while loop avoids overwrite with deleted particle at end of mlist
212 called every step from Comm::migrate_particles() when particles migrate
213 this is similar to compress_reactions(), but does not need
214 an auxiliary vector b/c indices are in ascending order
215 this does NOT preserve particle sorting
216 ------------------------------------------------------------------------- */
217
compress_migrate(int nmigrate,int * mlist)218 void Particle::compress_migrate(int nmigrate, int *mlist)
219 {
220 int i,j,k;
221 int nbytes = sizeof(OnePart);
222
223 if (!ncustom) {
224 for (i = 0; i < nmigrate; i++) {
225 j = mlist[i];
226 k = nlocal - 1;
227 while (k == mlist[nmigrate-1] && k > j) {
228 nmigrate--;
229 nlocal--;
230 k--;
231 }
232 nlocal--;
233 if (j == k) continue;
234 memcpy(&particles[j],&particles[k],nbytes);
235 }
236
237 } else {
238 for (i = 0; i < nmigrate; i++) {
239 j = mlist[i];
240 k = nlocal - 1;
241 while (k == mlist[nmigrate-1] && k > j) {
242 nmigrate--;
243 nlocal--;
244 k--;
245 }
246 nlocal--;
247 if (j == k) continue;
248 memcpy(&particles[j],&particles[k],nbytes);
249 copy_custom(j,k);
250 }
251 }
252
253 sorted = 0;
254 }
255
256 /* ----------------------------------------------------------------------
257 compress particle list to remove particles with icell < 0
258 all particles MUST be in owned cells
259 overwrite deleted particle with particle from end of nlocal list
260 called from Comm::migrate_cells() when cells+particles migrate on rebalance
261 called from AdaptGrid when coarsening occurs
262 called from ReadSurf to remove particles from cells with surfs
263 this does NOT preserve particle sorting
264 ------------------------------------------------------------------------- */
265
compress_rebalance()266 void Particle::compress_rebalance()
267 {
268 int nbytes = sizeof(OnePart);
269
270 if (!ncustom) {
271 int i = 0;
272 while (i < nlocal) {
273 if (particles[i].icell < 0) {
274 memcpy(&particles[i],&particles[nlocal-1],nbytes);
275 nlocal--;
276 } else i++;
277 }
278
279 } else {
280 int i = 0;
281 while (i < nlocal) {
282 if (particles[i].icell < 0) {
283 memcpy(&particles[i],&particles[nlocal-1],nbytes);
284 copy_custom(i,nlocal-1);
285 nlocal--;
286 } else i++;
287 }
288 }
289
290 sorted = 0;
291 }
292
293 /* ----------------------------------------------------------------------
294 compress particle list to remove particles with icell < 0
295 same as compress_rebalance() except this DOES preserve particle sorting
296 invoked by balance migrate_cells_less_memory()
297 ------------------------------------------------------------------------- */
298
compress_rebalance_sorted()299 void Particle::compress_rebalance_sorted()
300 {
301 int nbytes = sizeof(OnePart);
302
303 Grid::ChildInfo *cinfo = grid->cinfo;
304
305 if (!ncustom) {
306 int i = 0;
307 while (i < nlocal) {
308 if (particles[i].icell < 0) {
309 int icell = particles[nlocal-1].icell;
310 if (icell >= 0) {
311 if (cinfo[icell].first == nlocal-1) cinfo[icell].first = i;
312 else {
313 int ip = cinfo[icell].first;
314 while (ip >= 0) {
315 if (next[ip] == nlocal-1) {
316 next[ip] = i;
317 break;
318 }
319 ip = next[ip];
320 }
321 }
322 next[i] = next[nlocal-1];
323 }
324 memcpy(&particles[i],&particles[nlocal-1],nbytes);
325 nlocal--;
326 } else i++;
327 }
328
329 } else {
330 int i = 0;
331 while (i < nlocal) {
332 if (particles[i].icell < 0) {
333 int icell = particles[nlocal-1].icell;
334 if (icell >= 0) {
335 if (cinfo[icell].first == nlocal-1) cinfo[icell].first = i;
336 else {
337 int ip = cinfo[icell].first;
338 while (ip >= 0) {
339 if (next[ip] == nlocal-1) {
340 next[ip] = i;
341 break;
342 }
343 ip = next[ip];
344 }
345 }
346 next[i] = next[nlocal-1];
347 }
348 memcpy(&particles[i],&particles[nlocal-1],nbytes);
349 copy_custom(i,nlocal-1);
350 nlocal--;
351 } else i++;
352 }
353 }
354 }
355
356 /* ----------------------------------------------------------------------
357 compress particle list to remove particles with indices in dellist
358 dellist indices can be in ANY order
359 overwrite deleted particle with particle from end of nlocal list
360 use of next vector does bookkeeping for particles
361 that are moved from their original location before they are deleted
362 called from Collide::migrate_particles() each timestep
363 if any particles were deleted by gas-phase collision reactions
364 this is similar to compress_migrate(), but needs to use
365 an auxiliary vector b/c indices are in random order
366 ------------------------------------------------------------------------- */
367
compress_reactions(int ndelete,int * dellist)368 void Particle::compress_reactions(int ndelete, int *dellist)
369 {
370 int i;
371
372 int nbytes = sizeof(OnePart);
373
374 // reallocate next list as needed
375
376 if (maxsort < maxlocal) {
377 maxsort = maxlocal;
378 memory->destroy(next);
379 memory->create(next,maxsort,"particle:next");
380 }
381
382 // use next as a scratch vector
383 // is always defined when performing collisions and thus gas reactions
384 // next is only used for upper locs from nlocal-ndelete to nlocal
385 // next[i] = current index of atom originally at index i, when i >= nlocal
386 // next[i] = original index of atom currently at index i, when i <= nlocal
387
388 int upper = nlocal-ndelete;
389 for (i = upper; i < nlocal; i++) next[i] = i;
390
391 // i = current index of atom to remove, even if it previously moved
392
393 if (!ncustom) {
394 for (int m = 0; m < ndelete; m++) {
395 i = dellist[m];
396 if (i >= nlocal) i = next[i];
397 nlocal--;
398 if (i == nlocal) continue;
399 memcpy(&particles[i],&particles[nlocal],nbytes);
400 if (i >= upper) next[i] = next[nlocal];
401 next[next[nlocal]] = i;
402 }
403
404 } else {
405 for (int m = 0; m < ndelete; m++) {
406 i = dellist[m];
407 if (i >= nlocal) i = next[i];
408 nlocal--;
409 if (i == nlocal) continue;
410 memcpy(&particles[i],&particles[nlocal],nbytes);
411 copy_custom(i,nlocal);
412 if (i >= upper) next[i] = next[nlocal];
413 next[next[nlocal]] = i;
414 }
415 }
416 }
417
418 /* ----------------------------------------------------------------------
419 sort particles into grid cells
420 set cinfo.first = index of first particle in cell
421 set cinfo.count = # of particles in cell
422 next[] = index of next particle in same cell, -1 for no more
423 ------------------------------------------------------------------------- */
424
sort()425 void Particle::sort()
426 {
427 sorted = 1;
428
429 // reallocate next list as needed
430 // NOTE: why not just compare maxsort to nlocal?
431 // then could realloc less often?
432 // ditto for all reallocs of next in related methods
433
434 if (maxsort < maxlocal) {
435 maxsort = maxlocal;
436 memory->destroy(next);
437 memory->create(next,maxsort,"particle:next");
438 }
439
440 // initialize linked list of particles in cells I own
441
442 Grid::ChildInfo *cinfo = grid->cinfo;
443 int nglocal = grid->nlocal;
444
445 for (int icell = 0; icell < nglocal; icell++) {
446 cinfo[icell].first = -1;
447 cinfo[icell].count = 0;
448 }
449
450 // reverse loop over partlcles to store linked lists in forward order
451 // icell = global cell the particle is in
452
453 int icell;
454 for (int i = nlocal-1; i >= 0; i--) {
455 icell = particles[i].icell;
456 next[i] = cinfo[icell].first;
457 cinfo[icell].first = i;
458 cinfo[icell].count++;
459 }
460 }
461
462 /* ----------------------------------------------------------------------
463 reallocate next list if necessary
464 called before partial sort by FixEmit classes in subsonic case
465 ------------------------------------------------------------------------- */
466
sort_allocate()467 void Particle::sort_allocate()
468 {
469 if (maxsort < maxlocal) {
470 maxsort = maxlocal;
471 memory->destroy(next);
472 memory->create(next,maxsort,"particle:next");
473 }
474 }
475
476 /* ----------------------------------------------------------------------
477 mark all particles in a grid cell, stating with index ip, with icell = -1
478 so can be deleted by compress_rebalance()
479 assume particles are already sorted
480 called from ReadSurf to remove particles from a grid cell with surfaces
481 ------------------------------------------------------------------------- */
482
remove_all_from_cell(int ip)483 void Particle::remove_all_from_cell(int ip)
484 {
485 while (ip >= 0) {
486 particles[ip].icell = -1;
487 ip = next[ip];
488 }
489 }
490
491 /* ----------------------------------------------------------------------
492 set the initial weight of each particle
493 called by Update before particle move
494 only called if particle weighting is enabled
495 only grid-based weighting is currently implemented
496 ------------------------------------------------------------------------- */
497
pre_weight()498 void Particle::pre_weight()
499 {
500 int icell;
501 Grid::ChildInfo *cinfo = grid->cinfo;
502
503 for (int i = 0; i < nlocal; i++) {
504 icell = particles[i].icell;
505 particles[i].weight = cinfo[icell].weight;
506 }
507 }
508
509 /* ----------------------------------------------------------------------
510 clone/delete each particle based on ratio of its initial/final weights
511 called by Update after particle move and migration
512 only called if particle weighting is enabled
513 only grid-based weighting is currently implemented
514 ------------------------------------------------------------------------- */
515
post_weight()516 void Particle::post_weight()
517 {
518 int m,icell,nclone;
519 double ratio,fraction;
520
521 int nbytes = sizeof(OnePart);
522 Grid::ChildInfo *cinfo = grid->cinfo;
523
524 // nlocal_original-1 = index of last original particle
525
526 int nlocal_original = nlocal;
527 int i = 0;
528
529 while (i < nlocal_original) {
530 icell = particles[i].icell;
531
532 // next particle will be an original particle
533 // skip it if no weight change
534
535 if (particles[i].weight == cinfo[icell].weight) {
536 i++;
537 continue;
538 }
539
540 // ratio < 1.0 is candidate for deletion
541 // if deleted and particle that takes its place is cloned (Nloc > Norig)
542 // then skip it via i++, else will examine it on next iteration
543
544 ratio = particles[i].weight / cinfo[icell].weight;
545
546 if (ratio < 1.0) {
547 if (wrandom->uniform() > ratio) {
548 memcpy(&particles[i],&particles[nlocal-1],nbytes);
549 if (ncustom) copy_custom(i,nlocal-1);
550 if (nlocal > nlocal_original) i++;
551 else nlocal_original--;
552 nlocal--;
553 } else i++;
554 continue;
555 }
556
557 // ratio > 1.0 is candidate for cloning
558 // create Nclone new particles each with unique ID
559
560 nclone = static_cast<int> (ratio);
561 fraction = ratio - nclone;
562 nclone--;
563 if (wrandom->uniform() < fraction) nclone++;
564
565 for (m = 0; m < nclone; m++) {
566 clone_particle(i);
567 particles[nlocal-1].id = MAXSMALLINT*wrandom->uniform();
568 }
569 i++;
570 }
571 }
572
573 /* ----------------------------------------------------------------------
574 insure particle list can hold nextra new particles
575 if defined, also grow custom particle arrays and initialize with zeroes
576 ------------------------------------------------------------------------- */
577
grow(int nextra)578 void Particle::grow(int nextra)
579 {
580 bigint target = (bigint) nlocal + nextra;
581 if (target <= maxlocal) return;
582
583 int oldmax = maxlocal;
584 bigint newmax = maxlocal;
585 while (newmax < target) newmax += DELTA;
586
587 if (newmax > MAXSMALLINT)
588 error->one(FLERR,"Per-processor particle count is too big");
589
590 maxlocal = newmax;
591 particles = (OnePart *)
592 memory->srealloc(particles,maxlocal*sizeof(OnePart),
593 "particle:particles");
594 memset(&particles[oldmax],0,(maxlocal-oldmax)*sizeof(OnePart));
595
596 if (ncustom == 0) return;
597
598 for (int i = 0; i < ncustom; i++) {
599 if (ename[i] == NULL) continue;
600 grow_custom(i,oldmax,maxlocal);
601 }
602 }
603
604 /* ----------------------------------------------------------------------
605 insure species list can hold maxspecies species
606 assumes that maxspecies has already been increased
607 ------------------------------------------------------------------------- */
608
grow_species()609 void Particle::grow_species()
610 {
611 species = (Species *)
612 memory->srealloc(species,maxspecies*sizeof(Species),"particle:species");
613 }
614
615 /* ----------------------------------------------------------------------
616 grow next list if more particles now exist than there is room for
617 called from Grid::unpack_particles_adapt() when grid adaptation
618 takes place and acquire particles from other procs due to coarsening
619 unlike sort(), this requires next list be grown, not destroy/create
620 b/c sorted particle list is maintained during adaptation
621 ------------------------------------------------------------------------- */
622
grow_next()623 void Particle::grow_next()
624 {
625 // compare maxsort (length of next) to new particle count (nlocal)
626 // grow to maxlocal (max length of particles) to avoid frequent re-grow
627
628 if (maxsort < nlocal) {
629 maxsort = maxlocal;
630 memory->grow(next,maxsort,"particle:next");
631 }
632 }
633
634 /* ----------------------------------------------------------------------
635 add a particle to particle list
636 return 1 if particle array was reallocated, else 0
637 ------------------------------------------------------------------------- */
638
add_particle(int id,int ispecies,int icell,double * x,double * v,double erot,double evib)639 int Particle::add_particle(int id, int ispecies, int icell,
640 double *x, double *v, double erot, double evib)
641 {
642 int reallocflag = 0;
643 if (nlocal == maxlocal) {
644 grow(1);
645 reallocflag = 1;
646 }
647
648 OnePart *p = &particles[nlocal];
649
650 p->id = id;
651 p->ispecies = ispecies;
652 p->icell = icell;
653 p->x[0] = x[0];
654 p->x[1] = x[1];
655 p->x[2] = x[2];
656 p->v[0] = v[0];
657 p->v[1] = v[1];
658 p->v[2] = v[2];
659 p->erot = erot;
660 p->evib = evib;
661 p->flag = PKEEP;
662
663 //p->dtremain = 0.0; not needed due to memset in grow() ??
664 //p->weight = 1.0; not needed due to memset in grow() ??
665
666 nlocal++;
667 return reallocflag;
668 }
669
670 /* ----------------------------------------------------------------------
671 add an empty particle to particle list, caller will fill it
672 return 1 if particle array was reallocated, else 0
673 ------------------------------------------------------------------------- */
674
add_particle()675 int Particle::add_particle()
676 {
677 int reallocflag = 0;
678 if (nlocal == maxlocal) {
679 grow(1);
680 reallocflag = 1;
681 }
682
683 nlocal++;
684 return reallocflag;
685 }
686
687 /* ----------------------------------------------------------------------
688 clone particle Index and add to end of particle list
689 return 1 if particle array was reallocated, else 0
690 ------------------------------------------------------------------------- */
691
clone_particle(int index)692 int Particle::clone_particle(int index)
693 {
694 int reallocflag = 0;
695 if (nlocal == maxlocal) {
696 grow(1);
697 reallocflag = 1;
698 }
699
700 memcpy(&particles[nlocal],&particles[index],sizeof(OnePart));
701 if (ncustom) copy_custom(nlocal,index);
702
703 nlocal++;
704 return reallocflag;
705 }
706
707 /* ----------------------------------------------------------------------
708 add one or more species to species list
709 ------------------------------------------------------------------------- */
710
add_species(int narg,char ** arg)711 void Particle::add_species(int narg, char **arg)
712 {
713 int i,j,k,n;;
714
715 if (narg < 2) error->all(FLERR,"Illegal species command");
716
717 if (me == 0) {
718 fp = fopen(arg[0],"r");
719 if (fp == NULL) {
720 char str[128];
721 sprintf(str,"Cannot open species file %s",arg[0]);
722 error->one(FLERR,str);
723 }
724 }
725
726 // nfile = # of species defined in file
727 // filespecies = list of species defined in file
728
729 nfile = maxfile = 0;
730 filespecies = NULL;
731
732 if (me == 0) read_species_file();
733 MPI_Bcast(&nfile,1,MPI_INT,0,world);
734 if (comm->me) {
735 filespecies = (Species *)
736 memory->smalloc(nfile*sizeof(Species),"particle:filespecies");
737 }
738 MPI_Bcast(filespecies,nfile*sizeof(Species),MPI_BYTE,0,world);
739
740 // newspecies = # of new user-requested species,
741 // not including trailing optional keywords
742 // names = list of new species IDs
743 // customize abbreviations by adding new keyword in 2 places
744
745 char line[MAXLINE];
746
747 int newspecies = 0;
748 int iarg = 1;
749 while (iarg < narg) {
750 if (strcmp(arg[iarg],"air") == 0) {
751 strcpy(line,AIR);
752 newspecies += wordcount(line,NULL);
753 } else if (strcmp(arg[iarg],"rotfile") == 0) {
754 break;
755 } else if (strcmp(arg[iarg],"vibfile") == 0) {
756 break;
757 } else {
758 newspecies++;
759 }
760 iarg++;
761 }
762
763 char **names = new char*[newspecies];
764 newspecies = 0;
765
766 iarg = 1;
767 while (iarg < narg) {
768 if (strcmp(arg[iarg],"air") == 0) {
769 strcpy(line,AIR);
770 newspecies += wordcount(line,&names[newspecies]);
771 } else if (strcmp(arg[iarg],"rotfile") == 0) {
772 break;
773 } else if (strcmp(arg[iarg],"vibfile") == 0) {
774 break;
775 } else {
776 names[newspecies++] = arg[iarg];
777 }
778 iarg++;
779 }
780
781 // species ID must be all alphanumeric chars, underscore, plus/minus
782
783 for (i = 0; i < newspecies; i++) {
784 n = strlen(names[i]);
785 for (j = 0; j < n-1; j++)
786 if (!isalnum(names[i][j]) && names[i][j] != '_' &&
787 names[i][j] != '+' && names[i][j] != '-')
788 error->all(FLERR,"Invalid character in species ID");
789 }
790
791 // extend species list if necessary
792
793 if (nspecies + newspecies > maxspecies) {
794 while (nspecies+newspecies > maxspecies) maxspecies += DELTASPECIES;
795 grow_species();
796 }
797
798 // extract info on user-requested species from file species list
799 // add new species to default mixtures "all" and "species"
800
801 int nspecies_original = nspecies;
802
803 int imix_all = find_mixture((char *) "all");
804 int imix_species = find_mixture((char *) "species");
805
806 for (i = 0; i < newspecies; i++) {
807 for (j = 0; j < nspecies; j++)
808 if (strcmp(names[i],species[j].id) == 0) break;
809 if (j < nspecies) error->all(FLERR,"Species ID is already defined");
810 for (j = 0; j < nfile; j++)
811 if (strcmp(names[i],filespecies[j].id) == 0) break;
812 if (j == nfile)
813 error->all(FLERR,"Species ID does not appear in species file");
814 memcpy(&species[nspecies],&filespecies[j],sizeof(Species));
815 nspecies++;
816
817 mixture[imix_all]->add_species_default(species[nspecies-1].id);
818 mixture[imix_species]->add_species_default(species[nspecies-1].id);
819 }
820
821 memory->sfree(filespecies);
822
823 // process any optional keywords
824 // NOTE: rotfile is not yet supported
825
826 int rotindex = 0;
827 int vibindex = 0;
828
829 while (iarg < narg) {
830 if (strcmp(arg[iarg],"rotfile") == 0) {
831 // not yet supported
832 error->all(FLERR,"Illegal species command");
833 if (iarg+2 > narg) error->all(FLERR,"Illegal species command");
834 if (rotindex)
835 error->all(FLERR,"Species command can only use a single rotfile");
836 rotindex = iarg+1;
837 iarg += 2;
838 } else if (strcmp(arg[iarg],"vibfile") == 0) {
839 if (iarg+2 > narg) error->all(FLERR,"Illegal species command");
840 if (vibindex)
841 error->all(FLERR,"Species command can only use a single vibfile");
842 vibindex = iarg+1;
843 iarg += 2;
844 } else error->all(FLERR,"Illegal species command");
845 }
846
847 // read rotational species file and setup per-species params
848
849 if (rotindex) {
850 if (me == 0) {
851 fp = fopen(arg[rotindex],"r");
852 if (fp == NULL) {
853 char str[128];
854 sprintf(str,"Cannot open rotation file %s",arg[rotindex]);
855 error->one(FLERR,str);
856 }
857 }
858
859 nfile = maxfile = 0;
860 filerot = NULL;
861
862 if (me == 0) read_rotation_file();
863 MPI_Bcast(&nfile,1,MPI_INT,0,world);
864 if (comm->me) {
865 filerot = (RotFile *)
866 memory->smalloc(nfile*sizeof(RotFile),"particle:filerot");
867 }
868 MPI_Bcast(filerot,nfile*sizeof(RotFile),MPI_BYTE,0,world);
869
870 for (i = 0; i < newspecies; i++) {
871 int ii = nspecies_original + i;
872
873 for (j = 0; j < nfile; j++)
874 if (strcmp(names[i],filerot[j].id) == 0) break;
875 if (j == nfile) {
876 if (species[ii].rotdof == 0) continue;
877 error->all(FLERR,"Species ID does not appear in rotation file");
878 }
879
880 int ntemp = filerot[j].ntemp;
881 if ((species[ii].rotdof == 0) ||
882 (species[ii].rotdof == 2 && ntemp != 1) ||
883 (species[ii].rotdof == 3 && ntemp != 3))
884 error->all(FLERR,"Mismatch between species rotdof "
885 "and rotation file entry");
886
887 species[ii].nrottemp = ntemp;
888 for (k = 0; k < ntemp; k++) {
889 species[ii].rottemp[k] = filerot[j].rottemp[k];
890 }
891 }
892
893 memory->sfree(filerot);
894 }
895
896 // read vibrational species file and setup per-species params
897
898 if (vibindex) {
899 if (me == 0) {
900 fp = fopen(arg[vibindex],"r");
901 if (fp == NULL) {
902 char str[128];
903 sprintf(str,"Cannot open vibration file %s",arg[vibindex]);
904 error->one(FLERR,str);
905 }
906 }
907
908 nfile = maxfile = 0;
909 filevib = NULL;
910
911 if (me == 0) read_vibration_file();
912 MPI_Bcast(&nfile,1,MPI_INT,0,world);
913 if (comm->me) {
914 filevib = (VibFile *)
915 memory->smalloc(nfile*sizeof(VibFile),"particle:filevib");
916 }
917 MPI_Bcast(filevib,nfile*sizeof(VibFile),MPI_BYTE,0,world);
918
919 for (i = 0; i < newspecies; i++) {
920 int ii = nspecies_original + i;
921
922 for (j = 0; j < nfile; j++)
923 if (strcmp(names[i],filevib[j].id) == 0) break;
924 if (j == nfile) {
925 if (species[ii].vibdof <= 2) continue;
926 error->all(FLERR,"Species ID does not appear in vibration file");
927 }
928
929 int nmode = filevib[j].nmode;
930 if (species[ii].nvibmode != nmode)
931 error->all(FLERR,"Mismatch between species vibdof "
932 "and vibration file entry");
933
934 species[ii].nvibmode = nmode;
935 for (k = 0; k < nmode; k++) {
936 species[ii].vibtemp[k] = filevib[j].vibtemp[k];
937 species[ii].vibrel[k] = filevib[j].vibrel[k];
938 species[ii].vibdegen[k] = filevib[j].vibdegen[k];
939 }
940
941 maxvibmode = MAX(maxvibmode,species[ii].nvibmode);
942 species[ii].vibdiscrete_read = 1;
943 }
944
945 memory->sfree(filevib);
946 }
947
948 // clean up
949
950 delete [] names;
951 }
952
953 /* ----------------------------------------------------------------------
954 create or augment a mixture of species
955 ------------------------------------------------------------------------- */
956
add_mixture(int narg,char ** arg)957 void Particle::add_mixture(int narg, char **arg)
958 {
959 if (narg < 1) error->all(FLERR,"Illegal mixture command");
960
961 // imix = index if mixture ID already exists
962 // else instantiate a new mixture
963
964 int imix = find_mixture(arg[0]);
965
966 if (imix < 0) {
967 if (nmixture == maxmixture) {
968 maxmixture += DELTAMIXTURE;
969 mixture = (Mixture **) memory->srealloc(mixture,
970 maxmixture*sizeof(Mixture *),
971 "particle:mixture");
972 }
973 imix = nmixture;
974 mixture[nmixture++] = new Mixture(sparta,arg[0]);
975 }
976
977 mixture[imix]->command(narg,arg);
978 }
979
980 /* ----------------------------------------------------------------------
981 return index of ID in list of species IDs
982 return -1 if not found
983 ------------------------------------------------------------------------- */
984
find_species(char * id)985 int Particle::find_species(char *id)
986 {
987 for (int i = 0; i < nspecies; i++)
988 if (strcmp(id,species[i].id) == 0) return i;
989 return -1;
990 }
991
992 /* ----------------------------------------------------------------------
993 return index of ID in list of mixture IDs
994 return -1 if not found
995 ------------------------------------------------------------------------- */
996
find_mixture(char * id)997 int Particle::find_mixture(char *id)
998 {
999 for (int i = 0; i < nmixture; i++)
1000 if (strcmp(id,mixture[i]->id) == 0) return i;
1001 return -1;
1002 }
1003
1004 /* ----------------------------------------------------------------------
1005 generate random rotational energy for a particle
1006 only a function of species index and species properties
1007 ------------------------------------------------------------------------- */
1008
erot(int isp,double temp_thermal,RanKnuth * erandom)1009 double Particle::erot(int isp, double temp_thermal, RanKnuth *erandom)
1010 {
1011 double eng,a,erm,b;
1012 int rotstyle = NONE;
1013 if (collide) rotstyle = collide->rotstyle;
1014
1015 if (!collide || collide->rotstyle == NONE) return 0.0;
1016 if (species[isp].rotdof < 2) return 0.0;
1017
1018 if (rotstyle == DISCRETE && species[isp].rotdof == 2) {
1019 int irot = -log(erandom->uniform()) * temp_thermal /
1020 particle->species[isp].rottemp[0];
1021 eng = irot * update->boltz * particle->species[isp].rottemp[0];
1022 } else if (rotstyle == SMOOTH && species[isp].rotdof == 2) {
1023 eng = -log(erandom->uniform()) * update->boltz * temp_thermal;
1024 } else {
1025 a = 0.5*particle->species[isp].rotdof-1.0;
1026 while (1) {
1027 // energy cut-off at 10 kT
1028 erm = 10.0*erandom->uniform();
1029 b = pow(erm/a,a) * exp(a-erm);
1030 if (b > erandom->uniform()) break;
1031 }
1032 eng = erm * update->boltz * temp_thermal;
1033 }
1034
1035 return eng;
1036 }
1037
1038 /* ----------------------------------------------------------------------
1039 generate random vibrational energy for a particle
1040 only a function of species index and species properties
1041 index_vibmode = index of extra per-particle vibrational mode storage
1042 -1 if not defined for this model
1043 ------------------------------------------------------------------------- */
1044
evib(int isp,double temp_thermal,RanKnuth * erandom)1045 double Particle::evib(int isp, double temp_thermal, RanKnuth *erandom)
1046 {
1047 double eng,a,erm,b;
1048
1049 int vibstyle = NONE;
1050 if (collide) vibstyle = collide->vibstyle;
1051 if (vibstyle == NONE || species[isp].vibdof < 2) return 0.0;
1052
1053 // for DISCRETE, only need set evib for vibdof = 2
1054 // mode levels and evib will be set by FixVibmode::update_custom()
1055
1056 eng = 0.0;
1057
1058 if (vibstyle == DISCRETE && species[isp].vibdof == 2) {
1059 int ivib = -log(erandom->uniform()) * temp_thermal /
1060 particle->species[isp].vibtemp[0];
1061 eng = ivib * update->boltz * particle->species[isp].vibtemp[0];
1062 } else if (vibstyle == SMOOTH || species[isp].vibdof >= 2) {
1063 if (species[isp].vibdof == 2)
1064 eng = -log(erandom->uniform()) * update->boltz * temp_thermal;
1065 else if (species[isp].vibdof > 2) {
1066 a = 0.5*particle->species[isp].vibdof-1.;
1067 while (1) {
1068 // energy cut-off at 10 kT
1069 erm = 10.0*erandom->uniform();
1070 b = pow(erm/a,a) * exp(a-erm);
1071 if (b > erandom->uniform()) break;
1072 }
1073 eng = erm * update->boltz * temp_thermal;
1074 }
1075 }
1076
1077 return eng;
1078 }
1079
1080 /* ----------------------------------------------------------------------
1081 read list of species defined in species file
1082 store info in filespecies and nfile
1083 only invoked by proc 0
1084 ------------------------------------------------------------------------- */
1085
read_species_file()1086 void Particle::read_species_file()
1087 {
1088 // read file line by line
1089 // skip blank lines or comment lines starting with '#'
1090 // all other lines must have NWORDS
1091
1092 int NWORDS = 10;
1093 char **words = new char*[NWORDS];
1094 char line[MAXLINE],copy[MAXLINE];
1095
1096 while (fgets(line,MAXLINE,fp)) {
1097 int pre = strspn(line," \t\n\r");
1098 if (pre == strlen(line) || line[pre] == '#') continue;
1099
1100 strcpy(copy,line);
1101 int nwords = wordcount(copy,NULL);
1102 if (nwords != NWORDS)
1103 error->one(FLERR,"Incorrect line format in species file");
1104
1105 if (nfile == maxfile) {
1106 maxfile += DELTASPECIES;
1107 filespecies = (Species *)
1108 memory->srealloc(filespecies,maxfile*sizeof(Species),
1109 "particle:filespecies");
1110 memset(&filespecies[nfile],0,(maxfile-nfile)*sizeof(Species));
1111 }
1112
1113 nwords = wordcount(line,words);
1114 Species *fsp = &filespecies[nfile];
1115
1116 if (strlen(words[0]) + 1 > 16)
1117 error->one(FLERR,"Invalid species ID in species file");
1118 strcpy(fsp->id,words[0]);
1119
1120 fsp->molwt = atof(words[1]);
1121 fsp->mass = atof(words[2]);
1122 fsp->rotdof = atoi(words[3]);
1123 fsp->rotrel = atof(words[4]);
1124 fsp->vibdof = atoi(words[5]);
1125 fsp->vibrel[0] = atof(words[6]);
1126 fsp->vibtemp[0] = atof(words[7]);
1127 fsp->specwt = atof(words[8]);
1128 fsp->charge = atof(words[9]);
1129
1130 if (fsp->rotdof > 0 || fsp->vibdof > 0) fsp->internaldof = 1;
1131 else fsp->internaldof = 0;
1132
1133 // error checks
1134 // NOTE: allow rotdof = 3 when implement rotate = DISCRETE
1135
1136 if (fsp->rotdof != 0 && fsp->rotdof != 2)
1137 error->all(FLERR,"Invalid rotational DOF in species file");
1138 //if (fsp->rotdof != 0 && fsp->rotdof != 2 && fsp->rotdof != 3)
1139 // error->all(FLERR,"Invalid rotational DOF in species file");
1140
1141 if (fsp->vibdof < 0 || fsp->vibdof > 8 || fsp->vibdof % 2)
1142 error->all(FLERR,"Invalid vibrational DOF in species file");
1143
1144 // initialize additional rotation/vibration fields
1145 // may be overwritten by rotfile or vibfile
1146
1147 fsp->nrottemp = 0;
1148 fsp->nvibmode = fsp->vibdof / 2;
1149
1150 fsp->rottemp[0] = fsp->rottemp[1] = fsp->rottemp[2] = 0.0;
1151 fsp->vibtemp[1] = fsp->vibtemp[2] = fsp->vibtemp[3] = 0.0;
1152 fsp->vibrel[1] = fsp->vibrel[2] = fsp->vibrel[3] = 0.0;
1153 fsp->vibdegen[0] = fsp->vibdegen[1] =
1154 fsp->vibdegen[2] = fsp->vibdegen[3] = 0;
1155
1156 fsp->vibdiscrete_read = 0;
1157
1158 nfile++;
1159 }
1160
1161 delete [] words;
1162
1163 fclose(fp);
1164 }
1165
1166 /* ----------------------------------------------------------------------
1167 read list of extra rotation info in rotation file
1168 store info in filerot and nfile
1169 only invoked by proc 0
1170 ------------------------------------------------------------------------- */
1171
read_rotation_file()1172 void Particle::read_rotation_file()
1173 {
1174 // read file line by line
1175 // skip blank lines or comment lines starting with '#'
1176 // all other lines can have up to NWORDS
1177
1178 int NWORDS = 5;
1179 char **words = new char*[NWORDS];
1180 char line[MAXLINE],copy[MAXLINE];
1181
1182 while (fgets(line,MAXLINE,fp)) {
1183 int pre = strspn(line," \t\n\r");
1184 if (pre == strlen(line) || line[pre] == '#') continue;
1185
1186 strcpy(copy,line);
1187 int nwords = wordcount(copy,NULL);
1188 if (nwords > NWORDS)
1189 error->one(FLERR,"Incorrect line format in rotation file");
1190
1191 if (nfile == maxfile) {
1192 maxfile += DELTASPECIES;
1193 filerot = (RotFile *)
1194 memory->srealloc(filerot,maxfile*sizeof(RotFile),
1195 "particle:filerot");
1196 memset(&filerot[nfile],0,(maxfile-nfile)*sizeof(RotFile));
1197 }
1198
1199 nwords = wordcount(line,words);
1200 RotFile *rsp = &filerot[nfile];
1201
1202 if (strlen(words[0]) + 1 > 16)
1203 error->one(FLERR,"Invalid species ID in rotation file");
1204 strcpy(rsp->id,words[0]);
1205
1206 rsp->ntemp = atoi(words[1]);
1207 if (rsp->ntemp != 1 && rsp->ntemp != 3)
1208 error->one(FLERR,"Invalid N count in rotation file");
1209 if (nwords != 2 + rsp->ntemp)
1210 error->one(FLERR,"Incorrect line format in rotation file");
1211
1212 int j = 2;
1213 for (int i = 0; i < rsp->ntemp; i++)
1214 rsp->rottemp[i] = atof(words[j++]);
1215
1216 nfile++;
1217 }
1218
1219 delete [] words;
1220
1221 fclose(fp);
1222 }
1223
1224 /* ----------------------------------------------------------------------
1225 read list of extra rotation info in vibration file
1226 store info in filevib and nfile
1227 only invoked by proc 0
1228 ------------------------------------------------------------------------- */
1229
read_vibration_file()1230 void Particle::read_vibration_file()
1231 {
1232 // read file line by line
1233 // skip blank lines or comment lines starting with '#'
1234 // all other lines can have up to NWORDS
1235
1236 int NWORDS = 14;
1237 char **words = new char*[NWORDS];
1238 char line[MAXLINE],copy[MAXLINE];
1239
1240 while (fgets(line,MAXLINE,fp)) {
1241 int pre = strspn(line," \t\n\r");
1242 if (pre == strlen(line) || line[pre] == '#') continue;
1243
1244 strcpy(copy,line);
1245 int nwords = wordcount(copy,NULL);
1246 if (nwords > NWORDS)
1247 error->one(FLERR,"Incorrect line format in vibration file");
1248
1249 if (nfile == maxfile) {
1250 maxfile += DELTASPECIES;
1251 filevib = (VibFile *)
1252 memory->srealloc(filevib,maxfile*sizeof(VibFile),
1253 "particle:filevib");
1254 memset(&filevib[nfile],0,(maxfile-nfile)*sizeof(VibFile));
1255 }
1256
1257 nwords = wordcount(line,words);
1258 VibFile *vsp = &filevib[nfile];
1259
1260 if (strlen(words[0]) + 1 > 16)
1261 error->one(FLERR,"Invalid species ID in vibration file");
1262 strcpy(vsp->id,words[0]);
1263
1264 vsp->nmode = atoi(words[1]);
1265 if (vsp->nmode < 2 || vsp->nmode > 4)
1266 error->one(FLERR,"Invalid N count in vibration file");
1267 if (nwords != 2 + 3*vsp->nmode)
1268 error->one(FLERR,"Incorrect line format in vibration file");
1269
1270 int j = 2;
1271 for (int i = 0; i < vsp->nmode; i++) {
1272 vsp->vibtemp[i] = atof(words[j++]);
1273 vsp->vibrel[i] = atof(words[j++]);
1274 vsp->vibdegen[i] = atoi(words[j++]);
1275 }
1276 nfile++;
1277 }
1278
1279 delete [] words;
1280
1281 fclose(fp);
1282 }
1283
1284 /* ----------------------------------------------------------------------
1285 count whitespace-delimited words in line
1286 line will be modified, since strtok() inserts NULLs
1287 if words is non-NULL, store ptr to each word
1288 ------------------------------------------------------------------------- */
1289
wordcount(char * line,char ** words)1290 int Particle::wordcount(char *line, char **words)
1291 {
1292 int nwords = 0;
1293 char *word = strtok(line," \t");
1294
1295 while (word) {
1296 if (words) words[nwords] = word;
1297 nwords++;
1298 word = strtok(NULL," \t");
1299 }
1300
1301 return nwords;
1302 }
1303
1304 /* ----------------------------------------------------------------------
1305 proc 0 writes species info to restart file
1306 ------------------------------------------------------------------------- */
1307
write_restart_species(FILE * fp)1308 void Particle::write_restart_species(FILE *fp)
1309 {
1310 fwrite(&nspecies,sizeof(int),1,fp);
1311 fwrite(species,sizeof(Species),nspecies,fp);
1312 }
1313
1314 /* ----------------------------------------------------------------------
1315 proc 0 reads species info from restart file
1316 bcast to other procs
1317 ------------------------------------------------------------------------- */
1318
read_restart_species(FILE * fp)1319 void Particle::read_restart_species(FILE *fp)
1320 {
1321 if (me == 0) fread(&nspecies,sizeof(int),1,fp);
1322 MPI_Bcast(&nspecies,1,MPI_INT,0,world);
1323
1324 if (nspecies > maxspecies) {
1325 while (nspecies > maxspecies) maxspecies += DELTASPECIES;
1326 grow_species();
1327 }
1328
1329 if (me == 0) fread(species,sizeof(Species),nspecies,fp);
1330 MPI_Bcast(species,nspecies*sizeof(Species),MPI_CHAR,0,world);
1331
1332 maxvibmode = 0;
1333 for (int isp = 0; isp < nspecies; isp++)
1334 maxvibmode = MAX(maxvibmode,species[isp].nvibmode);
1335 }
1336
1337 /* ----------------------------------------------------------------------
1338 proc 0 writes mixture info to restart file
1339 ------------------------------------------------------------------------- */
1340
write_restart_mixture(FILE * fp)1341 void Particle::write_restart_mixture(FILE *fp)
1342 {
1343 fwrite(&nmixture,sizeof(int),1,fp);
1344 for (int i = 0; i < nmixture; i++) {
1345 int n = strlen(mixture[i]->id) + 1;
1346 fwrite(&n,sizeof(int),1,fp);
1347 fwrite(mixture[i]->id,sizeof(char),n,fp);
1348 mixture[i]->write_restart(fp);
1349 }
1350 }
1351
1352 /* ----------------------------------------------------------------------
1353 proc 0 reads mixture info from restart file
1354 bcast to other procs and all procs instantiate series of Mixtures
1355 ------------------------------------------------------------------------- */
1356
read_restart_mixture(FILE * fp)1357 void Particle::read_restart_mixture(FILE *fp)
1358 {
1359 // must first clear existing default mixtures
1360
1361 for (int i = 0; i < nmixture; i++) delete mixture[i];
1362 nmixture = 0;
1363
1364 // now process restart file data
1365
1366 if (me == 0) fread(&nmixture,sizeof(int),1,fp);
1367 MPI_Bcast(&nmixture,1,MPI_INT,0,world);
1368
1369 if (nmixture > maxmixture) {
1370 while (nmixture > maxmixture) maxmixture += DELTAMIXTURE;
1371 mixture = (Mixture **)
1372 memory->srealloc(mixture,maxmixture*sizeof(Mixture *),"particle:mixture");
1373 }
1374
1375 int n;
1376 char *id;
1377
1378 for (int i = 0; i < nmixture; i++) {
1379 if (me == 0) fread(&n,sizeof(int),1,fp);
1380 MPI_Bcast(&n,1,MPI_INT,0,world);
1381 id = new char[n];
1382 if (me == 0) fread(id,sizeof(char),n,fp);
1383 MPI_Bcast(id,n,MPI_CHAR,0,world);
1384 mixture[i] = new Mixture(sparta,id);
1385 mixture[i]->read_restart(fp);
1386 delete [] id;
1387 }
1388 }
1389
1390 /* ----------------------------------------------------------------------
1391 return size of particle restart info for this proc
1392 NOTE: worry about N overflowing int and IROUNDUP ???
1393 ------------------------------------------------------------------------- */
1394
size_restart()1395 int Particle::size_restart()
1396 {
1397 int n = sizeof(int);
1398 n = IROUNDUP(n);
1399 n += nlocal * sizeof(OnePartRestart);
1400 n += nlocal * sizeof_custom();
1401 n = IROUNDUP(n);
1402 return n;
1403 }
1404
1405 /* ----------------------------------------------------------------------
1406 return size of particle restart info for this proc
1407 ------------------------------------------------------------------------- */
1408
size_restart_big()1409 bigint Particle::size_restart_big()
1410 {
1411 bigint n = sizeof(int);
1412 n = BIROUNDUP(n);
1413 n += nlocal * sizeof(OnePartRestart);
1414 n += nlocal * sizeof_custom();
1415 n = BIROUNDUP(n);
1416 return n;
1417 }
1418
1419 /* ----------------------------------------------------------------------
1420 pack my particle info into buf
1421 use OnePartRestart data struct for permanent info and to encode cell ID
1422 include per-particle custom attributes if defined
1423 ------------------------------------------------------------------------- */
1424
pack_restart(char * buf)1425 int Particle::pack_restart(char *buf)
1426 {
1427 Grid::ChildCell *cells = grid->cells;
1428 OnePart *p;
1429 OnePartRestart *pr;
1430 int nbytes_custom = sizeof_custom();
1431
1432 char *ptr = buf;
1433 int *ibuf = (int *) ptr;
1434 *ibuf = nlocal;
1435 ptr += sizeof(int);
1436 ptr = ROUNDUP(ptr);
1437
1438 for (int i = 0; i < nlocal; i++) {
1439 p = &particles[i];
1440 pr = (OnePartRestart *) ptr;
1441 pr->id = p->id;
1442 pr->ispecies = p->ispecies;
1443 pr->icell = cells[p->icell].id;
1444 pr->nsplit = cells[p->icell].nsplit;
1445 pr->x[0] = p->x[0];
1446 pr->x[1] = p->x[1];
1447 pr->x[2] = p->x[2];
1448 pr->v[0] = p->v[0];
1449 pr->v[1] = p->v[1];
1450 pr->v[2] = p->v[2];
1451 pr->erot = p->erot;
1452 pr->evib = p->evib;
1453
1454 ptr += sizeof(OnePartRestart);
1455 if (!ncustom) continue;
1456
1457 pack_custom(i,ptr);
1458 ptr += nbytes_custom;
1459 }
1460
1461 ptr = ROUNDUP(ptr);
1462 return ptr - buf;
1463 }
1464
1465 /* ----------------------------------------------------------------------
1466 pack my particle info into buf
1467 use multiple passes to reduce memory use
1468 use OnePartRestart data struct for permanent info and to encode cell ID
1469 include per-particle custom attributes if defined
1470 NOTE: does not ROUNDUP(ptr) at the end, this is done by caller
1471 ------------------------------------------------------------------------- */
1472
pack_restart(char * buf,int step,int pass)1473 void Particle::pack_restart(char *buf, int step, int pass)
1474 {
1475 Grid::ChildCell *cells = grid->cells;
1476 OnePart *p;
1477 OnePartRestart *pr;
1478 int nbytes_custom = sizeof_custom();
1479
1480 char *ptr = buf;
1481 if (pass == 0) {
1482 int *ibuf = (int *) ptr;
1483 *ibuf = nlocal;
1484 ptr += sizeof(int);
1485 ptr = ROUNDUP(ptr);
1486 }
1487
1488 int start = step*pass;
1489 int end = start+step;
1490 end = MIN(nlocal,end);
1491 for (int i = start; i < end; i++) {
1492 p = &particles[i];
1493 pr = (OnePartRestart *) ptr;
1494 pr->id = p->id;
1495 pr->ispecies = p->ispecies;
1496 pr->icell = cells[p->icell].id;
1497 pr->nsplit = cells[p->icell].nsplit;
1498 pr->x[0] = p->x[0];
1499 pr->x[1] = p->x[1];
1500 pr->x[2] = p->x[2];
1501 pr->v[0] = p->v[0];
1502 pr->v[1] = p->v[1];
1503 pr->v[2] = p->v[2];
1504 pr->erot = p->erot;
1505 pr->evib = p->evib;
1506
1507 ptr += sizeof(OnePartRestart);
1508 if (!ncustom) continue;
1509
1510 pack_custom(i,ptr);
1511 ptr += nbytes_custom;
1512 }
1513 }
1514
1515 /* ----------------------------------------------------------------------
1516 unpack particle info into restart storage
1517 allocate data structure here, will be deallocated by ReadRestart
1518 include per-particle custom attributes if defined
1519 ------------------------------------------------------------------------- */
1520
unpack_restart(char * buf)1521 int Particle::unpack_restart(char *buf)
1522 {
1523 int nbytes_particle = sizeof(OnePartRestart);
1524 int nbytes_custom = sizeof_custom();
1525 int nbytes = nbytes_particle + nbytes_custom;
1526
1527 char *ptr = buf;
1528 int *ibuf = (int *) buf;
1529 nlocal_restart = *ibuf;
1530 ptr += sizeof(int);
1531 ptr = ROUNDUP(ptr);
1532
1533 particle_restart = (char *)
1534 memory->smalloc(nlocal_restart*nbytes,"particle:particle_restart");
1535
1536 memcpy(particle_restart,ptr,nlocal_restart*nbytes);
1537 ptr += nlocal_restart * sizeof(OnePartRestart);
1538 ptr = ROUNDUP(ptr);
1539
1540 return ptr - buf;
1541 }
1542
1543 /* ----------------------------------------------------------------------
1544 unpack particle info into restart storage
1545 use multiple passes to reduce memory use
1546 allocate data structure here, will be deallocated by ReadRestart
1547 include per-particle custom attributes if defined
1548 NOTE: does not ROUNDUP(ptr) at the end, this is done by caller
1549 ------------------------------------------------------------------------- */
1550
unpack_restart(char * buf,int & nlocal_restart,int step,int pass)1551 void Particle::unpack_restart(char *buf, int &nlocal_restart, int step, int pass)
1552 {
1553 int nbytes_particle = sizeof(OnePartRestart);
1554 int nbytes_custom = sizeof_custom();
1555 int nbytes = nbytes_particle + nbytes_custom;
1556
1557 char *ptr = buf;
1558 if (pass == 0) {
1559 int *ibuf = (int *) buf;
1560 nlocal_restart = *ibuf;
1561 ptr += sizeof(int);
1562 ptr = ROUNDUP(ptr);
1563 }
1564 int start = step*pass;
1565 int end = start+step;
1566 end = MIN(nlocal_restart,end);
1567 step = end - start;
1568
1569 particle_restart = (char *)
1570 memory->smalloc(step*nbytes,"particle:particle_restart");
1571
1572 memcpy(particle_restart,ptr,step*nbytes);
1573
1574 this->nlocal_restart = step;
1575 }
1576
1577 // ----------------------------------------------------------------------
1578 // methods for per-particle custom attributes
1579 // ----------------------------------------------------------------------
1580
1581 /* ----------------------------------------------------------------------
1582 find custom per-atom vector/array with name
1583 return index if found
1584 return -1 if not found
1585 ------------------------------------------------------------------------- */
1586
find_custom(char * name)1587 int Particle::find_custom(char *name)
1588 {
1589 for (int i = 0; i < ncustom; i++)
1590 if (ename[i] && strcmp(ename[i],name) == 0) return i;
1591 return -1;
1592 }
1593
1594 /* ----------------------------------------------------------------------
1595 error checks on existence of custom vectors/arrays
1596 ------------------------------------------------------------------------- */
1597
error_custom()1598 void Particle::error_custom()
1599 {
1600 if (collide && collide->vibstyle == DISCRETE && maxvibmode > 1) {
1601 int index = find_custom((char *) "vibmode");
1602 if (index < 0)
1603 error->all(FLERR,"No custom particle vibmode array defined");
1604 if (esize[index] != maxvibmode)
1605 error->all(FLERR,"Custom particle vibmode array is wrong size");
1606 }
1607 }
1608
1609 /* ----------------------------------------------------------------------
1610 add a custom attribute with name
1611 assumes name does not already exist, except in case of restart
1612 type = 0/1 for int/double
1613 size = 0 for vector, size > 0 for array with size columns
1614 allocate the vector or array to current maxlocal via grow_custom()
1615 return index of its location;
1616 ------------------------------------------------------------------------- */
1617
add_custom(char * name,int type,int size)1618 int Particle::add_custom(char *name, int type, int size)
1619 {
1620 int index;
1621
1622 // if name already exists
1623 // just return index if a restart script and re-defining the name
1624 // else error
1625
1626 index = find_custom(name);
1627 if (index >= 0) {
1628 if (custom_restart_flag == NULL || custom_restart_flag[index] == 1)
1629 error->all(FLERR,"Custom particle attribute name already exists");
1630 custom_restart_flag[index] = 1;
1631 return index;
1632 }
1633
1634 // use first available NULL entry or allocate a new one
1635
1636 for (index = 0; index < ncustom; index++)
1637 if (ename[index] == NULL) break;
1638
1639 if (index == ncustom) {
1640 ncustom++;
1641 ename = (char **) memory->srealloc(ename,ncustom*sizeof(char *),
1642 "particle:ename");
1643 memory->grow(etype,ncustom,"particle:etype");
1644 memory->grow(esize,ncustom,"particle:esize");
1645 memory->grow(ewhich,ncustom,"particle:ewhich");
1646 }
1647
1648 int n = strlen(name) + 1;
1649 ename[index] = new char[n];
1650 strcpy(ename[index],name);
1651 etype[index] = type;
1652 esize[index] = size;
1653
1654 if (type == INT) {
1655 if (size == 0) {
1656 ewhich[index] = ncustom_ivec++;
1657 eivec = (int **)
1658 memory->srealloc(eivec,ncustom_ivec*sizeof(int *),"particle:eivec");
1659 eivec[ncustom_ivec-1] = NULL;
1660 memory->grow(icustom_ivec,ncustom_ivec,"particle:icustom_ivec");
1661 icustom_ivec[ncustom_ivec-1] = index;
1662 } else {
1663 ewhich[index] = ncustom_iarray++;
1664 eiarray = (int ***)
1665 memory->srealloc(eiarray,ncustom_iarray*sizeof(int **),
1666 "particle:eiarray");
1667 eiarray[ncustom_iarray-1] = NULL;
1668 memory->grow(icustom_iarray,ncustom_iarray,"particle:icustom_iarray");
1669 icustom_iarray[ncustom_iarray-1] = index;
1670 memory->grow(eicol,ncustom_iarray,"particle:eicol");
1671 eicol[ncustom_iarray-1] = size;
1672 }
1673 } else if (type == DOUBLE) {
1674 if (size == 0) {
1675 ewhich[index] = ncustom_dvec++;
1676 edvec = (double **)
1677 memory->srealloc(edvec,ncustom_dvec*sizeof(double *),"particle:edvec");
1678 edvec[ncustom_dvec-1] = NULL;
1679 memory->grow(icustom_dvec,ncustom_dvec,"particle:icustom_dvec");
1680 icustom_dvec[ncustom_dvec-1] = index;
1681 } else {
1682 ewhich[index] = ncustom_darray++;
1683 edarray = (double ***)
1684 memory->srealloc(edarray,ncustom_darray*sizeof(double **),
1685 "particle:edarray");
1686 edarray[ncustom_darray-1] = NULL;
1687 memory->grow(icustom_darray,ncustom_darray,"particle:icustom_darray");
1688 icustom_darray[ncustom_darray-1] = index;
1689 memory->grow(edcol,ncustom_darray,"particle:edcol");
1690 edcol[ncustom_darray-1] = size;
1691 }
1692 }
1693
1694 grow_custom(index,0,maxlocal);
1695
1696 return index;
1697 }
1698
1699 /* ----------------------------------------------------------------------
1700 grow the vector/array associated with custom attribute with index
1701 nold = old length, nnew = new length (typically maxlocal)
1702 set new values to 0 via memset()
1703 ------------------------------------------------------------------------- */
1704
grow_custom(int index,int nold,int nnew)1705 void Particle::grow_custom(int index, int nold, int nnew)
1706 {
1707 if (etype[index] == INT) {
1708 if (esize[index] == 0) {
1709 int *ivector = eivec[ewhich[index]];
1710 memory->grow(ivector,nnew,"particle:eivec");
1711 if (ivector) memset(&ivector[nold],0,(nnew-nold)*sizeof(int));
1712 eivec[ewhich[index]] = ivector;
1713 } else {
1714 int **iarray = eiarray[ewhich[index]];
1715 memory->grow(iarray,nnew,esize[index],"particle:eiarray");
1716 if (iarray)
1717 memset(&iarray[nold][0],0,(nnew-nold)*esize[index]*sizeof(int));
1718 eiarray[ewhich[index]] = iarray;
1719 }
1720
1721 } else {
1722 if (esize[index] == 0) {
1723 double *dvector = edvec[ewhich[index]];
1724 memory->grow(dvector,nnew,"particle:edvec");
1725 if (dvector) memset(&dvector[nold],0,(nnew-nold)*sizeof(double));
1726 edvec[ewhich[index]] = dvector;
1727 } else {
1728 double **darray = edarray[ewhich[index]];
1729 memory->grow(darray,nnew,esize[index],"particle:edarray");
1730 if (darray)
1731 memset(&darray[nold][0],0,(nnew-nold)*esize[index]*sizeof(double));
1732 edarray[ewhich[index]] = darray;
1733 }
1734 }
1735 }
1736
1737 /* ----------------------------------------------------------------------
1738 remove a custom attribute at location index
1739 free memory for name and vector/array and set ptrs to NULL
1740 ncustom lists never shrink, but indices stored between
1741 the ncustom list and the dense vector/array lists must be reset
1742 ------------------------------------------------------------------------- */
1743
remove_custom(int index)1744 void Particle::remove_custom(int index)
1745 {
1746 delete [] ename[index];
1747 ename[index] = NULL;
1748
1749 if (etype[index] == INT) {
1750 if (esize[index] == 0) {
1751 memory->destroy(eivec[ewhich[index]]);
1752 ncustom_ivec--;
1753 for (int i = ewhich[index]; i < ncustom_ivec; i++) {
1754 icustom_ivec[i] = icustom_ivec[i+1];
1755 ewhich[icustom_ivec[i]] = i;
1756 eivec[i] = eivec[i+1];
1757 }
1758 } else{
1759 memory->destroy(eiarray[ewhich[index]]);
1760 ncustom_iarray--;
1761 for (int i = ewhich[index]; i < ncustom_iarray; i++) {
1762 icustom_iarray[i] = icustom_iarray[i+1];
1763 ewhich[icustom_iarray[i]] = i;
1764 eiarray[i] = eiarray[i+1];
1765 eicol[i] = eicol[i+1];
1766 }
1767 }
1768 } else if (etype[index] == DOUBLE) {
1769 if (esize[index] == 0) {
1770 memory->destroy(edvec[ewhich[index]]);
1771 ncustom_dvec--;
1772 for (int i = ewhich[index]; i < ncustom_dvec; i++) {
1773 icustom_dvec[i] = icustom_dvec[i+1];
1774 ewhich[icustom_dvec[i]] = i;
1775 edvec[i] = edvec[i+1];
1776 }
1777 } else{
1778 memory->destroy(edarray[ewhich[index]]);
1779 ncustom_darray--;
1780 for (int i = ewhich[index]; i < ncustom_darray; i++) {
1781 icustom_darray[i] = icustom_darray[i+1];
1782 ewhich[icustom_darray[i]] = i;
1783 edarray[i] = edarray[i+1];
1784 edcol[i] = edcol[i+1];
1785 }
1786 }
1787 }
1788
1789 // set ncustom = 0 if custom list is now entirely empty
1790
1791 int empty = 1;
1792 for (int i = 0; i < ncustom; i++)
1793 if (ename[i]) empty = 0;
1794 if (empty) ncustom = 0;
1795 }
1796
1797 /* ----------------------------------------------------------------------
1798 copy info for one particle in custom attribute vectors/arrays
1799 into location I from location J
1800 ------------------------------------------------------------------------- */
1801
copy_custom(int i,int j)1802 void Particle::copy_custom(int i, int j)
1803 {
1804 int m;
1805
1806 // caller does not always check this
1807 // shouldn't be a problem, but valgrind can complain if memcpy to self
1808 // oddly memcpy(&particles[i],&particles[j],sizeof(OnePart)) seems OK
1809
1810 if (i == j) return;
1811
1812 // 4 flavors of vectors/arrays
1813
1814 if (ncustom_ivec) {
1815 for (m = 0; m < ncustom_ivec; m++) eivec[m][i] = eivec[m][j];
1816 }
1817 if (ncustom_iarray) {
1818 for (m = 0; m < ncustom_iarray; m++)
1819 memcpy(eiarray[m][i],eiarray[m][j],eicol[m]*sizeof(int));
1820 }
1821 if (ncustom_dvec) {
1822 for (m = 0; m < ncustom_dvec; m++) edvec[m][i] = edvec[m][j];
1823 }
1824 if (ncustom_darray) {
1825 for (m = 0; m < ncustom_darray; m++)
1826 memcpy(edarray[m][i],edarray[m][j],edcol[m]*sizeof(double));
1827 }
1828 }
1829
1830 /* ----------------------------------------------------------------------
1831 return size of all custom attributes in bytes for one particle
1832 used by callers to allocate buffer memory for particles
1833 assume integer attributes can be put at start of buffer
1834 only alignment needed is between integers and doubles
1835 ------------------------------------------------------------------------- */
1836
sizeof_custom()1837 int Particle::sizeof_custom()
1838 {
1839 int n = 0;
1840
1841 n += ncustom_ivec*sizeof(int);
1842 if (ncustom_iarray)
1843 for (int i = 0; i < ncustom_iarray; i++)
1844 n += eicol[i]*sizeof(int);
1845
1846 n = IROUNDUP(n);
1847
1848 n += ncustom_dvec*sizeof(double);
1849 if (ncustom_darray)
1850 for (int i = 0; i < ncustom_darray; i++)
1851 n += edcol[i]*sizeof(double);
1852
1853 return n;
1854 }
1855
1856 /* ----------------------------------------------------------------------
1857 proc 0 writes custom attribute definition info to restart file
1858 ------------------------------------------------------------------------- */
1859
write_restart_custom(FILE * fp)1860 void Particle::write_restart_custom(FILE *fp)
1861 {
1862 int m,index;
1863
1864 // nactive = # of ncustom that have active vectors/arrays
1865
1866 int nactive = 0;
1867 for (int i = 0; i < ncustom; i++)
1868 if (ename[i]) nactive++;
1869
1870 fwrite(&nactive,sizeof(int),1,fp);
1871
1872 // must write custom info in same order
1873 // the per-particle custom values will be written into file
1874 // not necessarily the same as ncustom list, due to deletions & additions
1875
1876 for (m = 0; m < ncustom_ivec; m++) {
1877 index = icustom_ivec[m];
1878 int n = strlen(ename[index]) + 1;
1879 fwrite(&n,sizeof(int),1,fp);
1880 fwrite(ename[index],sizeof(char),n,fp);
1881 fwrite(&etype[index],sizeof(int),1,fp);
1882 fwrite(&esize[index],sizeof(int),1,fp);
1883 }
1884 for (m = 0; m < ncustom_iarray; m++) {
1885 index = icustom_iarray[m];
1886 int n = strlen(ename[index]) + 1;
1887 fwrite(&n,sizeof(int),1,fp);
1888 fwrite(ename[index],sizeof(char),n,fp);
1889 fwrite(&etype[index],sizeof(int),1,fp);
1890 fwrite(&esize[index],sizeof(int),1,fp);
1891 }
1892 for (m = 0; m < ncustom_dvec; m++) {
1893 index = icustom_dvec[m];
1894 int n = strlen(ename[index]) + 1;
1895 fwrite(&n,sizeof(int),1,fp);
1896 fwrite(ename[index],sizeof(char),n,fp);
1897 fwrite(&etype[index],sizeof(int),1,fp);
1898 fwrite(&esize[index],sizeof(int),1,fp);
1899 }
1900 for (m = 0; m < ncustom_darray; m++) {
1901 index = icustom_darray[m];
1902 int n = strlen(ename[index]) + 1;
1903 fwrite(&n,sizeof(int),1,fp);
1904 fwrite(ename[index],sizeof(char),n,fp);
1905 fwrite(&etype[index],sizeof(int),1,fp);
1906 fwrite(&esize[index],sizeof(int),1,fp);
1907 }
1908 }
1909
1910 /* ----------------------------------------------------------------------
1911 proc 0 reads custom attribute definition info from restart file
1912 bcast to other procs and all procs instantiate series of Mixtures
1913 ------------------------------------------------------------------------- */
1914
read_restart_custom(FILE * fp)1915 void Particle::read_restart_custom(FILE *fp)
1916 {
1917 // ncustom is 0 at time restart file is read
1918 // will be incremented as add_custom() for each nactive
1919
1920 int nactive;
1921 if (me == 0) fread(&nactive,sizeof(int),1,fp);
1922 MPI_Bcast(&nactive,1,MPI_INT,0,world);
1923 if (nactive == 0) return;
1924
1925 // order that custom vectors/arrays are in restart file
1926 // matches order the per-particle custom values will be read from file
1927
1928 int n,type,size;
1929 char *name;
1930
1931 for (int i = 0; i < nactive; i++) {
1932 if (me == 0) fread(&n,sizeof(int),1,fp);
1933 MPI_Bcast(&n,1,MPI_INT,0,world);
1934 name = new char[n];
1935 if (me == 0) fread(name,sizeof(char),n,fp);
1936 MPI_Bcast(name,n,MPI_CHAR,0,world);
1937 if (me == 0) fread(&type,sizeof(int),1,fp);
1938 MPI_Bcast(&type,n,MPI_CHAR,0,world);
1939 if (me == 0) fread(&size,sizeof(int),1,fp);
1940 MPI_Bcast(&size,n,MPI_CHAR,0,world);
1941
1942 // create the custom attribute
1943
1944 add_custom(name,type,size);
1945 delete [] name;
1946 }
1947
1948 // set flag for each newly created custom attribute to 0
1949 // will be reset to 1 if restart script redefines attribute with same name
1950
1951 custom_restart_flag = new int[ncustom];
1952 for (int i = 0; i < ncustom; i++) custom_restart_flag[i] = 0;
1953 }
1954
1955 /* ----------------------------------------------------------------------
1956 pack a custom attributes for a single particle N into buf
1957 this is done in order of 4 styles of vectors/arrays, not in ncustom order
1958 ------------------------------------------------------------------------- */
1959
pack_custom(int n,char * buf)1960 void Particle::pack_custom(int n, char *buf)
1961 {
1962 int i;
1963 char *ptr = buf;
1964
1965 if (ncustom_ivec) {
1966 for (i = 0; i < ncustom_ivec; i++) {
1967 memcpy(ptr,&eivec[i][n],sizeof(int));
1968 ptr += sizeof(int);
1969 }
1970 }
1971 if (ncustom_iarray) {
1972 for (i = 0; i < ncustom_iarray; i++) {
1973 memcpy(ptr,eiarray[i][n],eicol[i]*sizeof(int));
1974 ptr += eicol[i]*sizeof(int);
1975 }
1976 }
1977
1978 ptr = ROUNDUP(ptr);
1979
1980 if (ncustom_dvec) {
1981 for (i = 0; i < ncustom_dvec; i++) {
1982 memcpy(ptr,&edvec[i][n],sizeof(double));
1983 ptr += sizeof(double);
1984 }
1985 }
1986 if (ncustom_darray) {
1987 for (i = 0; i < ncustom_darray; i++) {
1988 memcpy(ptr,edarray[i][n],edcol[i]*sizeof(double));
1989 ptr += edcol[i]*sizeof(double);
1990 }
1991 }
1992 }
1993
1994 /* ----------------------------------------------------------------------
1995 unpack custom attributes for a single particle N from buf
1996 this is done in order of 4 styles of vectors/arrays, not in ncustom order
1997 ------------------------------------------------------------------------- */
1998
unpack_custom(char * buf,int n)1999 void Particle::unpack_custom(char *buf, int n)
2000 {
2001 int i;
2002 char *ptr = buf;
2003
2004 if (ncustom_ivec) {
2005 for (i = 0; i < ncustom_ivec; i++) {
2006 memcpy(&eivec[i][n],ptr,sizeof(int));
2007 ptr += sizeof(int);
2008 }
2009 }
2010 if (ncustom_iarray) {
2011 for (i = 0; i < ncustom_iarray; i++) {
2012 memcpy(eiarray[i][n],ptr,eicol[i]*sizeof(int));
2013 ptr += eicol[i]*sizeof(int);
2014 }
2015 }
2016
2017 ptr = ROUNDUP(ptr);
2018
2019 if (ncustom_dvec) {
2020 for (i = 0; i < ncustom_dvec; i++) {
2021 memcpy(&edvec[i][n],ptr,sizeof(double));
2022 ptr += sizeof(double);
2023 }
2024 }
2025 if (ncustom_darray) {
2026 for (i = 0; i < ncustom_darray; i++) {
2027 memcpy(edarray[i][n],ptr,edcol[i]*sizeof(double));
2028 ptr += edcol[i]*sizeof(double);
2029 }
2030 }
2031 }
2032
2033 /* ---------------------------------------------------------------------- */
2034
memory_usage()2035 bigint Particle::memory_usage()
2036 {
2037 bigint bytes = (bigint) maxlocal * sizeof(OnePart);
2038 bytes += (bigint) maxlocal * sizeof(int);
2039 for (int i = 0; i < ncustom_ivec; i++)
2040 bytes += (bigint) maxlocal * sizeof(int);
2041 for (int i = 0; i < ncustom_iarray; i++)
2042 bytes += (bigint) maxlocal*eicol[i] * sizeof(int);
2043 for (int i = 0; i < ncustom_dvec; i++)
2044 bytes += (bigint) maxlocal * sizeof(double);
2045 for (int i = 0; i < ncustom_darray; i++)
2046 bytes += (bigint) maxlocal*edcol[i] * sizeof(double);
2047 return bytes;
2048 }
2049