diff --git a/libNeonDomain/include/Neon/domain/details/mGrid/mGrid.h b/libNeonDomain/include/Neon/domain/details/mGrid/mGrid.h index 8f8796c7..8c2fa310 100644 --- a/libNeonDomain/include/Neon/domain/details/mGrid/mGrid.h +++ b/libNeonDomain/include/Neon/domain/details/mGrid/mGrid.h @@ -190,12 +190,19 @@ class mGrid auto levelBitMaskIsSet(int l, const Neon::index_3d& blockID, const Neon::index_3d& localChild) const -> bool; + // check if the bitmask is set, addressing the bitmask directly by its own + // coordinate (i.e. blockID * 2 + localChild, already combined by the caller) + auto levelBitMaskIsSetAt(int l, const Neon::index_3d& bxyz) const -> bool; + // set the bitmask assuming a dense domain auto setLevelBitMask(int l, const Neon::index_3d& blockID, const Neon::index_3d& localChild) -> void; // clear the bitmask assuming a dense domain auto clearLevelBitMask(int l, const Neon::index_3d& blockID, const Neon::index_3d& localChild) -> void; + // clear the bitmask, addressing it directly by its own coordinate + auto clearLevelBitMaskAt(int l, const Neon::index_3d& bxyz) -> void; + struct Data { Neon::index_3d domainSize; diff --git a/libNeonDomain/include/Neon/domain/details/mGrid/sparseBitMask.h b/libNeonDomain/include/Neon/domain/details/mGrid/sparseBitMask.h index b059fc63..fb1a44bd 100644 --- a/libNeonDomain/include/Neon/domain/details/mGrid/sparseBitMask.h +++ b/libNeonDomain/include/Neon/domain/details/mGrid/sparseBitMask.h @@ -1,4 +1,3 @@ - #pragma once #include "Neon/domain/tools/PointHashTable.h" @@ -237,17 +236,19 @@ struct BitBlock template auto setON(Neon::index_3d const& point) -> void { - auto blockLocalPoint = point % BitBlock::blockSize; - size_t pitch = blockLocalPoint.mPitch(blockSize); - auto brickID = pitch / widthBrick; - auto localBit = pitch % widthBrick; + auto blockLocalPoint = point % BitBlock::blockSize; + size_t pitch = blockLocalPoint.mPitch(blockSize); + auto brickID = pitch / widthBrick; + auto localBit = pitch % widthBrick; + Brick const mask = Brick(1) << localBit; if constexpr (ThreadSafe) { -#pragma omp critical(settingBitOp) - { - bits[brickID] |= (1 << localBit); - } + // An atomic read-modify-write on a single brick. This replaces a named + // critical section, which was a *global* lock shared by every BitBlock in + // the process and serialized all bitmask writes during mGrid construction. +#pragma omp atomic update + bits[brickID] |= mask; } else { - bits[brickID] |= (1 << localBit); + bits[brickID] |= mask; } } @@ -292,17 +293,17 @@ struct BitBlock template auto setOFF(Neon::index_3d const& point) -> void { - auto blockLocalPoint = point % BitBlock::blockSize; - size_t pitch = blockLocalPoint.mPitch(blockSize); - auto brickID = pitch / widthBrick; - auto localBit = pitch % widthBrick; + auto blockLocalPoint = point % BitBlock::blockSize; + size_t pitch = blockLocalPoint.mPitch(blockSize); + auto brickID = pitch / widthBrick; + auto localBit = pitch % widthBrick; + Brick const mask = ~(Brick(1) << localBit); if constexpr (ThreadSafe) { -#pragma omp critical(settingBitOp) - { - bits[brickID] &= ~(1 << localBit); - } + // See setON: an atomic update instead of a process-wide critical section. +#pragma omp atomic update + bits[brickID] &= mask; } else { - bits[brickID] &= ~(1 << localBit); + bits[brickID] &= mask; } } }; @@ -421,11 +422,20 @@ struct BitBlock template class SparseBitBlocks { - Neon::domain::tool::PointHashTable mHashTable; /**< Spatial hash table mapping block coordinates to BitBlock pointers */ + Neon::domain::tool::PointHashTable mHashTable; /**< Spatial hash table mapping block coordinates to BitBlock pointers (fallback for huge coordinate spaces) */ + std::vector mDenseLookup; /**< Flat block-coordinate -> BitBlock* table; empty when the fallback hash table is in use */ Neon::index_3d mBBox; /**< Bounding box defining the valid coordinate space */ using Pool = std::array; /**< Type alias for memory pool arrays */ std::vector memoryPool; /**< Dynamic memory pool for BitBlock allocation */ size_t firstFreeIndex; /**< Index of the first free BitBlock in the current memory pool */ + std::vector mAllocatedBlocks; /**< Coordinate of every BitBlock handed out, in allocation order */ + + /** + * Above this many block slots the flat lookup table is abandoned in favour of the + * hash table. 64M slots is 512MB of pointers, which is already well past the point + * where a domain of this shape is practical. + */ + static constexpr size_t denseLookupMaxEntries = 64ull * 1024ull * 1024ull; public: /** @@ -466,9 +476,23 @@ class SparseBitBlocks * @see isActivePoint() To query point states */ SparseBitBlocks(const Neon::index_3d& bbox) - : mBBox(bbox) + : mBBox( (bbox.x + BitBlock::blockSize.x - 1) / BitBlock::blockSize.x, + (bbox.y + BitBlock::blockSize.y - 1) / BitBlock::blockSize.y, + (bbox.z + BitBlock::blockSize.z - 1) / BitBlock::blockSize.z ) { - mHashTable = Neon::domain::tool::PointHashTable(bbox); + // The block-coordinate space is small enough (one slot per 8^3 region) that a + // flat table almost always fits: it costs 8 bytes per slot and turns every + // isActivePoint() from a hash probe into a single indexed load. Only fall back + // to the hash table for coordinate spaces where that would be unreasonable. + const size_t numBlockSlots = static_cast(mBBox.x) * + static_cast(mBBox.y) * + static_cast(mBBox.z); + if (numBlockSlots <= denseLookupMaxEntries) { + mDenseLookup.assign(numBlockSlots, nullptr); + } else { + mHashTable = Neon::domain::tool::PointHashTable(mBBox); + } + auto newPoolPtr = new Pool{}; memoryPool.emplace_back(newPoolPtr); firstFreeIndex = 0; @@ -534,50 +558,57 @@ class SparseBitBlocks template auto activatePoint(const Neon::index_3d& point) -> void { - if (!(point < mBBox)) { - std::cout << "Error -> point outside of valid range" << std::endl; + Neon::index_3d block_coord(point.x / BitBlock::blockSize.x, + point.y / BitBlock::blockSize.y, + point.z / BitBlock::blockSize.z); + Neon::index_3d local(point.x % BitBlock::blockSize.x, + point.y % BitBlock::blockSize.y, + point.z % BitBlock::blockSize.z); + + if (!(block_coord < mBBox)) { + std::cout << "Error -> block outside of valid range" << std::endl; std::exit(1); } - BitBlock* bitBlock = getBitBlockPrt(point); + + BitBlock* bitBlock = getBitBlockPrt(block_coord); if (bitBlock != nullptr) { - bitBlock->setON(point); + bitBlock->setON(local); + return; } + if constexpr (ThreadSafe == true) { #pragma omp critical(SparseBitBlocks_addPoint) { - bitBlock = getBitBlockPrt(point); + bitBlock = getBitBlockPrt(block_coord); if (bitBlock == nullptr) { - if (firstFreeIndex > memoryPoolGranularity) { - std::cout << "Error -> firstFreeIndex == memoryPoolGranularity" << std::endl; - std::exit(1); - } if (firstFreeIndex == memoryPoolGranularity) { - //auto new_pool = Pool{}; - memoryPool.emplace_back({}); + auto new_pool_ptr = new Pool{}; + memoryPool.emplace_back(new_pool_ptr); firstFreeIndex = 0; } - bitBlock = &memoryPool[memoryPool.size() - 1][firstFreeIndex]; + bitBlock = &((*memoryPool.back())[firstFreeIndex]); firstFreeIndex++; - mHashTable.addPoint(point, bitBlock); + // Set the bit *before* publishing the pointer. Readers do not take + // this lock, so a block must never become reachable in a state where + // its first bit is still being written non-atomically. + bitBlock->setON(local); + publishBitBlockPtr(block_coord, bitBlock); + } else { + bitBlock->setON(local); } - bitBlock->setON(point); } } else { // We are in a critical section managed by the calling - if (firstFreeIndex > memoryPoolGranularity) { - std::cout << "Error -> firstFreeIndex == memoryPoolGranularity" << std::endl; - std::exit(1); - } if (firstFreeIndex == memoryPoolGranularity) { auto new_pool_ptr = new Pool{}; memoryPool.emplace_back(new_pool_ptr); firstFreeIndex = 0; } - bitBlock = &(memoryPool[memoryPool.size() - 1]->at(firstFreeIndex)); + bitBlock = &((*memoryPool.back())[firstFreeIndex]); firstFreeIndex++; - mHashTable.addPoint(point, bitBlock); - bitBlock->setON(point); + bitBlock->setON(local); + publishBitBlockPtr(block_coord, bitBlock); return; } } @@ -630,20 +661,109 @@ class SparseBitBlocks template auto removePoint(const Neon::index_3d& point) -> void { - BitBlock* bitBlock = getBitBlockPrt(point); + Neon::index_3d block_coord(point.x / BitBlock::blockSize.x, + point.y / BitBlock::blockSize.y, + point.z / BitBlock::blockSize.z); + Neon::index_3d local(point.x % BitBlock::blockSize.x, + point.y % BitBlock::blockSize.y, + point.z % BitBlock::blockSize.z); + + BitBlock* bitBlock = getBitBlockPrt(block_coord); if (bitBlock == nullptr) { - std::cout << "Error -> " << std::endl; + return; } - bitBlock->setOFF(point); + bitBlock->setOFF(local); } - auto getBitBlockPrt(const Neon::index_3d& point) const -> BitBlock* + /** + * @brief Makes a freshly allocated BitBlock reachable from the given block coordinate. + * + * Must only be called while holding the SparseBitBlocks_addPoint critical section. + */ + auto publishBitBlockPtr(const Neon::index_3d& coord, BitBlock* bitBlock) -> void { - if (!(point < mBBox)) { - std::cout << "Error -> point outside of valid range" << std::endl; + mAllocatedBlocks.push_back(coord); + if (!mDenseLookup.empty()) { + mDenseLookup[coord.mPitch(mBBox)] = bitBlock; + return; + } + mHashTable.addPoint(coord, bitBlock); + } + + /** + * @brief Coordinates of every BitBlock that has been allocated. + * + * That is, every 8^3 region that has held at least one active point at some point. + * A caller that needs to visit the active set can iterate these instead of sweeping + * the whole coordinate space: for a sparse domain that is orders of magnitude less + * work, and any point outside an allocated block is inactive by construction. + * + * Blocks are never freed, so a block whose bits have all since been cleared still + * appears here. Returned by value because callers commonly mutate the collection + * while iterating. + */ + auto getAllocatedBlockCoords() const -> std::vector + { + return mAllocatedBlocks; + } + + /** + * @brief Invokes f(point) for every active point inside one allocated BitBlock. + * + * The point handed to the lambda is in this collection's coordinate space, not + * block-local. Bricks holding nothing are skipped, so an allocated block that is + * mostly empty costs far less than its 512 bits suggest. + * + * Bits are read as the traversal proceeds. A concurrent update may or may not be + * observed, so a caller that mutates while iterating must re-check the points it + * cares about with isActivePoint(). + */ + template + auto forEachActivePointInBlock(const Neon::index_3d& blockCoord, + const UserLambda& f) const -> void + { + BitBlock* bitBlock = getBitBlockPrt(blockCoord); + if (bitBlock == nullptr) { + return; + } + + constexpr int edge = BitBlock::blockEdge; + const int baseX = blockCoord.x * edge; + const int baseY = blockCoord.y * edge; + const int baseZ = blockCoord.z * edge; + + for (int brickID = 0; brickID < BitBlock::numBricks; ++brickID) { + const BitBlock::Brick brick = bitBlock->bits[brickID]; + if (brick == 0) { + continue; + } + for (int bit = 0; bit < BitBlock::widthBrick; ++bit) { + if ((brick & (BitBlock::Brick(1) << bit)) == 0) { + continue; + } + // Matches Integer_3d::mPitch: pitch = x + y * edge + z * edge * edge + const int pitch = brickID * BitBlock::widthBrick + bit; + f(Neon::index_3d(baseX + pitch % edge, + baseY + (pitch / edge) % edge, + baseZ + pitch / (edge * edge))); + } + } + } + + auto getBitBlockPrt(const Neon::index_3d& coord) const -> BitBlock* + { + if (!(coord < mBBox)) { + std::cout << "Error -> coord outside of valid range" << std::endl; std::exit(1); } - BitBlock* const* tmp = mHashTable.getMetadata(point); + if (!mDenseLookup.empty()) { + if (coord.x < 0 || coord.y < 0 || coord.z < 0) { + return nullptr; + } + return mDenseLookup[coord.mPitch(mBBox)]; + } + + BitBlock* const* tmp = mHashTable.getMetadata(coord); if (tmp == nullptr) { return nullptr; } @@ -704,11 +824,18 @@ class SparseBitBlocks */ auto isActivePoint(const Neon::index_3d& point) const -> bool { - BitBlock* bitBlock = getBitBlockPrt(point); + Neon::index_3d block_coord(point.x / BitBlock::blockSize.x, + point.y / BitBlock::blockSize.y, + point.z / BitBlock::blockSize.z); + Neon::index_3d local(point.x % BitBlock::blockSize.x, + point.y % BitBlock::blockSize.y, + point.z % BitBlock::blockSize.z); + + BitBlock* bitBlock = getBitBlockPrt(block_coord); if (bitBlock == nullptr) { return false; } - return bitBlock->isON(point); + return bitBlock->isON(local); }; }; } // namespace Neon::domain::details::mGrid diff --git a/libNeonDomain/include/Neon/domain/tools/Partitioner1D.h b/libNeonDomain/include/Neon/domain/tools/Partitioner1D.h index fe16898f..a7274b70 100644 --- a/libNeonDomain/include/Neon/domain/tools/Partitioner1D.h +++ b/libNeonDomain/include/Neon/domain/tools/Partitioner1D.h @@ -167,6 +167,9 @@ class Partitioner1D mData->spanDecomposition); timeMamager.stop_with_trace("SpanClassifier"); + // The classifier was the only consumer of the block activity mask. + mData->spanDecomposition->releaseBlockActiveMask(); + timeMamager.start_with_trace("SpanLayout"); mData->mSpanLayout = std::make_shared( backend, diff --git a/libNeonDomain/include/Neon/domain/tools/partitioning/SpanClassifier.h b/libNeonDomain/include/Neon/domain/tools/partitioning/SpanClassifier.h index 487ce9f9..45c939ea 100644 --- a/libNeonDomain/include/Neon/domain/tools/partitioning/SpanClassifier.h +++ b/libNeonDomain/include/Neon/domain/tools/partitioning/SpanClassifier.h @@ -193,6 +193,19 @@ SpanClassifier::SpanClassifier(const Neon::Backend& back return maxRadius; }(); + // SpanDecomposition has already swept this exact block space with this exact + // activation lambda. Reuse its verdict rather than paying for it twice: for a + // large finest level the sweep is billions of lambda evaluations. + auto const& blockActive = mSpanDecomposition->getBlockActiveMask(); + bool const hasBlockActive = !blockActive.empty(); + + auto const blockPitch = [block3DSpan](int bx, int by, int bz) -> size_t { + return static_cast(bx) + + static_cast(by) * static_cast(block3DSpan.x) + + static_cast(bz) * static_cast(block3DSpan.x) * + static_cast(block3DSpan.y); + }; + // For each Partition backend.devSet() .forEachSetIdxSeq( @@ -218,35 +231,52 @@ SpanClassifier::SpanClassifier(const Neon::Backend& back auto inspectBlock = [&](int bx, int by, int bz, ByPartition byPartition, ByDirection byDirection) { - Neon::int32_3d blockOrigin = block3dIdxToBlockOrigin({bx, by, bz}); - - bool doBreak = false; bool isActiveBlock = false; ByDomain byDomain = ByDomain::bulk; - for (int z = 0; (z < dataBlockSize3D.z && !doBreak); z++) { - for (int y = 0; (y < dataBlockSize3D.y && !doBreak); y++) { - for (int x = 0; (x < dataBlockSize3D.x && !doBreak); x++) { - - const Neon::int32_3d globalId = getVoxelAbsolute3DIdx(blockOrigin, - {x, y, z}); - if (globalId < domainSize * discreteVoxelSpacing) { - - if constexpr (std::is_same_v) { - if (activeCellLambda(globalId)) { - byDomain = ByDomain::bulk; - isActiveBlock = true; - doBreak = true; - break; - } - } else if constexpr (std::is_same_v::type, bool>) { - NEON_THROW_UNSUPPORTED_OPERATION("bool"); - } else if constexpr (std::is_same_v::type, ByDomain>) { - auto whatdomain = bcLambda(globalId); - if (activeCellLambda(globalId)) { - isActiveBlock = true; - if (whatdomain == ByDomain::bc) { - byDomain = ByDomain::bc; + + // Whether the block holds anything at all is already known. + bool needsVoxelSweep = !hasBlockActive; + if (hasBlockActive) { + if (blockActive[blockPitch(bx, by, bz)] == 0) { + return; + } + isActiveBlock = true; + } + // The bc classification, however, is not something the + // decomposition sweep computes, so it still needs the voxels. + if constexpr (!std::is_same_v) { + needsVoxelSweep = true; + } + + if (needsVoxelSweep) { + Neon::int32_3d blockOrigin = block3dIdxToBlockOrigin({bx, by, bz}); + + bool doBreak = false; + for (int z = 0; (z < dataBlockSize3D.z && !doBreak); z++) { + for (int y = 0; (y < dataBlockSize3D.y && !doBreak); y++) { + for (int x = 0; (x < dataBlockSize3D.x && !doBreak); x++) { + + const Neon::int32_3d globalId = getVoxelAbsolute3DIdx(blockOrigin, + {x, y, z}); + if (globalId < domainSize * discreteVoxelSpacing) { + + if constexpr (std::is_same_v) { + if (activeCellLambda(globalId)) { + byDomain = ByDomain::bulk; + isActiveBlock = true; doBreak = true; + break; + } + } else if constexpr (std::is_same_v::type, bool>) { + NEON_THROW_UNSUPPORTED_OPERATION("bool"); + } else if constexpr (std::is_same_v::type, ByDomain>) { + auto whatdomain = bcLambda(globalId); + if (activeCellLambda(globalId)) { + isActiveBlock = true; + if (whatdomain == ByDomain::bc) { + byDomain = ByDomain::bc; + doBreak = true; + } } } } diff --git a/libNeonDomain/include/Neon/domain/tools/partitioning/SpanDecomposition.h b/libNeonDomain/include/Neon/domain/tools/partitioning/SpanDecomposition.h index 977ebaf6..b45a9e74 100644 --- a/libNeonDomain/include/Neon/domain/tools/partitioning/SpanDecomposition.h +++ b/libNeonDomain/include/Neon/domain/tools/partitioning/SpanDecomposition.h @@ -1,4 +1,8 @@ #pragma once + +#include +#include + #include "Neon/core/core.h" #include "Neon/set/Containter.h" @@ -41,6 +45,27 @@ class SpanDecomposition auto getLastZSliceIdx() const -> const Neon::set::DataSet&; + /** + * Per-block activity flags over the full block3DSpan, laid out with bx varying + * fastest. A block is flagged when it holds at least one active voxel. + * + * Deciding this is the only thing the decomposition sweep does, and it is exactly + * the test SpanClassifier would otherwise run a second time over the same block + * space. Publishing it lets the classifier skip that repeat. + * + * An empty vector means the mask is not available (it has been released); callers + * must fall back to evaluating the activation lambda themselves. + */ + auto getBlockActiveMask() const + -> const std::vector&; + + /** + * Frees the block activity mask. For a large finest level it is tens of MB that + * would otherwise be retained for the lifetime of the grid. + */ + auto releaseBlockActiveMask() + -> void; + auto toString(Neon::Backend const&) const -> std::string; @@ -48,6 +73,7 @@ class SpanDecomposition Neon::set::DataSet mZFirstIdx; Neon::set::DataSet mZLastIdx; Neon::set::DataSet mNumBlocks; + std::vector mBlockActive; size_t mDomainBlocksCount; }; @@ -68,6 +94,13 @@ SpanDecomposition::SpanDecomposition(const Neon::Backend& backend, mDomainBlocksCount = 0; std::vector nBlockProjectedToZ(block3DSpan.z); + // Record which blocks turned out to be active so that SpanClassifier does not have + // to rediscover it. One byte per block rather than one bit, so that neighbouring + // blocks handled by different threads never touch the same location. + size_t const sliceStride = static_cast(block3DSpan.x) * + static_cast(block3DSpan.y); + mBlockActive.assign(sliceStride * static_cast(block3DSpan.z), uint8_t(0)); + for (int bz = 0; bz < block3DSpan.z; bz++) { size_t count_on_bz = 0; #pragma omp parallel for reduction(+ : count_on_bz) schedule(static) collapse(2) @@ -86,6 +119,9 @@ SpanDecomposition::SpanDecomposition(const Neon::Backend& backend, if (activeCellLambda(id)) { doBreak = true; count_on_bz++; + mBlockActive[bx64 + + by64 * static_cast(block3DSpan.x) + + static_cast(bz) * sliceStride] = 1; } } } diff --git a/libNeonDomain/src/domain/details/mGrid/mGrid.cpp b/libNeonDomain/src/domain/details/mGrid/mGrid.cpp index 5ea1141e..b453f816 100644 --- a/libNeonDomain/src/domain/details/mGrid/mGrid.cpp +++ b/libNeonDomain/src/domain/details/mGrid/mGrid.cpp @@ -260,11 +260,10 @@ mGrid::mGrid( // } else { if (activeCellLambda[l](voxel)) { containVoxels = true; -#pragma omp critical - { - // Set the bitmask for this voxel if it is active - setLevelBitMask(l, {bx, by, bz}, {x, y, z}); - } + // Set the bitmask for this voxel if it is active. + // setLevelBitMask is thread safe: the bit write is + // atomic and only block allocation takes a lock. + setLevelBitMask(l, {bx, by, bz}, {x, y, z}); } //} } @@ -281,10 +280,7 @@ mGrid::mGrid( const Neon::int32_3d voxel = mData->mDescriptor.parentToChild(blockOrigin, l, {x, y, z}); if (voxel < domainSize) { -#pragma omp critical - { - setLevelBitMask(l, {bx, by, bz}, {x, y, z}); - } + setLevelBitMask(l, {bx, by, bz}, {x, y, z}); } } } @@ -301,11 +297,9 @@ mGrid::mGrid( // Find local position within the parent block Neon::int32_3d indexInParentBlock = mData->mDescriptor.toLocalIndex(blockOrigin, l + 1); -#pragma omp critical - { - // Activate the corresponding voxel in the parent block - setLevelBitMask(l + 1, parentBlock, indexInParentBlock); - } + + // Activate the corresponding voxel in the parent block + setLevelBitMask(l + 1, parentBlock, indexInParentBlock); } } } @@ -396,71 +390,69 @@ mGrid::mGrid( // Process levels from coarsest to finest (skip level 0 which has no children) for (int l = mData->mDescriptor.getDepth() - 1; l > 0; --l) { - const int refFactor = mData->mDescriptor.getRefFactor(l); - // Process all blocks at this level in parallel -#pragma omp parallel for collapse(3) schedule(static) - for (size_t bzUint64 = 0; bzUint64 < static_cast(mData->mTotalNumBlocks[l].z); bzUint64++) { - for (size_t byUint64 = 0; byUint64 < static_cast(mData->mTotalNumBlocks[l].y); byUint64++) { - for (size_t bxUint64 = 0; bxUint64 < static_cast(mData->mTotalNumBlocks[l].x); bxUint64++) { - int const bz = static_cast(bzUint64); - int const by = static_cast(byUint64); - int const bx = static_cast(bxUint64); + // A voxel's position in the base index space is its bitmask coordinate + // scaled by the level's voxel spacing: + // blockOrigin + localChild * spacing(l-1) + // = (blockID * 2 + localChild) * spacing(l-1) + // = bxyz * spacing(l-1) + const int voxelSpacing = mData->mDescriptor.getSpacing(l - 1); + + // Only an allocated bitmask region can hold an active voxel, and this pass + // does nothing to inactive ones. Walking the allocated regions instead of + // the dense block space turns hundreds of millions of block visits into a + // few tens of thousands of region visits. + std::vector const activeRegions = + mData->sparseLevelsBitmask.at(l).getAllocatedBlockCoords(); + +#pragma omp parallel for schedule(dynamic, 8) + for (int64_t r = 0; r < static_cast(activeRegions.size()); r++) { + mData->sparseLevelsBitmask.at(l).forEachActivePointInBlock( + activeRegions[r], + [&](const Neon::index_3d& bxyz) { + // Another thread may have culled this voxel after its brick was read + if (!levelBitMaskIsSetAt(l, bxyz)) { + return; + } - const Neon::index_3d blockOrigin = mData->mDescriptor.toBaseIndexSpace({bx, by, bz}, l + 1); + const Neon::int32_3d voxel = bxyz * voxelSpacing; - // Check each voxel in this block for potential culling - for (int z = 0; z < refFactor; z++) { - for (int y = 0; y < refFactor; y++) { - for (int x = 0; x < refFactor; x++) { + // Only cull if voxel is within domain and is refined + if (!(voxel < domainSize)) { + return; + } + if (!isRefined(l, voxel)) { + return; + } - // Only consider active voxels - if (levelBitMaskIsSet(l, {bx, by, bz}, {x, y, z})) { - - const Neon::int32_3d voxel = mData->mDescriptor.parentToChild(blockOrigin, l, {x, y, z}); - - // Only cull if voxel is within domain and is refined - if (voxel < domainSize) { - if (isRefined(l, voxel)) { - - // Check all 26 neighbors in 3D - // Deactivate only if ALL neighbors are also refined - bool deactivate = true; - for (int k = -1; k < 2; k++) { - for (int j = -1; j < 2; j++) { - for (int i = -1; i < 2; i++) { - if (i == 0 && j == 0 && k == 0) { - continue; // Skip center voxel - } - - const Neon::int32_3d neighborVoxel = mData->mDescriptor.neighbourBlock(voxel, l, {i, j, k}); - - // Check if neighbor is within domain bounds - if (neighborVoxel.x >= 0 && neighborVoxel.y >= 0 && neighborVoxel.z >= 0 && neighborVoxel < domainSize) { - // If any neighbor is not refined, don't deactivate - if (!isRefined(l, neighborVoxel)) { - deactivate = false; - } - } - } - } - } + // Check all 26 neighbors in 3D + // Deactivate only if ALL neighbors are also refined + bool deactivate = true; + for (int k = -1; k < 2 && deactivate; k++) { + for (int j = -1; j < 2 && deactivate; j++) { + for (int i = -1; i < 2 && deactivate; i++) { + if (i == 0 && j == 0 && k == 0) { + continue; // Skip center voxel + } - // Deactivate voxel if it and all neighbors are refined - if (deactivate) { -#pragma omp critical - { - clearLevelBitMask(l, {bx, by, bz}, {x, y, z}); - } - } - } + const Neon::int32_3d neighborVoxel = mData->mDescriptor.neighbourBlock(voxel, l, {i, j, k}); + + // Check if neighbor is within domain bounds + if (neighborVoxel.x >= 0 && neighborVoxel.y >= 0 && neighborVoxel.z >= 0 && neighborVoxel < domainSize) { + // If any neighbor is not refined, don't deactivate + if (!isRefined(l, neighborVoxel)) { + deactivate = false; } } } } } - } - } + + // Deactivate voxel if it and all neighbors are refined + if (deactivate) { + clearLevelBitMaskAt(l, bxyz); + } + }); } } } @@ -504,106 +496,98 @@ mGrid::mGrid( * - Reduced aliasing in multi-scale computations */ if (mData->mStrongBalanced) { - // Iteratively refine grid until strong balance condition is satisfied - bool again = true; - while (again) { - again = false; + // Iteratively refine grid until strong balance condition is satisfied. + // An int rather than a bool so that the worker threads can flag "converged = no" + // with an omp atomic write instead of a critical section. + int again = 1; + while (again != 0) { + again = 0; // Check all levels for balance violations for (int l = 0; l < mData->mDescriptor.getDepth(); ++l) { - const int refFactor = mData->mDescriptor.getRefFactor(l); const int childSpacing = mData->mDescriptor.getSpacing(l - 1); -#pragma omp parallel for collapse(3) schedule(static) - for (size_t bzUint64 = 0; bzUint64 < static_cast(mData->mTotalNumBlocks[l].z); bzUint64++) { - for (size_t byUint64 = 0; byUint64 < static_cast(mData->mTotalNumBlocks[l].y); byUint64++) { - for (size_t bxUint64 = 0; bxUint64 < static_cast(mData->mTotalNumBlocks[l].x); bxUint64++) { - int const bz = static_cast(bzUint64); - int const by = static_cast(byUint64); - int const bx = static_cast(bxUint64); - - // Check each voxel in the current block - for (int z = 0; z < refFactor; z++) { - for (int y = 0; y < refFactor; y++) { - for (int x = 0; x < refFactor; x++) { - - // Only process active voxels - if (levelBitMaskIsSet(l, {bx, by, bz}, {x, y, z})) { - - // Calculate global position of this voxel - const Neon::int32_3d voxel(bx * refFactor + x, - by * refFactor + y, - bz * refFactor + z); - - // Check all 26 neighbors for balance violations - for (int k = -1; k < 2; k++) { - for (int j = -1; j < 2; j++) { - for (int i = -1; i < 2; i++) { - if (i == 0 && j == 0 && k == 0) { - continue; // Skip center voxel - } - - // Calculate neighbor position - Neon::int32_3d proxyVoxel(voxel.x + i, - voxel.y + j, - voxel.z + k); - - // Convert to physical coordinates - const Neon::int32_3d proxyVoxelLocation(proxyVoxel.x * childSpacing, - proxyVoxel.y * childSpacing, - proxyVoxel.z * childSpacing); - - if (proxyVoxelLocation < domainSize && proxyVoxelLocation >= 0) { - - // Store previous level information for potential activation - Neon::int32_3d prv_nVoxelBlockOrigin(0), prv_nVoxelLocalID(0); - - // Search through all coarser levels to find neighbor - for (int l_n = l; l_n < mData->mDescriptor.getDepth(); ++l_n) { - const int l_n_ref_factor = mData->mDescriptor.getRefFactor(l_n); - - // Calculate block and local indices at level l_n - const Neon::int32_3d nVoxelBlockOrigin(proxyVoxel.x / l_n_ref_factor, - proxyVoxel.y / l_n_ref_factor, - proxyVoxel.z / l_n_ref_factor); - - const Neon::int32_3d nVoxelLocalID(proxyVoxel.x % l_n_ref_factor, - proxyVoxel.y % l_n_ref_factor, - proxyVoxel.z % l_n_ref_factor); - - // Check if neighbor exists at this level - if (levelBitMaskIsSet(l_n, nVoxelBlockOrigin, nVoxelLocalID)) { - - // Strong balance: neighbors can differ by at most 1 level - if (l_n == l || l_n == l + 1) { - break; // Balance satisfied - } else { -#pragma omp critical - { - // Balance violation: activate intermediate level - setLevelBitMask(l_n - 1, prv_nVoxelBlockOrigin, prv_nVoxelLocalID); - again = true; // Need another iteration - } - } - } - - // Move to next coarser level - proxyVoxel = nVoxelBlockOrigin; - - // Cache current level info for potential activation - prv_nVoxelBlockOrigin = nVoxelBlockOrigin; - prv_nVoxelLocalID = nVoxelLocalID; - } - } + // This pass only ever reads from and writes around *active* voxels, so + // walk the allocated bitmask regions rather than the dense block space. + // + // Taking a snapshot of the coordinates is deliberate: a balance violation + // activates voxels at coarser levels, which can allocate new regions + // there. Those levels are visited later in this same sweep (l ascending, + // and a violation at level l only ever writes to levels > l), so the list + // being iterated right now is never the one being appended to. + std::vector const activeRegions = + mData->sparseLevelsBitmask.at(l).getAllocatedBlockCoords(); + +#pragma omp parallel for schedule(dynamic, 8) + for (int64_t r = 0; r < static_cast(activeRegions.size()); r++) { + mData->sparseLevelsBitmask.at(l).forEachActivePointInBlock( + activeRegions[r], + // The bitmask coordinate is exactly the voxel position this pass + // works in: blockID * refFactor + localChild. + [&](const Neon::index_3d& voxel) { + // Check all 26 neighbors for balance violations + for (int k = -1; k < 2; k++) { + for (int j = -1; j < 2; j++) { + for (int i = -1; i < 2; i++) { + if (i == 0 && j == 0 && k == 0) { + continue; // Skip center voxel + } + + // Calculate neighbor position + Neon::int32_3d proxyVoxel(voxel.x + i, + voxel.y + j, + voxel.z + k); + + // Convert to physical coordinates + const Neon::int32_3d proxyVoxelLocation(proxyVoxel.x * childSpacing, + proxyVoxel.y * childSpacing, + proxyVoxel.z * childSpacing); + + if (proxyVoxelLocation < domainSize && proxyVoxelLocation >= 0) { + + // Store previous level information for potential activation + Neon::int32_3d prv_nVoxelBlockOrigin(0), prv_nVoxelLocalID(0); + + // Search through all coarser levels to find neighbor + for (int l_n = l; l_n < mData->mDescriptor.getDepth(); ++l_n) { + const int l_n_ref_factor = mData->mDescriptor.getRefFactor(l_n); + + // Calculate block and local indices at level l_n + const Neon::int32_3d nVoxelBlockOrigin(proxyVoxel.x / l_n_ref_factor, + proxyVoxel.y / l_n_ref_factor, + proxyVoxel.z / l_n_ref_factor); + + const Neon::int32_3d nVoxelLocalID(proxyVoxel.x % l_n_ref_factor, + proxyVoxel.y % l_n_ref_factor, + proxyVoxel.z % l_n_ref_factor); + + // Check if neighbor exists at this level + if (levelBitMaskIsSet(l_n, nVoxelBlockOrigin, nVoxelLocalID)) { + + // Strong balance: neighbors can differ by at most 1 level + if (l_n == l || l_n == l + 1) { + break; // Balance satisfied + } else { + // Balance violation: activate intermediate level + setLevelBitMask(l_n - 1, prv_nVoxelBlockOrigin, prv_nVoxelLocalID); + // Need another iteration +#pragma omp atomic write + again = 1; } } + + // Move to next coarser level + proxyVoxel = nVoxelBlockOrigin; + + // Cache current level info for potential activation + prv_nVoxelBlockOrigin = nVoxelBlockOrigin; + prv_nVoxelLocalID = nVoxelLocalID; } } } } } - } - } + }); } } } @@ -626,15 +610,49 @@ mGrid::mGrid( Neon::int32_3d levelDomainSize(mData->mTotalNumBlocks[l].x * blockSize, mData->mTotalNumBlocks[l].y * blockSize, mData->mTotalNumBlocks[l].z * blockSize); + + /** + * The bitmask coordinate of a voxel reduces to a single shift: + * + * blockID = id / getSpacing(l) = id / (2 * voxelSpacing) + * localID = (id / getSpacing(l-1)) % refFactor = (id / voxelSpacing) % 2 + * bxyz = blockID * 2 + localID + * + * mGrid is always an octree (refFactor == 2 at every level, enforced by the + * static_assert and the descriptor validation at the top of this constructor), + * so writing q = id / voxelSpacing gives bxyz = (q / 2) * 2 + q % 2, which is + * just q. Since voxelSpacing is a power of two, that is id >> levelShift. + * + * This matters because the partitioner evaluates this lambda roughly 1.6 + * billion times for level 0, and the childToParent/toLocalIndex form costs + * nine integer divisions per call with a divisor the compiler cannot see. + */ + if ((voxelSpacing & (voxelSpacing - 1)) != 0) { + NeonException exp("mGrid::mGrid"); + exp << "Level spacing is expected to be a power of two. Level = " << l + << " spacing = " << voxelSpacing; + NEON_THROW(exp); + } + int const levelShift = [voxelSpacing] { + int s = 0; + while ((1 << s) < voxelSpacing) { + ++s; + } + return s; + }(); + mData->grids[l] = InternalGrid( backend, levelDomainSize, - [&](Neon::int32_3d id) { + [this, l, levelShift, domainSize](Neon::int32_3d id) { + if (id.x < 0 || id.y < 0 || id.z < 0) { + return false; + } if (id < domainSize) { - Neon::index_3d blockID = mData->mDescriptor.childToParent(id, l); - Neon::index_3d localID = mData->mDescriptor.toLocalIndex(id, l); - return levelBitMaskIsSet(l, blockID, localID); + return levelBitMaskIsSetAt(l, Neon::index_3d(id.x >> levelShift, + id.y >> levelShift, + id.z >> levelShift)); } else { return false; } @@ -954,6 +972,20 @@ auto mGrid::levelBitMaskIsSet(int l, const Neon::index_3d& blockID, cons return mData->sparseLevelsBitmask.at(l).isActivePoint(bxyz); }; +/** + * @brief Check if a voxel is active, given its bitmask coordinate directly. + * + * Same query as levelBitMaskIsSet(), but for callers that have already reduced + * (blockID, localChild) to the single bitmask coordinate blockID * 2 + localChild. + * See the InternalGrid activation lambda for why that reduction is worth doing + * outside this function. + */ +template +auto mGrid::levelBitMaskIsSetAt(int l, const Neon::index_3d& bxyz) const -> bool +{ + return mData->sparseLevelsBitmask.at(l).isActivePoint(bxyz); +}; + /** * @brief Activate a voxel at a specific resolution level. @@ -970,7 +1002,7 @@ auto mGrid:: setLevelBitMask(int l, const Neon::index_3d& blockID, co auto const bxyz = blockID * 2 + localChild; - return mData->sparseLevelsBitmask.at(l).template activatePoint(bxyz); + return mData->sparseLevelsBitmask.at(l).template activatePoint(bxyz); }; /** @@ -988,7 +1020,19 @@ auto mGrid::clearLevelBitMask(int l, const Neon::index_3d& blockID, cons auto const bxyz = blockID * 2 + localChild; - return mData->sparseLevelsBitmask.at(l).template removePoint(bxyz); + return mData->sparseLevelsBitmask.at(l).template removePoint(bxyz); +}; + +/** + * @brief Deactivate a voxel, given its bitmask coordinate directly. + * + * Counterpart of levelBitMaskIsSetAt() for the sweeps that already work in bitmask + * coordinates rather than in (blockID, localChild) pairs. + */ +template +auto mGrid::clearLevelBitMaskAt(int l, const Neon::index_3d& bxyz) -> void +{ + return mData->sparseLevelsBitmask.at(l).template removePoint(bxyz); }; /** diff --git a/libNeonDomain/src/domain/tools/partitioning/SpanDecomposition.cpp b/libNeonDomain/src/domain/tools/partitioning/SpanDecomposition.cpp index 78484998..efcdc397 100644 --- a/libNeonDomain/src/domain/tools/partitioning/SpanDecomposition.cpp +++ b/libNeonDomain/src/domain/tools/partitioning/SpanDecomposition.cpp @@ -16,6 +16,16 @@ auto SpanDecomposition::getLastZSliceIdx() const -> const Neon::set::DataSet const std::vector& +{ + return mBlockActive; +} + +auto SpanDecomposition::releaseBlockActiveMask() -> void +{ + std::vector().swap(mBlockActive); +} + auto SpanDecomposition::toString(Neon::Backend const& bk) const -> std::string { std::stringstream s;