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 #ifndef SPARTA_GRID_H
16 #define SPARTA_GRID_H
17 
18 #include "stdio.h"
19 #include "pointers.h"
20 #include "hash3.h"
21 #include "my_page.h"
22 #include "surf.h"
23 
24 namespace SPARTA_NS {
25 
26 class Grid : protected Pointers {
27  public:
28   int exist;            // 1 if grid is defined
29   int exist_ghost;      // 1 if ghost cells exist
30   int clumped;          // 1 if grid ownership is clumped, due to RCB
31                         // if not, some operations are not allowed
32 
33   bigint ncell;         // global count of child cells (unsplit+split, no sub)
34   bigint nunsplit;      // global count of unsplit cells
35   int nsplit;           // global count of split cells
36   int nsub;             // global count of split sub cells
37 
38   int maxlevel;         // max level of any child cell in grid, 0 = root
39   int plevel_limit;     // allocation bound of plevels
40 
41   int uniform;          // 1 if all child cells are at same level, else 0
42   int unx,uny,unz;      // if uniform, effective global Nx,Ny,Nz of finest grid
43   double cutoff;        // cutoff for ghost cells, -1.0 = infinite
44   double cell_epsilon;  // half of smallest cellside of any cell in any dim
45   int cellweightflag;   // 0/1+ for no/yes usage of cellwise fnum weighting
46 
47   int surfgrid_algorithm;  // algorithm for overlap of surfs & grid cells
48   int maxsurfpercell;   // max surf elements in one child cell
49   int maxsplitpercell;  // max split cells in one child cell
50 
51   int ngroup;           // # of defined groups
52   char **gnames;        // name of each group
53   int *bitmask;         // one-bit mask for each group
54   int *inversemask;     // inverse mask for each group
55 
56   double tmap,tsplit;   // timing breakdowns of both grid2surf() algs
57   double tcomm1,tcomm2,tcomm3,tcomm4;
58 
59   int copy,copymode;    // 1 if copy of class (prevents deallocation of
60                         //  base class when child copy is destroyed)
61 
62   // cell ID hash (owned + ghost, no sub-cells)
63 
64 #ifdef SPARTA_MAP
65   typedef std::map<cellint,int> MyHash;
66 #elif SPARTA_UNORDERED_MAP
67   typedef std::unordered_map<cellint,int> MyHash;
68 #else
69   typedef std::tr1::unordered_map<cellint,int> MyHash;
70 #endif
71 
72   MyHash *hash;
73   int hashfilled;             // 1 if hash is filled with cell IDs
74 
75   // list data structs
76 
77   MyPage<surfint> *csurfs;    // lists of surf indices for
78                               // owned + ghost child cells with surfs
79   MyPage<int> *csplits;       // lists of sub cell offsets for
80                               // owned + ghost split info
81   MyPage<int> *csubs;         // lists of sub cell indices for
82                               // owned + ghost split info
83 
84   // owned or ghost child cell
85   // includes unsplit cells, split cells, sub cells in any order
86   // ghost cells are appended to owned
87 
88   struct ChildCell {
89     cellint id;               // ID of child cell
90     int level;                // level of cell in hierarchical grid, 0 = root
91     int proc;                 // proc that owns this cell
92     int ilocal;               // index of this cell on owning proc
93                               // must be correct for all ghost cells
94 
95     cellint neigh[6];         // info on 6 neighbor cells that fully overlap faces
96                               // order = XLO,XHI,YLO,YHI,ZLO,ZHI
97                               // nmask has flags for what the values represent
98                               // for nmask 0,3: index into cells (own or ghost)
99                               // for nmask 1,4: index into pcells
100                               // for nmask 2,5: neigh is unknown, value ignored
101                               // if unknown or non-PBC boundary or 2d ZLO/ZHI, ignored
102                               // must be cellint, b/c sometimes stores cell IDs
103 
104     int nmask;                // 3 bits for each of 6 values in neigh
105                               // 0 = index of child neighbor I own
106                               // 1 = index of parent neighbor in pcells
107                               // 2 = unknown child neighbor
108                               // 3 = index of PBC child neighbor I own
109                               // 4 = index of PBC parent neighbor in pcells
110                               // 5 = unknown PBC child neighbor
111                               // 6 = non-PBC boundary or ZLO/ZHI in 2d
112 
113     double lo[3],hi[3];       // opposite corner pts of cell
114 
115     int nsurf;                // # of surf elements in cell
116                               // -1 = empty ghost cell
117     surfint *csurfs;          // indices of surf elements in cell
118                               // sometimes global surf IDs are stored
119                               // for sub cells, lo/hi/nsurf/csurfs
120                               //   are same as in split cell containing them
121 
122     int nsplit;               // 1, unsplit cell
123                               // N > 1, split cell with N sub cells
124                               // N <= 0, neg of sub cell index (0 to Nsplit-1)
125     int isplit;               // index into sinfo
126                               // set for split and sub cells, -1 if unsplit
127   };
128 
129   // info specific to owned child cell
130   // includes unsplit cells, split cells, sub cells in same order as cells
131 
132   struct ChildInfo {
133     int count;                // # of particles in this cell, 0 if split cell
134     int first;                // index of 1st particle in this cell, -1 if none
135 
136     int mask;                 // grid group mask
137     int type;                 // OUTSIDE,INSIDE,OVERLAP,UNKNOWN
138     int corner[8];            // corner flags, 4/8 in 2d/3d
139                               // OUTSIDE,INSIDE,UNKNOWN
140                               // no OVERLAP is used for this anymore I think
141                               // ordered x first, y next, z last
142                               // for sub cells, type/corner
143                               //   are same as in split cell containing them
144 
145     double volume;            // flow volume of cell or sub cell
146                               // entire cell volume for split cell
147     double weight;            // fnum weighting for this cell
148   };
149 
150   // additional info for owned or ghost split cell or sub cell
151   // ghost split cell info is appended to owned split cell info
152 
153   struct SplitInfo {
154     int icell;                // index of split cell this sub cell belongs to
155     int xsub;                 // which sub cell (0 to Nsplit-1) xsplit is in
156     double xsplit[3];         // coords of point in split cell
157     int *csplits;             // sub cell (0 to Nsplit-1) each Nsurf belongs to
158     int *csubs;               // indices in cells of Nsplit sub cells
159   };
160 
161   struct ParentLevel {
162     int nbits;                // nbits = # of bits to store parent ID at this level
163     int newbits;              // newbits = extra bits to store children of this parent
164     int nx,ny,nz;             // child grid below this parent level
165     bigint nxyz;              // # of child cells of this parent level
166   };
167 
168   struct ParentCell {
169     cellint id;               // ID of parent cell
170     double lo[3],hi[3];       // opposite corner pts of cell
171   };
172 
173   int nlocal;                 // # of child cells I own (all 3 kinds)
174   int nghost;                 // # of ghost child cells I store (all 3 kinds)
175   int nempty;                 // # of empty ghost cells I store
176   int nunsplitlocal;          // # of unsplit cells I own
177   int nunsplitghost;          // # of ghost unsplit cells I store
178   int nsplitlocal;            // # of split cells I own
179   int nsplitghost;            // # of ghost split cells I store
180   int nsublocal;              // # of sub cells I own
181   int nsubghost;              // # of ghost sub cells I store
182   int nparent;                // # of parent cell neighbors I store
183 
184   int maxlocal;               // size of cinfo
185   int maxparent;              // size of pcells
186 
187   ChildCell *cells;           // list of owned and ghost child cells
188   ChildInfo *cinfo;           // extra info for nlocal owned cells
189   SplitInfo *sinfo;           // extra info for owned and ghost split cells
190 
191   ParentLevel *plevels;       // list of parent levels, level = root = simulation box
192   ParentCell *pcells;         // list of parent cell neighbors
193 
194   // restart buffers, filled by read_restart
195 
196   int nlocal_restart;
197   cellint *id_restart;
198   int *level_restart,*nsplit_restart;
199 
200   // methods
201 
202   Grid(class SPARTA *);
203   ~Grid();
204   void remove();
205   void init();
206   void add_child_cell(cellint, int, double *, double *);
207   void add_split_cell(int);
208   void add_sub_cell(int, int);
209   void notify_changed();
210   int set_minlevel();
211   void set_maxlevel();
212   void setup_owned();
213   void remove_ghosts();
214   void acquire_ghosts(int surfflag=1);
215   void rehash();
216   void find_neighbors();
217   void unset_neighbors();
218   void reset_neighbors();
219   void set_inout();
220   void check_uniform();
221   void type_check(int flag=1);
222   void weight(int, char **);
223   void weight_one(int);
224 
225   void refine_cell(int, int *, class Cut2d *, class Cut3d *);
226   void coarsen_cell(cellint, int, double *, double *,
227 		    int, int *, int *, int *, void **, char **,
228 		    class Cut2d *, class Cut3d *);
229 
230   void group(int, char **);
231   int add_group(const char *);
232   int find_group(const char *);
233   int check_uniform_group(int, int *, double *, double *);
234 
235   void write_restart(FILE *);
236   void read_restart(FILE *);
237   int size_restart();
238   int size_restart(int);
239   int pack_restart(char *);
240   int unpack_restart(char *);
241 
242   bigint memory_usage();
243 
244   void debug();
245 
246   // grid_comm.cpp
247 
248   int pack_one(int, char *, int, int, int, int);
249   int unpack_one(char *, int, int, int, int sortflag=0);
250   int pack_one_adapt(char *, char *, int);
251   int pack_particles(int, char *, int);
252   int unpack_particles(char *, int, int);
253   void unpack_particles_adapt(int, char *);
254   void compress();
255 
256   // grid_surf.cpp
257 
258   void surf2grid(int, int outflag=1);
259   void surf2grid_implicit(int, int outflag=1);
260   void surf2grid_one(int, int, int, int, class Cut3d *, class Cut2d *);
261   void clear_surf();
262   void clear_surf_restart();
263   void combine_split_cell_particles(int, int);
264   void assign_split_cell_particles(int);
265   int outside_surfs(int, double *, class Cut3d *, class Cut2d *);
266   void allocate_surf_arrays();
267   int *csubs_request(int);
268 
269   // grid_id.cpp
270 
271   void id_point_child(double *, double *, double *, int, int, int,
272 		      int &, int &, int &);
273   cellint id_parent_of_child(cellint, int);
274   int id_find_child(cellint, int, double *, double *, double *);
275   cellint id_uniform_level(int, int, int, int);
276   void id_find_child_uniform_level(int, int, double *, double *, double *,
277 				   int &, int &, int &);
278   cellint id_neigh_same_parent(cellint, int, int);
279   cellint id_neigh_same_level(cellint, int, int);
280   cellint id_refine(cellint, int, int);
281   cellint id_coarsen(cellint, int);
282   cellint id_ichild(cellint, cellint, int);
283   int id_level(cellint);
284   void id_child_lohi(int, double *, double *, cellint, double *, double *);
285   void id_lohi(cellint, int, double *, double *, double *, double *);
286   int id_bits(int, int, int);
287   void id_num2str(cellint, char *);
288 
289   // extract/return neighbor flag for iface from per-cell nmask
290   // inlined for efficiency
291 
neigh_decode(int nmask,int iface)292   inline int neigh_decode(int nmask, int iface) {
293     return (nmask & neighmask[iface]) >> neighshift[iface];
294   }
295 
296   // overwrite neighbor flag for iface in per-cell nmask
297   // first line zeroes the iface bits via one's complement of mask
298   // inlined for efficiency
299   // return updated nmask
300 
neigh_encode(int flag,int nmask,int iface)301   inline int neigh_encode(int flag, int nmask, int iface) {
302     nmask &= ~neighmask[iface];
303     nmask |= flag << neighshift[iface];
304     return nmask;
305   }
306 
307  protected:
308   int me;
309   int maxcell;             // size of cells
310   int maxsplit;            // size of sinfo
311   int maxbits;             // max bits allowed in a cell ID
312 
313   int neighmask[6];        // bit-masks for each face in nmask
314   int neighshift[6];       // bit-shifts for each face in nmask
315 
316   class Cut2d *cut2d;
317   class Cut3d *cut3d;
318 
319   // connection between one of my cells and a neighbor cell on another proc
320 
321   struct Connect {
322     int itype;           // type of sending cell
323     int marktype;        // new type value (IN/OUT) for neighbor cell
324     int jcell;           // index of neighbor cell on receiving proc (owner)
325   };
326 
327   // bounding box for a clump of grid cells
328 
329   struct Box {
330     double lo[3],hi[3];    // opposite corners of extended bbox
331     int proc;              // proc that owns it
332   };
333 
334   // tree of RCB cuts for a partitioned uniform block of grid cells
335 
336   struct GridTree {
337     int dim,cut;
338   };
339 
340   // data structs for rendezvous comm
341 
342   struct InRvous {
343     int proc;
344     surfint surfID;
345   };
346 
347   struct OutRvousLine {
348     Surf::Line line;
349   };
350 
351   struct OutRvousTri {
352     Surf::Tri tri;
353   };
354 
355   // surf ID hashes
356 
357 #ifdef SPARTA_MAP
358     typedef std::map<surfint,int> MySurfHash;
359     typedef std::map<surfint,int>::iterator MyIterator;
360 #elif SPARTA_UNORDERED_MAP
361     typedef std::unordered_map<surfint,int> MySurfHash;
362     typedef std::unordered_map<surfint,int>::iterator MyIterator;
363 #else
364     typedef std::tr1::unordered_map<surfint,int> MySurfHash;
365     typedef std::tr1::unordered_map<surfint,int>::iterator MyIterator;
366 #endif
367 
368   // Particle class values used for packing/unpacking particles in grid comm
369 
370   int ncustom;
371   int nbytes_particle,nbytes_custom,nbytes_total;
372 
373   // private methods
374 
375   void surf2grid_cell_algorithm(int);
376   void surf2grid_surf_algorithm(int);
377   void surf2grid_split(int, int);
378   void recurse2d(cellint, int, double *, double *,
379 		 int, Surf::Line *, double *, double *,
380 		 int &, int &, int **&, MyHash *, MyHash *);
381   void recurse3d(cellint, int, double *, double *,
382 		 int, Surf::Tri *, double *, double *,
383 		 int &, int &, int **&, MyHash *, MyHash *);
384   void partition_grid(int, int, int, int, int, int, int, int, GridTree *);
385   void mybox(int, int, int, int &, int &, int &, int &, int &, int &,
386 	     GridTree *);
387   void box_drop(int *, int *, int, int, GridTree *, int &, int *);
388 
389   void acquire_ghosts_all(int);
390   void acquire_ghosts_near(int);
391   void acquire_ghosts_near_less_memory(int);
392 
393   void box_intersect(double *, double *, double *, double *,
394                      double *, double *);
395   int box_overlap(double *, double *, double *, double *);
396   int box_periodic(double *, double *, Box *);
397 
398   virtual void grow_cells(int, int);
399   virtual void grow_pcells();
400   virtual void grow_sinfo(int);
401 
402   void surf2grid_stats();
403   void flow_stats();
404   double flow_volume();
405 
406   // callback function for ring communication and class variable to access
407 
408   int unpack_ghosts_surfflag;
409   static void unpack_ghosts(int, char *, void *);
410 
411   // callback functions for rendezvous communication
412 
413   static int rendezvous_surfrequest(int, char *, int &, int *&, char *&, void *);
414 };
415 
416 }
417 
418 #endif
419 
420 /* ERROR/WARNING messages:
421 
422 E: Cell ID has too many bits
423 
424 Cell IDs must fit in 32 bits (SPARTA small integer) or 64 bits (SPARTA
425 big integer), as specified by the -DSPARTA_SMALL, -DSPARTA_BIG, or
426 -DSPARTA_BIGBIG options in the machine Makefile used to build
427 SPARTA.  See Section 2.2 of the manual for details.  And see Section
428 4.8 for details on how cell IDs are formatted.
429 
430 E: Owned cells with unknown neighbors = %d
431 
432 One or more grid cells have unknown neighbors which will prevent
433 particles from moving correctly.  Please report the issue to the
434 SPARTA developers.
435 
436 E: Grid in/out self-mark error %d for icell %d, icorner %d, connect %d %d, other cell %d, other corner %d, values %d %d\n
437 
438 A grid cell was incorrectly marked as inside, outside, or overlapping
439 with surface elements.  Please report the issue to the SPARTA
440 developers.
441 
442 E: Grid in/out other-mark error %d\n
443 
444 Grid cell marking as inside, outside, or overlapping with surface
445 elements failed.  Please report the issue to the SPARTA developers.
446 
447 E: Cell type mis-match when marking on self
448 
449 Grid cell marking as inside, outside, or overlapping with surface
450 elements failed.  Please report the issue to the SPARTA developers.
451 
452 E: Parent cell child missing
453 
454 Hierarchical grid traversal failed.  Please report the issue to the
455 SPARTA developers.
456 
457 E: Cell type mis-match when marking on neigh proc
458 
459 Grid cell marking as inside, outside, or overlapping with surface
460 elements failed.  Please report the issue to the SPARTA developers.
461 
462 E: Grid cells marked as unknown = %d
463 
464 Grid cell marking as inside, outside, or overlapping with surface
465 elements did not successfully mark all cells.  Please report the issue
466 to the SPARTA developers.
467 
468 W: Grid cell interior corner points marked as unknown = %d
469 
470 Corner points of grid cells interior to the simulation domain were not
471 all marked successfully as inside, outside, or overlapping with
472 surface elements.  This should normally not happen, but does
473 not affect simulations.
474 
475 E: Grid cell corner points on boundary marked as unknown = %d
476 
477 Corner points of grid cells on the boundary of the simulation domain
478 were not all marked successfully as inside, outside, or overlapping
479 with surface elements.  Please report the issue to the SPARTA
480 developers.
481 
482 E: Cannot weight cells before grid is defined
483 
484 Self-explanatory.
485 
486 E: Illegal ... command
487 
488 Self-explanatory.  Check the input script syntax and compare to the
489 documentation for the command.  You can use -echo screen as a
490 command-line option when running SPARTA to see the offending line.
491 
492 E: Cannot use weight cell radius unless axisymmetric
493 
494 An axisymmetric model is required for this style of cell weighting.
495 
496 */
497