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 "string.h"
16 #include "grid.h"
17 #include "geometry.h"
18 #include "domain.h"
19 #include "region.h"
20 #include "surf.h"
21 #include "comm.h"
22 #include "modify.h"
23 #include "fix.h"
24 #include "compute.h"
25 #include "output.h"
26 #include "dump.h"
27 #include "irregular.h"
28 #include "math_const.h"
29 #include "memory.h"
30 #include "error.h"
31
32 // DEBUG
33 #include "update.h"
34
35 using namespace SPARTA_NS;
36 using namespace MathConst;
37
38 #define DELTA 8192
39 #define DELTAPARENT 1024
40 #define LARGE 256000
41 #define BIG 1.0e20
42 #define MAXGROUP 32
43 #define MAXLEVEL 32
44
45 // default values, can be overridden by global command
46
47 #define MAXSURFPERCELL 100
48 #define MAXSPLITPERCELL 10
49
50 enum{XLO,XHI,YLO,YHI,ZLO,ZHI,INTERIOR}; // same as Domain
51 enum{PERIODIC,OUTFLOW,REFLECT,SURFACE,AXISYM}; // same as Domain
52 enum{REGION_ALL,REGION_ONE,REGION_CENTER}; // same as Surf
53 enum{PERAUTO,PERCELL,PERSURF}; // several files
54
55 // cell type = OUTSIDE/INSIDE/OVERLAP if entirely outside/inside surfs
56 // or has any overlap with surfs including grazing or touching
57 // corner point = OUTSIDE/INSIDE (not OVERLAP) if outside/inside
58 // if exactly on surf, is still marked OUTSIDE or INSIDE by cut2d/cut3d
59 // corner pts are only set if cell type = OVERLAP
60
61 enum{UNKNOWN,OUTSIDE,INSIDE,OVERLAP}; // several files
62 enum{NCHILD,NPARENT,NUNKNOWN,NPBCHILD,NPBPARENT,NPBUNKNOWN,NBOUND}; // Update
63 enum{NOWEIGHT,VOLWEIGHT,RADWEIGHT,RADONLYWEIGHT};
64
65 // corners[i][j] = J corner points of face I of a grid cell
66 // works for 2d quads and 3d hexes
67
68 int corners[6][4] = {{0,2,4,6}, {1,3,5,7}, {0,1,4,5}, {2,3,6,7},
69 {0,1,2,3}, {4,5,6,7}};
70
71 /* ---------------------------------------------------------------------- */
72
Grid(SPARTA * sparta)73 Grid::Grid(SPARTA *sparta) : Pointers(sparta)
74 {
75 exist = exist_ghost = clumped = 0;
76 MPI_Comm_rank(world,&me);
77
78 gnames = (char **) memory->smalloc(MAXGROUP*sizeof(char *),"grid:gnames");
79 bitmask = (int *) memory->smalloc(MAXGROUP*sizeof(int),"grid:bitmask");
80 inversemask = (int *) memory->smalloc(MAXGROUP*sizeof(int),
81 "grid:inversemask");
82
83 for (int i = 0; i < MAXGROUP; i++) bitmask[i] = 1 << i;
84 for (int i = 0; i < MAXGROUP; i++) inversemask[i] = bitmask[i] ^ ~0;
85
86 ngroup = 1;
87 int n = strlen("all") + 1;
88 gnames[0] = new char[n];
89 strcpy(gnames[0],"all");
90
91 ncell = nunsplit = nsplit = nsub = 0;
92
93 nlocal = nghost = maxlocal = maxcell = 0;
94 nsplitlocal = nsplitghost = maxsplit = 0;
95 nsublocal = nsubghost = 0;
96 nparent = maxparent = 0;
97 maxlevel = 0;
98 plevel_limit = MAXLEVEL;
99
100 cells = NULL;
101 cinfo = NULL;
102 sinfo = NULL;
103 pcells = NULL;
104
105 plevels = new ParentLevel[MAXLEVEL];
106 memset(plevels,0,MAXLEVEL*sizeof(ParentLevel));
107
108 surfgrid_algorithm = PERAUTO;
109 maxsurfpercell = MAXSURFPERCELL;
110 maxsplitpercell = MAXSPLITPERCELL;
111 csurfs = NULL; csplits = NULL; csubs = NULL;
112 allocate_surf_arrays();
113
114 neighshift[XLO] = 0;
115 neighshift[XHI] = 3;
116 neighshift[YLO] = 6;
117 neighshift[YHI] = 9;
118 neighshift[ZLO] = 12;
119 neighshift[ZHI] = 15;
120
121 neighmask[XLO] = 7 << neighshift[XLO];
122 neighmask[XHI] = 7 << neighshift[XHI];
123 neighmask[YLO] = 7 << neighshift[YLO];
124 neighmask[YHI] = 7 << neighshift[YHI];
125 neighmask[ZLO] = 7 << neighshift[ZLO];
126 neighmask[ZHI] = 7 << neighshift[ZHI];
127
128 cutoff = -1.0;
129 cellweightflag = NOWEIGHT;
130
131 // allocate hash for cell IDs
132
133 hash = new MyHash();
134 hashfilled = 0;
135
136 copy = copymode = 0;
137 }
138
139 /* ---------------------------------------------------------------------- */
140
~Grid()141 Grid::~Grid()
142 {
143 if (copy || copymode) return;
144
145 for (int i = 0; i < ngroup; i++) delete [] gnames[i];
146 memory->sfree(gnames);
147 memory->sfree(bitmask);
148 memory->sfree(inversemask);
149
150 memory->sfree(cells);
151 memory->sfree(cinfo);
152 memory->sfree(sinfo);
153 memory->sfree(pcells);
154
155 delete [] plevels;
156
157 delete csurfs;
158 delete csplits;
159 delete csubs;
160 delete hash;
161 }
162
163 /* ----------------------------------------------------------------------
164 remove existing grid
165 called by AdaptGrid when it reads a new grid from file
166 ------------------------------------------------------------------------- */
167
remove()168 void Grid::remove()
169 {
170 memory->sfree(cells);
171 memory->sfree(cinfo);
172 memory->sfree(sinfo);
173 memory->sfree(pcells);
174
175 delete csurfs;
176 delete csplits;
177 delete csubs;
178
179 exist_ghost = clumped = 0;
180 ncell = nunsplit = nsplit = nsub = 0;
181 nlocal = nghost = maxlocal = maxcell = 0;
182 nsplitlocal = nsplitghost = maxsplit = 0;
183 nsublocal = nsubghost = 0;
184 maxlevel = 0;
185
186 hash->clear();
187 hashfilled = 0;
188
189 cells = NULL;
190 cinfo = NULL;
191 sinfo = NULL;
192
193 csurfs = NULL; csplits = NULL; csubs = NULL;
194 allocate_surf_arrays();
195
196 // NOTE: what about cutoff and cellweightflag
197 }
198
199 /* ----------------------------------------------------------------------
200 store copy of Particle class settings
201 ------------------------------------------------------------------------- */
202
init()203 void Grid::init()
204 {
205 ncustom = particle->ncustom;
206 nbytes_particle = sizeof(Particle::OnePart);
207 nbytes_custom = particle->sizeof_custom();
208 nbytes_total = nbytes_particle + nbytes_custom;
209 }
210
211 /* ----------------------------------------------------------------------
212 add a single owned child cell to cells and cinfo
213 assume no surfs at this point, so cell is OUTSIDE, corners are ignored
214 neighs and nmask will be set later
215 ------------------------------------------------------------------------- */
216
add_child_cell(cellint id,int level,double * lo,double * hi)217 void Grid::add_child_cell(cellint id, int level, double *lo, double *hi)
218 {
219 grow_cells(1,1);
220
221 int ncorner;
222 if (domain->dimension == 3) ncorner = 8;
223 else ncorner = 4;
224
225 ChildCell *c = &cells[nlocal];
226 c->id = id;
227 c->proc = me;
228 c->ilocal = nlocal;
229 c->level = level;
230 c->lo[0] = lo[0];
231 c->lo[1] = lo[1];
232 c->lo[2] = lo[2];
233 c->hi[0] = hi[0];
234 c->hi[1] = hi[1];
235 c->hi[2] = hi[2];
236 c->nsurf = 0;
237 c->csurfs = NULL;
238 c->nsplit = 1;
239 c->isplit = -1;
240
241 ChildInfo *ci = &cinfo[nlocal];
242 ci->count = 0;
243 ci->first = -1;
244 ci->mask = 1;
245 ci->type = OUTSIDE;
246 for (int i = 0; i < ncorner; i++) ci->corner[i] = UNKNOWN;
247 ci->weight = 1.0;
248
249 if (domain->dimension == 3)
250 ci->volume = (hi[0]-lo[0]) * (hi[1]-lo[1]) * (hi[2]-lo[2]);
251 else if (domain->axisymmetric)
252 ci->volume = MY_PI * (hi[1]*hi[1]-lo[1]*lo[1]) * (hi[0]-lo[0]);
253 else
254 ci->volume = (hi[0]-lo[0]) * (hi[1]-lo[1]);
255
256 // increment both since are adding an unsplit cell
257
258 nunsplitlocal++;
259 nlocal++;
260 }
261
262 /* ----------------------------------------------------------------------
263 add a single split cell to sinfo
264 ownflag = 1/0 if split cell is owned or ghost
265 sinfo fields will be set by caller
266 ------------------------------------------------------------------------- */
267
add_split_cell(int ownflag)268 void Grid::add_split_cell(int ownflag)
269 {
270 grow_sinfo(1);
271 if (ownflag) nsplitlocal++;
272 else nsplitghost++;
273 }
274
275 /* ----------------------------------------------------------------------
276 add a sub cell to cells and cinfo (if owned)
277 ownflag = 1/0 if sub cell is owned or ghost
278 icell = index of split cell that contains sub cell
279 copy cells and cinfo for split cell to new sub cell
280 sub cell fields that should be different will be set by caller
281 ------------------------------------------------------------------------- */
282
add_sub_cell(int icell,int ownflag)283 void Grid::add_sub_cell(int icell, int ownflag)
284 {
285 grow_cells(1,1);
286
287 int inew;
288 if (ownflag) inew = nlocal;
289 else inew = nlocal + nghost;
290
291 memcpy(&cells[inew],&cells[icell],sizeof(ChildCell));
292 if (ownflag) memcpy(&cinfo[inew],&cinfo[icell],sizeof(ChildInfo));
293
294 if (ownflag) {
295 nsublocal++;
296 nlocal++;
297 } else {
298 nsubghost++;
299 nghost++;
300 }
301 }
302
303 /* ----------------------------------------------------------------------
304 called during a run when per-processor list of grid cells may have changed
305 trigger fixes, computes, dumps to change their allocated per-grid data
306 ------------------------------------------------------------------------- */
307
notify_changed()308 void Grid::notify_changed()
309 {
310 if (modify->n_pergrid) modify->grid_changed();
311
312 Compute **compute = modify->compute;
313 for (int i = 0; i < modify->ncompute; i++) {
314 if (compute[i]->per_grid_flag) compute[i]->reallocate();
315 if (compute[i]->per_surf_flag) compute[i]->reallocate();
316 }
317
318 for (int i = 0; i < output->ndump; i++)
319 output->dump[i]->reset_grid_count();
320 }
321
322 /* ----------------------------------------------------------------------
323 return minlevel = minlevel of any grid cell
324 ------------------------------------------------------------------------- */
325
set_minlevel()326 int Grid::set_minlevel()
327 {
328 int mylevel = MAXLEVEL;
329 for (int i = 0; i < nlocal; i++)
330 mylevel = MIN(mylevel,cells[i].level);
331
332 int minlevel;
333 MPI_Allreduce(&mylevel,&minlevel,1,MPI_INT,MPI_MIN,world);
334 return minlevel;
335 }
336
337 /* ----------------------------------------------------------------------
338 set maxlevel = maxlevel of any grid cell
339 ------------------------------------------------------------------------- */
340
set_maxlevel()341 void Grid::set_maxlevel()
342 {
343 int mylevel = 0;
344 for (int i = 0; i < nlocal; i++)
345 mylevel = MAX(mylevel,cells[i].level);
346
347 MPI_Allreduce(&mylevel,&maxlevel,1,MPI_INT,MPI_MAX,world);
348 }
349
350 /* ----------------------------------------------------------------------
351 set grid stats for owned cells
352 ------------------------------------------------------------------------- */
353
setup_owned()354 void Grid::setup_owned()
355 {
356 // global counts for ncell, nunsplit, nsplit, nsub
357
358 bigint one = nlocal - nsublocal;
359 MPI_Allreduce(&one,&ncell,1,MPI_SPARTA_BIGINT,MPI_SUM,world);
360 one = nunsplitlocal = nlocal - nsplitlocal - nsublocal;
361 MPI_Allreduce(&one,&nunsplit,1,MPI_SPARTA_BIGINT,MPI_SUM,world);
362 MPI_Allreduce(&nsplitlocal,&nsplit,1,MPI_INT,MPI_SUM,world);
363 MPI_Allreduce(&nsublocal,&nsub,1,MPI_INT,MPI_SUM,world);
364
365 // set cell_epsilon to 1/2 the smallest dimension of any grid cell
366
367 int dimension = domain->dimension;
368
369 double eps = BIG;
370 for (int i = 0; i < nlocal; i++) {
371 if (cells[i].nsplit <= 0) continue;
372 eps = MIN(eps,cells[i].hi[0]-cells[i].lo[0]);
373 eps = MIN(eps,cells[i].hi[1]-cells[i].lo[1]);
374 if (dimension == 3) eps = MIN(eps,cells[i].hi[2]-cells[i].lo[2]);
375 }
376
377 MPI_Allreduce(&eps,&cell_epsilon,1,MPI_DOUBLE,MPI_MIN,world);
378 cell_epsilon *= 0.5;
379 }
380
381 /* ----------------------------------------------------------------------
382 remove ghost grid cells and any allocated data they have
383 currently, cells just have ptrs into pages that are deallocated separately
384 also remove ghost surfaces, either explicit or implicit
385 ------------------------------------------------------------------------- */
386
remove_ghosts()387 void Grid::remove_ghosts()
388 {
389 hashfilled = 0;
390 exist_ghost = 0;
391 nghost = nunsplitghost = nsplitghost = nsubghost = 0;
392 surf->remove_ghosts();
393 }
394
395 /* ----------------------------------------------------------------------
396 acquire ghost cells from local cells of other procs
397 if surfs are distributed, also acquire ghost cell surfs
398 explicit distributed surfs require use of hash
399 method used depends on ghost cutoff
400 no-op if grid is not clumped and want to acquire only nearby ghosts
401 ------------------------------------------------------------------------- */
402
acquire_ghosts(int surfflag)403 void Grid::acquire_ghosts(int surfflag)
404 {
405 if (surf->distributed && !surf->implicit) surf->rehash();
406
407 if (cutoff < 0.0) acquire_ghosts_all(surfflag);
408 else if (clumped) acquire_ghosts_near(surfflag);
409 else if (comm->me == 0)
410 error->warning(FLERR,"Could not acquire nearby ghost cells b/c "
411 "grid partition is not clumped");
412
413 if (surf->distributed && !surf->implicit) {
414 surf->hash->clear();
415 surf->hashfilled = 0;
416 }
417 }
418
419 /* ----------------------------------------------------------------------
420 acquire ghost cells from local cells of other procs
421 use ring comm to get copy of all other cells in global system
422 ------------------------------------------------------------------------- */
423
acquire_ghosts_all(int surfflag)424 void Grid::acquire_ghosts_all(int surfflag)
425 {
426 exist_ghost = 1;
427 nempty = 0;
428
429 // compute total # of ghosts so can pre-allocate cells array
430 // issue a memory warning if grid cell count >= LARGE and
431 // user has not specified a grid cutoff
432
433 int nghost_new;
434 MPI_Allreduce(&nlocal,&nghost_new,1,MPI_INT,MPI_SUM,world);
435
436 if (nghost_new >= LARGE && comm->nprocs > 1 && comm->me == 0)
437 error->warning(FLERR,"Per-processor grid cell memory will be large "
438 "because global gridcut < 0.0");
439
440 nghost_new -= nlocal;
441 grow_cells(nghost_new,0);
442
443 // create buf for holding all of my cells, not including sub cells
444
445 bigint bsendsize = 0;
446 for (int icell = 0; icell < nlocal; icell++) {
447 if (cells[icell].nsplit <= 0) continue;
448 bsendsize += pack_one(icell,NULL,0,0,surfflag,0);
449 }
450
451 if (bsendsize > MAXSMALLINT)
452 error->one(FLERR,"Acquire ghosts all send buffer exceeds 2 GB");
453 int sendsize = bsendsize;
454
455 char *sbuf;
456 memory->create(sbuf,sendsize,"grid:sbuf");
457 memset(sbuf,0,sendsize);
458
459 // pack each unsplit or split cell
460 // subcells will be packed by split cell
461
462 sendsize = 0;
463 for (int icell = 0; icell < nlocal; icell++) {
464 if (cells[icell].nsplit <= 0) continue;
465 sendsize += pack_one(icell,&sbuf[sendsize],0,0,surfflag,1);
466 }
467
468 // circulate buf of my grid cells around to all procs
469 // unpack augments my ghost cells with info from other procs
470
471 unpack_ghosts_surfflag = surfflag;
472 comm->ring(sendsize,sizeof(char),sbuf,1,unpack_ghosts,NULL,0,(void *) this);
473
474 memory->destroy(sbuf);
475 }
476
477 /* ----------------------------------------------------------------------
478 acquire ghost cells from local cells of other procs
479 use irregular comm to only get copy of nearby cells
480 within extended bounding box = bounding box of owned cells + cutoff
481 ------------------------------------------------------------------------- */
482
acquire_ghosts_near(int surfflag)483 void Grid::acquire_ghosts_near(int surfflag)
484 {
485 if (update->have_mem_limit())
486 return acquire_ghosts_near_less_memory(surfflag);
487
488 exist_ghost = 1;
489
490 // bb lo/hi = bounding box for my owned cells
491
492 int i;
493 double bblo[3],bbhi[3];
494 double *lo,*hi;
495
496 for (i = 0; i < 3; i++) {
497 bblo[i] = BIG;
498 bbhi[i] = -BIG;
499 }
500
501 for (int icell = 0; icell < nlocal; icell++) {
502 if (cells[icell].nsplit <= 0) continue;
503 lo = cells[icell].lo;
504 hi = cells[icell].hi;
505 for (i = 0; i < 3; i++) {
506 bblo[i] = MIN(bblo[i],lo[i]);
507 bbhi[i] = MAX(bbhi[i],hi[i]);
508 }
509 }
510
511 // ebb lo/hi = bbox + grid cutoff
512 // trim to simulation box in non-periodic dims
513 // if bblo/hi is at periodic boundary and cutoff is 0.0,
514 // add cell_epsilon to insure ghosts across periodic boundary acquired,
515 // else may be UNKNOWN to owned cell
516
517 double *boxlo = domain->boxlo;
518 double *boxhi = domain->boxhi;
519 int *bflag = domain->bflag;
520
521 double ebblo[3],ebbhi[3];
522 for (i = 0; i < 3; i++) {
523 ebblo[i] = bblo[i] - cutoff;
524 ebbhi[i] = bbhi[i] + cutoff;
525 if (bflag[2*i] != PERIODIC) ebblo[i] = MAX(ebblo[i],boxlo[i]);
526 if (bflag[2*i] != PERIODIC) ebbhi[i] = MIN(ebbhi[i],boxhi[i]);
527 if (bflag[2*i] == PERIODIC && bblo[i] == boxlo[i] && cutoff == 0.0)
528 ebblo[i] -= cell_epsilon;
529 if (bflag[2*i] == PERIODIC && bbhi[i] == boxhi[i] && cutoff == 0.0)
530 ebbhi[i] += cell_epsilon;
531 }
532
533 // box = ebbox split across periodic BC
534 // 27 is max number of periodic images in 3d
535
536 Box box[27];
537 int nbox = box_periodic(ebblo,ebbhi,box);
538
539 // boxall = collection of boxes from all procs
540
541 int me = comm->me;
542 int nprocs = comm->nprocs;
543
544 int nboxall;
545 MPI_Allreduce(&nbox,&nboxall,1,MPI_INT,MPI_SUM,world);
546
547 int *recvcounts,*displs;
548 memory->create(recvcounts,nprocs,"grid:recvcounts");
549 memory->create(displs,nprocs,"grid:displs");
550
551 int nsend = nbox*sizeof(Box);
552 MPI_Allgather(&nsend,1,MPI_INT,recvcounts,1,MPI_INT,world);
553 displs[0] = 0;
554 for (i = 1; i < nprocs; i++) displs[i] = displs[i-1] + recvcounts[i-1];
555
556 Box *boxall = new Box[nboxall];
557 MPI_Allgatherv(box,nsend,MPI_CHAR,boxall,recvcounts,displs,MPI_CHAR,world);
558
559 memory->destroy(recvcounts);
560 memory->destroy(displs);
561
562 // nlist = # of boxes that overlap with my bbox, skipping self boxes
563 // list = indices into boxall of overlaps
564 // overlap = true overlap or just touching
565
566 int nlist = 0;
567 int *list;
568 memory->create(list,nboxall,"grid:list");
569
570 for (i = 0; i < nboxall; i++) {
571 if (boxall[i].proc == me) continue;
572 if (box_overlap(bblo,bbhi,boxall[i].lo,boxall[i].hi)) list[nlist++] = i;
573 }
574
575 // loop over my owned cells, not including sub cells
576 // each may overlap with multiple boxes in list
577 // on 1st pass, just tally memory to send copies of my cells
578 // use lastproc to insure a cell only overlaps once per other proc
579 // if oflag = 2 = my cell just touches box,
580 // so flag grid cell as EMPTY ghost by setting nsurf = -1
581
582 int j,oflag,lastproc,nsurf_hold;
583
584 nsend = 0;
585 bigint bsendsize = 0;
586
587 for (int icell = 0; icell < nlocal; icell++) {
588 if (cells[icell].nsplit <= 0) continue;
589 lo = cells[icell].lo;
590 hi = cells[icell].hi;
591 lastproc = -1;
592 for (i = 0; i < nlist; i++) {
593 j = list[i];
594 oflag = box_overlap(lo,hi,boxall[j].lo,boxall[j].hi);
595 if (!oflag) continue;
596 if (boxall[j].proc == lastproc) continue;
597 lastproc = boxall[j].proc;
598
599 if (oflag == 2) {
600 nsurf_hold = cells[icell].nsurf;
601 cells[icell].nsurf = -1;
602 }
603 bsendsize += pack_one(icell,NULL,0,0,surfflag,0);
604 if (oflag == 2) cells[icell].nsurf = nsurf_hold;
605 nsend++;
606 }
607 }
608
609 if (bsendsize > MAXSMALLINT)
610 error->one(FLERR,"Acquire ghosts near send buffer exceeds 2 GB");
611 int sendsize = bsendsize;
612
613 // create send buf and auxiliary irregular comm vectors
614
615 char *sbuf;
616 memory->create(sbuf,sendsize,"grid:sbuf");
617 memset(sbuf,0,sendsize);
618
619 int *proclist,*sizelist;
620 memory->create(proclist,nsend,"grid:proclist");
621 memory->create(sizelist,nsend,"grid:sizelist");
622
623 // on 2nd pass over local cells, fill the send buf
624 // use lastproc to insure a cell only overlaps once per other proc
625 // if oflag = 2 = my cell just touches box,
626 // so flag grid cell as EMPTY ghost by setting nsurf = -1
627
628 nsend = 0;
629 sendsize = 0;
630 for (int icell = 0; icell < nlocal; icell++) {
631 if (cells[icell].nsplit <= 0) continue;
632 lo = cells[icell].lo;
633 hi = cells[icell].hi;
634 lastproc = -1;
635 for (i = 0; i < nlist; i++) {
636 j = list[i];
637 oflag = box_overlap(lo,hi,boxall[j].lo,boxall[j].hi);
638 if (!oflag) continue;
639 if (boxall[j].proc == lastproc) continue;
640 lastproc = boxall[j].proc;
641
642 if (oflag == 2) {
643 nsurf_hold = cells[icell].nsurf;
644 cells[icell].nsurf = -1;
645 }
646 sizelist[nsend] = pack_one(icell,&sbuf[sendsize],0,0,surfflag,1);
647 if (oflag == 2) cells[icell].nsurf = nsurf_hold;
648 proclist[nsend] = lastproc;
649 sendsize += sizelist[nsend];
650 nsend++;
651 }
652 }
653
654 // clean up
655
656 memory->destroy(list);
657 delete [] boxall;
658
659 // perform irregular communication of list of ghost cells
660
661 Irregular *irregular = new Irregular(sparta);
662 int recvsize;
663 int nrecv = irregular->create_data_variable(nsend,proclist,sizelist,
664 recvsize,comm->commsortflag);
665
666 char *rbuf;
667 memory->create(rbuf,recvsize,"grid:rbuf");
668 memset(rbuf,0,recvsize);
669
670 irregular->exchange_variable(sbuf,sizelist,rbuf);
671 delete irregular;
672
673 // unpack received grid cells as ghost cells
674
675 int offset = 0;
676 for (i = 0; i < nrecv; i++)
677 offset += grid->unpack_one(&rbuf[offset],0,0,surfflag);
678
679 // more clean up
680
681 memory->destroy(proclist);
682 memory->destroy(sizelist);
683 memory->destroy(sbuf);
684 memory->destroy(rbuf);
685
686 // set nempty = # of EMPTY ghost cells I store
687
688 nempty = 0;
689 for (int icell = nlocal; icell < nlocal+nghost; icell++)
690 if (cells[icell].nsurf < 0) nempty++;
691 }
692
693 /* ----------------------------------------------------------------------
694 acquire ghost cells from local cells of other procs
695 use irregular comm to only get copy of nearby cells
696 within extended bounding box = bounding box of owned cells + cutoff
697 ------------------------------------------------------------------------- */
698
acquire_ghosts_near_less_memory(int surfflag)699 void Grid::acquire_ghosts_near_less_memory(int surfflag)
700 {
701 exist_ghost = 1;
702
703 // bb lo/hi = bounding box for my owned cells
704
705 int i;
706 double bblo[3],bbhi[3];
707 double *lo,*hi;
708
709 for (i = 0; i < 3; i++) {
710 bblo[i] = BIG;
711 bbhi[i] = -BIG;
712 }
713
714 for (int icell = 0; icell < nlocal; icell++) {
715 if (cells[icell].nsplit <= 0) continue;
716 lo = cells[icell].lo;
717 hi = cells[icell].hi;
718 for (i = 0; i < 3; i++) {
719 bblo[i] = MIN(bblo[i],lo[i]);
720 bbhi[i] = MAX(bbhi[i],hi[i]);
721 }
722 }
723
724 // ebb lo/hi = bbox + grid cutoff
725 // trim to simulation box in non-periodic dims
726 // if bblo/hi is at periodic boundary and cutoff is 0.0,
727 // add cell_epsilon to insure ghosts across periodic boundary acquired,
728 // else may be UNKNOWN to owned cell
729
730 double *boxlo = domain->boxlo;
731 double *boxhi = domain->boxhi;
732 int *bflag = domain->bflag;
733
734 double ebblo[3],ebbhi[3];
735 for (i = 0; i < 3; i++) {
736 ebblo[i] = bblo[i] - cutoff;
737 ebbhi[i] = bbhi[i] + cutoff;
738 if (bflag[2*i] != PERIODIC) ebblo[i] = MAX(ebblo[i],boxlo[i]);
739 if (bflag[2*i] != PERIODIC) ebbhi[i] = MIN(ebbhi[i],boxhi[i]);
740 if (bflag[2*i] == PERIODIC && bblo[i] == boxlo[i] && cutoff == 0.0)
741 ebblo[i] -= cell_epsilon;
742 if (bflag[2*i] == PERIODIC && bbhi[i] == boxhi[i] && cutoff == 0.0)
743 ebbhi[i] += cell_epsilon;
744 }
745
746 // box = ebbox split across periodic BC
747 // 27 is max number of periodic images in 3d
748
749 Box box[27];
750 int nbox = box_periodic(ebblo,ebbhi,box);
751
752 // boxall = collection of boxes from all procs
753
754 int me = comm->me;
755 int nprocs = comm->nprocs;
756
757 int nboxall;
758 MPI_Allreduce(&nbox,&nboxall,1,MPI_INT,MPI_SUM,world);
759
760 int *recvcounts,*displs;
761 memory->create(recvcounts,nprocs,"grid:recvcounts");
762 memory->create(displs,nprocs,"grid:displs");
763
764 int nsend = nbox*sizeof(Box);
765 MPI_Allgather(&nsend,1,MPI_INT,recvcounts,1,MPI_INT,world);
766 displs[0] = 0;
767 for (i = 1; i < nprocs; i++) displs[i] = displs[i-1] + recvcounts[i-1];
768
769 Box *boxall = new Box[nboxall];
770 MPI_Allgatherv(box,nsend,MPI_CHAR,boxall,recvcounts,displs,MPI_CHAR,world);
771
772 memory->destroy(recvcounts);
773 memory->destroy(displs);
774
775 // nlist = # of boxes that overlap with my bbox, skipping self boxes
776 // list = indices into boxall of overlaps
777 // overlap = true overlap or just touching
778
779 int nlist = 0;
780 int *list;
781 memory->create(list,nboxall,"grid:list");
782
783 for (i = 0; i < nboxall; i++) {
784 if (boxall[i].proc == me) continue;
785 if (box_overlap(bblo,bbhi,boxall[i].lo,boxall[i].hi)) list[nlist++] = i;
786 }
787
788 // loop over my owned cells, not including sub cells
789 // each may overlap with multiple boxes in list
790 // on 1st pass, just tally memory to send copies of my cells
791 // use lastproc to insure a cell only overlaps once per other proc
792 // if oflag = 2 = my cell just touches box,
793 // so flag grid cell as EMPTY ghost by setting nsurf = -1
794
795 int j,oflag,lastproc,nsurf_hold;
796
797 int icell_start = 0;
798 int icell_end = nlocal;
799 int i_start = 0;
800 int i_end = nlist;
801
802 int not_done = 1;
803
804 while (not_done) {
805 nsend = 0;
806 bigint bsendsize = 0;
807
808 for (int icell = icell_start; icell < nlocal; icell++) {
809 icell_end = icell+1;
810 if (cells[icell].nsplit <= 0) continue;
811 lo = cells[icell].lo;
812 hi = cells[icell].hi;
813 lastproc = -1;
814 int break_flag = 0;
815 int i_first = 0;
816 if (icell == icell_start) i_first = i_start;
817 for (i = i_first; i < nlist; i++) {
818 i_end = i+1;
819 j = list[i];
820 oflag = box_overlap(lo,hi,boxall[j].lo,boxall[j].hi);
821 if (!oflag) continue;
822 if (boxall[j].proc == lastproc) continue;
823 lastproc = boxall[j].proc;
824
825 int n = pack_one(icell,NULL,0,0,surfflag,0);
826 if (n > 0 && bsendsize > 0 && bsendsize+n > update->global_mem_limit) {
827 i_end = i;
828 break_flag = 1;
829 break;
830 }
831 bsendsize += n;
832 nsend++;
833 }
834 if (break_flag) break;
835 }
836
837 if (bsendsize > MAXSMALLINT)
838 error->one(FLERR,"Acquire ghosts near send buffer exceeds 2 GB");
839 int sendsize = bsendsize;
840
841 // create send buf and auxiliary irregular comm vectors
842
843 char *sbuf;
844 memory->create(sbuf,sendsize,"grid:sbuf");
845 memset(sbuf,0,sendsize);
846
847 int *proclist,*sizelist;
848 memory->create(proclist,nsend,"grid:proclist");
849 memory->create(sizelist,nsend,"grid:sizelist");
850
851 // on 2nd pass over local cells, fill the send buf
852 // use lastproc to insure a cell only overlaps once per other proc
853 // if oflag = 2 = my cell just touches box,
854 // so flag grid cell as EMPTY ghost by setting nsurf = -1
855
856 nsend = 0;
857 sendsize = 0;
858 for (int icell = icell_start; icell < icell_end; icell++) {
859 if (cells[icell].nsplit <= 0) continue;
860 lo = cells[icell].lo;
861 hi = cells[icell].hi;
862 lastproc = -1;
863 int i_first = 0;
864 if (icell == icell_start) i_first = i_start;
865 int i_last = nlist;
866 if (icell == icell_end-1) i_last = i_end;
867 for (i = i_first; i < i_last; i++) {
868 j = list[i];
869 oflag = box_overlap(lo,hi,boxall[j].lo,boxall[j].hi);
870 if (!oflag) continue;
871 if (boxall[j].proc == lastproc) continue;
872 lastproc = boxall[j].proc;
873
874 if (oflag == 2) {
875 nsurf_hold = cells[icell].nsurf;
876 cells[icell].nsurf = -1;
877 }
878 sizelist[nsend] = pack_one(icell,&sbuf[sendsize],0,0,surfflag,1);
879 if (oflag == 2) cells[icell].nsurf = nsurf_hold;
880 proclist[nsend] = lastproc;
881 sendsize += sizelist[nsend];
882 nsend++;
883 }
884 }
885
886 // perform irregular communication of list of ghost cells
887
888 Irregular *irregular = new Irregular(sparta);
889 int recvsize;
890 int nrecv = irregular->create_data_variable(nsend,proclist,sizelist,
891 recvsize,comm->commsortflag);
892
893 char *rbuf;
894 memory->create(rbuf,recvsize,"grid:rbuf");
895 memset(rbuf,0,recvsize);
896
897 irregular->exchange_variable(sbuf,sizelist,rbuf);
898 delete irregular;
899 irregular = NULL;
900
901 // unpack received grid cells as ghost cells
902
903 int offset = 0;
904 for (i = 0; i < nrecv; i++)
905 offset += grid->unpack_one(&rbuf[offset],0,0,surfflag);
906
907 // more clean up
908
909 memory->destroy(proclist);
910 memory->destroy(sizelist);
911 memory->destroy(sbuf);
912 memory->destroy(rbuf);
913
914 icell_start = icell_end-1;
915 i_start = i_end;
916 int not_done_local = (icell_end < nlocal || i_end < nlist);
917 if (not_done_local && i_end == nlist) i_start = 0;
918 MPI_Allreduce(¬_done_local,¬_done,1,MPI_INT,MPI_SUM,world);
919
920 } // end while loop
921
922 // clean up
923
924 memory->destroy(list);
925 delete [] boxall;
926
927 // set nempty = # of EMPTY ghost cells I store
928
929 nempty = 0;
930 for (int icell = nlocal; icell < nlocal+nghost; icell++)
931 if (cells[icell].nsurf < 0) nempty++;
932 }
933
934 /* ----------------------------------------------------------------------
935 check for overlap of 2 orthongal boxes, alo/ahi and blo/bhi
936 if no overlap, return 0
937 if non-grazing overlap, return 1
938 if grazing overlap, return 2
939 works for 2d and 3d, since in 2d both boxes have zhi > zlo
940 ------------------------------------------------------------------------- */
941
box_intersect(double * alo,double * ahi,double * blo,double * bhi,double * lo,double * hi)942 void Grid::box_intersect(double *alo, double *ahi, double *blo, double *bhi,
943 double *lo, double *hi)
944 {
945 lo[0] = MAX(alo[0],blo[0]);
946 hi[0] = MIN(ahi[0],bhi[0]);
947 lo[1] = MAX(alo[1],blo[1]);
948 hi[1] = MIN(ahi[1],bhi[1]);
949 lo[2] = MAX(alo[2],blo[2]);
950 hi[2] = MIN(ahi[2],bhi[2]);
951 }
952
953 /* ----------------------------------------------------------------------
954 check for overlap of 2 orthongal boxes, alo/ahi and blo/bhi
955 if no overlap, return 0
956 if non-grazing overlap, return 1
957 if grazing overlap, return 2
958 works for 2d and 3d, since in 2d both boxes have zhi > zlo
959 ------------------------------------------------------------------------- */
960
box_overlap(double * alo,double * ahi,double * blo,double * bhi)961 int Grid::box_overlap(double *alo, double *ahi, double *blo, double *bhi)
962 {
963 double lo[3],hi[3];
964 box_intersect(alo,ahi,blo,bhi,lo,hi);
965
966 if (lo[0] > hi[0]) return 0;
967 if (lo[1] > hi[1]) return 0;
968 if (lo[2] > hi[2]) return 0;
969
970 if (lo[0] == hi[0]) return 2;
971 if (lo[1] == hi[1]) return 2;
972 if (lo[2] == hi[2]) return 2;
973
974 return 1;
975 }
976
977 /* ----------------------------------------------------------------------
978 split lo/hi box into multilple boxes straddling periodic boundaries
979 return # of split boxes and list of new boxes in box
980 ------------------------------------------------------------------------- */
981
box_periodic(double * lo,double * hi,Box * box)982 int Grid::box_periodic(double *lo, double *hi, Box *box)
983 {
984 int ilo,ihi,jlo,jhi,klo,khi;
985 ilo = ihi = jlo = jhi = klo = khi = 0;
986
987 int *bflag = domain->bflag;
988 if (bflag[0] == PERIODIC) {
989 ilo = -1; ihi = 1;
990 }
991 if (bflag[2] == PERIODIC) {
992 jlo = -1; jhi = 1;
993 }
994 if (bflag[4] == PERIODIC && domain->dimension == 3) {
995 klo = -1; khi = 1;
996 }
997
998 double *boxlo = domain->boxlo;
999 double *boxhi = domain->boxhi;
1000 double *prd = domain->prd;
1001
1002 double plo[3],phi[3],olo[3],ohi[3];
1003
1004 int n = 0;
1005
1006 int i,j,k;
1007 for (k = klo; k <= khi; k++)
1008 for (j = jlo; j <= jhi; j++)
1009 for (i = ilo; i <= ihi; i++) {
1010 plo[0] = lo[0] + i*prd[0];
1011 phi[0] = hi[0] + i*prd[0];
1012 plo[1] = lo[1] + j*prd[1];
1013 phi[1] = hi[1] + j*prd[1];
1014 plo[2] = lo[2] + k*prd[2];
1015 phi[2] = hi[2] + k*prd[2];
1016 box_intersect(plo,phi,boxlo,boxhi,olo,ohi);
1017 if (olo[0] >= ohi[0] || olo[1] >= ohi[1] || olo[2] >= ohi[2]) continue;
1018
1019 box[n].proc = me;
1020 box[n].lo[0] = olo[0];
1021 box[n].lo[1] = olo[1];
1022 box[n].lo[2] = olo[2];
1023 box[n].hi[0] = ohi[0];
1024 box[n].hi[1] = ohi[1];
1025 box[n].hi[2] = ohi[2];
1026 n++;
1027 }
1028
1029 return n;
1030 }
1031
1032 /* ----------------------------------------------------------------------
1033 hash all my child IDs, owned + ghost
1034 ------------------------------------------------------------------------- */
1035
rehash()1036 void Grid::rehash()
1037 {
1038 // hash all owned/ghost child cell IDs
1039 // key = ID, value = index for child cells
1040 // skip sub cells
1041
1042 hash->clear();
1043
1044 for (int icell = 0; icell < nlocal+nghost; icell++) {
1045 if (cells[icell].nsplit <= 0) continue;
1046 (*hash)[cells[icell].id] = icell;
1047 }
1048
1049 hashfilled = 1;
1050 }
1051
1052 /* ----------------------------------------------------------------------
1053 find and set neighbor indices for all owned and ghost cells
1054 when done, hash table is valid for all owned and ghost cells
1055 no-op if ghosts don't exist
1056 ------------------------------------------------------------------------- */
1057
find_neighbors()1058 void Grid::find_neighbors()
1059 {
1060 int icell,idim,iface,ilevel,level,nmask,boundary,periodic;
1061 int found,face_touching,unknownflag;
1062 cellint id,neighID,refineID,coarsenID;
1063 cellint *neigh;
1064 double *lo,*hi;
1065
1066 if (!exist_ghost) return;
1067
1068 int dimension = domain->dimension;
1069 int *bflag = domain->bflag;
1070 double *boxlo = domain->boxlo;
1071 double *boxhi = domain->boxhi;
1072
1073 // insure cell IDs (owned + ghost) are hashed
1074
1075 rehash();
1076
1077 // clear parent cells data structure since will (re)build it here
1078
1079 nparent = 0;
1080
1081 // set neigh flags and nmask for each owned and ghost child cell
1082 // skip sub cells, loop below copies their info from original split cell
1083
1084 int nunknown = 0;
1085
1086 for (icell = 0; icell < nlocal+nghost; icell++) {
1087 if (cells[icell].nsplit <= 0) continue;
1088
1089 id = cells[icell].id;
1090 level = cells[icell].level;
1091 lo = cells[icell].lo;
1092 hi = cells[icell].hi;
1093 neigh = cells[icell].neigh;
1094 nmask = 0;
1095 unknownflag = 0;
1096
1097 for (iface = 0; iface < 6; iface++) {
1098 idim = iface/2;
1099
1100 // set boundary and periodic flags for this face
1101 // treat 2d Z boundaries as non-periodic
1102
1103 if (iface % 2 == 0 && lo[idim] == boxlo[idim]) boundary = 1;
1104 else if (iface % 2 == 1 && hi[idim] == boxhi[idim]) boundary = 1;
1105 else boundary = 0;
1106 if (bflag[iface] == PERIODIC) periodic = 1;
1107 else periodic = 0;
1108 if (dimension == 2 && (iface == ZLO || iface == ZHI))
1109 periodic = 0;
1110
1111 // face = non-periodic boundary, neighbor is BOUND
1112
1113 if (boundary && !periodic) {
1114 neigh[iface] = 0;
1115 nmask = neigh_encode(NBOUND,nmask,iface);
1116 continue;
1117 }
1118
1119 // neighID = ID of neighbor cell at same level as icell
1120
1121 neighID = id_neigh_same_parent(id,level,iface);
1122 if (neighID == 0) neighID = id_neigh_same_level(id,level,iface);
1123
1124 // if in hash, neighbor is CHILD
1125
1126 if (hash->find(neighID) != hash->end()) {
1127 neigh[iface] = (*hash)[neighID];
1128 if (!boundary) nmask = neigh_encode(NCHILD,nmask,iface);
1129 else nmask = neigh_encode(NPBCHILD,nmask,iface);
1130 continue;
1131 }
1132
1133 // refine from neighID until reach maxlevel
1134 // look for a child cell on touching face I do own (or ghost)
1135 // if find one, neighbor is PARENT, add its ID and lo/hi to local pcells
1136
1137 face_touching = iface % 2 ? iface-1 : iface+1;
1138 refineID = neighID;
1139 ilevel = level;
1140 found = 0;
1141
1142 while (ilevel < maxlevel) {
1143 refineID = id_refine(refineID,ilevel,face_touching);
1144 if (hash->find(refineID) != hash->end()) {
1145 neigh[iface] = nparent;
1146 if (!boundary) nmask = neigh_encode(NPARENT,nmask,iface);
1147 else nmask = neigh_encode(NPBPARENT,nmask,iface);
1148
1149 if (nparent == maxparent) grow_pcells();
1150 pcells[nparent].id = neighID;
1151 id_lohi(neighID,level,boxlo,boxhi,pcells[nparent].lo,pcells[nparent].hi);
1152 nparent++;
1153
1154 found = 1;
1155 break;
1156 } else ilevel++;
1157 }
1158 if (found) continue;
1159
1160 // coarsen from neighID until reach top level
1161 // if find one, neighbor is CHILD
1162
1163 coarsenID = neighID;
1164 ilevel = level;
1165 found = 0;
1166
1167 while (ilevel > 1) {
1168 coarsenID = id_coarsen(coarsenID,ilevel);
1169 if (hash->find(coarsenID) != hash->end()) {
1170 neigh[iface] = (*hash)[coarsenID];
1171 if (!boundary) nmask = neigh_encode(NCHILD,nmask,iface);
1172 else nmask = neigh_encode(NPBCHILD,nmask,iface);
1173 found = 1;
1174 break;
1175 } else ilevel--;
1176 }
1177 if (found) continue;
1178
1179 // found nothing, so UNKNOWN neighbor
1180 // should never happend for an owned cell (error check below)
1181 // can happen for a ghost cell
1182
1183 neigh[iface] = 0;
1184 if (!boundary) nmask = neigh_encode(NUNKNOWN,nmask,iface);
1185 else nmask = neigh_encode(NPBUNKNOWN,nmask,iface);
1186
1187 if (icell < nlocal) unknownflag = 1;
1188 }
1189
1190 cells[icell].nmask = nmask;
1191 if (unknownflag) nunknown++;
1192 }
1193
1194 // for sub cells, copy neighbor info from original split cell
1195
1196 int m,splitcell;
1197
1198 for (icell = 0; icell < nlocal+nghost; icell++) {
1199 if (cells[icell].nsplit >= 1) continue;
1200 splitcell = sinfo[cells[icell].isplit].icell;
1201 for (m = 0; m < 6; m++)
1202 cells[icell].neigh[m] = cells[splitcell].neigh[m];
1203 cells[icell].nmask = cells[splitcell].nmask;
1204 }
1205
1206 // error if any UNKNOWN neighbors for an owned cell
1207 // cannot move particle to new proc to continue move
1208
1209 int nall;
1210 MPI_Allreduce(&nunknown,&nall,1,MPI_INT,MPI_SUM,world);
1211
1212 if (nall) {
1213 char str[128];
1214 sprintf(str,"Owned cells with unknown neighbors = %d",nall);
1215 error->all(FLERR,str);
1216 }
1217 }
1218
1219 /* ----------------------------------------------------------------------
1220 change neigh[] in owned cells that are local indices to global cell IDs
1221 nmask is not changed
1222 called before cells data structure changes
1223 e.g. due to cell migration or new surfs inducing split cells
1224 can later use reset_neighbors() instead of find_neighbors() to
1225 change neigh[] back to local indices
1226 no-op if ghosts don't exist
1227 ------------------------------------------------------------------------- */
1228
unset_neighbors()1229 void Grid::unset_neighbors()
1230 {
1231 if (!exist_ghost) return;
1232
1233 // no change in neigh[] needed if nflag = NUNKNOWN, NPBUNKNOWN, or NBOUND
1234
1235 int i,index,nmask,nflag;
1236 cellint *neigh;
1237
1238 for (int icell = 0; icell < nlocal; icell++) {
1239 neigh = cells[icell].neigh;
1240 nmask = cells[icell].nmask;
1241
1242 for (i = 0; i < 6; i++) {
1243 index = neigh[i];
1244 nflag = neigh_decode(nmask,i);
1245 if (nflag == NCHILD || nflag == NPBCHILD)
1246 neigh[i] = cells[index].id;
1247 else if (nflag == NPARENT || nflag == NPBPARENT)
1248 neigh[i] = pcells[neigh[i]].id;
1249 }
1250 }
1251 }
1252
1253 /* ----------------------------------------------------------------------
1254 reset neigh[] and nmask for owned and ghost cells after call to unset()
1255 neigh[] currently stores cell IDs
1256 convert them to local indices using hash table
1257 when done, hash table is valid for all owned and ghost cells
1258 no-op if ghosts don't exist
1259 ------------------------------------------------------------------------- */
1260
reset_neighbors()1261 void Grid::reset_neighbors()
1262 {
1263 if (!exist_ghost) return;
1264
1265 // insure all cell IDs (owned + ghost) are hashed
1266
1267 rehash();
1268
1269 // clear parent cells data structure since will (re)build it here
1270
1271 nparent = 0;
1272
1273 // set neigh[] and nmask of each owned and ghost child cell
1274 // hash lookup can reset nmask to CHILD or UNKNOWN
1275 // no change in neigh[] or nmask needed if nflag = NBOUND
1276
1277 double *boxlo = domain->boxlo;
1278 double *boxhi = domain->boxhi;
1279
1280 int i,level,nmask,nflag;
1281 cellint *neigh;
1282
1283 for (int icell = 0; icell < nlocal+nghost; icell++) {
1284 level = cells[icell].level;
1285 neigh = cells[icell].neigh;
1286 nmask = cells[icell].nmask;
1287
1288 for (i = 0; i < 6; i++) {
1289 nflag = neigh_decode(nmask,i);
1290 if (nflag == NCHILD || nflag == NPBCHILD) {
1291 if (hash->find(neigh[i]) == hash->end()) {
1292 if (nflag == NCHILD) nmask = neigh_encode(NUNKNOWN,nmask,i);
1293 else nmask = neigh_encode(NPBUNKNOWN,nmask,i);
1294 } else neigh[i] = (*hash)[neigh[i]];
1295
1296 } else if (nflag == NPARENT || nflag == NPBPARENT) {
1297 if (nparent == maxparent) grow_pcells();
1298 pcells[nparent].id = neigh[i];
1299 id_lohi(neigh[i],level,boxlo,boxhi,pcells[nparent].lo,pcells[nparent].hi);
1300 neigh[i] = nparent;
1301 nparent++;
1302
1303 } else if (nflag == NUNKNOWN || nflag == NPBUNKNOWN) {
1304 if (hash->find(neigh[i]) != hash->end()) {
1305 neigh[i] = (*hash)[neigh[i]];
1306 if (nflag == NUNKNOWN) nmask = neigh_encode(NCHILD,nmask,i);
1307 else nmask = neigh_encode(NPBCHILD,nmask,i);
1308 }
1309 }
1310 }
1311
1312 cells[icell].nmask = nmask;
1313 }
1314 }
1315
1316 /* ----------------------------------------------------------------------
1317 set type = INSIDE/OUTSIDE for all owned cells if not already set
1318 input:
1319 all cells marked as OVERLAP or UNKNOWN
1320 most OVERLAP cells have their corner points marked as INSIDE/OUTSIDE
1321 a corner point is never marked as OVERLAP
1322 some OVERLAP cells may have all corner points UNKNOWN
1323 due to Cut2d/Cut3d split() failing to mark them
1324 b/c surf(s) only touch cell at point(s) on cell surface
1325 and not enough info to determine if cell otherwise inside vs outside
1326 for these cells, volume = 0.0
1327 output:
1328 all UNKNOWN cells marked as INSIDE/OUTSIDE
1329 mark corner pts of any OVERLAP cells with UNKNOWN corner pts
1330 to be all INSIDE or all OUTSIDE, set cell volume accordingly
1331 output status will be checked when type_check() is invoked
1332 errors or warnings will be triggered
1333 NOTE: if UNKNOWN corner pts of OVERLAP cells are not set,
1334 may just be flagged as warning in type_check(),
1335 but it means a non-zero volume for a cell that
1336 is effectively OUTSIDE will not be set correctly
1337 ------------------------------------------------------------------------- */
1338
set_inout()1339 void Grid::set_inout()
1340 {
1341 int i,j,m,icell,jcell,itype,jtype,marktype;
1342 int nflag,iface,icorner,ctype,ic;
1343 int iset,nset,nsetnew;
1344 int nsend,maxsend,nrecv,maxrecv;
1345 int *set,*setnew,*jcorner;
1346 int *cflags;
1347 int *proclist;
1348 double xcorner[3];
1349 ParentCell *pcell;
1350 Connect *sbuf,*rbuf;
1351
1352 if (!exist_ghost)
1353 error->all(FLERR,"Cannot mark grid cells as inside/outside surfs because "
1354 "ghost cells do not exist");
1355
1356 // if all surfs are transparent, just mark all cells as OUTSIDE
1357
1358 if (surf->all_transparent()) {
1359 for (icell = 0; icell < nlocal; icell++) cinfo[icell].type = OUTSIDE;
1360 return;
1361 }
1362
1363 // set dimensional dependent quantities
1364
1365 int faceflip[6] = {XHI,XLO,YHI,YLO,ZHI,ZLO};
1366
1367 int me = comm->me;
1368 int dimension = domain->dimension;
1369 int ncorner,nface,nface_pts;
1370 if (dimension == 3) {
1371 ncorner = 8;
1372 nface = 6;
1373 nface_pts = 4;
1374 } else {
1375 ncorner = 4;
1376 nface = 4;
1377 nface_pts = 2;
1378 }
1379
1380 // create irregular communicator for exchanging off-processor cell types
1381
1382 Irregular *irregular = new Irregular(sparta);
1383
1384 // create set1 and set2 lists so can swap between them
1385
1386 int *set1,*set2;
1387 memory->create(set1,nlocal*nface,"grid:set1");
1388 memory->create(set2,nlocal*nface,"grid:set2");
1389
1390 // initial set list = overlapped cells with corner values which are set
1391
1392 nset = nsetnew = 0;
1393 set = set1;
1394 setnew = set2;
1395
1396 for (icell = 0; icell < nlocal; icell++) {
1397 if (cells[icell].nsplit <= 0) continue;
1398 if (cinfo[icell].type == OVERLAP && cinfo[icell].corner[0] != UNKNOWN)
1399 set[nset++] = icell;
1400 }
1401
1402 // loop until no more cell types are set
1403
1404 maxsend = maxrecv = 0;
1405 proclist = NULL;
1406 sbuf = rbuf = NULL;
1407
1408 while (1) {
1409 nsend = 0;
1410
1411 // process local set list
1412 // for unmarked neighbor cells I own, mark them and add to new set list
1413 // for neighbor cells I would mark but don't own, add to comm list
1414
1415 while (nset) {
1416 nsetnew = 0;
1417 for (iset = 0; iset < nset; iset++) {
1418 icell = set[iset];
1419 itype = cinfo[icell].type;
1420 cflags = cinfo[icell].corner;
1421
1422 // loop over my cell's faces
1423
1424 for (iface = 0; iface < nface; iface++) {
1425
1426 // if I am an OUTSIDE/INSIDE cell: marktype = my itype
1427 // if I am an OVERLAP cell:
1428 // can only mark neighbor if all corner pts on face are same
1429 // marktype = value of those corner pts = OUTSIDE or INSIDE
1430
1431 if (itype != OVERLAP) marktype = itype;
1432 else {
1433 ctype = cflags[corners[iface][0]];
1434 for (m = 1; m < nface_pts; m++)
1435 if (cflags[corners[iface][m]] != ctype) break;
1436 if (m < nface_pts) continue;
1437 if (ctype == OUTSIDE) marktype = OUTSIDE;
1438 else if (ctype == INSIDE) marktype = INSIDE;
1439 else continue;
1440 }
1441
1442 jcell = cells[icell].neigh[iface];
1443 nflag = neigh_decode(cells[icell].nmask,iface);
1444
1445 if (nflag == NCHILD || nflag == NPBCHILD) {
1446
1447 // this proc owns neighbor cell
1448 // perform this marking logic (here and below):
1449 // (1) if jtype = UNKNOWN:
1450 // set jtype = marktype = OUTSIDE/INSIDE
1451 // add jcell to new set
1452 // (2) else if jtype = OVERLAP:
1453 // skip if jcell's corners are already marked
1454 // set jcell corner pts = marktype = OUTSIDE/INSIDE
1455 // also set its volume to full cell or zero
1456 // do not add jcell to new set, since is OVERLAP cell
1457 // maybe could add, but not sure need to,
1458 // and could be marked wrong ?
1459 // (3) else jcell = OUTSIDE/INSIDE:
1460 // if jcell mark is different than marktype,
1461 // error b/c markings are inconsistent
1462
1463 if (cells[jcell].proc == me) {
1464 jtype = cinfo[jcell].type;
1465
1466 if (jtype == UNKNOWN) {
1467 cinfo[jcell].type = marktype;
1468 setnew[nsetnew++] = jcell;
1469 } else if (jtype == OVERLAP) {
1470 // don't think this restriction needed (Mar 2017)
1471 //if (itype == OVERLAP) continue;
1472 jcorner = cinfo[jcell].corner;
1473 if (jcorner[0] != UNKNOWN) continue;
1474 for (icorner = 0; icorner < ncorner; icorner++)
1475 jcorner[icorner] = marktype;
1476 if (marktype == INSIDE) cinfo[jcell].volume = 0.0;
1477 else if (marktype == OUTSIDE) {
1478 double *lo = cells[jcell].lo;
1479 double *hi = cells[jcell].hi;
1480 if (dimension == 3)
1481 cinfo[jcell].volume =
1482 (hi[0]-lo[0]) * (hi[1]-lo[1]) * (hi[2]-lo[2]);
1483 else if (domain->axisymmetric)
1484 cinfo[jcell].volume =
1485 MY_PI * (hi[1]*hi[1]-lo[1]*lo[1]) * (hi[0]-lo[0]);
1486 else
1487 cinfo[jcell].volume = (hi[0]-lo[0]) * (hi[1]-lo[1]);
1488 }
1489 // do not add to setnew, see comment above
1490 //setnew[nsetnew++] = jcell;
1491 } else {
1492 if (jtype != marktype) {
1493 printf("CELL1 proc %d icell %d id " CELLINT_FORMAT " iface %d "
1494 "jcell %d id " CELLINT_FORMAT " marktype %d jtype %d\n",
1495 me,icell,cells[icell].id,iface,jcell,cells[jcell].id,
1496 marktype,jtype);
1497 error->one(FLERR,"Cell type mis-match when marking on self");
1498 }
1499 }
1500
1501 // this proc does not own neighbor cell
1502 // pack sbuf with info to send
1503
1504 } else {
1505 if (nsend == maxsend) {
1506 maxsend += DELTA;
1507 memory->grow(proclist,maxsend,"grid:proclist");
1508 sbuf = (Connect *)
1509 memory->srealloc(sbuf,maxsend*sizeof(Connect),"grid:sbuf");
1510 }
1511 proclist[nsend] = cells[jcell].proc;
1512 sbuf[nsend].itype = itype;
1513 sbuf[nsend].marktype = marktype;
1514 sbuf[nsend].jcell = cells[jcell].ilocal;
1515 nsend++;
1516 }
1517
1518 } else if (nflag == NPARENT || nflag == NPBPARENT) {
1519
1520 // neighbor is a parent cell, not a child cell
1521 // for each corner points of icell face:
1522 // find jcell = child cell owner of the face corner pt
1523 // if I own the child cell, mark it in same manner as above
1524
1525 pcell = &pcells[jcell];
1526 for (m = 0; m < nface_pts; m++) {
1527 ic = corners[iface][m];
1528 if (ic % 2) xcorner[0] = cells[icell].hi[0];
1529 else xcorner[0] = cells[icell].lo[0];
1530 if ((ic/2) % 2) xcorner[1] = cells[icell].hi[1];
1531 else xcorner[1] = cells[icell].lo[1];
1532 if (ic/4) xcorner[2] = cells[icell].hi[2];
1533 else xcorner[2] = cells[icell].lo[2];
1534
1535 if (nflag == NPBPARENT)
1536 domain->uncollide(faceflip[iface],xcorner);
1537
1538 jcell = id_find_child(pcell->id,cells[icell].level,
1539 pcell->lo,pcell->hi,xcorner);
1540 if (jcell < 0) error->one(FLERR,"Parent cell child missing");
1541
1542 // this proc owns neighbor cell
1543 // perform same marking logic as above
1544
1545 if (cells[jcell].proc == me) {
1546 jtype = cinfo[jcell].type;
1547
1548 if (jtype == UNKNOWN) {
1549 cinfo[jcell].type = marktype;
1550 setnew[nsetnew++] = jcell;
1551 } else if (jtype == OVERLAP) {
1552 // don't think this restriction needed (Mar 2017)
1553 //if (itype == OVERLAP) continue;
1554 jcorner = cinfo[jcell].corner;
1555 if (jcorner[0] != UNKNOWN) continue;
1556 for (icorner = 0; icorner < ncorner; icorner++)
1557 jcorner[icorner] = marktype;
1558 if (marktype == INSIDE) cinfo[jcell].volume = 0.0;
1559 else if (marktype == OUTSIDE) {
1560 double *lo = cells[jcell].lo;
1561 double *hi = cells[jcell].hi;
1562 if (dimension == 3)
1563 cinfo[jcell].volume =
1564 (hi[0]-lo[0]) * (hi[1]-lo[1]) * (hi[2]-lo[2]);
1565 else if (domain->axisymmetric)
1566 cinfo[jcell].volume =
1567 MY_PI * (hi[1]*hi[1]-lo[1]*lo[1]) * (hi[0]-lo[0]);
1568 else
1569 cinfo[jcell].volume = (hi[0]-lo[0]) * (hi[1]-lo[1]);
1570 }
1571 // do not add to setnew, see comment above
1572 //setnew[nsetnew++] = jcell;
1573 } else {
1574 if (jtype != marktype) {
1575 printf("CELL2 proc %d icell %d id " CELLINT_FORMAT
1576 " iface %d "
1577 "jcell %d id " CELLINT_FORMAT " marktype %d "
1578 "jtype %d\n",
1579 me,icell,cells[icell].id,iface,jcell,cells[jcell].id,
1580 marktype,jtype);
1581 error->one(FLERR,
1582 "Cell type mis-match when marking on self");
1583 }
1584 }
1585
1586 // this proc does not own neighbor cell
1587 // pack sbuf with info to send
1588
1589 } else {
1590 if (nsend == maxsend) {
1591 maxsend += DELTA;
1592 memory->grow(proclist,maxsend,"grid:proclist");
1593 sbuf = (Connect *)
1594 memory->srealloc(sbuf,maxsend*sizeof(Connect),"grid:sbuf");
1595 }
1596 proclist[nsend] = cells[jcell].proc;
1597 sbuf[nsend].itype = itype;
1598 sbuf[nsend].marktype = marktype;
1599 sbuf[nsend].jcell = cells[jcell].ilocal;
1600 nsend++;
1601 }
1602 }
1603 } else continue;
1604 }
1605 }
1606
1607 // swap set lists
1608
1609 nset = nsetnew;
1610 if (set == set1) {
1611 set = set2;
1612 setnew = set1;
1613 } else {
1614 set = set1;
1615 setnew = set2;
1616 }
1617 }
1618
1619 // if no proc has info to communicate, then done iterating
1620
1621 int anysend;
1622 MPI_Allreduce(&nsend,&anysend,1,MPI_INT,MPI_MAX,world);
1623 if (!anysend) break;
1624
1625 // perform irregular comm of each proc's comm list
1626 // realloc rbuf as needed
1627
1628 nrecv = irregular->create_data_uniform(nsend,proclist,comm->commsortflag);
1629 if (nrecv > maxrecv) {
1630 memory->sfree(rbuf);
1631 maxrecv = nrecv;
1632 rbuf = (Connect *) memory->smalloc(maxrecv*sizeof(Connect),"grid:rbuf");
1633 }
1634
1635 irregular->exchange_uniform((char *) sbuf,sizeof(Connect),(char *) rbuf);
1636
1637 // this proc received info to mark neighbor cell it owns
1638 // perform same marking logic as above
1639
1640 for (i = 0; i < nrecv; i++) {
1641 itype = rbuf[i].itype;
1642 marktype = rbuf[i].marktype;
1643 jcell = rbuf[i].jcell;
1644 jtype = cinfo[jcell].type;
1645
1646 if (jtype == UNKNOWN) {
1647 cinfo[jcell].type = marktype;
1648 set[nset++] = jcell;
1649 } else if (jtype == OVERLAP) {
1650 // don't think this restriction needed (Mar 2017)
1651 //if (itype == OVERLAP) continue;
1652 jcorner = cinfo[jcell].corner;
1653 if (jcorner[0] != UNKNOWN) continue;
1654 for (icorner = 0; icorner < ncorner; icorner++)
1655 jcorner[icorner] = marktype;
1656 if (marktype == INSIDE) cinfo[jcell].volume = 0.0;
1657 else if (marktype == OUTSIDE) {
1658 double *lo = cells[jcell].lo;
1659 double *hi = cells[jcell].hi;
1660 if (dimension == 3)
1661 cinfo[jcell].volume =
1662 (hi[0]-lo[0]) * (hi[1]-lo[1]) * (hi[2]-lo[2]);
1663 else if (domain->axisymmetric)
1664 cinfo[jcell].volume =
1665 MY_PI * (hi[1]*hi[1]-lo[1]*lo[1]) * (hi[0]-lo[0]);
1666 else
1667 cinfo[jcell].volume = (hi[0]-lo[0]) * (hi[1]-lo[1]);
1668 }
1669 // do not add to setnew, see comment above
1670 //set[nset++] = jcell;
1671 } else {
1672 if (marktype != jtype) {
1673 printf("CELL3 me %d jcell %d id " CELLINT_FORMAT
1674 " marktype %d jtype %d\n",
1675 me,jcell,cells[jcell].id,marktype,jtype);
1676 error->one(FLERR,"Cell type mis-match when marking on neigh proc");
1677 }
1678 }
1679 }
1680 }
1681
1682 // NOTE: at this point could make a final attempt to mark
1683 // any remaining UNKNOWN corner pts of an overlap cell
1684 // to avoid warnings and errors in type_check()
1685 // when doing grid adaptation, a new cell could possibly still be unmarked,
1686 // because its already marked neighbors which are non-OVERLAP cells
1687 // will not try to mark it
1688 // this is in contrast to first-time marking, when sweep entire grid
1689 // the final attempt logic could do this:
1690 // instead of having marked cells mark their unmarked neighbors
1691 // have unmarked cells look at their neighbors to acquire markings
1692 // not necessary when doing initial full-grid sweep,
1693 // but may be necessary when doing incremental adaptation
1694
1695 // all done with marking
1696 // set type and cflags for all sub cells from split cell it belongs to
1697
1698 int splitcell;
1699
1700 for (icell = 0; icell < nlocal; icell++) {
1701 if (cells[icell].nsplit > 0) continue;
1702 splitcell = sinfo[cells[icell].isplit].icell;
1703 cinfo[icell].type = cinfo[splitcell].type;
1704 for (j = 0; j < ncorner; j++)
1705 cinfo[icell].corner[j] = cinfo[splitcell].corner[j];
1706 }
1707
1708 // set volume of cells that are now INSIDE to 0.0
1709 // this allows error check in Collide and FixGridCheck for particles
1710 // in zero-volume cells
1711
1712 for (icell = 0; icell < nlocal; icell++)
1713 if (cinfo[icell].type == INSIDE) cinfo[icell].volume = 0.0;
1714
1715 /*
1716 printf("POST INOUT %d: %d\n",comm->me,grid->nlocal);
1717 for (int i = 0; i < grid->nlocal; i++) {
1718 Grid::ChildCell *g = &grid->cells[i];
1719 if (g->id == 52)
1720 printf("ICELL %d: %d id %d pid %d lo %g %g "
1721 "hi %g %g type %d corners %d %d %d %d vol %g\n",
1722 comm->me,i,g->id,pcells[g->iparent].id,
1723 g->lo[0],
1724 g->lo[1],
1725 g->hi[0],
1726 g->hi[1],
1727 grid->cinfo[i].type,
1728 grid->cinfo[i].corner[0],
1729 grid->cinfo[i].corner[1],
1730 grid->cinfo[i].corner[2],
1731 grid->cinfo[i].corner[3],grid->cinfo[i].volume);
1732 }
1733 */
1734
1735 // clean up
1736
1737 delete irregular;
1738 memory->destroy(set1);
1739 memory->destroy(set2);
1740 memory->destroy(proclist);
1741 memory->sfree(sbuf);
1742 memory->sfree(rbuf);
1743 }
1744
1745 /* ----------------------------------------------------------------------
1746 check if the hierarchical grid is uniform
1747 uniform means all child cells are at same level
1748 thus finest level grid is effectively uniform across entire domain
1749 if uniform, set unx,uny,unz
1750 ------------------------------------------------------------------------- */
1751
check_uniform()1752 void Grid::check_uniform()
1753 {
1754 // uniform = 1 only if all child cells are at maxlevel
1755
1756 int myuniform = 1;
1757 for (int i = 0; i < nlocal; i++)
1758 if (cells[i].level != maxlevel) myuniform = 0;
1759
1760 MPI_Allreduce(&myuniform,&uniform,1,MPI_INT,MPI_MIN,world);
1761
1762 // if grid is uniform, compute grid extent in each dimension
1763
1764 if (uniform) {
1765 unx = uny = unz = 1;
1766 for (int i = 0; i < maxlevel; i++) {
1767 unx *= plevels[i].nx;
1768 uny *= plevels[i].ny;
1769 unz *= plevels[i].nz;
1770 }
1771 }
1772 }
1773
1774 /* ----------------------------------------------------------------------
1775 require all cell types be set
1776 check if all corner pts of OVERLAP cells are set
1777 error if corner pts on global box boundary are not set
1778 warning if corner pts interior to simulation box are not set
1779 check that no OUTSIDE cell has zero volume
1780 flag = 1 to output flow stats (default)
1781 ------------------------------------------------------------------------- */
1782
type_check(int outflag)1783 void Grid::type_check(int outflag)
1784 {
1785 int i;
1786
1787 // check cell types
1788
1789 int unknown = 0;
1790 for (int icell = 0; icell < nlocal; icell++) {
1791 if (cells[icell].nsplit <= 0) continue;
1792 if (cinfo[icell].type == UNKNOWN) unknown++;
1793 }
1794 int unknownall;
1795 MPI_Allreduce(&unknown,&unknownall,1,MPI_INT,MPI_SUM,world);
1796
1797 if (unknownall) {
1798 char str[128];
1799 sprintf(str,"Grid cells marked as unknown = %d",unknownall);
1800 error->all(FLERR,str);
1801 }
1802
1803 // check corner flags of cells that are OVERLAP
1804 // warn if any interior corner flags are not set
1805 // error if any corner flags on global boundaries are unset
1806
1807 double *boxlo = domain->boxlo;
1808 double *boxhi = domain->boxhi;
1809 int dimension = domain->dimension;
1810
1811 int ncorner = 4;
1812 if (dimension == 3) ncorner = 8;
1813
1814 double x[3];
1815 int inside = 0;
1816 int outside = 0;
1817
1818 for (int icell = 0; icell < nlocal; icell++) {
1819 if (cells[icell].nsplit <= 0) continue;
1820 if (cinfo[icell].type != OVERLAP) continue;
1821 for (i = 0; i < ncorner; i++) {
1822 if (cinfo[icell].corner[i] != UNKNOWN) continue;
1823 if (i % 2 == 0) x[0] = cells[icell].lo[0];
1824 else x[0] = cells[icell].hi[0];
1825 if ((i/2) % 2 == 0) x[1] = cells[icell].lo[1];
1826 else x[1] = cells[icell].hi[1];
1827 if (dimension == 3) {
1828 if (i/4 == 0) x[2] = cells[icell].lo[2];
1829 else x[2] = cells[icell].hi[2];
1830 } else x[2] = 0.0;
1831
1832 if (Geometry::point_on_hex(x,boxlo,boxhi)) {
1833 printf("BAD CORNER icell %d id %d type %d "
1834 "icorner %d x %g %g %g cflags %d %d %d %d\n",
1835 icell,cells[icell].id,cinfo[icell].type,i,x[0],x[1],x[2],
1836 cinfo[icell].corner[0],
1837 cinfo[icell].corner[1],
1838 cinfo[icell].corner[2],
1839 cinfo[icell].corner[3]);
1840 outside++;
1841 }
1842 else inside++;
1843 }
1844 }
1845
1846 int insideall;
1847 MPI_Allreduce(&inside,&insideall,1,MPI_INT,MPI_SUM,world);
1848 if (insideall) {
1849 char str[128];
1850 sprintf(str,"Grid cell interior corner points marked as unknown "
1851 "(volume will be wrong if cell is effectively outside) = %d",
1852 insideall);
1853 if (comm->me == 0) error->warning(FLERR,str);
1854 }
1855
1856 int outsideall;
1857 MPI_Allreduce(&outside,&outsideall,1,MPI_INT,MPI_SUM,world);
1858 if (outsideall) {
1859 char str[128];
1860 sprintf(str,"Grid cell corner points on boundary marked as unknown = %d",
1861 outsideall);
1862 error->all(FLERR,str);
1863 }
1864
1865 int volzero = 0;
1866 for (int icell = 0; icell < nlocal; icell++) {
1867 if (cells[icell].nsplit <= 0) continue;
1868 if (cinfo[icell].type == OUTSIDE && cinfo[icell].volume == 0.0) volzero++;
1869 }
1870 int volzeroall;
1871 MPI_Allreduce(&volzero,&volzeroall,1,MPI_INT,MPI_SUM,world);
1872 if (outsideall) {
1873 char str[128];
1874 sprintf(str,"Grid cells marked outside, but with zero volume = %d",
1875 volzeroall);
1876 error->all(FLERR,str);
1877 }
1878
1879 if (outflag) flow_stats();
1880 }
1881
1882 /* ----------------------------------------------------------------------
1883 set static cellwise fnum weights based on uncut volumes
1884 for volume, use volume of cell, whether axisymmetric or not
1885 for radius, use radius of cell centroid from axisymmetric axis
1886 weight() called from input script and read_restart (with narg = -1)
1887 weight_one() is called only for adapted cells from grid_adapt
1888 ------------------------------------------------------------------------- */
1889
weight(int narg,char ** arg)1890 void Grid::weight(int narg, char **arg)
1891 {
1892 int i;
1893
1894 if (!exist) error->all(FLERR,"Cannot weight cells before grid is defined");
1895 if (narg > 0 && narg != 1) error->all(FLERR,"Illegal weight command");
1896
1897 // if called from read_restart with narg = -1, cellweightflag is already set
1898
1899 if (narg == 1) {
1900 if (strcmp(arg[0],"none") == 0) cellweightflag = NOWEIGHT;
1901 else if (strcmp(arg[0],"volume") == 0) cellweightflag = VOLWEIGHT;
1902 else if (strcmp(arg[0],"radius") == 0) cellweightflag = RADWEIGHT;
1903 else if (strcmp(arg[0],"radius/only") == 0) cellweightflag = RADONLYWEIGHT;
1904 else error->all(FLERR,"Illegal weight command");
1905 }
1906
1907 if (cellweightflag == RADWEIGHT && !domain->axisymmetric)
1908 error->all(FLERR,"Cannot use weight cell radius unless axisymmetric");
1909
1910 if (cellweightflag == RADONLYWEIGHT && !domain->axisymmetric)
1911 error->all(FLERR,"Cannot use weight cell radius/only unless axisymmetric");
1912
1913 // set per-cell weights
1914
1915 for (i = 0; i < nlocal; i++) weight_one(i);
1916 }
1917
weight_one(int icell)1918 void Grid::weight_one(int icell)
1919 {
1920 double *lo,*hi;
1921
1922 int dimension = domain->dimension;
1923 int axisymmetric = domain->axisymmetric;
1924
1925 if (cellweightflag == NOWEIGHT) {
1926 cinfo[icell].weight = 1.0;
1927 } else if (cellweightflag == VOLWEIGHT) {
1928 lo = cells[icell].lo;
1929 hi = cells[icell].hi;
1930 if (dimension == 3)
1931 cinfo[icell].weight = (hi[0]-lo[0]) * (hi[1]-lo[1]) * (hi[2]-lo[2]);
1932 else if (axisymmetric)
1933 cinfo[icell].weight = MY_PI * (hi[1]*hi[1]-lo[1]*lo[1]) * (hi[0]-lo[0]);
1934 else
1935 cinfo[icell].weight = (hi[0]-lo[0]) * (hi[1]-lo[1]);
1936 } else if (cellweightflag == RADWEIGHT) {
1937 lo = cells[icell].lo;
1938 hi = cells[icell].hi;
1939 cinfo[icell].weight = 0.5*(hi[1]+lo[1]) * (hi[0]-lo[0]);
1940 } else if (cellweightflag == RADONLYWEIGHT) {
1941 lo = cells[icell].lo;
1942 hi = cells[icell].hi;
1943 cinfo[icell].weight = 0.5*(hi[1]+lo[1]);
1944 }
1945 }
1946
1947 ///////////////////////////////////////////////////////////////////////////
1948 // grow cell list data structures
1949 ///////////////////////////////////////////////////////////////////////////
1950
1951 /* ----------------------------------------------------------------------
1952 insure cells and cinfo can hold N and M new cells respectively
1953 ------------------------------------------------------------------------- */
1954
grow_cells(int n,int m)1955 void Grid::grow_cells(int n, int m)
1956 {
1957 if (nlocal+nghost+n >= maxcell) {
1958 int oldmax = maxcell;
1959 while (maxcell < nlocal+nghost+n) maxcell += DELTA;
1960 cells = (ChildCell *)
1961 memory->srealloc(cells,maxcell*sizeof(ChildCell),"grid:cells");
1962 memset(&cells[oldmax],0,(maxcell-oldmax)*sizeof(ChildCell));
1963 }
1964
1965 if (nlocal+m >= maxlocal) {
1966 int oldmax = maxlocal;
1967 while (maxlocal < nlocal+m) maxlocal += DELTA;
1968 cinfo = (ChildInfo *)
1969 memory->srealloc(cinfo,maxlocal*sizeof(ChildInfo),"grid:cinfo");
1970 memset(&cinfo[oldmax],0,(maxlocal-oldmax)*sizeof(ChildInfo));
1971 }
1972 }
1973
1974 /* ----------------------------------------------------------------------
1975 grow pcells
1976 NOTE: need this function for Kokkos
1977 ------------------------------------------------------------------------- */
1978
grow_pcells()1979 void Grid::grow_pcells()
1980 {
1981 maxparent += DELTAPARENT;
1982 pcells = (ParentCell *)
1983 memory->srealloc(pcells,maxparent*sizeof(ParentCell),"grid:pcells");
1984 }
1985
1986 /* ----------------------------------------------------------------------
1987 insure sinfo can hold N new split cells
1988 ------------------------------------------------------------------------- */
1989
grow_sinfo(int n)1990 void Grid::grow_sinfo(int n)
1991 {
1992 if (nsplitlocal+nsplitghost+n >= maxsplit) {
1993 int oldmax = maxsplit;
1994 while (maxsplit < nsplitlocal+nsplitghost+n) maxsplit += DELTA;
1995 sinfo = (SplitInfo *)
1996 memory->srealloc(sinfo,maxsplit*sizeof(SplitInfo),"grid:sinfo");
1997 memset(&sinfo[oldmax],0,(maxsplit-oldmax)*sizeof(SplitInfo));
1998 }
1999 }
2000
2001 /* ----------------------------------------------------------------------
2002 group grid command called via input script
2003 ------------------------------------------------------------------------- */
2004
group(int narg,char ** arg)2005 void Grid::group(int narg, char **arg)
2006 {
2007 int i,flag;
2008 double x[3];
2009
2010 if (narg < 3) error->all(FLERR,"Illegal group command");
2011
2012 int dimension = domain->dimension;
2013
2014 int igroup = find_group(arg[0]);
2015 if (igroup < 0) igroup = add_group(arg[0]);
2016 int bit = bitmask[igroup];
2017
2018 // style = region
2019 // add grid cell to group if in region
2020
2021 if (strcmp(arg[2],"region") == 0) {
2022 if (narg != 5) error->all(FLERR,"Illegal group command");
2023 int iregion = domain->find_region(arg[3]);
2024 if (iregion == -1) error->all(FLERR,"Group region ID does not exist");
2025 Region *region = domain->regions[iregion];
2026
2027 int rstyle;
2028 if (strcmp(arg[4],"all") == 0) rstyle = REGION_ALL;
2029 else if (strcmp(arg[4],"one") == 0) rstyle = REGION_ONE;
2030 else if (strcmp(arg[4],"center") == 0) rstyle = REGION_CENTER;
2031 else error->all(FLERR,"Illegal group command");
2032
2033 if (dimension == 2) x[2] = 0.0;
2034
2035 if (rstyle == REGION_ALL) {
2036 for (i = 0; i < nlocal; i++) {
2037 flag = 1;
2038
2039 if (dimension == 3) x[2] = cells[i].lo[2];
2040 x[0] = cells[i].lo[0];
2041 x[1] = cells[i].lo[1];
2042 if (!region->match(x)) flag = 0;
2043 x[0] = cells[i].hi[0];
2044 x[1] = cells[i].lo[1];
2045 if (!region->match(x)) flag = 0;
2046 x[0] = cells[i].lo[0];
2047 x[1] = cells[i].hi[1];
2048 if (!region->match(x)) flag = 0;
2049 x[0] = cells[i].hi[0];
2050 x[1] = cells[i].hi[1];
2051 if (!region->match(x)) flag = 0;
2052
2053 if (dimension == 3) x[2] = cells[i].hi[2];
2054 if (dimension == 3) {
2055 x[0] = cells[i].lo[0];
2056 x[1] = cells[i].lo[1];
2057 if (!region->match(x)) flag = 0;
2058 x[0] = cells[i].hi[0];
2059 x[1] = cells[i].lo[1];
2060 if (!region->match(x)) flag = 0;
2061 x[0] = cells[i].lo[0];
2062 x[1] = cells[i].hi[1];
2063 if (!region->match(x)) flag = 0;
2064 x[0] = cells[i].hi[0];
2065 x[1] = cells[i].hi[1];
2066 if (!region->match(x)) flag = 0;
2067 }
2068
2069 if (flag) cinfo[i].mask |= bit;
2070 }
2071
2072 } else if (rstyle == REGION_ONE) {
2073 for (i = 0; i < nlocal; i++) {
2074 flag = 0;
2075
2076 if (dimension == 3) x[2] = cells[i].lo[2];
2077 x[0] = cells[i].lo[0];
2078 x[1] = cells[i].lo[1];
2079 if (region->match(x)) flag = 1;
2080
2081 x[0] = cells[i].hi[0];
2082 x[1] = cells[i].lo[1];
2083 if (region->match(x)) flag = 1;
2084 x[0] = cells[i].lo[0];
2085 x[1] = cells[i].hi[1];
2086 if (region->match(x)) flag = 1;
2087 x[0] = cells[i].hi[0];
2088 x[1] = cells[i].hi[1];
2089 if (region->match(x)) flag = 1;
2090
2091 if (dimension == 3) x[2] = cells[i].hi[2];
2092 if (dimension == 3) {
2093 x[0] = cells[i].lo[0];
2094 x[1] = cells[i].lo[1];
2095 if (region->match(x)) flag = 1;
2096 x[0] = cells[i].hi[0];
2097 x[1] = cells[i].lo[1];
2098 if (region->match(x)) flag = 1;
2099 x[0] = cells[i].lo[0];
2100 x[1] = cells[i].hi[1];
2101 if (region->match(x)) flag = 1;
2102 x[0] = cells[i].hi[0];
2103 x[1] = cells[i].hi[1];
2104 if (region->match(x)) flag = 1;
2105 }
2106
2107 if (flag) cinfo[i].mask |= bit;
2108 }
2109
2110 } else if (rstyle == REGION_CENTER) {
2111 for (i = 0; i < nlocal; i++) {
2112 x[0] = 0.5 * (cells[i].lo[0] + cells[i].hi[0]);
2113 x[1] = 0.5 * (cells[i].lo[1] + cells[i].hi[1]);
2114 if (dimension == 3) x[2] = 0.5 * (cells[i].lo[2] + cells[i].hi[2]);
2115 if (region->match(x)) cinfo[i].mask |= bit;
2116 }
2117 }
2118
2119 // style = subtract
2120
2121 } else if (strcmp(arg[2],"subtract") == 0) {
2122 if (narg < 5) error->all(FLERR,"Illegal group command");
2123
2124 int length = narg-3;
2125 int *list = new int[length];
2126
2127 int jgroup;
2128 for (int iarg = 3; iarg < narg; iarg++) {
2129 jgroup = find_group(arg[iarg]);
2130 if (jgroup == -1) error->all(FLERR,"Group ID does not exist");
2131 list[iarg-3] = jgroup;
2132 }
2133
2134 // add to group if in 1st group in list
2135
2136 int otherbit = bitmask[list[0]];
2137
2138 for (i = 0; i < nlocal; i++)
2139 if (cinfo[i].mask & otherbit) cinfo[i].mask |= bit;
2140
2141 // remove grid cells if they are in any of the other groups
2142 // AND with inverse mask removes the grid cell from group
2143
2144 int inverse = inversemask[igroup];
2145
2146 for (int ilist = 1; ilist < length; ilist++) {
2147 otherbit = bitmask[list[ilist]];
2148 for (i = 0; i < nlocal; i++)
2149 if (cinfo[i].mask & otherbit) cinfo[i].mask &= inverse;
2150 }
2151
2152 delete [] list;
2153
2154 // style = union
2155
2156 } else if (strcmp(arg[2],"union") == 0) {
2157 if (narg < 4) error->all(FLERR,"Illegal group command");
2158
2159 int length = narg-3;
2160 int *list = new int[length];
2161
2162 int jgroup;
2163 for (int iarg = 3; iarg < narg; iarg++) {
2164 jgroup = find_group(arg[iarg]);
2165 if (jgroup == -1) error->all(FLERR,"Group ID does not exist");
2166 list[iarg-3] = jgroup;
2167 }
2168
2169 // add to group if in any other group in list
2170
2171 int otherbit;
2172
2173 for (int ilist = 0; ilist < length; ilist++) {
2174 otherbit = bitmask[list[ilist]];
2175 for (i = 0; i < nlocal; i++)
2176 if (cinfo[i].mask & otherbit) cinfo[i].mask |= bit;
2177 }
2178
2179 delete [] list;
2180
2181 // style = intersect
2182
2183 } else if (strcmp(arg[2],"intersect") == 0) {
2184 if (narg < 5) error->all(FLERR,"Illegal group command");
2185
2186 int length = narg-3;
2187 int *list = new int[length];
2188
2189 int jgroup;
2190 for (int iarg = 3; iarg < narg; iarg++) {
2191 jgroup = find_group(arg[iarg]);
2192 if (jgroup == -1) error->all(FLERR,"Group ID does not exist");
2193 list[iarg-3] = jgroup;
2194 }
2195
2196 // add to group if in all groups in list
2197
2198 int otherbit,ok,ilist;
2199
2200 for (i = 0; i < nlocal; i++) {
2201 ok = 1;
2202 for (ilist = 0; ilist < length; ilist++) {
2203 otherbit = bitmask[list[ilist]];
2204 if ((cinfo[i].mask & otherbit) == 0) ok = 0;
2205 }
2206 if (ok) cinfo[i].mask |= bit;
2207 }
2208
2209 delete [] list;
2210
2211 // style = clear
2212 // remove all grid cells from group
2213
2214 } else if (strcmp(arg[2],"clear") == 0) {
2215 if (igroup == 0) error->all(FLERR,"Cannot clear group all");
2216 int inversebits = inversemask[igroup];
2217
2218 for (i = 0; i < nlocal; i++) cinfo[i].mask &= inversebits;
2219 }
2220
2221 // print stats for changed group
2222
2223 bigint n = 0;
2224 for (i = 0; i < nlocal; i++)
2225 if (cinfo[i].mask & bit) n++;
2226
2227 bigint nall;
2228 MPI_Allreduce(&n,&nall,1,MPI_SPARTA_BIGINT,MPI_SUM,world);
2229
2230 if (comm->me == 0) {
2231 if (screen)
2232 fprintf(screen,BIGINT_FORMAT " grid cells in group %s\n",
2233 nall,gnames[igroup]);
2234 if (logfile)
2235 fprintf(logfile,BIGINT_FORMAT " grid cells in group %s\n",
2236 nall,gnames[igroup]);
2237 }
2238 }
2239
2240 /* ----------------------------------------------------------------------
2241 add a new grid group ID, assumed to be unique
2242 ------------------------------------------------------------------------- */
2243
add_group(const char * id)2244 int Grid::add_group(const char *id)
2245 {
2246 if (ngroup == MAXGROUP)
2247 error->all(FLERR,"Cannot have more than 32 surface groups");
2248
2249 int n = strlen(id) + 1;
2250 gnames[ngroup] = new char[n];
2251 strcpy(gnames[ngroup],id);
2252 ngroup++;
2253 return ngroup-1;
2254 }
2255
2256 /* ----------------------------------------------------------------------
2257 find a grid group ID
2258 return index of group or -1 if not found
2259 ------------------------------------------------------------------------- */
2260
find_group(const char * id)2261 int Grid::find_group(const char *id)
2262 {
2263 int igroup;
2264 for (igroup = 0; igroup < ngroup; igroup++)
2265 if (strcmp(id,gnames[igroup]) == 0) break;
2266 if (igroup == ngroup) return -1;
2267 return igroup;
2268 }
2269
2270 /* ----------------------------------------------------------------------
2271 check if a grid igroup is a uniform grid
2272 no cells with surfs
2273 all cells are at same level
2274 group forms a contiguous 3d block of cells
2275 return count of my cells in the group
2276 set nxyz = extent of 3d block in each dim
2277 set corner = lower left corner of 3d block
2278 set xyzsize = size of one grid cell (all are same size)
2279 ------------------------------------------------------------------------- */
2280
check_uniform_group(int igroup,int * nxyz,double * corner,double * xyzsize)2281 int Grid::check_uniform_group(int igroup, int *nxyz,
2282 double *corner, double *xyzsize)
2283 {
2284 double lo[3],hi[3],onesize[3];
2285
2286 int sflag = 0;
2287 int minlev = maxlevel;
2288 int maxlev = 0;
2289
2290 int count = 0;
2291 lo[0] = domain->boxhi[0];
2292 lo[1] = domain->boxhi[1];
2293 lo[2] = domain->boxhi[2];
2294 hi[0] = domain->boxlo[0];
2295 hi[1] = domain->boxlo[1];
2296 hi[2] = domain->boxlo[2];
2297
2298 int groupbit = bitmask[igroup];
2299
2300 for (int icell = 0; icell < nlocal; icell++) {
2301 if (!(cinfo[icell].mask & groupbit)) continue;
2302 if (cells[icell].nsurf) sflag++;
2303 minlev = MIN(minlev,cells[icell].level);
2304 maxlev = MAX(maxlev,cells[icell].level);
2305 lo[0] = MIN(lo[0],cells[icell].lo[0]);
2306 lo[1] = MIN(lo[1],cells[icell].lo[1]);
2307 lo[2] = MIN(lo[2],cells[icell].lo[2]);
2308 hi[0] = MAX(hi[0],cells[icell].hi[0]);
2309 hi[1] = MAX(hi[1],cells[icell].hi[1]);
2310 hi[2] = MAX(hi[2],cells[icell].hi[2]);
2311 count++;
2312 }
2313
2314 // check that no cells already have surfs
2315
2316 int allsflag;
2317 MPI_Allreduce(&sflag,&allsflag,1,MPI_INT,MPI_SUM,world);
2318 if (allsflag) {
2319 char str[128];
2320 sprintf(str,
2321 "Read_isurfs adding surfs to %d cells which already have surfs",
2322 allsflag);
2323 error->all(FLERR,str);
2324 }
2325
2326 // check that all cells are at same level
2327
2328 int allminlev,allmaxlev;
2329 MPI_Allreduce(&minlev,&allminlev,1,MPI_INT,MPI_MIN,world);
2330 MPI_Allreduce(&maxlev,&allmaxlev,1,MPI_INT,MPI_MAX,world);
2331 if (allminlev != allmaxlev)
2332 error->all(FLERR,"Read_isurfs grid group is not all at uniform level");
2333
2334 // check that cell count matches a contiguous block of cells
2335 // xyzsize = size of one cell at allmaxlev
2336
2337 xyzsize[0] = domain->xprd;
2338 xyzsize[1] = domain->yprd;
2339 xyzsize[2] = domain->zprd;
2340
2341 for (int ilevel = 0; ilevel < allmaxlev; ilevel++) {
2342 xyzsize[0] /= plevels[ilevel].nx;
2343 xyzsize[1] /= plevels[ilevel].ny;
2344 xyzsize[2] /= plevels[ilevel].nz;
2345 }
2346
2347 double alllo[3],allhi[3];
2348 MPI_Allreduce(lo,&alllo,3,MPI_DOUBLE,MPI_MIN,world);
2349 MPI_Allreduce(hi,&allhi,3,MPI_DOUBLE,MPI_MAX,world);
2350
2351 corner[0] = alllo[0];
2352 corner[1] = alllo[1];
2353 corner[2] = alllo[2];
2354
2355 nxyz[0] = static_cast<int> ((allhi[0]-alllo[0])/xyzsize[0] + 0.5);
2356 nxyz[1] = static_cast<int> ((allhi[1]-alllo[1])/xyzsize[1] + 0.5);
2357 nxyz[2] = static_cast<int> ((allhi[2]-alllo[2])/xyzsize[2] + 0.5);
2358
2359 bigint allbcount;
2360 bigint bcount = count;
2361 MPI_Allreduce(&bcount,&allbcount,1,MPI_SPARTA_BIGINT,MPI_SUM,world);
2362
2363 if ((bigint) nxyz[0] * nxyz[1]*nxyz[2] != allbcount)
2364 error->all(FLERR,"Read_isurfs grid group is not a contiguous brick");
2365
2366 return count;
2367 }
2368
2369 /* ----------------------------------------------------------------------
2370 proc 0 writes grid info and group info to restart file
2371 ------------------------------------------------------------------------- */
2372
write_restart(FILE * fp)2373 void Grid::write_restart(FILE *fp)
2374 {
2375 fwrite(&maxlevel,sizeof(int),1,fp);
2376 fwrite(plevels,sizeof(ParentLevel),maxlevel,fp);
2377
2378 fwrite(&ngroup,sizeof(int),1,fp);
2379
2380 int n;
2381 for (int i = 0; i < ngroup; i++) {
2382 n = strlen(gnames[i]) + 1;
2383 fwrite(&n,sizeof(int),1,fp);
2384 fwrite(gnames[i],sizeof(char),n,fp);
2385 }
2386 }
2387
2388 /* ----------------------------------------------------------------------
2389 proc 0 reads grid info and group info from restart file
2390 bcast to other procs
2391 ------------------------------------------------------------------------- */
2392
read_restart(FILE * fp)2393 void Grid::read_restart(FILE *fp)
2394 {
2395 // read_restart may have reset maxsurfpercell in its header() method
2396 // if so, need to reallocate surf arrays to correct max length
2397
2398 if (maxsurfpercell != MAXSURFPERCELL) allocate_surf_arrays();
2399
2400 // read level info
2401
2402 if (me == 0) {
2403 fread(&maxlevel,sizeof(int),1,fp);
2404 fread(plevels,sizeof(ParentLevel),maxlevel,fp);
2405 }
2406 MPI_Bcast(&maxlevel,1,MPI_INT,0,world);
2407 MPI_Bcast(plevels,maxlevel*sizeof(ParentLevel),MPI_CHAR,0,world);
2408
2409 // if any exist, clear existing group names, before reading new ones
2410
2411 for (int i = 0; i < ngroup; i++) delete [] gnames[i];
2412
2413 if (me == 0) fread(&ngroup,sizeof(int),1,fp);
2414 MPI_Bcast(&ngroup,1,MPI_INT,0,world);
2415
2416 int n;
2417 for (int i = 0; i < ngroup; i++) {
2418 if (me == 0) fread(&n,sizeof(int),1,fp);
2419 MPI_Bcast(&n,1,MPI_INT,0,world);
2420 gnames[i] = new char[n];
2421 if (me == 0) fread(gnames[i],sizeof(char),n,fp);
2422 MPI_Bcast(gnames[i],n,MPI_CHAR,0,world);
2423 }
2424 }
2425
2426 /* ----------------------------------------------------------------------
2427 return size of child grid restart info for this proc
2428 using count of all owned cells
2429 // NOTE: worry about N overflowing int, and in IROUNDUP ???
2430 ------------------------------------------------------------------------- */
2431
size_restart()2432 int Grid::size_restart()
2433 {
2434 int n = 2*sizeof(int);
2435 n = IROUNDUP(n);
2436 n += nlocal * sizeof(cellint);
2437 n = IROUNDUP(n);
2438 n += nlocal * sizeof(int);
2439 n = IROUNDUP(n);
2440 n += nlocal * sizeof(int);
2441 n = IROUNDUP(n);
2442 return n;
2443 }
2444
2445 /* ----------------------------------------------------------------------
2446 return size of child grid restart info
2447 using nlocal_restart count of all owned cells
2448 ------------------------------------------------------------------------- */
2449
size_restart(int nlocal_restart)2450 int Grid::size_restart(int nlocal_restart)
2451 {
2452 int n = 2*sizeof(int);
2453 n = IROUNDUP(n);
2454 n += nlocal_restart * sizeof(cellint);
2455 n = IROUNDUP(n);
2456 n += nlocal_restart * sizeof(int);
2457 n = IROUNDUP(n);
2458 n += nlocal_restart * sizeof(int);
2459 n = IROUNDUP(n);
2460 return n;
2461 }
2462
2463 /* ----------------------------------------------------------------------
2464 pack my child grid info into buf
2465 nlocal, clumped as scalars
2466 ID, level, nsplit as vectors for all owned cells
2467 // NOTE: worry about N overflowing int, and in IROUNDUP ???
2468 ------------------------------------------------------------------------- */
2469
pack_restart(char * buf)2470 int Grid::pack_restart(char *buf)
2471 {
2472 int n;
2473
2474 int *ibuf = (int *) buf;
2475 ibuf[0] = nlocal;
2476 ibuf[1] = clumped;
2477 n = 2*sizeof(int);
2478 n = IROUNDUP(n);
2479
2480 cellint *cbuf = (cellint *) &buf[n];
2481 for (int i = 0; i < nlocal; i++)
2482 cbuf[i] = cells[i].id;
2483 n += nlocal * sizeof(cellint);
2484 n = IROUNDUP(n);
2485
2486 ibuf = (int *) &buf[n];
2487 for (int i = 0; i < nlocal; i++)
2488 ibuf[i] = cells[i].level;
2489 n += nlocal * sizeof(int);
2490 n = IROUNDUP(n);
2491
2492 ibuf = (int *) &buf[n];
2493 for (int i = 0; i < nlocal; i++)
2494 ibuf[i] = cells[i].nsplit;
2495 n += nlocal * sizeof(int);
2496 n = IROUNDUP(n);
2497
2498 return n;
2499 }
2500
2501 /* ----------------------------------------------------------------------
2502 unpack child grid info into restart storage
2503 nlocal_restart, clumped as scalars
2504 id_restart, nsplit_restart as vectors
2505 allocate vectors here, will be deallocated by ReadRestart
2506 ------------------------------------------------------------------------- */
2507
unpack_restart(char * buf)2508 int Grid::unpack_restart(char *buf)
2509 {
2510 int n;
2511
2512 int *ibuf = (int *) buf;
2513 nlocal_restart = ibuf[0];
2514 clumped = ibuf[1];
2515 n = 2*sizeof(int);
2516 n = IROUNDUP(n);
2517
2518 memory->create(id_restart,nlocal_restart,"grid:id_restart");
2519 memory->create(level_restart,nlocal_restart,"grid:nlevel_restart");
2520 memory->create(nsplit_restart,nlocal_restart,"grid:nsplit_restart");
2521
2522 cellint *cbuf = (cellint *) &buf[n];
2523 for (int i = 0; i < nlocal_restart; i++)
2524 id_restart[i] = cbuf[i];
2525 n += nlocal_restart * sizeof(cellint);
2526 n = IROUNDUP(n);
2527
2528 ibuf = (int *) &buf[n];
2529 for (int i = 0; i < nlocal_restart; i++)
2530 level_restart[i] = ibuf[i];
2531 n += nlocal_restart * sizeof(int);
2532 n = IROUNDUP(n);
2533
2534 ibuf = (int *) &buf[n];
2535 for (int i = 0; i < nlocal_restart; i++)
2536 nsplit_restart[i] = ibuf[i];
2537 n += nlocal_restart * sizeof(int);
2538 n = IROUNDUP(n);
2539
2540 return n;
2541 }
2542
2543 /* ---------------------------------------------------------------------- */
2544
memory_usage()2545 bigint Grid::memory_usage()
2546 {
2547 bigint bytes = maxcell * sizeof(ChildCell);
2548 bytes += maxlocal * sizeof(ChildInfo);
2549 bytes += maxsplit * sizeof(SplitInfo);
2550 bytes += csurfs->size();
2551 bytes += csplits->size();
2552
2553 return bytes;
2554 }
2555
2556 /* ---------------------------------------------------------------------- */
2557
debug()2558 void Grid::debug()
2559 {
2560 for (int i = 0; i < nlocal; i++) {
2561 printf("GRID %d " CELLINT_FORMAT ": \n",i,cells[i].id);
2562 printf(" neigh " CELLINT_FORMAT " " CELLINT_FORMAT " "
2563 CELLINT_FORMAT " " CELLINT_FORMAT " "
2564 CELLINT_FORMAT " " CELLINT_FORMAT "\n",
2565 cells[i].neigh[0],cells[i].neigh[1],cells[i].neigh[2],
2566 cells[i].neigh[3],cells[i].neigh[4],cells[i].neigh[5]);
2567 printf(" lohi %g %g %g: %g %g %g\n",
2568 cells[i].lo[0],cells[i].lo[1],cells[i].lo[2],
2569 cells[i].hi[0],cells[i].hi[1],cells[i].hi[2]);
2570 printf(" nsurf %d:",cells[i].nsurf);
2571 for (int j = 0; j < cells[i].nsurf; j++)
2572 printf(" %d",cells[i].csurfs[j]);
2573 printf("\n");
2574 printf(" nsplit %d isplit %d\n",cells[i].nsplit,cells[i].isplit);
2575 printf(" type %d corner %d %d %d %d %d %d %d %d\n",
2576 cinfo[i].type,
2577 cinfo[i].corner[0],cinfo[i].corner[1],cinfo[i].corner[2],
2578 cinfo[i].corner[3],cinfo[i].corner[4],cinfo[i].corner[5],
2579 cinfo[i].corner[6],cinfo[i].corner[7]);
2580 printf(" volume %g\n",cinfo[i].volume);
2581 }
2582 }
2583