Peano 4
Loading...
Searching...
No Matches
computeGradients.cpph
Go to the documentation of this file.
1
13#ifndef _EXAHYPE_KERNELS_ADERDG_GENERIC_COMPUTEGRADIENTES_H_
14#define _EXAHYPE_KERNELS_ADERDG_GENERIC_COMPUTEGRADIENTES_H_
15
16#include <algorithm> // fill_n
19
20
21namespace kernels {
22namespace aderdg {
23namespace generic {
24namespace c {
25
26/* A helper to determine basisSize^Dimensions */
27inline constexpr int basisSizeD(int basisSize) {
28 //const int basisX = basisSize;
29 //const int basisY = basisSize;
30 //const int basisZ = Dimensions == 3 ? basisSize : 1;
31 //return basisX * basisY * basisZ;
32 return basisSize * basisSize * (Dimensions == 3 ? basisSize : 1);
33}
34
68inline void computeGradQ(double* gradQ, const double* const u, double*** dudx, const tarch::la::Vector<Dimensions, double>& sizeOfPatch, const int numberOfVariables, const int order) {
69 const int basisSize = order + 1;
70
71 const int basisX = basisSize;
72 const int basisY = basisSize;
73 const int basisZ = Dimensions == 3 ? basisSize : 1;
74
75 // spatial gradient of q: gradQ(z,y,x,nDim,nVar)
76 index idx_gradQ(basisZ, basisY, basisX, Dimensions, numberOfVariables);
77 index idx_u(basisZ, basisY, basisX, numberOfVariables);
78 std::fill_n(gradQ, idx_gradQ.size, 0.0);
79
80 // x direction (independent from the y and z derivatives)
81 for (int j = 0; j < basisZ; j++) { // z
82 for (int k = 0; k < basisY; k++) { // y
83 // Matrix operation
84 for (int l = 0; l < basisX; l++) { // x
85 for (int m = 0; m < numberOfVariables; m++) {
86 for (int n = 0; n < basisSize; n++) { // matmul x
87 gradQ[idx_gradQ(j, k, l, /*x*/0, m)] += 1.0 / sizeOfPatch[0] *
88 u[idx_u(j, k, n, m)] * dudx[order][l][n];
89 }
90 }
91 }
92 }
93 }
94
95 // y direction (independent from the x and z derivatives)
96 for (int j = 0; j < basisZ; j++) { // z
97 for (int k = 0; k < basisX; k++) { // x
98 // Matrix operation
99 for (int l = 0; l < basisY; l++) { // y
100 for (int m = 0; m < numberOfVariables; m++) {
101 for (int n = 0; n < basisSize; n++) { // matmul y
102 gradQ[idx_gradQ(j, l, k, /*y*/1, m)] += 1.0 / sizeOfPatch[1] *
103 u[idx_u(j, n, k, m)] * dudx[order][l][n];
104 }
105 }
106 }
107 }
108 }
109
110 // z direction (independent from the x and y derivatives)
111 if(Dimensions == 3) {
112 for (int j = 0; j < basisY; j++) { // y
113 for (int k = 0; k < basisX; k++) { // x
114 // Matrix operation
115 for (int l = 0; l < basisZ; l++) { // z
116 for (int m = 0; m < numberOfVariables; m++) {
117 for (int n = 0; n < basisSize; n++) { // matmul z
118 gradQ[idx_gradQ(l, j, k, /*z*/2, m)] += 1.0 / sizeOfPatch[2] *
119 u[idx_u(n, j, k, m)] * dudx[order][l][n];
120 }
121 }
122 }
123 }
124 }
125 }
126} // computeGradQ
127
139template <typename SolverType>
140inline void computeGradQi(double* gradQi, const double* const u, double*** dudx, const int requestedDirection, const int requestedVariableIndex, const tarch::la::Vector<Dimensions, double>& sizeOfPatch, const int numberOfVariables, const int order) {
141 const int basisSize = order + 1;
142 const int basisX = basisSize;
143 const int basisY = basisSize;
144 const int basisZ = Dimensions == 3 ? basisSize : 1;
145
146 index idx_gradQi(basisZ, basisY, basisX);
147 index idx_u(basisZ, basisY, basisX, numberOfVariables);
148 std::fill_n(gradQi, idx_gradQi.size, 0.0);
149
150 // x direction (independent from the y and z derivatives)
151 if(requestedDirection == 0) {
152 for (int j = 0; j < basisZ; j++) { // z
153 for (int k = 0; k < basisY; k++) { // y
154 // Matrix operation
155 for (int l = 0; l < basisX; l++) { // x
156 for (int n = 0; n < basisSize; n++) { // matmul x
157 gradQi[idx_gradQi(j, k, l)] += 1.0 / sizeOfPatch[0] *
158 u[idx_u(j, k, n, requestedVariableIndex)] * SolverType::dudx[order][l][n];
159 }
160 }
161 }
162 }
163 return;
164 }
165
166 // y direction (independent from the x and z derivatives)
167 if(requestedDirection == 2) {
168 for (int j = 0; j < basisZ; j++) { // z
169 for (int k = 0; k < basisX; k++) { // x
170 // Matrix operation
171 for (int l = 0; l < basisY; l++) { // y
172 for (int n = 0; n < basisSize; n++) { // matmul y
173 gradQi[idx_gradQi(j, l, k)] += 1.0 / sizeOfPatch[1] *
174 u[idx_u(j, n, k, requestedVariableIndex)] * SolverType::dudx[order][l][n];
175 }
176 }
177 }
178 }
179 return;
180 }
181
182 // z direction (independent from the x and y derivatives)
183 if(requestedDirection == 3) {
184 for (int j = 0; j < basisY; j++) { // y
185 for (int k = 0; k < basisX; k++) { // x
186 // Matrix operation
187 for (int l = 0; l < basisZ; l++) { // z
188 for (int n = 0; n < basisSize; n++) { // matmul z
189 gradQi[idx_gradQi(l, j, k)] += 1.0 / sizeOfPatch[2] *
190 u[idx_u(n, j, k, requestedVariableIndex)] * SolverType::dudx[order][l][n];
191 }
192 }
193 }
194 }
195 return;
196 }
197} // computeGradQi
198
202template <typename SolverType>
203void computeGradQ(double* gradQ, const double* const u, const tarch::la::Vector<Dimensions, double>& sizeOfPatch) {
204 const int numberOfVariables = SolverType::NumberOfVariables;
205 const int order = SolverType::Order;
206 computeGradQ(gradQ, u, kernels::legendre::dudx, sizeOfPatch, numberOfVariables, order);
207}
208
210template <typename SolverType>
211void computeGradQi(double* gradQi, const double* const u, const int requestedDirection, const int requestedVariableIndex, const tarch::la::Vector<Dimensions, double>& sizeOfPatch) {
212 const int numberOfVariables = SolverType::NumberOfVariables;
213 const int order = SolverType::Order;
214 computeGradQi<SolverType>(gradQi, u, kernels::legendre::dudx, requestedDirection, requestedVariableIndex, sizeOfPatch, numberOfVariables, order);
215}
216
217
218} // c
219} // generic
220} // aderdg
221} // kernels
222
223#endif
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
-with-ExaHyPE2-benchmark
constexpr int basisSizeD(int basisSize)
void computeGradQ(double *gradQ, const double *const u, const tarch::la::Vector< DIMENSIONS, double > &sizeOfPatch, const int numberOfVariables, const int order)
Helper function to compute the gradients for all entries in Q with the ADER-DG scheme.
void computeGradQi(double *gradQi, const double *const u, const int requestedDirection, const int requestedVariableIndex, const tarch::la::Vector< DIMENSIONS, double > &sizeOfPatch, const int numberOfVariables, const int order)
Convenience function to compute the gradient of only a single variable in a given direction.
This file is part of the ExaHyPE project.
This is a single successor class for the idx2, idx3, idx4, idx5, idx6 classes.
Definition KernelUtils.h:62
const int size
Definition KernelUtils.h:66
Simple vector class.
Definition Vector.h:134