Peano 4
Loading...
Searching...
No Matches
SecondOrderAuxiliaryVariablesReconstruction.cpp
Go to the documentation of this file.
3
4
5namespace {
9 template <void (*recompute_LoopBody)(
10 double* __restrict__ QIn,
11 const ::exahype2::enumerator::AoSLexicographicEnumerator& QInEnumerator,
12 const tarch::la::Vector<Dimensions,double>& patchCentre,
14 int patchIndex,
16 int normal
17 )>
18 void recomputeAuxiliaryVariablesFD4(
19 ::exahype2::CellData& patchData,
20 int numberOfGridCellsPerPatchPerAxis,
21 int haloSize,
22 int unknowns,
23 int auxiliaryVariables
24 ) {
25 assertion( haloSize >= 3 );
26 ::exahype2::enumerator::AoSLexicographicEnumerator QInEnumerator (1,numberOfGridCellsPerPatchPerAxis,haloSize,unknowns,auxiliaryVariables);
27
28 for (int patchIndex=0; patchIndex<patchData.numberOfCells; patchIndex++) {
29 if (Dimensions==2) {
30 for (int x = 0; x < numberOfGridCellsPerPatchPerAxis; x++)
31 for (int y = 0; y < numberOfGridCellsPerPatchPerAxis; y++) {
32 recompute_LoopBody(
33 patchData.QIn[patchIndex],
34 QInEnumerator,
35 patchData.cellCentre[patchIndex],
36 patchData.cellSize[patchIndex],
37 patchIndex,
38 gridCellIndex2d(x,y),
39 0 // normal
40 );
41 }
42
43 for (int x = 0; x < numberOfGridCellsPerPatchPerAxis; x++)
44 for (int y = 0; y < numberOfGridCellsPerPatchPerAxis; y++) {
45 recompute_LoopBody(
46 patchData.QIn[patchIndex],
47 QInEnumerator,
48 patchData.cellCentre[patchIndex],
49 patchData.cellSize[patchIndex],
50 patchIndex,
51 gridCellIndex2d(x,y),
52 1 // normal
53 );
54 }
55 }
56 else {
57 for (int x = 0; x < numberOfGridCellsPerPatchPerAxis; x++)
58 for (int y = 0; y < numberOfGridCellsPerPatchPerAxis; y++)
59 for (int z = 0; z < numberOfGridCellsPerPatchPerAxis; z++) {
60 recompute_LoopBody(
61 patchData.QIn[patchIndex],
62 QInEnumerator,
63 patchData.cellCentre[patchIndex],
64 patchData.cellSize[patchIndex],
65 patchIndex,
66 gridCellIndex3d(x,y,z),
67 0 // normal
68 );
69 }
70
71 for (int x = 0; x < numberOfGridCellsPerPatchPerAxis; x++)
72 for (int y = 0; y < numberOfGridCellsPerPatchPerAxis; y++)
73 for (int z = 0; z < numberOfGridCellsPerPatchPerAxis; z++) {
74 recompute_LoopBody(
75 patchData.QIn[patchIndex],
76 QInEnumerator,
77 patchData.cellCentre[patchIndex],
78 patchData.cellSize[patchIndex],
79 patchIndex,
80 gridCellIndex3d(x,y,z),
81 1 // normal
82 );
83 }
84
85 for (int x = 0; x < numberOfGridCellsPerPatchPerAxis; x++)
86 for (int y = 0; y < numberOfGridCellsPerPatchPerAxis; y++)
87 for (int z = 0; z < numberOfGridCellsPerPatchPerAxis; z++) {
88 recompute_LoopBody(
89 patchData.QIn[patchIndex],
90 QInEnumerator,
91 patchData.cellCentre[patchIndex],
92 patchData.cellSize[patchIndex],
93 patchIndex,
94 gridCellIndex3d(x,y,z),
95 2 // normal
96 );
97 }
98 }
99 }
100 }
101}
102
103
105 ::exahype2::CellData& patchData,
106 int numberOfGridCellsPerPatchPerAxis,
107 int haloSize,
108 int unknowns,
109 int auxiliaryVariables
110) {
111 recomputeAuxiliaryVariablesFD4<internal::recomputeAuxiliaryVariablesFD4_centralDifferences_LoopBody>(
112 patchData,
113 numberOfGridCellsPerPatchPerAxis,
114 haloSize,
115 unknowns,
116 auxiliaryVariables
117 );
118}
119
120
122 ::exahype2::CellData& patchData,
123 int numberOfGridCellsPerPatchPerAxis,
124 int haloSize,
125 int unknowns,
126 int auxiliaryVariables
127) {
128 recomputeAuxiliaryVariablesFD4<internal::recomputeAuxiliaryVariablesFD4_leftDifferences_LoopBody>(
129 patchData,
130 numberOfGridCellsPerPatchPerAxis,
131 haloSize,
132 unknowns,
133 auxiliaryVariables
134 );
135}
136
137
139 ::exahype2::CellData& patchData,
140 int numberOfGridCellsPerPatchPerAxis,
141 int haloSize,
142 int unknowns,
143 int auxiliaryVariables
144) {
145 recomputeAuxiliaryVariablesFD4<internal::recomputeAuxiliaryVariablesFD4_rightDifferences_LoopBody>(
146 patchData,
147 numberOfGridCellsPerPatchPerAxis,
148 haloSize,
149 unknowns,
150 auxiliaryVariables
151 );
152}
153
154
156 ::exahype2::CellData& patchData,
157 int numberOfGridCellsPerPatchPerAxis,
158 int haloSize,
159 int unknowns,
160 int auxiliaryVariables
161) {
162 assertion( haloSize >= 3 );
163 ::exahype2::enumerator::AoSLexicographicEnumerator QInEnumerator (1,numberOfGridCellsPerPatchPerAxis,haloSize,unknowns,auxiliaryVariables);
164
165 for (int patchIndex=0; patchIndex<patchData.numberOfCells; patchIndex++) {
166 if (Dimensions==2) {
167 for (int x = -1; x < numberOfGridCellsPerPatchPerAxis+1; x++)
168 for (int y = 0; y < numberOfGridCellsPerPatchPerAxis; y++) {
169 internal::recomputeAuxiliaryVariablesFD4_4thOrder_LoopBody(
170 patchData.QIn[patchIndex],
171 QInEnumerator,
172 patchData.cellCentre[patchIndex],
173 patchData.cellSize[patchIndex],
174 patchIndex,
175 gridCellIndex2d(x,y),
176 0 // normal
177 );
178 }
179
180 for (int x = 0; x < numberOfGridCellsPerPatchPerAxis; x++)
181 for (int y = -1; y < numberOfGridCellsPerPatchPerAxis+1; y++) {
182 internal::recomputeAuxiliaryVariablesFD4_4thOrder_LoopBody(
183 patchData.QIn[patchIndex],
184 QInEnumerator,
185 patchData.cellCentre[patchIndex],
186 patchData.cellSize[patchIndex],
187 patchIndex,
188 gridCellIndex2d(x,y),
189 1 // normal
190 );
191 }
192 }
193 else {
194 for (int x = -1; x < numberOfGridCellsPerPatchPerAxis+1; x++)
195 for (int y = 0; y < numberOfGridCellsPerPatchPerAxis; y++)
196 for (int z = 0; z < numberOfGridCellsPerPatchPerAxis; z++) {
197 internal::recomputeAuxiliaryVariablesFD4_4thOrder_LoopBody(
198 patchData.QIn[patchIndex],
199 QInEnumerator,
200 patchData.cellCentre[patchIndex],
201 patchData.cellSize[patchIndex],
202 patchIndex,
203 gridCellIndex3d(x,y,z),
204 0 // normal
205 );
206 }
207
208 for (int x = 0; x < numberOfGridCellsPerPatchPerAxis; x++)
209 for (int y = -1; y < numberOfGridCellsPerPatchPerAxis+1; y++)
210 for (int z = 0; z < numberOfGridCellsPerPatchPerAxis; z++) {
211 internal::recomputeAuxiliaryVariablesFD4_4thOrder_LoopBody(
212 patchData.QIn[patchIndex],
213 QInEnumerator,
214 patchData.cellCentre[patchIndex],
215 patchData.cellSize[patchIndex],
216 patchIndex,
217 gridCellIndex3d(x,y,z),
218 1 // normal
219 );
220 }
221
222 for (int x = 0; x < numberOfGridCellsPerPatchPerAxis; x++)
223 for (int y = 0; y < numberOfGridCellsPerPatchPerAxis; y++)
224 for (int z = -1; z < numberOfGridCellsPerPatchPerAxis+1; z++) {
225 internal::recomputeAuxiliaryVariablesFD4_4thOrder_LoopBody(
226 patchData.QIn[patchIndex],
227 QInEnumerator,
228 patchData.cellCentre[patchIndex],
229 patchData.cellSize[patchIndex],
230 patchIndex,
231 gridCellIndex3d(x,y,z),
232 2 // normal
233 );
234 }
235 }
236 }
237}
238
239
241 double* __restrict__ QIn,
242 const ::exahype2::enumerator::AoSLexicographicEnumerator& QInEnumerator,
243 const tarch::la::Vector<Dimensions,double>& patchCentre,
245 int patchIndex,
246 const tarch::la::Vector<Dimensions,int>& volumeIndex,
247 int normal
248) {
249 tarch::la::Vector<Dimensions,int> centralVolume = volumeIndex;
250 tarch::la::Vector<Dimensions,int> leftAdjacentVolume = volumeIndex;
251 tarch::la::Vector<Dimensions,int> rightAdjacentVolume = volumeIndex;
252
253 rightAdjacentVolume(normal)++;
254 leftAdjacentVolume(normal)--;
255
256 double cellSize = patchSize(normal) / QInEnumerator._numberOfDoFsPerAxisInCell;
257
258 auto finiteDifferences = [&](int source, int target) {
259 double dx = QIn[ QInEnumerator(patchIndex, rightAdjacentVolume,source) ]
260 - QIn[ QInEnumerator(patchIndex, leftAdjacentVolume, source) ];
261 QIn[ QInEnumerator(patchIndex,centralVolume,target) ] = dx / 2.0 / cellSize;
262 };
263
264 finiteDifferences(16,23+normal);
265 finiteDifferences(54,55+normal);
266
267 for (int i=0; i<3; i++) {
268 finiteDifferences(17+i,26+3*normal+i);
269 }
270
271 for (int i=0; i<3; i++)
272 for (int j=i; j<3; j++) {
273 int k = (i==0?j:i+j+1);
274 finiteDifferences(k,35+6*normal+k);
275 }
276}
277
278
279
280
282 double* __restrict__ QIn,
283 const ::exahype2::enumerator::AoSLexicographicEnumerator& QInEnumerator,
284 const tarch::la::Vector<Dimensions,double>& patchCentre,
286 int patchIndex,
287 const tarch::la::Vector<Dimensions,int>& volumeIndex,
288 int normal
289) {
290 tarch::la::Vector<Dimensions,int> centralVolume = volumeIndex;
291 tarch::la::Vector<Dimensions,int> leftAdjacentVolume = volumeIndex;
292
293 leftAdjacentVolume(normal)--;
294
295 double cellSize = patchSize(normal) / QInEnumerator._numberOfDoFsPerAxisInCell;
296
297 auto finiteDifferences = [&](int source, int target) {
298 double dx = QIn[ QInEnumerator(patchIndex, centralVolume,source) ]
299 - QIn[ QInEnumerator(patchIndex, leftAdjacentVolume, source) ];
300 QIn[ QInEnumerator(patchIndex,centralVolume,target) ] = dx / cellSize;
301 };
302
303 finiteDifferences(16,23+normal);
304 finiteDifferences(54,55+normal);
305
306 for (int i=0; i<3; i++) {
307 finiteDifferences(17+i,26+3*normal+i);
308 }
309
310 for (int i=0; i<3; i++)
311 for (int j=i; j<3; j++) {
312 int k = (i==0?j:i+j+1);
313 finiteDifferences(k,35+6*normal+k);
314 }
315}
316
317
318
319
321 double* __restrict__ QIn,
322 const ::exahype2::enumerator::AoSLexicographicEnumerator& QInEnumerator,
323 const tarch::la::Vector<Dimensions,double>& patchCentre,
325 int patchIndex,
326 const tarch::la::Vector<Dimensions,int>& volumeIndex,
327 int normal
328) {
329 tarch::la::Vector<Dimensions,int> centralVolume = volumeIndex;
330 tarch::la::Vector<Dimensions,int> rightAdjacentVolume = volumeIndex;
331
332 rightAdjacentVolume(normal)++;
333
334 double cellSize = patchSize(normal) / QInEnumerator._numberOfDoFsPerAxisInCell;
335
336 auto finiteDifferences = [&](int source, int target) {
337 double dx = QIn[ QInEnumerator(patchIndex, rightAdjacentVolume,source) ]
338 - QIn[ QInEnumerator(patchIndex, centralVolume, source) ];
339 QIn[ QInEnumerator(patchIndex,centralVolume,target) ] = dx / cellSize;
340 };
341
342 finiteDifferences(16,23+normal);
343 finiteDifferences(54,55+normal);
344
345 for (int i=0; i<3; i++) {
346 finiteDifferences(17+i,26+3*normal+i);
347 }
348
349 for (int i=0; i<3; i++)
350 for (int j=i; j<3; j++) {
351 int k = (i==0?j:i+j+1);
352 finiteDifferences(k,35+6*normal+k);
353 }
354}
355
356
358 double* __restrict__ QIn,
359 const ::exahype2::enumerator::AoSLexicographicEnumerator& QInEnumerator,
360 const tarch::la::Vector<Dimensions,double>& patchCentre,
362 int patchIndex,
363 const tarch::la::Vector<Dimensions,int>& volumeIndex,
364 int normal
365) {
366 tarch::la::Vector<Dimensions,int> centralVolume = volumeIndex;
367 tarch::la::Vector<Dimensions,int> leftAdjacentVolume1 = volumeIndex;
368 tarch::la::Vector<Dimensions,int> leftAdjacentVolume2 = volumeIndex;
369 tarch::la::Vector<Dimensions,int> rightAdjacentVolume1 = volumeIndex;
370 tarch::la::Vector<Dimensions,int> rightAdjacentVolume2 = volumeIndex;
371
372 rightAdjacentVolume1(normal)++;
373 rightAdjacentVolume2(normal)++;
374 rightAdjacentVolume2(normal)++;
375 leftAdjacentVolume1(normal)--;
376 leftAdjacentVolume2(normal)--;
377 leftAdjacentVolume2(normal)--;
378
379 double cellsize = patchSize(normal) / QInEnumerator._numberOfDoFsPerAxisInCell;
380 double gradcoef = 1.0/(12.0*cellsize);
381
382 QIn[ QInEnumerator(patchIndex,centralVolume,23+normal) ] = gradcoef * QIn[ QInEnumerator(patchIndex, leftAdjacentVolume2,16) ] -
383 8.0*gradcoef * QIn[ QInEnumerator(patchIndex, leftAdjacentVolume1,16) ] +
384 8.0*gradcoef * QIn[ QInEnumerator(patchIndex,rightAdjacentVolume1,16) ] -
385 gradcoef * QIn[ QInEnumerator(patchIndex,rightAdjacentVolume2,16) ];
386
387 // I'm afraid I misunderstood the enumeration
388 QIn[ QInEnumerator(patchIndex,centralVolume,55+normal) ] = gradcoef * QIn[ QInEnumerator(patchIndex, leftAdjacentVolume2,54) ] -
389 8.0*gradcoef * QIn[ QInEnumerator(patchIndex, leftAdjacentVolume1,54) ] +
390 8.0*gradcoef * QIn[ QInEnumerator(patchIndex,rightAdjacentVolume1,54) ] -
391 gradcoef * QIn[ QInEnumerator(patchIndex,rightAdjacentVolume2,54) ];
392
393 for (int i=0; i<3; i++) {
394 QIn[ QInEnumerator(patchIndex,centralVolume,26+3*normal+i) ] = gradcoef * QIn[ QInEnumerator(patchIndex, leftAdjacentVolume2,17+i) ] -
395 8.0*gradcoef * QIn[ QInEnumerator(patchIndex, leftAdjacentVolume1,17+i) ] +
396 8.0*gradcoef * QIn[ QInEnumerator(patchIndex,rightAdjacentVolume1,17+i) ] -
397 gradcoef * QIn[ QInEnumerator(patchIndex,rightAdjacentVolume2,17+i) ];
398 }
399
400 for (int i=0; i<3; i++)
401 for (int j=i; j<3; j++) {
402 int k = (i==0?j:i+j+1);
403 QIn[ QInEnumerator(patchIndex,centralVolume,35+6*normal+k) ] = 0.5*gradcoef * QIn[ QInEnumerator(patchIndex, leftAdjacentVolume2,k) ] -
404 4.0*gradcoef * QIn[ QInEnumerator(patchIndex, leftAdjacentVolume1,k) ] +
405 4.0*gradcoef * QIn[ QInEnumerator(patchIndex,rightAdjacentVolume1,k) ] -
406 0.5*gradcoef * QIn[ QInEnumerator(patchIndex,rightAdjacentVolume2,k) ];
407 }
408}
#define assertion(expr)
And from this we can write down f$ nabla phi_i nabla phi_i dx but since we are constructing matrix let s investigate the f$ j
we integrate over each cell and then take the sum across each of the cells We also consider the terms that enter the f$ k
we integrate over each cell and then take the sum across each of the cells We also consider the terms that enter the f$ so we are integrating with f$ phi_k phi_k dx
void recomputeAuxiliaryVariablesFD4_rightDifferences_LoopBody(double *__restrict__ QIn, const ::exahype2::enumerator::AoSLexicographicEnumerator &QInEnumerator, const tarch::la::Vector< Dimensions, double > &patchCentre, const tarch::la::Vector< Dimensions, double > &patchSize, int patchIndex, const tarch::la::Vector< Dimensions, int > &volumeIndex, int normal)
void recomputeAuxiliaryVariablesFD4_leftDifferences_LoopBody(double *__restrict__ QIn, const ::exahype2::enumerator::AoSLexicographicEnumerator &QInEnumerator, const tarch::la::Vector< Dimensions, double > &patchCentre, const tarch::la::Vector< Dimensions, double > &patchSize, int patchIndex, const tarch::la::Vector< Dimensions, int > &volumeIndex, int normal)
void recomputeAuxiliaryVariablesFD4_4thOrder_LoopBody(double *__restrict__ QIn, const ::exahype2::enumerator::AoSLexicographicEnumerator &QInEnumerator, const tarch::la::Vector< Dimensions, double > &patchCentre, const tarch::la::Vector< Dimensions, double > &patchSize, int patchIndex, const tarch::la::Vector< Dimensions, int > &volumeIndex, int normal)
Recompute auxiliary variables.
void recomputeAuxiliaryVariablesFD4_centralDifferences_LoopBody(double *__restrict__ QIn, const ::exahype2::enumerator::AoSLexicographicEnumerator &QInEnumerator, const tarch::la::Vector< Dimensions, double > &patchCentre, const tarch::la::Vector< Dimensions, double > &patchSize, int patchIndex, const tarch::la::Vector< Dimensions, int > &volumeIndex, int normal)
void recomputeAuxiliaryVariablesFD4_4thOrder(::exahype2::CellData &patchData, int numberOfGridCellsPerPatchPerAxis, int haloSize, int unknowns, int auxiliaryVariables)
Recompute auxiliary variables for FD4 scheme with a 4th order scheme.
KeywordToAvoidDuplicateSymbolsForInlinedFunctions void source(double *S, const double *const Q, const int CCZ4LapseType, const double CCZ4ds, const double CCZ4c, const double CCZ4e, const double CCZ4f, const double CCZ4bs, const double CCZ4sk, const double CCZ4xi, const double CCZ4itau, const double CCZ4eta, const double CCZ4k1, const double CCZ4k2, const double CCZ4k3, const double CCZ4SO) InlineMethod
The source term is one out of two terms that we use in our CCZ4 formulation.
void recomputeAuxiliaryVariablesFD4_centralDifferences(::exahype2::CellData &patchData, int numberOfGridCellsPerPatchPerAxis, int haloSize, int unknowns, int auxiliaryVariables)
void recomputeAuxiliaryVariablesFD4_leftDifferences(::exahype2::CellData &patchData, int numberOfGridCellsPerPatchPerAxis, int haloSize, int unknowns, int auxiliaryVariables)
void recomputeAuxiliaryVariablesFD4_rightDifferences(::exahype2::CellData &patchData, int numberOfGridCellsPerPatchPerAxis, int haloSize, int unknowns, int auxiliaryVariables)
auto volumeIndex(Args... args)
Definition VolumeIndex.h:54
tuple z
Definition makeIC.py:53
Representation of a number of cells which contains all information that's required to process the sto...
Definition CellData.h:79
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
Simple vector class.
Definition Vector.h:134