xfygogo/phasicFlow
1
0
Fork 0
mirror of https://github.com/PhasicFlow/phasicFlow.git synced 2025-07-08 03:07:03 +00:00
Code Issues Packages Projects Releases Wiki Activity

Merge branch 'develop' into MPI

Browse Source 
- This is the merge from develop branch to get latest update from the this branch to start MPI developemet.
- This stage boundaries and data exchange between processors through boundaries are handled.
This commit is contained in:
Hamidreza Norouzi
2024-04-18 10:19:19 -07:00
parent 656e03de36 e4c387a801
commit 5864cea7e8
408 changed files with 22048 additions and 16233 deletions
Download Patch File Download Diff File Expand all files Collapse all files

7
.gitignore vendored
View File

@ -1,6 +1,12 @@
# files
.clang-format
.vscode
.dependencygraph
# Prerequisites
*.d
# Compiled Object files
*.slo
*.lo
@ -38,7 +44,6 @@ bin/**
lib/**
test*/**
**/**notnow
**/MPIParallel*/**
doc/code-documentation/
doc/DTAGS
# all possible time folders

6
CMakeLists.txt
View File

@ -6,7 +6,7 @@ project(phasicFlow VERSION 1.0 )
set(CMAKE_CXX_STANDARD 17 CACHE STRING "" FORCE)
set(CMAKE_CXX_STANDARD_REQUIRED True)
set(CMAKE_INSTALL_PREFIX ${phasicFlow_SOURCE_DIR} CACHE PATH "Install path of phasicFlow" FORCE)
set(CMAKE_BUILD_TYPE Release CACHE STRING "build type" FORCE)
set(CMAKE_BUILD_TYPE Debug CACHE STRING "build type" FORCE)
set(BUILD_SHARED_LIBS ON CACHE BOOL "Build using shared libraries" FORCE)
mark_as_advanced(FORCE var BUILD_SHARED_LIBS)
@ -71,9 +71,9 @@ include_directories(src/setHelpers src/demComponent "${PROJECT_BINARY_DIR}")
add_subdirectory(src)
#add_subdirectory(solvers)
add_subdirectory(solvers)
#add_subdirectory(utilities)
add_subdirectory(utilities)
#add_subdirectory(DEMSystems)
add_subdirectory(testIO)

2
solvers/CMakeLists.txt
View File

@ -1,6 +1,6 @@
#add_subdirectory(iterateSphereParticles)
add_subdirectory(iterateSphereParticles)
add_subdirectory(iterateGeometry)

10
solvers/iterateGeometry/iterateGeometry.cpp
View File

@ -29,12 +29,13 @@ Licence:
* folder.
*/
#include "vocabs.hpp"
#include "systemControl.hpp"
#include "geometryMotion.hpp"
#include "geometry.hpp"
#include "commandLine.hpp"
#include "readControlDict.hpp"
//#include "readControlDict.hpp"
using pFlow::commandLine;
using namespace pFlow;
int main( int argc, char* argv[] )
{
@ -54,6 +55,8 @@ commandLine cmds(
// this should be palced in each main
processors::initProcessors(argc, argv);
initialize_pFlowProcessors();
#include "initialize_Control.hpp"
#include "setProperty.hpp"
@ -68,6 +71,7 @@ commandLine cmds(
// this should be palced in each main
#include "finalize.hpp"
processors::finalizeProcessors();
}

7
solvers/iterateSphereParticles/CMakeLists.txt Normal file
View File

@ -0,0 +1,7 @@
set(source_files
iterateSphereParticles.cpp
)
set(link_lib Particles)
pFlow_make_executable_install(iterateSphereParticles source_files link_lib)

10
src/phasicFlow/structuredData/peakableRegion/peakableRegions.cpp → solvers/iterateSphereParticles/createDEMComponents.hpp Normal file → Executable file
View File

@ -18,10 +18,8 @@ Licence:
-----------------------------------------------------------------------------*/
//
REPORT(0)<<"\nReading sphere particles . . ."<<END_REPORT;
sphereParticles sphParticles(Control, proprties);
#include "peakableRegions.hpp"
#ifdef BUILD_SHARED_LIBS
#include "peakableRegionInstantiate.cpp"
#endif
WARNING<<"Particle insertion has not been set yet!"<<END_WARNING;

89
solvers/iterateSphereParticles/iterateSphereParticles.cpp Executable file
View File

@ -0,0 +1,89 @@
/*------------------------------- 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.
-----------------------------------------------------------------------------*/
/**
* \file sphereGranFlow.cpp
* \brief sphereGranFlow solver
*
* This solver simulate granular flow of cohesion-less, spherical particles.
* Particle insertion can be activated in the solver to gradually insert
* particles into the simulation. Geometry can be defined generally with
* built-in features of the code or through ASCII stl files or a combination
* of both. For more information refer to [tutorials/sphereGranFlow/]
* (https://github.com/PhasicFlow/phasicFlow/tree/main/tutorials/sphereGranFlow) folder.
*/
#include "vocabs.hpp"
#include "property.hpp"
#include "sphereParticles.hpp"
#include "systemControl.hpp"
#include "commandLine.hpp"
using namespace pFlow;
/**
* DEM solver for simulating granular flow of cohesion-less particles.
*
* In the root case directory just simply enter the following command to
* run the simulation. For command line options use flag -h.
*/
int main( int argc, char* argv[])
{
commandLine cmds(
"sphereGranFlow",
"DEM solver for non-cohesive spherical particles with particle insertion "
"mechanism and moving geometry");
bool isCoupling = false;
if(!cmds.parse(argc, argv)) return 0;
// this should be palced in each main
processors::initProcessors(argc, argv);
initialize_pFlowProcessors();
#include "initialize_Control.hpp"
#include "setProperty.hpp"
#include "createDEMComponents.hpp"
REPORT(0)<<"\nStart of time loop . . .\n"<<END_REPORT;
do
{
sphParticles.beforeIteration();
sphParticles.iterate();
sphParticles.afterIteration();
}while(Control++);
REPORT(0)<<"\nEnd of time loop.\n"<<END_REPORT;
// this should be palced in each main
#include "finalize.hpp"
processors::finalizeProcessors();
}

9
solvers/sphereGranFlow/createDEMComponents.hpp
View File

@ -19,11 +19,11 @@ Licence:
-----------------------------------------------------------------------------*/
//
REPORT(0)<<"\nReading sphere particles . . ."<<endREPORT;
REPORT(0)<<"\nReading sphere particles . . ."<<END_REPORT;
sphereParticles sphParticles(Control, proprties);
//
REPORT(0)<<"\nCreating particle insertion object . . ."<<endREPORT;
REPORT(0)<<"\nCreating particle insertion object . . ."<<END_REPORT;
/*auto& sphInsertion =
Control.caseSetup().emplaceObject<sphereInsertion>(
objectFile(
@ -37,11 +37,10 @@ REPORT(0)<<"\nCreating particle insertion object . . ."<<endREPORT;
);*/
auto sphInsertion = sphereInsertion(
Control.caseSetup().path()+insertionFile__,
sphParticles,
sphParticles.shapes());
sphParticles.spheres());
REPORT(0)<<"\nCreating interaction model for sphere-sphere contact and sphere-wall contact . . ."<<endREPORT;
REPORT(0)<<"\nCreating interaction model for sphere-sphere contact and sphere-wall contact . . ."<<END_REPORT;
auto interactionPtr = interaction::create(
Control,
sphParticles,

49
solvers/sphereGranFlow/sphereGranFlow.cpp
View File

@ -30,30 +30,22 @@ Licence:
* (https://github.com/PhasicFlow/phasicFlow/tree/main/tutorials/sphereGranFlow) folder.
*/
#include "vocabs.hpp"
#include "phasicFlowKokkos.hpp"
#include "systemControl.hpp"
#include "commandLine.hpp"
#include "property.hpp"
#include "geometry.hpp"
#include "sphereParticles.hpp"
#include "interaction.hpp"
#include "Insertions.hpp"
#include "systemControl.hpp"
#include "contactSearch.hpp"
#include "sphereInteraction.hpp"
#include "commandLine.hpp"
#include "readControlDict.hpp"
using pFlow::output;
using pFlow::endl;
using pFlow::property;
using pFlow::sphereParticles;
using pFlow::objectFile;
using pFlow::sphereInsertion;
using pFlow::insertionFile__;
using pFlow::interactionFile__;
using pFlow::contactSearch;
using pFlow::interaction;
using pFlow::commandLine;
//#include "readControlDict.hpp"
using namespace pFlow;
/**
* DEM solver for simulating granular flow of cohesion-less particles.
*
@ -73,6 +65,8 @@ bool isCoupling = false;
if(!cmds.parse(argc, argv)) return 0;
// this should be palced in each main
processors::initProcessors(argc, argv);
initialize_pFlowProcessors();
#include "initialize_Control.hpp"
#include "setProperty.hpp"
@ -81,12 +75,13 @@ if(!cmds.parse(argc, argv)) return 0;
#include "createDEMComponents.hpp"
REPORT(0)<<"\nStart of time loop . . .\n"<<endREPORT;
REPORT(0)<<"\nStart of time loop . . .\n"<<END_REPORT;
do
{
if(! sphInsertion.insertParticles(
Control.time().currentIter(),
Control.time().currentTime(),
Control.time().dt() ) )
{
@ -95,30 +90,34 @@ if(!cmds.parse(argc, argv)) return 0;
return 1;
}
// set force to zero
surfGeometry.beforeIteration();
// set force to zero, predict, particle deletion and etc.
sphParticles.beforeIteration();
sphInteraction.beforeIteration();
sphParticles.beforeIteration();
sphInteraction.iterate();
sphParticles.iterate();
surfGeometry.iterate();
sphParticles.afterIteration();
sphParticles.iterate();
sphInteraction.afterIteration();
surfGeometry.afterIteration();
sphParticles.afterIteration();
}while(Control++);
REPORT(0)<<"\nEnd of time loop.\n"<<endREPORT;
REPORT(0)<<"\nEnd of time loop.\n"<<END_REPORT;
// this should be palced in each main
#include "finalize.hpp"
processors::finalizeProcessors();
}

14
src/CMakeLists.txt
View File

@ -1,17 +1,15 @@
add_subdirectory(phasicFlow)
#add_subdirectory(Integration)
add_subdirectory(Integration)
#add_subdirectory(Property)
add_subdirectory(Property)
#add_subdirectory(Particles)
add_subdirectory(Particles)
#add_subdirectory(Interaction)
add_subdirectory(Geometry)
#add_subdirectory(MotionModel)
add_subdirectory(Interaction)
#add_subdirectory(Geometry)
add_subdirectory(MotionModel)
#add_subdirectory(MPIParallelization)

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

@ -19,32 +19,40 @@ Licence:
-----------------------------------------------------------------------------*/
#include "geometry.hpp"
#include "systemControl.hpp"
#include "vocabs.hpp"
bool pFlow::geometry::findPropertyId()
bool pFlow::geometry::createPropertyId()
{
int8Vector propId(0, surface().capacity(),RESERVE());
propId.clear();
uint32 pId;
ForAll(matI, materialName_)
{
if( !wallProperty_.nameToIndex( materialName_[matI], pId ) )
if(materialName_.size() != numSurfaces() )
{
fatalErrorInFunction<<
"material name for the geometry is invalid: "<< materialName_[matI]<<endl;
"number of subSurface and material names do not match"<<endl;
return false;
}
int32 surfSize = surface().surfNumTriangles(matI);
uint32Vector propId(
"propId",
capacity(),
size(),
RESERVE());
for(int32 i=0; i<surfSize; i++)
ForAll(i, materialName_)
{
propId.push_back(pId);
uint32 pIdx =0;
if( !wallProperty_.nameToIndex(materialName_[i], pIdx) )
{
fatalErrorInFunction<<
"Property/material name is invalid "<<materialName_[i]<<
". A list of valid names are \n"<< wallProperty_.materials()<<endl;
return false;
}
auto triRange = subSurfaceRange(i);
propId.fill(triRange.start(), triRange.end(), pIdx);
}
propertyId_.assign(propId);
@ -53,83 +61,361 @@ bool pFlow::geometry::findPropertyId()
}
void pFlow::geometry::zeroForce()
{
contactForceWall_.fill(zero3);
}
pFlow::geometry::geometry
(
systemControl& control,
const property& prop
)
:
demGeometry(control),
multiTriSurface
(
objectFile
(
triSurfaceFile__,
"",
objectFile::READ_ALWAYS,
objectFile::WRITE_ALWAYS
),
&control.geometry()
),
demComponent
(
"geometry",
control
),
wallProperty_(prop),
geometryRepository_(control.geometry()),
triSurface_(
control.geometry().emplaceObject<multiTriSurface>(
objectFile(
"triSurface",
"",
objectFile::READ_ALWAYS,
objectFile::WRITE_ALWAYS
)
)
),
motionComponentName_(
control.geometry().emplaceObject<wordField>(
objectFile(
"motionComponentName",
"",
objectFile::READ_ALWAYS,
objectFile::WRITE_ALWAYS
),
"motionNamesList"
)
),
materialName_(
control.geometry().emplaceObject<wordField>(
objectFile(
"materialName",
"",
objectFile::READ_ALWAYS,
objectFile::WRITE_ALWAYS
),
"materialNamesList"
)
),
propertyId_(
control.geometry().emplaceObject<int8TriSurfaceField_D>(
objectFile(
propertyId_
(
objectFile
(
"propertyId",
"",
objectFile::READ_NEVER,
objectFile::WRITE_NEVER),
surface(),
0 ) ),
contactForceWall_(
control.geometry().emplaceObject<realx3TriSurfaceField_D>(
objectFile(
"contactForceWall",
objectFile::WRITE_NEVER
),
*this,
0u
),
contactForceWall_
(
objectFile
(
"contactForcWall",
"",
objectFile::READ_IF_PRESENT,
objectFile::WRITE_ALWAYS),
surface(),
zero3) ),
stressWall_(
control.geometry().emplaceObject<realx3TriSurfaceField_D>(
objectFile(
"stressWall",
objectFile::WRITE_NEVER
),
*this,
zero3
),
normalStressWall_
(
objectFile
(
"normalStressWall",
"",
objectFile::READ_IF_PRESENT,
objectFile::WRITE_ALWAYS),
surface(),
zero3) )
objectFile::WRITE_NEVER
),
*this,
zero3
),
shearStressWall_
(
objectFile
(
"shearStressWall",
"",
objectFile::READ_IF_PRESENT,
objectFile::WRITE_NEVER
),
*this,
zero3
)
{
if(!findPropertyId())
readWholeObject_ = false;
if( !IOobject::readObject() )
{
fatalErrorInFunction<<
"Error in reading from file "<<IOobject::path()<<endl;
fatalExit;
}
readWholeObject_ = true;
if( this->numSurfaces() != motionComponentName_.size() )
{
fatalErrorInFunction<<
"Number of surfaces is not equal to number of motion component names"<<endl;
fatalExit;
}
if(!createPropertyId())
{
fatalExit;
}
}
pFlow::geometry::geometry
(
systemControl &control,
const property &prop,
multiTriSurface &surf,
const wordVector &motionCompName,
const wordVector &materialName,
const dictionary& motionDict
)
:
multiTriSurface
(
objectFile
(
triSurfaceFile__,
"",
objectFile::READ_NEVER,
objectFile::WRITE_ALWAYS
),
&control.geometry(),
surf
),
demComponent
(
"geometry",
control
),
wallProperty_
(
prop
),
propertyId_
(
objectFile
(
"propertyId",
"",
objectFile::READ_NEVER,
objectFile::WRITE_NEVER
),
*this,
0u
),
contactForceWall_
(
objectFile
(
"contactForcWall",
"",
objectFile::READ_IF_PRESENT,
objectFile::WRITE_NEVER
),
*this,
zero3
),
normalStressWall_
(
objectFile
(
"normalStressWall",
"",
objectFile::READ_IF_PRESENT,
objectFile::WRITE_NEVER
),
*this,
zero3
),
shearStressWall_
(
objectFile
(
"shearStressWall",
"",
objectFile::READ_IF_PRESENT,
objectFile::WRITE_NEVER
),
*this,
zero3
)
{
motionComponentName_.assign(motionCompName);
materialName_.assign(materialName);
if( this->numSurfaces() != motionComponentName_.size() )
{
fatalErrorInFunction<<
"Number of surfaces is not equal to number of motion component names"<<endl;
fatalExit;
}
if(!createPropertyId())
{
fatalExit;
}
}
bool pFlow::geometry::beforeIteration()
{
zeroForce();
return true;
}
bool pFlow::geometry::iterate()
{
return true;
}
bool pFlow::geometry::afterIteration()
{
auto numTris = size();
auto& normalsD = normals().deviceViewAll();
auto& areaD = area().deviceViewAll();
auto& cForceD = contactForceWall_.deviceViewAll();
auto& sStressD = shearStressWall_.deviceViewAll();
auto& nStressD = normalStressWall_.deviceViewAll();
Kokkos::parallel_for(
"pFlow::geometry::afterIteration",
deviceRPolicyStatic(0, numTris),
LAMBDA_HD(uint32 i)
{
realx3 n = normalsD[i];
real A = max(areaD[i],smallValue);
realx3 nF = dot(cForceD[i],n)*n;
realx3 sF = cForceD[i] - nF;
nStressD[i] = nF/A;
sStressD[i] = sF/A;
});
Kokkos::fence();
return true;
}
bool pFlow::geometry::read(iIstream &is, const IOPattern &iop)
{
motionComponentName_.read(is, iop);
materialName_.read(is, iop);
if( readWholeObject_ )
{
return multiTriSurface::read(is, iop);
}
return true;
}
bool pFlow::geometry::write(iOstream &os, const IOPattern &iop) const
{
if( !motionComponentName_.write(os, iop) )
{
fatalErrorInFunction;
return false;
}
if( !materialName_.write(os,iop))
{
fatalErrorInFunction;
return false;
}
return multiTriSurface::write(os,iop);
}
pFlow::uniquePtr<pFlow::geometry>
pFlow::geometry::create
(
systemControl& control,
const property& prop
)
{
//
fileDictionary dict( motionModelFile__, control.time().geometry().path());
word model = dict.getVal<word>("motionModel");
auto geomModel = angleBracketsNames("geometry", model);
REPORT(1)<< "Selecting geometry model . . ."<<END_REPORT;
if( systemControlvCtorSelector_.search(geomModel) )
{
auto objPtr = systemControlvCtorSelector_[geomModel] (control, prop);
REPORT(2)<<"Model "<< Green_Text(geomModel)<<" is created.\n"<<END_REPORT;
return objPtr;
}
else
{
printKeys
(
fatalError << "Ctor Selector "<< Yellow_Text(geomModel) << " dose not exist. \n"
<<"Avaiable ones are: \n\n"
,
systemControlvCtorSelector_
);
fatalExit;
}
return nullptr;
}
pFlow::uniquePtr<pFlow::geometry>
pFlow::geometry::create
(
systemControl& control,
const property& prop,
multiTriSurface& surf,
const wordVector& motionCompName,
const wordVector& materialName,
const dictionary& motionDic
)
{
word model = motionDic.getVal<word>("motionModel");
auto geomModel = angleBracketsNames("geometry", model);
REPORT(1)<< "Selecting geometry model . . ."<<END_REPORT;
if( dictionaryvCtorSelector_.search(geomModel) )
{
auto objPtr = dictionaryvCtorSelector_[geomModel]
(
control,
prop,
surf,
motionCompName,
materialName,
motionDic
);
REPORT(2)<<"Model "<< Green_Text(geomModel)<<" is created.\n"<<END_REPORT;
return objPtr;
}
else
{
printKeys
(
fatalError << "Ctor Selector "<< Yellow_Text(geomModel) << " dose not exist. \n"
<<"Avaiable ones are: \n\n"
,
dictionaryvCtorSelector_
);
fatalExit;
}
return nullptr;
}
/*pFlow::geometry::geometry
(
systemControl& control,
const property& prop,
@ -224,52 +510,7 @@ pFlow::geometry::geometry
{}
pFlow::uniquePtr<pFlow::geometry>
pFlow::geometry::create
(
systemControl& control,
const property& prop
)
{
//motionModelFile__
auto motionDictPtr = IOobject::make<dictionary>
(
objectFile
(
motionModelFile__,
control.geometry().path(),
objectFile::READ_ALWAYS,
objectFile::WRITE_NEVER
),
motionModelFile__,
true
);
word model = motionDictPtr().getObject<dictionary>().getVal<word>("motionModel");
auto geomModel = angleBracketsNames("geometry", model);
REPORT(1)<< "Selecting geometry model . . ."<<endREPORT;
if( systemControlvCtorSelector_.search(geomModel) )
{
auto objPtr = systemControlvCtorSelector_[geomModel] (control, prop);
REPORT(2)<<"Model "<< greenText(geomModel)<<" is created.\n"<<endREPORT;
return objPtr;
}
else
{
printKeys
(
fatalError << "Ctor Selector "<< yellowText(geomModel) << " dose not exist. \n"
<<"Avaiable ones are: \n\n"
,
systemControlvCtorSelector_
);
fatalExit;
}
return nullptr;
}
bool pFlow::geometry::beforeIteration()
{
@ -296,48 +537,4 @@ bool pFlow::geometry::afterIteration()
Kokkos::fence();
return true;
}
pFlow::uniquePtr<pFlow::geometry>
pFlow::geometry::create(
systemControl& control,
const property& prop,
const dictionary& dict,
const multiTriSurface& triSurface,
const wordVector& motionCompName,
const wordVector& propName)
{
word model = dict.getVal<word>("motionModel");
auto geomModel = angleBracketsNames("geometry", model);
REPORT(1)<< "Selecting geometry model . . ."<<endREPORT;
if( dictionaryvCtorSelector_.search(geomModel) )
{
auto objPtr = dictionaryvCtorSelector_[geomModel]
(
control,
prop,
dict,
triSurface,
motionCompName,
propName
);
REPORT(2)<<"Model "<< greenText(geomModel)<<" is created.\n"<<endREPORT;
return objPtr;
}
else
{
printKeys
(
fatalError << "Ctor Selector "<< yellowText(geomModel) << " dose not exist. \n"
<<"Avaiable ones are: \n\n"
,
dictionaryvCtorSelector_
);
fatalExit;
}
return nullptr;
}
}*/

191
src/Geometry/geometry/geometry.hpp
View File

@ -23,13 +23,14 @@ Licence:
#include "virtualConstructor.hpp"
#include "demGeometry.hpp"
#include "demComponent.hpp"
#include "property.hpp"
#include "Fields.hpp"
#include "Vectors.hpp"
#include "multiTriSurface.hpp"
#include "triSurfaceFields.hpp"
#include "dictionary.hpp"
//#include "Fields.hpp"
//#include "Vectors.hpp"
namespace pFlow
{
@ -42,48 +43,49 @@ namespace pFlow
*/
class geometry
:
public demGeometry
public multiTriSurface,
public demComponent
{
protected:
private:
// - Protected members
/// Const reference to physical property of materials
const property& wallProperty_;
/// Repository to store geometry data at each simulation moment
repository& geometryRepository_;
/// All triangles in the set of wall surfaces
multiTriSurface& triSurface_;
/// The name of motion component of each wall surface
wordField& motionComponentName_;
wordField_H motionComponentName_{
"motionComponentName",
"motionComponentName"};
/// Material name of each wall surface
wordField& materialName_;
wordField_H materialName_{
"materialName",
"materialName"};
/// Property id of each triangle in the set of wall surfaces
int8TriSurfaceField_D& propertyId_;
uint32TriSurfaceField_D propertyId_;
/// Contact force on each triangle in the set of wall surfaces
realx3TriSurfaceField_D& contactForceWall_;
realx3TriSurfaceField_D contactForceWall_;
/// Stress on ech triangle in the set of wall surfaces
realx3TriSurfaceField_D& stressWall_;
/// Stress on each triangle in the set of wall surfaces
realx3TriSurfaceField_D normalStressWall_;
/// Stress on each triangle in the set of wall surfaces
realx3TriSurfaceField_D shearStressWall_;
bool readWholeObject_ = true;
// - Protected member functions
/// Find property id of each triangle based on the supplied material name
/// and the surface wall that the triangle belongs to.
bool findPropertyId();
bool createPropertyId();
/// Initialize contact force to zero
void zeroForce()
{
contactForceWall_.fill(zero3);
}
void zeroForce();
public:
@ -95,8 +97,15 @@ public:
/// Construct from controlSystem and property, for reading from file
geometry(systemControl& control, const property& prop);
/// Construct from components
geometry(systemControl& control,
const property& prop,
multiTriSurface& surf,
const wordVector& motionCompName,
const wordVector& materialName,
const dictionary& motionDict);
/// Construct from components
/*geometry(systemControl& control,
const property& prop,
const multiTriSurface& triSurface,
const wordVector& motionCompName,
@ -110,10 +119,10 @@ public:
const dictionary& dict,
const multiTriSurface& triSurface,
const wordVector& motionCompName,
const wordVector& propName);
const wordVector& propName);*/
/// Destructor
virtual ~geometry() = default;
~geometry()override = default;
/// Virtual constructor
create_vCtor
@ -132,156 +141,96 @@ public:
(
geometry,
dictionary,
(systemControl& control,
(
systemControl& control,
const property& prop,
const dictionary& dict,
const multiTriSurface& triSurface,
multiTriSurface& surf,
const wordVector& motionCompName,
const wordVector& propName),
(control, prop, dict, triSurface, motionCompName, propName)
const wordVector& materialName,
const dictionary& motionDic),
(control, prop, surf, motionCompName, materialName, motionDic)
);
//- Methods
/// Size of tri-surface
inline
auto size()const
const auto& motionComponentName()const
{
return triSurface_.size();
}
/// Number of points in the set of surface walls
inline
auto numPoints()const
{
return triSurface_.numPoints();
}
/// Number of triangles in the set of surface walls
inline
auto numTriangles()const
{
return size();
}
/// Access to the points
inline
const auto& points()const
{
return triSurface_.points();
}
/// Access to the vertices
inline
const auto& vertices()const
{
return triSurface_.vertices();
}
/// Obtain an object for accessing triangles
inline auto getTriangleAccessor()const
{
return triSurface_.getTriangleAccessor();
}
/// Surface
inline auto& surface()
{
return triSurface_;
}
/// Surface
inline const auto& surface()const
{
return triSurface_;
return motionComponentName_;
}
/// Access to contact force
inline
realx3TriSurfaceField_D& contactForceWall()
auto& contactForceWall()
{
return contactForceWall_;
}
/// Access to contact force
inline
const realx3TriSurfaceField_D& contactForceWall() const
const auto& contactForceWall() const
{
return contactForceWall_;
}
/// Property ide of triangles
inline
const auto& propertyId()const
{
return propertyId_;
}
/// Access to property
inline const auto& wallProperty()const
{
return wallProperty_;
}
/// Owner repository
inline
const repository& owner()const
{
return geometryRepository_;
}
/// Owner repository
inline
repository& owner()
{
return geometryRepository_;
}
/// Path to the repository folder
inline auto path()
{
return owner().path();
}
/// The name of motion model
virtual
word motionModelTypeName()const = 0;
/// Motion model index of triangles
virtual
const int8Vector_HD& triMotionIndex() const =0;
const uint32Field_D& triMotionIndex() const =0;
/// Motion model index of points
virtual
const int8Vector_HD& pointMotionIndex()const = 0;
const uint32Field_D& pointMotionIndex()const = 0;
/// Property ide of triangles
const int8TriSurfaceField_D& propertyId() const
{
return propertyId_;
}
/// Operations before each iteration
bool beforeIteration() override;
/// Operations after each iteration
/// This is called in time loop. Perform the main calculations
/// when the component should evolve along time.
bool iterate() override;
/// This is called in time loop, after iterate.
bool afterIteration() override;
//- IO
bool read(iIstream& is, const IOPattern& iop) override;
/// write
bool write()const
{
return owner().write();
}
bool write( iOstream& os, const IOPattern& iop )const override;
//- Static members
static
uniquePtr<geometry> create(systemControl& control, const property& prop);
uniquePtr<geometry> create(
systemControl& control,
const property& prop);
static
uniquePtr<geometry> create(
systemControl& control,
const property& prop,
const dictionary& dict,
const multiTriSurface& triSurface,
multiTriSurface& surf,
const wordVector& motionCompName,
const wordVector& propName);
const wordVector& materialName,
const dictionary& motionDic);
};

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

@ -1,3 +1,4 @@
#include "geometryMotion.hpp"
/*------------------------------- phasicFlow ---------------------------------
O C enter of
O O E ngineering and
@ -21,42 +22,88 @@ Licence:
template<typename MotionModel>
bool pFlow::geometryMotion<MotionModel>::findMotionIndex()
{
motionIndex_.clear();
triMotionIndex_.reserve( this->surface().capacity() );
triMotionIndex_.clear();
ForAll( surfI, motionComponentName_)
if(motionComponentName().size() != numSurfaces() )
{
auto mName = motionComponentName_[surfI];
auto mInd = motionModel_.nameToIndex(mName);
motionIndex_.push_back(mInd);
// fill motionIndex for triangles of the surface
int32 surfSize = this->surface().surfNumTriangles(surfI);
for(int32 i=0; i<surfSize; i++)
{
triMotionIndex_.push_back(mInd);
fatalErrorInFunction<<
"size of motion component names in the triSurface is not"<<
" equal to size of number of sub-surfaces"<<endl;
return false;
}
}
motionIndex_.syncViews();
triMotionIndex_.syncViews();
pointMotionIndex_.reserve(triSurface_.numPoints());
pointMotionIndex_.clear();
uint32Vector surfMotionIndex("surfMotionIndex");
uint32Vector triMotionIndex("triMotionIndex");
uint32Vector pointMotionIndex("pointMotionIndex");
ForAll(surfI, motionIndex_)
ForAll( surfI, motionComponentName())
{
auto nP = triSurface_.surfNumPoints(surfI);
for(int32 i=0; i<nP; i++)
auto mName = motionComponentName()[surfI];
uint32 mInd=0;
if( !motionModel_.nameToIndex(mName, mInd) )
{
pointMotionIndex_.push_back(motionIndex_[surfI]);
fatalErrorInFunction<<
mName<< " does not exist in the list of motion names -> "<<
motionModel_.componentNames();
return false;
}
surfMotionIndex.push_back(mInd);
auto surfRange = subSurfaceRange(surfI);
for(uint32 i=0; i<surfRange.numElements(); i++)
{
triMotionIndex.push_back(mInd);
}
auto pointRange = subSurfacePointRange(surfI);
for(uint32 n=0; n<pointRange.numElements(); n++)
{
pointMotionIndex.push_back(mInd);
}
}
pointMotionIndex_.syncViews();
surfMotionIndex_.assign(surfMotionIndex);
triMotionIndex_.assign(triMotionIndex);
pointMotionIndex_.assign(pointMotionIndex);
return true;
}
template<typename MotionModel>
bool pFlow::geometryMotion<MotionModel>::moveGeometry()
{
uint32 iter = this->currentIter();
real t = this->currentTime();
real dt = this->dt();
auto mModel = motionModel_.getModelInterface(iter, t, dt);
auto& pointMIndexD= pointMotionIndex_.deviceViewAll();
auto& pointsD = points().deviceViewAll();
Kokkos::parallel_for(
"geometryMotion<MotionModel>::movePoints",
deviceRPolicyStatic(0, numPoints()),
LAMBDA_HD(uint32 i){
auto newPos = mModel.transferPoint(pointMIndexD[i], pointsD[i], dt);
pointsD[i] = newPos;
});
Kokkos::fence();
// move the motion components
motionModel_.move(iter, t,dt);
// end of calculations
return true;
}
template<typename MotionModel>
pFlow::geometryMotion<MotionModel>::geometryMotion
(
@ -65,96 +112,65 @@ pFlow::geometryMotion<MotionModel>::geometryMotion
)
:
geometry(control, prop),
motionModel_(
this->owner().template emplaceObject<MotionModel>(
objectFile(
motionModel_
(
objectFile
(
motionModelFile__,
"",
objectFile::READ_ALWAYS,
objectFile::WRITE_ALWAYS
)
)
),
owner()
),
moveGeomTimer_("move geometry", &this->timers())
{
findMotionIndex();
if(!findMotionIndex())
{
fatalExit;
}
}
template<typename MotionModel>
pFlow::geometryMotion<MotionModel>::geometryMotion
template <typename MotionModelType>
pFlow::geometryMotion<MotionModelType>::geometryMotion
(
systemControl &control,
const property &prop,
const multiTriSurface& triSurface,
multiTriSurface &surf,
const wordVector &motionCompName,
const wordVector& propName,
const MotionModel& motionModel
const wordVector &materialName,
const dictionary &motionDict
)
:
geometry(
geometry
(
control,
prop,
triSurface,
surf,
motionCompName,
propName
materialName,
motionDict
),
motionModel_(
this->owner().template emplaceObject<MotionModel>(
objectFile(
motionModel_
(
objectFile
(
motionModelFile__,
"",
objectFile::READ_NEVER,
objectFile::WRITE_ALWAYS
),
motionModel
)
motionDict,
owner()
),
moveGeomTimer_("move geometry", &this->timers())
{
findMotionIndex();
if(!findMotionIndex())
{
fatalExit;
}
}
template<typename MotionModel>
pFlow::geometryMotion<MotionModel>::geometryMotion
(
systemControl& control,
const property& prop,
const dictionary& dict,
const multiTriSurface& triSurface,
const wordVector& motionCompName,
const wordVector& propName
)
:
geometry(
control,
prop,
dict,
triSurface,
motionCompName,
propName
),
motionModel_(
this->owner().template emplaceObject<MotionModel>(
objectFile(
motionModelFile__,
"",
objectFile::READ_NEVER,
objectFile::WRITE_ALWAYS
),
dict
)
),
moveGeomTimer_("move geometry", &this->timers())
{
findMotionIndex();
}
template<typename MotionModel>
bool pFlow::geometryMotion<MotionModel>::beforeIteration()
{
geometry::beforeIteration();
return true;
}
template<typename MotionModel>
bool pFlow::geometryMotion<MotionModel>::iterate()
@ -163,46 +179,8 @@ bool pFlow::geometryMotion<MotionModel>::iterate()
{
moveGeomTimer_.start();
moveGeometry();
this->calculateNormals();
moveGeomTimer_.end();
}
return true;
}
template<typename MotionModel>
bool pFlow::geometryMotion<MotionModel>::afterIteration()
{
geometry::afterIteration();
return true;
}
template<typename MotionModel>
bool pFlow::geometryMotion<MotionModel>::moveGeometry()
{
real dt = this->dt();
real t = this->currentTime();
auto pointMIndex= pointMotionIndex_.deviceVector();
auto mModel = motionModel_.getModel(t);
realx3* points = triSurface_.pointsData_D();
auto numPoints = triSurface_.numPoints();
Kokkos::parallel_for(
"geometryMotion<MotionModel>::movePoints",
numPoints,
LAMBDA_HD(int32 i){
auto newPos = mModel.transferPoint(pointMIndex[i], points[i], dt);
points[i] = newPos;
});
Kokkos::fence();
// move the motion components
motionModel_.move(t,dt);
// end of calculations
moveGeomTimer_.end();
return true;
}

54
src/Geometry/geometryMotion/geometryMotion.hpp
View File

@ -20,9 +20,8 @@ Licence:
#ifndef __geometryMotion_hpp__
#define __geometryMotion_hpp__
#include "vocabs.hpp"
#include "geometry.hpp"
#include "VectorDuals.hpp"
namespace pFlow
{
@ -41,19 +40,21 @@ public:
using MotionModel = MotionModelType;
protected:
using ModelComponent = typename MotionModelType::ModelComponent;
private:
/// Ref to motion model
MotionModel& motionModel_;
MotionModelType motionModel_;
/// motion indext mapped on each surface
int32Vector_HD motionIndex_;
uint32Field_D surfMotionIndex_{"triMotionIndex"};
/// motion index mapped on each triangle
int8Vector_HD triMotionIndex_;
uint32Field_D triMotionIndex_ {"surfMotionIndex"};
/// motion index mapped on each point
int8Vector_HD pointMotionIndex_;
uint32Field_D pointMotionIndex_{"pointMotionIndex"};
/// timer for moveGeometry
Timer moveGeomTimer_;
@ -61,32 +62,25 @@ protected:
/// determine the motion index of each triangle
bool findMotionIndex();
/// Move geometry
bool moveGeometry();
public:
/// Type info
TypeInfoTemplate("geometry", MotionModel);
TypeInfoTemplate11("geometry", ModelComponent);
// - Constructors
geometryMotion(systemControl& control, const property& prop);
geometryMotion(
systemControl& control,
const property& prop,
const multiTriSurface& triSurface,
multiTriSurface& surf,
const wordVector& motionCompName,
const wordVector& propName,
const MotionModel& motionModel);
/// Construct from components and dictionary that contains
/// motionModel
geometryMotion(systemControl& control,
const property& prop,
const dictionary& dict,
const multiTriSurface& triSurface,
const wordVector& motionCompName,
const wordVector& propName);
const wordVector& materialName,
const dictionary& motionDict);
/// Add virtual constructor
add_vCtor
@ -107,9 +101,9 @@ public:
// - Methods
/// Obtain motion model at time t
auto getModel(real t)const
auto getModel(uint32 iter, real t, real dt)const
{
return motionModel_.getModel(t);
return motionModel_.getModelInterface(iter, t, dt);
}
/// TypeName / TypeInfo of motion model
@ -119,28 +113,21 @@ public:
}
/// Access to motion model index of triangles
const int8Vector_HD& triMotionIndex()const override
const uint32Field_D& triMotionIndex()const override
{
return triMotionIndex_;
}
/// Access to motion model index of points
const int8Vector_HD& pointMotionIndex()const override
const uint32Field_D& pointMotionIndex()const override
{
return pointMotionIndex_;
}
/// Operations before each iteration
bool beforeIteration() override;
/// Iterate geometry one time step
bool iterate() override ;
/// Operations after each iteration
bool afterIteration() override;
/// Move geometry
bool moveGeometry();
};
@ -148,9 +135,6 @@ public:
#include "geometryMotion.cpp"
#ifndef BUILD_SHARED_LIBS
#include "geometryMotionsInstantiate.cpp"
#endif
#endif //__geometryMotion_hpp__

12
src/Geometry/geometryMotion/geometryMotions.cpp
View File

@ -17,9 +17,13 @@ Licence:
implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
-----------------------------------------------------------------------------*/
#include "geometryMotions.hpp"
#ifdef BUILD_SHARED_LIBS
#include "geometryMotionsInstantiate.cpp"
#endif
template class pFlow::geometryMotion<pFlow::vibratingMotion>;
template class pFlow::geometryMotion<pFlow::rotatingAxisMotion>;
template class pFlow::geometryMotion<pFlow::stationaryWall>;
//template class pFlow::geometryMotion<pFlow::multiRotatingAxisMotion>;

13
src/Geometry/geometryMotion/geometryMotions.hpp
View File

@ -22,22 +22,23 @@ Licence:
#define __geometryMotions_hpp__
#include "geometryMotion.hpp"
#include "fixedWall.hpp"
#include "stationaryWall.hpp"
#include "rotatingAxisMotion.hpp"
#include "multiRotatingAxisMotion.hpp"
//#include "multiRotatingAxisMotion.hpp"
#include "vibratingMotion.hpp"
namespace pFlow
{
typedef geometryMotion<vibratingMotion> vibratingMotionGeometry;
using vibratingMotionGeometry = geometryMotion<vibratingMotion>;
typedef geometryMotion<rotatingAxisMotion> rotationAxisMotionGeometry;
using rotationAxisMotionGeometry = geometryMotion<rotatingAxisMotion>;
typedef geometryMotion<multiRotatingAxisMotion> multiRotationAxisMotionGeometry;
using stationaryGeometry = geometryMotion<stationaryWall>;
//typedef geometryMotion<multiRotatingAxisMotion> multiRotationAxisMotionGeometry;
typedef geometryMotion<fixedWall> fixedGeometry;

32
src/Geometry/geometryMotion/geometryMotionsInstantiate.cpp
View File

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

155
src/Integration/AdamsBashforth2/AdamsBashforth2.cpp
View File

@ -19,59 +19,146 @@ Licence:
-----------------------------------------------------------------------------*/
#include "AdamsBashforth2.hpp"
#include "pointStructure.hpp"
#include "Time.hpp"
#include "vocabs.hpp"
namespace pFlow
{
/// Range policy for integration kernel (alias)
using rpIntegration = Kokkos::RangePolicy<
DefaultExecutionSpace,
Kokkos::Schedule<Kokkos::Static>,
Kokkos::IndexType<uint32>
>;
bool intAllActive(
real dt,
realx3Field_D& y,
realx3PointField_D& dy,
realx3PointField_D& dy1)
{
auto d_dy = dy.deviceView();
auto d_y = y.deviceView();
auto d_dy1= dy1.deviceView();
auto activeRng = dy1.activeRange();
Kokkos::parallel_for(
"AdamsBashforth2::correct",
rpIntegration (activeRng.start(), activeRng.end()),
LAMBDA_HD(uint32 i){
d_y[i] += dt*(static_cast<real>(1.5) * d_dy[i] - static_cast<real>(0.5) * d_dy1[i]);
d_dy1[i] = d_dy[i];
});
Kokkos::fence();
return true;
}
bool intScattered
(
real dt,
realx3Field_D& y,
realx3PointField_D& dy,
realx3PointField_D& dy1
)
{
auto d_dy = dy.deviceView();
auto d_y = y.deviceView();
auto d_dy1 = dy1.deviceView();
auto activeRng = dy1.activeRange();
const auto& activeP = dy1.activePointsMaskDevice();
Kokkos::parallel_for(
"AdamsBashforth2::correct",
rpIntegration (activeRng.start(), activeRng.end()),
LAMBDA_HD(uint32 i){
if( activeP(i))
{
d_y[i] += dt*(static_cast<real>(1.5) * d_dy[i] - static_cast<real>(0.5) * d_dy1[i]);
d_dy1[i] = d_dy[i];
}
});
Kokkos::fence();
return true;
}
}
//const real AB2_coef[] = { 3.0 / 2.0, 1.0 / 2.0};
pFlow::AdamsBashforth2::AdamsBashforth2
(
const word& baseName,
repository& owner,
const pointStructure& pStruct,
const word& method
pointStructure& pStruct,
const word& method,
const realx3Field_D& initialValField
)
:
integration(baseName, owner, pStruct, method),
dy1_(
owner.emplaceObject<realx3PointField_D>(
objectFile(
integration(baseName, pStruct, method, initialValField),
realx3PointField_D
(
objectFile
(
groupNames(baseName,"dy1"),
"",
pStruct.time().integrationFolder(),
objectFile::READ_IF_PRESENT,
objectFile::WRITE_ALWAYS),
objectFile::WRITE_ALWAYS
),
pStruct,
zero3))
{
}
zero3,
zero3
)
{}
bool pFlow::AdamsBashforth2::predict
(
real UNUSED(dt),
realx3Vector_D& UNUSED(y),
realx3Vector_D& UNUSED(dy)
realx3PointField_D& UNUSED(y),
realx3PointField_D& UNUSED(dy)
)
{
return true;
}
bool pFlow::AdamsBashforth2::predict
(
real dt,
realx3Field_D &y,
realx3PointField_D &dy
)
{
return true;
}
bool pFlow::AdamsBashforth2::correct
(
real dt,
realx3Vector_D& y,
realx3Vector_D& dy
realx3PointField_D& y,
realx3PointField_D& dy
)
{
if(this->pStruct().allActive())
return correct(dt, y.field(), dy);
}
bool pFlow::AdamsBashforth2::correct(real dt, realx3Field_D &y, realx3PointField_D &dy)
{
return intAll(dt, y, dy, this->pStruct().activeRange());
auto& dy1l = dy1();
if(dy1l.isAllActive())
{
return intAllActive(dt, y, dy, dy1l);
}
else
{
return intRange(dt, y, dy, this->pStruct().activePointsMaskD());
return intScattered(dt, y, dy, dy1l);
}
return true;
return false;
}
bool pFlow::AdamsBashforth2::setInitialVals(
@ -81,25 +168,3 @@ bool pFlow::AdamsBashforth2::setInitialVals(
return true;
}
bool pFlow::AdamsBashforth2::intAll(
real dt,
realx3Vector_D& y,
realx3Vector_D& dy,
range activeRng)
{
auto d_dy = dy.deviceVectorAll();
auto d_y = y.deviceVectorAll();
auto d_dy1= dy1_.deviceVectorAll();
Kokkos::parallel_for(
"AdamsBashforth2::correct",
rpIntegration (activeRng.first, activeRng.second),
LAMBDA_HD(int32 i){
d_y[i] += dt*(static_cast<real>(3.0 / 2.0) * d_dy[i] - static_cast<real>(1.0 / 2.0) * d_dy1[i]);
d_dy1[i] = d_dy[i];
});
Kokkos::fence();
return true;
}

114
src/Integration/AdamsBashforth2/AdamsBashforth2.hpp
View File

@ -36,41 +36,32 @@ namespace pFlow
*/
class AdamsBashforth2
:
public integration
public integration,
public realx3PointField_D
{
protected:
private:
/// dy at t-dt
realx3PointField_D& dy1_;
/// Range policy for integration kernel (alias)
using rpIntegration = Kokkos::RangePolicy<
DefaultExecutionSpace,
Kokkos::Schedule<Kokkos::Static>,
Kokkos::IndexType<int32>
>;
auto& dy1()
{
return static_cast<realx3PointField_D&>(*this);
}
public:
/// Type info
TypeInfo("AdamsBashforth2");
/// Class info
ClassInfo("AdamsBashforth2");
// - Constructors
/// Construct from components
AdamsBashforth2(
const word& baseName,
repository& owner,
const pointStructure& pStruct,
const word& method);
uniquePtr<integration> clone()const override
{
return makeUnique<AdamsBashforth2>(*this);
}
pointStructure& pStruct,
const word& method,
const realx3Field_D& initialValField);
/// Destructor
virtual ~AdamsBashforth2()=default;
~AdamsBashforth2()final = default;
/// Add this to the virtual constructor table
add_vCtor(
@ -80,71 +71,52 @@ public:
// - Methods
/// return integration method
word method()const override
{
return "AdamsBashforth2";
}
bool predict(
real UNUSED(dt),
realx3Vector_D& UNUSED(y),
realx3Vector_D& UNUSED(dy)) override;
realx3PointField_D& UNUSED(y),
realx3PointField_D& UNUSED(dy)) final;
bool predict(
real dt,
realx3Field_D& y,
realx3PointField_D& dy) final;
bool correct(
real dt,
realx3Vector_D& y,
realx3Vector_D& dy) override;
realx3PointField_D& y,
realx3PointField_D& dy) final;
bool correct(
real dt,
realx3Field_D& y,
realx3PointField_D& dy) final;
/*bool hearChanges
(
real t,
real dt,
uint32 iter,
const message& msg,
const anyList& varList
) override;*/
bool setInitialVals(
const int32IndexContainer& newIndices,
const realx3Vector& y) override;
const realx3Vector& y) final;
bool needSetInitialVals()const override
bool needSetInitialVals()const final
{
return false;
}
/// Integrate on all points in the active range
bool intAll(
real dt,
realx3Vector_D& y,
realx3Vector_D& dy,
range activeRng);
/// Integrate on active points in the active range
template<typename activeFunctor>
bool intRange(
real dt,
realx3Vector_D& y,
realx3Vector_D& dy,
activeFunctor activeP );
};
template<typename activeFunctor>
bool pFlow::AdamsBashforth2::intRange(
real dt,
realx3Vector_D& y,
realx3Vector_D& dy,
activeFunctor activeP )
{
auto d_dy = dy.deviceVectorAll();
auto d_y = y.deviceVectorAll();
auto d_dy1= dy1_.deviceVectorAll();
auto activeRng = activeP.activeRange();
Kokkos::parallel_for(
"AdamsBashforth2::correct",
rpIntegration (activeRng.first, activeRng.second),
LAMBDA_HD(int32 i){
if( activeP(i))
{
d_y[i] += dt*(static_cast<real>(3.0 / 2.0) * d_dy[i] - static_cast<real>(1.0 / 2.0) * d_dy1[i]);
d_dy1[i] = d_dy[i];
}
});
Kokkos::fence();
return true;
}
} // pFlow

6
src/Integration/AdamsBashforth3/AdamsBashforth3.cpp
View File

@ -89,9 +89,9 @@ bool pFlow::AdamsBashforth3::intAll(
realx3Vector_D& dy,
range activeRng)
{
auto d_dy = dy.deviceVectorAll();
auto d_y = y.deviceVectorAll();
auto d_history = history_.deviceVectorAll();
auto d_dy = dy.deviceViewAll();
auto d_y = y.deviceViewAll();
auto d_history = history_.deviceViewAll();
Kokkos::parallel_for(
"AdamsBashforth3::correct",

6
src/Integration/AdamsBashforth3/AdamsBashforth3.hpp
View File

@ -160,9 +160,9 @@ bool pFlow::AdamsBashforth3::intRange(
realx3Vector_D& dy,
activeFunctor activeP )
{
auto d_dy = dy.deviceVectorAll();
auto d_y = y.deviceVectorAll();
auto d_history = history_.deviceVectorAll();
auto d_dy = dy.deviceViewAll();
auto d_y = y.deviceViewAll();
auto d_history = history_.deviceViewAll();
auto activeRng = activeP.activeRange();
Kokkos::parallel_for(

6
src/Integration/AdamsBashforth4/AdamsBashforth4.cpp
View File

@ -89,9 +89,9 @@ bool pFlow::AdamsBashforth4::intAll(
realx3Vector_D& dy,
range activeRng)
{
auto d_dy = dy.deviceVectorAll();
auto d_y = y.deviceVectorAll();
auto d_history = history_.deviceVectorAll();
auto d_dy = dy.deviceViewAll();
auto d_y = y.deviceViewAll();
auto d_history = history_.deviceViewAll();
Kokkos::parallel_for(
"AdamsBashforth4::correct",

6
src/Integration/AdamsBashforth4/AdamsBashforth4.hpp
View File

@ -165,9 +165,9 @@ bool pFlow::AdamsBashforth4::intRange(
realx3Vector_D& dy,
activeFunctor activeP )
{
auto d_dy = dy.deviceVectorAll();
auto d_y = y.deviceVectorAll();
auto d_history = history_.deviceVectorAll();
auto d_dy = dy.deviceViewAll();
auto d_y = y.deviceViewAll();
auto d_history = history_.deviceViewAll();
auto activeRng = activeP.activeRange();
Kokkos::parallel_for(

6
src/Integration/AdamsBashforth5/AdamsBashforth5.cpp
View File

@ -89,9 +89,9 @@ bool pFlow::AdamsBashforth5::intAll(
realx3Vector_D& dy,
range activeRng)
{
auto d_dy = dy.deviceVectorAll();
auto d_y = y.deviceVectorAll();
auto d_history = history_.deviceVectorAll();
auto d_dy = dy.deviceViewAll();
auto d_y = y.deviceViewAll();
auto d_history = history_.deviceViewAll();
Kokkos::parallel_for(
"AdamsBashforth5::correct",

6
src/Integration/AdamsBashforth5/AdamsBashforth5.hpp
View File

@ -165,9 +165,9 @@ bool pFlow::AdamsBashforth5::intRange(
realx3Vector_D& dy,
activeFunctor activeP )
{
auto d_dy = dy.deviceVectorAll();
auto d_y = y.deviceVectorAll();
auto d_history = history_.deviceVectorAll();
auto d_dy = dy.deviceViewAll();
auto d_y = y.deviceViewAll();
auto d_history = history_.deviceViewAll();
auto activeRng = activeP.activeRange();
Kokkos::parallel_for(

20
src/Integration/AdamsMoulton3/AdamsMoulton3.cpp
View File

@ -124,11 +124,11 @@ bool pFlow::AdamsMoulton3::predictAll(
range activeRng)
{
auto d_dy = dy.deviceVectorAll();
auto d_y = y.deviceVectorAll();
auto d_y0 = y0_.deviceVectorAll();
auto d_dy0 = dy0_.deviceVectorAll();
auto d_dy1= dy1_.deviceVectorAll();
auto d_dy = dy.deviceViewAll();
auto d_y = y.deviceViewAll();
auto d_y0 = y0_.deviceViewAll();
auto d_dy0 = dy0_.deviceViewAll();
auto d_dy1= dy1_.deviceViewAll();
Kokkos::parallel_for(
"AdamsMoulton3::predict",
@ -150,12 +150,12 @@ bool pFlow::AdamsMoulton3::intAll(
range activeRng)
{
auto d_dy = dy.deviceVectorAll();
auto d_y = y.deviceVectorAll();
auto d_dy = dy.deviceViewAll();
auto d_y = y.deviceViewAll();
auto d_dy0 = dy0_.deviceVectorAll();
auto d_y0 = y0_.deviceVectorAll();
auto d_dy1 = dy1_.deviceVectorAll();
auto d_dy0 = dy0_.deviceViewAll();
auto d_y0 = y0_.deviceViewAll();
auto d_dy1 = dy1_.deviceViewAll();
Kokkos::parallel_for(
"AdamsMoulton3::correct",

20
src/Integration/AdamsMoulton3/AdamsMoulton3.hpp
View File

@ -144,11 +144,11 @@ bool AdamsMoulton3::predictRange(
realx3Vector_D& dy,
activeFunctor activeP)
{
auto d_dy = dy.deviceVectorAll();
auto d_y = y.deviceVectorAll();
auto d_y0 = y0_.deviceVectorAll();
auto d_dy0 = dy0_.deviceVectorAll();
auto d_dy1= dy1_.deviceVectorAll();
auto d_dy = dy.deviceViewAll();
auto d_y = y.deviceViewAll();
auto d_y0 = y0_.deviceViewAll();
auto d_dy0 = dy0_.deviceViewAll();
auto d_dy1= dy1_.deviceViewAll();
auto activeRng = activeP.activeRange();
@ -182,12 +182,12 @@ bool pFlow::AdamsMoulton3::intRange(
activeFunctor activeP)
{
auto d_dy = dy.deviceVectorAll();
auto d_y = y.deviceVectorAll();
auto d_dy = dy.deviceViewAll();
auto d_y = y.deviceViewAll();
auto d_dy0 = dy0_.deviceVectorAll();
auto d_y0 = y0_.deviceVectorAll();
auto d_dy1 = dy1_.deviceVectorAll();
auto d_dy0 = dy0_.deviceViewAll();
auto d_y0 = y0_.deviceViewAll();
auto d_dy1 = dy1_.deviceViewAll();
auto activeRng = activeP.activeRange();

24
src/Integration/AdamsMoulton4/AdamsMoulton4.cpp
View File

@ -135,13 +135,13 @@ bool pFlow::AdamsMoulton4::predictAll(
range activeRng)
{
auto d_dy = dy.deviceVectorAll();
auto d_y = y.deviceVectorAll();
auto d_dy = dy.deviceViewAll();
auto d_y = y.deviceViewAll();
auto d_y0 = y0_.deviceVectorAll();
auto d_dy0 = dy0_.deviceVectorAll();
auto d_dy1 = dy1_.deviceVectorAll();
auto d_dy2 = dy2_.deviceVectorAll();
auto d_y0 = y0_.deviceViewAll();
auto d_dy0 = dy0_.deviceViewAll();
auto d_dy1 = dy1_.deviceViewAll();
auto d_dy2 = dy2_.deviceViewAll();
Kokkos::parallel_for(
"AdamsMoulton4::predict",
@ -165,13 +165,13 @@ bool pFlow::AdamsMoulton4::intAll(
range activeRng)
{
auto d_dy = dy.deviceVectorAll();
auto d_y = y.deviceVectorAll();
auto d_dy = dy.deviceViewAll();
auto d_y = y.deviceViewAll();
auto d_dy0 = dy0_.deviceVectorAll();
auto d_y0 = y0_.deviceVectorAll();
auto d_dy1 = dy1_.deviceVectorAll();
auto d_dy2 = dy2_.deviceVectorAll();
auto d_dy0 = dy0_.deviceViewAll();
auto d_y0 = y0_.deviceViewAll();
auto d_dy1 = dy1_.deviceViewAll();
auto d_dy2 = dy2_.deviceViewAll();
Kokkos::parallel_for(
"AdamsMoulton4::correct",

24
src/Integration/AdamsMoulton4/AdamsMoulton4.hpp
View File

@ -147,13 +147,13 @@ bool AdamsMoulton4::predictRange(
realx3Vector_D& dy,
activeFunctor activeP )
{
auto d_dy = dy.deviceVectorAll();
auto d_y = y.deviceVectorAll();
auto d_dy = dy.deviceViewAll();
auto d_y = y.deviceViewAll();
auto d_y0 = y0_.deviceVectorAll();
auto d_dy0 = dy0_.deviceVectorAll();
auto d_dy1 = dy1_.deviceVectorAll();
auto d_dy2 = dy2_.deviceVectorAll();
auto d_y0 = y0_.deviceViewAll();
auto d_dy0 = dy0_.deviceViewAll();
auto d_dy1 = dy1_.deviceViewAll();
auto d_dy2 = dy2_.deviceViewAll();
auto activeRng = activeP.activeRange();
@ -185,13 +185,13 @@ bool pFlow::AdamsMoulton4::intRange(
activeFunctor activeP )
{
auto d_dy = dy.deviceVectorAll();
auto d_y = y.deviceVectorAll();
auto d_dy = dy.deviceViewAll();
auto d_y = y.deviceViewAll();
auto d_dy0 = dy0_.deviceVectorAll();
auto d_y0 = y0_.deviceVectorAll();
auto d_dy1 = dy1_.deviceVectorAll();
auto d_dy2 = dy2_.deviceVectorAll();
auto d_dy0 = dy0_.deviceViewAll();
auto d_y0 = y0_.deviceViewAll();
auto d_dy1 = dy1_.deviceViewAll();
auto d_dy2 = dy2_.deviceViewAll();
auto activeRng = activeP.activeRange();

28
src/Integration/AdamsMoulton5/AdamsMoulton5.cpp
View File

@ -145,14 +145,14 @@ bool pFlow::AdamsMoulton5::predictAll(
range activeRng)
{
auto d_dy = dy.deviceVectorAll();
auto d_y = y.deviceVectorAll();
auto d_dy = dy.deviceViewAll();
auto d_y = y.deviceViewAll();
auto d_y0 = y0_.deviceVectorAll();
auto d_dy0 = dy0_.deviceVectorAll();
auto d_dy1 = dy1_.deviceVectorAll();
auto d_dy2 = dy2_.deviceVectorAll();
auto d_dy3 = dy3_.deviceVectorAll();
auto d_y0 = y0_.deviceViewAll();
auto d_dy0 = dy0_.deviceViewAll();
auto d_dy1 = dy1_.deviceViewAll();
auto d_dy2 = dy2_.deviceViewAll();
auto d_dy3 = dy3_.deviceViewAll();
Kokkos::parallel_for(
"AdamsMoulton5::predict",
@ -178,14 +178,14 @@ bool pFlow::AdamsMoulton5::intAll(
range activeRng)
{
auto d_dy = dy.deviceVectorAll();
auto d_y = y.deviceVectorAll();
auto d_dy = dy.deviceViewAll();
auto d_y = y.deviceViewAll();
auto d_dy0 = dy0_.deviceVectorAll();
auto d_y0 = y0_.deviceVectorAll();
auto d_dy1 = dy1_.deviceVectorAll();
auto d_dy2 = dy2_.deviceVectorAll();
auto d_dy3 = dy3_.deviceVectorAll();
auto d_dy0 = dy0_.deviceViewAll();
auto d_y0 = y0_.deviceViewAll();
auto d_dy1 = dy1_.deviceViewAll();
auto d_dy2 = dy2_.deviceViewAll();
auto d_dy3 = dy3_.deviceViewAll();
Kokkos::parallel_for(
"AdamsMoulton5::correct",

28
src/Integration/AdamsMoulton5/AdamsMoulton5.hpp
View File

@ -150,14 +150,14 @@ bool AdamsMoulton5::predictRange(
realx3Vector_D& dy,
activeFunctor activeP )
{
auto d_dy = dy.deviceVectorAll();
auto d_y = y.deviceVectorAll();
auto d_dy = dy.deviceViewAll();
auto d_y = y.deviceViewAll();
auto d_y0 = y0_.deviceVectorAll();
auto d_dy0 = dy0_.deviceVectorAll();
auto d_dy1 = dy1_.deviceVectorAll();
auto d_dy2 = dy2_.deviceVectorAll();
auto d_dy3 = dy3_.deviceVectorAll();
auto d_y0 = y0_.deviceViewAll();
auto d_dy0 = dy0_.deviceViewAll();
auto d_dy1 = dy1_.deviceViewAll();
auto d_dy2 = dy2_.deviceViewAll();
auto d_dy3 = dy3_.deviceViewAll();
auto activeRng = activeP.activeRange();
@ -189,14 +189,14 @@ bool pFlow::AdamsMoulton5::intRange(
activeFunctor activeP )
{
auto d_dy = dy.deviceVectorAll();
auto d_y = y.deviceVectorAll();
auto d_dy = dy.deviceViewAll();
auto d_y = y.deviceViewAll();
auto d_dy0 = dy0_.deviceVectorAll();
auto d_y0 = y0_.deviceVectorAll();
auto d_dy1 = dy1_.deviceVectorAll();
auto d_dy2 = dy2_.deviceVectorAll();
auto d_dy3 = dy3_.deviceVectorAll();
auto d_dy0 = dy0_.deviceViewAll();
auto d_y0 = y0_.deviceViewAll();
auto d_dy1 = dy1_.deviceViewAll();
auto d_dy2 = dy2_.deviceViewAll();
auto d_dy3 = dy3_.deviceViewAll();
auto activeRng = activeP.activeRange();

12
src/Integration/CMakeLists.txt
View File

@ -1,13 +1,13 @@
list(APPEND SourceFiles
integration/integration.cpp
AdamsBashforth5/AdamsBashforth5.cpp
AdamsBashforth4/AdamsBashforth4.cpp
AdamsBashforth3/AdamsBashforth3.cpp
AdamsBashforth2/AdamsBashforth2.cpp
AdamsMoulton3/AdamsMoulton3.cpp
AdamsMoulton4/AdamsMoulton4.cpp
AdamsMoulton5/AdamsMoulton5.cpp
#AdamsBashforth5/AdamsBashforth5.cpp
#AdamsBashforth4/AdamsBashforth4.cpp
#AdamsBashforth3/AdamsBashforth3.cpp
#AdamsMoulton3/AdamsMoulton3.cpp
#AdamsMoulton4/AdamsMoulton4.cpp
#AdamsMoulton5/AdamsMoulton5.cpp
)
set(link_libs Kokkos::kokkos phasicFlow)

30
src/Integration/integration/integration.cpp
View File

@ -19,33 +19,35 @@ Licence:
-----------------------------------------------------------------------------*/
#include "integration.hpp"
#include "pointStructure.hpp"
#include "repository.hpp"
pFlow::integration::integration
(
const word& baseName,
repository& owner,
const pointStructure& pStruct,
const word& method
pointStructure& pStruct,
const word&,
const realx3Field_D&
)
:
owner_(owner),
baseName_(baseName),
pStruct_(pStruct)
{
CONSUME_PARAM(method);
}
owner_(*pStruct.owner()),
pStruct_(pStruct),
baseName_(baseName)
{}
pFlow::uniquePtr<pFlow::integration>
pFlow::integration::create(
pFlow::integration::create
(
const word& baseName,
repository& owner,
const pointStructure& pStruct,
const word& method)
pointStructure& pStruct,
const word& method,
const realx3Field_D& initialValField
)
{
if( wordvCtorSelector_.search(method) )
{
return wordvCtorSelector_[method] (baseName, owner, pStruct, method);
return wordvCtorSelector_[method] (baseName, pStruct, method, initialValField);
}
else
{

62
src/Integration/integration/integration.hpp
View File

@ -23,14 +23,16 @@ Licence:
#include "virtualConstructor.hpp"
#include "Vectors.hpp"
#include "pointStructure.hpp"
#include "repository.hpp"
#include "pointFields.hpp"
namespace pFlow
{
// - Forward
class repository;
class pointStructure;
/**
* Base class for integrating the first order ODE (IVP)
*
@ -48,19 +50,19 @@ namespace pFlow
*/
class integration
{
protected:
private:
// - Protected data members
/// The owner repository that all fields are storred in
repository& owner_;
/// The base name for integration
const word baseName_;
/// A reference to pointStructure
const pointStructure& pStruct_;
/// The base name for integration
const word baseName_;
public:
/// Type info
@ -72,9 +74,9 @@ public:
/// Construct from components
integration(
const word& baseName,
repository& owner,
const pointStructure& pStruct,
const word& method);
pointStructure& pStruct,
const word& method,
const realx3Field_D& initialValField);
/// Copy constructor
integration(const integration&) = default;
@ -88,22 +90,22 @@ public:
/// Move assignment
integration& operator = (integration&&) = default;
/// Polymorphic copy/cloning
virtual
uniquePtr<integration> clone()const=0;
/// Destructor
virtual ~integration()=default;
/// Add a virtual constructor
create_vCtor(
create_vCtor
(
integration,
word,
(const word& baseName,
repository& owner,
const pointStructure& pStruct,
const word& method),
(baseName, owner, pStruct, method) );
(
const word& baseName,
pointStructure& pStruct,
const word& method,
const realx3Field_D& initialValField
),
(baseName, pStruct, method, initialValField)
);
// - Methods
@ -129,13 +131,23 @@ public:
return owner_;
}
/// return integration method
virtual
word method()const = 0 ;
/// Prediction step in integration
virtual
bool predict(real dt, realx3Vector_D& y, realx3Vector_D& dy) = 0;
bool predict(real dt, realx3PointField_D& y, realx3PointField_D& dy) = 0;
virtual
bool predict(real dt, realx3Field_D& y, realx3PointField_D& dy) = 0;
/// Correction/main integration step
virtual
bool correct(real dt, realx3Vector_D& y, realx3Vector_D& dy) = 0;
bool correct(real dt, realx3PointField_D& y, realx3PointField_D& dy) = 0;
virtual
bool correct(real dt, realx3Field_D& y, realx3PointField_D& dy) = 0;
/// Set the initial values for new indices
virtual
@ -152,9 +164,9 @@ public:
static
uniquePtr<integration> create(
const word& baseName,
repository& owner,
const pointStructure& pStruct,
const word& method);
pointStructure& pStruct,
const word& method,
const realx3Field_D& initialValField);
}; // integration

18
src/Interaction/CMakeLists.txt
View File

@ -1,9 +1,25 @@
set(SourceFiles
contactSearch/methods/cellBased/NBS/mapperNBS.cpp
contactSearch/methods/cellBased/NBS/mapperNBSKernels.cpp
contactSearch/methods/cellBased/NBS/NBSLevel0.cpp
contactSearch/methods/cellBased/NBS/NBS.cpp
contactSearch/methods/cellBased/NBS/cellsWallLevel0.cpp
contactSearch/boundaries/boundaryContactSearch/boundaryContactSearch.cpp
contactSearch/boundaries/periodicBoundaryContactSearch/ppwBndryContactSearchKernels.cpp
contactSearch/boundaries/periodicBoundaryContactSearch/ppwBndryContactSearch.cpp
contactSearch/boundaries/periodicBoundaryContactSearch/wallBoundaryContactSearch.cpp
contactSearch/boundaries/periodicBoundaryContactSearch/periodicBoundaryContactSearch.cpp
contactSearch/boundaries/boundaryContactSearchList.cpp
contactSearch/contactSearch/contactSearch.cpp
contactSearch/ContactSearch/ContactSearchs.cpp
interaction/interaction.cpp
sphereInteraction/sphereInteractions.cpp
sphereInteraction/sphereInteractionsLinearModels.cpp
sphereInteraction/sphereInteractionsNonLinearModels.cpp
sphereInteraction/sphereInteractionsNonLinearModModels.cpp
)
set(link_libs Kokkos::kokkos phasicFlow Property Particles Geometry)

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

@ -133,8 +133,9 @@ protected:
return false;
}
realVector etha_n(nElem);
realVector etha_t(nElem);
realVector etha_n("etha_n", nElem);
realVector etha_t("etha_t", nElem);
ForAll(i , kn)
{
etha_n[i] = -2.0*log(en[i])*sqrt(kn[i])/
@ -144,7 +145,7 @@ protected:
sqrt(pow(log(et[i]),2.0)+ pow(Pi,2.0));
}
Vector<linearProperties> prop(nElem);
Vector<linearProperties> prop("prop", nElem);
ForAll(i,kn)
{
prop[i] = {kn[i], kt[i], etha_n[i], etha_t[i], mu[i]};
@ -219,10 +220,10 @@ public:
void contactForce
(
const real dt,
const int32 i,
const int32 j,
const int32 propId_i,
const int32 propId_j,
const uint32 i,
const uint32 j,
const uint32 propId_i,
const uint32 propId_j,
const real Ri,
const real Rj,
const real ovrlp_n,

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

@ -121,7 +121,7 @@ protected:
}
realVector etha_n(nElem);
realVector etha_n("etha_n",nElem);
ForAll(i , en)
{
@ -137,7 +137,7 @@ protected:
}
Vector<nonLinearProperties> prop(nElem);
Vector<nonLinearProperties> prop("prop",nElem);
ForAll(i,Yeff)
{
prop[i] = {Yeff[i], Geff[i], etha_n[i], mu[i]};
@ -214,10 +214,10 @@ public:
void contactForce
(
const real dt,
const int32 i,
const int32 j,
const int32 propId_i,
const int32 propId_j,
const uint32 i,
const uint32 j,
const uint32 propId_i,
const uint32 propId_j,
const real Ri,
const real Rj,
const real ovrlp_n,

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

@ -121,7 +121,7 @@ protected:
}
Vector<nonLinearProperties> prop(nElem);
Vector<nonLinearProperties> prop("prop",nElem);
ForAll(i,Yeff)
{
prop[i] = {Yeff[i], Geff[i], etha_n[i], mu[i]};
@ -198,10 +198,10 @@ public:
void contactForce
(
const real dt,
const int32 i,
const int32 j,
const int32 propId_i,
const int32 propId_j,
const uint32 i,
const uint32 j,
const uint32 propId_i,
const uint32 propId_j,
const real Ri,
const real Rj,
const real ovrlp_n,

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

@ -76,10 +76,10 @@ public:
void rollingFriction
(
const real dt,
const int32 i,
const int32 j,
const int32 propId_i,
const int32 propId_j,
const uint32 i,
const uint32 j,
const uint32 propId_i,
const uint32 propId_j,
const real Ri,
const real Rj,
const realx3& wi,

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

@ -53,7 +53,7 @@ protected:
ViewType1D<ValueType,ExecutionSpace> values_;
int32 size0_ = 0;
uint32 size0_ = 0;
ViewType1D<PairType,ExecutionSpace> sortedPairs0_;
@ -73,7 +73,7 @@ protected:
using rpReFillPairs = Kokkos::RangePolicy<
ExecutionSpace,
Kokkos::Schedule<Kokkos::Static>,
Kokkos::IndexType<int32>,
Kokkos::IndexType<uint32>,
TagReFillPairs>;
public:
@ -81,7 +81,7 @@ public:
TypeInfoNV("sortedContactList");
sortedContactList(int32 initialSize =1)
explicit sortedContactList(uint32 initialSize =1)
:
SortedPairs(initialSize),
values_("values", SortedPairs::capacity()),
@ -114,31 +114,31 @@ public:
}
INLINE_FUNCTION_HD
ValueType getValue(int32 idx)const
ValueType getValue(uint32 idx)const
{
return values_[idx];
}
INLINE_FUNCTION_HD
void setValue(int32 idx, const ValueType& val)const
void setValue(uint32 idx, const ValueType& val)const
{
values_[idx] = val;
}
INLINE_FUNCTION_HD
void operator()(TagReFillPairs, int32 idx)const
void operator()(TagReFillPairs, uint32 idx)const
{
auto searchLen = max(size0_/1000,10);
auto start = max(0,idx-searchLen);
auto end = min(size0_,idx+searchLen);
uint32 searchLen = max(size0_/1000u,10u);
uint32 start = idx-min(searchLen,idx);
uint32 end = min(size0_,idx+searchLen);
auto newPair = this->sortedPairs_[idx];
if( auto idx0 = binarySearch(
sortedPairs0_,
start,
end,
newPair);
idx0>=0)
idx0!=-1)
{
values_[idx] = values0_[idx0];
}
@ -147,7 +147,7 @@ public:
start,
end,
newPair);
idx0>=0)
idx0!=-1)
{
values_[idx] = values0_[idx0];

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

@ -52,24 +52,24 @@ public:
{
using PairType = typename sortedPairs::PairType;
int32 size_;
uint32 size_;
ViewType1D<PairType,ExecutionSpace> sortedParis_;
INLINE_FUNCTION_HD
int32 size()const { return size_; }
uint32 size()const { return size_; }
INLINE_FUNCTION_HD
int32 loopCount()const { return size_; }
uint32 loopCount()const { return size_; }
INLINE_FUNCTION_HD
bool isValid(int32 i)const { return i<size_; }
bool isValid(uint32 i)const { return i<size_; }
INLINE_FUNCTION_HD
PairType getPair(int i)const { return sortedParis_[i]; }
PairType getPair(uint32 i)const { return sortedParis_[i]; }
INLINE_FUNCTION_HD
bool getPair(int32 i, PairType& pair)const {
bool getPair(uint32 i, PairType& pair)const {
if(i<size_) {
pair = sortedParis_[i];
return true;
@ -85,22 +85,22 @@ public:
protected:
/// size of pair list
int32 size_ = 0;
uint32 size_ = 0;
ViewType1D<int32,ExecutionSpace> flags_;
ViewType1D<uint32,ExecutionSpace> flags_;
ViewType1D<PairType,ExecutionSpace> sortedPairs_;
using rpFillFlag = Kokkos::RangePolicy<
ExecutionSpace,
Kokkos::Schedule<Kokkos::Static>,
Kokkos::IndexType<int32>,
Kokkos::IndexType<uint32>,
TagFillFlag >;
using rpFillPairs = Kokkos::RangePolicy<
ExecutionSpace,
Kokkos::Schedule<Kokkos::Static>,
Kokkos::IndexType<int32>,
Kokkos::IndexType<uint32>,
TagFillPairs>;
public:
@ -110,7 +110,7 @@ public:
// constructors
sortedPairs(int32 initialSize =1)
explicit sortedPairs(uint32 initialSize =1)
:
UnsortedPairs(initialSize),
flags_("flags_",UnsortedPairs::capacity()+1),
@ -134,7 +134,7 @@ public:
// return the pair at index idx
// perform no check for size and existance
INLINE_FUNCTION_HD
PairType getPair(int32 idx)const
PairType getPair(uint32 idx)const
{
return sortedPairs_[idx];
}
@ -142,7 +142,7 @@ public:
// - Device/host call
// return the pair at index idx
INLINE_FUNCTION_HD
bool getPair(int32 idx, PairType& p)const
bool getPair(uint32 idx, PairType& p)const
{
if(isValid(idx))
{
@ -156,7 +156,7 @@ public:
}
INLINE_FUNCTION_HD
bool isValid(int32 idx)const
bool isValid(uint32 idx)const
{
return idx < size_;
}
@ -164,12 +164,12 @@ public:
//use this when the value of size_ is updated
INLINE_FUNCTION_H
int32 size()const
uint32 size()const
{
return size_;
}
int32 loopCount()const
uint32 loopCount()const
{
return size_;
}
@ -189,7 +189,7 @@ public:
void prepareSorted()
{
// first check the size of flags_
int32 capacity = UnsortedPairs::capacity();
uint32 capacity = UnsortedPairs::capacity();
if( capacity+1 > flags_.size() )
{
@ -218,7 +218,7 @@ public:
if( size_>sortedPairs_.size() )
{
// get more space to prevent reallocations in next iterations
int32 len = size_*1.1+1;
uint32 len = size_*1.1+1;
reallocNoInit(sortedPairs_, len);
}
@ -235,7 +235,7 @@ public:
}
INLINE_FUNCTION_HD
void operator()(TagFillFlag, int32 i)const
void operator()(TagFillFlag, uint32 i)const
{
if(this->container_.valid_at(i) )
flags_[i] = 1;
@ -244,7 +244,7 @@ public:
}
INLINE_FUNCTION_HD
void operator()(TagFillPairs, int32 i)const
void operator()(TagFillPairs, uint32 i)const
{
auto fi = flags_[i];
if(fi!=flags_[i+1])

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

@ -17,10 +17,11 @@ Licence:
implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
-----------------------------------------------------------------------------*/
#ifndef __unsortedContactList_hpp__
#define __unsortedContactList_hpp__
#include "unsortedPairs.hpp"
namespace pFlow
{
@ -72,7 +73,7 @@ protected:
using rpFillPairs = Kokkos::RangePolicy<
ExecutionSpace,
Kokkos::Schedule<Kokkos::Static>,
Kokkos::IndexType<int32>,
Kokkos::IndexType<uint32>,
TagReFillPairs>;
@ -80,7 +81,7 @@ public:
TypeInfoNV("unsortedContactList");
unsortedContactList(int32 capacity=1)
explicit unsortedContactList(uint32 capacity=1)
:
UnsortedPairs(capacity),
values_("values", UnsortedPairs::capacity()),
@ -122,7 +123,7 @@ public:
INLINE_FUNCTION_HD
bool getValue(const PairType& p, ValueType& val)const
{
if(auto idx = this->find(p); idx>=0)
if(auto idx = this->find(p); idx!=-1)
{
val = getValue(idx);
return true;
@ -131,7 +132,7 @@ public:
}
INLINE_FUNCTION_HD
void setValue(int32 idx, const ValueType& val)const
void setValue(uint32 idx, const ValueType& val)const
{
values_[idx] = val;
}
@ -139,7 +140,7 @@ public:
INLINE_FUNCTION_HD
bool setValue(const PairType& p, const ValueType& val)const
{
if(auto idx = this->find(p); idx>=0)
if(uint32 idx = this->find(p); idx!=-1)
{
setValue(idx, val);
return true;;
@ -148,13 +149,13 @@ public:
}
INLINE_FUNCTION_HD
void operator()(TagReFillPairs, int32 idx)const
void operator()(TagReFillPairs, uint32 idx)const
{
if( this->isValid(idx) )
{
if( int32 idx0 =
if( uint32 idx0 =
container0_.find(this->getPair(idx));
idx0>=0 )
idx0!=-1 )
{
values_[idx] = values0_[idx0];
}

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

@ -41,7 +41,7 @@ public:
using memory_space = typename ExecutionSpace::memory_space;
using PairType = kPair<idType,idType>;
using PairType = Pair<idType,idType>;
using ContainerType = unorderedSet<PairType, ExecutionSpace>;
@ -52,19 +52,19 @@ public:
ContainerType Container_;
INLINE_FUNCTION_HD
int32 size()const { return Container_.size(); }
uint32 size()const { return Container_.size(); }
INLINE_FUNCTION_HD
int32 loopCount()const { return Container_.capacity(); }
uint32 loopCount()const { return Container_.capacity(); }
INLINE_FUNCTION_HD
bool isValid(int32 idx)const { return Container_.valid_at(idx); }
bool isValid(uint32 idx)const { return Container_.valid_at(idx); }
INLINE_FUNCTION_HD
PairType getPair(int idx)const { return Container_.key_at(idx); }
PairType getPair(uint32 idx)const { return Container_.key_at(idx); }
INLINE_FUNCTION_HD
bool getPair(int32 idx, PairType& pair)const {
bool getPair(uint32 idx, PairType& pair)const {
if(Container_.valid_at(idx)) {
pair = Container_.key_at(idx);
return true;
@ -84,7 +84,7 @@ public:
TypeInfoNV("unsorderedPairs");
// constructor
unsortedPairs(int32 capacity=1)
explicit unsortedPairs(uint32 capacity=1)
:
container_(capacity) // the minimum capacity would be 128
{}
@ -102,7 +102,7 @@ public:
// - Device call
INLINE_FUNCTION_HD
int32 insert(idType i, idType j)const
uint32 insert(idType i, idType j)const
{
if(auto insertResult = container_.insert(PairType(i,j)); insertResult.failed())
return -1;
@ -112,7 +112,7 @@ public:
}
INLINE_FUNCTION_HD
int32 insert(const PairType& p)const
uint32 insert(const PairType& p)const
{
if(auto insertResult = container_.insert(p); insertResult.failed())
return -1;
@ -125,7 +125,7 @@ public:
// return the pair at index idx
// perform no check for size and existance
INLINE_FUNCTION_HD
PairType getPair(int32 idx)const
PairType getPair(uint32 idx)const
{
return container_.key_at(idx);
}
@ -133,7 +133,7 @@ public:
// - Device call
// return the pair at index idx
INLINE_FUNCTION_HD
bool getPair(int32 idx, PairType& p)const
bool getPair(uint32 idx, PairType& p)const
{
if(container_.valid_at(idx))
{
@ -148,7 +148,7 @@ public:
}
INLINE_FUNCTION_HD
int32 find(const PairType & p)const
uint32 find(const PairType & p)const
{
if( auto idx = container_.find(p);
idx != Kokkos::UnorderedMapInvalidIndex )
@ -158,26 +158,26 @@ public:
}
INLINE_FUNCTION_HD
bool isValid(int32 idx)const
bool isValid(uint32 idx)const
{
return container_.valid_at(idx);
}
INLINE_FUNCTION_HD
int32 capacity() const
uint32 capacity() const
{
return container_.capacity();
}
int32 loopCount()const
uint32 loopCount()const
{
return container_.capacity();
}
//use this when the value of size_ is updated
INLINE_FUNCTION_H
int32 size()const
uint32 size()const
{
return container_.size();
}
@ -190,7 +190,7 @@ public:
/// increase the capacity of the container by at-least len
/// the content will be erased.
INLINE_FUNCTION_H
void increaseCapacityBy(int32 len)
void increaseCapacityBy(uint32 len)
{
uint newCap = container_.capacity()+len;
this->clear();

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

@ -22,16 +22,18 @@ Licence:
#ifndef __ContactSearch_hpp__
#define __ContactSearch_hpp__
#include "contactSearchGlobals.hpp"
#include "contactSearch.hpp"
#include "box.hpp"
#include "particles.hpp"
#include "geometry.hpp"
#include "boundaryContactSearchList.hpp"
namespace pFlow
{
template<
template<class> class BaseMethod,
template<class> class WallMapping
class searchMethod
>
class ContactSearch
:
@ -39,89 +41,77 @@ class ContactSearch
{
public:
using IdType = typename contactSearch::IdType;
using IdType = uint32;
using IndexType = typename contactSearch::IndexType;
using ExecutionSpace = DefaultExecutionSpace;
using ExecutionSpace = typename contactSearch::ExecutionSpace;
using SearchMethodType = searchMethod;
using PairContainerType = typename contactSearch::PairContainerType;
private:
using ParticleContactSearchType =
BaseMethod<
ExecutionSpace>;
uniquePtr<SearchMethodType> ppwContactSearch_ = nullptr;
using WallMappingType =
WallMapping<
ExecutionSpace>;
protected:
uniquePtr<ParticleContactSearchType> particleContactSearch_ = nullptr;
uniquePtr<WallMappingType> wallMapping_ = nullptr;
boundaryContactSearchList csBoundaries_;
public:
TypeInfoTemplate2("ContactSearch", ParticleContactSearchType, WallMappingType);
TypeInfoTemplate11("ContactSearch", SearchMethodType);
ContactSearch(
const dictionary& csDict,
const box& domain,
const box& extDomain,
const particles& prtcl,
const geometry& geom,
Timers& timers)
:
contactSearch(csDict, domain, prtcl, geom, timers)
contactSearch(
csDict,
extDomain,
prtcl,
geom,
timers),
csBoundaries_(
csDict,
Particles().pStruct().boundaries(),
*this)
{
auto method = dict().getVal<word>("method");
auto wmMethod = dict().getVal<word>("wallMapping");
/*auto method = dict().getVal<word>("method");
auto nbDict = dict().subDict(method+"Info");
auto nbDict = dict().subDict(method+"Info");*/
real minD, maxD;
real minD;
real maxD;
this->Particles().boundingSphereMinMax(minD, maxD);
const auto& position = this->Particles().pointPosition().deviceVectorAll();
const auto& diam = this->Particles().boundingSphere().deviceVectorAll();
const auto& position = this->Particles().pointPosition().deviceViewAll();
const auto& flags = this->Particles().dynPointStruct().activePointsMaskDevice();
const auto& diam = this->Particles().boundingSphere().deviceViewAll();
particleContactSearch_ =
makeUnique<ParticleContactSearchType>
uint32 wnPoints = this->Geometry().numPoints();
uint32 wnTri = this->Geometry().size();
const auto& wPoints = this->Geometry().points().deviceViewAll();
const auto& wVertices = this->Geometry().vertices().deviceViewAll();
const auto& wNormals = this->Geometry().normals().deviceViewAll();
ppwContactSearch_ =
makeUnique<SearchMethodType>
(
nbDict,
this->domain(),
dict(),
this->extendedDomainBox(),
minD,
maxD,
position,
diam
);
REPORT(2)<<"Contact search algorithm for particle-particle is "<<
greenText(particleContactSearch_().typeName())<<endREPORT;
auto wmDict = dict().subDict(wmMethod+"Info");
int32 wnPoints = this->Geometry().numPoints();
int32 wnTri = this->Geometry().size();
const auto& wPoints = this->Geometry().points().deviceVectorAll();
const auto& wVertices = this->Geometry().vertices().deviceVectorAll();
wallMapping_ =
makeUnique<WallMappingType>(
wmDict,
particleContactSearch_().numLevels(),
particleContactSearch_().getCellsLevels(),
flags,
diam,
wnPoints,
wnTri,
wPoints,
wVertices
wVertices,
wNormals
);
REPORT(2)<<"Wall mapping algorithm for particle-wall is "<<
greenText(wallMapping_().typeName())<< endREPORT;
REPORT(2)<<"Contact search algorithm for particle-particle is "<<
Green_Text(ppwContactSearch_().typeName())<<END_REPORT;
}
@ -131,114 +121,74 @@ public:
ContactSearch,
dictionary);
bool broadSearch(
PairContainerType& ppPairs,
PairContainerType& pwPairs,
uint32 iter,
real t,
real dt,
csPairContainerType& ppPairs,
csPairContainerType& pwPairs,
bool force = false) override
{
ppTimer().start();
if(particleContactSearch_)
const auto& position = Particles().pointPosition().deviceViewAll();
const auto& flags = Particles().dynPointStruct().activePointsMaskDevice();
const auto& diam = Particles().boundingSphere().deviceViewAll();
if( !ppwContactSearch_().broadSearch(
iter,
t,
dt,
ppPairs,
pwPairs,
position,
flags,
diam,
force) )
{
auto activeRange = this->Particles().activeRange();
sphereSphereTimer_.start();
if(this->Particles().allActive())
{
particleContactSearch_().broadSearch(ppPairs, activeRange, force);
}
else
{
particleContactSearch_().broadSearch(ppPairs, activeRange, this->Particles().activePointsMaskD(), force);
}
sphereSphereTimer_.end();
}
else
fatalErrorInFunction;
return false;
if(wallMapping_)
{
sphereWallTimer_.start();
wallMapping_().broadSearch(pwPairs, particleContactSearch_(), force);
sphereWallTimer_.end();
}
else
return false;
ppTimer().end();
return true;
}
bool boundaryBroadSearch(
uint32 i,
uint32 iter,
real t,
real dt,
csPairContainerType& ppPairs,
csPairContainerType& pwPairs,
bool force = false)override
{
return csBoundaries_[i].broadSearch(
iter,
t,
dt,
ppPairs,
pwPairs,
force);
}
bool ppEnterBroadSearch()const override
bool enterBroadSearch(uint32 iter, real t, real dt)const override
{
if(particleContactSearch_)
if(ppwContactSearch_)
{
return particleContactSearch_().enterBoadSearch();
return ppwContactSearch_().performSearch(iter);
}
return false;
}
bool pwEnterBroadSearch()const override
bool performedBroadSearch()const override
{
if(wallMapping_)
{
return wallMapping_().enterBoadSearch();
return ppwContactSearch_().performedSearch();
}
return false;
}
bool ppPerformedBroadSearch()const override
{
if(particleContactSearch_)
{
return particleContactSearch_().performedSearch();
}
return false;
}
bool pwPerformedBroadSearch()const override
{
if(wallMapping_)
{
return wallMapping_().performedSearch();
}
return false;
}
/*bool update(const eventMessage& msg)
{
if(msg.isSizeChanged() )
{
auto newSize = this->prtcl().size();
if(!particleContactSearch_().objectSizeChanged(newSize))
{
fatalErrorInFunction<<
"erro in changing the size for particleContactSearch_ \n";
return false;
}
}
if(msg.isCapacityChanged() )
{
auto newSize = this->prtcl().capacity();
if(!particleContactSearch_().objectSizeChanged(newSize))
{
fatalErrorInFunction<<
"erro in changing the capacity for particleContactSearch_ \n";
return false;
}
}
return true;
}*/
};

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

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

77
src/Interaction/contactSearch/boundaries/boundaryContactSearch/boundaryContactSearch.cpp Normal file
View File

@ -0,0 +1,77 @@
/*------------------------------- 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 "boundaryContactSearch.hpp"
#include "contactSearch.hpp"
pFlow::boundaryContactSearch::boundaryContactSearch(
const dictionary &dict,
const boundaryBase &boundary,
const contactSearch &cSearch)
: generalBoundary(
boundary,
cSearch.pStruct(),
"",
""),
contactSearch_(cSearch),
updateInterval_(dict.getVal<uint32>("updateInterval"))
{
}
pFlow::uniquePtr<pFlow::boundaryContactSearch>
pFlow::boundaryContactSearch::create(
const dictionary &dict,
const boundaryBase &boundary,
const contactSearch &cSearch)
{
word bType = angleBracketsNames2(
"boundaryContactSearch",
pFlowProcessors().localRunTypeName(),
boundary.type());
word altBType{"boundaryContactSearch<none>"};
if( boundaryBasevCtorSelector_.search(bType) )
{
REPORT(2)<<"Creating contact search boundary "<< Green_Text(bType)<<
" for "<<boundary.name()<<endl;
return boundaryBasevCtorSelector_[bType](dict, boundary, cSearch);
}
else if(boundaryBasevCtorSelector_.search(altBType))
{
REPORT(2)<<"Creating contact search boundary "<< Green_Text(altBType)<<
" for "<<boundary.name()<<endl;
return boundaryBasevCtorSelector_[altBType](dict, boundary, cSearch);
}
else
{
printKeys
(
fatalError << "Ctor Selector "<< bType<<
" and "<< altBType << " do not exist. \n"
<<"Avaiable ones are: \n"
,
boundaryBasevCtorSelector_
);
fatalExit;
}
return nullptr;
}

114
src/Interaction/contactSearch/boundaries/boundaryContactSearch/boundaryContactSearch.hpp Normal file
View File

@ -0,0 +1,114 @@
/*------------------------------- 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 __boundaryContactSearch_hpp__
#define __boundaryContactSearch_hpp__
#include "generalBoundary.hpp"
#include "contactSearchGlobals.hpp"
#include "virtualConstructor.hpp"
namespace pFlow
{
class contactSearch;
class boundaryContactSearch
: public generalBoundary
{
private:
const contactSearch &contactSearch_;
/// @brief update interval in terms of iteration numebr
uint32 updateInterval_;
/// @brief last iteration number which contact search has been performed
uint32 lastUpdated_ = 0;
/// @brief performed search?
bool performedSearch_ = false;
public:
// type info
TypeInfo("boundaryContactSearch<none>");
boundaryContactSearch(
const dictionary &dict,
const boundaryBase &boundary,
const contactSearch &cSearch);
create_vCtor(
boundaryContactSearch,
boundaryBase,
(
const dictionary &dict,
const boundaryBase &boundary,
const contactSearch &cSearch),
(dict, boundary, cSearch));
add_vCtor(
boundaryContactSearch,
boundaryContactSearch,
boundaryBase);
const contactSearch &cSearch() const
{
return contactSearch_;
}
void fill(const std::any &val) override
{
return;
}
bool hearChanges(
real t,
real dt,
uint32 iter,
const message &msg,
const anyList &varList) override
{
if (msg.equivalentTo(message::BNDR_RESET))
{
// do nothing
}
return true;
}
virtual bool broadSearch(
uint32 iter,
real t,
real dt,
csPairContainerType &ppPairs,
csPairContainerType &pwPairs,
bool force = false)
{
return true;
}
static uniquePtr<boundaryContactSearch> create(
const dictionary &dict,
const boundaryBase &boundary,
const contactSearch &cSearch);
};
}
#endif //__boundaryContactSearch_hpp__

29
src/Interaction/contactSearch/boundaries/boundaryContactSearchList.cpp Normal file
View File

@ -0,0 +1,29 @@
#include "boundaryContactSearchList.hpp"
#include "boundaryList.hpp"
void pFlow::boundaryContactSearchList::setList(
const dictionary &dict,
const contactSearch &cSearch)
{
for(auto i=0; i<boundaries_.size(); i++)
{
this->set
(
i,
boundaryContactSearch::create(dict, boundaries_[i], cSearch)
);
}
}
pFlow::boundaryContactSearchList::boundaryContactSearchList(
const dictionary &dict,
const boundaryList& bndrs,
const contactSearch &cSearch)
:
ListPtr(bndrs.size()),
boundaries_(bndrs)
{
setList(dict, cSearch);
}

35
src/Interaction/contactSearch/boundaries/boundaryContactSearchList.hpp Normal file
View File

@ -0,0 +1,35 @@
#include "ListPtr.hpp"
#include "boundaryContactSearch.hpp"
namespace pFlow
{
class boundaryList;
class contactSearch;
class boundaryContactSearchList
:
public ListPtr<boundaryContactSearch>
{
private:
const boundaryList& boundaries_;
void setList(
const dictionary& dict,
const contactSearch& cSearch);
public:
TypeInfoNV("boundaryContactSearchList");
boundaryContactSearchList(
const dictionary& dict,
const boundaryList& bndrs,
const contactSearch& cSearch);
~boundaryContactSearchList()=default;
};
}

131
src/Interaction/contactSearch/boundaries/periodicBoundaryContactSearch/periodicBoundaryContactSearch.cpp Normal file
View File

@ -0,0 +1,131 @@
/*------------------------------- 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 "periodicBoundaryContactSearch.hpp"
#include "contactSearch.hpp"
#include "particles.hpp"
#include "pointStructure.hpp"
#include "geometry.hpp"
void pFlow::periodicBoundaryContactSearch::setSearchBox()
{
auto db = pStruct().thisDomain().domainBox();
auto n1 = boundary().mirrorBoundary().boundaryPlane().normal();
auto l1 = boundary().mirrorBoundary().neighborLength();
auto n2 = boundary().boundaryPlane().normal();
auto l2 = boundary().neighborLength();
realx3 minP = db.minPoint() + (db.maxPoint()-db.minPoint())* n1+(n2*l2);
realx3 maxP = db.maxPoint() + (n1*l1);
searchBox_={minP, maxP};
}
pFlow::periodicBoundaryContactSearch::periodicBoundaryContactSearch(
const dictionary &dict,
const boundaryBase &boundary,
const contactSearch &cSearch)
:
boundaryContactSearch(dict, boundary, cSearch),
transferVec_(boundary.mirrorBoundary().displacementVectroToMirror()),
thisIndex_(thisBoundaryIndex()),
mirrorIndex_(mirrorBoundaryindex()),
diameter_(cSearch.Particles().boundingSphere())
{
if(thisIndex_%2==1)
{
masterSearch_ = true;
setSearchBox();
real minD;
real maxD;
cSearch.Particles().boundingSphereMinMax(minD, maxD);
ppContactSearch_ = makeUnique<ppwBndryContactSearch>(
searchBox_,
maxD);
const auto& geom = cSearch.Geometry();
pwContactSearch_ = makeUnique<wallBoundaryContactSearch>(
0.5,
geom.numPoints(),
geom.size(),
geom.points().deviceViewAll(),
geom.vertices().deviceViewAll(),
geom.normals().deviceViewAll());
}
else
{
masterSearch_ = false;
searchBox_={{0,0,0},{0,0,0}};
}
}
bool pFlow::periodicBoundaryContactSearch::broadSearch
(
uint32 iter,
real t,
real dt,
csPairContainerType &ppPairs,
csPairContainerType &pwPairs,
bool force
)
{
if(masterSearch_)
{
auto thisP = boundary().thisPoints();
auto thisDiams = diameter_.BoundaryField(thisIndex_).thisField();
auto mirrorP = mirrorBoundary().thisPoints();
auto mirrorDiams = diameter_.BoundaryField(mirrorIndex_).thisField();
ppContactSearch_().broadSearchPP(
ppPairs,
thisP,
thisDiams,
mirrorP,
mirrorDiams,
transferVec_);
/*pwContactSearch_().broadSearch(
pwPairs,
ppContactSearch_().searchCells(),
thisP,
thisDiams,
mirrorP,
mirrorDiams,
transferVec_,
ppContactSearch_().sizeRatio());*/
//output<<t<<" boundary pp size "<< ppPairs.size()<<endl;
//output<<t<<" boundary pw size "<< pwPairs.size()<<endl;
return true;
}else
{
return true;
}
}

81
src/Interaction/contactSearch/boundaries/periodicBoundaryContactSearch/periodicBoundaryContactSearch.hpp Normal file
View File

@ -0,0 +1,81 @@
/*------------------------------- 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 __periodicBoundaryContactSearch_hpp__
#define __periodicBoundaryContactSearch_hpp__
#include "boundaryContactSearch.hpp"
#include "box.hpp"
#include "ppwBndryContactSearch.hpp"
#include "pointFields.hpp"
#include "wallBoundaryContactSearch.hpp"
namespace pFlow
{
class periodicBoundaryContactSearch
: public boundaryContactSearch
{
private:
box searchBox_;
realx3 transferVec_;
uint32 thisIndex_;
uint32 mirrorIndex_;
uniquePtr<ppwBndryContactSearch> ppContactSearch_ = nullptr;
uniquePtr<wallBoundaryContactSearch> pwContactSearch_ = nullptr;
const realPointField_D &diameter_;
bool masterSearch_ = false;
void setSearchBox();
public:
TypeInfo("boundaryContactSearch<regular,periodic>")
periodicBoundaryContactSearch(
const dictionary &dict,
const boundaryBase &boundary,
const contactSearch &cSearch);
~periodicBoundaryContactSearch() override = default;
add_vCtor(
boundaryContactSearch,
periodicBoundaryContactSearch,
boundaryBase);
bool broadSearch(
uint32 iter,
real t,
real dt,
csPairContainerType &ppPairs,
csPairContainerType &pwPairs,
bool force = false) override;
};
}
#endif //__periodicBoundaryContactSearch_hpp__

109
src/Interaction/contactSearch/boundaries/periodicBoundaryContactSearch/ppwBndryContactSearch.cpp Normal file
View File

@ -0,0 +1,109 @@
#include "ppwBndryContactSearch.hpp"
#include "ppwBndryContactSearchKernels.hpp"
#include "phasicFlowKokkos.hpp"
#include "streams.hpp"
void pFlow::ppwBndryContactSearch::checkAllocateNext(uint32 n)
{
if( nextCapacity_ < n)
{
nextCapacity_ = n;
reallocNoInit(next_, n);
}
}
void pFlow::ppwBndryContactSearch::nullifyHead()
{
fill(head_, static_cast<uint32>(-1));
}
void pFlow::ppwBndryContactSearch::nullifyNext(uint32 n)
{
fill(next_, 0u, n, static_cast<uint32>(-1));
}
void pFlow::ppwBndryContactSearch::buildList(
const deviceScatteredFieldAccess<realx3> &points)
{
if(points.empty())return;
uint32 n = points.size();
checkAllocateNext(n);
nullifyNext(n);
nullifyHead();
pFlow::pweBndryContactSearchKernels::buildNextHead(
points,
searchCells_,
head_,
next_
);
}
pFlow::ppwBndryContactSearch::ppwBndryContactSearch
(
const box &domain,
real cellSize,
real sizeRatio
)
:
searchCells_(domain, cellSize),
head_("periodic:head",searchCells_.nx(), searchCells_.ny(), searchCells_.nz()),
sizeRatio_(sizeRatio)
{
}
bool pFlow::ppwBndryContactSearch::broadSearchPP
(
csPairContainerType &ppPairs,
const deviceScatteredFieldAccess<realx3> &points,
const deviceScatteredFieldAccess<real>& diams,
const deviceScatteredFieldAccess<realx3> &mirrorPoints,
const deviceScatteredFieldAccess<real>& mirrorDiams,
const realx3& transferVec
)
{
buildList(points);
uint32 nNotInserted = 1;
// loop until the container size fits the numebr of contact pairs
while (nNotInserted > 0)
{
nNotInserted = pFlow::pweBndryContactSearchKernels::broadSearchPP
(
ppPairs,
points,
diams,
mirrorPoints,
mirrorDiams,
transferVec,
head_,
next_,
searchCells_,
sizeRatio_
);
if(nNotInserted)
{
// - resize the container
// note that getFull now shows the number of failed insertions.
uint32 len = max(nNotInserted,100u) ;
auto oldCap = ppPairs.capacity();
ppPairs.increaseCapacityBy(len);
INFORMATION<< "Particle-particle contact pair container capacity increased from "<<
oldCap << " to "<<ppPairs.capacity()<<" in peiodicBoundaryContactSearch."<<END_INFO;
}
}
return true;
}

85
src/Interaction/contactSearch/boundaries/periodicBoundaryContactSearch/ppwBndryContactSearch.hpp Normal file
View File

@ -0,0 +1,85 @@
/*------------------------------- 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 __ppwBndryContactSearch_hpp__
#define __ppwBndryContactSearch_hpp__
#include "contactSearchGlobals.hpp"
#include "scatteredFieldAccess.hpp"
#include "cells.hpp"
namespace pFlow
{
class ppwBndryContactSearch
{
public:
using HeadType = deviceViewType3D<uint32>;
using NextType = deviceViewType1D<uint32>;
private:
cells searchCells_;
HeadType head_{"periodic::head", 1, 1, 1};
NextType next_{"periodic::next", 1};
real sizeRatio_ = 1.0;
uint32 nextCapacity_ = 0;
void checkAllocateNext(uint32 n);
void nullifyHead();
void nullifyNext(uint32 n);
void buildList(
const deviceScatteredFieldAccess<realx3> &points);
public:
ppwBndryContactSearch(
const box &domain,
real cellSize,
real sizeRatio = 1.0);
bool broadSearchPP(
csPairContainerType &ppPairs,
const deviceScatteredFieldAccess<realx3> &points,
const deviceScatteredFieldAccess<real> &diams,
const deviceScatteredFieldAccess<realx3> &mirrorPoints,
const deviceScatteredFieldAccess<real> &mirrorDiams,
const realx3 &transferVec);
const auto& searchCells()const
{
return searchCells_;
}
real sizeRatio()const
{
return sizeRatio_;
}
};
}
#endif //__ppwBndryContactSearch_hpp__

105
src/Interaction/contactSearch/boundaries/periodicBoundaryContactSearch/ppwBndryContactSearchKernels.cpp Normal file
View File

@ -0,0 +1,105 @@
#include "ppwBndryContactSearchKernels.hpp"
INLINE_FUNCTION_HD
bool sphereSphereCheckB(const pFlow::realx3& p1, const pFlow::realx3 p2, pFlow::real d1, pFlow::real d2)
{
return pFlow::length(p2-p1) < 0.5*(d2+d1);
}
void pFlow::pweBndryContactSearchKernels::buildNextHead
(
const deviceScatteredFieldAccess<realx3> &points,
const cells &searchCells,
deviceViewType3D<uint32> &head,
deviceViewType1D<uint32> &next
)
{
if(points.empty())return;
uint32 n= points.size();
Kokkos::parallel_for(
"pFlow::ppwBndryContactSearch::buildList",
deviceRPolicyStatic(0,n),
LAMBDA_HD(uint32 i)
{
int32x3 ind;
if( searchCells.pointIndexInDomain(points[i], ind) )
{
// discards points out of searchCell
uint32 old = Kokkos::atomic_exchange(&head(ind.x(), ind.y(), ind.z()), i);
next[i] = old;
}
}
);
Kokkos::fence();
}
pFlow::uint32 pFlow::pweBndryContactSearchKernels::broadSearchPP
(
csPairContainerType &ppPairs,
const deviceScatteredFieldAccess<realx3> &points,
const deviceScatteredFieldAccess<real> &diams,
const deviceScatteredFieldAccess<realx3> &mirrorPoints,
const deviceScatteredFieldAccess<real> &mirrorDiams,
const realx3 &transferVec,
const deviceViewType3D<uint32> &head,
const deviceViewType1D<uint32> &next,
const cells &searchCells,
const real sizeRatio
)
{
if(points.empty()) return 0;
if(mirrorPoints.empty())return 0;
auto nMirror = mirrorPoints.size();
uint32 getFull = 0;
Kokkos::parallel_reduce(
"pFlow::pweBndryContactSearchKernels::broadSearchPP",
deviceRPolicyStatic(0, nMirror),
LAMBDA_HD(const uint32 mrrI, uint32 &getFullUpdate)
{
realx3 p_m = mirrorPoints(mrrI) + transferVec;
int32x3 ind_m;
if( !searchCells.pointIndexInDomain(p_m, ind_m))return;
real d_m = sizeRatio*mirrorDiams[mrrI];
for(int ii=-1; ii<2; ii++)
{
for(int jj=-1; jj<2; jj++)
{
for(int kk =-1; kk<2; kk++)
{
auto ind = ind_m + int32x3{ii,jj,kk};
if(!searchCells.inCellRange(ind))continue;
uint32 thisI = head(ind.x(),ind.y(),ind.z());
while (thisI!=-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)
{
getFullUpdate++;
}
thisI = next(thisI);
}
}
}
}
},
getFull
);
return getFull;
}

31
src/Interaction/contactSearch/boundaries/periodicBoundaryContactSearch/ppwBndryContactSearchKernels.hpp Normal file
View File

@ -0,0 +1,31 @@
#include "contactSearchGlobals.hpp"
#include "cells.hpp"
#include "contactSearchFunctions.hpp"
#include "scatteredFieldAccess.hpp"
namespace pFlow::pweBndryContactSearchKernels
{
void buildNextHead(
const deviceScatteredFieldAccess<realx3> &points,
const cells &searchCells,
deviceViewType3D<uint32> &head,
deviceViewType1D<uint32> &next );
uint32 broadSearchPP
(
csPairContainerType &ppPairs,
const deviceScatteredFieldAccess<realx3> &points,
const deviceScatteredFieldAccess<real> &diams,
const deviceScatteredFieldAccess<realx3> &mirrorPoints,
const deviceScatteredFieldAccess<real> &mirrorDiams,
const realx3 &transferVec,
const deviceViewType3D<uint32> &head,
const deviceViewType1D<uint32> &next,
const cells &searchCells,
real sizeRatio
);
}

141
src/Interaction/contactSearch/boundaries/periodicBoundaryContactSearch/wallBoundaryContactSearch.cpp Normal file
View File

@ -0,0 +1,141 @@
#include "wallBoundaryContactSearch.hpp"
#include "streams.hpp"
pFlow::wallBoundaryContactSearch::wallBoundaryContactSearch
(
real cellExtent,
uint32 numPoints,
uint32 numElements,
const ViewType1D<realx3, memory_space> &points,
const ViewType1D<uint32x3, memory_space> &vertices,
const ViewType1D<realx3, memory_space> &normals
)
:
cellExtent_( max(cellExtent, 0.5 ) ),
numElements_(numElements),
numPoints_(numPoints),
vertices_(vertices),
points_(points),
normals_(normals)
{
allocateArrays();
}
bool pFlow::wallBoundaryContactSearch::build(const cells &searchBox, const realx3& transferVec)
{
Kokkos::parallel_for(
"pFlow::cellsWallLevel0::build",
deviceRPolicyStatic(0,numElements_),
CLASS_LAMBDA_HD(uint32 i)
{
auto v = vertices_[i];
auto p1 = points_[v.x()]+transferVec;
auto p2 = points_[v.y()]+transferVec;
auto p3 = points_[v.z()]+transferVec;
realx3 minP;
realx3 maxP;
searchBox.extendBox(p1, p2, p3, cellExtent_, minP, maxP);
elementBox_[i] = iBoxType(searchBox.pointIndex(minP), searchBox.pointIndex(maxP));
auto d = elementBox_[i].maxPoint()-elementBox_[i].minPoint();
validBox_[i] = (d.x()*d.y()*d.z())==0? 0:1;
});
Kokkos::fence();
return true;
}
bool pFlow::wallBoundaryContactSearch::broadSearch
(
csPairContainerType &pairs,
const cells &searchCells,
const deviceScatteredFieldAccess<realx3> &thisPoints,
const deviceScatteredFieldAccess<real> &thisDiams,
const deviceScatteredFieldAccess<realx3> &mirrorPoints,
const deviceScatteredFieldAccess<real> &mirroDiams,
const realx3 &transferVec,
real sizeRatio
)
{
uint32 nNotInserted = 1;
while (nNotInserted>0u)
{
build(searchCells,{0,0,0});
nNotInserted = findPairsElementRangeCount(
pairs,
searchCells,
thisPoints,
thisDiams,
{0,0,0},
0
);
build(searchCells, transferVec);
nNotInserted += findPairsElementRangeCount(
pairs,
searchCells,
mirrorPoints,
mirroDiams,
transferVec,
BASE_MIRROR_WALL_INDEX
);
if(nNotInserted>0u)
{
// note that getFull now shows the number of failed insertions.
uint32 incCap = max(nNotInserted,50u) ;
auto oldCap = pairs.capacity();
pairs.increaseCapacityBy(incCap);
INFORMATION<< "The contact pair container capacity increased from "<<
oldCap << " to "<<pairs.capacity()<<" in wallBoundaryContactSearch."<<END_INFO;
}
}
return true;
}
pFlow::uint32 pFlow::wallBoundaryContactSearch::findPairsElementRangeCount
(
csPairContainerType &pairs,
const cells &searchCells,
const deviceScatteredFieldAccess<realx3> &pPoints,
const deviceScatteredFieldAccess<real> &pDiams,
const realx3 &transferVec,
uint baseTriIndex
)
{
if(pPoints.empty())return 0u;
uint32 nNotInserted = 0;
uint32 nThis = pPoints.size();
Kokkos::parallel_reduce(
"pFlow::wallBoundaryContactSearch::findPairsElementRangeCount",
deviceRPolicyDynamic(0,nThis),
LAMBDA_HD(uint32 i, uint32 &notInsertedUpdate)
{
auto p = pPoints[i]+transferVec;
int32x3 ind;
if( searchCells.pointIndexInDomain(p, ind) )
{
for(uint32 nTri=0; nTri<numElements_; nTri++)
{
if( validBox_[nTri]== 0)continue;
if( elementBox_[nTri].isInside(ind)&&
pairs.insert(i,nTri+baseTriIndex) == -1)
{
notInsertedUpdate++;
}
}
}
},
nNotInserted
);
return nNotInserted;
}

131
src/Interaction/contactSearch/boundaries/periodicBoundaryContactSearch/wallBoundaryContactSearch.hpp Normal file
View File

@ -0,0 +1,131 @@
/*------------------------------- 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 __wallBoundaryContactSearch_hpp__
#define __wallBoundaryContactSearch_hpp__
#include "contactSearchGlobals.hpp"
#include "contactSearchFunctions.hpp"
#include "scatteredFieldAccess.hpp"
#include "iBox.hpp"
#include "cells.hpp"
namespace pFlow
{
class wallBoundaryContactSearch
{
public:
using execution_space = csExecutionSpace;
using memory_space = typename execution_space::memory_space;
using iBoxType = iBox<int32>;
private:
// - box extent
real cellExtent_ = 0.5;
// - number of triangle elements
uint32 numElements_ = 0;
// - number of points
uint32 numPoints_ = 0;
// - ref to vectices (borrowed)
ViewType1D<uint32x3, memory_space> vertices_;
// - ref to points in the trisurface (borrowed)
ViewType1D<realx3, memory_space> points_;
// - ref to normal vectors of triangles (borrowed)
ViewType1D<realx3, memory_space> normals_;
// cell range of element/triangle bounding box
ViewType1D<iBoxType, memory_space> elementBox_;
ViewType1D<uint8, memory_space> validBox_;
FUNCTION_H
void allocateArrays()
{
reallocNoInit( elementBox_, numElements_);
reallocNoInit( validBox_, numElements_);
}
public:
TypeInfoNV("wallBoundaryContactSearch");
INLINE_FUNCTION_HD
wallBoundaryContactSearch()=default;
FUNCTION_H
wallBoundaryContactSearch(
real cellExtent,
uint32 numPoints,
uint32 numElements,
const ViewType1D<realx3,memory_space>& points,
const ViewType1D<uint32x3,memory_space>& vertices,
const ViewType1D<realx3, memory_space>& normals);
INLINE_FUNCTION_HD
uint32 numElements()const
{
return numElements_;
}
bool build(const cells& searchBox, const realx3& transferVec);
bool broadSearch(
csPairContainerType& pairs,
const cells &searchCells,
const deviceScatteredFieldAccess<realx3>& thisPoints,
const deviceScatteredFieldAccess<real>& thisDiams,
const deviceScatteredFieldAccess<realx3>& mirrorPoints,
const deviceScatteredFieldAccess<real>& mirroDiams,
const realx3& transferVec,
real sizeRatio);
uint32 findPairsElementRangeCount(
csPairContainerType &pairs,
const cells &searchCells,
const deviceScatteredFieldAccess<realx3> &pPoints,
const deviceScatteredFieldAccess<real> &pDiams,
const realx3 &transferVec,
uint baseTriIndex);
}; // wallBoundaryContactSearch
} // pFlow
#endif // __wallBoundaryContactSearch_hpp__

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

@ -19,54 +19,56 @@ Licence:
-----------------------------------------------------------------------------*/
#include "contactSearch.hpp"
#include "streams.hpp"
#include "particles.hpp"
pFlow::contactSearch::contactSearch(
const dictionary& dict,
const box& domain,
const box& extDomain,
const particles& prtcl,
const geometry& geom,
Timers& timers)
:
interactionBase(prtcl, geom),
domain_(domain),
dict_(dict),
sphereSphereTimer_("particle-particle contact search", &timers),
sphereWallTimer_("particle-wall contact search", &timers)
extendedDomainBox_(extDomain),
particles_(prtcl),
geometry_(geom),
ppTimer_("particle-particle contact search", &timers),
pwTimer_("particle-wall contact search", &timers),
dict_(dict)
{
}
const pFlow::pointStructure &pFlow::contactSearch::pStruct() const
{
return particles_.pStruct();
}
pFlow::uniquePtr<pFlow::contactSearch> pFlow::contactSearch::create(
const dictionary &dict,
const box& domain,
const box &extDomain,
const particles &prtcl,
const geometry &geom,
Timers &timers)
{
word baseMethName = dict.getVal<word>("method");
word wallMethod = dict.getVal<word>("wallMapping");
auto model = angleBracketsNames2("ContactSearch", baseMethName, wallMethod);
REPORT(1)<<"Selecting contact search model . . ."<<endREPORT;
auto model = angleBracketsNames("ContactSearch", baseMethName);
REPORT(1)<<"Selecting contact search model . . ."<<END_REPORT;
if( dictionaryvCtorSelector_.search(model))
{
auto objPtr = dictionaryvCtorSelector_[model] (dict, domain, prtcl, geom, timers);
REPORT(2)<<"Model "<< greenText(model)<<" is created."<<endREPORT;
auto objPtr = dictionaryvCtorSelector_[model] (dict, extDomain, prtcl, geom, timers);
REPORT(2)<<"Model "<< Green_Text(model)<<" is created."<< END_REPORT;
return objPtr;
}
else
{
printKeys
(
fatalError << "Ctor Selector "<< model << " dose not exist. \n"
fatalError << "Ctor Selector "<< model << " does not exist. \n"
<<"Avaiable ones are: \n\n"
,
dictionaryvCtorSelector_

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

@ -23,50 +23,44 @@ Licence:
#define __contactSearch_hpp__
#include "interactionBase.hpp"
#include "unsortedPairs.hpp"
#include "box.hpp"
#include "contactSearchGlobals.hpp"
#include "dictionary.hpp"
#include "virtualConstructor.hpp"
#include "Timer.hpp"
namespace pFlow
{
// - forward
class box;
class particles;
class geometry;
class pointStructure;
class contactSearch
:
public interactionBase
{
public:
using IdType = typename interactionBase::IdType;
private:
using IndexType = typename interactionBase::IndexType;
const box& extendedDomainBox_;
using ExecutionSpace = typename interactionBase::ExecutionSpace;
const particles& particles_;
using PairContainerType = unsortedPairs<ExecutionSpace, IdType>;
const geometry& geometry_;
protected:
Timer ppTimer_;
const box& domain_;
Timer pwTimer_;
dictionary dict_;
Timer sphereSphereTimer_;
Timer sphereWallTimer_;
auto& dict()
{
return dict_;
}
public:
TypeInfo("contactSearch");
contactSearch(
const dictionary& dict,
const box& domain,
const box& extDomain,
const particles& prtcl,
const geometry& geom,
Timers& timers);
@ -88,40 +82,68 @@ public:
(dict, domain, prtcl, geom, timers)
);
const auto& domain()const
{
return domain_;
}
const auto& dict()const
{
return dict_;
}
const auto& extendedDomainBox()const
{
return extendedDomainBox_;
}
const auto& Particles()const
{
return particles_;
}
const pointStructure& pStruct()const;
const auto& Geometry()const
{
return geometry_;
}
auto& ppTimer()
{
return ppTimer_;
}
auto& pwTimer()
{
return pwTimer_;
}
virtual
bool broadSearch(
PairContainerType& ppPairs,
PairContainerType& pwPairs,
uint32 iter,
real t,
real dt,
csPairContainerType& ppPairs,
csPairContainerType& pwPairs,
bool force = false) = 0;
virtual
bool ppEnterBroadSearch()const = 0;
bool boundaryBroadSearch(
uint32 i,
uint32 iter,
real t,
real dt,
csPairContainerType& ppPairs,
csPairContainerType& pwPairs,
bool force = false)=0;
virtual
bool pwEnterBroadSearch()const = 0;
bool enterBroadSearch(uint32 iter, real t, real dt)const = 0;
virtual
bool ppPerformedBroadSearch()const = 0;
virtual
bool pwPerformedBroadSearch()const = 0;
bool performedBroadSearch()const = 0;
static
uniquePtr<contactSearch> create(
const dictionary& dict,
const box& domain,
const box& extDomain,
const particles& prtcl,
const geometry& geom,
Timers& timers);

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

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

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

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

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

@ -1,209 +0,0 @@
/*------------------------------- phasicFlow ---------------------------------
O C enter of
O O E ngineering and
O O M ultiscale modeling of
OOOOOOO F luid flow
------------------------------------------------------------------------------
Copyright (C): www.cemf.ir
email: hamid.r.norouzi AT gmail.com
------------------------------------------------------------------------------
Licence:
This file is part of phasicFlow code. It is a free software for simulating
granular and multiphase flows. You can redistribute it and/or modify it under
the terms of GNU General Public License v3 or any other later versions.
phasicFlow is distributed to help others in their research in the field of
granular and multiphase flows, but WITHOUT ANY WARRANTY; without even the
implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
-----------------------------------------------------------------------------*/
#ifndef __NBS_hpp__
#define __NBS_hpp__
#include "NBSLevel0.hpp"
namespace pFlow
{
template<typename executionSpace>
class NBS
{
public:
using NBSLevel0Type = NBSLevel0<executionSpace>;
using cellIterator = typename NBSLevel0Type::cellIterator;
using IdType = typename NBSLevel0Type::IdType;
using IndexType = typename NBSLevel0Type::IndexType;
using Cells = typename NBSLevel0Type::Cells;
using CellType = typename Cells::CellType;
using execution_space = typename NBSLevel0Type::execution_space;
using memory_space = typename NBSLevel0Type::memory_space;
protected:
real sizeRatio_ = 1.0;
int32 updateFrequency_= 1;
int32 currentIter_ = 0;
bool performedSearch_ = false;
NBSLevel0Type NBSLevel0_;
private:
bool performSearch()
{
if(currentIter_ % updateFrequency_ == 0)
{
currentIter_++;
return true;
}else
{
currentIter_++;
return false;
}
}
public:
TypeInfoNV("NBS");
NBS(
const dictionary& dict,
const box& domain,
real minSize,
real maxSize,
const ViewType1D<realx3, memory_space>& position,
const ViewType1D<real, memory_space>& diam)
:
sizeRatio_(
max(
dict.getValOrSet<real>("sizeRatio", 1.0),
1.0
)),
updateFrequency_(
max(
dict.getValOrSet<int32>("updateFrequency", 1),
1
)),
NBSLevel0_(
domain,
maxSize*sizeRatio_,
sizeRatio_,
position,
diam)
{}
INLINE_FUNCTION_HD
NBS(const NBS&) = default;
INLINE_FUNCTION_HD
NBS& operator = (const NBS&) = default;
INLINE_FUNCTION_HD
~NBS()=default;
//// - Methods
bool enterBoadSearch()const
{
return currentIter_%updateFrequency_==0;
}
bool performedSearch()const
{
return performedSearch_;
}
Vector<cellIterator> getCellIteratorLevels()
{
return Vector<cellIterator>("cellIterator", 1, NBSLevel0_.getCellIterator());
}
auto getCellIterator(int32 lvl)const
{
return NBSLevel0_.getCellIterator();
}
int32 numLevels()const
{
return 1;
}
Vector<Cells> getCellsLevels()const
{
return Vector<Cells>("Cells", 1, NBSLevel0_.getCells());
}
auto getCells()const
{
return NBSLevel0_.getCells();
}
bool objectSizeChanged(int32 newSize)
{
NBSLevel0_.checkAllocateNext(newSize);
return true;
}
// - Perform the broad search to find pairs
// with force = true, perform broad search regardless of
// updateFrequency_ value
// on all the points in the range of [0,numPoints_)
template<typename PairsContainer>
bool broadSearch(PairsContainer& pairs, range activeRange, bool force=false)
{
if(force) currentIter_ = 0;
performedSearch_ = false;
if( !performSearch() ) return true;
NBSLevel0_.build(activeRange);
NBSLevel0_.findPairs(pairs);
performedSearch_ = true;
return true;
}
// - Perform the broad search to find pairs,
// ignore particles with incld(i) = true,
// with force = true, perform broad search regardless of
// updateFrequency_ value
template<typename PairsContainer, typename IncludeFunction>
bool broadSearch(PairsContainer& pairs, range activeRange, IncludeFunction incld, bool force = false)
{
if(force) currentIter_ = 0;
performedSearch_ = false;
if( !performSearch() ) return true;
NBSLevel0_.build(activeRange, incld);
NBSLevel0_.findPairs(pairs);
performedSearch_ = true;
return true;
}
};
}
#endif

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

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

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

@ -1,127 +0,0 @@
#ifndef __NBSLevel_hpp__
#define __NBSLevel_hpp__
#include "NBSLevel0.hpp"
namespace pFlow
{
INLINE_FUNCTION_HD
int32x3 mapIndexLevels(const int32x3& ind, int32 lowerLevel, int32 upperLevel);
template<typename executionSpace>
class
NBSLevel
:
public NBSLevel0<executionSpace>
{
public:
using NBSLevel0Type = NBSLevel0<executionSpace>;
using cellIterator = typename NBSLevel0Type::cellIterator;
using IdType = typename NBSLevel0Type::IdType;
using IndexType = typename NBSLevel0Type::IndexType;
using Cells = typename NBSLevel0Type::Cells;
using CellType = typename Cells::CellType;
using execution_space = typename NBSLevel0Type::execution_space;
using memory_space = typename NBSLevel0Type::memory_space;
using mdrPolicyFindPairs = typename NBSLevel0Type::mdrPolicyFindPairs;
using HeadType = typename NBSLevel0Type::HeadType;
using NextType = typename NBSLevel0Type::NextType;
template<typename exeSpace>
friend class NBSLevels;
protected:
int32 level_ = 0;
public:
TypeInfoNV("NBSLevel0");
INLINE_FUNCTION_HD
NBSLevel(){}
NBSLevel(
int32 lvl,
const box& domain,
real cellSize,
real sizeRatio,
const ViewType1D<realx3, memory_space>& position,
const ViewType1D<real, memory_space>& diam)
:
NBSLevel0Type(
domain,
cellSize,
sizeRatio,
position,
diam,
lvl==0),
level_(lvl)
{}
INLINE_FUNCTION_HD
NBSLevel(const NBSLevel&) = default;
INLINE_FUNCTION_HD
NBSLevel& operator = (const NBSLevel&) = default;
INLINE_FUNCTION_HD
~NBSLevel() = default;
INLINE_FUNCTION_HD
auto level()const
{
return level_;
}
template<typename PairsContainer>
INLINE_FUNCTION_H
int32 findPairsCountCross(PairsContainer& pairs, NBSLevel& upperLevel)
{
mdrPolicyFindPairs
mdrPolicy(
{0,0,0},
{this->nx(),this->ny(),this->nz()} );
int32 notInsertedPairs;
Kokkos::parallel_reduce (
"NBSLevel::findPairsCountCross",
mdrPolicy,
CLASS_LAMBDA_HD(int32 i, int32 j, int32 k, int32& getFullUpdate){
#include "NBSCrossLoop.hpp"
}, notInsertedPairs);
return notInsertedPairs;
}
};
INLINE_FUNCTION_HD
int32x3 mapIndexLevels( const int32x3& ind, int32 lowerLevel, int32 upperLevel)
{
int32 a = pow(2, static_cast<int32>(upperLevel-lowerLevel));
return ind/a;
}
}
#endif

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

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

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

@ -1,438 +0,0 @@
#ifndef __NBSLevels_hpp__
#define __NBSLevels_hpp__
#include "NBSLevel.hpp"
#include "NBSLevel0.hpp"
#include "KokkosTypes.hpp"
namespace pFlow
{
template<typename executionSpace>
class NBSLevels
{
public:
using NBSLevelType = NBSLevel<executionSpace>;
using cellIterator = typename NBSLevelType::cellIterator;
using IdType = typename NBSLevelType::IdType;
using IndexType = typename NBSLevelType::IndexType;
using Cells = typename NBSLevelType::Cells;
using CellType = typename Cells::CellType;
using execution_space = typename NBSLevelType::execution_space;
using memory_space = typename NBSLevelType::memory_space;
using realRange = kPair<real,real>;
protected:
real minSize_;
real maxSize_;
int32 numLevels_=1;
Vector<NBSLevelType> nbsLevels_;
ViewType1D<realRange, memory_space> sizeRangeLevels_;
ViewType1D<realRange, HostSpace> sizeRangeLevelsHost_;
ViewType1D<real, memory_space> maxSizeLevels_;
ViewType1D<real, HostSpace> maxSizeLevelsHost_;
ViewType1D<int8, memory_space> particleLevel_;
range activeRange_{0,0};
using rangePolicyType =
Kokkos::RangePolicy<
Kokkos::IndexType<int32>,
Kokkos::Schedule<Kokkos::Static>,
execution_space>;
int32 setNumLevels()
{
int32 maxOvermin = static_cast<int32>(maxSize_/minSize_);
if (maxOvermin <=1)
return 1;
else if(maxOvermin<=3)
return 2;
else if(maxOvermin<=7)
return 3;
else if(maxOvermin<15)
return 4;
else if(maxOvermin<31)
return 5;
else if(maxOvermin<63)
return 6;
else if(maxOvermin <127)
return 7;
else
{
fatalErrorInFunction<<
"size ratio is not valid for multi-grid NBS "<< maxOvermin<<endl;
fatalExit;
}
return -1;
}
bool setDiameterRange(real sizeRatio)
{
real lvl_maxD = sizeRatio* maxSize_;
real lvl_minD = lvl_maxD/2.0;
for(int32 lvl=numLevels_-1; lvl>=0; lvl--)
{
if(lvl == 0 ) lvl_minD = 0.01*minSize_;
sizeRangeLevelsHost_[lvl] = {lvl_minD, lvl_maxD};
maxSizeLevelsHost_[lvl] = lvl_maxD;
lvl_maxD = lvl_minD;
lvl_minD /= 2.0;
}
copy(sizeRangeLevels_, sizeRangeLevelsHost_);
copy(maxSizeLevels_, maxSizeLevelsHost_);
REPORT(2)<<"Grids with "<< yellowText(numLevels_)<< " levels have been created."<<endREPORT;
for(int32 lvl=0; lvl<numLevels_; lvl++)
{
REPORT(3)<<"Cell gird No "<< yellowText(lvl)<<" with size range ("
<<sizeRangeLevelsHost_[lvl].first<<","<<sizeRangeLevelsHost_[lvl].second<<"]."<<endREPORT;
}
return true;
}
bool initLevels(
const box& domain,
real sizeRatio,
const ViewType1D<realx3, memory_space>& position,
const ViewType1D<real, memory_space>& diam)
{
for(int32 lvl = 0; lvl<numLevels_; lvl++)
{
nbsLevels_[lvl] = NBSLevelType(
lvl,
domain,
maxSizeLevelsHost_[lvl]*sizeRatio,
sizeRatio,
position,
diam );
}
auto next0 = nbsLevels_[0].next();
for(int32 lvl=1; lvl<numLevels_; lvl++)
{
nbsLevels_[lvl].setNext(next0);
}
return true;
}
void manageAllocateNext(range active)
{
activeRange_ = active;
if(activeRange_.second > nbsLevels_[0].capacity())
{
nbsLevels_[0].checkAllocateNext(activeRange_.second);
auto next0 = nbsLevels_[0].next();
for(int32 lvl=1; lvl<numLevels_; lvl++)
{
nbsLevels_[lvl].setNext(next0);
}
}
}
void nullify( range active)
{
for(int32 lvl=0; lvl<numLevels_; lvl++)
{
nbsLevels_[lvl].nullify(active);
}
}
public:
NBSLevels(
const box& domain,
real minSize,
real maxSize,
real sizeRatio,
const ViewType1D<realx3, memory_space>& position,
const ViewType1D<real, memory_space>& diam)
:
minSize_(minSize),
maxSize_(maxSize),
numLevels_(setNumLevels()),
nbsLevels_("nbsLevels", numLevels_, numLevels_, RESERVE()),
sizeRangeLevels_("sizeRangeLevels", numLevels_),
sizeRangeLevelsHost_("sizeRangeLevelsHost", numLevels_),
maxSizeLevels_("maxSizeLevels", numLevels_),
maxSizeLevelsHost_("maxSizeLevelsHost", numLevels_)
{
setDiameterRange(sizeRatio);
initLevels(domain, sizeRatio, position, diam);
}
auto getCellIterator(int32 lvl)const
{
return nbsLevels_[lvl].getCellIterator();
}
int32 numLevels()const
{
return numLevels_;
}
Cells getCells(int32 lvl)const
{
return nbsLevels_[lvl].getCells();
}
template<typename PairsContainer>
INLINE_FUNCTION_H
bool findPairs(PairsContainer& pairs)
{
int32 getFull = 1;
// loop until the container size fits the numebr of contact pairs
while (getFull > 0)
{
getFull = findPairsCount(pairs);
if(getFull)
{
// - resize the container
// note that getFull now shows the number of failed insertions.
uint32 len = max(getFull,100) ;
auto oldCap = pairs.capacity();
pairs.increaseCapacityBy(len);
INFORMATION<< "The contact pair container capacity increased from "<<
oldCap << " to "<<pairs.capacity()<<" in NBSLevels."<<endINFO;
}
Kokkos::fence();
}
return true;
}
template<typename PairsContainer>
INLINE_FUNCTION_H
int32 findPairsCount(PairsContainer& pairs)
{
int32 notInsertedCount = 0;
for(int32 lvl=0; lvl<numLevels_; lvl++)
{
// the same level
notInsertedCount+= nbsLevels_[lvl].findPairsCount(pairs);
for(int32 crsLvl = lvl+1; crsLvl<numLevels_; crsLvl++)
{
// cross levels
notInsertedCount+=
nbsLevels_[lvl].findPairsCountCross(pairs, nbsLevels_[crsLvl]);
}
}
return notInsertedCount;
}
INLINE_FUNCTION_H
void build(range activeRange)
{
// nullify next and heads
findParticleLevel(activeRange.first, activeRange.second);
manageAllocateNext(activeRange);
nullify(activeRange);
//
rangePolicyType rPolicy(activeRange.first, activeRange.second);
auto nbsLevel0 = nbsLevels_[0];
auto pointPosition = nbsLevel0.pointPosition();
auto particleLevel = particleLevel_;
auto next = nbsLevel0.next();
auto head0 = nbsLevel0.head();
typename NBSLevelType::HeadType head1, head2, head3, head4, head5, head6;
if(numLevels_>1) head1 = nbsLevels_[1].head();
if(numLevels_>2) head2 = nbsLevels_[2].head();
if(numLevels_>3) head3 = nbsLevels_[3].head();
if(numLevels_>4) head4 = nbsLevels_[4].head();
if(numLevels_>5) head5 = nbsLevels_[5].head();
if(numLevels_>6) head6 = nbsLevels_[6].head();
Kokkos::parallel_for(
"NBSLevels::build",
rPolicy,
LAMBDA_HD(int32 i){
int8 lvl = particleLevel[i];
auto ind = nbsLevel0.pointIndex(pointPosition[i]);
ind = mapIndexLevels(ind, 0, lvl);
int32 old;
if(lvl==0)
old =Kokkos::atomic_exchange(&head0(ind.x(), ind.y(), ind.z()),i);
else if(lvl==1)
old =Kokkos::atomic_exchange(&head1(ind.x(), ind.y(), ind.z()),i);
else if(lvl==2)
old =Kokkos::atomic_exchange(&head2(ind.x(), ind.y(), ind.z()),i);
else if(lvl==3)
old =Kokkos::atomic_exchange(&head3(ind.x(), ind.y(), ind.z()),i);
else if(lvl==4)
old =Kokkos::atomic_exchange(&head4(ind.x(), ind.y(), ind.z()),i);
else if(lvl==5)
old =Kokkos::atomic_exchange(&head5(ind.x(), ind.y(), ind.z()),i);
else if(lvl==6)
old =Kokkos::atomic_exchange(&head6(ind.x(), ind.y(), ind.z()),i);
next(i) = old;
});
Kokkos::fence();
}
template<typename IncludeFunction>
INLINE_FUNCTION_H
void build(range activeRange, IncludeFunction incld)
{
// nullify next and heads
findParticleLevel(activeRange.first, activeRange.second);
manageAllocateNext(activeRange);
nullify(activeRange);
rangePolicyType rPolicy(activeRange.first, activeRange.second);
auto nbsLevel0 = nbsLevels_[0];
auto pointPosition = nbsLevel0.pointPosition();
auto particleLevel = particleLevel_;
auto next = nbsLevel0.next();
auto head0 = nbsLevel0.head();
typename NBSLevelType::HeadType head1, head2, head3, head4, head5, head6;
if(numLevels_>1) head1 = nbsLevels_[1].head();
if(numLevels_>2) head2 = nbsLevels_[2].head();
if(numLevels_>3) head3 = nbsLevels_[3].head();
if(numLevels_>4) head4 = nbsLevels_[4].head();
if(numLevels_>5) head5 = nbsLevels_[5].head();
if(numLevels_>6) head6 = nbsLevels_[6].head();
Kokkos::parallel_for(
"NBSLevels::build",
rPolicy,
LAMBDA_HD(int32 i){
if(!incld(i)) return;
int8 lvl = particleLevel[i];
auto ind = nbsLevel0.pointIndex(pointPosition[i]);
ind = mapIndexLevels(ind, 0, lvl);
int32 old;
if(lvl==0)
old =Kokkos::atomic_exchange(&head0(ind.x(), ind.y(), ind.z()),i);
else if(lvl==1)
old =Kokkos::atomic_exchange(&head1(ind.x(), ind.y(), ind.z()),i);
else if(lvl==2)
old =Kokkos::atomic_exchange(&head2(ind.x(), ind.y(), ind.z()),i);
else if(lvl==3)
old =Kokkos::atomic_exchange(&head3(ind.x(), ind.y(), ind.z()),i);
else if(lvl==4)
old =Kokkos::atomic_exchange(&head4(ind.x(), ind.y(), ind.z()),i);
else if(lvl==5)
old =Kokkos::atomic_exchange(&head5(ind.x(), ind.y(), ind.z()),i);
else if(lvl==6)
old =Kokkos::atomic_exchange(&head6(ind.x(), ind.y(), ind.z()),i);
next(i) = old;
});
Kokkos::fence();
}
bool findParticleLevel(int32 first, int32 last)
{
if(last > particleLevel_.size())
{
reallocNoInit(particleLevel_,last);
}
auto diameter = nbsLevels_[0].diameter();
auto const maxSizes = maxSizeLevels_;
auto particleLevel = particleLevel_;
auto const sizeRatio = 0.999*nbsLevels_[0].sizeRatio();
int8 maxLvl = sizeRangeLevels_.size();
rangePolicyType rPolicy(first, last);
Kokkos::parallel_for(
"NBSLevels::findParticleLevel",
rPolicy,
LAMBDA_HD(int32 i)
{
for(int8 lvl = 0; lvl<maxLvl; lvl++)
{
if( sizeRatio*diameter[i]<= maxSizes[lvl] )
{
particleLevel[i] = lvl;
return;
}
}
particleLevel[i] = static_cast<int8>(-1);
});
Kokkos::fence();
return true;
}
}; //NBSLevels
}
#endif

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

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

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

@ -0,0 +1,69 @@
/*------------------------------- 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 "NBS.hpp"
pFlow::NBS::NBS
(
const dictionary& dict,
const box& domainBox,
real minBSSize,
real maxBSSize,
const deviceViewType1D<realx3> &position,
const pFlagTypeDevice &flags,
const deviceViewType1D<real> &diam,
uint32 nWallPoints,
uint32 nWallElements,
const ViewType1D<realx3,memory_space>& wallPoints,
const ViewType1D<uint32x3,memory_space>& wallVertices,
const ViewType1D<realx3,memory_space>& wallNormals
)
:
particleWallContactSearchs<NBS>(
dict,
domainBox,
minBSSize,
maxBSSize,
position,
flags,
diam),
sizeRatio_(max(dict.getVal<real>("sizeRatio"), 1.0)),
cellExtent_(max(dict.getVal<real>("cellExtent"), 0.5)),
adjustableBox_(dict.getVal<Logical>("adjustableBox")),
NBSLevel0_
(
this->domainBox_,
maxBSSize,
sizeRatio_,
position,
flags,
adjustableBox_()
),
cellsWallLevel0_
(
cellExtent_,
nWallPoints,
nWallElements,
wallPoints,
wallVertices,
wallNormals
)
{
}

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

@ -0,0 +1,152 @@
/*------------------------------- phasicFlow ---------------------------------
O C enter of
O O E ngineering and
O O M ultiscale modeling of
OOOOOOO F luid flow
------------------------------------------------------------------------------
Copyright (C): www.cemf.ir
email: hamid.r.norouzi AT gmail.com
------------------------------------------------------------------------------
Licence:
This file is part of phasicFlow code. It is a free software for simulating
granular and multiphase flows. You can redistribute it and/or modify it under
the terms of GNU General Public License v3 or any other later versions.
phasicFlow is distributed to help others in their research in the field of
granular and multiphase flows, but WITHOUT ANY WARRANTY; without even the
implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
-----------------------------------------------------------------------------*/
#ifndef __NBS_hpp__
#define __NBS_hpp__
#include "particleWallContactSearchs.hpp"
#include "NBSLevel0.hpp"
#include "cellsWallLevel0.hpp"
#include "Vector.hpp"
namespace pFlow
{
class NBS
:
public particleWallContactSearchs<NBS>
{
public:
using CellIterator = typename NBSLevel0::CellIterator;
private:
real sizeRatio_ = 1.0;
real cellExtent_ = 0.5;
Logical adjustableBox_;
NBSLevel0 NBSLevel0_;
cellsWallLevel0 cellsWallLevel0_;
protected:
friend particleWallContactSearchs<NBS>;
bool impl_broadSearch
(
csPairContainerType& ppPairs,
csPairContainerType& pwPairs,
const deviceViewType1D<realx3>& pointPos,
const pFlagTypeDevice& flags,
const deviceViewType1D<real>& diameter
)
{
bool searchBoxChanged;
if( !NBSLevel0_.broadSearch(
ppPairs,
pointPos,
flags,
diameter,
searchBoxChanged))
{
fatalErrorInFunction<<
"Error in broadSearch for NBS (particle-particle)"<<endl;
return false;
}
if(!cellsWallLevel0_.broadSearch(
pwPairs,
NBSLevel0_.getSearchCells(),
NBSLevel0_.getCellIterator(),
pointPos,
diameter,
sizeRatio_))
{
fatalErrorInFunction<<
"Error in broadSearch for NBS (particle-wall)"<<endl;
return false;
}
return true;
}
public:
TypeInfoNV("NBS");
NBS(
const dictionary& dict,
const box& domainBox,
real minBSSize,
real maxBSSize,
const deviceViewType1D<realx3> &position,
const pFlagTypeDevice &flags,
const deviceViewType1D<real> &diam,
uint32 nWallPoints,
uint32 nWallElements,
const ViewType1D<realx3,memory_space>& wallPoints,
const ViewType1D<uint32x3,memory_space>& wallVertices,
const ViewType1D<realx3,memory_space>& wallNormals );
INLINE_FUNCTION_HD
NBS(const NBS&) = default;
INLINE_FUNCTION_HD
NBS(NBS&&) = default;
INLINE_FUNCTION_HD
NBS& operator = (const NBS&) = default;
INLINE_FUNCTION_HD
NBS& operator = (NBS&&) = default;
INLINE_FUNCTION_HD
~NBS()=default;
//// - Methods
uint32 numLevels()const
{
return 1;
}
auto getCellIterator([[maybe_unused]] uint32 lvl)const
{
return NBSLevel0_.getCellIterator();
}
Vector<cells> getDomainCellsLevels()const
{
return Vector<cells>("Cells", 1, NBSLevel0_.getDomainCells());
}
};
}
#endif

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

@ -0,0 +1,117 @@
/*------------------------------- phasicFlow ---------------------------------
O C enter of
O O E ngineering and
O O M ultiscale modeling of
OOOOOOO F luid flow
------------------------------------------------------------------------------
Copyright (C): www.cemf.ir
email: hamid.r.norouzi AT gmail.com
------------------------------------------------------------------------------
Licence:
This file is part of phasicFlow code. It is a free software for simulating
granular and multiphase flows. You can redistribute it and/or modify it under
the terms of GNU General Public License v3 or any other later versions.
phasicFlow is distributed to help others in their research in the field of
granular and multiphase flows, but WITHOUT ANY WARRANTY; without even the
implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
-----------------------------------------------------------------------------*/
#include "NBSLevel0.hpp"
#include "streams.hpp"
bool pFlow::NBSLevel0::findPairs
(
csPairContainerType &pairs,
const deviceViewType1D<realx3> &pointPos,
const pFlagTypeDevice &flags,
const deviceViewType1D<real> &diameter
)
{
uint32 getFull = 1;
// loop until the container size fits the numebr of contact pairs
while (getFull > 0)
{
getFull = pFlow::NBSLevel0Kernels::findPairsCount
(
pairs,
sizeRatio_,
pointPos,
flags,
diameter,
getCellIterator()
);
if(getFull)
{
// - resize the container
// note that getFull now shows the number of failed insertions.
uint32 len = max(getFull,500u) ;
auto oldCap = pairs.capacity();
pairs.increaseCapacityBy(len);
INFORMATION<< "The contact pair container capacity increased from "<<
oldCap << " to "<<pairs.capacity()<<" in NBSLevel0."<<END_INFO;
}
}
return true;
}
pFlow::NBSLevel0::NBSLevel0
(
const box& domain,
real cellSize,
real sizeRatio,
const deviceViewType1D<realx3>& pointPos,
const pFlagTypeDevice& flags,
bool adjustableBox
)
:
mapperNBS
(
domain,
cellSize,
pointPos,
flags,
adjustableBox,
true
),
sizeRatio_(sizeRatio)
{
}
bool pFlow::NBSLevel0::broadSearch
(
csPairContainerType &pairs,
const deviceViewType1D<realx3> &pointPos,
const pFlagTypeDevice &flags,
const deviceViewType1D<real> &diameter,
bool& searchBoxChanged
)
{
if(!build(pointPos, flags, searchBoxChanged))
{
fatalErrorInFunction;
return false;
}
if(!findPairs(pairs, pointPos, flags, diameter))
{
fatalErrorInFunction;
return false;
}
return true;
}

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

@ -0,0 +1,98 @@
/*------------------------------- phasicFlow ---------------------------------
O C enter of
O O E ngineering and
O O M ultiscale modeling of
OOOOOOO F luid flow
------------------------------------------------------------------------------
Copyright (C): www.cemf.ir
email: hamid.r.norouzi AT gmail.com
------------------------------------------------------------------------------
Licence:
This file is part of phasicFlow code. It is a free software for simulating
granular and multiphase flows. You can redistribute it and/or modify it under
the terms of GNU General Public License v3 or any other later versions.
phasicFlow is distributed to help others in their research in the field of
granular and multiphase flows, but WITHOUT ANY WARRANTY; without even the
implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
-----------------------------------------------------------------------------*/
#ifndef __NBSLevel0_hpp__
#define __NBSLevel0_hpp__
#include "NBSLevel0Kernels.hpp"
#include "mapperNBS.hpp"
namespace pFlow
{
class NBSLevel0
:
public mapperNBS
{
public:
using MapperType = mapperNBS;
using CellIterator = typename MapperType::CellIterator;
using HeadType = typename MapperType::HeadType;
using NextType = typename MapperType::NextType;
private:
real sizeRatio_ = 1.0;
bool findPairs
(
csPairContainerType& pairs,
const deviceViewType1D<realx3>& pointPos,
const pFlagTypeDevice& flags,
const deviceViewType1D<real>& diameter
);
public:
TypeInfoNV("NBSLevel0");
INLINE_FUNCTION_HD
NBSLevel0() = default;
NBSLevel0(
const box& domain,
real cellSize,
real sizeRatio,
const deviceViewType1D<realx3>& pointPos,
const pFlagTypeDevice& flags,
bool adjustableBox);
INLINE_FUNCTION_HD
NBSLevel0(const NBSLevel0&) = default;
INLINE_FUNCTION_HD
NBSLevel0& operator = (const NBSLevel0&) = default;
INLINE_FUNCTION_HD
~NBSLevel0()=default;
bool broadSearch(
csPairContainerType& pairs,
const deviceViewType1D<realx3>& pointPos,
const pFlagTypeDevice& flags,
const deviceViewType1D<real>& diameter,
bool& searchBoxChanged);
};
} // pFlow
#endif // __NBSLevel0_hpp__

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

@ -0,0 +1,83 @@
/*------------------------------- phasicFlow ---------------------------------
O C enter of
O O E ngineering and
O O M ultiscale modeling of
OOOOOOO F luid flow
------------------------------------------------------------------------------
Copyright (C): www.cemf.ir
email: hamid.r.norouzi AT gmail.com
------------------------------------------------------------------------------
Licence:
This file is part of phasicFlow code. It is a free software for simulating
granular and multiphase flows. You can redistribute it and/or modify it under
the terms of GNU General Public License v3 or any other later versions.
phasicFlow is distributed to help others in their research in the field of
granular and multiphase flows, but WITHOUT ANY WARRANTY; without even the
implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
-----------------------------------------------------------------------------*/
#include "mapperNBS.hpp"
#include "contactSearchGlobals.hpp"
#include "streams.hpp"
namespace pFlow::NBSLevel0Kernels
{
using mdrPolicyFindPairs =
Kokkos::MDRangePolicy<
Kokkos::IndexType<uint32>,
Kokkos::Rank<3>,
Kokkos::Schedule<Kokkos::Dynamic>,
DefaultExecutionSpace>;
template<typename T>
INLINE_FUNCTION_HD
void Swap(T& x, T& y)
{
T tmp = x;
x = y;
y = tmp;
}
INLINE_FUNCTION_HD
bool sphereSphereCheck(const realx3& p1, const realx3 p2, real d1, real d2)
{
return length(p2-p1) < 0.5*(d2+d1);
}
inline
uint32 findPairsCount
(
csPairContainerType& pairs,
real sizeRatio,
const deviceViewType1D<realx3>& pointPos,
const pFlagTypeDevice& flags,
const deviceViewType1D<real>& diameter,
mapperNBS::CellIterator cellIter
)
{
auto nCells = cellIter.numCells();
mdrPolicyFindPairs
mdrPolicy(
{0,0,0},
{nCells.x(), nCells.y(), nCells.z()} );
uint32 notInsertedPairs = 0u;
Kokkos::parallel_reduce (
"pFlow::NBSLevel0Kernels::findPairsCount",
mdrPolicy,
LAMBDA_HD(uint32 i, uint32 j, uint32 k, uint32& getFullUpdate){
#include "NBSLoop.hpp"
}, notInsertedPairs);
return notInsertedPairs;
}
}

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

@ -0,0 +1,100 @@
/*------------------------------- phasicFlow ---------------------------------
O C enter of
O O E ngineering and
O O M ultiscale modeling of
OOOOOOO F luid flow
------------------------------------------------------------------------------
Copyright (C): www.cemf.ir
email: hamid.r.norouzi AT gmail.com
------------------------------------------------------------------------------
Licence:
This file is part of phasicFlow code. It is a free software for simulating
granular and multiphase flows. You can redistribute it and/or modify it under
the terms of GNU General Public License v3 or any other later versions.
phasicFlow is distributed to help others in their research in the field of
granular and multiphase flows, but WITHOUT ANY WARRANTY; without even the
implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
-----------------------------------------------------------------------------*/
uint32 m = cellIter.start(i,j,k);
int32x3 currentCell(i,j,k);
uint32 n = mapperNBS::NoPos;
while( m != mapperNBS::NoPos)
{
auto p_m = pointPos[m];
auto d_m = sizeRatio* diameter[m];
// the same cell
n = cellIter.getNext(m);
while(n != mapperNBS::NoPos)
{
auto p_n = pointPos[n];
auto d_n = sizeRatio*diameter[n];
if( sphereSphereCheck(p_m, p_n, d_m, d_n) )
{
auto ln = n;
auto lm = m;
if(lm>ln) Swap(lm,ln);
if( pairs.insert(lm,ln) == -1)
{
getFullUpdate++;
}
}
n = cellIter.getNext(n);;
}
// neighbor cells
int32x3 neighborCell;
for(uint32 ni=0u; ni<13; ni++)
{
if(ni==0) neighborCell = currentCell + int32x3( 0, 0,-1);
else if(ni==1) neighborCell = currentCell + int32x3(-1, 0,-1);
else if(ni==2) neighborCell = currentCell + int32x3(-1, 0, 0);
else if(ni==3) neighborCell = currentCell + int32x3(-1, 0, 1);
else if(ni==4) neighborCell = currentCell + int32x3( 0,-1,-1);
else if(ni==5) neighborCell = currentCell + int32x3( 0,-1, 0);
else if(ni==6) neighborCell = currentCell + int32x3( 0,-1, 1);
else if(ni==7) neighborCell = currentCell + int32x3(-1,-1,-1);
else if(ni==8) neighborCell = currentCell + int32x3(-1,-1, 0);
else if(ni==9) neighborCell = currentCell + int32x3(-1,-1, 1);
else if(ni==10) neighborCell = currentCell + int32x3( 1,-1,-1);
else if(ni==11) neighborCell = currentCell + int32x3( 1,-1, 0);
else if(ni==12) neighborCell = currentCell + int32x3( 1,-1, 1);
if( neighborCell.x()>=0 && neighborCell.y()>=0 && neighborCell.z()>=0 &&
neighborCell.x()<nCells.x() && neighborCell.y()<nCells.y() && neighborCell.z()<nCells.z() )
{
n = cellIter.start(neighborCell.x(), neighborCell.y(), neighborCell.z());
while(n != mapperNBS::NoPos)
{
auto p_n = pointPos[n];
auto d_n = sizeRatio*diameter[n];
if(sphereSphereCheck(p_m, p_n, d_m, d_n))
{
auto ln = n;
auto lm = m;
if(lm>ln) Swap(lm,ln);
if( pairs.insert(lm,ln) == -1)
{
getFullUpdate++;
}
}
n = cellIter.next(n);
}
}
}
m = cellIter.next(m);
}

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

@ -0,0 +1,209 @@
/*------------------------------- phasicFlow ---------------------------------
O C enter of
O O E ngineering and
O O M ultiscale modeling of
OOOOOOO F luid flow
------------------------------------------------------------------------------
Copyright (C): www.cemf.ir
email: hamid.r.norouzi AT gmail.com
------------------------------------------------------------------------------
Licence:
This file is part of phasicFlow code. It is a free software for simulating
granular and multiphase flows. You can redistribute it and/or modify it under
the terms of GNU General Public License v3 or any other later versions.
phasicFlow is distributed to help others in their research in the field of
granular and multiphase flows, but WITHOUT ANY WARRANTY; without even the
implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
-----------------------------------------------------------------------------*/
#include "cellsWallLevel0.hpp"
#include "streams.hpp"
pFlow::cellsWallLevel0::cellsWallLevel0
(
real cellExtent,
uint32 numPoints,
uint32 numElements,
const ViewType1D<realx3, memory_space> &points,
const ViewType1D<uint32x3, memory_space> &vertices,
const ViewType1D<realx3, memory_space>& normals
)
:
cellExtent_( max(cellExtent, 0.5 ) ),
numElements_(numElements),
numPoints_(numPoints),
vertices_(vertices),
points_(points),
normals_(normals)
{
allocateArrays();
}
bool pFlow::cellsWallLevel0::resetElements
(
uint32 numElements,
uint32 numPoints,
const ViewType1D<realx3, memory_space>& points,
const ViewType1D<uint32x3, memory_space>& vertices,
const ViewType1D<realx3, memory_space>& normals
)
{
numElements_ = numElements;
numPoints_ = numPoints;
points_ = points;
vertices_ = vertices;
normals_ = normals;
allocateArrays();
return true;
}
bool pFlow::cellsWallLevel0::broadSearch
(
csPairContainerType &pairs,
const cells& searchBox,
const mapperNBS::CellIterator &particleMap,
const deviceViewType1D<realx3>& pPoints,
const deviceViewType1D<real>& pDiams,
real sizeRatio
)
{
// map walls onto the cells
this->build(searchBox);
this->particleWallFindPairs(pairs, particleMap, pPoints, pDiams, sizeRatio);
return true;
}
bool pFlow::cellsWallLevel0::build(const cells & searchBox)
{
Kokkos::parallel_for(
"pFlow::cellsWallLevel0::build",
deviceRPolicyStatic(0,numElements_),
CLASS_LAMBDA_HD(uint32 i)
{
auto v = vertices_[i];
auto p1 = points_[v.x()];
auto p2 = points_[v.y()];
auto p3 = points_[v.z()];
realx3 minP;
realx3 maxP;
searchBox.extendBox(p1, p2, p3, cellExtent_, minP, maxP);
elementBox_[i] = iBoxType(searchBox.pointIndex(minP), searchBox.pointIndex(maxP));
});
Kokkos::fence();
return true;
}
bool pFlow::cellsWallLevel0::particleWallFindPairs
(
csPairContainerType &pairs,
const mapperNBS::CellIterator &particleMap,
const deviceViewType1D<realx3>& pPoints,
const deviceViewType1D<real>& pDiams,
real sizeRatio
)
{
uint32 getFull = 1;
while (getFull)
{
getFull = findPairsElementRangeCount(pairs, particleMap, pPoints, pDiams, sizeRatio);
if(getFull)
{
// - resize the container
// note that getFull now shows the number of failed insertions.
uint32 len = max(getFull, 50u);
auto oldCap = pairs.capacity();
pairs.increaseCapacityBy(len);
INFORMATION<<"Contact pair container capacity increased from "<<
oldCap << " to "
<< pairs.capacity() <<" in cellsWallLevel0."<<END_INFO;
Kokkos::fence();
}
}
return true;
}
pFlow::int32 pFlow::cellsWallLevel0::findPairsElementRangeCount
(
csPairContainerType &pairs,
const mapperNBS::CellIterator &particleMap,
const deviceViewType1D<realx3>& pPoints,
const deviceViewType1D<real>& pDiams,
real sizeRatio
)
{
uint32 getFull =0;
Kokkos::parallel_reduce(
"pFlow::cellsWallLevel0::findPairsElementRangeCount",
tpPWContactSearch(numElements_, Kokkos::AUTO),
LAMBDA_HD(
const typename tpPWContactSearch::member_type & teamMember,
uint32& valueToUpdate){
const uint32 iTri = teamMember.league_rank();
const auto triBox = elementBox_[iTri];
const auto triPlane = infinitePlane(
normals_[iTri],
points_[vertices_[iTri].x()]);
uint32 getFull2 = 0;
auto bExtent = boxExtent(triBox);
uint32 numCellBox = bExtent.x()*bExtent.y()*bExtent.z();
Kokkos::parallel_reduce(
Kokkos::TeamThreadRange( teamMember, numCellBox ),
[&] ( const uint32 linIndex, uint32 &innerUpdate )
{
int32x3 cell;
indexToCell(linIndex, triBox, cell);
uint32 n = particleMap.start(cell.x(),cell.y(),cell.z());
while( n != particleMap.NoPos)
{
// id is wall id the pair is (particle id, wall id)
if( abs(triPlane.pointFromPlane(pPoints[n]))< pDiams[n]*sizeRatio*cellExtent_)
{
if( pairs.insert(
static_cast<csIdType>(n),
static_cast<csIdType>(iTri) ) == -1 )
innerUpdate++;
}
n = particleMap.next(n);
}
},
getFull2
);
if ( teamMember.team_rank() == 0 ) valueToUpdate += getFull2;
},
getFull
);
return getFull;
}

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

@ -0,0 +1,153 @@
/*------------------------------- phasicFlow ---------------------------------
O C enter of
O O E ngineering and
O O M ultiscale modeling of
OOOOOOO F luid flow
------------------------------------------------------------------------------
Copyright (C): www.cemf.ir
email: hamid.r.norouzi AT gmail.com
------------------------------------------------------------------------------
Licence:
This file is part of phasicFlow code. It is a free software for simulating
granular and multiphase flows. You can redistribute it and/or modify it under
the terms of GNU General Public License v3 or any other later versions.
phasicFlow is distributed to help others in their research in the field of
granular and multiphase flows, but WITHOUT ANY WARRANTY; without even the
implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
-----------------------------------------------------------------------------*/
#ifndef __cellsWallLevel0_hpp__
#define __cellsWallLevel0_hpp__
#include "contactSearchGlobals.hpp"
#include "contactSearchFunctions.hpp"
#include "mapperNBS.hpp"
#include "iBox.hpp"
namespace pFlow
{
class cellsWallLevel0
{
public:
using execution_space = csExecutionSpace;
using memory_space = typename execution_space::memory_space;
using iBoxType = iBox<int32>;
class TagFindCellRange2{};
private:
// - box extent
real cellExtent_ = 0.5;
// - number of triangle elements
uint32 numElements_ = 0;
// - number of points
uint32 numPoints_ = 0;
// - ref to vectices (borrowed)
ViewType1D<uint32x3, memory_space> vertices_;
// - ref to points in the trisurface (borrowed)
ViewType1D<realx3, memory_space> points_;
// - ref to normal vectors of triangles (borrowed)
ViewType1D<realx3, memory_space> normals_;
// cell range of element/triangle bounding box
ViewType1D<iBoxType, memory_space> elementBox_;
using tpPWContactSearch = Kokkos::TeamPolicy<
execution_space,
Kokkos::Schedule<Kokkos::Dynamic>,
Kokkos::IndexType<uint32>>;
FUNCTION_H
void allocateArrays()
{
reallocNoInit( elementBox_, numElements_);
}
public:
TypeInfoNV("cellsWallLevel0");
INLINE_FUNCTION_HD
cellsWallLevel0(){}
FUNCTION_H
cellsWallLevel0(
real cellExtent,
uint32 numPoints,
uint32 numElements,
const ViewType1D<realx3,memory_space>& points,
const ViewType1D<uint32x3,memory_space>& vertices,
const ViewType1D<realx3, memory_space>& normals);
// - host call
// reset triangle elements if they have changed
bool resetElements(
uint32 numElements,
uint32 numPoints,
const ViewType1D<realx3, memory_space>& points,
const ViewType1D<uint32x3, memory_space>& vertices,
const ViewType1D<realx3, memory_space>& normals);
INLINE_FUNCTION_HD
iBoxType elementBox(uint32 i)const
{
return elementBox_[i];
}
INLINE_FUNCTION_HD
uint32 numElements()const
{
return numElements_;
}
bool broadSearch(
csPairContainerType& pairs,
const cells& searchBox,
const mapperNBS::CellIterator& particleMap,
const deviceViewType1D<realx3>& pPoints,
const deviceViewType1D<real>& pDiams,
real sizeRatio);
bool build(const cells& searchBox);
bool particleWallFindPairs(
csPairContainerType& pairs,
const mapperNBS::CellIterator& particleMap,
const deviceViewType1D<realx3>& pPoints,
const deviceViewType1D<real>& pDiams,
real sizeRatio);
int32 findPairsElementRangeCount(
csPairContainerType& pairs,
const mapperNBS::CellIterator& particleMap,
const deviceViewType1D<realx3>& pPoints,
const deviceViewType1D<real>& pDiams,
real sizeRatio);
}; // cellsWallLevel0
} // pFlow
#endif // __cellsWallLevel0_hpp__

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

@ -0,0 +1,194 @@
/*------------------------------- phasicFlow ---------------------------------
O C enter of
O O E ngineering and
O O M ultiscale modeling of
OOOOOOO F luid flow
------------------------------------------------------------------------------
Copyright (C): www.cemf.ir
email: hamid.r.norouzi AT gmail.com
------------------------------------------------------------------------------
Licence:
This file is part of phasicFlow code. It is a free software for simulating
granular and multiphase flows. You can redistribute it and/or modify it under
the terms of GNU General Public License v3 or any other later versions.
phasicFlow is distributed to help others in their research in the field of
granular and multiphase flows, but WITHOUT ANY WARRANTY; without even the
implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
-----------------------------------------------------------------------------*/
#include "mapperNBS.hpp"
#include "mapperNBSKernels.hpp"
#include "streams.hpp"
pFlow::uint32 pFlow::mapperNBS::checkInterval_ = 1000;
pFlow::real pFlow::mapperNBS::enlargementFactor_ = 1.1;
bool pFlow::mapperNBS::setSearchBox
(
const deviceViewType1D<realx3> &pointPos,
const pFlagTypeDevice &flags,
real cellSize
)
{
box domainBox = domainCells_.domainBox();
if(adjustableBox_)
{
lastCheckForBox_ = buildCount_;
realx3 minP;
realx3 maxP;
pFlow::mapperNBSKernels::findPointExtends
(
pointPos,
flags,
minP, maxP
);
minP = max( minP - enlargementFactor_*cellSize, domainBox.minPoint());
maxP = min( maxP + enlargementFactor_*cellSize, domainBox.maxPoint());
box searchBox = {minP, maxP};
searchCells_ = cells(searchBox, cellSize);
INFORMATION<<"Search box for contact search has changed: "<<
"search box is ["<<searchCells_.domainBox().minPoint()<<
" "<<searchCells_.domainBox().maxPoint()<<"]"<<END_INFO;
return true;
}
else
{
searchCells_ = cells(domainBox, cellSize);
INFORMATION<<"Search box for contact search is: ["
<< domainBox.minPoint()<<" "<<domainBox.maxPoint()<<"]"<<END_INFO;
return false;
}
}
void pFlow::mapperNBS::allocateArrays(rangeU32 nextRng)
{
checkAllocateNext(nextRng);
nullifyNext(nextRng);
reallocFill(head_, searchCells_.nx(), searchCells_.ny(), searchCells_.nz(), NoPos);
}
void pFlow::mapperNBS::checkAllocateNext(rangeU32 nextRng)
{
auto newCap = nextRng.end();
if( nextCapacity_ < newCap)
{
nextCapacity_ = newCap;
if(!nextOwner_)return;
reallocNoInit(next_, nextCapacity_);
}
}
void pFlow::mapperNBS::nullifyHead()
{
fill(head_, NoPos);
}
void pFlow::mapperNBS::nullifyNext(rangeU32 nextRng)
{
if(!nextOwner_)return;
fill(next_, nextRng, NoPos);
}
pFlow::mapperNBS::mapperNBS(
const box &domain,
real cellSize,
const deviceViewType1D<realx3> &pointPos,
const pFlagTypeDevice &flags,
bool adjustableBox,
bool nextOwner
)
:
domainCells_(domain, cellSize),
searchCells_(domain, cellSize),
adjustableBox_(adjustableBox),
nextOwner_(nextOwner)
{
setSearchBox(pointPos, flags, cellSize);
allocateArrays(flags.activeRange());
}
bool pFlow::mapperNBS::build
(
const deviceViewType1D<realx3>& pointPos,
const pFlagTypeDevice & flags,
bool& searchBoxChanged
)
{
auto aRange = flags.activeRange();
buildCount_++;
if(adjustableBox_ && buildCount_%checkInterval_ == 0)
{
if(searchBoxChanged =
setSearchBox(pointPos, flags, searchCells_.cellSize());searchBoxChanged)
{
allocateArrays(aRange);
}
}
else
{
checkAllocateNext(aRange);
nullifyHead();
nullifyNext(aRange);
}
if( adjustableBox_ )
{
if(!pFlow::mapperNBSKernels::buildListsReduce(
searchCells_,
head_,
next_,
pointPos,
flags) )
{
buildCount_++;
setSearchBox(pointPos, flags, searchCells_.cellSize());
searchBoxChanged = true;
allocateArrays(flags.activeRange());
if(!pFlow::mapperNBSKernels::buildListsReduce(
searchCells_,
head_,
next_,
pointPos,
flags))
{
fatalErrorInFunction<<"failed to build list in anjustable search box mode!"<<endl;
return false;
}
}
}
else
{
pFlow::mapperNBSKernels::buildLists(
searchCells_,
head_,
next_,
pointPos,
flags
);
searchBoxChanged = false;
}
return true;
}

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

@ -0,0 +1,212 @@
/*------------------------------- phasicFlow ---------------------------------
O C enter of
O O E ngineering and
O O M ultiscale modeling of
OOOOOOO F luid flow
------------------------------------------------------------------------------
Copyright (C): www.cemf.ir
email: hamid.r.norouzi AT gmail.com
------------------------------------------------------------------------------
Licence:
This file is part of phasicFlow code. It is a free software for simulating
granular and multiphase flows. You can redistribute it and/or modify it under
the terms of GNU General Public License v3 or any other later versions.
phasicFlow is distributed to help others in their research in the field of
granular and multiphase flows, but WITHOUT ANY WARRANTY; without even the
implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
-----------------------------------------------------------------------------*/
#ifndef __mapperNBS_hpp__
#define __mapperNBS_hpp__
#include "phasicFlowKokkos.hpp"
#include "pointFlag.hpp"
#include "cells.hpp"
//#include "contactSearchFunctions.hpp"
namespace pFlow
{
class mapperNBS
{
public:
using HeadType = deviceViewType3D<uint32>;
using NextType = deviceViewType1D<uint32>;
static constexpr uint32 NoPos = 0xFFFFFFFF;
class CellIterator
{
private:
HeadType head_;
NextType next_;
public:
CellIterator(const HeadType& head, const NextType& next)
:
head_(head),
next_(next)
{}
static constexpr uint32 NoPos = 0xFFFFFFFF;
INLINE_FUNCTION_HD
int32x3 numCells()const {
return int32x3(head_.extent(0), head_.extent(1), head_.extent(2));}
INLINE_FUNCTION_HD
uint32 start(int32 i, int32 j, int32 k)const {
return head_(i,j,k); }
INLINE_FUNCTION_HD
uint32 getNext(uint32 n)const {
if(n == NoPos ) return NoPos;
return next_(n); }
INLINE_FUNCTION_HD
uint32 next(uint32 n)const{
return next_(n);}
};
private:
cells domainCells_;
cells searchCells_;
HeadType head_{"NBS::head",1,1,1};
NextType next_{"NBS::next", 1};
uint32 nextCapacity_ = 0;
uint32 lastCheckForBox_ = 0;
uint32 buildCount_ = 0;
bool adjustableBox_ = false;
bool nextOwner_ = true;
static uint32 checkInterval_;
static real enlargementFactor_;
bool setSearchBox(
const deviceViewType1D<realx3>& pointPos,
const pFlagTypeDevice& flags,
real cellSize
);
void allocateArrays(rangeU32 nextRng);
void checkAllocateNext(rangeU32 nextRng);
void nullifyHead();
void nullifyNext(rangeU32 nextRng);
public:
TypeInfoNV("mapperNBS");
INLINE_FUNCTION_HD
mapperNBS() = default;
mapperNBS(
const box& domain,
real cellSize,
const deviceViewType1D<realx3>& pointPos,
const pFlagTypeDevice& flags,
bool adjustableBox,
bool nextOwner = true);
INLINE_FUNCTION_HD
mapperNBS(const mapperNBS&) = default;
INLINE_FUNCTION_HD
mapperNBS(mapperNBS&&) = default;
INLINE_FUNCTION_HD
mapperNBS& operator = (const mapperNBS&) = default;
INLINE_FUNCTION_HD
mapperNBS& operator = (mapperNBS&&) = default;
INLINE_FUNCTION_HD
~mapperNBS()=default;
//// - Methods
auto getCellIterator()const
{
return CellIterator(head_, next_);
}
const auto& getDomainCells()const
{
return domainCells_;
}
const auto& getSearchCells()const
{
return searchCells_;
}
bool build(
const deviceViewType1D<realx3>& pointPos,
const pFlagTypeDevice& flags,
bool& searchBoxChanged);
};
} // pFlow
#endif // __mapperNBS_hpp__
/*
INLINE_FUNCTION_HD
auto& head()
{
return head_;
}
INLINE_FUNCTION_HD
auto& next()
{
return next_;
}
INLINE_FUNCTION_HD
const auto& head()const
{
return head_;
}
INLINE_FUNCTION_HD
const auto& next()const
{
return next_;
}
INLINE_FUNCTION_HD
auto& pointPosition()
{
return pointPosition_;
}
*/

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

@ -0,0 +1,192 @@
/*------------------------------- phasicFlow ---------------------------------
O C enter of
O O E ngineering and
O O M ultiscale modeling of
OOOOOOO F luid flow
------------------------------------------------------------------------------
Copyright (C): www.cemf.ir
email: hamid.r.norouzi AT gmail.com
------------------------------------------------------------------------------
Licence:
This file is part of phasicFlow code. It is a free software for simulating
granular and multiphase flows. You can redistribute it and/or modify it under
the terms of GNU General Public License v3 or any other later versions.
phasicFlow is distributed to help others in their research in the field of
granular and multiphase flows, but WITHOUT ANY WARRANTY; without even the
implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
-----------------------------------------------------------------------------*/
#include "mapperNBSKernels.hpp"
void pFlow::mapperNBSKernels::findPointExtends
(
const deviceViewType1D<realx3>& points,
const pFlagTypeDevice& flags,
realx3& minPoint,
realx3& maxPoint
)
{
if(flags.numActive() == 0)
{
minPoint = {0,0,0};
maxPoint = {0,0,0};
return;
}
real minX;
real minY;
real minZ;
real maxX;
real maxY;
real maxZ;
auto aRange = flags.activeRange();
Kokkos::parallel_reduce(
"pFlow::mapperNBSKernels::findPointExtends",
deviceRPolicyStatic(aRange.start(), aRange.end()),
LAMBDA_HD(
uint32 i,
real& minXUpdate,
real& minYUpdate,
real& minZUpdate,
real& maxXUpdate,
real& maxYUpdate,
real& maxZUpdate)
{
if(flags(i))
{
auto p = points(i);
minXUpdate = min(p.x(), minXUpdate);
minYUpdate = min(p.y(), minYUpdate);
minZUpdate = min(p.z(), minZUpdate);
maxXUpdate = max(p.x(), maxXUpdate);
maxYUpdate = max(p.y(), maxYUpdate);
maxZUpdate = max(p.z(), maxZUpdate);
}
},
Kokkos::Min<real>(minX),
Kokkos::Min<real>(minY),
Kokkos::Min<real>(minZ),
Kokkos::Max<real>(maxX),
Kokkos::Max<real>(maxY),
Kokkos::Max<real>(maxZ)
);
minPoint = {minX, minY, minZ};
maxPoint = {maxX, maxY, maxZ};
}
bool pFlow::mapperNBSKernels::buildListsReduce
(
const cells &searchCell,
const deviceViewType3D<uint32> &head,
const deviceViewType1D<uint32> &next,
const deviceViewType1D<realx3> &points,
const pFlagTypeDevice &flags
)
{
uint32 numOut = 0;
auto aRange = flags.activeRange();
if(flags.isAllActive())
{
Kokkos::parallel_reduce
(
"pFlow::mapperNBSKernels::buildListsReduce",
deviceRPolicyStatic(aRange.start(), aRange.end()),
LAMBDA_HD(uint32 i, uint32& valToUpdate)
{
int32x3 ind;
if( searchCell.pointIndexInDomain(points[i], ind) )
{
uint32 old = Kokkos::atomic_exchange(&head(ind.x(), ind.y(), ind.z()), i);
next[i] = old;
}
else
{
valToUpdate++;
}
},
numOut
);
}
else
{
Kokkos::parallel_reduce
(
"pFlow::mapperNBSKernels::buildListsReduce",
deviceRPolicyStatic(aRange.start(), aRange.end()),
LAMBDA_HD(uint32 i, uint32& valToUpdate)
{
int32x3 ind;
if( flags(i) )
{
if( searchCell.pointIndexInDomain(points[i], ind) )
{
uint32 old = Kokkos::atomic_exchange(&head(ind.x(), ind.y(), ind.z()), i);
next[i] = old;
}
else
{
valToUpdate++;
}
}
},
numOut
);
}
return numOut == 0u ;
}
bool pFlow::mapperNBSKernels::buildLists
(
const cells &searchCell,
const deviceViewType3D<uint32> &head,
const deviceViewType1D<uint32> &next,
const deviceViewType1D<realx3> &points,
const pFlagTypeDevice &flags
)
{
auto aRange = flags.activeRange();
if(flags.isAllActive() )
{
Kokkos::parallel_for
(
"pFlow::mapperNBSKernels::buildLists",
deviceRPolicyStatic(aRange.start(), aRange.end()),
LAMBDA_HD(uint32 i)
{
auto ind = searchCell.pointIndex(points[i]);
uint32 old = Kokkos::atomic_exchange(&head(ind.x(), ind.y(), ind.z()), i);
next[i] = old;
}
);
Kokkos::fence();
}
else
{
Kokkos::parallel_for
(
"pFlow::mapperNBSKernels::buildLists",
deviceRPolicyStatic(aRange.start(), aRange.end()),
LAMBDA_HD(uint32 i)
{
if( flags(i) )
{
auto ind = searchCell.pointIndex(points[i]);
uint32 old = Kokkos::atomic_exchange(&head(ind.x(), ind.y(), ind.z()), i);
next[i] = old;
}
}
);
Kokkos::fence();
}
return true;
}

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

@ -0,0 +1,49 @@
/*------------------------------- 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 "phasicFlowKokkos.hpp"
#include "cells.hpp"
#include "pointFlag.hpp"
namespace pFlow::mapperNBSKernels
{
void findPointExtends(
const deviceViewType1D<realx3>& points,
const pFlagTypeDevice& flags,
realx3& minPoint,
realx3& maxPoint);
bool buildListsReduce(
const cells& searchCell,
const deviceViewType3D<uint32>& head,
const deviceViewType1D<uint32>& next,
const deviceViewType1D<realx3>& points,
const pFlagTypeDevice& flags);
bool buildLists(
const cells& searchCell,
const deviceViewType3D<uint32>& head,
const deviceViewType1D<uint32>& next,
const deviceViewType1D<realx3>& points,
const pFlagTypeDevice& flags);
}

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

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

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

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

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

@ -1,213 +0,0 @@
/*------------------------------- phasicFlow ---------------------------------
O C enter of
O O E ngineering and
O O M ultiscale modeling of
OOOOOOO F luid flow
------------------------------------------------------------------------------
Copyright (C): www.cemf.ir
email: hamid.r.norouzi AT gmail.com
------------------------------------------------------------------------------
Licence:
This file is part of phasicFlow code. It is a free software for simulating
granular and multiphase flows. You can redistribute it and/or modify it under
the terms of GNU General Public License v3 or any other later versions.
phasicFlow is distributed to help others in their research in the field of
granular and multiphase flows, but WITHOUT ANY WARRANTY; without even the
implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
-----------------------------------------------------------------------------*/
#ifndef __multiGridNBS_hpp__
#define __multiGridNBS_hpp__
#include "NBSLevels.hpp"
namespace pFlow
{
template<typename executionSpace>
class multiGridNBS
{
public:
using NBSLevelsType = NBSLevels<executionSpace>;
using cellIterator = typename NBSLevelsType::cellIterator;
using IdType = typename NBSLevelsType::IdType;
using IndexType = typename NBSLevelsType::IndexType;
using Cells = typename NBSLevelsType::Cells;
using CellType = typename Cells::CellType;
using execution_space = typename NBSLevelsType::execution_space;
using memory_space = typename NBSLevelsType::memory_space;
protected:
real sizeRatio_ = 1.0;
int32 updateFrequency_= 1;
int32 currentIter_ = 0;
bool performedSearch_ = false;
NBSLevelsType NBSLevels_;
private:
bool performSearch()
{
if(currentIter_ % updateFrequency_ == 0)
{
currentIter_++;
return true;
}else
{
currentIter_++;
return false;
}
}
public:
TypeInfoNV("multiGridNBS");
multiGridNBS(
const dictionary& dict,
const box& domain,
real minSize,
real maxSize,
const ViewType1D<realx3, memory_space>& position,
const ViewType1D<real, memory_space>& diam)
:
sizeRatio_(
max(
dict.getVal<real>("sizeRatio"),
1.0
)),
updateFrequency_(
max(
dict.getVal<int32>("updateFrequency"),
1
)),
NBSLevels_(
domain,
minSize,
maxSize,
sizeRatio_,
position,
diam)
{}
INLINE_FUNCTION_HD
multiGridNBS(const multiGridNBS&) = default;
INLINE_FUNCTION_HD
multiGridNBS& operator = (const multiGridNBS&) = default;
INLINE_FUNCTION_HD
~multiGridNBS()=default;
//// - Methods
bool enterBoadSearch()const
{
return currentIter_%updateFrequency_==0;
}
bool performedSearch()const
{
return performedSearch_;
}
int32 numLevels()const
{
return NBSLevels_.numLevels();
}
auto getCellsLevels()const
{
Vector<Cells> cellsLvl("cells", numLevels(), numLevels(), RESERVE());
for(int32 lvl=0; lvl<numLevels(); lvl++)
{
cellsLvl[lvl] = NBSLevels_.getCells(lvl) ;
}
return cellsLvl;
}
auto getCells(int32 lvl)const
{
return NBSLevels_.getCells(lvl);
}
auto getCellIterator(int32 lvl)const
{
return NBSLevels_.getCellIterator(lvl);
}
bool objectSizeChanged(int32 newSize)
{
NBSLevels_.checkAllocateNext(newSize);
return true;
}
// - Perform the broad search to find pairs
// with force = true, perform broad search regardless of
// updateFrequency_ value
// on all the points in the range of [0,numPoints_)
template<typename PairsContainer>
bool broadSearch(PairsContainer& pairs, range activeRange, bool force=false)
{
if(force) currentIter_ = 0;
performedSearch_ = false;
if( !performSearch() ) return true;
NBSLevels_.build(activeRange);
NBSLevels_.findPairs(pairs);
performedSearch_ = true;
return true;
}
// - Perform the broad search to find pairs,
// ignore particles with incld(i) = true,
// with force = true, perform broad search regardless of
// updateFrequency_ value
template<typename PairsContainer, typename IncludeFunction>
bool broadSearch(PairsContainer& pairs, range activeRange, IncludeFunction incld, bool force = false)
{
if(force) currentIter_ = 0;
performedSearch_ = false;
if( !performSearch() ) return true;
NBSLevels_.build(activeRange, incld);
NBSLevels_.findPairs(pairs);
performedSearch_ = true;
return true;
}
};
}
#endif

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

@ -0,0 +1,78 @@
/*------------------------------- phasicFlow ---------------------------------
O C enter of
O O E ngineering and
O O M ultiscale modeling of
OOOOOOO F luid flow
------------------------------------------------------------------------------
Copyright (C): www.cemf.ir
email: hamid.r.norouzi AT gmail.com
------------------------------------------------------------------------------
Licence:
This file is part of phasicFlow code. It is a free software for simulating
granular and multiphase flows. You can redistribute it and/or modify it under
the terms of GNU General Public License v3 or any other later versions.
phasicFlow is distributed to help others in their research in the field of
granular and multiphase flows, but WITHOUT ANY WARRANTY; without even the
implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
-----------------------------------------------------------------------------*/
template <typename method>
pFlow::particleWallContactSearchs<method>::particleWallContactSearchs
(
const dictionary &dict,
const box &domain,
real minSize,
real maxSize,
const ViewType1D<realx3, memory_space> &position,
const pFlagTypeDevice &flags,
const ViewType1D<real, memory_space> &diam
)
:
domainBox_(domain),
updateInterval_
(
max(dict.getValOrSet<uint32>("updateInterval", 1),1u)
)
{
}
template <typename method>
bool pFlow::particleWallContactSearchs<method>::broadSearch
(
uint32 iter,
real t,
real dt,
csPairContainerType& ppPairs,
csPairContainerType& pwPairs,
const deviceViewType1D<realx3>& pointPos,
const pFlagTypeDevice& flags,
const deviceViewType1D<real>& diameter,
bool force
)
{
performedSearch_ = false;
if( !performSearch(iter, force) ) return true;
if(!getMethod().impl_broadSearch(
ppPairs,
pwPairs,
pointPos,
flags,
diameter))
{
fatalErrorInFunction<<
"Error in performing particle-particle broadSearch in method"<<
getMethod().typeName()<<endl;
return false;
}
lastUpdated_ = iter;
performedSearch_ = true;
return true;
}

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

@ -0,0 +1,123 @@
/*------------------------------- phasicFlow ---------------------------------
O C enter of
O O E ngineering and
O O M ultiscale modeling of
OOOOOOO F luid flow
------------------------------------------------------------------------------
Copyright (C): www.cemf.ir
email: hamid.r.norouzi AT gmail.com
------------------------------------------------------------------------------
Licence:
This file is part of phasicFlow code. It is a free software for simulating
granular and multiphase flows. You can redistribute it and/or modify it under
the terms of GNU General Public License v3 or any other later versions.
phasicFlow is distributed to help others in their research in the field of
granular and multiphase flows, but WITHOUT ANY WARRANTY; without even the
implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
-----------------------------------------------------------------------------*/
#ifndef __particleWallContactSearchs_hpp__
#define __particleWallContactSearchs_hpp__
#include "pointFlag.hpp"
#include "contactSearchGlobals.hpp"
#include "box.hpp"
namespace pFlow
{
template<typename method>
class particleWallContactSearchs
{
public:
using MethodType = method;
using IndexType = uint32;
using execution_space = DefaultExecutionSpace;
using memory_space = typename execution_space::memory_space;
private:
// friend
friend MethodType;
/// @brief box enclosing the simulation domain (local to processor)
box domainBox_;
/*/// @brief box enclosing the area for contact search (region with points)
box searchBox_;*/
/// @brief update interval in terms of iteration numebr
uint32 updateInterval_= 1;
/// @brief last iteration number which contact search has been performed
uint32 lastUpdated_ = 0;
/// @brief performed search?
bool performedSearch_ = false;
protected:
inline
auto& getMethod()
{
return static_cast<MethodType&>(*this);
}
inline
const auto& getMethod()const
{
return static_cast<const MethodType&>(*this);
}
public:
particleWallContactSearchs(
const dictionary& dict,
const box& domain,
real minSize,
real maxSize,
const ViewType1D<realx3, memory_space>& position,
const pFlagTypeDevice &flags,
const ViewType1D<real, memory_space>& diam
);
bool broadSearch
(
uint32 iter,
real t,
real dt,
csPairContainerType& ppPairs,
csPairContainerType& pwPairs,
const deviceViewType1D<realx3>& pointPos,
const pFlagTypeDevice& flags,
const deviceViewType1D<real>& diameter,
bool force = false
);
bool performedSearch()const
{
return performedSearch_;
}
bool performSearch(uint32 iter, bool force = false)const
{
if((iter-lastUpdated_) % updateInterval_ == 0 || iter == 0 || force )
{
return true;
}
return false;
}
};
} // pFlow
#include "particleWallContactSearchs.cpp"
#endif //__particleWallContactSearchs_hpp__

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

@ -1,150 +0,0 @@
/*------------------------------- phasicFlow ---------------------------------
O C enter of
O O E ngineering and
O O M ultiscale modeling of
OOOOOOO F luid flow
------------------------------------------------------------------------------
Copyright (C): www.cemf.ir
email: hamid.r.norouzi AT gmail.com
------------------------------------------------------------------------------
Licence:
This file is part of phasicFlow code. It is a free software for simulating
granular and multiphase flows. You can redistribute it and/or modify it under
the terms of GNU General Public License v3 or any other later versions.
phasicFlow is distributed to help others in their research in the field of
granular and multiphase flows, but WITHOUT ANY WARRANTY; without even the
implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
-----------------------------------------------------------------------------*/
#ifndef __cellMapping_hpp__
#define __cellMapping_hpp__
#include "cellsWallLevel0.hpp"
#include "dictionary.hpp"
namespace pFlow
{
template<
typename executionSpace
>
class cellMapping
{
public:
using cellsWallLevel0Type = cellsWallLevel0<executionSpace>;
using IdType = typename cellsWallLevel0Type::IdType;
using IndexType = typename cellsWallLevel0Type::IndexType;
using Cells = typename cellsWallLevel0Type::Cells;
using CellType = typename Cells::CellType;
using execution_space = typename cellsWallLevel0Type::execution_space;
using memory_space = typename cellsWallLevel0Type::memory_space;
using iBoxType = iBox<IndexType>;
protected:
// - update frequency
int32 updateFrequency_=1;
real cellExtent_;
int32 currentIter_ = 0;
/// a broad search has been occured during last pass?
bool performedSearch_ = false;
cellsWallLevel0Type cellsWallLevle_;
private:
bool performSearch()
{
if(currentIter_ % updateFrequency_ == 0)
{
currentIter_++;
return true;
}else
{
currentIter_++;
return false;
}
}
public:
TypeInfoNV("cellMapping");
cellMapping(
const dictionary& dict,
int32 numLevels,
const Vector<Cells>& ppCells,
int32 numPoints,
int32 numElements,
const ViewType1D<realx3,memory_space>& points,
const ViewType1D<int32x3,memory_space>& vertices
)
:
updateFrequency_(
max(
dict.getValOrSet<int32>(
"updateFrequency",
1),
1)),
cellExtent_(
max(
dict.getValOrSet<real>(
"cellExtent",
0.5),
0.5)),
cellsWallLevle_(
ppCells[0],
cellExtent_,
numPoints,
numElements,
points,
vertices
)
{}
bool enterBoadSearch()const
{
return currentIter_%updateFrequency_==0;
}
bool performedSearch()const
{
return performedSearch_;
}
template<typename PairsContainer, typename particleMapType>
bool broadSearch(PairsContainer& pairs, particleMapType& particleMap, bool force=false)
{
if(force) currentIter_ = 0;
performedSearch_= false;
if(!performSearch())return true;
cellsWallLevle_.broadSearch(pairs, particleMap);
performedSearch_ = true;
return true;
}
}; // cellMapping
} // pFlow
#endif

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

@ -1,285 +0,0 @@
/*------------------------------- phasicFlow ---------------------------------
O C enter of
O O E ngineering and
O O M ultiscale modeling of
OOOOOOO F luid flow
------------------------------------------------------------------------------
Copyright (C): www.cemf.ir
email: hamid.r.norouzi AT gmail.com
------------------------------------------------------------------------------
Licence:
This file is part of phasicFlow code. It is a free software for simulating
granular and multiphase flows. You can redistribute it and/or modify it under
the terms of GNU General Public License v3 or any other later versions.
phasicFlow is distributed to help others in their research in the field of
granular and multiphase flows, but WITHOUT ANY WARRANTY; without even the
implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
-----------------------------------------------------------------------------*/
#ifndef __cellsWallLevel0_hpp__
#define __cellsWallLevel0_hpp__
#include "types.hpp"
#include "KokkosTypes.hpp"
#include "cells.hpp"
#include "iBox.hpp"
namespace pFlow
{
template<
typename executionSpace
>
class cellsWallLevel0
:
public cells<int32>
{
public:
using IdType = int32;
using IndexType = int32;
using Cells = cells<IndexType>;
using CellType = typename Cells::CellType;
using execution_space = executionSpace;
using memory_space = typename execution_space::memory_space;
using iBoxType = iBox<IndexType>;
class TagFindCellRange2{};
protected:
// - box extent
real cellExtent_ = 0.5;
// - number of triangle elements
int32 numElements_ = 0;
// - number of points
int32 numPoints_ = 0;
// - ref to vectices (borrowed)
ViewType1D<int32x3, memory_space> vertices_;
// - ref to points in the trisurface (borrowed)
ViewType1D<realx3, memory_space> points_;
// cell range of element/triangle bounding box
ViewType1D<iBoxType, memory_space> elementBox_;
using tpPWContactSearch = Kokkos::TeamPolicy<
execution_space,
Kokkos::Schedule<Kokkos::Dynamic>,
Kokkos::IndexType<int32>
>;
using rpFindCellRange2Type =
Kokkos::RangePolicy<TagFindCellRange2, execution_space, Kokkos::IndexType<int32>>;
FUNCTION_H
void allocateArrays()
{
reallocNoInit( elementBox_, numElements_);
}
public:
TypeInfoNV("cellsWallLevel0");
INLINE_FUNCTION_HD
cellsWallLevel0(){}
FUNCTION_H
cellsWallLevel0(
const Cells& ppCells,
real cellExtent,
int32 numPoints,
int32 numElements,
const ViewType1D<realx3,memory_space>& points,
const ViewType1D<int32x3,memory_space>& vertices
)
:
Cells(ppCells),
cellExtent_( max(cellExtent, 0.5 ) ),
numElements_(numElements),
numPoints_(numPoints),
vertices_(vertices),
points_(points)
{
allocateArrays();
}
// - host call
// reset triangle elements if they have changed
bool resetElements(
int32 numElements,
int32 numPoints,
ViewType1D<realx3, memory_space>& points,
ViewType1D<int32x3, memory_space>& vertices )
{
numElements_ = numElements;
numPoints_ = numPoints;
points_ = points;
vertices_ = vertices;
allocateArrays();
return true;
}
INLINE_FUNCTION_HD
iBoxType elementBox(int32 i)const
{
return elementBox_[i];
}
INLINE_FUNCTION_HD
int32 numElements()const
{
return numElements_;
}
template<typename PairsContainer, typename particleMapType>
bool broadSearch(PairsContainer& pairs, particleMapType& particleMap)
{
// map walls onto the cells
this->build();
this->particleWallFindPairs(pairs, particleMap);
return true;
}
bool build()
{
Kokkos::parallel_for(
"cellsSimple::findcellrange2",
rpFindCellRange2Type(0,numElements_),
*this);
Kokkos::fence();
return true;
}
template<typename PairsContainer, typename particleMapType>
bool particleWallFindPairs(PairsContainer& pairs, particleMapType& particleMap)
{
int32 getFull = 1;
while (getFull)
{
getFull = findPairsElementRangeCount(pairs, particleMap.getCellIterator(0));
if(getFull)
{
// - resize the container
// note that getFull now shows the number of failed insertions.
uint32 len = max(getFull, 50);
auto oldCap = pairs.capacity();
pairs.increaseCapacityBy(len);
INFORMATION<<"Contact pair container capacity increased from "<<
oldCap << " to "
<< pairs.capacity() <<" in cellsWallLevel0."<<endINFO;
Kokkos::fence();
}
}
return true;
}
template<typename PairsContainer, typename CellIteratorType>
int32 findPairsElementRangeCount(PairsContainer& pairs, CellIteratorType cellIter)
{
int32 getFull =0;
const auto pwPairs = pairs;
const auto elementBox = elementBox_;
Kokkos::parallel_reduce(
"cellsSimple::findPairsElementRangeModified2",
tpPWContactSearch(numElements_, Kokkos::AUTO),
LAMBDA_HD(
const typename tpPWContactSearch::member_type & teamMember,
int32& valueToUpdate){
const int32 iTri = teamMember.league_rank();
const auto triBox = elementBox[iTri];
int32 getFull2 = 0;
auto bExtent = boxExtent(triBox);
int32 numCellBox = bExtent.x()*bExtent.y()*bExtent.z();
Kokkos::parallel_reduce(
Kokkos::TeamThreadRange( teamMember, numCellBox ),
[&] ( const int32 linIndex, int32 &innerUpdate )
{
CellType cell;
indexToCell(linIndex, triBox, cell);
int32 n = cellIter.start(cell.x(),cell.y(),cell.z());
while( n>-1)
{
// id is wall id the pair is (particle id, wall id)
if( pairs.insert(static_cast<IdType>(n), iTri) < 0 )
innerUpdate++;
n = cellIter.getNext(n);
}
},
getFull2
);
if ( teamMember.team_rank() == 0 ) valueToUpdate += getFull2;
},
getFull
);
return getFull;
}
INLINE_FUNCTION_HD
void operator()(TagFindCellRange2, int32 i) const
{
auto v = vertices_[i];
auto p1 = points_[v.x()];
auto p2 = points_[v.y()];
auto p3 = points_[v.z()];
realx3 minP, maxP;
this->extendBox(p1, p2, p3, cellExtent_, minP, maxP);
elementBox_[i] = iBoxType(this->pointIndex(minP), this->pointIndex(maxP));
}
}; // cellsWallLevel0
} // pFlow
#endif // __cellsWallLevel0_hpp__

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

@ -1,152 +0,0 @@
/*------------------------------- phasicFlow ---------------------------------
O C enter of
O O E ngineering and
O O M ultiscale modeling of
OOOOOOO F luid flow
------------------------------------------------------------------------------
Copyright (C): www.cemf.ir
email: hamid.r.norouzi AT gmail.com
------------------------------------------------------------------------------
Licence:
This file is part of phasicFlow code. It is a free software for simulating
granular and multiphase flows. You can redistribute it and/or modify it under
the terms of GNU General Public License v3 or any other later versions.
phasicFlow is distributed to help others in their research in the field of
granular and multiphase flows, but WITHOUT ANY WARRANTY; without even the
implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
-----------------------------------------------------------------------------*/
#ifndef __cellsWallLevels_hpp__
#define __cellsWallLevels_hpp__
#include "cellsWallLevel0.hpp"
namespace pFlow
{
template<
typename executionSpace
>
class cellsWallLevels
{
public:
using cellsWallLevel0Type = cellsWallLevel0<executionSpace>;
using IdType = typename cellsWallLevel0Type::IdType;
using IndexType = typename cellsWallLevel0Type::IndexType;
using Cells = typename cellsWallLevel0Type::Cells;
using CellType = typename Cells::CellType;
using execution_space = typename cellsWallLevel0Type::execution_space;
using memory_space = typename cellsWallLevel0Type::memory_space;
using iBoxType = iBox<IndexType>;
protected:
int32 numLevles_=1;
Vector<cellsWallLevel0Type> cellsWallLevels_;
public:
TypeInfoNV("cellsWallLevels");
FUNCTION_H
cellsWallLevels(
int32 numLevels,
const Vector<Cells>& cellsLevels,
real cellExtent,
int32 numPoints,
int32 numElements,
const ViewType1D<realx3,memory_space>& points,
const ViewType1D<int32x3,memory_space>& vertices
)
:
numLevles_(numLevels),
cellsWallLevels_("cellsWallLevels",numLevels, numLevels, RESERVE())
{
for(int32 lvl=0; lvl<numLevles_; lvl++)
{
cellsWallLevels_[lvl] =
cellsWallLevel0Type(
cellsLevels[lvl],
cellExtent,
numPoints,
numElements,
points,
vertices);
}
}
template<typename PairsContainer, typename particleMapType>
bool broadSearch(PairsContainer& pairs, particleMapType& particleMap)
{
// map walls onto the cells
for(int32 lvl=0; lvl<numLevles_; lvl++)
{
cellsWallLevels_[lvl].build();
}
this->particleWallFindPairs(pairs, particleMap);
return true;
}
template<typename PairsContainer, typename particleMapType>
bool particleWallFindPairs(PairsContainer& pairs, particleMapType& particleMap)
{
int32 getFull = 1;
while (getFull)
{
getFull = 0;
for(int32 lvl=0; lvl<numLevles_; lvl++)
{
getFull +=
cellsWallLevels_[lvl].findPairsElementRangeCount(
pairs,
particleMap.getCellIterator(lvl));
}
if(getFull)
{
// - resize the container
// note that getFull now shows the number of failed insertions.
uint32 len = max(getFull, 5000);
auto oldCap = pairs.capacity();
pairs.increaseCapacityBy(len);
INFORMATION<<"Contact pair container capacity increased from "<<
oldCap << " to "
<< pairs.capacity() <<" in cellsWallLevels."<<endINFO;
Kokkos::fence();
}
}
return true;
}
}; // cellsWallLevels
} // pFlow
#endif // __cellsWallLevels_hpp__

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

@ -1,153 +0,0 @@
/*------------------------------- phasicFlow ---------------------------------
O C enter of
O O E ngineering and
O O M ultiscale modeling of
OOOOOOO F luid flow
------------------------------------------------------------------------------
Copyright (C): www.cemf.ir
email: hamid.r.norouzi AT gmail.com
------------------------------------------------------------------------------
Licence:
This file is part of phasicFlow code. It is a free software for simulating
granular and multiphase flows. You can redistribute it and/or modify it under
the terms of GNU General Public License v3 or any other later versions.
phasicFlow is distributed to help others in their research in the field of
granular and multiphase flows, but WITHOUT ANY WARRANTY; without even the
implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
-----------------------------------------------------------------------------*/
#ifndef __multiGridMapping_hpp__
#define __multiGridMapping_hpp__
#include "cellsWallLevels.hpp"
#include "dictionary.hpp"
namespace pFlow
{
template<
typename executionSpace
>
class multiGridMapping
{
public:
using CellsWallLevelType = cellsWallLevels<executionSpace>;
using IdType = typename CellsWallLevelType::IdType;
using IndexType = typename CellsWallLevelType::IndexType;
using Cells = typename CellsWallLevelType::Cells;
using CellType = typename Cells::CellType;
using execution_space = typename CellsWallLevelType::execution_space;
using memory_space = typename CellsWallLevelType::memory_space;
using iBoxType = iBox<IndexType>;
protected:
// - update frequency
int32 updateFrequency_=1;
real cellExtent_;
int32 currentIter_ = 0;
/// a broad search has been occured during last pass?
bool performedSearch_ = false;
CellsWallLevelType cellsWallLevle_;
private:
bool performSearch()
{
if(currentIter_ % updateFrequency_ == 0)
{
currentIter_++;
return true;
}else
{
currentIter_++;
return false;
}
}
public:
TypeInfoNV("multiGridMapping");
multiGridMapping(
const dictionary& dict,
int32 numLevels,
const Vector<Cells>& ppCells,
int32 numPoints,
int32 numElements,
const ViewType1D<realx3,memory_space>& points,
const ViewType1D<int32x3,memory_space>& vertices
)
:
updateFrequency_(
max(
dict.getVal<int32>("updateFrequency"),
1)),
cellExtent_(
max(
dict.getVal<real>("cellExtent"),
0.5)),
cellsWallLevle_(
numLevels,
ppCells,
cellExtent_,
numPoints,
numElements,
points,
vertices
)
{
REPORT(3)<<"Multi-grid wall mapping with "<<
yellowText(numLevels)<<" levels has been created."<<endREPORT;
}
bool enterBoadSearch()const
{
return currentIter_%updateFrequency_==0;
}
bool performedSearch()const
{
return performedSearch_;
}
template<typename PairsContainer, typename particleMapType>
bool broadSearch(PairsContainer& pairs, particleMapType& particleMap, bool force=false)
{
if(force) currentIter_ = 0;
performedSearch_= false;
if(!performSearch())return true;
cellsWallLevle_.broadSearch(pairs, particleMap);
performedSearch_ = true;
return true;
}
}; // multiGridMapping
} // pFlow
#endif

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

@ -19,7 +19,9 @@ Licence:
-----------------------------------------------------------------------------*/
#include "interaction.hpp"
#include "particles.hpp"
#include "vocabs.hpp"
#include "geometry.hpp"
pFlow::interaction::interaction
(
@ -28,26 +30,16 @@ pFlow::interaction::interaction
const geometry& geom
)
:
demInteraction(control, control.caseSetup().path()+interactionFile__),
interactionBase(prtcl, geom),
fileDict_(control.caseSetup().emplaceObject<dictionary>(
objectFile(
interactionFile__,
"",
objectFile::READ_ALWAYS,
objectFile::WRITE_NEVER),
interactionFile__,
true ))
property(interactionFile__, &control.caseSetup()),
observer
(
&prtcl.dynPointStruct(),
msg_
),
demComponent("interaction", control),
particles_(prtcl),
geometry_(geom)
{
this->subscribe(prtcl.pStruct());
contactSearch_ = contactSearch::create(
fileDict_.subDict("contactSearch"),
this->control().domain(),
prtcl,
geom,
timers()
);
}
@ -79,14 +71,15 @@ pFlow::uniquePtr<pFlow::interaction> pFlow::interaction::create
clType);
REPORT(1)<< "Selecting interaction model..."<<endREPORT;
REPORT(1)<<"Creating interaction "<<Green_Text(interactionModel)<<" . . ."<<END_REPORT;
if( systemControlvCtorSelector_.search(interactionModel) )
{
auto objPtr =
systemControlvCtorSelector_[interactionModel]
(control, prtcl, geom);
REPORT(2)<<"Model "<<greenText(interactionModel)<<" is created."<<endREPORT;
return objPtr;
}
else

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

@ -21,9 +21,10 @@ Licence:
#ifndef __interaction_hpp__
#define __interaction_hpp__
#include "demInteraction.hpp"
#include "eventObserver.hpp"
#include "interactionBase.hpp"
#include "demComponent.hpp"
#include "property.hpp"
#include "observer.hpp"
#include "systemControl.hpp"
#include "contactSearch.hpp"
@ -31,27 +32,26 @@ Licence:
namespace pFlow
{
class particles;
class geometry;
class interaction
:
public demInteraction,
public eventObserver,
public interactionBase
public property,
public observer,
public demComponent
{
public:
private:
using IdType = typename interactionBase::IdType;
/// reference to particles object
const particles& particles_;
using IndexType = typename interactionBase::IndexType;
/// reference to geometry object
const geometry& geometry_;
using ExecutionSpace = typename interactionBase::ExecutionSpace;
protected:
/// interaction file dictionary
dictionary& fileDict_;
/// contact search object for pp and pw interactions
uniquePtr<contactSearch> contactSearch_ = nullptr;
static inline
const message msg_ = message::ITEM_REARRANGE;
public:
@ -65,11 +65,10 @@ public:
const particles& prtcl,
const geometry& geom );
~interaction() override = default;
virtual ~interaction() = default;
create_vCtor(
create_vCtor
(
interaction,
systemControl,
(
@ -80,19 +79,16 @@ public:
(control, prtcl, geom)
);
auto& contactSearchPtr()
inline
const auto& Particles()const
{
return contactSearch_;
return particles_;
}
auto& contactSearchRef()
inline
const auto& Geometry()const
{
return contactSearch_();
}
const auto& fileDict()const
{
return fileDict_;
return geometry_;
}
static

91
src/Interaction/sphereInteraction/boundaries/boundarySphereInteraction/boundarySphereInteraction.cpp Normal file
View File

@ -0,0 +1,91 @@
/*------------------------------- phasicFlow ---------------------------------
O C enter of
O O E ngineering and
O O M ultiscale modeling of
OOOOOOO F luid flow
------------------------------------------------------------------------------
Copyright (C): www.cemf.ir
email: hamid.r.norouzi AT gmail.com
------------------------------------------------------------------------------
Licence:
This file is part of phasicFlow code. It is a free software for simulating
granular and multiphase flows. You can redistribute it and/or modify it under
the terms of GNU General Public License v3 or any other later versions.
phasicFlow is distributed to help others in their research in the field of
granular and multiphase flows, but WITHOUT ANY WARRANTY; without even the
implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
-----------------------------------------------------------------------------*/
template <typename cFM, typename gMM>
pFlow::boundarySphereInteraction<cFM, gMM>::boundarySphereInteraction(
const boundaryBase &boundary,
const sphereParticles &sphPrtcls,
const GeometryMotionModel &geomMotion)
: generalBoundary(boundary, sphPrtcls.pStruct(), "", ""),
geometryMotion_(geomMotion),
sphParticles_(sphPrtcls)
{
ppPairs_ = makeUnique<ContactListType>(1);
pwPairs_ = makeUnique<ContactListType>(1);
}
template <typename cFM, typename gMM>
pFlow::uniquePtr<pFlow::boundarySphereInteraction<cFM, gMM>>
pFlow::boundarySphereInteraction<cFM, gMM>::create(
const boundaryBase &boundary,
const sphereParticles &sphPrtcls,
const GeometryMotionModel &geomMotion)
{
word cfTypeName = ContactForceModel::TYPENAME();
word gmTypeName = MotionModel::TYPENAME();
word bType = boundary.type();
word boundaryTypeName = angleBracketsNames3(
"boundarySphereInteraction",
bType,
cfTypeName,
gmTypeName);
word altBTypeName = angleBracketsNames2(
"boundarySphereInteraction",
cfTypeName,
gmTypeName);
if (boundaryBasevCtorSelector_.search(boundaryTypeName))
{
REPORT(2) << "Creating boundry type " << Green_Text(boundaryTypeName) <<
" for boundary " << boundary.name() << " . . ." << END_REPORT;
return boundaryBasevCtorSelector_[boundaryTypeName](
boundary,
sphPrtcls,
geomMotion);
}
else if(boundaryBasevCtorSelector_[altBTypeName])
{
// if boundary condition is not implemented, the default is used
REPORT(2) << "Creating boundry type " << Green_Text(altBTypeName) <<
" for boundary " << boundary.name() << " . . ." << END_REPORT;
return boundaryBasevCtorSelector_[altBTypeName](
boundary,
sphPrtcls,
geomMotion);
}
else
{
printKeys
(
fatalError << "Ctor Selector "<< boundaryTypeName<<
" and "<< altBTypeName << " do not exist. \n"
<<"Avaiable ones are: \n\n"
,
boundaryBasevCtorSelector_
);
fatalExit;
}
return nullptr;
}

168
src/Interaction/sphereInteraction/boundaries/boundarySphereInteraction/boundarySphereInteraction.hpp Normal file
View File

@ -0,0 +1,168 @@
/*------------------------------- 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 __boundarySphereInteraction_hpp__
#define __boundarySphereInteraction_hpp__
#include "virtualConstructor.hpp"
#include "generalBoundary.hpp"
#include "unsortedContactList.hpp"
#include "sphereParticles.hpp"
namespace pFlow
{
template<typename contactForceModel,typename geometryMotionModel>
class boundarySphereInteraction
:
public generalBoundary
{
public:
using BoundarySphereInteractionType
= boundarySphereInteraction<contactForceModel, geometryMotionModel>;
using GeometryMotionModel = geometryMotionModel;
using ContactForceModel = contactForceModel;
using MotionModel = typename geometryMotionModel::MotionModel;
using ModelStorage = typename ContactForceModel::contactForceStorage;
using IdType = uint32;
using IndexType = uint32;
using ContactListType =
unsortedContactList<ModelStorage, DefaultExecutionSpace, IdType>;
private:
const GeometryMotionModel& geometryMotion_;
/// const reference to sphere particles
const sphereParticles& sphParticles_;
uniquePtr<ContactListType> ppPairs_;
uniquePtr<ContactListType> pwPairs_;
public:
TypeInfoTemplate12("boundarySphereInteraction", ContactForceModel, MotionModel);
boundarySphereInteraction(
const boundaryBase& boundary,
const sphereParticles& sphPrtcls,
const GeometryMotionModel& geomMotion);
create_vCtor
(
BoundarySphereInteractionType,
boundaryBase,
(
const boundaryBase& boundary,
const sphereParticles& sphPrtcls,
const GeometryMotionModel& geomMotion
),
(boundary, sphPrtcls, geomMotion)
);
add_vCtor
(
BoundarySphereInteractionType,
BoundarySphereInteractionType,
boundaryBase
);
~boundarySphereInteraction()override=default;
const auto& sphParticles()const
{
return sphParticles_;
}
const auto& geometryMotion()const
{
return geometryMotion_;
}
ContactListType& ppPairs()
{
return ppPairs_();
}
const ContactListType& ppPairs()const
{
return ppPairs_();
}
ContactListType& pwPairs()
{
return pwPairs_();
}
const ContactListType& pwPairs()const
{
return pwPairs_();
}
virtual
bool sphereSphereInteraction(
real dt,
const ContactForceModel& cfModel)
{
// for default boundary, no thing to be done
return true;
}
bool hearChanges
(
real t,
real dt,
uint32 iter,
const message& msg,
const anyList& varList
) override
{
notImplementedFunction;
return true;
}
void fill(const std::any& val)override
{
notImplementedFunction;
}
static
uniquePtr<BoundarySphereInteractionType> create(
const boundaryBase& boundary,
const sphereParticles& sphPrtcls,
const GeometryMotionModel& geomMotion
);
};
}
#include "boundarySphereInteraction.cpp"
#endif //__boundarySphereInteraction_hpp__

23
src/Interaction/sphereInteraction/boundaries/boundarySphereInteractionList.cpp Normal file
View File

@ -0,0 +1,23 @@
template <typename CFModel, typename gMModel>
pFlow::boundarySphereInteractionList<CFModel, gMModel>::boundarySphereInteractionList
(
const sphereParticles &sphPrtcls,
const gMModel &geomMotion
)
:
ListPtr<boundarySphereInteraction<CFModel,gMModel>>(6),
boundaries_(sphPrtcls.pStruct().boundaries())
{
for(uint32 i=0; i<6; i++)
{
this->set(
i,
boundarySphereInteraction<CFModel, gMModel>::create(
boundaries_[i],
sphPrtcls,
geomMotion));
}
}

40
src/Interaction/sphereInteraction/boundaries/boundarySphereInteractionList.hpp Normal file
View File

@ -0,0 +1,40 @@
#ifndef __boundarySphereInteractionList_hpp__
#define __boundarySphereInteractionList_hpp__
#include "boundaryList.hpp"
#include "ListPtr.hpp"
#include "boundarySphereInteraction.hpp"
namespace pFlow
{
template<typename contactForceModel,typename geometryMotionModel>
class boundarySphereInteractionList
:
public ListPtr<boundarySphereInteraction<contactForceModel,geometryMotionModel>>
{
private:
const boundaryList& boundaries_;
public:
boundarySphereInteractionList(
const sphereParticles& sphPrtcls,
const geometryMotionModel& geomMotion
);
~boundarySphereInteractionList()=default;
};
}
#include "boundarySphereInteractionList.cpp"
#endif //__boundarySphereInteractionList_hpp__

31
src/Interaction/sphereInteraction/boundaries/createBoundarySphereInteraction.hpp Normal file
View File

@ -0,0 +1,31 @@
/*------------------------------- 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 "boundarySphereInteraction.hpp"
#include "periodicBoundarySphereInteraction.hpp"
#define createBoundarySphereInteraction(ForceModel,GeomModel) \
template class pFlow::boundarySphereInteraction< \
ForceModel, \
GeomModel>; \
\
template class pFlow::periodicBoundarySphereInteraction< \
ForceModel, \
GeomModel>;

125
src/Interaction/sphereInteraction/boundaries/periodicBoundarySphereInteraction/periodicBoundarySIKernels.hpp Normal file
View File

@ -0,0 +1,125 @@
namespace pFlow::periodicBoundarySIKernels
{
template<typename ContactListType, typename ContactForceModel>
void sphereSphereInteraction
(
real dt,
const ContactListType& cntctList,
const ContactForceModel& forceModel,
const realx3& transferVec,
const deviceScatteredFieldAccess<realx3>& thisPoints,
const deviceScatteredFieldAccess<realx3>& mirrorPoints,
const deviceViewType1D<real>& diam,
const deviceViewType1D<uint32>& propId,
const deviceViewType1D<realx3>& vel,
const deviceViewType1D<realx3>& rVel,
const deviceViewType1D<realx3>& cForce,
const deviceViewType1D<realx3>& cTorque
)
{
using ValueType = typename ContactListType::ValueType;
uint32 ss = cntctList.size();
uint32 lastItem = cntctList.loopCount();
if(lastItem == 0u)return;
Kokkos::parallel_for(
"pFlow::periodicBoundarySIKernels::sphereSphereInteraction",
deviceRPolicyDynamic(0,lastItem),
LAMBDA_HD(uint32 n)
{
if(!cntctList.isValid(n))return;
auto [i,j] = cntctList.getPair(n);
uint32 ind_i = thisPoints.index(i);
uint32 ind_j = mirrorPoints.index(j);
real Ri = 0.5*diam[ind_i];
real Rj = 0.5*diam[ind_j];
realx3 xi = thisPoints.field()[ind_i];
realx3 xj = mirrorPoints.field()[ind_j]+transferVec;
real dist = length(xj-xi);
real ovrlp = (Ri+Rj) - dist;
if( ovrlp >0.0 )
{
/*auto Vi = thisVel[i];
auto Vj = mirrorVel[j];
auto wi = thisRVel[i];
auto wj = mirrorRVel[j];
auto Nij = (xj-xi)/dist;
auto Vr = Vi - Vj + cross((Ri*wi+Rj*wj), Nij);*/
auto Nij = (xj-xi)/dist;
auto wi = rVel[ind_i];
auto wj = rVel[ind_j];
auto Vr = vel[ind_i] - vel[ind_j] + cross((Ri*wi+Rj*wj), Nij);
auto history = cntctList.getValue(n);
int32 propId_i = propId[ind_i];
int32 propId_j = propId[ind_j];
realx3 FCn, FCt, Mri, Mrj, Mij, Mji;
// calculates contact force
forceModel.contactForce(
dt, i, j,
propId_i, propId_j,
Ri, Rj,
ovrlp,
Vr, Nij,
history,
FCn, FCt);
forceModel.rollingFriction(
dt, i, j,
propId_i, propId_j,
Ri, Rj,
wi, wj,
Nij,
FCn,
Mri, Mrj);
auto M = cross(Nij,FCt);
Mij = Ri*M+Mri;
Mji = Rj*M+Mrj;
auto FC = FCn + FCt;
Kokkos::atomic_add(&cForce[ind_i].x_,FC.x_);
Kokkos::atomic_add(&cForce[ind_i].y_,FC.y_);
Kokkos::atomic_add(&cForce[ind_i].z_,FC.z_);
Kokkos::atomic_add(&cForce[ind_j].x_,-FC.x_);
Kokkos::atomic_add(&cForce[ind_j].y_,-FC.y_);
Kokkos::atomic_add(&cForce[ind_j].z_,-FC.z_);
Kokkos::atomic_add(&cTorque[ind_i].x_, Mij.x_);
Kokkos::atomic_add(&cTorque[ind_i].y_, Mij.y_);
Kokkos::atomic_add(&cTorque[ind_i].z_, Mij.z_);
Kokkos::atomic_add(&cTorque[ind_j].x_, Mji.x_);
Kokkos::atomic_add(&cTorque[ind_j].y_, Mji.y_);
Kokkos::atomic_add(&cTorque[ind_j].z_, Mji.z_);
cntctList.setValue(n,history);
}
else
{
cntctList.setValue(n, ValueType());
}
});
Kokkos::fence();
}
} //pFlow::periodicBoundarySIKernels

65
src/Interaction/sphereInteraction/boundaries/periodicBoundarySphereInteraction/periodicBoundarySphereInteraction.cpp Normal file
View File

@ -0,0 +1,65 @@
/*------------------------------- 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 "periodicBoundarySIKernels.hpp"
template <typename cFM, typename gMM>
pFlow::periodicBoundarySphereInteraction<cFM, gMM>::periodicBoundarySphereInteraction(
const boundaryBase &boundary,
const sphereParticles &sphPrtcls,
const GeometryMotionModel &geomMotion)
: boundarySphereInteraction<cFM,gMM>(boundary, sphPrtcls, geomMotion),
transferVec_(boundary.mirrorBoundary().displacementVectroToMirror())
{
if(boundary.thisBoundaryIndex()%2==1)
{
masterInteraction_ = true;
}
else
{
masterInteraction_ = false;
}
}
template <typename cFM, typename gMM>
bool pFlow::periodicBoundarySphereInteraction<cFM, gMM>::sphereSphereInteraction
(
real dt,
const ContactForceModel &cfModel
)
{
if(!masterInteraction_) return true;
pFlow::periodicBoundarySIKernels::sphereSphereInteraction(
dt,
this->ppPairs(),
cfModel,
transferVec_,
this->boundary().thisPoints(),
this->mirrorBoundary().thisPoints(),
this->sphParticles().diameter().deviceViewAll(),
this->sphParticles().propertyId().deviceViewAll(),
this->sphParticles().velocity().deviceViewAll(),
this->sphParticles().rVelocity().deviceViewAll(),
this->sphParticles().contactForce().deviceViewAll(),
this->sphParticles().contactTorque().deviceViewAll());
return true;
}

95
src/Interaction/sphereInteraction/boundaries/periodicBoundarySphereInteraction/periodicBoundarySphereInteraction.hpp Normal file
View File

@ -0,0 +1,95 @@
/*------------------------------- 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 __periodicBoundarySphereInteraction_hpp__
#define __periodicBoundarySphereInteraction_hpp__
#include "boundarySphereInteraction.hpp"
namespace pFlow
{
template<typename contactForceModel,typename geometryMotionModel>
class periodicBoundarySphereInteraction
:
public boundarySphereInteraction<contactForceModel, geometryMotionModel>
{
public:
using PBSInteractionType =
periodicBoundarySphereInteraction<contactForceModel,geometryMotionModel>;
using BSInteractionType =
boundarySphereInteraction<contactForceModel, geometryMotionModel>;
using GeometryMotionModel = typename BSInteractionType::GeometryMotionModel;
using ContactForceModel = typename BSInteractionType::ContactForceModel;
using MotionModel = typename geometryMotionModel::MotionModel;
using ModelStorage = typename ContactForceModel::contactForceStorage;
using IdType = typename BSInteractionType::IdType;
using IndexType = typename BSInteractionType::IndexType;
using ContactListType = typename BSInteractionType::ContactListType;
private:
realx3 transferVec_;
bool masterInteraction_;
public:
TypeInfoTemplate22("boundarySphereInteraction", "periodic",ContactForceModel, MotionModel);
periodicBoundarySphereInteraction(
const boundaryBase& boundary,
const sphereParticles& sphPrtcls,
const GeometryMotionModel& geomMotion
);
add_vCtor
(
BSInteractionType,
PBSInteractionType,
boundaryBase
);
~periodicBoundarySphereInteraction()override = default;
bool sphereSphereInteraction(
real dt,
const ContactForceModel& cfModel)override;
};
}
#include "periodicBoundarySphereInteraction.cpp"
#endif //__periodicBoundarySphereInteraction_hpp__

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

@ -1,130 +0,0 @@
/*------------------------------- phasicFlow ---------------------------------
O C enter of
O O E ngineering and
O O M ultiscale modeling of
OOOOOOO F luid flow
------------------------------------------------------------------------------
Copyright (C): www.cemf.ir
email: hamid.r.norouzi AT gmail.com
------------------------------------------------------------------------------
Licence:
This file is part of phasicFlow code. It is a free software for simulating
granular and multiphase flows. You can redistribute it and/or modify it under
the terms of GNU General Public License v3 or any other later versions.
phasicFlow is distributed to help others in their research in the field of
granular and multiphase flows, but WITHOUT ANY WARRANTY; without even the
implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
-----------------------------------------------------------------------------*/
template<
typename contactForceModel,
typename geometryMotionModel,
template <class, class, class> class contactListType >
bool pFlow::sphereInteraction<contactForceModel,geometryMotionModel, contactListType>::
createSphereInteraction()
{
realVector_D rhoD(this->densities());
auto modelDict = this->fileDict().subDict("model");
REPORT(1)<<"Createing contact force model . . ."<<endREPORT;
forceModel_ = makeUnique<ContactForceModel>(
this->numMaterials(),
rhoD.deviceVector(),
modelDict );
uint32 nPrtcl = sphParticles_.size();
ppContactList_ = makeUnique<ContactListType>(nPrtcl+1);
pwContactList_ = makeUnique<ContactListType>(nPrtcl/4+1);
return true;
}
template<
typename contactForceModel,
typename geometryMotionModel,
template <class, class, class> class contactListType >
bool pFlow::sphereInteraction<contactForceModel,geometryMotionModel, contactListType>::
sphereSphereInteraction()
{
auto lastItem = ppContactList_().loopCount();
// create the kernel functor
pFlow::sphereInteractionKernels::ppInteractionFunctor
ppInteraction(
this->dt(),
this->forceModel_(),
ppContactList_(), // to be read
sphParticles_.diameter().deviceVectorAll(),
sphParticles_.propertyId().deviceVectorAll(),
sphParticles_.pointPosition().deviceVectorAll(),
sphParticles_.velocity().deviceVectorAll(),
sphParticles_.rVelocity().deviceVectorAll(),
sphParticles_.contactForce().deviceVectorAll(),
sphParticles_.contactTorque().deviceVectorAll()
);
Kokkos::parallel_for(
"",
rpPPInteraction(0,lastItem),
ppInteraction
);
Kokkos::fence();
return true;
}
template<
typename contactForceModel,
typename geometryMotionModel,
template <class, class, class> class contactListType >
bool pFlow::sphereInteraction<contactForceModel,geometryMotionModel, contactListType>::
sphereWallInteraction()
{
int32 lastItem = pwContactList_().loopCount();
real t = this->currentTime();
pFlow::sphereInteractionKernels::pwInteractionFunctor
pwInteraction(
this->dt(),
this->forceModel_(),
pwContactList_(),
geometryMotion_.getTriangleAccessor(),
geometryMotion_.getModel(t) ,
sphParticles_.diameter().deviceVectorAll() ,
sphParticles_.propertyId().deviceVectorAll(),
sphParticles_.pointPosition().deviceVectorAll(),
sphParticles_.velocity().deviceVectorAll(),
sphParticles_.rVelocity().deviceVectorAll() ,
sphParticles_.contactForce().deviceVectorAll(),
sphParticles_.contactTorque().deviceVectorAll() ,
geometryMotion_.triMotionIndex().deviceVectorAll(),
geometryMotion_.propertyId().deviceVectorAll(),
geometryMotion_.contactForceWall().deviceVectorAll()
);
Kokkos::parallel_for(
"",
rpPWInteraction(0,lastItem),
pwInteraction
);
Kokkos::fence();
return true;
}

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

@ -1,226 +0,0 @@
/*------------------------------- phasicFlow ---------------------------------
O C enter of
O O E ngineering and
O O M ultiscale modeling of
OOOOOOO F luid flow
------------------------------------------------------------------------------
Copyright (C): www.cemf.ir
email: hamid.r.norouzi AT gmail.com
------------------------------------------------------------------------------
Licence:
This file is part of phasicFlow code. It is a free software for simulating
granular and multiphase flows. You can redistribute it and/or modify it under
the terms of GNU General Public License v3 or any other later versions.
phasicFlow is distributed to help others in their research in the field of
granular and multiphase flows, but WITHOUT ANY WARRANTY; without even the
implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
-----------------------------------------------------------------------------*/
#ifndef __sphereInteraction_hpp__
#define __sphereInteraction_hpp__
#include "interaction.hpp"
#include "sphereParticles.hpp"
#include "sphereInteractionKernels.hpp"
namespace pFlow
{
template<
typename contactForceModel,
typename geometryMotionModel,
template <class, class, class> class contactListType>
class sphereInteraction
:
public interaction
{
public:
using GeometryMotionModel = geometryMotionModel;
using ContactForceModel = contactForceModel;
using MotionModel = typename geometryMotionModel::MotionModel;
using ModelStorage = typename ContactForceModel::contactForceStorage;
using IdType = typename interaction::IdType;
using IndexType = typename interaction::IndexType;
using ExecutionSpace = typename interaction::ExecutionSpace;
using ContactListType =
contactListType<ModelStorage, ExecutionSpace, IdType>;
using PairsContainerType= typename contactSearch::PairContainerType;
protected:
/// const reference to geometry
const GeometryMotionModel& geometryMotion_;
/// const reference to particles
const sphereParticles& sphParticles_;
/// contact force model
uniquePtr<ContactForceModel> forceModel_ = nullptr;
/// contact list for particle-particle interactoins (keeps the history)
uniquePtr<ContactListType> ppContactList_ = nullptr;
/// contact list for particle-wall interactions (keeps the history)
uniquePtr<ContactListType> pwContactList_ = nullptr;
/// timer for particle-particle interaction computations
Timer ppInteractionTimer_;
/// timer for particle-wall interaction computations
Timer pwInteractionTimer_;
bool createSphereInteraction();
bool managePPContactLists();
bool managePWContactLists();
/// range policy for p-p interaction execution
using rpPPInteraction =
Kokkos::RangePolicy<Kokkos::IndexType<int32>, Kokkos::Schedule<Kokkos::Dynamic>>;
/// range policy for p-w interaction execution
using rpPWInteraction = rpPPInteraction;
public:
TypeInfoTemplate3("sphereInteraction", ContactForceModel, MotionModel, ContactListType);
// constructor
sphereInteraction(
systemControl& control,
const particles& prtcl,
const geometry& geom)
:
interaction(control, prtcl, geom),
geometryMotion_(dynamic_cast<const GeometryMotionModel&>(geom)),
sphParticles_(dynamic_cast<const sphereParticles&>(prtcl)),
ppInteractionTimer_("sphere-sphere interaction", &this->timers()),
pwInteractionTimer_("sphere-wall interaction", &this->timers())
{
if(!createSphereInteraction())
{
fatalExit;
}
}
add_vCtor
(
interaction,
sphereInteraction,
systemControl
);
bool beforeIteration() override
{
return true;
}
bool iterate() override
{
//Info<<"before contact search"<<endInfo;
////Info<<"interaction iterrate start"<<endInfo;
if(this->contactSearch_)
{
if( this->contactSearch_().ppEnterBroadSearch())
{
//Info<<" before ppEnterBroadSearch"<<endInfo;
ppContactList_().beforeBroadSearch();
//Info<<" after ppEnterBroadSearch"<<endInfo;
}
if(this->contactSearch_().pwEnterBroadSearch())
{
//Info<<" before pwEnterBroadSearch"<<endInfo;
pwContactList_().beforeBroadSearch();
//Info<<" after pwEnterBroadSearch"<<endInfo;
}
//Info<<" before broadSearch"<<endInfo;
if( !contactSearch_().broadSearch(
ppContactList_(),
pwContactList_()) )
{
fatalErrorInFunction<<
"unable to perform broadSearch.\n";
fatalExit;
}
//Info<<" before broadSearch"<<endInfo;
if(this->contactSearch_().ppPerformedBroadSearch())
{
//Info<<" before afterBroadSearch"<<endInfo;
ppContactList_().afterBroadSearch();
//Info<<" after afterBroadSearch"<<endInfo;
}
if(this->contactSearch_().pwPerformedBroadSearch())
{
//Info<<" before pwContactList_().afterBroadSearch()"<<endInfo;
pwContactList_().afterBroadSearch();
//Info<<" after pwContactList_().afterBroadSearch()"<<endInfo;
}
}
//Info<<"after contact search"<<endInfo;
if( sphParticles_.numActive()<=0)return true;
//Info<<"before sphereSphereInteraction "<<endInfo;
ppInteractionTimer_.start();
sphereSphereInteraction();
ppInteractionTimer_.end();
//Info<<"after sphereSphereInteraction "<<endInfo;
//Info<<"before sphereWallInteraction "<<endInfo;
pwInteractionTimer_.start();
sphereWallInteraction();
pwInteractionTimer_.end();
//Info<<"after sphereWallInteraction "<<endInfo;
return true;
}
bool afterIteration() override
{
return true;
}
bool update(const eventMessage& msg)override
{
// it requires not action regarding any changes in the
// point structure
return true;
}
bool sphereSphereInteraction();
bool sphereWallInteraction();
};
}
#include "sphereInteraction.cpp"
#endif //__sphereInteraction_hpp__

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