Peano 4
Loading...
Searching...
No Matches
RiemannAoS.cpp
Go to the documentation of this file.
1#include "RiemannAoS.h"
2
3#include "Generic.h"
4
5#include "KernelUtils.h"
6
7/*
8#if defined(GPUOffloadingOMP)
9#pragma omp declare target
10#endif
11*/
13 std::function< void(
14 double * __restrict__ QInside,
15 double * __restrict__ OOutside,
17 double t,
18 int normal
19 ) > boundaryState,
20 double * __restrict__ QOut,
21 double * __restrict__ QIn,
22 const double * const __restrict__ nodes,
24 const double dx,
25 const double t,
26 const double dt,
27 const int nodesPerAxis,
28 const int unknowns,
29 const int strideQ,
30 const int direction,
31 const int scalarIndexFace) {
32 const tarch::la::Vector<Dimensions+1,int> indexQFace = // (t,y,z,face=0) , (t,x,z,face=0), (t,x,y,face=0)
33 delineariseIndex(scalarIndexFace,getStrides(nodesPerAxis));
35 getCoordinatesOnFace(indexQFace,faceCentre,direction,dx,t,dt,nodes);
36 const tarch::la::Vector<Dimensions,double> x(&coords[1]);
37 const double time = coords[0];
38
39 boundaryState( &QIn[ scalarIndexFace*strideQ ], &QOut[ scalarIndexFace*strideQ ], x, time, direction );
40}
41/*
42#if defined(GPUOffloadingOMP)
43#pragma omp end declare target
44#endif
45
46#if defined(GPUOffloadingOMP)
47#pragma omp declare target
48#endif
49*/
51 std::function< double(
52 const double * const __restrict__ Q,
54 double t,
55 const int normal
56 ) > maxAbsoluteEigenvalue,
57 const double * const __restrict__ QLR,
58 const double * const __restrict__ nodes,
60 const double dx,
61 const double t,
62 const double dt,
63 const int nodesPerAxis,
64 const int unknowns,
65 const int strideQ,
66 const int direction,
67 const int scalarIndexFace) {
68 const tarch::la::Vector<Dimensions+1,int> strides = getStrides(nodesPerAxis);
69 const tarch::la::Vector<Dimensions+1,int> indexQHull = // (t,y,z,block) , (t,x,z,block), (t,x,y,block)
70 delineariseIndex(scalarIndexFace,strides);
72 getCoordinatesOnFace(indexQHull,faceCentre,direction,dx,t,dt,nodes); // ignores indexHull[Dimensions]
73 const tarch::la::Vector<Dimensions,double> x(&coords[1]);
74 const double time = coords[0];
75
76 return maxAbsoluteEigenvalue( &QLR[ scalarIndexFace*strideQ ], x, time, direction );
77 // @todo: check this value is > 0.
78}
79/*
80#if defined(GPUOffloadingOMP)
81#pragma omp end declare target
82#endif
83*/
84
85// launchers
86
88 std::function< void(
89 const double * const __restrict__ QInside,
90 double * __restrict__ OOutside,
92 double t,
93 int normal
94 ) > boundaryState,
95 double * QHullOut[Dimensions*2],
96 const double * const __restrict__ nodes,
98 const double dx,
99 const double t,
100 const double dt,
101 const int order,
102 const int unknowns,
103 const int auxiliaryVariables,
104 const bool atBoundary[Dimensions*2]) {
105 const int nodesPerAxis = order + 1;
106
107 const int strideQ = unknowns + auxiliaryVariables;
108
109 const int spaceTimeNodesOnFace = getNodesPerCell(nodesPerAxis); // nodesPerAxis^d
110
111 const int strideQLR = spaceTimeNodesOnFace;
112
113 for (int orientation=0; orientation < 2; orientation++) {
114 for (int direction = 0; direction < Dimensions; direction++) {
115 const int face = Dimensions*orientation + direction;
116 tarch::la::Vector<Dimensions,double> faceCentre = cellCentre;
117 faceCentre[ direction ] += (-1+2*orientation)*0.5*dx;
118
119 // set Dirichlet boundary conditions
120 if ( atBoundary[ face ] ) {
121 double* QL = QHullOut[ face ];
122 double* QR = QHullOut[ face ] + strideQLR*strideQ;
123 double* QOut = ( orientation == 0 ) ? QL : QR;
124 double* QIn = ( orientation == 0 ) ? QR : QL;
125 for ( int scalarIndexFace = 0; scalarIndexFace < spaceTimeNodesOnFace; scalarIndexFace++ ) {
127 boundaryState,
128 QOut,
129 QIn,
130 nodes,
131 faceCentre,
132 dx,
133 t,
134 dt,
135 nodesPerAxis,
136 unknowns,
137 strideQ,
138 direction,
139 scalarIndexFace);
140 }
141 }
142 }
143 }
144}
145
147 std::function< double(
148 const double * const __restrict__ Q,
150 double t,
151 const int direction
152 ) > maxAbsoluteEigenvalue,
153 double maxEigenvaluePerFaceOut[Dimensions*2],
154 double * const QHullIn[Dimensions*2],
155 const double * const __restrict__ nodes,
156 const tarch::la::Vector<Dimensions,double>& cellCentre,
157 const double dx,
158 const double t,
159 const double dt,
160 const int order,
161 const int unknowns,
162 const int auxiliaryVariables) {
163 const int nodesPerAxis = order+1;
164
165 const int strideQ = unknowns+auxiliaryVariables;
166
167 // Each face stores first: the degrees of freedom of QL ("left") and then of QR ("right")
168 const int strideQLR = getNodesPerCell(nodesPerAxis);
169
170 // 6144 work items for order 7 in 3D
171 const int collocatedSpaceTimeNodesOnCellHull = Dimensions*2 * 2 * strideQLR; // #faces * 2 * nodesPerAxis^d,
172
173 // init output
174 for (int face=0; face<2*Dimensions; face++) {
175 maxEigenvaluePerFaceOut[ face ] = 0.0;
176 }
177 // compute
178 // scalarIndexHull = face*2*scalarIndexFace
179 for ( int scalarIndexHull = 0; scalarIndexHull < collocatedSpaceTimeNodesOnCellHull; scalarIndexHull++ ) {
180 const int face = scalarIndexHull / (2*strideQLR);
181 // face = Dimensions*orientation + direction, direction < Dimensions
182 const int orientation = face/Dimensions;
183 const int direction = face-Dimensions*orientation;
184 // strip off face index info to get DoF index on face
185 const int scalarIndexFace = scalarIndexHull - face*2*strideQLR;
186
187 tarch::la::Vector<Dimensions,double> faceCentre = cellCentre;
188 faceCentre[ direction ] += (-1+2*orientation)*0.5*dx;
189 const double smax = riemann_maxAbsoluteEigenvalue_body_AoS(
190 maxAbsoluteEigenvalue,
191 QHullIn[ face ],
192 nodes,
193 faceCentre,
194 dx,
195 t,
196 dt,
197 nodesPerAxis,
198 unknowns,
199 strideQ,
200 direction,
201 scalarIndexFace);
202 maxEigenvaluePerFaceOut[ face ] = std::max( maxEigenvaluePerFaceOut[ face ], smax );
203 }
204}
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
GPUCallableMethod tarch::la::Vector< Dimensions+1, int > delineariseIndex(int scalarIndex, const tarch::la::Vector< Dimensions+1, int > &strides)
void riemann_setBoundaryState_body_AoS(std::function< void(double *__restrict__ QInside, double *__restrict__ OOutside, const tarch::la::Vector< Dimensions, double > &x, double t, int normal) > boundaryState, double *__restrict__ QOut, double *__restrict__ QIn, const double *const __restrict__ nodes, const tarch::la::Vector< Dimensions, double > &faceCentre, const double dx, const double t, const double dt, const int nodesPerAxis, const int unknowns, const int strideQ, const int direction, const int scalarIndexFace)
Set boundary state at given coordinates and time.
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 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:
double riemann_maxAbsoluteEigenvalue_body_AoS(std::function< double(const double *const __restrict__ Q, const tarch::la::Vector< Dimensions, double > &x, double t, const int normal) > maxAbsoluteEigenvalue, const double *const __restrict__ QLR, const double *const __restrict__ nodes, const tarch::la::Vector< Dimensions, double > &faceCentre, const double dx, const double t, const double dt, const int nodesPerAxis, const int unknowns, const int strideQ, const int direction, const int scalarIndexFace)
Determines the maximum absolute value among the eigenvalues computed at a space-time quadrature point...
void riemann_setBoundaryState_loop_AoS(std::function< void(const double *const __restrict__ QInside, double *__restrict__ OOutside, const tarch::la::Vector< Dimensions, double > &x, double t, int normal) > boundaryState, double *QHullOut[Dimensions *2], const double *const __restrict__ nodes, const tarch::la::Vector< Dimensions, double > &cellCentre, const double dx, const double t, const double dt, const int order, const int unknowns, const int auxiliaryVariables, const bool atBoundary[Dimensions *2])
Prescribe boundary values (QOut) and overwrite interior values (QInside).
GPUCallableMethod int getNodesPerCell(int nodesPerAxis)
void riemann_maxAbsoluteEigenvalue_loop_AoS(std::function< double(const double *const __restrict__ Q, const tarch::la::Vector< Dimensions, double > &x, double t, const int direction) > maxAbsoluteEigenvalue, double maxEigenvaluePerFaceOut[Dimensions *2], double *const QHullIn[Dimensions *2], const double *const __restrict__ nodes, const tarch::la::Vector< Dimensions, double > &cellCentre, const double dx, const double t, const double dt, const int order, const int unknowns, const int auxiliaryVariables)
Per face, determine the eigenvalue with the largest absolute value among the boundary-extrapolated pr...
Simple vector class.
Definition Vector.h:134