From e6d7fbcda3264ff04e5e1cea968a66e9f1abb6b8 Mon Sep 17 00:00:00 2001 From: hamidrezanorouzi Date: Thu, 27 Oct 2022 14:19:53 +0330 Subject: [PATCH] contact search NBS refactored --- .../ContactSearch/ContactSearch.H | 19 +- src/Interaction/contactSearch/methods/NBS.H | 424 ++------------ .../contactSearch/methods/NBSLevel.H | 70 +++ .../contactSearch/methods/NBSLevel0.H | 211 +++++++ .../contactSearch/methods/NBSLoop.H | 20 +- .../contactSearch/methods/NBSold.H | 523 ++++++++++++++++++ .../contactSearch/methods/mapperNBS.H | 331 +++++++++++ .../sphereInteraction/sphereInteraction.H | 1 + src/phasicFlow/Kokkos/KokkosUtilities.H | 36 ++ .../positionRandom/positionRandom.C | 6 +- .../postprocessPhasicFlow/fieldOperations.H | 10 +- .../postprocessPhasicFlow/pointRectCell.H | 33 +- 12 files changed, 1256 insertions(+), 428 deletions(-) create mode 100644 src/Interaction/contactSearch/methods/NBSLevel.H create mode 100644 src/Interaction/contactSearch/methods/NBSLevel0.H create mode 100644 src/Interaction/contactSearch/methods/NBSold.H create mode 100644 src/Interaction/contactSearch/methods/mapperNBS.H diff --git a/src/Interaction/contactSearch/ContactSearch/ContactSearch.H b/src/Interaction/contactSearch/ContactSearch/ContactSearch.H index 276e1ae3..4220b820 100644 --- a/src/Interaction/contactSearch/ContactSearch/ContactSearch.H +++ b/src/Interaction/contactSearch/ContactSearch/ContactSearch.H @@ -30,7 +30,7 @@ namespace pFlow { template< - template class BaseMethod, + template class BaseMethod, template class WallMapping > class ContactSearch @@ -49,10 +49,7 @@ public: using ParticleContactSearchType = BaseMethod< - ExecutionSpace, - IdType, - IndexType - >; + ExecutionSpace>; using WallMappingType = WallMapping< @@ -141,7 +138,7 @@ public: bool force = false) override { - //output<<" ContactSearch::broadSearch::PP before.\n"; + if(particleContactSearch_) { auto activeRange = this->Particles().activeRange(); @@ -149,11 +146,11 @@ public: sphereSphereTimer_.start(); if(this->Particles().allActive()) - { + { particleContactSearch_().broadSearch(ppPairs, activeRange, force); } else - { + { particleContactSearch_().broadSearch(ppPairs, activeRange, this->Particles().activePointsMaskD(), force); } @@ -162,9 +159,7 @@ public: } else return false; - //output<<" ContactSearch::broadSearch::PP fater.\n"; - - //output<<" ContactSearch::broadSearch::PW before.\n"; + if(wallMapping_) { sphereWallTimer_.start(); @@ -174,7 +169,7 @@ public: else return false; - //output<<" ContactSearch::broadSearch::PW after.\n"; + return true; } diff --git a/src/Interaction/contactSearch/methods/NBS.H b/src/Interaction/contactSearch/methods/NBS.H index cf5aed45..a3610066 100644 --- a/src/Interaction/contactSearch/methods/NBS.H +++ b/src/Interaction/contactSearch/methods/NBS.H @@ -22,75 +22,35 @@ Licence: #ifndef __NBS_H__ #define __NBS_H__ -#include "cells.H" -#include "contactSearchFunctions.H" -#include "baseAlgorithms.H" +#include "NBSLevel0.H" namespace pFlow { -template< - typename executionSpace, - typename idType, - typename indexType=int32 - > +template class NBS -: - public cells { public: - using IdType = idType; - - using IndexType = indexType; + using NBSLevel0Type = NBSLevel0; - using Cells = cells; + using IdType = typename NBSLevel0Type::IdType; + + using IndexType = typename NBSLevel0Type::IndexType; + + using Cells = typename NBSLevel0Type::Cells; using CellType = typename Cells::CellType; - using ExecutionSpace= executionSpace; + using execution_space = typename NBSLevel0Type::execution_space; - using memory_space = typename ExecutionSpace::memory_space; + using memory_space = typename NBSLevel0Type::memory_space; - struct TagBuild{}; - - struct TagFindPairs{}; - - class cellIterator - { - private: - ViewType3D head_; - - ViewType1D next_; - - public: - - cellIterator(ViewType3D head, ViewType1D next) - : - head_(head), - next_(next) - {} - - INLINE_FUNCTION_HD - Cells cellsSize()const { - return Cells(head_.extent(0), head_.extent(1), head_.extent(2));} - - INLINE_FUNCTION_HD - int32 start(indexType i, indexType j, indexType k)const { - return head_(i,j,k); } - - INLINE_FUNCTION_HD - int32 getNext(int32 n)const { - if(n<0) return n; - return next_(n); } - }; protected: - int32 capacity_ = 1; - - real sizeRatio_ = 1.0; + real sizeRatio_ = 1.0; int32 updateFrequency_= 1; @@ -98,74 +58,10 @@ protected: bool performedSearch_ = false; - ViewType1D pointPosition_; - - ViewType1D diameter_; - - ViewType3D head_; - - ViewType1D next_; - - - INLINE_FUNCTION_H - void nullify() - { - fill( - head_, - range(0,this->nx()), - range(0,this->ny()), - range(0,this->nz()), - static_cast(-1) - ); - - fill( - next_, - range(0,capacity_), - static_cast(-1) - ); - } - - void nullify(range nextRng) - { - fill( - head_, - range(0,this->nx()), - range(0,this->ny()), - range(0,this->nz()), - static_cast(-1) - ); - - fill( - next_, - nextRng, - static_cast(-1) - ); - } - - using mdrPolicyFindPairs = - Kokkos::MDRangePolicy< - Kokkos::Rank<3>, - Kokkos::Schedule, - ExecutionSpace>; + NBSLevel0Type NBSLevel0_; private: - void checkAllocateNext(int newCap) - { - - if( capacity_ < newCap) - { - capacity_ = newCap; - Kokkos::realloc(next_, capacity_); - } - - } - - void allocateHead() - { - Kokkos::realloc(head_, this->nx(), this->ny(), this->nz()); - } - bool performSearch() { if(currentIter_ % updateFrequency_ == 0) @@ -180,79 +76,34 @@ private: } } - static INLINE_FUNCTION_HD - void Swap(int32& x, int32& y) - { - int32 tmp = x; - x = y; - y = tmp; - } - public: TypeNameNV("NBS"); - NBS( - const box& domain, - real cellSize, - const ViewType1D& position, - const ViewType1D& diam, - int32 initialContainerSize = 1) - : - Cells(domain, cellSize), - pointPosition_(position), - diameter_(diam), - head_("NBS::head_",1,1,1), //, this->nx(), this->ny(), this->nz()), - next_("NBS::next_",1) //,position.size()), - { - checkAllocateNext(pointPosition_.size()); - allocateHead(); - } - - NBS( - const box& domain, - int32 nx, - int32 ny, - int32 nz, - const ViewType1D& position, - const ViewType1D& diam, - int32 initialContainerSize = 1) - : - Cells(domain, nx, ny, nz), - pointPosition_(position), - diameter_(diam), - head_("NBS::head_",nx,ny,nz), //, this->nx(), this->ny(), this->nz()), - next_("NBS::next_",1) //,position.size()), - { - checkAllocateNext(pointPosition_.size()); - } - NBS( dictionary dict, const box& domain, - real cellSize, + real maxSize, const ViewType1D& position, - const ViewType1D& diam, - int32 initialContainerSize = 1 - ) + const ViewType1D& diam) : - Cells(domain, cellSize), - pointPosition_(position), - diameter_(diam), - head_("NBS::head_",1,1,1), //, this->nx(), this->ny(), this->nz()), - next_("NBS::next_",1) //,position.size()), - { - - updateFrequency_ = max( - dict.getVal("updateFrequency"),1 ); - - sizeRatio_ = max(dict.getVal( - "sizeRatio"),1.0); - - this->setCellSize(cellSize*sizeRatio_); - checkAllocateNext(pointPosition_.size()); - allocateHead(); - } + sizeRatio_( + max( + dict.getValOrSet("sizeRatio", 1.0), + 1.0 + )), + updateFrequency_( + max( + dict.getValOrSet("updateFrequency", 1), + 1 + )), + NBSLevel0_( + domain, + maxSize*sizeRatio_, + sizeRatio_, + position, + diam) + {} INLINE_FUNCTION_HD NBS(const NBS&) = default; @@ -275,19 +126,20 @@ public: return performedSearch_; } - const auto& Head()const + auto getCellIterator()const { - return head_; + return NBSLevel0_.getCellIterator(); } - const auto& Next()const + auto getCells()const { - return next_; + return NBSLevel0_.getCells(); } - cellIterator getCellIterator()const + bool objectSizeChanged(int32 newSize) { - return cellIterator(head_, next_); + NBSLevel0_.checkAllocateNext(newSize); + return true; } // - Perform the broad search to find pairs @@ -302,13 +154,12 @@ public: performedSearch_ = false; if( !performSearch() ) return true; - //Info<<"NBS::broadSearch(PairsContainer& pairs, range activeRange, bool force=false) before build"<(*this); - auto points = pointPosition_; - auto next = next_; - auto head = head_; - - Kokkos::RangePolicy< - Kokkos::IndexType, - Kokkos::Schedule, - ExecutionSpace> rPolicy(activeRange.first, activeRange.second); - - Kokkos::parallel_for( - "NBS::build", - rPolicy, - LAMBDA_HD(int32 i){ - - CellType ind = cellIndex.pointIndex(points[i]); - int32 old = Kokkos::atomic_exchange(&head(ind.x(), ind.y(), ind.z()), i); - next[i] = old; - }); - Kokkos::fence(); - - } - - - template - INLINE_FUNCTION_H - void build(range activeRange, IncludeFunction incld) - { - checkAllocateNext(activeRange.second); - nullify(activeRange); - - Cells cellIndex = static_cast(*this); - auto points = pointPosition_; - auto next = next_; - auto head = head_; - - Kokkos::RangePolicy< - Kokkos::IndexType, - Kokkos::Schedule, - ExecutionSpace> rPolicy(activeRange.first, activeRange.second); - - Kokkos::parallel_for( - "NBS::build", - rPolicy, - LAMBDA_HD(int32 i){ - - if( incld(i) ) - { - CellType ind = cellIndex.pointIndex(points[i]); - auto old = Kokkos::atomic_exchange(&head(ind.x(), ind.y(), ind.z()), i); - next[i] = old; - } - }); - Kokkos::fence(); - - } - - // - build based on all points in range [0, numPoints_) - INLINE_FUNCTION_H - void buildCheckInDomain(range activeRange) - { - checkAllocateNext(activeRange.second); - nullify(activeRange); - - Cells cellIndex = static_cast(*this); - auto points = pointPosition_; - auto next = next_; - auto head = head_; - - Kokkos::RangePolicy< - Kokkos::IndexType, - Kokkos::Schedule, - ExecutionSpace> rPolicy(activeRange.first, activeRange.second); - - Kokkos::parallel_for( - "NBS::buildCheckInDomain", - rPolicy, - LAMBDA_HD(int32 i){ - - CellType ind; - if( cellIndex.pointIndexInDomain(points[i], ind) ) - { - int32 old = Kokkos::atomic_exchange(&head(ind.x(), ind.y(), ind.z()), i); - next[i] = old; - } - - } - ); - Kokkos::fence(); - - } - - template - INLINE_FUNCTION_H - void buildCheckInDomain(range activeRange, IncludeFunction incld) - { - checkAllocateNext(activeRange.second); - nullify(activeRange); - - Cells cellIndex = static_cast(*this); - auto points = pointPosition_; - auto next = next_; - auto head = head_; - - Kokkos::RangePolicy< - Kokkos::IndexType, - Kokkos::Schedule, - ExecutionSpace> rPolicy(activeRange.first, activeRange.second); - - Kokkos::parallel_for( - "NBS::buildCheckInDomain", - rPolicy, - LAMBDA_HD(int32 i){ - - CellType ind; - if( incld(i) && cellIndex.pointIndexInDomain(points[i], ind) ) - { - auto old = Kokkos::atomic_exchange(&head(ind.x(), ind.y(), ind.z()), i); - next[i] = old; - } - }); - Kokkos::fence(); - } - - - template - INLINE_FUNCTION_H - bool findPairs(PairsContainer& pairs) - { - mdrPolicyFindPairs - mdrPolicy( - {0,0,0}, - {this->nx(),this->ny(),this->nz()} ); - - int32 getFull = 1; - - // loop until the container size fits the numebr of contact pairs - while (getFull > 0) - { - - getFull = 0; - - Kokkos::parallel_reduce ( - "NBS::broadSearch", - mdrPolicy, - CLASS_LAMBDA_HD(int32 i, int32 j, int32 k, int32& getFullUpdate){ - #include "NBSLoop.H" - }, getFull); - - if(getFull) - { - // - resize the container - // note that getFull now shows the number of failed insertions. - uint32 len = max(getFull,50) ; - - auto oldCap = pairs.capacity(); - - pairs.increaseCapacityBy(len); - - Info<< "The contact pair container capacity increased from "<< - oldCap << " to "< +class +NBSLevel +{ +public: + + + using execution_space= executionSpace; + + using memory_space = typename execution_space::memory_space; + + using rangePolicy = Kokkos::RangePolicy< + Kokkos::IndexType, + Kokkos::Schedule, + execution_space> + +protected: + + int32x3 gridExtent_; + + ViewType3D head_; + + int8 level_ = 1; + + +public: + + NBSLevel(int8 lvl, int32x3 gridExtent) + : + gridExtent_(gridExtent), + head_("NBSLevel::head", gridExtent.x(), gridExtent.y(), gridExtent.z()), + level_(lvl) + {} + + INLINE_FUNCION_HD + auto& head(int32 i, int32 j, int32 k) + { + return head_(i,j,k); + } + + INLINE_FUNCION_HD + auto& head() + { + return head_; + } + + INLINE_FUNCION_HD + auto level()const + { + return level_; + } + + INLINE_FUNCION_HD + const auto& gridExtent()const + { + return gridExtent_; + } + +}; + + + +#endif diff --git a/src/Interaction/contactSearch/methods/NBSLevel0.H b/src/Interaction/contactSearch/methods/NBSLevel0.H new file mode 100644 index 00000000..ea7dc275 --- /dev/null +++ b/src/Interaction/contactSearch/methods/NBSLevel0.H @@ -0,0 +1,211 @@ +/*------------------------------- phasicFlow --------------------------------- + O C enter of + O O E ngineering and + O O M ultiscale modeling of + OOOOOOO F luid flow +------------------------------------------------------------------------------ + Copyright (C): www.cemf.ir + email: hamid.r.norouzi AT gmail.com +------------------------------------------------------------------------------ +Licence: + This file is part of phasicFlow code. It is a free software for simulating + granular and multiphase flows. You can redistribute it and/or modify it under + the terms of GNU General Public License v3 or any other later versions. + + phasicFlow is distributed to help others in their research in the field of + granular and multiphase flows, but WITHOUT ANY WARRANTY; without even the + implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. + +-----------------------------------------------------------------------------*/ + + +#ifndef __NBSLevel0_H__ +#define __NBSLevel0_H__ + +#include "mapperNBS.H" + +namespace pFlow +{ + + +template +class NBSLevel0 +: + public mapperNBS +{ +public: + + using MapperType = mapperNBS; + + using IdType = typename MapperType::IdType; + + using IndexType = typename MapperType::IndexType; + + using Cells = typename MapperType::Cells; + + using CellType = typename Cells::CellType; + + using execution_space = typename MapperType::execution_space; + + using memory_space = typename MapperType::memory_space; + + struct TagFindPairs{}; + + +protected: + + real sizeRatio_ = 1.0; + + // borrowed ownership + ViewType1D diameter_; + + + using mdrPolicyFindPairs = + Kokkos::MDRangePolicy< + Kokkos::Rank<3>, + Kokkos::Schedule, + execution_space>; + +private: + + + static INLINE_FUNCTION_HD + void Swap(int32& x, int32& y) + { + int32 tmp = x; + x = y; + y = tmp; + } + +public: + + TypeNameNV("NBSLevel0"); + + NBSLevel0( + const box& domain, + real cellSize, + const ViewType1D& position, + const ViewType1D& diam) + : + MapperType(domain, cellSize, position), + diameter_(diam) + {} + + NBSLevel0( + const box& domain, + int32 nx, + int32 ny, + int32 nz, + const ViewType1D& position, + const ViewType1D& diam) + : + MapperType(domain, nx, ny, nz, position), + diameter_(diam) + { } + + NBSLevel0( + const box& domain, + real cellSize, + real sizeRatio, + const ViewType1D& position, + const ViewType1D& diam) + : + MapperType(domain, cellSize, position), + sizeRatio_(sizeRatio), + diameter_(diam) + {} + + INLINE_FUNCTION_HD + NBSLevel0(const NBSLevel0&) = default; + + INLINE_FUNCTION_HD + NBSLevel0& operator = (const NBSLevel0&) = default; + + INLINE_FUNCTION_HD + ~NBSLevel0()=default; + + // - Perform the broad search to find pairs + // with force = true, perform broad search regardless of + // updateFrequency_ value + // on all the points in the range of [0,numPoints_) + template + bool broadSearch(PairsContainer& pairs, range activeRange) + { + + + this->build(activeRange); + + findPairs(pairs); + + return true; + } + + // - Perform the broad search to find pairs, + // ignore particles with incld(i) = true, + // with force = true, perform broad search regardless of + // updateFrequency_ value + template + bool broadSearch(PairsContainer& pairs, range activeRange, IncludeFunction incld) + { + + this->build(activeRange, incld); + + findPairs(pairs); + + return true; + } + + template + INLINE_FUNCTION_H + bool findPairs(PairsContainer& pairs) + { + mdrPolicyFindPairs + mdrPolicy( + {0,0,0}, + {this->nx(),this->ny(),this->nz()} ); + + int32 getFull = 1; + + + // loop until the container size fits the numebr of contact pairs + while (getFull > 0) + { + + getFull = 0; + + Kokkos::parallel_reduce ( + "NBSLevel0::findPairs", + mdrPolicy, + CLASS_LAMBDA_HD(int32 i, int32 j, int32 k, int32& getFullUpdate){ + #include "NBSLoop.H" + }, getFull); + + + if(getFull) + { + // - resize the container + // note that getFull now shows the number of failed insertions. + uint32 len = max(getFull,500) ; + + auto oldCap = pairs.capacity(); + + pairs.increaseCapacityBy(len); + + Info<< "The contact pair container capacity increased from "<< + oldCap << " to "<head_(i,j,k); CellType currentCell(i,j,k); int32 n = -1; while( m > -1 ) { - auto p_m = pointPosition_[m]; + auto p_m = this->pointPosition_[m]; auto d_m = sizeRatio_* diameter_[m]; // the same cell - n = next_(m); + n = this->next_(m); while(n >-1) { - auto p_n = pointPosition_[n]; + auto p_n = this->pointPosition_[n]; auto d_n = sizeRatio_*diameter_[n]; if( sphereSphereCheck(p_m, p_n, d_m, d_n) ) { auto ln = n; auto lm = m; - //output<<"lm "<< lm <<" ln "<< ln << " i j k "<< int32x3(i,j,k)<ln) Swap(lm,ln); if( auto res = pairs.insert(lm,ln); res <0) { @@ -50,7 +50,7 @@ while( m > -1 ) } } - n = next_(n); + n = this->next_(n); } // neighbor cells @@ -74,11 +74,11 @@ while( m > -1 ) if( this->isInRange(neighborCell) ) { - n = head_(neighborCell.x(), neighborCell.y(), neighborCell.z()); + n = this->head_(neighborCell.x(), neighborCell.y(), neighborCell.z()); while( n>-1) { - auto p_n = pointPosition_[n]; + auto p_n = this->pointPosition_[n]; auto d_n = sizeRatio_*diameter_[n]; if(sphereSphereCheck(p_m, p_n, d_m, d_n)) @@ -91,11 +91,11 @@ while( m > -1 ) getFullUpdate++; } } - n = next_[n]; + n = this->next_[n]; } } } - m = next_[m]; + m = this->next_[m]; } diff --git a/src/Interaction/contactSearch/methods/NBSold.H b/src/Interaction/contactSearch/methods/NBSold.H new file mode 100644 index 00000000..cf5aed45 --- /dev/null +++ b/src/Interaction/contactSearch/methods/NBSold.H @@ -0,0 +1,523 @@ +/*------------------------------- phasicFlow --------------------------------- + O C enter of + O O E ngineering and + O O M ultiscale modeling of + OOOOOOO F luid flow +------------------------------------------------------------------------------ + Copyright (C): www.cemf.ir + email: hamid.r.norouzi AT gmail.com +------------------------------------------------------------------------------ +Licence: + This file is part of phasicFlow code. It is a free software for simulating + granular and multiphase flows. You can redistribute it and/or modify it under + the terms of GNU General Public License v3 or any other later versions. + + phasicFlow is distributed to help others in their research in the field of + granular and multiphase flows, but WITHOUT ANY WARRANTY; without even the + implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. + +-----------------------------------------------------------------------------*/ + + +#ifndef __NBS_H__ +#define __NBS_H__ + +#include "cells.H" +#include "contactSearchFunctions.H" +#include "baseAlgorithms.H" + +namespace pFlow +{ + + +template< + typename executionSpace, + typename idType, + typename indexType=int32 + > +class NBS +: + public cells +{ +public: + + using IdType = idType; + + using IndexType = indexType; + + using Cells = cells; + + using CellType = typename Cells::CellType; + + using ExecutionSpace= executionSpace; + + using memory_space = typename ExecutionSpace::memory_space; + + struct TagBuild{}; + + struct TagFindPairs{}; + + class cellIterator + { + private: + ViewType3D head_; + + ViewType1D next_; + + public: + + cellIterator(ViewType3D head, ViewType1D next) + : + head_(head), + next_(next) + {} + + INLINE_FUNCTION_HD + Cells cellsSize()const { + return Cells(head_.extent(0), head_.extent(1), head_.extent(2));} + + INLINE_FUNCTION_HD + int32 start(indexType i, indexType j, indexType k)const { + return head_(i,j,k); } + + INLINE_FUNCTION_HD + int32 getNext(int32 n)const { + if(n<0) return n; + return next_(n); } + }; + +protected: + + int32 capacity_ = 1; + + real sizeRatio_ = 1.0; + + int32 updateFrequency_= 1; + + int32 currentIter_ = 0; + + bool performedSearch_ = false; + + ViewType1D pointPosition_; + + ViewType1D diameter_; + + ViewType3D head_; + + ViewType1D next_; + + + INLINE_FUNCTION_H + void nullify() + { + fill( + head_, + range(0,this->nx()), + range(0,this->ny()), + range(0,this->nz()), + static_cast(-1) + ); + + fill( + next_, + range(0,capacity_), + static_cast(-1) + ); + } + + void nullify(range nextRng) + { + fill( + head_, + range(0,this->nx()), + range(0,this->ny()), + range(0,this->nz()), + static_cast(-1) + ); + + fill( + next_, + nextRng, + static_cast(-1) + ); + } + + using mdrPolicyFindPairs = + Kokkos::MDRangePolicy< + Kokkos::Rank<3>, + Kokkos::Schedule, + ExecutionSpace>; + +private: + + void checkAllocateNext(int newCap) + { + + if( capacity_ < newCap) + { + capacity_ = newCap; + Kokkos::realloc(next_, capacity_); + } + + } + + void allocateHead() + { + Kokkos::realloc(head_, this->nx(), this->ny(), this->nz()); + } + + bool performSearch() + { + if(currentIter_ % updateFrequency_ == 0) + { + currentIter_++; + return true; + + }else + { + currentIter_++; + return false; + } + } + + static INLINE_FUNCTION_HD + void Swap(int32& x, int32& y) + { + int32 tmp = x; + x = y; + y = tmp; + } + +public: + + TypeNameNV("NBS"); + + NBS( + const box& domain, + real cellSize, + const ViewType1D& position, + const ViewType1D& diam, + int32 initialContainerSize = 1) + : + Cells(domain, cellSize), + pointPosition_(position), + diameter_(diam), + head_("NBS::head_",1,1,1), //, this->nx(), this->ny(), this->nz()), + next_("NBS::next_",1) //,position.size()), + { + checkAllocateNext(pointPosition_.size()); + allocateHead(); + } + + NBS( + const box& domain, + int32 nx, + int32 ny, + int32 nz, + const ViewType1D& position, + const ViewType1D& diam, + int32 initialContainerSize = 1) + : + Cells(domain, nx, ny, nz), + pointPosition_(position), + diameter_(diam), + head_("NBS::head_",nx,ny,nz), //, this->nx(), this->ny(), this->nz()), + next_("NBS::next_",1) //,position.size()), + { + checkAllocateNext(pointPosition_.size()); + } + + NBS( + dictionary dict, + const box& domain, + real cellSize, + const ViewType1D& position, + const ViewType1D& diam, + int32 initialContainerSize = 1 + ) + : + Cells(domain, cellSize), + pointPosition_(position), + diameter_(diam), + head_("NBS::head_",1,1,1), //, this->nx(), this->ny(), this->nz()), + next_("NBS::next_",1) //,position.size()), + { + + updateFrequency_ = max( + dict.getVal("updateFrequency"),1 ); + + sizeRatio_ = max(dict.getVal( + "sizeRatio"),1.0); + + this->setCellSize(cellSize*sizeRatio_); + checkAllocateNext(pointPosition_.size()); + allocateHead(); + } + + INLINE_FUNCTION_HD + NBS(const NBS&) = default; + + INLINE_FUNCTION_HD + NBS& operator = (const NBS&) = default; + + INLINE_FUNCTION_HD + ~NBS()=default; + + //// - Methods + + bool enterBoadSearch()const + { + return currentIter_%updateFrequency_==0; + } + + bool performedSearch()const + { + return performedSearch_; + } + + const auto& Head()const + { + return head_; + } + + const auto& Next()const + { + return next_; + } + + cellIterator getCellIterator()const + { + return cellIterator(head_, next_); + } + + // - Perform the broad search to find pairs + // with force = true, perform broad search regardless of + // updateFrequency_ value + // on all the points in the range of [0,numPoints_) + template + bool broadSearch(PairsContainer& pairs, range activeRange, bool force=false) + { + + if(force) currentIter_ = 0; + performedSearch_ = false; + + if( !performSearch() ) return true; + //Info<<"NBS::broadSearch(PairsContainer& pairs, range activeRange, bool force=false) before build"< + bool broadSearch(PairsContainer& pairs, range activeRange, IncludeFunction incld, bool force = false) + { + if(force) currentIter_ = 0; + performedSearch_ = false; + + if( !performSearch() ) return true; + + build(activeRange, incld); + + findPairs(pairs); + + performedSearch_ = true; + return true; + } + + // - build based on all points in range [0, numPoints_) + INLINE_FUNCTION_H + void build(range activeRange) + { + checkAllocateNext(activeRange.second); + nullify(activeRange); + + Cells cellIndex = static_cast(*this); + auto points = pointPosition_; + auto next = next_; + auto head = head_; + + Kokkos::RangePolicy< + Kokkos::IndexType, + Kokkos::Schedule, + ExecutionSpace> rPolicy(activeRange.first, activeRange.second); + + Kokkos::parallel_for( + "NBS::build", + rPolicy, + LAMBDA_HD(int32 i){ + + CellType ind = cellIndex.pointIndex(points[i]); + int32 old = Kokkos::atomic_exchange(&head(ind.x(), ind.y(), ind.z()), i); + next[i] = old; + }); + Kokkos::fence(); + + } + + + template + INLINE_FUNCTION_H + void build(range activeRange, IncludeFunction incld) + { + checkAllocateNext(activeRange.second); + nullify(activeRange); + + Cells cellIndex = static_cast(*this); + auto points = pointPosition_; + auto next = next_; + auto head = head_; + + Kokkos::RangePolicy< + Kokkos::IndexType, + Kokkos::Schedule, + ExecutionSpace> rPolicy(activeRange.first, activeRange.second); + + Kokkos::parallel_for( + "NBS::build", + rPolicy, + LAMBDA_HD(int32 i){ + + if( incld(i) ) + { + CellType ind = cellIndex.pointIndex(points[i]); + auto old = Kokkos::atomic_exchange(&head(ind.x(), ind.y(), ind.z()), i); + next[i] = old; + } + }); + Kokkos::fence(); + + } + + // - build based on all points in range [0, numPoints_) + INLINE_FUNCTION_H + void buildCheckInDomain(range activeRange) + { + checkAllocateNext(activeRange.second); + nullify(activeRange); + + Cells cellIndex = static_cast(*this); + auto points = pointPosition_; + auto next = next_; + auto head = head_; + + Kokkos::RangePolicy< + Kokkos::IndexType, + Kokkos::Schedule, + ExecutionSpace> rPolicy(activeRange.first, activeRange.second); + + Kokkos::parallel_for( + "NBS::buildCheckInDomain", + rPolicy, + LAMBDA_HD(int32 i){ + + CellType ind; + if( cellIndex.pointIndexInDomain(points[i], ind) ) + { + int32 old = Kokkos::atomic_exchange(&head(ind.x(), ind.y(), ind.z()), i); + next[i] = old; + } + + } + ); + Kokkos::fence(); + + } + + template + INLINE_FUNCTION_H + void buildCheckInDomain(range activeRange, IncludeFunction incld) + { + checkAllocateNext(activeRange.second); + nullify(activeRange); + + Cells cellIndex = static_cast(*this); + auto points = pointPosition_; + auto next = next_; + auto head = head_; + + Kokkos::RangePolicy< + Kokkos::IndexType, + Kokkos::Schedule, + ExecutionSpace> rPolicy(activeRange.first, activeRange.second); + + Kokkos::parallel_for( + "NBS::buildCheckInDomain", + rPolicy, + LAMBDA_HD(int32 i){ + + CellType ind; + if( incld(i) && cellIndex.pointIndexInDomain(points[i], ind) ) + { + auto old = Kokkos::atomic_exchange(&head(ind.x(), ind.y(), ind.z()), i); + next[i] = old; + } + }); + Kokkos::fence(); + } + + + template + INLINE_FUNCTION_H + bool findPairs(PairsContainer& pairs) + { + mdrPolicyFindPairs + mdrPolicy( + {0,0,0}, + {this->nx(),this->ny(),this->nz()} ); + + int32 getFull = 1; + + // loop until the container size fits the numebr of contact pairs + while (getFull > 0) + { + + getFull = 0; + + Kokkos::parallel_reduce ( + "NBS::broadSearch", + mdrPolicy, + CLASS_LAMBDA_HD(int32 i, int32 j, int32 k, int32& getFullUpdate){ + #include "NBSLoop.H" + }, getFull); + + if(getFull) + { + // - resize the container + // note that getFull now shows the number of failed insertions. + uint32 len = max(getFull,50) ; + + auto oldCap = pairs.capacity(); + + pairs.increaseCapacityBy(len); + + Info<< "The contact pair container capacity increased from "<< + oldCap << " to "< +class mapperNBS +: + public cells +{ +public: + + using IdType = int32; + + using IndexType = int32; + + using Cells = cells; + + using CellType = typename Cells::CellType; + + using execution_space = executionSpace; + + using memory_space = typename execution_space::memory_space; + + class cellIterator + { + private: + ViewType3D head_; + + ViewType1D next_; + + public: + + cellIterator(ViewType3D head, ViewType1D next) + : + head_(head), + next_(next) + {} + + INLINE_FUNCTION_HD + Cells cellsSize()const { + return Cells(head_.extent(0), head_.extent(1), head_.extent(2));} + + INLINE_FUNCTION_HD + int32 start(IndexType i, IndexType j, IndexType k)const { + return head_(i,j,k); } + + INLINE_FUNCTION_HD + int32 getNext(int32 n)const { + if(n<0) return n; + return next_(n); } + }; + +protected: + + int32 capacity_ = 1; + + ViewType3D head_; + + ViewType1D next_; + + // borrowed ownership + ViewType1D pointPosition_; + + + INLINE_FUNCTION_H + void nullify() + { + fill( + head_, + range(0,this->nx()), + range(0,this->ny()), + range(0,this->nz()), + static_cast(-1) + ); + + fill( + next_, + range(0,capacity_), + static_cast(-1) + ); + } + + void nullify(range nextRng) + { + fill( + head_, + range(0,this->nx()), + range(0,this->ny()), + range(0,this->nz()), + static_cast(-1) + ); + + fill( + next_, + nextRng, + static_cast(-1) + ); + } + + using rangePolicyType = + Kokkos::RangePolicy< + Kokkos::IndexType, + Kokkos::Schedule, + execution_space>; + +private: + + void checkAllocateNext(int newCap) + { + if( capacity_ < newCap) + { + capacity_ = newCap; + reallocNoInit(next_, capacity_); + } + } + + void allocateHead() + { + reallocNoInit(head_, this->nx(), this->ny(), this->nz()); + } + +public: + + TypeNameNV("mapperNBS"); + + mapperNBS( + const box& domain, + real cellSize, + const ViewType1D& position) + : + Cells(domain, cellSize), + pointPosition_(position), + head_( + "mapperNBS::head_", + this->nx(), + this->ny(), + this->nz() + ), + next_("mapperNBS::next_",1) //,position.size()), + { + checkAllocateNext(pointPosition_.size()); + } + + mapperNBS( + const box& domain, + int32 nx, + int32 ny, + int32 nz, + const ViewType1D& position) + : + Cells(domain, nx, ny, nz), + pointPosition_(position), + head_("mapperNBS::head_",nx,ny,nz), //, this->nx(), this->ny(), this->nz()), + next_("mapperNBS::next_",1) //,position.size()), + { + checkAllocateNext(pointPosition_.size()); + } + + + INLINE_FUNCTION_HD + mapperNBS(const mapperNBS&) = default; + + INLINE_FUNCTION_HD + mapperNBS& operator = (const mapperNBS&) = default; + + INLINE_FUNCTION_HD + ~mapperNBS()=default; + + //// - Methods + INLINE_FUNCTION_HD + auto capacity()const + { + return capacity_; + } + + cellIterator getCellIterator()const + { + return cellIterator(head_, next_); + } + + bool particlesSizeChanged(int32 newSize) + { + checkAllocateNext(newSize); + return true; + } + + // - build based on all points in active range + INLINE_FUNCTION_H + void build(range activeRange) + { + checkAllocateNext(activeRange.second); + nullify(activeRange); + + Cells cellIndex = static_cast(*this); + auto points = pointPosition_; + auto next = next_; + auto head = head_; + + rangePolicyType rPolicy(activeRange.first, activeRange.second); + + Kokkos::parallel_for( + "mapperNBS::build", + rPolicy, + LAMBDA_HD(int32 i){ + CellType ind = cellIndex.pointIndex(points[i]); + int32 old = Kokkos::atomic_exchange(&head(ind.x(), ind.y(), ind.z()), i); + next[i] = old; + }); + Kokkos::fence(); + } + + + template + INLINE_FUNCTION_H + void build(range activeRange, IncludeFunction incld) + { + checkAllocateNext(activeRange.second); + nullify(activeRange); + + Cells cellIndex = static_cast(*this); + auto points = pointPosition_; + auto next = next_; + auto head = head_; + + rangePolicyType rPolicy(activeRange.first, activeRange.second); + + Kokkos::parallel_for( + "mapperNBS::build_Include", + rPolicy, + LAMBDA_HD(int32 i){ + if( incld(i) ) + { + CellType ind = cellIndex.pointIndex(points[i]); + auto old = Kokkos::atomic_exchange(&head(ind.x(), ind.y(), ind.z()), i); + next[i] = old; + } + }); + Kokkos::fence(); + + } + + + INLINE_FUNCTION_H + void buildCheckInDomain(range activeRange) + { + checkAllocateNext(activeRange.second); + nullify(activeRange); + + Cells cellIndex = static_cast(*this); + auto points = pointPosition_; + auto next = next_; + auto head = head_; + + rangePolicyType rPolicy(activeRange.first, activeRange.second); + + Kokkos::parallel_for( + "mapperNBS::buildCheckInDomain", + rPolicy, + LAMBDA_HD(int32 i){ + CellType ind; + if( cellIndex.pointIndexInDomain(points[i], ind) ) + { + int32 old = Kokkos::atomic_exchange(&head(ind.x(), ind.y(), ind.z()), i); + next[i] = old; + } + }); + + Kokkos::fence(); + + } + + template + INLINE_FUNCTION_H + void buildCheckInDomain(range activeRange, IncludeFunction incld) + { + checkAllocateNext(activeRange.second); + nullify(activeRange); + + Cells cellIndex = static_cast(*this); + auto points = pointPosition_; + auto next = next_; + auto head = head_; + + rangePolicyType rPolicy(activeRange.first, activeRange.second); + + Kokkos::parallel_for( + "mapperNBS::buildCheckInDomain_Include", + rPolicy, + LAMBDA_HD(int32 i){ + CellType ind; + if( incld(i) && cellIndex.pointIndexInDomain(points[i], ind) ) + { + auto old = Kokkos::atomic_exchange(&head(ind.x(), ind.y(), ind.z()), i); + next[i] = old; + } + }); + Kokkos::fence(); + } + +}; + +} // pFlow + +#endif // __mapperNBS_H__ \ No newline at end of file diff --git a/src/Interaction/sphereInteraction/sphereInteraction.H b/src/Interaction/sphereInteraction/sphereInteraction.H index 950276e2..db085645 100644 --- a/src/Interaction/sphereInteraction/sphereInteraction.H +++ b/src/Interaction/sphereInteraction/sphereInteraction.H @@ -136,6 +136,7 @@ public: { //Info<<"before contact search"<contactSearch_) { diff --git a/src/phasicFlow/Kokkos/KokkosUtilities.H b/src/phasicFlow/Kokkos/KokkosUtilities.H index 972565f5..05d11a9e 100644 --- a/src/phasicFlow/Kokkos/KokkosUtilities.H +++ b/src/phasicFlow/Kokkos/KokkosUtilities.H @@ -78,6 +78,42 @@ void reallocFill( ViewType1D& view, int32 len, Type val) } +template < + typename Type, + typename... Properties> +INLINE_FUNCTION_H +void realloc( ViewType3D& view, int32 len1, int32 len2, int32 len3) +{ + Kokkos::realloc(view, len1, len2, len3); +} + +template < + typename Type, + typename... Properties> +INLINE_FUNCTION_H +void reallocNoInit(ViewType3D& view, int32 len1, int32 len2, int32 len3) +{ + using ViewType = ViewType3D; + word vl = view.label(); + view = ViewType(); // Deallocate first + view = ViewType( + Kokkos::view_alloc( + Kokkos::WithoutInitializing, + vl), + len1, len2, len3); +} + +template < + typename Type, + typename... Properties> +INLINE_FUNCTION_H +void reallocFill( ViewType3D& view, int32 len1, int32 len2, int32 len3, Type val) +{ + reallocNoInit(view, len1, len2, len3); + Kokkos::deep_copy(view, val); +} + + template INLINE_FUNCTION_H void swapViews(ViewType& v1, ViewType &v2) diff --git a/utilities/particlesPhasicFlow/positionRandom/positionRandom.C b/utilities/particlesPhasicFlow/positionRandom/positionRandom.C index 11ea40dd..52cbccb4 100755 --- a/utilities/particlesPhasicFlow/positionRandom/positionRandom.C +++ b/utilities/particlesPhasicFlow/positionRandom/positionRandom.C @@ -21,7 +21,7 @@ Licence: #include "positionRandom.H" #include "uniformRandomReal.H" -#include "NBS.H" +#include "NBSLevel0.H" #include "unsortedPairs.H" #include "box.H" @@ -30,7 +30,7 @@ Licence: namespace pFlow { -using SearchType = NBS ; +using SearchType = NBSLevel0 ; using ContainerType = unsortedPairs; @@ -80,7 +80,7 @@ bool pFlow::positionRandom::positionOnePass(int32 pass, int32 startNum) fillPoints(startNum, positionHD, flagHD); - search.broadSearch(pairs, range(0, startNum), true); + search.broadSearch(pairs, range(0, startNum)); int32 numCollisions = findCollisions(pairs, flagHD); diff --git a/utilities/postprocessPhasicFlow/fieldOperations.H b/utilities/postprocessPhasicFlow/fieldOperations.H index ad4f26e4..8dbc1191 100644 --- a/utilities/postprocessPhasicFlow/fieldOperations.H +++ b/utilities/postprocessPhasicFlow/fieldOperations.H @@ -35,6 +35,7 @@ rectMeshField_H sumOp( const pointField_H field, const pointRectCell& poin { // create field const auto& mesh = pointToCell.mesh(); + auto iterator = pointToCell.getCellIterator(); rectMeshField_H results(mesh, T(0)); @@ -44,12 +45,12 @@ rectMeshField_H sumOp( const pointField_H field, const pointRectCell& poin { for(int32 k=0; k-1) { res += field[n]; - n = pointToCell.Next(n); + n = iterator.getNext(n); } results(i,j,k) = res; @@ -65,6 +66,7 @@ rectMeshField_H sumMaksOp( const pointField_H field, const pointRectCell& { // create field const auto& mesh = pointToCell.mesh(); + auto iterator = pointToCell.getCellIterator(); rectMeshField_H results(mesh, T(0)); @@ -75,7 +77,7 @@ rectMeshField_H sumMaksOp( const pointField_H field, const pointRectCell& for(int32 k=0; k-1) @@ -86,7 +88,7 @@ rectMeshField_H sumMaksOp( const pointField_H field, const pointRectCell& res += field[n]; } - n = pointToCell.Next(n); + n = iterator.getNext(n); } results(i,j,k) = res; diff --git a/utilities/postprocessPhasicFlow/pointRectCell.H b/utilities/postprocessPhasicFlow/pointRectCell.H index d4a9a324..dc709505 100644 --- a/utilities/postprocessPhasicFlow/pointRectCell.H +++ b/utilities/postprocessPhasicFlow/pointRectCell.H @@ -21,7 +21,7 @@ Licence: #ifndef __pointRectCell_H__ #define __pointRectCell_H__ -#include "NBS.H" +#include "mapperNBS.H" #include "rectMeshFields.H" #include "pointStructure.H" @@ -33,7 +33,7 @@ class pointRectCell { public: - using mapType = NBS; + using mapType = mapperNBS; using memory_space = typename mapType::memory_space; @@ -48,9 +48,6 @@ protected: ViewType1D pointPosition_; - ViewType1D diameter_; - - mapType map_; int32RectMeshField_H nPointInCell_; @@ -79,14 +76,12 @@ public: ), pStruct_(pStruct), pointPosition_(pStruct_.pointPosition().hostVectorAll()), - diameter_("diameter", pStruct_.capacity()), map_( mesh_.domain(), mesh_.nx(), mesh_.ny(), mesh_.nz(), - pointPosition_, - diameter_), + pointPosition_), nPointInCell_(mesh_, 0) { @@ -111,9 +106,9 @@ public: map_.buildCheckInDomain(activeRange, activeMask); + auto iterator = map_.getCellIterator(); - const auto& Next = map_.Next(); - const auto& Head = map_.Head(); + for(int32 i=0; i-1) { res+=1; - n = Next[n]; + n = iterator.getNext(n); } nPointInCell_(i,j,k) = res; @@ -136,20 +131,16 @@ public: } + auto getCellIterator()const + { + return map_.getCellIterator(); + } + int32 nPointInCell(int32 i, int32 j, int32 k)const { return nPointInCell_(i,j,k); } - auto inline Head(int32 i, int32 j, int32 k)const - { - return map_.Head()(i,j,k); - } - - auto inline Next(int32 n)const - { - return map_.Next()(n); - } //auto };