Peano 4
Loading...
Searching...
No Matches
CorrectorAoS.cpp
Go to the documentation of this file.
1#include "CorrectorAoS.h"
2
3#include "KernelUtils.h"
4
5#include "Generic.h"
6
7#include "tarch/la/Vector.h"
8
9#include <iostream>
10
11/*
12#if defined(GPUOffloadingOMP)
13#pragma omp declare target
14#endif
15*/
17 std::function< void(
18 double * __restrict__ Q,
20 double t
21 ) > adjustSolution,
22 std::function< double(
23 const double * const __restrict__ Q,
25 double t,
26 const int normal
27 ) > maxAbsoluteEigenvalue,
28 double * __restrict__ UOut,
29 const double * const __restrict__ nodes,
31 const double dx,
32 const double t,
33 const int nodesPerAxis,
34 const int unknowns,
35 const int strideQ,
36 const int callMaxEigenvalue,
37 const int scalarIndex) {
38 const tarch::la::Vector<Dimensions+1,int> index = delineariseIndex(scalarIndex,getStrides(nodesPerAxis,false));
39 const tarch::la::Vector<Dimensions+1, double> coords = getCoordinates(index,cellCentre,dx,t,0,nodes);
40 const tarch::la::Vector<Dimensions, double> x( ( coords.data() + 1 ) );
41
42 double * __restrict__ Q = &UOut[ scalarIndex * strideQ];
43 adjustSolution(Q, x, t);
44
45 if ( callMaxEigenvalue ) {
46 double lambdaMax = 0.0;
47 for ( int d = 0; d < Dimensions; d++ ) {
48 lambdaMax = std::max( lambdaMax, maxAbsoluteEigenvalue(Q,x,t,d) );
49 }
50 return lambdaMax;
51 } else {
52 return 0.0;
53 }
54}
55/*
56#if defined(GPUOffloadingOMP)
57#pragma omp end declare target
58#endif
59
60#if defined(GPUOffloadingOMP)
61#pragma omp declare target
62#endif
63*/
65 std::function< void(
66 const double * const __restrict__ Q,
68 double t,
69 int normal,
70 double * __restrict__ F
71 ) > flux,
72 double* __restrict__ UOut,
73 const double* __restrict__ QIn,
74 double* __restrict__ FAux, // must be allocated per thread as size is runtime parameter
75 const double* __restrict__ nodes,
76 const double* __restrict__ weights,
77 const double* __restrict__ Kxi,
79 const double dx,
80 const double t,
81 const double dt,
82 const int nodesPerAxis,
83 const int unknowns,
84 const int strideQ,
85 const int scalarIndex) {
86 tarch::la::Vector<Dimensions+1,int> strides = getStrides(nodesPerAxis,false);
87 tarch::la::Vector<Dimensions+1,int> index = delineariseIndex(scalarIndex,strides);
88 const tarch::la::Vector<Dimensions+1, double> coords = getCoordinates(index,cellCentre,dx,t,dt,nodes);
89 const tarch::la::Vector<Dimensions, double> x( ( coords.data() + 1 ) );
90
91 const double coeff0 = dt/dx/*[d]*/;
92 for ( int d = 0; d < Dimensions; d++) { // direction
93 const double coeff1 = coeff0 / weights[index[d+1]];
94
95 for ( int id = 0; id < nodesPerAxis; id++ ) { // loop over spatial neighbours
96 const double coeff2 = coeff1 * Kxi[ index[d+1]*nodesPerAxis + id ]; // note: transposed vs. predictor flux computation
97
98 const int scalarIndexNeighbour = scalarIndex + (id - index[d+1])*strides[d+1];
99 for ( int it = 0; it < nodesPerAxis; it++ ) { // time loop
100 const double coeff = coeff2 * weights[it];
101
102 const double* Q = &QIn[ (scalarIndexNeighbour*nodesPerAxis + it)*strideQ ];
103 const double time = t + nodes[it] * dt;
104 flux( Q, x, time, d, FAux );
105 for (int var=0; var < unknowns; var++) {
106 UOut[ scalarIndex*strideQ + var ] += coeff * FAux[ var ];
107 }
108 }
109 }
110 }
111}
112/*
113#if defined(GPUOffloadingOMP)
114#pragma omp end declare target
115#endif
116
117#if defined(GPUOffloadingOMP)
118#pragma omp declare target
119#endif
120*/
122 std::function< void(
123 const double * const __restrict__ Q,
125 double t,
126 double * __restrict__ S
127 ) > algebraicSource,
128 double * __restrict__ UOut,
129 double * __restrict__ SAux,
130 const double * const __restrict__ QIn,
131 const double * const __restrict__ nodes,
132 const double * const __restrict__ weights,
133 const tarch::la::Vector<Dimensions,double>& cellCentre,
134 const double dx,
135 const double t,
136 const double dt,
137 const int nodesPerAxis,
138 const int unknowns,
139 const int strideQ,
140 const int scalarIndex) {
141 const tarch::la::Vector<Dimensions+1,int> index = delineariseIndex(scalarIndex,getStrides(nodesPerAxis,false));
142 const tarch::la::Vector<Dimensions+1, double> coords = getCoordinates(index,cellCentre,dx,t,dt,nodes);
143 const tarch::la::Vector<Dimensions, double> x( ( coords.data() + 1 ) );
144
145 for ( int it=0; it < nodesPerAxis; it++ ) { // time-integral
146 const double coeff = dt * weights[it];
147
148 const double* Q = &QIn[ (scalarIndex*nodesPerAxis + it)*strideQ ];
149 const double time = t + nodes[it] * dt;
150 algebraicSource(Q, x, time, SAux);
151
152 for (int var = 0; var < unknowns; var++) {
153 UOut[ scalarIndex*strideQ + var ] += coeff * SAux[var];
154 }
155 }
156}
157/*
158#if defined(GPUOffloadingOMP)
159#pragma omp end declare target
160#endif
161
162#if defined(GPUOffloadingOMP)
163#pragma omp declare target
164#endif
165*/
167 std::function< void(
168 const double * const __restrict__ Q,
169 double * __restrict__ dQ_or_deltaQ,
171 double t,
172 int normal,
173 double * __restrict__ BgradQ
174 ) > nonconservativeProduct,
175 double * __restrict__ UOut,
176 double * __restrict__ gradQAux,
177 double * __restrict__ SAux,
178 const double * const __restrict__ QIn,
179 const double * const __restrict__ nodes,
180 const double * const __restrict__ weights,
181 const double * const __restrict__ dudx,
182 const tarch::la::Vector<Dimensions,double>& cellCentre,
183 const double dx,
184 const double t,
185 const double dt,
186 const int nodesPerAxis,
187 const int unknowns,
188 const int strideQ,
189 const int scalarIndex) {
190 const tarch::la::Vector<Dimensions+1,int> index = delineariseIndex(scalarIndex,getStrides(nodesPerAxis,false));
191 const tarch::la::Vector<Dimensions+1, double> coords = getCoordinates(index,cellCentre,dx,t,dt,nodes);
192 const tarch::la::Vector<Dimensions, double> x( ( coords.data() + 1 ) );
193
194 const double invDx = 1.0/dx;
195
196 for ( int it=0; it < nodesPerAxis; it++ ) { // time-integral
197 const double coeff = dt * weights[it];
198
199 gradient_AoS( QIn, dudx, invDx, nodesPerAxis, strideQ, scalarIndex*nodesPerAxis+it, gradQAux );
200 const double time = t + nodes[it] * dt;
201 const double* Q = &QIn[ (scalarIndex*nodesPerAxis + it)*strideQ ];
202 for ( int direction = 0; direction < Dimensions; direction++ ) {
203 nonconservativeProduct( Q, &gradQAux[ direction*strideQ ], x, time, direction, SAux );
204 for(int var=0; var<unknowns; var++) {
205 UOut[ scalarIndex*strideQ + var ] += coeff * SAux[var];
206 }
207 }
208 }
209}
210/*
211#if defined(GPUOffloadingOMP)
212#pragma omp end declare target
213#endif
214
215#if defined(GPUOffloadingOMP)
216#pragma omp declare target
217#endif
218*/
219
220#if defined(GPUOffloadingCPP)
222#else
224#endif
225 double * __restrict__ UOut,
226 const double * const __restrict__ riemannResultIn,
227 const double * const __restrict__ weights,
228 const double * const __restrict__ FLCoeff,
229 const double dx,
230 const double dt,
231 const int nodesPerAxis,
232 const int unknowns,
233 const int strideQ,
234 const int strideF,
235 const int scalarIndexCell) {
237 delineariseIndex(scalarIndexCell,getStrides(nodesPerAxis,false));
238
239 const double coeff0 = dt / dx/*[d]*/;
240 for (int d=0; d < Dimensions; d++) { // direction
241 const double coeff1 = coeff0 / weights[index[d+1]];
242
243 const double coeff_L = coeff1 * FLCoeff[index[d+1]];
244 const double coeff_R = coeff1 * FLCoeff[nodesPerAxis-1 - index[d+1]]; // FR is mirrored FL thanks to symmetric distrib. of nodes in [0,1]
245 //const double coeff_R = coeff1 * FRCoeff[index[d+1]];
246
247 const int scalarIndexHull_L = mapCellIndexToScalarHullIndex(index,d,0,nodesPerAxis);
248 const int scalarIndexHull_R = mapCellIndexToScalarHullIndex(index,d,1,nodesPerAxis);
249 // if the boundary data is symmetric, the contributions should cancel
250 // "left" minus "right" flux
251 for (int var=0; var < unknowns; var++) {
252 UOut[ scalarIndexCell*strideQ + var ] += coeff_L * riemannResultIn[ scalarIndexHull_L*strideF + var ];
253 }
254 for (int var=0; var < unknowns; var++) {
255 UOut[ scalarIndexCell*strideQ + var ] -= coeff_R * riemannResultIn[ scalarIndexHull_R*strideF + var ];
256 }
257 }
258}
259/*
260#if defined(GPUOffloadingOMP)
261#pragma omp end declare target
262#endif
263*/
264
265// CPU launchers
267 std::function< void(
268 const double * const __restrict__ Q,
270 double t,
271 int normal,
272 double * __restrict__ F
273 ) > flux,
274 std::function< void(
275 const double * const __restrict__ Q,
277 double t,
278 double * __restrict__ S
279 ) > algebraicSource,
280 std::function< void(
281 const double * const __restrict__ Q,
282 double * __restrict__ dQ_or_deltaQ,
284 double t,
285 int normal,
286 double * __restrict__ BgradQ
287 ) > nonconservativeProduct,
288 double * __restrict__ UOut,
289 const double * const __restrict__ QIn,
290 const double * const __restrict__ nodes,
291 const double * const __restrict__ weights,
292 const double * const __restrict__ Kxi,
293 const double * const __restrict__ dudx,
294 const tarch::la::Vector<Dimensions,double>& cellCentre,
295 const double dx,
296 const double t,
297 const double dt,
298 const int order,
299 const int unknowns,
300 const int auxiliaryVariables,
301 const bool callFlux,
302 const bool callSource,
303 const bool callNonconservativeProduct) {
304 const int nodesPerAxis = order+1;
305
306 const int nodesPerCell = getNodesPerCell(nodesPerAxis);
307
308 const int strideQ = unknowns+auxiliaryVariables;
309 const int strideRhs = unknowns;
310 const int strideS = unknowns;
311 const int strideF = unknowns;
312 const int strideGradQ = strideQ*Dimensions; // gradient of auxiliary variables needed for some apps
313
314 double* FAux = nullptr;
315 double* SAux = nullptr;
316 double* gradQAux = nullptr;
317 if ( callFlux ) {
318 FAux = new double[nodesPerCell*strideF]{0.0};
319 }
320 if ( callSource || callNonconservativeProduct ) {
321 SAux = new double[nodesPerCell*strideS]{0.0};
322 }
323 if ( callNonconservativeProduct ) {
324 gradQAux = new double[nodesPerCell*strideGradQ]{0.0};
325 }
326 for ( int scalarIndexCell = 0; scalarIndexCell < nodesPerCell; scalarIndexCell++ ) {
327 if ( callFlux ) {
329 flux,
330 UOut,
331 QIn,
332 FAux + scalarIndexCell*strideF,
333 nodes,
334 weights,
335 Kxi,
336 cellCentre,
337 dx,
338 t,
339 dt,
340 nodesPerAxis,
341 unknowns,
342 strideQ,
343 scalarIndexCell);
344 }
345 if ( callSource ) {
347 algebraicSource,
348 UOut,
349 SAux + scalarIndexCell*strideS,
350 QIn,
351 nodes,
352 weights,
353 cellCentre,
354 dx,
355 t,
356 dt,
357 nodesPerAxis,
358 unknowns,
359 strideQ,
360 scalarIndexCell);
361 }
362 if ( callNonconservativeProduct ) {
364 nonconservativeProduct,
365 UOut,
366 gradQAux + scalarIndexCell*strideGradQ,
367 SAux + scalarIndexCell*strideS,
368 QIn,
369 nodes,
370 weights,
371 dudx,
372 cellCentre,
373 dx,
374 t,
375 dt,
376 nodesPerAxis,
377 unknowns,
378 strideQ,
379 scalarIndexCell);
380 }
381 } // scalarIndexCell
382
383 if ( callFlux ) {
384 delete [] FAux;
385 }
386 if ( callSource || callNonconservativeProduct ) {
387 delete [] SAux;
388 }
389 if ( callNonconservativeProduct ) {
390 delete [] gradQAux;
391 }
392}
393
395 std::function< void(
396 double * __restrict__ Q,
398 double t
399 ) > adjustSolution,
400 std::function< double(
401 const double * const __restrict__ Q,
403 double t,
404 const int normal
405 ) > maxAbsoluteEigenvalue,
406 double * __restrict__ UOut,
407 const double * const __restrict__ riemannResultIn,
408 const double * const __restrict__ nodes,
409 const double * const __restrict__ weights,
410 const double * const __restrict__ FLCoeff,
411 const tarch::la::Vector<Dimensions,double>& cellCentre,
412 const double dx,
413 const double t,
414 const double dt,
415 const int order,
416 const int unknowns,
417 const int auxiliaryVariables,
418 const bool callMaxEigenvalue) {
419 const int nodesPerAxis = order+1;
420
421 const int nodesPerCell = getNodesPerCell(nodesPerAxis);
422
423 const int strideQ = unknowns+auxiliaryVariables;
424 const int strideF = unknowns;
425
426 double lambdaMax = 0.0;
427 for ( int scalarIndexCell = 0; scalarIndexCell < nodesPerCell; scalarIndexCell++ ) {
429 UOut,
430 riemannResultIn,
431 weights,
432 FLCoeff,
433 dx,
434 dt,
435 nodesPerAxis,
436 unknowns,
437 strideQ,
438 strideF,
439 scalarIndexCell);
440
441 const double result =
443 adjustSolution,
444 maxAbsoluteEigenvalue,
445 UOut,
446 nodes,
447 cellCentre,
448 dx,
449 t+dt,
450 nodesPerAxis,
451 unknowns,
452 strideQ,
453 callMaxEigenvalue,
454 scalarIndexCell);
455 lambdaMax = std::max( lambdaMax, result );
456 }
457 return lambdaMax;
458}
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$ our matrix elements will nabla phi_i dx f By this will be a sparse as these basis functions are chosen to not overlap with each other almost everywhere In other they have only local support We can read off the right hand side by taking our known right hand side f$ f f$ and integrating against an appropriate test phi_i dx f Please excuse the slight abuse of notation here There should probably be a clearer indication that we move from a continuous f$ f f$ to some discrete f$ f_i f$ We can demonstrate simply It s worth as when we discuss the discontinuous version of this it will no longer disappear We take our left hand side and discretise it
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
#define GPUCallableMethod
Definition accelerator.h:25
#define FLCoeff(i)
Definition KernelUtils.h:8
GPUCallableMethod tarch::la::Vector< Dimensions+1, int > delineariseIndex(int scalarIndex, const tarch::la::Vector< Dimensions+1, int > &strides)
GPUCallableMethod void gradient_AoS(const double *__restrict__ const QIn, const double *__restrict__ const dudx, const double invDx, const int nodesPerAxis, const int strideQ, const int scalarIndex, double *__restrict__ gradQAux)
Computes the gradient at a given (space-time) quadrature point.
double corrector_adjustSolution_computeMaxEigenvalue_body_AoS(std::function< void(double *__restrict__ Q, const tarch::la::Vector< Dimensions, double > &x, double t) > adjustSolution, std::function< double(const double *const __restrict__ Q, const tarch::la::Vector< Dimensions, double > &x, double t, const int normal) > maxAbsoluteEigenvalue, double *__restrict__ UOut, const double *const __restrict__ nodes, const tarch::la::Vector< Dimensions, double > &cellCentre, const double dx, const double t, const int nodesPerAxis, const int unknowns, const int strideQ, const int callMaxEigenvalue, const int scalarIndex)
Allow user to modify solution at the given coordinates and time and compute the max eigenvalue afterw...
void corrector_addCellContributions_loop_AoS(std::function< void(const double *const __restrict__ Q, const tarch::la::Vector< Dimensions, double > &x, double t, int normal, double *__restrict__ F) > flux, std::function< void(const double *const __restrict__ Q, const tarch::la::Vector< Dimensions, double > &x, double t, double *__restrict__ S) > algebraicSource, std::function< void(const double *const __restrict__ Q, double *__restrict__ dQ_or_deltaQ, const tarch::la::Vector< Dimensions, double > &x, double t, int normal, double *__restrict__ BgradQ) > nonconservativeProduct, double *__restrict__ UOut, const double *const __restrict__ QIn, const double *const __restrict__ nodes, const double *const __restrict__ weights, const double *const __restrict__ Kxi, const double *const __restrict__ dudx, 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 callFlux, const bool callSource, const bool callNonconservativeProduct)
Add cell-local contributions to new solution or update vector.
double corrector_addRiemannContributions_loop_AoS(std::function< void(double *__restrict__ Q, const tarch::la::Vector< Dimensions, double > &x, double t) > adjustSolution, std::function< double(const double *const __restrict__ Q, const tarch::la::Vector< Dimensions, double > &x, double t, const int normal) > maxAbsoluteEigenvalue, double *__restrict__ UOut, const double *const __restrict__ riemannResultIn, const double *const __restrict__ nodes, const double *const __restrict__ weights, const double *const __restrict__ FLCoeff, 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 callMaxEigenvalue)
Add cell-local contributions to new solution or update vector.
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:
GPUCallableMethod void corrector_addRiemannContributions_body_AoS(double *__restrict__ UOut, const double *const __restrict__ riemannResultIn, const double *const __restrict__ weights, const double *const __restrict__ FLCoeff, const double dx, const double dt, const int nodesPerAxis, const int unknowns, const int strideQ, const int strideF, const int scalarIndexCell)
Add Riemann flux contributions to the solution.
GPUCallableMethod int mapCellIndexToScalarHullIndex(const tarch::la::Vector< Dimensions+1, int > &indexCell, const int direction, const int orientation, const int nodesPerAxis)
Map cell index to an index for the whole hull, i.e.
void corrector_addSourceContributions_body_AoS(std::function< void(const double *const __restrict__ Q, const tarch::la::Vector< Dimensions, double > &x, double t, double *__restrict__ S) > algebraicSource, double *__restrict__ UOut, double *__restrict__ SAux, const double *__restrict__ QIn, const double *__restrict__ nodes, const double *__restrict__ weights, const tarch::la::Vector< Dimensions, double > &cellCentre, const double dx, const double t, const double dt, const int nodesPerAxis, const int unknowns, const int strideQ, const int scalarIndex)
Add source contributions to the solution.
void corrector_addFluxContributions_body_AoS(std::function< void(const double *const __restrict__ Q, const tarch::la::Vector< Dimensions, double > &x, double t, int normal, double *__restrict__ F) > flux, double *__restrict__ UOut, const double *__restrict__ QIn, double *__restrict__ FAux, const double *__restrict__ nodes, const double *__restrict__ weights, const double *__restrict__ Kxi, const tarch::la::Vector< Dimensions, double > &cellCentre, const double dx, const double t, const double dt, const int nodesPerAxis, const int unknowns, const int strideQ, const int scalarIndex)
Add space-time volume flux contributions to the solution.
GPUCallableMethod int getNodesPerCell(int nodesPerAxis)
GPUCallableMethod tarch::la::Vector< Dimensions+1, double > getCoordinates(const tarch::la::Vector< Dimensions+1, int > &index, const tarch::la::Vector< Dimensions, double > &centre, const double dx, const double t, const double dt, const double *__restrict__ nodes)
Compute coordinates from cell geometry and quadrature index.
void corrector_addNcpContributions_body_AoS(std::function< void(const double *const __restrict__ Q, double *__restrict__ dQ_or_deltaQ, const tarch::la::Vector< Dimensions, double > &x, double t, int normal, double *__restrict__ BgradQ) > nonconservativeProduct, double *__restrict__ UOut, double *__restrict__ gradQAux, double *__restrict__ SAux, const double *__restrict__ QIn, const double *__restrict__ nodes, const double *__restrict__ weights, const double *__restrict__ dudx, const tarch::la::Vector< Dimensions, double > &cellCentre, const double dx, const double t, const double dt, const int nodesPerAxis, const int unknowns, const int strideQ, const int scalarIndex)
Add nonconservative product contributions to the solution.
Simple vector class.
Definition Vector.h:134