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 "compute_lambda_grid.h"
17 #include "update.h"
18 #include "grid.h"
19 #include "domain.h"
20 #include "collide.h"
21 #include "modify.h"
22 #include "fix.h"
23 #include "compute.h"
24 #include "math_const.h"
25 #include "memory.h"
26 #include "error.h"
27 
28 using namespace SPARTA_NS;
29 using namespace MathConst;
30 
31 enum{NONE,COMPUTE,FIX};
32 enum{KNONE,KALL,KX,KY,KZ};
33 
34 #define INVOKED_PER_GRID 16
35 #define BIG 1.0e20
36 
37 /* ---------------------------------------------------------------------- */
38 
ComputeLambdaGrid(SPARTA * sparta,int narg,char ** arg)39 ComputeLambdaGrid::ComputeLambdaGrid(SPARTA *sparta, int narg, char **arg) :
40   Compute(sparta, narg, arg)
41 {
42   if (narg < 5 || narg > 6)
43     error->all(FLERR,"Illegal compute lambda/grid command");
44 
45   // parse three required input fields
46   // customize a new keyword by adding to if statement
47 
48   id_nrho = id_temp = NULL;
49 
50   if (strncmp(arg[2],"c_",2) == 0 || strncmp(arg[2],"f_",2) == 0) {
51     int n = strlen(arg[2]);
52     id_nrho = new char[n];
53     strcpy(id_nrho,&arg[2][2]);
54 
55     char *ptr = strchr(id_nrho,'[');
56     if (ptr) {
57       if (id_nrho[strlen(id_nrho)-1] != ']')
58         error->all(FLERR,"Invalid nrho in compute lambda/grid command");
59       nrhoindex = atoi(ptr+1);
60       *ptr = '\0';
61     } else nrhoindex = 0;
62 
63     if (strncmp(arg[2],"c_",2) == 0) nrhowhich = COMPUTE;
64     else nrhowhich = FIX;
65 
66     if (nrhowhich == COMPUTE) {
67       int n = modify->find_compute(id_nrho);
68       if (n < 0)
69         error->all(FLERR,"Could not find compute lambda/grid compute ID");
70       if (modify->compute[n]->per_grid_flag == 0)
71 	error->all(FLERR,"Compute lambda/grid compute does not "
72                    "compute per-grid info");
73       if (nrhoindex == 0 && modify->compute[n]->size_per_grid_cols > 0)
74 	error->all(FLERR,
75 		   "Compute lambda/grid compute does not "
76                    "compute per-grid vector");
77       if (nrhoindex > 0 && modify->compute[n]->size_per_grid_cols == 0)
78 	error->all(FLERR,
79 		   "Compute lambda/grid compute does not "
80                    "compute per-grid array");
81       if (nrhoindex > 0 && nrhoindex > modify->compute[n]->size_per_grid_cols)
82 	error->all(FLERR,"Compute lambda compute vector is "
83                    "accessed out-of-range");
84     } else {
85       int n = modify->find_fix(id_nrho);
86       if (n < 0) error->all(FLERR,"Could not find compute lambda/grid fix ID");
87       if (modify->fix[n]->per_grid_flag == 0)
88 	error->all(FLERR,"Compute lambda/grid fix does not "
89                    "compute per-grid info");
90       if (nrhoindex == 0 && modify->fix[n]->size_per_grid_cols > 0)
91 	error->all(FLERR,"Compute lambda/grid fix does not "
92                    "compute per-grid vector");
93       if (nrhoindex > 0 && modify->fix[n]->size_per_grid_cols == 0)
94 	error->all(FLERR,"Compute lambda/grid fix does not "
95                    "compute per-grid array");
96       if (nrhoindex > 0 && nrhoindex > modify->fix[n]->size_per_grid_cols)
97 	error->all(FLERR,"Compute lambda/grid fix array is "
98                    "accessed out-of-range");
99     }
100   } else error->all(FLERR,"Illegal compute lambda/grid command");
101 
102   if (strncmp(arg[3],"c_",2) == 0 || strncmp(arg[3],"f_",2) == 0) {
103     int n = strlen(arg[3]);
104     id_temp = new char[n];
105     strcpy(id_temp,&arg[3][2]);
106 
107     char *ptr = strchr(id_temp,'[');
108     if (ptr) {
109       if (id_temp[strlen(id_temp)-1] != ']')
110         error->all(FLERR,"Invalid temp in compute lambda/grid command");
111       tempindex = atoi(ptr+1);
112       *ptr = '\0';
113     } else tempindex = 0;
114 
115     if (strncmp(arg[3],"c_",2) == 0) tempwhich = COMPUTE;
116     else tempwhich = FIX;
117 
118     if (tempwhich == COMPUTE) {
119       int n = modify->find_compute(id_temp);
120       if (n < 0)
121         error->all(FLERR,"Could not find compute lambda/grid compute ID");
122       if (modify->compute[n]->per_grid_flag == 0)
123 	error->all(FLERR,"Compute lambda/grid compute does not "
124                    "compute per-grid info");
125       if (tempindex == 0 && modify->compute[n]->size_per_grid_cols > 0)
126 	error->all(FLERR,
127 		   "Compute lambda/grid compute does not "
128                    "compute per-grid vector");
129       if (tempindex > 0 && modify->compute[n]->size_per_grid_cols == 0)
130 	error->all(FLERR,
131 		   "Compute lambda/grid compute does not "
132                    "compute per-grid array");
133       if (tempindex > 0 && tempindex > modify->compute[n]->size_per_grid_cols)
134 	error->all(FLERR,"Compute lambda compute vector is "
135                    "accessed out-of-range");
136     } else {
137       int n = modify->find_fix(id_temp);
138       if (n < 0) error->all(FLERR,"Could not find compute lambda/grid fix ID");
139       if (modify->fix[n]->per_grid_flag == 0)
140 	error->all(FLERR,"Compute lambda/grid fix does not "
141                    "compute per-grid info");
142       if (tempindex == 0 && modify->fix[n]->size_per_grid_cols > 0)
143 	error->all(FLERR,"Compute lambda/grid fix does not "
144                    "compute per-grid vector");
145       if (tempindex > 0 && modify->fix[n]->size_per_grid_cols == 0)
146 	error->all(FLERR,"Compute lambda/grid fix does not "
147                    "compute per-grid array");
148       if (tempindex > 0 && tempindex > modify->fix[n]->size_per_grid_cols)
149 	error->all(FLERR,"Compute lambda/grid fix array is "
150                    "accessed out-of-range");
151     }
152   } else if (strcmp(arg[3],"NULL") == 0) tempwhich = NONE;
153   else error->all(FLERR,"Illegal compute lambda/grid command");
154 
155   int n = strlen(arg[4]) + 1;
156   species = new char[n];
157   strcpy(species,arg[4]);
158 
159   // optional arg
160 
161   kflag = KNONE;
162   if (narg == 6) {
163     if (strcmp(arg[5],"kall") == 0) kflag = KALL;
164     else if (strcmp(arg[5],"kx") == 0) kflag = KX;
165     else if (strcmp(arg[5],"ky") == 0) kflag = KY;
166     else if (strcmp(arg[5],"kz") == 0) kflag = KZ;
167     else error->all(FLERR,"Illegal compute lambda/grid command");
168   }
169 
170   if (kflag == KZ && domain->dimension == 2)
171     error->all(FLERR,"Cannot use compute lambda/grid kz for 2d simulation");
172 
173   // initialize data structures
174 
175   per_grid_flag = 1;
176   if (kflag == KNONE) size_per_grid_cols = 0;
177   else size_per_grid_cols = 2;
178 
179   nglocal = 0;
180   vector_grid = NULL;
181   array_grid = NULL;
182   nrho = temp = NULL;
183 }
184 
185 /* ---------------------------------------------------------------------- */
186 
~ComputeLambdaGrid()187 ComputeLambdaGrid::~ComputeLambdaGrid()
188 {
189   if (copymode) return;
190 
191   delete [] id_nrho;
192   delete [] id_temp;
193   delete [] species;
194   memory->destroy(vector_grid);
195   memory->destroy(array_grid);
196   memory->destroy(nrho);
197   memory->destroy(temp);
198 }
199 
200 /* ---------------------------------------------------------------------- */
201 
init()202 void ComputeLambdaGrid::init()
203 {
204   reallocate();
205 
206   // setup computes and fixes
207 
208   if (nrhowhich == COMPUTE) {
209     int icompute = modify->find_compute(id_nrho);
210     if (icompute < 0)
211       error->all(FLERR,"Could not find compute lambda/grid compute ID");
212     cnrho = modify->compute[icompute];
213   } else if (nrhowhich == FIX) {
214     int ifix = modify->find_fix(id_nrho);
215     if (ifix < 0)
216       error->all(FLERR,"Could not find compute lambda/grid fix ID");
217     fnrho = modify->fix[ifix];
218   }
219 
220   if (tempwhich == COMPUTE) {
221     int icompute = modify->find_compute(id_temp);
222     if (icompute < 0)
223       error->all(FLERR,"Could not find compute lambda/grid compute ID");
224     ctemp = modify->compute[icompute];
225   } else if (tempwhich == FIX) {
226     int ifix = modify->find_fix(id_temp);
227     if (ifix < 0)
228       error->all(FLERR,"Could not find compute lambda/grid fix ID");
229     ftemp = modify->fix[ifix];
230   }
231 
232   // setup collision parameter
233 
234   if (!collide)
235     error->all(FLERR,"Compute lambda/grid requires a "
236                "collision style be defined");
237 
238   int ispecies = particle->find_species(species);
239   if (ispecies < 0)
240     error->all(FLERR,"Compute lambda/grid species is not defined");
241 
242   dref = collide->extract(ispecies,ispecies,"diam");
243   tref = collide->extract(ispecies,ispecies,"tref");
244   omega = collide->extract(ispecies,ispecies,"omega");
245   prefactor = sqrt(2.0) * MY_PI * dref*dref;
246 }
247 
248 /* ---------------------------------------------------------------------- */
249 
compute_per_grid()250 void ComputeLambdaGrid::compute_per_grid()
251 {
252   invoked_per_grid = update->ntimestep;
253 
254   if (nrhowhich == FIX && update->ntimestep % fnrho->per_grid_freq)
255     error->all(FLERR,"Compute lambda/grid fix not computed at compatible time");
256   if (tempwhich == FIX && update->ntimestep % ftemp->per_grid_freq)
257     error->all(FLERR,"Compute lambda/grid fix not computed at compatible time");
258 
259   // grab nrho and temp values from compute or fix
260   // invoke nrho and temp computes as needed
261 
262   if (nrhowhich == COMPUTE) {
263     if (!(cnrho->invoked_flag & INVOKED_PER_GRID)) {
264       cnrho->compute_per_grid();
265       cnrho->invoked_flag |= INVOKED_PER_GRID;
266     }
267 
268     if (cnrho->post_process_grid_flag)
269       cnrho->post_process_grid(nrhoindex,1,NULL,NULL,NULL,1);
270 
271     if (nrhoindex == 0 || cnrho->post_process_grid_flag)
272       memcpy(nrho,cnrho->vector_grid,nglocal*sizeof(double));
273     else {
274       int index = nrhoindex-1;
275       double **array = cnrho->array_grid;
276       for (int i = 0; i < nglocal; i++)
277         nrho[i] = array[i][index];
278     }
279 
280   } else if (nrhowhich == FIX) {
281     if (nrhoindex == 0) {
282       memcpy(nrho,fnrho->vector_grid,nglocal*sizeof(double));
283     } else {
284       double **array = fnrho->array_grid;
285       int index = nrhoindex-1;
286       for (int i = 0; i < nglocal; i++)
287         nrho[i] = array[i][index];
288     }
289   }
290 
291   if (tempwhich == COMPUTE) {
292     if (!(ctemp->invoked_flag & INVOKED_PER_GRID)) {
293       ctemp->compute_per_grid();
294       ctemp->invoked_flag |= INVOKED_PER_GRID;
295     }
296 
297     if (ctemp->post_process_grid_flag)
298       ctemp->post_process_grid(tempindex,1,NULL,NULL,NULL,1);
299 
300     if (tempindex == 0 || ctemp->post_process_grid_flag)
301       memcpy(temp,ctemp->vector_grid,nglocal*sizeof(double));
302     else {
303       int index = tempindex-1;
304       double **array = ctemp->array_grid;
305       for (int i = 0; i < nglocal; i++)
306         temp[i] = array[i][index];
307     }
308 
309   } else if (tempwhich == FIX) {
310     if (tempindex == 0)
311       memcpy(temp,ftemp->vector_grid,nglocal*sizeof(double));
312     else {
313       double **array = ftemp->array_grid;
314       int index = tempindex-1;
315       for (int i = 0; i < nglocal; i++)
316         temp[i] = array[i][index];
317     }
318   }
319 
320   // compute mean free path for each grid cell
321   // formula from Bird, eq 4.65
322 
323   double lambda;
324 
325   for (int i = 0; i < nglocal; i++) {
326     if (nrho[i] == 0.0) lambda = BIG;
327     else if (tempwhich == NONE || temp[i] == 0.0)
328       lambda = 1.0 / (prefactor * nrho[i]);
329     else
330       lambda = 1.0 / (prefactor * nrho[i] * pow(tref/temp[i],omega-0.5));
331 
332     if (kflag == KNONE) vector_grid[i] = lambda;
333     else array_grid[i][0] = lambda;
334   }
335 
336   // calculate per-cell Knudsen number
337 
338   if (kflag == KNONE) return;
339 
340   Grid::ChildCell *cells = grid->cells;
341   int dimension = domain->dimension;
342 
343   if (kflag == KALL) {
344     double size;
345     for (int i = 0; i < nglocal; i++) {
346       size = (cells[i].hi[0] - cells[i].lo[0]);
347       size += (cells[i].hi[1] - cells[i].lo[1]);
348       if (dimension == 2) size *= 0.5;
349       else {
350         size += (cells[i].hi[2] - cells[i].lo[2]);
351         size /= 3.0;
352       }
353       array_grid[i][1] = array_grid[i][0] / size;
354     }
355   } else if (kflag == KX) {
356     for (int i = 0; i < nglocal; i++)
357       array_grid[i][1] = array_grid[i][0] / (cells[i].hi[0] - cells[i].lo[0]);
358   } else if (kflag == KY) {
359     for (int i = 0; i < nglocal; i++)
360       array_grid[i][1] = array_grid[i][0] / (cells[i].hi[1] - cells[i].lo[1]);
361   } else if (kflag == KZ) {
362     for (int i = 0; i < nglocal; i++)
363       array_grid[i][1] = array_grid[i][0] / (cells[i].hi[2] - cells[i].lo[2]);
364   }
365 }
366 
367 /* ----------------------------------------------------------------------
368    reallocate arrays if nglocal has changed
369    called by init() and load balancer
370 ------------------------------------------------------------------------- */
371 
reallocate()372 void ComputeLambdaGrid::reallocate()
373 {
374   if (grid->nlocal == nglocal) return;
375 
376   nglocal = grid->nlocal;
377   if (kflag == KNONE) {
378     memory->destroy(vector_grid);
379     memory->create(vector_grid,nglocal,"lambda/grid:vector_grid");
380   } else {
381     memory->destroy(array_grid);
382     memory->create(array_grid,nglocal,2,"lambda/grid:array_grid");
383   }
384 
385   memory->destroy(nrho);
386   memory->create(nrho,nglocal,"lambda/grid:nrho");
387   if (tempwhich != NONE) {
388     memory->destroy(temp);
389     memory->create(temp,nglocal,"lambda/grid:temp");
390   }
391 }
392 
393 /* ----------------------------------------------------------------------
394    memory usage of local grid-based array
395 ------------------------------------------------------------------------- */
396 
memory_usage()397 bigint ComputeLambdaGrid::memory_usage()
398 {
399   bigint bytes;
400   bytes = nglocal * sizeof(double);
401   if (kflag != KNONE) bytes += nglocal * sizeof(double);
402   bytes += nglocal * sizeof(double);                            // nrho
403   if (tempwhich != KNONE) bytes += nglocal * sizeof(double);    // temp
404   return bytes;
405 }
406