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