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_eflux_grid.h"
17 #include "particle.h"
18 #include "mixture.h"
19 #include "grid.h"
20 #include "update.h"
21 #include "modify.h"
22 #include "memory.h"
23 #include "error.h"
24
25 using namespace SPARTA_NS;
26
27 // user keywords
28
29 enum{HEATX,HEATY,HEATZ};
30
31 // internal accumulators
32
33 enum{MASSSUM,mVx,mVy,mVz,mVxVx,mVyVy,mVzVz,mVxVy,mVyVz,mVxVz,
34 mVxVxVx,mVyVyVy,mVzVzVz,
35 mVxVyVy,mVxVzVz,mVyVxVx,mVyVzVz,mVzVxVx,mVzVyVy,LASTSIZE};
36
37 // max # of quantities to accumulate for any user value
38
39 #define MAXACCUMULATE 12
40
41 /* ---------------------------------------------------------------------- */
42
ComputeEFluxGrid(SPARTA * sparta,int narg,char ** arg)43 ComputeEFluxGrid::ComputeEFluxGrid(SPARTA *sparta, int narg, char **arg) :
44 Compute(sparta, narg, arg)
45 {
46 if (narg < 5) error->all(FLERR,"Illegal compute eflux/grid command");
47
48 int igroup = grid->find_group(arg[2]);
49 if (igroup < 0)
50 error->all(FLERR,"Compute eflux/grid group ID does not exist");
51 groupbit = grid->bitmask[igroup];
52
53 imix = particle->find_mixture(arg[3]);
54 if (imix < 0) error->all(FLERR,"Compute eflux/grid mixture ID "
55 "does not exist");
56 ngroup = particle->mixture[imix]->ngroup;
57
58 nvalue = narg - 4;
59 value = new int[nvalue];
60
61 npergroup = 0;
62 unique = new int[LASTSIZE];
63 nmap = new int[nvalue];
64 memory->create(map,ngroup*nvalue,MAXACCUMULATE,"eflux/grid:map");
65 for (int i = 0; i < nvalue; i++) nmap[i] = 0;
66
67 int ivalue = 0;
68 int iarg = 4;
69 while (iarg < narg) {
70 if (strcmp(arg[iarg],"heatx") == 0) {
71 value[ivalue] = HEATX;
72 set_map(ivalue,MASSSUM);
73 set_map(ivalue,mVx);
74 set_map(ivalue,mVy);
75 set_map(ivalue,mVz);
76 set_map(ivalue,mVxVx);
77 set_map(ivalue,mVyVy);
78 set_map(ivalue,mVzVz);
79 set_map(ivalue,mVxVy);
80 set_map(ivalue,mVxVz);
81 set_map(ivalue,mVxVxVx);
82 set_map(ivalue,mVxVyVy);
83 set_map(ivalue,mVxVzVz);
84 } else if (strcmp(arg[iarg],"heaty") == 0) {
85 value[ivalue] = HEATX;
86 set_map(ivalue,MASSSUM);
87 set_map(ivalue,mVy);
88 set_map(ivalue,mVx);
89 set_map(ivalue,mVz);
90 set_map(ivalue,mVyVy);
91 set_map(ivalue,mVxVx);
92 set_map(ivalue,mVzVz);
93 set_map(ivalue,mVxVy);
94 set_map(ivalue,mVyVz);
95 set_map(ivalue,mVyVyVy);
96 set_map(ivalue,mVyVxVx);
97 set_map(ivalue,mVyVzVz);
98 } else if (strcmp(arg[iarg],"heatz") == 0) {
99 value[ivalue] = HEATX;
100 set_map(ivalue,MASSSUM);
101 set_map(ivalue,mVz);
102 set_map(ivalue,mVx);
103 set_map(ivalue,mVy);
104 set_map(ivalue,mVzVz);
105 set_map(ivalue,mVxVx);
106 set_map(ivalue,mVyVy);
107 set_map(ivalue,mVxVz);
108 set_map(ivalue,mVyVz);
109 set_map(ivalue,mVzVzVz);
110 set_map(ivalue,mVzVxVx);
111 set_map(ivalue,mVzVyVy);
112 } else error->all(FLERR,"Illegal compute eflux/grid command");
113
114 ivalue++;
115 iarg++;
116 }
117
118 // ntotal = total # of columns in tally array
119 // reset_map() adjusts indices in initial map() using final npergroup
120 // also adds columns to tally array for CELLCOUNT/CELLMASS
121
122 ntotal = ngroup*npergroup;
123 reset_map();
124
125 per_grid_flag = 1;
126 size_per_grid_cols = ngroup*nvalue;
127 post_process_grid_flag = 1;
128
129 nglocal = 0;
130 vector_grid = NULL;
131 tally = NULL;
132 }
133
134 /* ---------------------------------------------------------------------- */
135
~ComputeEFluxGrid()136 ComputeEFluxGrid::~ComputeEFluxGrid()
137 {
138 if (copymode) return;
139
140 delete [] value;
141 delete [] unique;
142
143 delete [] nmap;
144 memory->destroy(map);
145
146 memory->destroy(vector_grid);
147 memory->destroy(tally);
148 }
149
150 /* ---------------------------------------------------------------------- */
151
init()152 void ComputeEFluxGrid::init()
153 {
154 if (ngroup != particle->mixture[imix]->ngroup)
155 error->all(FLERR,"Number of groups in compute eflux/grid mixture "
156 "has changed");
157
158 reallocate();
159 }
160
161 /* ---------------------------------------------------------------------- */
162
compute_per_grid()163 void ComputeEFluxGrid::compute_per_grid()
164 {
165 invoked_per_grid = update->ntimestep;
166
167 Grid::ChildInfo *cinfo = grid->cinfo;
168 Particle::Species *species = particle->species;
169 Particle::OnePart *particles = particle->particles;
170 int *s2g = particle->mixture[imix]->species2group;
171 int nlocal = particle->nlocal;
172
173 int i,j,k,m,ispecies,igroup,icell;
174 double mass;
175 double *v,*vec;
176
177 // zero all accumulators - could do this with memset()
178
179 for (i = 0; i < nglocal; i++)
180 for (j = 0; j < ntotal; j++)
181 tally[i][j] = 0.0;
182
183 // loop over all particles, skip species not in mixture group
184 // perform all tallies needed for each particle
185 // depends on its species group and the user-requested values
186
187 for (i = 0; i < nlocal; i++) {
188 ispecies = particles[i].ispecies;
189 igroup = s2g[ispecies];
190 if (igroup < 0) continue;
191 icell = particles[i].icell;
192 if (!(cinfo[icell].mask & groupbit)) continue;
193
194 mass = species[ispecies].mass;
195 v = particles[i].v;
196
197 vec = tally[icell];
198
199 // loop has all possible values particle needs to accumulate
200 // subset defined by user values are indexed by accumulate vector
201 // NOTE: at some point may need prefactors v,v^2,v^3 converted to p,eng
202
203 k = igroup*npergroup;
204
205 for (m = 0; m < npergroup; m++) {
206 switch (unique[m]) {
207 case MASSSUM:
208 vec[k++] += mass;
209 break;
210 case mVx:
211 vec[k++] += mass*v[0];
212 break;
213 case mVy:
214 vec[k++] += mass*v[1];
215 break;
216 case mVz:
217 vec[k++] += mass*v[2];
218 break;
219 case mVxVx:
220 vec[k++] += mass*v[0]*v[0];
221 break;
222 case mVyVy:
223 vec[k++] += mass*v[1]*v[1];
224 break;
225 case mVzVz:
226 vec[k++] += mass*v[2]*v[2];
227 break;
228 case mVxVy:
229 vec[k++] += mass*v[0]*v[1];
230 break;
231 case mVyVz:
232 vec[k++] += mass*v[1]*v[2];
233 break;
234 case mVxVz:
235 vec[k++] += mass*v[0]*v[2];
236 break;
237 case mVxVxVx:
238 vec[k++] += mass*v[0]*v[0]*v[0];
239 break;
240 case mVyVyVy:
241 vec[k++] += mass*v[1]*v[1]*v[1];
242 break;
243 case mVzVzVz:
244 vec[k++] += mass*v[2]*v[2]*v[2];
245 break;
246 case mVxVyVy:
247 vec[k++] += mass*v[0]*v[1]*v[1];
248 break;
249 case mVxVzVz:
250 vec[k++] += mass*v[0]*v[2]*v[2];
251 break;
252 case mVyVxVx:
253 vec[k++] += mass*v[1]*v[0]*v[0];
254 break;
255 case mVyVzVz:
256 vec[k++] += mass*v[1]*v[2]*v[2];
257 break;
258 case mVzVxVx:
259 vec[k++] += mass*v[2]*v[0]*v[0];
260 break;
261 case mVzVyVy:
262 vec[k++] += mass*v[2]*v[1]*v[1];
263 break;
264 }
265 }
266 }
267 }
268
269 /* ----------------------------------------------------------------------
270 query info about internal tally array for this compute
271 index = which column of output (0 for vec, 1 to N for array)
272 return # of tally quantities for this index
273 also return array = ptr to tally array
274 also return cols = ptr to list of columns in tally for this index
275 ------------------------------------------------------------------------- */
276
query_tally_grid(int index,double ** & array,int * & cols)277 int ComputeEFluxGrid::query_tally_grid(int index, double **&array, int *&cols)
278 {
279 index--;
280 int ivalue = index % nvalue;
281 array = tally;
282 cols = map[index];
283 return nmap[ivalue];
284 }
285
286 /* ----------------------------------------------------------------------
287 tally accumulated info to compute final normalized values
288 index = which column of output (0 for vec, 1 to N for array)
289 for etally = NULL:
290 use internal tallied info for single timestep, set nsample = 1
291 compute values for all grid cells
292 store results in vector_grid with nstride = 1 (single col of array_grid)
293 for etally = ptr to caller array:
294 use external tallied info for many timesteps
295 nsample = additional normalization factor used by some values
296 emap = list of etally columns to use, # of columns determined by index
297 store results in caller's vec, spaced by nstride
298 if norm = 0.0, set result to 0.0 directly so do not divide by 0.0
299 ------------------------------------------------------------------------- */
300
post_process_grid(int index,int nsample,double ** etally,int * emap,double * vec,int nstride)301 void ComputeEFluxGrid::post_process_grid(int index, int nsample,
302 double **etally, int *emap,
303 double *vec, int nstride)
304 {
305 index--;
306 int ivalue = index % nvalue;
307
308 int lo = 0;
309 int hi = nglocal;
310 int k = 0;
311
312 if (!etally) {
313 nsample = 1;
314 etally = tally;
315 emap = map[index];
316 vec = vector_grid;
317 nstride = 1;
318 }
319
320 // compute normalized final value for each grid cell
321 // Vcm = Sum mv / Sum m = (Wx,Wy,Wz)
322 // Wi = Sum mVi / M
323 // heati = 0.5 * F/V Sum m (Vi - Wi) (V - W)^2
324 // (Vi - Wi) (V - W)^2 = (Vi - Wi)
325 // [ (Vi - Wi)^2 + (V1 - W1)^2 + (V2 - W2)^2 ]
326 // i = xyz and 1,2 = xyx indices different than i
327 // heati = 3 terms = h+h1+h2 where h2 is same as h1 with V1 replaced by V2
328 // h = F/V Sum m (Vi-Wi)^3
329 // (Vi-Wi)^3 = Vi^3 - 3WiVi^2 + 3Wi^2Vi - Wi^3
330 // Sum m (Vi-Wi)^3 = Sum(mVi^3) - 3 Sum(mVi^2) Sum(mVi) / M +
331 // 2 Sum(mVi)^3 / M^2
332 // h1 = F/V Sum m (Vi-Wi) (V1-W1)^2
333 // (Vi-Wi) (V1-W1)^2 = ViV1^2 - 2ViV1W1 + ViW1^2 - WiV1^2 + 2V1WiW1 - WiW1^2
334 // 3 terms in previous equation combine to 1 term in next equation
335 // Sum m (Vi-Wi) (V1-W1)^2 = Sum(mViV1^2) - 2 Sum(mViV1) Sum(mV1) / M -
336 // Sum(mVi) Sum(mV1^2) / M +
337 // 2 Sum(mVi) Sum(mV1)^2 / M^2
338
339 double summass,h,h1,h2,wt;
340 double *t;
341
342 double fnum = update->fnum;
343 Grid::ChildInfo *cinfo = grid->cinfo;
344
345 int mass = emap[0];
346 int mv = emap[1];
347 int mv1 = emap[2];
348 int mv2 = emap[3];
349 int mvv = emap[4];
350 int mv1v1 = emap[5];
351 int mv2v2 = emap[6];
352 int mvv1 = emap[7];
353 int mvv2 = emap[8];
354 int mvvv = emap[9];
355 int mvv1v1 = emap[10];
356 int mvv2v2 = emap[11];
357
358 for (int icell = lo; icell < hi; icell++) {
359 t = etally[icell];
360 summass = t[mass];
361 if (summass == 0.0) vec[k] = 0.0;
362 else {
363 h = t[mvvv] - 3.0*t[mv]*t[mvv]/summass +
364 2.0*t[mv]*t[mv]*t[mv]/summass/summass;
365 h1 = t[mvv1v1] - 2.0*t[mvv1]*t[mv1]/summass - t[mv]*t[mv1v1]/summass +
366 2.0*t[mv]*t[mv1]*t[mv1]/summass/summass;
367 h2 = t[mvv2v2] - 2.0*t[mvv2]*t[mv2]/summass - t[mv]*t[mv2v2]/summass +
368 2.0*t[mv]*t[mv2]*t[mv2]/summass/summass;
369 wt = 0.5 * fnum * cinfo[icell].weight / cinfo[icell].volume;
370 vec[k] = wt/nsample * (h + h1 + h2);
371 }
372 k += nstride;
373 }
374 }
375
376 /* ----------------------------------------------------------------------
377 add a tally quantity to all groups for ivalue
378 also add it to unique list if first time this name is used
379 name = name of tally quantity from enum{} at top of file
380 nmap[i] = # of tally quantities for user value I
381 map[i][k] = index of Kth tally quantity for output value I
382 npergroup = length of unique list
383 ------------------------------------------------------------------------- */
384
set_map(int ivalue,int name)385 void ComputeEFluxGrid::set_map(int ivalue, int name)
386 {
387 // index = loc of name in current unique list if there, else npergroup
388
389 int index = 0;
390 for (index = 0; index < npergroup; index++)
391 if (unique[index] == name) break;
392
393 // if name is not already in unique, add it and increment npergroup
394
395 if (index == npergroup) {
396 index = npergroup;
397 unique[npergroup++] = name;
398 }
399
400 // add index to map and nmap for all groups
401 // will add group offset in reset_map()
402
403 for (int igroup = 0; igroup < ngroup; igroup++)
404 map[igroup*nvalue+ivalue][nmap[ivalue]] = index;
405 nmap[ivalue]++;
406 }
407
408 /* ----------------------------------------------------------------------
409 reset map indices to reflect final npergroup = unique quantities/group
410 ------------------------------------------------------------------------- */
411
reset_map()412 void ComputeEFluxGrid::reset_map()
413 {
414 for (int i = 0; i < ngroup*nvalue; i++) {
415 int igroup = i / nvalue;
416 int ivalue = i % nvalue;
417 for (int k = 0; k < nmap[ivalue]; k++)
418 map[i][k] += igroup*npergroup;
419 }
420 }
421
422 /* ----------------------------------------------------------------------
423 reallocate data storage if nglocal has changed
424 called by init() and whenever grid changes
425 ------------------------------------------------------------------------- */
426
reallocate()427 void ComputeEFluxGrid::reallocate()
428 {
429 if (grid->nlocal == nglocal) return;
430
431 memory->destroy(vector_grid);
432 memory->destroy(tally);
433 nglocal = grid->nlocal;
434 memory->create(vector_grid,nglocal,"grid:vector_grid");
435 memory->create(tally,nglocal,ntotal,"grid:tally");
436 }
437
438 /* ----------------------------------------------------------------------
439 memory usage of local grid-based data
440 ------------------------------------------------------------------------- */
441
memory_usage()442 bigint ComputeEFluxGrid::memory_usage()
443 {
444 bigint bytes;
445 bytes = nglocal * sizeof(double);
446 bytes = ntotal*nglocal * sizeof(double);
447 return bytes;
448 }
449