Peano
Loading...
Searching...
No Matches
MusclHancock_patchwise_functors.cpph
Go to the documentation of this file.
2#include "MusclHancock.h"
3
4
7 int numberOfVolumesPerAxisInPatch,
8 int haloSize,
9 int unknowns,
10 int auxiliaryVariables,
11 bool evaluateFlux,
12 bool evaluateNonconservativeProduct,
13 bool evaluateSource,
14 bool evaluateMaximumEigenvalueAfterTimeStep,
15 Flux flux,
16 NonconservativeProduct nonconservativeProduct,
17 Source source,
18 MaxEigenvalue maxEigenvalue
19) {
20 static tarch::logging::Log _log( "exahype2::fv::musclhancock" );
21
26 constexpr int solversHaloSize = 1;
27
28 assertionEquals(haloSize,1);
29 assertionEquals(3*unknowns,auxiliaryVariables);
30
31 exahype2::enumerator::AoSLexicographicEnumerator QInEnumerator ( 1, numberOfVolumesPerAxisInPatch, haloSize,unknowns,auxiliaryVariables);
32 exahype2::enumerator::AoSLexicographicEnumerator QOutEnumerator ( 1, numberOfVolumesPerAxisInPatch, 0,unknowns,auxiliaryVariables);
33 exahype2::enumerator::AoSLexicographicEnumerator fluxEnumerator (patchData.numberOfCells, numberOfVolumesPerAxisInPatch, solversHaloSize,unknowns, 0);
34 exahype2::enumerator::AoSLexicographicEnumerator ncpEnumerator (patchData.numberOfCells, numberOfVolumesPerAxisInPatch+1, solversHaloSize,unknowns, 0);
35 exahype2::enumerator::AoSLexicographicEnumerator eigenvalueEnumerator(patchData.numberOfCells, numberOfVolumesPerAxisInPatch, solversHaloSize, 1, 0);
36
37 double* tempFluxXL = evaluateFlux ? tarch::allocateMemory( fluxEnumerator.size(), tarch::MemoryLocation::Heap ) : nullptr;
38 double* tempFluxYL = evaluateFlux ? tarch::allocateMemory( fluxEnumerator.size(), tarch::MemoryLocation::Heap ) : nullptr;
39 double* tempFluxZL = (evaluateFlux and Dimensions==3) ? tarch::allocateMemory( fluxEnumerator.size(), tarch::MemoryLocation::Heap ) : nullptr;
40 double* tempFluxXR = evaluateFlux ? tarch::allocateMemory( fluxEnumerator.size(), tarch::MemoryLocation::Heap ) : nullptr;
41 double* tempFluxYR = evaluateFlux ? tarch::allocateMemory( fluxEnumerator.size(), tarch::MemoryLocation::Heap ) : nullptr;
42 double* tempFluxZR = (evaluateFlux and Dimensions==3) ? tarch::allocateMemory( fluxEnumerator.size(), tarch::MemoryLocation::Heap ) : nullptr;
43 double* tempDX = evaluateNonconservativeProduct ? tarch::allocateMemory( ncpEnumerator.size(), tarch::MemoryLocation::Heap ) : nullptr;
44 double* tempDY = evaluateNonconservativeProduct ? tarch::allocateMemory( ncpEnumerator.size(), tarch::MemoryLocation::Heap ) : nullptr;
45 double* tempDZ = (evaluateNonconservativeProduct and Dimensions==3) ? tarch::allocateMemory( ncpEnumerator.size(), tarch::MemoryLocation::Heap ) : nullptr;
46 double* tempEigenvalueX = tarch::allocateMemory( eigenvalueEnumerator.size(), tarch::MemoryLocation::Heap );
47 double* tempEigenvalueY = tarch::allocateMemory( eigenvalueEnumerator.size(), tarch::MemoryLocation::Heap );
48 double* tempEigenvalueZ = (Dimensions==3) ? tarch::allocateMemory( eigenvalueEnumerator.size(), tarch::MemoryLocation::Heap ) : nullptr;
49
50 exahype2::enumerator::AoSLexicographicEnumerator QInterEnumerator ( 1, numberOfVolumesPerAxisInPatch, haloSize,unknowns, 0);
51 //used to sotre the time derivatives
52 double* timeDerivative = tarch::allocateMemory( QInterEnumerator.size(), tarch::MemoryLocation::Heap );
53 //face value
54 double* QfaceXneg = tarch::allocateMemory( QInterEnumerator.size(), tarch::MemoryLocation::Heap );
55 double* QfaceXpos = tarch::allocateMemory( QInterEnumerator.size(), tarch::MemoryLocation::Heap );
56 double* QfaceYneg = tarch::allocateMemory( QInterEnumerator.size(), tarch::MemoryLocation::Heap );
57 double* QfaceYpos = tarch::allocateMemory( QInterEnumerator.size(), tarch::MemoryLocation::Heap );
58 double* QfaceZneg = tarch::allocateMemory( QInterEnumerator.size(), tarch::MemoryLocation::Heap );
59 double* QfaceZpos = tarch::allocateMemory( QInterEnumerator.size(), tarch::MemoryLocation::Heap );
60
61
62 for (int patchIndex=0; patchIndex<patchData.numberOfCells; patchIndex++) {
63 //logDebug( "timeStepWithMusclHancock_sequential(...)", "evaluate source terms and/or copy solution over" );
64 //step1:copy solution
65 if (Dimensions==3) {
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],
72 QInEnumerator,
73 patchIndex,
74 volumeIndex3d(x,y,z),
75 unknown,
76 patchData.QOut[patchIndex],
77 QOutEnumerator
78 );
79 }
80 }
81
82 //step2: calculate the time derivatives
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++)
86 {
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)) )
90 {}
91 //for now, we assume that we use a full setup (include flux, ncp and source)
92 else {
93 internal::computeTimeDerivative_LoopBody(
94 patchData.QIn[patchIndex],
95 QInEnumerator,
96 flux, nonconservativeProduct, source,
97 patchData.cellCentre[patchIndex],
98 patchData.cellSize[patchIndex],
99 patchIndex,
100 volumeIndex3d(x,y,z),
101 patchData.t[patchIndex],
102 patchData.dt[patchIndex],
103 timeDerivative,
104 QInterEnumerator
105 );
106 }
107 }
108
109 //step2.1: calculate the quantites on the face
110 //Q_i, dt, dx_i -> Q^n+1/2_i-1/2, so on...
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)) )
117 {}
118 else {
119 internal::computeQonFace_LoopBody(
120 patchData.QIn[patchIndex],
121 QInEnumerator,
122 patchData.cellCentre[patchIndex],
123 patchData.cellSize[patchIndex],
124 patchIndex,
125 volumeIndex3d(x,y,z),
126 patchData.t[patchIndex],
127 patchData.dt[patchIndex],
128 timeDerivative,
129 QfaceXneg, QfaceXpos,
130 QfaceYneg, QfaceYpos,
131 QfaceZneg, QfaceZpos,
132 QInterEnumerator
133 );
134 }
135
136 //step3: calculate max eigenvalue and update the solution with eigenvalue damping
137 if (Dimensions==3) {
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],
143 QInEnumerator,
144 maxEigenvalue,
145 patchData.cellCentre[patchIndex],
146 patchData.cellSize[patchIndex],
147 patchIndex,
148 volumeIndex3d(x,y,z),
149 patchData.t[patchIndex], patchData.dt[patchIndex],
150 0,
151 tempEigenvalueX,
152 eigenvalueEnumerator
153 );
154 }
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],
160 QInEnumerator,
161 maxEigenvalue,
162 patchData.cellCentre[patchIndex],
163 patchData.cellSize[patchIndex],
164 patchIndex,
165 volumeIndex3d(x,y,z),
166 patchData.t[patchIndex], patchData.dt[patchIndex],
167 1,
168 tempEigenvalueY,
169 eigenvalueEnumerator
170 );
171 }
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],
177 QInEnumerator,
178 maxEigenvalue,
179 patchData.cellCentre[patchIndex],
180 patchData.cellSize[patchIndex],
181 patchIndex,
182 volumeIndex3d(x,y,z),
183 patchData.t[patchIndex], patchData.dt[patchIndex],
184 2,
185 tempEigenvalueZ,
186 eigenvalueEnumerator
187 );
188 }
189
190
192//todo everyting below ask for q on the face rather than q at volume center!
194
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],
201 QInEnumerator,
202 tempEigenvalueX,
203 tempEigenvalueY,
204 tempEigenvalueZ,
205 eigenvalueEnumerator,
206 patchData.cellCentre[patchIndex],
207 patchData.cellSize[patchIndex],
208 patchIndex,
209 volumeIndex3d(x,y,z),
210 unknown,
211 patchData.dt[patchIndex],
212 patchData.QOut[patchIndex],
213 QOutEnumerator,
214 QfaceXneg, QfaceXpos,
215 QfaceYneg, QfaceYpos,
216 QfaceZneg, QfaceZpos,
217 QInterEnumerator
218 );
219 }
220 }
221
222
223//now we have both space and time derivative
224
225 //step4: compute and update flux part
226 //std::cout<<"flux check: "<<std::to_string(1)<<" fluxRight: "<<std::to_string(2)<<std::endl;
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,
233 flux,
234 patchData.cellCentre[patchIndex],
235 patchData.cellSize[patchIndex],
236 patchIndex,
237 volumeIndex3d(x,y,z),
238 patchData.t[patchIndex],
239 patchData.dt[patchIndex],
240 0, // normal
241 tempFluxXL, tempFluxXR,
242 fluxEnumerator
243 );
244 }
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,
250 flux,
251 patchData.cellCentre[patchIndex],
252 patchData.cellSize[patchIndex],
253 patchIndex,
254 volumeIndex3d(x,y,z),
255 patchData.t[patchIndex],
256 patchData.dt[patchIndex],
257 1, // normal
258 tempFluxYL, tempFluxYR,
259 fluxEnumerator
260 );
261 }
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,
267 flux,
268 patchData.cellCentre[patchIndex],
269 patchData.cellSize[patchIndex],
270 patchIndex,
271 volumeIndex3d(x,y,z),
272 patchData.t[patchIndex],
273 patchData.dt[patchIndex],
274 2, // normal
275 tempFluxZL, tempFluxZR,
276 fluxEnumerator
277 );
278 }
279
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,
287 fluxEnumerator,
288 patchData.cellCentre[patchIndex],
289 patchData.cellSize[patchIndex],
290 patchIndex,
291 volumeIndex3d(x,y,z),
292 unknown,
293 patchData.dt[patchIndex],
294 *(patchData.QOut + patchIndex),
295 QOutEnumerator
296 );
297 }
298 }
299
300
301 //step5: compute and update D part
302 //tricky thing: it is actually identical what we use for ncp now
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,
310 patchData.cellCentre[patchIndex],
311 patchData.cellSize[patchIndex],
312 patchIndex,
313 volumeIndex3d(x,y,z),
314 patchData.t[patchIndex],
315 patchData.dt[patchIndex],
316 0, // normal
317 tempDX,
318 ncpEnumerator
319 );
320 }
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,
327 patchData.cellCentre[patchIndex],
328 patchData.cellSize[patchIndex],
329 patchIndex,
330 volumeIndex3d(x,y,z),
331 patchData.t[patchIndex],
332 patchData.dt[patchIndex],
333 1, // normal
334 tempDY,
335 ncpEnumerator
336 );
337 }
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,
344 patchData.cellCentre[patchIndex],
345 patchData.cellSize[patchIndex],
346 patchIndex,
347 volumeIndex3d(x,y,z),
348 patchData.t[patchIndex],
349 patchData.dt[patchIndex],
350 2, // normal
351 tempDZ,
352 ncpEnumerator
353 );
354 }
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,
361 ncpEnumerator,
362 patchData.cellCentre[patchIndex],
363 patchData.cellSize[patchIndex],
364 patchIndex,
365 volumeIndex3d(x,y,z),
366 unknown,
367 patchData.dt[patchIndex],
368 *(patchData.QOut + patchIndex),
369 QOutEnumerator
370 );
371 }
372 }
373
374 //step6: update the source term and ncp term
375 //this function now cause crash in foccz4
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++) {
380 //for now, we assume that we use a full setup (include ncp and source)
381 internal::updateSolutionwithNCPandSource_LoopBody(
382 patchData.QIn[patchIndex],
383 QInterEnumerator, QInEnumerator,
384 nonconservativeProduct,
385 source,
386 patchData.cellCentre[patchIndex],
387 patchData.cellSize[patchIndex],
388 patchIndex,
389 volumeIndex3d(x,y,z),
390 patchData.t[patchIndex],
391 patchData.dt[patchIndex],
392 timeDerivative,
393 patchData.QOut[patchIndex],
394 QOutEnumerator,
395 evaluateNonconservativeProduct,
396 evaluateSource
397 );
398 }
399
400 //step7: compute maximum eigenvalue for subsequent time step
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(
407 newMaxEigenvalue,
408 internal::reduceMaxEigenvalue_LoopBody(
409 *(patchData.QOut + patchIndex),
410 QOutEnumerator,
411 maxEigenvalue,
412 patchData.cellCentre[patchIndex],
413 patchData.cellSize[patchIndex],
414 patchIndex,
415 volumeIndex3d(x,y,z),
416 patchData.t[patchIndex],
417 patchData.dt[patchIndex]
418 )
419 );
420 }
421 patchData.maxEigenvalue[patchIndex] = newMaxEigenvalue;
422 }
423 }
424
425 if (tempFluxXL!=nullptr) tarch::freeMemory(tempFluxXL, tarch::MemoryLocation::Heap);
426 if (tempFluxYL!=nullptr) tarch::freeMemory(tempFluxYL, tarch::MemoryLocation::Heap);
427 if (tempFluxZL!=nullptr) tarch::freeMemory(tempFluxZL, tarch::MemoryLocation::Heap);
428 if (tempFluxXR!=nullptr) tarch::freeMemory(tempFluxXR, tarch::MemoryLocation::Heap);
429 if (tempFluxYR!=nullptr) tarch::freeMemory(tempFluxYR, tarch::MemoryLocation::Heap);
430 if (tempFluxZR!=nullptr) tarch::freeMemory(tempFluxZR, tarch::MemoryLocation::Heap);
431 if (tempDX!=nullptr) tarch::freeMemory(tempDX, tarch::MemoryLocation::Heap);
432 if (tempDY!=nullptr) tarch::freeMemory(tempDY, tarch::MemoryLocation::Heap);
433 if (tempDZ!=nullptr) tarch::freeMemory(tempDZ, tarch::MemoryLocation::Heap);
434 if (tempEigenvalueX!=nullptr) tarch::freeMemory(tempEigenvalueX, tarch::MemoryLocation::Heap);
435 if (tempEigenvalueY!=nullptr) tarch::freeMemory(tempEigenvalueY, tarch::MemoryLocation::Heap);
436 if (tempEigenvalueZ!=nullptr) tarch::freeMemory(tempEigenvalueZ, tarch::MemoryLocation::Heap);
444}
445
#define assertionEquals(lhs, rhs)
Log Device.
Definition Log.h:516
tarch::logging::Log _log("exahype2::fv")
static void timeStepWithMusclHancock_patchwise_heap_functors(::exahype2::CellData< double, double > &patchData, int numberOfVolumesPerAxisInPatch, int haloSize, int unknowns, int auxiliaryVariables, bool evaluateFlux, bool evaluateNonconservativeProduct, bool evaluateSource, bool evaluateMaximumEigenvalueAfterTimeStep, Flux flux, NonconservativeProduct nonconservativeProduct, Source source, MaxEigenvalue maxEigenvalue)
void freeMemory(void *data, MemoryLocation location, int device=accelerator::Device::HostDevice)
Free memory.
@ Heap
Create data on the heap of the local device.
T * allocateMemory(std::size_t count, MemoryLocation location, int device=accelerator::Device::HostDevice)
Allocates memory on the specified memory location.
Definition accelerator.h:99
Representation of a number of cells which contains all information that's required to process the sto...
Definition CellData.h:77
outType ** QOut
Out values.
Definition CellData.h:116
inType ** QIn
QIn may not be const, as some kernels delete it straightaway once the input data has been handled.
Definition CellData.h:82
const int numberOfCells
As we store data as SoA, we have to know how big the actual arrays are.
Definition CellData.h:99
double * maxEigenvalue
Out values.
Definition CellData.h:121
tarch::la::Vector< Dimensions, double > * cellCentre
Definition CellData.h:83
tarch::la::Vector< Dimensions, double > * cellSize
Definition CellData.h:84
GPUCallableInlineMethod int size() const InlineMethod