Peano 4
Loading...
Searching...
No Matches
CellIntegral.cpp
Go to the documentation of this file.
1#include "CellIntegral.h"
2
3#include "peano4/utils/Loop.h"
4
6
8 ::exahype2::CellData& cellData,
9 const int order,
10 const int unknowns,
11 const int auxiliaryVariables,
12 Flux flux,
13 NonConservativeProduct nonconservativeProduct,
14 Source source,
15 PointSources pointSources,
16 const double* __restrict__ quadratureNodes,
17 const double* __restrict__ massMatrixDiagonal,
18 const double* __restrict__ stiffnessMatrix,
19 const double* __restrict__ derivativeOperator,
20 bool evaluateFlux,
21 bool evaluateNonconservativeProduct,
22 bool evaluateSource,
23 bool evaluatePointSources
24) {
26 1, // number of cells in one memory chunk
27 order+1, // numberOfDoFsPerAxisInCell
28 0, // no halo
30 auxiliaryVariables
31 );
32
34 1, // number of cells in one memory chunk
35 order+1, // numberOfDoFsPerAxisInCell
36 0, // no halo
38 0 // no auxiliary variables
39 );
40
41 double* fluxValues = evaluateFlux ? new double[unknowns] : nullptr;
42 double* sourceValues = evaluateSource ? new double[unknowns] : nullptr;
43 double* gradientValues = evaluateNonconservativeProduct ? new double[(unknowns+auxiliaryVariables)] : nullptr;
44 double* ncpValues = evaluateNonconservativeProduct ? new double[unknowns] : nullptr;
45
46 for (int currentCell=0; currentCell<cellData.numberOfCells; currentCell++) {
47 const double* __restrict__ cellQin = cellData.QIn[currentCell];
48 double* __restrict__ cellQout = cellData.QOut[currentCell];
49 const tarch::la::Vector<Dimensions,double> x = cellData.cellCentre[currentCell];
50 const tarch::la::Vector<Dimensions,double> h = cellData.cellSize[currentCell];
51 const double t = cellData.t[currentCell];
52 const double dt = cellData.dt[currentCell];
53
54 const int nodesPerAxis = order+1;
55 const int strideQ=unknowns+auxiliaryVariables;
56 const int nodesPerCell = getNodesPerCell(nodesPerAxis);
57
58 dfor( dof, order+1 ) {
59 for(int var=0; var<unknowns; var++){
60 cellQout[enumeratorOut(currentCell,dof,var)] = 0.0;
61 }
62
63 if(evaluateFlux){
64 // direction in which flux contributions are being computed
65 for(int dim=0; dim<Dimensions; dim++)
66 for(int i=0; i<nodesPerAxis; i++){
67 // compute contributing node real position
68 tarch::la::Vector<Dimensions,int> contributingNodeIndex = dof;
69 contributingNodeIndex[dim] = i;
70 tarch::la::Vector<Dimensions,double> position = getQuadraturePoint(
71 x,
72 h,
73 contributingNodeIndex,
74 order,
75 quadratureNodes);
76
77 flux( &cellQin[enumeratorIn(currentCell, contributingNodeIndex, 0)],
78 position,
79 t,
80 dt,
81 dim, //normal
82 fluxValues);
83
84 // Coefficient within matrix is tensor product over dimensions, too.
85 // So we take the stiffness matrix in the direction that we are
86 // evaluating and multiply it for the d-1 remaining directions with
87 // the mass matrix. We integrate over the whole cell, but the h-terms
88 // are eliminated along the derivative direction.
89 double coeff = stiffnessMatrix[i+dof(dim)*nodesPerAxis];
90 for (int dd=0; dd<Dimensions; dd++) {
91 coeff *= (dd==dim) ? 1.0 : massMatrixDiagonal[ dof(dd) ] * h(dd);
92 }
93
94 for(int var=0; var<unknowns; var++){
95 cellQout[enumeratorOut(currentCell, dof, var)] += coeff * fluxValues[var];
96 } //var
97 } // dim/i combination
98 }//if useFlux
99
100 if(evaluateSource){
101 //compute real node position
102 tarch::la::Vector<Dimensions,double> position = getQuadraturePoint(
103 x,
104 h,
105 dof,
106 order,
107 quadratureNodes
108 );
109
110 //compute source term contributions
111 source(
112 &cellQin[enumeratorIn(currentCell, dof, 0)],
113 position,
114 t,
115 dt,
116 sourceValues
117 );
118
119 #if Dimensions==2
120 double coeff = massMatrixDiagonal[ dof(0) ] * massMatrixDiagonal[ dof(1) ] * tarch::la::volume(h);
121 #else
122 double coeff = massMatrixDiagonal[ dof(0) ] * massMatrixDiagonal[ dof(1) ] * massMatrixDiagonal[ dof(2) ] * tarch::la::volume(h);
123 #endif
124
125 for(int var=0; var<unknowns; var++){
126 cellQout[enumeratorOut(currentCell, dof, var)] += coeff * sourceValues[var];
127 } //var
128 } //if useSource
129
130 // @todo Point Source missing
131
132 if (evaluateNonconservativeProduct) {
133 //compute real node position
134 tarch::la::Vector<Dimensions,double> position = getQuadraturePoint(
135 x,
136 h,
137 dof,
138 order,
139 quadratureNodes
140 );
141
142 // Mass matrix.
143 double coeff = 1.0;
144 for (int dd=0; dd<Dimensions; dd++) {
145 coeff *= massMatrixDiagonal[ dof(dd) ] * h(dd);
146 }
147
148 for(int dim=0; dim<Dimensions; dim++) {
149 //computes gradient values in all directions for the given node
150 std::fill(gradientValues, &gradientValues[strideQ-1], 0.0);
151
152 const int invDx = 1 / h(dim);
153
154 // computing gradients in direction dim
155 for(int node=0; node<nodesPerAxis; node++){
156 tarch::la::Vector<Dimensions,int> contributingNodeIndex = dof;
157 contributingNodeIndex[dim] = node;
158 const double coeffDerivative = invDx * derivativeOperator[node + nodesPerAxis*dof(dim)];
159
160 for(int var=0; var<strideQ; var++){
161 gradientValues[var] += coeffDerivative * cellQin[enumeratorIn(currentCell, contributingNodeIndex, var)];
162 assertion(gradientValues[var]==gradientValues[var]);
163 } // var
164 } //node
165
166 nonconservativeProduct(
167 &cellQin[enumeratorIn(currentCell, dof, 0)],
168 gradientValues,
169 position,
170 t,
171 dt,
172 dim,
173 ncpValues);
174
175 for(int var=0; var<unknowns; var++){
176 assertion(ncpValues[var]==ncpValues[var]);
177 cellQout[enumeratorOut(currentCell, dof, var)] -= coeff * ncpValues[var];
178 } //var
179 } //dim
180 } //if useNCP
181 } // dof
182 } //currentCell
183
184 if (evaluateFlux) {
185 delete [] fluxValues;
186 }
187 if (evaluateSource) {
188 delete [] sourceValues;
189 }
190 if (evaluateNonconservativeProduct) {
191 delete [] gradientValues;
192 delete [] ncpValues;
193 }
194}
195
196
198 ::exahype2::CellData& cellData,
199 const int order,
200 const int unknowns,
201 const int auxiliaryVariables,
202 MaxEigenvalue maxEigenvalue,
203 const double* __restrict__ QuadratureNodes1d
204) {
206 1, // number of patches stored per continuous memory segment
207 order+1,
208 0, // halo
209 unknowns,
210 auxiliaryVariables
211 );
212
213 for (int cellIndex=0; cellIndex<cellData.numberOfCells; cellIndex++) {
214 double eigenvalue = std::numeric_limits<double>::min();
215
216 dfor(k,order+1) {
217 for(int direction=0; direction<Dimensions; direction++) {
218 double nodeEigenvalue = maxEigenvalue(
219 cellData.QOut[cellIndex] + enumerator(0,k,0),
220 getQuadraturePoint(
221 cellData.cellCentre[cellIndex],
222 cellData.cellSize[cellIndex],
223 k,
224 order,
225 QuadratureNodes1d
226 ),
227 cellData.t[cellIndex],
228 cellData.dt[cellIndex],
229 direction
230 );
231
232 assertion( nodeEigenvalue>=0.0 );
233
234 eigenvalue = std::max(eigenvalue, nodeEigenvalue);
235 }
236 }
237
238 cellData.maxEigenvalue[cellIndex] = eigenvalue;
239 }
240}
241
242
244 const double* __restrict__ cellQin,
245 const int order,
246 const int unknowns,
247 const int auxiliaryVariables,
249 double* __restrict__ cellQout
250) {
251 const int nodeSerialised = peano4::utils::dLinearised( node, order+1 );
252 for(int var=0; var<unknowns+auxiliaryVariables; var++) {
253 cellQout[nodeSerialised+var] = cellQin[nodeSerialised+var];
254 }
255}
256
257
259 ::exahype2::CellData& cellData,
260 const int order,
261 const int unknowns,
262 const int auxiliaryVariables,
263 const double* __restrict__ massMatrixDiagonal1d
264) {
266 1, // number of cells in one memory chunk
267 order+1, // numberOfDoFsPerAxisInCell
268 0, // no halo
269 unknowns,
270 0 // no auxiliary variables here
271 );
272
273 for (int currentCell=0; currentCell<cellData.numberOfCells; currentCell++) {
274 dfor( dof, order+1 ) {
275 double* __restrict__ cellQInOut = cellData.QOut[currentCell];
276
277 #if Dimensions==2
278 double diagonal = massMatrixDiagonal1d[ dof(0) ] * massMatrixDiagonal1d[ dof(1) ] * (cellData.cellSize[currentCell])(0) * (cellData.cellSize[currentCell])(1);
279 #else
280 double diagonal = massMatrixDiagonal1d[ dof(0) ]
281 * massMatrixDiagonal1d[ dof(1) ]
282 * massMatrixDiagonal1d[ dof(2) ]
283 * (cellData.cellSize[currentCell])(0)
284 * (cellData.cellSize[currentCell])(1)
285 * (cellData.cellSize[currentCell])(2);
286 #endif
287
288 assertion( diagonal>0.0 );
289
290 for(int var=0; var<unknowns; var++) {
291 cellQInOut[ enumerator(currentCell,dof,var) ] = 1.0 / diagonal * cellQInOut[enumerator(currentCell,dof,var)];
292 }
293 }
294 }
295}
#define assertion(expr)
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
#define dfor(counter, max)
d-dimensional Loop
Definition Loop.h:308
void copySolution(const double *__restrict__ cellQin, const int order, const int unknowns, const int auxiliaryVariables, const tarch::la::Vector< Dimensions, int > &node, double *__restrict__ cellQout)
Copy solution over for one node.
void cellIntegral_patchwise_in_situ_GaussLegendre_functors(::exahype2::CellData &cellData, const int order, const int unknowns, const int auxiliaryVariables, Flux flux, NonConservativeProduct nonconservativeProduct, Source source, PointSources pointSources, const double *__restrict__ QuadratureNodes1d, const double *__restrict__ MassMatrixDiagonal1d, const double *__restrict__ StiffnessMatrix1d, const double *__restrict__ DerivativeOperator1d, bool evaluateFlux, bool evaluateNonconservativeProduct, bool evaluateSource, bool evaluatePointSources)
void reduceMaxEigenvalue_patchwise_functors(::exahype2::CellData &cellData, const int order, const int unknowns, const int auxiliaryVariables, MaxEigenvalue maxEigenvalue, const double *__restrict__ QuadratureNodes1d)
Compute the maximum eigenvalues over a sequence of cells and store the result in the respective CellD...
void multiplyWithInvertedMassMatrix_GaussLegendre(::exahype2::CellData &cellData, const int order, const int unknowns, const int auxiliaryVariables, const double *__restrict__ MassMatrixDiagonal1d)
Final step of DG algorithm.
CPUGPUMethod int dLinearised(const tarch::la::Vector< Dimensions, int > &counter, int max)
Map d-dimensional vector onto integer index.
Definition Loop.cpp:106
Scalar volume(const Vector< Size, Scalar > &vector)
Computes the volume of the tetrahedron spanned by the Cartesian unit vectors scaled by the correspond...
Representation of a number of cells which contains all information that's required to process the sto...
Definition CellData.h:79
double * maxEigenvalue
Out values.
Definition CellData.h:111
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
Simple vector class.
Definition Vector.h:134