Peano 4
Loading...
Searching...
No Matches
CellIntegral.cpph
Go to the documentation of this file.
1#include "peano4/utils/Loop.h"
3
4
5template <
6 typename Solver,
7 int order,
8 int unknowns,
9 int auxiliaryVariables
10>
12 ::exahype2::CellData& cellData,
13 bool evaluateFlux,
14 bool evaluateNonconservativeProduct,
15 bool evaluateSource,
16 bool evaluatePointSources
17) {
19 1, // number of cells in one memory chunk
20 order+1, // numberOfDoFsPerAxisInCell
21 0, // no halo
23 auxiliaryVariables
24 );
25
27 1, // number of cells in one memory chunk
28 order+1, // numberOfDoFsPerAxisInCell
29 0, // no halo
31 0 // no auxiliary variables
32 );
33
34 for (int currentCell=0; currentCell<cellData.numberOfCells; currentCell++) {
35 const double* __restrict__ cellQin = cellData.QIn[currentCell];
36 double* __restrict__ cellQout = cellData.QOut[currentCell];
37 const tarch::la::Vector<Dimensions,double> x = cellData.cellCentre[currentCell];
38 const tarch::la::Vector<Dimensions,double> h = cellData.cellSize[currentCell];
39 const double t = cellData.t[currentCell];
40 const double dt = cellData.dt[currentCell];
41
42 const int nodesPerAxis = order+1;
43 const int strideQ=unknowns+auxiliaryVariables;
44 const int nodesPerCell = getNodesPerCell(nodesPerAxis);
45
46 dfor( dof, order+1 ) {
47 double fluxValues[unknowns];
48 double sourceValues[unknowns];
49 double gradientValues[(unknowns+auxiliaryVariables)];
50 double ncpValues[unknowns];
51
52 for(int var=0; var<unknowns; var++){
53 cellQout[enumeratorOut(currentCell,dof,var)] = 0.0;
54 }
55
56 if(evaluateFlux){
57 // direction in which flux contributions are being computed
58 for(int dim=0; dim<Dimensions; dim++)
59 for(int i=0; i<nodesPerAxis; i++){
60 // compute contributing node real position
61 tarch::la::Vector<Dimensions,int> contributingNodeIndex = dof;
62 contributingNodeIndex[dim] = i;
63 tarch::la::Vector<Dimensions,double> position = getQuadraturePoint(
64 x,
65 h,
66 contributingNodeIndex,
67 order,
68 Solver::QuadraturePoints1d);
69
70 Solver::flux(
71 &cellQin[enumeratorIn(currentCell, contributingNodeIndex, 0)],
72 position,
73 t,
74 dt,
75 dim, //normal
76 fluxValues,
77 Solver::Offloadable::Yes
78 );
79
80 // Coefficient within matrix is tensor product over dimensions, too
81 double coeff = Solver::StiffnessOperator1d[i+dof(dim)*nodesPerAxis];
82 for (int dd=0; dd<Dimensions; dd++) {
83 coeff *= (dd==dim) ? 1.0 : Solver::MassMatrixDiagonal1d[ dof(dd) ] * h(dd);
84 }
85
86 for(int var=0; var<unknowns; var++){
87 cellQout[enumeratorOut(currentCell, dof, var)] += coeff * fluxValues[var];
88 } //var
89 } // dim/i combination
90 }//if useFlux
91
92 if(evaluateSource){
93 //compute real node position
94 tarch::la::Vector<Dimensions,double> position = getQuadraturePoint(
95 x,
96 h,
97 dof,
98 order,
99 Solver::QuadraturePoints1d);
100
101 //compute source term contributions
102 Solver::sourceTerm(
103 &cellQin[enumeratorIn(currentCell, dof, 0)],
104 position,
105 t,
106 dt,
107 sourceValues,
108 Solver::Offloadable::Yes
109 );
110
111 #if Dimensions==2
112 double coeff = Solver::MassMatrixDiagonal1d[ dof(0) ] * Solver::MassMatrixDiagonal1d[ dof(1) ] * tarch::la::volume(h);
113 #else
114 double coeff = Solver::MassMatrixDiagonal1d[ dof(0) ] * Solver::MassMatrixDiagonal1d[ dof(1) ] * Solver::MassMatrixDiagonal1d[ dof(2) ] * tarch::la::volume(h);
115 #endif
116
117 for(int var=0; var<unknowns; var++){
118 cellQout[enumeratorOut(currentCell, dof, var)] += coeff * sourceValues[var];
119 } //var
120 } //if useSource
121
122 // @todo Point Source missing
123
124 if (evaluateNonconservativeProduct) {
125
126 //compute real node position
127 tarch::la::Vector<Dimensions,double> position = getQuadraturePoint( x,
128 h,
129 dof,
130 order,
131 Solver::QuadraturePoints1d);
132
133 double coeff = 1.0;
134 for (int dd=0; dd<Dimensions; dd++) {
135 coeff *= Solver::MassMatrixDiagonal1d[ dof(dd) ] * h(dd);
136 }
137
138 for(int dim=0; dim<Dimensions; dim++) {
139 //computes gradient values in all directions for the given node
140 std::fill(gradientValues, &gradientValues[strideQ-1], 0.0);
141
142 const int invDx = 1 / h(dim);
143
144 // computing gradients in direction dim
145 for(int node=0; node<nodesPerAxis; node++){
146 tarch::la::Vector<Dimensions,int> contributingNodeIndex = dof;
147 contributingNodeIndex[dim] = node;
148 const double coeffDerivative = invDx * Solver::DerivativeOperator1d[node + nodesPerAxis*dof(dim)];
149
150 for(int var=0; var<strideQ; var++){
151 gradientValues[var] += coeffDerivative * cellQin[enumeratorIn(currentCell, contributingNodeIndex, var)];
152 assertion(gradientValues[var]==gradientValues[var]);
153 } // var
154
155 } //node
156
157 Solver::nonconservativeProduct(
158 &cellQin[enumeratorIn(currentCell, dof, 0)],
159 gradientValues,
160 position,
161 t,
162 dt,
163 dim,
164 ncpValues,
165 Solver::Offloadable::Yes
166 );
167
168 for(int var=0; var<unknowns; var++){
169 assertion(ncpValues[var]==ncpValues[var]);
170 cellQout[enumeratorOut(currentCell, dof, var)] -= coeff * ncpValues[var];
171 } //var
172 } //dim
173 } //if useNCP
174 } // dof
175 } //currentCell
176}
#define assertion(expr)
#define dfor(counter, max)
d-dimensional Loop
Definition Loop.h:308
void cellIntegral_patchwise_in_situ_GaussLegendre(::exahype2::CellData &cellData, bool evaluateFlux, bool evaluateNonconservativeProduct, bool evaluateSource, bool evaluatePointSources)
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 ** 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