xfygogo
/
phasicFlow
mirror of https://github.com/PhasicFlow/phasicFlow.git
1
0
Fork 0
Code Issues Packages Projects Releases Wiki Activity

Interaction folder is finalized with adjustable box and tested for rotating drum

This commit is contained in:
Hamidreza Norouzi 2024-03-24 01:46:19 -07:00
parent be56d8ee2e
commit d21b7276e1
48 changed files with 2252 additions and 3260 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

5
src/Interaction/CMakeLists.txt
View File

@ -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

15
src/Interaction/Models/contactForce/linearCF.hpp
View File

@ -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<linearProperties> prop(nElem);
Vector<linearProperties> 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,

12
src/Interaction/Models/contactForce/nonLinearCF.hpp
View File

@ -121,7 +121,7 @@ protected:
}
realVector etha_n(nElem);
realVector etha_n("etha_n",nElem);
ForAll(i , en)
{
@ -137,7 +137,7 @@ protected:
}
Vector<nonLinearProperties> prop(nElem);
Vector<nonLinearProperties> 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,

10
src/Interaction/Models/contactForce/nonLinearMod.hpp
View File

@ -121,7 +121,7 @@ protected:
}
Vector<nonLinearProperties> prop(nElem);
Vector<nonLinearProperties> 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,

8
src/Interaction/Models/rolling/normalRolling.hpp
View File

@ -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,

22
src/Interaction/contactLists/sortedContactList.hpp
View File

@ -53,7 +53,7 @@ protected:
ViewType1D<ValueType,ExecutionSpace> values_;
int32 size0_ = 0;
uint32 size0_ = 0;
ViewType1D<PairType,ExecutionSpace> sortedPairs0_;
@ -73,7 +73,7 @@ protected:
using rpReFillPairs = Kokkos::RangePolicy<
ExecutionSpace,
Kokkos::Schedule<Kokkos::Static>,
Kokkos::IndexType<int32>,
Kokkos::IndexType<uint32>,
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];

40
src/Interaction/contactLists/sortedPairs.hpp
View File

@ -52,24 +52,24 @@ public:
{
using PairType = typename sortedPairs::PairType;
int32 size_;
uint32 size_;
ViewType1D<PairType,ExecutionSpace> 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<size_; }
bool isValid(uint32 i)const { return i<size_; }
INLINE_FUNCTION_HD
PairType getPair(int i)const { return sortedParis_[i]; }
PairType getPair(uint32 i)const { return sortedParis_[i]; }
INLINE_FUNCTION_HD
bool getPair(int32 i, PairType& pair)const {
bool getPair(uint32 i, PairType& pair)const {
if(i<size_) {
pair = sortedParis_[i];
return true;
@ -85,22 +85,22 @@ public:
protected:
/// size of pair list
int32 size_ = 0;
uint32 size_ = 0;
ViewType1D<int32,ExecutionSpace> flags_;
ViewType1D<uint32,ExecutionSpace> flags_;
ViewType1D<PairType,ExecutionSpace> sortedPairs_;
using rpFillFlag = Kokkos::RangePolicy<
ExecutionSpace,
Kokkos::Schedule<Kokkos::Static>,
Kokkos::IndexType<int32>,
Kokkos::IndexType<uint32>,
TagFillFlag >;
using rpFillPairs = Kokkos::RangePolicy<
ExecutionSpace,
Kokkos::Schedule<Kokkos::Static>,
Kokkos::IndexType<int32>,
Kokkos::IndexType<uint32>,
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])

16
src/Interaction/contactLists/unsortedContactList.hpp
View File

@ -72,7 +72,7 @@ protected:
using rpFillPairs = Kokkos::RangePolicy<
ExecutionSpace,
Kokkos::Schedule<Kokkos::Static>,
Kokkos::IndexType<int32>,
Kokkos::IndexType<uint32>,
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];
}

34
src/Interaction/contactLists/unsortedPairs.hpp
View File

@ -41,7 +41,7 @@ public:
using memory_space = typename ExecutionSpace::memory_space;
using PairType = kPair<idType,idType>;
using PairType = Pair<idType,idType>;
using ContainerType = unorderedSet<PairType, ExecutionSpace>;
@ -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();

199
src/Interaction/contactSearch/ContactSearch/ContactSearch.hpp
View File

@ -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> class BaseMethod,
template<class> 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<ParticleContactSearchType> particleContactSearch_ = nullptr;
uniquePtr<WallMappingType> wallMapping_ = nullptr;
uniquePtr<SearchMethodType> 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<word>("method");
auto wmMethod = dict().getVal<word>("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<ParticleContactSearchType>
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<SearchMethodType>
(
nbDict,
this->domain(),
this->domainBox(),
minD,
maxD,
position,
diam
);
REPORT(2)<<"Contact search algorithm for particle-particle is "<<
greenText(particleContactSearch_().typeName())<<endREPORT;
auto wmDict = dict().subDict(wmMethod+"Info");
int32 wnPoints = this->Geometry().numPoints();
int32 wnTri = this->Geometry().size();
const auto& wPoints = this->Geometry().points().deviceVectorAll();
const auto& wVertices = this->Geometry().vertices().deviceVectorAll();
wallMapping_ =
makeUnique<WallMappingType>(
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())<<END_REPORT;
}
@ -130,116 +105,60 @@ public:
contactSearch,
ContactSearch,
dictionary);
bool broadSearch(
PairContainerType& ppPairs,
PairContainerType& pwPairs,
uint32 iter,
real t,
real dt,
csPairContainerType& ppPairs,
csPairContainerType& pwPairs,
bool force = false) override
{
ppTimer().start();
const auto& position = Particles().pointPosition().deviceViewAll();
const auto& flags = Particles().dynPointStruct().activePointsMaskDevice();
const auto& diam = Particles().boundingSphere().deviceViewAll();
if(particleContactSearch_)
if( !ppwContactSearch_().broadSearch(
iter,
t,
dt,
ppPairs,
pwPairs,
position,
flags,
diam,
force) )
{
auto activeRange = this->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;
}*/
};
}

10
src/Interaction/contactSearch/ContactSearch/ContactSearchs.cpp
View File

@ -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<pFlow::NBS, pFlow::cellMapping>;
template class pFlow::ContactSearch<pFlow::multiGridNBS, pFlow::multiGridMapping>;
template class pFlow::ContactSearch<pFlow::NBS>;
//template class pFlow::ContactSearch<pFlow::multiGridNBS, pFlow::multiGridMapping>;

25
src/Interaction/contactSearch/contactSearch/contactSearch.cpp
View File

@ -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> pFlow::contactSearch::create(
const geometry& geom,
Timers& timers)
{
word baseMethName = dict.getVal<word>("method");
word baseMethName = dict.getVal<word>("method");
word wallMethod = dict.getVal<word>("wallMapping");
auto model = angleBracketsNames("ContactSearch", baseMethName);
auto model = angleBracketsNames2("ContactSearch", baseMethName, wallMethod);
REPORT(1)<<"Selecting contact search model . . ."<<endREPORT;
REPORT(1)<<"Selecting contact search model . . ."<<END_REPORT;
if( dictionaryvCtorSelector_.search(model))
{
auto objPtr = dictionaryvCtorSelector_[model] (dict, domain, prtcl, geom, timers);
REPORT(2)<<"Model "<< greenText(model)<<" is created."<<endREPORT;
REPORT(2)<<"Model "<< Green_Text(model)<<" is created."<< END_REPORT;
return objPtr;
}
else

83
src/Interaction/contactSearch/contactSearch/contactSearch.hpp
View File

@ -23,43 +23,36 @@ Licence:
#define __contactSearch_hpp__
#include "interactionBase.hpp"
#include "unsortedPairs.hpp"
#include "box.hpp"
#include "contactSearchGlobals.hpp"
#include "dictionary.hpp"
#include "virtualConstructor.hpp"
#include "Timer.hpp"
namespace pFlow
{
// - forward
class box;
class particles;
class geometry;
class contactSearch
:
public interactionBase
{
public:
using IdType = typename interactionBase::IdType;
private:
using IndexType = typename interactionBase::IndexType;
const box& domainBox_;
using ExecutionSpace = typename interactionBase::ExecutionSpace;
const particles& particles_;
using PairContainerType = unsortedPairs<ExecutionSpace, IdType>;
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<contactSearch> create(

11
src/Interaction/contactSearch/contactSearchFunctions.hpp
View File

@ -96,17 +96,6 @@ void indexToCell(const indexType idx, const iBox<cellIndexType>& box, triple<cel
cell+= box.minPoint();
}
INLINE_FUNCTION_HD
bool sphereSphereCheck(const realx3& p1, const realx3 p2, real d1, real d2)
{
return length(p2-p1) < 0.5*(d2+d1);
}
}
#endif //__broadSearchFunctions_hpp__

21
src/Interaction/interaction/interactionTypes.hpp → src/Interaction/contactSearch/contactSearchGlobals.hpp
View File

@ -18,28 +18,21 @@ Licence:
-----------------------------------------------------------------------------*/
#ifndef __interactionTypes_hpp__
#define __interactionTypes_hpp__
#include "types.hpp"
#include "unsortedPairs.hpp"
#ifndef __contactSearchGlobals_hpp__
#define __contactSearchGlobals_hpp__
namespace pFlow
{
using csExecutionSpace = DefaultExecutionSpace;
// always use signed types
using CELL_INDEX_TYPE = int32;
using csIdType = uint32;
using ID_TYPE = int32;
//constexpr int32 minCellIndex = largestNegative<CELL_INDEX_TYPE>();
//constexpr int32 maxCellIndex = largestPositive<CELL_INDEX_TYPE>();
using csPairContainerType = unsortedPairs<DefaultExecutionSpace, uint32>;
}
#endif //__interactionTypes_hpp__
#endif

209
src/Interaction/contactSearch/methods/NBS.hpp
View File

@ -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<typename executionSpace>
class NBS
{
public:
using NBSLevel0Type = NBSLevel0<executionSpace>;
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<realx3, memory_space>& position,
const ViewType1D<real, memory_space>& diam)
:
sizeRatio_(
max(
dict.getValOrSet<real>("sizeRatio", 1.0),
1.0
)),
updateFrequency_(
max(
dict.getValOrSet<int32>("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<cellIterator> getCellIteratorLevels()
{
return Vector<cellIterator>("cellIterator", 1, NBSLevel0_.getCellIterator());
}
auto getCellIterator(int32 lvl)const
{
return NBSLevel0_.getCellIterator();
}
int32 numLevels()const
{
return 1;
}
Vector<Cells> getCellsLevels()const
{
return Vector<Cells>("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<typename PairsContainer>
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<typename PairsContainer, typename IncludeFunction>
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

82
src/Interaction/contactSearch/methods/NBSCrossLoop.hpp
View File

@ -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];
}

127
src/Interaction/contactSearch/methods/NBSLevel.hpp
View File

@ -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<typename executionSpace>
class
NBSLevel
:
public NBSLevel0<executionSpace>
{
public:
using NBSLevel0Type = NBSLevel0<executionSpace>;
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<typename exeSpace>
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<realx3, memory_space>& position,
const ViewType1D<real, memory_space>& 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<typename PairsContainer>
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<int32>(upperLevel-lowerLevel));
return ind/a;
}
}
#endif

240
src/Interaction/contactSearch/methods/NBSLevel0.hpp
View File

@ -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<typename executionSpace>
class NBSLevel0
:
public mapperNBS<executionSpace>
{
public:
using MapperType = mapperNBS<executionSpace>;
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<real, memory_space> diameter_;
using mdrPolicyFindPairs =
Kokkos::MDRangePolicy<
Kokkos::Rank<3>,
Kokkos::Schedule<Kokkos::Dynamic>,
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<realx3, memory_space>& position,
const ViewType1D<real, memory_space>& diam)
:
MapperType(domain, cellSize, position),
diameter_(diam)
{}
NBSLevel0(
const box& domain,
int32 nx,
int32 ny,
int32 nz,
const ViewType1D<realx3, memory_space>& position,
const ViewType1D<real, memory_space>& diam)
:
MapperType(domain, nx, ny, nz, position),
diameter_(diam)
{ }
NBSLevel0(
const box& domain,
real cellSize,
real sizeRatio,
const ViewType1D<realx3, memory_space>& position,
const ViewType1D<real, memory_space>& 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<typename PairsContainer>
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<typename PairsContainer, typename IncludeFunction>
bool broadSearch(PairsContainer& pairs, range activeRange, IncludeFunction incld)
{
this->build(activeRange, incld);
findPairs(pairs);
return true;
}
template<typename PairsContainer>
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 "<<pairs.capacity()<<" in NBSLevel0."<<endINFO;
}
Kokkos::fence();
}
return true;
}
template<typename PairsContainer>
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__

438
src/Interaction/contactSearch/methods/NBSLevels.hpp
View File

@ -1,438 +0,0 @@
#ifndef __NBSLevels_hpp__
#define __NBSLevels_hpp__
#include "NBSLevel.hpp"
#include "NBSLevel0.hpp"
#include "KokkosTypes.hpp"
namespace pFlow
{
template<typename executionSpace>
class NBSLevels
{
public:
using NBSLevelType = NBSLevel<executionSpace>;
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<real,real>;
protected:
real minSize_;
real maxSize_;
int32 numLevels_=1;
Vector<NBSLevelType> nbsLevels_;
ViewType1D<realRange, memory_space> sizeRangeLevels_;
ViewType1D<realRange, HostSpace> sizeRangeLevelsHost_;
ViewType1D<real, memory_space> maxSizeLevels_;
ViewType1D<real, HostSpace> maxSizeLevelsHost_;
ViewType1D<int8, memory_space> particleLevel_;
range activeRange_{0,0};
using rangePolicyType =
Kokkos::RangePolicy<
Kokkos::IndexType<int32>,
Kokkos::Schedule<Kokkos::Static>,
execution_space>;
int32 setNumLevels()
{
int32 maxOvermin = static_cast<int32>(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<<endl;
fatalExit;
}
return -1;
}
bool setDiameterRange(real sizeRatio)
{
real lvl_maxD = sizeRatio* maxSize_;
real lvl_minD = lvl_maxD/2.0;
for(int32 lvl=numLevels_-1; lvl>=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."<<endREPORT;
for(int32 lvl=0; lvl<numLevels_; lvl++)
{
REPORT(3)<<"Cell gird No "<< yellowText(lvl)<<" with size range ("
<<sizeRangeLevelsHost_[lvl].first<<","<<sizeRangeLevelsHost_[lvl].second<<"]."<<endREPORT;
}
return true;
}
bool initLevels(
const box& domain,
real sizeRatio,
const ViewType1D<realx3, memory_space>& position,
const ViewType1D<real, memory_space>& diam)
{
for(int32 lvl = 0; lvl<numLevels_; lvl++)
{
nbsLevels_[lvl] = NBSLevelType(
lvl,
domain,
maxSizeLevelsHost_[lvl]*sizeRatio,
sizeRatio,
position,
diam );
}
auto next0 = nbsLevels_[0].next();
for(int32 lvl=1; lvl<numLevels_; lvl++)
{
nbsLevels_[lvl].setNext(next0);
}
return true;
}
void manageAllocateNext(range active)
{
activeRange_ = active;
if(activeRange_.second > nbsLevels_[0].capacity())
{
nbsLevels_[0].checkAllocateNext(activeRange_.second);
auto next0 = nbsLevels_[0].next();
for(int32 lvl=1; lvl<numLevels_; lvl++)
{
nbsLevels_[lvl].setNext(next0);
}
}
}
void nullify( range active)
{
for(int32 lvl=0; lvl<numLevels_; lvl++)
{
nbsLevels_[lvl].nullify(active);
}
}
public:
NBSLevels(
const box& domain,
real minSize,
real maxSize,
real sizeRatio,
const ViewType1D<realx3, memory_space>& position,
const ViewType1D<real, memory_space>& 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<typename PairsContainer>
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 "<<pairs.capacity()<<" in NBSLevels."<<endINFO;
}
Kokkos::fence();
}
return true;
}
template<typename PairsContainer>
INLINE_FUNCTION_H
int32 findPairsCount(PairsContainer& pairs)
{
int32 notInsertedCount = 0;
for(int32 lvl=0; lvl<numLevels_; lvl++)
{
// the same level
notInsertedCount+= nbsLevels_[lvl].findPairsCount(pairs);
for(int32 crsLvl = lvl+1; crsLvl<numLevels_; crsLvl++)
{
// cross levels
notInsertedCount+=
nbsLevels_[lvl].findPairsCountCross(pairs, nbsLevels_[crsLvl]);
}
}
return notInsertedCount;
}
INLINE_FUNCTION_H
void build(range activeRange)
{
// 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){
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<typename IncludeFunction>
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<maxLvl; lvl++)
{
if( sizeRatio*diameter[i]<= maxSizes[lvl] )
{
particleLevel[i] = lvl;
return;
}
}
particleLevel[i] = static_cast<int8>(-1);
});
Kokkos::fence();
return true;
}
}; //NBSLevels
}
#endif

101
src/Interaction/contactSearch/methods/NBSLoop.hpp
View File

@ -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];
}

105
src/Interaction/interaction/interactionBase.hpp → src/Interaction/contactSearch/methods/cellBased/NBS/NBS.cpp
View File

@ -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<realx3> &position,
const pFlagTypeDevice &flags,
const deviceViewType1D<real> &diam,
uint32 nWallPoints,
uint32 nWallElements,
const ViewType1D<realx3,memory_space>& wallPoints,
const ViewType1D<uint32x3,memory_space>& wallVertices
)
:
particleWallContactSearchs<NBS>(
dict,
domainBox,
minBSSize,
maxBSSize,
position,
flags,
diam),
sizeRatio_(max(dict.getVal<real>("sizeRatio"), 1.0)),
cellExtent_(max(dict.getVal<real>("cellExtent"), 0.5)),
adjustableBox_(dict.getVal<Logical>("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__

148
src/Interaction/contactSearch/methods/cellBased/NBS/NBS.hpp Normal file
View File

@ -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<NBS>
{
public:
using CellIterator = typename NBSLevel0::CellIterator;
private:
real sizeRatio_ = 1.0;
real cellExtent_ = 0.5;
Logical adjustableBox_;
NBSLevel0 NBSLevel0_;
cellsWallLevel0 cellsWallLevel0_;
protected:
friend particleWallContactSearchs<NBS>;
bool impl_broadSearch
(
csPairContainerType& ppPairs,
csPairContainerType& pwPairs,
const deviceViewType1D<realx3>& pointPos,
const pFlagTypeDevice& flags,
const deviceViewType1D<real>& diameter
)
{
bool searchBoxChanged;
if( !NBSLevel0_.broadSearch(
ppPairs,
pointPos,
flags,
diameter,
searchBoxChanged))
{
fatalErrorInFunction<<
"Error in broadSearch for NBS (particle-particle)"<<endl;
return false;
}
if(!cellsWallLevel0_.broadSearch(
pwPairs,
NBSLevel0_.getSearchCells(),
NBSLevel0_.getCellIterator()))
{
fatalErrorInFunction<<
"Error in broadSearch for NBS (particle-wall)"<<endl;
return false;
}
return true;
}
public:
TypeInfoNV("NBS");
NBS(
const dictionary& dict,
const box& domainBox,
real minBSSize,
real maxBSSize,
const deviceViewType1D<realx3> &position,
const pFlagTypeDevice &flags,
const deviceViewType1D<real> &diam,
uint32 nWallPoints,
uint32 nWallElements,
const ViewType1D<realx3,memory_space>& wallPoints,
const ViewType1D<uint32x3,memory_space>& 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<cells> getDomainCellsLevels()const
{
return Vector<cells>("Cells", 1, NBSLevel0_.getDomainCells());
}
};
}
#endif

118
src/Interaction/contactSearch/methods/cellBased/NBS/NBSLevel0.cpp Normal file
View File

@ -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<realx3> &pointPos,
const pFlagTypeDevice &flags,
const deviceViewType1D<real> &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 "<<pairs.capacity()<<" in NBSLevel0."<<END_INFO;
}
Kokkos::fence();
}
return true;
}
pFlow::NBSLevel0::NBSLevel0
(
const box& domain,
real cellSize,
real sizeRatio,
const deviceViewType1D<realx3>& pointPos,
const pFlagTypeDevice& flags,
bool adjustableBox
)
:
mapperNBS
(
domain,
cellSize,
pointPos,
flags,
adjustableBox,
true
),
sizeRatio_(sizeRatio)
{
}
bool pFlow::NBSLevel0::broadSearch
(
csPairContainerType &pairs,
const deviceViewType1D<realx3> &pointPos,
const pFlagTypeDevice &flags,
const deviceViewType1D<real> &diameter,
bool& searchBoxChanged
)
{
if(!build(pointPos, flags, searchBoxChanged))
{
fatalErrorInFunction;
return false;
}
if(!findPairs(pairs, pointPos, flags, diameter))
{
fatalErrorInFunction;
return false;
}
return true;
}

98
src/Interaction/contactSearch/methods/cellBased/NBS/NBSLevel0.hpp Normal file
View File

@ -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<realx3>& pointPos,
const pFlagTypeDevice& flags,
const deviceViewType1D<real>& diameter
);
public:
TypeInfoNV("NBSLevel0");
INLINE_FUNCTION_HD
NBSLevel0() = default;
NBSLevel0(
const box& domain,
real cellSize,
real sizeRatio,
const deviceViewType1D<realx3>& 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<realx3>& pointPos,
const pFlagTypeDevice& flags,
const deviceViewType1D<real>& diameter,
bool& searchBoxChanged);
};
} // pFlow
#endif // __NBSLevel0_hpp__

83
src/Interaction/contactSearch/methods/cellBased/NBS/NBSLevel0Kernels.hpp Normal file
View File

@ -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<uint32>,
Kokkos::Rank<3>,
Kokkos::Schedule<Kokkos::Dynamic>,
DefaultExecutionSpace>;
template<typename T>
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<realx3>& pointPos,
const pFlagTypeDevice& flags,
const deviceViewType1D<real>& 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;
}
}

100
src/Interaction/contactSearch/methods/cellBased/NBS/NBSLoop.hpp Normal file
View File

@ -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()<nCells.x() && neighborCell.y()<nCells.y() && neighborCell.z()<nCells.z() )
{
n = cellIter.start(neighborCell.x(), neighborCell.y(), neighborCell.z());
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.next(n);
}
}
}
m = cellIter.next(m);
}

195
src/Interaction/contactSearch/methods/cellBased/NBS/cellsWallLevel0.cpp Normal file
View File

@ -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<realx3,
memory_space> &points,
const ViewType1D<uint32x3, memory_space> &vertices
)
:
cellExtent_( max(cellExtent, 0.5 ) ),
numElements_(numElements),
numPoints_(numPoints),
vertices_(vertices),
points_(points)
{
allocateArrays();
}
bool pFlow::cellsWallLevel0::resetElements
(
uint32 numElements,
uint32 numPoints,
ViewType1D<realx3, memory_space> &points,
ViewType1D<uint32x3, memory_space> &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."<<END_INFO;
Kokkos::fence();
}
}
return true;
}
pFlow::int32 pFlow::cellsWallLevel0::findPairsElementRangeCount
(
csPairContainerType &pairs,
const mapperNBS::CellIterator &particleMap
)
{
uint32 getFull =0;
//const auto pwPairs = pairs;
const auto elementBox = elementBox_;
Kokkos::parallel_reduce(
"pFlow::cellsWallLevel0::findPairsElementRangeCount",
tpPWContactSearch(numElements_, Kokkos::AUTO),
LAMBDA_HD(
const typename tpPWContactSearch::member_type & teamMember,
uint32& valueToUpdate){
const uint32 iTri = teamMember.league_rank();
const auto triBox = elementBox[iTri];
uint32 getFull2 = 0;
auto bExtent = boxExtent(triBox);
uint32 numCellBox = bExtent.x()*bExtent.y()*bExtent.z();
Kokkos::parallel_reduce(
Kokkos::TeamThreadRange( teamMember, numCellBox ),
[&] ( const uint32 linIndex, uint32 &innerUpdate )
{
int32x3 cell;
indexToCell(linIndex, triBox, cell);
uint32 n = particleMap.start(cell.x(),cell.y(),cell.z());
while( n != particleMap.NoPos)
{
// id is wall id the pair is (particle id, wall id)
if( pairs.insert(
static_cast<csIdType>(n),
static_cast<csIdType>(iTri) ) == -1 )
innerUpdate++;
n = particleMap.next(n);
}
},
getFull2
);
if ( teamMember.team_rank() == 0 ) valueToUpdate += getFull2;
},
getFull
);
return getFull;
}

139
src/Interaction/contactSearch/methods/cellBased/NBS/cellsWallLevel0.hpp Normal file
View File

@ -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<int32>;
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<uint32x3, memory_space> vertices_;
// - ref to points in the trisurface (borrowed)
ViewType1D<realx3, memory_space> points_;
// cell range of element/triangle bounding box
ViewType1D<iBoxType, memory_space> elementBox_;
using tpPWContactSearch = Kokkos::TeamPolicy<
execution_space,
Kokkos::Schedule<Kokkos::Dynamic>,
Kokkos::IndexType<uint32>>;
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<realx3,memory_space>& points,
const ViewType1D<uint32x3,memory_space>& vertices);
// - host call
// reset triangle elements if they have changed
bool resetElements(
uint32 numElements,
uint32 numPoints,
ViewType1D<realx3, memory_space>& points,
ViewType1D<uint32x3, memory_space>& 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__

194
src/Interaction/contactSearch/methods/cellBased/NBS/mapperNBS.cpp Normal file
View File

@ -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<realx3> &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 ["<<searchCells_.domainBox().minPoint()<<
" "<<searchCells_.domainBox().maxPoint()<<"]"<<END_INFO;
return true;
}
else
{
searchCells_ = cells(domainBox, cellSize);
INFORMATION<<"Search box for contact search is: ["
<< domainBox.minPoint()<<" "<<domainBox.maxPoint()<<"]"<<END_INFO;
return false;
}
}
void pFlow::mapperNBS::allocateArrays(rangeU32 nextRng)
{
checkAllocateNext(nextRng);
nullifyNext(nextRng);
reallocFill(head_, searchCells_.nx(), searchCells_.ny(), searchCells_.nz(), NoPos);
}
void pFlow::mapperNBS::checkAllocateNext(rangeU32 nextRng)
{
auto newCap = nextRng.end();
if( nextCapacity_ < newCap)
{
nextCapacity_ = newCap;
if(!nextOwner_)return;
reallocNoInit(next_, nextCapacity_);
}
}
void pFlow::mapperNBS::nullifyHead()
{
fill(head_, NoPos);
}
void pFlow::mapperNBS::nullifyNext(rangeU32 nextRng)
{
if(!nextOwner_)return;
fill(next_, nextRng, NoPos);
}
pFlow::mapperNBS::mapperNBS(
const box &domain,
real cellSize,
const deviceViewType1D<realx3> &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<realx3>& 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!"<<endl;
return false;
}
}
}
else
{
pFlow::mapperNBSKernels::buildLists(
searchCells_,
head_,
next_,
pointPos,
flags
);
searchBoxChanged = false;
}
return true;
}

212
src/Interaction/contactSearch/methods/cellBased/NBS/mapperNBS.hpp Normal file
View File

@ -0,0 +1,212 @@
/*------------------------------- 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 "phasicFlowKokkos.hpp"
#include "pointFlag.hpp"
#include "cells.hpp"
//#include "contactSearchFunctions.hpp"
namespace pFlow
{
class mapperNBS
{
public:
using HeadType = deviceViewType3D<uint32>;
using NextType = deviceViewType1D<uint32>;
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<realx3>& 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<realx3>& 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<realx3>& 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_;
}
*/

192
src/Interaction/contactSearch/methods/cellBased/NBS/mapperNBSKernels.cpp Normal file
View File

@ -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<realx3>& 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<real>(minX),
Kokkos::Min<real>(minY),
Kokkos::Min<real>(minZ),
Kokkos::Max<real>(maxX),
Kokkos::Max<real>(maxY),
Kokkos::Max<real>(maxZ)
);
minPoint = {minX, minY, minZ};
maxPoint = {maxX, maxY, maxZ};
}
bool pFlow::mapperNBSKernels::buildListsReduce
(
const cells &searchCell,
const deviceViewType3D<uint32> &head,
const deviceViewType1D<uint32> &next,
const deviceViewType1D<realx3> &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<uint32> &head,
const deviceViewType1D<uint32> &next,
const deviceViewType1D<realx3> &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;
}

54
src/Interaction/interaction/demInteraction.hpp → src/Interaction/contactSearch/methods/cellBased/NBS/mapperNBSKernels.hpp
View File

@ -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<realx3>& points,
const pFlagTypeDevice& flags,
realx3& minPoint,
realx3& maxPoint);
public:
demInteraction(systemControl& control)
:
property(),
demComponent("interaction", control)
{}
bool buildListsReduce(
const cells& searchCell,
const deviceViewType3D<uint32>& head,
const deviceViewType1D<uint32>& next,
const deviceViewType1D<realx3>& points,
const pFlagTypeDevice& flags);
demInteraction(systemControl& control, const fileSystem& file)
:
property(file),
demComponent("interaction", control)
{}
bool buildLists(
const cells& searchCell,
const deviceViewType3D<uint32>& head,
const deviceViewType1D<uint32>& next,
const deviceViewType1D<realx3>& points,
const pFlagTypeDevice& flags);
};
}
#endif //__interaction_hpp__
}

107
src/Interaction/contactSearch/cells.hpp → src/Interaction/contactSearch/methods/cellBased/cells.hpp
View File

@ -28,58 +28,40 @@ Licence:
namespace pFlow
{
template<typename indexType>
class cells
{
public:
using CellType = triple<indexType>;
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<indexType>(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<int64>(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()))
);
}
};

389
src/Interaction/contactSearch/methods/mapperNBS.hpp
View File

@ -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<typename executionSpace>
class mapperNBS
:
public cells<int32>
{
public:
using IdType = int32;
using IndexType = int32;
using Cells = cells<IndexType>;
using CellType = typename Cells::CellType;
using execution_space = executionSpace;
using memory_space = typename execution_space::memory_space;
using HeadType = ViewType3D<int32, memory_space>;
using NextType = ViewType1D<int32, memory_space>;
class cellIterator
{
private:
HeadType head_;
NextType next_;
public:
cellIterator(ViewType3D<int32, memory_space> head, ViewType1D<int32, memory_space> 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<int32, memory_space> head_;
ViewType1D<int32, memory_space> next_;
bool nextOwner_ = true;
// borrowed ownership
ViewType1D<realx3, memory_space> pointPosition_;
using rangePolicyType =
Kokkos::RangePolicy<
Kokkos::IndexType<int32>,
Kokkos::Schedule<Kokkos::Static>,
execution_space>;
INLINE_FUNCTION_H
void nullifyHead()
{
fill(
head_,
range(0,this->nx()),
range(0,this->ny()),
range(0,this->nz()),
static_cast<int32>(-1)
);
}
void nullifyNext(range nextRng)
{
if(!nextOwner_)return;
fill(
next_,
nextRng,
static_cast<int32>(-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<realx3, memory_space>& 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<realx3, memory_space>& 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<int32, memory_space>& 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<Cells>(*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<typename IncludeFunction>
INLINE_FUNCTION_H
void build(range activeRange, IncludeFunction incld)
{
checkAllocateNext(activeRange.second);
nullify(activeRange);
Cells cellIndex = static_cast<Cells>(*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<Cells>(*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<typename IncludeFunction>
INLINE_FUNCTION_H
void buildCheckInDomain(range activeRange, IncludeFunction incld)
{
checkAllocateNext(activeRange.second);
nullify(activeRange);
Cells cellIndex = static_cast<Cells>(*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__

213
src/Interaction/contactSearch/methods/multiGridNBS.hpp
View File

@ -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<typename executionSpace>
class multiGridNBS
{
public:
using NBSLevelsType = NBSLevels<executionSpace>;
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<realx3, memory_space>& position,
const ViewType1D<real, memory_space>& diam)
:
sizeRatio_(
max(
dict.getVal<real>("sizeRatio"),
1.0
)),
updateFrequency_(
max(
dict.getVal<int32>("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<Cells> cellsLvl("cells", numLevels(), numLevels(), RESERVE());
for(int32 lvl=0; lvl<numLevels(); lvl++)
{
cellsLvl[lvl] = NBSLevels_.getCells(lvl) ;
}
return cellsLvl;
}
auto getCells(int32 lvl)const
{
return NBSLevels_.getCells(lvl);
}
auto getCellIterator(int32 lvl)const
{
return NBSLevels_.getCellIterator(lvl);
}
bool objectSizeChanged(int32 newSize)
{
NBSLevels_.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<typename PairsContainer>
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<typename PairsContainer, typename IncludeFunction>
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

78
src/Interaction/contactSearch/methods/particleWallContactSearchs.cpp Normal file
View File

@ -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 <typename method>
pFlow::particleWallContactSearchs<method>::particleWallContactSearchs
(
const dictionary &dict,
const box &domain,
real minSize,
real maxSize,
const ViewType1D<realx3, memory_space> &position,
const pFlagTypeDevice &flags,
const ViewType1D<real, memory_space> &diam
)
:
domainBox_(domain),
updateInterval_
(
dict.getValOrSet<uint32>("updateInterval", 1)
)
{
}
template <typename method>
bool pFlow::particleWallContactSearchs<method>::broadSearch
(
uint32 iter,
real t,
real dt,
csPairContainerType& ppPairs,
csPairContainerType& pwPairs,
const deviceViewType1D<realx3>& pointPos,
const pFlagTypeDevice& flags,
const deviceViewType1D<real>& 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()<<endl;
return false;
}
lastUpdated_ = iter;
performedSearch_ = true;
return true;
}

123
src/Interaction/contactSearch/methods/particleWallContactSearchs.hpp Normal file
View File

@ -0,0 +1,123 @@
/*------------------------------- 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 __particleWallContactSearchs_hpp__
#define __particleWallContactSearchs_hpp__
#include "pointFlag.hpp"
#include "contactSearchGlobals.hpp"
#include "box.hpp"
namespace pFlow
{
template<typename method>
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<MethodType&>(*this);
}
inline
const auto& getMethod()const
{
return static_cast<const MethodType&>(*this);
}
public:
particleWallContactSearchs(
const dictionary& dict,
const box& domain,
real minSize,
real maxSize,
const ViewType1D<realx3, memory_space>& position,
const pFlagTypeDevice &flags,
const ViewType1D<real, memory_space>& diam
);
bool broadSearch
(
uint32 iter,
real t,
real dt,
csPairContainerType& ppPairs,
csPairContainerType& pwPairs,
const deviceViewType1D<realx3>& pointPos,
const pFlagTypeDevice& flags,
const deviceViewType1D<real>& 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__

150
src/Interaction/contactSearch/wallMappings/cellMapping.hpp
View File

@ -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<executionSpace>;
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<IndexType>;
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<Cells>& ppCells,
int32 numPoints,
int32 numElements,
const ViewType1D<realx3,memory_space>& points,
const ViewType1D<int32x3,memory_space>& vertices
)
:
updateFrequency_(
max(
dict.getValOrSet<int32>(
"updateFrequency",
1),
1)),
cellExtent_(
max(
dict.getValOrSet<real>(
"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<typename PairsContainer, typename particleMapType>
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

285
src/Interaction/contactSearch/wallMappings/cellsWallLevel0.hpp
View File

@ -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<int32>
{
public:
using IdType = int32;
using IndexType = int32;
using Cells = cells<IndexType>;
using CellType = typename Cells::CellType;
using execution_space = executionSpace;
using memory_space = typename execution_space::memory_space;
using iBoxType = iBox<IndexType>;
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<int32x3, memory_space> vertices_;
// - ref to points in the trisurface (borrowed)
ViewType1D<realx3, memory_space> points_;
// cell range of element/triangle bounding box
ViewType1D<iBoxType, memory_space> elementBox_;
using tpPWContactSearch = Kokkos::TeamPolicy<
execution_space,
Kokkos::Schedule<Kokkos::Dynamic>,
Kokkos::IndexType<int32>
>;
using rpFindCellRange2Type =
Kokkos::RangePolicy<TagFindCellRange2, execution_space, Kokkos::IndexType<int32>>;
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<realx3,memory_space>& points,
const ViewType1D<int32x3,memory_space>& 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<realx3, memory_space>& points,
ViewType1D<int32x3, memory_space>& 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<typename PairsContainer, typename particleMapType>
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<typename PairsContainer, typename particleMapType>
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."<<endINFO;
Kokkos::fence();
}
}
return true;
}
template<typename PairsContainer, typename CellIteratorType>
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<IdType>(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__

152
src/Interaction/contactSearch/wallMappings/cellsWallLevels.hpp
View File

@ -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<executionSpace>;
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<IndexType>;
protected:
int32 numLevles_=1;
Vector<cellsWallLevel0Type> cellsWallLevels_;
public:
TypeInfoNV("cellsWallLevels");
FUNCTION_H
cellsWallLevels(
int32 numLevels,
const Vector<Cells>& cellsLevels,
real cellExtent,
int32 numPoints,
int32 numElements,
const ViewType1D<realx3,memory_space>& points,
const ViewType1D<int32x3,memory_space>& vertices
)
:
numLevles_(numLevels),
cellsWallLevels_("cellsWallLevels",numLevels, numLevels, RESERVE())
{
for(int32 lvl=0; lvl<numLevles_; lvl++)
{
cellsWallLevels_[lvl] =
cellsWallLevel0Type(
cellsLevels[lvl],
cellExtent,
numPoints,
numElements,
points,
vertices);
}
}
template<typename PairsContainer, typename particleMapType>
bool broadSearch(PairsContainer& pairs, particleMapType& particleMap)
{
// map walls onto the cells
for(int32 lvl=0; lvl<numLevles_; lvl++)
{
cellsWallLevels_[lvl].build();
}
this->particleWallFindPairs(pairs, particleMap);
return true;
}
template<typename PairsContainer, typename particleMapType>
bool particleWallFindPairs(PairsContainer& pairs, particleMapType& particleMap)
{
int32 getFull = 1;
while (getFull)
{
getFull = 0;
for(int32 lvl=0; lvl<numLevles_; lvl++)
{
getFull +=
cellsWallLevels_[lvl].findPairsElementRangeCount(
pairs,
particleMap.getCellIterator(lvl));
}
if(getFull)
{
// - resize the container
// note that getFull now shows the number of failed insertions.
uint32 len = max(getFull, 5000);
auto oldCap = pairs.capacity();
pairs.increaseCapacityBy(len);
INFORMATION<<"Contact pair container capacity increased from "<<
oldCap << " to "
<< pairs.capacity() <<" in cellsWallLevels."<<endINFO;
Kokkos::fence();
}
}
return true;
}
}; // cellsWallLevels
} // pFlow
#endif // __cellsWallLevels_hpp__

153
src/Interaction/contactSearch/wallMappings/multiGridMapping.hpp
View File

@ -1,153 +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 __multiGridMapping_hpp__
#define __multiGridMapping_hpp__
#include "cellsWallLevels.hpp"
#include "dictionary.hpp"
namespace pFlow
{
template<
typename executionSpace
>
class multiGridMapping
{
public:
using CellsWallLevelType = cellsWallLevels<executionSpace>;
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<IndexType>;
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<Cells>& ppCells,
int32 numPoints,
int32 numElements,
const ViewType1D<realx3,memory_space>& points,
const ViewType1D<int32x3,memory_space>& vertices
)
:
updateFrequency_(
max(
dict.getVal<int32>("updateFrequency"),
1)),
cellExtent_(
max(
dict.getVal<real>("cellExtent"),
0.5)),
cellsWallLevle_(
numLevels,
ppCells,
cellExtent_,
numPoints,
numElements,
points,
vertices
)
{
REPORT(3)<<"Multi-grid wall mapping with "<<
yellowText(numLevels)<<" levels has been created."<<endREPORT;
}
bool enterBoadSearch()const
{
return currentIter_%updateFrequency_==0;
}
bool performedSearch()const
{
return performedSearch_;
}
template<typename PairsContainer, typename particleMapType>
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

36
src/Interaction/interaction/interaction.cpp
View File

@ -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<dictionary>(
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> pFlow::interaction::create
clType);
REPORT(1)<< "Selecting interaction model..."<<endREPORT;
REPORT(1)<< "Selecting interaction model..."<<END_REPORT;
if( systemControlvCtorSelector_.search(interactionModel) )
{
auto objPtr =
systemControlvCtorSelector_[interactionModel]
(control, prtcl, geom);
REPORT(2)<<"Model "<<greenText(interactionModel)<<" is created."<<endREPORT;
REPORT(2)<<"Model "<<Green_Text(interactionModel)<<" is created."<<END_REPORT;
return objPtr;
}
else

59
src/Interaction/interaction/interaction.hpp
View File

@ -21,9 +21,10 @@ Licence:
#ifndef __interaction_hpp__
#define __interaction_hpp__
#include "demInteraction.hpp"
#include "eventObserver.hpp"
#include "interactionBase.hpp"
#include "demComponent.hpp"
#include "property.hpp"
#include "observer.hpp"
#include "systemControl.hpp"
#include "contactSearch.hpp"
@ -31,27 +32,26 @@ Licence:
namespace pFlow
{
class particles;
class geometry;
class interaction
:
public demInteraction,
public eventObserver,
public interactionBase
{
public:
public property,
public observer,
public demComponent
using IdType = typename interactionBase::IdType;
{
private:
using IndexType = typename interactionBase::IndexType;
/// reference to particles object
const particles& particles_;
using ExecutionSpace = typename interactionBase::ExecutionSpace;
/// reference to geometry object
const geometry& geometry_;
protected:
/// interaction file dictionary
dictionary& fileDict_;
/// contact search object for pp and pw interactions
uniquePtr<contactSearch> 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<interaction> create(

201
src/Interaction/sphereInteraction/sphereInteraction.cpp
View File

@ -18,22 +18,18 @@ Licence:
-----------------------------------------------------------------------------*/
template<
typename contactForceModel,
typename geometryMotionModel,
template <class, class, class> class contactListType >
bool pFlow::sphereInteraction<contactForceModel,geometryMotionModel, contactListType>::
createSphereInteraction()
template<typename cFM,typename gMM,template <class, class, class> class cLT>
bool pFlow::sphereInteraction<cFM,gMM, cLT>::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 . . ."<<endREPORT;
REPORT(1)<<"Createing contact force model . . ."<<END_REPORT;
forceModel_ = makeUnique<ContactForceModel>(
this->numMaterials(),
rhoD.deviceVector(),
rhoD.deviceView(),
modelDict );
@ -47,17 +43,9 @@ bool pFlow::sphereInteraction<contactForceModel,geometryMotionModel, contactList
}
template<
typename contactForceModel,
typename geometryMotionModel,
template <class, class, class> class contactListType >
bool pFlow::sphereInteraction<contactForceModel,geometryMotionModel, contactListType>::
sphereSphereInteraction()
template<typename cFM,typename gMM,template <class, class, class> class cLT>
bool pFlow::sphereInteraction<cFM,gMM, cLT>::sphereSphereInteraction()
{
auto lastItem = ppContactList_().loopCount();
// create the kernel functor
@ -66,17 +54,17 @@ bool pFlow::sphereInteraction<contactForceModel,geometryMotionModel, contactList
this->dt(),
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<contactForceModel,geometryMotionModel, contactList
}
template<
typename contactForceModel,
typename geometryMotionModel,
template <class, class, class> class contactListType >
bool pFlow::sphereInteraction<contactForceModel,geometryMotionModel, contactListType>::
sphereWallInteraction()
template<typename cFM,typename gMM,template <class, class, class> class cLT>
bool pFlow::sphereInteraction<cFM,gMM, cLT>::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<contactForceModel,geometryMotionModel, contactList
Kokkos::fence();
return true;
}
}
template<typename cFM,typename gMM,template <class, class, class> class cLT>
pFlow::sphereInteraction<cFM,gMM, cLT>::sphereInteraction
(
systemControl& control,
const particles& prtcl,
const geometry& geom
)
:
interaction(control, prtcl, geom),
geometryMotion_(dynamic_cast<const GeometryMotionModel&>(geom)),
sphParticles_(dynamic_cast<const sphereParticles&>(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<typename cFM,typename gMM,template <class, class, class> class cLT>
bool pFlow::sphereInteraction<cFM,gMM, cLT>::beforeIteration()
{
return true;
}
template<typename cFM,typename gMM,template <class, class, class> class cLT>
bool pFlow::sphereInteraction<cFM,gMM, cLT>::iterate()
{
auto iter = this->currentIter();
auto t = this->currentTime();
auto dt = this->dt();
//output<<"iter, t, dt "<< iter<<" "<< t << " "<<dt<<endl;
bool broadSearch = contactSearch_().enterBroadSearch(iter, t, dt);
if(broadSearch)
{
contactListTimer0_.start();
ppContactList_().beforeBroadSearch();
pwContactList_().beforeBroadSearch();
contactListTimer0_.end();
}
if( sphParticles_.numActive()<=0)return true;
if( !contactSearch_().broadSearch(
iter,
t,
dt,
ppContactList_(),
pwContactList_()) )
{
fatalErrorInFunction<<
"unable to perform broadSearch.\n";
fatalExit;
}
if(broadSearch && contactSearch_().performedBroadSearch())
{
contactListTimer_.start();
ppContactList_().afterBroadSearch();
pwContactList_().afterBroadSearch();
contactListTimer_.end();
}
ppInteractionTimer_.start();
sphereSphereInteraction();
ppInteractionTimer_.end();
pwInteractionTimer_.start();
sphereWallInteraction();
pwInteractionTimer_.end();
return true;
}
template<typename cFM,typename gMM,template <class, class, class> class cLT>
bool pFlow::sphereInteraction<cFM,gMM, cLT>::afterIteration()
{
return true;
}
template<typename cFM,typename gMM,template <class, class, class> class cLT>
bool pFlow::sphereInteraction<cFM,gMM, cLT>::hearChanges
(
real t,
real dt,
uint32 iter,
const message& msg,
const anyList& varList
)
{
if(msg.equivalentTo(message::ITEM_REARRANGE))
{
notImplementedFunction;
}
return true;
}

167
src/Interaction/sphereInteraction/sphereInteraction.hpp
View File

@ -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<ModelStorage, ExecutionSpace, IdType>;
contactListType<ModelStorage, DefaultExecutionSpace, IdType>;
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> contactSearch_ = nullptr;
/// contact force model
uniquePtr<ContactForceModel> forceModel_ = nullptr;
uniquePtr<ContactForceModel> forceModel_ = nullptr;
/// contact list for particle-particle interactoins (keeps the history)
uniquePtr<ContactListType> ppContactList_ = nullptr;
uniquePtr<ContactListType> ppContactList_ = nullptr;
/// contact list for particle-wall interactions (keeps the history)
uniquePtr<ContactListType> pwContactList_ = nullptr;
uniquePtr<ContactListType> 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::IndexType<int32>, Kokkos::Schedule<Kokkos::Dynamic>>;
Kokkos::RangePolicy<Kokkos::IndexType<uint32>, Kokkos::Schedule<Kokkos::Dynamic>>;
/// 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<const GeometryMotionModel&>(geom)),
sphParticles_(dynamic_cast<const sphereParticles&>(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"<<endInfo;
////Info<<"interaction iterrate start"<<endInfo;
if(this->contactSearch_)
{
if( this->contactSearch_().ppEnterBroadSearch())
{
//Info<<" before ppEnterBroadSearch"<<endInfo;
ppContactList_().beforeBroadSearch();
//Info<<" after ppEnterBroadSearch"<<endInfo;
}
if(this->contactSearch_().pwEnterBroadSearch())
{
//Info<<" before pwEnterBroadSearch"<<endInfo;
pwContactList_().beforeBroadSearch();
//Info<<" after pwEnterBroadSearch"<<endInfo;
}
//Info<<" before broadSearch"<<endInfo;
if( !contactSearch_().broadSearch(
ppContactList_(),
pwContactList_()) )
{
fatalErrorInFunction<<
"unable to perform broadSearch.\n";
fatalExit;
}
//Info<<" before broadSearch"<<endInfo;
if(this->contactSearch_().ppPerformedBroadSearch())
{
//Info<<" before afterBroadSearch"<<endInfo;
ppContactList_().afterBroadSearch();
//Info<<" after afterBroadSearch"<<endInfo;
}
if(this->contactSearch_().pwPerformedBroadSearch())
{
//Info<<" before pwContactList_().afterBroadSearch()"<<endInfo;
pwContactList_().afterBroadSearch();
//Info<<" after pwContactList_().afterBroadSearch()"<<endInfo;
}
}
//Info<<"after contact search"<<endInfo;
if( sphParticles_.numActive()<=0)return true;
//Info<<"before sphereSphereInteraction "<<endInfo;
ppInteractionTimer_.start();
sphereSphereInteraction();
ppInteractionTimer_.end();
//Info<<"after sphereSphereInteraction "<<endInfo;
//Info<<"before sphereWallInteraction "<<endInfo;
pwInteractionTimer_.start();
sphereWallInteraction();
pwInteractionTimer_.end();
//Info<<"after sphereWallInteraction "<<endInfo;
return true;
}
bool afterIteration() override
{
return true;
}
bool update(const eventMessage& msg)override
{
// it requires not action regarding any changes in the
// point structure
return true;
}
bool sphereSphereInteraction();
bool sphereWallInteraction();
/// Check for changes in the point structures. (overriden from observer)
bool hearChanges(
real t,
real dt,
uint32 iter,
const message& msg,
const anyList& varList)override;
};

18
src/Interaction/sphereInteraction/sphereInteractionKernels.hpp
View File

@ -42,7 +42,7 @@ struct ppInteractionFunctor
ContactListType tobeFilled_;
deviceViewType1D<real> diam_;
deviceViewType1D<int8> propId_;
deviceViewType1D<uint32> propId_;
deviceViewType1D<realx3> pos_;
deviceViewType1D<realx3> lVel_;
deviceViewType1D<realx3> rVel_;
@ -55,7 +55,7 @@ struct ppInteractionFunctor
ContactForceModel forceModel,
ContactListType tobeFilled,
deviceViewType1D<real> diam,
deviceViewType1D<int8> propId,
deviceViewType1D<uint32> propId,
deviceViewType1D<realx3> pos,
deviceViewType1D<realx3> lVel,
deviceViewType1D<realx3> 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<real> diam_;
deviceViewType1D<int8> propId_;
deviceViewType1D<uint32> propId_;
deviceViewType1D<realx3> pos_;
deviceViewType1D<realx3> lVel_;
deviceViewType1D<realx3> rVel_;
deviceViewType1D<realx3> cForce_;
deviceViewType1D<realx3> cTorque_;
deviceViewType1D<int8> wTriMotionIndex_;
deviceViewType1D<int8> wPropId_;
deviceViewType1D<uint32> wTriMotionIndex_;
deviceViewType1D<uint32> wPropId_;
deviceViewType1D<realx3> wCForce_;
@ -199,14 +199,14 @@ struct pwInteractionFunctor
TraingleAccessor triangles,
MotionModel motionModel ,
deviceViewType1D<real> diam ,
deviceViewType1D<int8> propId,
deviceViewType1D<uint32> propId,
deviceViewType1D<realx3> pos ,
deviceViewType1D<realx3> lVel,
deviceViewType1D<realx3> rVel,
deviceViewType1D<realx3> cForce,
deviceViewType1D<realx3> cTorque ,
deviceViewType1D<int8> wTriMotionIndex,
deviceViewType1D<int8> wPropId,
deviceViewType1D<uint32> wTriMotionIndex,
deviceViewType1D<uint32> wPropId,
deviceViewType1D<realx3> wCForce)
:
dt_(dt),

35
src/Interaction/sphereInteraction/sphereInteractions.cpp
View File

@ -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>;
*/