Vlasiator ebf0dd394 on dev (v5.4.0 + 1054 commits)
Loading...
Searching...
No Matches
gpu_moments.cpp
Go to the documentation of this file.
1/*
2 * This file is part of Vlasiator.
3 * Copyright 2010-2024 Finnish Meteorological Institute and University of Helsinki
4 *
5 * For details of usage, see the COPYING file and read the "Rules of the Road"
6 * at http://www.physics.helsinki.fi/vlasiator/
7 *
8 * This program is free software; you can redistribute it and/or modify
9 * it under the terms of the GNU General Public License as published by
10 * the Free Software Foundation; either version 2 of the License, or
11 * (at your option) any later version.
12 *
13 * This program is distributed in the hope that it will be useful,
14 * but WITHOUT ANY WARRANTY; without even the implied warranty of
15 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
16 * GNU General Public License for more details.
17 *
18 * You should have received a copy of the GNU General Public License along
19 * with this program; if not, write to the Free Software Foundation, Inc.,
20 * 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
21 */
22
23#include <phiprof.hpp>
24#include "gpu_moments.h"
25#include "vlasovmover.h"
26#include "../object_wrapper.h"
27#include "../fieldsolver/fs_common.h" // divideIfNonZero(
28#include "../arch/gpu_base.hpp"
29
30using namespace std;
31
36__global__ void __launch_bounds__(WID3) first_moments_kernel (
37 const vmesh::VelocityBlockContainer* __restrict__ const *dev_VBC,
38 Real* dev_moments1,
39 const uint nAllCells)
40{
41 const uint celli = blockIdx.x; // used for pointer to cell
42 const uint stride = gridDim.y; // used for faster looping over contents
43 const uint strideOffset = blockIdx.y; // used for faster looping over contents
44 const uint ti = threadIdx.z*blockDim.x*blockDim.y + threadIdx.y*blockDim.x + threadIdx.x;
45 const int i = threadIdx.x;
46 const int j = threadIdx.y;
47 const int k = threadIdx.z;
48 const int blockSize = blockDim.x*blockDim.y*blockDim.z;
49
50 extern __shared__ Real smom[];
51 Real myMom[nMom1] = {0};
52
53 const vmesh::VelocityBlockContainer* __restrict__ blockContainer = dev_VBC[celli];
54 if (blockContainer==0) {
55 return;
56 }
57 const uint thisVBCSize = blockContainer->size();
58 if (thisVBCSize==0) {
59 return;
60 }
61 const Realf* __restrict__ data = blockContainer->getData();
62 const Real* __restrict__ blockParameters = blockContainer->getParameters();
63 const Real HALF = 0.5;
64 const int indx = ti%WID3;
65 for (uint blockIndex = strideOffset; blockIndex < thisVBCSize; blockIndex += stride) {
66 const Realf cellValue = data[blockIndex*WID3+indx];
67 const Real* __restrict__ blockParamsZ = &blockParameters[blockIndex*BlockParams::N_VELOCITY_BLOCK_PARAMS];
68
69 const Real paramX = blockParamsZ[BlockParams::DVX];
70 const Real paramY = blockParamsZ[BlockParams::DVY];
71 const Real paramZ = blockParamsZ[BlockParams::DVZ];
72
73 const Real DV3 = paramX*paramY*paramZ;
74
75 const Real VX = blockParamsZ[BlockParams::VXCRD] + (i+HALF)*paramX;
76 const Real VY = blockParamsZ[BlockParams::VYCRD] + (j+HALF)*paramY;
77 const Real VZ = blockParamsZ[BlockParams::VZCRD] + (k+HALF)*paramZ;
78
79 const Real f = (Real)cellValue;
80 const Real fDV3 = f * DV3;
81
82 myMom[0] += fDV3;
83 myMom[1] += fDV3 * VX;
84 myMom[2] += fDV3 * VY;
85 myMom[3] += fDV3 * VZ;
86 }
87
88 const int indexInsideWarp = ti % GPUTHREADS;
89 const int warpIndex = ti / GPUTHREADS;
90
91 //Now reduce one-by-one for cell
92 for (int offset = GPUTHREADS/2; offset > 0; offset /= 2) {
93 for (uint imom=0; imom<nMom1; imom++) {
94 myMom[imom] += gpuKernelShflDown(myMom[imom], offset);
95 }
96 }
97
98 if (indexInsideWarp == 0) {
99 for (uint imom=0; imom<nMom1; imom++) {
100 smom[warpIndex*nMom1+imom] = myMom[imom];
101 }
102 }
103
104 __syncthreads();
105
106 const int warpsPerBlock = blockSize/GPUTHREADS;
107
108 for (uint imom=warpIndex; imom<nMom1; imom +=warpsPerBlock ) {
109 myMom[imom] = (indexInsideWarp < warpsPerBlock) ? smom[indexInsideWarp*nMom1+imom] : 0.0;
110 for (int offset = (warpsPerBlock)/2; offset > 0; offset /= 2) {
111 myMom[imom] += gpuKernelShflDown(myMom[imom], offset);
112 }
113
114 if (indexInsideWarp == 0) {
115 if (gridDim.y == 1) {
116 dev_moments1[celli*nMom1 + imom] = myMom[imom];
117 } else {
118 atomicAdd(&dev_moments1[celli*nMom1 + imom],myMom[imom]);
119 }
120 }
121 }
122 }
123
128__global__ void __launch_bounds__(WID3) second_moments_kernel (
129 const vmesh::VelocityBlockContainer* __restrict__ const *dev_VBC,
130 const Real* __restrict__ dev_moments1,
131 Real* dev_moments2,
132 const uint nAllCells)
133{
134 const uint celli = blockIdx.x; // used for pointer to cell
135 const uint stride = gridDim.y; // used for faster looping over contents
136 const uint strideOffset = blockIdx.y; // used for faster looping over contents
137 const uint ti = threadIdx.z*blockDim.x*blockDim.y + threadIdx.y*blockDim.x + threadIdx.x;
138 const int i = threadIdx.x;
139 const int j = threadIdx.y;
140 const int k = threadIdx.z;
141 const int blockSize = blockDim.x*blockDim.y*blockDim.z;
142
143 Real myMom[nMom2] = {0};
144 extern __shared__ Real smom[];
145
146 const vmesh::VelocityBlockContainer* __restrict__ blockContainer = dev_VBC[celli];
147 if (blockContainer==0) {
148 return;
149 }
150 const uint thisVBCSize = blockContainer->size();
151 if (thisVBCSize==0) {
152 return;
153 }
154 const Realf* __restrict__ data = blockContainer->getData();
155 const Real* __restrict__ blockParameters = blockContainer->getParameters();
156
157 const Real HALF = 0.5;
158 const int indx = ti%WID3;
159
160 // index +0 is number density
161 const Real averageVX = dev_moments1[celli*nMom1 + 1];
162 const Real averageVY = dev_moments1[celli*nMom1 + 2];
163 const Real averageVZ = dev_moments1[celli*nMom1 + 3];
164
165 for (uint blockIndex = strideOffset; blockIndex < thisVBCSize; blockIndex += stride) {
166 const Realf cellValue = data[blockIndex*WID3+indx];
167 const Real* __restrict__ blockParamsZ = &blockParameters[blockIndex*BlockParams::N_VELOCITY_BLOCK_PARAMS];
168
169 const Real paramX = blockParamsZ[BlockParams::DVX];
170 const Real paramY = blockParamsZ[BlockParams::DVY];
171 const Real paramZ = blockParamsZ[BlockParams::DVZ];
172
173 const Real DV3 = paramX*paramY*paramZ;
174
175 const Real VX = blockParamsZ[BlockParams::VXCRD] + (i+HALF)*paramX;
176 const Real VY = blockParamsZ[BlockParams::VYCRD] + (j+HALF)*paramY;
177 const Real VZ = blockParamsZ[BlockParams::VZCRD] + (k+HALF)*paramZ;
178
179 const Real VXDifference = VX - averageVX;
180 const Real VYDifference = VY - averageVY;
181 const Real VZDifference = VZ - averageVZ;
182
183 const Real f = (Real)cellValue;
184 const Real fDV3 = f * DV3;
185
186 const Real fDV3VXDifference = fDV3 * VXDifference;
187 const Real fDV3VYDifference = fDV3 * VYDifference;
188
189 myMom[0] += fDV3VXDifference * VXDifference;
190 myMom[1] += fDV3VYDifference * VYDifference;
191 myMom[2] += fDV3 * VZDifference * VZDifference;
192 myMom[3] += fDV3VYDifference * VZDifference;
193 myMom[4] += fDV3VXDifference * VZDifference;
194 myMom[5] += fDV3VXDifference * VYDifference;
195 }
196
197 const int indexInsideWarp = ti % GPUTHREADS;
198 const int warpIndex = ti / GPUTHREADS;
199
200 //Now reduce one-by-one for cell
201 for (int offset = GPUTHREADS/2; offset > 0; offset /= 2) {
202 for (uint imom=0; imom<nMom2; imom++) {
203 myMom[imom] += gpuKernelShflDown(myMom[imom], offset);
204 }
205 }
206
207 if (indexInsideWarp == 0) {
208 for (uint imom=0; imom<nMom2; imom++) {
209 smom[warpIndex*nMom2+imom] = myMom[imom];
210 }
211 }
212
213 __syncthreads();
214
215 const int warpsPerBlock = blockSize/GPUTHREADS;
216
217 for (uint imom=warpIndex; imom<nMom2; imom +=warpsPerBlock ) {
218 myMom[imom] = (indexInsideWarp < warpsPerBlock) ? smom[indexInsideWarp*nMom2+imom] : 0.0;
219 for (int offset = (warpsPerBlock)/2; offset > 0; offset /= 2) {
220 myMom[imom] += gpuKernelShflDown(myMom[imom], offset);
221 }
222
223 if (indexInsideWarp == 0) {
224 if (gridDim.y == 1) {
225 dev_moments2[celli*nMom2 + imom] = myMom[imom];
226 } else {
227 atomicAdd(&dev_moments2[celli*nMom2 + imom],myMom[imom]);
228 }
229 }
230 }
231}
232
238
249void gpu_calculateMoments_R(
250 dccrg::Dccrg<SpatialCell,dccrg::Cartesian_Geometry>& mpiGrid,
251 const std::vector<CellID>& cells_in,
252 const bool computeSecond,
253 const bool initialCompute) {
254
255 phiprof::Timer computeMomentsTimer {"Compute _R moments"};
256
257 // Ensure unique cells
258 std::vector<CellID> cells = cells_in;
259 std::sort( cells.begin(), cells.end() );
260 cells.erase( std::unique( cells.begin(), cells.end() ), cells.end() );
261
262 const uint nAllCells = cells.size();
263 if (nAllCells==0) {
264 return;
265 }
266
267 gpuMemoryManager.startSession(0,0);
268
269 SESSION_HOST_ALLOCATE(gpuMemoryManager, vmesh::VelocityBlockContainer*, host_VBC, nAllCells*sizeof(vmesh::VelocityBlockContainer*));
270 SESSION_HOST_ALLOCATE(gpuMemoryManager, Real, host_moments1, nAllCells*nMom1*sizeof(Real));
271 SESSION_HOST_ALLOCATE(gpuMemoryManager, Real, host_moments2, nAllCells*nMom2*sizeof(Real));
272 SESSION_ALLOCATE(gpuMemoryManager, vmesh::VelocityBlockContainer*, dev_VBC, nAllCells*sizeof(vmesh::VelocityBlockContainer*));
273 SESSION_ALLOCATE(gpuMemoryManager, Real, dev_moments1, nAllCells*nMom1*sizeof(Real));
274 SESSION_ALLOCATE(gpuMemoryManager, Real, dev_moments2, nAllCells*nMom2*sizeof(Real));
275
276 vmesh::VelocityBlockContainer** host_VBC = GET_SESSION_HOST_POINTER(gpuMemoryManager, vmesh::VelocityBlockContainer*, host_VBC);
277 Real* host_moments1 = GET_SESSION_HOST_POINTER(gpuMemoryManager, Real, host_moments1);
278 Real* host_moments2 = GET_SESSION_HOST_POINTER(gpuMemoryManager, Real, host_moments2);
279 vmesh::VelocityBlockContainer** dev_VBC = GET_SESSION_POINTER(gpuMemoryManager, vmesh::VelocityBlockContainer*, dev_VBC);
280 Real* dev_moments1 = GET_SESSION_POINTER(gpuMemoryManager, Real, dev_moments1);
281 Real* dev_moments2 = GET_SESSION_POINTER(gpuMemoryManager, Real, dev_moments2);
282
283 std::vector<vmesh::LocalID> maxVmeshSizes;
284
285 for (uint popID=0; popID<getObjectWrapper().particleSpecies.size(); ++popID) {
286 // Gather VBCs
287 maxVmeshSizes.push_back(0);
288 #pragma omp parallel
289 {
290 vmesh::LocalID threadMaxVmeshSize = 0;
291 #pragma omp for schedule(static)
292 for(uint celli = 0; celli < nAllCells; celli++){
293 SpatialCell* cell = mpiGrid[cells[celli]];
294 if (cell->sysBoundaryFlag == sysboundarytype::DO_NOT_COMPUTE) {
295 host_VBC[celli] = 0;
296 continue;
297 }
298 if (cell->sysBoundaryFlag == sysboundarytype::OUTFLOW && cell->sysBoundaryLayer != 1 && !initialCompute) {
299 // these should have been handled by the boundary code
300 host_VBC[celli] = 0;
301 continue;
302 }
303 host_VBC[celli] = cell->dev_get_velocity_blocks(popID); // GPU-side VBC
304 // Evaluate cached vmesh size
305 const vmesh::LocalID meshSize = cell->get_velocity_mesh(popID)->size();
306 threadMaxVmeshSize = meshSize > threadMaxVmeshSize ? meshSize : threadMaxVmeshSize;
307
308 // Clear old moments to zero value
309 if (popID == 0) {
310 cell->parameters[CellParams::RHOM_R ] = 0.0;
311 cell->parameters[CellParams::VX_R] = 0.0;
312 cell->parameters[CellParams::VY_R] = 0.0;
313 cell->parameters[CellParams::VZ_R] = 0.0;
314 cell->parameters[CellParams::RHOQ_R ] = 0.0;
315 cell->parameters[CellParams::P_11_R] = 0.0;
316 cell->parameters[CellParams::P_22_R] = 0.0;
317 cell->parameters[CellParams::P_33_R] = 0.0;
318 cell->parameters[CellParams::P_23_R] = 0.0;
319 cell->parameters[CellParams::P_13_R] = 0.0;
320 cell->parameters[CellParams::P_12_R] = 0.0;
321 }
322 }
323 #pragma omp critical
324 {
325 maxVmeshSizes.at(popID) = maxVmeshSizes.at(popID) > threadMaxVmeshSize ? maxVmeshSizes.at(popID) : threadMaxVmeshSize;
326 }
327 }
328 if (maxVmeshSizes.at(popID) == 0) {
329 continue;
330 }
331
332 vmesh::LocalID maxVmeshLaunch = (gpuMultiProcessorCount*min(threadsPerMP/WID3, blocksPerMP))/nAllCells;
333 maxVmeshLaunch = maxVmeshLaunch < 1 ? 1 : maxVmeshLaunch;
334 // Send pointers, set initial data to zero
335 CHK_ERR( gpuMemcpy(dev_VBC, host_VBC, nAllCells*sizeof(vmesh::VelocityBlockContainer*), gpuMemcpyHostToDevice) );
336 CHK_ERR( gpuMemset(dev_moments1, 0, nAllCells*nMom1*sizeof(Real)) );
337 // Launch kernel calculating this species' contribution to first velocity moments
338 dim3 blockSize(WID,WID,WID);
339 dim3 gridSize(nAllCells,maxVmeshLaunch,1);
340 int sharedMemorySizeFirstMoments = nMom1 * (WID3 / GPUTHREADS) * sizeof(Real);
341 first_moments_kernel<<<gridSize, blockSize, sharedMemorySizeFirstMoments, 0>>> (
342 dev_VBC,
343 dev_moments1,
344 nAllCells
345 );
346 CHK_ERR( gpuPeekAtLastError() );
347 CHK_ERR( gpuDeviceSynchronize());
348 CHK_ERR( gpuMemcpy(host_moments1, dev_moments1, nAllCells*nMom1*sizeof(Real), gpuMemcpyDeviceToHost) );
349
350 const Real mass = getObjectWrapper().particleSpecies[popID].mass;
351 const Real charge = getObjectWrapper().particleSpecies[popID].charge;
352 #pragma omp parallel for schedule(static)
353 for (uint celli = 0; celli < nAllCells; celli++){
354 SpatialCell* cell = mpiGrid[cells[celli]];
355 if (cell->sysBoundaryFlag == sysboundarytype::DO_NOT_COMPUTE) {
356 continue;
357 }
358 if (cell->sysBoundaryFlag == sysboundarytype::OUTFLOW && cell->sysBoundaryLayer != 1 && !initialCompute) {
359 // these should have been handled by the boundary code
360 continue;
361 }
362
363 // Store species' contribution to bulk velocity moments
364 Population &pop = cell->get_population(popID);
365 pop.RHO_R = host_moments1[nMom1*celli];
366 pop.V_R[0] = divideIfNonZero(host_moments1[nMom1*celli + 1], host_moments1[nMom1*celli]);
367 pop.V_R[1] = divideIfNonZero(host_moments1[nMom1*celli + 2], host_moments1[nMom1*celli]);
368 pop.V_R[2] = divideIfNonZero(host_moments1[nMom1*celli + 3], host_moments1[nMom1*celli]);
369
370 cell->parameters[CellParams::RHOM_R ] += host_moments1[nMom1*celli]*mass;
371 cell->parameters[CellParams::VX_R] += host_moments1[nMom1*celli + 1]*mass;
372 cell->parameters[CellParams::VY_R] += host_moments1[nMom1*celli + 2]*mass;
373 cell->parameters[CellParams::VZ_R] += host_moments1[nMom1*celli + 3]*mass;
374 cell->parameters[CellParams::RHOQ_R ] += host_moments1[nMom1*celli]*charge;
375 } // for-loop over spatial cells
376 } // for-loop over particle species
377
378 #pragma omp parallel for schedule(static)
379 for (size_t celli=0; celli<nAllCells; ++celli) {
380 SpatialCell* cell = mpiGrid[cells[celli]];
381 if (cell->sysBoundaryFlag == sysboundarytype::DO_NOT_COMPUTE) {
382 continue;
383 }
384 if (cell->sysBoundaryFlag == sysboundarytype::OUTFLOW && cell->sysBoundaryLayer != 1 && !initialCompute) {
385 // these should have been handled by the boundary code
386 continue;
387 }
388 cell->parameters[CellParams::VX_R] = divideIfNonZero(cell->parameters[CellParams::VX_R], cell->parameters[CellParams::RHOM_R]);
389 cell->parameters[CellParams::VY_R] = divideIfNonZero(cell->parameters[CellParams::VY_R], cell->parameters[CellParams::RHOM_R]);
390 cell->parameters[CellParams::VZ_R] = divideIfNonZero(cell->parameters[CellParams::VZ_R], cell->parameters[CellParams::RHOM_R]);
391
392 // copy the bulk flow frame back to device
393 //host_moments1[nMom1*celli + 0] = cell->parameters[CellParams::RHOM_R];
394 host_moments1[nMom1*celli + 1] = cell->parameters[CellParams::VX_R];
395 host_moments1[nMom1*celli + 2] = cell->parameters[CellParams::VY_R];
396 host_moments1[nMom1*celli + 3] = cell->parameters[CellParams::VZ_R];
397 }
398
399 // Compute second moments only if requested.
400 if (computeSecond == false) {
401 gpuMemoryManager.endSession();
402 return;
403 }
404
405 CHK_ERR( gpuMemcpy(dev_moments1, host_moments1, nAllCells*nMom1*sizeof(Real), gpuMemcpyHostToDevice) );
406 CHK_ERR( gpuMemset(dev_moments2, 0, nAllCells*nMom2*sizeof(Real)) );
407 for (uint popID=0; popID<getObjectWrapper().particleSpecies.size(); ++popID) {
408 if (maxVmeshSizes.at(popID) == 0) {
409 continue;
410 }
411 // Launch kernel calculating this species' contribution to second velocity moments
412
413 vmesh::LocalID maxVmeshLaunch = (gpuMultiProcessorCount*min(threadsPerMP/WID3, blocksPerMP))/nAllCells;
414 maxVmeshLaunch = maxVmeshLaunch < 1 ? 1 : maxVmeshLaunch;
415
416 dim3 blockSize(WID,WID,WID);
417 dim3 gridSize(nAllCells,maxVmeshLaunch,1);
418 int sharedMemorySizeSecondMoments = nMom2 * (WID3 / GPUTHREADS) * sizeof(Real);
419 second_moments_kernel<<<gridSize, blockSize, sharedMemorySizeSecondMoments, 0>>> (
420 dev_VBC,
421 dev_moments1,
422 dev_moments2,
423 nAllCells
424 );
425 CHK_ERR( gpuPeekAtLastError() );
426 CHK_ERR( gpuDeviceSynchronize());
427 CHK_ERR( gpuMemcpy(host_moments2, dev_moments2, nAllCells*nMom2*sizeof(Real), gpuMemcpyDeviceToHost) );
428
429 const Real mass = getObjectWrapper().particleSpecies[popID].mass;
430 #pragma omp parallel for schedule(static)
431 for (uint celli = 0; celli < nAllCells; celli++){
432 SpatialCell* cell = mpiGrid[cells[celli]];
433 if (cell->sysBoundaryFlag == sysboundarytype::DO_NOT_COMPUTE) {
434 continue;
435 }
436 if (cell->sysBoundaryFlag == sysboundarytype::OUTFLOW && cell->sysBoundaryLayer != 1 && !initialCompute) {
437 // these should have been handled by the boundary code
438 continue;
439 }
440
441 // Store species' contribution to bulk velocity moments
442 Population &pop = cell->get_population(popID);
443 for (size_t i = 0; i < nMom2; ++i) {
444 pop.P_R[i] = mass*host_moments2[nMom2*celli + i];
445 }
446
447 cell->parameters[CellParams::P_11_R] += pop.P_R[0];
448 cell->parameters[CellParams::P_22_R] += pop.P_R[1];
449 cell->parameters[CellParams::P_33_R] += pop.P_R[2];
450 cell->parameters[CellParams::P_23_R] += pop.P_R[3];
451 cell->parameters[CellParams::P_13_R] += pop.P_R[4];
452 cell->parameters[CellParams::P_12_R] += pop.P_R[5];
453 } // for-loop over spatial cells
454 } // for-loop over particle species
455
456 gpuMemoryManager.endSession();
457}
458
469void gpu_calculateMoments_V(
470 dccrg::Dccrg<SpatialCell,dccrg::Cartesian_Geometry>& mpiGrid,
471 const std::vector<CellID>& cells_in,
472 const bool computeSecond,
473 const bool initialCompute) {
474
475 phiprof::Timer computeMomentsTimer {"Compute _V moments"};
476
477 // Ensure unique cells
478 std::vector<CellID> cells = cells_in;
479 std::sort( cells.begin(), cells.end() );
480 cells.erase( std::unique( cells.begin(), cells.end() ), cells.end() );
481
482 const uint nAllCells = cells.size();
483 if (nAllCells==0) {
484 return;
485 }
486
487 gpuMemoryManager.startSession(0,0);
488
489 SESSION_HOST_ALLOCATE(gpuMemoryManager, vmesh::VelocityBlockContainer*, host_VBC, nAllCells*sizeof(vmesh::VelocityBlockContainer*));
490 SESSION_HOST_ALLOCATE(gpuMemoryManager, Real, host_moments1, nAllCells*nMom1*sizeof(Real));
491 SESSION_HOST_ALLOCATE(gpuMemoryManager, Real, host_moments2, nAllCells*nMom2*sizeof(Real));
492 SESSION_ALLOCATE(gpuMemoryManager, vmesh::VelocityBlockContainer*, dev_VBC, nAllCells*sizeof(vmesh::VelocityBlockContainer*));
493 SESSION_ALLOCATE(gpuMemoryManager, Real, dev_moments1, nAllCells*nMom1*sizeof(Real));
494 SESSION_ALLOCATE(gpuMemoryManager, Real, dev_moments2, nAllCells*nMom2*sizeof(Real));
495
496 vmesh::VelocityBlockContainer** host_VBC = GET_SESSION_HOST_POINTER(gpuMemoryManager, vmesh::VelocityBlockContainer*, host_VBC);
497 Real* host_moments1 = GET_SESSION_HOST_POINTER(gpuMemoryManager, Real, host_moments1);
498 Real* host_moments2 = GET_SESSION_HOST_POINTER(gpuMemoryManager, Real, host_moments2);
499 vmesh::VelocityBlockContainer** dev_VBC = GET_SESSION_POINTER(gpuMemoryManager, vmesh::VelocityBlockContainer*, dev_VBC);
500 Real* dev_moments1 = GET_SESSION_POINTER(gpuMemoryManager, Real, dev_moments1);
501 Real* dev_moments2 = GET_SESSION_POINTER(gpuMemoryManager, Real, dev_moments2);
502
503 std::vector<vmesh::LocalID> maxVmeshSizes;
504
505 for (uint popID=0; popID<getObjectWrapper().particleSpecies.size(); ++popID) {
506 // Gather VBCs
507 maxVmeshSizes.push_back(0);
508 #pragma omp parallel
509 {
510 vmesh::LocalID threadMaxVmeshSize = 0;
511 #pragma omp for schedule(static)
512 for(uint celli = 0; celli < nAllCells; celli++){
513 SpatialCell* cell = mpiGrid[cells[celli]];
514 if (cell->sysBoundaryFlag == sysboundarytype::DO_NOT_COMPUTE) {
515 host_VBC[celli] = 0;
516 continue;
517 }
518 if (cell->sysBoundaryFlag == sysboundarytype::OUTFLOW && cell->sysBoundaryLayer != 1 && !initialCompute) {
519 // these should have been handled by the boundary code
520 host_VBC[celli] = 0;
521 continue;
522 }
523 host_VBC[celli] = cell->dev_get_velocity_blocks(popID); // GPU-side VBC
524 // Evaluate cached vmesh size
525 const vmesh::LocalID meshSize = cell->get_velocity_mesh(popID)->size();
526 threadMaxVmeshSize = meshSize > threadMaxVmeshSize ? meshSize : threadMaxVmeshSize;
527
528 // Clear old moments to zero value
529 if (popID == 0) {
530 cell->parameters[CellParams::RHOM_V ] = 0.0;
531 cell->parameters[CellParams::VX_V] = 0.0;
532 cell->parameters[CellParams::VY_V] = 0.0;
533 cell->parameters[CellParams::VZ_V] = 0.0;
534 cell->parameters[CellParams::RHOQ_V ] = 0.0;
535 cell->parameters[CellParams::P_11_V] = 0.0;
536 cell->parameters[CellParams::P_22_V] = 0.0;
537 cell->parameters[CellParams::P_33_V] = 0.0;
538 cell->parameters[CellParams::P_23_V] = 0.0;
539 cell->parameters[CellParams::P_13_V] = 0.0;
540 cell->parameters[CellParams::P_12_V] = 0.0;
541 }
542 }
543 #pragma omp critical
544 {
545 maxVmeshSizes.at(popID) = maxVmeshSizes.at(popID) > threadMaxVmeshSize ? maxVmeshSizes.at(popID) : threadMaxVmeshSize;
546 }
547 }
548 if (maxVmeshSizes.at(popID) == 0) {
549 continue;
550 }
551
552 vmesh::LocalID maxVmeshLaunch = (gpuMultiProcessorCount*min(threadsPerMP/WID3, blocksPerMP))/nAllCells;
553 maxVmeshLaunch = maxVmeshLaunch < 1 ? 1 : maxVmeshLaunch;
554 // Send pointers, set initial data to zero
555 CHK_ERR( gpuMemcpy(dev_VBC, host_VBC, nAllCells*sizeof(vmesh::VelocityBlockContainer*), gpuMemcpyHostToDevice) );
556 CHK_ERR( gpuMemset(dev_moments1, 0, nAllCells*nMom1*sizeof(Real)) );
557 // Launch kernel calculating this species' contribution to first velocity moments
558 dim3 blockSize(WID,WID,WID);
559 dim3 gridSize(nAllCells,maxVmeshLaunch,1);
560 int sharedMemorySizeFirstMoments = nMom1 * (WID3 / GPUTHREADS) * sizeof(Real);
561 first_moments_kernel<<<gridSize, blockSize, sharedMemorySizeFirstMoments, 0>>> (
562 dev_VBC,
563 dev_moments1,
564 nAllCells
565 );
566 CHK_ERR( gpuPeekAtLastError() );
567 CHK_ERR( gpuDeviceSynchronize());
568 CHK_ERR( gpuMemcpy(host_moments1, dev_moments1, nAllCells*nMom1*sizeof(Real), gpuMemcpyDeviceToHost) );
569
570 const Real mass = getObjectWrapper().particleSpecies[popID].mass;
571 const Real charge = getObjectWrapper().particleSpecies[popID].charge;
572 #pragma omp parallel for schedule(static)
573 for (uint celli = 0; celli < nAllCells; celli++){
574 SpatialCell* cell = mpiGrid[cells[celli]];
575 if (cell->sysBoundaryFlag == sysboundarytype::DO_NOT_COMPUTE) {
576 continue;
577 }
578 if (cell->sysBoundaryFlag == sysboundarytype::OUTFLOW && cell->sysBoundaryLayer != 1 && !initialCompute) {
579 // these should have been handled by the boundary code
580 continue;
581 }
582
583 // Store species' contribution to bulk velocity moments
584 Population &pop = cell->get_population(popID);
585 pop.RHO_V = host_moments1[nMom1*celli];
586 pop.V_V[0] = divideIfNonZero(host_moments1[nMom1*celli + 1], host_moments1[nMom1*celli]);
587 pop.V_V[1] = divideIfNonZero(host_moments1[nMom1*celli + 2], host_moments1[nMom1*celli]);
588 pop.V_V[2] = divideIfNonZero(host_moments1[nMom1*celli + 3], host_moments1[nMom1*celli]);
589
590 cell->parameters[CellParams::RHOM_V ] += host_moments1[nMom1*celli]*mass;
591 cell->parameters[CellParams::VX_V] += host_moments1[nMom1*celli + 1]*mass;
592 cell->parameters[CellParams::VY_V] += host_moments1[nMom1*celli + 2]*mass;
593 cell->parameters[CellParams::VZ_V] += host_moments1[nMom1*celli + 3]*mass;
594 cell->parameters[CellParams::RHOQ_V ] += host_moments1[nMom1*celli]*charge;
595 } // for-loop over spatial cells
596 } // for-loop over particle species
597
598 #pragma omp parallel for schedule(static)
599 for (size_t celli=0; celli<nAllCells; ++celli) {
600 SpatialCell* cell = mpiGrid[cells[celli]];
601 if (cell->sysBoundaryFlag == sysboundarytype::DO_NOT_COMPUTE) {
602 continue;
603 }
604 if (cell->sysBoundaryFlag == sysboundarytype::OUTFLOW && cell->sysBoundaryLayer != 1 && !initialCompute) {
605 // these should have been handled by the boundary code
606 continue;
607 }
608 cell->parameters[CellParams::VX_V] = divideIfNonZero(cell->parameters[CellParams::VX_V], cell->parameters[CellParams::RHOM_V]);
609 cell->parameters[CellParams::VY_V] = divideIfNonZero(cell->parameters[CellParams::VY_V], cell->parameters[CellParams::RHOM_V]);
610 cell->parameters[CellParams::VZ_V] = divideIfNonZero(cell->parameters[CellParams::VZ_V], cell->parameters[CellParams::RHOM_V]);
611
612 // copy the bulk flow frame back to device
613 //host_moments1[nMom1*celli + 0] = cell->parameters[CellParams::RHOM_V];
614 host_moments1[nMom1*celli + 1] = cell->parameters[CellParams::VX_V];
615 host_moments1[nMom1*celli + 2] = cell->parameters[CellParams::VY_V];
616 host_moments1[nMom1*celli + 3] = cell->parameters[CellParams::VZ_V];
617 }
618
619 // Compute second moments only if requested.
620 if (computeSecond == false) {
621 gpuMemoryManager.endSession();
622 return;
623 }
624
625 CHK_ERR( gpuMemcpy(dev_moments1, host_moments1, nAllCells*nMom1*sizeof(Real), gpuMemcpyHostToDevice) );
626 CHK_ERR( gpuMemset(dev_moments2, 0, nAllCells*nMom2*sizeof(Real)) );
627 for (uint popID=0; popID<getObjectWrapper().particleSpecies.size(); ++popID) {
628 if (maxVmeshSizes.at(popID) == 0) {
629 continue;
630 }
631
632 vmesh::LocalID maxVmeshLaunch = (gpuMultiProcessorCount*min(threadsPerMP/WID3, blocksPerMP))/nAllCells;
633 maxVmeshLaunch = maxVmeshLaunch < 1 ? 1 : maxVmeshLaunch;
634
635 dim3 blockSize(WID,WID,WID);
636 dim3 gridSize(nAllCells,maxVmeshLaunch,1);
637 int sharedMemorySizeSecondMoments = nMom2 * (WID3 / GPUTHREADS) * sizeof(Real);
638 second_moments_kernel<<<gridSize, blockSize, sharedMemorySizeSecondMoments, 0>>> (
639 dev_VBC,
640 dev_moments1,
641 dev_moments2,
642 nAllCells
643 );
644 CHK_ERR( gpuPeekAtLastError() );
645 CHK_ERR( gpuDeviceSynchronize());
646 CHK_ERR( gpuMemcpy(host_moments2, dev_moments2, nAllCells*nMom2*sizeof(Real), gpuMemcpyDeviceToHost) );
647
648 const Real mass = getObjectWrapper().particleSpecies[popID].mass;
649 #pragma omp parallel for schedule(static)
650 for (uint celli = 0; celli < nAllCells; celli++){
651 SpatialCell* cell = mpiGrid[cells[celli]];
652 if (cell->sysBoundaryFlag == sysboundarytype::DO_NOT_COMPUTE) {
653 continue;
654 }
655 if (cell->sysBoundaryFlag == sysboundarytype::OUTFLOW && cell->sysBoundaryLayer != 1 && !initialCompute) {
656 // these should have been handled by the boundary code
657 continue;
658 }
659
660 // Store species' contribution to bulk velocity moments
661 Population &pop = cell->get_population(popID);
662 for (size_t i = 0; i < nMom2; ++i) {
663 pop.P_V[i] = mass*host_moments2[nMom2*celli + i];
664 }
665
666 cell->parameters[CellParams::P_11_V] += pop.P_V[0];
667 cell->parameters[CellParams::P_22_V] += pop.P_V[1];
668 cell->parameters[CellParams::P_33_V] += pop.P_V[2];
669 cell->parameters[CellParams::P_23_V] += pop.P_V[3];
670 cell->parameters[CellParams::P_13_V] += pop.P_V[4];
671 cell->parameters[CellParams::P_12_V] += pop.P_V[5];
672 } // for-loop over spatial cells
673 } // for-loop over particle species
674
675 gpuMemoryManager.endSession();
676}
i
for i
Definition Dispersion.m:24
gpuPeekAtLastError
#define gpuPeekAtLastError
Definition arch_device_cuda.h:18
gpuMemcpyHostToDevice
#define gpuMemcpyHostToDevice
Definition arch_device_cuda.h:87
CHK_ERR
#define CHK_ERR(err)
Definition arch_device_cuda.h:123
gpuMemcpy
#define gpuMemcpy
Definition arch_device_cuda.h:42
gpuMemcpyDeviceToHost
#define gpuMemcpyDeviceToHost
Definition arch_device_cuda.h:86
gpuKernelShflDown
#define gpuKernelShflDown(val, offset)
Definition arch_device_cuda.h:95
gpuMemset
#define gpuMemset
Definition arch_device_cuda.h:44
gpuDeviceSynchronize
#define gpuDeviceSynchronize
Definition arch_device_cuda.h:25
GPUTHREADS
#define GPUTHREADS
Definition arch_device_cuda.h:113
nMom1
#define nMom1
Definition arch_moments.h:37
nMom2
#define nMom2
Definition arch_moments.h:38
spatial_cell::SpatialCell
Definition spatial_cell_cpu.hpp:235
spatial_cell::SpatialCell::get_velocity_mesh
vmesh::VelocityMesh * get_velocity_mesh(const size_t &popID)
Definition spatial_cell_cpu.hpp:673
spatial_cell::SpatialCell::sysBoundaryLayer
uint sysBoundaryLayer
Definition spatial_cell_cpu.hpp:372
spatial_cell::SpatialCell::get_population
Population & get_population(const uint popID)
Definition spatial_cell_cpu.hpp:495
spatial_cell::SpatialCell::parameters
std::array< Real, CellParams::N_SPATIAL_CELL_PARAMS > parameters
Definition spatial_cell_cpu.hpp:361
spatial_cell::SpatialCell::sysBoundaryFlag
uint sysBoundaryFlag
Definition spatial_cell_cpu.hpp:370
spatial_cell::SpatialCell::dev_get_velocity_blocks
vmesh::VelocityBlockContainer * dev_get_velocity_blocks(const size_t &popID)
Definition spatial_cell_gpu.hpp:1074
vmesh::VelocityMesh::size
size_t size(bool dummy=0) const
Definition velocity_mesh_cpu.h:497
WID
#define WID
Definition common.h:514
WID3
const int WID3
Definition common.h:517
Real
float Real
Definition definitions.h:41
Realf
float Realf
Definition definitions.h:33
divideIfNonZero
Real divideIfNonZero(creal numerator, creal denominator)
Helper function.
Definition fs_common.cpp:33
fs_common.h
HALF
const Real HALF
Definition fs_common.h:49
ti
const uint ti
Definition gpu_acc_map.cpp:764
blocksPerMP
int blocksPerMP
Definition gpu_base.cpp:43
gpuMemoryManager
GPUMemoryManager gpuMemoryManager
Definition gpu_base.cpp:64
threadsPerMP
int threadsPerMP
Definition gpu_base.cpp:44
gpuMultiProcessorCount
int gpuMultiProcessorCount
Definition gpu_base.cpp:42
gpu_base.hpp
SESSION_HOST_ALLOCATE
#define SESSION_HOST_ALLOCATE(object, type, member, bytes)
Definition gpu_base.hpp:600
SESSION_ALLOCATE
#define SESSION_ALLOCATE(object, type, member, bytes)
Definition gpu_base.hpp:557
GET_SESSION_POINTER
#define GET_SESSION_POINTER(object, type, member)
Definition gpu_base.hpp:833
GET_SESSION_HOST_POINTER
#define GET_SESSION_HOST_POINTER(object, type, member)
Definition gpu_base.hpp:853
gpu_calculateMoments_V
void gpu_calculateMoments_V(dccrg::Dccrg< SpatialCell, dccrg::Cartesian_Geometry > &mpiGrid, const std::vector< CellID > &cells_in, const bool computeSecond, const bool initialCompute)
Definition gpu_moments.cpp:469
gpu_calculateMoments_R
void gpu_calculateMoments_R(dccrg::Dccrg< SpatialCell, dccrg::Cartesian_Geometry > &mpiGrid, const std::vector< CellID > &cells_in, const bool computeSecond, const bool initialCompute)
Definition gpu_moments.cpp:249
gpu_moments.h
VY
const Real VY
Definition gpu_pitch_angle_diffusion.cpp:69
cellValue
const Realf cellValue
Definition gpu_pitch_angle_diffusion.cpp:84
blockSize
const int blockSize
Definition gpu_pitch_angle_diffusion.cpp:392
VZ
const Real VZ
Definition gpu_pitch_angle_diffusion.cpp:70
j
const int j
Definition gpu_pitch_angle_diffusion.cpp:60
blockParameters
const Real *__restrict__ blockParameters
Definition gpu_pitch_angle_diffusion.cpp:65
__syncthreads
__syncthreads()
VX
const Real VX
Definition gpu_pitch_angle_diffusion.cpp:68
k
const int k
Definition gpu_pitch_angle_diffusion.cpp:61
getObjectWrapper
ObjectWrapper & getObjectWrapper()
Definition main.cpp:33
BlockParams::N_VELOCITY_BLOCK_PARAMS
@ N_VELOCITY_BLOCK_PARAMS
Definition common.h:115
arch::uint
uint32_t uint
Definition arch_device_api.h:17
faceestimates_amr.f
f
Definition faceestimates_amr.py:51
sort_refined_ids.offset
int offset
Definition sort_refined_ids.py:260
spatial_cell::__launch_bounds__
static __global__ void __launch_bounds__(WID3, 4) population_scale_kernel(vmesh
Definition spatial_cell_gpu.hpp:69
vmesh::LocalID
uint32_t LocalID
Definition definitions.h:60
object_wrapper.h
ObjectWrapper::particleSpecies
std::vector< species::Species > particleSpecies
Definition object_wrapper.h:48
spatial_cell::Population
Definition spatial_cell_cpu.hpp:66
spatial_cell::Population::V_V
Real V_V[3]
Definition spatial_cell_cpu.hpp:72
spatial_cell::Population::P_R
Real P_R[6]
Definition spatial_cell_cpu.hpp:74
spatial_cell::Population::RHO_V
Real RHO_V
Definition spatial_cell_cpu.hpp:71
spatial_cell::Population::RHO_R
Real RHO_R
Definition spatial_cell_cpu.hpp:69
spatial_cell::Population::V_R
Real V_R[3]
Definition spatial_cell_cpu.hpp:70
spatial_cell::Population::P_V
Real P_V[6]
Definition spatial_cell_cpu.hpp:75
min
static ARCH_HOSTDEV VecSimple< T > min(VecSimple< T > const &l, VecSimple< T > const &r)
Definition vectorclass_fallback.h:463