PlasCom2  1.0
XPACC Multi-physics simluation application
ViscidUtil.C
Go to the documentation of this file.
1 #include "PCPPTypes.H"
2 #include "EulerUtil.H"
3 
4 namespace viscid {
5  namespace util {
6 
14  int ComputeTVBufferPower(const pcpp::IndexIntervalType &regionInterval,
15  const std::vector<size_t> &bufferSize,
16  const double * temperatureBuffer,
17  std::vector<double *> &tvBuffer,
18  const double beta,
19  const double power,
20  const double bulkViscFac,
21  const double specificHeat,
22  const double prandtlNumber)
23  {
24  int numDim = bufferSize.size();
25 
26  if(numDim < 1)
27  return(1);
28 
29  double *lmuPtr = tvBuffer[0];
30  double *llambdaPtr = tvBuffer[1];
31  double *lkappaPtr = tvBuffer[2];
32 
33  if(numDim == 1){
34 
35  size_t xStart = regionInterval[0].first;
36  size_t xEnd = regionInterval[0].second;
37  size_t xSize = bufferSize[0];
38  size_t nPoints = xSize;
39 
40  for(size_t iX = xStart;iX <= xEnd;iX++){
41  size_t xyzIndex = iX;
42  lmuPtr[xyzIndex] = beta*pow(temperatureBuffer[xyzIndex],power);
43  llambdaPtr[xyzIndex] = (bulkViscFac - 2.0/3.0)*lmuPtr[xyzIndex];
44  lkappaPtr[xyzIndex] = specificHeat*lmuPtr[xyzIndex]/prandtlNumber;
45  }
46  } else if (numDim == 2){
47 
48  size_t xStart = regionInterval[0].first;
49  size_t xEnd = regionInterval[0].second;
50  size_t yStart = regionInterval[1].first;
51  size_t yEnd = regionInterval[1].second;
52 
53  size_t xSize = bufferSize[0];
54  size_t ySize = bufferSize[1];
55  size_t nPoints = xSize*ySize;
56 
57  for(size_t iY = yStart;iY <= yEnd;iY++){
58  size_t yzIndex = iY*xSize;
59  for(size_t iX = xStart;iX <= xEnd;iX++){
60  size_t xyzIndex = yzIndex + iX;
61  lmuPtr[xyzIndex] = beta*pow(temperatureBuffer[xyzIndex],power);
62  llambdaPtr[xyzIndex] = (bulkViscFac - 2.0/3.0)*lmuPtr[xyzIndex];
63  lkappaPtr[xyzIndex] = specificHeat*lmuPtr[xyzIndex]/prandtlNumber;
64  }
65  }
66  } else if(numDim == 3) {
67  size_t xSize = bufferSize[0];
68  size_t ySize = bufferSize[1];
69  size_t zSize = bufferSize[2];
70 
71  size_t xStart = regionInterval[0].first;
72  size_t xEnd = regionInterval[0].second;
73  size_t yStart = regionInterval[1].first;
74  size_t yEnd = regionInterval[1].second;
75  size_t zStart = regionInterval[2].first;
76  size_t zEnd = regionInterval[2].second;
77 
78  size_t nPlane = xSize*ySize;
79  size_t nPoints = xSize*ySize*zSize;
80 
81  for(size_t iZ = zStart;iZ <= zEnd;iZ++){
82  size_t zIndex = iZ*nPlane;
83  for(size_t iY = yStart;iY <= yEnd;iY++){
84  size_t yzIndex = iY*xSize+zIndex;
85  for(size_t iX = xStart;iX <= xEnd;iX++){
86  size_t xyzIndex = yzIndex + iX;
87  lmuPtr[xyzIndex] = beta*pow(temperatureBuffer[xyzIndex],power);
88  llambdaPtr[xyzIndex] = (bulkViscFac - 2.0/3.0)*lmuPtr[xyzIndex];
89  lkappaPtr[xyzIndex] = specificHeat*lmuPtr[xyzIndex]/prandtlNumber;
90  }
91  }
92  }
93  }
94  return(0);
95  };
96 
97  int ComputeTVBuffer(const int *numDimPtr,
98  const size_t *bufferSize,
99  const size_t *numPointsPtr,
100  const size_t *bufferInterval,
101  const double *lrhoPtr,
102  const double *lrhoVPtr,
103  const double *lrhoEPtr,
104  double *lpressurePtr,
105  double *temperatureBuffer,
106  double *lrhom1Ptr,
107  double *lVPtr)
108  {
109 
110  return(-1);
111  };
112 
113  int ComputeTauBuffer(const pcpp::IndexIntervalType &regionInterval,
114  const std::vector<size_t> &bufferSize,
115  int gridType,
116  const std::vector<double> &gridMetrics,
117  const std::vector<double> &gridJacobian,
118  const double &Re,
119  const std::vector<double *> &tvBuffers,
120  const std::vector<double *> &velGradBuffers,
121  std::vector<double *> &tauBuffers)
122  {
123 
124  int numDim = bufferSize.size();
125 
126  if(numDim < 1)
127  return(1);
128 
129  double *lmuPtr = tvBuffers[0];
130  double *llambdaPtr = tvBuffers[1];
131  double ReM1=1/Re;
132 
133  if(gridType < simulation::grid::RECTILINEAR){ // uniform rectangular or Cartesian
134 
135  // Let ReM1 = J/Re for this type of grid
136  ReM1 *= gridJacobian[0];
137 
138  if(numDim == 1){
139 
140  size_t xStart = regionInterval[0].first;
141  size_t xEnd = regionInterval[0].second;
142  size_t xSize = bufferSize[0];
143  size_t nPoints = xSize;
144 
145  double* ltau11Ptr = tauBuffers[0];
146  double* velGradData0 = velGradBuffers[0];
147 
148  const double *ldudxiPtr = &velGradData0[0];
149  double dxidx = gridMetrics[0];
150 
151  for(size_t iX = xStart;iX <= xEnd;iX++){
152  size_t xyzIndex = iX;
153  double bulkTau = dxidx*ldudxiPtr[xyzIndex];
154 
155  ltau11Ptr[xyzIndex] = ReM1*(2.0*lmuPtr[xyzIndex]*dxidx*ldudxiPtr[xyzIndex] +
156  llambdaPtr[xyzIndex]*bulkTau);
157 
158  }
159  } else if (numDim == 2){
160 
161  size_t xStart = regionInterval[0].first;
162  size_t xEnd = regionInterval[0].second;
163  size_t yStart = regionInterval[1].first;
164  size_t yEnd = regionInterval[1].second;
165 
166  size_t xSize = bufferSize[0];
167  size_t ySize = bufferSize[1];
168  size_t nPoints = xSize*ySize;
169 
170  // brute force here, probably can be refactored to be more elegant
171  double* ltau11Ptr = tauBuffers[0];
172  double* ltau12Ptr = tauBuffers[1];
173  double* ltau22Ptr = tauBuffers[3];
174 
175  double* velGradData0 = velGradBuffers[0];
176  double* velGradData1 = velGradBuffers[1];
177 
178  const double *ldudxiPtr = &velGradData0[0];
179  const double *ldvdxiPtr = &velGradData0[nPoints];
180  const double *ldudetaPtr = &velGradData1[0];
181  const double *ldvdetaPtr = &velGradData1[nPoints];
182 
183  double dxidx = gridMetrics[0];
184  double detady = gridMetrics[1];
185 
186  for(size_t iY = yStart;iY <= yEnd;iY++){
187  size_t yzIndex = iY*xSize;
188  for(size_t iX = xStart;iX <= xEnd;iX++){
189 
190  size_t xyzIndex = yzIndex + iX;
191 
192  double dudx = ldudxiPtr[xyzIndex]*dxidx;
193  double dudy = ldudetaPtr[xyzIndex]*detady;
194 
195  double dvdx = ldvdxiPtr[xyzIndex]*dxidx;
196  double dvdy = ldvdetaPtr[xyzIndex]*detady;
197 
198  double bulkTau = dudx + dvdy;
199 
200  ltau11Ptr[xyzIndex] = ReM1*(2.0*lmuPtr[xyzIndex]*dudx +
201  llambdaPtr[xyzIndex]*bulkTau);
202 
203  ltau22Ptr[xyzIndex] = ReM1*(2.0*lmuPtr[xyzIndex]*dvdy +
204  llambdaPtr[xyzIndex]*bulkTau);
205 
206  ltau12Ptr[xyzIndex] = ReM1*lmuPtr[xyzIndex]*(dudy + dvdx);
207 
208 
209  }
210  }
211  } else if(numDim == 3) {
212 
213  size_t xSize = bufferSize[0];
214  size_t ySize = bufferSize[1];
215  size_t zSize = bufferSize[2];
216 
217  size_t xStart = regionInterval[0].first;
218  size_t xEnd = regionInterval[0].second;
219  size_t yStart = regionInterval[1].first;
220  size_t yEnd = regionInterval[1].second;
221  size_t zStart = regionInterval[2].first;
222  size_t zEnd = regionInterval[2].second;
223 
224  size_t nPlane = xSize*ySize;
225  size_t nPoints = xSize*ySize*zSize;
226 
227  // brute force here, probably can be refactored to be more elegant
228  double* ltau11Ptr = tauBuffers[0];
229  double* ltau12Ptr = tauBuffers[1];
230  double* ltau13Ptr = tauBuffers[2];
231  double* ltau22Ptr = tauBuffers[4];
232  double* ltau23Ptr = tauBuffers[5];
233  double* ltau33Ptr = tauBuffers[8];
234 
235  double* velGradData0 = velGradBuffers[0];
236  double* velGradData1 = velGradBuffers[1];
237  double* velGradData2 = velGradBuffers[2];
238 
239  const double* ldudxiPtr = &velGradData0[0];
240  const double* ldvdxiPtr = &velGradData0[nPoints];
241  const double* ldwdxiPtr = &velGradData0[2*nPoints];
242  const double* ldudetaPtr = &velGradData1[0];
243  const double* ldvdetaPtr = &velGradData1[nPoints];
244  const double* ldwdetaPtr = &velGradData1[2*nPoints];
245  const double* ldudzetaPtr = &velGradData2[0];
246  const double* ldvdzetaPtr = &velGradData2[nPoints];
247  const double* ldwdzetaPtr = &velGradData2[2*nPoints];
248 
249  double dxidx = gridMetrics[0];
250  double detady = gridMetrics[1];
251  double dzetadz = gridMetrics[2];
252 
253  for(size_t iZ = zStart;iZ <= zEnd;iZ++){
254  size_t zIndex = iZ*nPlane;
255  for(size_t iY = yStart;iY <= yEnd;iY++){
256  size_t yzIndex = iY*xSize+zIndex;
257  for(size_t iX = xStart;iX <= xEnd;iX++){
258  size_t xyzIndex = yzIndex + iX;
259 
260  double dudx = ldudxiPtr[xyzIndex]*dxidx;
261  double dvdx = ldvdxiPtr[xyzIndex]*dxidx;
262  double dwdx = ldwdxiPtr[xyzIndex]*dxidx;
263 
264  double dudy = ldudetaPtr[xyzIndex]*detady;
265  double dvdy = ldvdetaPtr[xyzIndex]*detady;
266  double dwdy = ldwdetaPtr[xyzIndex]*detady;
267 
268  double dudz = ldudzetaPtr[xyzIndex]*dzetadz;
269  double dvdz = ldvdzetaPtr[xyzIndex]*dzetadz;
270  double dwdz = ldwdzetaPtr[xyzIndex]*dzetadz;
271 
272  double bulkTau = dudx + dvdy + dwdz;
273 
274  ltau11Ptr[xyzIndex] = ReM1*(2.0*lmuPtr[xyzIndex]*dudx +
275  llambdaPtr[xyzIndex]*bulkTau);
276  ltau22Ptr[xyzIndex] = ReM1*(2.0*lmuPtr[xyzIndex]*dvdy +
277  llambdaPtr[xyzIndex]*bulkTau);
278  ltau33Ptr[xyzIndex] = ReM1*(2.0*lmuPtr[xyzIndex]*dwdz +
279  llambdaPtr[xyzIndex]*bulkTau);
280 
281  ltau12Ptr[xyzIndex] = ReM1*lmuPtr[xyzIndex]*(dudy + dvdx);
282  ltau13Ptr[xyzIndex] = ReM1*lmuPtr[xyzIndex]*(dudz + dwdx);
283  ltau23Ptr[xyzIndex] = ReM1*lmuPtr[xyzIndex]*(dvdz + dwdy);
284 
285  }
286  }
287  }
288  }
289  } else if (gridType == simulation::grid::RECTILINEAR || numDim == 1) {
290 
291  if(numDim == 1){
292 
293  size_t xStart = regionInterval[0].first;
294  size_t xEnd = regionInterval[0].second;
295  size_t xSize = bufferSize[0];
296  size_t nPoints = xSize;
297 
298  double* ltau11Ptr = tauBuffers[0];
299  double* velGradData0 = velGradBuffers[0];
300  const double *ldudxiPtr = &velGradData0[0];
301 
302  for(size_t iX = xStart;iX <= xEnd;iX++){
303  size_t xyzIndex = iX;
304  double dxidx = gridMetrics[xyzIndex];
305  double dudx = ldudxiPtr[xyzIndex]*dxidx;
306  double bulkTau = dudx;
307  double jacobianDeterminant = gridJacobian[xyzIndex];
308 
309  ltau11Ptr[xyzIndex] = gridJacobian[xyzIndex]*ReM1*(2.0*lmuPtr[xyzIndex]*dudx +
310  llambdaPtr[xyzIndex]*bulkTau);
311 
312  }
313  } else if (numDim == 2){
314 
315  size_t xStart = regionInterval[0].first;
316  size_t xEnd = regionInterval[0].second;
317  size_t yStart = regionInterval[1].first;
318  size_t yEnd = regionInterval[1].second;
319 
320  size_t xSize = bufferSize[0];
321  size_t ySize = bufferSize[1];
322  size_t nPoints = xSize*ySize;
323 
324  // brute force here, probably can be refactored to be more elegant
325  double* ltau11Ptr = tauBuffers[0];
326  double* ltau12Ptr = tauBuffers[1];
327  double* ltau22Ptr = tauBuffers[3];
328 
329  double* velGradData0 = velGradBuffers[0];
330  double* velGradData1 = velGradBuffers[1];
331 
332  const double *ldudxiPtr = &velGradData0[0];
333  const double *ldvdxiPtr = &velGradData0[nPoints];
334  const double *ldudetaPtr = &velGradData1[0];
335  const double *ldvdetaPtr = &velGradData1[nPoints];
336 
337  for(size_t iY = yStart;iY <= yEnd;iY++){
338  size_t yzIndex = iY*xSize;
339  for(size_t iX = xStart;iX <= xEnd;iX++){
340  size_t xyzIndex = yzIndex + iX;
341 
342  double dxidx = gridMetrics[xyzIndex];
343  double detady = gridMetrics[nPoints+xyzIndex];
344 
345 
346  double dudx = ldudxiPtr[xyzIndex]*dxidx;
347  double dudy = ldudetaPtr[xyzIndex]*detady;
348  double dvdx = ldvdxiPtr[xyzIndex]*dxidx;
349  double dvdy = ldvdetaPtr[xyzIndex]*detady;
350 
351  double bulkTau = dudx + dvdy;
352 
353  double fac = gridJacobian[xyzIndex]*ReM1;
354 
355  ltau11Ptr[xyzIndex] = fac*(2.0*lmuPtr[xyzIndex]*dudx +
356  llambdaPtr[xyzIndex]*bulkTau);
357  ltau22Ptr[xyzIndex] = fac*(2.0*lmuPtr[xyzIndex]*dvdy +
358  llambdaPtr[xyzIndex]*bulkTau);
359 
360  ltau12Ptr[xyzIndex] = fac*lmuPtr[xyzIndex]*(dudy + dvdx);
361 
362  }
363  }
364  } else if(numDim == 3) {
365 
366  size_t xSize = bufferSize[0];
367  size_t ySize = bufferSize[1];
368  size_t zSize = bufferSize[2];
369 
370  size_t xStart = regionInterval[0].first;
371  size_t xEnd = regionInterval[0].second;
372  size_t yStart = regionInterval[1].first;
373  size_t yEnd = regionInterval[1].second;
374  size_t zStart = regionInterval[2].first;
375  size_t zEnd = regionInterval[2].second;
376 
377  size_t nPlane = xSize*ySize;
378  size_t nPoints = xSize*ySize*zSize;
379 
380  // brute force here, probably can be refactored to be more elegant
381  double* ltau11Ptr = tauBuffers[0];
382  double* ltau12Ptr = tauBuffers[1];
383  double* ltau13Ptr = tauBuffers[2];
384  double* ltau22Ptr = tauBuffers[4];
385  double* ltau23Ptr = tauBuffers[5];
386  double* ltau33Ptr = tauBuffers[8];
387 
388  double* velGradData0 = velGradBuffers[0];
389  double* velGradData1 = velGradBuffers[1];
390  double* velGradData2 = velGradBuffers[2];
391 
392  const double* ldudxiPtr = &velGradData0[0];
393  const double* ldvdxiPtr = &velGradData0[nPoints];
394  const double* ldwdxiPtr = &velGradData0[2*nPoints];
395  const double* ldudetaPtr = &velGradData1[0];
396  const double* ldvdetaPtr = &velGradData1[nPoints];
397  const double* ldwdetaPtr = &velGradData1[2*nPoints];
398  const double* ldudzetaPtr = &velGradData2[0];
399  const double* ldvdzetaPtr = &velGradData2[nPoints];
400  const double* ldwdzetaPtr = &velGradData2[2*nPoints];
401 
402  for(size_t iZ = zStart;iZ <= zEnd;iZ++){
403  size_t zIndex = iZ*nPlane;
404  for(size_t iY = yStart;iY <= yEnd;iY++){
405  size_t yzIndex = iY*xSize+zIndex;
406  for(size_t iX = xStart;iX <= xEnd;iX++){
407  size_t xyzIndex = yzIndex + iX;
408 
409  double dxidx = gridMetrics[xyzIndex];
410  double detady = gridMetrics[nPoints+xyzIndex];
411  double dzetadz = gridMetrics[2*nPoints+xyzIndex];
412 
413 
414  double dudx = ldudxiPtr[xyzIndex]*dxidx;
415  double dudy = ldudetaPtr[xyzIndex]*detady;
416  double dudz = ldudzetaPtr[xyzIndex]*dzetadz;
417 
418  double dvdx = ldvdxiPtr[xyzIndex]*dxidx;
419  double dvdy = ldvdetaPtr[xyzIndex]*detady;
420  double dvdz = ldvdzetaPtr[xyzIndex]*dzetadz;
421 
422  double dwdx = ldwdxiPtr[xyzIndex]*dxidx;
423  double dwdy = ldwdetaPtr[xyzIndex]*detady;
424  double dwdz = ldwdzetaPtr[xyzIndex]*dzetadz;
425 
426  double bulkTau = dudx + dvdy + dwdz;
427 
428  double fac = ReM1*gridJacobian[xyzIndex];
429 
430  ltau11Ptr[xyzIndex] = fac*(2.0*lmuPtr[xyzIndex]*dudx +
431  llambdaPtr[xyzIndex]*bulkTau);
432  ltau22Ptr[xyzIndex] = fac*(2.0*lmuPtr[xyzIndex]*dvdy +
433  llambdaPtr[xyzIndex]*bulkTau);
434  ltau33Ptr[xyzIndex] = fac*(2.0*lmuPtr[xyzIndex]*dwdz +
435  llambdaPtr[xyzIndex]*bulkTau);
436 
437  ltau12Ptr[xyzIndex] = fac*lmuPtr[xyzIndex]*(dudy + dvdx);
438  ltau13Ptr[xyzIndex] = fac*lmuPtr[xyzIndex]*(dudz + dwdx);
439  ltau23Ptr[xyzIndex] = fac*lmuPtr[xyzIndex]*(dvdz + dwdy);
440 
441  }
442  }
443  }
444  }
445  } else { // must be curvilinear
446 
447  if (numDim == 2){
448 
449  size_t xStart = regionInterval[0].first;
450  size_t xEnd = regionInterval[0].second;
451  size_t yStart = regionInterval[1].first;
452  size_t yEnd = regionInterval[1].second;
453 
454  size_t xSize = bufferSize[0];
455  size_t ySize = bufferSize[1];
456  size_t nPoints = xSize*ySize;
457 
458  // brute force here, probably can be refactored to be more elegant
459  double* ltau11Ptr = tauBuffers[0];
460  double* ltau12Ptr = tauBuffers[1];
461  double* ltau22Ptr = tauBuffers[3];
462 
463  double* velGradData0 = velGradBuffers[0];
464  double* velGradData1 = velGradBuffers[1];
465 
466  const double *ldudxiPtr = &velGradData0[0];
467  const double *ldvdxiPtr = &velGradData0[nPoints];
468  const double *ldudetaPtr = &velGradData1[0];
469  const double *ldvdetaPtr = &velGradData1[nPoints];
470 
471  const size_t offset1 = nPoints;
472  const size_t offset2 = 2*nPoints;
473  const size_t offset3 = 3*nPoints;
474 
475  for(size_t iY = yStart;iY <= yEnd;iY++){
476  size_t yzIndex = iY*xSize;
477  for(size_t iX = xStart;iX <= xEnd;iX++){
478  size_t xyzIndex = yzIndex + iX;
479 
480  double dxidx = gridMetrics[xyzIndex];
481  double dxidy = gridMetrics[offset1+xyzIndex];
482 
483  double detadx = gridMetrics[offset2+xyzIndex];
484  double detady = gridMetrics[offset3+xyzIndex];
485 
486  double dudx = dxidx*ldudxiPtr[xyzIndex] + detadx*ldudetaPtr[xyzIndex];
487  double dudy = dxidy*ldudxiPtr[xyzIndex] + detady*ldudetaPtr[xyzIndex];
488 
489  double dvdx = dxidx*ldvdxiPtr[xyzIndex] + detadx*ldvdetaPtr[xyzIndex];
490  double dvdy = dxidy*ldvdxiPtr[xyzIndex] + detady*ldvdetaPtr[xyzIndex];
491 
492  double bulkTau = dudx + dvdy;
493 
494  double fac = ReM1*gridJacobian[xyzIndex];
495 
496  ltau11Ptr[xyzIndex] = fac*(2.0*lmuPtr[xyzIndex]*dudx + llambdaPtr[xyzIndex]*bulkTau);
497  ltau22Ptr[xyzIndex] = fac*(2.0*lmuPtr[xyzIndex]*dvdy + llambdaPtr[xyzIndex]*bulkTau);
498  ltau12Ptr[xyzIndex] = fac*lmuPtr[xyzIndex]*(dudy+dvdx);
499 
500  }
501  }
502  } else if(numDim == 3) {
503 
504  size_t xSize = bufferSize[0];
505  size_t ySize = bufferSize[1];
506  size_t zSize = bufferSize[2];
507 
508  size_t xStart = regionInterval[0].first;
509  size_t xEnd = regionInterval[0].second;
510  size_t yStart = regionInterval[1].first;
511  size_t yEnd = regionInterval[1].second;
512  size_t zStart = regionInterval[2].first;
513  size_t zEnd = regionInterval[2].second;
514 
515  size_t nPlane = xSize*ySize;
516  size_t nPoints = xSize*ySize*zSize;
517 
518  size_t offset1 = nPoints;
519  size_t offset2 = 2*nPoints;
520  size_t offset3 = 3*nPoints;
521  size_t offset4 = 4*nPoints;
522  size_t offset5 = 5*nPoints;
523  size_t offset6 = 6*nPoints;
524  size_t offset7 = 7*nPoints;
525  size_t offset8 = 8*nPoints;
526 
527  // brute force here, probably can be refactored to be more elegant
528  double* ltau11Ptr = tauBuffers[0];
529  double* ltau12Ptr = tauBuffers[1];
530  double* ltau13Ptr = tauBuffers[2];
531  double* ltau22Ptr = tauBuffers[4];
532  double* ltau23Ptr = tauBuffers[5];
533  double* ltau33Ptr = tauBuffers[8];
534 
535  double* velGradData0 = velGradBuffers[0];
536  double* velGradData1 = velGradBuffers[1];
537  double* velGradData2 = velGradBuffers[2];
538 
539  const double* ldudxiPtr = &velGradData0[0];
540  const double* ldvdxiPtr = &velGradData0[offset1];
541  const double* ldwdxiPtr = &velGradData0[offset2];
542  const double* ldudetaPtr = &velGradData1[0];
543  const double* ldvdetaPtr = &velGradData1[offset1];
544  const double* ldwdetaPtr = &velGradData1[offset2];
545  const double* ldudzetaPtr = &velGradData2[0];
546  const double* ldvdzetaPtr = &velGradData2[offset1];
547  const double* ldwdzetaPtr = &velGradData2[offset2];
548 
549  for(size_t iZ = zStart;iZ <= zEnd;iZ++){
550  size_t zIndex = iZ*nPlane;
551  for(size_t iY = yStart;iY <= yEnd;iY++){
552  size_t yzIndex = iY*xSize+zIndex;
553  for(size_t iX = xStart;iX <= xEnd;iX++){
554  size_t xyzIndex = yzIndex + iX;
555 
556  double dxidx = gridMetrics[xyzIndex];
557  double dxidy = gridMetrics[offset1+xyzIndex];
558  double dxidz = gridMetrics[offset2+xyzIndex];
559 
560  double detadx = gridMetrics[offset3+xyzIndex];
561  double detady = gridMetrics[offset4+xyzIndex];
562  double detadz = gridMetrics[offset5+xyzIndex];
563 
564  double dzetadx = gridMetrics[offset6+xyzIndex];
565  double dzetady = gridMetrics[offset7+xyzIndex];
566  double dzetadz = gridMetrics[offset8+xyzIndex];
567 
568  double jacobianDeterminant = gridJacobian[xyzIndex];
569 
570  double dudx = (dxidx*ldudxiPtr[xyzIndex] +
571  detadx*ldudetaPtr[xyzIndex] +
572  dzetadx*ldudzetaPtr[xyzIndex]);
573  double dudy = (dxidy*ldudxiPtr[xyzIndex] +
574  detady*ldudetaPtr[xyzIndex] +
575  dzetady*ldudzetaPtr[xyzIndex]);
576  double dudz = (dxidz*ldudxiPtr[xyzIndex] +
577  detadz*ldudetaPtr[xyzIndex] +
578  dzetadz*ldudzetaPtr[xyzIndex]);
579 
580  double dvdx = (dxidx*ldvdxiPtr[xyzIndex] +
581  detadx*ldvdetaPtr[xyzIndex] +
582  dzetadx*ldvdzetaPtr[xyzIndex]);
583  double dvdy = (dxidy*ldvdxiPtr[xyzIndex] +
584  detady*ldvdetaPtr[xyzIndex] +
585  dzetady*ldvdzetaPtr[xyzIndex]);
586  double dvdz = (dxidz*ldvdxiPtr[xyzIndex] +
587  detadz*ldvdetaPtr[xyzIndex] +
588  dzetadz*ldvdzetaPtr[xyzIndex]);
589 
590  double dwdx = (dxidx*ldwdxiPtr[xyzIndex] +
591  detadx*ldwdetaPtr[xyzIndex] +
592  dzetadx*ldwdzetaPtr[xyzIndex]);
593  double dwdy = (dxidy*ldwdxiPtr[xyzIndex] +
594  detady*ldwdetaPtr[xyzIndex] +
595  dzetady*ldwdzetaPtr[xyzIndex]);
596  double dwdz = (dxidz*ldwdxiPtr[xyzIndex] +
597  detadz*ldwdetaPtr[xyzIndex] +
598  dzetadz*ldwdzetaPtr[xyzIndex]);
599 
600  double bulkTau = dudx + dvdy + dwdz;
601 
602  double fac = gridJacobian[xyzIndex] * ReM1;
603 
604  ltau11Ptr[xyzIndex] = fac*(2.0*lmuPtr[xyzIndex]*dudx +
605  llambdaPtr[xyzIndex]*bulkTau);
606 
607  ltau22Ptr[xyzIndex] = fac*(2.0*lmuPtr[xyzIndex]*dvdy +
608  llambdaPtr[xyzIndex]*bulkTau);
609 
610  ltau33Ptr[xyzIndex] = fac*(2.0*lmuPtr[xyzIndex]*dwdz +
611  llambdaPtr[xyzIndex]*bulkTau);
612 
613  ltau12Ptr[xyzIndex] = fac*lmuPtr[xyzIndex]*(dudy + dvdx);
614  ltau13Ptr[xyzIndex] = fac*lmuPtr[xyzIndex]*(dudz + dwdx);
615  ltau23Ptr[xyzIndex] = fac*lmuPtr[xyzIndex]*(dvdz + dwdy);
616 
617  }
618  }
619  }
620  }
621 
622  }
623  return(0);
624  };
625 
627  int ComputeHeatFluxBuffer(const pcpp::IndexIntervalType &regionInterval,
628  const std::vector<size_t> &bufferSize,
629  int gridType,
630  const std::vector<double> &gridMetrics,
631  const std::vector<double> &gridJacobian,
632  const double &Re,
633  const std::vector<double *> &tvBuffers,
634  const std::vector<double *> &temperatureGradBuffers,
635  std::vector<double *> &heatFluxBuffers)
636  {
637  int numDim = bufferSize.size();
638 
639  if(numDim < 1)
640  return(1);
641 
642  double *lmuPtr = tvBuffers[0];
643  double *llambdaPtr = tvBuffers[1];
644  double *lkappaPtr = tvBuffers[2];
645 
646  double ReM1=1.0/Re;
647 
648  if(gridType < simulation::grid::RECTILINEAR){
649 
650  double jacobianDeterminant = gridJacobian[0];
651 
652  if(numDim == 1){
653 
654  const double *gradTxi = temperatureGradBuffers[0];
655 
656  double dXidX = gridMetrics[0];
657 
658  ReM1 *= jacobianDeterminant*dXidX;
659 
660  size_t xStart = regionInterval[0].first;
661  size_t xEnd = regionInterval[0].second;
662  size_t xSize = bufferSize[0];
663  size_t nPoints = xSize;
664 
665  for(size_t iX = xStart;iX <= xEnd;iX++){
666  size_t xyzIndex = iX;
667  heatFluxBuffers[0][xyzIndex] = lkappaPtr[xyzIndex]*ReM1*gradTxi[xyzIndex];
668  }
669 
670  } else if (numDim == 2){
671 
672  double dXidX = gridMetrics[0];
673  double dEtadY = gridMetrics[1];
674 
675  const double *gradTxi = temperatureGradBuffers[0];
676  const double *gradTeta = temperatureGradBuffers[1];
677 
678  ReM1 *= jacobianDeterminant;
679 
680  size_t xStart = regionInterval[0].first;
681  size_t xEnd = regionInterval[0].second;
682  size_t yStart = regionInterval[1].first;
683  size_t yEnd = regionInterval[1].second;
684 
685  size_t xSize = bufferSize[0];
686  size_t ySize = bufferSize[1];
687  size_t nPoints = xSize*ySize;
688 
689  for(size_t iY = yStart;iY <= yEnd;iY++){
690  size_t yzIndex = iY*xSize;
691  for(size_t iX = xStart;iX <= xEnd;iX++){
692  size_t xyzIndex = yzIndex + iX;
693  double fac = lkappaPtr[xyzIndex]*ReM1;
694 
695  heatFluxBuffers[0][xyzIndex] = fac*dXidX*gradTxi[xyzIndex];
696  heatFluxBuffers[1][xyzIndex] = fac*dEtadY*gradTeta[xyzIndex];
697 
698  }
699  }
700  } else if(numDim == 3) {
701 
702  double dXidX = gridMetrics[0];
703  double dEtadY = gridMetrics[1];
704  double dZetadZ = gridMetrics[2];
705 
706  const double *gradTxi = temperatureGradBuffers[0];
707  const double *gradTeta = temperatureGradBuffers[1];
708  const double *gradTzeta = temperatureGradBuffers[2];
709 
710  ReM1 *= jacobianDeterminant;
711 
712  size_t xSize = bufferSize[0];
713  size_t ySize = bufferSize[1];
714  size_t zSize = bufferSize[2];
715 
716  size_t xStart = regionInterval[0].first;
717  size_t xEnd = regionInterval[0].second;
718  size_t yStart = regionInterval[1].first;
719  size_t yEnd = regionInterval[1].second;
720  size_t zStart = regionInterval[2].first;
721  size_t zEnd = regionInterval[2].second;
722 
723  size_t nPlane = xSize*ySize;
724  size_t nPoints = xSize*ySize*zSize;
725 
726  for(size_t iZ = zStart;iZ <= zEnd;iZ++){
727  size_t zIndex = iZ*nPlane;
728  for(size_t iY = yStart;iY <= yEnd;iY++){
729  size_t yzIndex = iY*xSize+zIndex;
730  for(size_t iX = xStart;iX <= xEnd;iX++){
731  size_t xyzIndex = yzIndex + iX;
732  double fac = lkappaPtr[xyzIndex]*ReM1;
733  heatFluxBuffers[0][xyzIndex] = fac*dXidX*gradTxi[xyzIndex];
734  heatFluxBuffers[1][xyzIndex] = fac*dEtadY*gradTeta[xyzIndex];
735  heatFluxBuffers[2][xyzIndex] = fac*dZetadZ*gradTzeta[xyzIndex];
736 
737  }
738  }
739  }
740  }
741  } else if (gridType == simulation::grid::RECTILINEAR || (numDim == 1)) {
742 
743  // numDim == 1 is the same in rectilinear or curvilinear
744  if(numDim == 1){
745 
746  const double *gradTxi = temperatureGradBuffers[0];
747 
748  size_t xStart = regionInterval[0].first;
749  size_t xEnd = regionInterval[0].second;
750  size_t xSize = bufferSize[0];
751  size_t nPoints = xSize;
752 
753 
754  for(size_t iX = xStart;iX <= xEnd;iX++){
755  size_t xyzIndex = iX;
756  heatFluxBuffers[0][xyzIndex] = lkappaPtr[xyzIndex]*ReM1*
757  gridJacobian[xyzIndex]*gridMetrics[xyzIndex]*gradTxi[xyzIndex];
758  }
759 
760  } else if (numDim == 2){
761 
762  const double *gradTxi = temperatureGradBuffers[0];
763  const double *gradTeta = temperatureGradBuffers[1];
764 
765  size_t xStart = regionInterval[0].first;
766  size_t xEnd = regionInterval[0].second;
767  size_t yStart = regionInterval[1].first;
768  size_t yEnd = regionInterval[1].second;
769 
770  size_t xSize = bufferSize[0];
771  size_t ySize = bufferSize[1];
772  size_t nPoints = xSize*ySize;
773 
774  for(size_t iY = yStart;iY <= yEnd;iY++){
775  size_t yzIndex = iY*xSize;
776  for(size_t iX = xStart;iX <= xEnd;iX++){
777  size_t xyzIndex = yzIndex + iX;
778 
779  double dXidX = gridMetrics[xyzIndex];
780  double dEtadY = gridMetrics[xyzIndex+nPoints];
781 
782  double fac = gridJacobian[xyzIndex]*lkappaPtr[xyzIndex]*ReM1;
783 
784  heatFluxBuffers[0][xyzIndex] = fac*dXidX*gradTxi[xyzIndex];
785  heatFluxBuffers[1][xyzIndex] = fac*dEtadY*gradTeta[xyzIndex];
786 
787  }
788  }
789  } else if(numDim == 3) {
790 
791  const double *gradTxi = temperatureGradBuffers[0];
792  const double *gradTeta = temperatureGradBuffers[1];
793  const double *gradTzeta = temperatureGradBuffers[2];
794 
795  size_t xSize = bufferSize[0];
796  size_t ySize = bufferSize[1];
797  size_t zSize = bufferSize[2];
798 
799  size_t xStart = regionInterval[0].first;
800  size_t xEnd = regionInterval[0].second;
801  size_t yStart = regionInterval[1].first;
802  size_t yEnd = regionInterval[1].second;
803  size_t zStart = regionInterval[2].first;
804  size_t zEnd = regionInterval[2].second;
805 
806  size_t nPlane = xSize*ySize;
807  size_t nPoints = xSize*ySize*zSize;
808 
809  for(size_t iZ = zStart;iZ <= zEnd;iZ++){
810  size_t zIndex = iZ*nPlane;
811  for(size_t iY = yStart;iY <= yEnd;iY++){
812  size_t yzIndex = iY*xSize+zIndex;
813  for(size_t iX = xStart;iX <= xEnd;iX++){
814  size_t xyzIndex = yzIndex + iX;
815 
816  double dXidX = gridMetrics[xyzIndex];
817  double dEtadY = gridMetrics[xyzIndex+nPoints];
818  double dZetadZ = gridMetrics[xyzIndex+2*nPoints];
819 
820  double fac = gridJacobian[xyzIndex]*lkappaPtr[xyzIndex]*ReM1;
821 
822  heatFluxBuffers[0][xyzIndex] = fac*dXidX*gradTxi[xyzIndex];
823  heatFluxBuffers[1][xyzIndex] = fac*dEtadY*gradTeta[xyzIndex];
824  heatFluxBuffers[2][xyzIndex] = fac*dZetadZ*gradTzeta[xyzIndex];
825 
826  }
827  }
828  }
829  }
830  } else { // must be curvilinear
831 
832  if (numDim == 2){
833 
834  const double *gradTxi = temperatureGradBuffers[0];
835  const double *gradTeta = temperatureGradBuffers[1];
836 
837  size_t xStart = regionInterval[0].first;
838  size_t xEnd = regionInterval[0].second;
839  size_t yStart = regionInterval[1].first;
840  size_t yEnd = regionInterval[1].second;
841 
842  size_t xSize = bufferSize[0];
843  size_t ySize = bufferSize[1];
844  size_t nPoints = xSize*ySize;
845 
846  const double *dXidX = &gridMetrics[0];
847  const double *dXidY = &gridMetrics[nPoints];
848  const double *dEtadX = &gridMetrics[2*nPoints];
849  const double *dEtadY = &gridMetrics[3*nPoints];
850 
851  for(size_t iY = yStart;iY <= yEnd;iY++){
852  size_t yzIndex = iY*xSize;
853  for(size_t iX = xStart;iX <= xEnd;iX++){
854  size_t xyzIndex = yzIndex + iX;
855 
856  double dxidx = dXidX[xyzIndex];
857  double dxidy = dXidY[xyzIndex];
858  double detadx = dEtadX[xyzIndex];
859  double detady = dEtadY[xyzIndex];
860 
861  double fac = lkappaPtr[xyzIndex]*ReM1*gridJacobian[xyzIndex];
862 
863  heatFluxBuffers[0][xyzIndex] = fac*(dxidx*gradTxi[xyzIndex] +
864  detadx*gradTeta[xyzIndex]);
865  heatFluxBuffers[1][xyzIndex] = fac*(detady*gradTeta[xyzIndex] +
866  dxidy*gradTxi[xyzIndex]);
867 
868  }
869  }
870  } else if(numDim == 3) {
871 
872  const double *gradTxi = temperatureGradBuffers[0];
873  const double *gradTeta = temperatureGradBuffers[1];
874  const double *gradTzeta = temperatureGradBuffers[2];
875 
876  size_t xSize = bufferSize[0];
877  size_t ySize = bufferSize[1];
878  size_t zSize = bufferSize[2];
879 
880  size_t xStart = regionInterval[0].first;
881  size_t xEnd = regionInterval[0].second;
882  size_t yStart = regionInterval[1].first;
883  size_t yEnd = regionInterval[1].second;
884  size_t zStart = regionInterval[2].first;
885  size_t zEnd = regionInterval[2].second;
886 
887  size_t nPlane = xSize*ySize;
888  size_t nPoints = xSize*ySize*zSize;
889 
890  const double *dXidX = &gridMetrics[0];
891  const double *dXidY = &gridMetrics[nPoints];
892  const double *dXidZ = &gridMetrics[2*nPoints];
893  const double *dEtadX = &gridMetrics[3*nPoints];
894  const double *dEtadY = &gridMetrics[4*nPoints];
895  const double *dEtadZ = &gridMetrics[5*nPoints];
896  const double *dZetadX = &gridMetrics[6*nPoints];
897  const double *dZetadY = &gridMetrics[7*nPoints];
898  const double *dZetadZ = &gridMetrics[8*nPoints];
899 
900  for(size_t iZ = zStart;iZ <= zEnd;iZ++){
901  size_t zIndex = iZ*nPlane;
902  for(size_t iY = yStart;iY <= yEnd;iY++){
903  size_t yzIndex = iY*xSize+zIndex;
904  for(size_t iX = xStart;iX <= xEnd;iX++){
905  size_t xyzIndex = yzIndex + iX;
906 
907 
908  double dxidx = dXidX[xyzIndex];
909  double dxidy = dXidY[xyzIndex];
910  double dxidz = dXidZ[xyzIndex];
911 
912  double detadx = dEtadX[xyzIndex];
913  double detady = dEtadY[xyzIndex];
914  double detadz = dEtadZ[xyzIndex];
915 
916  double dzetadx = dZetadX[xyzIndex];
917  double dzetady = dZetadY[xyzIndex];
918  double dzetadz = dZetadZ[xyzIndex];
919 
920  double fac = lkappaPtr[xyzIndex]*ReM1*gridJacobian[xyzIndex];
921 
922 
923  heatFluxBuffers[0][xyzIndex] = fac*(dxidx*gradTxi[xyzIndex] +
924  detadx*gradTeta[xyzIndex] +
925  dzetadx*gradTzeta[xyzIndex]);
926  heatFluxBuffers[1][xyzIndex] = fac*(dxidy*gradTxi[xyzIndex] +
927  detady*gradTeta[xyzIndex] +
928  dzetady*gradTzeta[xyzIndex]);
929  heatFluxBuffers[2][xyzIndex] = fac*(dxidz*gradTxi[xyzIndex] +
930  detadz*gradTeta[xyzIndex] +
931  dzetadz*gradTzeta[xyzIndex]);
932 
933  }
934  }
935  }
936  }
937  }
938  return(0);
939  };
940 
941  } // namespace util
942 } // namespace viscid
viscid::util::ComputeTauBuffer
int ComputeTauBuffer(const pcpp::IndexIntervalType &regionInterval, const std::vector< size_t > &bufferSize, int gridType, const std::vector< double > &gridMetrics, const std::vector< double > &gridJacobian, const double &Re, const std::vector< double *> &tvBuffers, const std::vector< double *> &velGradBuffers, std::vector< double *> &tauBuffers)
Definition: ViscidUtil.C:113
bufferSize
void size_t int size_t int size_t int int int int double int int double double *void size_t int size_t int int int int int double int size_t size_t size_t double double *void size_t int size_t int size_t size_t int double int double double *void size_t size_t * bufferSize
Definition: OperatorKernels.H:29
viscid::util::ComputeHeatFluxBuffer
int ComputeHeatFluxBuffer(const pcpp::IndexIntervalType &regionInterval, const std::vector< size_t > &bufferSize, int gridType, const std::vector< double > &gridMetrics, const std::vector< double > &gridJacobian, const double &Re, const std::vector< double *> &tvBuffers, const std::vector< double *> &temperatureGradBuffers, std::vector< double *> &heatFluxBuffers)
Definition: ViscidUtil.C:627
gridType
void const size_t const size_t const size_t const int const int * gridType
Definition: EulerKernels.H:10
jacobianDeterminant
void size_t size_t int size_t size_t size_t size_t int double double * jacobianDeterminant
Definition: SATKernels.H:21
viscid::util::ComputeTVBuffer
int ComputeTVBuffer(const int *numDimPtr, const size_t *bufferSize, const size_t *numPointsPtr, const size_t *bufferInterval, const double *rhoBuffer, const double *rhoVBuffer, const double *rhoEBuffer, double *pressureBuffer, double *tempBuffer, double *rhom1Buffer, double *velBuffer)
Definition: ViscidUtil.C:97
viscid
Definition: Viscid.f90:1
viscid::util::ComputeTVBufferPower
int ComputeTVBufferPower(const pcpp::IndexIntervalType &regionInterval, const std::vector< size_t > &bufferSize, const double *temperatureBuffer, std::vector< double *> &tvBuffer, const double beta, const double power, const double bulkViscFac, const double specificHeat, const double prandtlNumber)
Compute transport coefficients using the power law.
Definition: ViscidUtil.C:14
gridJacobian
void const size_t const size_t const size_t const int const double * gridJacobian
Definition: MetricKernels.H:9
bufferInterval
void size_t int size_t int size_t int int int int double int int double double *void size_t int size_t int int int int int double int size_t size_t size_t double double *void size_t int size_t int size_t size_t int double int double double *void size_t size_t size_t * bufferInterval
Definition: OperatorKernels.H:29
ix::util::sizeextent
Simple Block Structured Mesh object.
Definition: IndexUtil.H:21