diff --git a/src/Interaction/CMakeLists.txt b/src/Interaction/CMakeLists.txt index 25d43014..7e9212da 100644 --- a/src/Interaction/CMakeLists.txt +++ b/src/Interaction/CMakeLists.txt @@ -1,5 +1,10 @@ set(SourceFiles +contactSearch/methods/cellBased/NBS/mapperNBS.cpp +contactSearch/methods/cellBased/NBS/mapperNBSKernels.cpp +contactSearch/methods/cellBased/NBS/NBSLevel0.cpp +contactSearch/methods/cellBased/NBS/NBS.cpp +contactSearch/methods/cellBased/NBS/cellsWallLevel0.cpp contactSearch/contactSearch/contactSearch.cpp contactSearch/ContactSearch/ContactSearchs.cpp interaction/interaction.cpp diff --git a/src/Interaction/Models/contactForce/linearCF.hpp b/src/Interaction/Models/contactForce/linearCF.hpp index 5ef6da3e..88bf2108 100644 --- a/src/Interaction/Models/contactForce/linearCF.hpp +++ b/src/Interaction/Models/contactForce/linearCF.hpp @@ -133,8 +133,9 @@ protected: return false; } - realVector etha_n(nElem); - realVector etha_t(nElem); + realVector etha_n("etha_n", nElem); + realVector etha_t("etha_t", nElem); + ForAll(i , kn) { etha_n[i] = -2.0*log(en[i])*sqrt(kn[i])/ @@ -144,7 +145,7 @@ protected: sqrt(pow(log(et[i]),2.0)+ pow(Pi,2.0)); } - Vector prop(nElem); + Vector prop("prop", nElem); ForAll(i,kn) { prop[i] = {kn[i], kt[i], etha_n[i], etha_t[i], mu[i]}; @@ -219,10 +220,10 @@ public: void contactForce ( const real dt, - const int32 i, - const int32 j, - const int32 propId_i, - const int32 propId_j, + const uint32 i, + const uint32 j, + const uint32 propId_i, + const uint32 propId_j, const real Ri, const real Rj, const real ovrlp_n, diff --git a/src/Interaction/Models/contactForce/nonLinearCF.hpp b/src/Interaction/Models/contactForce/nonLinearCF.hpp index b9f42d69..3ba1c5c5 100644 --- a/src/Interaction/Models/contactForce/nonLinearCF.hpp +++ b/src/Interaction/Models/contactForce/nonLinearCF.hpp @@ -121,7 +121,7 @@ protected: } - realVector etha_n(nElem); + realVector etha_n("etha_n",nElem); ForAll(i , en) { @@ -137,7 +137,7 @@ protected: } - Vector prop(nElem); + Vector prop("prop",nElem); ForAll(i,Yeff) { prop[i] = {Yeff[i], Geff[i], etha_n[i], mu[i]}; @@ -214,10 +214,10 @@ public: void contactForce ( const real dt, - const int32 i, - const int32 j, - const int32 propId_i, - const int32 propId_j, + const uint32 i, + const uint32 j, + const uint32 propId_i, + const uint32 propId_j, const real Ri, const real Rj, const real ovrlp_n, diff --git a/src/Interaction/Models/contactForce/nonLinearMod.hpp b/src/Interaction/Models/contactForce/nonLinearMod.hpp index 3bfbd29b..9a771970 100644 --- a/src/Interaction/Models/contactForce/nonLinearMod.hpp +++ b/src/Interaction/Models/contactForce/nonLinearMod.hpp @@ -121,7 +121,7 @@ protected: } - Vector prop(nElem); + Vector prop("prop",nElem); ForAll(i,Yeff) { prop[i] = {Yeff[i], Geff[i], etha_n[i], mu[i]}; @@ -198,10 +198,10 @@ public: void contactForce ( const real dt, - const int32 i, - const int32 j, - const int32 propId_i, - const int32 propId_j, + const uint32 i, + const uint32 j, + const uint32 propId_i, + const uint32 propId_j, const real Ri, const real Rj, const real ovrlp_n, diff --git a/src/Interaction/Models/rolling/normalRolling.hpp b/src/Interaction/Models/rolling/normalRolling.hpp index de850e3c..80b8e321 100644 --- a/src/Interaction/Models/rolling/normalRolling.hpp +++ b/src/Interaction/Models/rolling/normalRolling.hpp @@ -76,10 +76,10 @@ public: void rollingFriction ( const real dt, - const int32 i, - const int32 j, - const int32 propId_i, - const int32 propId_j, + const uint32 i, + const uint32 j, + const uint32 propId_i, + const uint32 propId_j, const real Ri, const real Rj, const realx3& wi, diff --git a/src/Interaction/contactLists/sortedContactList.hpp b/src/Interaction/contactLists/sortedContactList.hpp index 45c8e4b0..73f0ce4b 100644 --- a/src/Interaction/contactLists/sortedContactList.hpp +++ b/src/Interaction/contactLists/sortedContactList.hpp @@ -53,7 +53,7 @@ protected: ViewType1D values_; - int32 size0_ = 0; + uint32 size0_ = 0; ViewType1D sortedPairs0_; @@ -73,7 +73,7 @@ protected: using rpReFillPairs = Kokkos::RangePolicy< ExecutionSpace, Kokkos::Schedule, - Kokkos::IndexType, + Kokkos::IndexType, TagReFillPairs>; public: @@ -81,7 +81,7 @@ public: TypeInfoNV("sortedContactList"); - sortedContactList(int32 initialSize =1) + explicit sortedContactList(uint32 initialSize =1) : SortedPairs(initialSize), values_("values", SortedPairs::capacity()), @@ -114,31 +114,31 @@ public: } INLINE_FUNCTION_HD - ValueType getValue(int32 idx)const + ValueType getValue(uint32 idx)const { return values_[idx]; } INLINE_FUNCTION_HD - void setValue(int32 idx, const ValueType& val)const + void setValue(uint32 idx, const ValueType& val)const { values_[idx] = val; } INLINE_FUNCTION_HD - void operator()(TagReFillPairs, int32 idx)const + void operator()(TagReFillPairs, uint32 idx)const { - auto searchLen = max(size0_/1000,10); - auto start = max(0,idx-searchLen); - auto end = min(size0_,idx+searchLen); + uint32 searchLen = max(size0_/1000u,10u); + uint32 start = idx-min(searchLen,idx); + uint32 end = min(size0_,idx+searchLen); auto newPair = this->sortedPairs_[idx]; if( auto idx0 = binarySearch( sortedPairs0_, start, end, newPair); - idx0>=0) + idx0!=-1) { values_[idx] = values0_[idx0]; } @@ -147,7 +147,7 @@ public: start, end, newPair); - idx0>=0) + idx0!=-1) { values_[idx] = values0_[idx0]; diff --git a/src/Interaction/contactLists/sortedPairs.hpp b/src/Interaction/contactLists/sortedPairs.hpp index d8b74a50..12593810 100644 --- a/src/Interaction/contactLists/sortedPairs.hpp +++ b/src/Interaction/contactLists/sortedPairs.hpp @@ -52,24 +52,24 @@ public: { using PairType = typename sortedPairs::PairType; - int32 size_; + uint32 size_; ViewType1D sortedParis_; INLINE_FUNCTION_HD - int32 size()const { return size_; } + uint32 size()const { return size_; } INLINE_FUNCTION_HD - int32 loopCount()const { return size_; } + uint32 loopCount()const { return size_; } INLINE_FUNCTION_HD - bool isValid(int32 i)const { return i flags_; + ViewType1D flags_; ViewType1D sortedPairs_; using rpFillFlag = Kokkos::RangePolicy< ExecutionSpace, Kokkos::Schedule, - Kokkos::IndexType, + Kokkos::IndexType, TagFillFlag >; using rpFillPairs = Kokkos::RangePolicy< ExecutionSpace, Kokkos::Schedule, - Kokkos::IndexType, + Kokkos::IndexType, TagFillPairs>; public: @@ -110,7 +110,7 @@ public: // constructors - sortedPairs(int32 initialSize =1) + explicit sortedPairs(uint32 initialSize =1) : UnsortedPairs(initialSize), flags_("flags_",UnsortedPairs::capacity()+1), @@ -134,7 +134,7 @@ public: // return the pair at index idx // perform no check for size and existance INLINE_FUNCTION_HD - PairType getPair(int32 idx)const + PairType getPair(uint32 idx)const { return sortedPairs_[idx]; } @@ -142,7 +142,7 @@ public: // - Device/host call // return the pair at index idx INLINE_FUNCTION_HD - bool getPair(int32 idx, PairType& p)const + bool getPair(uint32 idx, PairType& p)const { if(isValid(idx)) { @@ -156,7 +156,7 @@ public: } INLINE_FUNCTION_HD - bool isValid(int32 idx)const + bool isValid(uint32 idx)const { return idx < size_; } @@ -164,12 +164,12 @@ public: //use this when the value of size_ is updated INLINE_FUNCTION_H - int32 size()const + uint32 size()const { return size_; } - int32 loopCount()const + uint32 loopCount()const { return size_; } @@ -189,7 +189,7 @@ public: void prepareSorted() { // first check the size of flags_ - int32 capacity = UnsortedPairs::capacity(); + uint32 capacity = UnsortedPairs::capacity(); if( capacity+1 > flags_.size() ) { @@ -218,7 +218,7 @@ public: if( size_>sortedPairs_.size() ) { // get more space to prevent reallocations in next iterations - int32 len = size_*1.1+1; + uint32 len = size_*1.1+1; reallocNoInit(sortedPairs_, len); } @@ -235,7 +235,7 @@ public: } INLINE_FUNCTION_HD - void operator()(TagFillFlag, int32 i)const + void operator()(TagFillFlag, uint32 i)const { if(this->container_.valid_at(i) ) flags_[i] = 1; @@ -244,7 +244,7 @@ public: } INLINE_FUNCTION_HD - void operator()(TagFillPairs, int32 i)const + void operator()(TagFillPairs, uint32 i)const { auto fi = flags_[i]; if(fi!=flags_[i+1]) diff --git a/src/Interaction/contactLists/unsortedContactList.hpp b/src/Interaction/contactLists/unsortedContactList.hpp index f2200fd9..2a3150b1 100644 --- a/src/Interaction/contactLists/unsortedContactList.hpp +++ b/src/Interaction/contactLists/unsortedContactList.hpp @@ -72,7 +72,7 @@ protected: using rpFillPairs = Kokkos::RangePolicy< ExecutionSpace, Kokkos::Schedule, - Kokkos::IndexType, + Kokkos::IndexType, TagReFillPairs>; @@ -80,7 +80,7 @@ public: TypeInfoNV("unsortedContactList"); - unsortedContactList(int32 capacity=1) + explicit unsortedContactList(uint32 capacity=1) : UnsortedPairs(capacity), values_("values", UnsortedPairs::capacity()), @@ -122,7 +122,7 @@ public: INLINE_FUNCTION_HD bool getValue(const PairType& p, ValueType& val)const { - if(auto idx = this->find(p); idx>=0) + if(auto idx = this->find(p); idx!=-1) { val = getValue(idx); return true; @@ -131,7 +131,7 @@ public: } INLINE_FUNCTION_HD - void setValue(int32 idx, const ValueType& val)const + void setValue(uint32 idx, const ValueType& val)const { values_[idx] = val; } @@ -139,7 +139,7 @@ public: INLINE_FUNCTION_HD bool setValue(const PairType& p, const ValueType& val)const { - if(auto idx = this->find(p); idx>=0) + if(uint32 idx = this->find(p); idx!=-1) { setValue(idx, val); return true;; @@ -148,13 +148,13 @@ public: } INLINE_FUNCTION_HD - void operator()(TagReFillPairs, int32 idx)const + void operator()(TagReFillPairs, uint32 idx)const { if( this->isValid(idx) ) { - if( int32 idx0 = + if( uint32 idx0 = container0_.find(this->getPair(idx)); - idx0>=0 ) + idx0!=-1 ) { values_[idx] = values0_[idx0]; } diff --git a/src/Interaction/contactLists/unsortedPairs.hpp b/src/Interaction/contactLists/unsortedPairs.hpp index e9434751..fd4f117d 100644 --- a/src/Interaction/contactLists/unsortedPairs.hpp +++ b/src/Interaction/contactLists/unsortedPairs.hpp @@ -41,7 +41,7 @@ public: using memory_space = typename ExecutionSpace::memory_space; - using PairType = kPair; + using PairType = Pair; using ContainerType = unorderedSet; @@ -52,19 +52,19 @@ public: ContainerType Container_; INLINE_FUNCTION_HD - int32 size()const { return Container_.size(); } + uint32 size()const { return Container_.size(); } INLINE_FUNCTION_HD - int32 loopCount()const { return Container_.capacity(); } + uint32 loopCount()const { return Container_.capacity(); } INLINE_FUNCTION_HD - bool isValid(int32 idx)const { return Container_.valid_at(idx); } + bool isValid(uint32 idx)const { return Container_.valid_at(idx); } INLINE_FUNCTION_HD - PairType getPair(int idx)const { return Container_.key_at(idx); } + PairType getPair(uint32 idx)const { return Container_.key_at(idx); } INLINE_FUNCTION_HD - bool getPair(int32 idx, PairType& pair)const { + bool getPair(uint32 idx, PairType& pair)const { if(Container_.valid_at(idx)) { pair = Container_.key_at(idx); return true; @@ -84,7 +84,7 @@ public: TypeInfoNV("unsorderedPairs"); // constructor - unsortedPairs(int32 capacity=1) + explicit unsortedPairs(uint32 capacity=1) : container_(capacity) // the minimum capacity would be 128 {} @@ -102,7 +102,7 @@ public: // - Device call INLINE_FUNCTION_HD - int32 insert(idType i, idType j)const + uint32 insert(idType i, idType j)const { if(auto insertResult = container_.insert(PairType(i,j)); insertResult.failed()) return -1; @@ -112,7 +112,7 @@ public: } INLINE_FUNCTION_HD - int32 insert(const PairType& p)const + uint32 insert(const PairType& p)const { if(auto insertResult = container_.insert(p); insertResult.failed()) return -1; @@ -125,7 +125,7 @@ public: // return the pair at index idx // perform no check for size and existance INLINE_FUNCTION_HD - PairType getPair(int32 idx)const + PairType getPair(uint32 idx)const { return container_.key_at(idx); } @@ -133,7 +133,7 @@ public: // - Device call // return the pair at index idx INLINE_FUNCTION_HD - bool getPair(int32 idx, PairType& p)const + bool getPair(uint32 idx, PairType& p)const { if(container_.valid_at(idx)) { @@ -148,7 +148,7 @@ public: } INLINE_FUNCTION_HD - int32 find(const PairType & p)const + uint32 find(const PairType & p)const { if( auto idx = container_.find(p); idx != Kokkos::UnorderedMapInvalidIndex ) @@ -158,26 +158,26 @@ public: } INLINE_FUNCTION_HD - bool isValid(int32 idx)const + bool isValid(uint32 idx)const { return container_.valid_at(idx); } INLINE_FUNCTION_HD - int32 capacity() const + uint32 capacity() const { return container_.capacity(); } - int32 loopCount()const + uint32 loopCount()const { return container_.capacity(); } //use this when the value of size_ is updated INLINE_FUNCTION_H - int32 size()const + uint32 size()const { return container_.size(); } @@ -190,7 +190,7 @@ public: /// increase the capacity of the container by at-least len /// the content will be erased. INLINE_FUNCTION_H - void increaseCapacityBy(int32 len) + void increaseCapacityBy(uint32 len) { uint newCap = container_.capacity()+len; this->clear(); diff --git a/src/Interaction/contactSearch/ContactSearch/ContactSearch.hpp b/src/Interaction/contactSearch/ContactSearch/ContactSearch.hpp index d6428064..84fe8c95 100644 --- a/src/Interaction/contactSearch/ContactSearch/ContactSearch.hpp +++ b/src/Interaction/contactSearch/ContactSearch/ContactSearch.hpp @@ -22,16 +22,17 @@ Licence: #ifndef __ContactSearch_hpp__ #define __ContactSearch_hpp__ - +#include "contactSearchGlobals.hpp" #include "contactSearch.hpp" #include "box.hpp" +#include "particles.hpp" +#include "geometry.hpp" namespace pFlow { template< - template class BaseMethod, - template class WallMapping + class searchMethod > class ContactSearch : @@ -39,32 +40,19 @@ class ContactSearch { public: - using IdType = typename contactSearch::IdType; + using IdType = uint32; - using IndexType = typename contactSearch::IndexType; + using ExecutionSpace = DefaultExecutionSpace; - using ExecutionSpace = typename contactSearch::ExecutionSpace; + using SearchMethodType = searchMethod; - using PairContainerType = typename contactSearch::PairContainerType; - - using ParticleContactSearchType = - BaseMethod< - ExecutionSpace>; - - using WallMappingType = - WallMapping< - ExecutionSpace>; - -protected: +private: - - uniquePtr particleContactSearch_ = nullptr; - - uniquePtr wallMapping_ = nullptr; + uniquePtr ppwContactSearch_ = nullptr; public: - TypeInfoTemplate2("ContactSearch", ParticleContactSearchType, WallMappingType); + TypeInfoTemplate11("ContactSearch", SearchMethodType); ContactSearch( const dictionary& csDict, @@ -74,54 +62,41 @@ public: Timers& timers) : contactSearch(csDict, domain, prtcl, geom, timers) - { auto method = dict().getVal("method"); - auto wmMethod = dict().getVal("wallMapping"); - + auto nbDict = dict().subDict(method+"Info"); real minD, maxD; this->Particles().boundingSphereMinMax(minD, maxD); - const auto& position = this->Particles().pointPosition().deviceVectorAll(); - const auto& diam = this->Particles().boundingSphere().deviceVectorAll(); + const auto& position = this->Particles().pointPosition().deviceViewAll(); + const auto& flags = this->Particles().dynPointStruct().activePointsMaskDevice(); + const auto& diam = this->Particles().boundingSphere().deviceViewAll(); - particleContactSearch_ = - makeUnique + uint32 wnPoints = this->Geometry().numPoints(); + uint32 wnTri = this->Geometry().size(); + const auto& wPoints = this->Geometry().points().deviceViewAll(); + const auto& wVertices = this->Geometry().vertices().deviceViewAll(); + + ppwContactSearch_ = + makeUnique ( nbDict, - this->domain(), + this->domainBox(), minD, maxD, position, - diam - ); - REPORT(2)<<"Contact search algorithm for particle-particle is "<< - greenText(particleContactSearch_().typeName())<Geometry().numPoints(); - int32 wnTri = this->Geometry().size(); - - const auto& wPoints = this->Geometry().points().deviceVectorAll(); - const auto& wVertices = this->Geometry().vertices().deviceVectorAll(); - - wallMapping_ = - makeUnique( - wmDict, - particleContactSearch_().numLevels(), - particleContactSearch_().getCellsLevels(), + flags, + diam, wnPoints, - wnTri, + wnTri, wPoints, wVertices - ); - REPORT(2)<<"Wall mapping algorithm for particle-wall is "<< - greenText(wallMapping_().typeName())<< endREPORT; + ); + REPORT(2)<<"Contact search algorithm for particle-particle is "<< + Green_Text(ppwContactSearch_().typeName())<Particles().activeRange(); - - sphereSphereTimer_.start(); - - if(this->Particles().allActive()) - { - particleContactSearch_().broadSearch(ppPairs, activeRange, force); - } - else - { - particleContactSearch_().broadSearch(ppPairs, activeRange, this->Particles().activePointsMaskD(), force); - } - - sphereSphereTimer_.end(); - - } - else + fatalErrorInFunction; return false; - - if(wallMapping_) - { - sphereWallTimer_.start(); - wallMapping_().broadSearch(pwPairs, particleContactSearch_(), force); - sphereWallTimer_.end(); } - else - return false; + ppTimer().end(); - - + + return true; } - bool ppEnterBroadSearch()const override + bool enterBroadSearch(uint32 iter, real t, real dt)const override { - if(particleContactSearch_) + if(ppwContactSearch_) { - return particleContactSearch_().enterBoadSearch(); + return ppwContactSearch_().performSearch(iter); } return false; } - bool pwEnterBroadSearch()const override + bool performedBroadSearch()const override { - if(wallMapping_) - { - return wallMapping_().enterBoadSearch(); - } - return false; + return ppwContactSearch_().performedSearch(); } - - bool ppPerformedBroadSearch()const override - { - if(particleContactSearch_) - { - return particleContactSearch_().performedSearch(); - } - return false; - } - - - bool pwPerformedBroadSearch()const override - { - if(wallMapping_) - { - return wallMapping_().performedSearch(); - } - return false; - } - - /*bool update(const eventMessage& msg) - { - if(msg.isSizeChanged() ) - { - auto newSize = this->prtcl().size(); - if(!particleContactSearch_().objectSizeChanged(newSize)) - { - fatalErrorInFunction<< - "erro in changing the size for particleContactSearch_ \n"; - return false; - } - } - - if(msg.isCapacityChanged() ) - { - auto newSize = this->prtcl().capacity(); - if(!particleContactSearch_().objectSizeChanged(newSize)) - { - fatalErrorInFunction<< - "erro in changing the capacity for particleContactSearch_ \n"; - return false; - } - } - - return true; - }*/ - - }; } diff --git a/src/Interaction/contactSearch/ContactSearch/ContactSearchs.cpp b/src/Interaction/contactSearch/ContactSearch/ContactSearchs.cpp index 043e9b78..61024220 100644 --- a/src/Interaction/contactSearch/ContactSearch/ContactSearchs.cpp +++ b/src/Interaction/contactSearch/ContactSearch/ContactSearchs.cpp @@ -20,11 +20,11 @@ Licence: #include "ContactSearch.hpp" -#include "cellMapping.hpp" +//#include "cellMapping.hpp" #include "NBS.hpp" -#include "multiGridNBS.hpp" -#include "multiGridMapping.hpp" +//#include "multiGridNBS.hpp" +//#include "multiGridMapping.hpp" -template class pFlow::ContactSearch; -template class pFlow::ContactSearch; +template class pFlow::ContactSearch; +//template class pFlow::ContactSearch; diff --git a/src/Interaction/contactSearch/contactSearch/contactSearch.cpp b/src/Interaction/contactSearch/contactSearch/contactSearch.cpp index 17c1dbd3..5753d4a2 100644 --- a/src/Interaction/contactSearch/contactSearch/contactSearch.cpp +++ b/src/Interaction/contactSearch/contactSearch/contactSearch.cpp @@ -19,7 +19,7 @@ Licence: -----------------------------------------------------------------------------*/ #include "contactSearch.hpp" - +#include "streams.hpp" pFlow::contactSearch::contactSearch( @@ -29,11 +29,12 @@ pFlow::contactSearch::contactSearch( const geometry& geom, Timers& timers) : - interactionBase(prtcl, geom), - domain_(domain), - dict_(dict), - sphereSphereTimer_("particle-particle contact search", &timers), - sphereWallTimer_("particle-wall contact search", &timers) + domainBox_(domain), + particles_(prtcl), + geometry_(geom), + ppTimer_("particle-particle contact search", &timers), + pwTimer_("particle-wall contact search", &timers), + dict_(dict) { } @@ -46,20 +47,16 @@ pFlow::uniquePtr pFlow::contactSearch::create( const geometry& geom, Timers& timers) { + word baseMethName = dict.getVal("method"); - word baseMethName = dict.getVal("method"); - word wallMethod = dict.getVal("wallMapping"); + auto model = angleBracketsNames("ContactSearch", baseMethName); - auto model = angleBracketsNames2("ContactSearch", baseMethName, wallMethod); - - - REPORT(1)<<"Selecting contact search model . . ."<; + const geometry& geometry_; -protected: + Timer ppTimer_; - const box& domain_; + Timer pwTimer_; dictionary dict_; - Timer sphereSphereTimer_; - - Timer sphereWallTimer_; - - auto& dict() - { - return dict_; - } - public: TypeInfo("contactSearch"); @@ -88,35 +81,51 @@ public: (dict, domain, prtcl, geom, timers) ); - const auto& domain()const - { - return domain_; - } - + const auto& dict()const { return dict_; } + const auto& domainBox()const + { + return domainBox_; + } + + const auto& Particles()const + { + return particles_; + } + + const auto& Geometry()const + { + return geometry_; + } + + auto& ppTimer() + { + return ppTimer_; + } + + auto& pwTimer() + { + return pwTimer_; + } virtual bool broadSearch( - PairContainerType& ppPairs, - PairContainerType& pwPairs, + uint32 iter, + real t, + real dt, + csPairContainerType& ppPairs, + csPairContainerType& pwPairs, bool force = false) = 0; virtual - bool ppEnterBroadSearch()const = 0; + bool enterBroadSearch(uint32 iter, real t, real dt)const = 0; virtual - bool pwEnterBroadSearch()const = 0; - - virtual - bool ppPerformedBroadSearch()const = 0; - - virtual - bool pwPerformedBroadSearch()const = 0; - + bool performedBroadSearch()const = 0; static uniquePtr create( diff --git a/src/Interaction/contactSearch/contactSearchFunctions.hpp b/src/Interaction/contactSearch/contactSearchFunctions.hpp index 55879c95..96ee9596 100644 --- a/src/Interaction/contactSearch/contactSearchFunctions.hpp +++ b/src/Interaction/contactSearch/contactSearchFunctions.hpp @@ -96,17 +96,6 @@ void indexToCell(const indexType idx, const iBox& box, triple(); - -//constexpr int32 maxCellIndex = largestPositive(); +using csPairContainerType = unsortedPairs; } - -#endif //__interactionTypes_hpp__ +#endif diff --git a/src/Interaction/contactSearch/methods/NBS.hpp b/src/Interaction/contactSearch/methods/NBS.hpp deleted file mode 100644 index fb3cc4f8..00000000 --- a/src/Interaction/contactSearch/methods/NBS.hpp +++ /dev/null @@ -1,209 +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_hpp__ -#define __NBS_hpp__ - -#include "NBSLevel0.hpp" - -namespace pFlow -{ - - -template -class NBS -{ -public: - - using NBSLevel0Type = NBSLevel0; - - using cellIterator = typename NBSLevel0Type::cellIterator; - - 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; - - -protected: - - real sizeRatio_ = 1.0; - - int32 updateFrequency_= 1; - - int32 currentIter_ = 0; - - bool performedSearch_ = false; - - NBSLevel0Type NBSLevel0_; - -private: - - bool performSearch() - { - if(currentIter_ % updateFrequency_ == 0) - { - currentIter_++; - return true; - - }else - { - currentIter_++; - return false; - } - } - -public: - - TypeInfoNV("NBS"); - - NBS( - const dictionary& dict, - const box& domain, - real minSize, - real maxSize, - const ViewType1D& position, - const ViewType1D& diam) - : - 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; - - 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_; - } - - 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(); - } - - bool objectSizeChanged(int32 newSize) - { - NBSLevel0_.checkAllocateNext(newSize); - return true; - } - - // - 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; - - - NBSLevel0_.build(activeRange); - - NBSLevel0_.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; - - NBSLevel0_.build(activeRange, incld); - - NBSLevel0_.findPairs(pairs); - - performedSearch_ = true; - return true; - } - -}; - -} - -#endif diff --git a/src/Interaction/contactSearch/methods/NBSCrossLoop.hpp b/src/Interaction/contactSearch/methods/NBSCrossLoop.hpp deleted file mode 100644 index adce7409..00000000 --- a/src/Interaction/contactSearch/methods/NBSCrossLoop.hpp +++ /dev/null @@ -1,82 +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. - ------------------------------------------------------------------------------*/ - - -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.hpp b/src/Interaction/contactSearch/methods/NBSLevel.hpp deleted file mode 100644 index f355d566..00000000 --- a/src/Interaction/contactSearch/methods/NBSLevel.hpp +++ /dev/null @@ -1,127 +0,0 @@ -#ifndef __NBSLevel_hpp__ -#define __NBSLevel_hpp__ - - -#include "NBSLevel0.hpp" - - -namespace pFlow -{ - -INLINE_FUNCTION_HD -int32x3 mapIndexLevels(const int32x3& ind, int32 lowerLevel, int32 upperLevel); - -template -class -NBSLevel -: - public NBSLevel0 -{ -public: - - using NBSLevel0Type = NBSLevel0; - - using cellIterator = typename NBSLevel0Type::cellIterator; - - 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 mdrPolicyFindPairs = typename NBSLevel0Type::mdrPolicyFindPairs; - - using HeadType = typename NBSLevel0Type::HeadType; - - using NextType = typename NBSLevel0Type::NextType; - - template - friend class NBSLevels; - -protected: - - int32 level_ = 0; - - -public: - - TypeInfoNV("NBSLevel0"); - - INLINE_FUNCTION_HD - NBSLevel(){} - - NBSLevel( - int32 lvl, - const box& domain, - real cellSize, - real sizeRatio, - const ViewType1D& position, - const ViewType1D& diam) - : - NBSLevel0Type( - domain, - cellSize, - sizeRatio, - position, - diam, - lvl==0), - level_(lvl) - {} - - INLINE_FUNCTION_HD - NBSLevel(const NBSLevel&) = default; - - INLINE_FUNCTION_HD - NBSLevel& operator = (const NBSLevel&) = default; - - INLINE_FUNCTION_HD - ~NBSLevel() = default; - - INLINE_FUNCTION_HD - auto level()const - { - return level_; - } - - template - INLINE_FUNCTION_H - int32 findPairsCountCross(PairsContainer& pairs, NBSLevel& upperLevel) - { - - 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.hpp" - }, notInsertedPairs); - - return notInsertedPairs; - } - - -}; - -INLINE_FUNCTION_HD -int32x3 mapIndexLevels( const 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.hpp b/src/Interaction/contactSearch/methods/NBSLevel0.hpp deleted file mode 100644 index d9ae2f83..00000000 --- a/src/Interaction/contactSearch/methods/NBSLevel0.hpp +++ /dev/null @@ -1,240 +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 __NBSLevel0_hpp__ -#define __NBSLevel0_hpp__ - -#include "mapperNBS.hpp" - -namespace pFlow -{ - - -template -class NBSLevel0 -: - public mapperNBS -{ -public: - - using MapperType = mapperNBS; - - using cellIterator = typename MapperType::cellIterator; - - 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; - - using HeadType = typename MapperType::HeadType; - - using NextType = typename MapperType::NextType; - - struct TagFindPairs{}; - - -protected: - - real sizeRatio_ = 1.0; - - // borrowed ownership - ViewType1D diameter_; - - - using mdrPolicyFindPairs = - Kokkos::MDRangePolicy< - Kokkos::Rank<3>, - Kokkos::Schedule, - execution_space>; - - static INLINE_FUNCTION_HD - void Swap(int32& x, int32& y) - { - int32 tmp = x; - x = y; - y = tmp; - } - -public: - - TypeInfoNV("NBSLevel0"); - - INLINE_FUNCTION_HD - 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, - bool nextOwner = true) - : - MapperType(domain, cellSize, position, nextOwner), - sizeRatio_(sizeRatio), - diameter_(diam) - {} - - INLINE_FUNCTION_HD - NBSLevel0(const NBSLevel0&) = default; - - INLINE_FUNCTION_HD - NBSLevel0& operator = (const NBSLevel0&) = default; - - INLINE_FUNCTION_HD - ~NBSLevel0()=default; - - - INLINE_FUNCTION_HD - auto sizeRatio()const - { - 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 - // 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) - { - - - 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,500) ; - - auto oldCap = pairs.capacity(); - - pairs.increaseCapacityBy(len); - - INFORMATION<< "The contact pair container capacity increased from "<< - oldCap << " to "< - 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.hpp" - }, notInsertedPairs); - - return notInsertedPairs; - - } - -}; - -} // pFlow - -#endif // __NBSLevel0_hpp__ diff --git a/src/Interaction/contactSearch/methods/NBSLevels.hpp b/src/Interaction/contactSearch/methods/NBSLevels.hpp deleted file mode 100644 index 30603bd3..00000000 --- a/src/Interaction/contactSearch/methods/NBSLevels.hpp +++ /dev/null @@ -1,438 +0,0 @@ -#ifndef __NBSLevels_hpp__ -#define __NBSLevels_hpp__ - -#include "NBSLevel.hpp" -#include "NBSLevel0.hpp" -#include "KokkosTypes.hpp" - -namespace pFlow -{ - -template -class NBSLevels -{ - -public: - - using NBSLevelType = NBSLevel; - - using cellIterator = typename NBSLevelType::cellIterator; - - using IdType = typename NBSLevelType::IdType; - - using IndexType = typename NBSLevelType::IndexType; - - using Cells = typename NBSLevelType::Cells; - - using CellType = typename Cells::CellType; - - using execution_space = typename NBSLevelType::execution_space; - - using memory_space = typename NBSLevelType::memory_space; - - using realRange = kPair; - -protected: - - real minSize_; - - real maxSize_; - - int32 numLevels_=1; - - Vector nbsLevels_; - - ViewType1D sizeRangeLevels_; - - ViewType1D sizeRangeLevelsHost_; - - ViewType1D maxSizeLevels_; - - ViewType1D maxSizeLevelsHost_; - - ViewType1D particleLevel_; - - range activeRange_{0,0}; - - using rangePolicyType = - Kokkos::RangePolicy< - Kokkos::IndexType, - Kokkos::Schedule, - execution_space>; - - int32 setNumLevels() - { - - int32 maxOvermin = static_cast(maxSize_/minSize_); - - if (maxOvermin <=1) - return 1; - else if(maxOvermin<=3) - return 2; - else if(maxOvermin<=7) - return 3; - else if(maxOvermin<15) - return 4; - else if(maxOvermin<31) - return 5; - else if(maxOvermin<63) - return 6; - else if(maxOvermin <127) - return 7; - else - { - fatalErrorInFunction<< - "size ratio is not valid for multi-grid NBS "<< maxOvermin<=0; lvl--) - { - - 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; - } - - copy(sizeRangeLevels_, sizeRangeLevelsHost_); - copy(maxSizeLevels_, maxSizeLevelsHost_); - - REPORT(2)<<"Grids with "<< yellowText(numLevels_)<< " levels have been created."<& 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) - : - minSize_(minSize), - maxSize_(maxSize), - numLevels_(setNumLevels()), - nbsLevels_("nbsLevels", numLevels_, numLevels_, RESERVE()), - sizeRangeLevels_("sizeRangeLevels", numLevels_), - sizeRangeLevelsHost_("sizeRangeLevelsHost", numLevels_), - maxSizeLevels_("maxSizeLevels", numLevels_), - maxSizeLevelsHost_("maxSizeLevelsHost", numLevels_) - { - - setDiameterRange(sizeRatio); - - 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); - - INFORMATION<< "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, - LAMBDA_HD(int32 i){ - - 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(); - } - - template - INLINE_FUNCTION_H - void build(range activeRange, IncludeFunction incld) - { - // 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, - 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) - { - - if(last > particleLevel_.size()) - { - reallocNoInit(particleLevel_,last); - } - - auto diameter = nbsLevels_[0].diameter(); - auto const maxSizes = maxSizeLevels_; - auto particleLevel = particleLevel_; - auto const sizeRatio = 0.999*nbsLevels_[0].sizeRatio(); - - int8 maxLvl = sizeRangeLevels_.size(); - - rangePolicyType rPolicy(first, last); - - Kokkos::parallel_for( - "NBSLevels::findParticleLevel", - rPolicy, - LAMBDA_HD(int32 i) - { - for(int8 lvl = 0; lvl(-1); - }); - Kokkos::fence(); - - return true; - } - - -}; //NBSLevels - -} - -#endif diff --git a/src/Interaction/contactSearch/methods/NBSLoop.hpp b/src/Interaction/contactSearch/methods/NBSLoop.hpp deleted file mode 100644 index 8aed3e26..00000000 --- a/src/Interaction/contactSearch/methods/NBSLoop.hpp +++ /dev/null @@ -1,101 +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. - ------------------------------------------------------------------------------*/ - - -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 = sizeRatio_* diameter_[m]; - - // the same cell - n = this->next_(m); - - while(n >-1) - { - - 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; - - if(lm>ln) Swap(lm,ln); - if( auto res = pairs.insert(lm,ln); res <0) - { - getFullUpdate++; - } - } - - n = this->next_(n); - } - - // neighbor cells - CellType neighborCell; - for(int32 ni=0; ni<13; ni++) - { - if(ni==0) neighborCell = currentCell + CellType( 0, 0,-1); - else if(ni==1) neighborCell = currentCell + CellType(-1, 0,-1); - else if(ni==2) neighborCell = currentCell + CellType(-1, 0, 0); - else if(ni==3) neighborCell = currentCell + CellType(-1, 0, 1); - else if(ni==4) neighborCell = currentCell + CellType( 0,-1,-1); - else if(ni==5) neighborCell = currentCell + CellType( 0,-1, 0); - else if(ni==6) neighborCell = currentCell + CellType( 0,-1, 1); - else if(ni==7) neighborCell = currentCell + CellType(-1,-1,-1); - else if(ni==8) neighborCell = currentCell + CellType(-1,-1, 0); - else if(ni==9) neighborCell = currentCell + CellType(-1,-1, 1); - else if(ni==10) neighborCell = currentCell + CellType( 1,-1,-1); - else if(ni==11) neighborCell = currentCell + CellType( 1,-1, 0); - else if(ni==12) neighborCell = currentCell + CellType( 1,-1, 1); - - if( this->isInRange(neighborCell) ) - { - - n = this->head_(neighborCell.x(), neighborCell.y(), neighborCell.z()); - while( n>-1) - { - - 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; - if(lm>ln) Swap(lm,ln); - if( auto res = pairs.insert(lm,ln); res <0) - { - getFullUpdate++; - } - } - n = this->next_[n]; - } - } - - } - m = this->next_[m]; -} - diff --git a/src/Interaction/interaction/interactionBase.hpp b/src/Interaction/contactSearch/methods/cellBased/NBS/NBS.cpp similarity index 50% rename from src/Interaction/interaction/interactionBase.hpp rename to src/Interaction/contactSearch/methods/cellBased/NBS/NBS.cpp index b5c9651e..9c8b3a29 100644 --- a/src/Interaction/interaction/interactionBase.hpp +++ b/src/Interaction/contactSearch/methods/cellBased/NBS/NBS.cpp @@ -18,67 +18,50 @@ Licence: -----------------------------------------------------------------------------*/ -#ifndef __interactionBase_hpp__ -#define __interactionBase_hpp__ +#include "NBS.hpp" -#include "interactionTypes.hpp" -#include "particles.hpp" -#include "geometry.hpp" - -namespace pFlow +pFlow::NBS::NBS +( + const dictionary& dict, + const box& domainBox, + real minBSSize, + real maxBSSize, + const deviceViewType1D &position, + const pFlagTypeDevice &flags, + const deviceViewType1D &diam, + uint32 nWallPoints, + uint32 nWallElements, + const ViewType1D& wallPoints, + const ViewType1D& wallVertices +) +: + particleWallContactSearchs( + dict, + domainBox, + minBSSize, + maxBSSize, + position, + flags, + diam), + sizeRatio_(max(dict.getVal("sizeRatio"), 1.0)), + cellExtent_(max(dict.getVal("cellExtent"), 0.5)), + adjustableBox_(dict.getVal("adjustableBox")), + NBSLevel0_ + ( + this->domainBox_, + maxBSSize, + sizeRatio_, + position, + flags, + adjustableBox_() + ), + cellsWallLevel0_ + ( + cellExtent_, + nWallPoints, + nWallElements, + wallPoints, + wallVertices + ) { - -class interactionBase -{ -public: - - using IndexType = CELL_INDEX_TYPE; - - using IdType = ID_TYPE; - - using ExecutionSpace = DefaultExecutionSpace; - -protected: - - const particles& particles_; - - const geometry& geometry_; - -public: - - interactionBase( - const particles& prtcl, - const geometry& geom) - : - particles_(prtcl), - geometry_(geom) - {} - - inline - const auto& pStruct()const - { - return particles_.pStruct(); - } - - inline - const auto& surface()const - { - return geometry_.surface(); - } - - inline - const auto& Particles()const - { - return particles_; - } - - inline auto& Geometry()const - { - return geometry_; - } -}; - } - - -#endif //__interactionBase_hpp__ diff --git a/src/Interaction/contactSearch/methods/cellBased/NBS/NBS.hpp b/src/Interaction/contactSearch/methods/cellBased/NBS/NBS.hpp new file mode 100644 index 00000000..f0f83436 --- /dev/null +++ b/src/Interaction/contactSearch/methods/cellBased/NBS/NBS.hpp @@ -0,0 +1,148 @@ +/*------------------------------- 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_hpp__ +#define __NBS_hpp__ + +#include "particleWallContactSearchs.hpp" +#include "NBSLevel0.hpp" +#include "cellsWallLevel0.hpp" +#include "Vector.hpp" + +namespace pFlow +{ + + +class NBS +: + public particleWallContactSearchs +{ +public: + + using CellIterator = typename NBSLevel0::CellIterator; + +private: + + real sizeRatio_ = 1.0; + + real cellExtent_ = 0.5; + + Logical adjustableBox_; + + NBSLevel0 NBSLevel0_; + + cellsWallLevel0 cellsWallLevel0_; + +protected: + + friend particleWallContactSearchs; + + bool impl_broadSearch + ( + csPairContainerType& ppPairs, + csPairContainerType& pwPairs, + const deviceViewType1D& pointPos, + const pFlagTypeDevice& flags, + const deviceViewType1D& diameter + ) + { + bool searchBoxChanged; + if( !NBSLevel0_.broadSearch( + ppPairs, + pointPos, + flags, + diameter, + searchBoxChanged)) + { + fatalErrorInFunction<< + "Error in broadSearch for NBS (particle-particle)"< &position, + const pFlagTypeDevice &flags, + const deviceViewType1D &diam, + uint32 nWallPoints, + uint32 nWallElements, + const ViewType1D& wallPoints, + const ViewType1D& wallVertices); + + + INLINE_FUNCTION_HD + NBS(const NBS&) = default; + + INLINE_FUNCTION_HD + NBS(NBS&&) = default; + + INLINE_FUNCTION_HD + NBS& operator = (const NBS&) = default; + + INLINE_FUNCTION_HD + NBS& operator = (NBS&&) = default; + + INLINE_FUNCTION_HD + ~NBS()=default; + + + //// - Methods + + uint32 numLevels()const + { + return 1; + } + + auto getCellIterator([[maybe_unused]] uint32 lvl)const + { + return NBSLevel0_.getCellIterator(); + } + + Vector getDomainCellsLevels()const + { + return Vector("Cells", 1, NBSLevel0_.getDomainCells()); + } + +}; + +} + +#endif diff --git a/src/Interaction/contactSearch/methods/cellBased/NBS/NBSLevel0.cpp b/src/Interaction/contactSearch/methods/cellBased/NBS/NBSLevel0.cpp new file mode 100644 index 00000000..96e47917 --- /dev/null +++ b/src/Interaction/contactSearch/methods/cellBased/NBS/NBSLevel0.cpp @@ -0,0 +1,118 @@ +/*------------------------------- 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. + +-----------------------------------------------------------------------------*/ + +#include "NBSLevel0.hpp" +#include "streams.hpp" + +bool pFlow::NBSLevel0::findPairs +( + csPairContainerType &pairs, + const deviceViewType1D &pointPos, + const pFlagTypeDevice &flags, + const deviceViewType1D &diameter +) +{ + uint32 getFull = 1; + + // loop until the container size fits the numebr of contact pairs + while (getFull > 0) + { + + getFull = pFlow::NBSLevel0Kernels::findPairsCount + ( + pairs, + sizeRatio_, + pointPos, + flags, + diameter, + getCellIterator() + ); + + if(getFull) + { + // - resize the container + // note that getFull now shows the number of failed insertions. + uint32 len = max(getFull,500u) ; + + auto oldCap = pairs.capacity(); + + pairs.increaseCapacityBy(len); + + INFORMATION<< "The contact pair container capacity increased from "<< + oldCap << " to "<& pointPos, + const pFlagTypeDevice& flags, + bool adjustableBox +) +: + mapperNBS + ( + domain, + cellSize, + pointPos, + flags, + adjustableBox, + true + ), + sizeRatio_(sizeRatio) +{ + +} + +bool pFlow::NBSLevel0::broadSearch +( + csPairContainerType &pairs, + const deviceViewType1D &pointPos, + const pFlagTypeDevice &flags, + const deviceViewType1D &diameter, + bool& searchBoxChanged +) +{ + if(!build(pointPos, flags, searchBoxChanged)) + { + fatalErrorInFunction; + return false; + } + + if(!findPairs(pairs, pointPos, flags, diameter)) + { + fatalErrorInFunction; + return false; + } + + return true; +} + + diff --git a/src/Interaction/contactSearch/methods/cellBased/NBS/NBSLevel0.hpp b/src/Interaction/contactSearch/methods/cellBased/NBS/NBSLevel0.hpp new file mode 100644 index 00000000..2f87c151 --- /dev/null +++ b/src/Interaction/contactSearch/methods/cellBased/NBS/NBSLevel0.hpp @@ -0,0 +1,98 @@ +/*------------------------------- 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_hpp__ +#define __NBSLevel0_hpp__ + +#include "NBSLevel0Kernels.hpp" +#include "mapperNBS.hpp" + +namespace pFlow +{ + + +class NBSLevel0 +: + public mapperNBS +{ +public: + + using MapperType = mapperNBS; + + using CellIterator = typename MapperType::CellIterator; + + using HeadType = typename MapperType::HeadType; + + using NextType = typename MapperType::NextType; + + +private: + + real sizeRatio_ = 1.0; + + bool findPairs + ( + csPairContainerType& pairs, + const deviceViewType1D& pointPos, + const pFlagTypeDevice& flags, + const deviceViewType1D& diameter + ); + +public: + + TypeInfoNV("NBSLevel0"); + + INLINE_FUNCTION_HD + NBSLevel0() = default; + + NBSLevel0( + const box& domain, + real cellSize, + real sizeRatio, + const deviceViewType1D& pointPos, + const pFlagTypeDevice& flags, + bool adjustableBox); + + INLINE_FUNCTION_HD + NBSLevel0(const NBSLevel0&) = default; + + INLINE_FUNCTION_HD + NBSLevel0& operator = (const NBSLevel0&) = default; + + INLINE_FUNCTION_HD + ~NBSLevel0()=default; + + + bool broadSearch( + csPairContainerType& pairs, + const deviceViewType1D& pointPos, + const pFlagTypeDevice& flags, + const deviceViewType1D& diameter, + bool& searchBoxChanged); + + + + +}; + +} // pFlow + +#endif // __NBSLevel0_hpp__ diff --git a/src/Interaction/contactSearch/methods/cellBased/NBS/NBSLevel0Kernels.hpp b/src/Interaction/contactSearch/methods/cellBased/NBS/NBSLevel0Kernels.hpp new file mode 100644 index 00000000..31df063c --- /dev/null +++ b/src/Interaction/contactSearch/methods/cellBased/NBS/NBSLevel0Kernels.hpp @@ -0,0 +1,83 @@ +/*------------------------------- 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. + +-----------------------------------------------------------------------------*/ + +#include "mapperNBS.hpp" +#include "contactSearchGlobals.hpp" + +#include "streams.hpp" + +namespace pFlow::NBSLevel0Kernels +{ + +using mdrPolicyFindPairs = + Kokkos::MDRangePolicy< + Kokkos::IndexType, + Kokkos::Rank<3>, + Kokkos::Schedule, + DefaultExecutionSpace>; + +template +INLINE_FUNCTION_HD +void Swap(T& x, T& y) +{ + T tmp = x; + x = y; + y = tmp; +} + +INLINE_FUNCTION_HD +bool sphereSphereCheck(const realx3& p1, const realx3 p2, real d1, real d2) +{ + return length(p2-p1) < 0.5*(d2+d1); +} + +inline +uint32 findPairsCount +( + csPairContainerType& pairs, + real sizeRatio, + const deviceViewType1D& pointPos, + const pFlagTypeDevice& flags, + const deviceViewType1D& diameter, + mapperNBS::CellIterator cellIter +) +{ + + auto nCells = cellIter.numCells(); + + mdrPolicyFindPairs + mdrPolicy( + {0,0,0}, + {nCells.x(), nCells.y(), nCells.z()} ); + + uint32 notInsertedPairs = 0u; + + + Kokkos::parallel_reduce ( + "pFlow::NBSLevel0Kernels::findPairsCount", + mdrPolicy, + LAMBDA_HD(uint32 i, uint32 j, uint32 k, uint32& getFullUpdate){ + #include "NBSLoop.hpp" + }, notInsertedPairs); + + return notInsertedPairs; +} + +} \ No newline at end of file diff --git a/src/Interaction/contactSearch/methods/cellBased/NBS/NBSLoop.hpp b/src/Interaction/contactSearch/methods/cellBased/NBS/NBSLoop.hpp new file mode 100644 index 00000000..ab3eeea7 --- /dev/null +++ b/src/Interaction/contactSearch/methods/cellBased/NBS/NBSLoop.hpp @@ -0,0 +1,100 @@ +/*------------------------------- 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. + +-----------------------------------------------------------------------------*/ + + +uint32 m = cellIter.start(i,j,k); +int32x3 currentCell(i,j,k); +uint32 n = mapperNBS::NoPos; + +while( m != mapperNBS::NoPos) +{ + + auto p_m = pointPos[m]; + auto d_m = sizeRatio* diameter[m]; + + // the same cell + n = cellIter.getNext(m); + + while(n != mapperNBS::NoPos) + { + auto p_n = pointPos[n]; + auto d_n = sizeRatio*diameter[n]; + + if( sphereSphereCheck(p_m, p_n, d_m, d_n) ) + { + auto ln = n; + auto lm = m; + + if(lm>ln) Swap(lm,ln); + if( pairs.insert(lm,ln) == -1) + { + getFullUpdate++; + } + } + + n = cellIter.getNext(n);; + } + + // neighbor cells + int32x3 neighborCell; + for(uint32 ni=0u; ni<13; ni++) + { + if(ni==0) neighborCell = currentCell + int32x3( 0, 0,-1); + else if(ni==1) neighborCell = currentCell + int32x3(-1, 0,-1); + else if(ni==2) neighborCell = currentCell + int32x3(-1, 0, 0); + else if(ni==3) neighborCell = currentCell + int32x3(-1, 0, 1); + else if(ni==4) neighborCell = currentCell + int32x3( 0,-1,-1); + else if(ni==5) neighborCell = currentCell + int32x3( 0,-1, 0); + else if(ni==6) neighborCell = currentCell + int32x3( 0,-1, 1); + else if(ni==7) neighborCell = currentCell + int32x3(-1,-1,-1); + else if(ni==8) neighborCell = currentCell + int32x3(-1,-1, 0); + else if(ni==9) neighborCell = currentCell + int32x3(-1,-1, 1); + else if(ni==10) neighborCell = currentCell + int32x3( 1,-1,-1); + else if(ni==11) neighborCell = currentCell + int32x3( 1,-1, 0); + else if(ni==12) neighborCell = currentCell + int32x3( 1,-1, 1); + + if( neighborCell.x()>=0 && neighborCell.y()>=0 && neighborCell.z()>=0 && + neighborCell.x()ln) Swap(lm,ln); + if( pairs.insert(lm,ln) == -1) + { + getFullUpdate++; + } + } + n = cellIter.next(n); + } + } + } + m = cellIter.next(m); +} + diff --git a/src/Interaction/contactSearch/methods/cellBased/NBS/cellsWallLevel0.cpp b/src/Interaction/contactSearch/methods/cellBased/NBS/cellsWallLevel0.cpp new file mode 100644 index 00000000..5d85a118 --- /dev/null +++ b/src/Interaction/contactSearch/methods/cellBased/NBS/cellsWallLevel0.cpp @@ -0,0 +1,195 @@ +/*------------------------------- 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. + +-----------------------------------------------------------------------------*/ + +#include "cellsWallLevel0.hpp" +#include "streams.hpp" + +pFlow::cellsWallLevel0::cellsWallLevel0 +( + real cellExtent, + uint32 numPoints, + uint32 numElements, + const ViewType1D &points, + const ViewType1D &vertices +) +: + cellExtent_( max(cellExtent, 0.5 ) ), + numElements_(numElements), + numPoints_(numPoints), + vertices_(vertices), + points_(points) +{ + + allocateArrays(); +} + +bool pFlow::cellsWallLevel0::resetElements +( + uint32 numElements, + uint32 numPoints, + ViewType1D &points, + ViewType1D &vertices +) +{ + + numElements_ = numElements; + numPoints_ = numPoints; + points_ = points; + vertices_ = vertices; + + allocateArrays(); + + return true; +} + +bool pFlow::cellsWallLevel0::broadSearch +( + csPairContainerType &pairs, + const cells& searchBox, + const mapperNBS::CellIterator &particleMap +) +{ + + // map walls onto the cells + + this->build(searchBox); + + this->particleWallFindPairs(pairs, particleMap); + + return true; +} + +bool pFlow::cellsWallLevel0::build(const cells & searchBox) +{ + + Kokkos::parallel_for( + "pFlow::cellsWallLevel0::build", + deviceRPolicyStatic(0,numElements_), + CLASS_LAMBDA_HD(uint32 i) + { + auto v = vertices_[i]; + auto p1 = points_[v.x()]; + auto p2 = points_[v.y()]; + auto p3 = points_[v.z()]; + + realx3 minP; + realx3 maxP; + + searchBox.extendBox(p1, p2, p3, cellExtent_, minP, maxP); + elementBox_[i] = iBoxType(searchBox.pointIndex(minP), searchBox.pointIndex(maxP)); + }); + Kokkos::fence(); + + return true; +} + +bool pFlow::cellsWallLevel0::particleWallFindPairs +( + csPairContainerType &pairs, + const mapperNBS::CellIterator &particleMap +) +{ + + uint32 getFull = 1; + + while (getFull) + { + + getFull = findPairsElementRangeCount(pairs, particleMap); + + if(getFull) + { + // - resize the container + // note that getFull now shows the number of failed insertions. + uint32 len = max(getFull, 50u); + auto oldCap = pairs.capacity(); + pairs.increaseCapacityBy(len); + + INFORMATION<<"Contact pair container capacity increased from "<< + oldCap << " to " + << pairs.capacity() <<" in cellsWallLevel0."<(n), + static_cast(iTri) ) == -1 ) + innerUpdate++; + n = particleMap.next(n); + + } + + }, + getFull2 + ); + + if ( teamMember.team_rank() == 0 ) valueToUpdate += getFull2; + }, + getFull + ); + return getFull; +} \ No newline at end of file diff --git a/src/Interaction/contactSearch/methods/cellBased/NBS/cellsWallLevel0.hpp b/src/Interaction/contactSearch/methods/cellBased/NBS/cellsWallLevel0.hpp new file mode 100644 index 00000000..2ffb4414 --- /dev/null +++ b/src/Interaction/contactSearch/methods/cellBased/NBS/cellsWallLevel0.hpp @@ -0,0 +1,139 @@ +/*------------------------------- 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 __cellsWallLevel0_hpp__ +#define __cellsWallLevel0_hpp__ + +#include "contactSearchGlobals.hpp" +#include "contactSearchFunctions.hpp" +#include "mapperNBS.hpp" +#include "iBox.hpp" + + + +namespace pFlow +{ + +class cellsWallLevel0 +{ +public: + + using execution_space = csExecutionSpace; + + using memory_space = typename execution_space::memory_space; + + using iBoxType = iBox; + + class TagFindCellRange2{}; + +private: + + // - box extent + real cellExtent_ = 0.5; + + // - number of triangle elements + uint32 numElements_ = 0; + + // - number of points + uint32 numPoints_ = 0; + + // - ref to vectices (borrowed) + ViewType1D vertices_; + + // - ref to points in the trisurface (borrowed) + ViewType1D points_; + + // cell range of element/triangle bounding box + ViewType1D elementBox_; + + + using tpPWContactSearch = Kokkos::TeamPolicy< + execution_space, + Kokkos::Schedule, + Kokkos::IndexType>; + + FUNCTION_H + void allocateArrays() + { + reallocNoInit( elementBox_, numElements_); + } + +public: + + TypeInfoNV("cellsWallLevel0"); + + INLINE_FUNCTION_HD + cellsWallLevel0(){} + + FUNCTION_H + cellsWallLevel0( + real cellExtent, + uint32 numPoints, + uint32 numElements, + const ViewType1D& points, + const ViewType1D& vertices); + + + + // - host call + // reset triangle elements if they have changed + bool resetElements( + uint32 numElements, + uint32 numPoints, + ViewType1D& points, + ViewType1D& vertices); + + + INLINE_FUNCTION_HD + iBoxType elementBox(uint32 i)const + { + return elementBox_[i]; + } + + INLINE_FUNCTION_HD + uint32 numElements()const + { + return numElements_; + } + + bool broadSearch( + csPairContainerType& pairs, + const cells& searchBox, + const mapperNBS::CellIterator& particleMap); + + + bool build(const cells& searchBox); + + bool particleWallFindPairs( + csPairContainerType& pairs, + const mapperNBS::CellIterator& particleMap); + + int32 findPairsElementRangeCount( + csPairContainerType& pairs, + const mapperNBS::CellIterator& particleMap); + + + +}; // cellsWallLevel0 + +} // pFlow + + +#endif // __cellsWallLevel0_hpp__ diff --git a/src/Interaction/contactSearch/methods/cellBased/NBS/mapperNBS.cpp b/src/Interaction/contactSearch/methods/cellBased/NBS/mapperNBS.cpp new file mode 100644 index 00000000..416ec57d --- /dev/null +++ b/src/Interaction/contactSearch/methods/cellBased/NBS/mapperNBS.cpp @@ -0,0 +1,194 @@ +/*------------------------------- 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. + +-----------------------------------------------------------------------------*/ + +#include "mapperNBS.hpp" +#include "mapperNBSKernels.hpp" +#include "streams.hpp" + +pFlow::uint32 pFlow::mapperNBS::checkInterval_ = 1000; +pFlow::real pFlow::mapperNBS::enlargementFactor_ = 1; + +bool pFlow::mapperNBS::setSearchBox +( + const deviceViewType1D &pointPos, + const pFlagTypeDevice &flags, + real cellSize +) +{ + box domainBox = domainCells_.domainBox(); + + + if(adjustableBox_) + { + lastCheckForBox_ = buildCount_; + + realx3 minP; + realx3 maxP; + pFlow::mapperNBSKernels::findPointExtends + ( + pointPos, + flags, + minP, maxP + ); + + minP = max( minP - enlargementFactor_*cellSize, domainBox.minPoint()); + maxP = min( maxP + enlargementFactor_*cellSize, domainBox.maxPoint()); + + box searchBox = {minP, maxP}; + searchCells_ = cells(searchBox, cellSize); + INFORMATION<<"Search box for contact search has changed: "<< + "search box is ["< &pointPos, + const pFlagTypeDevice &flags, + bool adjustableBox, + bool nextOwner +) +: + domainCells_(domain, cellSize), + searchCells_(domain, cellSize), + adjustableBox_(adjustableBox), + nextOwner_(nextOwner) +{ + setSearchBox(pointPos, flags, cellSize); + + allocateArrays(flags.activeRange()); +} + +bool pFlow::mapperNBS::build +( + const deviceViewType1D& pointPos, + const pFlagTypeDevice & flags, + bool& searchBoxChanged +) +{ + auto aRange = flags.activeRange(); + buildCount_++; + if(adjustableBox_ && buildCount_%checkInterval_ == 0) + { + + if(searchBoxChanged = + setSearchBox(pointPos, flags, searchCells_.cellSize());searchBoxChanged) + { + allocateArrays(aRange); + } + + } + else + { + checkAllocateNext(aRange); + nullifyHead(); + nullifyNext(aRange); + } + + if( adjustableBox_ ) + { + if(!pFlow::mapperNBSKernels::buildListsReduce( + searchCells_, + head_, + next_, + pointPos, + flags) ) + { + + buildCount_++; + setSearchBox(pointPos, flags, searchCells_.cellSize()); + + searchBoxChanged = true; + + allocateArrays(flags.activeRange()); + + if(!pFlow::mapperNBSKernels::buildListsReduce( + searchCells_, + head_, + next_, + pointPos, + flags)) + { + fatalErrorInFunction<<"failed to build list in anjustable search box mode!"<; + + using NextType = deviceViewType1D; + + + static constexpr uint32 NoPos = 0xFFFFFFFF; + + class CellIterator + { + private: + HeadType head_; + + NextType next_; + + public: + + CellIterator(const HeadType& head, const NextType& next) + : + head_(head), + next_(next) + {} + + static constexpr uint32 NoPos = 0xFFFFFFFF; + + INLINE_FUNCTION_HD + int32x3 numCells()const { + return int32x3(head_.extent(0), head_.extent(1), head_.extent(2));} + + INLINE_FUNCTION_HD + uint32 start(int32 i, int32 j, int32 k)const { + return head_(i,j,k); } + + INLINE_FUNCTION_HD + uint32 getNext(uint32 n)const { + if(n == NoPos ) return NoPos; + return next_(n); } + + INLINE_FUNCTION_HD + uint32 next(uint32 n)const{ + return next_(n);} + }; + +private: + + cells domainCells_; + + cells searchCells_; + + HeadType head_{"NBS::head",1,1,1}; + + NextType next_{"NBS::next", 1}; + + uint32 nextCapacity_ = 0; + + uint32 lastCheckForBox_ = 0; + + uint32 buildCount_ = 0; + + bool adjustableBox_ = false; + + bool nextOwner_ = true; + + static uint32 checkInterval_; + + static real enlargementFactor_; + + bool setSearchBox( + const deviceViewType1D& pointPos, + const pFlagTypeDevice& flags, + real cellSize + ); + + void allocateArrays(rangeU32 nextRng); + + void checkAllocateNext(rangeU32 nextRng); + + void nullifyHead(); + + void nullifyNext(rangeU32 nextRng); + +public: + + TypeInfoNV("mapperNBS"); + + INLINE_FUNCTION_HD + mapperNBS() = default; + + mapperNBS( + const box& domain, + real cellSize, + const deviceViewType1D& pointPos, + const pFlagTypeDevice& flags, + bool adjustableBox, + bool nextOwner = true); + + INLINE_FUNCTION_HD + mapperNBS(const mapperNBS&) = default; + + INLINE_FUNCTION_HD + mapperNBS(mapperNBS&&) = default; + + INLINE_FUNCTION_HD + mapperNBS& operator = (const mapperNBS&) = default; + + INLINE_FUNCTION_HD + mapperNBS& operator = (mapperNBS&&) = default; + + INLINE_FUNCTION_HD + ~mapperNBS()=default; + + //// - Methods + + auto getCellIterator()const + { + return CellIterator(head_, next_); + } + + const auto& getDomainCells()const + { + return domainCells_; + } + + const auto& getSearchCells()const + { + return searchCells_; + } + + bool build( + const deviceViewType1D& pointPos, + const pFlagTypeDevice& flags, + bool& searchBoxChanged); + + + + +}; + +} // pFlow + +#endif // __mapperNBS_hpp__ + + + + +/* +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_; + } +*/ \ No newline at end of file diff --git a/src/Interaction/contactSearch/methods/cellBased/NBS/mapperNBSKernels.cpp b/src/Interaction/contactSearch/methods/cellBased/NBS/mapperNBSKernels.cpp new file mode 100644 index 00000000..01141eac --- /dev/null +++ b/src/Interaction/contactSearch/methods/cellBased/NBS/mapperNBSKernels.cpp @@ -0,0 +1,192 @@ + + +/*------------------------------- 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. + +-----------------------------------------------------------------------------*/ + +#include "mapperNBSKernels.hpp" + + +void pFlow::mapperNBSKernels::findPointExtends +( + const deviceViewType1D& points, + const pFlagTypeDevice& flags, + realx3& minPoint, + realx3& maxPoint +) +{ + if(flags.numActive() == 0) + { + minPoint = {0,0,0}; + maxPoint = {0,0,0}; + return; + } + + real minX; + real minY; + real minZ; + real maxX; + real maxY; + real maxZ; + auto aRange = flags.activeRange(); + Kokkos::parallel_reduce( + "pFlow::mapperNBSKernels::findPointExtends", + deviceRPolicyStatic(aRange.start(), aRange.end()), + LAMBDA_HD( + uint32 i, + real& minXUpdate, + real& minYUpdate, + real& minZUpdate, + real& maxXUpdate, + real& maxYUpdate, + real& maxZUpdate) + { + if(flags(i)) + { + auto p = points(i); + minXUpdate = min(p.x(), minXUpdate); + minYUpdate = min(p.y(), minYUpdate); + minZUpdate = min(p.z(), minZUpdate); + maxXUpdate = max(p.x(), maxXUpdate); + maxYUpdate = max(p.y(), maxYUpdate); + maxZUpdate = max(p.z(), maxZUpdate); + } + + }, + Kokkos::Min(minX), + Kokkos::Min(minY), + Kokkos::Min(minZ), + Kokkos::Max(maxX), + Kokkos::Max(maxY), + Kokkos::Max(maxZ) + ); + + minPoint = {minX, minY, minZ}; + maxPoint = {maxX, maxY, maxZ}; +} + +bool pFlow::mapperNBSKernels::buildListsReduce +( + const cells &searchCell, + const deviceViewType3D &head, + const deviceViewType1D &next, + const deviceViewType1D &points, + const pFlagTypeDevice &flags +) +{ + uint32 numOut = 0; + auto aRange = flags.activeRange(); + + if(flags.isAllActive()) + { + Kokkos::parallel_reduce + ( + "pFlow::mapperNBSKernels::buildListsReduce", + deviceRPolicyStatic(aRange.start(), aRange.end()), + LAMBDA_HD(uint32 i, uint32& valToUpdate) + { + int32x3 ind; + if( searchCell.pointIndexInDomain(points[i], ind) ) + { + uint32 old = Kokkos::atomic_exchange(&head(ind.x(), ind.y(), ind.z()), i); + next[i] = old; + } + else + { + valToUpdate++; + } + }, + numOut + ); + } + else + { + Kokkos::parallel_reduce + ( + "pFlow::mapperNBSKernels::buildListsReduce", + deviceRPolicyStatic(aRange.start(), aRange.end()), + LAMBDA_HD(uint32 i, uint32& valToUpdate) + { + int32x3 ind; + if( flags(i) ) + { + if( searchCell.pointIndexInDomain(points[i], ind) ) + { + uint32 old = Kokkos::atomic_exchange(&head(ind.x(), ind.y(), ind.z()), i); + next[i] = old; + } + else + { + valToUpdate++; + } + } + }, + numOut + ); + } + + return numOut == 0u ; +} + +bool pFlow::mapperNBSKernels::buildLists +( + const cells &searchCell, + const deviceViewType3D &head, + const deviceViewType1D &next, + const deviceViewType1D &points, + const pFlagTypeDevice &flags +) +{ + auto aRange = flags.activeRange(); + if(flags.isAllActive() ) + { + Kokkos::parallel_for + ( + "pFlow::mapperNBSKernels::buildLists", + deviceRPolicyStatic(aRange.start(), aRange.end()), + LAMBDA_HD(uint32 i) + { + auto ind = searchCell.pointIndex(points[i]); + uint32 old = Kokkos::atomic_exchange(&head(ind.x(), ind.y(), ind.z()), i); + next[i] = old; + } + ); + Kokkos::fence(); + } + else + { + Kokkos::parallel_for + ( + "pFlow::mapperNBSKernels::buildLists", + deviceRPolicyStatic(aRange.start(), aRange.end()), + LAMBDA_HD(uint32 i) + { + if( flags(i) ) + { + auto ind = searchCell.pointIndex(points[i]); + uint32 old = Kokkos::atomic_exchange(&head(ind.x(), ind.y(), ind.z()), i); + next[i] = old; + } + } + ); + Kokkos::fence(); + } + + return true; +} diff --git a/src/Interaction/interaction/demInteraction.hpp b/src/Interaction/contactSearch/methods/cellBased/NBS/mapperNBSKernels.hpp similarity index 59% rename from src/Interaction/interaction/demInteraction.hpp rename to src/Interaction/contactSearch/methods/cellBased/NBS/mapperNBSKernels.hpp index 8c74639c..6730ff22 100644 --- a/src/Interaction/interaction/demInteraction.hpp +++ b/src/Interaction/contactSearch/methods/cellBased/NBS/mapperNBSKernels.hpp @@ -18,42 +18,32 @@ Licence: -----------------------------------------------------------------------------*/ -#ifndef __demInteraction_hpp__ -#define __demInteraction_hpp__ +#include "phasicFlowKokkos.hpp" +#include "cells.hpp" +#include "pointFlag.hpp" -#include "property.hpp" -#include "demComponent.hpp" -#include "pointFields.hpp" -#include "triSurfaceFields.hpp" -namespace pFlow +namespace pFlow::mapperNBSKernels { -class demInteraction -: - public property, - public demComponent -{ -protected: +void findPointExtends( + const deviceViewType1D& points, + const pFlagTypeDevice& flags, + realx3& minPoint, + realx3& maxPoint); -public: - demInteraction(systemControl& control) - : - property(), - demComponent("interaction", control) - {} +bool buildListsReduce( + const cells& searchCell, + const deviceViewType3D& head, + const deviceViewType1D& next, + const deviceViewType1D& points, + const pFlagTypeDevice& flags); - demInteraction(systemControl& control, const fileSystem& file) - : - property(file), - demComponent("interaction", control) - {} +bool buildLists( + const cells& searchCell, + const deviceViewType3D& head, + const deviceViewType1D& next, + const deviceViewType1D& points, + const pFlagTypeDevice& flags); - - -}; - - -} - -#endif //__interaction_hpp__ +} \ No newline at end of file diff --git a/src/Interaction/contactSearch/cells.hpp b/src/Interaction/contactSearch/methods/cellBased/cells.hpp similarity index 62% rename from src/Interaction/contactSearch/cells.hpp rename to src/Interaction/contactSearch/methods/cellBased/cells.hpp index c77b1a69..c4b50813 100644 --- a/src/Interaction/contactSearch/cells.hpp +++ b/src/Interaction/contactSearch/methods/cellBased/cells.hpp @@ -28,58 +28,40 @@ Licence: namespace pFlow { -template class cells { -public: - - using CellType = triple; - -protected: +private: // - domain - box domain_{realx3(0.0), realx3(1.0)}; + box domainBox_{realx3(0.0), realx3(1.0)}; // - cell size - realx3 cellSize_{1,1,1}; - - CellType numCells_{1,1,1}; + real celldx_{1}; + int32x3 numCells_{1,1,1}; // - protected methods INLINE_FUNCTION_H void calculate() { - numCells_ = (domain_.maxPoint()-domain_.minPoint())/cellSize_ + realx3(1.0); - numCells_ = max( numCells_ , CellType(static_cast(1)) ); + numCells_ = (domainBox_.maxPoint()-domainBox_.minPoint())/celldx_ + realx3(1.0); + numCells_ = max( numCells_ , int32x3(1) ); } public: INLINE_FUNCTION_HD - cells() - {} + cells() = default; INLINE_FUNCTION_H cells(const box& domain, real cellSize) : - domain_(domain), - cellSize_(cellSize) + domainBox_(domain), + celldx_(cellSize) { calculate(); } - - INLINE_FUNCTION_H - cells(const box& domain, int32 nx, int32 ny, int32 nz) - : - domain_(domain), - cellSize_( - (domain_.maxPoint() - domain_.minPoint())/realx3(nx, ny, nz) - ), - numCells_(nx, ny, nz) - {} - INLINE_FUNCTION_HD cells(const cells&) = default; @@ -100,43 +82,36 @@ public: INLINE_FUNCTION_H void setCellSize(real cellSize) { - cellSize_ = cellSize; - calculate(); - } - - INLINE_FUNCTION_H - void setCellSize(realx3 cellSize) - { - cellSize_ = cellSize; + celldx_ = cellSize; calculate(); } INLINE_FUNCTION_HD - realx3 cellSize()const + real cellSize()const { - return cellSize_; + return celldx_; } INLINE_FUNCTION_HD - const CellType& numCells()const + const int32x3& numCells()const { return numCells_; } INLINE_FUNCTION_HD - indexType nx()const + int32 nx()const { return numCells_.x(); } INLINE_FUNCTION_HD - indexType ny()const + int32 ny()const { return numCells_.y(); } INLINE_FUNCTION_HD - indexType nz()const + int32 nz()const { return numCells_.z(); } @@ -149,22 +124,21 @@ public: static_cast(numCells_.z()); } - const auto& domain()const + const auto& domainBox()const { - return domain_; + return domainBox_; } INLINE_FUNCTION_HD - CellType pointIndex(const realx3& p)const + int32x3 pointIndex(const realx3& p)const { - return CellType( (p - domain_.minPoint())/cellSize_ ); + return int32x3( (p - domainBox_.minPoint())/celldx_ ); } INLINE_FUNCTION_HD - bool pointIndexInDomain(const realx3 p, CellType& index)const + bool pointIndexInDomain(const realx3 p, int32x3& index)const { - if( !domain_.isInside(p) ) return false; - + if(!inDomain(p))return false; index = this->pointIndex(p); return true; } @@ -172,11 +146,11 @@ public: INLINE_FUNCTION_HD bool inDomain(const realx3& p)const { - return domain_.isInside(p); + return domainBox_.isInside(p); } INLINE_FUNCTION_HD - bool isInRange(const CellType& cell)const + bool inCellRange(const int32x3& cell)const { if(cell.x()<0)return false; if(cell.y()<0)return false; @@ -188,7 +162,7 @@ public: } INLINE_FUNCTION_HD - bool isInRange(indexType i, indexType j, indexType k)const + bool inCellRange(int32 i, int32 j, int32 k)const { if(i<0)return false; if(j<0)return false; @@ -199,22 +173,7 @@ public: return true; } - INLINE_FUNCTION_HD - void extendBox( - const CellType& p1, - const CellType& p2, - const CellType& p3, - indexType extent, - CellType& minP, - CellType& maxP)const - { - minP = min( min( p1, p2), p3)-extent; - maxP = max( max( p1, p2), p3)+extent; - - minP = bound(minP); - maxP = bound(maxP); - } - + INLINE_FUNCTION_HD void extendBox( const realx3& p1, @@ -224,17 +183,17 @@ public: realx3& minP, realx3& maxP)const { - minP = min(min(p1,p2),p3) - extent*cellSize_ ; - maxP = max(max(p1,p2),p3) + extent*cellSize_ ; + minP = min(min(p1,p2),p3) - extent*celldx_ ; + maxP = max(max(p1,p2),p3) + extent*celldx_ ; minP = bound(minP); maxP = bound(maxP); } INLINE_FUNCTION_HD - CellType bound(CellType p)const + int32x3 bound(int32x3 p)const { - return CellType( + return int32x3( min( numCells_.x()-1, max(0,p.x())), min( numCells_.y()-1, max(0,p.y())), min( numCells_.z()-1, max(0,p.z())) @@ -245,9 +204,9 @@ public: realx3 bound(realx3 p)const { return realx3( - min( domain_.maxPoint().x(), max(domain_.minPoint().x(),p.x())), - min( domain_.maxPoint().y(), max(domain_.minPoint().y(),p.y())), - min( domain_.maxPoint().z(), max(domain_.minPoint().z(),p.z())) + min( domainBox_.maxPoint().x(), max(domainBox_.minPoint().x(),p.x())), + min( domainBox_.maxPoint().y(), max(domainBox_.minPoint().y(),p.y())), + min( domainBox_.maxPoint().z(), max(domainBox_.minPoint().z(),p.z())) ); } }; diff --git a/src/Interaction/contactSearch/methods/mapperNBS.hpp b/src/Interaction/contactSearch/methods/mapperNBS.hpp deleted file mode 100644 index 41b47068..00000000 --- a/src/Interaction/contactSearch/methods/mapperNBS.hpp +++ /dev/null @@ -1,389 +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 __mapperNBS_hpp__ -#define __mapperNBS_hpp__ - -#include "cells.hpp" -#include "contactSearchFunctions.hpp" -#include "baseAlgorithms.hpp" -#include "ViewAlgorithms.hpp" - -namespace pFlow -{ - -template -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; - - using HeadType = ViewType3D; - - using NextType = ViewType1D; - - class cellIterator - { - private: - HeadType head_; - - NextType 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_; - - bool nextOwner_ = true; - - // borrowed ownership - ViewType1D pointPosition_; - - using rangePolicyType = - Kokkos::RangePolicy< - Kokkos::IndexType, - Kokkos::Schedule, - execution_space>; - - INLINE_FUNCTION_H - void nullifyHead() - { - fill( - head_, - range(0,this->nx()), - range(0,this->ny()), - range(0,this->nz()), - static_cast(-1) - ); - } - - void nullifyNext(range nextRng) - { - if(!nextOwner_)return; - fill( - next_, - nextRng, - static_cast(-1) - ); - } - - void nullify() - { - nullifyHead(); - - 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_); - } - } - - void allocateHead() - { - reallocNoInit(head_, this->nx(), this->ny(), this->nz()); - } - - - -public: - - TypeInfoNV("mapperNBS"); - - INLINE_FUNCTION_HD - mapperNBS(){} - - mapperNBS( - const box& domain, - real cellSize, - const ViewType1D& position, - bool nextOwner = true) - : - Cells(domain, cellSize), - pointPosition_(position), - head_( - "mapperNBS::head_", - this->nx(), - this->ny(), - this->nz() - ), - next_("mapperNBS::next_",1), //,position.size()), - nextOwner_(nextOwner) - { - checkAllocateNext(pointPosition_.size()); - } - - mapperNBS( - const box& domain, - int32 nx, - int32 ny, - int32 nz, - const ViewType1D& position, - bool nextOwner = true) - : - Cells(domain, nx, ny, nz), - pointPosition_(position), - head_("mapperNBS::head_",nx,ny,nz), - next_("mapperNBS::next_",1), - nextOwner_(nextOwner) - { - 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 particlesCapcityChanged(int32 newCap) - { - 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) - { - 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_hpp__ \ No newline at end of file diff --git a/src/Interaction/contactSearch/methods/multiGridNBS.hpp b/src/Interaction/contactSearch/methods/multiGridNBS.hpp deleted file mode 100644 index dbce390b..00000000 --- a/src/Interaction/contactSearch/methods/multiGridNBS.hpp +++ /dev/null @@ -1,213 +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 __multiGridNBS_hpp__ -#define __multiGridNBS_hpp__ - -#include "NBSLevels.hpp" - -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: - - TypeInfoNV("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/methods/particleWallContactSearchs.cpp b/src/Interaction/contactSearch/methods/particleWallContactSearchs.cpp new file mode 100644 index 00000000..25064346 --- /dev/null +++ b/src/Interaction/contactSearch/methods/particleWallContactSearchs.cpp @@ -0,0 +1,78 @@ +/*------------------------------- 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. + +-----------------------------------------------------------------------------*/ + + +template +pFlow::particleWallContactSearchs::particleWallContactSearchs +( + const dictionary &dict, + const box &domain, + real minSize, + real maxSize, + const ViewType1D &position, + const pFlagTypeDevice &flags, + const ViewType1D &diam +) +: + domainBox_(domain), + updateInterval_ + ( + dict.getValOrSet("updateInterval", 1) + ) +{ + +} + +template +bool pFlow::particleWallContactSearchs::broadSearch + ( + uint32 iter, + real t, + real dt, + csPairContainerType& ppPairs, + csPairContainerType& pwPairs, + const deviceViewType1D& pointPos, + const pFlagTypeDevice& flags, + const deviceViewType1D& diameter, + bool force + ) + { + + performedSearch_ = false; + if( !performSearch(iter, force) ) return true; + + if(!getMethod().impl_broadSearch( + ppPairs, + pwPairs, + pointPos, + flags, + diameter)) + { + fatalErrorInFunction<< + "Error in performing particle-particle broadSearch in method"<< + getMethod().typeName()< +class particleWallContactSearchs +{ +public: + + using MethodType = method; + + using IndexType = uint32; + + using execution_space = DefaultExecutionSpace; + + using memory_space = typename execution_space::memory_space; + +private: + + // friend + friend MethodType; + + /// @brief box enclosing the simulation domain (local to processor) + box domainBox_; + + /*/// @brief box enclosing the area for contact search (region with points) + box searchBox_;*/ + + /// @brief update interval in terms of iteration numebr + uint32 updateInterval_= 1; + + /// @brief last iteration number which contact search has been performed + uint32 lastUpdated_ = 0; + + /// @brief performed search? + bool performedSearch_ = false; + +protected: + + inline + auto& getMethod() + { + return static_cast(*this); + } + + inline + const auto& getMethod()const + { + return static_cast(*this); + } + +public: + + particleWallContactSearchs( + const dictionary& dict, + const box& domain, + real minSize, + real maxSize, + const ViewType1D& position, + const pFlagTypeDevice &flags, + const ViewType1D& diam + ); + + bool broadSearch + ( + uint32 iter, + real t, + real dt, + csPairContainerType& ppPairs, + csPairContainerType& pwPairs, + const deviceViewType1D& pointPos, + const pFlagTypeDevice& flags, + const deviceViewType1D& diameter, + bool force = false + ); + + bool performedSearch()const + { + return performedSearch_; + } + + bool performSearch(uint32 iter, bool force = false)const + { + if((iter-lastUpdated_) % updateInterval_ == 0 || iter == 0 || force ) + { + return true; + } + return false; + } + +}; + + +} // pFlow + +#include "particleWallContactSearchs.cpp" + +#endif //__particleWallContactSearchs_hpp__ \ No newline at end of file diff --git a/src/Interaction/contactSearch/wallMappings/cellMapping.hpp b/src/Interaction/contactSearch/wallMappings/cellMapping.hpp deleted file mode 100644 index 5a10c248..00000000 --- a/src/Interaction/contactSearch/wallMappings/cellMapping.hpp +++ /dev/null @@ -1,150 +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 __cellMapping_hpp__ -#define __cellMapping_hpp__ - -#include "cellsWallLevel0.hpp" -#include "dictionary.hpp" - - -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: - - TypeInfoNV("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/cellsWallLevel0.hpp b/src/Interaction/contactSearch/wallMappings/cellsWallLevel0.hpp deleted file mode 100644 index 65034507..00000000 --- a/src/Interaction/contactSearch/wallMappings/cellsWallLevel0.hpp +++ /dev/null @@ -1,285 +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 __cellsWallLevel0_hpp__ -#define __cellsWallLevel0_hpp__ - -#include "types.hpp" -#include "KokkosTypes.hpp" -#include "cells.hpp" -#include "iBox.hpp" - - - -namespace pFlow -{ - -template< - typename executionSpace - > -class cellsWallLevel0 -: - 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; - - using iBoxType = iBox; - - class TagFindCellRange2{}; - -protected: - - // - box extent - real cellExtent_ = 0.5; - - // - number of triangle elements - int32 numElements_ = 0; - - // - number of points - int32 numPoints_ = 0; - - // - ref to vectices (borrowed) - ViewType1D vertices_; - - // - ref to points in the trisurface (borrowed) - ViewType1D points_; - - // cell range of element/triangle bounding box - ViewType1D elementBox_; - - - using tpPWContactSearch = Kokkos::TeamPolicy< - execution_space, - Kokkos::Schedule, - Kokkos::IndexType - >; - - using rpFindCellRange2Type = - Kokkos::RangePolicy>; - - - FUNCTION_H - void allocateArrays() - { - reallocNoInit( elementBox_, numElements_); - } - -public: - - TypeInfoNV("cellsWallLevel0"); - - INLINE_FUNCTION_HD - cellsWallLevel0(){} - - FUNCTION_H - cellsWallLevel0( - const Cells& ppCells, - real cellExtent, - int32 numPoints, - int32 numElements, - const ViewType1D& points, - const ViewType1D& vertices - ) - : - Cells(ppCells), - cellExtent_( max(cellExtent, 0.5 ) ), - numElements_(numElements), - numPoints_(numPoints), - vertices_(vertices), - points_(points) - { - - allocateArrays(); - } - - - // - host call - // reset triangle elements if they have changed - bool resetElements( - int32 numElements, - int32 numPoints, - ViewType1D& points, - ViewType1D& vertices ) - { - - numElements_ = numElements; - numPoints_ = numPoints; - points_ = points; - vertices_ = vertices; - - allocateArrays(); - - return true; - } - - INLINE_FUNCTION_HD - iBoxType elementBox(int32 i)const - { - return elementBox_[i]; - } - - INLINE_FUNCTION_HD - int32 numElements()const - { - return numElements_; - } - - - - template - bool broadSearch(PairsContainer& pairs, particleMapType& particleMap) - { - - // map walls onto the cells - this->build(); - - this->particleWallFindPairs(pairs, particleMap); - - return true; - } - - bool build() - { - Kokkos::parallel_for( - "cellsSimple::findcellrange2", - rpFindCellRange2Type(0,numElements_), - *this); - Kokkos::fence(); - return true; - } - - template - bool particleWallFindPairs(PairsContainer& pairs, particleMapType& particleMap) - { - - int32 getFull = 1; - - while (getFull) - { - - getFull = findPairsElementRangeCount(pairs, particleMap.getCellIterator(0)); - - 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); - - INFORMATION<<"Contact pair container capacity increased from "<< - oldCap << " to " - << pairs.capacity() <<" in cellsWallLevel0."< - int32 findPairsElementRangeCount(PairsContainer& pairs, CellIteratorType cellIter) - { - int32 getFull =0; - - const auto pwPairs = pairs; - const auto elementBox = elementBox_; - - Kokkos::parallel_reduce( - "cellsSimple::findPairsElementRangeModified2", - tpPWContactSearch(numElements_, Kokkos::AUTO), - LAMBDA_HD( - const typename tpPWContactSearch::member_type & teamMember, - int32& valueToUpdate){ - - const int32 iTri = teamMember.league_rank(); - - const auto triBox = elementBox[iTri]; - - int32 getFull2 = 0; - - auto bExtent = boxExtent(triBox); - int32 numCellBox = bExtent.x()*bExtent.y()*bExtent.z(); - - Kokkos::parallel_reduce( - Kokkos::TeamThreadRange( teamMember, numCellBox ), - [&] ( const int32 linIndex, int32 &innerUpdate ) - { - - CellType cell; - indexToCell(linIndex, triBox, cell); - - int32 n = cellIter.start(cell.x(),cell.y(),cell.z()); - - while( n>-1) - { - // id is wall id the pair is (particle id, wall id) - if( pairs.insert(static_cast(n), iTri) < 0 ) - innerUpdate++; - n = cellIter.getNext(n); - } - - }, - getFull2 - ); - - if ( teamMember.team_rank() == 0 ) valueToUpdate += getFull2; - }, - getFull - ); - - return getFull; - } - - INLINE_FUNCTION_HD - void operator()(TagFindCellRange2, int32 i) const - { - auto v = vertices_[i]; - auto p1 = points_[v.x()]; - auto p2 = points_[v.y()]; - auto p3 = points_[v.z()]; - - realx3 minP, maxP; - - this->extendBox(p1, p2, p3, cellExtent_, minP, maxP); - elementBox_[i] = iBoxType(this->pointIndex(minP), this->pointIndex(maxP)); - - } - -}; // cellsWallLevel0 - -} // pFlow - - -#endif // __cellsWallLevel0_hpp__ diff --git a/src/Interaction/contactSearch/wallMappings/cellsWallLevels.hpp b/src/Interaction/contactSearch/wallMappings/cellsWallLevels.hpp deleted file mode 100644 index 9f2c80ec..00000000 --- a/src/Interaction/contactSearch/wallMappings/cellsWallLevels.hpp +++ /dev/null @@ -1,152 +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 __cellsWallLevels_hpp__ -#define __cellsWallLevels_hpp__ - -#include "cellsWallLevel0.hpp" - -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: - - TypeInfoNV("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: - - TypeInfoNV("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/Interaction/interaction/interaction.cpp b/src/Interaction/interaction/interaction.cpp index 8d368c29..ec8938d7 100644 --- a/src/Interaction/interaction/interaction.cpp +++ b/src/Interaction/interaction/interaction.cpp @@ -19,7 +19,9 @@ Licence: -----------------------------------------------------------------------------*/ #include "interaction.hpp" - +#include "particles.hpp" +#include "vocabs.hpp" +#include "geometry.hpp" pFlow::interaction::interaction ( @@ -28,26 +30,16 @@ pFlow::interaction::interaction const geometry& geom ) : - demInteraction(control, control.caseSetup().path()+interactionFile__), - interactionBase(prtcl, geom), - fileDict_(control.caseSetup().emplaceObject( - objectFile( - interactionFile__, - "", - objectFile::READ_ALWAYS, - objectFile::WRITE_NEVER), - interactionFile__, - true )) + property(interactionFile__, &control.caseSetup()), + observer + ( + &prtcl.dynPointStruct(), + msg_ + ), + demComponent("interaction", control), + particles_(prtcl), + geometry_(geom) { - this->subscribe(prtcl.pStruct()); - - contactSearch_ = contactSearch::create( - fileDict_.subDict("contactSearch"), - this->control().domain(), - prtcl, - geom, - timers() - ); } @@ -79,14 +71,14 @@ pFlow::uniquePtr pFlow::interaction::create clType); - REPORT(1)<< "Selecting interaction model..."< contactSearch_ = nullptr; + static inline + const message msg_ = message::ITEM_REARRANGE; public: @@ -65,11 +65,10 @@ public: const particles& prtcl, const geometry& geom ); + ~interaction() override = default; - - virtual ~interaction() = default; - - create_vCtor( + create_vCtor + ( interaction, systemControl, ( @@ -78,22 +77,8 @@ public: const geometry& geom ), (control, prtcl, geom) - ); + ); - auto& contactSearchPtr() - { - return contactSearch_; - } - - auto& contactSearchRef() - { - return contactSearch_(); - } - - const auto& fileDict()const - { - return fileDict_; - } static uniquePtr create( diff --git a/src/Interaction/sphereInteraction/sphereInteraction.cpp b/src/Interaction/sphereInteraction/sphereInteraction.cpp index 0dacee74..1cb1472e 100644 --- a/src/Interaction/sphereInteraction/sphereInteraction.cpp +++ b/src/Interaction/sphereInteraction/sphereInteraction.cpp @@ -18,22 +18,18 @@ Licence: -----------------------------------------------------------------------------*/ -template< - typename contactForceModel, - typename geometryMotionModel, - template class contactListType > -bool pFlow::sphereInteraction:: - createSphereInteraction() +template class cLT> +bool pFlow::sphereInteraction::createSphereInteraction() { - realVector_D rhoD(this->densities()); + realVector_D rhoD("densities", this->densities()); - auto modelDict = this->fileDict().subDict("model"); + auto modelDict = this->subDict("model"); - REPORT(1)<<"Createing contact force model . . ."<( this->numMaterials(), - rhoD.deviceVector(), + rhoD.deviceView(), modelDict ); @@ -47,17 +43,9 @@ bool pFlow::sphereInteraction class contactListType > -bool pFlow::sphereInteraction:: - sphereSphereInteraction() +template class cLT> +bool pFlow::sphereInteraction::sphereSphereInteraction() { - - - auto lastItem = ppContactList_().loopCount(); // create the kernel functor @@ -66,17 +54,17 @@ bool pFlow::sphereInteractiondt(), this->forceModel_(), ppContactList_(), // to be read - sphParticles_.diameter().deviceVectorAll(), - sphParticles_.propertyId().deviceVectorAll(), - sphParticles_.pointPosition().deviceVectorAll(), - sphParticles_.velocity().deviceVectorAll(), - sphParticles_.rVelocity().deviceVectorAll(), - sphParticles_.contactForce().deviceVectorAll(), - sphParticles_.contactTorque().deviceVectorAll() + sphParticles_.diameter().deviceViewAll(), + sphParticles_.propertyId().deviceViewAll(), + sphParticles_.pointPosition().deviceViewAll(), + sphParticles_.velocity().deviceViewAll(), + sphParticles_.rVelocity().deviceViewAll(), + sphParticles_.contactForce().deviceViewAll(), + sphParticles_.contactTorque().deviceViewAll() ); Kokkos::parallel_for( - "", + "ppInteraction", rpPPInteraction(0,lastItem), ppInteraction ); @@ -87,34 +75,32 @@ bool pFlow::sphereInteraction class contactListType > -bool pFlow::sphereInteraction:: - sphereWallInteraction() +template class cLT> +bool pFlow::sphereInteraction::sphereWallInteraction() { - int32 lastItem = pwContactList_().loopCount(); + uint32 lastItem = pwContactList_().loopCount(); + uint32 iter = this->currentIter(); real t = this->currentTime(); + real dt = this->dt(); pFlow::sphereInteractionKernels::pwInteractionFunctor pwInteraction( - this->dt(), + dt, this->forceModel_(), pwContactList_(), geometryMotion_.getTriangleAccessor(), - geometryMotion_.getModel(t) , - sphParticles_.diameter().deviceVectorAll() , - sphParticles_.propertyId().deviceVectorAll(), - sphParticles_.pointPosition().deviceVectorAll(), - sphParticles_.velocity().deviceVectorAll(), - sphParticles_.rVelocity().deviceVectorAll() , - sphParticles_.contactForce().deviceVectorAll(), - sphParticles_.contactTorque().deviceVectorAll() , - geometryMotion_.triMotionIndex().deviceVectorAll(), - geometryMotion_.propertyId().deviceVectorAll(), - geometryMotion_.contactForceWall().deviceVectorAll() + geometryMotion_.getModel(iter, t, dt) , + sphParticles_.diameter().deviceViewAll() , + sphParticles_.propertyId().deviceViewAll(), + sphParticles_.pointPosition().deviceViewAll(), + sphParticles_.velocity().deviceViewAll(), + sphParticles_.rVelocity().deviceViewAll() , + sphParticles_.contactForce().deviceViewAll(), + sphParticles_.contactTorque().deviceViewAll() , + geometryMotion_.triMotionIndex().deviceViewAll(), + geometryMotion_.propertyId().deviceViewAll(), + geometryMotion_.contactForceWall().deviceViewAll() ); Kokkos::parallel_for( @@ -127,4 +113,123 @@ bool pFlow::sphereInteraction class cLT> +pFlow::sphereInteraction::sphereInteraction +( + systemControl& control, + const particles& prtcl, + const geometry& geom +) +: + interaction(control, prtcl, geom), + geometryMotion_(dynamic_cast(geom)), + sphParticles_(dynamic_cast(prtcl)), + ppInteractionTimer_("sphere-sphere interaction", &this->timers()), + pwInteractionTimer_("sphere-wall interaction", &this->timers()), + contactListTimer_("contact list management", &this->timers()), + contactListTimer0_("contact list clear", &this->timers()) +{ + contactSearch_ = contactSearch::create( + subDict("contactSearch"), + prtcl.thisDomain().domainBox(), + prtcl, + geom, + timers()); + + + if(!createSphereInteraction()) + { + fatalExit; + } +} + +template class cLT> +bool pFlow::sphereInteraction::beforeIteration() +{ + return true; +} + +template class cLT> +bool pFlow::sphereInteraction::iterate() +{ + + auto iter = this->currentIter(); + auto t = this->currentTime(); + auto dt = this->dt(); + + //output<<"iter, t, dt "<< iter<<" "<< t << " "< class cLT> +bool pFlow::sphereInteraction::afterIteration() +{ + return true; +} + +template class cLT> +bool pFlow::sphereInteraction::hearChanges +( + real t, + real dt, + uint32 iter, + const message& msg, + const anyList& varList +) +{ + if(msg.equivalentTo(message::ITEM_REARRANGE)) + { + notImplementedFunction; + } + return true; +} + diff --git a/src/Interaction/sphereInteraction/sphereInteraction.hpp b/src/Interaction/sphereInteraction/sphereInteraction.hpp index b256c452..56c00047 100644 --- a/src/Interaction/sphereInteraction/sphereInteraction.hpp +++ b/src/Interaction/sphereInteraction/sphereInteraction.hpp @@ -23,7 +23,6 @@ Licence: #include "interaction.hpp" #include "sphereParticles.hpp" - #include "sphereInteractionKernels.hpp" namespace pFlow @@ -43,22 +42,19 @@ public: using ContactForceModel = contactForceModel; - using MotionModel = typename geometryMotionModel::MotionModel; + using MotionModel = typename geometryMotionModel::MotionModel; - using ModelStorage = typename ContactForceModel::contactForceStorage; + using ModelStorage = typename ContactForceModel::contactForceStorage; - using IdType = typename interaction::IdType; + using IdType = uint32; - using IndexType = typename interaction::IndexType; - - using ExecutionSpace = typename interaction::ExecutionSpace; + using IndexType = uint32; using ContactListType = - contactListType; + contactListType; - using PairsContainerType= typename contactSearch::PairContainerType; -protected: +private: /// const reference to geometry const GeometryMotionModel& geometryMotion_; @@ -66,15 +62,17 @@ protected: /// const reference to particles const sphereParticles& sphParticles_; + /// contact search object for pp and pw interactions + uniquePtr contactSearch_ = nullptr; /// contact force model - uniquePtr forceModel_ = nullptr; + uniquePtr forceModel_ = nullptr; /// contact list for particle-particle interactoins (keeps the history) - uniquePtr ppContactList_ = nullptr; + uniquePtr ppContactList_ = nullptr; /// contact list for particle-wall interactions (keeps the history) - uniquePtr pwContactList_ = nullptr; + uniquePtr pwContactList_ = nullptr; /// timer for particle-particle interaction computations Timer ppInteractionTimer_; @@ -82,42 +80,39 @@ protected: /// timer for particle-wall interaction computations Timer pwInteractionTimer_; + Timer contactListTimer_; + + Timer contactListTimer0_; + bool createSphereInteraction(); - bool managePPContactLists(); + bool sphereSphereInteraction(); - bool managePWContactLists(); + bool sphereWallInteraction(); + + //bool managePPContactLists(); + + //bool managePWContactLists(); /// range policy for p-p interaction execution using rpPPInteraction = - Kokkos::RangePolicy, Kokkos::Schedule>; + Kokkos::RangePolicy, Kokkos::Schedule>; /// range policy for p-w interaction execution using rpPWInteraction = rpPPInteraction; public: - TypeInfoTemplate3("sphereInteraction", ContactForceModel, MotionModel, ContactListType); - - // constructor + TypeInfoTemplate13("sphereInteraction", ContactForceModel, MotionModel, ContactListType); + /// Constructor from components sphereInteraction( systemControl& control, const particles& prtcl, - const geometry& geom) - : - interaction(control, prtcl, geom), - geometryMotion_(dynamic_cast(geom)), - sphParticles_(dynamic_cast(prtcl)), - ppInteractionTimer_("sphere-sphere interaction", &this->timers()), - pwInteractionTimer_("sphere-wall interaction", &this->timers()) - { - if(!createSphereInteraction()) - { - fatalExit; - } - } + const geometry& geom); + + /// Add virtual constructor add_vCtor ( interaction, @@ -125,97 +120,25 @@ public: systemControl ); + /// This is called in time loop, before iterate. (overriden from demComponent) + bool beforeIteration() override; + + /// This is called in time loop. Perform the main calculations + /// when the component should evolve along time. (overriden from demComponent) + bool iterate() override; + + /// This is called in time loop, after iterate. (overriden from demComponent) + bool afterIteration() override; - bool beforeIteration() override - { - return true; - } - - - bool iterate() override - { - -//Info<<"before contact search"<contactSearch_) - { - - if( this->contactSearch_().ppEnterBroadSearch()) - { - //Info<<" before ppEnterBroadSearch"<contactSearch_().pwEnterBroadSearch()) - { - //Info<<" before pwEnterBroadSearch"<contactSearch_().ppPerformedBroadSearch()) - { - //Info<<" before afterBroadSearch"<contactSearch_().pwPerformedBroadSearch()) - { - //Info<<" before pwContactList_().afterBroadSearch()"< diam_; - deviceViewType1D propId_; + deviceViewType1D propId_; deviceViewType1D pos_; deviceViewType1D lVel_; deviceViewType1D rVel_; @@ -55,7 +55,7 @@ struct ppInteractionFunctor ContactForceModel forceModel, ContactListType tobeFilled, deviceViewType1D diam, - deviceViewType1D propId, + deviceViewType1D propId, deviceViewType1D pos, deviceViewType1D lVel, deviceViewType1D rVel, @@ -75,7 +75,7 @@ struct ppInteractionFunctor {} INLINE_FUNCTION_HD - void operator()(const int32 n)const + void operator()(const uint32 n)const { if(!tobeFilled_.isValid(n))return; @@ -181,14 +181,14 @@ struct pwInteractionFunctor MotionModel motionModel_; deviceViewType1D diam_; - deviceViewType1D propId_; + deviceViewType1D propId_; deviceViewType1D pos_; deviceViewType1D lVel_; deviceViewType1D rVel_; deviceViewType1D cForce_; deviceViewType1D cTorque_; - deviceViewType1D wTriMotionIndex_; - deviceViewType1D wPropId_; + deviceViewType1D wTriMotionIndex_; + deviceViewType1D wPropId_; deviceViewType1D wCForce_; @@ -199,14 +199,14 @@ struct pwInteractionFunctor TraingleAccessor triangles, MotionModel motionModel , deviceViewType1D diam , - deviceViewType1D propId, + deviceViewType1D propId, deviceViewType1D pos , deviceViewType1D lVel, deviceViewType1D rVel, deviceViewType1D cForce, deviceViewType1D cTorque , - deviceViewType1D wTriMotionIndex, - deviceViewType1D wPropId, + deviceViewType1D wTriMotionIndex, + deviceViewType1D wPropId, deviceViewType1D wCForce) : dt_(dt), diff --git a/src/Interaction/sphereInteraction/sphereInteractions.cpp b/src/Interaction/sphereInteraction/sphereInteractions.cpp index 569fa7ad..1bd0da9c 100644 --- a/src/Interaction/sphereInteraction/sphereInteractions.cpp +++ b/src/Interaction/sphereInteraction/sphereInteractions.cpp @@ -28,22 +28,22 @@ Licence: template class pFlow::sphereInteraction< pFlow::cfModels::limitedLinearNormalRolling, - pFlow::fixedGeometry, + pFlow::stationaryGeometry, pFlow::unsortedContactList>; template class pFlow::sphereInteraction< pFlow::cfModels::limitedLinearNormalRolling, - pFlow::fixedGeometry, + pFlow::stationaryGeometry, pFlow::sortedContactList>; template class pFlow::sphereInteraction< pFlow::cfModels::nonLimitedLinearNormalRolling, - pFlow::fixedGeometry, + pFlow::stationaryGeometry, pFlow::unsortedContactList>; template class pFlow::sphereInteraction< pFlow::cfModels::nonLimitedLinearNormalRolling, - pFlow::fixedGeometry, + pFlow::stationaryGeometry, pFlow::sortedContactList>; template class pFlow::sphereInteraction< @@ -88,7 +88,7 @@ template class pFlow::sphereInteraction< pFlow::sortedContactList>; -template class pFlow::sphereInteraction< +/*template class pFlow::sphereInteraction< pFlow::cfModels::limitedLinearNormalRolling, pFlow::multiRotationAxisMotionGeometry, pFlow::unsortedContactList>; @@ -106,28 +106,28 @@ template class pFlow::sphereInteraction< template class pFlow::sphereInteraction< pFlow::cfModels::nonLimitedLinearNormalRolling, pFlow::multiRotationAxisMotionGeometry, - pFlow::sortedContactList>; + pFlow::sortedContactList>;*/ /// non-linear models template class pFlow::sphereInteraction< pFlow::cfModels::limitedNonLinearNormalRolling, - pFlow::fixedGeometry, + pFlow::stationaryGeometry, pFlow::unsortedContactList>; template class pFlow::sphereInteraction< pFlow::cfModels::limitedNonLinearNormalRolling, - pFlow::fixedGeometry, + pFlow::stationaryGeometry, pFlow::sortedContactList>; template class pFlow::sphereInteraction< pFlow::cfModels::nonLimitedNonLinearNormalRolling, - pFlow::fixedGeometry, + pFlow::stationaryGeometry, pFlow::unsortedContactList>; template class pFlow::sphereInteraction< pFlow::cfModels::nonLimitedNonLinearNormalRolling, - pFlow::fixedGeometry, + pFlow::stationaryGeometry, pFlow::sortedContactList>; template class pFlow::sphereInteraction< @@ -171,7 +171,7 @@ template class pFlow::sphereInteraction< pFlow::vibratingMotionGeometry, pFlow::sortedContactList>; -template class pFlow::sphereInteraction< +/*template class pFlow::sphereInteraction< pFlow::cfModels::limitedNonLinearNormalRolling, pFlow::multiRotationAxisMotionGeometry, pFlow::unsortedContactList>; @@ -189,28 +189,28 @@ template class pFlow::sphereInteraction< template class pFlow::sphereInteraction< pFlow::cfModels::nonLimitedNonLinearNormalRolling, pFlow::multiRotationAxisMotionGeometry, - pFlow::sortedContactList>; + pFlow::sortedContactList>;*/ // - nonLinearMod models template class pFlow::sphereInteraction< pFlow::cfModels::limitedNonLinearModNormalRolling, - pFlow::fixedGeometry, + pFlow::stationaryGeometry, pFlow::unsortedContactList>; template class pFlow::sphereInteraction< pFlow::cfModels::limitedNonLinearModNormalRolling, - pFlow::fixedGeometry, + pFlow::stationaryGeometry, pFlow::sortedContactList>; template class pFlow::sphereInteraction< pFlow::cfModels::nonLimitedNonLinearModNormalRolling, - pFlow::fixedGeometry, + pFlow::stationaryGeometry, pFlow::unsortedContactList>; template class pFlow::sphereInteraction< pFlow::cfModels::nonLimitedNonLinearModNormalRolling, - pFlow::fixedGeometry, + pFlow::stationaryGeometry, pFlow::sortedContactList>; template class pFlow::sphereInteraction< @@ -255,7 +255,7 @@ template class pFlow::sphereInteraction< pFlow::sortedContactList>; -template class pFlow::sphereInteraction< +/*template class pFlow::sphereInteraction< pFlow::cfModels::limitedNonLinearModNormalRolling, pFlow::multiRotationAxisMotionGeometry, pFlow::unsortedContactList>; @@ -274,3 +274,4 @@ template class pFlow::sphereInteraction< pFlow::cfModels::nonLimitedNonLinearModNormalRolling, pFlow::multiRotationAxisMotionGeometry, pFlow::sortedContactList>; +*/ \ No newline at end of file