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