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

Merge pull request #49 from PhasicFlow/multigridNBS 

Multigrid nbs
This commit is contained in:
PhasicFlow 2022-10-30 18:12:04 +03:30 committed by GitHub
parent 61c832a27c 39534d7a6a
commit 9fff4638ed
No known key found for this signature in database
GPG Key ID: 4AEE18F83AFDEB23
17 changed files with 1506 additions and 750 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

15
src/Interaction/contactSearch/ContactSearch/ContactSearch.H
View File

@ -31,7 +31,7 @@ namespace pFlow
template< template<
template<class> class BaseMethod, template<class> class BaseMethod,
template<class, class, class> class WallMapping template<class> class WallMapping
> >
class ContactSearch class ContactSearch
: :
@ -53,10 +53,7 @@ public:
using WallMappingType = using WallMappingType =
WallMapping< WallMapping<
ExecutionSpace, ExecutionSpace>;
IdType,
IndexType
>;
protected: protected:
@ -83,7 +80,7 @@ public:
auto method = dict().getVal<word>("method"); auto method = dict().getVal<word>("method");
auto wmMethod = dict().getVal<word>("wallMapping"); auto wmMethod = dict().getVal<word>("wallMapping");
auto nbDict = dict().subDictOrCreate(method+"Info"); auto nbDict = dict().subDict(method+"Info");
real minD, maxD; real minD, maxD;
this->Particles().boundingSphereMinMax(minD, maxD); this->Particles().boundingSphereMinMax(minD, maxD);
@ -96,6 +93,7 @@ public:
( (
nbDict, nbDict,
this->domain(), this->domain(),
minD,
maxD, maxD,
position, position,
diam diam
@ -104,7 +102,7 @@ public:
greenText(particleContactSearch_().typeName())<<endReport; greenText(particleContactSearch_().typeName())<<endReport;
auto wmDict = dict().subDictOrCreate(wmMethod+"Info"); auto wmDict = dict().subDict(wmMethod+"Info");
int32 wnPoints = this->Geometry().numPoints(); int32 wnPoints = this->Geometry().numPoints();
int32 wnTri = this->Geometry().size(); int32 wnTri = this->Geometry().size();
@ -115,7 +113,8 @@ public:
wallMapping_ = wallMapping_ =
makeUnique<WallMappingType>( makeUnique<WallMappingType>(
wmDict, wmDict,
particleContactSearch_().getCells(), particleContactSearch_().numLevels(),
particleContactSearch_().getCellsLevels(),
wnPoints, wnPoints,
wnTri, wnTri,
wPoints, wPoints,

8
src/Interaction/contactSearch/ContactSearch/ContactSearchs.C
View File

@ -20,7 +20,11 @@ Licence:
#include "ContactSearch.H" #include "ContactSearch.H"
#include "cellsSimple.H" #include "cellMapping.H"
#include "NBS.H" #include "NBS.H"
#include "multiGridNBS.H"
#include "multiGridMapping.H"
template class pFlow::ContactSearch<pFlow::NBS, pFlow::cellsSimple>;
template class pFlow::ContactSearch<pFlow::NBS, pFlow::cellMapping>;
template class pFlow::ContactSearch<pFlow::multiGridNBS, pFlow::multiGridMapping>;

22
src/Interaction/contactSearch/methods/NBS.H
View File

@ -35,6 +35,8 @@ public:
using NBSLevel0Type = NBSLevel0<executionSpace>; using NBSLevel0Type = NBSLevel0<executionSpace>;
using cellIterator = typename NBSLevel0Type::cellIterator;
using IdType = typename NBSLevel0Type::IdType; using IdType = typename NBSLevel0Type::IdType;
using IndexType = typename NBSLevel0Type::IndexType; using IndexType = typename NBSLevel0Type::IndexType;
@ -81,8 +83,9 @@ public:
TypeNameNV("NBS"); TypeNameNV("NBS");
NBS( NBS(
dictionary dict, const dictionary& dict,
const box& domain, const box& domain,
real minSize,
real maxSize, real maxSize,
const ViewType1D<realx3, memory_space>& position, const ViewType1D<realx3, memory_space>& position,
const ViewType1D<real, memory_space>& diam) const ViewType1D<real, memory_space>& diam)
@ -126,11 +129,26 @@ public:
return performedSearch_; return performedSearch_;
} }
auto getCellIterator()const Vector<cellIterator> getCellIteratorLevels()
{
return Vector<cellIterator>("cellIterator", 1, NBSLevel0_.getCellIterator());
}
auto getCellIterator(int32 lvl)const
{ {
return NBSLevel0_.getCellIterator(); return NBSLevel0_.getCellIterator();
} }
int32 numLevels()const
{
return 1;
}
Vector<Cells> getCellsLevels()const
{
return Vector<Cells>("Cells", 1, NBSLevel0_.getCells());
}
auto getCells()const auto getCells()const
{ {
return NBSLevel0_.getCells(); return NBSLevel0_.getCells();

82
src/Interaction/contactSearch/methods/NBSCrossLoop.H Normal file
View File

@ -0,0 +1,82 @@
/*------------------------------- phasicFlow ---------------------------------
O C enter of
O O E ngineering and
O O M ultiscale modeling of
OOOOOOO F luid flow
------------------------------------------------------------------------------
Copyright (C): www.cemf.ir
email: hamid.r.norouzi AT gmail.com
------------------------------------------------------------------------------
Licence:
This file is part of phasicFlow code. It is a free software for simulating
granular and multiphase flows. You can redistribute it and/or modify it under
the terms of GNU General Public License v3 or any other later versions.
phasicFlow is distributed to help others in their research in the field of
granular and multiphase flows, but WITHOUT ANY WARRANTY; without even the
implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
-----------------------------------------------------------------------------*/
int32 m = this->head_(i,j,k);
CellType currentCell(i,j,k);
int32 n = -1;
while( m > -1 )
{
auto p_m = this->pointPosition_[m];
auto d_m = this->sizeRatio_* this->diameter_[m];
int32x3 crossIndex = mapIndexLevels(
int32x3(i,j,k),
level_,
upperLevel.level());
for(int32 ii = -1; ii<2; ii++)
{
for(int32 jj=-1; jj<2; jj++)
{
int32 kk=-1;
while( kk < 2)
{
int32x3 nghbrCI = crossIndex + int32x3(ii,jj,kk);
if( upperLevel.isInRange(nghbrCI) )
{
n = upperLevel.head_(
nghbrCI.x(),
nghbrCI.y(),
nghbrCI.z());
while( n >-1)
{
auto p_n = this->pointPosition_[n];
auto d_n = this->sizeRatio_*this->diameter_[n];
if( sphereSphereCheck(p_m, p_n, d_m, d_n) )
{
auto ln = n;
auto lm = m;
if(lm>ln) this->Swap(lm,ln);
if( auto res = pairs.insert(lm,ln); res <0)
{
getFullUpdate++;
}
}
n = this->next_[n];
}
}
kk++;
}
}
}
m = this->next_[m];
}

117
src/Interaction/contactSearch/methods/NBSLevel.H
View File

@ -2,69 +2,126 @@
#define __NBSLevel_H__ #define __NBSLevel_H__
#include "cells.H" #include "NBSLevel0.H"
#include "contactSearchFunctions.H"
namespace pFlow
{
INLINE_FUNCTION_HD
int32x3 mapIndexLevels(const int32x3& ind, int32 lowerLevel, int32 upperLevel);
template<typename executionSpace> template<typename executionSpace>
class class
NBSLevel NBSLevel
:
public NBSLevel0<executionSpace>
{ {
public: public:
using NBSLevel0Type = NBSLevel0<executionSpace>;
using execution_space= executionSpace; using cellIterator = typename NBSLevel0Type::cellIterator;
using memory_space = typename execution_space::memory_space; using IdType = typename NBSLevel0Type::IdType;
using rangePolicy = Kokkos::RangePolicy< using IndexType = typename NBSLevel0Type::IndexType;
Kokkos::IndexType<int32>,
Kokkos::Schedule<Kokkos::Static>, using Cells = typename NBSLevel0Type::Cells;
execution_space>
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: protected:
int32x3 gridExtent_; int32 level_ = 0;
ViewType3D<int32, memory_space> head_;
int8 level_ = 1;
public: public:
NBSLevel(int8 lvl, int32x3 gridExtent) TypeNameNV("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)
: :
gridExtent_(gridExtent), NBSLevel0Type(
head_("NBSLevel::head", gridExtent.x(), gridExtent.y(), gridExtent.z()), domain,
cellSize,
sizeRatio,
position,
diam,
lvl==0),
level_(lvl) level_(lvl)
{} {}
INLINE_FUNCION_HD INLINE_FUNCTION_HD
auto& head(int32 i, int32 j, int32 k) NBSLevel(const NBSLevel&) = default;
{
return head_(i,j,k);
}
INLINE_FUNCION_HD INLINE_FUNCTION_HD
auto& head() NBSLevel& operator = (const NBSLevel&) = default;
{
return head_;
}
INLINE_FUNCION_HD INLINE_FUNCTION_HD
~NBSLevel() = default;
INLINE_FUNCTION_HD
auto level()const auto level()const
{ {
return level_; return level_;
} }
INLINE_FUNCION_HD template<typename PairsContainer>
const auto& gridExtent()const INLINE_FUNCTION_H
int32 findPairsCountCross(PairsContainer& pairs, NBSLevel& upperLevel)
{ {
return gridExtent_;
mdrPolicyFindPairs
mdrPolicy(
{0,0,0},
{this->nx(),this->ny(),this->nz()} );
int32 notInsertedPairs;
Kokkos::parallel_reduce (
"NBSLevel::findPairsCountCross",
mdrPolicy,
CLASS_LAMBDA_HD(int32 i, int32 j, int32 k, int32& getFullUpdate){
#include "NBSCrossLoop.H"
}, notInsertedPairs);
return notInsertedPairs;
} }
}; };
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 #endif

65
src/Interaction/contactSearch/methods/NBSLevel0.H
View File

@ -37,6 +37,8 @@ public:
using MapperType = mapperNBS<executionSpace>; using MapperType = mapperNBS<executionSpace>;
using cellIterator = typename MapperType::cellIterator;
using IdType = typename MapperType::IdType; using IdType = typename MapperType::IdType;
using IndexType = typename MapperType::IndexType; using IndexType = typename MapperType::IndexType;
@ -49,6 +51,10 @@ public:
using memory_space = typename MapperType::memory_space; using memory_space = typename MapperType::memory_space;
using HeadType = typename MapperType::HeadType;
using NextType = typename MapperType::NextType;
struct TagFindPairs{}; struct TagFindPairs{};
@ -66,9 +72,6 @@ protected:
Kokkos::Schedule<Kokkos::Dynamic>, Kokkos::Schedule<Kokkos::Dynamic>,
execution_space>; execution_space>;
private:
static INLINE_FUNCTION_HD static INLINE_FUNCTION_HD
void Swap(int32& x, int32& y) void Swap(int32& x, int32& y)
{ {
@ -81,6 +84,9 @@ public:
TypeNameNV("NBSLevel0"); TypeNameNV("NBSLevel0");
INLINE_FUNCTION_HD
NBSLevel0(){}
NBSLevel0( NBSLevel0(
const box& domain, const box& domain,
real cellSize, real cellSize,
@ -108,9 +114,10 @@ public:
real cellSize, real cellSize,
real sizeRatio, real sizeRatio,
const ViewType1D<realx3, memory_space>& position, const ViewType1D<realx3, memory_space>& position,
const ViewType1D<real, memory_space>& diam) const ViewType1D<real, memory_space>& diam,
bool nextOwner = true)
: :
MapperType(domain, cellSize, position), MapperType(domain, cellSize, position, nextOwner),
sizeRatio_(sizeRatio), sizeRatio_(sizeRatio),
diameter_(diam) diameter_(diam)
{} {}
@ -124,6 +131,19 @@ public:
INLINE_FUNCTION_HD INLINE_FUNCTION_HD
~NBSLevel0()=default; ~NBSLevel0()=default;
INLINE_FUNCTION_HD
auto sizeRatio()const
{
return sizeRatio_;
}
INLINE_FUNCTION_HD
auto& diameter()
{
return diameter_;
}
// - Perform the broad search to find pairs // - Perform the broad search to find pairs
// with force = true, perform broad search regardless of // with force = true, perform broad search regardless of
// updateFrequency_ value // updateFrequency_ value
@ -159,10 +179,7 @@ public:
INLINE_FUNCTION_H INLINE_FUNCTION_H
bool findPairs(PairsContainer& pairs) bool findPairs(PairsContainer& pairs)
{ {
mdrPolicyFindPairs
mdrPolicy(
{0,0,0},
{this->nx(),this->ny(),this->nz()} );
int32 getFull = 1; int32 getFull = 1;
@ -171,15 +188,7 @@ public:
while (getFull > 0) while (getFull > 0)
{ {
getFull = 0; getFull = findPairsCount(pairs);
Kokkos::parallel_reduce (
"NBSLevel0::findPairs",
mdrPolicy,
CLASS_LAMBDA_HD(int32 i, int32 j, int32 k, int32& getFullUpdate){
#include "NBSLoop.H"
}, getFull);
if(getFull) if(getFull)
{ {
@ -202,7 +211,27 @@ public:
return true; 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.H"
}, notInsertedPairs);
return notInsertedPairs;
}
}; };

436
src/Interaction/contactSearch/methods/NBSLevels.H Normal file
View File

@ -0,0 +1,436 @@
#ifndef __NBSLevels_H__
#define __NBSLevels_H__
#include "NBSLevel.H"
#include "NBSLevel0.H"
#include "KokkosTypes.H"
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);
Info<< "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 maxSizes = maxSizeLevels_;
auto particleLevel = particleLevel_;
int8 maxLvl = sizeRangeLevels_.size();
rangePolicyType rPolicy(first, last);
Kokkos::parallel_for(
"NBSLevels::findParticleLevel",
rPolicy,
LAMBDA_HD(int32 i)
{
for(int8 lvl = 0; lvl<maxLvl; lvl++)
{
if( diameter[i]<= maxSizes[lvl] )
{
particleLevel[i] = lvl;
return;
}
}
particleLevel[i] = static_cast<int8>(-1);
});
Kokkos::fence();
return true;
}
}; //NBSLevels
}
#endif

523
src/Interaction/contactSearch/methods/NBSold.H
View File

@ -1,523 +0,0 @@
/*------------------------------- phasicFlow ---------------------------------
O C enter of
O O E ngineering and
O O M ultiscale modeling of
OOOOOOO F luid flow
------------------------------------------------------------------------------
Copyright (C): www.cemf.ir
email: hamid.r.norouzi AT gmail.com
------------------------------------------------------------------------------
Licence:
This file is part of phasicFlow code. It is a free software for simulating
granular and multiphase flows. You can redistribute it and/or modify it under
the terms of GNU General Public License v3 or any other later versions.
phasicFlow is distributed to help others in their research in the field of
granular and multiphase flows, but WITHOUT ANY WARRANTY; without even the
implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
-----------------------------------------------------------------------------*/
#ifndef __NBS_H__
#define __NBS_H__
#include "cells.H"
#include "contactSearchFunctions.H"
#include "baseAlgorithms.H"
namespace pFlow
{
template<
typename executionSpace,
typename idType,
typename indexType=int32
>
class NBS
:
public cells<indexType>
{
public:
using IdType = idType;
using IndexType = indexType;
using Cells = cells<IndexType>;
using CellType = typename Cells::CellType;
using ExecutionSpace= executionSpace;
using memory_space = typename ExecutionSpace::memory_space;
struct TagBuild{};
struct TagFindPairs{};
class cellIterator
{
private:
ViewType3D<int32, memory_space> head_;
ViewType1D<int32, memory_space> 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;
real sizeRatio_ = 1.0;
int32 updateFrequency_= 1;
int32 currentIter_ = 0;
bool performedSearch_ = false;
ViewType1D<realx3, memory_space> pointPosition_;
ViewType1D<real, memory_space> diameter_;
ViewType3D<int32, memory_space> head_;
ViewType1D<int32, memory_space> next_;
INLINE_FUNCTION_H
void nullify()
{
fill(
head_,
range(0,this->nx()),
range(0,this->ny()),
range(0,this->nz()),
static_cast<int32>(-1)
);
fill(
next_,
range(0,capacity_),
static_cast<int32>(-1)
);
}
void nullify(range nextRng)
{
fill(
head_,
range(0,this->nx()),
range(0,this->ny()),
range(0,this->nz()),
static_cast<int32>(-1)
);
fill(
next_,
nextRng,
static_cast<int32>(-1)
);
}
using mdrPolicyFindPairs =
Kokkos::MDRangePolicy<
Kokkos::Rank<3>,
Kokkos::Schedule<Kokkos::Dynamic>,
ExecutionSpace>;
private:
void checkAllocateNext(int newCap)
{
if( capacity_ < newCap)
{
capacity_ = newCap;
Kokkos::realloc(next_, capacity_);
}
}
void allocateHead()
{
Kokkos::realloc(head_, this->nx(), this->ny(), this->nz());
}
bool performSearch()
{
if(currentIter_ % updateFrequency_ == 0)
{
currentIter_++;
return true;
}else
{
currentIter_++;
return false;
}
}
static INLINE_FUNCTION_HD
void Swap(int32& x, int32& y)
{
int32 tmp = x;
x = y;
y = tmp;
}
public:
TypeNameNV("NBS");
NBS(
const box& domain,
real cellSize,
const ViewType1D<realx3, memory_space>& position,
const ViewType1D<real, memory_space>& diam,
int32 initialContainerSize = 1)
:
Cells(domain, cellSize),
pointPosition_(position),
diameter_(diam),
head_("NBS::head_",1,1,1), //, this->nx(), this->ny(), this->nz()),
next_("NBS::next_",1) //,position.size()),
{
checkAllocateNext(pointPosition_.size());
allocateHead();
}
NBS(
const box& domain,
int32 nx,
int32 ny,
int32 nz,
const ViewType1D<realx3, memory_space>& position,
const ViewType1D<real, memory_space>& diam,
int32 initialContainerSize = 1)
:
Cells(domain, nx, ny, nz),
pointPosition_(position),
diameter_(diam),
head_("NBS::head_",nx,ny,nz), //, this->nx(), this->ny(), this->nz()),
next_("NBS::next_",1) //,position.size()),
{
checkAllocateNext(pointPosition_.size());
}
NBS(
dictionary dict,
const box& domain,
real cellSize,
const ViewType1D<realx3, memory_space>& position,
const ViewType1D<real, memory_space>& diam,
int32 initialContainerSize = 1
)
:
Cells(domain, cellSize),
pointPosition_(position),
diameter_(diam),
head_("NBS::head_",1,1,1), //, this->nx(), this->ny(), this->nz()),
next_("NBS::next_",1) //,position.size()),
{
updateFrequency_ = max(
dict.getVal<int32>("updateFrequency"),1 );
sizeRatio_ = max(dict.getVal<real>(
"sizeRatio"),1.0);
this->setCellSize(cellSize*sizeRatio_);
checkAllocateNext(pointPosition_.size());
allocateHead();
}
INLINE_FUNCTION_HD
NBS(const NBS&) = default;
INLINE_FUNCTION_HD
NBS& operator = (const NBS&) = default;
INLINE_FUNCTION_HD
~NBS()=default;
//// - Methods
bool enterBoadSearch()const
{
return currentIter_%updateFrequency_==0;
}
bool performedSearch()const
{
return performedSearch_;
}
const auto& Head()const
{
return head_;
}
const auto& Next()const
{
return next_;
}
cellIterator getCellIterator()const
{
return cellIterator(head_, next_);
}
// - Perform the broad search to find pairs
// with force = true, perform broad search regardless of
// updateFrequency_ value
// on all the points in the range of [0,numPoints_)
template<typename PairsContainer>
bool broadSearch(PairsContainer& pairs, range activeRange, bool force=false)
{
if(force) currentIter_ = 0;
performedSearch_ = false;
if( !performSearch() ) return true;
//Info<<"NBS::broadSearch(PairsContainer& pairs, range activeRange, bool force=false) before build"<<endInfo;
build(activeRange);
//Info<<"NBS::broadSearch(PairsContainer& pairs, range activeRange, bool force=false) after build"<<endInfo;
//Info<<"NBS::broadSearch(PairsContainer& pairs, range activeRange, bool force=false) before findPairs"<<endInfo;
findPairs(pairs);
//Info<<"NBS::broadSearch(PairsContainer& pairs, range activeRange, bool force=false) after findPairs"<<endInfo;
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;
build(activeRange, incld);
findPairs(pairs);
performedSearch_ = true;
return true;
}
// - build based on all points in range [0, numPoints_)
INLINE_FUNCTION_H
void build(range activeRange)
{
checkAllocateNext(activeRange.second);
nullify(activeRange);
Cells cellIndex = static_cast<Cells>(*this);
auto points = pointPosition_;
auto next = next_;
auto head = head_;
Kokkos::RangePolicy<
Kokkos::IndexType<int32>,
Kokkos::Schedule<Kokkos::Static>,
ExecutionSpace> rPolicy(activeRange.first, activeRange.second);
Kokkos::parallel_for(
"NBS::build",
rPolicy,
LAMBDA_HD(int32 i){
CellType ind = cellIndex.pointIndex(points[i]);
int32 old = Kokkos::atomic_exchange(&head(ind.x(), ind.y(), ind.z()), i);
next[i] = old;
});
Kokkos::fence();
}
template<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_;
Kokkos::RangePolicy<
Kokkos::IndexType<int32>,
Kokkos::Schedule<Kokkos::Dynamic>,
ExecutionSpace> rPolicy(activeRange.first, activeRange.second);
Kokkos::parallel_for(
"NBS::build",
rPolicy,
LAMBDA_HD(int32 i){
if( incld(i) )
{
CellType ind = cellIndex.pointIndex(points[i]);
auto old = Kokkos::atomic_exchange(&head(ind.x(), ind.y(), ind.z()), i);
next[i] = old;
}
});
Kokkos::fence();
}
// - build based on all points in range [0, numPoints_)
INLINE_FUNCTION_H
void buildCheckInDomain(range activeRange)
{
checkAllocateNext(activeRange.second);
nullify(activeRange);
Cells cellIndex = static_cast<Cells>(*this);
auto points = pointPosition_;
auto next = next_;
auto head = head_;
Kokkos::RangePolicy<
Kokkos::IndexType<int32>,
Kokkos::Schedule<Kokkos::Static>,
ExecutionSpace> rPolicy(activeRange.first, activeRange.second);
Kokkos::parallel_for(
"NBS::buildCheckInDomain",
rPolicy,
LAMBDA_HD(int32 i){
CellType ind;
if( cellIndex.pointIndexInDomain(points[i], ind) )
{
int32 old = Kokkos::atomic_exchange(&head(ind.x(), ind.y(), ind.z()), i);
next[i] = old;
}
}
);
Kokkos::fence();
}
template<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_;
Kokkos::RangePolicy<
Kokkos::IndexType<int32>,
Kokkos::Schedule<Kokkos::Dynamic>,
ExecutionSpace> rPolicy(activeRange.first, activeRange.second);
Kokkos::parallel_for(
"NBS::buildCheckInDomain",
rPolicy,
LAMBDA_HD(int32 i){
CellType ind;
if( incld(i) && cellIndex.pointIndexInDomain(points[i], ind) )
{
auto old = Kokkos::atomic_exchange(&head(ind.x(), ind.y(), ind.z()), i);
next[i] = old;
}
});
Kokkos::fence();
}
template<typename PairsContainer>
INLINE_FUNCTION_H
bool findPairs(PairsContainer& pairs)
{
mdrPolicyFindPairs
mdrPolicy(
{0,0,0},
{this->nx(),this->ny(),this->nz()} );
int32 getFull = 1;
// loop until the container size fits the numebr of contact pairs
while (getFull > 0)
{
getFull = 0;
Kokkos::parallel_reduce (
"NBS::broadSearch",
mdrPolicy,
CLASS_LAMBDA_HD(int32 i, int32 j, int32 k, int32& getFullUpdate){
#include "NBSLoop.H"
}, getFull);
if(getFull)
{
// - resize the container
// note that getFull now shows the number of failed insertions.
uint32 len = max(getFull,50) ;
auto oldCap = pairs.capacity();
pairs.increaseCapacityBy(len);
Info<< "The contact pair container capacity increased from "<<
oldCap << " to "<<pairs.capacity()<<" in NBS."<<endInfo;
}
Kokkos::fence();
}
return true;
}
bool objectSizeChanged(int32 newSize)
{
checkAllocateNext(newSize);
return true;
}
};
}
#endif

149
src/Interaction/contactSearch/methods/mapperNBS.H
View File

@ -48,12 +48,16 @@ public:
using memory_space = typename execution_space::memory_space; using memory_space = typename execution_space::memory_space;
using HeadType = ViewType3D<int32, memory_space>;
using NextType = ViewType1D<int32, memory_space>;
class cellIterator class cellIterator
{ {
private: private:
ViewType3D<int32, memory_space> head_; HeadType head_;
ViewType1D<int32, memory_space> next_; NextType next_;
public: public:
@ -85,58 +89,60 @@ protected:
ViewType1D<int32, memory_space> next_; ViewType1D<int32, memory_space> next_;
bool nextOwner_ = true;
// borrowed ownership // borrowed ownership
ViewType1D<realx3, memory_space> pointPosition_; ViewType1D<realx3, memory_space> pointPosition_;
INLINE_FUNCTION_H
void nullify()
{
fill(
head_,
range(0,this->nx()),
range(0,this->ny()),
range(0,this->nz()),
static_cast<int32>(-1)
);
fill(
next_,
range(0,capacity_),
static_cast<int32>(-1)
);
}
void nullify(range nextRng)
{
fill(
head_,
range(0,this->nx()),
range(0,this->ny()),
range(0,this->nz()),
static_cast<int32>(-1)
);
fill(
next_,
nextRng,
static_cast<int32>(-1)
);
}
using rangePolicyType = using rangePolicyType =
Kokkos::RangePolicy< Kokkos::RangePolicy<
Kokkos::IndexType<int32>, Kokkos::IndexType<int32>,
Kokkos::Schedule<Kokkos::Static>, Kokkos::Schedule<Kokkos::Static>,
execution_space>; execution_space>;
private: 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) void checkAllocateNext(int newCap)
{ {
if( capacity_ < newCap) if( capacity_ < newCap)
{ {
capacity_ = newCap; capacity_ = newCap;
if(!nextOwner_)return;
reallocNoInit(next_, capacity_); reallocNoInit(next_, capacity_);
} }
} }
@ -146,14 +152,20 @@ private:
reallocNoInit(head_, this->nx(), this->ny(), this->nz()); reallocNoInit(head_, this->nx(), this->ny(), this->nz());
} }
public: public:
TypeNameNV("mapperNBS"); TypeNameNV("mapperNBS");
INLINE_FUNCTION_HD
mapperNBS(){}
mapperNBS( mapperNBS(
const box& domain, const box& domain,
real cellSize, real cellSize,
const ViewType1D<realx3, memory_space>& position) const ViewType1D<realx3, memory_space>& position,
bool nextOwner = true)
: :
Cells(domain, cellSize), Cells(domain, cellSize),
pointPosition_(position), pointPosition_(position),
@ -163,7 +175,8 @@ public:
this->ny(), this->ny(),
this->nz() this->nz()
), ),
next_("mapperNBS::next_",1) //,position.size()), next_("mapperNBS::next_",1), //,position.size()),
nextOwner_(nextOwner)
{ {
checkAllocateNext(pointPosition_.size()); checkAllocateNext(pointPosition_.size());
} }
@ -173,12 +186,14 @@ public:
int32 nx, int32 nx,
int32 ny, int32 ny,
int32 nz, int32 nz,
const ViewType1D<realx3, memory_space>& position) const ViewType1D<realx3, memory_space>& position,
bool nextOwner = true)
: :
Cells(domain, nx, ny, nz), Cells(domain, nx, ny, nz),
pointPosition_(position), pointPosition_(position),
head_("mapperNBS::head_",nx,ny,nz), //, this->nx(), this->ny(), this->nz()), head_("mapperNBS::head_",nx,ny,nz),
next_("mapperNBS::next_",1) //,position.size()), next_("mapperNBS::next_",1),
nextOwner_(nextOwner)
{ {
checkAllocateNext(pointPosition_.size()); checkAllocateNext(pointPosition_.size());
} }
@ -205,12 +220,54 @@ public:
return cellIterator(head_, next_); return cellIterator(head_, next_);
} }
bool particlesSizeChanged(int32 newSize) bool particlesCapcityChanged(int32 newCap)
{ {
checkAllocateNext(newSize); checkAllocateNext(newCap);
return true; 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 // - build based on all points in active range
INLINE_FUNCTION_H INLINE_FUNCTION_H
void build(range activeRange) void build(range activeRange)

213
src/Interaction/contactSearch/methods/multiGridNBS.H Normal file
View File

@ -0,0 +1,213 @@
/*------------------------------- phasicFlow ---------------------------------
O C enter of
O O E ngineering and
O O M ultiscale modeling of
OOOOOOO F luid flow
------------------------------------------------------------------------------
Copyright (C): www.cemf.ir
email: hamid.r.norouzi AT gmail.com
------------------------------------------------------------------------------
Licence:
This file is part of phasicFlow code. It is a free software for simulating
granular and multiphase flows. You can redistribute it and/or modify it under
the terms of GNU General Public License v3 or any other later versions.
phasicFlow is distributed to help others in their research in the field of
granular and multiphase flows, but WITHOUT ANY WARRANTY; without even the
implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
-----------------------------------------------------------------------------*/
#ifndef __multiGridNBS_H__
#define __multiGridNBS_H__
#include "NBSLevels.H"
namespace pFlow
{
template<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:
TypeNameNV("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

150
src/Interaction/contactSearch/wallMappings/cellMapping.H Normal file
View File

@ -0,0 +1,150 @@
/*------------------------------- phasicFlow ---------------------------------
O C enter of
O O E ngineering and
O O M ultiscale modeling of
OOOOOOO F luid flow
------------------------------------------------------------------------------
Copyright (C): www.cemf.ir
email: hamid.r.norouzi AT gmail.com
------------------------------------------------------------------------------
Licence:
This file is part of phasicFlow code. It is a free software for simulating
granular and multiphase flows. You can redistribute it and/or modify it under
the terms of GNU General Public License v3 or any other later versions.
phasicFlow is distributed to help others in their research in the field of
granular and multiphase flows, but WITHOUT ANY WARRANTY; without even the
implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
-----------------------------------------------------------------------------*/
#ifndef __cellMapping_H__
#define __cellMapping_H__
#include "cellsWallLevel0.H"
#include "dictionary.H"
namespace pFlow
{
template<
typename executionSpace
>
class cellMapping
{
public:
using cellsWallLevel0Type = cellsWallLevel0<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:
TypeNameNV("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

140
src/Interaction/contactSearch/wallMappings/cellsSimple.H → src/Interaction/contactSearch/wallMappings/cellsWallLevel0.H
View File

@ -18,57 +18,48 @@ Licence:
-----------------------------------------------------------------------------*/ -----------------------------------------------------------------------------*/
#ifndef __cellsSimple_H__ #ifndef __cellsWallLevel0_H__
#define __cellsSimple_H__ #define __cellsWallLevel0_H__
#include "types.H" #include "types.H"
#include "KokkosTypes.H" #include "KokkosTypes.H"
#include "cells.H" #include "cells.H"
#include "iBox.H" #include "iBox.H"
#include "dictionary.H"
namespace pFlow namespace pFlow
{ {
template< template<
typename executionSpace, typename executionSpace
typename idType,
typename indexType = int32
> >
class cellsSimple class cellsWallLevel0
: :
public cells<indexType> public cells<int32>
{ {
public: public:
using IdType = idType; using IdType = int32;
using IndexType = indexType; using IndexType = int32;
using Cells = cells<IndexType>; using Cells = cells<IndexType>;
using CellType = typename Cells::CellType; using CellType = typename Cells::CellType;
using ExecutionSpace = executionSpace; using execution_space = executionSpace;
using memory_space = typename ExecutionSpace::memory_space; using memory_space = typename execution_space::memory_space;
using iBoxType = iBox<IndexType>; using iBoxType = iBox<IndexType>;
bool constexpr static LOOP_ELEMENT_RANGE = true;
class TagFindCellRange2{}; class TagFindCellRange2{};
protected: protected:
// - box extent // - box extent
real cellExtent_ = 0.6; real cellExtent_ = 0.5;
// - update frequency
int32 updateFrequency_=1;
int32 currentIter_ = 0;
// - number of triangle elements // - number of triangle elements
int32 numElements_ = 0; int32 numElements_ = 0;
@ -76,13 +67,10 @@ protected:
// - number of points // - number of points
int32 numPoints_ = 0; int32 numPoints_ = 0;
/// a broad search has been occured during last pass? // - ref to vectices (borrowed)
bool performedSearch_ = false;
// - ref to vectices
ViewType1D<int32x3, memory_space> vertices_; ViewType1D<int32x3, memory_space> vertices_;
// - ref to points in the trisurface // - ref to points in the trisurface (borrowed)
ViewType1D<realx3, memory_space> points_; ViewType1D<realx3, memory_space> points_;
// cell range of element/triangle bounding box // cell range of element/triangle bounding box
@ -90,44 +78,31 @@ protected:
using tpPWContactSearch = Kokkos::TeamPolicy< using tpPWContactSearch = Kokkos::TeamPolicy<
ExecutionSpace, execution_space,
Kokkos::Schedule<Kokkos::Dynamic>, Kokkos::Schedule<Kokkos::Dynamic>,
Kokkos::IndexType<int32> Kokkos::IndexType<int32>
>; >;
using rpFindCellRange2Type = using rpFindCellRange2Type =
Kokkos::RangePolicy<TagFindCellRange2, ExecutionSpace, Kokkos::IndexType<int32>>; Kokkos::RangePolicy<TagFindCellRange2, execution_space, Kokkos::IndexType<int32>>;
FUNCTION_H FUNCTION_H
void allocateArrays() void allocateArrays()
{ {
Kokkos::realloc( elementBox_, numElements_); reallocNoInit( elementBox_, numElements_);
}
private:
bool performSearch()
{
if(currentIter_ % updateFrequency_ == 0)
{
currentIter_++;
return true;
}else
{
currentIter_++;
return false;
}
} }
public: public:
TypeNameNV("cellsSimple"); TypeNameNV("cellsWallLevel0");
INLINE_FUNCTION_HD
cellsWallLevel0(){}
FUNCTION_H FUNCTION_H
cellsSimple( cellsWallLevel0(
Cells ppCells, const Cells& ppCells,
real cellExtent, real cellExtent,
int32 numPoints, int32 numPoints,
int32 numElements, int32 numElements,
@ -136,7 +111,7 @@ public:
) )
: :
Cells(ppCells), Cells(ppCells),
cellExtent_( max(cellExtent, 0.6 ) ), cellExtent_( max(cellExtent, 0.5 ) ),
numElements_(numElements), numElements_(numElements),
numPoints_(numPoints), numPoints_(numPoints),
vertices_(vertices), vertices_(vertices),
@ -146,43 +121,9 @@ public:
allocateArrays(); allocateArrays();
} }
cellsSimple(
dictionary dict,
Cells ppCells,
int32 numPoints,
int32 numElements,
const ViewType1D<realx3,memory_space>& points,
const ViewType1D<int32x3,memory_space>& vertices
)
:
Cells(ppCells),
numElements_(numElements),
numPoints_(numPoints),
vertices_(vertices),
points_(points)
{
updateFrequency_ = dict.getVal<int32>(
"updateFrequency" );
updateFrequency_ = max(updateFrequency_,1);
cellExtent_ = dict.getVal<real>(
"cellExtent");
cellExtent_ = max(cellExtent_,0.6);
allocateArrays();
}
constexpr bool loopElementRange()const
{
return LOOP_ELEMENT_RANGE;
}
// - host call // - host call
// reset triangle elements if they have changed // reset triangle elements if they have changed
FUNCTION_H
bool resetElements( bool resetElements(
int32 numElements, int32 numElements,
int32 numPoints, int32 numPoints,
@ -212,29 +153,17 @@ public:
return numElements_; return numElements_;
} }
bool enterBoadSearch()const
{
return currentIter_%updateFrequency_==0;
}
bool performedSearch()const
{
return performedSearch_;
}
template<typename PairsContainer, typename particleMapType> template<typename PairsContainer, typename particleMapType>
bool broadSearch(PairsContainer& pairs, particleMapType& particleMap, bool force=false) bool broadSearch(PairsContainer& pairs, particleMapType& particleMap)
{ {
if(force) currentIter_ = 0;
performedSearch_= false;
if(!performSearch())return true;
// map walls onto the cells // map walls onto the cells
this->build(); this->build();
this->particleWallFindPairs(pairs, particleMap); this->particleWallFindPairs(pairs, particleMap);
performedSearch_ = true;
return true; return true;
} }
@ -257,7 +186,7 @@ public:
while (getFull) while (getFull)
{ {
getFull = findPairsElementRange(pairs, particleMap); getFull = findPairsElementRangeCount(pairs, particleMap.getCellIterator(0));
if(getFull) if(getFull)
{ {
@ -269,7 +198,7 @@ public:
Info<<"Contact pair container capacity increased from "<< Info<<"Contact pair container capacity increased from "<<
oldCap << " to " oldCap << " to "
<< pairs.capacity() <<" in cellsSimple."<<endInfo; << pairs.capacity() <<" in cellsWallLevel0."<<endInfo;
Kokkos::fence(); Kokkos::fence();
} }
@ -279,21 +208,20 @@ public:
} }
template<typename PairsContainer, typename particleMapType> template<typename PairsContainer, typename CellIteratorType>
int32 findPairsElementRange(PairsContainer& pairs, particleMapType& particleMap) int32 findPairsElementRangeCount(PairsContainer& pairs, CellIteratorType cellIter)
{ {
int32 getFull =0; int32 getFull =0;
const auto pwPairs = pairs; const auto pwPairs = pairs;
const auto elementBox = elementBox_; const auto elementBox = elementBox_;
auto cellIter = particleMap.getCellIterator();
Kokkos::parallel_reduce( Kokkos::parallel_reduce(
"cellsSimple::findPairsElementRangeModified2", "cellsSimple::findPairsElementRangeModified2",
tpPWContactSearch(numElements_, Kokkos::AUTO), tpPWContactSearch(numElements_, Kokkos::AUTO),
LAMBDA_HD( LAMBDA_HD(
const typename tpPWContactSearch::member_type & teamMember, int32& valueToUpdate){ const typename tpPWContactSearch::member_type & teamMember,
int32& valueToUpdate){
const int32 iTri = teamMember.league_rank(); const int32 iTri = teamMember.league_rank();
@ -349,9 +277,9 @@ public:
} }
}; }; // cellsWallLevel0
} } // pFlow
#endif #endif // __cellsWallLevel0_H__

152
src/Interaction/contactSearch/wallMappings/cellsWallLevels.H Normal file
View File

@ -0,0 +1,152 @@
/*------------------------------- phasicFlow ---------------------------------
O C enter of
O O E ngineering and
O O M ultiscale modeling of
OOOOOOO F luid flow
------------------------------------------------------------------------------
Copyright (C): www.cemf.ir
email: hamid.r.norouzi AT gmail.com
------------------------------------------------------------------------------
Licence:
This file is part of phasicFlow code. It is a free software for simulating
granular and multiphase flows. You can redistribute it and/or modify it under
the terms of GNU General Public License v3 or any other later versions.
phasicFlow is distributed to help others in their research in the field of
granular and multiphase flows, but WITHOUT ANY WARRANTY; without even the
implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
-----------------------------------------------------------------------------*/
#ifndef __cellsWallLevels_H__
#define __cellsWallLevels_H__
#include "cellsWallLevel0.H"
namespace pFlow
{
template<
typename executionSpace
>
class cellsWallLevels
{
public:
using cellsWallLevel0Type = cellsWallLevel0<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:
TypeNameNV("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);
Info<<"Contact pair container capacity increased from "<<
oldCap << " to "
<< pairs.capacity() <<" in cellsWallLevels."<<endInfo;
Kokkos::fence();
}
}
return true;
}
}; // cellsWallLevels
} // pFlow
#endif // __cellsWallLevels_H__

153
src/Interaction/contactSearch/wallMappings/multiGridMapping.H Normal file
View File

@ -0,0 +1,153 @@
/*------------------------------- 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_H__
#define __multiGridMapping_H__
#include "cellsWallLevels.H"
#include "dictionary.H"
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:
TypeNameNV("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

3
src/phasicFlow/Kokkos/KokkosTypes.H
View File

@ -61,6 +61,9 @@ template<typename T, typename... properties>
template<typename T, typename... properties> template<typename T, typename... properties>
using ViewType1D = Kokkos::View<T*,properties...>; using ViewType1D = Kokkos::View<T*,properties...>;
template<typename T, typename... properties>
using DualViewType1D = Kokkos::View<T*,properties...>;
template<typename T, typename... properties> template<typename T, typename... properties>
using ViewType3D = Kokkos::View<T***,properties...>; using ViewType3D = Kokkos::View<T***,properties...>;

10
src/phasicFlow/structuredData/box/box.C
View File

@ -21,16 +21,6 @@ Licence:
#include "box.H" #include "box.H"
FUNCTION_HD
pFlow::box::box
(
const realx3& minP,
const realx3& maxP
)
:
min_(minP),
max_(maxP)
{}
FUNCTION_H FUNCTION_H
pFlow::box::box pFlow::box::box

16
src/phasicFlow/structuredData/box/box.H
View File

@ -34,10 +34,10 @@ class box
protected: protected:
// - min point // - min point
realx3 min_; realx3 min_{0,0,0};
// - max point // - max point
realx3 max_; realx3 max_{1,1,1};
public: public:
@ -45,8 +45,16 @@ public:
TypeNameNV("box"); TypeNameNV("box");
//// - Constructors //// - Constructors
FUNCTION_HD INLINE_FUNCTION_HD
box(const realx3& minP, const realx3& maxP); box(){}
INLINE_FUNCTION_HD
box(const realx3& minP, const realx3& maxP)
:
min_(minP),
max_(maxP)
{}
FUNCTION_H FUNCTION_H
box(const dictionary& dict); box(const dictionary& dict);