7 int numberOfVolumesPerAxisInPatch,
10 int auxiliaryVariables,
12 bool evaluateNonconservativeProduct,
14 bool evaluateMaximumEigenvalueAfterTimeStep,
16 NonconservativeProduct nonconservativeProduct,
18 MaxEigenvalue maxEigenvalue
26 constexpr int solversHaloSize = 1;
62 for (
int patchIndex=0; patchIndex<patchData.
numberOfCells; patchIndex++) {
66 for (
int x = 0; x < numberOfVolumesPerAxisInPatch; x++)
67 for (
int y = 0; y < numberOfVolumesPerAxisInPatch; y++)
68 for (
int z = 0; z < numberOfVolumesPerAxisInPatch; z++)
69 for (
int unknown=0; unknown<unknowns+auxiliaryVariables; unknown++) {
70 internal::copySolution_LoopBody(
71 patchData.
QIn[patchIndex],
76 patchData.
QOut[patchIndex],
83 for (
int x = -solversHaloSize; x < numberOfVolumesPerAxisInPatch+solversHaloSize; x++)
84 for (
int y = -solversHaloSize; y < numberOfVolumesPerAxisInPatch+solversHaloSize; y++)
85 for (
int z = -solversHaloSize; z < numberOfVolumesPerAxisInPatch+solversHaloSize; z++)
87 if ( ((x<0 or x>=numberOfVolumesPerAxisInPatch) and (y<0 or y>=numberOfVolumesPerAxisInPatch)) or
88 ((x<0 or x>=numberOfVolumesPerAxisInPatch) and (z<0 or z>=numberOfVolumesPerAxisInPatch)) or
89 ((y<0 or y>=numberOfVolumesPerAxisInPatch) and (z<0 or z>=numberOfVolumesPerAxisInPatch)) )
93 internal::computeTimeDerivative_LoopBody(
94 patchData.
QIn[patchIndex],
96 flux, nonconservativeProduct, source,
100 volumeIndex3d(x,y,z),
101 patchData.
t[patchIndex],
102 patchData.
dt[patchIndex],
111 for (
int x = -solversHaloSize; x < numberOfVolumesPerAxisInPatch+solversHaloSize; x++)
112 for (
int y = -solversHaloSize; y < numberOfVolumesPerAxisInPatch+solversHaloSize; y++)
113 for (
int z = -solversHaloSize; z < numberOfVolumesPerAxisInPatch+solversHaloSize; z++)
114 if ( ((x<0 or x>=numberOfVolumesPerAxisInPatch) and (y<0 or y>=numberOfVolumesPerAxisInPatch)) or
115 ((x<0 or x>=numberOfVolumesPerAxisInPatch) and (z<0 or z>=numberOfVolumesPerAxisInPatch)) or
116 ((y<0 or y>=numberOfVolumesPerAxisInPatch) and (z<0 or z>=numberOfVolumesPerAxisInPatch)) )
119 internal::computeQonFace_LoopBody(
120 patchData.
QIn[patchIndex],
125 volumeIndex3d(x,y,z),
126 patchData.
t[patchIndex],
127 patchData.
dt[patchIndex],
129 QfaceXneg, QfaceXpos,
130 QfaceYneg, QfaceYpos,
131 QfaceZneg, QfaceZpos,
138 for (
int x = -solversHaloSize; x < numberOfVolumesPerAxisInPatch+solversHaloSize; x++)
139 for (
int y = 0; y < numberOfVolumesPerAxisInPatch; y++)
140 for (
int z = 0; z < numberOfVolumesPerAxisInPatch; z++) {
141 internal::computeMaxEigenvalue_LoopBody(
142 patchData.
QIn[patchIndex],
148 volumeIndex3d(x,y,z),
149 patchData.
t[patchIndex], patchData.
dt[patchIndex],
155 for (
int x = 0; x < numberOfVolumesPerAxisInPatch; x++)
156 for (
int y = -solversHaloSize; y < numberOfVolumesPerAxisInPatch+solversHaloSize; y++)
157 for (
int z = 0; z < numberOfVolumesPerAxisInPatch; z++) {
158 internal::computeMaxEigenvalue_LoopBody(
159 patchData.
QIn[patchIndex],
165 volumeIndex3d(x,y,z),
166 patchData.
t[patchIndex], patchData.
dt[patchIndex],
172 for (
int x = 0; x < numberOfVolumesPerAxisInPatch; x++)
173 for (
int y = 0; y < numberOfVolumesPerAxisInPatch; y++)
174 for (
int z = -solversHaloSize; z < numberOfVolumesPerAxisInPatch+solversHaloSize; z++) {
175 internal::computeMaxEigenvalue_LoopBody(
176 patchData.
QIn[patchIndex],
182 volumeIndex3d(x,y,z),
183 patchData.
t[patchIndex], patchData.
dt[patchIndex],
195 for (
int x = 0; x < numberOfVolumesPerAxisInPatch; x++)
196 for (
int y = 0; y < numberOfVolumesPerAxisInPatch; y++)
197 for (
int z = 0; z < numberOfVolumesPerAxisInPatch; z++)
198 for (
int unknown = 0; unknown < unknowns; unknown++) {
199 internal::updateSolutionWithEigenvalueDamping_LoopBody(
200 patchData.
QIn[patchIndex],
205 eigenvalueEnumerator,
209 volumeIndex3d(x,y,z),
211 patchData.
dt[patchIndex],
212 patchData.
QOut[patchIndex],
214 QfaceXneg, QfaceXpos,
215 QfaceYneg, QfaceYpos,
216 QfaceZneg, QfaceZpos,
227 if (evaluateFlux and Dimensions==3) {
228 for (
int x = -solversHaloSize; x < numberOfVolumesPerAxisInPatch+solversHaloSize; x++)
229 for (
int y = 0; y < numberOfVolumesPerAxisInPatch; y++)
230 for (
int z = 0; z < numberOfVolumesPerAxisInPatch; z++) {
231 internal::computeFlux_LoopBody(
232 QfaceXneg, QfaceXpos, QInterEnumerator, QInEnumerator,
237 volumeIndex3d(x,y,z),
238 patchData.
t[patchIndex],
239 patchData.
dt[patchIndex],
241 tempFluxXL, tempFluxXR,
245 for (
int x = 0; x < numberOfVolumesPerAxisInPatch; x++)
246 for (
int y = -solversHaloSize; y < numberOfVolumesPerAxisInPatch+solversHaloSize; y++)
247 for (
int z = 0; z < numberOfVolumesPerAxisInPatch; z++) {
248 internal::computeFlux_LoopBody(
249 QfaceYneg, QfaceYpos, QInterEnumerator, QInEnumerator,
254 volumeIndex3d(x,y,z),
255 patchData.
t[patchIndex],
256 patchData.
dt[patchIndex],
258 tempFluxYL, tempFluxYR,
262 for (
int x = 0; x < numberOfVolumesPerAxisInPatch; x++)
263 for (
int y = 0; y < numberOfVolumesPerAxisInPatch; y++)
264 for (
int z = -solversHaloSize; z < numberOfVolumesPerAxisInPatch+solversHaloSize; z++) {
265 internal::computeFlux_LoopBody(
266 QfaceZneg, QfaceZpos, QInterEnumerator, QInEnumerator,
271 volumeIndex3d(x,y,z),
272 patchData.
t[patchIndex],
273 patchData.
dt[patchIndex],
275 tempFluxZL, tempFluxZR,
280 for (
int x = 0; x < numberOfVolumesPerAxisInPatch; x++)
281 for (
int y = 0; y < numberOfVolumesPerAxisInPatch; y++)
282 for (
int z = 0; z < numberOfVolumesPerAxisInPatch; z++)
283 for (
int unknown=0; unknown<unknowns; unknown++) {
284 internal::updateSolutionWithFlux_LoopBody(
285 tempFluxXL, tempFluxYL, tempFluxZL,
286 tempFluxXR, tempFluxYR, tempFluxZR,
291 volumeIndex3d(x,y,z),
293 patchData.
dt[patchIndex],
294 *(patchData.
QOut + patchIndex),
303 if (evaluateNonconservativeProduct and Dimensions==3) {
304 for (
int z = 0; z < numberOfVolumesPerAxisInPatch; z++)
305 for (
int y = 0; y < numberOfVolumesPerAxisInPatch; y++)
306 for (
int x = -solversHaloSize; x < numberOfVolumesPerAxisInPatch+solversHaloSize-1; x++) {
307 internal::computeDTerm_LoopBody(
308 QfaceXneg, QfaceXpos, QInterEnumerator, QInEnumerator,
309 nonconservativeProduct,
313 volumeIndex3d(x,y,z),
314 patchData.
t[patchIndex],
315 patchData.
dt[patchIndex],
321 for (
int x = 0; x < numberOfVolumesPerAxisInPatch; x++)
322 for (
int y = -solversHaloSize; y < numberOfVolumesPerAxisInPatch+solversHaloSize-1; y++)
323 for (
int z = 0; z < numberOfVolumesPerAxisInPatch; z++) {
324 internal::computeDTerm_LoopBody(
325 QfaceYneg, QfaceYpos, QInterEnumerator, QInEnumerator,
326 nonconservativeProduct,
330 volumeIndex3d(x,y,z),
331 patchData.
t[patchIndex],
332 patchData.
dt[patchIndex],
338 for (
int x = 0; x < numberOfVolumesPerAxisInPatch; x++)
339 for (
int y = 0; y < numberOfVolumesPerAxisInPatch; y++)
340 for (
int z = -solversHaloSize; z < numberOfVolumesPerAxisInPatch+solversHaloSize-1; z++) {
341 internal::computeDTerm_LoopBody(
342 QfaceZneg, QfaceZpos, QInterEnumerator, QInEnumerator,
343 nonconservativeProduct,
347 volumeIndex3d(x,y,z),
348 patchData.
t[patchIndex],
349 patchData.
dt[patchIndex],
355 for (
int x = 0; x < numberOfVolumesPerAxisInPatch; x++)
356 for (
int y = 0; y < numberOfVolumesPerAxisInPatch; y++)
357 for (
int z = 0; z < numberOfVolumesPerAxisInPatch; z++)
358 for (
int unknown=0; unknown<unknowns; unknown++) {
359 internal::updateSolutionWithDTerm_LoopBody(
360 tempDX, tempDY, tempDZ,
365 volumeIndex3d(x,y,z),
367 patchData.
dt[patchIndex],
368 *(patchData.
QOut + patchIndex),
376 for (
int x = 0; x < numberOfVolumesPerAxisInPatch; x++)
377 for (
int y = 0; y < numberOfVolumesPerAxisInPatch; y++)
378 for (
int z = 0; z < numberOfVolumesPerAxisInPatch; z++)
379 for (
int unknown = 0; unknown < unknowns; unknown++) {
381 internal::updateSolutionwithNCPandSource_LoopBody(
382 patchData.
QIn[patchIndex],
383 QInterEnumerator, QInEnumerator,
384 nonconservativeProduct,
389 volumeIndex3d(x,y,z),
390 patchData.
t[patchIndex],
391 patchData.
dt[patchIndex],
393 patchData.
QOut[patchIndex],
395 evaluateNonconservativeProduct,
401 if (evaluateMaximumEigenvalueAfterTimeStep) {
402 double newMaxEigenvalue = 0.0;
403 for (
int x = 0; x < numberOfVolumesPerAxisInPatch; x++)
404 for (
int y = 0; y < numberOfVolumesPerAxisInPatch; y++)
405 for (
int z = 0; z < numberOfVolumesPerAxisInPatch; z++) {
406 newMaxEigenvalue = std::max(
408 internal::reduceMaxEigenvalue_LoopBody(
409 *(patchData.
QOut + patchIndex),
415 volumeIndex3d(x,y,z),
416 patchData.
t[patchIndex],
417 patchData.
dt[patchIndex]