Peano
Loading...
Searching...
No Matches
FD4_batched_static_calls.cpph
Go to the documentation of this file.
1#include "FD4.h"
2#include "tarch/accelerator/accelerator.h"
3#include "tarch/multicore/otter.h"
4#include "exahype2/fd/LoopBody.h"
5#include "exahype2/fd/fd4/LoopBody.h"
6
7
8namespace exahype2 {
9 namespace fd {
10 namespace fd4 {
11 namespace internal{
12 template <
13 typename Solver,
14 int numberOfGridCellsPerPatchPerAxis,
15 int haloSize,
16 int unknowns,
17 int auxiliaryVariables,
18 typename TempDataEnumerator
19 >
20 static void timeStep_batched_static_calls(
21 ::exahype2::CellData<double, double>& patchData,
22 double KOSigma,
23 bool evaluateFlux,
24 bool evaluateDifferentialSource, //for ncp
25 bool evaluateAlgebraicSource, //real source
26 bool copyOldTimeStepAndScaleWithTimeStepSize,
27 DifferentialSourceTermVariant variant,
28 const exahype2::enumerator::AoSLexicographicEnumerator& QInEnumerator,
29 const exahype2::enumerator::AoSLexicographicEnumerator& QOutEnumerator,
30 const TempDataEnumerator& fluxEnumerator,
31 const TempDataEnumerator& differentialSourceEnumerator,
32 const TempDataEnumerator& KODissipationEnumerator,
33 double* tempFluxX,
34 double* tempFluxY,
35 double* tempFluxZ,
36 double* tempDifferentialSourceX,
37 double* tempDifferentialSourceY,
38 double* tempDifferentialSourceZ,
39 double* tempKODissipationX,
40 double* tempKODissipationY,
41 double* tempKODissipationZ
42 ) {
43 const bool evaluateKODissipation = KOSigma>0.0;
44
45 // STEP 1: copy solution from old array to new array, or set all to zero
46 if ( Dimensions==3 and copyOldTimeStepAndScaleWithTimeStepSize ) {
47 for (int x = 0; x < numberOfGridCellsPerPatchPerAxis; x++)
48 for (int y = 0; y < numberOfGridCellsPerPatchPerAxis; y++)
49 for (int z = 0; z < numberOfGridCellsPerPatchPerAxis; z++)
50 for (int unknown=0; unknown<unknowns; unknown++)
51 for (int patchIndex=0; patchIndex<patchData.numberOfCells; patchIndex++) {
52 const double timeStepSize = copyOldTimeStepAndScaleWithTimeStepSize ? patchData.dt[patchIndex] : 1.0;
53 ::exahype2::fd::internal::copySolution_LoopBody(
54 patchData.QIn[patchIndex],
55 QInEnumerator,
56 patchIndex,
57 gridCellIndex3d(x,y,z),
58 unknown,
59 patchData.QOut[patchIndex],
60 QOutEnumerator
61 );
62 }
63 }
64 else if (Dimensions==3 and not copyOldTimeStepAndScaleWithTimeStepSize){
65 for (int x = 0; x < numberOfGridCellsPerPatchPerAxis; x++)
66 for (int y = 0; y < numberOfGridCellsPerPatchPerAxis; y++)
67 for (int z = 0; z < numberOfGridCellsPerPatchPerAxis; z++)
68 for (int patchIndex=0; patchIndex<patchData.numberOfCells; patchIndex++) {
69 const double timeStepSize = copyOldTimeStepAndScaleWithTimeStepSize ? patchData.dt[patchIndex] : 1.0;
70 for (int unknown=0; unknown<unknowns; unknown++) {
71 ::exahype2::fd::internal::clearSolution_LoopBody(
72 patchIndex,
73 gridCellIndex3d(x,y,z),
74 unknown,
75 patchData.QOut[patchIndex],
76 QOutEnumerator
77 );
78 }
79 }
80 }
81 else if (Dimensions==2 and copyOldTimeStepAndScaleWithTimeStepSize) {
82 for (int x = 0; x < numberOfGridCellsPerPatchPerAxis; x++)
83 for (int y = 0; y < numberOfGridCellsPerPatchPerAxis; y++)
84 for (int unknown=0; unknown<unknowns; unknown++)
85 for (int patchIndex=0; patchIndex<patchData.numberOfCells; patchIndex++) {
86 const double timeStepSize = copyOldTimeStepAndScaleWithTimeStepSize ? patchData.dt[patchIndex] : 1.0;
87 ::exahype2::fd::internal::copySolution_LoopBody(
88 patchData.QIn[patchIndex],
89 QInEnumerator,
90 patchIndex,
91 gridCellIndex2d(x,y),
92 unknown,
93 patchData.QOut[patchIndex],
94 QOutEnumerator
95 );
96 }
97 }
98 else {
99 for (int x = 0; x < numberOfGridCellsPerPatchPerAxis; x++)
100 for (int y = 0; y < numberOfGridCellsPerPatchPerAxis; y++)
101 for (int unknown=0; unknown<unknowns; unknown++)
102 for (int patchIndex=0; patchIndex<patchData.numberOfCells; patchIndex++) {
103 const double timeStepSize = copyOldTimeStepAndScaleWithTimeStepSize ? patchData.dt[patchIndex] : 1.0;
104 ::exahype2::fd::internal::clearSolution_LoopBody(
105 patchIndex,
106 gridCellIndex2d(x,y),
107 unknown,
108 patchData.QOut[patchIndex],
109 QOutEnumerator
110 );
111 }
112 }
113
114
115 // Step 2: Add source terms
116 if (Dimensions==2 and evaluateAlgebraicSource) {
117 for (int x = 0; x < numberOfGridCellsPerPatchPerAxis; x++)
118 for (int y = 0; y < numberOfGridCellsPerPatchPerAxis; y++)
119 for (int patchIndex=0; patchIndex<patchData.numberOfCells; patchIndex++) {
120 const double timeStepSize = copyOldTimeStepAndScaleWithTimeStepSize ? patchData.dt[patchIndex] : 1.0;
121 ::exahype2::fd::internal::addAlgebraicSourceTerm_LoopBody<Solver>(
122 patchData.QIn[patchIndex],
123 QInEnumerator,
124 patchData.cellCentre[patchIndex],
125 patchData.cellSize[patchIndex],
126 patchIndex,
127 gridCellIndex2d(x,y),
128 patchData.t[patchIndex],
129 timeStepSize,
130 patchData.QOut[patchIndex],
131 QOutEnumerator
132 );
133 }
134 }
135 else if (Dimensions==3 and evaluateAlgebraicSource) {
136 for (int x = 0; x < numberOfGridCellsPerPatchPerAxis; x++)
137 for (int y = 0; y < numberOfGridCellsPerPatchPerAxis; y++)
138 for (int z = 0; z < numberOfGridCellsPerPatchPerAxis; z++)
139 for (int patchIndex=0; patchIndex<patchData.numberOfCells; patchIndex++) {
140 const double timeStepSize = copyOldTimeStepAndScaleWithTimeStepSize ? patchData.dt[patchIndex] : 1.0;
141 ::exahype2::fd::internal::addAlgebraicSourceTerm_LoopBody<Solver>(
142 patchData.QIn[patchIndex],
143 QInEnumerator,
144 patchData.cellCentre[patchIndex],
145 patchData.cellSize[patchIndex],
146 patchIndex,
147 gridCellIndex3d(x,y,z),
148 patchData.t[patchIndex],
149 timeStepSize,
150 patchData.QOut[patchIndex],
151 QOutEnumerator
152 );
153 }
154 }
155
156 // STEP 3: Evaluate flux and update solution accordingly
157 if (evaluateFlux and Dimensions==2) {
158 for (int x = 0; x < numberOfGridCellsPerPatchPerAxis; x++)
159 for (int y = 0; y < numberOfGridCellsPerPatchPerAxis; y++)
160 for (int patchIndex=0; patchIndex<patchData.numberOfCells; patchIndex++) {
161 const double timeStepSize = copyOldTimeStepAndScaleWithTimeStepSize ? patchData.dt[patchIndex] : 1.0;
162 internal::computeFlux_LoopBody<Solver>(
163 patchData.QIn[patchIndex],
164 QInEnumerator,
165 patchData.cellCentre[patchIndex],
166 patchData.cellSize[patchIndex],
167 patchIndex,
168 gridCellIndex2d(x,y),
169 patchData.t[patchIndex],
170 timeStepSize,
171 0, // normal
172 tempFluxX,
173 QOutEnumerator
174 );
175 }
176 for (int x = 0; x < numberOfGridCellsPerPatchPerAxis; x++)
177 for (int y = 0; y < numberOfGridCellsPerPatchPerAxis; y++)
178 for (int patchIndex=0; patchIndex<patchData.numberOfCells; patchIndex++) {
179 const double timeStepSize = copyOldTimeStepAndScaleWithTimeStepSize ? patchData.dt[patchIndex] : 1.0;
180 internal::computeFlux_LoopBody<Solver>(
181 patchData.QIn[patchIndex],
182 QInEnumerator,
183 patchData.cellCentre[patchIndex],
184 patchData.cellSize[patchIndex],
185 patchIndex,
186 gridCellIndex2d(x,y),
187 patchData.t[patchIndex],
188 timeStepSize,
189 1, // normal
190 tempFluxY,
191 QOutEnumerator
192 );
193 }
194 for (int x = 0; x < numberOfGridCellsPerPatchPerAxis; x++)
195 for (int y = 0; y < numberOfGridCellsPerPatchPerAxis; y++)
196 for (int unknown=0; unknown<unknowns; unknown++)
197 for (int patchIndex=0; patchIndex<patchData.numberOfCells; patchIndex++) {
198 const double timeStepSize = copyOldTimeStepAndScaleWithTimeStepSize ? patchData.dt[patchIndex] : 1.0;
199 // @todo 2d missing
200 internal::updateSolutionWithFlux_LoopBody(
201 tempFluxX, tempFluxY, tempFluxZ,
202 QOutEnumerator,
203 patchData.cellCentre[patchIndex],
204 patchData.cellSize[patchIndex],
205 patchIndex,
206 gridCellIndex2d(x,y),
207 unknown,
208 timeStepSize,
209 *(patchData.QOut + patchIndex),
210 QOutEnumerator
211 );
212 }
213 }
214 else if (evaluateFlux and Dimensions==3) {
215 for (int x = 0; x < numberOfGridCellsPerPatchPerAxis; x++)
216 for (int y = 0; y < numberOfGridCellsPerPatchPerAxis; y++)
217 for (int z = 0; z < numberOfGridCellsPerPatchPerAxis; z++)
218 for (int patchIndex=0; patchIndex<patchData.numberOfCells; patchIndex++) {
219 const double timeStepSize = copyOldTimeStepAndScaleWithTimeStepSize ? patchData.dt[patchIndex] : 1.0;
220 internal::computeFlux_LoopBody<Solver>(
221 patchData.QIn[patchIndex],
222 QInEnumerator,
223 patchData.cellCentre[patchIndex],
224 patchData.cellSize[patchIndex],
225 patchIndex,
226 gridCellIndex3d(x,y,z),
227 patchData.t[patchIndex],
228 timeStepSize,
229 0, // normal
230 tempFluxX,
231 QOutEnumerator
232 );
233 }
234 for (int x = 0; x < numberOfGridCellsPerPatchPerAxis; x++)
235 for (int y = 0; y < numberOfGridCellsPerPatchPerAxis; y++)
236 for (int z = 0; z < numberOfGridCellsPerPatchPerAxis; z++)
237 for (int patchIndex=0; patchIndex<patchData.numberOfCells; patchIndex++) {
238 const double timeStepSize = copyOldTimeStepAndScaleWithTimeStepSize ? patchData.dt[patchIndex] : 1.0;
239 internal::computeFlux_LoopBody<Solver>(
240 patchData.QIn[patchIndex],
241 QInEnumerator,
242 patchData.cellCentre[patchIndex],
243 patchData.cellSize[patchIndex],
244 patchIndex,
245 gridCellIndex3d(x,y,z),
246 patchData.t[patchIndex],
247 timeStepSize,
248 1, // normal
249 tempFluxY,
250 QOutEnumerator
251 );
252 }
253 for (int x = 0; x < numberOfGridCellsPerPatchPerAxis; x++)
254 for (int y = 0; y < numberOfGridCellsPerPatchPerAxis; y++)
255 for (int z = 0; z < numberOfGridCellsPerPatchPerAxis; z++)
256 for (int patchIndex=0; patchIndex<patchData.numberOfCells; patchIndex++) {
257 const double timeStepSize = copyOldTimeStepAndScaleWithTimeStepSize ? patchData.dt[patchIndex] : 1.0;
258 internal::computeFlux_LoopBody<Solver>(
259 patchData.QIn[patchIndex],
260 QInEnumerator,
261 patchData.cellCentre[patchIndex],
262 patchData.cellSize[patchIndex],
263 patchIndex,
264 gridCellIndex3d(x,y,z),
265 patchData.t[patchIndex],
266 timeStepSize,
267 2, // normal
268 tempFluxZ,
269 QOutEnumerator
270 );
271 }
272
273 for (int x = 0; x < numberOfGridCellsPerPatchPerAxis; x++)
274 for (int y = 0; y < numberOfGridCellsPerPatchPerAxis; y++)
275 for (int z = 0; z < numberOfGridCellsPerPatchPerAxis; z++)
276 for (int unknown=0; unknown<unknowns; unknown++)
277 for (int patchIndex=0; patchIndex<patchData.numberOfCells; patchIndex++) {
278 const double timeStepSize = copyOldTimeStepAndScaleWithTimeStepSize ? patchData.dt[patchIndex] : 1.0;
279 internal::updateSolutionWithFlux_LoopBody(
280 tempFluxX, tempFluxY, tempFluxZ,
281 QOutEnumerator,
282 patchData.cellCentre[patchIndex],
283 patchData.cellSize[patchIndex],
284 patchIndex,
285 gridCellIndex3d(x,y,z),
286 unknown,
287 timeStepSize,
288 *(patchData.QOut + patchIndex),
289 QOutEnumerator
290 );
291 }
292 }
293
294 // STEP 4: Evaluate differential source and update solution accordingly
295 if (evaluateDifferentialSource and Dimensions==2) {
296 for (int x = 0; x < numberOfGridCellsPerPatchPerAxis; x++)
297 for (int y = 0; y < numberOfGridCellsPerPatchPerAxis; y++)
298 for (int patchIndex=0; patchIndex<patchData.numberOfCells; patchIndex++) {
299 const double timeStepSize = copyOldTimeStepAndScaleWithTimeStepSize ? patchData.dt[patchIndex] : 1.0;
300 internal::computeDifferentialSourceTerm_LoopBody<Solver>(
301 patchData.QIn[patchIndex],
302 QInEnumerator,
303 patchData.cellCentre[patchIndex],
304 patchData.cellSize[patchIndex],
305 patchIndex,
306 gridCellIndex2d(x,y),
307 patchData.t[patchIndex],
308 timeStepSize,
309 0, // normal
310 tempDifferentialSourceX,
311 QOutEnumerator,
312 variant
313 );
314 }
315 for (int x = 0; x < numberOfGridCellsPerPatchPerAxis; x++)
316 for (int y = 0; y < numberOfGridCellsPerPatchPerAxis; y++)
317 for (int patchIndex=0; patchIndex<patchData.numberOfCells; patchIndex++) {
318 const double timeStepSize = copyOldTimeStepAndScaleWithTimeStepSize ? patchData.dt[patchIndex] : 1.0;
319 internal::computeDifferentialSourceTerm_LoopBody<Solver>(
320 patchData.QIn[patchIndex],
321 QInEnumerator,
322 patchData.cellCentre[patchIndex],
323 patchData.cellSize[patchIndex],
324 patchIndex,
325 gridCellIndex2d(x,y),
326 patchData.t[patchIndex],
327 timeStepSize,
328 1, // normal
329 tempDifferentialSourceY,
330 QOutEnumerator,
331 variant
332 );
333 }
334 for (int x = 0; x < numberOfGridCellsPerPatchPerAxis; x++)
335 for (int y = 0; y < numberOfGridCellsPerPatchPerAxis; y++)
336 for (int unknown=0; unknown<unknowns; unknown++)
337 for (int patchIndex=0; patchIndex<patchData.numberOfCells; patchIndex++) {
338 const double timeStepSize = copyOldTimeStepAndScaleWithTimeStepSize ? patchData.dt[patchIndex] : 1.0;
339 internal::updateSolutionWithDifferentialSourceTerm_LoopBody(
340 tempDifferentialSourceX, tempDifferentialSourceY, tempDifferentialSourceZ,
341 QOutEnumerator,
342 patchData.cellCentre[patchIndex],
343 patchData.cellSize[patchIndex],
344 patchIndex,
345 gridCellIndex2d(x,y),
346 unknown,
347 timeStepSize,
348 *(patchData.QOut + patchIndex),
349 QOutEnumerator
350 );
351 }
352 }
353 else if (evaluateDifferentialSource and Dimensions==3) {
354 for (int x = 0; x < numberOfGridCellsPerPatchPerAxis; x++)
355 for (int y = 0; y < numberOfGridCellsPerPatchPerAxis; y++)
356 for (int z = 0; z < numberOfGridCellsPerPatchPerAxis; z++)
357 for (int patchIndex=0; patchIndex<patchData.numberOfCells; patchIndex++) {
358 const double timeStepSize = copyOldTimeStepAndScaleWithTimeStepSize ? patchData.dt[patchIndex] : 1.0;
359 internal::computeDifferentialSourceTerm_LoopBody<Solver>(
360 patchData.QIn[patchIndex],
361 QInEnumerator,
362 patchData.cellCentre[patchIndex],
363 patchData.cellSize[patchIndex],
364 patchIndex,
365 gridCellIndex3d(x,y,z),
366 patchData.t[patchIndex],
367 timeStepSize,
368 0, // normal
369 tempDifferentialSourceX,
370 QOutEnumerator,
371 variant
372 );
373 }
374 for (int x = 0; x < numberOfGridCellsPerPatchPerAxis; x++)
375 for (int y = 0; y < numberOfGridCellsPerPatchPerAxis; y++)
376 for (int z = 0; z < numberOfGridCellsPerPatchPerAxis; z++)
377 for (int patchIndex=0; patchIndex<patchData.numberOfCells; patchIndex++) {
378 const double timeStepSize = copyOldTimeStepAndScaleWithTimeStepSize ? patchData.dt[patchIndex] : 1.0;
379 internal::computeDifferentialSourceTerm_LoopBody<Solver>(
380 patchData.QIn[patchIndex],
381 QInEnumerator,
382 patchData.cellCentre[patchIndex],
383 patchData.cellSize[patchIndex],
384 patchIndex,
385 gridCellIndex3d(x,y,z),
386 patchData.t[patchIndex],
387 timeStepSize,
388 1, // normal
389 tempDifferentialSourceY,
390 QOutEnumerator,
391 variant
392 );
393 }
394 for (int x = 0; x < numberOfGridCellsPerPatchPerAxis; x++)
395 for (int y = 0; y < numberOfGridCellsPerPatchPerAxis; y++)
396 for (int z = 0; z < numberOfGridCellsPerPatchPerAxis; z++)
397 for (int patchIndex=0; patchIndex<patchData.numberOfCells; patchIndex++) {
398 const double timeStepSize = copyOldTimeStepAndScaleWithTimeStepSize ? patchData.dt[patchIndex] : 1.0;
399 internal::computeDifferentialSourceTerm_LoopBody<Solver>(
400 patchData.QIn[patchIndex],
401 QInEnumerator,
402 patchData.cellCentre[patchIndex],
403 patchData.cellSize[patchIndex],
404 patchIndex,
405 gridCellIndex3d(x,y,z),
406 patchData.t[patchIndex],
407 timeStepSize,
408 2, // normal
409 tempDifferentialSourceZ,
410 QOutEnumerator,
411 variant
412 );
413 }
414 for (int x = 0; x < numberOfGridCellsPerPatchPerAxis; x++)
415 for (int y = 0; y < numberOfGridCellsPerPatchPerAxis; y++)
416 for (int z = 0; z < numberOfGridCellsPerPatchPerAxis; z++)
417 for (int unknown=0; unknown<unknowns; unknown++)
418 for (int patchIndex=0; patchIndex<patchData.numberOfCells; patchIndex++) {
419 const double timeStepSize = copyOldTimeStepAndScaleWithTimeStepSize ? patchData.dt[patchIndex] : 1.0;
420 internal::updateSolutionWithDifferentialSourceTerm_LoopBody(
421 tempDifferentialSourceX, tempDifferentialSourceY, tempDifferentialSourceZ,
422 QOutEnumerator,
423 patchData.cellCentre[patchIndex],
424 patchData.cellSize[patchIndex],
425 patchIndex,
426 gridCellIndex3d(x,y,z),
427 unknown,
428 timeStepSize,
429 *(patchData.QOut + patchIndex),
430 QOutEnumerator
431 );
432 }
433 }
434
435 // STEP 5: compute the Kreiss Oliger dissipation and add it to the solution
436 if (evaluateKODissipation and Dimensions==2) {
437 for (int x = 0; x < numberOfGridCellsPerPatchPerAxis; x++)
438 for (int y = 0; y < numberOfGridCellsPerPatchPerAxis; y++)
439 for (int patchIndex=0; patchIndex<patchData.numberOfCells; patchIndex++) {
440 const double timeStepSize = copyOldTimeStepAndScaleWithTimeStepSize ? patchData.dt[patchIndex] : 1.0;
441 internal::computeKreissOligerDissipationTerm_LoopBody(
442 patchData.QIn[patchIndex],
443 QInEnumerator,
444 patchData.cellCentre[patchIndex],
445 patchData.cellSize[patchIndex],
446 patchIndex,
447 gridCellIndex2d(x,y),
448 patchData.t[patchIndex],
449 timeStepSize,
450 0, // normal
451 tempKODissipationX,
452 QOutEnumerator
453 );
454 }
455 for (int x = 0; x < numberOfGridCellsPerPatchPerAxis; x++)
456 for (int y = 0; y < numberOfGridCellsPerPatchPerAxis; y++)
457 for (int patchIndex=0; patchIndex<patchData.numberOfCells; patchIndex++) {
458 const double timeStepSize = copyOldTimeStepAndScaleWithTimeStepSize ? patchData.dt[patchIndex] : 1.0;
459 internal::computeKreissOligerDissipationTerm_LoopBody(
460 patchData.QIn[patchIndex],
461 QInEnumerator,
462 patchData.cellCentre[patchIndex],
463 patchData.cellSize[patchIndex],
464 patchIndex,
465 gridCellIndex2d(x,y),
466 patchData.t[patchIndex],
467 timeStepSize,
468 1, // normal
469 tempKODissipationY,
470 QOutEnumerator
471 );
472 }
473 for (int x = 0; x < numberOfGridCellsPerPatchPerAxis; x++)
474 for (int y = 0; y < numberOfGridCellsPerPatchPerAxis; y++)
475 for (int unknown=0; unknown<unknowns; unknown++)
476 for (int patchIndex=0; patchIndex<patchData.numberOfCells; patchIndex++) {
477 const double timeStepSize = copyOldTimeStepAndScaleWithTimeStepSize ? patchData.dt[patchIndex] : 1.0;
478 internal::updateSolutionWithKODissipationTerm_LoopBody(
479 KOSigma, tempKODissipationX, tempKODissipationY, tempKODissipationZ,
480 QOutEnumerator,
481 patchData.cellCentre[patchIndex],
482 patchData.cellSize[patchIndex],
483 patchIndex,
484 gridCellIndex2d(x,y),
485 unknown,
486 timeStepSize,
487 *(patchData.QOut + patchIndex),
488 QOutEnumerator
489 );
490 }
491 }
492 else if (evaluateKODissipation and Dimensions==3) {
493 for (int x = 0; x < numberOfGridCellsPerPatchPerAxis; x++)
494 for (int y = 0; y < numberOfGridCellsPerPatchPerAxis; y++)
495 for (int z = 0; z < numberOfGridCellsPerPatchPerAxis; z++)
496 for (int patchIndex=0; patchIndex<patchData.numberOfCells; patchIndex++) {
497 const double timeStepSize = copyOldTimeStepAndScaleWithTimeStepSize ? patchData.dt[patchIndex] : 1.0;
498 internal::computeKreissOligerDissipationTerm_LoopBody(
499 patchData.QIn[patchIndex],
500 QInEnumerator,
501 patchData.cellCentre[patchIndex],
502 patchData.cellSize[patchIndex],
503 patchIndex,
504 gridCellIndex3d(x,y,z),
505 patchData.t[patchIndex],
506 timeStepSize,
507 0, // normal
508 tempKODissipationX,
509 QOutEnumerator
510 );
511 }
512 for (int x = 0; x < numberOfGridCellsPerPatchPerAxis; x++)
513 for (int y = 0; y < numberOfGridCellsPerPatchPerAxis; y++)
514 for (int z = 0; z < numberOfGridCellsPerPatchPerAxis; z++)
515 for (int patchIndex=0; patchIndex<patchData.numberOfCells; patchIndex++) {
516 const double timeStepSize = copyOldTimeStepAndScaleWithTimeStepSize ? patchData.dt[patchIndex] : 1.0;
517 internal::computeKreissOligerDissipationTerm_LoopBody(
518 patchData.QIn[patchIndex],
519 QInEnumerator,
520 patchData.cellCentre[patchIndex],
521 patchData.cellSize[patchIndex],
522 patchIndex,
523 gridCellIndex3d(x,y,z),
524 patchData.t[patchIndex],
525 timeStepSize,
526 1, // normal
527 tempKODissipationY,
528 QOutEnumerator
529 );
530 }
531 for (int x = 0; x < numberOfGridCellsPerPatchPerAxis; x++)
532 for (int y = 0; y < numberOfGridCellsPerPatchPerAxis; y++)
533 for (int z = 0; z < numberOfGridCellsPerPatchPerAxis; z++)
534 for (int patchIndex=0; patchIndex<patchData.numberOfCells; patchIndex++) {
535 const double timeStepSize = copyOldTimeStepAndScaleWithTimeStepSize ? patchData.dt[patchIndex] : 1.0;
536 internal::computeKreissOligerDissipationTerm_LoopBody(
537 patchData.QIn[patchIndex],
538 QInEnumerator,
539 patchData.cellCentre[patchIndex],
540 patchData.cellSize[patchIndex],
541 patchIndex,
542 gridCellIndex3d(x,y,z),
543 patchData.t[patchIndex],
544 timeStepSize,
545 2, // normal
546 tempKODissipationZ,
547 QOutEnumerator
548 );
549 }
550 for (int x = 0; x < numberOfGridCellsPerPatchPerAxis; x++)
551 for (int y = 0; y < numberOfGridCellsPerPatchPerAxis; y++)
552 for (int z = 0; z < numberOfGridCellsPerPatchPerAxis; z++)
553 for (int unknown=0; unknown<unknowns; unknown++)
554 for (int patchIndex=0; patchIndex<patchData.numberOfCells; patchIndex++) {
555 const double timeStepSize = copyOldTimeStepAndScaleWithTimeStepSize ? patchData.dt[patchIndex] : 1.0;
556 internal::updateSolutionWithKODissipationTerm_LoopBody(
557 KOSigma, tempKODissipationX, tempKODissipationY, tempKODissipationZ,
558 QOutEnumerator,
559 patchData.cellCentre[patchIndex],
560 patchData.cellSize[patchIndex],
561 patchIndex,
562 gridCellIndex3d(x,y,z),
563 unknown,
564 timeStepSize,
565 *(patchData.QOut + patchIndex),
566 QOutEnumerator
567 );
568 }
569 }
570 }
571 }
572 }
573 }
574}
575
576
577template <
578 typename Solver,
579 int numberOfGridCellsPerPatchPerAxis,
580 int haloSize,
581 int unknowns,
582 int auxiliaryVariables,
583 typename TempDataEnumerator
584>
585void exahype2::fd::fd4::timeStep_batched_heap_static_calls(
586 ::exahype2::CellData<double, double>& patchData,
587 double KOSigma,
588 bool evaluateFlux,
589 bool evaluateDifferentialSource, //for ncp
590 bool evaluateAlgebraicSource, //real source
591 bool copyOldTimeStepAndScaleWithTimeStepSize,
592 DifferentialSourceTermVariant variant
593) {
594 assertionEquals(haloSize,3);
595
596 const bool evaluateKODissipation = KOSigma>0.0;
597
598 exahype2::enumerator::AoSLexicographicEnumerator QInEnumerator (1,numberOfGridCellsPerPatchPerAxis,haloSize,unknowns,auxiliaryVariables);
599 exahype2::enumerator::AoSLexicographicEnumerator QOutEnumerator(1,numberOfGridCellsPerPatchPerAxis,0, unknowns,0);
600
601 TempDataEnumerator fluxEnumerator (patchData.numberOfCells, numberOfGridCellsPerPatchPerAxis, 0, unknowns, 0);
602 TempDataEnumerator differentialSourceEnumerator(patchData.numberOfCells, numberOfGridCellsPerPatchPerAxis, 0, unknowns, 0);
603 TempDataEnumerator KODissipationEnumerator (patchData.numberOfCells, numberOfGridCellsPerPatchPerAxis, 0, unknowns, 0);
604
605 double* tempFluxX = evaluateFlux ? tarch::allocateMemory( fluxEnumerator.size(),tarch::MemoryLocation::Heap ) : nullptr;
606 double* tempFluxY = evaluateFlux ? tarch::allocateMemory( fluxEnumerator.size(),tarch::MemoryLocation::Heap ) : nullptr;
607 double* tempFluxZ = (evaluateFlux and Dimensions==3) ? tarch::allocateMemory( fluxEnumerator.size(),tarch::MemoryLocation::Heap ) : nullptr;
608 double* tempDifferentialSourceX = evaluateDifferentialSource ? tarch::allocateMemory( differentialSourceEnumerator.size(),tarch::MemoryLocation::Heap ) : nullptr;
609 double* tempDifferentialSourceY = evaluateDifferentialSource ? tarch::allocateMemory( differentialSourceEnumerator.size(),tarch::MemoryLocation::Heap ) : nullptr;
610 double* tempDifferentialSourceZ = (evaluateDifferentialSource and Dimensions==3) ? tarch::allocateMemory( differentialSourceEnumerator.size(),tarch::MemoryLocation::Heap ) : nullptr;
611 double* tempKODissipationX = evaluateKODissipation ? tarch::allocateMemory( KODissipationEnumerator.size(),tarch::MemoryLocation::Heap ) : nullptr;
612 double* tempKODissipationY = evaluateKODissipation ? tarch::allocateMemory( KODissipationEnumerator.size(),tarch::MemoryLocation::Heap ) : nullptr;
613 double* tempKODissipationZ = (evaluateKODissipation and Dimensions==3) ? tarch::allocateMemory( KODissipationEnumerator.size(),tarch::MemoryLocation::Heap ) : nullptr;
614
615 internal::timeStep_batched_static_calls<Solver, numberOfGridCellsPerPatchPerAxis, haloSize, unknowns, auxiliaryVariables, TempDataEnumerator>(
616 patchData,
617 KOSigma,
618 evaluateFlux,
619 evaluateDifferentialSource,
620 evaluateAlgebraicSource,
621 copyOldTimeStepAndScaleWithTimeStepSize,
622 variant,
623 QInEnumerator,
624 QOutEnumerator,
625 fluxEnumerator,
626 differentialSourceEnumerator,
627 KODissipationEnumerator,
628 tempFluxX,
629 tempFluxY,
630 tempFluxZ,
631 tempDifferentialSourceX,
632 tempDifferentialSourceY,
633 tempDifferentialSourceZ,
634 tempKODissipationX,
635 tempKODissipationY,
636 tempKODissipationZ
637 );
638
639 if (tempFluxX!=nullptr) tarch::freeMemory(tempFluxX, tarch::MemoryLocation::Heap);
640 if (tempFluxY!=nullptr) tarch::freeMemory(tempFluxY, tarch::MemoryLocation::Heap);
641 if (tempFluxZ!=nullptr) tarch::freeMemory(tempFluxZ, tarch::MemoryLocation::Heap);
642 if (tempDifferentialSourceX!=nullptr) tarch::freeMemory(tempDifferentialSourceX, tarch::MemoryLocation::Heap);
643 if (tempDifferentialSourceY!=nullptr) tarch::freeMemory(tempDifferentialSourceY, tarch::MemoryLocation::Heap);
644 if (tempDifferentialSourceZ!=nullptr) tarch::freeMemory(tempDifferentialSourceZ, tarch::MemoryLocation::Heap);
645 if (tempKODissipationX!=nullptr) tarch::freeMemory(tempKODissipationX, tarch::MemoryLocation::Heap);
646 if (tempKODissipationY!=nullptr) tarch::freeMemory(tempKODissipationY, tarch::MemoryLocation::Heap);
647 if (tempKODissipationZ!=nullptr) tarch::freeMemory(tempKODissipationZ, tarch::MemoryLocation::Heap);
648}
649
650
651
652template <
653 typename Solver,
654 int numberOfGridCellsPerPatchPerAxis,
655 int haloSize,
656 int unknowns,
657 int auxiliaryVariables,
658 typename TempDataEnumerator
659>
660static void exahype2::fd::fd4::timeStep_batched_heap_multicore_static_calls(
661 ::exahype2::CellData<double, double>& patchData,
662 double KOSigma,
663 bool evaluateFlux,
664 bool evaluateDifferentialSource, //for ncp
665 bool evaluateAlgebraicSource, //real source
666 bool copyOldTimeStepAndScaleWithTimeStepSize,
667 DifferentialSourceTermVariant variant
668) {
669 #if defined(SharedOMP)
670 ::exahype2::fd::fd4::omp::timeStep_batched_heap_multicore_static_calls<
671 Solver,
672 numberOfGridCellsPerPatchPerAxis,
673 haloSize,
674 unknowns,
675 auxiliaryVariables,
676 TempDataEnumerator
677 >(
678 patchData,
679 KOSigma,
680 evaluateFlux,
681 evaluateDifferentialSource, //for ncp
682 evaluateAlgebraicSource, //real source
683 copyOldTimeStepAndScaleWithTimeStepSize,
684 variant
685 );
686 #else
687 assertionMsg(false, "not implemented yet");
688 #endif
689}
assertionEquals
#define assertionEquals(lhs, rhs)
Definition Assertions.h:1034
assertionMsg
#define assertionMsg(expr, message)
Definition Assertions.h:1076
exahype2::fd::fd4::internal::updateSolutionWithFlux_LoopBody
static void updateSolutionWithFlux_LoopBody(const double *__restrict__ tempFluxX, const double *__restrict__ tempFluxY, const double *__restrict__ tempFluxZ, const QOutEnumeratorType &fluxEnumerator, const tarch::la::Vector< Dimensions, double > &patchCentre, const tarch::la::Vector< Dimensions, double > &patchSize, int patchIndex, const tarch::la::Vector< Dimensions, int > &volumeIndex, int unknown, double dt, double *__restrict__ QOut, const QOutEnumeratorType &QOutEnumerator) InlineMethod
Plain update of flux in a finite differences scheme.
exahype2::fd::fd4::internal::updateSolutionWithDifferentialSourceTerm_LoopBody
static void updateSolutionWithDifferentialSourceTerm_LoopBody(const double *__restrict__ QDiffSrcX, const double *__restrict__ QDiffSrcY, const double *__restrict__ QDiffSrcZ, const QOutEnumeratorType &QDiffSrcEnumerator, const tarch::la::Vector< Dimensions, double > &patchCentre, const tarch::la::Vector< Dimensions, double > &patchSize, int patchIndex, const tarch::la::Vector< Dimensions, int > &volumeIndex, int unknown, double dt, double *__restrict__ QOut, const QOutEnumeratorType &QOutEnumerator) InlineMethod
exahype2::fd::fd4::internal::updateSolutionWithKODissipationTerm_LoopBody
static void updateSolutionWithKODissipationTerm_LoopBody(const double KOSigma, const double *__restrict__ QKODspX, const double *__restrict__ QKODspY, const double *__restrict__ QKODspZ, const QOutEnumeratorType &QKODspEnumerator, const tarch::la::Vector< Dimensions, double > &patchCentre, const tarch::la::Vector< Dimensions, double > &patchSize, int patchIndex, const tarch::la::Vector< Dimensions, int > &volumeIndex, int unknown, double dt, double *__restrict__ QOut, const QOutEnumeratorType &QOutEnumerator) InlineMethod
exahype2::fd::fd4::internal::computeKreissOligerDissipationTerm_LoopBody
static void computeKreissOligerDissipationTerm_LoopBody(const double *__restrict__ QIn, const QInEnumeratorType &QInEnumerator, const tarch::la::Vector< Dimensions, double > &patchCentre, const tarch::la::Vector< Dimensions, double > &patchSize, int patchIndex, const tarch::la::Vector< Dimensions, int > &volumeIndex, double t, double dt, int normal, double *__restrict__ QKODsp, const QOutEnumeratorType &QKODspEnumerator) InlineMethod
exahype2::fd::fd4::internal::timeStep_batched_static_calls
static void timeStep_batched_static_calls(::exahype2::CellData< double, double > &patchData, double KOSigma, bool evaluateFlux, bool evaluateDifferentialSource, bool evaluateAlgebraicSource, bool copyOldTimeStepAndScaleWithTimeStepSize, DifferentialSourceTermVariant variant, const exahype2::enumerator::AoSLexicographicEnumerator &QInEnumerator, const exahype2::enumerator::AoSLexicographicEnumerator &QOutEnumerator, const TempDataEnumerator &fluxEnumerator, const TempDataEnumerator &differentialSourceEnumerator, const TempDataEnumerator &KODissipationEnumerator, double *tempFluxX, double *tempFluxY, double *tempFluxZ, double *tempDifferentialSourceX, double *tempDifferentialSourceY, double *tempDifferentialSourceZ, double *tempKODissipationX, double *tempKODissipationY, double *tempKODissipationZ)
Definition FD4_batched_static_calls.cpph:20
exahype2::fd::fd4::omp::timeStep_batched_heap_multicore_static_calls
static void timeStep_batched_heap_multicore_static_calls(::exahype2::CellData< double, double > &patchData, double KOSigma, bool evaluateFlux, bool evaluateDifferentialSource, bool evaluateAlgebraicSource, bool copyOldTimeStepAndScaleWithTimeStepSize, DifferentialSourceTermVariant variant)
exahype2::fd::fd4::timeStep_batched_heap_static_calls
static void timeStep_batched_heap_static_calls(::exahype2::CellData< double, double > &patchData, double KOSigma, bool evaluateFlux, bool evaluateDifferentialSource, bool evaluateAlgebraicSource, bool copyOldTimeStepAndScaleWithTimeStepSize, DifferentialSourceTermVariant variant)
exahype2::fd::fd4::timeStep_batched_heap_multicore_static_calls
static void timeStep_batched_heap_multicore_static_calls(::exahype2::CellData< double, double > &patchData, double KOSigma, bool evaluateFlux, bool evaluateDifferentialSource, bool evaluateAlgebraicSource, bool copyOldTimeStepAndScaleWithTimeStepSize, DifferentialSourceTermVariant variant)
exahype2::fd::internal::copySolution_LoopBody
static void copySolution_LoopBody(const double *__restrict__ QIn, const QInEnumeratorType &QInEnumerator, int patchIndex, const tarch::la::Vector< Dimensions, int > &volumeIndex, int unknown, double *__restrict__ QOut, const QOutEnumeratorType &QOutEnumerator) InlineMethod
exahype2::fd::internal::clearSolution_LoopBody
static void clearSolution_LoopBody(int patchIndex, const tarch::la::Vector< Dimensions, int > &volumeIndex, int unknown, double *__restrict__ QOut, const QOutEnumeratorType &QOutEnumerator) InlineMethod
The complicated way to write =0.
exahype2
For the generic kernels that I use here most of the time.
Definition CellAccess.h:13
tarch::freeMemory
void freeMemory(void *data, MemoryLocation location, int device=accelerator::Device::HostDevice)
Free memory.
Definition accelerator.cpp:30
tarch::MemoryLocation::Heap
@ Heap
Create data on the heap of the local device.
tarch::allocateMemory
T * allocateMemory(std::size_t count, MemoryLocation location, int device=accelerator::Device::HostDevice)
Allocates memory on the specified memory location.
Definition accelerator.h:99
exahype2::CellData
Representation of a number of cells which contains all information that's required to process the sto...
Definition CellData.h:77
exahype2::CellData::QOut
outType ** QOut
Out values.
Definition CellData.h:116
exahype2::CellData::QIn
inType ** QIn
QIn may not be const, as some kernels delete it straightaway once the input data has been handled.
Definition CellData.h:82
exahype2::CellData::numberOfCells
const int numberOfCells
As we store data as SoA, we have to know how big the actual arrays are.
Definition CellData.h:99
exahype2::CellData::cellCentre
tarch::la::Vector< Dimensions, double > * cellCentre
Definition CellData.h:83
exahype2::CellData::cellSize
tarch::la::Vector< Dimensions, double > * cellSize
Definition CellData.h:84