Vlasiator ebf0dd394 on dev (v5.4.0 + 1054 commits)
Loading...
Searching...
No Matches
velocity_mesh_gpu.h
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// Non-AMR implementation of the velocity space, still in use despite the filename
24
25#ifndef VELOCITY_MESH_GPU_H
26#define VELOCITY_MESH_GPU_H
27
28#include <iostream>
29#include <algorithm>
30#include <sstream>
31#include <stdint.h>
32#include <vector>
33#include <cmath>
34
35#include "../velocity_mesh_parameters.h"
36
37#include "include/hashinator/hashinator.h"
38#include "include/splitvector/splitvec.h"
39
40#include "../arch/gpu_base.hpp"
41//#include "../arch/arch_device_api.h" // included in above
42
43#if defined(DEBUG_VLASIATOR) || defined(DEBUG_SPATIAL_CELL)
44 #ifndef DEBUG_VMESH
45 #define DEBUG_VMESH
46 #endif
47#endif
48
49namespace vmesh {
50
51 class VelocityMesh {
52 public:
53 VelocityMesh();
54 ~VelocityMesh();
55 VelocityMesh(const VelocityMesh& other);
56 const VelocityMesh& operator=(const VelocityMesh& other);
57
58 ARCH_HOSTDEV size_t capacity() const;
59 size_t capacityInBytes() const;
60 ARCH_HOSTDEV bool check(); // These are no longer const as they prefetch back and forth
61 ARCH_HOSTDEV void print();
62 void clear(bool shrink);
63 ARCH_HOSTDEV bool move(const vmesh::LocalID sourceLocalID,const vmesh::LocalID targetLocalID);
64 ARCH_DEV bool warpMove(const vmesh::LocalID sourceLocalID,const vmesh::LocalID targetLocalID, const size_t b_tid);
65 ARCH_HOSTDEV size_t count(const vmesh::GlobalID globalID) const;
66 ARCH_DEV size_t warpCount(const vmesh::GlobalID globalID, const size_t b_tid) const;
67 ARCH_HOSTDEV vmesh::GlobalID findBlock(vmesh::GlobalID cellIndices[3]) const;
68 ARCH_DEV vmesh::GlobalID warpFindBlock(vmesh::GlobalID cellIndices[3], const size_t b_tid) const;
69 ARCH_HOSTDEV bool getBlockCoordinates(const vmesh::GlobalID globalID,Real coords[3]) const;
70 ARCH_HOSTDEV void getBlockInfo(const vmesh::GlobalID globalID,Real* array) const;
71 ARCH_HOSTDEV const Real* getBlockSize() const;
72 ARCH_HOSTDEV bool getBlockSize(const vmesh::GlobalID globalID,Real size[3]) const;
73 ARCH_HOSTDEV const Real* getCellSize() const;
74 ARCH_HOSTDEV bool getCellSize(const vmesh::GlobalID globalID,Real size[3]) const;
75 ARCH_HOSTDEV vmesh::GlobalID getGlobalID(const vmesh::LocalID localID) const;
76 ARCH_HOSTDEV vmesh::GlobalID getGlobalID(const Real* coords) const;
77 ARCH_HOSTDEV vmesh::GlobalID getGlobalID(const vmesh::LocalID indices[3]) const;
78 ARCH_HOSTDEV vmesh::GlobalID getGlobalID(const vmesh::LocalID i,const vmesh::LocalID j,const vmesh::LocalID k) const;
79 ARCH_HOSTDEV split::SplitVector<vmesh::GlobalID>* getGrid();
80 ARCH_HOSTDEV const vmesh::LocalID* getGridLength() const;
81 ARCH_HOSTDEV void getIndices(const vmesh::GlobalID globalID,vmesh::LocalID& i,vmesh::LocalID& j,vmesh::LocalID& k) const;
82 ARCH_HOSTDEV void getIndicesX(const vmesh::GlobalID globalID,vmesh::LocalID& i) const;
83 ARCH_HOSTDEV void getIndicesY(const vmesh::GlobalID globalID,vmesh::LocalID& j) const;
84 ARCH_HOSTDEV void getIndicesZ(const vmesh::GlobalID globalID,vmesh::LocalID& k) const;
85 ARCH_HOSTDEV size_t getMesh() const;
86 ARCH_HOSTDEV vmesh::LocalID getLocalID(const vmesh::GlobalID globalID) const;
87 ARCH_DEV vmesh::LocalID warpGetLocalID(const vmesh::GlobalID globalID, const size_t b_tid) const;
88 ARCH_HOSTDEV const Real* getMeshMaxLimits() const;
89 ARCH_HOSTDEV const Real* getMeshMinLimits() const;
90 ARCH_HOSTDEV bool initialize(const size_t meshID);
91 ARCH_HOSTDEV static vmesh::LocalID invalidBlockIndex();
92 ARCH_HOSTDEV static vmesh::GlobalID invalidGlobalID();
93 ARCH_HOSTDEV static vmesh::LocalID invalidLocalID();
94 ARCH_HOSTDEV bool isInitialized() const;
95 ARCH_HOSTDEV void pop();
96 ARCH_DEV void warpPop(const size_t b_tid);
97 ARCH_HOSTDEV bool push_back(const vmesh::GlobalID globalID);
98 ARCH_DEV bool warpPush_back(const vmesh::GlobalID globalID, const size_t b_tid);
99 vmesh::LocalID push_back(const std::vector<vmesh::GlobalID>& blocks);
100 ARCH_HOSTDEV vmesh::LocalID push_back(split::SplitVector<vmesh::GlobalID>* blocks);
101 ARCH_DEV vmesh::LocalID warpPush_back(const split::SplitVector<vmesh::GlobalID>& blocks, const size_t b_tid);
102 ARCH_DEV void replaceBlock(const vmesh::GlobalID GIDold,const vmesh::LocalID LID,const vmesh::GlobalID GIDnew);
103 ARCH_DEV void warpReplaceBlock(const vmesh::GlobalID GIDold,const vmesh::LocalID LID,const vmesh::GlobalID GIDnew, const size_t b_tid);
104 ARCH_DEV void placeBlock(const vmesh::GlobalID GID,const vmesh::LocalID LID);
105 ARCH_DEV void warpPlaceBlock(const vmesh::GlobalID GID,const vmesh::LocalID LID, const size_t b_tid);
106 ARCH_DEV void deleteBlock(const vmesh::GlobalID GID,const vmesh::LocalID LID);
107 ARCH_DEV void warpDeleteBlock(const vmesh::GlobalID GID,const vmesh::LocalID LID, const size_t b_tid);
108 void setGrid();
109 bool setGrid(const std::vector<vmesh::GlobalID>& globalIDs);
110 bool setGrid(const split::SplitVector<vmesh::GlobalID>& globalIDs);
111 bool setMesh(const size_t meshID);
112 size_t getMeshID();
113 void setNewSize(const vmesh::LocalID newSize);
114 // bool setNewSizeClear(const vmesh::LocalID& newSize, gpuStream_t stream);
115 ARCH_DEV void device_setNewSize(const vmesh::LocalID newSize);
116 void setNewCachedSize(const vmesh::LocalID newSize);
117 void updateCachedSize();
118 void updateCachedCapacity();
119 bool setNewCapacity(const vmesh::LocalID newCapacity);
120 ARCH_HOSTDEV size_t size() const;
121 ARCH_HOSTDEV size_t sizeInBytes() const;
122 // void swap(VelocityMesh& vm);
123
124 void gpu_prefetchHost(gpuStream_t stream);
125 void gpu_prefetchDevice(gpuStream_t stream);
126 void gpu_cleanHashMap(gpuStream_t stream);
127 void print_addresses();
128 void print_sizes();
129
130 ARCH_HOSTDEV Hashinator::Hashmap<vmesh::GlobalID,vmesh::LocalID>* gpu_expose_map();
131
132 private:
133 size_t meshID;
134 size_t ltg_size=0, ltg_capacity=0, gtl_sizepower=0; // host-cached values
135
136 Hashinator::Hashmap<vmesh::GlobalID,vmesh::LocalID> globalToLocalMap;
137 split::SplitVector<vmesh::GlobalID> localToGlobalMap;
138 };
139
140 // ***** DEFINITIONS OF MEMBER FUNCTIONS ***** //
141
142
143 inline VelocityMesh::VelocityMesh() {
144 meshID = std::numeric_limits<size_t>::max();
145 // Set sizepower to 10 (1024 blocks) straight away so there's enough room to grow?
146 globalToLocalMap = Hashinator::Hashmap<vmesh::GlobalID,vmesh::LocalID>(INIT_MAP_SIZE);
147 localToGlobalMap = split::SplitVector<vmesh::GlobalID>(INIT_VMESH_SIZE);
148 localToGlobalMap.clear();
149 ltg_size = 0;
150 ltg_capacity = INIT_VMESH_SIZE;
151 gtl_sizepower = INIT_MAP_SIZE;
152 }
153
154 inline VelocityMesh::~VelocityMesh() {}
155
156 inline VelocityMesh::VelocityMesh(const VelocityMesh& other) {
157 gpuStream_t stream = gpu_getStream();
158 meshID = other.meshID;
159 if (other.localToGlobalMap.size() > 0) {
160 globalToLocalMap = Hashinator::Hashmap<vmesh::GlobalID,vmesh::LocalID>(other.globalToLocalMap);
161 localToGlobalMap = split::SplitVector<vmesh::GlobalID>(other.localToGlobalMap.capacity());
162 // Overwrite is like a copy assign but takes a stream
163 localToGlobalMap.overwrite(other.localToGlobalMap,stream);
164 ltg_size = other.ltg_size;
165 ltg_capacity = localToGlobalMap.capacity();
166 gtl_sizepower = globalToLocalMap.getSizePower();
167 } else {
168 globalToLocalMap = Hashinator::Hashmap<vmesh::GlobalID,vmesh::LocalID>(INIT_MAP_SIZE);
169 localToGlobalMap = split::SplitVector<vmesh::GlobalID>(INIT_VMESH_SIZE);
170 localToGlobalMap.clear();
171 ltg_size = 0;
172 ltg_capacity = INIT_VMESH_SIZE;
173 gtl_sizepower = INIT_MAP_SIZE;
174 }
175 }
176
177 inline const VelocityMesh& VelocityMesh::operator=(const VelocityMesh& other) {
178 gpuStream_t stream = gpu_getStream();
179 meshID = other.meshID;
180 setNewCapacity(other.ltg_size);
181 // Overwrite is like a copy assign but takes a stream
182 globalToLocalMap.overwrite(other.globalToLocalMap,stream);
183 // This constructor is used e.g. for boundary cells, where we in fact
184 // don't need any extra capacity, so let's not call this reserve.
185 // localToGlobalMap->reserve((other.localToGlobalMap)->capacity(),true);
186 localToGlobalMap.overwrite(other.localToGlobalMap,stream);
187 ltg_size = other.ltg_size;
188 return *this;
189 }
190
191 inline size_t VelocityMesh::capacityInBytes() const {
192 #ifdef DEBUG_VMESH
193 const size_t cap1 = localToGlobalMap.capacity();
194 if (ltg_capacity != cap1) {
195 printf("VMESH CAPACITY ERROR: LTG capacity %lu vs cached value %lu in %s : %d\n",cap1,ltg_capacity,__FILE__,__LINE__);
196 }
197 const size_t cap2 = globalToLocalMap.getSizePower();
198 if (gtl_sizepower != cap2) {
199 printf("VMESH CAPACITY ERROR: GTL sizePower %lu vs cached value %lu in %s : %d\n",cap2,gtl_sizepower,__FILE__,__LINE__);
200 }
201 #endif
202 // *** Using cached values
203 const size_t currentCapacity = ltg_capacity;
204 const size_t currentBucketCount = std::pow(2,gtl_sizepower);
205
206 const size_t capacityInBytes = currentCapacity * sizeof(vmesh::GlobalID)
207 + currentBucketCount * sizeof(Hashinator::hash_pair<vmesh::GlobalID,vmesh::LocalID>);
208 return capacityInBytes;
209 }
210
211 ARCH_HOSTDEV inline size_t VelocityMesh::capacity() const {
212 #if defined(__CUDA_ARCH__) || defined(__HIP_DEVICE_COMPILE__)
213 return localToGlobalMap.capacity();
214 #else
215 #ifdef DEBUG_VMESH
216 const size_t cap1 = localToGlobalMap.capacity();
217 if (ltg_capacity != cap1) {
218 printf("VMESH CAPACITY ERROR: capacity %lu vs cached value %lu in %s : %d\n",cap1,ltg_capacity,__FILE__,__LINE__);
219 }
220 #endif
221 return ltg_capacity;
222 #endif
223 }
224
225 ARCH_HOSTDEV inline bool VelocityMesh::check() {
226 bool ok = true;
227 #if !defined(__CUDA_ARCH__) && !defined(__HIP_DEVICE_COMPILE__)
228 // Prefetch to host for the in-detail check
229 gpuStream_t stream = gpu_getStream();
230 localToGlobalMap.optimizeCPU(stream);
231 globalToLocalMap.optimizeCPU(stream);
232 //CHK_ERR( gpuStreamSynchronize(stream) );
233 CHK_ERR( gpuDeviceSynchronize() );
234 #endif
235 const size_t size1 = localToGlobalMap.size();
236 const size_t size2 = globalToLocalMap.size();
237 #if !defined(__CUDA_ARCH__) && !defined(__HIP_DEVICE_COMPILE__)
238 // Only check cached values on host
239 const size_t cap1 = localToGlobalMap.capacity();
240 const size_t cap2 = globalToLocalMap.getSizePower();
241 if (cap1 != ltg_capacity) {
242 printf("VMESH CHECK ERROR: capacity %lu vs cached value %lu in %s : %d\n",cap1,ltg_capacity,__FILE__,__LINE__);
243 return false;
244 }
245 if (cap2 != gtl_sizepower) {
246 printf("VMESH CHECK ERROR: sizepower %lu vs cached value %lu in %s : %d\n",cap2,gtl_sizepower,__FILE__,__LINE__);
247 return false;
248 }
249 if (size1 != ltg_size) {
250 printf("VMESH CHECK ERROR: size %lu vs cached value %lu in %s : %d\n",size1,ltg_size,__FILE__,__LINE__);
251 return false;
252 }
253 #endif
254 if (size1 != size2) {
255 printf("VMESH CHECK ERROR: LTG size %lu vs GTL size %lu in %s : %d\n",size1,size2,__FILE__,__LINE__);
256 return false;
257 //assert(0 && "VM check ERROR: sizes differ");
258 }
259 size_t fail = 0;
260 for (size_t b=0; b<size1; ++b) {
261 const vmesh::LocalID globalID = localToGlobalMap.at(b);
262 #if defined(__CUDA_ARCH__) || defined(__HIP_DEVICE_COMPILE__)
263 auto it = globalToLocalMap.device_find(globalID);
264 if (it != globalToLocalMap.device_end()) {
265 #else
266 auto it = globalToLocalMap.find(globalID);
267 if (it != globalToLocalMap.end()) {
268 #endif
269 const vmesh::GlobalID localID = it->second;
270 if (localID != b) {
271 ok = false;
272 printf("VMESH CHECK ERROR: localToGlobalMap[%lu] = %u but globalToLocalMap[%u] = %u\n",b,globalID,globalID,localID);
273 //assert(0 && "VM check ERROR");
274 fail++;
275 }
276 } else {
277 ok = false;
278 printf("VMESH CHECK ERROR: localToGlobalMap[%lu] = %u but could not find in globalToLocalMap ",b,globalID);
279 fail++;
280 }
281 }
282 if (fail>0) {
283 printf("VMESH CHECK ERROR Encountered %lu failures.\n",fail);
284 return false;
285 }
286 #if !defined(__CUDA_ARCH__) && !defined(__HIP_DEVICE_COMPILE__)
287 localToGlobalMap.optimizeGPU(stream);
288 globalToLocalMap.optimizeGPU(stream);
289 #endif
290 return ok;
291 }
292
293 ARCH_HOSTDEV inline void VelocityMesh::print() {
294 #if !defined(__CUDA_ARCH__) && !defined(__HIP_DEVICE_COMPILE__)
295 gpuStream_t stream = gpu_getStream();
296 localToGlobalMap.optimizeCPU(stream);
297 globalToLocalMap.optimizeCPU(stream);
298 #endif
299
300 if (localToGlobalMap.size() != globalToLocalMap.size()) {
301 printf("VMO ERROR: sizes differ, %lu vs %lu\n",localToGlobalMap.size(),globalToLocalMap.size());
302 assert(0 && "VM check ERROR: sizes differ");
303 }
304 vmesh::LocalID thisSize = size();
305 printf("VM Size: %u \n",thisSize);
306 for (vmesh::LocalID b=0; b<thisSize; ++b) {
307 const vmesh::LocalID globalID = localToGlobalMap.at(b);
308 #if defined(__CUDA_ARCH__) || defined(__HIP_DEVICE_COMPILE__)
309 auto it = globalToLocalMap.device_find(globalID);
310 #else
311 auto it = globalToLocalMap.find(globalID);
312 #endif
313 const vmesh::GlobalID localID = it->second;
314 printf("vmesh LID [%6u] => GID [%6u] => [%6u]\n",b,globalID,localID);
315 }
316 #if !defined(__CUDA_ARCH__) && !defined(__HIP_DEVICE_COMPILE__)
317 localToGlobalMap.optimizeGPU(stream);
318 globalToLocalMap.optimizeGPU(stream);
319 #endif
320 }
321
322 inline void VelocityMesh::clear(bool shrink=true) {
323 ltg_size = 0;
324 if (shrink) {
325 ltg_capacity = 1;
326 gtl_sizepower = 4;
327 globalToLocalMap = Hashinator::Hashmap<vmesh::GlobalID,vmesh::LocalID>(ltg_capacity);
328 localToGlobalMap = split::SplitVector<vmesh::GlobalID>(gtl_sizepower);
329 localToGlobalMap.clear();
330 } else {
331 gpuStream_t stream = gpu_getStream();
332 localToGlobalMap.clear();
333 globalToLocalMap.clear<false>(Hashinator::targets::device,stream, std::pow(2,gtl_sizepower));
334 CHK_ERR( gpuStreamSynchronize(stream) );
335 }
336
337 #ifdef DEBUG_VMESH
338 if ((localToGlobalMap.size() != 0) || (globalToLocalMap.size() != 0)) {
339 std::cerr<<"VMESH CLEAR FAILED"<<std::endl;
340 }
341 #endif
342 }
343
344 ARCH_HOSTDEV inline bool VelocityMesh::move(const vmesh::LocalID sourceLID,const vmesh::LocalID targetLID) {
345 const vmesh::GlobalID moveGID = localToGlobalMap.at(sourceLID); // block to move (at the end of list)
346 const vmesh::GlobalID removeGID = localToGlobalMap.at(targetLID); // removed block
347 #ifdef DEBUG_VMESH
348 if (sourceLID != size()-1) {
349 assert( 0 && "Error! Moving velocity mesh entry from position which is not last LID!");
350 }
351 #endif
352
353 // at-function will throw out_of_range exception for non-existing global ID:
354 #if defined(__CUDA_ARCH__) || defined(__HIP_DEVICE_COMPILE__)
355 globalToLocalMap.set_element(moveGID,targetLID);
356 globalToLocalMap.device_erase(removeGID);
357 #else
358 globalToLocalMap.at(moveGID) = targetLID;
359 globalToLocalMap.erase(removeGID);
360 #endif
361 localToGlobalMap.at(targetLID) = moveGID;
362 localToGlobalMap.pop_back();
363 // Update cached value
364 ltg_size--;
365 return true;
366 }
367 ARCH_HOSTDEV inline size_t VelocityMesh::count(const vmesh::GlobalID globalID) const {
368 #if defined(__CUDA_ARCH__) || defined(__HIP_DEVICE_COMPILE__)
369 return globalToLocalMap.device_count(globalID);
370 #else
371 return globalToLocalMap.count(globalID);
372 #endif
373 }
374 ARCH_HOSTDEV inline vmesh::GlobalID VelocityMesh::findBlock(vmesh::GlobalID cellIndices[3]) const {
375 // Calculate i/j/k indices of the block that would own the cell:
376 vmesh::GlobalID i_block = cellIndices[0] / (*(vmesh::getMeshWrapper()->velocityMeshes))[meshID].blockLength[0];
377 vmesh::GlobalID j_block = cellIndices[1] / (*(vmesh::getMeshWrapper()->velocityMeshes))[meshID].blockLength[1];
378 vmesh::GlobalID k_block = cellIndices[2] / (*(vmesh::getMeshWrapper()->velocityMeshes))[meshID].blockLength[2];
379
380 // Calculate block global ID:
381 vmesh::GlobalID blockGID = getGlobalID(i_block,j_block,k_block);
382
383 // If the block exists, return it:
384 #if defined(__CUDA_ARCH__) || defined(__HIP_DEVICE_COMPILE__)
385 if (globalToLocalMap.device_count(blockGID) != 0) {
386 #else
387 if (globalToLocalMap.count(blockGID) != 0) {
388 #endif
389 return blockGID;
390 } else {
391 return invalidGlobalID();
392 }
393 }
394
395 ARCH_HOSTDEV inline bool VelocityMesh::getBlockCoordinates(const vmesh::GlobalID globalID,Real coords[3]) const {
396 if (globalID == invalidGlobalID()) {
397 for (int i=0; i<3; ++i) coords[i] = std::numeric_limits<Real>::quiet_NaN();
398 return false;
399 }
400
401 vmesh::LocalID indices[3];
402 getIndices(globalID,indices[0],indices[1],indices[2]);
403 if (indices[0] == invalidBlockIndex()) {
404 for (int i=0; i<3; ++i) coords[i] = std::numeric_limits<Real>::quiet_NaN();
405 return false;
406 }
407
408 coords[0] = (*(vmesh::getMeshWrapper()->velocityMeshes))[meshID].meshMinLimits[0] + indices[0]*(*(vmesh::getMeshWrapper()->velocityMeshes))[meshID].blockSize[0];
409 coords[1] = (*(vmesh::getMeshWrapper()->velocityMeshes))[meshID].meshMinLimits[1] + indices[1]*(*(vmesh::getMeshWrapper()->velocityMeshes))[meshID].blockSize[1];
410 coords[2] = (*(vmesh::getMeshWrapper()->velocityMeshes))[meshID].meshMinLimits[2] + indices[2]*(*(vmesh::getMeshWrapper()->velocityMeshes))[meshID].blockSize[2];
411 return true;
412 }
413 ARCH_HOSTDEV inline void VelocityMesh::getBlockInfo(const vmesh::GlobalID globalID, Real* array) const {
414 #ifdef DEBUG_VMESH
415 if (globalID == invalidGlobalID()) {
416 for (int i=0; i<6; ++i) array[i] = std::numeric_limits<Real>::infinity();
417 }
418 #endif
419
420 vmesh::LocalID indices[3];
421 indices[0] = globalID % (*(vmesh::getMeshWrapper()->velocityMeshes))[meshID].gridLength[0];
422 indices[1] = (globalID / (*(vmesh::getMeshWrapper()->velocityMeshes))[meshID].gridLength[0]) % (*(vmesh::getMeshWrapper()->velocityMeshes))[meshID].gridLength[1];
423 indices[2] = globalID / ((*(vmesh::getMeshWrapper()->velocityMeshes))[meshID].gridLength[0] * (*(vmesh::getMeshWrapper()->velocityMeshes))[meshID].gridLength[1]);
424
425 // Indices 0-2 contain coordinates of the lower left corner.
426 // The values are the same as if getBlockCoordinates(globalID,&(array[0])) was called
427 array[0] = (*(vmesh::getMeshWrapper()->velocityMeshes))[meshID].meshMinLimits[0] + indices[0]*(*(vmesh::getMeshWrapper()->velocityMeshes))[meshID].blockSize[0];
428 array[1] = (*(vmesh::getMeshWrapper()->velocityMeshes))[meshID].meshMinLimits[1] + indices[1]*(*(vmesh::getMeshWrapper()->velocityMeshes))[meshID].blockSize[1];
429 array[2] = (*(vmesh::getMeshWrapper()->velocityMeshes))[meshID].meshMinLimits[2] + indices[2]*(*(vmesh::getMeshWrapper()->velocityMeshes))[meshID].blockSize[2];
430
431 // Indices 3-5 contain the cell size.
432 // The values are the same as if getCellSize(globalID,&(array[3])) was called
433 array[3] = (*(vmesh::getMeshWrapper()->velocityMeshes))[meshID].cellSize[0];
434 array[4] = (*(vmesh::getMeshWrapper()->velocityMeshes))[meshID].cellSize[1];
435 array[5] = (*(vmesh::getMeshWrapper()->velocityMeshes))[meshID].cellSize[2];
436 }
437
438 ARCH_HOSTDEV inline const Real* VelocityMesh::getBlockSize() const {
439 return (*(vmesh::getMeshWrapper()->velocityMeshes))[meshID].blockSize;
440 }
441
448
449 ARCH_HOSTDEV inline const Real* VelocityMesh::getCellSize() const {
450 return (*(vmesh::getMeshWrapper()->velocityMeshes))[meshID].cellSize;
451 }
452
453 ARCH_HOSTDEV inline bool VelocityMesh::getCellSize(const vmesh::GlobalID globalID, Real size[3]) const {
454 size[0] = (*(vmesh::getMeshWrapper()->velocityMeshes))[meshID].cellSize[0];
455 size[1] = (*(vmesh::getMeshWrapper()->velocityMeshes))[meshID].cellSize[1];
456 size[2] = (*(vmesh::getMeshWrapper()->velocityMeshes))[meshID].cellSize[2];
457 return true;
458 }
459
460 ARCH_HOSTDEV inline vmesh::GlobalID VelocityMesh::getGlobalID(const vmesh::LocalID localID) const {
461 #ifdef DEBUG_VMESH
462 if (localID >= localToGlobalMap.size()) {
463 assert (0 && "ERROR invalid local id");
464 }
465 #endif
466
467 return localToGlobalMap.at(localID);
468 }
469
470 ARCH_HOSTDEV inline vmesh::GlobalID VelocityMesh::getGlobalID(const Real* coords) const {
471 if (coords[0] < (*(vmesh::getMeshWrapper()->velocityMeshes))[meshID].meshMinLimits[0] || coords[0] >= (*(vmesh::getMeshWrapper()->velocityMeshes))[meshID].meshMaxLimits[0] ||
472 (coords[1] < (*(vmesh::getMeshWrapper()->velocityMeshes))[meshID].meshMinLimits[1] || coords[1] >= (*(vmesh::getMeshWrapper()->velocityMeshes))[meshID].meshMaxLimits[1] ||
473 coords[2] < (*(vmesh::getMeshWrapper()->velocityMeshes))[meshID].meshMinLimits[2] || coords[2] >= (*(vmesh::getMeshWrapper()->velocityMeshes))[meshID].meshMaxLimits[2])) {
474 return invalidGlobalID();
475 }
476
477 const vmesh::LocalID indices[3] = {
478 static_cast<vmesh::LocalID>(floor((coords[0] - (*(vmesh::getMeshWrapper()->velocityMeshes))[meshID].meshMinLimits[0]) / (*(vmesh::getMeshWrapper()->velocityMeshes))[meshID].blockSize[0])),
479 static_cast<vmesh::LocalID>(floor((coords[1] - (*(vmesh::getMeshWrapper()->velocityMeshes))[meshID].meshMinLimits[1]) / (*(vmesh::getMeshWrapper()->velocityMeshes))[meshID].blockSize[1])),
480 static_cast<vmesh::LocalID>(floor((coords[2] - (*(vmesh::getMeshWrapper()->velocityMeshes))[meshID].meshMinLimits[2]) / (*(vmesh::getMeshWrapper()->velocityMeshes))[meshID].blockSize[2]))
481 };
482
483 return indices[2]*(*(vmesh::getMeshWrapper()->velocityMeshes))[meshID].gridLength[1]*(*(vmesh::getMeshWrapper()->velocityMeshes))[meshID].gridLength[0]
484 + indices[1]*(*(vmesh::getMeshWrapper()->velocityMeshes))[meshID].gridLength[0] + indices[0];
485 }
486
487 ARCH_HOSTDEV inline vmesh::GlobalID VelocityMesh::getGlobalID(const vmesh::LocalID indices[3]) const {
488 if (indices[0] >= (*(vmesh::getMeshWrapper()->velocityMeshes))[meshID].gridLength[0]) {
489 return invalidGlobalID();
490 }
491 if (indices[1] >= (*(vmesh::getMeshWrapper()->velocityMeshes))[meshID].gridLength[1]) {
492 return invalidGlobalID();
493 }
494 if (indices[2] >= (*(vmesh::getMeshWrapper()->velocityMeshes))[meshID].gridLength[2]) {
495 return invalidGlobalID();
496 }
497 return indices[0] + indices[1] * (*(vmesh::getMeshWrapper()->velocityMeshes))[meshID].gridLength[0]
498 + indices[2] * (*(vmesh::getMeshWrapper()->velocityMeshes))[meshID].gridLength[1]
499 * (*(vmesh::getMeshWrapper()->velocityMeshes))[meshID].gridLength[0];
500 }
501
502 ARCH_HOSTDEV inline vmesh::GlobalID VelocityMesh::getGlobalID(const vmesh::LocalID i,const vmesh::LocalID j,const vmesh::LocalID k) const {
503 if (i >= (*(vmesh::getMeshWrapper()->velocityMeshes))[meshID].gridLength[0]) {
504 return invalidGlobalID();
505 }
506 if (j >= (*(vmesh::getMeshWrapper()->velocityMeshes))[meshID].gridLength[1]) {
507 return invalidGlobalID();
508 }
509 if (k >= (*(vmesh::getMeshWrapper()->velocityMeshes))[meshID].gridLength[2]) {
510 return invalidGlobalID();
511 }
512 return i + j * (*(vmesh::getMeshWrapper()->velocityMeshes))[meshID].gridLength[0]
513 + k * (*(vmesh::getMeshWrapper()->velocityMeshes))[meshID].gridLength[1]
514 * (*(vmesh::getMeshWrapper()->velocityMeshes))[meshID].gridLength[0];
515 }
516
517 ARCH_HOSTDEV inline split::SplitVector<vmesh::GlobalID>* VelocityMesh::getGrid() {
518 return &localToGlobalMap;
519 }
520
521 ARCH_HOSTDEV inline const vmesh::LocalID* VelocityMesh::getGridLength() const {
522 return (*(vmesh::getMeshWrapper()->velocityMeshes))[meshID].gridLength;
523 }
524
525 ARCH_HOSTDEV inline void VelocityMesh::getIndices(const vmesh::GlobalID globalID,vmesh::LocalID& i,vmesh::LocalID& j,vmesh::LocalID& k) const {
526 if (globalID >= invalidGlobalID()) {
527 i = j = k = invalidBlockIndex();
528 } else {
529 i = globalID % (*(vmesh::getMeshWrapper()->velocityMeshes))[meshID].gridLength[0];
530 j = (globalID / (*(vmesh::getMeshWrapper()->velocityMeshes))[meshID].gridLength[0]) % (*(vmesh::getMeshWrapper()->velocityMeshes))[meshID].gridLength[1];
531 k = globalID / ((*(vmesh::getMeshWrapper()->velocityMeshes))[meshID].gridLength[0] * (*(vmesh::getMeshWrapper()->velocityMeshes))[meshID].gridLength[1]);
532 }
533 }
534
535 ARCH_HOSTDEV inline void VelocityMesh::getIndicesX(const vmesh::GlobalID globalID,vmesh::LocalID& i) const {
536 if (globalID >= invalidGlobalID()) {
537 i = invalidBlockIndex();
538 } else {
539 i = globalID % (*(vmesh::getMeshWrapper()->velocityMeshes))[meshID].gridLength[0];
540 }
541 }
542 ARCH_HOSTDEV inline void VelocityMesh::getIndicesY(const vmesh::GlobalID globalID,vmesh::LocalID& j) const {
543 if (globalID >= invalidGlobalID()) {
544 j = invalidBlockIndex();
545 } else {
546 j = (globalID / (*(vmesh::getMeshWrapper()->velocityMeshes))[meshID].gridLength[0]) % (*(vmesh::getMeshWrapper()->velocityMeshes))[meshID].gridLength[1];
547 }
548 }
549 ARCH_HOSTDEV inline void VelocityMesh::getIndicesZ(const vmesh::GlobalID globalID,vmesh::LocalID& k) const {
550 if (globalID >= invalidGlobalID()) {
551 k = invalidBlockIndex();
552 } else {
553 k = globalID / ((*(vmesh::getMeshWrapper()->velocityMeshes))[meshID].gridLength[0] * (*(vmesh::getMeshWrapper()->velocityMeshes))[meshID].gridLength[1]);
554 }
555 }
556
557 ARCH_HOSTDEV inline vmesh::LocalID VelocityMesh::getLocalID(const vmesh::GlobalID globalID) const {
558 #if defined(__CUDA_ARCH__) || defined(__HIP_DEVICE_COMPILE__)
559 auto it = globalToLocalMap.device_find(globalID);
560 if (it != globalToLocalMap.device_end()) {
561 return it->second;
562 }
563 #else
564 auto it = globalToLocalMap.find(globalID);
565 if (it != globalToLocalMap.end()) {
566 return it->second;
567 }
568 #endif
569 return invalidLocalID();
570 }
571 ARCH_HOSTDEV inline size_t VelocityMesh::getMesh() const {
572 return meshID;
573 }
574
575 ARCH_HOSTDEV inline const Real* VelocityMesh::getMeshMaxLimits() const {
576 return (*(vmesh::getMeshWrapper()->velocityMeshes))[meshID].meshMaxLimits;
577 }
578
579 ARCH_HOSTDEV inline const Real* VelocityMesh::getMeshMinLimits() const {
580 return (*(vmesh::getMeshWrapper()->velocityMeshes))[meshID].meshMinLimits;
581 }
582
583 ARCH_HOSTDEV inline bool VelocityMesh::initialize(const size_t meshID) {
584 this->meshID = meshID;
585 return true;
586 }
587
588 ARCH_HOSTDEV inline vmesh::LocalID VelocityMesh::invalidBlockIndex() {
589 return INVALID_VEL_BLOCK_INDEX;
590 }
591
592 ARCH_HOSTDEV inline vmesh::GlobalID VelocityMesh::invalidGlobalID() {
593 return INVALID_GLOBALID;
594 }
595
596 ARCH_HOSTDEV inline vmesh::LocalID VelocityMesh::invalidLocalID() {
597 return INVALID_LOCALID;
598 }
599
600 ARCH_HOSTDEV inline bool VelocityMesh::isInitialized() const {
601 return (*(vmesh::getMeshWrapper()->velocityMeshes))[meshID].initialized;
602 }
603
604 ARCH_HOSTDEV inline void VelocityMesh::pop() {
605 const size_t mySize = size();
606 if (mySize == 0) {
607 return;
608 }
609 const vmesh::LocalID lastLID = mySize-1;
610
611 #ifdef DEBUG_VMESH
612 const vmesh::GlobalID lastGID = localToGlobalMap.at(lastLID);
613 #else
614 const vmesh::GlobalID lastGID = localToGlobalMap[lastLID];
615 #endif
616 #if defined(__CUDA_ARCH__) || defined(__HIP_DEVICE_COMPILE__)
617 auto last = globalToLocalMap.device_find(lastGID);
618 globalToLocalMap.device_erase(last);
619 #else
620 auto last = globalToLocalMap.find(lastGID);
621 globalToLocalMap.erase(last);
622 #endif
623 localToGlobalMap.pop_back();
624 // Update cached value
625 ltg_size--;
626 }
629 ARCH_HOSTDEV inline bool VelocityMesh::push_back(const vmesh::GlobalID globalID) {
630 const size_t mySize = size();
631 if (mySize >= (*(vmesh::getMeshWrapper()->velocityMeshes))[meshID].max_velocity_blocks) {
632 return false;
633 }
634 if (globalID == invalidGlobalID()) {
635 return false;
636 }
637
638 #if defined(__CUDA_ARCH__) || defined(__HIP_DEVICE_COMPILE__)
639 // device_insert is slower, returns iterator and true or false for whether inserted key was new
640 const bool newEntry = globalToLocalMap.set_element<true>(globalID,(vmesh::LocalID)mySize);
641 if (newEntry) {
642 localToGlobalMap.device_push_back(globalID);
643 ltg_size++; // Note: called from inside kernel, cached size must be updated separately
644 //ltg_capacity = localToGlobalMap.capacity(); // on-device, no recapacitate
645 }
646 return newEntry;
647 #else
648 auto position
649 = globalToLocalMap.insert(Hashinator::make_pair(globalID,(vmesh::LocalID)mySize));
650 if (position.second == true) {
651 localToGlobalMap.push_back(globalID); // May increase capacity
652 ltg_size++;
653 if (ltg_size > ltg_capacity) {
654 ltg_capacity = localToGlobalMap.capacity();
655 }
656 gtl_sizepower = globalToLocalMap.getSizePower();
657 }
658 return position.second;
659 #endif
660 }
661 inline vmesh::LocalID VelocityMesh::push_back(const std::vector<vmesh::GlobalID>& blocks) {
662 gpuStream_t stream = gpu_getStream();
663 const size_t blocksSize = blocks.size();
664 if (ltg_size+blocksSize > (*(vmesh::getMeshWrapper()->velocityMeshes))[meshID].max_velocity_blocks) {
665 printf("vmesh: too many blocks, current size is %lu",ltg_size);
666 printf(", adding %lu blocks", blocksSize);
667 printf(", max is %u\n",(*(vmesh::getMeshWrapper()->velocityMeshes))[meshID].max_velocity_blocks);
668 return false;
669 }
670
671 if (ltg_size==0) {
672 // Fast insertion into empty mesh
673 if (blocksSize > ltg_capacity) {
674 ltg_capacity = blocksSize*BLOCK_ALLOCATION_FACTOR;
675 localToGlobalMap.reserve(ltg_capacity,true,stream);
676 }
677 localToGlobalMap.insert(localToGlobalMap.end(),blocks.begin(),blocks.end());
678 vmesh::GlobalID* _localToGlobalMapData = localToGlobalMap.data();
679 localToGlobalMap.optimizeGPU(stream);
680 globalToLocalMap.insertIndex<false>(_localToGlobalMapData,blocksSize,0.5,stream);
681 ltg_size = blocksSize;
682 gtl_sizepower = globalToLocalMap.getSizePower();
683 return blocksSize;
684 } else {
685 // GPUTODO: do inside kernel? This function is used by SpatialCell::add_velocity_blocks.
686 if (ltg_size+blocksSize > ltg_capacity) {
687 ltg_capacity = (ltg_size+blocksSize)*BLOCK_ALLOCATION_FACTOR;
688 localToGlobalMap.reserve(ltg_capacity,true,stream);
689 }
690 localToGlobalMap.resize(ltg_size+blocksSize,true,stream);
691 size_t newElements = 0;
692 for (size_t b=0; b<blocksSize; ++b) {
693 auto position
694 = globalToLocalMap.insert(Hashinator::make_pair(blocks[b],(vmesh::LocalID)(ltg_size+b)));
695 // Verify insertion into map and update vector
696 if (position.second) { // this is true if the element did not previously exist in the map
697 localToGlobalMap.at(ltg_size+newElements) = blocks[b];
698 newElements++;
699 }
700 }
701 localToGlobalMap.resize(ltg_size+newElements,true,stream);
702 ltg_size += newElements;
703 gtl_sizepower = globalToLocalMap.getSizePower();
704 localToGlobalMap.optimizeGPU(stream);
705 return newElements;
706 }
707 }
708
709 ARCH_HOSTDEV inline vmesh::LocalID VelocityMesh::push_back(split::SplitVector<vmesh::GlobalID>* blocks) {
710 #if !(defined(__CUDA_ARCH__) || defined(__HIP_DEVICE_COMPILE__))
711 gpuStream_t stream = gpu_getStream();
712 localToGlobalMap.optimizeCPU(stream); // insert one-by-one on CPU
713 #endif
714 const size_t blocksSize = blocks->size();
715 if (ltg_size+blocksSize > (*(vmesh::getMeshWrapper()->velocityMeshes))[meshID].max_velocity_blocks) {
716 printf("vmesh: too many blocks, current size is %lu",ltg_size);
717 printf(", adding %lu blocks", blocksSize);
718 printf(", max is %u\n",(*(vmesh::getMeshWrapper()->velocityMeshes))[meshID].max_velocity_blocks);
719 return false;
720 }
721
722 #if defined(__CUDA_ARCH__) || defined(__HIP_DEVICE_COMPILE__)
723 localToGlobalMap.device_resize(ltg_size+blocksSize, false); //construct=false don't construct or set to zero
724 size_t newElements = 0;
725 for (size_t b=0; b<blocksSize; ++b) {
726 // device_insert is slower than set_element, returns iterator and true or false for whether inserted key was new
727 const bool newEntry = globalToLocalMap.set_element<true>((*blocks)[b],(vmesh::LocalID)(ltg_size+b));
728 // Verify insertion into map and update vector
729 if (newEntry) { // this is true if the element did not previously exist in the map
730 localToGlobalMap[ltg_size+newElements] = (*blocks)[b];
731 newElements++;
732 }
733 }
734 localToGlobalMap.device_resize(ltg_size+newElements); //only make smaller so no construct
735 ltg_size += newElements; // Note: called from inside kernel, cached size must be updated separately
736 //ltg_capacity = localToGlobalMap.capacity(); // on-device, no recapacitate
737 return newElements;
738 #else
739 if (ltg_size==0) {
740 // Fast insertion into empty mesh
741 if (blocksSize > ltg_capacity) {
742 ltg_capacity = blocksSize*BLOCK_ALLOCATION_FACTOR;
743 localToGlobalMap.reserve(ltg_capacity,true,stream);
744 }
745 localToGlobalMap.insert(localToGlobalMap.end(),blocks->begin(),blocks->end());
746 vmesh::GlobalID* _localToGlobalMapData = localToGlobalMap.data();
747 localToGlobalMap.optimizeGPU(stream);
748 globalToLocalMap.insertIndex<false>(_localToGlobalMapData,blocksSize,0.5,stream);
749 ltg_size = blocksSize;
750 ltg_capacity = localToGlobalMap.capacity();
751 gtl_sizepower = globalToLocalMap.getSizePower();
752 return blocksSize;
753 } else {
754 // GPUTODO: do inside kernel?
755 if (ltg_size+blocksSize > ltg_capacity) {
756 ltg_capacity = (ltg_size+blocksSize)*BLOCK_ALLOCATION_FACTOR;
757 localToGlobalMap.reserve(ltg_capacity,true,stream);
758 }
759 localToGlobalMap.resize(ltg_size+blocksSize,true,stream);
760 size_t newElements = 0;
761 for (size_t b=0; b<blocksSize; ++b) {
762 auto position
763 = globalToLocalMap.insert(Hashinator::make_pair((*blocks)[b],(vmesh::LocalID)(ltg_size+b)));
764 // Verify insertion into map and update vector
765 if (position.second) { // this is true if the element did not previously exist in the map
766 localToGlobalMap.at(ltg_size+newElements) = (*blocks)[b];
767 newElements++;
768 }
769 }
770 localToGlobalMap.resize(ltg_size+newElements,true,stream);
771 localToGlobalMap.optimizeGPU(stream);
772 ltg_size += newElements;
773 ltg_capacity = localToGlobalMap.capacity();
774 gtl_sizepower = globalToLocalMap.getSizePower();
775 return newElements;
776 }
777 #endif
778 }
779
780#if defined(__CUDA_ARCH__) || defined(__HIP_DEVICE_COMPILE__)
781 /****
782 Device-only accessors to be called from a single GPU thread
783 **/
784 ARCH_DEV inline void VelocityMesh::replaceBlock(const vmesh::GlobalID GIDold,const vmesh::LocalID LID,const vmesh::GlobalID GIDnew) {
785 #ifdef DEBUG_VMESH
786 if (LID > size()-1) {
787 printf("vmesh replaceBlock error: LID is too large!\n");
788 }
789 vmesh::LocalID LIDold = invalidLocalID();
790 auto it = globalToLocalMap.device_find(GIDold);
791 if (it != globalToLocalMap.device_end()) {
792 LIDold = it->second;
793 }
794 if (localToGlobalMap.at(LIDold) != GIDold) {
795 printf("vmesh replaceBlock error: oldGID and oldLID don't match!\n");
796 }
797 #endif
798 globalToLocalMap.device_erase(GIDold);
799 globalToLocalMap.set_element(GIDnew,LID);
800 localToGlobalMap.at(LID) = GIDnew;
801 }
802 // Note: this function does not adjust the vmesh size, as it is used from within a parallel kernel.
803 ARCH_DEV inline void VelocityMesh::placeBlock(const vmesh::GlobalID GID,const vmesh::LocalID LID) {
804 #ifdef DEBUG_VMESH
805 if (LID > size()-1) {
806 assert(0 && "vmesh placeBlock error: LID is too large!");
807 }
808 bool newEntry = globalToLocalMap.set_element(GID,LID);
809 if (!newEntry) {
810 assert(0 && "vmesh placeBlock error: set_element not a new entry!");
811 }
812 localToGlobalMap.at(LID) = GID;
813 #else
814 globalToLocalMap.set_element(GID,LID);
815 localToGlobalMap[LID] = GID;
816 #endif
817 }
818 ARCH_DEV inline void VelocityMesh::deleteBlock(const vmesh::GlobalID GID,const vmesh::LocalID LID) {
819 #ifdef DEBUG_VMESH
820 // Verify that GID and LID match
821 if (GID==invalidGlobalID()) {
822 printf("vmesh deleteBlock error: GID is invalidGlobalID! (LID %ul)\n",LID);
823 }
824 if (LID==invalidLocalID()) {
825 printf("vmesh deleteBlock error: LID is invalidLocalID! (GID %ul)\n",GID);
826 }
827 auto it = globalToLocalMap.device_find(GID);
828 if (it == globalToLocalMap.device_end()) {
829 printf("vmesh deleteBlock error: GID %ul does not exist! (LID %ul)\n",GID,LID);
830 } else {
831 if (it->second != LID) {
832 printf("vmesh deleteBlock error: LID %ul found with GID %ul does not match provided LID %ul!\n",it->second,GID,LID);
833 }
834 }
835 if (localToGlobalMap.at(LID) != GID) {
836 printf("vmesh deleteBlock error: GID %ul found with LID %ul does not match provided GID %ul!\n",localToGlobalMap.at(LID),LID,GID);
837 }
838 #endif
839 globalToLocalMap.device_erase(GID);
840 #ifdef DEBUG_VMESH
841 localToGlobalMap.at(LID) = invalidGlobalID();
842 #else
843 localToGlobalMap[LID] = invalidGlobalID();
844 #endif
845 ltg_size--;
846 }
847
848 /****
849 Warp accessor functions to be called from within GPU kernels over several threads
850 **/
851 ARCH_DEV inline void VelocityMesh::warpPop(const size_t b_tid) {
852 const size_t mySize = size();
853 if (mySize == 0) {
854 return;
855 }
856 const vmesh::LocalID lastLID = mySize-1;
857 #ifdef DEBUG_VMESH
858 const vmesh::GlobalID lastGID = localToGlobalMap.at(lastLID);
859 const vmesh::LocalID mapSize = globalToLocalMap.size();
860 #else
861 const vmesh::GlobalID lastGID = localToGlobalMap[lastLID];
862 #endif
863 if (b_tid < GPUTHREADS) {
864 globalToLocalMap.warpErase(lastGID, b_tid);
865 }
866 if (b_tid==0) {
867 localToGlobalMap.pop_back();
868 }
869 #ifdef DEBUG_VMESH
870 __syncthreads();
871 const vmesh::LocalID newMapSize = globalToLocalMap.size();
872 const vmesh::LocalID newVecSize = localToGlobalMap.size();
873 if (newMapSize != mapSize-1) {
874 printf("warpError in VelocityMesh::warpPop: map size %u is not expected %u! (thread %u)\n",newMapSize,(vmesh::LocalID)(mapSize-1),(vmesh::LocalID)b_tid);
875 assert(0);
876 }
877 if (newVecSize != mySize-1) {
878 printf("warpError in VelocityMesh::warpPop: vector size %u is not expected %u! (thread %u)\n",newVecSize,(vmesh::LocalID)(mySize-1),(vmesh::LocalID)b_tid);
879 assert(0);
880 }
881 #endif
882 if (b_tid==0) {
883 ltg_size--;
884 }
885 __syncthreads();
886 }
887 ARCH_DEV inline vmesh::LocalID VelocityMesh::warpGetLocalID(const vmesh::GlobalID globalID, const size_t b_tid) const {
888 vmesh::LocalID retval = invalidLocalID();
889 globalToLocalMap.warpFind(globalID, retval, b_tid % GPUTHREADS);
890 #ifdef DEBUG_VMESH
891 auto it = globalToLocalMap.device_find(globalID);
892 if (it == globalToLocalMap.device_end()) {
893 if (retval != invalidLocalID()) {
894 printf("Warp error in VelocityMesh::warpGetLocalID: thread %u search did not find entry, warp search found LID %u for GID %u\n",(vmesh::LocalID)b_tid,retval,globalID);
895 }
896 } else if (retval != it->second) {
897 printf("Warp error in VelocityMesh::warpGetLocalID: LID %u (warp) != %u (thread %u) for GID %u\n",
898 retval,it->second,(vmesh::LocalID)b_tid,globalID);
899 }
900 #endif
901 __syncthreads();
902 return retval;
903 }
904 ARCH_DEV inline bool VelocityMesh::warpMove(const vmesh::LocalID sourceLID,const vmesh::LocalID targetLID, const size_t b_tid) {
905 #ifdef DEBUG_VMESH
906 const vmesh::GlobalID moveGID = localToGlobalMap.at(sourceLID); // block to move (must at the end of list)
907 const vmesh::GlobalID removeGID = localToGlobalMap.at(targetLID); // removed block
908 if (sourceLID != size()-1) {
909 assert( 0 && "Error! Moving velocity mesh entry from position which is not last LID!");
910 }
911 const vmesh::LocalID preMapSize = globalToLocalMap.size();
912 const vmesh::LocalID preVecSize = localToGlobalMap.size();
913 #else
914 const vmesh::GlobalID moveGID = localToGlobalMap[sourceLID]; // block to move (must at the end of list)
915 const vmesh::GlobalID removeGID = localToGlobalMap[targetLID]; // removed block
916 #endif
917
918 // at-function will throw out_of_range exception for non-existing global ID:
919 if (b_tid < GPUTHREADS) {
920 globalToLocalMap.warpInsert(moveGID,targetLID,b_tid); // will overwrite
921 globalToLocalMap.warpErase(removeGID,b_tid);
922 }
923 if (b_tid==0) {
924 localToGlobalMap.at(targetLID) = moveGID;
925 localToGlobalMap.pop_back();
926 }
927 #ifdef DEBUG_VMESH
928 __syncthreads();
929 const vmesh::LocalID postMapSize = globalToLocalMap.size();
930 const vmesh::LocalID postVecSize = localToGlobalMap.size();
931 const vmesh::LocalID postLID = getLocalID(moveGID);
932 if (preMapSize-1 != postMapSize) {
933 printf("warpError in VelocityMesh::warpMove: map size %u is not expected %u! (thread %u)\n",postMapSize,preMapSize,(vmesh::LocalID)b_tid);
934 assert(0);
935 }
936 if (preVecSize-1 != postVecSize) {
937 printf("warpError in VelocityMesh::warpMove: vector size %u is not expected %u! (thread %u)\n",postVecSize,preVecSize,(vmesh::LocalID)b_tid);
938 assert(0);
939 }
940 if (targetLID != postLID) {
941 printf("warpError in VelocityMesh::warpMove: Entry at GID %u is %u instead of expected %u! (thread %u)\n",moveGID,postLID,targetLID,(vmesh::LocalID)b_tid);
942 assert(0);
943 }
944 if (count(removeGID) != 0) {
945 const vmesh::LocalID removedLIDfound = getLocalID(removeGID);
946 printf("warpError in VelocityMesh::warpMove: Deleted GID %u still found in map with LID %u! (used to be LID %u) (thread %u)\n",removeGID,removedLIDfound,targetLID,(vmesh::LocalID)b_tid);
947 assert(0);
948 }
949 #endif
950 if (b_tid==0) {
951 ltg_size--;
952 }
953 __syncthreads();
954 return true;
955 }
956 ARCH_DEV inline size_t VelocityMesh::warpCount(const vmesh::GlobalID globalID, const size_t b_tid) const {
957 vmesh::LocalID retval = invalidLocalID();
958 globalToLocalMap.warpFind(globalID, retval, b_tid % GPUTHREADS);
959 #ifdef DEBUG_VMESH
960 __syncthreads();
961 vmesh::LocalID verify = getLocalID(globalID);
962 if (verify != retval) {
963 printf("warpError in VelocityMesh::warpCount: Searched GID %u found in map with LID %u for thread %u, but warpFind returned %u!\n",globalID,verify,(vmesh::LocalID)b_tid,retval);
964 assert(0);
965 }
966 #endif
967 __syncthreads();
968 if (retval == invalidLocalID()) {
969 return 0;
970 } else {
971 return 1;
972 }
973 }
974 ARCH_DEV inline vmesh::GlobalID VelocityMesh::warpFindBlock(vmesh::GlobalID cellIndices[3], const size_t b_tid) const {
975 // Calculate i/j/k indices of the block that would own the cell:
976 vmesh::GlobalID i_block = cellIndices[0] / (*(vmesh::getMeshWrapper()->velocityMeshes))[meshID].blockLength[0];
977 vmesh::GlobalID j_block = cellIndices[1] / (*(vmesh::getMeshWrapper()->velocityMeshes))[meshID].blockLength[1];
978 vmesh::GlobalID k_block = cellIndices[2] / (*(vmesh::getMeshWrapper()->velocityMeshes))[meshID].blockLength[2];
979
980 // Calculate block global ID:
981 vmesh::GlobalID blockGID = getGlobalID(i_block,j_block,k_block);
982
983 // If the block exists, return it:
984 vmesh::LocalID retval = invalidLocalID();
985 globalToLocalMap.warpFind(blockGID, retval, b_tid % GPUTHREADS);
986
987 #ifdef DEBUG_VMESH
988 __syncthreads();
989 auto it = globalToLocalMap.device_find(blockGID);
990 if (it == globalToLocalMap.device_end()) {
991 printf("Warp error in VelocityMesh::warpFindBlock: single-thread %u search did not find entry, warp search found LID %u for GID %u\n",(vmesh::LocalID)b_tid,retval,blockGID);
992 } else if (retval != it->second) {
993 printf("Warp error in VelocityMesh::warpFindBlock: LID %u (warp) != %u (thread %u) for GID %u\n",
994 retval,it->second,(vmesh::LocalID)b_tid,blockGID);
995 }
996 #endif
997
998 __syncthreads();
999 if (retval == invalidLocalID()) {
1000 return invalidGlobalID();
1001 } else {
1002 return blockGID;
1003 }
1004 }
1005 ARCH_DEV inline bool VelocityMesh::warpPush_back(const vmesh::GlobalID globalID, const size_t b_tid) {
1006 __shared__ bool inserted;
1007 const vmesh::LocalID mySize = size();
1008 #ifdef DEBUG_VMESH
1009 if (mySize >= (*(vmesh::getMeshWrapper()->velocityMeshes))[meshID].max_velocity_blocks) {
1010 return false;
1011 }
1012 if (globalID == invalidGlobalID()) {
1013 return false;
1014 }
1015 const vmesh::LocalID mapSize = globalToLocalMap.size();
1016 #endif
1017
1018 __syncthreads();
1019 if (b_tid < GPUTHREADS) {
1020 // If exists, do not overwrite
1021 inserted = globalToLocalMap.warpInsert_V<true>(globalID,(vmesh::LocalID)mySize, b_tid);
1022 if (inserted == true && b_tid==0) {
1023 localToGlobalMap.device_push_back(globalID);
1024 ltg_size++; // Note: called from inside kernel, cached size must be updated separately
1025 }
1026 }
1027 #ifdef DEBUG_VMESH
1028 __syncthreads();
1029 const vmesh::LocalID newMapSize = globalToLocalMap.size();
1030 const vmesh::LocalID newVecSize = size();
1031 const vmesh::LocalID postLID = getLocalID(globalID);
1032 const vmesh::GlobalID postGID = getGlobalID(mySize);
1033 if (newMapSize != mapSize+1) {
1034 printf("warpError in VelocityMesh::warpPush_back: map size %u is not expected %u! (thread %u)\n",newMapSize,mapSize+1,(vmesh::LocalID)b_tid);
1035 assert(0);
1036 }
1037 if (newVecSize != mySize+1) {
1038 printf("warpError in VelocityMesh::warpPush_back: vector size %u is not expected %u! (thread %u)\n",newVecSize,mySize+1,(vmesh::LocalID)b_tid);
1039 assert(0);
1040 }
1041 if (postLID != mySize) {
1042 printf("warpError in VelocityMesh::warpPush_back: hashmap returns LID %u but expected %u! (thread %u)\n",postLID,mySize,(vmesh::LocalID)b_tid);
1043 assert(0);
1044 }
1045 if (postLID != mySize) {
1046 printf("warpError in VelocityMesh::warpPush_back: vector entry is GID %u but expected %u! (thread %u)\n",postGID,globalID,(vmesh::LocalID)b_tid);
1047 assert(0);
1048 }
1049 #endif
1050 __syncthreads();
1051 return inserted;
1052 }
1053 ARCH_DEV inline vmesh::LocalID VelocityMesh::warpPush_back(const split::SplitVector<vmesh::GlobalID>& blocks, const size_t b_tid) {
1054 //GPUTODO: ADD debugs
1055 const vmesh::LocalID mySize = size();
1056 const vmesh::LocalID blocksSize = blocks.size();
1057 if (mySize+blocksSize > (*(vmesh::getMeshWrapper()->velocityMeshes))[meshID].max_velocity_blocks) {
1058 if (b_tid==0) {
1059 printf("vmesh: too many blocks, current size is %u",mySize);
1060 printf(", adding %u blocks", blocksSize);
1061 printf(", max is %u\n",(*(vmesh::getMeshWrapper()->velocityMeshes))[meshID].max_velocity_blocks);
1062 }
1063 return false;
1064 }
1065 __shared__ vmesh::LocalID nInserted;
1066 __shared__ bool inserted;
1067 if (b_tid==0) {
1068 inserted = true;
1069 }
1070 __syncthreads();
1071 if (b_tid < GPUTHREADS) {
1072 for (vmesh::LocalID b=0; b<blocksSize; ++b) { // GPUTODO parallelize
1073 // Do not overwrite
1074 inserted = inserted && globalToLocalMap.warpInsert_V<true>(blocks[b],(vmesh::LocalID)(mySize+b), b_tid);
1075 if (inserted == true && b_tid == 0) {
1076 nInserted = b;
1077 }
1078 }
1079 }
1080 __syncthreads();
1081 if (inserted == false) {
1082 if (b_tid == 0) {
1083 printf("vmesh: failed to push_back new block! %d of %d\n",(uint)(nInserted+1),(uint)blocksSize);
1084 }
1085 return nInserted+1;
1086 }
1087 if (b_tid==0) {
1088 localToGlobalMap.device_insert(localToGlobalMap.end(),blocks.begin(),blocks.end());
1089 }
1090 if (b_tid==0) {
1091 ltg_size += nInserted; // Note: called from inside kernel, cached size must be updated separately
1092 }
1093 __syncthreads();
1094 return blocksSize;
1095 }
1096 ARCH_DEV inline void VelocityMesh::warpReplaceBlock(const vmesh::GlobalID GIDold,
1097 const vmesh::LocalID LID,
1098 const vmesh::GlobalID GIDnew,
1099 const size_t b_tid) {
1100 // Inserts a (possibly new) block into the vmesh at a given position, and removes the existing block from there.
1101 #ifdef DEBUG_VMESH
1102 if (b_tid==0) {
1103 if (LID > size()-1) {
1104 printf("vmesh replaceBlock error: LID %u is too large for size %u!\n",LID,(vmesh::LocalID)size());
1105 assert(0);
1106 }
1107 }
1108 __syncthreads();
1109 vmesh::LocalID LIDold = invalidLocalID();
1110 globalToLocalMap.warpFind(GIDold, LIDold, b_tid % GPUTHREADS);
1111 auto it = globalToLocalMap.device_find(GIDold);
1112 if (it == globalToLocalMap.device_end()) {
1113 if (LIDold != invalidLocalID()) {
1114 printf("Warp error in VelocityMesh::warpReplaceBlock: thread %u search did not find entry, warp search found LID %u for GID %u\n",(vmesh::LocalID)b_tid,LIDold,GIDold);
1115 assert(0);
1116 }
1117 } else if (LIDold != it->second) {
1118 printf("Warp error in VelocityMesh::warpReplaceBlock: LID %u (warp) != %u (thread %u) for GID %u\n",
1119 LIDold,it->second,(vmesh::LocalID)b_tid,GIDold);
1120 assert(0);
1121 }
1122 if (b_tid==0) {
1123 if (localToGlobalMap.at(LIDold) != GIDold) {
1124 printf("vmesh replaceBlock error: oldGID and oldLID don't match!\n");
1125 assert(0);
1126 }
1127 }
1128 __syncthreads();
1129 const vmesh::LocalID preVecSize = size();
1130 const vmesh::LocalID preMapSize = globalToLocalMap.size();
1131 #endif
1132 if (b_tid < GPUTHREADS) { // GPUTODO these in parallel?
1133 globalToLocalMap.warpErase(GIDold, b_tid);
1134
1135 #ifdef DEBUG_VMESH
1136 // if (globalToLocalMap.size() != preMapSize-1) {
1137 // printf("Warp error in VelocityMesh::warpReplaceBlock: map size %u does not match expected %u for thread %u!\n",(vmesh::LocalID)globalToLocalMap.size(),(vmesh::LocalID)(preMapSize-1),(vmesh::LocalID)b_tid);
1138 // assert(0);
1139 // }
1140 auto it2 = globalToLocalMap.device_find(GIDold);
1141 if (it2 != globalToLocalMap.device_end()) {
1142 printf("Warp error in VelocityMesh::warpReplaceBlock: warp-erased GID %u LID %u but thread %u still finds LID %u associated with it!\n",GIDold,LID,(vmesh::LocalID)b_tid,it2->second);
1143 assert(0);
1144 }
1145 bool newlyadded = false;
1146 // Do not overwrite
1147 newlyadded = globalToLocalMap.warpInsert_V(GIDnew,LID, b_tid);
1148 vmesh::LocalID postMapSize = globalToLocalMap.size();
1149 // int change = 0;
1150 // if (!newlyadded) {
1151 // change = -1;
1152 // }
1153 // if (postMapSize != preMapSize + change) {
1154 // printf("Warp error in VelocityMesh::warpReplaceBlock: map size %u does not match expected %u for thread %u!\n",postMapSize,(vmesh::LocalID)(preMapSize+change),(vmesh::LocalID)b_tid);
1155 // if (b_tid==0) {
1156 // globalToLocalMap.stats();
1157 // }
1158 // assert(0);
1159 // }
1160 auto it3 = globalToLocalMap.device_find(GIDnew);
1161 if (it3 == globalToLocalMap.device_end()) {
1162 int sizePower = globalToLocalMap.getSizePower();
1163 printf("Warp error in VelocityMesh::warpReplaceBlock: warp-inserted GID %u LID %u but thread %u cannot find it!\n",GIDnew,LID,(vmesh::LocalID)b_tid);
1164 if (b_tid==0) {
1165 if (newlyadded) {
1166 printf("warpAccessor reported true for insertion!\n");
1167 } else {
1168 printf("warpAccessor reported false for insertion!\n");
1169 }
1170 globalToLocalMap.stats();
1171 //globalToLocalMap.dump_buckets();
1172 }
1173 assert(0);
1174 } else if (it3->second != LID) {
1175 printf("Warp error in VelocityMesh::warpReplaceBlock: warp-inserted GID %u LID %u but thread %u instead finds LID %u!\n",GIDnew,LID,(vmesh::LocalID)b_tid,it3->second);
1176 if (b_tid==0) {
1177 globalToLocalMap.stats();
1178 }
1179 assert(0);
1180 }
1181 #else
1182 bool newlyadded = globalToLocalMap.warpInsert_V(GIDnew,LID, b_tid);
1183 //globalToLocalMap.warpInsert(GIDnew,LID,b_tid);
1184 #endif
1185 }
1186 //__syncthreads(); // not needed
1187 if (b_tid==0) {
1188 localToGlobalMap.at(LID) = GIDnew;
1189 }
1190 __syncthreads();
1191 }
1192 ARCH_DEV inline void VelocityMesh::warpPlaceBlock(const vmesh::GlobalID GID,const vmesh::LocalID LID, const size_t b_tid) {
1193 // Places block GID into the mesh with value LID. Assumes localToGlobalMap has already been grown sufficiently.
1194
1195 #ifdef DEBUG_VMESH
1196 if (b_tid==0) {
1197 if (LID > size()-1) printf("vmesh placeBlock error: LID is too large!\n");
1198 }
1199 __syncthreads();
1200 {
1201 auto it = globalToLocalMap.device_find(GID);
1202 if (it != globalToLocalMap.device_end()) {
1203 printf("Warp error in VelocityMesh::warpPlaceBlock: single-thread %u search found GID %u=%u LID %u before it was inserted!\n",(vmesh::LocalID)b_tid,GID,it->first,it->second);
1204 if (b_tid==0) {
1205 globalToLocalMap.stats();
1206 }
1207 __syncthreads();
1208 assert(0);
1209 }
1210 }
1211 __syncthreads();
1212 #endif
1213 if (b_tid < GPUTHREADS) {
1214 bool newlyadded = false;
1215 newlyadded = globalToLocalMap.warpInsert_V(GID,LID, b_tid);
1216 if (!newlyadded) {
1217 if (b_tid==0) {
1218 globalToLocalMap.stats();
1219 printf("warpPlaceBlock error GID %u LID %u reported as not newly added! Size %zu.\n",GID,LID,localToGlobalMap.size());
1220 //globalToLocalMap.dump_buckets();
1221 }
1222 assert(newlyadded && "newlyAdded warpPlaceBlock");
1223 }
1224 localToGlobalMap.at(LID) = GID;
1225 }
1226 #ifdef DEBUG_VMESH
1227 __syncthreads();
1228 // Note: no size check possible.
1229 auto it = globalToLocalMap.device_find(GID);
1230 if (it == globalToLocalMap.device_end()) {
1231 printf("Warp error in VelocityMesh::warpPlaceBlock: single-thread %u search did not find inserted GID %u LID %u\n",(vmesh::LocalID)b_tid,GID,LID);
1232 if (b_tid==0) {
1233 globalToLocalMap.stats();
1234 }
1235 assert(0);
1236 } else if (LID != it->second) {
1237 printf("Warp error in VelocityMesh::warpPlaceBlock: LID %u (warp) != %u (thread %u) for GID %u\n",LID,it->second,(vmesh::LocalID)b_tid,GID);
1238 if (b_tid==0) {
1239 globalToLocalMap.stats();
1240 }
1241 assert(0);
1242 }
1243 #endif
1244 __syncthreads();
1245 }
1246 ARCH_DEV inline void VelocityMesh::warpDeleteBlock(const vmesh::GlobalID GID,const vmesh::LocalID LID, const size_t b_tid) {
1247 #ifdef DEBUG_VMESH
1248 // Verify that GID and LID match
1249 if (b_tid==0) {
1250 if (GID==invalidGlobalID()) {
1251 printf("vmesh warpDeleteBlock error: GID is invalidGlobalID! (LID %ul)\n",LID);
1252 }
1253 if (LID==invalidLocalID()) {
1254 printf("vmesh warpDeleteBlock error: LID is invalidLocalID! (GID %ul)\n",GID);
1255 }
1256 }
1257 __syncthreads();
1258 vmesh::LocalID retval = invalidLocalID();
1259 globalToLocalMap.warpFind(GID, retval, b_tid % GPUTHREADS);
1260 if (b_tid==0) {
1261 if (retval == invalidLocalID()) {
1262 printf("vmesh warpDeleteBlock error: GID %ul does not exist! (LID %ul)\n",GID,LID);
1263 } else {
1264 if (retval != LID) {
1265 printf("vmesh warpDeleteBlock error: LID %ul warpFound with GID %ul does not match provided LID %ul!\n",
1266 retval,GID,LID);
1267 }
1268 }
1269 if (localToGlobalMap.at(LID) != GID) {
1270 printf("vmesh warpDeleteBlock error: GID %ul warpFound with LID %ul does not match provided GID %ul!\n",
1271 localToGlobalMap.at(LID),LID,GID);
1272 }
1273 }
1274 const vmesh::LocalID preMapSize = globalToLocalMap.size();
1275 __syncthreads();
1276 #endif
1277 if (b_tid < GPUTHREADS) {
1278 globalToLocalMap.warpErase(GID, b_tid);
1279 if (b_tid==0) {
1280 localToGlobalMap.at(LID) = invalidGlobalID();
1281 }
1282 }
1283 __syncthreads();
1284 #ifdef DEBUG_VMESH
1285 // const vmesh::LocalID postMapSize = globalToLocalMap.size();
1286 // if (postMapSize != preMapSize-1) {
1287 // printf("Warp error in VelocityMesh::warpDeleteBlock: map size %u does not match expected %u for thread %u!\n",postMapSize,preMapSize-1,(vmesh::LocalID)b_tid);
1288 // assert(0);
1289 // }
1290 auto it = globalToLocalMap.device_find(GID);
1291 if (it != globalToLocalMap.device_end()) {
1292 printf("Warp error in VelocityMesh::warpDeleteBlock: GID %u still found with LID %u for thread %u!\n",GID,it->second,(vmesh::LocalID)b_tid);
1293 assert(0);
1294 }
1295 __syncthreads();
1296 #endif
1297 }
1298#endif
1299
1300 inline void VelocityMesh::setGrid() {
1301 // Assumes we have a valid localToGlobalMap from e.g. MPI communication,
1302 // populates globalToLocalMap based on it.
1303 gpuStream_t stream = gpu_getStream();
1304 globalToLocalMap.clear<false>(Hashinator::targets::device,stream,std::pow(2,gtl_sizepower));
1305 CHK_ERR( gpuStreamSynchronize(stream) );
1306 size_t nBlocks = localToGlobalMap.size();
1307 globalToLocalMap.insertIndex<false>(localToGlobalMap.data(),nBlocks,0.5,stream);
1308 CHK_ERR( gpuStreamSynchronize(stream) );
1309 ltg_size = nBlocks;
1310 ltg_capacity = localToGlobalMap.capacity();
1311 gtl_sizepower = globalToLocalMap.getSizePower();
1312 }
1313
1314 // GPUTODO: These accessors are still slow, but we don't actually use them at all.
1315 inline bool VelocityMesh::setGrid(const std::vector<vmesh::GlobalID>& globalIDs) {
1316 printf("Warning! Slow version of VelocityMesh::setGrid.\n");
1317 gpuStream_t stream = gpu_getStream();
1318 globalToLocalMap.clear<false>(Hashinator::targets::device,stream,std::pow(2,gtl_sizepower));
1319 CHK_ERR( gpuStreamSynchronize(stream) );
1320 for (vmesh::LocalID i=0; i<globalIDs.size(); ++i) {
1321 globalToLocalMap.insert(Hashinator::make_pair(globalIDs[i],(vmesh::LocalID)i));
1322 }
1323 localToGlobalMap.clear();
1324 localToGlobalMap.insert(localToGlobalMap.end(),globalIDs.begin(),globalIDs.end());
1325 ltg_size = globalIDs.size();
1326 ltg_capacity = localToGlobalMap.capacity();
1327 gtl_sizepower = globalToLocalMap.getSizePower();
1328 return true;
1329 }
1330 inline bool VelocityMesh::setGrid(const split::SplitVector<vmesh::GlobalID>& globalIDs) {
1331 printf("Warning! Slow version of VelocityMesh::setGrid.\n");
1332 gpuStream_t stream = gpu_getStream();
1333 globalToLocalMap.clear<false>(Hashinator::targets::device,stream,std::pow(2,gtl_sizepower));
1334 CHK_ERR( gpuStreamSynchronize(stream) );
1335 for (vmesh::LocalID i=0; i<globalIDs.size(); ++i) {
1336 globalToLocalMap.insert(Hashinator::make_pair(globalIDs[i],(vmesh::LocalID)i));
1337 }
1338 localToGlobalMap.clear();
1339 localToGlobalMap.insert(localToGlobalMap.end(),globalIDs.begin(),globalIDs.end());
1340 ltg_size = globalIDs.size();
1341 ltg_capacity = localToGlobalMap.capacity();
1342 gtl_sizepower = globalToLocalMap.getSizePower();
1343 return true;
1344 }
1345
1346 inline bool VelocityMesh::setMesh(const size_t meshID) {
1347 if (meshID >= vmesh::getMeshWrapper()->velocityMeshes->size()) {
1348 return false;
1349 }
1350 this->meshID = meshID;
1351 return true;
1352 }
1353 inline size_t VelocityMesh::getMeshID() {
1354 return this->meshID;
1355 }
1356
1357 inline void VelocityMesh::setNewSize(const vmesh::LocalID newSize) {
1358 // Needed by GPU block adjustment
1359 setNewCapacity(newSize);
1360 gpuStream_t stream = gpu_getStream();
1361 localToGlobalMap.resize(newSize,true,stream);
1362 ltg_size = newSize;
1363 }
1364
1365 ARCH_DEV inline void VelocityMesh::device_setNewSize(const vmesh::LocalID newSize) {
1366 if (newSize>0) {
1367 const vmesh::LocalID currentCapacity = localToGlobalMap.capacity();
1368 assert(newSize <= currentCapacity && "insufficient vector capacity in vmesh::device_setNewSize");
1369 const int currentSizePower = globalToLocalMap.getSizePower();
1370 assert(ceil(log2(newSize)) <= currentSizePower && "insufficient map capacity in vmesh::device_setNewSize");
1371 }
1372 localToGlobalMap.device_resize(newSize,false); //construct=false don't construct or set to zero
1373 ltg_size = newSize; // Remember to update size on host as well
1374 }
1375 inline void VelocityMesh::setNewCachedSize(const vmesh::LocalID newSize) {
1376 // Should only be used to update host-side size if resizing on device
1377 ltg_size = newSize;
1378 }
1379 inline void VelocityMesh::updateCachedSize() {
1380 // More secure page-faulting way to update cached size
1381 ltg_size = localToGlobalMap.size();
1382 }
1383 inline void VelocityMesh::updateCachedCapacity() {
1384 // Should not be needed, added as an optional safeguard
1385 ltg_capacity = localToGlobalMap.capacity();
1386 gtl_sizepower = globalToLocalMap.getSizePower();
1387 }
1388
1389 // Used in initialization
1390 inline bool VelocityMesh::setNewCapacity(const vmesh::LocalID newCapacity) {
1391 // Ensure also that the map is large enough (newCapacity always greater than zero here)
1392 uint HashmapReqSize = (uint)ceil(log2(newCapacity)) +1;
1393
1394 if (ltg_capacity >= newCapacity && gtl_sizepower >= HashmapReqSize) {
1395 return false; // Still have enough buffer
1396 }
1397 gpuStream_t stream = gpu_getStream();
1398 ltg_capacity = newCapacity*BLOCK_ALLOCATION_FACTOR;
1399 // Passing eco flag = true to resize tells splitvector we manage padding manually.
1400 localToGlobalMap.reserve(ltg_capacity,true, stream);
1401 if (gtl_sizepower < HashmapReqSize) {
1402 gtl_sizepower = HashmapReqSize;
1403 globalToLocalMap.resize(gtl_sizepower, Hashinator::targets::device, stream);
1404 }
1405 CHK_ERR( gpuStreamSynchronize(stream) );
1406 return true;
1407 }
1408
1409 ARCH_HOSTDEV inline size_t VelocityMesh::size() const {
1410 #if defined(__CUDA_ARCH__) || defined(__HIP_DEVICE_COMPILE__)
1411 return localToGlobalMap.size();
1412 #else
1413 #ifdef DEBUG_VMESH
1414 const size_t size1 = localToGlobalMap.size();
1415 if (ltg_size != size1) {
1416 printf("VMESH SIZE ERROR: size %lu vs cached value %lu in %s : %d\n",size1,ltg_size,__FILE__,__LINE__);
1417 }
1418 #endif
1419 return ltg_size;
1420 #endif
1421 }
1422
1423 ARCH_HOSTDEV inline size_t VelocityMesh::sizeInBytes() const {
1424 #ifdef DEBUG_VMESH
1425 const size_t size1 = localToGlobalMap.size();
1426 const size_t size2 = globalToLocalMap.size();
1427 if (ltg_size != size1) {
1428 printf("VMESH SIZE ERROR: LTG size %lu vs cached value %lu in %s : %d\n",size1,ltg_size,__FILE__,__LINE__);
1429 }
1430 if (ltg_size != size2) {
1431 printf("VMESH SIZE ERROR: GTL size %lu vs cached value %lu in %s : %d\n",size1,ltg_size,__FILE__,__LINE__);
1432 }
1433 #endif
1434 const size_t sizeInBytes = ltg_size * sizeof(vmesh::GlobalID)
1435 + ltg_size * sizeof(Hashinator::hash_pair<vmesh::GlobalID,vmesh::LocalID>);
1436 return sizeInBytes;
1437 }
1438
1439 // inline void VelocityMesh::swap(VelocityMesh& vm) {
1440 // gpuStream_t stream = gpu_getStream();
1441 // globalToLocalMap.swap(vm.globalToLocalMap);
1442 // localToGlobalMap.swap(vm.localToGlobalMap);
1443 // }
1444 inline void VelocityMesh::gpu_prefetchHost(gpuStream_t stream=0) {
1445 if (stream==0) {
1446 stream = gpu_getStream();
1447 }
1448 localToGlobalMap.optimizeCPU(stream);
1449 globalToLocalMap.optimizeCPU(stream);
1450 return;
1451 }
1452
1453 inline void VelocityMesh::gpu_prefetchDevice(gpuStream_t stream=0) {
1454 if (stream==0) {
1455 stream = gpu_getStream();
1456 }
1457 //phiprof::Timer vmeshPrefetchTimer {"prefetch Vmesh"};
1458 localToGlobalMap.optimizeGPU(stream);
1459 globalToLocalMap.optimizeGPU(stream);
1460 //CHK_ERR( gpuStreamSynchronize(stream) );
1461 return;
1462 }
1463
1464 inline void VelocityMesh::gpu_cleanHashMap(gpuStream_t stream = 0) {
1465 if (stream==0) {
1466 stream = gpu_getStream();
1467 }
1468
1469 // phiprof::Timer resizeTimer {"Hashinator resize"};
1470 // //globalToLocalMap.performCleanupTasks<false>(stream);
1471 // globalToLocalMap.resize_to_lf(0.5, Hashinator::targets::device, stream);
1472 // CHK_ERR( gpuStreamSynchronize(stream) );
1473 // resizeTimer.stop();
1474
1475 phiprof::Timer cleanupTimer {"Hashinator tombstones"};
1476 // globalToLocalMap.clean_tombstones<false>(stream);
1477 globalToLocalMap.performCleanupTasks<false>(stream);
1478 CHK_ERR( gpuStreamSynchronize(stream) );
1479 cleanupTimer.stop();
1480 return;
1481 }
1482
1483 ARCH_HOSTDEV inline Hashinator::Hashmap<vmesh::GlobalID,vmesh::LocalID>* VelocityMesh::gpu_expose_map() {
1484 return &globalToLocalMap;
1485 }
1486
1487 inline void VelocityMesh::print_addresses() {
1488 printf("GPU localToGlobalMap %p\n GPU globalToLocalMap %p\n",&localToGlobalMap,&globalToLocalMap);
1489 printf("GPU localToGlobalMap capacity %zu size %zu \n GPU globalToLocalMap size %zu bucket count %zu\n",localToGlobalMap.capacity(),localToGlobalMap.size(),globalToLocalMap.size(),globalToLocalMap.bucket_count());
1490 printf("GPU localToGlobalMap data %p\n",localToGlobalMap.data());
1491 //printf("GPU localToGlobalMap iterators begin %p end %p\n",localToGlobalMap.begin(),localToGlobalMap.end());
1492 }
1493 inline void VelocityMesh::print_sizes() {
1494 printf("GPU localToGlobalMap size %lu capacity %lu cached size %lu cached capacity %lu\nGPU globalToLocalMap fill %lu sizePower %lu cached sizePower %lu\n",
1495 localToGlobalMap.size(),localToGlobalMap.capacity(),ltg_size,ltg_capacity,globalToLocalMap.size(),(size_t)globalToLocalMap.getSizePower(),gtl_sizepower);
1496 }
1497
1498} // namespace vmesh
1499
1500#endif
i
for i
Definition Dispersion.m:24
ARCH_DEV
#define ARCH_DEV
Definition arch_device_api.h:11
ARCH_HOSTDEV
#define ARCH_HOSTDEV
Definition arch_device_api.h:10
gpuStream_t
#define gpuStream_t
Definition arch_device_cuda.h:77
gpuStreamSynchronize
#define gpuStreamSynchronize
Definition arch_device_cuda.h:59
CHK_ERR
#define CHK_ERR(err)
Definition arch_device_cuda.h:123
gpuDeviceSynchronize
#define gpuDeviceSynchronize
Definition arch_device_cuda.h:25
GPUTHREADS
#define GPUTHREADS
Definition arch_device_cuda.h:113
OpenBucketHashtable::count
size_t count(const GID &key) const
Definition open_bucket_hashtable.h:163
vmesh::VelocityMesh::warpCount
ARCH_DEV size_t warpCount(const vmesh::GlobalID globalID, const size_t b_tid) const
vmesh::VelocityMesh::setNewCapacity
void setNewCapacity(const vmesh::LocalID &newCapacity)
Definition velocity_mesh_cpu.h:493
vmesh::VelocityMesh::getIndicesZ
ARCH_HOSTDEV void getIndicesZ(const vmesh::GlobalID globalID, vmesh::LocalID &k) const
Definition velocity_mesh_gpu.h:549
vmesh::VelocityMesh::capacity
ARCH_HOSTDEV size_t capacity() const
Definition velocity_mesh_gpu.h:211
vmesh::VelocityMesh::gpu_cleanHashMap
void gpu_cleanHashMap(gpuStream_t stream)
Definition velocity_mesh_gpu.h:1464
vmesh::VelocityMesh::getGridLength
ARCH_HOSTDEV const vmesh::LocalID * getGridLength() const
vmesh::VelocityMesh::push_back
bool push_back(const vmesh::GlobalID &globalID)
Definition velocity_mesh_cpu.h:429
vmesh::VelocityMesh::gpu_expose_map
ARCH_HOSTDEV Hashinator::Hashmap< vmesh::GlobalID, vmesh::LocalID > * gpu_expose_map()
Definition velocity_mesh_gpu.h:1483
vmesh::VelocityMesh::setNewCachedSize
void setNewCachedSize(const vmesh::LocalID newSize)
Definition velocity_mesh_gpu.h:1375
vmesh::VelocityMesh::setGrid
bool setGrid(const std::vector< vmesh::GlobalID > &globalIDs)
vmesh::VelocityMesh::operator=
const VelocityMesh & operator=(const VelocityMesh &other)
vmesh::VelocityMesh::pop
ARCH_HOSTDEV void pop()
vmesh::VelocityMesh::globalToLocalMap
OpenBucketHashtable< vmesh::GlobalID, vmesh::LocalID > globalToLocalMap
Definition velocity_mesh_cpu.h:102
vmesh::VelocityMesh::placeBlock
ARCH_DEV void placeBlock(const vmesh::GlobalID GID, const vmesh::LocalID LID)
vmesh::VelocityMesh::warpMove
ARCH_DEV bool warpMove(const vmesh::LocalID sourceLocalID, const vmesh::LocalID targetLocalID, const size_t b_tid)
vmesh::VelocityMesh::invalidBlockIndex
static vmesh::LocalID invalidBlockIndex()
Definition velocity_mesh_cpu.h:400
vmesh::VelocityMesh::warpPush_back
ARCH_DEV bool warpPush_back(const vmesh::GlobalID globalID, const size_t b_tid)
vmesh::VelocityMesh::getGlobalID
ARCH_HOSTDEV vmesh::GlobalID getGlobalID(const Real *coords) const
vmesh::VelocityMesh::size
ARCH_HOSTDEV size_t size() const
Definition velocity_mesh_gpu.h:1409
vmesh::VelocityMesh::getBlockCoordinates
bool getBlockCoordinates(const vmesh::GlobalID &globalID, Real coords[3]) const
Definition velocity_mesh_cpu.h:234
vmesh::VelocityMesh::getGrid
ARCH_HOSTDEV split::SplitVector< vmesh::GlobalID > * getGrid()
vmesh::VelocityMesh::warpPop
ARCH_DEV void warpPop(const size_t b_tid)
vmesh::VelocityMesh::warpPush_back
ARCH_DEV vmesh::LocalID warpPush_back(const split::SplitVector< vmesh::GlobalID > &blocks, const size_t b_tid)
vmesh::VelocityMesh::invalidGlobalID
static ARCH_HOSTDEV vmesh::GlobalID invalidGlobalID()
vmesh::VelocityMesh::push_back
vmesh::LocalID push_back(const std::vector< vmesh::GlobalID > &blocks)
vmesh::VelocityMesh::invalidLocalID
static vmesh::LocalID invalidLocalID()
Definition velocity_mesh_cpu.h:408
vmesh::VelocityMesh::getMeshMaxLimits
ARCH_HOSTDEV const Real * getMeshMaxLimits() const
vmesh::VelocityMesh::device_setNewSize
ARCH_DEV void device_setNewSize(const vmesh::LocalID newSize)
Definition velocity_mesh_gpu.h:1365
vmesh::VelocityMesh::move
bool move(const vmesh::LocalID &sourceLocalID, const vmesh::LocalID &targetLocalID)
Definition velocity_mesh_cpu.h:184
vmesh::VelocityMesh::getBlockSize
ARCH_HOSTDEV const Real * getBlockSize() const
vmesh::VelocityMesh::setMesh
bool setMesh(const size_t &meshID)
Definition velocity_mesh_cpu.h:480
vmesh::VelocityMesh::setNewSize
void setNewSize(const vmesh::LocalID &newSize)
Definition velocity_mesh_cpu.h:488
vmesh::VelocityMesh::count
size_t count(const vmesh::GlobalID &globalID) const
Definition velocity_mesh_cpu.h:200
vmesh::VelocityMesh::getMeshMinLimits
ARCH_HOSTDEV const Real * getMeshMinLimits() const
vmesh::VelocityMesh::print
ARCH_HOSTDEV void print()
Definition velocity_mesh_gpu.h:293
vmesh::VelocityMesh::getLocalID
vmesh::LocalID getLocalID(const vmesh::GlobalID &globalID) const
Definition velocity_mesh_cpu.h:371
vmesh::VelocityMesh::getGlobalID
ARCH_HOSTDEV vmesh::GlobalID getGlobalID(const vmesh::LocalID indices[3]) const
vmesh::VelocityMesh::deleteBlock
ARCH_DEV void deleteBlock(const vmesh::GlobalID GID, const vmesh::LocalID LID)
vmesh::VelocityMesh::clear
void clear(bool shrink=false)
Definition velocity_mesh_cpu.h:169
vmesh::VelocityMesh::findBlock
ARCH_HOSTDEV vmesh::GlobalID findBlock(vmesh::GlobalID cellIndices[3]) const
vmesh::VelocityMesh::warpDeleteBlock
ARCH_DEV void warpDeleteBlock(const vmesh::GlobalID GID, const vmesh::LocalID LID, const size_t b_tid)
vmesh::VelocityMesh::sizeInBytes
ARCH_HOSTDEV size_t sizeInBytes() const
vmesh::VelocityMesh::warpFindBlock
ARCH_DEV vmesh::GlobalID warpFindBlock(vmesh::GlobalID cellIndices[3], const size_t b_tid) const
vmesh::VelocityMesh::VelocityMesh
VelocityMesh(const VelocityMesh &other)
vmesh::VelocityMesh::getIndicesX
ARCH_HOSTDEV void getIndicesX(const vmesh::GlobalID globalID, vmesh::LocalID &i) const
Definition velocity_mesh_gpu.h:535
vmesh::VelocityMesh::replaceBlock
ARCH_DEV void replaceBlock(const vmesh::GlobalID GIDold, const vmesh::LocalID LID, const vmesh::GlobalID GIDnew)
vmesh::VelocityMesh::getCellSize
ARCH_HOSTDEV const Real * getCellSize() const
vmesh::VelocityMesh::warpReplaceBlock
ARCH_DEV void warpReplaceBlock(const vmesh::GlobalID GIDold, const vmesh::LocalID LID, const vmesh::GlobalID GIDnew, const size_t b_tid)
vmesh::VelocityMesh::gpu_prefetchDevice
void gpu_prefetchDevice(gpuStream_t stream)
Definition velocity_mesh_gpu.h:1453
vmesh::VelocityMesh::isInitialized
ARCH_HOSTDEV bool isInitialized() const
vmesh::VelocityMesh::localToGlobalMap
std::vector< vmesh::GlobalID > localToGlobalMap
Definition velocity_mesh_cpu.h:101
vmesh::VelocityMesh::getGlobalID
vmesh::GlobalID getGlobalID(const vmesh::LocalID &localID) const
Definition velocity_mesh_cpu.h:300
vmesh::VelocityMesh::clear
void clear(bool shrink)
vmesh::VelocityMesh::warpPlaceBlock
ARCH_DEV void warpPlaceBlock(const vmesh::GlobalID GID, const vmesh::LocalID LID, const size_t b_tid)
vmesh::VelocityMesh::invalidGlobalID
static vmesh::GlobalID invalidGlobalID()
Definition velocity_mesh_cpu.h:404
vmesh::VelocityMesh::size
size_t size(bool dummy=0) const
Definition velocity_mesh_cpu.h:497
vmesh::VelocityMesh::initialize
bool initialize(const size_t &meshID)
Definition velocity_mesh_cpu.h:395
vmesh::VelocityMesh::invalidBlockIndex
static ARCH_HOSTDEV vmesh::LocalID invalidBlockIndex()
vmesh::VelocityMesh::gpu_prefetchHost
void gpu_prefetchHost(gpuStream_t stream)
Definition velocity_mesh_gpu.h:1444
vmesh::VelocityMesh::getIndices
void getIndices(const vmesh::GlobalID &globalID, vmesh::LocalID &i, vmesh::LocalID &j, vmesh::LocalID &k) const
Definition velocity_mesh_cpu.h:361
vmesh::VelocityMesh::invalidLocalID
static ARCH_HOSTDEV vmesh::LocalID invalidLocalID()
vmesh::VelocityMesh::capacityInBytes
size_t capacityInBytes() const
vmesh::VelocityMesh::getMesh
ARCH_HOSTDEV size_t getMesh() const
vmesh::VelocityMesh::getBlockInfo
void getBlockInfo(const vmesh::GlobalID &globalID, Real *array) const
Definition velocity_mesh_cpu.h:253
vmesh::VelocityMesh::warpGetLocalID
ARCH_DEV vmesh::LocalID warpGetLocalID(const vmesh::GlobalID globalID, const size_t b_tid) const
vmesh::VelocityMesh::getIndicesY
ARCH_HOSTDEV void getIndicesY(const vmesh::GlobalID globalID, vmesh::LocalID &j) const
Definition velocity_mesh_gpu.h:542
Real
float Real
Definition definitions.h:41
gpu_getStream
__host__ gpuStream_t gpu_getStream()
Definition gpu_base.cpp:244
gpu_base.hpp
INIT_MAP_SIZE
static const uint INIT_MAP_SIZE(16 - WID)
INIT_VMESH_SIZE
static const uint INIT_VMESH_SIZE(32768/WID3)
BLOCK_ALLOCATION_FACTOR
static const double BLOCK_ALLOCATION_FACTOR
Definition gpu_base.hpp:61
blockSize
const int blockSize
Definition gpu_pitch_angle_diffusion.cpp:392
j
const int j
Definition gpu_pitch_angle_diffusion.cpp:60
__syncthreads
__syncthreads()
k
const int k
Definition gpu_pitch_angle_diffusion.cpp:61
faceestimates_amr.b
b
Definition faceestimates_amr.py:51
mosaic_worker.position
position
Definition mosaic_worker.py:43
vmesh::INVALID_LOCALID
static const LocalID INVALID_LOCALID
Definition definitions.h:66
vmesh::INVALID_GLOBALID
static const GlobalID INVALID_GLOBALID
Definition definitions.h:63
vmesh::LocalID
uint32_t LocalID
Definition definitions.h:60
vmesh::GlobalID
uint32_t GlobalID
Definition definitions.h:59
vmesh::getMeshWrapper
ARCH_HOSTDEV MeshWrapper * getMeshWrapper()
Definition velocity_mesh_parameters.h:120
vmesh::INVALID_VEL_BLOCK_INDEX
static const LocalID INVALID_VEL_BLOCK_INDEX
Definition definitions.h:69
vmesh::MeshWrapper::velocityMeshes
std::array< vmesh::MeshParameters, MAX_VMESH_PARAMETERS_COUNT > * velocityMeshes
Definition velocity_mesh_parameters.h:95
floor
static ARCH_HOSTDEV VecSimple< T > floor(VecSimple< T > const &a)
Definition vectorclass_fallback.h:598
velocity_mesh_parameters.h