PlasCom2  1.0
XPACC Multi-physics simluation application
EulerUtil.C
Go to the documentation of this file.
1 #include "PCPPTypes.H"
2 #include "EulerUtil.H"
3 
4 namespace euler {
5  namespace util {
6 
7 /*
8  // MJA, to be deleted
12  int ComputeDVBufferIdeal(const pcpp::IndexIntervalType &regionInterval,
13  const std::vector<size_t> &bufferSizes,
14  const std::vector<double *> &stateBuffers,
15  std::vector<double *> &dvBuffers,
16  const double gamma,
17  const double specGasConst)
18  {
19  double gammam1 = gamma-1.0;
20  double invR = 1.0/specGasConst;
21 
22  int numDim = bufferSizes.size();
23 
24  if(numDim < 1)
25  return(1);
26 
27  double *lpressurePtr = dvBuffers[0];
28  double *ltemperaturePtr = dvBuffers[1];
29  double *lrhom1Ptr = dvBuffers[2];
30  double *lV1Ptr = dvBuffers[3];
31 
32 
33  const double *lrhoPtr = stateBuffers[0];
34  const double *lrhoVPtr = stateBuffers[1];
35  const double *lrhoEPtr = stateBuffers[numDim+1];
36 
37  if(numDim == 1){
38 
39  size_t xStart = regionInterval[0].first;
40  size_t xEnd = regionInterval[0].second;
41  size_t xSize = bufferSizes[0];
42  size_t nPoints = xSize;
43 
44  for(size_t iX = xStart;iX <= xEnd;iX++){
45 
46  size_t xyzIndex = iX;
47 
48  double rho = lrhoPtr[xyzIndex];
49  double rhom1 = 1.0/rho;
50  double rhoV1 = lrhoVPtr[xyzIndex];
51  double kineticTerm = 0.0;
52 
53  lrhom1Ptr[xyzIndex] = rhom1;
54  lV1Ptr[xyzIndex] = rhoV1*rhom1;
55  kineticTerm += rhoV1*rhoV1;
56 
57  lpressurePtr[xyzIndex] = gammam1*(lrhoEPtr[xyzIndex] - 0.5*kineticTerm*rhom1);
58  ltemperaturePtr[xyzIndex] = lpressurePtr[xyzIndex]*rhom1*invR;
59 
60  }
61  } else if (numDim == 2){
62 
63  double *lV2Ptr = dvBuffers[4];
64 
65  size_t xStart = regionInterval[0].first;
66  size_t xEnd = regionInterval[0].second;
67  size_t yStart = regionInterval[1].first;
68  size_t yEnd = regionInterval[1].second;
69 
70  size_t xSize = bufferSizes[0];
71  size_t ySize = bufferSizes[1];
72  size_t nPoints = xSize*ySize;
73 
74  for(size_t iY = yStart;iY <= yEnd;iY++){
75  size_t yzIndex = iY*xSize;
76  for(size_t iX = xStart;iX <= xEnd;iX++){
77 
78  size_t xyzIndex = yzIndex + iX;
79  size_t xyzIndex2 = xyzIndex + nPoints;
80 
81  double rho = lrhoPtr[xyzIndex];
82  double rhom1 = 1.0/rho;
83  double rhoV1 = lrhoVPtr[xyzIndex];
84  double rhoV2 = lrhoVPtr[xyzIndex2];
85 
86  double kineticTerm = 0.0;
87 
88  lrhom1Ptr[xyzIndex] = rhom1;
89  lV1Ptr[xyzIndex] = rhoV1*rhom1;
90  kineticTerm += rhoV1*rhoV1;
91  lV2Ptr[xyzIndex] = rhoV2*rhom1;
92  kineticTerm += rhoV2*rhoV2;
93 
94  lpressurePtr[xyzIndex] = gammam1*(lrhoEPtr[xyzIndex] - 0.5*kineticTerm*rhom1);
95  ltemperaturePtr[xyzIndex] = lpressurePtr[xyzIndex]*rhom1*invR;
96  }
97  }
98  } else if(numDim == 3) {
99 
100  double *lV2Ptr = dvBuffers[4];
101  double *lV3Ptr = dvBuffers[5];
102 
103  size_t xSize = bufferSizes[0];
104  size_t ySize = bufferSizes[1];
105  size_t zSize = bufferSizes[2];
106 
107  size_t xStart = regionInterval[0].first;
108  size_t xEnd = regionInterval[0].second;
109  size_t yStart = regionInterval[1].first;
110  size_t yEnd = regionInterval[1].second;
111  size_t zStart = regionInterval[2].first;
112  size_t zEnd = regionInterval[2].second;
113 
114  size_t nPlane = xSize*ySize;
115  size_t nPoints = xSize*ySize*zSize;
116 
117  for(size_t iZ = zStart;iZ <= zEnd;iZ++){
118  size_t zIndex = iZ*nPlane;
119  for(size_t iY = yStart;iY <= yEnd;iY++){
120  size_t yzIndex = iY*xSize+zIndex;
121  for(size_t iX = xStart;iX <= xEnd;iX++){
122 
123  size_t xyzIndex = yzIndex + iX;
124  size_t xyzIndex2 = xyzIndex + nPoints;
125  size_t xyzIndex3 = xyzIndex + 2*nPoints;
126 
127  double rho = lrhoPtr[xyzIndex];
128  double rhom1 = 1.0/rho;
129  double rhoV1 = lrhoVPtr[xyzIndex];
130  double rhoV2 = lrhoVPtr[xyzIndex2];
131  double rhoV3 = lrhoVPtr[xyzIndex3];
132  double kineticTerm = 0.0;
133 
134  lrhom1Ptr[xyzIndex] = rhom1;
135  lV1Ptr[xyzIndex] = rhoV1*rhom1;
136  kineticTerm += rhoV1*rhoV1;
137  lV2Ptr[xyzIndex] = rhoV2*rhom1;
138  kineticTerm += rhoV2*rhoV2;
139  lV3Ptr[xyzIndex] = rhoV3*rhom1;
140  kineticTerm += rhoV3*rhoV3;
141 
142  lpressurePtr[xyzIndex] = gammam1*(lrhoEPtr[xyzIndex] - 0.5*kineticTerm*rhom1);
143  ltemperaturePtr[xyzIndex] = lpressurePtr[xyzIndex]*rhom1*invR;
144  }
145  }
146  }
147  }
148  return(0);
149  };
150 */
151 
152 
153  int ComputeDVBuffer(const pcpp::IndexIntervalType &regionInterval,
154  const std::vector<size_t> &bufferSizes,
155  const std::vector<double *> &stateBuffers,
156  std::vector<double *> &dvBuffers)
157  {
158 
159  int numDim = bufferSizes.size();
160 
161  if(numDim < 1)
162  return(1);
163 
164  double *lpressurePtr = dvBuffers[0];
165  double *ltemperaturePtr = dvBuffers[1];
166  double *lrhom1Ptr = dvBuffers[2];
167  double *lV1Ptr = dvBuffers[3];
168  double *linternalEnergyPtr = dvBuffers[numDim+3];
169 
170  const double *lrhoPtr = stateBuffers[0];
171  const double *lrhoVPtr = stateBuffers[1];
172  const double *lrhoEPtr = stateBuffers[numDim+1];
173 
174  if(numDim == 1){
175 
176  size_t xStart = regionInterval[0].first;
177  size_t xEnd = regionInterval[0].second;
178  size_t xSize = bufferSizes[0];
179  size_t nPoints = xSize;
180 
181  for(size_t iX = xStart;iX <= xEnd;iX++){
182 
183  size_t xyzIndex = iX;
184 
185  double rho = lrhoPtr[xyzIndex];
186  double rhom1 = 1.0/rho;
187  double rhoV1 = lrhoVPtr[xyzIndex];
188  double rhoE = lrhoEPtr[xyzIndex];
189  double kineticTerm = 0.0;
190 
191 
192  lrhom1Ptr[xyzIndex] = rhom1;
193  lV1Ptr[xyzIndex] = rhoV1*rhom1;
194  kineticTerm += rhoV1*rhoV1;
195  linternalEnergyPtr[xyzIndex] = rhoE - 0.5*kineticTerm*rhom1;
196  }
197  } else if (numDim == 2){
198 
199  double *lV2Ptr = dvBuffers[4];
200 
201  size_t xStart = regionInterval[0].first;
202  size_t xEnd = regionInterval[0].second;
203  size_t yStart = regionInterval[1].first;
204  size_t yEnd = regionInterval[1].second;
205 
206  size_t xSize = bufferSizes[0];
207  size_t ySize = bufferSizes[1];
208  size_t nPoints = xSize*ySize;
209 
210  for(size_t iY = yStart;iY <= yEnd;iY++){
211  size_t yzIndex = iY*xSize;
212  for(size_t iX = xStart;iX <= xEnd;iX++){
213 
214  size_t xyzIndex = yzIndex + iX;
215  size_t xyzIndex2 = xyzIndex + nPoints;
216 
217  double rho = lrhoPtr[xyzIndex];
218  double rhom1 = 1.0/rho;
219  double rhoV1 = lrhoVPtr[xyzIndex];
220  double rhoV2 = lrhoVPtr[xyzIndex2];
221  double rhoE = lrhoEPtr[xyzIndex];
222 
223  double kineticTerm = 0.0;
224 
225  lrhom1Ptr[xyzIndex] = rhom1;
226  lV1Ptr[xyzIndex] = rhoV1*rhom1;
227  kineticTerm += rhoV1*rhoV1;
228  lV2Ptr[xyzIndex] = rhoV2*rhom1;
229  kineticTerm += rhoV2*rhoV2;
230  linternalEnergyPtr[xyzIndex] = rhoE - 0.5*kineticTerm*rhom1;
231  }
232  }
233  } else if(numDim == 3) {
234 
235  double *lV2Ptr = dvBuffers[4];
236  double *lV3Ptr = dvBuffers[5];
237 
238  size_t xSize = bufferSizes[0];
239  size_t ySize = bufferSizes[1];
240  size_t zSize = bufferSizes[2];
241 
242  size_t xStart = regionInterval[0].first;
243  size_t xEnd = regionInterval[0].second;
244  size_t yStart = regionInterval[1].first;
245  size_t yEnd = regionInterval[1].second;
246  size_t zStart = regionInterval[2].first;
247  size_t zEnd = regionInterval[2].second;
248 
249  size_t nPlane = xSize*ySize;
250  size_t nPoints = xSize*ySize*zSize;
251 
252  for(size_t iZ = zStart;iZ <= zEnd;iZ++){
253  size_t zIndex = iZ*nPlane;
254  for(size_t iY = yStart;iY <= yEnd;iY++){
255  size_t yzIndex = iY*xSize+zIndex;
256  for(size_t iX = xStart;iX <= xEnd;iX++){
257 
258  size_t xyzIndex = yzIndex + iX;
259  size_t xyzIndex2 = xyzIndex + nPoints;
260  size_t xyzIndex3 = xyzIndex + 2*nPoints;
261 
262  double rho = lrhoPtr[xyzIndex];
263  double rhom1 = 1.0/rho;
264  double rhoV1 = lrhoVPtr[xyzIndex];
265  double rhoV2 = lrhoVPtr[xyzIndex2];
266  double rhoV3 = lrhoVPtr[xyzIndex3];
267  double rhoE = lrhoEPtr[xyzIndex];
268  double kineticTerm = 0.0;
269 
270  lrhom1Ptr[xyzIndex] = rhom1;
271  lV1Ptr[xyzIndex] = rhoV1*rhom1;
272  kineticTerm += rhoV1*rhoV1;
273  lV2Ptr[xyzIndex] = rhoV2*rhom1;
274  kineticTerm += rhoV2*rhoV2;
275  lV3Ptr[xyzIndex] = rhoV3*rhom1;
276  kineticTerm += rhoV3*rhoV3;
277  linternalEnergyPtr[xyzIndex] = rhoE - 0.5*kineticTerm*rhom1;
278  }
279  }
280  }
281  }
282  return(0);
283  };
284 
285  int ComputeDVBuffer2(const int *numDimPtr,
286  const size_t *bufferSize,
287  const size_t *numPointsPtr,
288  const size_t *bufferInterval,
289  const double *lrhoPtr,
290  const double *lrhoVPtr,
291  const double *lrhoEPtr,
292  double *lpressurePtr,
293  double *ltemperaturePtr,
294  double *lrhom1Ptr,
295  double *lVPtr)
296  {
297 
298  double gamma = 1.4;
299  double gammam1 = gamma-1.0;
300  double gammaR = gamma/gammam1;
301 
302  int numDim = *numDimPtr;
303  size_t numPoints = *numPointsPtr;
304 
305  if(numDim < 1)
306  return(1);
307 
308  if(numDim == 1){
309 
310  size_t xStart = bufferInterval[0];
311  size_t xEnd = bufferInterval[1];
312 
313  size_t xSize = bufferSize[0];
314 
315  for(size_t iX = xStart;iX <= xEnd;iX++){
316 
317  size_t xyzIndex = iX;
318 
319  double rho = lrhoPtr[xyzIndex];
320  double rhom1 = 1.0/rho;
321  double rhoV = lrhoVPtr[xyzIndex];
322  double kineticTerm = 0.0;
323 
324  lrhom1Ptr[xyzIndex] = rhom1;
325  lVPtr[xyzIndex] = rhoV*rhom1;
326  kineticTerm += rhoV*rhoV;
327 
328  lpressurePtr[xyzIndex] = gammam1*(lrhoEPtr[xyzIndex]
329  - 0.5*kineticTerm*rhom1);
330  ltemperaturePtr[xyzIndex] = gammaR*lpressurePtr[xyzIndex]*rhom1;
331 
332  }
333  } else if (numDim == 2){
334 
335  size_t xStart = bufferInterval[0];
336  size_t xEnd = bufferInterval[1];
337  size_t yStart = bufferInterval[2];
338  size_t yEnd = bufferInterval[3];
339 
340  size_t xSize = bufferSize[0];
341  size_t ySize = bufferSize[1];
342 
343  for(size_t iY = yStart;iY <= yEnd;iY++){
344  size_t yzIndex = iY*xSize;
345  for(size_t iX = xStart;iX <= xEnd;iX++){
346 
347  size_t xyzIndex = yzIndex + iX;
348  size_t xyzIndex2 = xyzIndex + numPoints;
349 
350  double rho = lrhoPtr[xyzIndex];
351  double rhom1 = 1.0/rho;
352  double rhoV1 = lrhoVPtr[xyzIndex];
353  double rhoV2 = lrhoVPtr[xyzIndex2];
354  double kineticTerm = 0.0;
355 
356  lrhom1Ptr[xyzIndex] = rhom1;
357  lVPtr[xyzIndex] = rhoV1*rhom1;
358  kineticTerm += rhoV1*rhoV1;
359  lVPtr[xyzIndex2] = rhoV2*rhom1;
360  kineticTerm += rhoV2*rhoV2;
361 
362  lpressurePtr[xyzIndex] = gammam1*(lrhoEPtr[xyzIndex]
363  - 0.5*kineticTerm*rhom1);
364  ltemperaturePtr[xyzIndex] = gammaR*lpressurePtr[xyzIndex]*rhom1;
365 
366  }
367 
368  }
369 
370  } else if(numDim == 3) {
371 
372  size_t xSize = bufferSize[0];
373  size_t ySize = bufferSize[1];
374  size_t zSize = bufferSize[2];
375 
376  size_t xStart = bufferInterval[0];
377  size_t xEnd = bufferInterval[1];
378  size_t yStart = bufferInterval[2];
379  size_t yEnd = bufferInterval[3];
380  size_t zStart = bufferInterval[4];
381  size_t zEnd = bufferInterval[5];
382 
383  size_t nPlane = xSize*ySize;
384 
385  for(size_t iZ = zStart;iZ <= zEnd;iZ++){
386  size_t zIndex = iZ*nPlane;
387  for(size_t iY = yStart;iY <= yEnd;iY++){
388  size_t yzIndex = iY*xSize+zIndex;
389  for(size_t iX = xStart;iX <= xEnd;iX++){
390 
391  size_t xyzIndex = yzIndex + iX;
392  size_t xyzIndex2 = xyzIndex + numPoints;
393  size_t xyzIndex3 = xyzIndex + 2*numPoints;
394 
395  double rho = lrhoPtr[xyzIndex];
396  double rhom1 = 1.0/rho;
397  double rhoV1 = lrhoVPtr[xyzIndex];
398  double rhoV2 = lrhoVPtr[xyzIndex2];
399  double rhoV3 = lrhoVPtr[xyzIndex3];
400  double kineticTerm = 0.0;
401 
402  lrhom1Ptr[xyzIndex] = rhom1;
403  lVPtr[xyzIndex] = rhoV1*rhom1;
404  kineticTerm += rhoV1*rhoV1;
405  lVPtr[xyzIndex2] = rhoV2*rhom1;
406  kineticTerm += rhoV2*rhoV2;
407  lVPtr[xyzIndex3] = rhoV3*rhom1;
408  kineticTerm += rhoV3*rhoV3;
409 
410  lpressurePtr[xyzIndex] = gammam1*(lrhoEPtr[xyzIndex]
411  - 0.5*kineticTerm*rhom1);
412  ltemperaturePtr[xyzIndex] = gammaR*lpressurePtr[xyzIndex]*rhom1;
413 
414  }
415 
416  }
417 
418  }
419 
420  }
421  return(0);
422  };
423 
424 
425  }
426 }
euler::util::ComputeDVBuffer
int ComputeDVBuffer(const pcpp::IndexIntervalType &regionInterval, const std::vector< size_t > &bufferSizes, const std::vector< double *> &stateBuffers, std::vector< double *> &dvBuffers)
Definition: EulerUtil.C:153
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
numPoints
void const size_t * numPoints
Definition: EulerKernels.H:10
euler
Definition: Euler.f90:1
euler::util::ComputeDVBuffer2
int ComputeDVBuffer2(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: EulerUtil.C:285
bufferSizes
void const size_t const size_t * bufferSizes
Definition: MetricKernels.H:19
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