xfygogo/phasicFlow
1
0
Fork 0
mirror of https://github.com/PhasicFlow/phasicFlow.git synced 2025-06-22 16:28:30 +00:00
Code Issues Packages Projects Releases Wiki Activity

Merge branch 'develop' into MPIdev

Browse Source
This commit is contained in:
Hamidreza Norouzi
2024-05-18 19:14:01 +03:30
parent 665879f8ca 4e8b921514
commit e05bd2c350
111 changed files with 2006 additions and 703 deletions
Download Patch File Download Diff File Expand all files Collapse all files

10
solvers/iterateGeometry/iterateGeometry.cpp
View File

@ -35,12 +35,10 @@ Licence:
#include "commandLine.hpp"
//#include "readControlDict.hpp"
using namespace pFlow;
int main( int argc, char* argv[] )
{
commandLine cmds(
pFlow::commandLine cmds(
"iterateGeometry",
"Performs simulation without particles, only geometry is solved");
@ -55,8 +53,8 @@ commandLine cmds(
// this should be palced in each main
processors::initProcessors(argc, argv);
initialize_pFlowProcessors();
pFlow::processors::initProcessors(argc, argv);
pFlow::initialize_pFlowProcessors();
#include "initialize_Control.hpp"
#include "setProperty.hpp"
@ -71,7 +69,7 @@ initialize_pFlowProcessors();
// this should be palced in each main
#include "finalize.hpp"
processors::finalizeProcessors();
pFlow::processors::finalizeProcessors();
}

2
solvers/iterateSphereParticles/createDEMComponents.hpp
View File

@ -20,6 +20,6 @@ Licence:
//
REPORT(0)<<"\nReading sphere particles . . ."<<END_REPORT;
sphereParticles sphParticles(Control, proprties);
pFlow::sphereParticles sphParticles(Control, proprties);
WARNING<<"Particle insertion has not been set yet!"<<END_WARNING;

10
solvers/iterateSphereParticles/iterateSphereParticles.cpp
View File

@ -36,7 +36,7 @@ Licence:
#include "systemControl.hpp"
#include "commandLine.hpp"
using namespace pFlow;
/**
* DEM solver for simulating granular flow of cohesion-less particles.
@ -47,7 +47,7 @@ using namespace pFlow;
int main( int argc, char* argv[])
{
commandLine cmds(
pFlow::commandLine cmds(
"sphereGranFlow",
"DEM solver for non-cohesive spherical particles with particle insertion "
"mechanism and moving geometry");
@ -57,8 +57,8 @@ bool isCoupling = false;
if(!cmds.parse(argc, argv)) return 0;
// this should be palced in each main
processors::initProcessors(argc, argv);
initialize_pFlowProcessors();
pFlow::processors::initProcessors(argc, argv);
pFlow::initialize_pFlowProcessors();
#include "initialize_Control.hpp"
#include "setProperty.hpp"
@ -82,7 +82,7 @@ initialize_pFlowProcessors();
// this should be palced in each main
#include "finalize.hpp"
processors::finalizeProcessors();
pFlow::processors::finalizeProcessors();
}

6
solvers/sphereGranFlow/createDEMComponents.hpp
View File

@ -20,7 +20,7 @@ Licence:
//
REPORT(0)<<"\nReading sphere particles . . ."<<END_REPORT;
sphereParticles sphParticles(Control, proprties);
pFlow::sphereParticles sphParticles(Control, proprties);
//
REPORT(0)<<"\nCreating particle insertion object . . ."<<END_REPORT;
@ -36,12 +36,12 @@ REPORT(0)<<"\nCreating particle insertion object . . ."<<END_REPORT;
sphParticles.shapes()
);*/
auto sphInsertion = sphereInsertion(
auto sphInsertion = pFlow::sphereInsertion(
sphParticles,
sphParticles.spheres());
REPORT(0)<<"\nCreating interaction model for sphere-sphere contact and sphere-wall contact . . ."<<END_REPORT;
auto interactionPtr = interaction::create(
auto interactionPtr = pFlow::interaction::create(
Control,
sphParticles,
surfGeometry

13
solvers/sphereGranFlow/sphereGranFlow.cpp
View File

@ -40,12 +40,8 @@ Licence:
#include "interaction.hpp"
#include "Insertions.hpp"
//#include "readControlDict.hpp"
using namespace pFlow;
/**
* DEM solver for simulating granular flow of cohesion-less particles.
*
@ -55,7 +51,7 @@ using namespace pFlow;
int main( int argc, char* argv[])
{
commandLine cmds(
pFlow::commandLine cmds(
"sphereGranFlow",
"DEM solver for non-cohesive spherical particles with particle insertion "
"mechanism and moving geometry");
@ -65,14 +61,13 @@ bool isCoupling = false;
if(!cmds.parse(argc, argv)) return 0;
// this should be palced in each main
processors::initProcessors(argc, argv);
initialize_pFlowProcessors();
pFlow::processors::initProcessors(argc, argv);
pFlow::initialize_pFlowProcessors();
#include "initialize_Control.hpp"
#include "setProperty.hpp"
#include "setSurfaceGeometry.hpp"
#include "createDEMComponents.hpp"
REPORT(0)<<"\nStart of time loop . . .\n"<<END_REPORT;
@ -117,7 +112,7 @@ initialize_pFlowProcessors();
// this should be palced in each main
#include "finalize.hpp"
processors::finalizeProcessors();
pFlow::processors::finalizeProcessors();
}

4
src/Geometry/geometry/geometry.cpp
View File

@ -249,7 +249,9 @@ pFlow::geometry::geometry
if( this->numSurfaces() != motionComponentName_.size() )
{
fatalErrorInFunction<<
"Number of surfaces is not equal to number of motion component names"<<endl;
"Number of surfaces ("<< this->numSurfaces() <<
") is not equal to number of motion component names("<<
motionComponentName_.size()<<")"<<endl;
fatalExit;
}

34
src/Geometry/geometryMotion/geometryMotion.cpp
View File

@ -1,4 +1,3 @@
#include "geometryMotion.hpp"
/*------------------------------- phasicFlow ---------------------------------
O C enter of
O O E ngineering and
@ -70,6 +69,28 @@ bool pFlow::geometryMotion<MotionModel>::findMotionIndex()
return true;
}
namespace pFlow::GMotion
{
template<typename ModelInterface>
void moveGeometry(
real dt,
uint32 nPoints,
const ModelInterface& mModel,
const deviceViewType1D<uint32>& pointMIndexD,
const deviceViewType1D<realx3>& pointsD
)
{
Kokkos::parallel_for(
"geometryMotion<MotionModel>::movePoints",
deviceRPolicyStatic(0, nPoints),
LAMBDA_HD(uint32 i){
auto newPos = mModel.transferPoint(pointMIndexD[i], pointsD[i], dt);
pointsD[i] = newPos;
});
Kokkos::fence();
}
}
template<typename MotionModel>
bool pFlow::geometryMotion<MotionModel>::moveGeometry()
@ -84,8 +105,15 @@ template<typename MotionModel>
auto& pointMIndexD= pointMotionIndex_.deviceViewAll();
auto& pointsD = points().deviceViewAll();
pFlow::GMotion::moveGeometry(
dt,
numPoints(),
motionModel_.getModelInterface(iter, t, dt),
pointMotionIndex_.deviceViewAll(),
points().deviceViewAll()
);
Kokkos::parallel_for(
/*Kokkos::parallel_for(
"geometryMotion<MotionModel>::movePoints",
deviceRPolicyStatic(0, numPoints()),
LAMBDA_HD(uint32 i){
@ -93,7 +121,7 @@ template<typename MotionModel>
pointsD[i] = newPos;
});
Kokkos::fence();
Kokkos::fence();*/
// move the motion components
motionModel_.move(iter, t,dt);

4
src/Interaction/CMakeLists.txt
View File

@ -7,8 +7,8 @@ contactSearch/methods/cellBased/NBS/NBS.cpp
contactSearch/methods/cellBased/NBS/cellsWallLevel0.cpp
contactSearch/boundaries/boundaryContactSearch/boundaryContactSearch.cpp
contactSearch/boundaries/twoPartContactSearch/twoPartContactSearchKernels.cpp
contactSearch/boundaries/twoPartContactSearch/twoPartContactSearch.cpp
#contactSearch/boundaries/twoPartContactSearch/twoPartContactSearchKernels.cpp
#contactSearch/boundaries/twoPartContactSearch/twoPartContactSearch.cpp
contactSearch/boundaries/periodicBoundaryContactSearch/ppwBndryContactSearchKernels.cpp
contactSearch/boundaries/periodicBoundaryContactSearch/ppwBndryContactSearch.cpp
contactSearch/boundaries/periodicBoundaryContactSearch/wallBoundaryContactSearch.cpp

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

@ -138,7 +138,7 @@ public:
start,
end,
newPair);
idx0!=-1)
idx0!=static_cast<uint32>(-1))
{
values_[idx] = values0_[idx0];
}
@ -147,7 +147,7 @@ public:
start,
end,
newPair);
idx0!=-1)
idx0!= static_cast<uint32>(-1) )
{
values_[idx] = values0_[idx0];

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

@ -95,7 +95,8 @@ public:
// swap conainer and values
swapViews(values0_, values_);
swapViews(container0_, this->container_);
return UnsortedPairs::beforeBroadSearch();
UnsortedPairs::beforeBroadSearch();
return true;
}
bool afterBroadSearch()
@ -123,7 +124,7 @@ public:
INLINE_FUNCTION_HD
bool getValue(const PairType& p, ValueType& val)const
{
if(auto idx = this->find(p); idx!=-1)
if(auto idx = this->find(p); idx!=static_cast<uint32>(-1))
{
val = getValue(idx);
return true;
@ -140,7 +141,7 @@ public:
INLINE_FUNCTION_HD
bool setValue(const PairType& p, const ValueType& val)const
{
if(uint32 idx = this->find(p); idx!=-1)
if(uint32 idx = this->find(p); idx!=static_cast<uint32>(-1))
{
setValue(idx, val);
return true;;
@ -155,7 +156,7 @@ public:
{
if( uint32 idx0 =
container0_.find(this->getPair(idx));
idx0!=-1 )
idx0!= static_cast<uint32>(-1) )
{
values_[idx] = values0_[idx0];
}

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

@ -105,7 +105,7 @@ public:
uint32 insert(idType i, idType j)const
{
if(auto insertResult = container_.insert(PairType(i,j)); insertResult.failed())
return -1;
return static_cast<uint32>(-1);
else
return insertResult.index();
@ -115,7 +115,7 @@ public:
uint32 insert(const PairType& p)const
{
if(auto insertResult = container_.insert(p); insertResult.failed())
return -1;
return static_cast<uint32>(-1);
else
return insertResult.index();
@ -154,7 +154,7 @@ public:
idx != Kokkos::UnorderedMapInvalidIndex )
return idx;
else
return -1;
return static_cast<uint32>(-1);
}
INLINE_FUNCTION_HD

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

@ -160,6 +160,7 @@ public:
csPairContainerType& pwPairs,
bool force = false)override
{
if(i==0u)
Particles().boundingSphere().updateBoundaries(DataDirection::SlaveToMaster);
return csBoundaries_[i].broadSearch(
iter,
@ -184,6 +185,23 @@ public:
return ppwContactSearch_().performedSearch();
}
uint32 updateInterval()const override
{
return ppwContactSearch_().updateInterval();
}
real sizeRatio()const override
{
return ppwContactSearch_().sizeRatio();
}
real cellExtent()const override
{
return ppwContactSearch_().cellExtent();
}
};
}

5
src/Interaction/contactSearch/boundaries/boundaryContactSearch/boundaryContactSearch.hpp
View File

@ -72,11 +72,6 @@ public:
return contactSearch_;
}
void fill(const std::any &val) override
{
return;
}
bool hearChanges(
real t,
real dt,

4
src/Interaction/contactSearch/boundaries/periodicBoundaryContactSearch/ppwBndryContactSearchKernels.cpp
View File

@ -80,14 +80,14 @@ pFlow::uint32 pFlow::pweBndryContactSearchKernels::broadSearchPP
if(!searchCells.inCellRange(ind))continue;
uint32 thisI = head(ind.x(),ind.y(),ind.z());
while (thisI!=-1)
while (thisI!=static_cast<uint32>(-1))
{
auto d_n = sizeRatio*diams[thisI];
// first item is for this boundary and second itme, for mirror
if(sphereSphereCheckB(p_m, points[thisI], d_m, d_n)&&
ppPairs.insert(thisI,mrrI) == -1)
ppPairs.insert(thisI,mrrI) == static_cast<uint32>(-1))
{
getFullUpdate++;
}

11
src/Interaction/contactSearch/boundaries/periodicBoundaryContactSearch/wallBoundaryContactSearch.cpp
View File

@ -113,6 +113,9 @@ pFlow::uint32 pFlow::wallBoundaryContactSearch::findPairsElementRangeCount
uint32 nNotInserted = 0;
uint32 nThis = pPoints.size();
const auto& numElements = numElements_;
const auto& elementBox = elementBox_;
const auto& validBox = validBox_;
Kokkos::parallel_reduce(
"pFlow::wallBoundaryContactSearch::findPairsElementRangeCount",
@ -123,11 +126,11 @@ pFlow::uint32 pFlow::wallBoundaryContactSearch::findPairsElementRangeCount
int32x3 ind;
if( searchCells.pointIndexInDomain(p, ind) )
{
for(uint32 nTri=0; nTri<numElements_; nTri++)
for(uint32 nTri=0; nTri<numElements; nTri++)
{
if( validBox_[nTri]== 0)continue;
if( elementBox_[nTri].isInside(ind)&&
pairs.insert(i,nTri+baseTriIndex) == -1)
if( validBox[nTri]== 0)continue;
if( elementBox[nTri].isInside(ind)&&
pairs.insert(i,nTri+baseTriIndex) == static_cast<uint32>(-1))
{
notInsertedUpdate++;
}

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

@ -140,6 +140,15 @@ public:
virtual
bool performedBroadSearch()const = 0;
virtual
uint32 updateInterval()const = 0;
virtual
real sizeRatio()const = 0;
virtual
real cellExtent()const = 0;
static
uniquePtr<contactSearch> create(
const dictionary& dict,

10
src/Interaction/contactSearch/methods/cellBased/NBS/NBS.hpp
View File

@ -135,6 +135,16 @@ public:
return 1;
}
real sizeRatio()const
{
return sizeRatio_;
}
real cellExtent()const
{
return cellExtent_;
}
auto getCellIterator([[maybe_unused]] uint32 lvl)const
{
return NBSLevel0_.getCellIterator();

4
src/Interaction/contactSearch/methods/cellBased/NBS/NBSLoop.hpp
View File

@ -43,7 +43,7 @@ while( m != mapperNBS::NoPos)
auto lm = m;
if(lm>ln) Swap(lm,ln);
if( pairs.insert(lm,ln) == -1)
if( pairs.insert(lm,ln) == static_cast<uint32>(-1))
{
getFullUpdate++;
}
@ -86,7 +86,7 @@ while( m != mapperNBS::NoPos)
auto ln = n;
auto lm = m;
if(lm>ln) Swap(lm,ln);
if( pairs.insert(lm,ln) == -1)
if( pairs.insert(lm,ln) == static_cast<uint32>(-1))
{
getFullUpdate++;
}

35
src/Interaction/contactSearch/methods/cellBased/NBS/cellsWallLevel0.cpp
View File

@ -85,21 +85,26 @@ bool pFlow::cellsWallLevel0::broadSearch
bool pFlow::cellsWallLevel0::build(const cells & searchBox)
{
const auto& points = points_;
const auto& vertices = vertices_;
const auto& elementBox = elementBox_;
const auto cellExtent = cellExtent_;
Kokkos::parallel_for(
"pFlow::cellsWallLevel0::build",
deviceRPolicyStatic(0,numElements_),
CLASS_LAMBDA_HD(uint32 i)
LAMBDA_HD(uint32 i)
{
auto v = vertices_[i];
auto p1 = points_[v.x()];
auto p2 = points_[v.y()];
auto p3 = points_[v.z()];
auto v = vertices[i];
auto p1 = points[v.x()];
auto p2 = points[v.y()];
auto p3 = points[v.z()];
realx3 minP;
realx3 maxP;
searchBox.extendBox(p1, p2, p3, cellExtent_, minP, maxP);
elementBox_[i] = iBoxType(searchBox.pointIndex(minP), searchBox.pointIndex(maxP));
searchBox.extendBox(p1, p2, p3, cellExtent, minP, maxP);
elementBox[i] = iBoxType(searchBox.pointIndex(minP), searchBox.pointIndex(maxP));
});
Kokkos::fence();
@ -153,6 +158,11 @@ pFlow::int32 pFlow::cellsWallLevel0::findPairsElementRangeCount
{
uint32 getFull =0;
const auto& elementBox = elementBox_;
const auto& normals = normals_;
const auto& points = points_;
const auto& vertices = vertices_;
const auto cellExtent = cellExtent_;
Kokkos::parallel_reduce(
"pFlow::cellsWallLevel0::findPairsElementRangeCount",
@ -163,10 +173,10 @@ pFlow::int32 pFlow::cellsWallLevel0::findPairsElementRangeCount
const uint32 iTri = teamMember.league_rank();
const auto triBox = elementBox_[iTri];
const auto triBox = elementBox[iTri];
const auto triPlane = infinitePlane(
normals_[iTri],
points_[vertices_[iTri].x()]);
normals[iTri],
points[vertices[iTri].x()]);
uint32 getFull2 = 0;
@ -186,11 +196,12 @@ pFlow::int32 pFlow::cellsWallLevel0::findPairsElementRangeCount
while( n != particleMap.NoPos)
{
// id is wall id the pair is (particle id, wall id)
if( abs(triPlane.pointFromPlane(pPoints[n]))< pDiams[n]*sizeRatio*cellExtent_)
if( abs(triPlane.pointFromPlane(pPoints[n]))< pDiams[n]*sizeRatio*cellExtent)
{
if( pairs.insert(
static_cast<csIdType>(n),
static_cast<csIdType>(iTri) ) == -1 )
static_cast<csIdType>(iTri) ) == static_cast<csIdType>(-1)
)
innerUpdate++;
}
n = particleMap.next(n);

15
src/Interaction/contactSearch/methods/particleWallContactSearchs.hpp
View File

@ -113,6 +113,21 @@ public:
return false;
}
uint32 updateInterval()const
{
return updateInterval_;
}
real sizeRatio()const
{
return getMethod().sizeRatio();
}
real cellExtent()const
{
return getMethod().cellExtent();
}
};

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

@ -72,7 +72,7 @@ pFlow::uniquePtr<pFlow::interaction> pFlow::interaction::create
gSettings::sleepMiliSeconds(
1000*(pFlowProcessors().localSize()-pFlowProcessors().localRank()-1));
100*(pFlowProcessors().localSize()-pFlowProcessors().localRank()-1));
pOutput.space(2)<<"Creating interaction "<<Green_Text(interactionModel)<<" . . ."<<END_REPORT;
if( systemControlvCtorSelector_.search(interactionModel) )
{

17
src/Interaction/sphereInteraction/boundaries/boundarySphereInteraction/boundarySphereInteraction.cpp
View File

@ -1,3 +1,4 @@
#include "boundarySphereInteraction.hpp"
/*------------------------------- phasicFlow ---------------------------------
O C enter of
O O E ngineering and
@ -18,6 +19,20 @@ Licence:
-----------------------------------------------------------------------------*/
template <typename cFM, typename gMM>
void pFlow::boundarySphereInteraction<cFM, gMM>::allocatePPPairs()
{
ppPairs_.reset(nullptr);
ppPairs_ = makeUnique<ContactListType>(1);
}
template <typename cFM, typename gMM>
void pFlow::boundarySphereInteraction<cFM, gMM>::allocatePWPairs()
{
pwPairs_.reset(nullptr);
pwPairs_ = makeUnique<ContactListType>(1);
}
template <typename cFM, typename gMM>
pFlow::boundarySphereInteraction<cFM, gMM>::boundarySphereInteraction(
@ -28,8 +43,6 @@ pFlow::boundarySphereInteraction<cFM, gMM>::boundarySphereInteraction(
geometryMotion_(geomMotion),
sphParticles_(sphPrtcls)
{
ppPairs_ = makeUnique<ContactListType>(1);
pwPairs_ = makeUnique<ContactListType>(1);
}
template <typename cFM, typename gMM>

38
src/Interaction/sphereInteraction/boundaries/boundarySphereInteraction/boundarySphereInteraction.hpp
View File

@ -22,7 +22,7 @@ Licence:
#include "virtualConstructor.hpp"
#include "generalBoundary.hpp"
#include "unsortedContactList.hpp"
#include "sortedContactList.hpp"
#include "sphereParticles.hpp"
namespace pFlow
@ -51,7 +51,7 @@ public:
using IndexType = uint32;
using ContactListType =
unsortedContactList<ModelStorage, DefaultExecutionSpace, IdType>;
sortedContactList<ModelStorage, DefaultExecutionSpace, IdType>;
private:
@ -60,9 +60,15 @@ private:
/// const reference to sphere particles
const sphereParticles& sphParticles_;
uniquePtr<ContactListType> ppPairs_;
uniquePtr<ContactListType> ppPairs_ = nullptr;
uniquePtr<ContactListType> pwPairs_;
uniquePtr<ContactListType> pwPairs_ = nullptr;
protected:
void allocatePPPairs();
void allocatePWPairs();
public:
@ -124,15 +130,30 @@ public:
return pwPairs_();
}
bool ppPairsAllocated()const
{
if( ppPairs_)return true;
return false;
}
bool pwPairsAllocated()const
{
if( pwPairs_)return true;
return false;
}
virtual
bool sphereSphereInteraction(
real dt,
const ContactForceModel& cfModel)
const ContactForceModel& cfModel,
uint32 step)
{
// for default boundary, no thing to be done
return true;
return false;
}
bool hearChanges
(
real t,
@ -149,11 +170,6 @@ public:
return true;
}
void fill(const std::any& val)override
{
notImplementedFunction;
}
static
uniquePtr<BoundarySphereInteractionType> create(
const boundaryBase& boundary,

10
src/Interaction/sphereInteraction/boundaries/periodicBoundarySphereInteraction/periodicBoundarySphereInteraction.cpp
View File

@ -31,6 +31,9 @@ pFlow::periodicBoundarySphereInteraction<cFM, gMM>::periodicBoundarySphereIntera
if(boundary.thisBoundaryIndex()%2==1)
{
masterInteraction_ = true;
this->allocatePPPairs();
this->allocatePWPairs();
}
else
{
@ -42,10 +45,11 @@ template <typename cFM, typename gMM>
bool pFlow::periodicBoundarySphereInteraction<cFM, gMM>::sphereSphereInteraction
(
real dt,
const ContactForceModel &cfModel
const ContactForceModel &cfModel,
uint32 step
)
{
if(!masterInteraction_) return true;
if(!masterInteraction_) return false;
pFlow::periodicBoundarySIKernels::sphereSphereInteraction(
dt,
@ -61,5 +65,5 @@ bool pFlow::periodicBoundarySphereInteraction<cFM, gMM>::sphereSphereInteraction
this->sphParticles().contactForce().deviceViewAll(),
this->sphParticles().contactTorque().deviceViewAll());
return true;
return false;
}

3
src/Interaction/sphereInteraction/boundaries/periodicBoundarySphereInteraction/periodicBoundarySphereInteraction.hpp
View File

@ -83,7 +83,8 @@ public:
bool sphereSphereInteraction(
real dt,
const ContactForceModel& cfModel)override;
const ContactForceModel& cfModel,
uint32 step)override;
};

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

@ -137,7 +137,9 @@ pFlow::sphereInteraction<cFM,gMM, cLT>::sphereInteraction
ppInteractionTimer_("sphere-sphere interaction", &this->timers()),
pwInteractionTimer_("sphere-wall interaction", &this->timers()),
contactListMangementTimer_("contact-list management", &this->timers()),
boundaryInteractionTimer_("interaction for boundary", &this->timers())
boundaryContactSearchTimer_("contact search for boundary", &this->timers()),
boundaryInteractionTimer_("interaction for boundary", &this->timers()),
contactListBoundaryTimer_("contact-list management for boundary", &this->timers())
{
if(!createSphereInteraction())
@ -179,11 +181,14 @@ bool pFlow::sphereInteraction<cFM,gMM, cLT>::iterate()
contactListMangementTimer_.pause();
}
contactListBoundaryTimer_.start();
for(uint32 i=0; i<6u; i++)
{
boundaryInteraction_[i].ppPairs().beforeBroadSearch();
boundaryInteraction_[i].pwPairs().beforeBroadSearch();
auto& BI = boundaryInteraction_[i];
if(BI.ppPairsAllocated()) BI.ppPairs().beforeBroadSearch();
if(BI.pwPairsAllocated()) BI.pwPairs().beforeBroadSearch();
}
contactListBoundaryTimer_.pause();
if( sphParticles_.numActive()<=0)return true;
@ -200,21 +205,27 @@ bool pFlow::sphereInteraction<cFM,gMM, cLT>::iterate()
fatalExit;
}
boundaryContactSearchTimer_.start();
for(uint32 i=0; i<6u; i++)
{
auto& BI =boundaryInteraction_[i];
if(BI.ppPairsAllocated())
{
if( !contactSearch_().boundaryBroadSearch(
i,
iter,
t,
dt,
boundaryInteraction_[i].ppPairs(),
boundaryInteraction_[i].pwPairs()))
BI.ppPairs(),
BI.pwPairs()))
{
fatalErrorInFunction<<
"failed to perform broadSearch for boundary index "<<i<<endl;
return false;
}
}
}
boundaryContactSearchTimer_.end();
if(broadSearch && contactSearch_().performedBroadSearch())
{
@ -224,12 +235,44 @@ bool pFlow::sphereInteraction<cFM,gMM, cLT>::iterate()
contactListMangementTimer_.end();
}
contactListBoundaryTimer_.resume();
for(uint32 i=0; i<6u; i++)
{
boundaryInteraction_[i].ppPairs().afterBroadSearch();
boundaryInteraction_[i].pwPairs().afterBroadSearch();
auto& BI = boundaryInteraction_[i];
if(BI.ppPairsAllocated()) BI.ppPairs().afterBroadSearch();
if(BI.pwPairsAllocated()) BI.pwPairs().afterBroadSearch();
}
contactListBoundaryTimer_.end();
{
boundaryInteractionTimer_.start();
std::array<bool,6> requireStep{
boundaryInteraction_[0].isBoundaryMaster(),
boundaryInteraction_[1].isBoundaryMaster(),
boundaryInteraction_[2].isBoundaryMaster(),
boundaryInteraction_[3].isBoundaryMaster(),
boundaryInteraction_[4].isBoundaryMaster(),
boundaryInteraction_[5].isBoundaryMaster()};
int step = 1;
const auto& cfModel = this->forceModel_();
while( std::any_of(requireStep.begin(), requireStep.end(), [](bool v) { return v==true; }))
{
for(uint32 i=0; i<6u; i++)
{
if(step==1u || requireStep[i] )
{
requireStep[i] = boundaryInteraction_[i].sphereSphereInteraction(
dt,
this->forceModel_(),
step
);
}
}
step++;
}
boundaryInteractionTimer_.pause();
}
ppInteractionTimer_.start();
sphereSphereInteraction();
ppInteractionTimer_.end();
@ -239,14 +282,36 @@ bool pFlow::sphereInteraction<cFM,gMM, cLT>::iterate()
sphereWallInteraction();
pwInteractionTimer_.end();
boundaryInteractionTimer_.start();
{
boundaryInteractionTimer_.resume();
std::array<bool,6> requireStep{
!boundaryInteraction_[0].isBoundaryMaster(),
!boundaryInteraction_[1].isBoundaryMaster(),
!boundaryInteraction_[2].isBoundaryMaster(),
!boundaryInteraction_[3].isBoundaryMaster(),
!boundaryInteraction_[4].isBoundaryMaster(),
!boundaryInteraction_[5].isBoundaryMaster()};
int step = 2;
const auto& cfModel = this->forceModel_();
while(std::any_of(requireStep.begin(), requireStep.end(), [](bool v) { return v==true; }))
{
for(uint32 i=0; i<6u; i++)
{
boundaryInteraction_[i].sphereSphereInteraction(
if(requireStep[i])
{
requireStep[i] = boundaryInteraction_[i].sphereSphereInteraction(
dt,
this->forceModel_());
this->forceModel_(),
step
);
}
}
step++;
}
boundaryInteractionTimer_.end();
}
return true;
}

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

@ -96,9 +96,12 @@ private:
/// timer for managing contact lists (only inernal points)
Timer contactListMangementTimer_;
Timer boundaryContactSearchTimer_;
/// timer for boundary interaction time
Timer boundaryInteractionTimer_;
Timer contactListBoundaryTimer_;
bool createSphereInteraction();

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

@ -261,7 +261,7 @@ struct pwInteractionFunctor
int32 propId_i = propId_[i];
int32 wPropId_j = wPropId_[tj];
realx3 FCn, FCt, Mri, Mrj, Mij, Mji;
realx3 FCn, FCt, Mri, Mrj, Mij;
//output<< "before "<<history.overlap_t_<<endl;
// calculates contact force
forceModel_.contactForce(

2
src/MotionModel/entities/stationary/stationary.hpp
View File

@ -62,7 +62,7 @@ public:
INLINE_FUNCTION_HD
realx3 linVelocityPoint(const realx3 &)const
{
return zero3;
return realx3(0);
}
INLINE_FUNCTION_HD

6
src/MotionModel/rotatingAxisMotion/rotatingAxisMotion.hpp
View File

@ -77,8 +77,6 @@ protected:
return true;
}
void impl_setTime(uint32 iter, real t, real dt)const;
public:
TypeInfo("rotatingAxisMotion");
@ -90,6 +88,7 @@ public:
const dictionary& dict,
repository* owner);
using fileDictionary::write;
bool write(iOstream& os, const IOPattern& iop)const override;
@ -98,6 +97,9 @@ public:
{
return rotatingAxis({0,0,0}, {1,0,0}, 0.0);
}
// TODO: make this method protected
void impl_setTime(uint32 iter, real t, real dt)const;
};
} // pFlow

2
src/MotionModel/stationaryWall/stationaryWall.hpp
View File

@ -84,6 +84,8 @@ public:
repository* owner);
using fileDictionary::write;
bool write(iOstream& os, const IOPattern& iop)const override;
static

6
src/MotionModel/vibratingMotion/vibratingMotion.hpp
View File

@ -82,7 +82,6 @@ protected:
return true;
}
void impl_setTime(uint32 iter, real t, real dt)const;
public:
@ -99,15 +98,20 @@ public:
/// Destructor
~vibratingMotion()override = default;
using fileDictionary::write;
bool write(iOstream& os, const IOPattern& iop)const override;
static
auto noneComponent()
{
return vibrating();
}
// TODO: make this protected
void impl_setTime(uint32 iter, real t, real dt)const;
};
} // pFlow

4
src/Particles/Insertion/collisionCheck/collisionCheck.cpp
View File

@ -60,7 +60,7 @@ pFlow::collisionCheck::checkPoint(const realx3& p, const real d) const
{
uint32 n = head_(i, j, k);
while( n != -1)
while( n != static_cast<uint32>(-1))
{
if( ((position_[n]-p).length() - 0.5*(diameters_[n]+d )) <= 0.0 )
{
@ -85,7 +85,7 @@ pFlow::collisionCheck::mapLastAddedParticle()
"size mismatch of next and position"<<endl;
return false;
}
next_.push_back(-1);
next_.push_back(static_cast<uint32>(-1));
const auto& p = position_[n];
if(!searchBox_.isInside(p))

3
src/Particles/Insertion/insertion/insertion.cpp
View File

@ -74,9 +74,6 @@ pFlow::insertion::pStruct() const
return particles_.pStruct();
}
bool
pFlow::insertion::readInsertionDict()
{

6
src/Particles/Insertion/insertion/insertion.hpp
View File

@ -32,7 +32,9 @@ class pointStructure;
/**
* Base class for particle insertion
*/
class insertion : public fileDictionary
class insertion
:
public fileDictionary
{
private:
@ -118,6 +120,8 @@ public:
/*/// read from iIstream
virtual bool read(iIstream& is) = 0;*/
using fileDictionary::write;
/// write to iOstream
bool write(iOstream& os, const IOPattern& iop)const override ;
};

7
src/Particles/SphereParticles/sphereParticles/sphereParticles.cpp
View File

@ -389,7 +389,9 @@ pFlow::sphereParticles::sphereParticles(
intPredictTimer_(
"Integration-predict", &this->timers() ),
intCorrectTimer_(
"Integration-correct", &this->timers() )
"Integration-correct", &this->timers() ),
fieldUpdateTimer_(
"fieldUpdate", &this->timers() )
{
auto intMethod = control.settingsDict().getVal<word>("integrationMethod");
@ -514,13 +516,14 @@ bool pFlow::sphereParticles::beforeIteration()
rVelIntegration_().predict(dt,rVelocity_, rAcceleration_);
intPredictTimer_.end();
fieldUpdateTimer_.start();
propertyId_.updateBoundariesSlaveToMasterIfRequested();
diameter_.updateBoundariesSlaveToMasterIfRequested();
mass_.updateBoundariesSlaveToMasterIfRequested();
I_.updateBoundariesSlaveToMasterIfRequested();
rVelocity_.updateBoundariesSlaveToMasterIfRequested();
rAcceleration_.updateBoundariesSlaveToMasterIfRequested();
fieldUpdateTimer_.end();
return true;
}

7
src/Particles/SphereParticles/sphereParticles/sphereParticles.hpp
View File

@ -79,10 +79,10 @@ private:
/// timer for integration computations (correction step)
Timer intCorrectTimer_;
Timer fieldUpdateTimer_;
private:
bool initializeParticles();
bool getParticlesInfoFromShape(
const wordVector& shapeNames,
@ -124,6 +124,9 @@ public:
const setFieldList& setField
) override ;*/
// TODO: make this method private later
bool initializeParticles();
/// const reference to shapes object
const auto& spheres() const
{

11
src/Particles/dynamicPointStructure/dynamicPointStructure.cpp
View File

@ -38,6 +38,7 @@ pFlow::dynamicPointStructure::dynamicPointStructure
*this,
zero3
),
velocityUpdateTimer_("velocity boundary update", &timers()),
integrationMethod_
(
control.settingsDict().getVal<word>("integrationMethod")
@ -81,11 +82,11 @@ pFlow::dynamicPointStructure::dynamicPointStructure
bool pFlow::dynamicPointStructure::beforeIteration()
{
return pointStructure::beforeIteration();
/*real dt = this->dt();
auto& acc = time().lookupObject<realx3PointField_D>("acceleration");
return predict(dt, acc);*/
if(!pointStructure::beforeIteration())return false;
velocityUpdateTimer_.start();
velocity_.updateBoundariesSlaveToMasterIfRequested();
velocityUpdateTimer_.end();
return true;
}
bool pFlow::dynamicPointStructure::iterate()

2
src/Particles/dynamicPointStructure/dynamicPointStructure.hpp
View File

@ -44,6 +44,8 @@ private:
uniquePtr<integration> integrationVel_ = nullptr;
Timer velocityUpdateTimer_;
/// @brief integration method for velocity and position
word integrationMethod_;

11
src/Particles/particles/particles.cpp
View File

@ -74,7 +74,8 @@ pFlow::particles::particles(systemControl& control)
dynPointStruct_,
zero3
),
idHandler_(particleIdHandler::create(dynPointStruct_))
idHandler_(particleIdHandler::create(dynPointStruct_)),
baseFieldBoundaryUpdateTimer_("baseFieldBoundaryUpdate",&timers())
{
this->addToSubscriber(dynPointStruct_);
@ -84,18 +85,18 @@ pFlow::particles::particles(systemControl& control)
bool
pFlow::particles::beforeIteration()
{
zeroForce();
zeroTorque();
if( !dynPointStruct_.beforeIteration())
{
return false;
}
zeroForce();
zeroTorque();
baseFieldBoundaryUpdateTimer_.start();
shapeIndex_.updateBoundariesSlaveToMasterIfRequested();
accelertion_.updateBoundariesSlaveToMasterIfRequested();
idHandler_().updateBoundariesSlaveToMasterIfRequested();
baseFieldBoundaryUpdateTimer_.end();
return true;
}

2
src/Particles/particles/particles.hpp
View File

@ -54,6 +54,8 @@ private:
/// handling new ids for new particles
uniquePtr<particleIdHandler> idHandler_ = nullptr;
Timer baseFieldBoundaryUpdateTimer_;
/// messages for this objects
static inline const message defaultMessage_{ message::DEFAULT };

7
src/Particles/particles/regularParticleIdHandler/regularParticleIdHandler.cpp
View File

@ -15,8 +15,11 @@ pFlow::regularParticleIdHandler::regularParticleIdHandler
pFlow::Pair<pFlow::uint32, pFlow::uint32>
pFlow::regularParticleIdHandler::getIdRange(uint32 nNewParticles)
{
if(nNewParticles==0) return {0,0};
uint32 startId;
if(maxId_==-1)
if(maxId_== static_cast<uint32>(-1))
{
startId = 0;
}
@ -37,7 +40,7 @@ bool pFlow::regularParticleIdHandler::initialIdCheck()
uint32 maxId = max( *this );
/// particles should get ids from 0 to size-1
if(maxId == -1)
if(maxId == static_cast<uint32>(-1))
{
fillSequence(*this,0u);
maxId_ = size()-1;

2
src/Particles/particles/regularParticleIdHandler/regularParticleIdHandler.hpp
View File

@ -31,7 +31,7 @@ class regularParticleIdHandler
{
private:
uint32 maxId_ = -1;
uint32 maxId_ = static_cast<uint32>(-1);
bool initialIdCheck()override;
public:

4
src/Particles/particles/shape/baseShapeNames.hpp
View File

@ -132,7 +132,7 @@ public:
return true;
}
}
idx = -1;
idx = static_cast<uint32>(-1);
return false;
}
@ -144,6 +144,8 @@ public:
// - IO
using fileDictionary::write;
bool write(iOstream& os)const override;
};

3
src/phasicFlow/CMakeLists.txt
View File

@ -59,6 +59,7 @@ Timer/Timers.cpp
containers/Vector/Vectors.cpp
containers/VectorHD/wordVectorHost.cpp
containers/VectorHD/VectorSingles.cpp
containers/Field/Fields.cpp
containers/symArrayHD/symArrays.cpp
@ -94,6 +95,8 @@ structuredData/pointStructure/selectors/selectorStridedRange/selectorStridedRang
structuredData/pointStructure/selectors/selectorRandomPoints/selectorRandomPoints.cpp
structuredData/pointStructure/selectors/selectBox/selectBox.cpp
structuredData/pointStructure/selectors/selectorGeometric/selectorGeometrics.cpp
structuredData/pointStructure/pointStructure/pointSorting/mortonIndexing.cpp
structuredData/pointStructure/pointStructure/pointSorting/pointSorting.cpp
triSurface/subSurface.cpp
triSurface/triSurface.cpp

4
src/phasicFlow/Kokkos/ViewAlgorithms.hpp
View File

@ -421,10 +421,10 @@ binarySearch(
)
{
if (end <= start)
return -1;
return static_cast<uint32>(-1);
if (auto res = binarySearch_(view.data() + start, end - start, val);
res != -1)
res != static_cast<uint32>(-1))
{
return res + start;
}

4
src/phasicFlow/commandLine/CLI/TypeTools.hpp
View File

@ -146,11 +146,11 @@ template <typename T, typename C> class is_direct_constructible {
static auto test(int, std::true_type) -> decltype(
// NVCC warns about narrowing conversions here
#ifdef __CUDACC__
#pragma diag_suppress 2361
#pragma nv_diag_suppress 2361
#endif
TT { std::declval<CC>() }
#ifdef __CUDACC__
#pragma diag_default 2361
#pragma nv_diag_default 2361
#endif
,
std::is_move_assignable<TT>());

1
src/phasicFlow/containers/Field/Field.hpp
View File

@ -24,6 +24,7 @@ Licence:
#include "types.hpp"
#include "VectorSingle.hpp"
#include "wordVectorHost.hpp"
#include "Vector.hpp"
#include "streams.hpp"

3
src/phasicFlow/containers/Field/Fields.cpp
View File

@ -31,5 +31,4 @@ template class pFlow::Field<pFlow::real>;
template class pFlow::Field<pFlow::realx3>;
template class pFlow::Field<pFlow::word, pFlow::HostSpace>;

23
src/phasicFlow/containers/List/ListPtr/ListPtrI.hpp
View File

@ -39,6 +39,7 @@ inline bool pFlow::ListPtr<T>::copy(const ListPtrType& src)
}
template<typename T>
inline
T* pFlow::ListPtr<T>::ptr(size_t i)
{
@ -51,6 +52,7 @@ T* pFlow::ListPtr<T>::ptr(size_t i)
}
template<typename T>
inline
const T* pFlow::ListPtr<T>::ptr
(
size_t i
@ -66,6 +68,7 @@ const T* pFlow::ListPtr<T>::ptr
}
template<typename T>
inline
auto pFlow::ListPtr<T>::pos
(
size_t i
@ -84,6 +87,7 @@ auto pFlow::ListPtr<T>::pos
}
template<typename T>
inline
auto pFlow::ListPtr<T>::pos
(
size_t i
@ -102,6 +106,7 @@ auto pFlow::ListPtr<T>::pos
}
template<typename T>
inline
pFlow::ListPtr<T>::ListPtr
(
const ListPtrType& src
@ -119,6 +124,7 @@ pFlow::ListPtr<T>::ListPtr
template<typename T>
inline
pFlow::ListPtr<T>& pFlow::ListPtr<T>::operator=
(
const ListPtrType& rhs
@ -144,6 +150,7 @@ pFlow::ListPtr<T>& pFlow::ListPtr<T>::operator=
}
template<typename T>
inline
pFlow::uniquePtr<T> pFlow::ListPtr<T>::set
(
size_t i, T* ptr
@ -155,6 +162,7 @@ pFlow::uniquePtr<T> pFlow::ListPtr<T>::set
template<typename T>
inline
pFlow::uniquePtr<T> pFlow::ListPtr<T>::set
(
size_t i,
@ -179,6 +187,7 @@ pFlow::uniquePtr<T> pFlow::ListPtr<T>::set
template<typename T>
template<typename... Args>
inline
pFlow::uniquePtr<T> pFlow::ListPtr<T>::setSafe
(
size_t i,
@ -190,6 +199,7 @@ pFlow::uniquePtr<T> pFlow::ListPtr<T>::setSafe
}
template<typename T>
inline
void pFlow::ListPtr<T>::push_back
(
T* ptr
@ -199,6 +209,7 @@ void pFlow::ListPtr<T>::push_back
}
template<typename T>
inline
void pFlow::ListPtr<T>::push_back(uniquePtr<T>&& ptr)
{
list_.push_back( ptr.release() );
@ -206,13 +217,15 @@ void pFlow::ListPtr<T>::push_back(uniquePtr<T>&& ptr)
template<typename T>
template<typename... Args>
inline
void pFlow::ListPtr<T>::push_backSafe(Args&&... args)
{
auto ptr=makeUnique<T>(std::forward<Args>(args)...) ;
uniquePtr<T> ptr = makeUnique<T>(std::forward<Args>(args)...) ;
push_back(ptr);
}
template<typename T>
inline
T& pFlow::ListPtr<T>::operator[]
(
size_t i
@ -231,6 +244,7 @@ T& pFlow::ListPtr<T>::operator[]
}
template<typename T>
inline
const T& pFlow::ListPtr<T>::operator[]
(
size_t i
@ -248,18 +262,21 @@ const T& pFlow::ListPtr<T>::operator[]
}
template<typename T>
inline
size_t pFlow::ListPtr<T>::size()const
{
return list_.size();
}
template<typename T>
inline
auto pFlow::ListPtr<T>::empty() const
{
return list_.emtpy();
}
template<typename T>
inline
pFlow::uniquePtr<T> pFlow::ListPtr<T>::release
(
size_t i
@ -273,15 +290,14 @@ pFlow::uniquePtr<T> pFlow::ListPtr<T>::release
template<typename T>
inline
void pFlow::ListPtr<T>::clear()
{
int i =0;
for( auto iter = list_.begin(); iter != list_.end(); ++iter )
{
if(*iter != nullptr)
{
delete *iter;
*iter = nullptr;
}
@ -291,6 +307,7 @@ void pFlow::ListPtr<T>::clear()
}
template<typename T>
inline
void pFlow::ListPtr<T>::clear
(
size_t i

195
src/phasicFlow/containers/VectorHD/VectorSingle.cpp
View File

@ -55,8 +55,7 @@ void pFlow::VectorSingle<T,MemorySpace>::changeCapacity
bool withInit
)
{
if constexpr( isTriviallyCopyable_ )
{
if(withInit)
{
resizeInit(view_, actualCap);
@ -65,21 +64,6 @@ void pFlow::VectorSingle<T,MemorySpace>::changeCapacity
{
resizeNoInit(view_, actualCap);
}
}
else if constexpr( isHostAccessible_ )
{
viewType newView(view_.label(), actualCap);
for(auto i=0u; i<min(size_,actualCap); i++)
{
newView(i) = view_(i);
}
view_ = newView;
}
else
{
static_assert("changeCapacity is not a valid operation for non-trivially-copyable data type on device memory");
}
}
@ -91,16 +75,7 @@ pFlow::uint32 pFlow::VectorSingle<T,MemorySpace>::reallocateCapacitySize
uint32 s
)
{
if constexpr (isTriviallyCopyable_)
{
reallocNoInit(view_, cap);
}
else
{
viewType newView(view_.label(), cap);
view_ = newView;
}
return setSize(s);
}
@ -209,21 +184,19 @@ pFlow::VectorSingle<T,MemorySpace>::VectorSingle
:
VectorSingle(name, src.capacity(), src.size(), RESERVE())
{
if constexpr(isTriviallyCopyable_)
{
copy(deviceView(), src.deviceView());
}
else if constexpr( isHostAccessible_)
{
for(auto i=0u; i<src.size(); i++)
{
view_[i] = src.view_[i];
}
}
else
{
static_assert("This constructor is not valid for non-trivially copyable data type on device memory");
}
}
template<typename T, typename MemorySpace>
pFlow::VectorSingle<T,MemorySpace>::VectorSingle
(
const word& name,
const ViewType1D<T, MemorySpace>& src
)
:
VectorSingle(name, src.size(), src.size(), RESERVE())
{
copy(deviceView(), src);
}
template<typename T, typename MemorySpace>
@ -299,18 +272,7 @@ INLINE_FUNCTION_H
auto pFlow::VectorSingle<T,MemorySpace>::hostViewAll()const
{
auto hView = Kokkos::create_mirror_view(view_);
if constexpr(isTriviallyCopyable_)
{
copy(hView, view_);
}
else if constexpr( isHostAccessible_ )
{
// nothing to be done, since it is already a host memory
}
else
{
static_assert("hostViewAll is not valid for non-trivially copyable data type on device memory");
}
return hView;
}
@ -319,19 +281,7 @@ INLINE_FUNCTION_H
auto pFlow::VectorSingle<T,MemorySpace>::hostView()const
{
auto hView = Kokkos::create_mirror_view(deviceView());
if constexpr(isTriviallyCopyable_)
{
copy(hView, deviceView());
}
else if constexpr( isHostAccessible_ )
{
// nothing to be done, since it is already a host memory
}
else
{
static_assert("hostView is not valid for non-trivially copyable data type on device memory");
}
return hView;
}
@ -367,6 +317,7 @@ template<typename T, typename MemorySpace>
INLINE_FUNCTION_H
void pFlow::VectorSingle<T,MemorySpace>::reserve(uint32 cap)
{
if(cap == capacity() ) return;
changeCapacity(cap);
}
@ -470,23 +421,9 @@ void pFlow::VectorSingle<T,MemorySpace>::assign
changeSize(srcSize);
}
if constexpr( isTriviallyCopyable_ )
{
// - unmanaged view in the host
hostViewType1D<const T> temp(src.data(), srcSize );
copy(deviceView(), temp);
}
else if constexpr( isHostAccessible_)
{
for(auto i=0u; i<srcSize; i++)
{
view_[i] = src[i];
}
}
else
{
static_assert("Not a valid operation for this data type on memory device");
}
}
@ -516,21 +453,8 @@ void pFlow::VectorSingle<T,MemorySpace>::assignFromHost(const VectorTypeHost& sr
changeSize(srcSize);
}
if constexpr(isTriviallyCopyable_)
{
copy(deviceView(), src.hostView());
}
else if constexpr( isHostAccessible_)
{
for(auto i=0u; i<src.size(); i++)
{
view_[i] = src.view_[i];
}
}
else
{
static_assert("Not a valid operation for this data type on device memory");
}
}
template<typename T, typename MemorySpace>
@ -553,22 +477,48 @@ void pFlow::VectorSingle<T,MemorySpace>::assign
changeSize(srcSize);
}
if constexpr(isTriviallyCopyable_)
{
copy(deviceView(), src.deviceView());
}
template<typename T, typename MemorySpace>
template<typename MSpace>
INLINE_FUNCTION_H
void pFlow::VectorSingle<T,MemorySpace>::assignFromDevice(
const VectorSingle<T, MSpace>& src,
bool srcCapacity
)
{
uint32 srcSize = src.size();
uint32 srcCap = src.capacity();
if(srcCapacity && srcCap != capacity()){
reallocateCapacitySize(srcCap, srcSize);
}
else if constexpr( isHostAccessible_)
{
for(auto i=0u; i<src.size(); i++)
{
view_[i] = src.view_[i];
}
}
else
{
static_assert("Not a valid operation for this data type on device memory");
else {
changeSize(srcSize);
}
copy(deviceView(), src.deviceView());
}
template <typename T, typename MemorySpace>
INLINE_FUNCTION_H
void pFlow::VectorSingle<T, MemorySpace>::append(const ViewType1D<T,MemorySpace>& appVec)
{
uint32 appSize = appVec.size();
if(appSize == 0) return;
uint32 oldS = size();
uint32 newSize = oldS + appSize;
changeSize(newSize);
auto appendView = Kokkos::subview(
view_,
Kokkos::make_pair<uint32>(oldS, newSize));
copy(appendView, appVec);
}
template <typename T, typename MemorySpace>
@ -589,22 +539,7 @@ void pFlow::VectorSingle<T, MemorySpace>::append
hostViewType1D<const T> temp(appVec.data(), srcSize );
auto dest = Kokkos::subview(view_, Kokkos::make_pair<uint32>(oldSize,newSize));
if constexpr( isTriviallyCopyable_)
{
copy(dest, temp);
}
else if constexpr( isHostAccessible_)
{
for(auto i=0u; i<appVec.size(); i++)
{
dest[i] = appVec[i];
}
}
else
{
static_assert("not a valid operation for this data type on device memory");
}
}
template <typename T, typename MemorySpace>
@ -626,21 +561,8 @@ void pFlow::VectorSingle<T, MemorySpace>::append
view_,
Kokkos::make_pair<uint32>(oldS, newSize));
if constexpr( isTriviallyCopyable_)
{
copy(appendView, appVec.deviceView());
}
else if constexpr( isHostAccessible_)
{
for(auto i=0u; i<appVec.size(); i++)
{
appendView[i] = appVec.view_[i];
}
}
else
{
static_assert("not a valid operation for this data type on device memory");
}
}
template <typename T, typename MemorySpace>
@ -813,7 +735,7 @@ bool pFlow::VectorSingle<T,MemorySpace>::insertSetElement
template<typename T, typename MemorySpace>
INLINE_FUNCTION_H
bool pFlow::VectorSingle<T,MemorySpace>::reorderItems(uint32IndexContainer indices)
bool pFlow::VectorSingle<T,MemorySpace>::reorderItems(const uint32IndexContainer& indices)
{
if(indices.size() == 0)
{
@ -833,8 +755,6 @@ bool pFlow::VectorSingle<T,MemorySpace>::reorderItems(uint32IndexContainer indic
uint32 newSize = indices.size();
setSize(newSize);
viewType sortedView(this->name(), newSize);
using policy = Kokkos::RangePolicy< execution_space,Kokkos::IndexType<uint32>>;
@ -875,7 +795,10 @@ bool pFlow::VectorSingle<T,MemorySpace>::reorderItems(uint32IndexContainer indic
Kokkos::fence();
}
setSize(newSize);
copy(deviceView(), sortedView);
return true;
}

31
src/phasicFlow/containers/VectorHD/VectorSingle.hpp
View File

@ -39,9 +39,6 @@ Licence:
namespace pFlow
{
//- Forward
template<typename T, typename MemorySpace>
class VectorSingle;
template<typename T, typename MemorySpace=void>
class VectorSingle
@ -101,6 +98,8 @@ private:
static constexpr
bool isTriviallyCopyable_ = std::is_trivially_copyable_v<T>;
static_assert(isTriviallyCopyable_, "This type is not trivially copyable");
/// Evaluate capacity based on the input size
static INLINE_FUNCTION_H uint32 evalCapacity(uint32 n)
{
@ -159,6 +158,9 @@ public:
/// Copy construct with a new name (perform deep copy)
VectorSingle(const word& name, const VectorSingle& src);
/// Copy construct with a new name (perform deep copy)
VectorSingle(const word& name, const ViewType1D<T, MemorySpace>& src);
/// Copy assignment (perform deep copy from rhs to *this)
VectorSingle& operator = (const VectorSingle& rhs) ;
@ -287,25 +289,10 @@ public:
template<typename MSpace>
INLINE_FUNCTION_H
void assignFromDevice(const VectorSingle<T, MSpace>& src, bool srcCapacity = true)
{
uint32 srcSize = src.size();
uint32 srcCap = src.capacity();
void assignFromDevice(const VectorSingle<T, MSpace>& src, bool srcCapacity = true);
if(srcCapacity && srcCap != capacity()){
reallocateCapacitySize(srcCap, srcSize);
}
else {
changeSize(srcSize);
}
if constexpr(isTriviallyCopyable_){
copy(deviceView(), src.deviceView());
}
else{
static_assert("Not a valid operation for this data type ");
}
}
INLINE_FUNCTION_H
void append(const ViewType1D<T,MemorySpace>& appVec);
INLINE_FUNCTION_H
void append(const std::vector<T>& appVec);
@ -331,7 +318,7 @@ public:
const ViewType1D<T, memory_space> vals);
INLINE_FUNCTION_H
bool reorderItems(uint32IndexContainer indices);
bool reorderItems(const uint32IndexContainer& indices);
/// @brief push a new element at the end (host call only)
/// resize if necessary and works on host accessible vector.

3
src/phasicFlow/containers/VectorHD/VectorSingles.hpp
View File

@ -25,6 +25,7 @@ Licence:
#include "types.hpp"
#include "VectorSingle.hpp"
#include "wordVectorHost.hpp"
namespace pFlow
{
@ -77,6 +78,8 @@ typedef VectorSingle<realx3x3> realx3x3Vector_D;
typedef VectorSingle<realx3x3, HostSpace> realx3x3Vector_H;
typedef VectorSingle<word, HostSpace> wordVector_H;
}
#endif

23
src/phasicFlow/containers/VectorHD/wordVectorHost.cpp Normal file
View File

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

558
src/phasicFlow/containers/VectorHD/wordVectorHost.hpp Normal file
View File

@ -0,0 +1,558 @@
/*------------------------------- 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 __wordVectorHost_hpp__
#define __wordVectorHost_hpp__
#include "VectorSingle.hpp"
#include "Vector.hpp"
namespace pFlow
{
template<>
class VectorSingle<word, HostSpace>
{
public:
//- typedefs for accessing data
using VectorType = VectorSingle<word, HostSpace>;
using VectorTypeHost = VectorSingle<word, HostSpace>;
using iterator = word*;
using const_iterator = const word*;
using reference = word&;
using const_reference = const word&;
using value_type = word;
using pointer = word*;
using const_pointer = const word*;
//- typedefs related to memory management
using viewType = ViewType1D<word, HostSpace>;
using device_type = typename viewType::device_type;
using memory_space = typename viewType::memory_space;
using execution_space = typename viewType::execution_space;
private:
// - Data members
Vector<word> container_;
mutable viewType unManagedView_;
// - protected members and methods
/// Is the memory of this vector accessible from Host
static constexpr
bool isHostAccessible_ = true;
/// Is the memory of this vector accessiple from Divce
static constexpr
bool isDeviceAccessible_ = false;
static constexpr
bool isTriviallyCopyable_ = false;
public:
/// Type info
TypeInfoTemplateNV111("VectorSingle", word , memoerySpaceName());
//// - Constructors
/// Empty vector
VectorSingle() = default;
/// Empty vector with a name (capacity = 2)
explicit VectorSingle(const word& name)
:
container_(name)
{}
/// Vector with name and size n
VectorSingle(const word& name, uint32 n)
:
container_(name, n)
{}
/// Vector with name, size and value
VectorSingle(const word& name, uint32 n, const word& val)
:
container_(name, n, val)
{}
/// Vector with name, size (n) and reserved capacity
VectorSingle(const word& name, uint32 cap, uint32 n, const RESERVE& r )
:
container_(name, cap, n, r)
{}
/// Construct with a name and form std::vector (host memory)
VectorSingle(const word& name, const std::vector<word> & src)
:
container_(name, src)
{}
/// Construct with a name and form std::vector (host memory) and with a desired capacity.
VectorSingle(const word& name, const std::vector<word> & src, uint32 cap)
:
container_(name, src, cap)
{}
/// Copy construct (performs deep copy)
VectorSingle(const VectorSingle& src)=default;
/// Copy construct with a new name (perform deep copy)
VectorSingle(const word& name, const VectorSingle& src)
:
container_(name, src.container_)
{}
/// Copy assignment (perform deep copy from rhs to *this)
VectorSingle& operator = (const VectorSingle& rhs)=default;
/// Move construct
VectorSingle(VectorSingle&&) = default;
/// Move assignment
VectorSingle& operator= (VectorSingle&&) = default;
/// @brief Descructor
/// This may not destroy the underlying memory, sice view is
/// shared_ptr and maybe referenced by another object too
~VectorSingle() = default;
/// Clone as a uniquePtr (perform deep copy)
INLINE_FUNCTION_H
uniquePtr<VectorSingle> clone() const
{
return makeUnique<VectorSingle>(*this);
}
//// - Methods
/// Return *this
INLINE_FUNCTION_H
VectorType& VectorField()
{
return *this;
}
/// Return *this
INLINE_FUNCTION_H
const VectorType& VectorField()const
{
return *this;
}
/// Device view range [0,capcity)
INLINE_FUNCTION_H
auto& deviceViewAll()
{
// return un-managed view
unManagedView_ = viewType(container_.data(), container_.capacity());
return unManagedView_;
}
/// Device view range [0,capcity)
INLINE_FUNCTION_H
const auto& deviceViewAll() const
{
// return un-managed view
unManagedView_ = viewType(const_cast<word*>(container_.data()), container_.capacity());
return unManagedView_;
}
/// Device view range [0, size)
INLINE_FUNCTION_H
auto deviceView()const
{
return viewType(const_cast<word*>(container_.data()), container_.size());
}
/// Return a view accessible on Host in range [0,capacity)
INLINE_FUNCTION_H
auto hostViewAll()const
{
return viewType(const_cast<word*>(container_.data()), container_.capacity());
}
/// Return a view accessible on Host in range [0,size)
INLINE_FUNCTION_H
auto hostView()const
{
return viewType(const_cast<word*>(container_.data()), container_.size());
}
/// Name of the vector
INLINE_FUNCTION_H
word name()const
{
return container_.name();
}
/// Size of the vector
INLINE_FUNCTION_H
uint32 size()const
{
return container_.size();
}
// Capcity of the vector
INLINE_FUNCTION_H
uint32 capacity()const
{
return container_.capacity();
}
/// If vector is empty
INLINE_FUNCTION_H
bool empty()const
{
return container_.size()==0uL;
}
/// Reserve capacity for vector
/// Preserve the content.
INLINE_FUNCTION_H
void reserve(uint32 cap)
{
container_.reserve(cap);
}
/// Reallocate memory to new cap and set size to 0.
/*INLINE_FUNCTION_H
void reallocate(uint32 cap);
/// Reallocate memory to new cap and set size to newSize.
/// Do not preserve the content
INLINE_FUNCTION_H
void reallocate(uint32 cap, uint32 newSize);*/
/// Resize the vector and preserve the content
INLINE_FUNCTION_H
void resize(uint32 n)
{
container_.resize(n);
}
/// Resize the vector and assign the value to it.
INLINE_FUNCTION_H
void resize(uint32 n, const word& val)
{
container_.resize(n, val);
}
/// Clear the vector, but keep the allocated memory unchanged
INLINE_FUNCTION_H
void clear()
{
container_.clear();
}
/// Fill the range [0,size) with val
INLINE_FUNCTION_H
void fill(const word& val)
{
container_.fill(val);
}
INLINE_FUNCTION_H
void fill(rangeU32 r, const word& val)
{
container_.fill(r.start(), r.end(), val);
}
/// Change size of the vector and assign val to vector and
INLINE_FUNCTION_H
void assign(size_t n, const word& val)
{
container_.assign(n, val);
}
/// Assign source vector with specified capacity.
/// The size of *this becomes the size of src.
INLINE_FUNCTION_H
void assign(const std::vector<word>& src, uint32 cap)
{
container_.reserve(cap);
this->assign(src);
}
/// Assign source vector.
/// The size of *this becomes the size of src.
/// The capacity of *this becomes the capacity of src.
INLINE_FUNCTION_H
void assign(const std::vector<word>& src)
{
container_.assign(src.begin(), src.end());
}
/// Assign source vector from host side.
/// The size of *this becomes the size of src.
/// The capacity of *this becomes the capacity of src.
INLINE_FUNCTION_H
void assignFromHost(const VectorTypeHost& src)
{
notImplementedFunction;
fatalExit;
}
INLINE_FUNCTION_H
void assign(const VectorType& src, bool srcCapacity = true)
{
notImplementedFunction;
fatalExit;
}
template<typename MSpace>
INLINE_FUNCTION_H
void assignFromDevice(const VectorSingle<word, MSpace>& src, bool srcCapacity = true)
{
notImplementedFunction;
fatalExit;
}
/*INLINE_FUNCTION_H
void append(const ViewType1D<T,MemorySpace>& appVec);
INLINE_FUNCTION_H
void append(const std::vector<T>& appVec);
INLINE_FUNCTION_H
void append(const VectorType& appVec);*/
INLINE_FUNCTION_H
auto getSpan()
{
return span<word>(container_.data(), container_.size());
}
INLINE_FUNCTION_H
auto getSpan()const
{
return span<word>(const_cast<word*>(container_.data()), container_.size());
}
INLINE_FUNCTION_H
bool insertSetElement(const uint32IndexContainer& indices, const word& val)
{
notImplementedFunction;
return false;
}
INLINE_FUNCTION_H
bool insertSetElement(const uint32IndexContainer& indices, const std::vector<word>& vals)
{
notImplementedFunction;
return false;
}
INLINE_FUNCTION_H
bool insertSetElement(
const uint32IndexContainer& indices,
const ViewType1D<word, memory_space> vals)
{
notImplementedFunction;
return false;
}
INLINE_FUNCTION_H
bool reorderItems(const uint32IndexContainer& indices)
{
notImplementedFunction;
return false;
}
/// @brief push a new element at the end (host call only)
/// resize if necessary and works on host accessible vector.
void push_back(const word& val)
{
container_.push_back(val);
}
INLINE_FUNCTION_H
pointer data(){
return container_.data();
}
INLINE_FUNCTION_H
const_pointer data()const{
return container_.data();
}
/// Return begin iterator. It works when devices is host accessible.
auto
begin(){
return container_.begin();
}
/// Return begin iterator. it works when host is accessible.
const auto
begin()const {
return container_.begin();
}
auto end(){
return container_.end();
}
/// Return end iterator. it works when host is accessible.
const auto end()const{
return container_.end();
}
/// Return reference to element i. it works when host is accessible.
word& operator[](size_t i){
return container_[i];
}
const word& operator[](size_t i)const{
return container_[i];
}
//// - IO operations
/// Read vector from stream
FUNCTION_H
bool read(iIstream& is)
{
return container_.read(is);
}
/// Read vector from stream
FUNCTION_H
bool read(iIstream& is, const IOPattern& iop)
{
return container_.read(is, iop);
}
/// Write the vector to os
FUNCTION_H
bool write(iOstream& os, const IOPattern& iop)const
{
return container_.write(os, iop);
}
FUNCTION_H
bool write(iOstream& os)const
{
return container_.write(os);
}
template<typename HostMask>
FUNCTION_H
bool write(iOstream& os, const IOPattern& iop, const HostMask& mask)const
{
notImplementedFunction;
return false;
}
/// Name of the memory space
static
constexpr const char* memoerySpaceName()
{
return memory_space::name();
}
}; // class wordVectorHost
inline iIstream& operator >> (iIstream & is, VectorSingle<word, HostSpace> & ivec )
{
if( !ivec.read(is) )
{
ioErrorInFile (is.name(), is.lineNumber());
fatalExit;
}
return is;
}
inline iOstream& operator << (iOstream& os, const VectorSingle<word, HostSpace>& ovec )
{
if( !ovec.write(os) )
{
ioErrorInFile(os.name(), os.lineNumber());
fatalExit;
}
return os;
}
} // - pFlow
#endif
/*
INLINE_FUNCTION_H
bool append(const deviceViewType1D<T>& dVec, size_t numElems)
{
if(numElems == 0 )return true;
auto oldSize = size_;
auto newSize = size_ + numElems;
if(this->empty() || newSize > capacity() )
{
resize(newSize);
}
else
{
size_ = size_+numElems;
}
auto dSubView = Kokkos::subview(view_, Kokkos::make_pair(oldSize, newSize));
copy(dSubView, dVec);
return true;
}
INLINE_FUNCTION_H
bool append(const VectorSingle& Vec)
{
return append(Vec.deviceView(), Vec.size());
}*/

20
src/phasicFlow/containers/indexContainer/indexContainer.hpp
View File

@ -143,6 +143,16 @@ public:
indexContainer(ind.data(), ind.size())
{}
indexContainer(const DeviceViewType& ind):
size_(ind.size()),
views_("indexContainer", size_)
{
copy(views_.h_view, ind);
copy(views_.d_view, ind);
min_ = pFlow::min(views_.d_view, 0, size_);
max_ = pFlow::max(views_.d_view, 0, size_);
}
/// Copy
indexContainer(const indexContainer&) = default;
@ -240,13 +250,13 @@ public:
void syncViews()
{
bool findMinMax = false;
if(views_.template need_sync<HostType>())
if(views_.template need_sync<DeviceType>())
{
Kokkos::deep_copy(views_.d_view, views_.h_view);
views_.clear_sync_state();
findMinMax = true;
}
else if(views_.template need_sync<DeviceType>())
else if(views_.template need_sync<HostType>())
{
Kokkos::deep_copy(views_.h_view, views_.d_view);
views_.clear_sync_state();
@ -260,6 +270,12 @@ public:
}
}
void syncViews(uint32 newSize)
{
size_ = newSize;
syncViews();
}
};

29
src/phasicFlow/containers/pointField/boundaryField/boundaryField/boundaryField.hpp
View File

@ -121,21 +121,45 @@ public:
return true;
}
inline
InternalFieldType& internal()
{
return internal_;
}
inline
const InternalFieldType& internal()const
{
return internal_;
}
FieldAccessType thisField()const
{
if constexpr(isDeviceAccessible<execution_space>())
return FieldAccessType(
this->size(),
this->indexList().deviceViewAll(),
internal_.deviceViewAll());
else
return FieldAccessType(
this->size(),
this->boundary().indexListHost().deviceViewAll(),
internal_.deviceViewAll());
}
FieldAccessType mirrorField()const
{
if constexpr(isDeviceAccessible<execution_space>())
return FieldAccessType(
this->mirrorBoundary().size(),
this->mirrorBoundary().indexList().deviceViewAll(),
internal_.deviceViewAll());
else
return FieldAccessType(
this->mirrorBoundary().size(),
this->mirrorBoundary().indexListHost().deviceViewAll(),
internal_.deviceViewAll());
}
virtual
@ -144,11 +168,6 @@ public:
virtual
const ProcVectorType& neighborProcField()const;
void fill(const std::any& val)override
{
return;
}
virtual
void fill(const T& val)
{

4
src/phasicFlow/containers/pointField/boundaryField/boundaryFieldList.hpp
View File

@ -43,7 +43,7 @@ protected:
const boundaryList& boundaries_;
uint32 slaveToMasterUpdateIter_ = -1;
uint32 slaveToMasterUpdateIter_ = static_cast<uint32>(-1);
public:
@ -95,7 +95,7 @@ public:
bool slaveToMasterUpdateRequested()const
{
return slaveToMasterUpdateIter_ != -1;
return slaveToMasterUpdateIter_ != static_cast<uint32>(-1);
}

4
src/phasicFlow/containers/pointField/boundaryField/generalBoundary/generalBoundary.hpp
View File

@ -148,10 +148,6 @@ public:
const Time& time()const;
virtual
void fill(const std::any& val)=0;
};
} // pFlow

32
src/phasicFlow/containers/pointField/internalField/internalField.cpp
View File

@ -52,9 +52,29 @@ bool pFlow::internalField<T, MemorySpace>::insert(const anyList& varList)
}
return true;
}
template<class T, class MemorySpace>
bool pFlow::internalField<T, MemorySpace>::rearrange(const anyList& varList)
{
const word eventName = message::eventName(message::ITEM_REARRANGE);
const auto& indices = varList.getObject<uint32IndexContainer>(
eventName);
field_.reserve( internalPoints_.capacity());
field_.resize(internalPoints_.size());
if(!field_.reorderItems(indices))
{
fatalErrorInFunction<<
"cannot reorder items in field "<< name()<<endl;
return false;
}
return true;
}
template<class T, class MemorySpace>
pFlow::internalField<T, MemorySpace>::internalField
(
@ -166,15 +186,15 @@ bool pFlow::internalField<T, MemorySpace>:: hearChanges
{
// do nothing
}
if(msg.equivalentTo(message::ITEM_REARRANGE))
{
notImplementedFunction;
return false;
}
if(msg.equivalentTo(message::ITEM_INSERT))
{
return insert(varList);
}
if(msg.equivalentTo(message::ITEM_REARRANGE))
{
return rearrange(varList);
}
return true;
}

2
src/phasicFlow/containers/pointField/internalField/internalField.hpp
View File

@ -73,6 +73,8 @@ protected:
bool insert(const anyList& varList);
bool rearrange(const anyList& varList);
public:
internalField(

4
src/phasicFlow/dictionary/dictionary.hpp
View File

@ -302,10 +302,10 @@ public:
//// - IO operations
/// read from stream
virtual bool read(iIstream& is);
bool read(iIstream& is) override;
/// write to stream
virtual bool write(iOstream& os) const;
bool write(iOstream& os) const override;
};

3
src/phasicFlow/dictionary/fileDictionary.hpp
View File

@ -46,6 +46,9 @@ public:
const dictionary& dict,
repository* owner=nullptr);
using dictionary::read;
using dictionary::write;
/// read from stream
bool read(iIstream& is, const IOPattern& iop) override;

12
src/phasicFlow/eventManagement/message.hpp
View File

@ -48,14 +48,16 @@ public:
BNDR_RESET = 10, // boundary indices reset entirely
BNDR_DELETE = 11, // boundary indices deleted
BNDR_APPEND = 12, //
BNDR_PROCTRANS1 = 13, // transfer of data between processors step 1
BNDR_PROCTRANS2 = 14 // transfer of data between processors step 2
BNDR_PROCTRANSFER_SEND = 13, // transfer of data between processors step 1
BNDR_PROCTRANSFER_RECIEVE = 14, // transfer of data between processors step 2
BNDR_PROCTRANSFER_WAITFILL = 15, // wait for data and fill the field
BNDR_PROC_SIZE_CHANGED = 16
};
private:
static constexpr size_t numberOfEvents_ = 15;
static constexpr size_t numberOfEvents_ = 17;
std::bitset<numberOfEvents_> events_{0x0000};
@ -76,7 +78,9 @@ private:
"deletedIndices",
"appendedIndices",
"transferredIndices",
"numberRecieved"
"numToRecieve",
"insertedIndices",
"size"
};
public:

4
src/phasicFlow/globals/pFlowMacros.hpp
View File

@ -37,6 +37,10 @@ Licence:
#define CONSUME_PARAM(x) (void)(x);
#if defined(pFlow_Build_Cuda) && !defined(__CUDACC__)
#define __CUDACC__
#endif
#ifdef __CUDACC__
#define INLINE_FUNCTION_HD inline __host__ __device__
#define INLINE_FUNCTION_D inline __device__

9
src/phasicFlow/processors/localProcessors.hpp
View File

@ -145,6 +145,15 @@ public:
}
#endif
static
bool builtForMPI()
{
#ifdef pFlow_Build_MPI
return true;
#else
return false;
#endif
}
};
}

1
src/phasicFlow/processors/processors.cpp
View File

@ -29,7 +29,6 @@ Licence:
#include "processors.hpp"
#include "streams.hpp"
static int numVarsInitialized__ = 0;
#ifdef pFlow_Build_MPI

22
src/phasicFlow/repository/IOobject/IOfileHeader.cpp
View File

@ -42,12 +42,22 @@ pFlow::uniquePtr<pFlow::oFstream> pFlow::IOfileHeader::outStream()const
return osPtr;
}
pFlow::IOfileHeader::IOfileHeader
(
const objectFile& objf
)
:
objectFile(objf)
pFlow::uniquePtr<pFlow::oFstream> pFlow::IOfileHeader::dummyOutStream() const
{
auto osPtr = makeUnique<oFstream>( CWD()+word("dummyFile") , outFileBinary());
if(osPtr && owner())
{
auto outPrecision = owner()->outFilePrecision();
osPtr->precision(static_cast<int>(outPrecision));
}
return osPtr;
}
pFlow::IOfileHeader::IOfileHeader(
const objectFile &objf)
: objectFile(objf)
{}
pFlow::fileSystem pFlow::IOfileHeader::path() const

2
src/phasicFlow/repository/IOobject/IOfileHeader.hpp
View File

@ -58,6 +58,8 @@ protected:
// - ouput file stream
uniquePtr<oFstream> outStream()const;
uniquePtr<oFstream> dummyOutStream()const;
public:
// with owner

19
src/phasicFlow/repository/IOobject/IOobject.cpp
View File

@ -122,7 +122,8 @@ bool pFlow::IOobject::writeObject() const
if(ioPattern().thisCallWrite())
{
if( ioPattern().thisProcWriteData())
{
if(auto ptrOS = outStream(); ptrOS )
{
return writeObject(ptrOS());
@ -133,6 +134,22 @@ bool pFlow::IOobject::writeObject() const
"error in opening file "<< path() <<endl;
return false;
}
}
else
{
if(auto ptrOS = dummyOutStream(); ptrOS )
{
return writeObject(ptrOS());
}
else
{
warningInFunction<<
"error in opening file "<< path() <<endl;
return false;
}
}
}

10
src/phasicFlow/repository/Time/baseTimeControl.cpp
View File

@ -65,6 +65,16 @@ pFlow::baseTimeControl::baseTimeControl
}
pFlow::baseTimeControl::baseTimeControl(int32 start, int32 end, int32 stride, const word &intervalPrefix)
:
isTimeStep_(true),
iRange_(start, end, std::max(stride,1)),
intervalPrefix_(
intervalPrefix.size()==0uL? word("interval"): intervalPrefix+"Interval"
)
{
}
bool pFlow::baseTimeControl::timeEvent(uint32 iter, real t, real dt) const
{
if(isTimeStep_)

7
src/phasicFlow/repository/Time/baseTimeControl.hpp
View File

@ -46,6 +46,13 @@ public:
real defStartTime = 0.0
);
baseTimeControl(
int32 start,
int32 end,
int32 stride,
const word& intervalPrefix = ""
);
inline bool isTimeStep() const
{
return isTimeStep_;

12
src/phasicFlow/repository/Time/timeControl.cpp
View File

@ -18,7 +18,7 @@ Licence:
-----------------------------------------------------------------------------*/
#include "math.hpp"
#include "timeControl.hpp"
#include "dictionary.hpp"
@ -137,14 +137,14 @@ pFlow::word pFlow::timeControl::timeName()const
bool pFlow::timeControl::finalTime()const
{
if( currentTime_ >= endTime_ ) return true;
if( abs(currentTime_-endTime_) < 0.5*dt_ )return true;
if( std::abs(currentTime_-endTime_) < 0.5*dt_ )return true;
return false;
}
bool pFlow::timeControl::reachedStopAt()const
{
if( currentTime_ >= stopAt_ ) return true;
if( abs(currentTime_-stopAt_) < 0.5*dt_ )return true;
if( std::abs(currentTime_-stopAt_) < 0.5*dt_ )return true;
return false;
}
@ -154,7 +154,7 @@ void pFlow::timeControl::checkForOutputToFile()
bool save = false;
if(managedExternaly_)
{
if( abs(currentTime_-writeTime_) < 0.5*dt_)
if( std::abs(currentTime_-writeTime_) < 0.5*dt_)
{
save = true;
lastSaved_ = currentTime_;
@ -162,12 +162,12 @@ void pFlow::timeControl::checkForOutputToFile()
}
else
{
if ( abs(currentTime_ - lastSaved_ - saveInterval_) < 0.5 * dt_ )
if ( std::abs(currentTime_ - lastSaved_ - saveInterval_) < 0.5 * dt_ )
{
lastSaved_ = currentTime_;
save = true;
}
else if( abs(currentTime_ - lastSaved_) < min( pow(10.0,-1.0*timePrecision_), 0.5 *dt_) )
else if( std::abs(currentTime_ - lastSaved_) < std::min( pow(10.0,-1.0*timePrecision_), 0.5 *dt_) )
{
lastSaved_ = currentTime_;
save = true;

7
src/phasicFlow/smartPointers/uniquePtr.hpp
View File

@ -37,16 +37,15 @@ namespace pFlow
template<
typename T,
typename Deleter = std::default_delete<T>
typename T
>
class uniquePtr
:
public std::unique_ptr<T, Deleter>
public std::unique_ptr<T>
{
public:
using uniquePtrType = std::unique_ptr<T, Deleter>;
using uniquePtrType = std::unique_ptr<T>;
// using base constructors
using uniquePtrType::unique_ptr;

2
src/phasicFlow/streams/TStream/iTstream.cpp
View File

@ -325,7 +325,7 @@ void pFlow::iTstream::reset()
size_t pFlow::iTstream::tell()
{
notImplementedFunction;
return -1;
return static_cast<size_t>(-1);
}
const pFlow::tokenList& pFlow::iTstream::tokens()const

4
src/phasicFlow/streams/dataIO/dataIO.cpp
View File

@ -22,7 +22,7 @@ template<typename T>
bool pFlow::writeDataAsciiBinary(iOstream& os, span<T> data)
{
if( os.isBinary() )
if( std::is_trivially_copyable_v<T> && os.isBinary() )
{
// first write the number of data
uint64 numData = data.size();
@ -83,7 +83,7 @@ bool pFlow::readDataAsciiBinary
)
{
if(is.isBinary())
if(std::is_trivially_copyable_v<T> && is.isBinary())
{
data.clear();
// read length of data

5
src/phasicFlow/structuredData/boundaries/boundaryBase/boundaryBase.cpp
View File

@ -42,6 +42,7 @@ void pFlow::boundaryBase::setNewIndices
unSyncLists();
}
bool pFlow::boundaryBase::appendNewIndices
(
const uint32Vector_D& newIndices
@ -144,7 +145,7 @@ bool pFlow::boundaryBase::setRemoveKeepIndices
return true;
}
bool pFlow::boundaryBase::transferPoints
bool pFlow::boundaryBase::transferPointsToMirror
(
uint32 numTransfer,
const uint32Vector_D& transferMask,
@ -203,6 +204,8 @@ pFlow::boundaryBase::boundaryBase
indexList_(groupNames("indexList", dict.name())),
indexListHost_(groupNames("hostIndexList",dict.name())),
neighborLength_(dict.getVal<real>("neighborLength")),
updateInetrval_(dict.getVal<uint32>("updateInterval")),
boundaryExtntionLengthRatio_(dict.getVal<real>("boundaryExtntionLengthRatio")),
internal_(internal),
boundaries_(bndrs),
thisBoundaryIndex_(thisIndex),

39
src/phasicFlow/structuredData/boundaries/boundaryBase/boundaryBase.hpp
View File

@ -68,6 +68,12 @@ private:
/// The length defined for creating neighbor list
real neighborLength_;
/// time steps between successive update of boundary lists
uint32 updateInetrval_;
/// the extra boundary extension beyound actual limits of boundary
real boundaryExtntionLengthRatio_;
/// a reference to internal points
internalPoints& internal_;
@ -106,7 +112,7 @@ protected:
const uint32Vector_D& removeIndices,
const uint32Vector_D& keepIndices);
bool transferPoints(
bool transferPointsToMirror(
uint32 numTransfer,
const uint32Vector_D& transferMask,
uint32 transferBoundaryIndex,
@ -126,6 +132,12 @@ protected:
}
}
/// Is this iter the right time for updating bounday list
bool boundaryListUpdate(uint32 iter)const
{
return iter%updateInetrval_ == 0u || iter == 0u;
}
/// Update this boundary data in two steps (1 and 2).
/// This is called after calling beforeIteration for
/// all boundaries, so any particle addition, deletion,
@ -133,7 +145,7 @@ protected:
/// This two-step update help to have a flexible mechanism
/// for data transfer, mostly for MPI related jobs.
virtual
bool updataBoundary(int step)
bool updataBoundaryData(int step)
{
return true;
}
@ -145,11 +157,13 @@ protected:
/// to complete. false: if the operation is complete and no need for
/// additional step in operation.
virtual
bool transferData(int step)
bool transferData(uint32 iter, int step)
{
return false;
}
public:
TypeInfo("boundaryBase");
@ -189,20 +203,21 @@ public:
/// The length from boundary plane into the domain
/// where beyond that distance internal points exist.
/// By conventions is it always equal to neighborLength_
inline
real neighborLengthIntoInternal()const
{
return neighborLength_;
}
/// The distance length from boundary plane
/// where neighbor particles exist in that distance.
/// where neighbor particles still exist in that distance.
/// This length may be modified in each boundary type
/// as required. In this case the boundaryExtensionLength
/// method should also be modified accordingly.
virtual
real neighborLength()const
{
return neighborLength_;
return (1+boundaryExtntionLengthRatio_)*neighborLength_;
}
/// The extention length (in vector form) for the boundary
@ -213,7 +228,7 @@ public:
virtual
realx3 boundaryExtensionLength()const
{
return {0,0,0};
return -boundaryExtntionLengthRatio_*neighborLength_ * boundaryPlane_.normal();
}
inline
@ -351,6 +366,18 @@ public:
virtual
const realx3Vector_D& neighborProcPoints()const;
virtual
uint32 numToTransfer()const
{
return 0u;
}
virtual
uint32 numToRecieve()const
{
return 0u;
}
/// - static create
static
uniquePtr<boundaryBase> create

11
src/phasicFlow/structuredData/boundaries/boundaryBase/scatteredFieldAccess.hpp
View File

@ -114,6 +114,17 @@ public:
return fieldVals_;
}
auto& indices()
{
return indices_;
}
const auto& indices()const
{
return indices_;
}
INLINE_FUNCTION_HD
uint32 index(uint32 i)const
{
return indices_[i];

3
src/phasicFlow/structuredData/boundaries/boundaryExit/boundaryExit.cpp
View File

@ -47,6 +47,9 @@ bool pFlow::boundaryExit::beforeIteration
real dt
)
{
if( !boundaryListUpdate(iterNum))return true;
// nothing have to be done
if(empty())
{

34
src/phasicFlow/structuredData/boundaries/boundaryList.cpp
View File

@ -91,18 +91,18 @@ pFlow::boundaryList::updateNeighborLists()
pFlow::boundaryList::boundaryList(pointStructure& pStruct)
: ListPtr<boundaryBase>(pStruct.simDomain().sizeOfBoundaries()),
pStruct_(pStruct),
timeControl_(
pStruct.simDomain().subDict("boundaries"),
"update",
pStruct.currentTime()
neighborListUpdateInterval_(
pStruct.simDomain().subDict("boundaries").getVal<uint32>(
"neighborListUpdateInterval"
)
)
{
}
bool
pFlow::boundaryList::updateNeighborLists(uint32 iter, real t, real dt)
pFlow::boundaryList::updateNeighborLists(uint32 iter, bool force)
{
if (timeControl_.timeEvent(iter, t, dt))
if(iter%neighborListUpdateInterval_==0u || iter == 0u || force)
{
return updateNeighborLists();
}
@ -155,15 +155,18 @@ pFlow::boundaryList::internalDomainBox() const
}
bool
pFlow::boundaryList::beforeIteration(uint32 iter, real t, real dt)
pFlow::boundaryList::beforeIteration(uint32 iter, real t, real dt, bool force)
{
// it is time to update lists
if (timeControl_.timeEvent(iter, t, dt) && !updateNeighborLists())
if(!updateNeighborLists(iter , force) )
{
fatalErrorInFunction;
return false;
}
// this forces performing updating the list on each boundary
if(force) iter = 0;
for (auto bdry : *this)
{
if (!bdry->beforeIteration(iter, t, dt))
@ -176,12 +179,12 @@ pFlow::boundaryList::beforeIteration(uint32 iter, real t, real dt)
for (auto bdry : *this)
{
bdry->updataBoundary(1);
bdry->updataBoundaryData(1);
}
for (auto bdry : *this)
{
bdry->updataBoundary(2);
bdry->updataBoundaryData(2);
}
return true;
@ -215,20 +218,11 @@ pFlow::boundaryList::afterIteration(uint32 iter, real t, real dt)
{
if(requireStep[i])
{
requireStep[i] = this->operator[](i).transferData(step);
requireStep[i] = this->operator[](i).transferData(iter, step);
}
}
step++;
}
/*for (auto& bdry : *this)
{
if (!bdry->afterIteration(iter, t, dt))
{
fatalErrorInFunction << "Error in afterIteration in boundary "
<< bdry->name() << endl;
return false;
}
}*/
return true;
}

14
src/phasicFlow/structuredData/boundaries/boundaryList.hpp
View File

@ -17,15 +17,12 @@ Licence:
implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
-----------------------------------------------------------------------------*/
#ifndef __boundaryList_hpp__
#define __boundaryList_hpp__
#include "domain.hpp"
#include "boundaryBase.hpp"
#include "ListPtr.hpp"
#include "baseTimeControl.hpp"
namespace pFlow
@ -42,7 +39,7 @@ private:
//// - data members
pointStructure& pStruct_;
baseTimeControl timeControl_;
uint32 neighborListUpdateInterval_;
domain extendedDomain_;
@ -72,7 +69,7 @@ public:
/// @brief update neighbor list of boundaries based on
/// the time intervals
bool updateNeighborLists(uint32 iter, real t, real dt);
bool updateNeighborLists(uint32 iter, bool force = false);
bool createBoundaries();
@ -94,11 +91,6 @@ public:
return ListPtr<boundaryBase>::operator[](i);
}
inline
const baseTimeControl& timeControl()const
{
return timeControl_;
}
inline
const auto& extendedDomain()const
@ -114,7 +106,7 @@ public:
box internalDomainBox()const;
bool beforeIteration(uint32 iter, real t, real dt);
bool beforeIteration(uint32 iter, real t, real dt, bool force = false);
bool iterate(uint32 iter, real t, real dt);

8
src/phasicFlow/structuredData/boundaries/boundaryNone/boundaryNone.hpp
View File

@ -43,7 +43,7 @@ public:
boundaryList &bndrs,
uint32 thisIndex);
~boundaryNone() override= default;
~boundaryNone() final= default;
add_vCtor
(
@ -52,11 +52,11 @@ public:
dictionary
);
bool beforeIteration(uint32 iterNum, real t, real dt) override;
bool beforeIteration(uint32 iterNum, real t, real dt) final;
bool iterate(uint32 iterNum, real t, real dt) override;
bool iterate(uint32 iterNum, real t, real dt) final;
bool afterIteration(uint32 iterNum, real t, real dt) override;
bool afterIteration(uint32 iterNum, real t, real dt) final;
};

15
src/phasicFlow/structuredData/boundaries/boundaryPeriodic/boundaryPeriodic.cpp
View File

@ -33,17 +33,20 @@ pFlow::boundaryPeriodic::boundaryPeriodic
)
:
boundaryBase(dict, bplane, internal, bndrs, thisIndex),
mirrorBoundaryIndex_(dict.getVal<uint32>("mirrorBoundaryIndex"))
{}
mirrorBoundaryIndex_(dict.getVal<uint32>("mirrorBoundaryIndex")),
extensionLength_(dict.getVal<real>("boundaryExtntionLengthRatio"))
{
extensionLength_ = max(extensionLength_, static_cast<real>(0.1));
}
pFlow::real pFlow::boundaryPeriodic::neighborLength() const
{
return (1+extensionLength_)*boundaryBase::neighborLength();
return (1+extensionLength_)*neighborLengthIntoInternal();
}
pFlow::realx3 pFlow::boundaryPeriodic::boundaryExtensionLength() const
{
return -extensionLength_*neighborLength()*boundaryBase::boundaryPlane().normal();
return -extensionLength_*neighborLengthIntoInternal()*boundaryBase::boundaryPlane().normal();
}
@ -58,6 +61,8 @@ bool pFlow::boundaryPeriodic::beforeIteration(
return true;
}
if( !boundaryListUpdate(iterNum))return true;
uint32 s = size();
uint32Vector_D transferFlags("transferFlags",s+1, s+1, RESERVE());
transferFlags.fill(0u);
@ -92,7 +97,7 @@ bool pFlow::boundaryPeriodic::beforeIteration(
// to obtain the transfer vector
realx3 transferVec = displacementVectroToMirror();
return transferPoints
return transferPointsToMirror
(
numTransfered,
transferFlags,

3
src/phasicFlow/structuredData/boundaries/boundaryPeriodic/boundaryPeriodic.hpp
View File

@ -35,8 +35,7 @@ private:
uint32 mirrorBoundaryIndex_;
static
inline const real extensionLength_ = 0.1;
real extensionLength_ = 0.1;
public:

4
src/phasicFlow/structuredData/boundaries/boundaryReflective/boundaryReflective.cpp
View File

@ -126,6 +126,8 @@ bool pFlow::boundaryReflective::afterIteration
const auto& velocity = time().lookupObject<realx3PointField_D>(velocityName_);
const auto& velocityD = velocity.deviceViewAll();
const auto restitution = restitution_;
Kokkos::parallel_for(
"pFlow::boundaryReflective::velocityChange",
deviceRPolicyStatic(0,numInContact),
@ -137,7 +139,7 @@ bool pFlow::boundaryReflective::afterIteration
if(vn < 0)
{
realx3 vt = vel - vn*p.normal();
vel = restitution_*(vt - vn*p.normal());
vel = restitution*(vt - vn*p.normal());
}
}
);

105
src/phasicFlow/structuredData/cells/cells.hpp
View File

@ -28,58 +28,40 @@ Licence:
namespace pFlow
{
template<typename indexType>
class cells
{
public:
using CellType = triple<indexType>;
protected:
private:
// - domain
box domain_{realx3(0.0), realx3(1.0)};
box domainBox_{realx3(0.0), realx3(1.0)};
// - cell size
realx3 cellSize_{1,1,1};
CellType numCells_{1,1,1};
real celldx_{1};
int32x3 numCells_{1,1,1};
// - protected methods
INLINE_FUNCTION_H
void calculate()
{
numCells_ = (domain_.maxPoint()-domain_.minPoint())/cellSize_ + realx3(1.0);
numCells_ = max( numCells_ , CellType(static_cast<indexType>(1)) );
numCells_ = (domainBox_.maxPoint()-domainBox_.minPoint())/celldx_ + realx3(1.0);
numCells_ = max( numCells_ , int32x3(1) );
}
public:
INLINE_FUNCTION_HD
cells()
{}
cells() = default;
INLINE_FUNCTION_H
cells(const box& domain, real cellSize)
:
domain_(domain),
cellSize_(cellSize)
domainBox_(domain),
celldx_(cellSize)
{
calculate();
}
INLINE_FUNCTION_H
cells(const box& domain, int32 nx, int32 ny, int32 nz)
:
domain_(domain),
cellSize_(
(domain_.maxPoint() - domain_.minPoint())/realx3(nx, ny, nz)
),
numCells_(nx, ny, nz)
{}
INLINE_FUNCTION_HD
cells(const cells&) = default;
@ -100,43 +82,36 @@ public:
INLINE_FUNCTION_H
void setCellSize(real cellSize)
{
cellSize_ = cellSize;
calculate();
}
INLINE_FUNCTION_H
void setCellSize(realx3 cellSize)
{
cellSize_ = cellSize;
celldx_ = cellSize;
calculate();
}
INLINE_FUNCTION_HD
realx3 cellSize()const
real cellSize()const
{
return cellSize_;
return celldx_;
}
INLINE_FUNCTION_HD
const CellType& numCells()const
const int32x3& numCells()const
{
return numCells_;
}
INLINE_FUNCTION_HD
indexType nx()const
int32 nx()const
{
return numCells_.x();
}
INLINE_FUNCTION_HD
indexType ny()const
int32 ny()const
{
return numCells_.y();
}
INLINE_FUNCTION_HD
indexType nz()const
int32 nz()const
{
return numCells_.z();
}
@ -149,22 +124,21 @@ public:
static_cast<int64>(numCells_.z());
}
const auto& domain()const
const auto& domainBox()const
{
return domain_;
return domainBox_;
}
INLINE_FUNCTION_HD
CellType pointIndex(const realx3& p)const
int32x3 pointIndex(const realx3& p)const
{
return CellType( (p - domain_.minPoint())/cellSize_ );
return int32x3( (p - domainBox_.minPoint())/celldx_ );
}
INLINE_FUNCTION_HD
bool pointIndexInDomain(const realx3 p, CellType& index)const
bool pointIndexInDomain(const realx3 p, int32x3& index)const
{
if( !domain_.isInside(p) ) return false;
if(!inDomain(p))return false;
index = this->pointIndex(p);
return true;
}
@ -172,11 +146,11 @@ public:
INLINE_FUNCTION_HD
bool inDomain(const realx3& p)const
{
return domain_.isInside(p);
return domainBox_.isInside(p);
}
INLINE_FUNCTION_HD
bool isInRange(const CellType& cell)const
bool inCellRange(const int32x3& cell)const
{
if(cell.x()<0)return false;
if(cell.y()<0)return false;
@ -188,7 +162,7 @@ public:
}
INLINE_FUNCTION_HD
bool isInRange(indexType i, indexType j, indexType k)const
bool inCellRange(int32 i, int32 j, int32 k)const
{
if(i<0)return false;
if(j<0)return false;
@ -199,21 +173,6 @@ public:
return true;
}
INLINE_FUNCTION_HD
void extendBox(
const CellType& p1,
const CellType& p2,
const CellType& p3,
indexType extent,
CellType& minP,
CellType& maxP)const
{
minP = min( min( p1, p2), p3)-extent;
maxP = max( max( p1, p2), p3)+extent;
minP = bound(minP);
maxP = bound(maxP);
}
INLINE_FUNCTION_HD
void extendBox(
@ -224,17 +183,17 @@ public:
realx3& minP,
realx3& maxP)const
{
minP = min(min(p1,p2),p3) - extent*cellSize_ ;
maxP = max(max(p1,p2),p3) + extent*cellSize_ ;
minP = min(min(p1,p2),p3) - extent*celldx_ ;
maxP = max(max(p1,p2),p3) + extent*celldx_ ;
minP = bound(minP);
maxP = bound(maxP);
}
INLINE_FUNCTION_HD
CellType bound(CellType p)const
int32x3 bound(int32x3 p)const
{
return CellType(
return int32x3(
min( numCells_.x()-1, max(0,p.x())),
min( numCells_.y()-1, max(0,p.y())),
min( numCells_.z()-1, max(0,p.z()))
@ -245,9 +204,9 @@ public:
realx3 bound(realx3 p)const
{
return realx3(
min( domain_.maxPoint().x(), max(domain_.minPoint().x(),p.x())),
min( domain_.maxPoint().y(), max(domain_.minPoint().y(),p.y())),
min( domain_.maxPoint().z(), max(domain_.minPoint().z(),p.z()))
min( domainBox_.maxPoint().x(), max(domainBox_.minPoint().x(),p.x())),
min( domainBox_.maxPoint().y(), max(domainBox_.minPoint().y(),p.y())),
min( domainBox_.maxPoint().z(), max(domainBox_.minPoint().z(),p.z()))
);
}
};

2
src/phasicFlow/structuredData/cylinder/cylinder.cpp
View File

@ -18,7 +18,7 @@ Licence:
-----------------------------------------------------------------------------*/
#include "math.hpp"
#include "cylinder.hpp"
#include "zAxis.hpp"
#include "streams.hpp"

13
src/phasicFlow/structuredData/domain/regularSimulationDomain.cpp
View File

@ -31,7 +31,10 @@ bool pFlow::regularSimulationDomain::createBoundaryDicts()
this->addDict("regularBoundaries", boundaries);
auto& rbBoundaries = this->subDict("regularBoundaries");
real neighborLength = boundaries.getVal<real>("neighborLength");
auto neighborLength = boundaries.getVal<real>("neighborLength");
auto boundaryExtntionLengthRatio =
boundaries.getValOrSet<real>("boundaryExtntionLengthRatio", 0.1);
auto updateIntercal = boundaries.getValOrSet<uint32>("updateInterval", 1u);
for(uint32 i=0; i<sizeOfBoundaries(); i++)
{
@ -51,6 +54,14 @@ bool pFlow::regularSimulationDomain::createBoundaryDicts()
return false;
}
if(!bDict.addOrReplace("updateInterval", updateIntercal))
{
fatalErrorInFunction<<"error in adding updateIntercal to "<< bName <<
"in dictionary "<< boundaries.globalName()<<endl;
}
bDict.addOrReplace("boundaryExtntionLengthRatio", boundaryExtntionLengthRatio);
if(!bDict.add("neighborProcessorNo",(uint32) processors::globalRank()))
{
fatalErrorInFunction<<"error in adding neighborProcessorNo to "<< bName <<

136
src/phasicFlow/structuredData/pointStructure/internalPoints/internalPoints.cpp
View File

@ -142,13 +142,24 @@ bool pFlow::internalPoints::changePointsFlagPosition
return true;
}
bool pFlow::internalPoints::sortPoints(const uint32IndexContainer &sortedIndices)
{
if(!pointPosition_.reorderItems(sortedIndices))
{
fatalErrorInFunction;
return false;
}
pFlagsD_.resetFlags(pFlagsD_.capacity(), 0u, sortedIndices.size());
unSyncFlag();
return true;
}
pFlow::internalPoints::internalPoints()
:
subscriber("internalPoints"),
: subscriber("internalPoints"),
pointPosition_("internalPoints", "internalPoints", initialCapacity_, 0, RESERVE()),
pFlagsD_(initialCapacity_, 0 , 0)
pFlagsD_(initialCapacity_, 0, 0)
{}
@ -306,6 +317,9 @@ pFlow::internalPoints::insertPoints(
auto aRange = pFlagsD_.activeRange();
uint32 emptySpots = pFlagsD_.capacity() - pFlagsD_.numActive();
if(emptySpots!= 0) emptySpots--;
message msg;
if( numNew > emptySpots )
@ -337,10 +351,22 @@ pFlow::internalPoints::insertPoints(
// we should fill the scattered empty spots
else
{
notImplementedFunction;
auto newIndices = pFlagsD_.getEmptyPoints(numNew);
if(numNew != newIndices.size())
{
fatalErrorInFunction<<"not enough empty points in pointFlag"<<
numNew<< " "<<newIndices.size() <<endl;
pOutput<< pFlagsD_.capacity()<<endl;
pOutput<< pFlagsD_.numActive()<<endl;
return false;
}
varList.emplaceBack<uint32IndexContainer>(
msg.addAndName(message::ITEM_INSERT),
newIndices
);
}
const auto& indices = varList.getObject<uint32IndexContainer>(
message::eventName(message::ITEM_INSERT)
);
@ -398,6 +424,108 @@ pFlow::internalPoints::insertPoints(
return true;
}
bool pFlow::internalPoints::insertPointsOnly(
const realx3Vector_D &points,
message& msg,
anyList &varList
)
{
uint32 numNew = static_cast<uint32>(points.size());
auto aRange = pFlagsD_.activeRange();
uint32 emptySpots = pFlagsD_.capacity() - pFlagsD_.numActive();
if(emptySpots!= 0) emptySpots--;
if( numNew > emptySpots )
{
// increase the capacity to hold new points
aRange = pFlagsD_.activeRange();
uint32 newCap = pFlagsD_.changeCapacity(numNew);
unSyncFlag();
varList.emplaceBack(
msg.addAndName(message::CAP_CHANGED),
newCap);
}
// first check if it is possible to add to the beggining of list
if(numNew <= aRange.start())
{
varList.emplaceBack<uint32IndexContainer>(
msg.addAndName(message::ITEM_INSERT),
0u, numNew);
}// check if it is possible to add to the end of the list
else if( numNew <= pFlagsD_.capacity() - aRange.end() )
{
varList.emplaceBack<uint32IndexContainer>(
msg.addAndName(message::ITEM_INSERT),
aRange.end(), aRange.end()+numNew);
}// we should fill the scattered empty spots
else
{
auto newIndices = pFlagsD_.getEmptyPoints(numNew);
if(numNew != newIndices.size())
{
fatalErrorInFunction<<"not enough empty points in pointFlag"<<
numNew<< " "<<newIndices.size() <<endl;
pOutput<< pFlagsD_.capacity()<<endl;
pOutput<< pFlagsD_.numActive()<<endl;
return false;
}
varList.emplaceBack<uint32IndexContainer>(
msg.addAndName(message::ITEM_INSERT),
newIndices
);
}
const auto& indices = varList.getObject<uint32IndexContainer>(
message::eventName(message::ITEM_INSERT)
);
auto nAdded = pFlagsD_.addInternalPoints(indices.deviceView());
unSyncFlag();
if(nAdded != numNew )
{
fatalErrorInFunction;
return false;
}
pointPosition_.reserve( pFlagsD_.capacity() );
if(!pointPosition_.insertSetElement(indices, points.deviceView()))
{
fatalErrorInFunction<<
"Error in inserting new positions into pointPosition"<<
" internalPoints field"<<endl;
return false;
}
auto newARange = pFlagsD_.activeRange();
if( aRange.end() != newARange.end() )
{
varList.emplaceBack(
msg.addAndName(message::RANGE_CHANGED),
newARange);
varList.emplaceBack(
msg.addAndName(message::SIZE_CHANGED),
newARange.end());
}
else if(aRange.start() != newARange.start())
{
varList.emplaceBack(
msg.addAndName(message::RANGE_CHANGED),
newARange);
}
return true;
}
bool pFlow::internalPoints::read
(
iIstream& is

4
src/phasicFlow/structuredData/pointStructure/internalPoints/internalPoints.hpp
View File

@ -101,6 +101,8 @@ protected:
pFlagsD_ = pFlagTypeDevice(cap, start, end);
}
bool sortPoints(const uint32IndexContainer& sortedIndices);
public:
//friend class dynamicinternalPoints;
@ -236,6 +238,8 @@ public:
bool insertPoints(const realx3Vector& points, anyList& varList);
bool insertPointsOnly(const realx3Vector_D& points, message& msg, anyList& varList);
//// - IO operations

27
src/phasicFlow/structuredData/pointStructure/internalPoints/pointFlag.hpp
View File

@ -57,6 +57,8 @@ class pointFlag
protected:
friend class internalPoints;
viewType flags_;
uint32 numActive_ = 0;
@ -108,9 +110,29 @@ protected:
return flg;
}
void resetFlags(uint32 cap, uint32 start, uint32 end)
{
if(cap != capacity() )
{
reallocFill(flags_, cap, static_cast<uint8>(Flag::DELETED));
}
else
{
fill(flags_, end, cap, static_cast<uint8>(Flag::DELETED));
}
activeRange_ = {start, end};
numActive_ = activeRange_.numElements();
isAllActive_ = true;
fill(flags_, activeRange_, static_cast<uint8>(Flag::INTERNAL));
nLeft_ = nRight_ = 0;
nBottom_ = nTop_ = 0;
nRear_ = nFront_ = 0;
}
public:
friend class internalPoints;
pointFlag() = default;
@ -285,6 +307,9 @@ public:
ViewType1D<uint32, memory_space> getActivePoints();
ViewType1D<uint32, memory_space> getEmptyPoints(uint32 numToGet)const;
/// @brief Loop over the active points and mark those out of the box
/// and return number of deleted points
/// @param validBox the box whose inside is valid

60
src/phasicFlow/structuredData/pointStructure/internalPoints/pointFlagKernels.hpp
View File

@ -50,6 +50,52 @@ pFlow::ViewType1D<pFlow::uint32, typename pFlow::pointFlag<ExecutionSpace>::memo
return aPoints;
}
template<typename ExecutionSpace>
pFlow::ViewType1D<pFlow::uint32, typename pFlow::pointFlag<ExecutionSpace>::memory_space>
pFlow::pointFlag<ExecutionSpace>::getEmptyPoints(uint32 numToGet)const
{
uint32 cap = capacity();
using rpAPoints = Kokkos::RangePolicy<execution_space,
Kokkos::IndexType<uint32>>;
ViewType1D<uint32,memory_space> indices("indices", cap+1);
ViewType1D<uint32,memory_space> emptyPoints("emptyPoints", numToGet);
auto flags = flags_;
Kokkos::parallel_for(
"getEmptyPoints",
rPolicy(0, cap),
LAMBDA_HD(uint32 i){
indices(i) = flags[i] == DELETED;
});
Kokkos::fence();
exclusiveScan(indices, 0u, cap, indices, 0u);
uint32 numFound = 0;
Kokkos::parallel_reduce(
"fillEmptyPoints",
rpAPoints(0, cap-1),
LAMBDA_HD(uint32 i, uint32& numFoundUpdate){
if(indices(i)!= indices(i+1) && indices(i)< numToGet)
{
emptyPoints(indices(i)) = i;
numFoundUpdate++;
}
},
numFound
);
if(numFound < numToGet)
{
return Kokkos::subview(emptyPoints, Kokkos::make_pair(0u, numFound));
}
return emptyPoints;
}
template<typename ExecutionSpace>
pFlow::uint32 pFlow::pointFlag<ExecutionSpace>::markOutOfBoxDelete
@ -106,7 +152,7 @@ pFlow::uint32 pFlow::pointFlag<ExecutionSpace>::markOutOfBoxDelete
{
minRange = 0;
maxRange = 0;
numDeleted == numActive_;
numDeleted = numActive_;
}
else
{
@ -135,6 +181,7 @@ pFlow::pointFlag<ExecutionSpace>::addInternalPoints
uint32 maxRange;
uint32 numAdded = 0;
const auto& flag = flags_;
Kokkos::parallel_reduce(
"pointFlagKernels::addInternalPoints",
@ -146,10 +193,10 @@ pFlow::pointFlag<ExecutionSpace>::addInternalPoints
uint32& addToUpdate)
{
uint32 idx = points(i);
if( flags_[idx] <= DELETED) addToUpdate ++;
if( flag[idx] <= DELETED) addToUpdate ++;
minUpdate = min(minUpdate,idx);
maxUpdate = max(maxUpdate,idx);
flags_[idx] = INTERNAL;
flag[idx] = INTERNAL;
},
Kokkos::Min<uint32>(minRange),
Kokkos::Max<uint32>(maxRange),
@ -202,7 +249,7 @@ bool pFlow::pointFlag<ExecutionSpace>::deletePoints
if(numDeleted >= numActive_)
{
activeRange_ = {0, 0};
numDeleted == numActive_;
numDeleted = numActive_;
}
numActive_ = numActive_ - numDeleted;
@ -245,7 +292,7 @@ bool pFlow::pointFlag<ExecutionSpace>::deletePoints
if(numDeleted >= numActive_)
{
activeRange_ = {0, 0};
numDeleted == numActive_;
numDeleted = numActive_;
}
numActive_ = numActive_ - numDeleted;
@ -264,13 +311,14 @@ bool pFlow::pointFlag<ExecutionSpace>::changeFlags
)
{
auto flg = getBoundaryFlag(boundaryIndex);
const auto& flags = flags_;
Kokkos::parallel_for
(
"pointFlag::changeFlags",
rPolicy(0, changePoints.size()),
LAMBDA_HD(uint32 i)
{
flags_[changePoints(i)] = flg;
flags[changePoints(i)] = flg;
}
);
Kokkos::fence();

58
src/phasicFlow/structuredData/pointStructure/pointStructure/mortonIndexing.cpp → src/phasicFlow/structuredData/pointStructure/pointStructure/pointSorting/mortonIndexing.cpp
View File

@ -17,66 +17,60 @@ Licence:
implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
-----------------------------------------------------------------------------*/
#include "mortonIndexing.hpp"
#include "cells.hpp"
#include "streams.hpp"
bool pFlow::getSortedIndex(
pFlow::uint32IndexContainer pFlow::getSortedIndices(
box boundingBox,
real dx,
range activeRange,
ViewType1D<realx3> pos,
ViewType1D<int8> flag,
int32IndexContainer& sortedIndex)
const ViewType1D<realx3>& pos,
const pFlagTypeDevice& flag
)
{
if(flag.numActive() == 0u)return uint32IndexContainer();
// obtain the morton code of the particles
cells<size_t> allCells( boundingBox, dx);
int32IndexContainer index(activeRange.first, activeRange.second);
cells allCells( boundingBox, dx);
auto aRange = flag.activeRange();
ViewType1D<uint64_t> mortonCode("mortonCode", activeRange.second);
uint32IndexContainer sortedIndex(aRange.start(), aRange.end());
output<<"before first kernel"<<endl;;
ViewType1D<uint64_t> mortonCode("mortonCode", aRange.end());
using rpMorton =
Kokkos::RangePolicy<Kokkos::IndexType<int32>>;
int32 numActive = 0;
Kokkos::parallel_reduce
Kokkos::parallel_for
(
"mortonIndexing::getIndex::morton",
rpMorton(activeRange.first, activeRange.second),
LAMBDA_HD(int32 i, int32& sumToUpdate){
if( flag[i] == 1 ) // active point
deviceRPolicyStatic(aRange.start(), aRange.end()),
LAMBDA_HD(uint32 i){
if( flag.isActive(i)) // active point
{
auto cellInd = allCells.pointIndex(pos[i]);
mortonCode[i] = xyzToMortonCode64(cellInd.x(), cellInd.y(), cellInd.z());
sumToUpdate++;
}else
{
mortonCode[i] = xyzToMortonCode64
(
static_cast<uint64_t>(-1),
static_cast<uint64_t>(-1),
static_cast<uint64_t>(-1)
largestPosInt32,
largestPosInt32,
largestPosInt32
);
}
},
numActive
}
);
Kokkos::fence();
permuteSort(
mortonCode,
activeRange.first,
activeRange.second,
index.deviceView(),
aRange.start(),
aRange.end(),
sortedIndex.deviceView(),
0 );
index.modifyOnDevice();
index.setSize(numActive);
index.syncViews();
sortedIndex = index;
sortedIndex.modifyOnDevice();
sortedIndex.syncViews(flag.numActive());
return true;
return sortedIndex;
}

10
src/phasicFlow/structuredData/pointStructure/pointStructure/mortonIndexing.hpp → src/phasicFlow/structuredData/pointStructure/pointStructure/pointSorting/mortonIndexing.hpp
View File

@ -17,24 +17,22 @@ Licence:
implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
-----------------------------------------------------------------------------*/
#ifndef __mortonIndexing_hpp__
#define __mortonIndexing_hpp__
#include "types.hpp"
#include "box.hpp"
#include "indexContainer.hpp"
#include "pointFlag.hpp"
namespace pFlow
{
bool getSortedIndex(
uint32IndexContainer getSortedIndices(
box boundingBox,
real dx,
range activeRange,
ViewType1D<realx3> pos,
ViewType1D<int8> flag,
int32IndexContainer& sortedIndex);
const ViewType1D<realx3>& pos,
const pFlagTypeDevice& flag);
INLINE_FUNCTION_HD

52
src/phasicFlow/structuredData/pointStructure/pointStructure/pointSorting/pointSorting.cpp Normal file
View File

@ -0,0 +1,52 @@
/*------------------------------- phasicFlow ---------------------------------
O C enter of
O O E ngineering and
O O M ultiscale modeling of
OOOOOOO F luid flow
------------------------------------------------------------------------------
Copyright (C): www.cemf.ir
email: hamid.r.norouzi AT gmail.com
------------------------------------------------------------------------------
Licence:
This file is part of phasicFlow code. It is a free software for simulating
granular and multiphase flows. You can redistribute it and/or modify it under
the terms of GNU General Public License v3 or any other later versions.
phasicFlow is distributed to help others in their research in the field of
granular and multiphase flows, but WITHOUT ANY WARRANTY; without even the
implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
-----------------------------------------------------------------------------*/
#include "pointSorting.hpp"
#include "mortonIndexing.hpp"
pFlow::pointSorting::pointSorting(const dictionary & dict)
:
performSorting_(dict.getValOrSet("active", Logical(false))),
timeControl_(
performSorting_()?
baseTimeControl(dict, "sorting"):
baseTimeControl(0,1,1, "sorting")
),
dx_(
performSorting_()?
dict.getVal<real>("dx"):
1.0
)
{
if( performSorting_() )
REPORT(1)<<"Point sorting is "<<Yellow_Text("active")<<" in simulation"<<END_REPORT;
else
REPORT(1)<<"Point sorting is "<<Yellow_Text("inactive")<<" in simulation"<<END_REPORT;
}
pFlow::uint32IndexContainer
pFlow::pointSorting::getSortedIndices(
const box& boundingBox,
const ViewType1D<realx3> &pos,
const pFlagTypeDevice &flag
) const
{
return pFlow::getSortedIndices(boundingBox, dx_, pos, flag);
}

68
src/phasicFlow/structuredData/pointStructure/pointStructure/pointSorting/pointSorting.hpp Normal file
View File

@ -0,0 +1,68 @@
/*------------------------------- 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 __pointSorting_hpp__
#define __pointSorting_hpp__
#include "baseTimeControl.hpp"
#include "indexContainer.hpp"
#include "pointFlag.hpp"
namespace pFlow
{
class box;
class pointSorting
{
private:
Logical performSorting_;
baseTimeControl timeControl_;
real dx_;
public:
explicit pointSorting(const dictionary& dict);
bool performSorting()const
{
return performSorting_();
}
bool sortTime(uint32 iter, real t, real dt)const
{
return performSorting_() && timeControl_.timeEvent(iter, t, dt);
}
uint32IndexContainer getSortedIndices(
const box& boundingBox,
const ViewType1D<realx3>& pos,
const pFlagTypeDevice& flag
)const;
};
}
#endif

58
src/phasicFlow/structuredData/pointStructure/pointStructure/pointStructure.cpp
View File

@ -22,6 +22,7 @@ Licence:
#include "pointStructure.hpp"
#include "systemControl.hpp"
#include "vocabs.hpp"
#include "anyList.hpp"
bool pFlow::pointStructure::setupPointStructure(const realx3Vector& points)
{
@ -73,6 +74,7 @@ bool pFlow::pointStructure::setupPointStructure(const PointsTypeHost &points)
}
boundaries_.createBoundaries();
boundaries_.updateNeighborLists(0u, true);
return true;
}
@ -111,10 +113,13 @@ pFlow::pointStructure::pointStructure
(
simulationDomain::create(control)
),
pointSorting_(simulationDomain_->subDictOrCreate("pointSorting")),
boundaries_
(
*this
)
),
boundaryUpdateTimer_("boundaryUpdate", &timers()),
boundaryDataTransferTimer_("boundaryDataTransferTimer", &timers())
{
REPORT(1)<< "Reading "<< Green_Text("pointStructure")<<
" from "<<IOobject::path()<<END_REPORT;
@ -150,6 +155,7 @@ pFlow::pointStructure::pointStructure(
(
simulationDomain::create(control)
),
pointSorting_(simulationDomain_->subDictOrCreate("pointSorting")),
boundaries_
(
*this
@ -165,12 +171,58 @@ pFlow::pointStructure::pointStructure(
bool pFlow::pointStructure::beforeIteration()
{
if( !boundaries_.beforeIteration(currentIter(), currentTime(), dt()) )
uint32 iter = currentIter();
real t = currentTime();
real deltat = dt();
if(pointSorting_.sortTime(iter, t, deltat))
{
auto sortedIndices = pointSorting_.getSortedIndices(
simulationDomain_().globalBox(),
pointPositionDevice(),
activePointsMaskDevice()
);
if( !sortPoints(sortedIndices) )
{
fatalErrorInFunction;
return false;
}
boundaryUpdateTimer_.start();
boundaries_.beforeIteration(iter, t, deltat, true);
boundaryUpdateTimer_.end();
INFORMATION<<"Reordering of particles has been done. New active range for particles is "<<
activeRange()<<END_INFO;
message msg;
anyList varList;
varList.emplaceBack(
msg.addAndName(message::ITEM_REARRANGE),
sortedIndices);
if(!notify(iter, t, deltat, msg, varList))
{
fatalErrorInFunction<<
"cannot notify for reordering items."<<endl;
return false;
}
return true;
// notify others about this change
}
else
{
boundaryUpdateTimer_.start();
if( !boundaries_.beforeIteration(iter, t, deltat) )
{
fatalErrorInFunction<<
"Unable to perform beforeIteration for boundaries"<<endl;
return false;
}
boundaryUpdateTimer_.end();
}
return true;
}
@ -188,12 +240,14 @@ bool pFlow::pointStructure::iterate()
bool pFlow::pointStructure::afterIteration()
{
boundaryDataTransferTimer_.start();
if( !boundaries_.afterIteration(currentIter(), currentTime(), dt()) )
{
fatalErrorInFunction<<
"Unable to perform afterIteration for boundaries"<<endl;
return false;
}
boundaryDataTransferTimer_.end();
return true;
}

9
src/phasicFlow/structuredData/pointStructure/pointStructure/pointStructure.hpp
View File

@ -24,6 +24,7 @@ Licence:
#include "demComponent.hpp"
#include "internalPoints.hpp"
#include "simulationDomain.hpp"
#include "pointSorting.hpp"
#include "boundaryList.hpp"
#include "streams.hpp"
@ -53,8 +54,16 @@ private:
//// - data members
uniquePtr<simulationDomain> simulationDomain_ = nullptr;
pointSorting pointSorting_;
boundaryList boundaries_;
Timer boundaryUpdateTimer_;
Timer boundaryDataTransferTimer_;
bool setupPointStructure(const realx3Vector& points);
bool setupPointStructure(const PointsTypeHost& points);

4
src/phasicFlow/structuredData/zAxis/array2D.hpp
View File

@ -33,12 +33,12 @@ struct array2D
constexpr size_t nCols()const noexcept
{
return nCols;
return nCol;
}
constexpr size_t nRows()const noexcept
{
return nRows;
return nRow;
}
};

Some files were not shown because too many files have changed in this diff Show More