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

multigrid NBS is added, not tested with particle insertion

This commit is contained in:
hamidrezanorouzi 2022-10-30 18:08:32 +03:30
parent c47d5f5cfa
commit 39534d7a6a
17 changed files with 1344 additions and 896 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<class> class BaseMethod,
template<class, class, class> class WallMapping
template<class> class WallMapping
>
class ContactSearch
:
@ -53,10 +53,7 @@ public:
using WallMappingType =
WallMapping<
ExecutionSpace,
IdType,
IndexType
>;
ExecutionSpace>;
protected:
@ -83,7 +80,7 @@ public:
auto method = dict().getVal<word>("method");
auto wmMethod = dict().getVal<word>("wallMapping");
auto nbDict = dict().subDictOrCreate(method+"Info");
auto nbDict = dict().subDict(method+"Info");
real minD, maxD;
this->Particles().boundingSphereMinMax(minD, maxD);
@ -96,6 +93,7 @@ public:
(
nbDict,
this->domain(),
minD,
maxD,
position,
diam
@ -104,7 +102,7 @@ public:
greenText(particleContactSearch_().typeName())<<endReport;
auto wmDict = dict().subDictOrCreate(wmMethod+"Info");
auto wmDict = dict().subDict(wmMethod+"Info");
int32 wnPoints = this->Geometry().numPoints();
int32 wnTri = this->Geometry().size();
@ -115,7 +113,8 @@ public:
wallMapping_ =
makeUnique<WallMappingType>(
wmDict,
particleContactSearch_().getCells(),
particleContactSearch_().numLevels(),
particleContactSearch_().getCellsLevels(),
wnPoints,
wnTri,
wPoints,

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

@ -20,7 +20,11 @@ Licence:
#include "ContactSearch.H"
#include "cellsSimple.H"
#include "cellMapping.H"
#include "NBS.H"
#include "multiGridNBS.H"
#include "multiGridMapping.H"
template class pFlow::ContactSearch<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 cellIterator = typename NBSLevel0Type::cellIterator;
using IdType = typename NBSLevel0Type::IdType;
using IndexType = typename NBSLevel0Type::IndexType;
@ -81,8 +83,9 @@ public:
TypeNameNV("NBS");
NBS(
dictionary dict,
const dictionary& dict,
const box& domain,
real minSize,
real maxSize,
const ViewType1D<realx3, memory_space>& position,
const ViewType1D<real, memory_space>& diam)
@ -126,11 +129,26 @@ public:
return performedSearch_;
}
auto getCellIterator()const
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();

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

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

@ -2,84 +2,120 @@
#define __NBSLevel_H__
#include "cells.H"
#include "contactSearchFunctions.H"
#include "NBSLevel0.H"
namespace pFlow
{
INLINE_FUNCTION_HD
int32x3 mapIndexLevels(const int32x3& ind, int32 lowerLevel, int32 upperLevel);
template<typename executionSpace>
class
NBSLevel
:
public NBSLevel0<executionSpace>
{
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<
Kokkos::IndexType<int32>,
Kokkos::Schedule<Kokkos::Static>,
execution_space>
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:
int32x3 numCells_{1,1,1};
ViewType3D<int32, memory_space> head_;
int8 level_ = 0;
int32 level_ = 0;
public:
NBSLevel()
{}
TypeNameNV("NBSLevel0");
NBSLevel(int8 lvl, int32x3 numCells)
INLINE_FUNCTION_HD
NBSLevel(){}
NBSLevel(
int32 lvl,
const box& domain,
real cellSize,
real sizeRatio,
const ViewType1D<realx3, memory_space>& position,
const ViewType1D<real, memory_space>& diam)
:
numCells_(gridExtent),
head_("NBSLevel::head", numCells_.x(), numCells_.y(), numCells_.z()),
NBSLevel0Type(
domain,
cellSize,
sizeRatio,
position,
diam,
lvl==0),
level_(lvl)
{}
INLINE_FUNCION_HD
auto& head(int32 i, int32 j, int32 k)
{
return head_(i,j,k);
}
INLINE_FUNCTION_HD
NBSLevel(const NBSLevel&) = default;
INLINE_FUNCION_HD
auto& head()
{
return head_;
}
INLINE_FUNCTION_HD
NBSLevel& operator = (const NBSLevel&) = default;
INLINE_FUNCION_HD
INLINE_FUNCTION_HD
~NBSLevel() = default;
INLINE_FUNCTION_HD
auto level()const
{
return level_;
}
INLINE_FUNCION_HD
const auto& numCells()const
template<typename PairsContainer>
INLINE_FUNCTION_H
int32 findPairsCountCross(PairsContainer& pairs, NBSLevel& upperLevel)
{
return gridExtent_;
mdrPolicyFindPairs
mdrPolicy(
{0,0,0},
{this->nx(),this->ny(),this->nz()} );
int32 notInsertedPairs;
Kokkos::parallel_reduce (
"NBSLevel::findPairsCountCross",
mdrPolicy,
CLASS_LAMBDA_HD(int32 i, int32 j, int32 k, int32& getFullUpdate){
#include "NBSCrossLoop.H"
}, notInsertedPairs);
return notInsertedPairs;
}
void nullify()
{
fill(
head_,
range(0,numCells_.x()),
range(0,numCells_.y()),
range(0,numCells_.z()),
static_cast<int32>(-1)
);
}
};
INLINE_FUNCION_HD
int32x3 mapIndexLevels( int32x3 ind, int32 lowerLevel, int32 upperLevel)
INLINE_FUNCTION_HD
int32x3 mapIndexLevels( const int32x3& ind, int32 lowerLevel, int32 upperLevel)
{
int32 a = pow(2, static_cast<int32>(upperLevel-lowerLevel));
return ind/a;

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

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

398
src/Interaction/contactSearch/methods/NBSLevels.H
View File

@ -1,51 +1,58 @@
#ifndef __NBSLevels_H__
#define __NBSLevels_H__
#include "NBSLevel.H"
#include "NBSLevel0.H"
#include "KokkosTypes.H"
namespace pFlow
{
template<typename executionSpace>
class NBSLevels
:
public NBSLevel0<executionSpace>
{
public:
using NBSLevel0Type = NBSLevel0<executionSpace>;
using NBSLevelType = NBSLevel<executionSpace>;
using IdType = typename NBSLevel0Type::IdType;
using cellIterator = typename NBSLevelType::cellIterator;
using IndexType = typename NBSLevel0Type::IndexType;
using IdType = typename NBSLevelType::IdType;
using Cells = typename NBSLevel0Type::Cells;
using IndexType = typename NBSLevelType::IndexType;
using Cells = typename NBSLevelType::Cells;
using CellType = typename Cells::CellType;
using execution_space = typename NBSLevel0Type::execution_space;
using execution_space = typename NBSLevelType::execution_space;
using memory_space = typename NBSLevel0Type::memory_space;
using memory_space = typename NBSLevelType::memory_space;
using realRange = kPair<real,real>;
protected:
real minSize_;
real maxSize_;
int32 numLevles_=1;
int32 numLevels_=1;
ViewType1D<NBSLevel, memory_space> nbsSuperLevels_;
Vector<NBSLevelType> nbsLevels_;
ViewType1D<realRange, memory_space> sizeRangeLevels_;
ViewTypeID<int8, memory_space> particleLevel_;
ViewType1D<realRange, HostSpace> sizeRangeLevelsHost_;
range activeRange_;
ViewType1D<real, memory_space> maxSizeLevels_;
ViewType1D<real, HostSpace> maxSizeLevelsHost_;
ViewType1D<int8, memory_space> particleLevel_;
range activeRange_{0,0};
using rangePolicyType =
Kokkos::RangePolicy<
@ -53,76 +60,115 @@ protected:
Kokkos::Schedule<Kokkos::Static>,
execution_space>;
bool setNumLevels()
int32 setNumLevels()
{
int32 maxOvermin = static_cast<int32>(maxSize_/minSize_);
if (maxOvermin <=1)
numLevles_ = 1;
return 1;
else if(maxOvermin<=3)
numLevles_ = 2;
return 2;
else if(maxOvermin<=7)
numLevles_ = 3;
return 3;
else if(maxOvermin<15)
numLevles_ = 4;
return 4;
else if(maxOvermin<31)
numLevles_ = 5;
return 5;
else if(maxOvermin<63)
numLevles_ = 6;
return 6;
else if(maxOvermin <127)
numLevles_ = 7;
return 7;
else
{
fatalErrorInFunction<<
"size ratio is not valid for multi-grid NBS "<< maxOvermin<<endl;
return false;
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 setDiameterRange()
bool initLevels(
const box& domain,
real sizeRatio,
const ViewType1D<realx3, memory_space>& position,
const ViewType1D<real, memory_space>& diam)
{
auto sizeRatio = this->sizeRatio();
real lvl_minD = 0.1*minSize_;
real lvl_maxD = sizeRatio* minSize_;
reallocNoInit(sizeRangeLevels_, numLevles_);
for(int32 i=0; i<numLevles_; i++)
for(int32 lvl = 0; lvl<numLevels_; lvl++)
{
sizeRangeLevels_[i] = {lvl_minD, lvl_maxD};
lvl_minD = lvl_maxD;
lvl_maxD *= 2.0;
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;
}
bool initSuperLevels()
{
if(numLevles_== 1) return true;
reallocNoInit(nbsSuperLevels_, numLevles_);
if(int lvl = 1; lvl<numLevles_-1; lvl++)
void manageAllocateNext(range active)
{
Cells lvlCells(this->domain(), sizeRangeLevels_[lvl].second);
nbsSuperLevels_[lvl] = NBSLevle(lvl, lvlCells.numCells());
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);
}
}
}
return true;
}
INLINE_FUNCTION_H
void nullify()
void nullify( range active)
{
// level 0
NBSLevel0Type::nullify(activeRange_);
// super levels
for(int32 lvl=1; lvl<numLevles_; lvl++)
for(int32 lvl=0; lvl<numLevels_; lvl++)
{
nbsSuperLevels_[lvl].nullify();
nbsLevels_[lvl].nullify(active);
}
}
@ -134,66 +180,152 @@ public:
real maxSize,
real sizeRatio,
const ViewType1D<realx3, memory_space>& position,
const ViewType1D<real, memory_space>& diam),
range active
const ViewType1D<real, memory_space>& diam)
:
NBSLevel0Type(
domain,
sizeRatio*minSize,
position,
diam),
minSize_(minSize),
maxSize_(maxSize),
activeRange_(active)
numLevels_(setNumLevels()),
nbsLevels_("nbsLevels", numLevels_, numLevels_, RESERVE()),
sizeRangeLevels_("sizeRangeLevels", numLevels_),
sizeRangeLevelsHost_("sizeRangeLevelsHost", numLevels_),
maxSizeLevels_("maxSizeLevels", numLevels_),
maxSizeLevelsHost_("maxSizeLevelsHost", numLevels_)
{
if(!setNumLevels())
{
fatalExit;
setDiameterRange(sizeRatio);
initLevels(domain, sizeRatio, position, diam);
}
setDiameterRange();
initSuperLevels();
findParticleLevel(activeRange_.first, activeRange_.last);
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)
{
// check for active ragne************************************************************************
checkAllocateNext(activeRange.second);
nullify();
// nullify next and heads
findParticleLevel(activeRange.first, activeRange.second);
manageAllocateNext(activeRange);
nullify(activeRange);
//
rangePolicyType rPolicy(activeRange.first, activeRange.second);
auto nbsLevel0 = nbsLevels_[0];
auto pointPosition = nbsLevel0.pointPosition();
auto particleLevel = particleLevel_;
auto next = nbsLevel0.next();
auto head0 = nbsLevel0.head();
typename NBSLevelType::HeadType head1, head2, head3, head4, head5, head6;
if(numLevels_>1) head1 = nbsLevels_[1].head();
if(numLevels_>2) head2 = nbsLevels_[2].head();
if(numLevels_>3) head3 = nbsLevels_[3].head();
if(numLevels_>4) head4 = nbsLevels_[4].head();
if(numLevels_>5) head5 = nbsLevels_[5].head();
if(numLevels_>6) head6 = nbsLevels_[6].head();
Kokkos::parallel_for(
"NBSLevels::build",
rPolicy,
CLASS_LAMBDA_HD(int32 i){
CellType ind = this->pointIndex(points[i]);
int8 lvl = particleLevel_[i];
LAMBDA_HD(int32 i){
int8 lvl = particleLevel[i];
auto ind = nbsLevel0.pointIndex(pointPosition[i]);
ind = mapIndexLevels(ind, 0, lvl);
int32 old;
if( lvl == 0)
{
int32 old =
Kokkos::atomic_exchange(
&this->head_(ind.x(), ind.y(), ind.z()),
i);
this->next_[i] = old;
}
else
{
ind = upper.mapIndexUpperLevel(ind, 0, lvl);
int32 old =
Kokkos::atomic_exchange(
&nbsSuperLevels_[lvl].head()(ind.x(), ind.y(), ind.z()),
i);
this->next_[i] = old;
}
if(lvl==0)
old =Kokkos::atomic_exchange(&head0(ind.x(), ind.y(), ind.z()),i);
else if(lvl==1)
old =Kokkos::atomic_exchange(&head1(ind.x(), ind.y(), ind.z()),i);
else if(lvl==2)
old =Kokkos::atomic_exchange(&head2(ind.x(), ind.y(), ind.z()),i);
else if(lvl==3)
old =Kokkos::atomic_exchange(&head3(ind.x(), ind.y(), ind.z()),i);
else if(lvl==4)
old =Kokkos::atomic_exchange(&head4(ind.x(), ind.y(), ind.z()),i);
else if(lvl==5)
old =Kokkos::atomic_exchange(&head5(ind.x(), ind.y(), ind.z()),i);
else if(lvl==6)
old =Kokkos::atomic_exchange(&head6(ind.x(), ind.y(), ind.z()),i);
next(i) = old;
});
Kokkos::fence();
@ -203,51 +335,73 @@ public:
INLINE_FUNCTION_H
void build(range activeRange, IncludeFunction incld)
{
// check for active ragne************************************************************************
checkAllocateNext(activeRange.second);
nullify();
//
// nullify next and heads
findParticleLevel(activeRange.first, activeRange.second);
manageAllocateNext(activeRange);
nullify(activeRange);
rangePolicyType rPolicy(activeRange.first, activeRange.second);
auto nbsLevel0 = nbsLevels_[0];
auto pointPosition = nbsLevel0.pointPosition();
auto particleLevel = particleLevel_;
auto next = nbsLevel0.next();
auto head0 = nbsLevel0.head();
typename NBSLevelType::HeadType head1, head2, head3, head4, head5, head6;
if(numLevels_>1) head1 = nbsLevels_[1].head();
if(numLevels_>2) head2 = nbsLevels_[2].head();
if(numLevels_>3) head3 = nbsLevels_[3].head();
if(numLevels_>4) head4 = nbsLevels_[4].head();
if(numLevels_>5) head5 = nbsLevels_[5].head();
if(numLevels_>6) head6 = nbsLevels_[6].head();
Kokkos::parallel_for(
"NBSLevels::build",
rPolicy,
CLASS_LAMBDA_HD(int32 i){
if(incld(i))
{
CellType ind = this->pointIndex(points[i]);
int8 lvl = particleLevel_[i];
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)
{
int32 old =
Kokkos::atomic_exchange(
&this->head_(ind.x(), ind.y(), ind.z()),
i);
this->next_[i] = old;
}
else
{
ind = upper.mapIndexUpperLevel(ind, 0, lvl);
int32 old =
Kokkos::atomic_exchange(
&nbsSuperLevels_[lvl].head()(ind.x(), ind.y(), ind.z()),
i);
this->next_[i] = old;
}
}
if(lvl==0)
old =Kokkos::atomic_exchange(&head0(ind.x(), ind.y(), ind.z()),i);
else if(lvl==1)
old =Kokkos::atomic_exchange(&head1(ind.x(), ind.y(), ind.z()),i);
else if(lvl==2)
old =Kokkos::atomic_exchange(&head2(ind.x(), ind.y(), ind.z()),i);
else if(lvl==3)
old =Kokkos::atomic_exchange(&head3(ind.x(), ind.y(), ind.z()),i);
else if(lvl==4)
old =Kokkos::atomic_exchange(&head4(ind.x(), ind.y(), ind.z()),i);
else if(lvl==5)
old =Kokkos::atomic_exchange(&head5(ind.x(), ind.y(), ind.z()),i);
else if(lvl==6)
old =Kokkos::atomic_exchange(&head6(ind.x(), ind.y(), ind.z()),i);
next(i) = old;
});
Kokkos::fence();
}
bool findParticleLevel(int32 first, int32 last)
{
auto diameter = this->diameter_;
auto sizeRange = sizeRangeLevels_;
if(last > particleLevel_.size())
{
reallocNoInit(particleLevel_,last);
}
auto diameter = nbsLevels_[0].diameter();
auto maxSizes = maxSizeLevels_;
auto particleLevel = particleLevel_;
int8 maxLvl = sizeRangeLevels_.size();
@ -261,13 +415,15 @@ public:
{
for(int8 lvl = 0; lvl<maxLvl; lvl++)
{
if( diameter[i]<= sizeRange[lvl].second )
if( diameter[i]<= maxSizes[lvl] )
{
particleLevel[i] = lvl;
break;
return;
}
}
particleLevel[i] = static_cast<int8>(-1);
});
Kokkos::fence();
return true;
}

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 HeadType = ViewType3D<int32, memory_space>;
using NextType = ViewType1D<int32, memory_space>;
class cellIterator
{
private:
ViewType3D<int32, memory_space> head_;
HeadType head_;
ViewType1D<int32, memory_space> next_;
NextType next_;
public:
@ -85,58 +89,60 @@ protected:
ViewType1D<int32, memory_space> next_;
bool nextOwner_ = true;
// borrowed ownership
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 =
Kokkos::RangePolicy<
Kokkos::IndexType<int32>,
Kokkos::Schedule<Kokkos::Static>,
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)
{
if( capacity_ < newCap)
{
capacity_ = newCap;
if(!nextOwner_)return;
reallocNoInit(next_, capacity_);
}
}
@ -146,14 +152,20 @@ private:
reallocNoInit(head_, this->nx(), this->ny(), this->nz());
}
public:
TypeNameNV("mapperNBS");
INLINE_FUNCTION_HD
mapperNBS(){}
mapperNBS(
const box& domain,
real cellSize,
const ViewType1D<realx3, memory_space>& position)
const ViewType1D<realx3, memory_space>& position,
bool nextOwner = true)
:
Cells(domain, cellSize),
pointPosition_(position),
@ -163,7 +175,8 @@ public:
this->ny(),
this->nz()
),
next_("mapperNBS::next_",1) //,position.size()),
next_("mapperNBS::next_",1), //,position.size()),
nextOwner_(nextOwner)
{
checkAllocateNext(pointPosition_.size());
}
@ -173,12 +186,14 @@ public:
int32 nx,
int32 ny,
int32 nz,
const ViewType1D<realx3, memory_space>& position)
const ViewType1D<realx3, memory_space>& position,
bool nextOwner = true)
:
Cells(domain, nx, ny, nz),
pointPosition_(position),
head_("mapperNBS::head_",nx,ny,nz), //, this->nx(), this->ny(), this->nz()),
next_("mapperNBS::next_",1) //,position.size()),
head_("mapperNBS::head_",nx,ny,nz),
next_("mapperNBS::next_",1),
nextOwner_(nextOwner)
{
checkAllocateNext(pointPosition_.size());
}
@ -205,12 +220,54 @@ public:
return cellIterator(head_, next_);
}
bool particlesSizeChanged(int32 newSize)
bool particlesCapcityChanged(int32 newCap)
{
checkAllocateNext(newSize);
checkAllocateNext(newCap);
return true;
}
INLINE_FUNCTION_HD
auto& head()
{
return head_;
}
INLINE_FUNCTION_HD
auto& next()
{
return next_;
}
INLINE_FUNCTION_HD
const auto& head()const
{
return head_;
}
INLINE_FUNCTION_HD
const auto& next()const
{
return next_;
}
INLINE_FUNCTION_HD
auto& pointPosition()
{
return pointPosition_;
}
INLINE_FUNCTION_H
void setNext(ViewType1D<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)

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

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>
using ViewType1D = Kokkos::View<T*,properties...>;
template<typename T, typename... properties>
using DualViewType1D = Kokkos::View<T*,properties...>;
template<typename T, typename... properties>
using ViewType3D = Kokkos::View<T***,properties...>;

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

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

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

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