Peano 4
Loading...
Searching...
No Matches
KernelUtils.cpp
Go to the documentation of this file.
1#include "KernelUtils.h"
2
3#if defined(GPUOffloadingOMP)
4#pragma omp declare target
5#endif
6#if defined(GPUOffloadingCPP)
7int exahype2::aderdg::getNodesPerCell(int nodesPerAxis) {
8#else
10#endif
11 int nodesPerCell = 1;
12 for ( int d = 0; d < Dimensions; d++ ) { nodesPerCell *= nodesPerAxis; }
13 return nodesPerCell;
14}
15#if defined(GPUOffloadingOMP)
16#pragma omp end declare target
17#endif
18
19#if defined(GPUOffloadingOMP)
20#pragma omp declare target
21#endif
22#if defined(GPUOffloadingCPP)
24#else
26#endif
27 int nodesPerCell = 1;
28 for ( int d = 0; d < Dimensions+1; d++ ) { nodesPerCell *= nodesPerAxis; }
29 return nodesPerCell;
30}
31#if defined(GPUOffloadingOMP)
32#pragma omp end declare target
33#endif
34
35#if defined(GPUOffloadingOMP)
36#pragma omp declare target
37#endif
38#if defined(GPUOffloadingCPP)
40#else
42#endif
43 int nodesPerAxis,
44 bool strides4SpaceTimeQuantity) {
46 for ( int d = 0; d < Dimensions; d++ ) {
47 strides[d+1] *= strides[d]*nodesPerAxis;
48 }
49 if ( !strides4SpaceTimeQuantity ) {
50 for ( int d = 0; d < Dimensions+1; d++ ) {
51 strides[d] /= nodesPerAxis;
52 }
53 }
54 return strides;
55}
56#if defined(GPUOffloadingOMP)
57#pragma omp end declare target
58#endif
59
60#if defined(GPUOffloadingOMP)
61#pragma omp declare target
62#endif
63#if defined(GPUOffloadingCPP)
65#else
67#endif
70) {
71 int scalarIndex = 0;
72 for ( int d = 0; d < Dimensions+1; d++ ) { // t -> x -> y -> z
73 scalarIndex += strides[d]*index[d];
74 }
75 return scalarIndex;
76}
77#if defined(GPUOffloadingOMP)
78#pragma omp end declare target
79#endif
80
81#if defined(GPUOffloadingOMP)
82#pragma omp declare target
83#endif
84#if defined(GPUOffloadingCPP)
86#else
88#endif
89 int scalarIndex,
91 tarch::la::Vector<Dimensions+1,int> index(-1); // valid index component values are non-negative
92 int tmp = 0;
93 if ( strides[0] > 0 ) {
94 for ( int d = Dimensions; d >= 0; d-- ) { // z -> y -> x -> t
95 index[d] = (scalarIndex - tmp) / strides[d];
96 tmp += strides[d]*index[d]; // iz -> iy -> ix -> it
97 }
98 } else {
99 for ( int d = Dimensions; d >= 1; d-- ) { // z -> y -> x
100 index[d] = (scalarIndex - tmp) / strides[d];
101 tmp += strides[d]*index[d]; // iz -> iy -> ix
102 }
103 }
104 return index;
105}
106#if defined(GPUOffloadingOMP)
107#pragma omp end declare target
108#endif
109
110#if defined(GPUOffloadingOMP)
111#pragma omp declare target
112#endif
113#if defined(GPUOffloadingCPP)
115#else
117#endif
120 const double dx,
121 const double t,
122 const double dt,
123 const double* __restrict__ nodes) {
125
126 coords[0]= t;
127 if ( index[0] >= 0 ) {
128 coords[0] += nodes[index[0]] * dt;
129 }
130 for ( int d = 1; d < Dimensions+1; d++ ) { // x -> y -> z
131 coords[d] = centre[d-1] + dx/*[d-1]*/ * (nodes[index[d]] - 0.5);
132 }
133 return coords;
134}
135#if defined(GPUOffloadingOMP)
136#pragma omp end declare target
137#endif
138
139#if defined(GPUOffloadingOMP)
140#pragma omp declare target
141#endif
142#if defined(GPUOffloadingCPP)
144#else
146#endif
147 const tarch::la::Vector<Dimensions+1,int>& indexOnFace,
148 const tarch::la::Vector<Dimensions,double>& faceCentre,
149 const int direction,
150 const double dx,
151 const double t,
152 const double dt,
153 const double* __restrict__ nodes) {
155
156 coords[0]= t;
157 if ( indexOnFace[0] >= 0 ) {
158 coords[0] += nodes[indexOnFace[0]] * dt;
159 }
160 int e = 1;
161 for ( int d = 1; d < Dimensions+1; d++ ) { // x -> y -> z
162 if ( d-1 == direction ) {
163 coords[d] = faceCentre[d-1];
164 } else {
165 coords[d] = faceCentre[d-1] + dx/*[d-1]*/ * (nodes[indexOnFace[e]] - 0.5);
166 }
167 e++;
168 }
169 return coords;
170}
171#if defined(GPUOffloadingOMP)
172#pragma omp end declare target
173#endif
174
175#if defined(GPUOffloadingOMP)
176#pragma omp declare target
177#endif
178#if defined(GPUOffloadingCPP)
180#else
182#endif
184 const int direction,
185 const int nodesPerAxis
186) {
187 // freeze spatial dimension direction (indexCell[direction+1])
188 int scalarIndexFace = 0;
189 int stride = 1;
190 for ( int e=0; e < Dimensions; e++ ) { // ordering (fastest running left): (y,z), (x,z), (x,y)
191 if ( e != direction ) {
192 scalarIndexFace += stride*indexCell[e+1];
193 stride *= nodesPerAxis;
194 }
195 }
196 return scalarIndexFace;
197}
198#if defined(GPUOffloadingOMP)
199#pragma omp end declare target
200#endif
201
202#if defined(GPUOffloadingOMP)
203#pragma omp declare target
204#endif
205#if defined(GPUOffloadingCPP)
207#else
209#endif
211 const int direction,
212 const int orientation,
213 const int nodesPerAxis
214) {
215 // freeze spatial dimension direction (indexCell[direction+1])
216 int scalarIndexHull = 0;
217 int stride = 1;
218 for ( int e=0; e < Dimensions; e++ ) { // ordering (fastest running left): (y,z), (x,z), (x,y)
219 if ( e != direction ) {
220 scalarIndexHull += stride*indexCell[e+1];
221 stride *= nodesPerAxis;
222 }
223 }
224 //const int faceOffset = ( direction*2 + orientation ) * stride; // stride = nodesPerAxis^{d-1}
225 const int faceOffset = ( direction + Dimensions*orientation ) * stride; // stride = nodesPerAxis^{d-1}
226 scalarIndexHull += faceOffset;
227 return scalarIndexHull;
228}
229#if defined(GPUOffloadingOMP)
230#pragma omp end declare target
231#endif
232
233#if defined(GPUOffloadingOMP)
234#pragma omp declare target
235#endif
236#if defined(GPUOffloadingCPP)
238#else
240#endif
242 const int direction,
243 const int id,
244 const int nodesPerAxis
245) {
247 indexCell[0] = indexFace[0]; // time
248 indexCell[direction+1] = id; // node in interior
249 int i = 1;
250 for ( int e = 0; e < Dimensions; e++) { // take from face
251 if ( e != direction ) {
252 indexCell[e+1] = indexFace[i++];
253 }
254 }
255 return lineariseIndex(indexCell,getStrides(nodesPerAxis));
256}
257#if defined(GPUOffloadingOMP)
258#pragma omp end declare target
259#endif
260
261#if defined(GPUOffloadingOMP)
262#pragma omp declare target
263#endif
264#if defined(GPUOffloadingCPP)
266#else
268#endif
269 const double* __restrict__ const QIn,
270 const double* __restrict__ const dudx,
271 const double invDx,
272 const int nodesPerAxis,
273 const int strideQ,
274 const int scalarIndex,
275 double* __restrict__ gradQ) {
276 // Always works with space time input data
277 const tarch::la::Vector<Dimensions+1,int> strides = getStrides(nodesPerAxis);
278 const tarch::la::Vector<Dimensions+1,int> index = delineariseIndex(scalarIndex, strides);
279
280 for ( int d = 0; d < Dimensions; d++ ) { // x -> y -> z
281 // zero
282 for (int var=0; var < strideQ; var++) {
283 gradQ[ d*strideQ + var ] = 0.0;
284 //gradQ[ d + var*Dimensions ] = 0.0;
285 }
286 // sum up
287 for (int a=0; a < nodesPerAxis; a++) {
288 const double coeff = invDx/*[d]*/ * dudx[ a*nodesPerAxis + index[d+1] ];
289 for (int var=0; var < strideQ; var++) {
290 gradQ[ d*strideQ + var ] += coeff * QIn[ ( scalarIndex + (a-index[d+1])*strides[d+1] )*strideQ + var ];
291 //gradQ[ d + var*Dimensions ] += coeff * QIn[ ( scalarIndex + (a-index[d+1])*strides[d+1] )*strideQ + var ];
292 }
293 }
294 }
295}
296#if defined(GPUOffloadingOMP)
297#pragma omp end declare target
298#endif
299
300#if defined(GPUOffloadingOMP)
301#pragma omp declare target
302#endif
303#if defined(GPUOffloadingCPP)
305#else
313#if defined(GPUOffloadingOMP)
314#pragma omp end declare target
315#endif
316
323#if defined(GPUOffloadingOMP)
324#pragma omp declare target
325#endif
326#if defined(GPUOffloadingCPP)
328#else
330#endif
331 const double* __restrict__ const UIn,
332 const double* __restrict__ const nodes,
333 const double* __restrict__ const barycentricWeights,
334 const tarch::la::Vector<Dimensions,double>& referenceCoodinates,
335 const int nodesPerAxis,
336 const int strideQ,
337 double* __restrict__ pointwiseQOut) {
338 double l = 1.0;
339 // clear
340 for ( int var = 0; var < strideQ; var++ ) {
341 pointwiseQOut[ var ] = 0.0;
342 }
343 // compute
344 for ( int d = 0; d < Dimensions; d++ ) {
345 for ( int i = 0; i < nodesPerAxis; i++ ) {
346 l *= (referenceCoodinates[d] - nodes[i]);
347 }
348 }
349 const int nodesPerCell = getNodesPerCell(nodesPerAxis);
350 for ( int scalarIndex = 0; scalarIndex < nodesPerCell; scalarIndex++ ) {
351 const tarch::la::Vector<Dimensions+1,int> strides = getStrides(nodesPerAxis,false);
352 const tarch::la::Vector<Dimensions+1,int> index = delineariseIndex(scalarIndex, strides);
353
354 double coeff = l;
355 for ( int d = 0; d < Dimensions; d++ ){
356 coeff *= barycentricWeights[index[1+d]] / (referenceCoodinates[d] - nodes[index[1+d]] );
357 }
358 for ( int var = 0; var < strideQ; var++ ) {
359 pointwiseQOut[ var ] += coeff * UIn[ scalarIndex*strideQ + var ];
360 }
361 }
362}
363#if defined(GPUOffloadingOMP)
364#pragma omp end declare target
365#endif
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
#define GPUCallableMethod
Definition accelerator.h:25
GPUCallableMethod tarch::la::Vector< Dimensions+1, int > delineariseIndex(int scalarIndex, const tarch::la::Vector< Dimensions+1, int > &strides)
GPUCallableMethod int mapSpaceTimeFaceIndexToScalarCellIndex(const tarch::la::Vector< Dimensions+1, int > &indexFace, const int direction, const int id, const int nodesPerAxis)
GPUCallableMethod void gradient_AoS(const double *__restrict__ const QIn, const double *__restrict__ const dudx, const double invDx, const int nodesPerAxis, const int strideQ, const int scalarIndex, double *__restrict__ gradQAux)
Computes the gradient at a given (space-time) quadrature point.
GPUCallableMethod int lineariseIndex(const tarch::la::Vector< Dimensions+1, int > &index, const tarch::la::Vector< Dimensions+1, int > &strides)
GPUCallableMethod tarch::la::Vector< Dimensions+1, double > getCoordinatesOnFace(const tarch::la::Vector< Dimensions+1, int > &indexOnFace, const tarch::la::Vector< Dimensions, double > &faceCentre, const int direction, const double dx, const double t, const double dt, const double *__restrict__ nodes)
Compute coordinates from cell geometry and quadrature index.
GPUCallableMethod void interpolate_AoS(const double *__restrict__ const UIn, const double *__restrict__ const nodes, const double *__restrict__ const barycentricWeights, const tarch::la::Vector< Dimensions, double > &referenceCoodinates, const int nodesPerAxis, const int strideQ, double *__restrict__ pointwiseQOut)
Evalutes DG polynomial at arbitrary reference space coordinate in reference domain [0,...
GPUCallableMethod tarch::la::Vector< Dimensions+1, int > getStrides(int nodesPerAxis, bool strides4SpaceTimeQuantity=true)
If the stride is used for a space-time quantity, then this function generates:
GPUCallableMethod int mapCellIndexToScalarHullIndex(const tarch::la::Vector< Dimensions+1, int > &indexCell, const int direction, const int orientation, const int nodesPerAxis)
Map cell index to an index for the whole hull, i.e.
GPUCallableMethod tarch::la::Vector< Dimensions, double > mapToReferenceCoordinates(const tarch::la::Vector< Dimensions, double > &x, const tarch::la::Vector< Dimensions, double > &cellCentre, const double dx)
Map point in physical domain to corresponding coordinates in reference domain.
GPUCallableMethod int getNodesPerCell(int nodesPerAxis)
GPUCallableMethod int mapCellIndexToScalarFaceIndex(const tarch::la::Vector< Dimensions+1, int > &indexCell, const int direction, const int nodesPerAxis)
Map cell index to an index for a single face.
GPUCallableMethod tarch::la::Vector< Dimensions+1, double > getCoordinates(const tarch::la::Vector< Dimensions+1, int > &index, const tarch::la::Vector< Dimensions, double > &centre, const double dx, const double t, const double dt, const double *__restrict__ nodes)
Compute coordinates from cell geometry and quadrature index.
GPUCallableMethod int getSpaceTimeNodesPerCell(int nodesPerAxis)
Simple vector class.
Definition Vector.h:134