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 "spatype.h"
16 #include "mpi.h"
17 #include "math.h"
18 #include "stdlib.h"
19 #include "string.h"
20 #include "update.h"
21 #include "math_const.h"
22 #include "particle.h"
23 #include "modify.h"
24 #include "fix.h"
25 #include "compute.h"
26 #include "domain.h"
27 #include "comm.h"
28 #include "collide.h"
29 #include "grid.h"
30 #include "surf.h"
31 #include "surf_collide.h"
32 #include "surf_react.h"
33 #include "input.h"
34 #include "output.h"
35 #include "geometry.h"
36 #include "random_mars.h"
37 #include "timer.h"
38 #include "math_extra.h"
39 #include "memory.h"
40 #include "error.h"
41
42 using namespace SPARTA_NS;
43
44 enum{XLO,XHI,YLO,YHI,ZLO,ZHI,INTERIOR}; // same as Domain
45 enum{PERIODIC,OUTFLOW,REFLECT,SURFACE,AXISYM}; // same as Domain
46 enum{OUTSIDE,INSIDE,ONSURF2OUT,ONSURF2IN}; // several files
47 enum{PKEEP,PINSERT,PDONE,PDISCARD,PENTRY,PEXIT,PSURF}; // several files
48 enum{NCHILD,NPARENT,NUNKNOWN,NPBCHILD,NPBPARENT,NPBUNKNOWN,NBOUND}; // Grid
49 enum{TALLYAUTO,TALLYREDUCE,TALLYRVOUS}; // same as Surf
50 enum{PERAUTO,PERCELL,PERSURF}; // several files
51
52 #define MAXSTUCK 20
53 #define EPSPARAM 1.0e-7
54
55 // either set ID or PROC/INDEX, set other to -1
56
57 //#define MOVE_DEBUG 1 // un-comment to debug one particle
58 #define MOVE_DEBUG_ID 308143534 // particle ID
59 #define MOVE_DEBUG_PROC -1 // owning proc
60 #define MOVE_DEBUG_INDEX -1 // particle index on owning proc
61 #define MOVE_DEBUG_STEP 4107 // timestep
62
63 /* ---------------------------------------------------------------------- */
64
Update(SPARTA * sparta)65 Update::Update(SPARTA *sparta) : Pointers(sparta)
66 {
67 MPI_Comm_rank(world,&me);
68 MPI_Comm_size(world,&nprocs);
69
70 ntimestep = 0;
71 firststep = laststep = 0;
72 beginstep = endstep = 0;
73 runflag = 0;
74
75 unit_style = NULL;
76 set_units("si");
77
78 fnum = 1.0;
79 nrho = 1.0;
80 vstream[0] = vstream[1] = vstream[2] = 0.0;
81 temp_thermal = 273.15;
82 gravity[0] = gravity[1] = gravity[2] = 0.0;
83
84 maxmigrate = 0;
85 mlist = NULL;
86
87 nslist_compute = nblist_compute = 0;
88 slist_compute = blist_compute = NULL;
89 slist_active = blist_active = NULL;
90
91 nulist_surfcollide = 0;
92 ulist_surfcollide = NULL;
93
94 ranmaster = new RanMars(sparta);
95
96 reorder_period = 0;
97 global_mem_limit = 0;
98 mem_limit_grid_flag = 0;
99
100 copymode = 0;
101 }
102
103 /* ---------------------------------------------------------------------- */
104
~Update()105 Update::~Update()
106 {
107 if (copymode) return;
108
109 delete [] unit_style;
110 memory->destroy(mlist);
111 delete [] slist_compute;
112 delete [] blist_compute;
113 delete [] slist_active;
114 delete [] blist_active;
115 delete [] ulist_surfcollide;
116 delete ranmaster;
117 }
118
119 /* ---------------------------------------------------------------------- */
120
set_units(const char * style)121 void Update::set_units(const char *style)
122 {
123 // physical constants from:
124 // http://physics.nist.gov/cuu/Constants/Table/allascii.txt
125
126 if (strcmp(style,"cgs") == 0) {
127 boltz = 1.3806488e-16;
128 mvv2e = 1.0;
129 dt = 1.0;
130
131 } else if (strcmp(style,"si") == 0) {
132 boltz = 1.3806488e-23;
133 mvv2e = 1.0;
134 dt = 1.0;
135
136 } else error->all(FLERR,"Illegal units command");
137
138 delete [] unit_style;
139 int n = strlen(style) + 1;
140 unit_style = new char[n];
141 strcpy(unit_style,style);
142 }
143
144 /* ---------------------------------------------------------------------- */
145
init()146 void Update::init()
147 {
148 // init the Update class if performing a run, else just return
149 // only set first_update if a run is being performed
150
151 if (runflag == 0) return;
152 first_update = 1;
153
154 // choose the appropriate move method
155
156 if (domain->dimension == 3) {
157 if (surf->exist) moveptr = &Update::move<3,1>;
158 else moveptr = &Update::move<3,0>;
159 } else if (domain->axisymmetric) {
160 if (surf->exist) moveptr = &Update::move<1,1>;
161 else moveptr = &Update::move<1,0>;
162 } else if (domain->dimension == 2) {
163 if (surf->exist) moveptr = &Update::move<2,1>;
164 else moveptr = &Update::move<2,0>;
165 }
166
167 // check gravity vector
168
169 if (domain->dimension == 2 && gravity[2] != 0.0)
170 error->all(FLERR,"Gravity in z not allowed for 2d");
171 if (domain->axisymmetric && gravity[1] != 0.0)
172 error->all(FLERR,"Gravity in y not allowed for axi-symmetric model");
173
174 // moveperturb method is set if particle motion is perturbed
175
176 moveperturb = NULL;
177 if (gravity[0] != 0.0 || gravity[1] != 0.0 || gravity[2] != 0.0) {
178 if (domain->dimension == 3) moveperturb = &Update::gravity3d;
179 if (domain->dimension == 2) moveperturb = &Update::gravity2d;
180 }
181
182 if (moveperturb) perturbflag = 1;
183 else perturbflag = 0;
184 }
185
186 /* ---------------------------------------------------------------------- */
187
setup()188 void Update::setup()
189 {
190 // initialize counters in case stats outputs them
191 // initialize running stats before each run
192
193 ntouch_one = ncomm_one = 0;
194 nboundary_one = nexit_one = 0;
195 nscheck_one = nscollide_one = 0;
196 surf->nreact_one = 0;
197
198 first_running_step = update->ntimestep;
199 niterate_running = 0;
200 nmove_running = ntouch_running = ncomm_running = 0;
201 nboundary_running = nexit_running = 0;
202 nscheck_running = nscollide_running = 0;
203 surf->nreact_running = 0;
204 nstuck = naxibad = 0;
205
206 collide_react = collide_react_setup();
207 bounce_tally = bounce_setup();
208
209 dynamic = 0;
210 dynamic_setup();
211
212 modify->setup();
213 output->setup(1);
214 }
215
216 /* ---------------------------------------------------------------------- */
217
run(int nsteps)218 void Update::run(int nsteps)
219 {
220 int n_start_of_step = modify->n_start_of_step;
221 int n_end_of_step = modify->n_end_of_step;
222
223 // cellweightflag = 1 if grid-based particle weighting is ON
224
225 int cellweightflag = 0;
226 if (grid->cellweightflag) cellweightflag = 1;
227
228 // loop over timesteps
229
230 for (int i = 0; i < nsteps; i++) {
231
232 ntimestep++;
233
234 if (collide_react) collide_react_reset();
235 if (bounce_tally) bounce_set(ntimestep);
236
237 timer->stamp();
238
239 // dynamic parameter updates
240
241 if (dynamic) dynamic_update();
242
243 // start of step fixes
244
245 if (n_start_of_step) {
246 modify->start_of_step();
247 timer->stamp(TIME_MODIFY);
248 }
249
250 // move particles
251
252 if (cellweightflag) particle->pre_weight();
253 (this->*moveptr)();
254 timer->stamp(TIME_MOVE);
255
256 // communicate particles
257
258 comm->migrate_particles(nmigrate,mlist);
259 if (cellweightflag) particle->post_weight();
260 timer->stamp(TIME_COMM);
261
262 if (collide) {
263 particle->sort();
264 timer->stamp(TIME_SORT);
265
266 collide->collisions();
267 timer->stamp(TIME_COLLIDE);
268 }
269
270 if (collide_react) collide_react_update();
271
272 // diagnostic fixes
273
274 if (n_end_of_step) {
275 modify->end_of_step();
276 timer->stamp(TIME_MODIFY);
277 }
278
279 // all output
280
281 if (ntimestep == output->next) {
282 output->write(ntimestep);
283 timer->stamp(TIME_OUTPUT);
284 }
285 }
286 }
287
288 /* ----------------------------------------------------------------------
289 advect particles thru grid
290 DIM = 2/3 for 2d/3d, 1 for 2d axisymmetric
291 SURF = 0/1 for no surfs or surfs
292 use multiple iterations of move/comm if necessary
293 ------------------------------------------------------------------------- */
294
move()295 template < int DIM, int SURF > void Update::move()
296 {
297 bool hitflag;
298 int m,icell,icell_original,nmask,outface,bflag,nflag,pflag,itmp;
299 int side,minside,minsurf,nsurf,cflag,isurf,exclude,stuck_iterate;
300 int pstart,pstop,entryexit,any_entryexit,reaction;
301 surfint *csurfs;
302 cellint *neigh;
303 double dtremain,frac,newfrac,param,minparam,rnew,dtsurf,tc,tmp;
304 double xnew[3],xhold[3],xc[3],vc[3],minxc[3],minvc[3];
305 double *x,*v,*lo,*hi;
306 Grid::ParentCell *pcell;
307 Surf::Tri *tri;
308 Surf::Line *line;
309 Particle::OnePart iorig;
310 Particle::OnePart *particles;
311 Particle::OnePart *ipart,*jpart;
312
313 // for 2d and axisymmetry only
314 // xnew,xc passed to geometry routines which use or set z component
315
316 if (DIM < 3) xnew[2] = xc[2] = 0.0;
317
318 // extend migration list if necessary
319
320 int nlocal = particle->nlocal;
321 int maxlocal = particle->maxlocal;
322
323 if (nlocal > maxmigrate) {
324 maxmigrate = maxlocal;
325 memory->destroy(mlist);
326 memory->create(mlist,maxmigrate,"particle:mlist");
327 }
328
329 // counters
330
331 niterate = 0;
332 ntouch_one = ncomm_one = 0;
333 nboundary_one = nexit_one = 0;
334 nscheck_one = nscollide_one = 0;
335 surf->nreact_one = 0;
336
337 // move/migrate iterations
338
339 Grid::ChildCell *cells = grid->cells;
340 Grid::ParentCell *pcells = grid->pcells;
341 Surf::Tri *tris = surf->tris;
342 Surf::Line *lines = surf->lines;
343 double dt = update->dt;
344
345 // DEBUG
346
347 while (1) {
348
349 // loop over particles
350 // first iteration = all my particles
351 // subsequent iterations = received particles
352
353 niterate++;
354 particles = particle->particles;
355 nmigrate = 0;
356 entryexit = 0;
357
358 if (niterate == 1) {
359 pstart = 0;
360 pstop = nlocal;
361 }
362
363 for (int i = pstart; i < pstop; i++) {
364 pflag = particles[i].flag;
365
366 // received from another proc and move is done
367 // if first iteration, PDONE is from a previous step,
368 // set pflag to PKEEP so move the particle on this step
369 // else do nothing
370
371 if (pflag == PDONE) {
372 pflag = particles[i].flag = PKEEP;
373 if (niterate > 1) continue;
374 }
375
376 x = particles[i].x;
377 v = particles[i].v;
378 exclude = -1;
379
380 // apply moveperturb() to PKEEP and PINSERT since are computing xnew
381 // not to PENTRY,PEXIT since are just re-computing xnew of sender
382 // set xnew[2] to linear move for axisymmetry, will be remapped later
383 // let pflag = PEXIT persist to check during axisymmetric cell crossing
384
385 if (pflag == PKEEP) {
386 dtremain = dt;
387 xnew[0] = x[0] + dtremain*v[0];
388 xnew[1] = x[1] + dtremain*v[1];
389 if (DIM != 2) xnew[2] = x[2] + dtremain*v[2];
390 if (perturbflag) (this->*moveperturb)(dtremain,xnew,v);
391 } else if (pflag == PINSERT) {
392 dtremain = particles[i].dtremain;
393 xnew[0] = x[0] + dtremain*v[0];
394 xnew[1] = x[1] + dtremain*v[1];
395 if (DIM != 2) xnew[2] = x[2] + dtremain*v[2];
396 if (perturbflag) (this->*moveperturb)(dtremain,xnew,v);
397 } else if (pflag == PENTRY) {
398 icell = particles[i].icell;
399 if (cells[icell].nsplit > 1) {
400 if (DIM == 3 && SURF) icell = split3d(icell,x);
401 if (DIM < 3 && SURF) icell = split2d(icell,x);
402 particles[i].icell = icell;
403 }
404 dtremain = particles[i].dtremain;
405 xnew[0] = x[0] + dtremain*v[0];
406 xnew[1] = x[1] + dtremain*v[1];
407 if (DIM != 2) xnew[2] = x[2] + dtremain*v[2];
408 } else if (pflag == PEXIT) {
409 dtremain = particles[i].dtremain;
410 xnew[0] = x[0] + dtremain*v[0];
411 xnew[1] = x[1] + dtremain*v[1];
412 if (DIM != 2) xnew[2] = x[2] + dtremain*v[2];
413 } else if (pflag >= PSURF) {
414 dtremain = particles[i].dtremain;
415 xnew[0] = x[0] + dtremain*v[0];
416 xnew[1] = x[1] + dtremain*v[1];
417 if (DIM != 2) xnew[2] = x[2] + dtremain*v[2];
418 if (pflag > PSURF) exclude = pflag - PSURF - 1;
419 }
420
421 particles[i].flag = PKEEP;
422 icell = particles[i].icell;
423 lo = cells[icell].lo;
424 hi = cells[icell].hi;
425 neigh = cells[icell].neigh;
426 nmask = cells[icell].nmask;
427 stuck_iterate = 0;
428 ntouch_one++;
429
430 // advect one particle from cell to cell and thru surf collides til done
431
432 //int iterate = 0;
433
434 while (1) {
435
436 #ifdef MOVE_DEBUG
437 if (DIM == 3) {
438 if (ntimestep == MOVE_DEBUG_STEP &&
439 (MOVE_DEBUG_ID == particles[i].id ||
440 (me == MOVE_DEBUG_PROC && i == MOVE_DEBUG_INDEX)))
441 printf("PARTICLE %d %ld: %d %d: %d: x %g %g %g: xnew %g %g %g: %d "
442 CELLINT_FORMAT ": lo %g %g %g: hi %g %g %g: DTR %g\n",
443 me,update->ntimestep,i,particles[i].id,
444 cells[icell].nsurf,
445 x[0],x[1],x[2],xnew[0],xnew[1],xnew[2],
446 icell,cells[icell].id,
447 lo[0],lo[1],lo[2],hi[0],hi[1],hi[2],dtremain);
448 }
449 if (DIM == 2) {
450 if (ntimestep == MOVE_DEBUG_STEP &&
451 (MOVE_DEBUG_ID == particles[i].id ||
452 (me == MOVE_DEBUG_PROC && i == MOVE_DEBUG_INDEX)))
453 printf("PARTICLE %d %ld: %d %d: %d: x %g %g: xnew %g %g: %d "
454 CELLINT_FORMAT ": lo %g %g: hi %g %g: DTR: %g\n",
455 me,update->ntimestep,i,particles[i].id,
456 cells[icell].nsurf,
457 x[0],x[1],xnew[0],xnew[1],
458 icell,cells[icell].id,
459 lo[0],lo[1],hi[0],hi[1],dtremain);
460 }
461 if (DIM == 1) {
462 if (ntimestep == MOVE_DEBUG_STEP &&
463 (MOVE_DEBUG_ID == particles[i].id ||
464 (me == MOVE_DEBUG_PROC && i == MOVE_DEBUG_INDEX)))
465 printf("PARTICLE %d %ld: %d %d: %d: x %g %g: xnew %g %g: %d "
466 CELLINT_FORMAT ": lo %g %g: hi %g %g: DTR: %g\n",
467 me,update->ntimestep,i,particles[i].id,
468 cells[icell].nsurf,
469 x[0],x[1],xnew[0],sqrt(xnew[1]*xnew[1]+xnew[2]*xnew[2]),
470 icell,cells[icell].id,
471 lo[0],lo[1],hi[0],hi[1],dtremain);
472 }
473 #endif
474
475 // check if particle crosses any cell face
476 // frac = fraction of move completed before hitting cell face
477 // this section should be as efficient as possible,
478 // since most particles won't do anything else
479 // axisymmetric y cell face crossings:
480 // these faces are curved cylindrical shells
481 // axi_horizontal_line() checks for intersection of
482 // straight-line y,z move with circle in y,z
483 // always check move against lower y face
484 // except when particle starts on face and
485 // PEXIT is set (just received) or particle is moving downward in y
486 // only check move against upper y face
487 // if remapped final y position (rnew) is within cell,
488 // or except when particle starts on face and
489 // PEXIT is set (just received) or particle is moving upward in y
490 // unset pflag so not checked again for this particle
491
492 outface = INTERIOR;
493 frac = 1.0;
494
495 if (xnew[0] < lo[0]) {
496 frac = (lo[0]-x[0]) / (xnew[0]-x[0]);
497 outface = XLO;
498 } else if (xnew[0] >= hi[0]) {
499 frac = (hi[0]-x[0]) / (xnew[0]-x[0]);
500 outface = XHI;
501 }
502
503 if (DIM != 1) {
504 if (xnew[1] < lo[1]) {
505 newfrac = (lo[1]-x[1]) / (xnew[1]-x[1]);
506 if (newfrac < frac) {
507 frac = newfrac;
508 outface = YLO;
509 }
510 } else if (xnew[1] >= hi[1]) {
511 newfrac = (hi[1]-x[1]) / (xnew[1]-x[1]);
512 if (newfrac < frac) {
513 frac = newfrac;
514 outface = YHI;
515 }
516 }
517 }
518
519 if (DIM == 1) {
520 if (x[1] == lo[1] && (pflag == PEXIT || v[1] < 0.0)) {
521 frac = 0.0;
522 outface = YLO;
523 } else if (Geometry::
524 axi_horizontal_line(dtremain,x,v,lo[1],itmp,tc,tmp)) {
525 newfrac = tc/dtremain;
526 if (newfrac < frac) {
527 frac = newfrac;
528 outface = YLO;
529 }
530 }
531
532 if (x[1] == hi[1] && (pflag == PEXIT || v[1] > 0.0)) {
533 frac = 0.0;
534 outface = YHI;
535 } else {
536 rnew = sqrt(xnew[1]*xnew[1] + xnew[2]*xnew[2]);
537 if (rnew >= hi[1]) {
538 if (Geometry::
539 axi_horizontal_line(dtremain,x,v,hi[1],itmp,tc,tmp)) {
540 newfrac = tc/dtremain;
541 if (newfrac < frac) {
542 frac = newfrac;
543 outface = YHI;
544 }
545 }
546 }
547 }
548
549 pflag = 0;
550 }
551
552 if (DIM == 3) {
553 if (xnew[2] < lo[2]) {
554 newfrac = (lo[2]-x[2]) / (xnew[2]-x[2]);
555 if (newfrac < frac) {
556 frac = newfrac;
557 outface = ZLO;
558 }
559 } else if (xnew[2] >= hi[2]) {
560 newfrac = (hi[2]-x[2]) / (xnew[2]-x[2]);
561 if (newfrac < frac) {
562 frac = newfrac;
563 outface = ZHI;
564 }
565 }
566 }
567
568 //if (iterate == 10) exit(1);
569 //iterate++;
570
571 #ifdef MOVE_DEBUG
572 if (ntimestep == MOVE_DEBUG_STEP &&
573 (MOVE_DEBUG_ID == particles[i].id ||
574 (me == MOVE_DEBUG_PROC && i == MOVE_DEBUG_INDEX))) {
575 if (outface != INTERIOR)
576 printf(" OUTFACE %d out: %d %d, frac %g\n",
577 outface,grid->neigh_decode(nmask,outface),
578 neigh[outface],frac);
579 else
580 printf(" INTERIOR %d %d\n",outface,INTERIOR);
581 }
582 #endif
583
584 // START of code specific to surfaces
585
586 if (SURF) {
587
588 // skip surf checks if particle flagged as EXITing this cell
589 // then unset pflag so not checked again for this particle
590
591 nsurf = cells[icell].nsurf;
592 if (pflag == PEXIT) {
593 nsurf = 0;
594 pflag = 0;
595 }
596 nscheck_one += nsurf;
597
598 if (nsurf) {
599
600 // particle crosses cell face, reset xnew exactly on face of cell
601 // so surface check occurs only for particle path within grid cell
602 // xhold = saved xnew so can restore below if no surf collision
603
604 if (outface != INTERIOR) {
605 xhold[0] = xnew[0];
606 xhold[1] = xnew[1];
607 if (DIM != 2) xhold[2] = xnew[2];
608
609 xnew[0] = x[0] + frac*(xnew[0]-x[0]);
610 xnew[1] = x[1] + frac*(xnew[1]-x[1]);
611 if (DIM != 2) xnew[2] = x[2] + frac*(xnew[2]-x[2]);
612
613 if (outface == XLO) xnew[0] = lo[0];
614 else if (outface == XHI) xnew[0] = hi[0];
615 else if (outface == YLO) xnew[1] = lo[1];
616 else if (outface == YHI) xnew[1] = hi[1];
617 else if (outface == ZLO) xnew[2] = lo[2];
618 else if (outface == ZHI) xnew[2] = hi[2];
619 }
620
621 // for axisymmetric, dtsurf = time that particle stays in cell
622 // used as arg to axi_line_intersect()
623
624 if (DIM == 1) {
625 if (outface == INTERIOR) dtsurf = dtremain;
626 else dtsurf = dtremain * frac;
627 }
628
629 // check for collisions with triangles or lines in cell
630 // find 1st surface hit via minparam
631 // skip collisions with previous surf, but not for axisymmetric
632 // not considered collision if 2 params are tied and one INSIDE surf
633 // if collision occurs, perform collision with surface model
634 // reset x,v,xnew,dtremain and continue single particle trajectory
635
636 cflag = 0;
637 minparam = 2.0;
638 csurfs = cells[icell].csurfs;
639
640 for (m = 0; m < nsurf; m++) {
641 isurf = csurfs[m];
642
643 if (DIM > 1) {
644 if (isurf == exclude) continue;
645 }
646 if (DIM == 3) {
647 tri = &tris[isurf];
648 hitflag = Geometry::
649 line_tri_intersect(x,xnew,tri->p1,tri->p2,tri->p3,
650 tri->norm,xc,param,side);
651 }
652 if (DIM == 2) {
653 line = &lines[isurf];
654 hitflag = Geometry::
655 line_line_intersect(x,xnew,line->p1,line->p2,
656 line->norm,xc,param,side);
657 }
658 if (DIM == 1) {
659 line = &lines[isurf];
660 hitflag = Geometry::
661 axi_line_intersect(dtsurf,x,v,outface,lo,hi,line->p1,line->p2,
662 line->norm,exclude == isurf,
663 xc,vc,param,side);
664 }
665
666 #ifdef MOVE_DEBUG
667 if (DIM == 3) {
668 if (hitflag && ntimestep == MOVE_DEBUG_STEP &&
669 (MOVE_DEBUG_ID == particles[i].id ||
670 (me == MOVE_DEBUG_PROC && i == MOVE_DEBUG_INDEX)))
671 printf("SURF COLLIDE: %d %d %d %d: "
672 "P1 %g %g %g: P2 %g %g %g: "
673 "T1 %g %g %g: T2 %g %g %g: T3 %g %g %g: "
674 "TN %g %g %g: XC %g %g %g: "
675 "Param %g: Side %d\n",
676 MOVE_DEBUG_INDEX,icell,nsurf,isurf,
677 x[0],x[1],x[2],xnew[0],xnew[1],xnew[2],
678 tri->p1[0],tri->p1[1],tri->p1[2],
679 tri->p2[0],tri->p2[1],tri->p2[2],
680 tri->p3[0],tri->p3[1],tri->p3[2],
681 tri->norm[0],tri->norm[1],tri->norm[2],
682 xc[0],xc[1],xc[2],param,side);
683 }
684 if (DIM == 2) {
685 if (hitflag && ntimestep == MOVE_DEBUG_STEP &&
686 (MOVE_DEBUG_ID == particles[i].id ||
687 (me == MOVE_DEBUG_PROC && i == MOVE_DEBUG_INDEX)))
688 printf("SURF COLLIDE: %d %d %d %d: P1 %g %g: P2 %g %g: "
689 "L1 %g %g: L2 %g %g: LN %g %g: XC %g %g: "
690 "Param %g: Side %d\n",
691 MOVE_DEBUG_INDEX,icell,nsurf,isurf,
692 x[0],x[1],xnew[0],xnew[1],
693 line->p1[0],line->p1[1],line->p2[0],line->p2[1],
694 line->norm[0],line->norm[1],
695 xc[0],xc[1],param,side);
696 }
697 if (DIM == 1) {
698 if (hitflag && ntimestep == MOVE_DEBUG_STEP &&
699 (MOVE_DEBUG_ID == particles[i].id ||
700 (me == MOVE_DEBUG_PROC && i == MOVE_DEBUG_INDEX)))
701 printf("SURF COLLIDE %d %ld: %d %d %d %d: P1 %g %g: P2 %g %g: "
702 "L1 %g %g: L2 %g %g: LN %g %g: XC %g %g: "
703 "VC %g %g %g: Param %g: Side %d\n",
704 hitflag,ntimestep,MOVE_DEBUG_INDEX,icell,nsurf,isurf,
705 x[0],x[1],
706 xnew[0],sqrt(xnew[1]*xnew[1]+xnew[2]*xnew[2]),
707 line->p1[0],line->p1[1],line->p2[0],line->p2[1],
708 line->norm[0],line->norm[1],
709 xc[0],xc[1],vc[0],vc[1],vc[2],param,side);
710 double edge1[3],edge2[3],xfinal[3],cross[3];
711 MathExtra::sub3(line->p2,line->p1,edge1);
712 MathExtra::sub3(x,line->p1,edge2);
713 MathExtra::cross3(edge2,edge1,cross);
714 if (hitflag && ntimestep == MOVE_DEBUG_STEP &&
715 MOVE_DEBUG_ID == particles[i].id)
716 printf("CROSSSTART %g %g %g\n",cross[0],cross[1],cross[2]);
717 xfinal[0] = xnew[0];
718 xfinal[1] = sqrt(xnew[1]*xnew[1]+xnew[2]*xnew[2]);
719 xfinal[2] = 0.0;
720 MathExtra::sub3(xfinal,line->p1,edge2);
721 MathExtra::cross3(edge2,edge1,cross);
722 if (hitflag && ntimestep == MOVE_DEBUG_STEP &&
723 MOVE_DEBUG_ID == particles[i].id)
724 printf("CROSSFINAL %g %g %g\n",cross[0],cross[1],cross[2]);
725 }
726 #endif
727
728 if (hitflag && param < minparam && side == OUTSIDE) {
729 cflag = 1;
730 minparam = param;
731 minside = side;
732 minsurf = isurf;
733 minxc[0] = xc[0];
734 minxc[1] = xc[1];
735 if (DIM == 3) minxc[2] = xc[2];
736 if (DIM == 1) {
737 minvc[1] = vc[1];
738 minvc[2] = vc[2];
739 }
740 }
741
742 } // END of for loop over surfs
743
744 // tri/line = surf that particle hit first
745
746 if (cflag) {
747 if (DIM == 3) tri = &tris[minsurf];
748 if (DIM != 3) line = &lines[minsurf];
749
750 // set x to collision point
751 // if axisymmetric, set v to remapped velocity at collision pt
752
753 x[0] = minxc[0];
754 x[1] = minxc[1];
755 if (DIM == 3) x[2] = minxc[2];
756 if (DIM == 1) {
757 v[1] = minvc[1];
758 v[2] = minvc[2];
759 }
760
761 // perform surface collision using surface collision model
762 // surface chemistry may destroy particle or create new one
763 // must update particle's icell to current icell so that
764 // if jpart is created, it will be added to correct cell
765 // if jpart, add new particle to this iteration via pstop++
766 // tally surface statistics if requested using iorig
767
768 ipart = &particles[i];
769 ipart->icell = icell;
770 dtremain *= 1.0 - minparam*frac;
771
772 if (nsurf_tally)
773 memcpy(&iorig,&particles[i],sizeof(Particle::OnePart));
774
775 if (DIM == 3)
776 jpart = surf->sc[tri->isc]->
777 collide(ipart,dtremain,minsurf,tri->norm,tri->isr,reaction);
778 if (DIM != 3)
779 jpart = surf->sc[line->isc]->
780 collide(ipart,dtremain,minsurf,line->norm,line->isr,reaction);
781
782 if (jpart) {
783 particles = particle->particles;
784 x = particles[i].x;
785 v = particles[i].v;
786 jpart->flag = PSURF + 1 + minsurf;
787 jpart->dtremain = dtremain;
788 jpart->weight = particles[i].weight;
789 pstop++;
790 }
791
792 if (nsurf_tally)
793 for (m = 0; m < nsurf_tally; m++)
794 slist_active[m]->surf_tally(minsurf,icell,reaction,
795 &iorig,ipart,jpart);
796
797 // stuck_iterate = consecutive iterations particle is immobile
798
799 if (minparam == 0.0) stuck_iterate++;
800 else stuck_iterate = 0;
801
802 // reset post-bounce xnew
803
804 xnew[0] = x[0] + dtremain*v[0];
805 xnew[1] = x[1] + dtremain*v[1];
806 if (DIM != 2) xnew[2] = x[2] + dtremain*v[2];
807
808 exclude = minsurf;
809 nscollide_one++;
810
811 #ifdef MOVE_DEBUG
812 if (DIM == 3) {
813 if (ntimestep == MOVE_DEBUG_STEP &&
814 (MOVE_DEBUG_ID == particles[i].id ||
815 (me == MOVE_DEBUG_PROC && i == MOVE_DEBUG_INDEX)))
816 printf("POST COLLISION %d: %g %g %g: %g %g %g: %g %g %g\n",
817 MOVE_DEBUG_INDEX,
818 x[0],x[1],x[2],xnew[0],xnew[1],xnew[2],
819 minparam,frac,dtremain);
820 }
821 if (DIM == 2) {
822 if (ntimestep == MOVE_DEBUG_STEP &&
823 (MOVE_DEBUG_ID == particles[i].id ||
824 (me == MOVE_DEBUG_PROC && i == MOVE_DEBUG_INDEX)))
825 printf("POST COLLISION %d: %g %g: %g %g: %g %g %g\n",
826 MOVE_DEBUG_INDEX,
827 x[0],x[1],xnew[0],xnew[1],
828 minparam,frac,dtremain);
829 }
830 if (DIM == 1) {
831 if (ntimestep == MOVE_DEBUG_STEP &&
832 (MOVE_DEBUG_ID == particles[i].id ||
833 (me == MOVE_DEBUG_PROC && i == MOVE_DEBUG_INDEX)))
834 printf("POST COLLISION %d: %g %g: %g %g: vel %g %g %g: %g %g %g\n",
835 MOVE_DEBUG_INDEX,
836 x[0],x[1],
837 xnew[0],sqrt(xnew[1]*xnew[1]+xnew[2]*xnew[2]),
838 v[0],v[1],v[2],
839 minparam,frac,dtremain);
840 }
841 #endif
842
843 // if ipart = NULL, particle discarded due to surface chem
844 // else if particle not stuck, continue advection while loop
845 // if stuck, mark for DISCARD, and drop out of SURF code
846
847 if (ipart == NULL) particles[i].flag = PDISCARD;
848 else if (stuck_iterate < MAXSTUCK) continue;
849 else {
850 particles[i].flag = PDISCARD;
851 nstuck++;
852 }
853
854 } // END of cflag if section that performed collision
855
856 // no collision, so restore saved xnew if changed it above
857
858 if (outface != INTERIOR) {
859 xnew[0] = xhold[0];
860 xnew[1] = xhold[1];
861 if (DIM != 2) xnew[2] = xhold[2];
862 }
863
864 } // END of if test for any surfs in this cell
865 } // END of code specific to surfaces
866
867 // break from advection loop if discarding particle
868
869 if (particles[i].flag == PDISCARD) break;
870
871 // no cell crossing and no surface collision
872 // set final particle position to xnew, then break from advection loop
873 // for axisymmetry, must first remap linear xnew and v
874 // for axisymmetry, check if final particle position is within cell
875 // can be rare epsilon round-off cases where particle ends up outside
876 // of final cell curved surf when move logic thinks it is inside
877 // example is when Geom::axi_horizontal_line() says no crossing of cell edge
878 // but axi_remap() puts particle outside the cell
879 // in this case, just DISCARD particle and tally it to naxibad
880 // if migrating to another proc,
881 // flag as PDONE so new proc won't move it more on this step
882
883 if (outface == INTERIOR) {
884 if (DIM == 1) axi_remap(xnew,v);
885 x[0] = xnew[0];
886 x[1] = xnew[1];
887 if (DIM == 3) x[2] = xnew[2];
888 if (DIM == 1) {
889 if (x[1] < lo[1] || x[1] > hi[1]) {
890 particles[i].flag = PDISCARD;
891 naxibad++;
892 break;
893 }
894 }
895 if (cells[icell].proc != me) particles[i].flag = PDONE;
896 break;
897 }
898
899 // particle crosses cell face
900 // decrement dtremain in case particle is passed to another proc
901 // for axisymmetry, must then remap linear x and v
902 // reset particle x to be exactly on cell face
903 // for axisymmetry, must reset xnew for next iteration since v changed
904
905 dtremain *= 1.0-frac;
906 exclude = -1;
907
908 x[0] += frac * (xnew[0]-x[0]);
909 x[1] += frac * (xnew[1]-x[1]);
910 if (DIM != 2) x[2] += frac * (xnew[2]-x[2]);
911 if (DIM == 1) axi_remap(x,v);
912
913 if (outface == XLO) x[0] = lo[0];
914 else if (outface == XHI) x[0] = hi[0];
915 else if (outface == YLO) x[1] = lo[1];
916 else if (outface == YHI) x[1] = hi[1];
917 else if (outface == ZLO) x[2] = lo[2];
918 else if (outface == ZHI) x[2] = hi[2];
919
920 if (DIM == 1) {
921 xnew[0] = x[0] + dtremain*v[0];
922 xnew[1] = x[1] + dtremain*v[1];
923 xnew[2] = x[2] + dtremain*v[2];
924 }
925
926 // nflag = type of neighbor cell: child, parent, unknown, boundary
927 // if parent, use id_find_child to identify child cell
928 // result can be -1 for unknown cell, occurs when:
929 // (a) particle hits face of ghost child cell
930 // (b) the ghost cell extends beyond ghost halo
931 // (c) cell on other side of face is a parent
932 // (d) its child, which the particle is in, is entirely beyond my halo
933 // if new cell is child and surfs exist, check if a split cell
934
935 nflag = grid->neigh_decode(nmask,outface);
936 icell_original = icell;
937
938 if (nflag == NCHILD) {
939 icell = neigh[outface];
940 if (DIM == 3 && SURF) {
941 if (cells[icell].nsplit > 1 && cells[icell].nsurf >= 0)
942 icell = split3d(icell,x);
943 }
944 if (DIM < 3 && SURF) {
945 if (cells[icell].nsplit > 1 && cells[icell].nsurf >= 0)
946 icell = split2d(icell,x);
947 }
948 } else if (nflag == NPARENT) {
949 pcell = &pcells[neigh[outface]];
950 icell = grid->id_find_child(pcell->id,cells[icell].level,
951 pcell->lo,pcell->hi,x);
952 if (icell >= 0) {
953 if (DIM == 3 && SURF) {
954 if (cells[icell].nsplit > 1 && cells[icell].nsurf >= 0)
955 icell = split3d(icell,x);
956 }
957 if (DIM < 3 && SURF) {
958 if (cells[icell].nsplit > 1 && cells[icell].nsurf >= 0)
959 icell = split2d(icell,x);
960 }
961 }
962 } else if (nflag == NUNKNOWN) icell = -1;
963
964 // neighbor cell is global boundary
965 // tally boundary stats if requested using iorig
966 // collide() updates x,v,xnew as needed due to boundary interaction
967 // may also update dtremain (piston BC)
968 // for axisymmetric, must recalculate xnew since v may have changed
969 // surface chemistry may destroy particle or create new one
970 // if jpart, add new particle to this iteration via pstop++
971 // OUTFLOW: exit with particle flag = PDISCARD
972 // PERIODIC: new cell via same logic as above for child/parent/unknown
973 // OTHER: reflected particle stays in same grid cell
974
975 else {
976 ipart = &particles[i];
977
978 if (nboundary_tally)
979 memcpy(&iorig,&particles[i],sizeof(Particle::OnePart));
980
981 bflag = domain->collide(ipart,outface,icell,xnew,dtremain,
982 jpart,reaction);
983
984 if (jpart) {
985 particles = particle->particles;
986 x = particles[i].x;
987 v = particles[i].v;
988 }
989
990 if (nboundary_tally)
991 for (m = 0; m < nboundary_tally; m++)
992 blist_active[m]->
993 boundary_tally(outface,bflag,reaction,&iorig,ipart,jpart);
994
995 if (DIM == 1) {
996 xnew[0] = x[0] + dtremain*v[0];
997 xnew[1] = x[1] + dtremain*v[1];
998 xnew[2] = x[2] + dtremain*v[2];
999 }
1000
1001 if (bflag == OUTFLOW) {
1002 particles[i].flag = PDISCARD;
1003 nexit_one++;
1004 break;
1005
1006 } else if (bflag == PERIODIC) {
1007 if (nflag == NPBCHILD) {
1008 icell = neigh[outface];
1009 if (DIM == 3 && SURF) {
1010 if (cells[icell].nsplit > 1 && cells[icell].nsurf >= 0)
1011 icell = split3d(icell,x);
1012 }
1013 if (DIM < 3 && SURF) {
1014 if (cells[icell].nsplit > 1 && cells[icell].nsurf >= 0)
1015 icell = split2d(icell,x);
1016 }
1017 } else if (nflag == NPBPARENT) {
1018 pcell = &pcells[neigh[outface]];
1019 icell = grid->id_find_child(pcell->id,cells[icell].level,
1020 pcell->lo,pcell->hi,x);
1021 if (icell >= 0) {
1022 if (DIM == 3 && SURF) {
1023 if (cells[icell].nsplit > 1 && cells[icell].nsurf >= 0)
1024 icell = split3d(icell,x);
1025 }
1026 if (DIM < 3 && SURF) {
1027 if (cells[icell].nsplit > 1 && cells[icell].nsurf >= 0)
1028 icell = split2d(icell,x);
1029 }
1030 } else domain->uncollide(outface,x);
1031 } else if (nflag == NPBUNKNOWN) {
1032 icell = -1;
1033 domain->uncollide(outface,x);
1034 }
1035
1036 } else if (bflag == SURFACE) {
1037 if (ipart == NULL) {
1038 particles[i].flag = PDISCARD;
1039 break;
1040 } else if (jpart) {
1041 jpart->flag = PSURF;
1042 jpart->dtremain = dtremain;
1043 jpart->weight = particles[i].weight;
1044 pstop++;
1045 }
1046 nboundary_one++;
1047 ntouch_one--; // decrement here since will increment below
1048
1049 } else {
1050 nboundary_one++;
1051 ntouch_one--; // decrement here since will increment below
1052 }
1053 }
1054
1055 // neighbor cell is unknown
1056 // reset icell to original icell which must be a ghost cell
1057 // exit with particle flag = PEXIT, so receiver can identify neighbor
1058
1059 if (icell < 0) {
1060 icell = icell_original;
1061 particles[i].flag = PEXIT;
1062 particles[i].dtremain = dtremain;
1063 entryexit = 1;
1064 break;
1065 }
1066
1067 // if nsurf < 0, new cell is EMPTY ghost
1068 // exit with particle flag = PENTRY, so receiver can continue move
1069
1070 if (cells[icell].nsurf < 0) {
1071 particles[i].flag = PENTRY;
1072 particles[i].dtremain = dtremain;
1073 entryexit = 1;
1074 break;
1075 }
1076
1077 // move particle into new grid cell for next stage of move
1078
1079 lo = cells[icell].lo;
1080 hi = cells[icell].hi;
1081 neigh = cells[icell].neigh;
1082 nmask = cells[icell].nmask;
1083 ntouch_one++;
1084 }
1085
1086 // END of while loop over advection of single particle
1087
1088 #ifdef MOVE_DEBUG
1089 if (ntimestep == MOVE_DEBUG_STEP &&
1090 (MOVE_DEBUG_ID == particles[i].id ||
1091 (me == MOVE_DEBUG_PROC && i == MOVE_DEBUG_INDEX)))
1092 printf("MOVE DONE %d %d %d: %g %g %g: DTR %g\n",
1093 MOVE_DEBUG_INDEX,particles[i].flag,icell,
1094 x[0],x[1],x[2],dtremain);
1095 #endif
1096
1097 // move is complete, or as much as can be done on this proc
1098 // update particle's grid cell
1099 // if particle flag set, add particle to migrate list
1100 // if discarding, migration will delete particle
1101
1102 particles[i].icell = icell;
1103
1104 if (particles[i].flag != PKEEP) {
1105 mlist[nmigrate++] = i;
1106 if (particles[i].flag != PDISCARD) {
1107 if (cells[icell].proc == me) {
1108 char str[128];
1109 sprintf(str,
1110 "Particle %d on proc %d being sent to self "
1111 "on step " BIGINT_FORMAT,
1112 i,me,update->ntimestep);
1113 error->one(FLERR,str);
1114 }
1115 ncomm_one++;
1116 }
1117 }
1118 }
1119
1120 // END of pstart/pstop loop advecting all particles
1121
1122 // if gridcut >= 0.0, check if another iteration of move is required
1123 // only the case if some particle flag = PENTRY/PEXIT
1124 // in which case perform particle migration
1125 // if not, move is done and final particle comm will occur in run()
1126 // if iterating, reset pstart/pstop and extend migration list if necessary
1127
1128 if (grid->cutoff < 0.0) break;
1129
1130 timer->stamp(TIME_MOVE);
1131 MPI_Allreduce(&entryexit,&any_entryexit,1,MPI_INT,MPI_MAX,world);
1132 timer->stamp();
1133
1134 if (any_entryexit) {
1135 timer->stamp(TIME_MOVE);
1136 pstart = comm->migrate_particles(nmigrate,mlist);
1137 timer->stamp(TIME_COMM);
1138 pstop = particle->nlocal;
1139 if (pstop-pstart > maxmigrate) {
1140 maxmigrate = pstop-pstart;
1141 memory->destroy(mlist);
1142 memory->create(mlist,maxmigrate,"particle:mlist");
1143 }
1144 } else break;
1145
1146 // END of single move/migrate iteration
1147
1148 }
1149
1150 // END of all move/migrate iterations
1151
1152 particle->sorted = 0;
1153
1154 // accumulate running totals
1155
1156 niterate_running += niterate;
1157 nmove_running += nlocal;
1158 ntouch_running += ntouch_one;
1159 ncomm_running += ncomm_one;
1160 nboundary_running += nboundary_one;
1161 nexit_running += nexit_one;
1162 nscheck_running += nscheck_one;
1163 nscollide_running += nscollide_one;
1164 surf->nreact_running += surf->nreact_one;
1165 }
1166
1167 /* ----------------------------------------------------------------------
1168 particle is entering split parent icell at x
1169 determine which split child cell it is in
1170 return index of sub-cell in ChildCell
1171 ------------------------------------------------------------------------- */
1172
split3d(int icell,double * x)1173 int Update::split3d(int icell, double *x)
1174 {
1175 int m,cflag,isurf,hitflag,side,minsurfindex;
1176 double param,minparam;
1177 double xc[3];
1178 Surf::Tri *tri;
1179
1180 Grid::ChildCell *cells = grid->cells;
1181 Grid::SplitInfo *sinfo = grid->sinfo;
1182 Surf::Tri *tris = surf->tris;
1183
1184 // check for collisions with lines in cell
1185 // find 1st surface hit via minparam
1186 // only consider tris that are mapped via csplits to a split cell
1187 // unmapped tris only touch cell surf at xnew
1188 // another mapped tri should include same xnew
1189 // NOTE: these next 2 lines do not seem correct compared to code
1190 // not considered a collision if particles starts on surf, moving out
1191 // not considered a collision if 2 params are tied and one is INSIDE surf
1192
1193 int nsurf = cells[icell].nsurf;
1194 surfint *csurfs = cells[icell].csurfs;
1195 int isplit = cells[icell].isplit;
1196 int *csplits = sinfo[isplit].csplits;
1197 double *xnew = sinfo[isplit].xsplit;
1198
1199 cflag = 0;
1200 minparam = 2.0;
1201
1202 for (m = 0; m < nsurf; m++) {
1203 if (csplits[m] < 0) continue;
1204 isurf = csurfs[m];
1205 tri = &tris[isurf];
1206 hitflag = Geometry::
1207 line_tri_intersect(x,xnew,tri->p1,tri->p2,tri->p3,
1208 tri->norm,xc,param,side);
1209
1210 if (hitflag && side != INSIDE && param < minparam) {
1211 cflag = 1;
1212 minparam = param;
1213 minsurfindex = m;
1214 }
1215 }
1216
1217 if (!cflag) return sinfo[isplit].csubs[sinfo[isplit].xsub];
1218 int index = csplits[minsurfindex];
1219 return sinfo[isplit].csubs[index];
1220 }
1221
1222 /* ----------------------------------------------------------------------
1223 particle is entering split ICELL at X
1224 determine which split sub-cell it is in
1225 return index of sub-cell in ChildCell
1226 ------------------------------------------------------------------------- */
1227
split2d(int icell,double * x)1228 int Update::split2d(int icell, double *x)
1229 {
1230 int m,cflag,isurf,hitflag,side,minsurfindex;
1231 double param,minparam;
1232 double xc[3];
1233 Surf::Line *line;
1234
1235 Grid::ChildCell *cells = grid->cells;
1236 Grid::SplitInfo *sinfo = grid->sinfo;
1237 Surf::Line *lines = surf->lines;
1238
1239 // check for collisions with lines in cell
1240 // find 1st surface hit via minparam
1241 // only consider lines that are mapped via csplits to a split cell
1242 // unmapped lines only touch cell surf at xnew
1243 // another mapped line should include same xnew
1244 // NOTE: these next 2 lines do not seem correct compared to code
1245 // not considered a collision if particle starts on surf, moving out
1246 // not considered a collision if 2 params are tied and one is INSIDE surf
1247
1248 int nsurf = cells[icell].nsurf;
1249 surfint *csurfs = cells[icell].csurfs;
1250 int isplit = cells[icell].isplit;
1251 int *csplits = sinfo[isplit].csplits;
1252 double *xnew = sinfo[isplit].xsplit;
1253
1254 cflag = 0;
1255 minparam = 2.0;
1256 for (m = 0; m < nsurf; m++) {
1257 if (csplits[m] < 0) continue;
1258 isurf = csurfs[m];
1259 line = &lines[isurf];
1260 hitflag = Geometry::
1261 line_line_intersect(x,xnew,line->p1,line->p2,line->norm,xc,param,side);
1262
1263 if (hitflag && side != INSIDE && param < minparam) {
1264 cflag = 1;
1265 minparam = param;
1266 minsurfindex = m;
1267 }
1268 }
1269
1270 if (!cflag) return sinfo[isplit].csubs[sinfo[isplit].xsub];
1271 int index = csplits[minsurfindex];
1272 return sinfo[isplit].csubs[index];
1273 }
1274
1275 /* ----------------------------------------------------------------------
1276 check if any surface collision or reaction models are defined
1277 return 1 if there are any, 0 if not
1278 ------------------------------------------------------------------------- */
1279
collide_react_setup()1280 int Update::collide_react_setup()
1281 {
1282 nsc = surf->nsc;
1283 sc = surf->sc;
1284 nsr = surf->nsr;
1285 sr = surf->sr;
1286
1287 if (nsc || nsr) return 1;
1288 return 0;
1289 }
1290
1291 /* ----------------------------------------------------------------------
1292 zero counters for tallying surface collisions/reactions
1293 done at start of each timestep
1294 done within individual SurfCollide and SurfReact instances
1295 ------------------------------------------------------------------------- */
1296
collide_react_reset()1297 void Update::collide_react_reset()
1298 {
1299 for (int i = 0; i < nsc; i++) sc[i]->tally_reset();
1300 for (int i = 0; i < nsr; i++) sr[i]->tally_reset();
1301 }
1302
1303 /* ----------------------------------------------------------------------
1304 update cummulative counters for tallying surface collisions/reactions
1305 done at end of each timestep
1306 this is done within individual SurfCollide and SurfReact instances
1307 ------------------------------------------------------------------------- */
1308
collide_react_update()1309 void Update::collide_react_update()
1310 {
1311 for (int i = 0; i < nsc; i++) sc[i]->tally_update();
1312 for (int i = 0; i < nsr; i++) sr[i]->tally_update();
1313 }
1314
1315 /* ----------------------------------------------------------------------
1316 setup lists of all computes that tally surface and boundary bounce info
1317 return 1 if there are any, 0 if not
1318 ------------------------------------------------------------------------- */
1319
bounce_setup()1320 int Update::bounce_setup()
1321 {
1322 delete [] slist_compute;
1323 delete [] blist_compute;
1324 delete [] slist_active;
1325 delete [] blist_active;
1326 slist_compute = blist_compute = NULL;
1327
1328 nslist_compute = nblist_compute = 0;
1329 for (int i = 0; i < modify->ncompute; i++) {
1330 if (modify->compute[i]->surf_tally_flag) nslist_compute++;
1331 if (modify->compute[i]->boundary_tally_flag) nblist_compute++;
1332 }
1333
1334 if (nslist_compute) slist_compute = new Compute*[nslist_compute];
1335 if (nblist_compute) blist_compute = new Compute*[nblist_compute];
1336 if (nslist_compute) slist_active = new Compute*[nslist_compute];
1337 if (nblist_compute) blist_active = new Compute*[nblist_compute];
1338
1339 nslist_compute = nblist_compute = 0;
1340 for (int i = 0; i < modify->ncompute; i++) {
1341 if (modify->compute[i]->surf_tally_flag)
1342 slist_compute[nslist_compute++] = modify->compute[i];
1343 if (modify->compute[i]->boundary_tally_flag)
1344 blist_compute[nblist_compute++] = modify->compute[i];
1345 }
1346
1347 if (nslist_compute || nblist_compute) return 1;
1348 nsurf_tally = nboundary_tally = 0;
1349 return 0;
1350 }
1351
1352 /* ----------------------------------------------------------------------
1353 set bounce tally flags for current timestep
1354 nsurf_tally = # of surface computes needing bounce info on this step
1355 nboundary_tally = # of boundary computes needing bounce info on this step
1356 clear accumulators in computes that will be invoked this step
1357 ------------------------------------------------------------------------- */
1358
bounce_set(bigint ntimestep)1359 void Update::bounce_set(bigint ntimestep)
1360 {
1361 int i;
1362
1363 nsurf_tally = 0;
1364 if (nslist_compute) {
1365 for (i = 0; i < nslist_compute; i++)
1366 if (slist_compute[i]->matchstep(ntimestep)) {
1367 slist_active[nsurf_tally++] = slist_compute[i];
1368 slist_compute[i]->clear();
1369 }
1370 }
1371
1372 nboundary_tally = 0;
1373 if (nblist_compute) {
1374 for (i = 0; i < nblist_compute; i++)
1375 if (blist_compute[i]->matchstep(ntimestep)) {
1376 blist_active[nboundary_tally++] = blist_compute[i];
1377 blist_compute[i]->clear();
1378 }
1379 }
1380 }
1381
1382 /* ----------------------------------------------------------------------
1383 make list of classes that reset dynamic parameters
1384 currently only surf collision models
1385 ------------------------------------------------------------------------- */
1386
dynamic_setup()1387 void Update::dynamic_setup()
1388 {
1389 delete [] ulist_surfcollide;
1390 ulist_surfcollide = NULL;
1391
1392 nulist_surfcollide = 0;
1393 for (int i = 0; i < surf->nsc; i++)
1394 if (surf->sc[i]->dynamicflag) nulist_surfcollide++;
1395
1396 if (nulist_surfcollide)
1397 ulist_surfcollide = new SurfCollide*[nulist_surfcollide];
1398
1399 nulist_surfcollide = 0;
1400 for (int i = 0; i < surf->nsc; i++)
1401 if (surf->sc[i]->dynamicflag)
1402 ulist_surfcollide[nulist_surfcollide++] = surf->sc[i];
1403
1404 if (nulist_surfcollide) dynamic = 1;
1405 }
1406
1407 /* ----------------------------------------------------------------------
1408 invoke class methods that reset dynamic parameters
1409 ------------------------------------------------------------------------- */
1410
dynamic_update()1411 void Update::dynamic_update()
1412 {
1413 if (nulist_surfcollide) {
1414 for (int i = 0; i < nulist_surfcollide; i++)
1415 ulist_surfcollide[i]->dynamic();
1416 }
1417 }
1418
1419 /* ----------------------------------------------------------------------
1420 set global properites via global command in input script
1421 ------------------------------------------------------------------------- */
1422
global(int narg,char ** arg)1423 void Update::global(int narg, char **arg)
1424 {
1425 if (narg < 1) error->all(FLERR,"Illegal global command");
1426
1427 int iarg = 0;
1428 while (iarg < narg) {
1429 if (strcmp(arg[iarg],"fnum") == 0) {
1430 if (iarg+2 > narg) error->all(FLERR,"Illegal global command");
1431 fnum = input->numeric(FLERR,arg[iarg+1]);
1432 if (fnum <= 0.0) error->all(FLERR,"Illegal global command");
1433 iarg += 2;
1434 } else if (strcmp(arg[iarg],"nrho") == 0) {
1435 if (iarg+2 > narg) error->all(FLERR,"Illegal global command");
1436 nrho = input->numeric(FLERR,arg[iarg+1]);
1437 if (nrho <= 0.0) error->all(FLERR,"Illegal global command");
1438 iarg += 2;
1439 } else if (strcmp(arg[iarg],"vstream") == 0) {
1440 if (iarg+4 > narg) error->all(FLERR,"Illegal global command");
1441 vstream[0] = input->numeric(FLERR,arg[iarg+1]);
1442 vstream[1] = input->numeric(FLERR,arg[iarg+2]);
1443 vstream[2] = input->numeric(FLERR,arg[iarg+3]);
1444 iarg += 4;
1445 } else if (strcmp(arg[iarg],"temp") == 0) {
1446 if (iarg+2 > narg) error->all(FLERR,"Illegal global command");
1447 temp_thermal = input->numeric(FLERR,arg[iarg+1]);
1448 if (temp_thermal <= 0.0) error->all(FLERR,"Illegal global command");
1449 iarg += 2;
1450 } else if (strcmp(arg[iarg],"gravity") == 0) {
1451 if (iarg+5 > narg) error->all(FLERR,"Illegal global command");
1452 double gmag = input->numeric(FLERR,arg[iarg+1]);
1453 gravity[0] = input->numeric(FLERR,arg[iarg+2]);
1454 gravity[1] = input->numeric(FLERR,arg[iarg+3]);
1455 gravity[2] = input->numeric(FLERR,arg[iarg+4]);
1456 if (gmag < 0.0) error->all(FLERR,"Illegal global command");
1457 if (gmag > 0.0 &&
1458 gravity[0] == 0.0 && gravity[1] == 0.0 && gravity[2] == 0.0)
1459 error->all(FLERR,"Illegal global command");
1460 if (gmag > 0.0) MathExtra::snorm3(gmag,gravity);
1461 else gravity[0] = gravity[1] = gravity[2] = 0.0;
1462 iarg += 5;
1463
1464 } else if (strcmp(arg[iarg],"surfs") == 0) {
1465 if (iarg+2 > narg) error->all(FLERR,"Illegal global command");
1466 surf->global(arg[iarg+1]);
1467 iarg += 2;
1468 } else if (strcmp(arg[iarg],"surfgrid") == 0) {
1469 if (iarg+2 > narg) error->all(FLERR,"Illegal global command");
1470 if (surf->exist)
1471 error->all(FLERR,
1472 "Cannot set global surfgrid when surfaces already exist");
1473 if (strcmp(arg[iarg+1],"auto") == 0) grid->surfgrid_algorithm = PERAUTO;
1474 else if (strcmp(arg[iarg+1],"percell") == 0)
1475 grid->surfgrid_algorithm = PERCELL;
1476 else if (strcmp(arg[iarg+1],"persurf") == 0)
1477 grid->surfgrid_algorithm = PERSURF;
1478 else error->all(FLERR,"Illegal global command");
1479 iarg += 2;
1480 } else if (strcmp(arg[iarg],"surfmax") == 0) {
1481 if (iarg+2 > narg) error->all(FLERR,"Illegal global command");
1482 if (surf->exist)
1483 error->all(FLERR,
1484 "Cannot set global surfmax when surfaces already exist");
1485 grid->maxsurfpercell = atoi(arg[iarg+1]);
1486 if (grid->maxsurfpercell <= 0) error->all(FLERR,"Illegal global command");
1487 // reallocate paged data structs for variable-length surf info
1488 grid->allocate_surf_arrays();
1489 iarg += 2;
1490 } else if (strcmp(arg[iarg],"splitmax") == 0) {
1491 if (iarg+2 > narg) error->all(FLERR,"Illegal global command");
1492 if (surf->exist)
1493 error->all(FLERR,
1494 "Cannot set global splitmax when surfaces already exist");
1495 grid->maxsplitpercell = atoi(arg[iarg+1]);
1496 if (grid->maxsplitpercell <= 0) error->all(FLERR,"Illegal global command");
1497 // reallocate paged data structs for variable-length cell info
1498 grid->allocate_surf_arrays();
1499 iarg += 2;
1500
1501 } else if (strcmp(arg[iarg],"gridcut") == 0) {
1502 if (iarg+2 > narg) error->all(FLERR,"Illegal global command");
1503 grid->cutoff = input->numeric(FLERR,arg[iarg+1]);
1504 if (grid->cutoff < 0.0 && grid->cutoff != -1.0)
1505 error->all(FLERR,"Illegal global command");
1506 // force ghost info to be regenerated with new cutoff
1507 grid->remove_ghosts();
1508 iarg += 2;
1509
1510 } else if (strcmp(arg[iarg],"comm/sort") == 0) {
1511 if (iarg+2 > narg) error->all(FLERR,"Illegal global command");
1512 if (strcmp(arg[iarg+1],"yes") == 0) comm->commsortflag = 1;
1513 else if (strcmp(arg[iarg+1],"no") == 0) comm->commsortflag = 0;
1514 else error->all(FLERR,"Illegal global command");
1515 iarg += 2;
1516 } else if (strcmp(arg[iarg],"comm/style") == 0) {
1517 if (iarg+2 > narg) error->all(FLERR,"Illegal global command");
1518 if (strcmp(arg[iarg+1],"neigh") == 0) comm->commpartstyle = 1;
1519 else if (strcmp(arg[iarg+1],"all") == 0) comm->commpartstyle = 0;
1520 else error->all(FLERR,"Illegal global command");
1521 iarg += 2;
1522 } else if (strcmp(arg[iarg],"surftally") == 0) {
1523 if (iarg+2 > narg) error->all(FLERR,"Illegal global command");
1524 if (strcmp(arg[iarg+1],"auto") == 0) surf->tally_comm = TALLYAUTO;
1525 else if (strcmp(arg[iarg+1],"reduce") == 0) surf->tally_comm = TALLYREDUCE;
1526 else if (strcmp(arg[iarg+1],"rvous") == 0) surf->tally_comm = TALLYRVOUS;
1527 else error->all(FLERR,"Illegal global command");
1528 iarg += 2;
1529
1530 } else if (strcmp(arg[iarg],"weight") == 0) {
1531 // for now assume just one arg after "cell"
1532 // may need to generalize later
1533 if (iarg+3 > narg) error->all(FLERR,"Illegal global command");
1534 if (strcmp(arg[iarg+1],"cell") == 0) grid->weight(1,&arg[iarg+2]);
1535 else error->all(FLERR,"Illegal weight command");
1536 iarg += 3;
1537 } else if (strcmp(arg[iarg],"particle/reorder") == 0) {
1538 reorder_period = input->inumeric(FLERR,arg[iarg+1]);
1539 if (reorder_period < 0) error->all(FLERR,"Illegal global command");
1540 iarg += 2;
1541 } else if (strcmp(arg[iarg],"mem/limit") == 0) {
1542 if (iarg+2 > narg) error->all(FLERR,"Illegal global command");
1543 if (strcmp(arg[iarg+1],"grid") == 0) mem_limit_grid_flag = 1;
1544 else {
1545 double factor = input->numeric(FLERR,arg[iarg+1]);
1546 bigint global_mem_limit_big = static_cast<bigint> (factor * 1024*1024);
1547 if (global_mem_limit_big < 0) error->all(FLERR,"Illegal global command");
1548 if (global_mem_limit_big > MAXSMALLINT)
1549 error->all(FLERR,"Global mem/limit setting cannot exceed 2GB");
1550 global_mem_limit = global_mem_limit_big;
1551 }
1552 iarg += 2;
1553 } else error->all(FLERR,"Illegal global command");
1554 }
1555 }
1556
1557 /* ----------------------------------------------------------------------
1558 reset timestep as called from input script
1559 ------------------------------------------------------------------------- */
1560
reset_timestep(int narg,char ** arg)1561 void Update::reset_timestep(int narg, char **arg)
1562 {
1563 if (narg != 1) error->all(FLERR,"Illegal reset_timestep command");
1564 bigint newstep = ATOBIGINT(arg[0]);
1565 reset_timestep(newstep);
1566 }
1567
1568 /* ----------------------------------------------------------------------
1569 reset timestep
1570 set atimestep to new timestep, so future update_time() calls will be correct
1571 trigger reset of timestep for output and for fixes that require it
1572 do not allow any timestep-dependent fixes to be defined
1573 reset eflag/vflag global so nothing will think eng/virial are current
1574 reset invoked flags of computes,
1575 so nothing will think they are current between runs
1576 clear timestep list of computes that store future invocation times
1577 called from rerun command and input script (indirectly)
1578 ------------------------------------------------------------------------- */
1579
reset_timestep(bigint newstep)1580 void Update::reset_timestep(bigint newstep)
1581 {
1582 ntimestep = newstep;
1583 if (ntimestep < 0) error->all(FLERR,"Timestep must be >= 0");
1584 if (ntimestep > MAXBIGINT) error->all(FLERR,"Too big a timestep");
1585
1586 output->reset_timestep(ntimestep);
1587
1588 for (int i = 0; i < modify->nfix; i++) {
1589 if (modify->fix[i]->time_depend)
1590 error->all(FLERR,
1591 "Cannot reset timestep with a time-dependent fix defined");
1592 //modify->fix[i]->reset_timestep(ntimestep);
1593 }
1594
1595 for (int i = 0; i < modify->ncompute; i++) {
1596 modify->compute[i]->invoked_scalar = -1;
1597 modify->compute[i]->invoked_vector = -1;
1598 modify->compute[i]->invoked_array = -1;
1599 modify->compute[i]->invoked_per_particle = -1;
1600 modify->compute[i]->invoked_per_grid = -1;
1601 modify->compute[i]->invoked_per_surf = -1;
1602 }
1603
1604 for (int i = 0; i < modify->ncompute; i++)
1605 if (modify->compute[i]->timeflag) modify->compute[i]->clearstep();
1606 }
1607
1608 /* ----------------------------------------------------------------------
1609 get mem/limit based on grid memory
1610 ------------------------------------------------------------------------- */
1611
set_mem_limit_grid(int gnlocal)1612 void Update::set_mem_limit_grid(int gnlocal)
1613 {
1614 if (gnlocal == 0) gnlocal = grid->nlocal;
1615
1616 bigint global_mem_limit_big = static_cast<bigint> (gnlocal*sizeof(Grid::ChildCell));
1617
1618 if (global_mem_limit_big > MAXSMALLINT)
1619 error->one(FLERR,"Global mem/limit setting cannot exceed 2GB");
1620
1621 global_mem_limit = global_mem_limit_big;
1622 }
1623
1624 /* ----------------------------------------------------------------------
1625 get mem/limit based on grid memory
1626 ------------------------------------------------------------------------- */
1627
have_mem_limit()1628 int Update::have_mem_limit()
1629 {
1630 if (mem_limit_grid_flag)
1631 set_mem_limit_grid();
1632
1633 int mem_limit_flag = 0;
1634
1635 if (global_mem_limit > 0 || (mem_limit_grid_flag && !grid->nlocal))
1636 mem_limit_flag = 1;
1637
1638 return mem_limit_flag;
1639 }
1640