Peano 4
Loading...
Searching...
No Matches
computeGradients.cpph
Go to the documentation of this file.
1
14#include <algorithm> // fill_n
15#include "kernels/GaussLegendreBasis.h"
16#include "kernels/KernelUtils.h"
17
18
19namespace kernels {
20namespace aderdg {
21namespace generic {
22namespace c {
23
24/* A helper to determine basisSize^DIMENSIONS */
25inline constexpr int basisSizeD(int basisSize) {
26 //const int basisX = basisSize;
27 //const int basisY = basisSize;
28 //const int basisZ = DIMENSIONS == 3 ? basisSize : 1;
29 //return basisX * basisY * basisZ;
30 return basisSize * basisSize * (DIMENSIONS == 3 ? basisSize : 1);
31}
32
65inline void computeGradQ(double* gradQ, const double* const u, const tarch::la::Vector<DIMENSIONS, double>& sizeOfPatch, const int numberOfVariables, const int order) {
66 const int basisSize = order + 1;
67
68 const int basisX = basisSize;
69 const int basisY = basisSize;
70 const int basisZ = DIMENSIONS == 3 ? basisSize : 1;
71
72 // spatial gradient of q: gradQ(z,y,x,nDim,nVar)
73 index idx_gradQ(basisZ, basisY, basisX, DIMENSIONS, numberOfVariables);
74 index idx_u(basisZ, basisY, basisX, numberOfVariables);
75 std::fill_n(gradQ, idx_gradQ.size, 0.0);
76
77 // x direction (independent from the y and z derivatives)
78 for (int j = 0; j < basisZ; j++) { // z
79 for (int k = 0; k < basisY; k++) { // y
80 // Matrix operation
81 for (int l = 0; l < basisX; l++) { // x
82 for (int m = 0; m < numberOfVariables; m++) {
83 for (int n = 0; n < basisSize; n++) { // matmul x
84 gradQ[idx_gradQ(j, k, l, /*x*/0, m)] += 1.0 / sizeOfPatch[0] *
85 u[idx_u(j, k, n, m)] * exahype2::solvers::aderdg::ADERDGSolver::dudx[order][l][n];
86 }
87 }
88 }
89 }
90 }
91
92 // y direction (independent from the x and z derivatives)
93 for (int j = 0; j < basisZ; j++) { // z
94 for (int k = 0; k < basisX; k++) { // x
95 // Matrix operation
96 for (int l = 0; l < basisY; l++) { // y
97 for (int m = 0; m < numberOfVariables; m++) {
98 for (int n = 0; n < basisSize; n++) { // matmul y
99 gradQ[idx_gradQ(j, l, k, /*y*/1, m)] += 1.0 / sizeOfPatch[1] *
100 u[idx_u(j, n, k, m)] * exahype2::solvers::aderdg::ADERDGSolver::dudx[order][l][n];
101 }
102 }
103 }
104 }
105 }
106
107 // z direction (independent from the x and y derivatives)
108 if(DIMENSIONS == 3) {
109 for (int j = 0; j < basisY; j++) { // y
110 for (int k = 0; k < basisX; k++) { // x
111 // Matrix operation
112 for (int l = 0; l < basisZ; l++) { // z
113 for (int m = 0; m < numberOfVariables; m++) {
114 for (int n = 0; n < basisSize; n++) { // matmul z
115 gradQ[idx_gradQ(l, j, k, /*z*/2, m)] += 1.0 / sizeOfPatch[2] *
116 u[idx_u(n, j, k, m)] * exahype2::solvers::aderdg::ADERDGSolver::dudx[order][l][n];
117 }
118 }
119 }
120 }
121 }
122 }
123} // computeGradQ
124
136inline 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) {
137 const int basisSize = order + 1;
138 const int basisX = basisSize;
139 const int basisY = basisSize;
140 const int basisZ = DIMENSIONS == 3 ? basisSize : 1;
141
142 index idx_gradQi(basisZ, basisY, basisX);
143 index idx_u(basisZ, basisY, basisX, numberOfVariables);
144 std::fill_n(gradQi, idx_gradQi.size, 0.0);
145
146 // x direction (independent from the y and z derivatives)
147 if(requestedDirection == 0) {
148 for (int j = 0; j < basisZ; j++) { // z
149 for (int k = 0; k < basisY; k++) { // y
150 // Matrix operation
151 for (int l = 0; l < basisX; l++) { // x
152 for (int n = 0; n < basisSize; n++) { // matmul x
153 gradQi[idx_gradQi(j, k, l)] += 1.0 / sizeOfPatch[0] *
154 u[idx_u(j, k, n, requestedVariableIndex)] * exahype2::solvers::aderdg::ADERDGSolver::dudx[order][l][n];
155 }
156 }
157 }
158 }
159 return;
160 }
161
162 // y direction (independent from the x and z derivatives)
163 if(requestedDirection == 2) {
164 for (int j = 0; j < basisZ; j++) { // z
165 for (int k = 0; k < basisX; k++) { // x
166 // Matrix operation
167 for (int l = 0; l < basisY; l++) { // y
168 for (int n = 0; n < basisSize; n++) { // matmul y
169 gradQi[idx_gradQi(j, l, k)] += 1.0 / sizeOfPatch[1] *
170 u[idx_u(j, n, k, requestedVariableIndex)] * exahype2::solvers::aderdg::ADERDGSolver::dudx[order][l][n];
171 }
172 }
173 }
174 }
175 return;
176 }
177
178 // z direction (independent from the x and y derivatives)
179 if(requestedDirection == 3) {
180 for (int j = 0; j < basisY; j++) { // y
181 for (int k = 0; k < basisX; k++) { // x
182 // Matrix operation
183 for (int l = 0; l < basisZ; l++) { // z
184 for (int n = 0; n < basisSize; n++) { // matmul z
185 gradQi[idx_gradQi(l, j, k)] += 1.0 / sizeOfPatch[2] *
186 u[idx_u(n, j, k, requestedVariableIndex)] * exahype2::solvers::aderdg::ADERDGSolver::dudx[order][l][n];
187 }
188 }
189 }
190 }
191 return;
192 }
193} // computeGradQi
194
198template <typename SolverType>
199void computeGradQ(double* gradQ, const double* const u, const tarch::la::Vector<DIMENSIONS, double>& sizeOfPatch) {
200 const int numberOfVariables = SolverType::NumberOfVariables;
201 const int order = SolverType::Order;
202 computeGradQ(gradQ, u, sizeOfPatch, numberOfVariables, order);
203}
204
206template <typename SolverType>
207void computeGradQi(double* gradQi, const double* const u, const int requestedDirection, const int requestedVariableIndex, const tarch::la::Vector<DIMENSIONS, double>& sizeOfPatch) {
208 const int numberOfVariables = SolverType::NumberOfVariables;
209 const int order = SolverType::Order;
210 computeGradQi(gradQi, u, requestedDirection, requestedVariableIndex, sizeOfPatch, numberOfVariables, order);
211}
212
213
214} // c
215} // generic
216} // aderdg
217} // kernels
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