From c47d5f5cfa9708dc4f0a05210d392bf30932fd23 Mon Sep 17 00:00:00 2001 From: hamidrezanorouzi Date: Thu, 27 Oct 2022 20:10:12 +0330 Subject: [PATCH 1/2] multigrid-step1 --- .../contactSearch/methods/NBSLevel.H | 35 ++- .../contactSearch/methods/NBSLevel0.H | 7 + .../contactSearch/methods/NBSLevels.H | 280 ++++++++++++++++++ 3 files changed, 315 insertions(+), 7 deletions(-) create mode 100644 src/Interaction/contactSearch/methods/NBSLevels.H diff --git a/src/Interaction/contactSearch/methods/NBSLevel.H b/src/Interaction/contactSearch/methods/NBSLevel.H index 3f3689c2..657acf29 100644 --- a/src/Interaction/contactSearch/methods/NBSLevel.H +++ b/src/Interaction/contactSearch/methods/NBSLevel.H @@ -23,19 +23,22 @@ public: protected: - int32x3 gridExtent_; + int32x3 numCells_{1,1,1}; ViewType3D head_; - int8 level_ = 1; + int8 level_ = 0; public: - NBSLevel(int8 lvl, int32x3 gridExtent) + NBSLevel() + {} + + NBSLevel(int8 lvl, int32x3 numCells) : - gridExtent_(gridExtent), - head_("NBSLevel::head", gridExtent.x(), gridExtent.y(), gridExtent.z()), + numCells_(gridExtent), + head_("NBSLevel::head", numCells_.x(), numCells_.y(), numCells_.z()), level_(lvl) {} @@ -58,13 +61,31 @@ public: } INLINE_FUNCION_HD - const auto& gridExtent()const + const auto& numCells()const { return gridExtent_; } - + + void nullify() + { + fill( + head_, + range(0,numCells_.x()), + range(0,numCells_.y()), + range(0,numCells_.z()), + static_cast(-1) + ); + } }; +INLINE_FUNCION_HD +int32x3 mapIndexLevels( int32x3 ind, int32 lowerLevel, int32 upperLevel) +{ + int32 a = pow(2, static_cast(upperLevel-lowerLevel)); + return ind/a; +} +} + #endif diff --git a/src/Interaction/contactSearch/methods/NBSLevel0.H b/src/Interaction/contactSearch/methods/NBSLevel0.H index ea7dc275..43f1a970 100644 --- a/src/Interaction/contactSearch/methods/NBSLevel0.H +++ b/src/Interaction/contactSearch/methods/NBSLevel0.H @@ -124,6 +124,13 @@ public: INLINE_FUNCTION_HD ~NBSLevel0()=default; + + INLINE_FUNCTION_HD + auto sizeRatio()const + { + return sizeRatio_; + } + // - Perform the broad search to find pairs // with force = true, perform broad search regardless of // updateFrequency_ value diff --git a/src/Interaction/contactSearch/methods/NBSLevels.H b/src/Interaction/contactSearch/methods/NBSLevels.H new file mode 100644 index 00000000..c1174808 --- /dev/null +++ b/src/Interaction/contactSearch/methods/NBSLevels.H @@ -0,0 +1,280 @@ +#ifndef __NBSLevels_H__ +#define __NBSLevels_H__ + + + +namespace pFlow +{ + +template +class NBSLevels +: + public NBSLevel0 +{ + +public: + + using NBSLevel0Type = NBSLevel0; + + using IdType = typename NBSLevel0Type::IdType; + + using IndexType = typename NBSLevel0Type::IndexType; + + using Cells = typename NBSLevel0Type::Cells; + + using CellType = typename Cells::CellType; + + using execution_space = typename NBSLevel0Type::execution_space; + + using memory_space = typename NBSLevel0Type::memory_space; + + using realRange = kPair; + +protected: + + + real minSize_; + + real maxSize_; + + int32 numLevles_=1; + + ViewType1D nbsSuperLevels_; + + ViewType1D sizeRangeLevels_; + + ViewTypeID particleLevel_; + + range activeRange_; + + using rangePolicyType = + Kokkos::RangePolicy< + Kokkos::IndexType, + Kokkos::Schedule, + execution_space>; + + bool setNumLevels() + { + + int32 maxOvermin = static_cast(maxSize_/minSize_); + + if (maxOvermin <=1) + numLevles_ = 1; + else if(maxOvermin<=3) + numLevles_ = 2; + else if(maxOvermin<=7) + numLevles_ = 3; + else if(maxOvermin<15) + numLevles_ = 4; + else if(maxOvermin<31) + numLevles_ = 5; + else if(maxOvermin<63) + numLevles_ = 6; + else if(maxOvermin <127) + numLevles_ = 7; + else + { + fatalErrorInFunction<< + "size ratio is not valid for multi-grid NBS "<< maxOvermin<sizeRatio(); + real lvl_minD = 0.1*minSize_; + real lvl_maxD = sizeRatio* minSize_; + reallocNoInit(sizeRangeLevels_, numLevles_); + + for(int32 i=0; idomain(), sizeRangeLevels_[lvl].second); + nbsSuperLevels_[lvl] = NBSLevle(lvl, lvlCells.numCells()); + } + + return true; + } + + INLINE_FUNCTION_H + void nullify() + { + // level 0 + NBSLevel0Type::nullify(activeRange_); + + // super levels + for(int32 lvl=1; lvl& position, + const ViewType1D& diam), + range active + : + NBSLevel0Type( + domain, + sizeRatio*minSize, + position, + diam), + minSize_(minSize), + maxSize_(maxSize), + activeRange_(active) + { + if(!setNumLevels()) + { + fatalExit; + } + + setDiameterRange(); + + initSuperLevels(); + + findParticleLevel(activeRange_.first, activeRange_.last); + } + + + INLINE_FUNCTION_H + void build(range activeRange) + { + + // check for active ragne************************************************************************ + checkAllocateNext(activeRange.second); + nullify(); + // + + rangePolicyType rPolicy(activeRange.first, activeRange.second); + + Kokkos::parallel_for( + "NBSLevels::build", + rPolicy, + CLASS_LAMBDA_HD(int32 i){ + CellType ind = this->pointIndex(points[i]); + int8 lvl = particleLevel_[i]; + int32 old; + if( lvl == 0) + { + int32 old = + Kokkos::atomic_exchange( + &this->head_(ind.x(), ind.y(), ind.z()), + i); + this->next_[i] = old; + } + else + { + ind = upper.mapIndexUpperLevel(ind, 0, lvl); + int32 old = + Kokkos::atomic_exchange( + &nbsSuperLevels_[lvl].head()(ind.x(), ind.y(), ind.z()), + i); + this->next_[i] = old; + } + }); + + Kokkos::fence(); + } + + template + INLINE_FUNCTION_H + void build(range activeRange, IncludeFunction incld) + { + // check for active ragne************************************************************************ + checkAllocateNext(activeRange.second); + nullify(); + // + + rangePolicyType rPolicy(activeRange.first, activeRange.second); + + Kokkos::parallel_for( + "NBSLevels::build", + rPolicy, + CLASS_LAMBDA_HD(int32 i){ + if(incld(i)) + { + CellType ind = this->pointIndex(points[i]); + int8 lvl = particleLevel_[i]; + int32 old; + if( lvl == 0) + { + int32 old = + Kokkos::atomic_exchange( + &this->head_(ind.x(), ind.y(), ind.z()), + i); + this->next_[i] = old; + } + else + { + ind = upper.mapIndexUpperLevel(ind, 0, lvl); + int32 old = + Kokkos::atomic_exchange( + &nbsSuperLevels_[lvl].head()(ind.x(), ind.y(), ind.z()), + i); + this->next_[i] = old; + } + } + }); + + Kokkos::fence(); + + } + + + bool findParticleLevel(int32 first, int32 last) + { + auto diameter = this->diameter_; + auto sizeRange = sizeRangeLevels_; + auto particleLevel = particleLevel_; + + int8 maxLvl = sizeRangeLevels_.size(); + + rangePolicyType rPolicy(first, last); + + Kokkos::parallel_for( + "NBSLevels::findParticleLevel", + rPolicy, + LAMBDA_HD(int32 i) + { + for(int8 lvl = 0; lvl Date: Sun, 30 Oct 2022 18:08:32 +0330 Subject: [PATCH 2/2] multigrid NBS is added, not tested with particle insertion --- .../ContactSearch/ContactSearch.H | 15 +- .../ContactSearch/ContactSearchs.C | 8 +- src/Interaction/contactSearch/methods/NBS.H | 22 +- .../contactSearch/methods/NBSCrossLoop.H | 82 +++ .../contactSearch/methods/NBSLevel.H | 124 +++-- .../contactSearch/methods/NBSLevel0.H | 68 ++- .../contactSearch/methods/NBSLevels.H | 422 +++++++++----- .../contactSearch/methods/NBSold.H | 523 ------------------ .../contactSearch/methods/mapperNBS.H | 125 +++-- .../contactSearch/methods/multiGridNBS.H | 213 +++++++ .../contactSearch/wallMappings/cellMapping.H | 150 +++++ .../{cellsSimple.H => cellsWallLevel0.H} | 154 ++---- .../wallMappings/cellsWallLevels.H | 152 +++++ .../wallMappings/multiGridMapping.H | 153 +++++ src/phasicFlow/Kokkos/KokkosTypes.H | 3 + src/phasicFlow/structuredData/box/box.C | 10 - src/phasicFlow/structuredData/box/box.H | 16 +- 17 files changed, 1344 insertions(+), 896 deletions(-) create mode 100644 src/Interaction/contactSearch/methods/NBSCrossLoop.H delete mode 100644 src/Interaction/contactSearch/methods/NBSold.H create mode 100644 src/Interaction/contactSearch/methods/multiGridNBS.H create mode 100644 src/Interaction/contactSearch/wallMappings/cellMapping.H rename src/Interaction/contactSearch/wallMappings/{cellsSimple.H => cellsWallLevel0.H} (66%) create mode 100644 src/Interaction/contactSearch/wallMappings/cellsWallLevels.H create mode 100644 src/Interaction/contactSearch/wallMappings/multiGridMapping.H diff --git a/src/Interaction/contactSearch/ContactSearch/ContactSearch.H b/src/Interaction/contactSearch/ContactSearch/ContactSearch.H index 4220b820..21f03934 100644 --- a/src/Interaction/contactSearch/ContactSearch/ContactSearch.H +++ b/src/Interaction/contactSearch/ContactSearch/ContactSearch.H @@ -31,7 +31,7 @@ namespace pFlow template< template class BaseMethod, - template class WallMapping + template class WallMapping > class ContactSearch : @@ -53,10 +53,7 @@ public: using WallMappingType = WallMapping< - ExecutionSpace, - IdType, - IndexType - >; + ExecutionSpace>; protected: @@ -83,7 +80,7 @@ public: auto method = dict().getVal("method"); auto wmMethod = dict().getVal("wallMapping"); - auto nbDict = dict().subDictOrCreate(method+"Info"); + auto nbDict = dict().subDict(method+"Info"); real minD, maxD; this->Particles().boundingSphereMinMax(minD, maxD); @@ -96,6 +93,7 @@ public: ( nbDict, this->domain(), + minD, maxD, position, diam @@ -104,7 +102,7 @@ public: greenText(particleContactSearch_().typeName())<Geometry().numPoints(); int32 wnTri = this->Geometry().size(); @@ -115,7 +113,8 @@ public: wallMapping_ = makeUnique( wmDict, - particleContactSearch_().getCells(), + particleContactSearch_().numLevels(), + particleContactSearch_().getCellsLevels(), wnPoints, wnTri, wPoints, diff --git a/src/Interaction/contactSearch/ContactSearch/ContactSearchs.C b/src/Interaction/contactSearch/ContactSearch/ContactSearchs.C index 75343dfa..7d617155 100644 --- a/src/Interaction/contactSearch/ContactSearch/ContactSearchs.C +++ b/src/Interaction/contactSearch/ContactSearch/ContactSearchs.C @@ -20,7 +20,11 @@ Licence: #include "ContactSearch.H" -#include "cellsSimple.H" +#include "cellMapping.H" #include "NBS.H" +#include "multiGridNBS.H" +#include "multiGridMapping.H" -template class pFlow::ContactSearch; + +template class pFlow::ContactSearch; +template class pFlow::ContactSearch; diff --git a/src/Interaction/contactSearch/methods/NBS.H b/src/Interaction/contactSearch/methods/NBS.H index a3610066..4fe46272 100644 --- a/src/Interaction/contactSearch/methods/NBS.H +++ b/src/Interaction/contactSearch/methods/NBS.H @@ -35,6 +35,8 @@ public: using NBSLevel0Type = NBSLevel0; + using cellIterator = typename NBSLevel0Type::cellIterator; + using IdType = typename NBSLevel0Type::IdType; using IndexType = typename NBSLevel0Type::IndexType; @@ -81,8 +83,9 @@ public: TypeNameNV("NBS"); NBS( - dictionary dict, + const dictionary& dict, const box& domain, + real minSize, real maxSize, const ViewType1D& position, const ViewType1D& diam) @@ -126,11 +129,26 @@ public: return performedSearch_; } - auto getCellIterator()const + Vector getCellIteratorLevels() + { + return Vector("cellIterator", 1, NBSLevel0_.getCellIterator()); + } + + auto getCellIterator(int32 lvl)const { return NBSLevel0_.getCellIterator(); } + int32 numLevels()const + { + return 1; + } + + Vector getCellsLevels()const + { + return Vector("Cells", 1, NBSLevel0_.getCells()); + } + auto getCells()const { return NBSLevel0_.getCells(); diff --git a/src/Interaction/contactSearch/methods/NBSCrossLoop.H b/src/Interaction/contactSearch/methods/NBSCrossLoop.H new file mode 100644 index 00000000..adce7409 --- /dev/null +++ b/src/Interaction/contactSearch/methods/NBSCrossLoop.H @@ -0,0 +1,82 @@ +/*------------------------------- 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. + +-----------------------------------------------------------------------------*/ + + +int32 m = this->head_(i,j,k); +CellType currentCell(i,j,k); +int32 n = -1; + +while( m > -1 ) +{ + + auto p_m = this->pointPosition_[m]; + auto d_m = this->sizeRatio_* this->diameter_[m]; + + int32x3 crossIndex = mapIndexLevels( + int32x3(i,j,k), + level_, + upperLevel.level()); + + + for(int32 ii = -1; ii<2; ii++) + { + for(int32 jj=-1; jj<2; jj++) + { + int32 kk=-1; + while( kk < 2) + { + int32x3 nghbrCI = crossIndex + int32x3(ii,jj,kk); + + if( upperLevel.isInRange(nghbrCI) ) + { + n = upperLevel.head_( + nghbrCI.x(), + nghbrCI.y(), + nghbrCI.z()); + + while( n >-1) + { + auto p_n = this->pointPosition_[n]; + auto d_n = this->sizeRatio_*this->diameter_[n]; + + if( sphereSphereCheck(p_m, p_n, d_m, d_n) ) + { + auto ln = n; + auto lm = m; + if(lm>ln) this->Swap(lm,ln); + if( auto res = pairs.insert(lm,ln); res <0) + { + getFullUpdate++; + } + } + + n = this->next_[n]; + } + + } + + kk++; + } + } + } + + m = this->next_[m]; +} + diff --git a/src/Interaction/contactSearch/methods/NBSLevel.H b/src/Interaction/contactSearch/methods/NBSLevel.H index 657acf29..2a1c6aaa 100644 --- a/src/Interaction/contactSearch/methods/NBSLevel.H +++ b/src/Interaction/contactSearch/methods/NBSLevel.H @@ -2,84 +2,120 @@ #define __NBSLevel_H__ -#include "cells.H" -#include "contactSearchFunctions.H" +#include "NBSLevel0.H" + + +namespace pFlow +{ + +INLINE_FUNCTION_HD +int32x3 mapIndexLevels(const int32x3& ind, int32 lowerLevel, int32 upperLevel); template class NBSLevel +: + public NBSLevel0 { public: + using NBSLevel0Type = NBSLevel0; - using execution_space= executionSpace; + using cellIterator = typename NBSLevel0Type::cellIterator; - using memory_space = typename execution_space::memory_space; + using IdType = typename NBSLevel0Type::IdType; + + using IndexType = typename NBSLevel0Type::IndexType; - using rangePolicy = Kokkos::RangePolicy< - Kokkos::IndexType, - Kokkos::Schedule, - execution_space> + using Cells = typename NBSLevel0Type::Cells; + + using CellType = typename Cells::CellType; + + using execution_space = typename NBSLevel0Type::execution_space; + + using memory_space = typename NBSLevel0Type::memory_space; + + using mdrPolicyFindPairs = typename NBSLevel0Type::mdrPolicyFindPairs; + + using HeadType = typename NBSLevel0Type::HeadType; + + using NextType = typename NBSLevel0Type::NextType; + + template + friend class NBSLevels; protected: - int32x3 numCells_{1,1,1}; - - ViewType3D head_; - - int8 level_ = 0; + int32 level_ = 0; public: - NBSLevel() - {} + TypeNameNV("NBSLevel0"); - NBSLevel(int8 lvl, int32x3 numCells) + INLINE_FUNCTION_HD + NBSLevel(){} + + NBSLevel( + int32 lvl, + const box& domain, + real cellSize, + real sizeRatio, + const ViewType1D& position, + const ViewType1D& diam) : - numCells_(gridExtent), - head_("NBSLevel::head", numCells_.x(), numCells_.y(), numCells_.z()), + NBSLevel0Type( + domain, + cellSize, + sizeRatio, + position, + diam, + lvl==0), level_(lvl) {} - INLINE_FUNCION_HD - auto& head(int32 i, int32 j, int32 k) - { - return head_(i,j,k); - } + INLINE_FUNCTION_HD + NBSLevel(const NBSLevel&) = default; + + INLINE_FUNCTION_HD + NBSLevel& operator = (const NBSLevel&) = default; - INLINE_FUNCION_HD - auto& head() - { - return head_; - } + INLINE_FUNCTION_HD + ~NBSLevel() = default; - INLINE_FUNCION_HD + INLINE_FUNCTION_HD auto level()const { return level_; } - INLINE_FUNCION_HD - const auto& numCells()const + template + INLINE_FUNCTION_H + int32 findPairsCountCross(PairsContainer& pairs, NBSLevel& upperLevel) { - return gridExtent_; + + mdrPolicyFindPairs + mdrPolicy( + {0,0,0}, + {this->nx(),this->ny(),this->nz()} ); + + int32 notInsertedPairs; + + Kokkos::parallel_reduce ( + "NBSLevel::findPairsCountCross", + mdrPolicy, + CLASS_LAMBDA_HD(int32 i, int32 j, int32 k, int32& getFullUpdate){ + #include "NBSCrossLoop.H" + }, notInsertedPairs); + + return notInsertedPairs; } - void nullify() - { - fill( - head_, - range(0,numCells_.x()), - range(0,numCells_.y()), - range(0,numCells_.z()), - static_cast(-1) - ); - } + }; -INLINE_FUNCION_HD -int32x3 mapIndexLevels( int32x3 ind, int32 lowerLevel, int32 upperLevel) +INLINE_FUNCTION_HD +int32x3 mapIndexLevels( const int32x3& ind, int32 lowerLevel, int32 upperLevel) { int32 a = pow(2, static_cast(upperLevel-lowerLevel)); return ind/a; diff --git a/src/Interaction/contactSearch/methods/NBSLevel0.H b/src/Interaction/contactSearch/methods/NBSLevel0.H index 43f1a970..af313402 100644 --- a/src/Interaction/contactSearch/methods/NBSLevel0.H +++ b/src/Interaction/contactSearch/methods/NBSLevel0.H @@ -35,20 +35,26 @@ class NBSLevel0 { public: - using MapperType = mapperNBS; + using MapperType = mapperNBS; + + using cellIterator = typename MapperType::cellIterator; using IdType = typename MapperType::IdType; - using IndexType = typename MapperType::IndexType; + using IndexType = typename MapperType::IndexType; using Cells = typename MapperType::Cells; - using CellType = typename Cells::CellType; + using CellType = typename Cells::CellType; - using execution_space = typename MapperType::execution_space; + using execution_space = typename MapperType::execution_space; using memory_space = typename MapperType::memory_space; + using HeadType = typename MapperType::HeadType; + + using NextType = typename MapperType::NextType; + struct TagFindPairs{}; @@ -66,9 +72,6 @@ protected: Kokkos::Schedule, execution_space>; -private: - - static INLINE_FUNCTION_HD void Swap(int32& x, int32& y) { @@ -81,6 +84,9 @@ public: TypeNameNV("NBSLevel0"); + INLINE_FUNCTION_HD + NBSLevel0(){} + NBSLevel0( const box& domain, real cellSize, @@ -108,9 +114,10 @@ public: real cellSize, real sizeRatio, const ViewType1D& position, - const ViewType1D& diam) + const ViewType1D& diam, + bool nextOwner = true) : - MapperType(domain, cellSize, position), + MapperType(domain, cellSize, position, nextOwner), sizeRatio_(sizeRatio), diameter_(diam) {} @@ -131,6 +138,12 @@ public: return sizeRatio_; } + INLINE_FUNCTION_HD + auto& diameter() + { + return diameter_; + } + // - Perform the broad search to find pairs // with force = true, perform broad search regardless of // updateFrequency_ value @@ -166,10 +179,7 @@ public: INLINE_FUNCTION_H bool findPairs(PairsContainer& pairs) { - mdrPolicyFindPairs - mdrPolicy( - {0,0,0}, - {this->nx(),this->ny(),this->nz()} ); + int32 getFull = 1; @@ -178,15 +188,7 @@ public: 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); - + getFull = findPairsCount(pairs); if(getFull) { @@ -209,7 +211,27 @@ public: return true; } - + template + INLINE_FUNCTION_H + int32 findPairsCount(PairsContainer& pairs) + { + mdrPolicyFindPairs + mdrPolicy( + {0,0,0}, + {this->nx(),this->ny(),this->nz()} ); + + int32 notInsertedPairs; + + Kokkos::parallel_reduce ( + "NBSLevel0::findPairs", + mdrPolicy, + CLASS_LAMBDA_HD(int32 i, int32 j, int32 k, int32& getFullUpdate){ + #include "NBSLoop.H" + }, notInsertedPairs); + + return notInsertedPairs; + + } }; diff --git a/src/Interaction/contactSearch/methods/NBSLevels.H b/src/Interaction/contactSearch/methods/NBSLevels.H index c1174808..d945ced4 100644 --- a/src/Interaction/contactSearch/methods/NBSLevels.H +++ b/src/Interaction/contactSearch/methods/NBSLevels.H @@ -1,51 +1,58 @@ #ifndef __NBSLevels_H__ #define __NBSLevels_H__ - +#include "NBSLevel.H" +#include "NBSLevel0.H" +#include "KokkosTypes.H" namespace pFlow { template class NBSLevels -: - public NBSLevel0 { public: - using NBSLevel0Type = NBSLevel0; + using NBSLevelType = NBSLevel; - using IdType = typename NBSLevel0Type::IdType; + using cellIterator = typename NBSLevelType::cellIterator; + + using IdType = typename NBSLevelType::IdType; - using IndexType = typename NBSLevel0Type::IndexType; + using IndexType = typename NBSLevelType::IndexType; - using Cells = typename NBSLevel0Type::Cells; + using Cells = typename NBSLevelType::Cells; using CellType = typename Cells::CellType; - using execution_space = typename NBSLevel0Type::execution_space; + using execution_space = typename NBSLevelType::execution_space; - using memory_space = typename NBSLevel0Type::memory_space; + using memory_space = typename NBSLevelType::memory_space; using realRange = kPair; protected: - real minSize_; real maxSize_; - int32 numLevles_=1; + int32 numLevels_=1; - ViewType1D nbsSuperLevels_; + Vector nbsLevels_; - ViewType1D sizeRangeLevels_; + ViewType1D sizeRangeLevels_; - ViewTypeID particleLevel_; + ViewType1D sizeRangeLevelsHost_; - range activeRange_; + ViewType1D maxSizeLevels_; + + ViewType1D maxSizeLevelsHost_; + + ViewType1D particleLevel_; + + range activeRange_{0,0}; using rangePolicyType = Kokkos::RangePolicy< @@ -53,76 +60,115 @@ protected: Kokkos::Schedule, execution_space>; - bool setNumLevels() + int32 setNumLevels() { int32 maxOvermin = static_cast(maxSize_/minSize_); if (maxOvermin <=1) - numLevles_ = 1; + return 1; else if(maxOvermin<=3) - numLevles_ = 2; + return 2; else if(maxOvermin<=7) - numLevles_ = 3; + return 3; else if(maxOvermin<15) - numLevles_ = 4; + return 4; else if(maxOvermin<31) - numLevles_ = 5; + return 5; else if(maxOvermin<63) - numLevles_ = 6; + return 6; else if(maxOvermin <127) - numLevles_ = 7; + return 7; else { fatalErrorInFunction<< "size ratio is not valid for multi-grid NBS "<< maxOvermin<sizeRatio(); - real lvl_minD = 0.1*minSize_; - real lvl_maxD = sizeRatio* minSize_; - reallocNoInit(sizeRangeLevels_, numLevles_); + real lvl_maxD = sizeRatio* maxSize_; + real lvl_minD = lvl_maxD/2.0; - for(int32 i=0; i=0; lvl--) { - sizeRangeLevels_[i] = {lvl_minD, lvl_maxD}; - lvl_minD = lvl_maxD; - lvl_maxD *= 2.0; + + if(lvl == 0 ) lvl_minD = 0.01*minSize_; + + sizeRangeLevelsHost_[lvl] = {lvl_minD, lvl_maxD}; + maxSizeLevelsHost_[lvl] = lvl_maxD; + lvl_maxD = lvl_minD; + lvl_minD /= 2.0; } - - return true; - } - - bool initSuperLevels() - { - if(numLevles_== 1) return true; - reallocNoInit(nbsSuperLevels_, numLevles_); - - if(int lvl = 1; lvldomain(), sizeRangeLevels_[lvl].second); - nbsSuperLevels_[lvl] = NBSLevle(lvl, lvlCells.numCells()); - } - - return true; - } - - INLINE_FUNCTION_H - void nullify() - { - // level 0 - NBSLevel0Type::nullify(activeRange_); - // super levels - for(int32 lvl=1; lvl& position, + const ViewType1D& diam) + { + + + for(int32 lvl = 0; lvl nbsLevels_[0].capacity()) + { + nbsLevels_[0].checkAllocateNext(activeRange_.second); + + auto next0 = nbsLevels_[0].next(); + + for(int32 lvl=1; lvl& position, - const ViewType1D& diam), - range active + const ViewType1D& diam) : - NBSLevel0Type( - domain, - sizeRatio*minSize, - position, - diam), minSize_(minSize), maxSize_(maxSize), - activeRange_(active) + numLevels_(setNumLevels()), + nbsLevels_("nbsLevels", numLevels_, numLevels_, RESERVE()), + sizeRangeLevels_("sizeRangeLevels", numLevels_), + sizeRangeLevelsHost_("sizeRangeLevelsHost", numLevels_), + maxSizeLevels_("maxSizeLevels", numLevels_), + maxSizeLevelsHost_("maxSizeLevelsHost", numLevels_) { - if(!setNumLevels()) - { - fatalExit; - } - setDiameterRange(); + setDiameterRange(sizeRatio); - initSuperLevels(); - - findParticleLevel(activeRange_.first, activeRange_.last); + initLevels(domain, sizeRatio, position, diam); } + auto getCellIterator(int32 lvl)const + { + return nbsLevels_[lvl].getCellIterator(); + } + + int32 numLevels()const + { + return numLevels_; + } + + Cells getCells(int32 lvl)const + { + return nbsLevels_[lvl].getCells(); + } + + template + INLINE_FUNCTION_H + bool findPairs(PairsContainer& pairs) + { + + int32 getFull = 1; + + // loop until the container size fits the numebr of contact pairs + while (getFull > 0) + { + + getFull = findPairsCount(pairs); + + if(getFull) + { + // - resize the container + // note that getFull now shows the number of failed insertions. + uint32 len = max(getFull,100) ; + + auto oldCap = pairs.capacity(); + + pairs.increaseCapacityBy(len); + + Info<< "The contact pair container capacity increased from "<< + oldCap << " to "< + INLINE_FUNCTION_H + int32 findPairsCount(PairsContainer& pairs) + { + + int32 notInsertedCount = 0; + + for(int32 lvl=0; lvl1) head1 = nbsLevels_[1].head(); + if(numLevels_>2) head2 = nbsLevels_[2].head(); + if(numLevels_>3) head3 = nbsLevels_[3].head(); + if(numLevels_>4) head4 = nbsLevels_[4].head(); + if(numLevels_>5) head5 = nbsLevels_[5].head(); + if(numLevels_>6) head6 = nbsLevels_[6].head(); + + + Kokkos::parallel_for( "NBSLevels::build", rPolicy, - CLASS_LAMBDA_HD(int32 i){ - CellType ind = this->pointIndex(points[i]); - int8 lvl = particleLevel_[i]; + LAMBDA_HD(int32 i){ + + int8 lvl = particleLevel[i]; + auto ind = nbsLevel0.pointIndex(pointPosition[i]); + ind = mapIndexLevels(ind, 0, lvl); int32 old; - if( lvl == 0) - { - int32 old = - Kokkos::atomic_exchange( - &this->head_(ind.x(), ind.y(), ind.z()), - i); - this->next_[i] = old; - } - else - { - ind = upper.mapIndexUpperLevel(ind, 0, lvl); - int32 old = - Kokkos::atomic_exchange( - &nbsSuperLevels_[lvl].head()(ind.x(), ind.y(), ind.z()), - i); - this->next_[i] = old; - } + if(lvl==0) + old =Kokkos::atomic_exchange(&head0(ind.x(), ind.y(), ind.z()),i); + else if(lvl==1) + old =Kokkos::atomic_exchange(&head1(ind.x(), ind.y(), ind.z()),i); + else if(lvl==2) + old =Kokkos::atomic_exchange(&head2(ind.x(), ind.y(), ind.z()),i); + else if(lvl==3) + old =Kokkos::atomic_exchange(&head3(ind.x(), ind.y(), ind.z()),i); + else if(lvl==4) + old =Kokkos::atomic_exchange(&head4(ind.x(), ind.y(), ind.z()),i); + else if(lvl==5) + old =Kokkos::atomic_exchange(&head5(ind.x(), ind.y(), ind.z()),i); + else if(lvl==6) + old =Kokkos::atomic_exchange(&head6(ind.x(), ind.y(), ind.z()),i); + + next(i) = old; }); Kokkos::fence(); @@ -203,57 +335,79 @@ public: INLINE_FUNCTION_H void build(range activeRange, IncludeFunction incld) { - // check for active ragne************************************************************************ - checkAllocateNext(activeRange.second); - nullify(); - // + // nullify next and heads + findParticleLevel(activeRange.first, activeRange.second); + manageAllocateNext(activeRange); + nullify(activeRange); + rangePolicyType rPolicy(activeRange.first, activeRange.second); + auto nbsLevel0 = nbsLevels_[0]; + auto pointPosition = nbsLevel0.pointPosition(); + auto particleLevel = particleLevel_; + auto next = nbsLevel0.next(); + auto head0 = nbsLevel0.head(); + + typename NBSLevelType::HeadType head1, head2, head3, head4, head5, head6; + + if(numLevels_>1) head1 = nbsLevels_[1].head(); + if(numLevels_>2) head2 = nbsLevels_[2].head(); + if(numLevels_>3) head3 = nbsLevels_[3].head(); + if(numLevels_>4) head4 = nbsLevels_[4].head(); + if(numLevels_>5) head5 = nbsLevels_[5].head(); + if(numLevels_>6) head6 = nbsLevels_[6].head(); + Kokkos::parallel_for( "NBSLevels::build", rPolicy, - CLASS_LAMBDA_HD(int32 i){ - if(incld(i)) - { - CellType ind = this->pointIndex(points[i]); - int8 lvl = particleLevel_[i]; - int32 old; - if( lvl == 0) - { - int32 old = - Kokkos::atomic_exchange( - &this->head_(ind.x(), ind.y(), ind.z()), - i); - this->next_[i] = old; - } - else - { - ind = upper.mapIndexUpperLevel(ind, 0, lvl); - int32 old = - Kokkos::atomic_exchange( - &nbsSuperLevels_[lvl].head()(ind.x(), ind.y(), ind.z()), - i); - this->next_[i] = old; - } - } + LAMBDA_HD(int32 i){ + if(!incld(i)) return; + + int8 lvl = particleLevel[i]; + auto ind = nbsLevel0.pointIndex(pointPosition[i]); + + ind = mapIndexLevels(ind, 0, lvl); + int32 old; + if(lvl==0) + old =Kokkos::atomic_exchange(&head0(ind.x(), ind.y(), ind.z()),i); + else if(lvl==1) + old =Kokkos::atomic_exchange(&head1(ind.x(), ind.y(), ind.z()),i); + else if(lvl==2) + old =Kokkos::atomic_exchange(&head2(ind.x(), ind.y(), ind.z()),i); + else if(lvl==3) + old =Kokkos::atomic_exchange(&head3(ind.x(), ind.y(), ind.z()),i); + else if(lvl==4) + old =Kokkos::atomic_exchange(&head4(ind.x(), ind.y(), ind.z()),i); + else if(lvl==5) + old =Kokkos::atomic_exchange(&head5(ind.x(), ind.y(), ind.z()),i); + else if(lvl==6) + old =Kokkos::atomic_exchange(&head6(ind.x(), ind.y(), ind.z()),i); + + next(i) = old; + }); Kokkos::fence(); } - bool findParticleLevel(int32 first, int32 last) { - auto diameter = this->diameter_; - auto sizeRange = sizeRangeLevels_; + + if(last > particleLevel_.size()) + { + reallocNoInit(particleLevel_,last); + } + + auto diameter = nbsLevels_[0].diameter(); + auto maxSizes = maxSizeLevels_; auto particleLevel = particleLevel_; int8 maxLvl = sizeRangeLevels_.size(); rangePolicyType rPolicy(first, last); - + Kokkos::parallel_for( "NBSLevels::findParticleLevel", rPolicy, @@ -261,13 +415,15 @@ public: { for(int8 lvl = 0; lvl(-1); }); + Kokkos::fence(); return true; } diff --git a/src/Interaction/contactSearch/methods/NBSold.H b/src/Interaction/contactSearch/methods/NBSold.H deleted file mode 100644 index cf5aed45..00000000 --- a/src/Interaction/contactSearch/methods/NBSold.H +++ /dev/null @@ -1,523 +0,0 @@ -/*------------------------------- 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 "<; + + using NextType = ViewType1D; + class cellIterator { private: - ViewType3D head_; + HeadType head_; - ViewType1D next_; + NextType next_; public: @@ -85,12 +89,19 @@ protected: ViewType1D next_; + bool nextOwner_ = true; + // borrowed ownership ViewType1D pointPosition_; - + + using rangePolicyType = + Kokkos::RangePolicy< + Kokkos::IndexType, + Kokkos::Schedule, + execution_space>; INLINE_FUNCTION_H - void nullify() + void nullifyHead() { fill( head_, @@ -99,24 +110,11 @@ protected: range(0,this->nz()), static_cast(-1) ); - - fill( - next_, - range(0,capacity_), - static_cast(-1) - ); } - void nullify(range nextRng) + void nullifyNext(range nextRng) { - fill( - head_, - range(0,this->nx()), - range(0,this->ny()), - range(0,this->nz()), - static_cast(-1) - ); - + if(!nextOwner_)return; fill( next_, nextRng, @@ -124,19 +122,27 @@ protected: ); } - using rangePolicyType = - Kokkos::RangePolicy< - Kokkos::IndexType, - Kokkos::Schedule, - execution_space>; + void nullify() + { + nullifyHead(); -private: + nullifyNext(range(0,capacity_)); + } + + void nullify(range nextRng) + { + nullifyHead(); + + nullifyNext(nextRng); + } + void checkAllocateNext(int newCap) - { + { if( capacity_ < newCap) { capacity_ = newCap; + if(!nextOwner_)return; reallocNoInit(next_, capacity_); } } @@ -145,15 +151,21 @@ private: { reallocNoInit(head_, this->nx(), this->ny(), this->nz()); } + + public: TypeNameNV("mapperNBS"); + INLINE_FUNCTION_HD + mapperNBS(){} + mapperNBS( const box& domain, - real cellSize, - const ViewType1D& position) + real cellSize, + const ViewType1D& position, + bool nextOwner = true) : Cells(domain, cellSize), pointPosition_(position), @@ -163,7 +175,8 @@ public: this->ny(), this->nz() ), - next_("mapperNBS::next_",1) //,position.size()), + next_("mapperNBS::next_",1), //,position.size()), + nextOwner_(nextOwner) { checkAllocateNext(pointPosition_.size()); } @@ -173,12 +186,14 @@ public: int32 nx, int32 ny, int32 nz, - const ViewType1D& position) + const ViewType1D& position, + bool nextOwner = true) : 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()), + head_("mapperNBS::head_",nx,ny,nz), + next_("mapperNBS::next_",1), + nextOwner_(nextOwner) { checkAllocateNext(pointPosition_.size()); } @@ -205,12 +220,54 @@ public: return cellIterator(head_, next_); } - bool particlesSizeChanged(int32 newSize) + bool particlesCapcityChanged(int32 newCap) { - checkAllocateNext(newSize); + checkAllocateNext(newCap); return true; } + INLINE_FUNCTION_HD + auto& head() + { + return head_; + } + + INLINE_FUNCTION_HD + auto& next() + { + return next_; + } + + INLINE_FUNCTION_HD + const auto& head()const + { + return head_; + } + + INLINE_FUNCTION_HD + const auto& next()const + { + return next_; + } + + INLINE_FUNCTION_HD + auto& pointPosition() + { + return pointPosition_; + } + + INLINE_FUNCTION_H + void setNext(ViewType1D& next) + { + if(!nextOwner_) + { + next_ = next; + capacity_ = next.size(); + } + } + + + // - build based on all points in active range INLINE_FUNCTION_H void build(range activeRange) diff --git a/src/Interaction/contactSearch/methods/multiGridNBS.H b/src/Interaction/contactSearch/methods/multiGridNBS.H new file mode 100644 index 00000000..bdbf60dd --- /dev/null +++ b/src/Interaction/contactSearch/methods/multiGridNBS.H @@ -0,0 +1,213 @@ +/*------------------------------- 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 __multiGridNBS_H__ +#define __multiGridNBS_H__ + +#include "NBSLevels.H" + +namespace pFlow +{ + + +template +class multiGridNBS +{ +public: + + using NBSLevelsType = NBSLevels; + + using cellIterator = typename NBSLevelsType::cellIterator; + + using IdType = typename NBSLevelsType::IdType; + + using IndexType = typename NBSLevelsType::IndexType; + + using Cells = typename NBSLevelsType::Cells; + + using CellType = typename Cells::CellType; + + using execution_space = typename NBSLevelsType::execution_space; + + using memory_space = typename NBSLevelsType::memory_space; + + +protected: + + real sizeRatio_ = 1.0; + + int32 updateFrequency_= 1; + + int32 currentIter_ = 0; + + bool performedSearch_ = false; + + NBSLevelsType NBSLevels_; + +private: + + bool performSearch() + { + if(currentIter_ % updateFrequency_ == 0) + { + currentIter_++; + return true; + + }else + { + currentIter_++; + return false; + } + } + +public: + + TypeNameNV("multiGridNBS"); + + multiGridNBS( + const dictionary& dict, + const box& domain, + real minSize, + real maxSize, + const ViewType1D& position, + const ViewType1D& diam) + : + sizeRatio_( + max( + dict.getVal("sizeRatio"), + 1.0 + )), + updateFrequency_( + max( + dict.getVal("updateFrequency"), + 1 + )), + NBSLevels_( + domain, + minSize, + maxSize, + sizeRatio_, + position, + diam) + {} + + INLINE_FUNCTION_HD + multiGridNBS(const multiGridNBS&) = default; + + INLINE_FUNCTION_HD + multiGridNBS& operator = (const multiGridNBS&) = default; + + INLINE_FUNCTION_HD + ~multiGridNBS()=default; + + //// - Methods + + bool enterBoadSearch()const + { + return currentIter_%updateFrequency_==0; + } + + bool performedSearch()const + { + return performedSearch_; + } + + int32 numLevels()const + { + return NBSLevels_.numLevels(); + } + + auto getCellsLevels()const + { + Vector cellsLvl("cells", numLevels(), numLevels(), RESERVE()); + + for(int32 lvl=0; lvl + bool broadSearch(PairsContainer& pairs, range activeRange, bool force=false) + { + + if(force) currentIter_ = 0; + performedSearch_ = false; + + if( !performSearch() ) return true; + + + NBSLevels_.build(activeRange); + + NBSLevels_.findPairs(pairs); + + + performedSearch_ = true; + 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, bool force = false) + { + if(force) currentIter_ = 0; + performedSearch_ = false; + + if( !performSearch() ) return true; + + NBSLevels_.build(activeRange, incld); + + NBSLevels_.findPairs(pairs); + + performedSearch_ = true; + return true; + } + +}; + +} + +#endif diff --git a/src/Interaction/contactSearch/wallMappings/cellMapping.H b/src/Interaction/contactSearch/wallMappings/cellMapping.H new file mode 100644 index 00000000..baaaaf9d --- /dev/null +++ b/src/Interaction/contactSearch/wallMappings/cellMapping.H @@ -0,0 +1,150 @@ +/*------------------------------- 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 __cellMapping_H__ +#define __cellMapping_H__ + +#include "cellsWallLevel0.H" +#include "dictionary.H" + + +namespace pFlow +{ + +template< + typename executionSpace + > +class cellMapping +{ +public: + + using cellsWallLevel0Type = cellsWallLevel0; + + using IdType = typename cellsWallLevel0Type::IdType; + + using IndexType = typename cellsWallLevel0Type::IndexType; + + using Cells = typename cellsWallLevel0Type::Cells; + + using CellType = typename Cells::CellType; + + using execution_space = typename cellsWallLevel0Type::execution_space; + + using memory_space = typename cellsWallLevel0Type::memory_space; + + using iBoxType = iBox; + + +protected: + + // - update frequency + int32 updateFrequency_=1; + + real cellExtent_; + + int32 currentIter_ = 0; + + /// a broad search has been occured during last pass? + bool performedSearch_ = false; + + cellsWallLevel0Type cellsWallLevle_; + +private: + + bool performSearch() + { + if(currentIter_ % updateFrequency_ == 0) + { + currentIter_++; + return true; + + }else + { + currentIter_++; + return false; + } + } + +public: + + TypeNameNV("cellMapping"); + + cellMapping( + const dictionary& dict, + int32 numLevels, + const Vector& ppCells, + int32 numPoints, + int32 numElements, + const ViewType1D& points, + const ViewType1D& vertices + ) + : + updateFrequency_( + max( + dict.getValOrSet( + "updateFrequency", + 1), + 1)), + cellExtent_( + max( + dict.getValOrSet( + "cellExtent", + 0.5), + 0.5)), + cellsWallLevle_( + ppCells[0], + cellExtent_, + numPoints, + numElements, + points, + vertices + ) + {} + + + bool enterBoadSearch()const + { + return currentIter_%updateFrequency_==0; + } + + bool performedSearch()const + { + return performedSearch_; + } + + template + bool broadSearch(PairsContainer& pairs, particleMapType& particleMap, bool force=false) + { + if(force) currentIter_ = 0; + performedSearch_= false; + if(!performSearch())return true; + + cellsWallLevle_.broadSearch(pairs, particleMap); + + performedSearch_ = true; + return true; + } + +}; // cellMapping + +} // pFlow + + +#endif diff --git a/src/Interaction/contactSearch/wallMappings/cellsSimple.H b/src/Interaction/contactSearch/wallMappings/cellsWallLevel0.H similarity index 66% rename from src/Interaction/contactSearch/wallMappings/cellsSimple.H rename to src/Interaction/contactSearch/wallMappings/cellsWallLevel0.H index 6f14ae14..309fe149 100644 --- a/src/Interaction/contactSearch/wallMappings/cellsSimple.H +++ b/src/Interaction/contactSearch/wallMappings/cellsWallLevel0.H @@ -18,57 +18,48 @@ Licence: -----------------------------------------------------------------------------*/ -#ifndef __cellsSimple_H__ -#define __cellsSimple_H__ +#ifndef __cellsWallLevel0_H__ +#define __cellsWallLevel0_H__ #include "types.H" #include "KokkosTypes.H" #include "cells.H" #include "iBox.H" -#include "dictionary.H" + namespace pFlow { template< - typename executionSpace, - typename idType, - typename indexType = int32 + typename executionSpace > -class cellsSimple +class cellsWallLevel0 : - public cells + public cells { public: - using IdType = idType; + using IdType = int32; - using IndexType = indexType; + using IndexType = int32; using Cells = cells; using CellType = typename Cells::CellType; - using ExecutionSpace = executionSpace; + using execution_space = executionSpace; - using memory_space = typename ExecutionSpace::memory_space; + using memory_space = typename execution_space::memory_space; - using iBoxType = iBox; + using iBoxType = iBox; - bool constexpr static LOOP_ELEMENT_RANGE = true; - class TagFindCellRange2{}; + protected: // - box extent - real cellExtent_ = 0.6; - - // - update frequency - int32 updateFrequency_=1; - - - int32 currentIter_ = 0; + real cellExtent_ = 0.5; // - number of triangle elements int32 numElements_ = 0; @@ -76,13 +67,10 @@ protected: // - number of points int32 numPoints_ = 0; - /// a broad search has been occured during last pass? - bool performedSearch_ = false; - - // - ref to vectices + // - ref to vectices (borrowed) ViewType1D vertices_; - // - ref to points in the trisurface + // - ref to points in the trisurface (borrowed) ViewType1D points_; // cell range of element/triangle bounding box @@ -90,53 +78,40 @@ protected: using tpPWContactSearch = Kokkos::TeamPolicy< - ExecutionSpace, + execution_space, Kokkos::Schedule, Kokkos::IndexType >; using rpFindCellRange2Type = - Kokkos::RangePolicy>; + Kokkos::RangePolicy>; FUNCTION_H void allocateArrays() { - Kokkos::realloc( elementBox_, numElements_); - } - -private: - - bool performSearch() - { - if(currentIter_ % updateFrequency_ == 0) - { - currentIter_++; - return true; - - }else - { - currentIter_++; - return false; - } + reallocNoInit( elementBox_, numElements_); } public: - TypeNameNV("cellsSimple"); + TypeNameNV("cellsWallLevel0"); + + INLINE_FUNCTION_HD + cellsWallLevel0(){} FUNCTION_H - cellsSimple( - Cells ppCells, - real cellExtent, - int32 numPoints, - int32 numElements, + cellsWallLevel0( + const Cells& ppCells, + real cellExtent, + int32 numPoints, + int32 numElements, const ViewType1D& points, const ViewType1D& vertices ) : Cells(ppCells), - cellExtent_( max(cellExtent, 0.6 ) ), + cellExtent_( max(cellExtent, 0.5 ) ), numElements_(numElements), numPoints_(numPoints), vertices_(vertices), @@ -146,45 +121,11 @@ public: allocateArrays(); } - cellsSimple( - dictionary dict, - Cells ppCells, - int32 numPoints, - int32 numElements, - const ViewType1D& points, - const ViewType1D& vertices - ) - : - Cells(ppCells), - numElements_(numElements), - numPoints_(numPoints), - vertices_(vertices), - points_(points) - { - - updateFrequency_ = dict.getVal( - "updateFrequency" ); - updateFrequency_ = max(updateFrequency_,1); - - cellExtent_ = dict.getVal( - "cellExtent"); - - cellExtent_ = max(cellExtent_,0.6); - - allocateArrays(); - - } - - constexpr bool loopElementRange()const - { - return LOOP_ELEMENT_RANGE; - } - + // - host call // reset triangle elements if they have changed - FUNCTION_H bool resetElements( - int32 numElements, + int32 numElements, int32 numPoints, ViewType1D& points, ViewType1D& vertices ) @@ -212,29 +153,17 @@ public: return numElements_; } - bool enterBoadSearch()const - { - return currentIter_%updateFrequency_==0; - } - - bool performedSearch()const - { - return performedSearch_; - } + template - bool broadSearch(PairsContainer& pairs, particleMapType& particleMap, bool force=false) + bool broadSearch(PairsContainer& pairs, particleMapType& particleMap) { - if(force) currentIter_ = 0; - performedSearch_= false; - if(!performSearch())return true; // map walls onto the cells this->build(); this->particleWallFindPairs(pairs, particleMap); - performedSearch_ = true; return true; } @@ -257,7 +186,7 @@ public: while (getFull) { - getFull = findPairsElementRange(pairs, particleMap); + getFull = findPairsElementRangeCount(pairs, particleMap.getCellIterator(0)); if(getFull) { @@ -269,7 +198,7 @@ public: Info<<"Contact pair container capacity increased from "<< oldCap << " to " - << pairs.capacity() <<" in cellsSimple."< - int32 findPairsElementRange(PairsContainer& pairs, particleMapType& particleMap) + template + int32 findPairsElementRangeCount(PairsContainer& pairs, CellIteratorType cellIter) { int32 getFull =0; const auto pwPairs = pairs; const auto elementBox = elementBox_; - - auto cellIter = particleMap.getCellIterator(); Kokkos::parallel_reduce( "cellsSimple::findPairsElementRangeModified2", tpPWContactSearch(numElements_, Kokkos::AUTO), LAMBDA_HD( - const typename tpPWContactSearch::member_type & teamMember, int32& valueToUpdate){ + const typename tpPWContactSearch::member_type & teamMember, + int32& valueToUpdate){ const int32 iTri = teamMember.league_rank(); @@ -349,9 +277,9 @@ public: } -}; +}; // cellsWallLevel0 -} +} // pFlow -#endif +#endif // __cellsWallLevel0_H__ diff --git a/src/Interaction/contactSearch/wallMappings/cellsWallLevels.H b/src/Interaction/contactSearch/wallMappings/cellsWallLevels.H new file mode 100644 index 00000000..276ed058 --- /dev/null +++ b/src/Interaction/contactSearch/wallMappings/cellsWallLevels.H @@ -0,0 +1,152 @@ +/*------------------------------- 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 __cellsWallLevels_H__ +#define __cellsWallLevels_H__ + +#include "cellsWallLevel0.H" + +namespace pFlow +{ + +template< + typename executionSpace + > +class cellsWallLevels +{ +public: + + using cellsWallLevel0Type = cellsWallLevel0; + + using IdType = typename cellsWallLevel0Type::IdType; + + using IndexType = typename cellsWallLevel0Type::IndexType; + + using Cells = typename cellsWallLevel0Type::Cells; + + using CellType = typename Cells::CellType; + + using execution_space = typename cellsWallLevel0Type::execution_space; + + using memory_space = typename cellsWallLevel0Type::memory_space; + + using iBoxType = iBox; + +protected: + + int32 numLevles_=1; + + + Vector cellsWallLevels_; + +public: + + TypeNameNV("cellsWallLevels"); + + FUNCTION_H + cellsWallLevels( + int32 numLevels, + const Vector& cellsLevels, + real cellExtent, + int32 numPoints, + int32 numElements, + const ViewType1D& points, + const ViewType1D& vertices + ) + : + numLevles_(numLevels), + cellsWallLevels_("cellsWallLevels",numLevels, numLevels, RESERVE()) + { + + + + for(int32 lvl=0; lvl + bool broadSearch(PairsContainer& pairs, particleMapType& particleMap) + { + + // map walls onto the cells + for(int32 lvl=0; lvlparticleWallFindPairs(pairs, particleMap); + + return true; + } + + template + bool particleWallFindPairs(PairsContainer& pairs, particleMapType& particleMap) + { + + int32 getFull = 1; + + while (getFull) + { + getFull = 0; + for(int32 lvl=0; lvl +class multiGridMapping +{ +public: + + using CellsWallLevelType = cellsWallLevels; + + using IdType = typename CellsWallLevelType::IdType; + + using IndexType = typename CellsWallLevelType::IndexType; + + using Cells = typename CellsWallLevelType::Cells; + + using CellType = typename Cells::CellType; + + using execution_space = typename CellsWallLevelType::execution_space; + + using memory_space = typename CellsWallLevelType::memory_space; + + using iBoxType = iBox; + + +protected: + + // - update frequency + int32 updateFrequency_=1; + + real cellExtent_; + + int32 currentIter_ = 0; + + /// a broad search has been occured during last pass? + bool performedSearch_ = false; + + CellsWallLevelType cellsWallLevle_; + +private: + + bool performSearch() + { + if(currentIter_ % updateFrequency_ == 0) + { + currentIter_++; + return true; + + }else + { + currentIter_++; + return false; + } + } + +public: + + TypeNameNV("multiGridMapping"); + + multiGridMapping( + const dictionary& dict, + int32 numLevels, + const Vector& ppCells, + int32 numPoints, + int32 numElements, + const ViewType1D& points, + const ViewType1D& vertices + ) + : + updateFrequency_( + max( + dict.getVal("updateFrequency"), + 1)), + cellExtent_( + max( + dict.getVal("cellExtent"), + 0.5)), + cellsWallLevle_( + numLevels, + ppCells, + cellExtent_, + numPoints, + numElements, + points, + vertices + ) + { + + Report(3)<<"Multi-grid wall mapping with "<< + yellowText(numLevels)<<" levels has been created."< + bool broadSearch(PairsContainer& pairs, particleMapType& particleMap, bool force=false) + { + if(force) currentIter_ = 0; + performedSearch_= false; + if(!performSearch())return true; + + + cellsWallLevle_.broadSearch(pairs, particleMap); + + + performedSearch_ = true; + return true; + } + +}; // multiGridMapping + +} // pFlow + + +#endif diff --git a/src/phasicFlow/Kokkos/KokkosTypes.H b/src/phasicFlow/Kokkos/KokkosTypes.H index 63ff7d40..a2742609 100644 --- a/src/phasicFlow/Kokkos/KokkosTypes.H +++ b/src/phasicFlow/Kokkos/KokkosTypes.H @@ -61,6 +61,9 @@ template template using ViewType1D = Kokkos::View; +template + using DualViewType1D = Kokkos::View; + template using ViewType3D = Kokkos::View; diff --git a/src/phasicFlow/structuredData/box/box.C b/src/phasicFlow/structuredData/box/box.C index 731c8106..652d76d1 100644 --- a/src/phasicFlow/structuredData/box/box.C +++ b/src/phasicFlow/structuredData/box/box.C @@ -21,16 +21,6 @@ Licence: #include "box.H" -FUNCTION_HD -pFlow::box::box -( - const realx3& minP, - const realx3& maxP -) -: - min_(minP), - max_(maxP) -{} FUNCTION_H pFlow::box::box diff --git a/src/phasicFlow/structuredData/box/box.H b/src/phasicFlow/structuredData/box/box.H index 85b25df3..b8d790f6 100644 --- a/src/phasicFlow/structuredData/box/box.H +++ b/src/phasicFlow/structuredData/box/box.H @@ -34,10 +34,10 @@ class box protected: // - min point - realx3 min_; + realx3 min_{0,0,0}; // - max point - realx3 max_; + realx3 max_{1,1,1}; public: @@ -45,8 +45,16 @@ public: TypeNameNV("box"); //// - Constructors - FUNCTION_HD - box(const realx3& minP, const realx3& maxP); + INLINE_FUNCTION_HD + box(){} + + INLINE_FUNCTION_HD + box(const realx3& minP, const realx3& maxP) + : + min_(minP), + max_(maxP) + {} + FUNCTION_H box(const dictionary& dict);