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

contact search NBS refactored

This commit is contained in:
hamidrezanorouzi 2022-10-27 14:19:53 +03:30
parent 490577dcd2
commit e6d7fbcda3
12 changed files with 1256 additions and 428 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

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

@ -30,7 +30,7 @@ namespace pFlow
{
template<
template<class, class, class> class BaseMethod,
template<class> class BaseMethod,
template<class, class, class> class WallMapping
>
class ContactSearch
@ -49,10 +49,7 @@ public:
using ParticleContactSearchType =
BaseMethod<
ExecutionSpace,
IdType,
IndexType
>;
ExecutionSpace>;
using WallMappingType =
WallMapping<
@ -141,7 +138,7 @@ public:
bool force = false) override
{
//output<<" ContactSearch::broadSearch::PP before.\n";
if(particleContactSearch_)
{
auto activeRange = this->Particles().activeRange();
@ -162,9 +159,7 @@ public:
}
else
return false;
//output<<" ContactSearch::broadSearch::PP fater.\n";
//output<<" ContactSearch::broadSearch::PW before.\n";
if(wallMapping_)
{
sphereWallTimer_.start();
@ -174,7 +169,7 @@ public:
else
return false;
//output<<" ContactSearch::broadSearch::PW after.\n";
return true;
}

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

@ -22,75 +22,35 @@ Licence:
#ifndef __NBS_H__
#define __NBS_H__
#include "cells.H"
#include "contactSearchFunctions.H"
#include "baseAlgorithms.H"
#include "NBSLevel0.H"
namespace pFlow
{
template<
typename executionSpace,
typename idType,
typename indexType=int32
>
template<typename executionSpace>
class NBS
:
public cells<indexType>
{
public:
using IdType = idType;
using NBSLevel0Type = NBSLevel0<executionSpace>;
using IndexType = indexType;
using IdType = typename NBSLevel0Type::IdType;
using Cells = cells<IndexType>;
using IndexType = typename NBSLevel0Type::IndexType;
using Cells = typename NBSLevel0Type::Cells;
using CellType = typename Cells::CellType;
using ExecutionSpace= executionSpace;
using execution_space = typename NBSLevel0Type::execution_space;
using memory_space = typename ExecutionSpace::memory_space;
using memory_space = typename NBSLevel0Type::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;
real sizeRatio_ = 1.0;
int32 updateFrequency_= 1;
@ -98,74 +58,10 @@ protected:
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>;
NBSLevel0Type NBSLevel0_;
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)
@ -180,79 +76,34 @@ private:
}
}
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,
real maxSize,
const ViewType1D<realx3, memory_space>& position,
const ViewType1D<real, memory_space>& diam,
int32 initialContainerSize = 1
)
const ViewType1D<real, memory_space>& diam)
:
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();
}
sizeRatio_(
max(
dict.getValOrSet<real>("sizeRatio", 1.0),
1.0
)),
updateFrequency_(
max(
dict.getValOrSet<int32>("updateFrequency", 1),
1
)),
NBSLevel0_(
domain,
maxSize*sizeRatio_,
sizeRatio_,
position,
diam)
{}
INLINE_FUNCTION_HD
NBS(const NBS&) = default;
@ -275,19 +126,20 @@ public:
return performedSearch_;
}
const auto& Head()const
auto getCellIterator()const
{
return head_;
return NBSLevel0_.getCellIterator();
}
const auto& Next()const
auto getCells()const
{
return next_;
return NBSLevel0_.getCells();
}
cellIterator getCellIterator()const
bool objectSizeChanged(int32 newSize)
{
return cellIterator(head_, next_);
NBSLevel0_.checkAllocateNext(newSize);
return true;
}
// - Perform the broad search to find pairs
@ -302,13 +154,12 @@ public:
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;
NBSLevel0_.build(activeRange);
NBSLevel0_.findPairs(pairs);
performedSearch_ = true;
return true;
}
@ -325,197 +176,14 @@ public:
if( !performSearch() ) return true;
build(activeRange, incld);
NBSLevel0_.build(activeRange, incld);
findPairs(pairs);
NBSLevel0_.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;
}
};
}

70
src/Interaction/contactSearch/methods/NBSLevel.H Normal file
View File

@ -0,0 +1,70 @@
#ifndef __NBSLevel_H__
#define __NBSLevel_H__
#include "cells.H"
#include "contactSearchFunctions.H"
template<typename executionSpace>
class
NBSLevel
{
public:
using execution_space= executionSpace;
using memory_space = typename execution_space::memory_space;
using rangePolicy = Kokkos::RangePolicy<
Kokkos::IndexType<int32>,
Kokkos::Schedule<Kokkos::Static>,
execution_space>
protected:
int32x3 gridExtent_;
ViewType3D<int32, memory_space> head_;
int8 level_ = 1;
public:
NBSLevel(int8 lvl, int32x3 gridExtent)
:
gridExtent_(gridExtent),
head_("NBSLevel::head", gridExtent.x(), gridExtent.y(), gridExtent.z()),
level_(lvl)
{}
INLINE_FUNCION_HD
auto& head(int32 i, int32 j, int32 k)
{
return head_(i,j,k);
}
INLINE_FUNCION_HD
auto& head()
{
return head_;
}
INLINE_FUNCION_HD
auto level()const
{
return level_;
}
INLINE_FUNCION_HD
const auto& gridExtent()const
{
return gridExtent_;
}
};
#endif

211
src/Interaction/contactSearch/methods/NBSLevel0.H Normal file
View File

@ -0,0 +1,211 @@
/*------------------------------- phasicFlow ---------------------------------
O C enter of
O O E ngineering and
O O M ultiscale modeling of
OOOOOOO F luid flow
------------------------------------------------------------------------------
Copyright (C): www.cemf.ir
email: hamid.r.norouzi AT gmail.com
------------------------------------------------------------------------------
Licence:
This file is part of phasicFlow code. It is a free software for simulating
granular and multiphase flows. You can redistribute it and/or modify it under
the terms of GNU General Public License v3 or any other later versions.
phasicFlow is distributed to help others in their research in the field of
granular and multiphase flows, but WITHOUT ANY WARRANTY; without even the
implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
-----------------------------------------------------------------------------*/
#ifndef __NBSLevel0_H__
#define __NBSLevel0_H__
#include "mapperNBS.H"
namespace pFlow
{
template<typename executionSpace>
class NBSLevel0
:
public mapperNBS<executionSpace>
{
public:
using MapperType = mapperNBS<executionSpace>;
using IdType = typename MapperType::IdType;
using IndexType = typename MapperType::IndexType;
using Cells = typename MapperType::Cells;
using CellType = typename Cells::CellType;
using execution_space = typename MapperType::execution_space;
using memory_space = typename MapperType::memory_space;
struct TagFindPairs{};
protected:
real sizeRatio_ = 1.0;
// borrowed ownership
ViewType1D<real, memory_space> diameter_;
using mdrPolicyFindPairs =
Kokkos::MDRangePolicy<
Kokkos::Rank<3>,
Kokkos::Schedule<Kokkos::Dynamic>,
execution_space>;
private:
static INLINE_FUNCTION_HD
void Swap(int32& x, int32& y)
{
int32 tmp = x;
x = y;
y = tmp;
}
public:
TypeNameNV("NBSLevel0");
NBSLevel0(
const box& domain,
real cellSize,
const ViewType1D<realx3, memory_space>& position,
const ViewType1D<real, memory_space>& diam)
:
MapperType(domain, cellSize, position),
diameter_(diam)
{}
NBSLevel0(
const box& domain,
int32 nx,
int32 ny,
int32 nz,
const ViewType1D<realx3, memory_space>& position,
const ViewType1D<real, memory_space>& diam)
:
MapperType(domain, nx, ny, nz, position),
diameter_(diam)
{ }
NBSLevel0(
const box& domain,
real cellSize,
real sizeRatio,
const ViewType1D<realx3, memory_space>& position,
const ViewType1D<real, memory_space>& diam)
:
MapperType(domain, cellSize, position),
sizeRatio_(sizeRatio),
diameter_(diam)
{}
INLINE_FUNCTION_HD
NBSLevel0(const NBSLevel0&) = default;
INLINE_FUNCTION_HD
NBSLevel0& operator = (const NBSLevel0&) = default;
INLINE_FUNCTION_HD
~NBSLevel0()=default;
// - Perform the broad search to find pairs
// with force = true, perform broad search regardless of
// updateFrequency_ value
// on all the points in the range of [0,numPoints_)
template<typename PairsContainer>
bool broadSearch(PairsContainer& pairs, range activeRange)
{
this->build(activeRange);
findPairs(pairs);
return true;
}
// - Perform the broad search to find pairs,
// ignore particles with incld(i) = true,
// with force = true, perform broad search regardless of
// updateFrequency_ value
template<typename PairsContainer, typename IncludeFunction>
bool broadSearch(PairsContainer& pairs, range activeRange, IncludeFunction incld)
{
this->build(activeRange, incld);
findPairs(pairs);
return true;
}
template<typename PairsContainer>
INLINE_FUNCTION_H
bool findPairs(PairsContainer& pairs)
{
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 (
"NBSLevel0::findPairs",
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,500) ;
auto oldCap = pairs.capacity();
pairs.increaseCapacityBy(len);
Info<< "The contact pair container capacity increased from "<<
oldCap << " to "<<pairs.capacity()<<" in NBSLevel0."<<endInfo;
}
Kokkos::fence();
}
return true;
}
};
} // pFlow
#endif // __NBSLevel0_H__

20
src/Interaction/contactSearch/methods/NBSLoop.H
View File

@ -19,30 +19,30 @@ Licence:
-----------------------------------------------------------------------------*/
int32 m = head_(i,j,k);
int32 m = this->head_(i,j,k);
CellType currentCell(i,j,k);
int32 n = -1;
while( m > -1 )
{
auto p_m = pointPosition_[m];
auto p_m = this->pointPosition_[m];
auto d_m = sizeRatio_* diameter_[m];
// the same cell
n = next_(m);
n = this->next_(m);
while(n >-1)
{
auto p_n = pointPosition_[n];
auto p_n = this->pointPosition_[n];
auto d_n = sizeRatio_*diameter_[n];
if( sphereSphereCheck(p_m, p_n, d_m, d_n) )
{
auto ln = n;
auto lm = m;
//output<<"lm "<< lm <<" ln "<< ln << " i j k "<< int32x3(i,j,k)<<endl;
if(lm>ln) Swap(lm,ln);
if( auto res = pairs.insert(lm,ln); res <0)
{
@ -50,7 +50,7 @@ while( m > -1 )
}
}
n = next_(n);
n = this->next_(n);
}
// neighbor cells
@ -74,11 +74,11 @@ while( m > -1 )
if( this->isInRange(neighborCell) )
{
n = head_(neighborCell.x(), neighborCell.y(), neighborCell.z());
n = this->head_(neighborCell.x(), neighborCell.y(), neighborCell.z());
while( n>-1)
{
auto p_n = pointPosition_[n];
auto p_n = this->pointPosition_[n];
auto d_n = sizeRatio_*diameter_[n];
if(sphereSphereCheck(p_m, p_n, d_m, d_n))
@ -91,11 +91,11 @@ while( m > -1 )
getFullUpdate++;
}
}
n = next_[n];
n = this->next_[n];
}
}
}
m = next_[m];
m = this->next_[m];
}

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

@ -0,0 +1,523 @@
/*------------------------------- 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

331
src/Interaction/contactSearch/methods/mapperNBS.H Normal file
View File

@ -0,0 +1,331 @@
/*------------------------------- phasicFlow ---------------------------------
O C enter of
O O E ngineering and
O O M ultiscale modeling of
OOOOOOO F luid flow
------------------------------------------------------------------------------
Copyright (C): www.cemf.ir
email: hamid.r.norouzi AT gmail.com
------------------------------------------------------------------------------
Licence:
This file is part of phasicFlow code. It is a free software for simulating
granular and multiphase flows. You can redistribute it and/or modify it under
the terms of GNU General Public License v3 or any other later versions.
phasicFlow is distributed to help others in their research in the field of
granular and multiphase flows, but WITHOUT ANY WARRANTY; without even the
implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
-----------------------------------------------------------------------------*/
#ifndef __mapperNBS_H__
#define __mapperNBS_H__
#include "cells.H"
#include "contactSearchFunctions.H"
#include "baseAlgorithms.H"
namespace pFlow
{
template<typename executionSpace>
class mapperNBS
:
public cells<int32>
{
public:
using IdType = int32;
using IndexType = int32;
using Cells = cells<IndexType>;
using CellType = typename Cells::CellType;
using execution_space = executionSpace;
using memory_space = typename execution_space::memory_space;
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;
ViewType3D<int32, memory_space> head_;
ViewType1D<int32, memory_space> next_;
// 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:
void checkAllocateNext(int newCap)
{
if( capacity_ < newCap)
{
capacity_ = newCap;
reallocNoInit(next_, capacity_);
}
}
void allocateHead()
{
reallocNoInit(head_, this->nx(), this->ny(), this->nz());
}
public:
TypeNameNV("mapperNBS");
mapperNBS(
const box& domain,
real cellSize,
const ViewType1D<realx3, memory_space>& position)
:
Cells(domain, cellSize),
pointPosition_(position),
head_(
"mapperNBS::head_",
this->nx(),
this->ny(),
this->nz()
),
next_("mapperNBS::next_",1) //,position.size()),
{
checkAllocateNext(pointPosition_.size());
}
mapperNBS(
const box& domain,
int32 nx,
int32 ny,
int32 nz,
const ViewType1D<realx3, memory_space>& position)
:
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()),
{
checkAllocateNext(pointPosition_.size());
}
INLINE_FUNCTION_HD
mapperNBS(const mapperNBS&) = default;
INLINE_FUNCTION_HD
mapperNBS& operator = (const mapperNBS&) = default;
INLINE_FUNCTION_HD
~mapperNBS()=default;
//// - Methods
INLINE_FUNCTION_HD
auto capacity()const
{
return capacity_;
}
cellIterator getCellIterator()const
{
return cellIterator(head_, next_);
}
bool particlesSizeChanged(int32 newSize)
{
checkAllocateNext(newSize);
return true;
}
// - build based on all points in active range
INLINE_FUNCTION_H
void build(range activeRange)
{
checkAllocateNext(activeRange.second);
nullify(activeRange);
Cells cellIndex = static_cast<Cells>(*this);
auto points = pointPosition_;
auto next = next_;
auto head = head_;
rangePolicyType rPolicy(activeRange.first, activeRange.second);
Kokkos::parallel_for(
"mapperNBS::build",
rPolicy,
LAMBDA_HD(int32 i){
CellType ind = cellIndex.pointIndex(points[i]);
int32 old = Kokkos::atomic_exchange(&head(ind.x(), ind.y(), ind.z()), i);
next[i] = old;
});
Kokkos::fence();
}
template<typename IncludeFunction>
INLINE_FUNCTION_H
void build(range activeRange, IncludeFunction incld)
{
checkAllocateNext(activeRange.second);
nullify(activeRange);
Cells cellIndex = static_cast<Cells>(*this);
auto points = pointPosition_;
auto next = next_;
auto head = head_;
rangePolicyType rPolicy(activeRange.first, activeRange.second);
Kokkos::parallel_for(
"mapperNBS::build_Include",
rPolicy,
LAMBDA_HD(int32 i){
if( incld(i) )
{
CellType ind = cellIndex.pointIndex(points[i]);
auto old = Kokkos::atomic_exchange(&head(ind.x(), ind.y(), ind.z()), i);
next[i] = old;
}
});
Kokkos::fence();
}
INLINE_FUNCTION_H
void buildCheckInDomain(range activeRange)
{
checkAllocateNext(activeRange.second);
nullify(activeRange);
Cells cellIndex = static_cast<Cells>(*this);
auto points = pointPosition_;
auto next = next_;
auto head = head_;
rangePolicyType rPolicy(activeRange.first, activeRange.second);
Kokkos::parallel_for(
"mapperNBS::buildCheckInDomain",
rPolicy,
LAMBDA_HD(int32 i){
CellType ind;
if( cellIndex.pointIndexInDomain(points[i], ind) )
{
int32 old = Kokkos::atomic_exchange(&head(ind.x(), ind.y(), ind.z()), i);
next[i] = old;
}
});
Kokkos::fence();
}
template<typename IncludeFunction>
INLINE_FUNCTION_H
void buildCheckInDomain(range activeRange, IncludeFunction incld)
{
checkAllocateNext(activeRange.second);
nullify(activeRange);
Cells cellIndex = static_cast<Cells>(*this);
auto points = pointPosition_;
auto next = next_;
auto head = head_;
rangePolicyType rPolicy(activeRange.first, activeRange.second);
Kokkos::parallel_for(
"mapperNBS::buildCheckInDomain_Include",
rPolicy,
LAMBDA_HD(int32 i){
CellType ind;
if( incld(i) && cellIndex.pointIndexInDomain(points[i], ind) )
{
auto old = Kokkos::atomic_exchange(&head(ind.x(), ind.y(), ind.z()), i);
next[i] = old;
}
});
Kokkos::fence();
}
};
} // pFlow
#endif // __mapperNBS_H__

1
src/Interaction/sphereInteraction/sphereInteraction.H
View File

@ -136,6 +136,7 @@ public:
{
//Info<<"before contact search"<<endInfo;
////Info<<"interaction iterrate start"<<endInfo;
if(this->contactSearch_)
{

36
src/phasicFlow/Kokkos/KokkosUtilities.H
View File

@ -78,6 +78,42 @@ void reallocFill( ViewType1D<Type,Properties...>& view, int32 len, Type val)
}
template <
typename Type,
typename... Properties>
INLINE_FUNCTION_H
void realloc( ViewType3D<Type,Properties...>& view, int32 len1, int32 len2, int32 len3)
{
Kokkos::realloc(view, len1, len2, len3);
}
template <
typename Type,
typename... Properties>
INLINE_FUNCTION_H
void reallocNoInit(ViewType3D<Type,Properties...>& view, int32 len1, int32 len2, int32 len3)
{
using ViewType = ViewType3D<Type,Properties...>;
word vl = view.label();
view = ViewType(); // Deallocate first
view = ViewType(
Kokkos::view_alloc(
Kokkos::WithoutInitializing,
vl),
len1, len2, len3);
}
template <
typename Type,
typename... Properties>
INLINE_FUNCTION_H
void reallocFill( ViewType3D<Type,Properties...>& view, int32 len1, int32 len2, int32 len3, Type val)
{
reallocNoInit(view, len1, len2, len3);
Kokkos::deep_copy(view, val);
}
template<typename ViewType>
INLINE_FUNCTION_H
void swapViews(ViewType& v1, ViewType &v2)

6
utilities/particlesPhasicFlow/positionRandom/positionRandom.C
View File

@ -21,7 +21,7 @@ Licence:
#include "positionRandom.H"
#include "uniformRandomReal.H"
#include "NBS.H"
#include "NBSLevel0.H"
#include "unsortedPairs.H"
#include "box.H"
@ -30,7 +30,7 @@ Licence:
namespace pFlow
{
using SearchType = NBS<DefaultExecutionSpace, int32> ;
using SearchType = NBSLevel0<DefaultExecutionSpace> ;
using ContainerType = unsortedPairs<DefaultExecutionSpace, int32>;
@ -80,7 +80,7 @@ bool pFlow::positionRandom::positionOnePass(int32 pass, int32 startNum)
fillPoints(startNum, positionHD, flagHD);
search.broadSearch(pairs, range(0, startNum), true);
search.broadSearch(pairs, range(0, startNum));
int32 numCollisions = findCollisions(pairs, flagHD);

10
utilities/postprocessPhasicFlow/fieldOperations.H
View File

@ -35,6 +35,7 @@ rectMeshField_H<T> sumOp( const pointField_H<T> field, const pointRectCell& poin
{
// create field
const auto& mesh = pointToCell.mesh();
auto iterator = pointToCell.getCellIterator();
rectMeshField_H<T> results(mesh, T(0));
@ -44,12 +45,12 @@ rectMeshField_H<T> sumOp( const pointField_H<T> field, const pointRectCell& poin
{
for(int32 k=0; k<mesh.nz(); k++)
{
auto n = pointToCell.Head(i,j,k);
auto n = iterator.start(i,j,k);
T res (0);
while(n>-1)
{
res += field[n];
n = pointToCell.Next(n);
n = iterator.getNext(n);
}
results(i,j,k) = res;
@ -65,6 +66,7 @@ rectMeshField_H<T> sumMaksOp( const pointField_H<T> field, const pointRectCell&
{
// create field
const auto& mesh = pointToCell.mesh();
auto iterator = pointToCell.getCellIterator();
rectMeshField_H<T> results(mesh, T(0));
@ -75,7 +77,7 @@ rectMeshField_H<T> sumMaksOp( const pointField_H<T> field, const pointRectCell&
for(int32 k=0; k<mesh.nz(); k++)
{
//auto [loop, n] = pointToCell.startLoop(i,j,k);
auto n = pointToCell.Head(i,j,k);
auto n = iterator.start(i,j,k);
T res (0);
while(n>-1)
@ -86,7 +88,7 @@ rectMeshField_H<T> sumMaksOp( const pointField_H<T> field, const pointRectCell&
res += field[n];
}
n = pointToCell.Next(n);
n = iterator.getNext(n);
}
results(i,j,k) = res;

33
utilities/postprocessPhasicFlow/pointRectCell.H
View File

@ -21,7 +21,7 @@ Licence:
#ifndef __pointRectCell_H__
#define __pointRectCell_H__
#include "NBS.H"
#include "mapperNBS.H"
#include "rectMeshFields.H"
#include "pointStructure.H"
@ -33,7 +33,7 @@ class pointRectCell
{
public:
using mapType = NBS<DefaultHostExecutionSpace, int32>;
using mapType = mapperNBS<DefaultHostExecutionSpace>;
using memory_space = typename mapType::memory_space;
@ -48,9 +48,6 @@ protected:
ViewType1D<realx3, memory_space> pointPosition_;
ViewType1D<real, memory_space> diameter_;
mapType map_;
int32RectMeshField_H nPointInCell_;
@ -79,14 +76,12 @@ public:
),
pStruct_(pStruct),
pointPosition_(pStruct_.pointPosition().hostVectorAll()),
diameter_("diameter", pStruct_.capacity()),
map_(
mesh_.domain(),
mesh_.nx(),
mesh_.ny(),
mesh_.nz(),
pointPosition_,
diameter_),
pointPosition_),
nPointInCell_(mesh_, 0)
{
@ -111,9 +106,9 @@ public:
map_.buildCheckInDomain(activeRange, activeMask);
auto iterator = map_.getCellIterator();
const auto& Next = map_.Next();
const auto& Head = map_.Head();
for(int32 i=0; i<map_.nx(); i++)
{
for(int32 j=0; j<map_.ny(); j++)
@ -122,11 +117,11 @@ public:
{
int32 res = 0;
int32 n = Head(i,j,k);
int32 n = iterator.start(i,j,k);
while( n>-1)
{
res+=1;
n = Next[n];
n = iterator.getNext(n);
}
nPointInCell_(i,j,k) = res;
@ -136,20 +131,16 @@ public:
}
auto getCellIterator()const
{
return map_.getCellIterator();
}
int32 nPointInCell(int32 i, int32 j, int32 k)const
{
return nPointInCell_(i,j,k);
}
auto inline Head(int32 i, int32 j, int32 k)const
{
return map_.Head()(i,j,k);
}
auto inline Next(int32 n)const
{
return map_.Next()(n);
}
//auto
};