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_UPDATE_H
16 #define SPARTA_UPDATE_H
17 
18 #include "math.h"
19 #include "pointers.h"
20 
21 namespace SPARTA_NS {
22 
23 class Update : protected Pointers {
24  public:
25   bigint ntimestep;               // current timestep
26   int nsteps;                     // # of steps to run
27   int runflag;                    // 0 for unset, 1 for run
28   bigint firststep,laststep;      // 1st & last step of this run
29   bigint beginstep,endstep;       // 1st and last step of multiple runs
30   int first_update;               // 0 before initial update, 1 after
31   double dt;                      // timestep size
32 
33   char *unit_style;      // style of units used throughout simulation
34   double boltz;          // Boltzmann constant (eng/degree K)
35   double mvv2e;          // conversion of mv^2 to energy
36 
37   double fnum;           // ratio of real particles to simulation particles
38   double nrho;           // number density of background gas
39   double vstream[3];     // streaming velocity of background gas
40   double temp_thermal;   // thermal temperature of background gas
41   double gravity[3];     // acceleration vector of gravity
42 
43   int nmigrate;          // # of particles to migrate to new procs
44   int *mlist;            // indices of particles to migrate
45 
46                          // current step counters
47   int niterate;          // iterations of move/comm
48   int ntouch_one;        // particle-cell touches
49   int ncomm_one;         // particles migrating to new procs
50   int nboundary_one;     // particles colliding with global boundary
51   int nexit_one;         // particles exiting outflow boundary
52   int nscheck_one;       // surface elements checked for collisions
53   int nscollide_one;     // particle/surface collisions
54 
55   bigint first_running_step; // timestep running counts start on
56   int niterate_running;      // running count of move/comm interations
57   bigint nmove_running;      // running count of total particle moves
58   bigint ntouch_running;     // running count of current step counters
59   bigint ncomm_running;
60   bigint nboundary_running;
61   bigint nexit_running;
62   bigint nscheck_running;
63   bigint nscollide_running;
64 
65   int nstuck;                // # of particles stuck on surfs and deleted
66   int naxibad;               // # of particles where axisymm move was bad
67                              // in this case, bad means particle ended up
68                              // outside of final cell curved surf by epsilon
69                              // when move logic thinks it is inside cell
70 
71   int reorder_period;        // # of timesteps between particle reordering
72   int global_mem_limit;      // max # of bytes in arrays for rebalance and reordering
73   int mem_limit_grid_flag;   // 1 if using size of grid as memory limit
74   void set_mem_limit_grid(int gnlocal = 0);
75   int have_mem_limit();      // 1 if have memory limit
76 
77   int copymode;          // 1 if copy of class (prevents deallocation of
78                          //  base class when child copy is destroyed)
79 
80   class RanMars *ranmaster;   // master random number generator
81 
82   double rcblo[3],rcbhi[3];    // debug info from RCB for dump image
83 
84   Update(class SPARTA *);
85   ~Update();
86   void set_units(const char *);
87   virtual void init();
88   virtual void setup();
89   virtual void run(int);
90   void global(int, char **);
91   void reset_timestep(int, char **);
92 
93   int split3d(int, double *);
94   int split2d(int, double *);
95 
96  protected:
97   int me,nprocs;
98   int maxmigrate;            // max # of particles in mlist
99   class RanKnuth *random;     // RNG for particle timestep moves
100 
101   int collide_react;         // 1 if any SurfCollide or React classes defined
102   int nsc,nsr;               // copy of Collide/React data in Surf class
103   class SurfCollide **sc;
104   class SurfReact **sr;
105 
106   int bounce_tally;               // 1 if any bounces are ever tallied
107   int nslist_compute;             // # of computes that tally surf bounces
108   int nblist_compute;             // # of computes that tally boundary bounces
109   class Compute **slist_compute;  // list of all surf bounce Computes
110   class Compute **blist_compute;  // list of all boundary bounce Computes
111 
112   int nsurf_tally;         // # of Cmp tallying surf bounce info this step
113   int nboundary_tally;     // # of Cmp tallying boundary bounce info this step
114   class Compute **slist_active;   // list of active surf Computes this step
115   class Compute **blist_active;   // list of active boundary Computes this step
116 
117   int surf_pre_tally;       // 1 to log particle stats before surf collide
118   int boundary_pre_tally;   // 1 to log particle stats before boundary collide
119 
120   int collide_react_setup();
121   void collide_react_reset();
122   void collide_react_update();
123 
124   int bounce_setup();
125   virtual void bounce_set(bigint);
126 
127   int nulist_surfcollide;
128   SurfCollide **ulist_surfcollide;
129 
130   int dynamic;              // 1 if any classes do dynamic updates of params
131   void dynamic_setup();
132   void dynamic_update();
133 
134   void reset_timestep(bigint);
135 
136   //int axi_vertical_line(double, double *, double *, double, double, double,
137   //                     double &);
138 
139   // remap x and v components into axisymmetric plane
140   // input x at end of linear move (x = xold + dt*v)
141   // change x[1] = sqrt(x[1]^2 + x[2]^2), x[2] = 0.0
142   // change vy,vz by rotation into axisymmetric plane
143 
axi_remap(double * x,double * v)144   inline void axi_remap(double *x, double *v) {
145     double ynew = x[1];
146     double znew = x[2];
147     x[1] = sqrt(ynew*ynew + znew*znew);
148     x[2] = 0.0;
149     double rn = ynew / x[1];
150     double wn = znew / x[1];
151     double vy = v[1];
152     double vz = v[2];
153     v[1] = vy*rn + vz*wn;
154     v[2] = -vy*wn + vz*rn;
155   };
156 
157   typedef void (Update::*FnPtr)();
158   FnPtr moveptr;             // ptr to move method
159   template < int, int > void move();
160 
161   int perturbflag;
162   typedef void (Update::*FnPtr2)(double, double *, double *);
163   FnPtr2 moveperturb;        // ptr to moveperturb method
164 
165   // variants of moveperturb method
166   // adjust end-of-move x,v due to perturbation on straight-line advection
167 
gravity2d(double dt,double * x,double * v)168   inline void gravity2d(double dt, double *x, double *v) {
169     double dtsq = 0.5*dt*dt;
170     x[0] += dtsq*gravity[0];
171     x[1] += dtsq*gravity[1];
172     v[0] += dt*gravity[0];
173     v[1] += dt*gravity[1];
174   };
175 
gravity3d(double dt,double * x,double * v)176   inline void gravity3d(double dt, double *x, double *v) {
177     double dtsq = 0.5*dt*dt;
178     x[0] += dtsq*gravity[0];
179     x[1] += dtsq*gravity[1];
180     x[2] += dtsq*gravity[2];
181     v[0] += dt*gravity[0];
182     v[1] += dt*gravity[1];
183     v[2] += dt*gravity[2];
184   };
185 };
186 
187 }
188 
189 #endif
190 
191 /* ERROR/WARNING messages:
192 
193 E: Illegal ... command
194 
195 Self-explanatory.  Check the input script syntax and compare to the
196 documentation for the command.  You can use -echo screen as a
197 command-line option when running SPARTA to see the offending line.
198 
199 E: Gravity in z not allowed for 2d
200 
201 Self-explanatory.
202 
203 E: Gravity in y not allowed for axi-symmetric model
204 
205 Self-explanatory.
206 
207 E: Particle %d on proc %d hit inside of surf %d on step %ld
208 
209 This error should not happen if particles start outside of physical
210 objects.  Please report the issue to the SPARTA developers.
211 
212 E: Sending particle to self
213 
214 This error should not occur.  Please report the issue to the SPARTA
215 developers.
216 
217 E: Cannot set global surfmax when surfaces already exist
218 
219 This setting must be made before any surfac elements are
220 read via the read_surf command.
221 
222 E: Global mem/limit setting cannot exceed 2GB
223 
224 Self-expanatory, prevents 32-bit interger overflow
225 
226 E: Timestep must be >= 0
227 
228 Reset_timestep cannot be used to set a negative timestep.
229 
230 E: Too big a timestep
231 
232 Reset_timestep timestep value must fit in a SPARTA big integer, as
233 specified by the -DSPARTA_SMALL, -DSPARTA_BIG, or -DSPARTA_BIGBIG
234 options in the low-level Makefile used to build SPARTA.  See
235 Section 2.2 of the manual for details.
236 
237 E: Cannot reset timestep with a time-dependent fix defined
238 
239 The timestep cannot be reset when a fix that keeps track of elapsed
240 time is in place.
241 
242 */
243