Peano 4
Loading...
Searching...
No Matches
RusanovNonlinearAoS.cpp
Go to the documentation of this file.
2
3#include <iostream>
4
5#include "Generic.h"
6
7#include "KernelUtils.h"
8
9 // geht net durch NVC++ muessen wir angucken, sobald es geht
10/*
11#if defined(GPUOffloadingOMP)
12#pragma omp declare target
13#endif
14*/
16 std::function< void(
17 const double * const __restrict__ Q,
19 double t,
20 int normal,
21 double * __restrict__ F
22 ) > flux,
23 std::function< void(
24 double * __restrict__ Q,
25 double * __restrict__ dQ_or_deltaQ,
27 double t,
28 int normal,
29 double * __restrict__ BgradQ
30 ) > nonconservativeProduct,
31 double * __restrict__ riemannResultOut,
32 double * __restrict__ FLAux,
33 double * __restrict__ FRAux,
34 double * __restrict__ QAvgAux,
35 double * __restrict__ dQAux,
36 double * __restrict__ SAux,
37 const double * const __restrict__ QLIn,
38 const double * const __restrict__ QRIn,
39 const double smax,
40 const double * const __restrict__ nodes,
41 const double * const __restrict__ weights,
43 const double dx,
44 const double t,
45 const double dt,
46 const int nodesPerAxis,
47 const int unknowns,
48 const int strideQ,
49 const int strideF,
50 const int direction,
51 const bool orientationToCell,
52 const bool callFlux,
53 const bool callNonconservativeProduct,
54 const int scalarIndexHull) {
55 const tarch::la::Vector<Dimensions+1,int> strides = getStrides(nodesPerAxis,false);
56 const tarch::la::Vector<Dimensions+1,int> indexUHull = // (y,z,face) , (x,z,face), (x,y,face)
57 delineariseIndex(scalarIndexHull,strides);
59 getCoordinatesOnFace(indexUHull,faceCentre,direction,dx,t,dt,nodes);
60 const tarch::la::Vector<Dimensions,double> x(&coords[1]);
61
62 const int scalarIndexFace = scalarIndexHull - indexUHull[Dimensions]*strides[Dimensions];
63
64 //std::cout << " scalarIndexHull=" << scalarIndexHull << std::endl;
65 //std::cout << " scalarIndexFace=" << scalarIndexFace << std::endl;
66
67 // zero result vector
68 for (int var = 0; var < unknowns; var++) {
69 riemannResultOut[ scalarIndexHull*strideF + var ] = 0.0;
70 }
71 for (int it = 0; it < nodesPerAxis; it++) { // time integration
72 const double time = t + nodes[it] * dt;
73 const double* qL = &QLIn[ (scalarIndexFace*nodesPerAxis + it)*strideQ ];
74 const double* qR = &QRIn[ (scalarIndexFace*nodesPerAxis + it)*strideQ ];
75
76 //std::cout << " qL[0]=" << qL[0] << std::endl;
77 //std::cout << " qR[0]=" << qR[0] << std::endl;
78 //std::cout << " direction=" << direction << std::endl;
79
80 // time averaged flux contributions
81 if ( callFlux ) {
82 const double coeff1 = weights[it] * 0.5;
83 const double coeff2 = coeff1 * smax;
84
85 flux( qL, x, time, direction, FLAux );
86 flux( qR, x, time, direction, FRAux );
87 for (int var = 0; var < unknowns; var++) {
88 riemannResultOut[ scalarIndexHull*strideF + var ] += coeff1 * (FRAux[ var ] + FLAux[ var ]);
89 }
90 for (int var = 0; var < unknowns; var++) {
91 riemannResultOut[ scalarIndexHull*strideF + var ] -= coeff2 * (qR[ var ] - qL[ var ]);
92 }
93 }
94 // nonconservative product contirbutions
95 if ( callNonconservativeProduct ) {
96 // orientationToCell = 0 => face is left to cell => cell is right to face => cell gets FR => sign is negative
97 // orientationToCell = 1 => face is right to cell => cell is left to face => cell gets FL => sign is positive
98 const double coeff0 = (-1.0 + 2.0*orientationToCell) * 0.5;
99 for (int var=0; var < strideQ; var++) {
100 dQAux[ var ] = qR[ var ] - qL[ var ];
101 }
102 for (int var=0; var < strideQ; var++) {
103 QAvgAux[var] = 0.5 * (qR[ var ] + qL[ var ]);
104 }
105
106 nonconservativeProduct(QAvgAux, dQAux, x, time, direction, SAux);
107
108 const double coeff = coeff0*weights[it];
109 for (int var = 0; var < unknowns; var++) {
110 riemannResultOut[ scalarIndexHull*strideF + var ] += coeff * SAux[var];
111 }
112 }
113 }
114}
115/*
116#if defined(GPUOffloadingOMP)
117#pragma omp end declare target
118#endif
119*/
120
121// launchers
123 std::function< void(
124 const double * const __restrict__ Q,
126 double t,
127 int normal,
128 double * __restrict__ F
129 ) > flux,
130 std::function< void(
131 const double * const __restrict__ Q,
133 double t,
134 int normal,
135 double * __restrict__ F
136 ) > boundaryFlux,
137 std::function< void(
138 double * __restrict__ Q,
139 double * __restrict__ dQ_or_deltaQ,
141 double t,
142 int normal,
143 double * __restrict__ BgradQ
144 ) > nonconservativeProduct,
145 std::function< void(
146 double * __restrict__ Q,
147 double * __restrict__ dQ_or_deltaQ,
149 double t,
150 int normal,
151 double * __restrict__ BgradQ
152 ) > boundaryNonconservativeProduct,
153 double * __restrict__ riemannResultOut,
154 const double * const __restrict__ QHullIn[Dimensions*2],
155 const double maxEigenvaluePerFace[Dimensions*2],
156 const double * const __restrict__ nodes,
157 const double * const __restrict__ weights,
158 const tarch::la::Vector<Dimensions,double>& cellCentre,
159 const double dx,
160 const double t,
161 const double dt,
162 const int order,
163 const int unknowns,
164 const int auxiliaryVariables,
165 const bool atBoundary[Dimensions*2],
166 const bool callFlux,
167 const bool callNonconservativeProduct) {
168 const int nodesPerAxis = order + 1;
169
170 const int strideQ = unknowns+auxiliaryVariables;
171 const int strideS = unknowns;
172 const int strideF = unknowns;
173 const int strideDQ = strideQ; // gradient of auxiliary variables needed for some apps
174
175 // 384 work items for order 7 in 3D
176 const int nodesOnFace = getNodesPerCell(nodesPerAxis)/nodesPerAxis;
177 const int nodesOnHull = Dimensions*2 * nodesOnFace; // 2*d*nodesPerAxis^(d-1)
178
179 const int strideQLR = nodesOnFace*nodesPerAxis;
180
181 double * FAux = nullptr;
182 double * QAvgAux = nullptr;
183 double * SAux = nullptr;
184 double * dQAux = nullptr;
185 if ( callFlux ) {
186 FAux = new double[2*nodesOnHull*strideF]{0.0};
187 }
188 if ( callNonconservativeProduct ) {
189 QAvgAux = new double[nodesOnHull*strideQ]{0.0};
190 SAux = new double[nodesOnHull*strideS]{0.0};
191 dQAux = new double[nodesOnHull*strideDQ]{0.0};
192 }
193
194 for ( int scalarIndexHull = 0; scalarIndexHull < nodesOnHull; scalarIndexHull++ ) {
195 const int face = scalarIndexHull / nodesOnFace;
196 const int orientationToCell = face/Dimensions;
197 const int direction = face-Dimensions*orientationToCell;
198 //std::cout << "face=" << face << std::endl;
199 //std::cout << " orientationToCell=" << orientationToCell << std::endl;
200 //std::cout << " direction=" << direction << std::endl;
201 //std::cout << " strideQLR=" << strideQLR << std::endl;
202 const int scalarIndexFace = scalarIndexHull - face*nodesOnFace;
203
204 const bool leftCellIsOutside = atBoundary[ face ] && orientationToCell == 0;
205 const bool rightCellIsOutside = atBoundary[ face ] && orientationToCell == 1;
206
207 tarch::la::Vector<Dimensions,double> faceCentre = cellCentre;
208 faceCentre[ direction ] += (-1+2*orientationToCell)*0.5*dx;
209
211 ( leftCellIsOutside || rightCellIsOutside ) ? boundaryFlux : flux,
212 ( leftCellIsOutside || rightCellIsOutside ) ? boundaryNonconservativeProduct : nonconservativeProduct,
213 riemannResultOut, // we index the whole hull
214 FAux + scalarIndexHull*2*strideF,
215 FAux + scalarIndexHull*2*strideF + strideF,
216 QAvgAux + scalarIndexHull*strideQ,
217 dQAux + scalarIndexHull*strideDQ,
218 SAux + scalarIndexHull*strideS,
219 QHullIn[ face ],
220 QHullIn[ face ] + strideQLR*strideQ,
221 maxEigenvaluePerFace[ face ],
222 nodes,
223 weights,
224 faceCentre,
225 dx,
226 t,
227 dt,
228 nodesPerAxis,
229 unknowns,
230 strideQ,
231 strideF,
232 direction,
233 orientationToCell,
234 callFlux,
235 callNonconservativeProduct,
236 scalarIndexHull);
237 }
238
239 if ( callFlux ) {
240 delete [] FAux;
241 }
242 if ( callNonconservativeProduct ) {
243 delete [] QAvgAux;
244 delete [] SAux;
245 delete [] dQAux;
246 }
247}
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
void rusanovNonlinear_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, int normal, double *__restrict__ F) > boundaryFlux, std::function< void(double *__restrict__ Q, double *__restrict__ dQ_or_deltaQ, const tarch::la::Vector< Dimensions, double > &x, double t, int normal, double *__restrict__ BgradQ) > nonconservativeProduct, std::function< void(double *__restrict__ Q, double *__restrict__ dQ_or_deltaQ, const tarch::la::Vector< Dimensions, double > &x, double t, int normal, double *__restrict__ BgradQ) > boundaryNonconservativeProduct, double *__restrict__ riemannResultOut, const double *const __restrict__ QHullIn[Dimensions *2], const double maxEigenvaluePerFace[Dimensions *2], const double *const __restrict__ nodes, const double *const __restrict__ weights, 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 atBoundary[Dimensions *2], const bool callFlux, const bool callNonconservativeProduct)
Apply a Rusanov Riemann solver to the pair of extrapolated predictor values from two neighbouring cel...
GPUCallableMethod tarch::la::Vector< Dimensions+1, int > delineariseIndex(int scalarIndex, const tarch::la::Vector< Dimensions+1, int > &strides)
void rusanovNonlinear_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, std::function< void(double *__restrict__ Q, double *__restrict__ dQ_or_deltaQ, const tarch::la::Vector< Dimensions, double > &x, double t, int normal, double *__restrict__ BgradQ) > nonconservativeProduct, double *__restrict__ riemannResultOut, double *__restrict__ FLAux, double *__restrict__ FRAux, double *__restrict__ QAvgAux, double *__restrict__ dQAux, double *__restrict__ SAux, const double *const __restrict__ QLIn, const double *const __restrict__ QRIn, const double smax, const double *const __restrict__ nodes, const double *const __restrict__ weights, const tarch::la::Vector< Dimensions, double > &faceCentre, const double dx, const double t, const double dt, const int nodesPerAxis, const int unknowns, const int strideQ, const int strideF, const int direction, const bool orientationToCell, const bool callFlux, const bool callNonconservativeProduct, const int scalarIndexFace)
GPUCallableMethod tarch::la::Vector< Dimensions+1, double > getCoordinatesOnFace(const tarch::la::Vector< Dimensions+1, int > &indexOnFace, const tarch::la::Vector< Dimensions, double > &faceCentre, const int direction, const double dx, const double t, const double dt, const double *__restrict__ nodes)
Compute coordinates from cell geometry and quadrature index.
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 int getNodesPerCell(int nodesPerAxis)
Simple vector class.
Definition Vector.h:134