diff --git a/CMakeLists.txt b/CMakeLists.txt index a3aa0978..3f9af446 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -6,6 +6,7 @@ project(phasicFlow VERSION 0.1 ) set(CMAKE_CXX_STANDARD 17) 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) message(STATUS ${CMAKE_INSTALL_PREFIX}) @@ -13,24 +14,26 @@ mark_as_advanced(FORCE var Kokkos_ENABLE_CUDA_LAMBDA) mark_as_advanced(FORCE var Kokkos_ENABLE_OPENMP) mark_as_advanced(FORCE var Kokkos_ENABLE_SERIAL) mark_as_advanced(FORCE var Kokkos_ENABLE_CUDA_LAMBDA) +mark_as_advanced(FORCE var BUILD_SHARED_LIBS) -option(BUILD_SHARED_LIBS "Build using shared libraries" ON) option(USE_STD_PARALLEL_ALG "Use TTB std parallel algorithms" ON) option(pFlow_Build_Serial "Build phasicFlow and backends for serial execution" ON) option(pFlow_Build_OpenMP "Build phasicFlow and backends for OpenMP execution" OFF) option(pFlow_Build_Cuda "Build phasicFlow and backends for Cuda execution" OFF) -if(dFlow_Build_Serial) +set(BUILD_SHARED_LIBS ON CACHE BOOL "Build using shared libraries" FORCE) + +if(pFlow_Build_Serial) set(Kokkos_ENABLE_SERIAL ON CACHE BOOL "Serial execution" FORCE) set(Kokkos_ENABLE_OPENMP OFF CACHE BOOL "OpenMP execution" FORCE) set(Kokkos_ENABLE_CUDA OFF CACHE BOOL "Cuda execution" FORCE) set(Kokkos_ENABLE_CUDA_LAMBDA OFF CACHE BOOL "Cuda execution" FORCE) -elseif(dFlow_Build_OpenMP ) +elseif(pFlow_Build_OpenMP ) set(Kokkos_ENABLE_SERIAL ON CACHE BOOL "Serial execution" FORCE) set(Kokkos_ENABLE_OPENMP ON CACHE BOOL "OpenMP execution" FORCE) set(Kokkos_ENABLE_CUDA OFF CACHE BOOL "Cuda execution" FORCE) set(Kokkos_ENABLE_CUDA_LAMBDA OFF CACHE BOOL "Cuda execution" FORCE) -elseif(dFlow_Build_Cuda) +elseif(pFlow_Build_Cuda) set(Kokkos_ENABLE_SERIAL ON CACHE BOOL "Serial execution" FORCE) set(Kokkos_ENABLE_OPENMP OFF CACHE BOOL "OpenMP execution" FORCE) set(Kokkos_ENABLE_CUDA ON CACHE BOOL "Cuda execution" FORCE) diff --git a/cmake/bashrc b/cmake/bashrc index 8e923ee9..50f8fa60 100644 --- a/cmake/bashrc +++ b/cmake/bashrc @@ -1,7 +1,7 @@ export pFlow_PROJECT_VERSION=v0.1 -export pFlow_PROJECT=phasicFlow-$pFlow_PROJECT_VERSION +export pFlow_PROJECT=phasicFlow projectDir="$HOME/PhasicFlow" diff --git a/cmake/makeExecutableGlobals.cmake b/cmake/makeExecutableGlobals.cmake index 5aef5261..24eb6e0f 100644 --- a/cmake/makeExecutableGlobals.cmake +++ b/cmake/makeExecutableGlobals.cmake @@ -1,4 +1,4 @@ -#add a library to dFlow with source files and target_link_libs (thouse under the main CMakeLists.txt) +#add a library to pFlow with source files and target_link_libs (thouse under the main CMakeLists.txt) macro(pFlow_make_executable_install target_name source_files target_link_libs) # add library diff --git a/cmake/makeLibraryGlobals.cmake b/cmake/makeLibraryGlobals.cmake index 26d567c7..b5e30156 100644 --- a/cmake/makeLibraryGlobals.cmake +++ b/cmake/makeLibraryGlobals.cmake @@ -1,4 +1,4 @@ -#add a library to dFlow with source files and target_link_libs (thouse under the main CMakeLists.txt) +#add a library to pFlow with source files and target_link_libs (thouse under the main CMakeLists.txt) macro(pFlow_add_library_install target_name src_files target_link_libs) set(source_files ${${src_files}}) diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index ff2582b5..b6a61171 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -1,15 +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(Interaction) -#add_subdirectory(MotionModel) +add_subdirectory(MotionModel) -#add_subdirectory(Geometry) +add_subdirectory(Geometry) diff --git a/src/Geometry/CMakeLists.txt b/src/Geometry/CMakeLists.txt new file mode 100644 index 00000000..e077dc46 --- /dev/null +++ b/src/Geometry/CMakeLists.txt @@ -0,0 +1,12 @@ + +list(APPEND SourceFiles +geometry/geometry.C +geometryMotion/geometryMotions.C +) + +set(link_libs Kokkos::kokkos phasicFlow MotionModel Property) + +pFlow_add_library_install(Geometry SourceFiles link_libs) + + + diff --git a/src/Geometry/geometry/demGeometry.H b/src/Geometry/geometry/demGeometry.H new file mode 100644 index 00000000..e44f4dd2 --- /dev/null +++ b/src/Geometry/geometry/demGeometry.H @@ -0,0 +1,47 @@ +/*------------------------------- 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 __demGeometry_H__ +#define __demGeometry_H__ + +#include "demComponent.H" + +namespace pFlow +{ + + +class demGeometry +: + public demComponent +{ + +public: + demGeometry(systemControl& control) + : + demComponent("geometry", control) + {} + + +}; + +} + +#endif diff --git a/src/Geometry/geometry/geometry.C b/src/Geometry/geometry/geometry.C new file mode 100644 index 00000000..723acd0f --- /dev/null +++ b/src/Geometry/geometry/geometry.C @@ -0,0 +1,316 @@ +/*------------------------------- 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 "geometry.H" +#include "vocabs.H" + + +bool pFlow::geometry::findPropertyId() +{ + + int8Vector propId(0, surface().capacity(),RESERVE()); + propId.clear(); + + uint32 pId; + forAll(matI, materialName_) + { + + if( !wallProperty_.nameToIndex( materialName_[matI], pId ) ) + { + fatalErrorInFunction<< + "material name for the geometry is invalid: "<< materialName_[matI]<( + objectFile( + "triSurface", + "", + objectFile::READ_ALWAYS, + objectFile::WRITE_ALWAYS + ) + ) + ), + motionComponentName_( + control.geometry().emplaceObject( + objectFile( + "motionComponentName", + "", + objectFile::READ_ALWAYS, + objectFile::WRITE_ALWAYS + ), + "motionNamesList" + ) + ), + materialName_( + control.geometry().emplaceObject( + objectFile( + "materialName", + "", + objectFile::READ_ALWAYS, + objectFile::WRITE_ALWAYS + ), + "materialNamesList" + ) + ), + propertyId_( + control.geometry().emplaceObject( + objectFile( + "propertyId", + "", + objectFile::READ_NEVER, + objectFile::WRITE_NEVER), + surface(), + 0 ) ), + contactForceWall_( + control.geometry().emplaceObject( + objectFile( + "contactForceWall", + "", + objectFile::READ_IF_PRESENT, + objectFile::WRITE_ALWAYS), + surface(), + zero3) ), + stressWall_( + control.geometry().emplaceObject( + objectFile( + "stressWall", + "", + objectFile::READ_IF_PRESENT, + objectFile::WRITE_ALWAYS), + surface(), + zero3) ) +{ + + if(!findPropertyId()) + { + fatalExit; + } +} + +pFlow::geometry::geometry +( + systemControl& control, + const property& prop, + const multiTriSurface& triSurface, + const wordVector& motionCompName, + const wordVector& matName +) +: + demGeometry(control), + wallProperty_(prop), + geometryRepository_(control.geometry()), + triSurface_( + control.geometry().emplaceObject( + objectFile( + "triSurface", + "", + objectFile::READ_NEVER, + objectFile::WRITE_ALWAYS + ), + triSurface + ) + ), + motionComponentName_( + control.geometry().emplaceObject( + objectFile( + "motionComponentName", + "", + objectFile::READ_NEVER, + objectFile::WRITE_ALWAYS + ), + "motionNamesList", + motionCompName + ) + ), + materialName_( + control.geometry().emplaceObject( + objectFile( + "materialName", + "", + objectFile::READ_NEVER, + objectFile::WRITE_ALWAYS + ), + "materialNamesList", + matName + ) + ), + propertyId_( + control.geometry().emplaceObject( + objectFile( + "propertyId", + "", + objectFile::READ_NEVER, + objectFile::WRITE_NEVER), + surface(), + 0 ) ), + contactForceWall_( + control.geometry().emplaceObject( + objectFile( + "contactForceWall", + "", + objectFile::READ_IF_PRESENT, + objectFile::WRITE_ALWAYS), + surface(), + zero3) ), + stressWall_( + control.geometry().emplaceObject( + objectFile( + "stressWall", + "", + objectFile::READ_IF_PRESENT, + objectFile::WRITE_ALWAYS), + surface(), + zero3) ) +{ + if(!findPropertyId()) + { + fatalExit; + } +} + +pFlow::geometry::geometry +( + systemControl& control, + const property& prop, + const dictionary& dict, + const multiTriSurface& triSurface, + const wordVector& motionCompName, + const wordVector& matName +) +: + geometry(control, prop, triSurface, motionCompName, matName) +{} + + +pFlow::uniquePtr + pFlow::geometry::create +( + systemControl& control, + const property& prop +) +{ + //motionModelFile__ + auto motionDictPtr = IOobject::make + ( + objectFile + ( + motionModelFile__, + control.geometry().path(), + objectFile::READ_ALWAYS, + objectFile::WRITE_NEVER + ), + motionModelFile__, + true + ); + + word model = motionDictPtr().getObject().getVal("motionModel"); + + auto geomModel = angleBracketsNames("geometry", model); + + Report(1)<< "Selecting geometry model . . ."< + pFlow::geometry::create( + systemControl& control, + const property& prop, + const dictionary& dict, + const multiTriSurface& triSurface, + const wordVector& motionCompName, + const wordVector& propName) +{ + + word model = dict.getVal("motionModel"); + + auto geomModel = angleBracketsNames("geometry", model); + + Report(1)<< "Selecting geometry model . . ."<zeroForce(); + return true; + + } + + bool afterIteration() override { + + auto Force = contactForceWall_.deviceVectorAll(); + auto area = triSurface_.area().deviceVectorAll(); + auto stress = stressWall_.deviceVectorAll(); + auto numTri =triSurface_.size(); + + Kokkos::parallel_for( + "geometry::calculateStress", + numTri, + LAMBDA_HD(int32 i){ + stress[i] = Force[i]/area[i]; + }); + Kokkos::fence(); + return true; + + } + + bool write()const + { + return owner().write(); + } + + + // static + + static + uniquePtr create(systemControl& control, const property& prop); + + static + uniquePtr create( + systemControl& control, + const property& prop, + const dictionary& dict, + const multiTriSurface& triSurface, + const wordVector& motionCompName, + const wordVector& propName); + +}; + +} + +#endif diff --git a/src/Geometry/geometryMotion/geometryMotion.C b/src/Geometry/geometryMotion/geometryMotion.C new file mode 100644 index 00000000..5e6136f5 --- /dev/null +++ b/src/Geometry/geometryMotion/geometryMotion.C @@ -0,0 +1,189 @@ +/*------------------------------- 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 +bool pFlow::geometryMotion::moveGeometry() +{ + + real dt = this->dt(); + + auto pointMIndex= pointMotionIndex_.deviceVector(); + auto mModel = motionModel_.getModel(); + realx3* points = triSurface_.pointsData_D(); + auto numPoints = triSurface_.numPoints(); + + + Kokkos::parallel_for( + "geometryMotion::movePoints", + numPoints, + LAMBDA_HD(int32 i){ + auto newPos = mModel.transferPoint(pointMIndex[i], points[i], dt); + points[i] = newPos; + }); + + Kokkos::fence(); + + // end of calculations + moveGeomTimer_.end(); + + return true; +} + +template +bool pFlow::geometryMotion::findMotionIndex() +{ + motionIndex_.clear(); + triMotionIndex_.reserve( this->surface().capacity() ); + triMotionIndex_.clear(); + + forAll( surfI, motionComponentName_) + { + 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 +pFlow::geometryMotion::geometryMotion +( + systemControl& control, + const property& prop +) +: + geometry(control, prop), + motionModel_( + this->owner().template emplaceObject( + objectFile( + motionModelFile__, + "", + objectFile::READ_ALWAYS, + objectFile::WRITE_ALWAYS + ) + ) + ), + moveGeomTimer_("move geometry", &this->timers()) +{ + findMotionIndex(); +} + +template +pFlow::geometryMotion::geometryMotion +( + systemControl& control, + const property& prop, + const multiTriSurface& triSurface, + const wordVector& motionCompName, + const wordVector& propName, + const MotionModel& motionModel +) +: + geometry( + control, + prop, + triSurface, + motionCompName, + propName + ), + motionModel_( + this->owner().template emplaceObject( + objectFile( + motionModelFile__, + "", + objectFile::READ_NEVER, + objectFile::WRITE_ALWAYS + ), + motionModel + ) + ), + moveGeomTimer_("move geometry", &this->timers()) +{ + findMotionIndex(); +} + +template +pFlow::geometryMotion::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( + objectFile( + motionModelFile__, + "", + objectFile::READ_NEVER, + objectFile::WRITE_ALWAYS + ), + dict + ) + ), + moveGeomTimer_("move geometry", &this->timers()) +{ + findMotionIndex(); +} + +template +bool pFlow::geometryMotion::iterate() +{ + if( motionModel_.isMoving() ) + { + moveGeomTimer_.start(); + moveGeometry(); + moveGeomTimer_.end(); + } + return true; +} \ No newline at end of file diff --git a/src/Geometry/geometryMotion/geometryMotion.H b/src/Geometry/geometryMotion/geometryMotion.H new file mode 100644 index 00000000..cdfad213 --- /dev/null +++ b/src/Geometry/geometryMotion/geometryMotion.H @@ -0,0 +1,155 @@ +/*------------------------------- 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 __geometryMotion_H__ +#define __geometryMotion_H__ + + +#include "geometry.H" +#include "VectorDuals.H" + +namespace pFlow +{ + + +template +class geometryMotion +: + public geometry +{ +public: + + using MotionModel = MotionModelType; + +protected: + + + MotionModel& motionModel_; + + // motion indext mapped on each surface + int32Vector_HD motionIndex_; + + // motion index mapped on each triangle + int8Vector_HD triMotionIndex_; + + /// motion index mapped on each point + int8Vector_HD pointMotionIndex_; + + // timer for moveGeometry + Timer moveGeomTimer_; + + + bool findMotionIndex(); + +public: + + // type info + TypeNameTemplate("geometry", MotionModel); + + //// - Constructors + + geometryMotion(systemControl& control, const property& prop); + + // construct from components + geometryMotion( + systemControl& control, + const property& prop, + const multiTriSurface& triSurface, + 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); + + + + add_vCtor + ( + geometry, + geometryMotion, + systemControl + ); + + add_vCtor + ( + geometry, + geometryMotion, + dictionary + ); + + //// - Methods + + auto getModel()const + { + return motionModel_.getModel(); + } + + word motionModelTypeName()const override + { + return motionModel_.typeName(); + } + + const int8Vector_HD& triMotionIndex()const override + { + return triMotionIndex_; + } + + const int8Vector_HD& pointMotionIndex()const override + { + return pointMotionIndex_; + } + + // - iterate + bool beforeIteration() override { + geometry::beforeIteration(); + return true; + } + + bool iterate() override ; + + bool afterIteration() override { + geometry::afterIteration(); + return true; + } + + + bool moveGeometry(); + +}; + + +} + +#include "geometryMotion.C" + +#ifndef BUILD_SHARED_LIBS + #include "geometryMotionsInstantiate.C" +#endif + + +#endif //__geometryMotion_H__ diff --git a/src/Geometry/geometryMotion/geometryMotions.C b/src/Geometry/geometryMotion/geometryMotions.C new file mode 100644 index 00000000..c79e8565 --- /dev/null +++ b/src/Geometry/geometryMotion/geometryMotions.C @@ -0,0 +1,25 @@ +/*------------------------------- 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 "geometryMotions.H" + +#ifdef BUILD_SHARED_LIBS +#include "geometryMotionsInstantiate.C" +#endif diff --git a/src/Geometry/geometryMotion/geometryMotions.H b/src/Geometry/geometryMotion/geometryMotions.H new file mode 100644 index 00000000..d7e6d704 --- /dev/null +++ b/src/Geometry/geometryMotion/geometryMotions.H @@ -0,0 +1,40 @@ +/*------------------------------- 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 __geometryMotions_H__ +#define __geometryMotions_H__ + +#include "geometryMotion.H" +#include "fixedWall.H" +#include "rotatingAxisMotion.H" + + +namespace pFlow +{ + +typedef geometryMotion rotationAxisMotionGeometry; + +typedef geometryMotion fixedGeometry; + + +} + + +#endif \ No newline at end of file diff --git a/src/Geometry/geometryMotion/geometryMotionsInstantiate.C b/src/Geometry/geometryMotion/geometryMotionsInstantiate.C new file mode 100644 index 00000000..2e700db0 --- /dev/null +++ b/src/Geometry/geometryMotion/geometryMotionsInstantiate.C @@ -0,0 +1,28 @@ +/*------------------------------- 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.H" +#include "rotatingAxisMotion.H" + +template class pFlow::geometryMotion; + +template class pFlow::geometryMotion; + + diff --git a/src/Integration/AdamsBashforth2/AdamsBashforth2.C b/src/Integration/AdamsBashforth2/AdamsBashforth2.C new file mode 100644 index 00000000..baf76a54 --- /dev/null +++ b/src/Integration/AdamsBashforth2/AdamsBashforth2.C @@ -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. + +-----------------------------------------------------------------------------*/ + +#include "AdamsBashforth2.H" + +//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 +) +: + integration(baseName, owner, pStruct, method), + dy1_( + owner.emplaceObject( + objectFile( + groupNames(baseName,"dy1"), + "", + objectFile::READ_IF_PRESENT, + objectFile::WRITE_ALWAYS), + pStruct, + zero3)) +{ + +} + +bool pFlow::AdamsBashforth2::predict +( + real UNUSED(dt), + realx3Vector_D& UNUSED(y), + realx3Vector_D& UNUSED(dy) +) +{ + + return true; +} + +bool pFlow::AdamsBashforth2::correct +( + real dt, + realx3Vector_D& y, + realx3Vector_D& dy +) +{ + if(this->pStruct().allActive()) + { + return intAll(dt, y, dy, this->pStruct().activeRange()); + } + else + { + return intRange(dt, y, dy, this->pStruct().activePointsMaskD()); + } + + 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(3.0 / 2.0) * d_dy[i] - static_cast(1.0 / 2.0) * d_dy1[i]); + d_dy1[i] = d_dy[i]; + }); + Kokkos::fence(); + + return true; +} \ No newline at end of file diff --git a/src/Integration/AdamsBashforth2/AdamsBashforth2.H b/src/Integration/AdamsBashforth2/AdamsBashforth2.H new file mode 100644 index 00000000..22dd792c --- /dev/null +++ b/src/Integration/AdamsBashforth2/AdamsBashforth2.H @@ -0,0 +1,109 @@ +/*------------------------------- 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 __AdamsBashforth2_H__ +#define __AdamsBashforth2_H__ + + +#include "integration.H" +#include "pointFields.H" + +namespace pFlow +{ + + +class AdamsBashforth2 +: + public integration +{ +protected: + + realx3PointField_D& dy1_; + + using rpIntegration = Kokkos::RangePolicy< + DefaultExecutionSpace, + Kokkos::Schedule, + Kokkos::IndexType + >; +public: + + // type info + TypeName("AdamsBashforth2"); + + //// - Constructors + AdamsBashforth2( + const word& baseName, + repository& owner, + const pointStructure& pStruct, + const word& method); + + virtual ~AdamsBashforth2()=default; + + // - add a virtual constructor + add_vCtor( + integration, + AdamsBashforth2, + word); + + + //// - Methods + bool predict(real UNUSED(dt), realx3Vector_D& UNUSED(y), realx3Vector_D& UNUSED(dy)) override; + + bool correct(real dt, realx3Vector_D& y, realx3Vector_D& dy) override; + + bool intAll(real dt, realx3Vector_D& y, realx3Vector_D& dy, range activeRng); + + template + bool intRange(real dt, realx3Vector_D& y, realx3Vector_D& dy, activeFunctor activeP ); + +}; + +template +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(3.0 / 2.0) * d_dy[i] - static_cast(1.0 / 2.0) * d_dy1[i]); + d_dy1[i] = d_dy[i]; + } + }); + Kokkos::fence(); + + + return true; +} + +} // pFlow + +#endif //__integration_H__ diff --git a/src/Integration/AdamsBashforth3/AdamsBashforth3.C b/src/Integration/AdamsBashforth3/AdamsBashforth3.C new file mode 100644 index 00000000..c772ad13 --- /dev/null +++ b/src/Integration/AdamsBashforth3/AdamsBashforth3.C @@ -0,0 +1,120 @@ +/*------------------------------- 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 "AdamsBashforth3.H" + +//const real AB3_coef[] = { 23.0 / 12.0, 16.0 / 12.0, 5.0 / 12.0 }; + +pFlow::AdamsBashforth3::AdamsBashforth3 +( + const word& baseName, + repository& owner, + const pointStructure& pStruct, + const word& method +) +: + integration(baseName, owner, pStruct, method), + /*dy1_( + owner.emplaceObject( + objectFile( + groupNames(baseName,"dy1"), + "", + objectFile::READ_IF_PRESENT, + objectFile::WRITE_ALWAYS), + pStruct, + zero3)), + dy2_( + owner.emplaceObject( + objectFile( + groupNames(baseName,"dy2"), + "", + objectFile::READ_IF_PRESENT, + objectFile::WRITE_ALWAYS), + pStruct, + zero3))*/ + history_( + owner.emplaceObject( + objectFile( + groupNames(baseName,"AB3History"), + "", + objectFile::READ_IF_PRESENT, + objectFile::WRITE_ALWAYS), + pStruct, + AB3History({zero3,zero3}))) + +{ + +} + +bool pFlow::AdamsBashforth3::predict +( + real UNUSED(dt), + realx3Vector_D& UNUSED(y), + realx3Vector_D& UNUSED(dy) +) +{ + + return true; +} + +bool pFlow::AdamsBashforth3::correct +( + real dt, + realx3Vector_D& y, + realx3Vector_D& dy +) +{ + + if(this->pStruct().allActive()) + { + return intAll(dt, y, dy, this->pStruct().activeRange()); + } + else + { + return intRange(dt, y, dy, this->pStruct().activePointsMaskD()); + } + + return true; +} + +bool pFlow::AdamsBashforth3::intAll( + real dt, + realx3Vector_D& y, + realx3Vector_D& dy, + range activeRng) +{ + auto d_dy = dy.deviceVectorAll(); + auto d_y = y.deviceVectorAll(); + auto d_history = history_.deviceVectorAll(); + + Kokkos::parallel_for( + "AdamsBashforth3::correct", + rpIntegration (activeRng.first, activeRng.second), + LAMBDA_HD(int32 i){ + auto ldy = d_dy[i]; + d_y[i] += dt*( static_cast(23.0 / 12.0) * ldy + - static_cast(16.0 / 12.0) * d_history[i].dy1_ + + static_cast(5.0 / 12.0) * d_history[i].dy2_); + d_history[i] = {ldy ,d_history[i].dy1_}; + }); + Kokkos::fence(); + + return true; +} \ No newline at end of file diff --git a/src/Integration/AdamsBashforth3/AdamsBashforth3.H b/src/Integration/AdamsBashforth3/AdamsBashforth3.H new file mode 100644 index 00000000..c432926e --- /dev/null +++ b/src/Integration/AdamsBashforth3/AdamsBashforth3.H @@ -0,0 +1,166 @@ +/*------------------------------- 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 __AdamsBashforth3_H__ +#define __AdamsBashforth3_H__ + + +#include "integration.H" +#include "pointFields.H" + +namespace pFlow +{ + +struct AB3History +{ + realx3 dy1_={0,0,0}; + realx3 dy2_={0,0,0}; + + TypeNameNV("AB3History"); +}; + + +INLINE_FUNCTION +iIstream& operator>>(iIstream& str, AB3History& ab3) +{ + str.readBegin("AB3History"); + + str >> ab3.dy1_; + str >> ab3.dy2_; + + str.readEnd("AB3History"); + + str.check(FUNCTION_NAME); + + return str; + +} + +INLINE_FUNCTION +iOstream& operator<<(iOstream& str, const AB3History& ab3) +{ + str << token::BEGIN_LIST << ab3.dy1_ + << token::SPACE << ab3.dy2_ + << token::END_LIST; + + str.check(FUNCTION_NAME); + + return str; +} + + +class AdamsBashforth3 +: + public integration +{ +protected: + + using HistoryFieldType = pointField; + //realx3PointField_D& dy1_; + + //realx3PointField_D& dy2_; + + // this is a device + HistoryFieldType history_; + + using rpIntegration = Kokkos::RangePolicy< + DefaultExecutionSpace, + Kokkos::Schedule, + Kokkos::IndexType + >; + +public: + + // type info + TypeName("AdamsBashforth3"); + + //// - Constructors + AdamsBashforth3( + const word& baseName, + repository& owner, + const pointStructure& pStruct, + const word& method); + + virtual ~AdamsBashforth3()=default; + + // - add a virtual constructor + add_vCtor( + integration, + AdamsBashforth3, + word); + + + //// - Methods + bool predict( + real UNUSED(dt), + realx3Vector_D & UNUSED(y), + realx3Vector_D& UNUSED(dy)) override; + + bool correct(real dt, + realx3Vector_D & y, + realx3Vector_D& dy) override; + + bool intAll(real dt, + realx3Vector_D& y, + realx3Vector_D& dy, + range activeRng); + + template + bool intRange(real dt, + realx3Vector_D& y, + realx3Vector_D& dy, + activeFunctor activeP ); + +}; + + +template +bool pFlow::AdamsBashforth3::intRange( + real dt, + realx3Vector_D& y, + realx3Vector_D& dy, + activeFunctor activeP ) +{ + auto d_dy = dy.deviceVectorAll(); + auto d_y = y.deviceVectorAll(); + auto d_history = history_.deviceVectorAll(); + auto activeRng = activeP.activeRange(); + + Kokkos::parallel_for( + "AdamsBashforth3::correct", + rpIntegration (activeRng.first, activeRng.second), + LAMBDA_HD(int32 i){ + if( activeP(i)) + { + auto ldy = d_dy[i]; + d_y[i] += dt*( static_cast(23.0 / 12.0) * ldy + - static_cast(16.0 / 12.0) * d_history[i].dy1_ + + static_cast(5.0 / 12.0) * d_history[i].dy2_); + d_history[i] = {ldy ,d_history[i].dy1_}; + } + }); + Kokkos::fence(); + + return true; +} + +} // pFlow + +#endif //__integration_H__ diff --git a/src/Integration/CMakeLists.txt b/src/Integration/CMakeLists.txt new file mode 100644 index 00000000..6c74df58 --- /dev/null +++ b/src/Integration/CMakeLists.txt @@ -0,0 +1,13 @@ + +list(APPEND SourceFiles +integration/integration.C +AdamsBashforth3/AdamsBashforth3.C +AdamsBashforth2/AdamsBashforth2.C +) + +set(link_libs Kokkos::kokkos phasicFlow) + +pFlow_add_library_install(Integration SourceFiles link_libs) + + + diff --git a/src/Integration/integration/integration.C b/src/Integration/integration/integration.C new file mode 100644 index 00000000..7a704a0c --- /dev/null +++ b/src/Integration/integration/integration.C @@ -0,0 +1,62 @@ +/*------------------------------- 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 "integration.H" + +pFlow::integration::integration +( + const word& baseName, + repository& owner, + const pointStructure& pStruct, + const word& method +) +: + owner_(owner), + baseName_(baseName), + pStruct_(pStruct) +{ + CONSUME_PARAM(method); +} + + +pFlow::uniquePtr + pFlow::integration::create( + const word& baseName, + repository& owner, + const pointStructure& pStruct, + const word& method) +{ + if( wordvCtorSelector_.search(method) ) + { + return wordvCtorSelector_[method] (baseName, owner, pStruct, method); + } + else + { + printKeys + ( + fatalError << "Ctor Selector "<< method << " dose not exist. \n" + <<"Avaiable ones are: \n\n" + , + wordvCtorSelector_ + ); + fatalExit; + } + return nullptr; +} diff --git a/src/Integration/integration/integration.H b/src/Integration/integration/integration.H new file mode 100644 index 00000000..7050499d --- /dev/null +++ b/src/Integration/integration/integration.H @@ -0,0 +1,103 @@ +/*------------------------------- 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 __integration_H__ +#define __integration_H__ + + +#include "virtualConstructor.H" +#include "Vectors.H" +#include "pointStructure.H" +#include "repository.H" + + +namespace pFlow +{ + + +class integration +{ +protected: + + repository& owner_; + + const word baseName_; + + const pointStructure& pStruct_; + +public: + + // type info + TypeName("integration"); + + //// - Constructors + integration( + const word& baseName, + repository& owner, + const pointStructure& pStruct, + const word& method); + + virtual ~integration()=default; + + // - add a virtual constructor + create_vCtor( + integration, + word, + (const word& baseName, + repository& owner, + const pointStructure& pStruct, + const word& method), + (baseName, owner, pStruct, method) ); + + + //// - Methods + + const auto& pStruct()const + { + return pStruct_; + } + + virtual bool predict(real dt, realx3Vector_D& y, realx3Vector_D& dy) = 0; + + virtual bool correct(real dt, realx3Vector_D& y, realx3Vector_D& dy) = 0; + + const word& baseName()const + { + return baseName_; + } + + repository& owner() + { + return owner_; + } + + static + uniquePtr create( + const word& baseName, + repository& owner, + const pointStructure& pStruct, + const word& method); + +}; + +} // pFlow + + +#endif //__integration_H__ diff --git a/src/Integration/integration/integrations.H b/src/Integration/integration/integrations.H new file mode 100644 index 00000000..ba777359 --- /dev/null +++ b/src/Integration/integration/integrations.H @@ -0,0 +1,28 @@ +/*------------------------------- 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 __integrations_H__ +#define __integrations_H__ + +#include "integration.H" +#include "AdamsBashforth2.H" +#include "AdamsBashforth3.H" + +#endif \ No newline at end of file diff --git a/src/Interaction/CMakeLists.txt b/src/Interaction/CMakeLists.txt new file mode 100644 index 00000000..122f2091 --- /dev/null +++ b/src/Interaction/CMakeLists.txt @@ -0,0 +1,16 @@ + +set(SourceFiles +contactSearch/contactSearch/contactSearch.C +contactSearch/ContactSearch/ContactSearchs.C +interaction/interaction.C +sphereInteraction/sphereInteractions.C +) + +set(link_libs Kokkos::kokkos phasicFlow Property Particles Geometry) + +pFlow_add_library_install(Interaction SourceFiles link_libs) + +#additional locations +target_include_directories(Interaction PUBLIC "./contactSearch") + + diff --git a/src/Interaction/Models/contactForce/linearCF.H b/src/Interaction/Models/contactForce/linearCF.H new file mode 100644 index 00000000..57b1b46a --- /dev/null +++ b/src/Interaction/Models/contactForce/linearCF.H @@ -0,0 +1,284 @@ +/*------------------------------- 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 __linearCF_H__ +#define __linearCF_H__ + +#include "types.H" +#include "symArrays.H" + +namespace pFlow::cfModels +{ + +template +class linear +{ +public: + + struct contactForceStorage + { + realx3 overlap_t_ = 0.0; + }; + + struct linearProperties + { + real kn_ = 1000.0; + real kt_ = 800.0; + real ethan_ = 0.0; + real ethat_ = 0.0; + real mu_ = 0.00001; + + INLINE_FUNCTION_HD + linearProperties(){} + + INLINE_FUNCTION_HD + linearProperties(real kn, real kt, real etha_n, real etha_t, real mu ): + kn_(kn), kt_(kt), ethan_(etha_n),ethat_(etha_t), mu_(mu) + {} + + INLINE_FUNCTION_HD + linearProperties(const linearProperties&)=default; + + INLINE_FUNCTION_HD + linearProperties& operator=(const linearProperties&)=default; + + INLINE_FUNCTION_HD + ~linearProperties() = default; + }; + +protected: + + using LinearArrayType = symArray; + + int32 numMaterial_ = 0; + + ViewType1D rho_; + + LinearArrayType linearProperties_; + + /**/ + + bool readLinearDictionary(const dictionary& dict) + { + auto kn = dict.getVal("kn"); + auto kt = dict.getVal("kt"); + auto en = dict.getVal("en"); + auto et = dict.getVal("et"); + auto mu = dict.getVal("mu"); + + auto nElem = kn.size(); + + + if(nElem != kt.size()) + { + fatalErrorInFunction<< + "sizes of kn("< prop(nElem); + forAll(i,kn) + { + prop[i] = {kn[i], kt[i], etha_n[i], etha_t[i], mu[i]}; + } + + linearProperties_.assign(prop); + + return true; + + } + + static const char* modelName() + { + if constexpr (limited) + { + return "linearLimited"; + } + else + { + return "linearNonLimited"; + } + return ""; + } + +public: + + + TypeNameNV(modelName()); + + INLINE_FUNCTION_HD + linear(){} + + linear(int32 nMaterial, const ViewType1D& rho, const dictionary& dict) + : + numMaterial_(nMaterial), + rho_("rho",nMaterial), + linearProperties_("linearProperties",nMaterial) + { + + Kokkos::deep_copy(rho_,rho); + if(!readLinearDictionary(dict)) + { + fatalExit; + } + } + + INLINE_FUNCTION_HD + linear(const linear&) = default; + + INLINE_FUNCTION_HD + linear(linear&&) = default; + + INLINE_FUNCTION_HD + linear& operator=(const linear&) = default; + + INLINE_FUNCTION_HD + linear& operator=(linear&&) = default; + + + INLINE_FUNCTION_HD + ~linear()=default; + + INLINE_FUNCTION_HD + int32 numMaterial()const + { + return numMaterial_; + } + + //// - Methods + + INLINE_FUNCTION_HD + void contactForce + ( + const real dt, + const int32 i, + const int32 j, + const int32 propId_i, + const int32 propId_j, + const real Ri, + const real Rj, + const real ovrlp_n, + const realx3& Vr, + const realx3& Nij, + contactForceStorage& history, + realx3& FCn, + realx3& FCt + )const + { + + auto prop = linearProperties_(propId_i,propId_j); + + + real vrn = dot(Vr, Nij); + realx3 Vt = Vr - vrn*Nij; + + history.overlap_t_ += Vt*dt; + + real mi = 3*Pi/4*pow(Ri,3.0)*rho_[propId_i]; + real mj = 3*Pi/4*pow(Rj,3.0)*rho_[propId_j]; + + real sqrt_meff = sqrt((mi*mj)/(mi+mj)); + + FCn = (-prop.kn_ * ovrlp_n - sqrt_meff * prop.ethan_ * vrn)*Nij; + FCt = -prop.kt_ * history.overlap_t_ - sqrt_meff * prop.ethat_*Vt; + + real ft = length(FCt); + real ft_fric = prop.mu_ * length(FCn); + + if(ft > ft_fric) + { + if( length(history.overlap_t_) >static_cast(0.0)) + { + if constexpr (limited) + { + FCt *= (ft_fric/ft); + history.overlap_t_ = - (FCt/prop.kt_); + } + else + { + FCt = (FCt/ft)*ft_fric; + } + //cout<<"friction is applied here \n"; + + } + else + { + FCt = 0.0; + } + } + + } + +}; + +} //pFlow::cfModels + +#endif diff --git a/src/Interaction/Models/contactForce/nonLinearCF.H b/src/Interaction/Models/contactForce/nonLinearCF.H new file mode 100644 index 00000000..fbd53152 --- /dev/null +++ b/src/Interaction/Models/contactForce/nonLinearCF.H @@ -0,0 +1,300 @@ +/*------------------------------- 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 __nonLinearCF_H__ +#define __nonLinearCF_H__ + +#include "types.H" + +namespace pFlow::cfModels +{ + +template +class nonLinear +{ +public: + + struct contactForceStorage + { + realx3 overlap_t_ = 0.0; + }; + + struct nonLinearProperties + { + real Yeff_ = 1000000.0; + real Geff_ = 8000000.0; + real ethan_ = 0.0; + real ethat_ = 0.0; + real mu_ = 0.00001; + + INLINE_FUNCTION_HD + nonLinearProperties(){} + + INLINE_FUNCTION_HD + nonLinearProperties(real Yeff, real Geff, real etha_n, real etha_t, real mu ): + Yeff_(Yeff), Geff_(Geff), ethan_(etha_n),ethat_(etha_t), mu_(mu) + {} + + INLINE_FUNCTION_HD + nonLinearProperties(const nonLinearProperties&)=default; + + INLINE_FUNCTION_HD + nonLinearProperties& operator=(const nonLinearProperties&)=default; + + INLINE_FUNCTION_HD + ~nonLinearProperties() = default; + }; + +protected: + + using NonLinearArrayType = symArray; + + int32 numMaterial_ = 0; + + ViewType1D rho_; + + NonLinearArrayType nonlinearProperties_; + + bool readNonLinearDictionary(const dictionary& dict) + { + auto Yeff = dict.getVal("Yeff"); + auto Geff = dict.getVal("Geff"); + auto nu = dict.getVal("nu"); + auto en = dict.getVal("en"); + auto et = dict.getVal("et"); + auto mu = dict.getVal("mu"); + + auto nElem = Yeff.size(); + + if(nElem != nu.size()) + { + fatalErrorInFunction<< + "sizes of Yeff("< prop(nElem); + forAll(i,Yeff) + { + prop[i] = {Yeff[i], Geff[i], etha_n[i], etha_t[i], mu[i]}; + } + + nonlinearProperties_.assign(prop); + + return true; + + } + + static const char* modelName() + { + if constexpr (limited) + { + return "nonLinearLimited"; + } + else + { + return "nonLinearNonLimited"; + } + return ""; + } + +public: + + TypeNameNV(modelName()); + + INLINE_FUNCTION_HD + nonLinear(){} + + nonLinear( + int32 nMaterial, + const ViewType1D& rho, + const dictionary& dict) + : + numMaterial_(nMaterial), + rho_("rho",nMaterial), + nonlinearProperties_("nonLinearProperties",nMaterial) + { + + Kokkos::deep_copy(rho_,rho); + if(!readNonLinearDictionary(dict)) + { + fatalExit; + } + } + + INLINE_FUNCTION_HD + nonLinear(const nonLinear&) = default; + + INLINE_FUNCTION_HD + nonLinear(nonLinear&&) = default; + + INLINE_FUNCTION_HD + nonLinear& operator=(const nonLinear&) = default; + + INLINE_FUNCTION_HD + nonLinear& operator=(nonLinear&&) = default; + + + INLINE_FUNCTION_HD + ~nonLinear()=default; + + INLINE_FUNCTION_HD + int32 numMaterial()const + { + return numMaterial_; + } + + //// - Methods + + INLINE_FUNCTION_HD + void contactForce + ( + const real dt, + const int32 i, + const int32 j, + const int32 propId_i, + const int32 propId_j, + const real Ri, + const real Rj, + const real ovrlp_n, + const realx3& Vr, + const realx3& Nij, + contactForceStorage& history, + realx3& FCn, + realx3& FCt + )const + { + + auto prop = nonlinearProperties_(propId_i,propId_j); + + + real vrn = dot(Vr, Nij); + realx3 Vt = Vr - vrn*Nij; + + history.overlap_t_ += Vt*dt; + + real mi = 3*Pi/4*pow(Ri,static_cast(3))*rho_[propId_i]; + real mj = 3*Pi/4*pow(Rj,static_cast(3))*rho_[propId_j]; + real Reff = 1.0/(1/Ri + 1/Rj); + + real K_hertz = 4.0/3.0*prop.Yeff_*sqrt(Reff); + real sqrt_meff_K_hertz = sqrt((mi*mj)/(mi+mj) * K_hertz); + + + //FCn = (-prop.kn_ * ovrlp_n - sqrt_meff * prop.ethan_ * vrn)*Nij; + //FCt = -prop.kt_ * history.overlap_t_ - sqrt_meff * prop.ethat_*Vt; + + + + FCn = (static_cast(-4.0/3.0) * prop.Yeff_ * sqrt(Reff)* pow(ovrlp_n,static_cast(1.5)) - + sqrt_meff_K_hertz*prop.ethan_*pow(ovrlp_n,static_cast(0.25))*vrn)*Nij; + + FCt = (- static_cast(16.0/3.0) * prop.Geff_ * sqrt(Reff*ovrlp_n) ) * history.overlap_t_; + + real ft = length(FCt); + real ft_fric = prop.mu_ * length(FCn); + + // apply friction + if(ft > ft_fric) + { + if( length(history.overlap_t_) >0.0) + { + if constexpr (limited) + { + real kt = static_cast(16.0/3.0) * prop.Geff_ * sqrt(Reff*ovrlp_n); + FCt *= (ft_fric/ft); + history.overlap_t_ = - (FCt/kt); + } + else + { + FCt = (FCt/ft)*ft_fric; + } + + } + else + { + FCt = 0.0; + } + } + + } + + +}; + +} //pFlow::CFModels + +#endif diff --git a/src/Interaction/Models/contactForceModels.H b/src/Interaction/Models/contactForceModels.H new file mode 100644 index 00000000..b359cfa2 --- /dev/null +++ b/src/Interaction/Models/contactForceModels.H @@ -0,0 +1,44 @@ +/*------------------------------- 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 __contactForceModels_H__ +#define __contactForceModels_H__ + + +#include "linearCF.H" +#include "nonLinearCF.H" +#include "normalRolling.H" + + +namespace pFlow::cfModels +{ + + +using limitedLinearNormalRolling = normalRolling>; +using nonLimitedLinearNormalRolling = normalRolling>; + +using limitedNonLinearNormalRolling = normalRolling>; +using nonLimitedNonLinearNormalRolling = normalRolling>; + +} + + + +#endif //__contactForceModels_H__ \ No newline at end of file diff --git a/src/Interaction/Models/rolling/normalRolling.H b/src/Interaction/Models/rolling/normalRolling.H new file mode 100644 index 00000000..ce602838 --- /dev/null +++ b/src/Interaction/Models/rolling/normalRolling.H @@ -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. + +-----------------------------------------------------------------------------*/ + +#ifndef __normalRolling_H__ +#define __normalRolling_H__ + + +namespace pFlow::cfModels +{ + +template +class normalRolling +: + public contactForceModel +{ +public: + + using contactForceStorage = + typename contactForceModel::contactForceStorage; + + + realSymArray_D mur_; + + bool readNormalDict(const dictionary& dict) + { + auto mur = dict.getVal("mur"); + + uint32 nMat; + + if(!realSymArray_D::getN(mur.size(),nMat) || nMat != this->numMaterial()) + { + fatalErrorInFunction<< + "wrong number of values supplied in mur. \n"; + return false; + } + + mur_.assign(mur); + + return true; + } + +public: + + TypeNameNV(word("normal<"+contactForceModel::TYPENAME()+">")); + + + normalRolling(int32 nMaterial, const ViewType1D& rho, const dictionary& dict) + : + contactForceModel(nMaterial, rho, dict), + mur_("mur", nMaterial) + { + if(!readNormalDict(dict)) + { + fatalExit; + } + } + + INLINE_FUNCTION_HD + void rollingFriction + ( + const real dt, + const int32 i, + const int32 j, + const int32 propId_i, + const int32 propId_j, + const real Ri, + const real Rj, + const realx3& wi, + const realx3& wj, + const realx3& Nij, + const realx3& FCn, + realx3& Mri, + realx3& Mrj + )const + { + + realx3 w_hat = wi-wj; + real w_hat_mag = length(w_hat); + + if( !equal(w_hat_mag,0.0) ) + w_hat /= w_hat_mag; + else + w_hat = 0.0; + + auto Reff = (Ri*Rj)/(Ri+Rj); + + Mri = ( -mur_(propId_i,propId_j) *length(FCn) * Reff ) * w_hat ; + + //removing the normal part + // Mri = Mri - ( (Mri .dot. nij)*nij ) + + Mrj = -Mri; + } + + +}; + +} + +#endif diff --git a/src/Interaction/contactLists/sortedContactList.H b/src/Interaction/contactLists/sortedContactList.H new file mode 100644 index 00000000..0c9a44fb --- /dev/null +++ b/src/Interaction/contactLists/sortedContactList.H @@ -0,0 +1,171 @@ +/*------------------------------- 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 __sortedContactList_H__ +#define __sortedContactList_H__ + +#include "sortedPairs.H" + +namespace pFlow +{ + +template +class sortedContactList +: + public sortedPairs +{ +public: + + using SortedPairs = sortedPairs; + + using ValueType = valueType; + + using IdType = typename SortedPairs::IdType; + + using ExecutionSpace = typename SortedPairs::ExecutionSpace; + + using memory_space = typename SortedPairs::memory_space; + + using PairType = typename SortedPairs::PairType; + + using ContainerType = typename SortedPairs::ContainerType; + + class TagReFillPairs{}; + +protected: + + ViewType1D values_; + + int32 size0_ = 0; + + ViewType1D sortedPairs0_; + + ViewType1D values0_; + + + void adjustCapacity() + { + if(auto s = this->size(); s > values_.size()) + { + reallocNoInit(values_, s); + } + } + + + + using rpReFillPairs = Kokkos::RangePolicy< + ExecutionSpace, + Kokkos::Schedule, + Kokkos::IndexType, + TagReFillPairs>; + +public: + + TypeNameNV("sortedContactList"); + + + sortedContactList(int32 initialSize =1) + : + SortedPairs(initialSize), + values_("values", SortedPairs::capacity()), + sortedPairs0_("sortedPairs0", SortedPairs::capacity()), + values0_("values0", SortedPairs::capacity()) + {} + + bool beforeBroadSearch() + { + swapViews(values0_, values_); + swapViews(sortedPairs0_, this->sortedPairs_); + size0_ = this->size(); + return SortedPairs::beforeBroadSearch(); + } + + bool afterBroadSearch() + { + SortedPairs::afterBroadSearch(); + + adjustCapacity(); + + Kokkos::parallel_for( + "sortedContactList::reFillPairs", + rpReFillPairs(0, this->size() ), + *this + ); + Kokkos::fence(); + + return true; + } + + INLINE_FUNCTION_HD + ValueType getValue(int32 idx)const + { + return values_[idx]; + } + + INLINE_FUNCTION_HD + void setValue(int32 idx, const ValueType& val)const + { + values_[idx] = val; + } + + INLINE_FUNCTION_HD + void operator()(TagReFillPairs, int32 idx)const + { + + auto searchLen = max(size0_/1000,10); + auto start = max(0,idx-searchLen); + auto end = min(size0_,idx+searchLen); + auto newPair = this->sortedPairs_[idx]; + if( auto idx0 = binarySearch( + sortedPairs0_, + start, + end, + newPair); + idx0>=0) + { + values_[idx] = values0_[idx0]; + } + else if(auto idx0 = binarySearch( + sortedPairs0_, + start, + end, + newPair); + idx0>=0) + { + values_[idx] = values0_[idx0]; + + } + else + { + values_[idx] = ValueType(); + } + } + + + +}; // sortedContactList + + + +} // pFlow + + + +#endif //__sortedContactList_H__ diff --git a/src/Interaction/contactLists/sortedPairs.H b/src/Interaction/contactLists/sortedPairs.H new file mode 100644 index 00000000..e127fcfe --- /dev/null +++ b/src/Interaction/contactLists/sortedPairs.H @@ -0,0 +1,258 @@ +/*------------------------------- 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 __sortedPairs_H__ +#define __sortedPairs_H__ + +#include "unsortedPairs.H" +#include "KokkosUtilities.H" + +namespace pFlow +{ + + +template +class sortedPairs +: + public unsortedPairs +{ +public: + + using UnsortedPairs = unsortedPairs; + + using IdType = typename UnsortedPairs::IdType; + + using ExecutionSpace = typename UnsortedPairs::ExecutionSpace; + + using memory_space = typename ExecutionSpace::memory_space; + + using PairType = typename UnsortedPairs::PairType; + + using ContainerType = typename UnsortedPairs::ContainerType; + + struct pairAccessor + { + using PairType = typename sortedPairs::PairType; + + int32 size_; + + ViewType1D sortedParis_; + + INLINE_FUNCTION_HD + int32 size()const { return size_; } + + INLINE_FUNCTION_HD + int32 loopCount()const { return size_; } + + INLINE_FUNCTION_HD + bool isValid(int32 i)const { return i flags_; + + ViewType1D sortedPairs_; + + using rpFillFlag = Kokkos::RangePolicy< + ExecutionSpace, + Kokkos::Schedule, + Kokkos::IndexType, + TagFillFlag >; + + using rpFillPairs = Kokkos::RangePolicy< + ExecutionSpace, + Kokkos::Schedule, + Kokkos::IndexType, + TagFillPairs>; + +public: + + // - type info + TypeNameNV("sortedPairs"); + + + // constructors + sortedPairs(int32 initialSize =1) + : + UnsortedPairs(initialSize), + flags_("flags_",UnsortedPairs::capacity()+1), + sortedPairs_("sortedPairs_",UnsortedPairs::capacity()) + {} + + + bool beforeBroadSearch() + { + this->clear(); + return true; + } + + bool afterBroadSearch() + { + prepareSorted(); + return true; + } + + // - Device call + // return the pair at index idx + // perform no check for size and existance + INLINE_FUNCTION_HD + PairType getPair(int32 idx)const + { + return sortedPairs_[idx]; + } + + // - Device/host call + // return the pair at index idx + INLINE_FUNCTION_HD + bool getPair(int32 idx, PairType& p)const + { + if(isValid(idx)) + { + p = getPair(idx); + return true; + } + else + { + return false; + } + } + + INLINE_FUNCTION_HD + bool isValid(int32 idx)const + { + return idx < size_; + } + + + //use this when the value of size_ is updated + INLINE_FUNCTION_H + int32 size()const + { + return size_; + } + + int32 loopCount()const + { + return size_; + } + + pairAccessor getPairs()const + { + return {size_, sortedPairs_}; + } + + INLINE_FUNCTION_H + void clear() + { + UnsortedPairs::clear(); + size_ = 0; + } + + void prepareSorted() + { + // first check the size of flags_ + int32 capacity = UnsortedPairs::capacity(); + + if( capacity+1 > flags_.size() ) + { + reallocNoInit(flags_, capacity+1); + } + + // fill the flags + Kokkos::parallel_for( + "contactPairsSorted::fillFlag", + rpFillFlag(0,capacity+1), + *this); + Kokkos::fence(); + + + // inclusive scan on flags_ + exclusiveScan(flags_, 0, capacity+1,flags_,0); + + Kokkos::fence(); // Kokkos's scan is blocking I guess. So, this could be removed. + + // get the last value of flags_ to obtain the size of sortedPairs_ + getNth(size_, flags_, capacity); + + if(size_ == 0 )return; + + // resize sortedPairs_ if necessary; + if( size_>sortedPairs_.size() ) + { + // get more space to prevent reallocations in next iterations + int32 len = size_*1.1+1; + reallocNoInit(sortedPairs_, len); + } + + Kokkos::parallel_for( + "contactPairsSorted::fillPairs", + rpFillPairs(0,this->capacity()), + *this); + Kokkos::fence(); + + // - sort paris based on the first and second + sort(sortedPairs_, 0, size_ ); + + + } + + INLINE_FUNCTION_HD + void operator()(TagFillFlag, int32 i)const + { + if(this->container_.valid_at(i) ) + flags_[i] = 1; + else + flags_[i] = 0; + } + + INLINE_FUNCTION_HD + void operator()(TagFillPairs, int32 i)const + { + auto fi = flags_[i]; + if(fi!=flags_[i+1]) + sortedPairs_[fi] = this->container_.key_at(i); + } + +}; + +} + +#endif //__sortedPairs_H__ diff --git a/src/Interaction/contactLists/unsortedContactList.H b/src/Interaction/contactLists/unsortedContactList.H new file mode 100644 index 00000000..f83960c6 --- /dev/null +++ b/src/Interaction/contactLists/unsortedContactList.H @@ -0,0 +1,175 @@ +/*------------------------------- 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 __unsortedContactList_H__ +#define __unsortedContactList_H__ + +namespace pFlow +{ + + +template +class unsortedContactList +: + public unsortedPairs +{ +public: + + using ValueType = valueType; + + using UnsortedPairs = unsortedPairs; + + using IdType = typename UnsortedPairs::IdType; + + using ExecutionSpace = typename UnsortedPairs::ExecutionSpace; + + using memory_space = typename ExecutionSpace::memory_space; + + using PairType = typename UnsortedPairs::PairType; + + using ContainerType = typename UnsortedPairs::ContainerType; + + class TagReFillPairs{}; + +protected: + + /// storage for keeping the values of the current list + ViewType1D values_; + + /// storage for keeping pairs from the previous list + ContainerType container0_; + + /// storage for keeping values from the previous list + ViewType1D values0_; + + void adjustCapacity() + { + auto cap = this->capacity(); + if(cap > values_.size()) + { + reallocNoInit(values_, cap); + } + + } + + using rpFillPairs = Kokkos::RangePolicy< + ExecutionSpace, + Kokkos::Schedule, + Kokkos::IndexType, + TagReFillPairs>; + + +public: + + TypeNameNV("unsortedContactList"); + + unsortedContactList(int32 capacity=1) + : + UnsortedPairs(capacity), + values_("values", UnsortedPairs::capacity()), + container0_(capacity), + values0_("values0",container0_.capacity()) + {} + + + bool beforeBroadSearch() + { + // swap conainer and values + swapViews(values0_, values_); + swapViews(container0_, this->container_); + return UnsortedPairs::beforeBroadSearch(); + } + + bool afterBroadSearch() + { + + UnsortedPairs::afterBroadSearch(); + + adjustCapacity(); + + Kokkos::parallel_for( + "reFillPairs", + rpFillPairs(0,this->capacity()), + *this); + Kokkos::fence(); + + return true; + } + + INLINE_FUNCTION_HD + ValueType getValue(int32 idx)const + { + return values_[idx]; + } + + INLINE_FUNCTION_HD + bool getValue(const PairType& p, ValueType& val)const + { + if(auto idx = this->find(p); idx>=0) + { + val = getValue(idx); + return true; + } + return false; + } + + INLINE_FUNCTION_HD + void setValue(int32 idx, const ValueType& val)const + { + values_[idx] = val; + } + + INLINE_FUNCTION_HD + bool setValue(const PairType& p, const ValueType& val)const + { + if(auto idx = this->find(p); idx>=0) + { + setValue(idx, val); + return true;; + } + return false; + } + + INLINE_FUNCTION_HD + void operator()(TagReFillPairs, int32 idx)const + { + if( this->isValid(idx) ) + { + if( int32 idx0 = + container0_.find(this->getPair(idx)); + idx0>=0 ) + { + values_[idx] = values0_[idx0]; + } + else + { + values_[idx] = ValueType(); + } + } + // invalid locations should not be filled. + } + +}; //unsortedContactList + + +} // pFlow + + +#endif //__unsortedContactList_H__ diff --git a/src/Interaction/contactLists/unsortedPairs.H b/src/Interaction/contactLists/unsortedPairs.H new file mode 100644 index 00000000..4905f1c3 --- /dev/null +++ b/src/Interaction/contactLists/unsortedPairs.H @@ -0,0 +1,216 @@ +/*------------------------------- 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 __unsortedPairs_H__ +#define __unsortedPairs_H__ + +#include "KokkosTypes.H" +#include "types.H" + +namespace pFlow +{ + + +template +class unsortedPairs +{ +public: + + using UnsortedPairs = unsortedPairs; + + using IdType = idType; + + using ExecutionSpace = executionSpace; + + using memory_space = typename ExecutionSpace::memory_space; + + using PairType = kPair; + + using ContainerType = unorderedSet; + + struct pairAccessor + { + using PairType = typename UnsortedPairs::PairType; + + ContainerType Container_; + + INLINE_FUNCTION_HD + int32 size()const { return Container_.size(); } + + INLINE_FUNCTION_HD + int32 loopCount()const { return Container_.capacity(); } + + INLINE_FUNCTION_HD + bool isValid(int32 idx)const { return Container_.valid_at(idx); } + + INLINE_FUNCTION_HD + PairType getPair(int idx)const { return Container_.key_at(idx); } + + INLINE_FUNCTION_HD + bool getPair(int32 idx, PairType& pair)const { + if(Container_.valid_at(idx)) { + pair = Container_.key_at(idx); + return true; + } + return false; + } + }; + +protected: + + ContainerType container_; + + +public: + + // - type info + TypeNameNV("unsorderedPairs"); + + // constructor + unsortedPairs(int32 capacity=1) + : + container_(capacity) // the minimum capacity would be 128 + {} + + bool beforeBroadSearch() + { + container_.clear(); + return true; + } + + bool afterBroadSearch() + { + return true; + } + + // - Device call + INLINE_FUNCTION_HD + int32 insert(idType i, idType j)const + { + if(auto insertResult = container_.insert(PairType(i,j)); insertResult.failed()) + return -1; + else + return insertResult.index(); + + } + + INLINE_FUNCTION_HD + int32 insert(const PairType& p)const + { + if(auto insertResult = container_.insert(p); insertResult.failed()) + return -1; + else + return insertResult.index(); + + } + + // - Device call + // return the pair at index idx + // perform no check for size and existance + INLINE_FUNCTION_HD + PairType getPair(int32 idx)const + { + return container_.key_at(idx); + } + + // - Device call + // return the pair at index idx + INLINE_FUNCTION_HD + bool getPair(int32 idx, PairType& p)const + { + if(container_.valid_at(idx)) + { + p = container_.key_at(idx); + return true; + } + else + { + + return false; + } + } + + INLINE_FUNCTION_HD + int32 find(const PairType & p)const + { + if( auto idx = container_.find(p); + idx != Kokkos::UnorderedMapInvalidIndex ) + return idx; + else + return -1; + } + + INLINE_FUNCTION_HD + bool isValid(int32 idx)const + { + return container_.valid_at(idx); + } + + + INLINE_FUNCTION_HD + int32 capacity() const + { + return container_.capacity(); + } + + int32 loopCount()const + { + return container_.capacity(); + } + + //use this when the value of size_ is updated + INLINE_FUNCTION_H + int32 size()const + { + return container_.size(); + } + + pairAccessor getPairs()const + { + return {container_}; + } + + /// increase the capacity of the container by at-least len + /// the content will be erased. + INLINE_FUNCTION_H + void increaseCapacityBy(int32 len) + { + uint newCap = container_.capacity()+len; + this->clear(); + container_.rehash(newCap); + } + + INLINE_FUNCTION_H + void clear() + { + container_.clear(); + } + + const ContainerType& container()const + { + return container_; + } + +}; + + +} + +#endif //__unsortedPairs_H__ \ No newline at end of file diff --git a/src/Interaction/contactSearch/ContactSearch/ContactSearch.H b/src/Interaction/contactSearch/ContactSearch/ContactSearch.H new file mode 100644 index 00000000..11a7f90c --- /dev/null +++ b/src/Interaction/contactSearch/ContactSearch/ContactSearch.H @@ -0,0 +1,250 @@ +/*------------------------------- 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 __ContactSearch_H__ +#define __ContactSearch_H__ + + +#include "contactSearch.H" +#include "box.H" + +namespace pFlow +{ + +template< + template class BaseMethod, + template class WallMapping +> +class ContactSearch +: + public contactSearch +{ +public: + + using IdType = typename contactSearch::IdType; + + using IndexType = typename contactSearch::IndexType; + + using ExecutionSpace = typename contactSearch::ExecutionSpace; + + using PairContainerType = typename contactSearch::PairContainerType; + + using ParticleContactSearchType = + BaseMethod< + ExecutionSpace, + IdType, + IndexType + >; + + using WallMappingType = + WallMapping< + ExecutionSpace, + IdType, + IndexType + >; + +protected: + + + uniquePtr particleContactSearch_ = nullptr; + + uniquePtr wallMapping_ = nullptr; + +public: + + TypeNameTemplate2("ContactSearch", ParticleContactSearchType, WallMappingType); + + ContactSearch( + const dictionary& csDict, + const box& domain, + const particles& prtcl, + const geometry& geom, + Timers& timers) + : + contactSearch(csDict, domain, prtcl, geom, timers) + + { + + auto method = dict().getVal("method"); + auto wmMethod = dict().getVal("wallMapping"); + + auto nbDict = dict().subDictOrCreate(method+"Info"); + + real minD, maxD; + this->Particles().boundingSphereMinMax(minD, maxD); + + const auto& position = this->Particles().pointPosition().deviceVectorAll(); + const auto& diam = this->Particles().boundingSphere().deviceVectorAll(); + + particleContactSearch_ = + makeUnique + ( + nbDict, + this->domain(), + maxD, + position, + diam + ); + Report(2)<<"Contact search algorithm for particle-particle is "<< + greenText(particleContactSearch_().typeName())<Geometry().numPoints(); + int32 wnTri = this->Geometry().size(); + + const auto& wPoints = this->Geometry().points().deviceVectorAll(); + const auto& wVertices = this->Geometry().vertices().deviceVectorAll(); + + wallMapping_ = + makeUnique( + wmDict, + particleContactSearch_().getCells(), + wnPoints, + wnTri, + wPoints, + wVertices + ); + Report(2)<<"Wall mapping algorithm for particle-wall is "<< + greenText(wallMapping_().typeName())<< endReport; + + } + + + add_vCtor( + contactSearch, + ContactSearch, + dictionary); + + bool broadSearch( + PairContainerType& ppPairs, + PairContainerType& pwPairs, + bool force = false) override + { + + //output<<" ContactSearch::broadSearch::PP before.\n"; + if(particleContactSearch_) + { + 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 + return false; + + if(wallMapping_) + { + sphereWallTimer_.start(); + wallMapping_().broadSearch(pwPairs, particleContactSearch_(), force); + sphereWallTimer_.end(); + } + else + return false; + + return true; + } + + + bool ppEnterBroadSearch()const override + { + if(particleContactSearch_) + { + return particleContactSearch_().enterBoadSearch(); + } + return false; + } + + bool pwEnterBroadSearch()const override + { + if(wallMapping_) + { + return wallMapping_().enterBoadSearch(); + } + 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; + }*/ + + +}; + +} + + +#endif //__ContactSearch_H__ diff --git a/src/Interaction/contactSearch/ContactSearch/ContactSearchs.C b/src/Interaction/contactSearch/ContactSearch/ContactSearchs.C new file mode 100644 index 00000000..75343dfa --- /dev/null +++ b/src/Interaction/contactSearch/ContactSearch/ContactSearchs.C @@ -0,0 +1,26 @@ +/*------------------------------- 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 "ContactSearch.H" + +#include "cellsSimple.H" +#include "NBS.H" + +template class pFlow::ContactSearch; diff --git a/src/Interaction/contactSearch/cells.H b/src/Interaction/contactSearch/cells.H new file mode 100644 index 00000000..f0513d20 --- /dev/null +++ b/src/Interaction/contactSearch/cells.H @@ -0,0 +1,233 @@ +/*------------------------------- 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 __cells_H__ +#define __cells_H__ + + +#include "types.H" +#include "box.H" + +namespace pFlow +{ + +template +class cells +{ +public: + + using CellType = triple; + +protected: + + // - domain + box domain_; + + // - cell size + real cellSize_; + + CellType numCells_; + + + // - protected methods + INLINE_FUNCTION_H + void calculate() + { + numCells_ = (domain_.maxPoint()-domain_.minPoint())/cellSize_ + realx3(1.0); + numCells_ = max( numCells_ , CellType(static_cast(1)) ); + } + +public: + + INLINE_FUNCTION_HD + cells() + {} + + INLINE_FUNCTION_H + cells(const box& domain, real cellSize) + : + domain_(domain), + cellSize_(cellSize) + { + calculate(); + } + + INLINE_FUNCTION_HD + cells(const cells&) = default; + + INLINE_FUNCTION_HD + cells& operator = (const cells&) = default; + + cells getCells()const + { + return *this; + } + + INLINE_FUNCTION_H + void setCellSize(real cellSize) + { + cellSize_ = cellSize; + calculate(); + } + + INLINE_FUNCTION_HD + real cellSize()const + { + return cellSize_; + } + + INLINE_FUNCTION_HD + const CellType& numCells()const + { + return numCells_; + } + + INLINE_FUNCTION_HD + indexType nx()const + { + return numCells_.x(); + } + + INLINE_FUNCTION_HD + indexType ny()const + { + return numCells_.y(); + } + + INLINE_FUNCTION_HD + indexType nz()const + { + return numCells_.z(); + } + + INLINE_FUNCTION_HD + int64 totalCells()const + { + return static_cast(numCells_.x())* + static_cast(numCells_.y())* + static_cast(numCells_.z()); + } + + INLINE_FUNCTION_HD + CellType pointIndex(const realx3& p)const + { + return CellType( (p - domain_.minPoint())/cellSize_ ); + } + + INLINE_FUNCTION_HD + bool pointIndexInDomain(const realx3 p, CellType& index)const + { + if( !domain_.isInside(p) ) return false; + + index = this->pointIndex(p); + return true; + } + + INLINE_FUNCTION_HD + bool inDomain(const realx3& p)const + { + return domain_.isInside(p); + } + + INLINE_FUNCTION_HD + bool isInRange(const CellType& cell)const + { + if(cell.x()<0)return false; + if(cell.y()<0)return false; + if(cell.z()<0)return false; + if(cell.x()>numCells_.x()-1) return false; + if(cell.y()>numCells_.y()-1) return false; + if(cell.z()>numCells_.z()-1) return false; + return true; + } + + INLINE_FUNCTION_HD + bool isInRange(indexType i, indexType j, indexType k)const + { + if(i<0)return false; + if(j<0)return false; + if(k<0)return false; + if(i>numCells_.x()-1) return false; + if(j>numCells_.y()-1) return false; + if(k>numCells_.z()-1) return false; + 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( + const realx3& p1, + const realx3& p2, + const realx3& p3, + real extent, + realx3& minP, + realx3& maxP)const + { + minP = min(min(p1,p2),p3) - extent*cellSize_ ; + maxP = max(max(p1,p2),p3) + extent*cellSize_ ; + + minP = bound(minP); + maxP = bound(maxP); + } + + INLINE_FUNCTION_HD + CellType bound(CellType p)const + { + return CellType( + min( numCells_.x()-1, max(0,p.x())), + min( numCells_.y()-1, max(0,p.y())), + min( numCells_.z()-1, max(0,p.z())) + ); + } + + INLINE_FUNCTION_HD + 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())) + ); + } + + + +}; + + +} + + +#endif diff --git a/src/Interaction/contactSearch/contactSearch/contactSearch.C b/src/Interaction/contactSearch/contactSearch/contactSearch.C new file mode 100644 index 00000000..8c05a03e --- /dev/null +++ b/src/Interaction/contactSearch/contactSearch/contactSearch.C @@ -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. + +-----------------------------------------------------------------------------*/ + +#include "contactSearch.H" + + + +pFlow::contactSearch::contactSearch( + const dictionary& dict, + const box& domain, + 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) +{ + +} + + +pFlow::uniquePtr pFlow::contactSearch::create( + const dictionary& dict, + const box& domain, + const particles& prtcl, + const geometry& geom, + Timers& timers) +{ + + word baseMethName = dict.getVal("method"); + word wallMethod = dict.getVal("wallMapping"); + + auto model = angleBracketsNames2("ContactSearch", baseMethName, wallMethod); + + + Report(1)<<"Selecting contact search model . . ."<; + +protected: + + const box& domain_; + + dictionary dict_; + + Timer sphereSphereTimer_; + + Timer sphereWallTimer_; + + auto& dict() + { + return dict_; + } + +public: + + TypeName("contactSearch"); + + contactSearch( + const dictionary& dict, + const box& domain, + const particles& prtcl, + const geometry& geom, + Timers& timers); + + virtual ~contactSearch()=default; + + + create_vCtor + ( + contactSearch, + dictionary, + ( + const dictionary& dict, + const box& domain, + const particles& prtcl, + const geometry& geom, + Timers& timers + ), + (dict, domain, prtcl, geom, timers) + ); + + const auto& domain()const + { + return domain_; + } + + const auto& dict()const + { + return dict_; + } + + + virtual + bool broadSearch( + PairContainerType& ppPairs, + PairContainerType& pwPairs, + bool force = false) = 0; + + virtual + bool ppEnterBroadSearch()const = 0; + + virtual + bool pwEnterBroadSearch()const = 0; + + virtual + bool ppPerformedBroadSearch()const = 0; + + virtual + bool pwPerformedBroadSearch()const = 0; + + + static + uniquePtr create( + const dictionary& dict, + const box& domain, + const particles& prtcl, + const geometry& geom, + Timers& timers); + +}; + + +} + + +#endif //__ContactSearch_H__ diff --git a/src/Interaction/contactSearch/contactSearchFunctions.H b/src/Interaction/contactSearch/contactSearchFunctions.H new file mode 100644 index 00000000..10eadce1 --- /dev/null +++ b/src/Interaction/contactSearch/contactSearchFunctions.H @@ -0,0 +1,112 @@ +/*------------------------------- 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 __broadSearchFunctions_H__ +#define __broadSearchFunctions_H__ + +#include "types.H" +#include "iBox.H" + +namespace pFlow +{ + +INLINE_FUNCTION_HD +uint64_t splitBy3(const uint64_t val){ + uint64_t x = val; + x = (x | x << 32) & 0x1f00000000ffff; + x = (x | x << 16) & 0x1f0000ff0000ff; + x = (x | x << 8) & 0x100f00f00f00f00f; + x = (x | x << 4) & 0x10c30c30c30c30c3; + x = (x | x << 2) & 0x1249249249249249; + return x; +} + +INLINE_FUNCTION_HD +uint64_t xyzToMortonCode64(uint64_t x, uint64_t y, uint64_t z) +{ + return splitBy3(x) | (splitBy3(y) << 1) | (splitBy3(z) << 2); +} + + +INLINE_FUNCTION_HD +uint64_t getThirdBits(uint64_t x) +{ + x = x & 0x9249249249249249; + x = (x | (x >> 2)) & 0x30c30c30c30c30c3; + x = (x | (x >> 4)) & 0xf00f00f00f00f00f; + x = (x | (x >> 8)) & 0x00ff0000ff0000ff; + x = (x | (x >> 16)) & 0xffff00000000ffff; + x = (x | (x >> 32)) & 0x00000000ffffffff; + return x; + +} + +INLINE_FUNCTION_HD +void mortonCode64Toxyz(uint64_t morton, uint64_t& x, uint64_t& y, uint64_t& z) +{ + x = getThirdBits(morton); + y = getThirdBits(morton >> 1); + z = getThirdBits(morton >> 2); +} + +template +INLINE_FUNCTION_HD +void indexToCell(const indexType idx, const triple& extent, triple& cell) +{ + indexType nxny = extent.x()*extent.y(); + cell.z() = static_cast(idx / nxny); + auto rem = idx % nxny; + cell.y() = static_cast(rem / extent.x()); + cell.x() = static_cast(rem % extent.x()); +} + +template +INLINE_FUNCTION_HD +triple boxExtent( const iBox& box) +{ + return triple( + box.maxPoint().x() - box.minPoint().x() + 1, + box.maxPoint().y() - box.minPoint().y() + 1, + box.maxPoint().z() - box.minPoint().z() + 1); +} + +template +INLINE_FUNCTION_HD +void indexToCell(const indexType idx, const iBox& box, triple& cell) +{ + auto extent = boxExtent(box); + indexToCell(idx, extent, cell); + 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_H__ diff --git a/src/Interaction/contactSearch/methods/NBS.H b/src/Interaction/contactSearch/methods/NBS.H new file mode 100644 index 00000000..a454819f --- /dev/null +++ b/src/Interaction/contactSearch/methods/NBS.H @@ -0,0 +1,522 @@ +/*------------------------------- phasicFlow --------------------------------- + O C enter of + O O E ngineering and + O O M ultiscale modeling of + OOOOOOO F luid flow +------------------------------------------------------------------------------ + Copyright (C): www.cemf.ir + email: hamid.r.norouzi AT gmail.com +------------------------------------------------------------------------------ +Licence: + This file is part of phasicFlow code. It is a free software for simulating + granular and multiphase flows. You can redistribute it and/or modify it under + the terms of GNU General Public License v3 or any other later versions. + + phasicFlow is distributed to help others in their research in the field of + granular and multiphase flows, but WITHOUT ANY WARRANTY; without even the + implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. + +-----------------------------------------------------------------------------*/ + + +#ifndef __NBS_H__ +#define __NBS_H__ + +#include "cells.H" +#include "contactSearchFunctions.H" + +namespace pFlow +{ + + +template< + typename executionSpace, + typename idType, + typename indexType=int32 + > +class NBS +: + public cells +{ +public: + + using IdType = idType; + + using IndexType = indexType; + + using Cells = cells; + + using CellType = typename Cells::CellType; + + using ExecutionSpace= executionSpace; + + using memory_space = typename ExecutionSpace::memory_space; + + struct TagBuild{}; + + struct TagFindPairs{}; + +protected: + + int32 capacity_ = 1; + + real sizeRatio_ = 1.0; + + int32 updateFrequency_= 1; + + int32 currentIter_ = 0; + + bool performedSearch_ = false; + + ViewType1D pointPosition_; + + ViewType1D diameter_; + + ViewType3D head_; + + ViewType1D next_; + + + INLINE_FUNCTION_H + void nullify() + { + fill( + head_, + range(0,this->nx()), + range(0,this->ny()), + range(0,this->nz()), + -1 + ); + + fill( + next_, + range(0,capacity_), + -1 + ); + } + + void nullify(range nextRng) + { + fill( + head_, + range(0,this->nx()), + range(0,this->ny()), + range(0,this->nz()), + -1 + ); + + fill( + next_, + nextRng, + -1 + ); + } + + using mdrPolicyFindPairs = + Kokkos::MDRangePolicy< + Kokkos::Rank<3>, + Kokkos::Schedule, + ExecutionSpace>; + +private: + + void checkAllocateNext(int newCap) + { + + if( capacity_ < newCap) + { + capacity_ = newCap; + Kokkos::realloc(next_, capacity_); + } + + } + + void allocateHead() + { + Kokkos::realloc(head_, this->nx(), this->ny(), this->nz()); + } + + bool performSearch() + { + if(currentIter_ % updateFrequency_ == 0) + { + currentIter_++; + return true; + + }else + { + currentIter_++; + return false; + } + } + + static INLINE_FUNCTION_HD + void Swap(int32& x, int32& y) + { + int32 tmp = x; + x = y; + y = tmp; + } + +public: + + TypeNameNV("NBS"); + + NBS( + const box& domain, + real cellSize, + const ViewType1D& position, + const ViewType1D& diam, + int32 initialContainerSize = 1) + : + Cells(domain, cellSize), + pointPosition_(position), + diameter_(diam), + head_("NBS::head_",1,1,1), //, this->nx(), this->ny(), this->nz()), + next_("NBS::next_",1) //,position.size()), + { + checkAllocateNext(pointPosition_.size()); + allocateHead(); + } + + NBS( + dictionary dict, + const box& domain, + real cellSize, + const ViewType1D& position, + const ViewType1D& diam, + int32 initialContainerSize = 1 + ) + : + Cells(domain, cellSize), + pointPosition_(position), + diameter_(diam), + head_("NBS::head_",1,1,1), //, this->nx(), this->ny(), this->nz()), + next_("NBS::next_",1) //,position.size()), + { + + updateFrequency_ = max( + dict.getVal("updateFrequency"),1 ); + + sizeRatio_ = max(dict.getVal( + "sizeRatio"),1.0); + + this->setCellSize(cellSize*sizeRatio_); + checkAllocateNext(pointPosition_.size()); + allocateHead(); + } + + INLINE_FUNCTION_HD + NBS(const NBS&) = default; + + INLINE_FUNCTION_HD + NBS& operator = (const NBS&) = default; + + INLINE_FUNCTION_HD + ~NBS()=default; + + //// - Methods + + bool enterBoadSearch()const + { + return currentIter_%updateFrequency_==0; + } + + bool performedSearch()const + { + return performedSearch_; + } + + // - 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 + bool broadSearch(PairsContainer& pairs, range activeRange, bool force=false) + { + + if(force) currentIter_ = 0; + performedSearch_ = false; + + if( !performSearch() ) return true; + + build(activeRange); + + 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 + bool broadSearch(PairsContainer& pairs, range activeRange, IncludeFunction incld, bool force = false) + { + if(force) currentIter_ = 0; + performedSearch_ = false; + + if( !performSearch() ) return true; + + build(activeRange, incld); + + findPairs(pairs); + + performedSearch_ = true; + return true; + } + + // - build based on all points in range [0, numPoints_) + INLINE_FUNCTION_H + void build(range activeRange) + { + checkAllocateNext(activeRange.second); + nullify(activeRange); + + Cells cellIndex = static_cast(*this); + auto points = pointPosition_; + auto next = next_; + auto head = head_; + + Kokkos::RangePolicy< + Kokkos::IndexType, + Kokkos::Schedule, + ExecutionSpace> rPolicy(activeRange.first, activeRange.second); + + Kokkos::parallel_for( + "NBS::build", + rPolicy, + LAMBDA_HD(int32 i){ + + CellType ind = cellIndex.pointIndex(points[i]); + int32 old = Kokkos::atomic_exchange(&head(ind.x(), ind.y(), ind.z()), i); + next[i] = old; + }); + Kokkos::fence(); + + } + + + template + INLINE_FUNCTION_H + void build(range activeRange, IncludeFunction incld) + { + checkAllocateNext(activeRange.second); + nullify(activeRange); + + Cells cellIndex = static_cast(*this); + auto points = pointPosition_; + auto next = next_; + auto head = head_; + + Kokkos::RangePolicy< + Kokkos::IndexType, + Kokkos::Schedule, + ExecutionSpace> rPolicy(activeRange.first, activeRange.second); + + Kokkos::parallel_for( + "NBS::build", + rPolicy, + LAMBDA_HD(int32 i){ + + if( incld(i) ) + { + CellType ind = cellIndex.pointIndex(points[i]); + auto old = Kokkos::atomic_exchange(&head(ind.x(), ind.y(), ind.z()), i); + next[i] = old; + } + }); + Kokkos::fence(); + + } + + template + INLINE_FUNCTION_H + bool findPairs(PairsContainer& pairs) + { + mdrPolicyFindPairs + mdrPolicy( + {0,0,0}, + {this->nx(),this->ny(),this->nz()} ); + + int32 getFull = 1; + + // loop until the container size fits the numebr of contact pairs + while (getFull > 0) + { + + getFull = 0; + + Kokkos::parallel_reduce ( + "NBS::broadSearch", + mdrPolicy, + CLASS_LAMBDA_HD(int32 i, int32 j, int32 k, int32& getFullUpdate){ + #include "NBSLoop.H" + }, getFull); + + if(getFull) + { + // - resize the container + // note that getFull now shows the number of failed insertions. + uint32 len = max(getFull,50) ; + + auto oldCap = pairs.capacity(); + + pairs.increaseCapacityBy(len); + + Info<< "The contact pair container capacity increased from "<< + oldCap << " to "< -1 ) + { + + auto p_m = pointPosition_[m]; + auto d_m = sizeRatio_* diameter_[m]; + + // the same cell + n = next_(m); + + while(n >-1) + { + + auto p_n = 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 = 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 = head_(neighborCell.x(), neighborCell.y(), neighborCell.z()); + while( n>-1) + { + + auto p_n = 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 = next_[n]; + } + } + + } + m = next_[m]; + } + }*/ + + template + INLINE_FUNCTION_HD + int32 addPointsInBoxToList( + const teamMemberType& teamMember, + IdType id, + const iBox& triBox, + const PairsContainer& pairs)const + { + int32 getFull = 0; + + auto bExtent = boxExtent(triBox); + int32 numCellBox = bExtent.x()*bExtent.y()*bExtent.z(); + + const auto head = head_; + const auto next = next_; + + // perform a loop over all cells in the triBox + Kokkos::parallel_reduce( + Kokkos::TeamThreadRange(teamMember,numCellBox), + [=](int32 linIndex, int32& valToUpdate){ + + CellType cell; + indexToCell(linIndex, triBox, cell); + + int32 n = head(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(n), id) < 0 ) + valToUpdate++; + n = next(n); + } + + }, + getFull + ); + + return getFull; + + } + +}; + +} + +#endif diff --git a/src/Interaction/contactSearch/methods/NBSLoop.H b/src/Interaction/contactSearch/methods/NBSLoop.H new file mode 100644 index 00000000..6033ffa5 --- /dev/null +++ b/src/Interaction/contactSearch/methods/NBSLoop.H @@ -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. + +-----------------------------------------------------------------------------*/ + + +int32 m = head_(i,j,k); +CellType currentCell(i,j,k); +int32 n = -1; + +while( m > -1 ) +{ + + auto p_m = pointPosition_[m]; + auto d_m = sizeRatio_* diameter_[m]; + + // the same cell + n = next_(m); + + while(n >-1) + { + + auto p_n = 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 = 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 = head_(neighborCell.x(), neighborCell.y(), neighborCell.z()); + while( n>-1) + { + + auto p_n = 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 = next_[n]; + } + } + + } + m = next_[m]; +} + diff --git a/src/Interaction/contactSearch/wallMappings/cellsSimple.H b/src/Interaction/contactSearch/wallMappings/cellsSimple.H new file mode 100644 index 00000000..d0436a80 --- /dev/null +++ b/src/Interaction/contactSearch/wallMappings/cellsSimple.H @@ -0,0 +1,338 @@ +/*------------------------------- 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 __cellsSimple_H__ +#define __cellsSimple_H__ + +#include "types.H" +#include "KokkosTypes.H" +#include "cells.H" +#include "iBox.H" +#include "dictionary.H" + + +namespace pFlow +{ + +template< + typename executionSpace, + typename idType, + typename indexType = int32 + > +class cellsSimple +: + public cells +{ +public: + + using IdType = idType; + + using IndexType = indexType; + + using Cells = cells; + + using CellType = typename Cells::CellType; + + using ExecutionSpace = executionSpace; + + using memory_space = typename ExecutionSpace::memory_space; + + using iBoxType = iBox; + + bool constexpr static LOOP_ELEMENT_RANGE = true; + + class TagFindCellRange2{}; +protected: + + // - box extent + real cellExtent_ = 0.6; + + // - update frequency + int32 updateFrequency_=1; + + + int32 currentIter_ = 0; + + // - number of triangle elements + int32 numElements_ = 0; + + // - number of points + int32 numPoints_ = 0; + + /// a broad search has been occured during last pass? + bool performedSearch_ = false; + + // - ref to vectices + ViewType1D vertices_; + + // - ref to points in the trisurface + ViewType1D points_; + + // cell range of element/triangle bounding box + ViewType1D elementBox_; + + + using tpPWContactSearch = Kokkos::TeamPolicy< + ExecutionSpace, + Kokkos::Schedule, + Kokkos::IndexType + >; + + using rpFindCellRange2Type = + Kokkos::RangePolicy>; + + + FUNCTION_H + void allocateArrays() + { + Kokkos::realloc( elementBox_, numElements_); + } + +private: + + bool performSearch() + { + if(currentIter_ % updateFrequency_ == 0) + { + currentIter_++; + return true; + + }else + { + currentIter_++; + return false; + } + } + +public: + + TypeNameNV("cellsSimple"); + + FUNCTION_H + cellsSimple( + Cells ppCells, + real cellExtent, + int32 numPoints, + int32 numElements, + const ViewType1D& points, + const ViewType1D& vertices + ) + : + Cells(ppCells), + cellExtent_( max(cellExtent, 0.6 ) ), + numElements_(numElements), + numPoints_(numPoints), + vertices_(vertices), + points_(points) + { + + allocateArrays(); + } + + cellsSimple( + dictionary dict, + Cells ppCells, + int32 numPoints, + int32 numElements, + const ViewType1D& points, + const ViewType1D& vertices + ) + : + Cells(ppCells), + numElements_(numElements), + numPoints_(numPoints), + vertices_(vertices), + points_(points) + { + + updateFrequency_ = dict.getVal( + "updateFrequency" ); + updateFrequency_ = max(updateFrequency_,1); + + cellExtent_ = dict.getVal( + "cellExtent"); + + cellExtent_ = max(cellExtent_,0.6); + + allocateArrays(); + + } + + constexpr bool loopElementRange()const + { + return LOOP_ELEMENT_RANGE; + } + + // - host call + // reset triangle elements if they have changed + FUNCTION_H + bool resetElements( + int32 numElements, + int32 numPoints, + ViewType1D& points, + ViewType1D& 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_; + } + + bool enterBoadSearch()const + { + return currentIter_%updateFrequency_==0; + } + + bool performedSearch()const + { + return performedSearch_; + } + + template + bool broadSearch(PairsContainer& pairs, particleMapType& particleMap, bool force=false) + { + if(force) currentIter_ = 0; + performedSearch_= false; + if(!performSearch())return true; + + // map walls onto the cells + this->build(); + + this->particleWallFindPairs(pairs, particleMap); + + performedSearch_ = true; + return true; + } + + bool build() + { + Kokkos::parallel_for( + "cellsSimple::findcellrange2", + rpFindCellRange2Type(0,numElements_), + *this); + Kokkos::fence(); + return true; + } + + template + bool particleWallFindPairs(PairsContainer& pairs, particleMapType& particleMap) + { + + int32 getFull = 1; + + while (getFull) + { + + getFull = findPairsElementRange(pairs, particleMap); + + if(getFull) + { + // - resize the container + // note that getFull now shows the number of failed insertions. + uint32 len = max(getFull, 50); + auto oldCap = pairs.capacity(); + pairs.increaseCapacityBy(len); + + Info<<"Contact pair container capacity increased from "<< + oldCap << " to " + << pairs.capacity() <<" in cellsSimple."< + int32 findPairsElementRange(PairsContainer& pairs, particleMapType& particleMap) + { + int32 getFull =0; + + const auto pwPairs = pairs; + const auto elementBox = elementBox_; + + + Kokkos::parallel_reduce( + "cellsSimple::findPairsElementRange", + tpPWContactSearch(numElements_, Kokkos::AUTO), + LAMBDA_HD( + const typename tpPWContactSearch::member_type & teamMember, + int32& valueToUpdate){ + int32 i = teamMember.league_rank(); + IdType id = i; + const auto triBox = elementBox[i]; + + valueToUpdate += + particleMap.addPointsInBoxToList( + teamMember, + id, + triBox, + pwPairs + ); + }, + 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); + //output<< minP << " maxP "<< maxP<pointIndex(minP), this->pointIndex(maxP)); + + } + +}; + +} + + +#endif diff --git a/src/Interaction/interaction/demInteraction.H b/src/Interaction/interaction/demInteraction.H new file mode 100644 index 00000000..b6b3efaf --- /dev/null +++ b/src/Interaction/interaction/demInteraction.H @@ -0,0 +1,59 @@ +/*------------------------------- 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 __demInteraction_H__ +#define __demInteraction_H__ + +#include "property.H" +#include "demComponent.H" +#include "pointFields.H" +#include "triSurfaceFields.H" + +namespace pFlow +{ + +class demInteraction +: + public property, + public demComponent +{ +protected: + +public: + demInteraction(systemControl& control) + : + property(), + demComponent("interaction", control) + {} + + demInteraction(systemControl& control, const fileSystem& file) + : + property(file), + demComponent("interaction", control) + {} + + + +}; + + +} + +#endif //__interaction_H__ diff --git a/src/Interaction/interaction/interaction.C b/src/Interaction/interaction/interaction.C new file mode 100644 index 00000000..4675cd28 --- /dev/null +++ b/src/Interaction/interaction/interaction.C @@ -0,0 +1,107 @@ +/*------------------------------- 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 "interaction.H" + + +pFlow::interaction::interaction +( + systemControl& control, + const particles& prtcl, + const geometry& geom +) +: + demInteraction(control, control.caseSetup().path()+interactionFile__), + interactionBase(prtcl, geom), + fileDict_(control.caseSetup().emplaceObject( + objectFile( + interactionFile__, + "", + objectFile::READ_ALWAYS, + objectFile::WRITE_NEVER), + interactionFile__, + true )) +{ + this->subscribe(prtcl.pStruct()); + + contactSearch_ = contactSearch::create( + fileDict_.subDict("contactSearch"), + this->control().domain(), + prtcl, + geom, + timers() + ); + +} + + +pFlow::uniquePtr pFlow::interaction::create +( + systemControl& control, + const particles& prtcl, + const geometry& geom +) +{ + word shapeTypeName = prtcl.shapeTypeName(); + word motionTypeName = geom.motionModelTypeName(); + + fileSystem file = control.caseSetup().path()+interactionFile__; + + dictionary dict(interactionFile__, file); + auto interactionDict= dict.subDict("model"); + + word clType = dict.getVal("contactListType"); + word cfModel = interactionDict.getVal("contactForceModel"); + word rfModel = interactionDict.getVal("rollingFrictionModel"); + + + auto interactionModel = angleBracketsNames3( + shapeTypeName+"Interaction", + angleBracketsNames(rfModel,cfModel), + motionTypeName, + clType); + + + Report(1)<< "Selecting interaction model..."< contactSearch_ = nullptr; + +public: + + TypeName("interaction"); + + + //// - constructors + + interaction( + systemControl& control, + const particles& prtcl, + const geometry& geom ); + + + + virtual ~interaction() = default; + + create_vCtor( + interaction, + systemControl, + ( + systemControl& control, + const particles& prtcl, + const geometry& geom + ), + (control, prtcl, geom) + ); + + auto& contactSearchPtr() + { + return contactSearch_; + } + + auto& contactSearchRef() + { + return contactSearch_(); + } + + const auto& fileDict()const + { + return fileDict_; + } + + static + uniquePtr create( + systemControl& control, + const particles& prtcl, + const geometry& geom); + +}; + + +} + +#endif //__interaction_H__ diff --git a/src/Interaction/interaction/interactionBase.H b/src/Interaction/interaction/interactionBase.H new file mode 100644 index 00000000..dce5ece9 --- /dev/null +++ b/src/Interaction/interaction/interactionBase.H @@ -0,0 +1,84 @@ +/*------------------------------- 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 __interactionBase_H__ +#define __interactionBase_H__ + +#include "interactionTypes.H" +#include "particles.H" +#include "geometry.H" + +namespace pFlow +{ + +class interactionBase +{ +public: + + using IndexType = CELL_INDEX_TYPE; + + using IdType = ID_TYPE; + + using ExecutionSpace = DefaultExecutionSpace; + +protected: + + const particles& particles_; + + const geometry& geometry_; + +public: + + interactionBase( + const particles& prtcl, + const geometry& geom) + : + particles_(prtcl), + geometry_(geom) + {} + + inline + const auto& pStruct()const + { + return particles_.pStruct(); + } + + inline + const auto& surface()const + { + return geometry_.surface(); + } + + inline + const auto& Particles()const + { + return particles_; + } + + inline auto& Geometry()const + { + return geometry_; + } +}; + +} + + +#endif //__interactionBase_H__ diff --git a/src/Interaction/interaction/interactionTypes.H b/src/Interaction/interaction/interactionTypes.H new file mode 100644 index 00000000..a907018d --- /dev/null +++ b/src/Interaction/interaction/interactionTypes.H @@ -0,0 +1,45 @@ +/*------------------------------- 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 __interactionTypes_H__ +#define __interactionTypes_H__ + + + +#include "types.H" + + +namespace pFlow +{ + + +// always use signed types +using CELL_INDEX_TYPE = int32; + +using ID_TYPE = int32; + +//constexpr int32 minCellIndex = largestNegative(); + +//constexpr int32 maxCellIndex = largestPositive(); + +} + + +#endif //__interactionTypes_H__ diff --git a/src/Interaction/sphereInteraction/pLine.H b/src/Interaction/sphereInteraction/pLine.H new file mode 100644 index 00000000..2b4ad9aa --- /dev/null +++ b/src/Interaction/sphereInteraction/pLine.H @@ -0,0 +1,107 @@ +/*------------------------------- 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 __pLine_H__ +#define __pLine_H__ + +#include "types.H" + +namespace pFlow::sphTriInteraction +{ + +struct pLine +{ + + realx3 p1_; // point 1 + realx3 p2_; // piont 2 + realx3 v_; // direction vector + real L_; // line lenght + + INLINE_FUNCTION_HD + pLine(){} + + INLINE_FUNCTION_HD + pLine(const realx3 &p1, const realx3 &p2) + : + p1_(p1), + p2_(p2), + v_(p2-p1), + L_(length(v_)) + {} + + // get a point on the line based on input 0<= t <= 1 + INLINE_FUNCTION_HD + realx3 point(real t)const { + return v_ * t + p1_; + } + + // return the projected point of point p on line + INLINE_FUNCTION_HD + realx3 projectPoint(const realx3 &p) const + { + return point(projectNormLength(p)); + } + + // calculates the normalized distance between projected p and p1 + INLINE_FUNCTION_HD + real projectNormLength(realx3 p) const + { + realx3 w = p - p1_; + return dot(w,v_) / (L_*L_); + } + + INLINE_FUNCTION_HD + bool lineSphereCheck( + const realx3 pos, + real Rad, + realx3 &nv, + realx3 &cp, + real &ovrlp)const + { + + + real t = projectNormLength(pos); + + if(t >= 0.0 && t <= 1.0) cp = point(t); + else if(t >= (-Rad / L_) && t < 0.0) cp = point(0.0); + else if(t>1.0 && t >= (1.0 + Rad / L_)) cp = point(1.0); + else return false; + + realx3 vec = pos - cp; // from cp to pos + + real dist = length(vec); + ovrlp = Rad - dist; + + if (ovrlp >= 0.0) + { + if (dist > 0) + nv = vec / dist; + else + nv = v_; + return true; + } + + return false; + } +}; + +} //pFlow::sphTriInteractio + +#endif diff --git a/src/Interaction/sphereInteraction/sphereInteraction.C b/src/Interaction/sphereInteraction/sphereInteraction.C new file mode 100644 index 00000000..eee24bda --- /dev/null +++ b/src/Interaction/sphereInteraction/sphereInteraction.C @@ -0,0 +1,128 @@ +/*------------------------------- 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 contactListType > +bool pFlow::sphereInteraction:: + createSphereInteraction() +{ + + realVector_D rhoD(this->densities()); + + auto modelDict = this->fileDict().subDict("model"); + + Report(1)<<"Createing contact force model . . ."<( + this->numMaterials(), + rhoD.deviceVector(), + modelDict ); + + + uint32 nPrtcl = sphParticles_.size(); + + ppContactList_ = makeUnique(nPrtcl+1); + + pwContactList_ = makeUnique(nPrtcl/4+1); + + return true; +} + + + +template< + typename contactForceModel, + typename geometryMotionModel, + template class contactListType > +bool pFlow::sphereInteraction:: + 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 contactListType > +bool pFlow::sphereInteraction:: + sphereWallInteraction() +{ + + int32 lastItem = pwContactList_().loopCount(); + + pFlow::sphereInteractionKernels::pwInteractionFunctor + pwInteraction( + this->dt(), + this->forceModel_(), + pwContactList_(), + geometryMotion_.getTriangleAccessor(), + geometryMotion_.getModel() , + 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; +} \ No newline at end of file diff --git a/src/Interaction/sphereInteraction/sphereInteraction.H b/src/Interaction/sphereInteraction/sphereInteraction.H new file mode 100644 index 00000000..89ca8cef --- /dev/null +++ b/src/Interaction/sphereInteraction/sphereInteraction.H @@ -0,0 +1,207 @@ +/*------------------------------- 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_H__ +#define __sphereInteraction_H__ + +#include "interaction.H" +#include "sphereParticles.H" + +#include "sphereInteractionKernels.H" + +namespace pFlow +{ + +template< + typename contactForceModel, + typename geometryMotionModel, + template 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; + + using PairsContainerType= typename contactSearch::PairContainerType; + +protected: + + /// const reference to geometry + const GeometryMotionModel& geometryMotion_; + + /// const reference to particles + const sphereParticles& sphParticles_; + + + /// contact force model + uniquePtr forceModel_ = nullptr; + + /// contact list for particle-particle interactoins (keeps the history) + uniquePtr ppContactList_ = nullptr; + + /// contact list for particle-wall interactions (keeps the history) + uniquePtr 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::Schedule>; + + /// range policy for p-w interaction execution + using rpPWInteraction = rpPPInteraction; + +public: + + TypeNameTemplate3("sphereInteraction", ContactForceModel, MotionModel, ContactListType); + + // constructor + + sphereInteraction( + systemControl& control, + const particles& prtcl, + const geometry& geom) + : + interaction(control, prtcl, geom), + geometryMotion_(dynamic_cast(geom)), + sphParticles_(dynamic_cast(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 + { + + if(this->contactSearch_) + { + if( this->contactSearch_().ppEnterBroadSearch()) + { + ppContactList_().beforeBroadSearch(); + } + + if(this->contactSearch_().pwEnterBroadSearch()) + { + pwContactList_().beforeBroadSearch(); + } + + if( !contactSearch_().broadSearch( + ppContactList_(), + pwContactList_()) ) + { + fatalErrorInFunction<< + "unable to perform broadSearch.\n"; + fatalExit; + } + + if(this->contactSearch_().ppPerformedBroadSearch()) + { + ppContactList_().afterBroadSearch(); + } + + if(this->contactSearch_().pwPerformedBroadSearch()) + { + pwContactList_().afterBroadSearch(); + } + } + + if( sphParticles_.numActive()<=0)return true; + + ppInteractionTimer_.start(); + sphereSphereInteraction(); + ppInteractionTimer_.end(); + + pwInteractionTimer_.start(); + sphereWallInteraction(); + pwInteractionTimer_.end(); + + 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.C" + +#endif //__sphereInteraction_H__ diff --git a/src/Interaction/sphereInteraction/sphereInteractionKernels.H b/src/Interaction/sphereInteraction/sphereInteractionKernels.H new file mode 100644 index 00000000..d6fdc553 --- /dev/null +++ b/src/Interaction/sphereInteraction/sphereInteractionKernels.H @@ -0,0 +1,325 @@ +/*------------------------------- 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 __sphereInteractionKernels_H__ +#define __sphereInteractionKernels_H__ + + +#include "sphereTriSurfaceContact.H" + +namespace pFlow::sphereInteractionKernels +{ + +template< + typename ContactForceModel, + typename ContactListType> +struct ppInteractionFunctor +{ + + using PairType = typename ContactListType::PairType; + using ValueType = typename ContactListType::ValueType; + + real dt_; + + ContactForceModel forceModel_; + ContactListType tobeFilled_; + + deviceViewType1D diam_; + deviceViewType1D propId_; + deviceViewType1D pos_; + deviceViewType1D lVel_; + deviceViewType1D rVel_; + deviceViewType1D cForce_; + deviceViewType1D cTorque_; + + + ppInteractionFunctor( + real dt, + ContactForceModel forceModel, + ContactListType tobeFilled, + deviceViewType1D diam, + deviceViewType1D propId, + deviceViewType1D pos, + deviceViewType1D lVel, + deviceViewType1D rVel, + deviceViewType1D cForce, + deviceViewType1D cTorque ) + : + dt_(dt), + forceModel_(forceModel), + tobeFilled_(tobeFilled), + diam_(diam), + propId_(propId), + pos_(pos), + lVel_(lVel), + rVel_(rVel), + cForce_(cForce), // this is converted to an atomic vector + cTorque_(cTorque) // this is converted to an atomic vector + {} + + INLINE_FUNCTION_HD + void operator()(const int32 n)const + { + + if(!tobeFilled_.isValid(n))return; + + auto [i,j] = tobeFilled_.getPair(n); + + real Ri = 0.5*diam_[i]; + real Rj = 0.5*diam_[j]; + realx3 xi = pos_[i]; + realx3 xj = pos_[j]; + real dist = length(xj-xi); + real ovrlp = (Ri+Rj) - dist; + + if( ovrlp >0.0 ) + { + + auto Vi = lVel_[i]; + auto Vj = lVel_[j]; + auto wi = rVel_[i]; + auto wj = rVel_[j]; + auto Nij = (xj-xi)/dist; + auto Vr = Vi - Vj + cross((Ri*wi+Rj*wj), Nij); + + auto history = tobeFilled_.getValue(n); + + int32 propId_i = propId_[i]; + int32 propId_j = propId_[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_[i].x_,FC.x_); + Kokkos::atomic_add(&cForce_[i].y_,FC.y_); + Kokkos::atomic_add(&cForce_[i].z_,FC.z_); + + Kokkos::atomic_add(&cForce_[j].x_,-FC.x_); + Kokkos::atomic_add(&cForce_[j].y_,-FC.y_); + Kokkos::atomic_add(&cForce_[j].z_,-FC.z_); + + Kokkos::atomic_add(&cTorque_[i].x_, Mij.x_); + Kokkos::atomic_add(&cTorque_[i].y_, Mij.y_); + Kokkos::atomic_add(&cTorque_[i].z_, Mij.z_); + + Kokkos::atomic_add(&cTorque_[j].x_, Mji.x_); + Kokkos::atomic_add(&cTorque_[j].y_, Mji.y_); + Kokkos::atomic_add(&cTorque_[j].z_, Mji.z_); + + + tobeFilled_.setValue(n,history); + + } + else + { + tobeFilled_.setValue(n, ValueType()); + } + + } +}; + + +template< + typename ContactForceModel, + typename ContactListType, + typename TraingleAccessor, + typename MotionModel> +struct pwInteractionFunctor +{ + using PairType = typename ContactListType::PairType; + using ValueType = typename ContactListType::ValueType; + + real dt_; + + ContactForceModel forceModel_; + ContactListType tobeFilled_; + + TraingleAccessor triangles_; + MotionModel motionModel_; + + deviceViewType1D diam_; + deviceViewType1D propId_; + deviceViewType1D pos_; + deviceViewType1D lVel_; + deviceViewType1D rVel_; + deviceViewType1D cForce_; + deviceViewType1D cTorque_; + deviceViewType1D wTriMotionIndex_; + deviceViewType1D wPropId_; + deviceViewType1D wCForce_; + + + pwInteractionFunctor( + real dt, + ContactForceModel forceModel, + ContactListType tobeFilled, + TraingleAccessor triangles, + MotionModel motionModel , + deviceViewType1D diam , + deviceViewType1D propId, + deviceViewType1D pos , + deviceViewType1D lVel, + deviceViewType1D rVel, + deviceViewType1D cForce, + deviceViewType1D cTorque , + deviceViewType1D wTriMotionIndex, + deviceViewType1D wPropId, + deviceViewType1D wCForce) + : + dt_(dt), + forceModel_(forceModel), + tobeFilled_(tobeFilled), + triangles_(triangles) , + motionModel_(motionModel) , + diam_(diam) , + propId_(propId), + pos_(pos) , + lVel_(lVel), + rVel_(rVel) , + cForce_(cForce), + cTorque_(cTorque) , + wTriMotionIndex_(wTriMotionIndex) , + wPropId_(wPropId), + wCForce_(wCForce) + {} + + INLINE_FUNCTION_HD + void operator()(const int32 n)const + { + + if(!tobeFilled_.isValid(n))return; + + auto [i,tj] = tobeFilled_.getPair(n); + + real Ri = 0.5*diam_[i]; + real Rj = 10000.0; + realx3 xi = pos_[i]; + + realx3x3 tri = triangles_(tj); + real ovrlp; + realx3 Nij, cp; + + if( pFlow::sphTriInteraction::isSphereInContactBothSides( + tri, xi, Ri, ovrlp, Nij, cp) ) + { + auto Vi = lVel_[i]; + auto wi = rVel_[i]; + + int32 mInd = wTriMotionIndex_[tj]; + + auto Vw = motionModel_(mInd, cp); + + //output<< "par-wall index "<< i<<" - "<< tj<; + +template class pFlow::sphereInteraction< + pFlow::cfModels::limitedLinearNormalRolling, + pFlow::fixedGeometry, + pFlow::sortedContactList>; + +template class pFlow::sphereInteraction< + pFlow::cfModels::nonLimitedLinearNormalRolling, + pFlow::fixedGeometry, + pFlow::unsortedContactList>; + +template class pFlow::sphereInteraction< + pFlow::cfModels::nonLimitedLinearNormalRolling, + pFlow::fixedGeometry, + pFlow::sortedContactList>; + +template class pFlow::sphereInteraction< + pFlow::cfModels::limitedLinearNormalRolling, + pFlow::rotationAxisMotionGeometry, + pFlow::unsortedContactList>; + +template class pFlow::sphereInteraction< + pFlow::cfModels::limitedLinearNormalRolling, + pFlow::rotationAxisMotionGeometry, + pFlow::sortedContactList>; + +template class pFlow::sphereInteraction< + pFlow::cfModels::nonLimitedLinearNormalRolling, + pFlow::rotationAxisMotionGeometry, + pFlow::unsortedContactList>; + +template class pFlow::sphereInteraction< + pFlow::cfModels::nonLimitedLinearNormalRolling, + pFlow::rotationAxisMotionGeometry, + pFlow::sortedContactList>; + + + +/// non-linear models + +template class pFlow::sphereInteraction< + pFlow::cfModels::limitedNonLinearNormalRolling, + pFlow::fixedGeometry, + pFlow::unsortedContactList>; + +template class pFlow::sphereInteraction< + pFlow::cfModels::limitedNonLinearNormalRolling, + pFlow::fixedGeometry, + pFlow::sortedContactList>; + +template class pFlow::sphereInteraction< + pFlow::cfModels::nonLimitedNonLinearNormalRolling, + pFlow::fixedGeometry, + pFlow::unsortedContactList>; + +template class pFlow::sphereInteraction< + pFlow::cfModels::nonLimitedNonLinearNormalRolling, + pFlow::fixedGeometry, + pFlow::sortedContactList>; + +template class pFlow::sphereInteraction< + pFlow::cfModels::limitedNonLinearNormalRolling, + pFlow::rotationAxisMotionGeometry, + pFlow::unsortedContactList>; + +template class pFlow::sphereInteraction< + pFlow::cfModels::limitedNonLinearNormalRolling, + pFlow::rotationAxisMotionGeometry, + pFlow::sortedContactList>; + +template class pFlow::sphereInteraction< + pFlow::cfModels::nonLimitedNonLinearNormalRolling, + pFlow::rotationAxisMotionGeometry, + pFlow::unsortedContactList>; + +template class pFlow::sphereInteraction< + pFlow::cfModels::nonLimitedNonLinearNormalRolling, + pFlow::rotationAxisMotionGeometry, + pFlow::sortedContactList>; + + + + + diff --git a/src/Interaction/sphereInteraction/sphereTriSurfaceContact.H b/src/Interaction/sphereInteraction/sphereTriSurfaceContact.H new file mode 100644 index 00000000..7a5b611a --- /dev/null +++ b/src/Interaction/sphereInteraction/sphereTriSurfaceContact.H @@ -0,0 +1,231 @@ +/*------------------------------- 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 __sphereTriSurfaceContact_H__ +#define __sphereTriSurfaceContact_H__ + +#include "triWall.H" +#include "pLine.H" + +namespace pFlow::sphTriInteraction +{ + +INLINE_FUNCTION_HD +bool pointInPlane( + const realx3& p1, + const realx3& p2, + const realx3& p3, + const realx3& p ) +{ + + realx3 p1p = p1 - p; + realx3 p2p = p2 - p; + realx3 p3p = p3 - p; + + real p1p2 = dot(p1p, p2p); + real p2p3 = dot(p2p, p3p); + real p2p2 = dot(p2p, p2p); + real p1p3 = dot(p1p, p3p); + + // first condition u.v < 0 + // u.v = [(p1-p)x(p2-p)].[(p2-p)x(p3-p)] = (p1p.p2p)(p2p.p3p) - (p2p.p2p)(p1p.p3p) + if (p1p2*p2p3 - p2p2*p1p3 < 0.0) return false; + + + // second condition u.w < 0 + // u.w = [(p1-p)x(p2-p)].[(p3-p)x(p1-p)] = (p1p.p3p)(p2p.p1p) - (p2p.p3p)(p1p.p1p) + real p1p1 = dot(p1p, p1p); + if (p1p3*p1p2 - p2p3*p1p1 < (0.0)) return false; + + return true; +} + +INLINE_FUNCTION_HD +void cramerRule2(real A[2][2], real B[2], real & x1, real &x2) +{ + real det = (A[0][0] * A[1][1] - A[1][0]*A[0][1]); + x1 = (B[0]*A[1][1] - B[1]*A[0][1]) / det; + x2 = (A[0][0] * B[1] - A[1][0] * B[0])/ det; +} + +INLINE_FUNCTION_HD +bool pointInPlane( + const realx3& p1, + const realx3& p2, + const realx3& p3, + const realx3 &p, + int32 &Ln) +{ + realx3 v0 = p2 - p1; + realx3 v1 = p3 - p1; + realx3 v2 = p - p1; + + real A[2][2] = { dot(v0, v0), dot(v0, v1), dot(v1, v0), dot(v1, v1) }; + real B[2] = { dot(v0, v2), dot(v1, v2) }; + real nu, w; + + cramerRule2(A, B, nu, w); + real nuW = nu + w; + + + if (nuW > 1) + { + Ln = 2; return true; + } + else if (nuW >= 0) + { + if (nu >= 0 && w >= 0) + { + Ln = 0; return true; + } + if (nu > 0 && w < 0) + { + Ln = 1; return true; + } + if (nu < 0 && w > 0) + { + Ln = 3; return true; + } + } + else + { + Ln = 1; return true; + } + + return false; +} + +INLINE_FUNCTION_HD +bool isSphereInContactActiveSide( + const realx3& p1, + const realx3& p2, + const realx3& p3, + const realx3& cntr, + real rad, + real& ovrlp, + realx3& norm, + realx3& cp) +{ + + triWall wall(true, p1,p2,p3); + + real dist = wall.normalDistFromWall(cntr); + + if(dist < 0.0 )return false; + + ovrlp = rad - dist; + + if (ovrlp > 0) + { + realx3 ptOnPlane = wall.nearestPointOnWall(cntr); + + if (pointInPlane(p1, p2, p3, ptOnPlane)) + { + cp = ptOnPlane; + norm = -wall.n_; + return true; + } + + realx3 lnv; + + if (pLine(p1,p2).lineSphereCheck(cntr, rad, lnv, cp, ovrlp)) + { + norm = -lnv; + return true; + } + + if ( pLine(p2,p3).lineSphereCheck(cntr, rad, lnv, cp, ovrlp)) + { + norm = -lnv; + return true; + } + + if ( pLine(p3,p1).lineSphereCheck(cntr, rad, lnv, cp, ovrlp)) + { + norm = -lnv; + return true; + } + } + + return false; +} + +INLINE_FUNCTION_HD +bool isSphereInContactBothSides( + const realx3x3& tri, + const realx3& cntr, + real Rad, + real& ovrlp, + realx3& norm, + realx3& cp) +{ + + triWall wall(true, tri.x_,tri.y_,tri.z_); + + real dist = wall.normalDistFromWall(cntr); + + + ovrlp = Rad - abs(dist); + + if (ovrlp > 0) + { + realx3 ptOnPlane = wall.nearestPointOnWall(cntr); + + if (pointInPlane(tri.x_,tri.y_,tri.z_,ptOnPlane)) + { + cp = ptOnPlane; + + if(dist >= 0.0) + norm = -wall.n_; + else + norm = wall.n_; + return true; + } + + realx3 lnv; + + if (pLine(tri.x_, tri.y_).lineSphereCheck(cntr, Rad, lnv, cp, ovrlp)) + { + norm = -lnv; + return true; + } + + if ( pLine(tri.y_, tri.z_).lineSphereCheck(cntr, Rad, lnv, cp, ovrlp)) + { + norm = -lnv; + return true; + } + + if ( pLine(tri.z_, tri.x_).lineSphereCheck(cntr, Rad, lnv, cp, ovrlp)) + { + norm = -lnv; + return true; + } + } + + return false; +} + + +} // pFlow::sphTriInteraction + + +#endif //__sphereTriSurfaceContact_H__ \ No newline at end of file diff --git a/src/Interaction/sphereInteraction/triWall.H b/src/Interaction/sphereInteraction/triWall.H new file mode 100644 index 00000000..36ca12f9 --- /dev/null +++ b/src/Interaction/sphereInteraction/triWall.H @@ -0,0 +1,108 @@ +/*------------------------------- 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 __triWall_H__ +#define __triWall_H__ + +#include "types.H" + +namespace pFlow::sphTriInteraction +{ + +struct triWall +{ + + realx3 n_; + real offset_; + + INLINE_FUNCTION_H + triWall(const realx3& p1, const realx3& p2, const realx3& p3) + { + if(!makeWall(p1,p2,p3, n_, offset_)) + { + fatalErrorInFunction<< + "bad input for the wall.\n"; + fatalExit; + } + } + + INLINE_FUNCTION_HD + triWall(bool, const realx3& p1, const realx3& p2, const realx3& p3) + { + makeWall(p1,p2,p3,n_,offset_); + } + + INLINE_FUNCTION_HD + triWall(const realx3x3& tri) + { + makeWall(tri.x_, tri.y_, tri.z_, n_, offset_); + } + + INLINE_FUNCTION_HD + triWall(const triWall&) = default; + + INLINE_FUNCTION_HD + triWall& operator=(const triWall&) = default; + + INLINE_FUNCTION_HD + triWall(triWall&&) = default; + + INLINE_FUNCTION_HD + triWall& operator=(triWall&&) = default; + + INLINE_FUNCTION_HD + ~triWall()=default; + + + INLINE_FUNCTION_HD + real normalDistFromWall(const realx3 &p) const + { + return dot(n_, p) + offset_; + } + + INLINE_FUNCTION_HD + realx3 nearestPointOnWall(const realx3 &p) const + { + real t = -(dot(n_, p) + offset_); + return realx3(n_.x_*t + p.x_, n_.y_*t + p.y_, n_.z_*t + p.z_); + } + + INLINE_FUNCTION_HD static + bool makeWall( + const realx3& p1, + const realx3& p2, + const realx3& p3, + realx3& n, real& offset) + { + n = cross(p2 - p1, p3 - p1); + real len = length(n); + + if (len < 0.00000000001) return false; + n /= len; + offset = -dot(n, p1); + return true; + } + + +}; + +} + +#endif diff --git a/src/MotionModel/CMakeLists.txt b/src/MotionModel/CMakeLists.txt new file mode 100644 index 00000000..fdc69ecb --- /dev/null +++ b/src/MotionModel/CMakeLists.txt @@ -0,0 +1,13 @@ + +list(APPEND SourceFiles +entities/rotatingAxis/rotatingAxis.C +fixedWall/fixedWall.C +rotatingAxisMotion/rotatingAxisMotion.C +) + +set(link_libs Kokkos::kokkos phasicFlow) + +pFlow_add_library_install(MotionModel SourceFiles link_libs) + + + diff --git a/src/MotionModel/entities/rotatingAxis/rotatingAxis.C b/src/MotionModel/entities/rotatingAxis/rotatingAxis.C new file mode 100644 index 00000000..edf14ae5 --- /dev/null +++ b/src/MotionModel/entities/rotatingAxis/rotatingAxis.C @@ -0,0 +1,144 @@ +/*------------------------------- 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 "rotatingAxis.H" +#include "dictionary.H" + +/*FUNCTION_HD +pFlow::rotatingAxis::rotatingAxis() +: + line(), + omega_(0.0), + rotating_(false) +{}*/ + +FUNCTION_H +pFlow::rotatingAxis::rotatingAxis +( + const dictionary& dict +) +{ + + if(!read(dict)) + { + fatalErrorInFunction<< + " error in reading rotatingAxis from dictionary "<< dict.globalName()<(0.0)); + + setOmega(omega); + return true; +} + +FUNCTION_H +bool pFlow::rotatingAxis::write +( + dictionary& dict +) const +{ + if( !line::write(dict) ) return false; + + if( !dict.add("omega", omega_) ) + { + fatalErrorInFunction<< + " error in writing omega to dictionary "<< dict.globalName()<> omegaw >> omega; + if( !is.check(FUNCTION_NAME)) + { + ioErrorInFile(is.name(), is.lineNumber()); + return false; + } + if(omegaw != "omega") + { + ioErrorInFile(is.name(), is.lineNumber())<< + " expected omega but found "<< omegaw <>(iIstream& is, rotatingAxis& ax) +{ + if( !ax.read(is) ) + { + fatalExit; + } + return is; +} + +} + + +#include "rotatingAxisI.H" + +#endif diff --git a/src/MotionModel/entities/rotatingAxis/rotatingAxisFwd.H b/src/MotionModel/entities/rotatingAxis/rotatingAxisFwd.H new file mode 100755 index 00000000..8cbfe593 --- /dev/null +++ b/src/MotionModel/entities/rotatingAxis/rotatingAxisFwd.H @@ -0,0 +1,34 @@ +/*------------------------------- 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. + +-----------------------------------------------------------------------------*/ + +INLINE_FUNCTION_HD +realx3 rotate(const realx3 &p, const line& ln, real theta); + +INLINE_FUNCTION_HD +realx3 rotate(const realx3& p, const rotatingAxis& ax, real dt); + +INLINE_FUNCTION_HD +void rotate(realx3* p, size_t n, const line& ln, real theta); + +INLINE_FUNCTION_HD +void rotate(realx3* p, size_t n, const rotatingAxis& ax, real dt); + + + diff --git a/src/MotionModel/entities/rotatingAxis/rotatingAxisI.H b/src/MotionModel/entities/rotatingAxis/rotatingAxisI.H new file mode 100755 index 00000000..a46a634a --- /dev/null +++ b/src/MotionModel/entities/rotatingAxis/rotatingAxisI.H @@ -0,0 +1,164 @@ +/*------------------------------- 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. + +-----------------------------------------------------------------------------*/ + +INLINE_FUNCTION_HD +pFlow::realx3 pFlow::rotatingAxis::linTangentialVelocityPoint(const realx3 &p)const +{ + realx3 L = p - projectPoint(p); + return cross(omega_*unitVector(),L); +} + +INLINE_FUNCTION_HD +pFlow::realx3 pFlow::rotate(const realx3& p, const rotatingAxis& ax, real dt) +{ + realx3 nv = ax.unitVector(); + real cos_tet = cos(ax.omega()*dt); + real sin_tet = sin(ax.omega()*dt); + real u2 = nv.x()*nv.x(); + real v2 = nv.y()*nv.y(); + real w2 = nv.z()*nv.z(); + realx3 lp1 = ax.point1(); + + // (a(v2+w2) - u( bv + cw - ux - vy - wz)) (1-cos_tet) + x cos_tet + (- cv + bw - wy + vz) sin_tet + realx3 res; + + res.x_ = (lp1.x_*(v2 + w2) - (nv.x_*(lp1.y_*nv.y_ + lp1.z_*nv.z_ - nv.x_*p.x_ - nv.y_*p.y_ - nv.z_*p.z_)))*(1 - cos_tet) + + p.x_ * cos_tet + + (-lp1.z_*nv.y_ + lp1.y_*nv.z_ - nv.z_*p.y_ + nv.y_*p.z_) * sin_tet; + + + // ( b(u2+w2) - v( au + cw - ux - vy - wz))(1-cos_tet) + y cos_tet + ( cu - aw + wx - uz ) sin_tet + res.y_ = (lp1.y_*(u2 + w2) - (nv.y_*(lp1.x_*nv.x_ + lp1.z_*nv.z_ - nv.x_*p.x_ - nv.y_*p.y_ - nv.z_*p.z_)))*(1 - cos_tet) + + p.y_ * cos_tet + + (lp1.z_*nv.x_ - lp1.x_*nv.z_ + nv.z_*p.x_ - nv.x_*p.z_) * sin_tet; + + // (c(u2+v2) - w( au + bv - ux - vy - wz ))(1-cos_tet) + z cos_tet + (-bu + av - vx + uy) sin_tet + res.z_ = (lp1.z_*(u2 + v2) - (nv.z_*(lp1.x_*nv.x_ + lp1.y_*nv.y_ - nv.x_*p.x_ - nv.y_*p.y_ - nv.z_*p.z_)))*(1 - cos_tet) + + p.z_ * cos_tet + + (-lp1.y_*nv.x_ + lp1.x_*nv.y_ - nv.y_*p.x_ + nv.x_*p.y_) * sin_tet; + + return res; +} + +INLINE_FUNCTION_HD +pFlow::realx3 pFlow::rotate(const realx3 &p, const line& ln, real theta) +{ + + + realx3 nv = ln.unitVector(); + real cos_tet = cos(theta); + real sin_tet = sin(theta); + real u2 = nv.x()*nv.x(); + real v2 = nv.y()*nv.y(); + real w2 = nv.z()*nv.z(); + realx3 lp1 = ln.point1(); + + // (a(v2+w2) - u( bv + cw - ux - vy - wz)) (1-cos_tet) + x cos_tet + (- cv + bw - wy + vz) sin_tet + realx3 res; + + res.x_ = (lp1.x_*(v2 + w2) - (nv.x_*(lp1.y_*nv.y_ + lp1.z_*nv.z_ - nv.x_*p.x_ - nv.y_*p.y_ - nv.z_*p.z_)))*(1 - cos_tet) + + p.x_ * cos_tet + + (-lp1.z_*nv.y_ + lp1.y_*nv.z_ - nv.z_*p.y_ + nv.y_*p.z_) * sin_tet; + + + // ( b(u2+w2) - v( au + cw - ux - vy - wz))(1-cos_tet) + y cos_tet + ( cu - aw + wx - uz ) sin_tet + res.y_ = (lp1.y_*(u2 + w2) - (nv.y_*(lp1.x_*nv.x_ + lp1.z_*nv.z_ - nv.x_*p.x_ - nv.y_*p.y_ - nv.z_*p.z_)))*(1 - cos_tet) + + p.y_ * cos_tet + + (lp1.z_*nv.x_ - lp1.x_*nv.z_ + nv.z_*p.x_ - nv.x_*p.z_) * sin_tet; + + // (c(u2+v2) - w( au + bv - ux - vy - wz ))(1-cos_tet) + z cos_tet + (-bu + av - vx + uy) sin_tet + res.z_ = (lp1.z_*(u2 + v2) - (nv.z_*(lp1.x_*nv.x_ + lp1.y_*nv.y_ - nv.x_*p.x_ - nv.y_*p.y_ - nv.z_*p.z_)))*(1 - cos_tet) + + p.z_ * cos_tet + + (-lp1.y_*nv.x_ + lp1.x_*nv.y_ - nv.y_*p.x_ + nv.x_*p.y_) * sin_tet; + + return res; +} + +INLINE_FUNCTION_HD +void pFlow::rotate(realx3* p, size_t n, const line& ln, real theta) +{ + realx3 nv = ln.unitVector(); + real cos_tet = cos(theta); + real sin_tet = sin(theta); + real u2 = nv.x()*nv.x(); + real v2 = nv.y()*nv.y(); + real w2 = nv.z()*nv.z(); + realx3 lp1 = ln.point1(); + + // (a(v2+w2) - u( bv + cw - ux - vy - wz)) (1-cos_tet) + x cos_tet + (- cv + bw - wy + vz) sin_tet + realx3 res; + + for(label i=0; i("motionModel"); + + if(motionModel != "fixedWall") + { + fatalErrorInFunction<< + " motionModel should be fixedWall, but found " << motionModel <("motionModel"); + + if(motionModel != "rotatingAxisMotion") + { + fatalErrorInFunction<< + " motionModel should be rotatingAxisMotion, but found " + << motionModel <(axDict); axPtr) + { + axis_.push_back(axPtr()); + axisName_.push_back(aName); + } + else + { + fatalErrorInFunction<< + "could not read rotating axis from "<< axDict.globalName()< axis_; + int32 numAxis_=0; + + public: + + INLINE_FUNCTION_HD + Model(deviceViewType1D axis, int32 numAxis): + axis_(axis), + numAxis_(numAxis) + {} + + INLINE_FUNCTION_HD + Model(const Model&) = default; + + + INLINE_FUNCTION_HD + Model& operator=(const Model&) = default; + + + INLINE_FUNCTION_HD + realx3 pointVelocity(int32 n, const realx3& p)const + { + return axis_[n].linTangentialVelocityPoint(p); + } + + INLINE_FUNCTION_HD + realx3 operator()(int32 n, const realx3& p)const + { + return pointVelocity(n,p); + } + + INLINE_FUNCTION_HD + realx3 transferPoint(int32 n, const realx3 p, real dt)const + { + return rotate(p, axis_[n], dt); + } + + INLINE_FUNCTION_HD int32 numComponents()const + { + return numAxis_; + } + }; + +protected: + + using axisVector_HD = VectorDual; + + axisVector_HD axis_; + + wordList axisName_; + + label numAxis_= 0; + + bool readDictionary(const dictionary& dict); + + bool writeDictionary(dictionary& dict)const; + +public: + + TypeNameNV("rotatingAxisMotion"); + + // empty + FUNCTION_H + rotatingAxisMotion(); + + // construct with dictionary + FUNCTION_H + rotatingAxisMotion(const dictionary& dict); + + // copy + FUNCTION_H + rotatingAxisMotion(const rotatingAxisMotion&) = default; + + rotatingAxisMotion(rotatingAxisMotion&&) = delete; + + FUNCTION_H + rotatingAxisMotion& operator=(const rotatingAxisMotion&) = default; + + rotatingAxisMotion& operator=(rotatingAxisMotion&&) = delete; + + FUNCTION_H + ~rotatingAxisMotion() = default; + + + Model getModel()const + { + return Model(axis_.deviceVector(), numAxis_); + } + + INLINE_FUNCTION_H + int32 nameToIndex(const word& name)const + { + if( auto i = axisName_.findi(name); i == -1) + { + fatalErrorInFunction<< + "axis name " << name << " does not exist. \n"; + fatalExit; + return i; + } + else + { + return i; + } + + } + + INLINE_FUNCTION_H + word indexToName(label i)const + { + if(i < numAxis_ ) + return axisName_[i]; + else + { + fatalErrorInFunction<< + "out of range access to the list of axes " << i <=numAxis_)return false; + + rotate( pVec, numP, axis_.deviceVectorAll()[n], dt); + + return true; + } + + INLINE_FUNCTION_HD + bool isMoving()const + { + return true; + } + + INLINE_FUNCTION_H + bool move(real dt) + { + return true; + } + + FUNCTION_H + bool read(iIstream& is); + + FUNCTION_H + bool write(iOstream& os)const; + + +}; + +} // pFlow + +#endif //__rotatingAxisMotion_H__ diff --git a/src/Particles/CMakeLists.txt b/src/Particles/CMakeLists.txt new file mode 100644 index 00000000..725c425f --- /dev/null +++ b/src/Particles/CMakeLists.txt @@ -0,0 +1,19 @@ + +set(SourceFiles +dynamicPointStructure/dynamicPointStructure.C +particles/particles.C +SphereParticles/sphereShape/sphereShape.C +SphereParticles/sphereParticles/sphereParticles.C +Insertion/shapeMixture/shapeMixture.C +Insertion/insertion/insertion.C +Insertion/Insertion/Insertions.C +Insertion/insertionRegion/insertionRegion.C +Insertion/insertionRegion/timeFlowControl.C +) + +set(link_libs Kokkos::kokkos phasicFlow Integration Property) + +pFlow_add_library_install(Particles SourceFiles link_libs) + + + diff --git a/src/Particles/Insertion/Insertion/Insertion.C b/src/Particles/Insertion/Insertion/Insertion.C new file mode 100644 index 00000000..08a7f4b9 --- /dev/null +++ b/src/Particles/Insertion/Insertion/Insertion.C @@ -0,0 +1,204 @@ +/*------------------------------- 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 +bool pFlow::Insertion::readInsertionDict +( + const dictionary& dict +) +{ + if(!insertion::readInsertionDict(dict)) return false; + + regions_.clear(); + + if( !this->isActive() ) + { + return true; + } + + wordList regionDicNames = dict.dictionaryKeywords(); + + for(auto& name:regionDicNames) + { + Report(2)<<"reading insertion region "<< greenText(name)< +bool pFlow::Insertion::writeInsertionDict +( + dictionary& dict +)const +{ + if( !insertion::writeInsertionDict(dict) ) return false; + + if( !this->isActive() ) return true; + + forAll(i,regions_) + { + auto& rgnDict = dict.subDictOrCreate(regions_[i].name()); + + if( !regions_[i].write(rgnDict) ) + { + return false; + } + } + + return true; +} + +template +pFlow::Insertion::Insertion +( + particles& prtcl, + const ShapeType& shapes +) +: + insertion(prtcl), + shapes_(shapes) +{ + + +} + + +template +bool pFlow::Insertion::insertParticles +( + real currentTime, + real dt +) +{ + if(!isActive()) return true; + + + forAll(i,regions_) + { + bool insertionOccured = false; + auto& rgn = regions_[i]; + if( rgn.insertionTime(currentTime, dt) ) + { + + realx3Vector pos; + wordVector shapes; + if( rgn.insertParticles(currentTime, dt, shapes, pos, insertionOccured) ) + { + + if(insertionOccured) + { + Report(0)<<"\nParticle insertion from "<< greenText(rgn.name())< +bool pFlow::Insertion::read +( + iIstream& is +) +{ + + // create an empty dictionary + dictionary dict(is.name(), true); + + if(!dict.read(is)) + { + ioErrorInFile( is.name(), is.lineNumber() )<< + " error in reading "<< insertionFile__ << "dictionary from file."< +bool pFlow::Insertion::write +( + iOstream& os +)const +{ + + dictionary dict(insertionFile__,true); + + if(! writeInsertionDict(dict) ) + { + fatalErrorInFunction<< + " error in writing to " << dict.globalName()< +class Insertion +: + public insertion +{ +protected: + + const ShapeType& shapes_; + + // - insertion regions + ListPtr> regions_; + + + bool readInsertionDict(const dictionary& dict); + + bool writeInsertionDict(dictionary& dict)const; + +public: + + TypeNameTemplateNV("Insertion",ShapeType); + + Insertion(particles& prtcl, const ShapeType& shapes); + + + bool insertParticles(real currentTime, real dt); + + virtual bool read(iIstream& is) override; + + virtual bool write(iOstream& os)const override; + +}; + +} + +#include "Insertion.C" + +#endif diff --git a/src/Particles/Insertion/Insertion/Insertions.C b/src/Particles/Insertion/Insertion/Insertions.C new file mode 100644 index 00000000..5951a090 --- /dev/null +++ b/src/Particles/Insertion/Insertion/Insertions.C @@ -0,0 +1,24 @@ +/*------------------------------- 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 "Insertions.H" + +template class pFlow::Insertion; \ No newline at end of file diff --git a/src/Particles/Insertion/Insertion/Insertions.H b/src/Particles/Insertion/Insertion/Insertions.H new file mode 100644 index 00000000..299df6c9 --- /dev/null +++ b/src/Particles/Insertion/Insertion/Insertions.H @@ -0,0 +1,35 @@ +/*------------------------------- 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 __Insertions_H__ +#define __Insertions_H__ + + +#include "Insertion.H" +#include "sphereShape.H" + +namespace pFlow +{ + +using sphereInsertion = Insertion ; + +} + +#endif diff --git a/src/Particles/Insertion/InsertionRegion/InsertionRegion.C b/src/Particles/Insertion/InsertionRegion/InsertionRegion.C new file mode 100644 index 00000000..5af830c9 --- /dev/null +++ b/src/Particles/Insertion/InsertionRegion/InsertionRegion.C @@ -0,0 +1,115 @@ +/*------------------------------- 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 +bool pFlow::InsertionRegion::checkForContact +( + const realx3Vector& pos, + const realVector& diams, + const realx3& p, + const real& d +) +{ + + forAll(i, pos) + { + if( length(pos[i]-p) < 0.5*(diams[i]+d) ) return true; + } + + return false; +} + +template +pFlow::InsertionRegion::InsertionRegion +( + const dictionary& dict, + const ShapeType& shapes +) +: + insertionRegion(dict), + shapes_(shapes) +{ + +} + + +template +bool pFlow::InsertionRegion::insertParticles +( + real currentTime, + real dt, + wordVector& names, + realx3Vector& pos, + bool& insertionOccured +) +{ + insertionOccured = false; + + if(!insertionTime( currentTime, dt)) return true; + + size_t newNum = numberToBeInserted(currentTime); + + if(newNum == 0) return true; + + names.reserve(max(newNum,names.capacity())); + pos.reserve(max(newNum,pos.capacity())); + names.clear(); + pos.clear(); + + realVector diams(newNum, RESERVE()); + + mixture_->getNextShapeNameN(newNum, names); + + if(!shapes_.shapeToDiameter(names,diams)) + { + fatalErrorInFunction<< + " error occured in insertion region "<< name() << + " while converting shapes to diameter. \n"; + return false; + } + + size_t n = 0; + + for(label iter=0; iter< 10*newNum ; ++iter) + { + if( !(n < newNum) ) + { + addToNumInserted(newNum); + insertionOccured = true; + return true; + } + realx3 p = pRegion_().peek(); + real d = diams[pos.size()]; + if( !checkForContact(pos, diams, p, d) ) + { + pos.push_back(p); + n++; + } + + } + + fatalErrorInFunction<< + " Cannot insert "<< newNum << " new particles from region "<< name()<<". \n" + " pFlow could position only "<< n<< " particles in this region. \n"; + addToNumInserted(n); + insertionOccured = false; + return false; + +} \ No newline at end of file diff --git a/src/Particles/Insertion/InsertionRegion/InsertionRegion.H b/src/Particles/Insertion/InsertionRegion/InsertionRegion.H new file mode 100644 index 00000000..bae48646 --- /dev/null +++ b/src/Particles/Insertion/InsertionRegion/InsertionRegion.H @@ -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 __InsertionRegion_H__ +#define __InsertionRegion_H__ + + +#include "insertionRegion.H" +#include "dictionary.H" + +namespace pFlow +{ + +template +class InsertionRegion +: + public insertionRegion +{ +protected: + // - type of particle shapes + const ShapeType& shapes_; + + static bool checkForContact( + const realx3Vector& pos, + const realVector& diams, + const realx3& p, + const real& d); + +public: + + // - type info + TypeNameTemplateNV("insertionRegion", ShapeType); + + InsertionRegion(const dictionary& dict, const ShapeType& shapes); + + InsertionRegion(const InsertionRegion& ) = default; + + InsertionRegion(InsertionRegion&&) = default; + + InsertionRegion& operator=(const InsertionRegion& ) = default; + + InsertionRegion& operator=(InsertionRegion&&) = default; + + + auto clone()const + { + return makeUnique>(*this); + } + + auto clonePtr()const + { + return new InsertionRegion(*this); + } + + + bool insertParticles + ( + real currentTime, + real dt, + wordVector& names, + realx3Vector& pos, + bool& insertionOccured + ); + + //bool read(const dictionary& dict); + + //bool write(dictionary& dict)const; + +}; + + + +} // pFlow + + +#include "InsertionRegion.C" + +#endif diff --git a/src/Particles/Insertion/insertion/insertion.C b/src/Particles/Insertion/insertion/insertion.C new file mode 100644 index 00000000..b5c4e12d --- /dev/null +++ b/src/Particles/Insertion/insertion/insertion.C @@ -0,0 +1,74 @@ +/*------------------------------- 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 "particles.H" +#include "dictionary.H" +#include "insertion.H" +#include "streams.H" + +bool pFlow::insertion::readInsertionDict +( + const dictionary& dict +) +{ + + active_ = dict.getVal("active"); + + if(active_) + Report(1)<< "Particle insertion mechanism is "<< + yellowText("active")<<" in the simulation."<("type"); + + pRegion_ = peakableRegion::create(type_, dict.subDict(type_+"Info")); + + mixture_ = makeUnique(dict.subDict("mixture")); + + addToNumInserted(mixture_().totalInserted()); + + if( !dict.containsDictionay("setFields")) + { + output<<"\n insertion region "<< name_ << " does not contain setFields dictionary." + " An empty dictoiinary is created for it. \n"; + setFields_ = makeUnique( dictionary("setFields") ); + } + else + { + setFields_ = makeUnique( dict.subDict("setFields") ); + } + + for(auto& sfEntry:setFields_()) + { + if(!sfEntry.checkForTypeAndValueAll()) + { + fatalErrorInFunction<< + " error in setFields dictionary "<< dict.globalName()<clone(): nullptr), + mixture_( src.mixture_? src.mixture_->clone(): nullptr), + setFields_( src.setFields_? src.setFields_->clone(): nullptr) +{} + +pFlow::insertionRegion& pFlow::insertionRegion::operator= +( + const insertionRegion& src +) +{ + + if(&src == this)return *this; + timeFlowControl::operator=(src); + + name_ = src.name_; + type_ = src.type_; + pRegion_ = (src.pRegion_? src.pRegion_->clone(): nullptr); + mixture_ = (src.mixture_? src.mixture_->clone(): nullptr); + setFields_ = (src.setFields_? src.setFields_->clone(): nullptr); + + return *this; +} + + diff --git a/src/Particles/Insertion/insertionRegion/insertionRegion.H b/src/Particles/Insertion/insertionRegion/insertionRegion.H new file mode 100644 index 00000000..b20b3116 --- /dev/null +++ b/src/Particles/Insertion/insertionRegion/insertionRegion.H @@ -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 __insertionRegion_H__ +#define __insertionRegion_H__ + +#include "timeFlowControl.H" +#include "shapeMixture.H" +#include "peakableRegions.H" +#include "setFieldList.H" + +namespace pFlow +{ + +class dictionary; + +class insertionRegion +: + public timeFlowControl +{ +protected: + + // - name of the region + word name_; + + // - type of insertion region + word type_; + + // peakable region of points + uniquePtr pRegion_ = nullptr; + + // mixture of shapes + uniquePtr mixture_ = nullptr; + + // setFields for insertion region + uniquePtr setFields_ = nullptr; + + + bool readInsertionRegion(const dictionary& dict); + + bool writeInsertionRegion(dictionary& dict) const; + + +public: + + TypeNameNV("insertionRegion"); + + //// - Constructors + + insertionRegion(const dictionary& dict); + + insertionRegion(const insertionRegion& src); + + insertionRegion(insertionRegion&&) = default; + + insertionRegion& operator=(const insertionRegion&); + + insertionRegion& operator=(insertionRegion&&) = default; + + + ~insertionRegion() = default; + + + //// - Methods + const auto& setFields()const + { + return setFields_(); + } + + const auto& name()const + { + return name_; + } + + + //// - IO operation + + bool read(const dictionary& dict) + { + if(!timeFlowControl::read(dict))return false; + + return readInsertionRegion(dict); + } + + bool write(dictionary& dict)const + { + if(!timeFlowControl::write(dict)) return false; + + return writeInsertionRegion(dict); + } + + +}; + +} //pFlow + +#endif //__insertionRegion_H__ diff --git a/src/Particles/Insertion/insertionRegion/timeFlowControl.C b/src/Particles/Insertion/insertionRegion/timeFlowControl.C new file mode 100644 index 00000000..98ce99d6 --- /dev/null +++ b/src/Particles/Insertion/insertionRegion/timeFlowControl.C @@ -0,0 +1,62 @@ +/*------------------------------- 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 "timeFlowControl.H" +#include "dictionary.H" + +bool pFlow::timeFlowControl::readTimeFlowControl +( + const dictionary& dict +) +{ + rate_ = dict.getVal("rate"); + startTime_ = dict.getVal("startTime"); + endTime_ = dict.getVal("endTime"); + interval_ = dict.getVal("interval"); + numInserted_=0; + return true; +} + +bool pFlow::timeFlowControl::writeTimeFlowControl +( + dictionary& dict +) const +{ + if(!dict.add("rate", rate_)) return false; + if(!dict.add("startTime", startTime_)) return false; + if(!dict.add("endTime", endTime_)) return false; + if(!dict.add("interval", interval_)) return false; + + return true; +} + + +pFlow::timeFlowControl::timeFlowControl +( + const dictionary& dict +) +{ + + if(!readTimeFlowControl(dict)) + { + fatalExit; + } + +} \ No newline at end of file diff --git a/src/Particles/Insertion/insertionRegion/timeFlowControl.H b/src/Particles/Insertion/insertionRegion/timeFlowControl.H new file mode 100644 index 00000000..e29e150b --- /dev/null +++ b/src/Particles/Insertion/insertionRegion/timeFlowControl.H @@ -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 __timeFlowControl_H__ +#define __timeFlowControl_H__ + +#include "types.H" +#include "streams.H" + +namespace pFlow +{ + +class dictionary; + +class timeFlowControl +{ +protected: + + real startTime_; + + real endTime_; + + real interval_; + + real rate_; + + size_t numInserted_ = 0; + + bool readTimeFlowControl(const dictionary& dict); + + bool writeTimeFlowControl(dictionary& dict) const; + + size_t numberToBeInserted(real currentTime) + { + if(currentTimeendTime_) return 0; + + return static_cast( (currentTime - startTime_ + interval_)*rate_ - numInserted_ ); + } + + size_t addToNumInserted(size_t newInserted) + { + return numInserted_ += newInserted; + } + +public: + + + timeFlowControl(const dictionary& dict); + + + bool insertionTime( real currentTime, real dt) + { + if(currentTime < startTime_) return false; + + if(currentTime > endTime_) return false; + if( mod(abs(currentTime-startTime_),interval_)/dt < 1 ) return true; + + return false; + } + + size_t totalInserted()const + { + return numInserted_; + } + + bool read(const dictionary& dict) + { + return readTimeFlowControl(dict); + } + + bool write(dictionary& dict)const + { + return writeTimeFlowControl(dict); + } + +}; + +} + +#endif //__timeFlowControl_H__ diff --git a/src/Particles/Insertion/shapeMixture/shapeMixture.C b/src/Particles/Insertion/shapeMixture/shapeMixture.C new file mode 100755 index 00000000..1a79285a --- /dev/null +++ b/src/Particles/Insertion/shapeMixture/shapeMixture.C @@ -0,0 +1,151 @@ +/*------------------------------- 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 "shapeMixture.H" +#include "dictionary.H" + + +pFlow::shapeMixture::shapeMixture +( + const dictionary & dict +) +{ + if( !read(dict)) + { + fatalExit; + } +} + +pFlow::word pFlow::shapeMixture::getNextShapeName() +{ + + forAll(i, names_) + { + if(current_[i]< number_[i]) + { + current_[i]++; + numberInserted_[i]++; + return names_[i]; + } + } + + fill(current_, static_cast(0)); + return getNextShapeName(); +} + +void pFlow::shapeMixture::getNextShapeNameN +( + size_t n, + wordVector& names +) +{ + names.clear(); + + + for(label i=0; i(nm); + if( num <= 0 ) + { + fatalErrorInFunction<< + " number inserte in front of "<< nm << " is invalid: "<< num<("numberInserted"); + } + else + { + numberInserted_ = uint32Vector(size(), static_cast(0)); + } + + if(numberInserted_.size() != names_.size() ) + { + fatalErrorInFunction<< + " number of elements in numberInserted ("< clone()const + { + return makeUnique(*this); + } + + shapeMixture* clonePtr()const + { + return new shapeMixture(*this); + } + + + // + ~shapeMixture() = default; + + //// - Methods + word getNextShapeName(); + + + void getNextShapeNameN(size_t n, wordVector& names); + + + auto size()const { + return names_.size(); + } + + auto totalInserted()const { + return sum(numberInserted_); + } + + //// - IO operatoins + bool read(const dictionary& dict); + + bool write(dictionary& dict) const; + +}; + +} // pFlow + +#endif //__shapeMixture_H__ \ No newline at end of file diff --git a/src/Particles/SphereParticles/sphereParticles/sphereParticles.C b/src/Particles/SphereParticles/sphereParticles/sphereParticles.C new file mode 100644 index 00000000..e3c0ba48 --- /dev/null +++ b/src/Particles/SphereParticles/sphereParticles/sphereParticles.C @@ -0,0 +1,353 @@ +/*------------------------------- 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 "sphereParticles.H" +#include "setFieldList.H" +#include "sphereParticlesKernels.H" + +pFlow::uniquePtr> +pFlow::sphereParticles::getFieldObjectList()const +{ + auto objListPtr = particles::getFieldObjectList(); + + objListPtr().push_back( + static_cast(&I_) ); + + return objListPtr; +} + +bool pFlow::sphereParticles::diameterMassInertiaPropId +( + const word& shName, + real& diam, + real& mass, + real& I, + int8& propIdx +) +{ + uint32 idx; + if( !shapes_.nameToIndex(shName, idx) ) + { + printKeys(fatalErrorInFunction<< + " wrong shape name in the input: "<< shName<(pIdx); + return true; +} + +bool pFlow::sphereParticles::initializeParticles() +{ + + int32IndexContainer indices(0, shapeName_.size()); + return insertSphereParticles(shapeName_, indices); +} + + +bool pFlow::sphereParticles::beforeIteration() +{ + particles::beforeIteration(); + + intPredictTimer_.start(); + dynPointStruct_.predict(this->dt(), accelertion_); + rVelIntegration_().predict(this->dt(),rVelocity_, rAcceleration_); + intPredictTimer_.end(); + + return true; +} + + +bool pFlow::sphereParticles::iterate() +{ + accelerationTimer_.start(); + pFlow::sphereParticlesKernels::acceleration( + control().g(), + mass().deviceVectorAll(), + contactForce().deviceVectorAll(), + I().deviceVectorAll(), + contactTorque().deviceVectorAll(), + pStruct().activePointsMaskD(), + accelertion().deviceVectorAll(), + rAcceleration().deviceVectorAll() + ); + accelerationTimer_.end(); + + intCorrectTimer_.start(); + dynPointStruct_.correct(this->dt(), accelertion_); + rVelIntegration_().correct(this->dt(), rVelocity_, rAcceleration_); + intCorrectTimer_.end(); + + return true; +} + + +bool pFlow::sphereParticles::afterIteration() +{ + return true; +} + + +bool pFlow::sphereParticles::insertSphereParticles( + const wordVector& names, + const int32IndexContainer& indices + ) +{ + + if(names.size()!= indices.size()) + { + fatalErrorInFunction << + "sizes of names ("<( + objectFile( + sphereShapeFile__, + "", + objectFile::READ_ALWAYS, + objectFile::WRITE_NEVER + ) + ) + ), + I_( + this->time().emplaceObject( + objectFile( + "I", + "", + objectFile::READ_NEVER, + objectFile::WRITE_ALWAYS + ), + pStruct(), + static_cast(0.0000000001) + ) + ), + rVelocity_( + this->time().emplaceObject( + objectFile( + "rVelocity", + "", + objectFile::READ_IF_PRESENT, + objectFile::WRITE_ALWAYS + ), + pStruct(), + zero3 + ) + ), + rAcceleration_( + this->time().emplaceObject( + objectFile( + "rAcceleration", + "", + objectFile::READ_IF_PRESENT, + objectFile::WRITE_ALWAYS + ), + pStruct(), + zero3 + ) + ), + accelerationTimer_( + "Acceleration", &this->timers() ), + intPredictTimer_( + "Integration-predict", &this->timers() ), + intCorrectTimer_( + "Integration-correct", &this->timers() ) + +{ + + Report(1)<<"Creating integration method "<integrationMethod()) + << " for rotational velocity."<time().integration(), + this->pStruct(), + this->integrationMethod()); + + if( !rVelIntegration_ ) + { + fatalErrorInFunction<< + " error in creating integration object for rVelocity. \n"; + fatalExit; + } + + + if(!initializeParticles()) + { + fatalExit; + } + +} + + +bool pFlow::sphereParticles::insertParticles +( + const realx3Vector& position, + const wordVector& shapes, + const setFieldList& setField + ) +{ + + if( position.size() != shapes.size() ) + { + fatalErrorInFunction<< + " size of vectors position ("<activeRange(); + + Report(1)<< "Active range is "<numActive())<< + " and pointStructure capacity is "<capacity())< rVelIntegration_ = nullptr; + + /// timer for acceleration computations + Timer accelerationTimer_; + + /// timer for integration computations (prediction step) + Timer intPredictTimer_; + + /// timer for integration computations (correction step) + Timer intCorrectTimer_; + + bool diameterMassInertiaPropId(const word& shName, real& diam, real& mass, real& I, int8& propIdx); + + bool initializeParticles(); + + bool insertSphereParticles(const wordVector& names, const int32IndexContainer& indices); + + virtual uniquePtr> getFieldObjectList()const override; + +public: + + /// construct from systemControl and property + sphereParticles(systemControl &control, const property& prop); + + /// no copy constructor + sphereParticles(const sphereParticles&) = delete; + + /// move constructor + sphereParticles(sphereParticles&&) = default; + + /// no copy-assignement + sphereParticles& operator=(const sphereParticles&) = delete; + + /// move-assignement + sphereParticles& operator=(sphereParticles&&) = default; + + virtual ~sphereParticles()=default; + + /** + * Insert new particles in position with specified shapes + * + * This function is involked by inserted object to insert new set of particles + * into the simulation. + * \param position position of new particles + * \param shape shape of new particles + * \param setField initial value of the selected fields for new particles + */ + bool insertParticles + ( + const realx3Vector& position, + const wordVector& shapes, + const setFieldList& setField + ) override ; + + /// const reference to shapes object + const auto& shapes()const + { + return shapes_; + } + + /// const reference to inertia pointField + const auto& I()const + { + return I_; + } + + /// reference to inertia pointField + auto& I() + { + return I_; + } + + const realx3Vector_D rVelocity()const + { + return rVelocity_; + } + + const realVector_D& boundingSphere()const override + { + return this->diameter(); + } + + word shapeTypeName() const override + { + return "sphere"; + } + + void boundingSphereMinMax( + real& minDiam, + real& maxDiam )const override + { + shapes_.diameterMinMax(minDiam, maxDiam); + } + + bool update(const eventMessage& msg) override + { + CONSUME_PARAM(msg); + return true; + } + + realx3PointField_D& rAcceleration() override + { + return rAcceleration_; + } + + const realx3PointField_D& rAcceleration() const override + { + return rAcceleration_; + } + + /// before iteration step + bool beforeIteration() override; + + /// iterate particles + bool iterate() override; + + /// after iteration step + bool afterIteration() override; + + + +}; //sphereParticles + +} // pFlow + +#endif //__sphereParticles_H__ diff --git a/src/Particles/SphereParticles/sphereParticles/sphereParticlesKernels.H b/src/Particles/SphereParticles/sphereParticles/sphereParticlesKernels.H new file mode 100644 index 00000000..0bab1397 --- /dev/null +++ b/src/Particles/SphereParticles/sphereParticles/sphereParticlesKernels.H @@ -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. + +-----------------------------------------------------------------------------*/ + +#ifndef __sphereParticlesKernels_H__ +#define __sphereParticlesKernels_H__ + +namespace pFlow::sphereParticlesKernels +{ + +using rpAcceleration = Kokkos::RangePolicy< + DefaultExecutionSpace, + Kokkos::Schedule, + Kokkos::IndexType + >; + +template +void acceleration( + realx3 g, + deviceViewType1D mass, + deviceViewType1D force, + deviceViewType1D I, + deviceViewType1D torque, + IncludeFunctionType incld, + deviceViewType1D lAcc, + deviceViewType1D rAcc + ) +{ + + auto activeRange = incld.activeRange(); + if(incld.allActive()) + { + Kokkos::parallel_for( + "pFlow::sphereParticlesKernels::acceleration", + rpAcceleration(activeRange.first, activeRange.second), + LAMBDA_HD(int32 i){ + lAcc[i] = force[i]/mass[i] + g; + rAcc[i] = torque[i]/I[i]; + }); + } + else + { + Kokkos::parallel_for( + "pFlow::sphereParticlesKernels::acceleration", + rpAcceleration(activeRange.first, activeRange.second), + LAMBDA_HD(int32 i){ + if(incld(i)) + { + lAcc[i] = force[i]/mass[i] + g; + rAcc[i] = torque[i]/I[i]; + } + }); + + } + + Kokkos::fence(); +} + +} + +#endif diff --git a/src/Particles/SphereParticles/sphereShape/sphereShape.C b/src/Particles/SphereParticles/sphereShape/sphereShape.C new file mode 100644 index 00000000..833d4748 --- /dev/null +++ b/src/Particles/SphereParticles/sphereShape/sphereShape.C @@ -0,0 +1,204 @@ +/*------------------------------- 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 "sphereShape.H" +#include "dictionary.H" +#include "vocabs.H" +#include "streams.H" + + +bool pFlow::sphereShape::insertNames +( + const wordVector& names +) +{ + names_.clear(); + uint32 i=0; + for(const auto& nm:names) + { + if(!names_.insertIf(nm, i)) + { + fatalErrorInFunction<< + " repeated name in the list of sphere names: " << names; + return false; + } + i++; + } + names_.rehash(0); + + numShapes_ = names_.size(); + + return true; +} + + +bool pFlow::sphereShape::readDictionary +( + const dictionary& dict +) +{ + diameters_ = dict.getVal("diameters"); + materials_ = dict.getVal("materials"); + auto names = dict.getVal("names"); + + if(diameters_.size() != materials_.size() ) + { + fatalErrorInFunction<< + " number of elements in diameters and properties are not the same in "<< dict.globalName()< names_; + + size_t numShapes_; + + + bool insertNames(const wordVector& names); + + bool readDictionary(const dictionary& dict); + + bool writeDictionary(dictionary& dict)const; + +public: + + // - type info + TypeNameNV("sphereShape"); + + + sphereShape(){} + + sphereShape( + const realVector& diameter, + const wordVector& property, + const wordVector& name + ); + + sphereShape(const sphereShape&) = default; + + sphereShape(sphereShape&&) = default; + + sphereShape& operator=(const sphereShape&) = default; + + sphereShape& operator=(sphereShape&&) = default; + + auto clone()const + { + return makeUnique(*this); + } + + sphereShape* clonePtr()const + { + return new sphereShape(*this); + } + + ~sphereShape() = default; + + //// - Methods + const auto& names()const{ + return names_; + } + + const auto& diameters()const{ + return diameters_; + } + + const auto& materials()const{ + return materials_; + } + + const auto diameter(label i)const{ + return diameters_[i]; + } + + const auto material(label i)const{ + return materials_[i]; + } + + + // name to index + bool nameToIndex(const word& name, uint32& index)const + { + if(auto[iter, found] = names_.findIf(name); found ) + { + index = iter->second; + return true; + } + else + { + index = 0; + return false; + } + } + + uint32 nameToIndex(const word& name)const + { + return names_.at(name); + } + + bool indexToName(uint32 i, word& name)const + { + for(auto& nm: names_) + { + if(nm.second == i) + { + name = nm.first; + return true; + } + } + name = ""; + return false; + } + + bool shapeToDiameter(wordVector& names, realVector& diams)const; + + void diameterMinMax(real& minD, real& maxD)const + { + minD = min(diameters_); + maxD = max(diameters_); + } + + //// - IO operatoin + + // - read from stream/file + bool read(iIstream& is); + + // - write to stream/file + bool write(iOstream& os)const; + + // - read from dictionary + bool read(const dictionary& dict) + { + return readDictionary(dict); + } + + // - write to dictionary + bool write(dictionary& dict)const + { + return writeDictionary(dict); + } + + +}; + +} // pFlow + +#endif //__sphereShape_H__ diff --git a/src/Particles/dynamicPointStructure/dynamicPointStructure.C b/src/Particles/dynamicPointStructure/dynamicPointStructure.C new file mode 100644 index 00000000..1828e87a --- /dev/null +++ b/src/Particles/dynamicPointStructure/dynamicPointStructure.C @@ -0,0 +1,115 @@ +/*------------------------------- 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 "dynamicPointStructure.H" + + +pFlow::dynamicPointStructure::dynamicPointStructure +( + Time& time, + const word& integrationMethod +) +: + time_(time), + integrationMethod_(integrationMethod), + pStruct_( + time_.emplaceObject( + objectFile( + pointStructureFile__, + "", + objectFile::READ_ALWAYS, + objectFile::WRITE_ALWAYS + ) + ) + ), + velocity_( + time_.emplaceObject( + objectFile( + "velocity", + "", + objectFile::READ_ALWAYS, + objectFile::WRITE_ALWAYS + ), + pStruct_, + zero3 + ) + ) + +{ + + Report(1)<< "Creating integration method "<< + greenText(integrationMethod_)<<" for dynamicPointStructure."< integrationPos_; + + uniquePtr integrationVel_; + +public: + + TypeName("dynamicPointStructure"); + + dynamicPointStructure(Time& time, const word& integrationMethod); + + virtual ~dynamicPointStructure() = default; + + inline pointStructure& pStruct() + { + return pStruct_; + } + + inline const pointStructure& pStruct() const + { + return pStruct_; + } + + inline const realx3PointField_D& velocity()const + { + return velocity_; + } + + auto markDeleteOutOfBox(const box& domain) + { + return pStruct_.markDeleteOutOfBox(domain); + } + + bool predict(real dt, realx3PointField_D& acceleration); + + bool correct(real dt, realx3PointField_D& acceleration); + + +}; + +} + + +#endif diff --git a/src/Particles/particles/demParticles.H b/src/Particles/particles/demParticles.H new file mode 100644 index 00000000..4e5f0191 --- /dev/null +++ b/src/Particles/particles/demParticles.H @@ -0,0 +1,47 @@ +/*------------------------------- 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 __demParticles_H__ +#define __demParticles_H__ + +#include "demComponent.H" + +namespace pFlow +{ + + +class demParticles +: + public demComponent +{ + +public: + + demParticles(systemControl& control): + demComponent("particles", control) + {} + + +}; + +} + +#endif diff --git a/src/Particles/particles/particleIdHandler.H b/src/Particles/particles/particleIdHandler.H new file mode 100644 index 00000000..15edda24 --- /dev/null +++ b/src/Particles/particles/particleIdHandler.H @@ -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 __particleIdHandler_H__ +#define __particleIdHandler_H__ + + +#include "pointFields.H" + +namespace pFlow +{ + +class particleIdHandler +{ +protected: + int32 nextId_=0; +public: + particleIdHandler(int32PointField_HD & id) + { + int32 maxID = maxActive(id); + + if( maxID != -1 && id.size() == 0 ) + { + nextId_ = 0; + } + else if( maxID == -1 && id.size()>0 ) + { + nextId_ = 0; + id.modifyOnHost(); + + forAll(i,id) + { + if(id.isActive(i)) + { + id[i] = getNextId(); + } + } + + id.syncViews(); + } + else if( maxID >= static_cast(id.size()) ) + { + nextId_ = maxID + 1; + } + else + { + nextId_ = id.size(); + } + } + + int32 getNextId() + { + return nextId_++; + } + + int32 nextId() const + { + return nextId_; + } +}; + +} + +#endif diff --git a/src/Particles/particles/particles.C b/src/Particles/particles/particles.C new file mode 100644 index 00000000..391c8797 --- /dev/null +++ b/src/Particles/particles/particles.C @@ -0,0 +1,162 @@ +/*------------------------------- 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 "particles.H" + + +pFlow::particles::particles +( + systemControl& control, + const word& integrationMethod +) +: + demParticles(control), + time_(control.time()), + integrationMethod_(integrationMethod), + dynPointStruct_( + control.time(), + integrationMethod + ), + shapeName_( + control.time().emplaceObject( + objectFile( + "shapeName", + "", + objectFile::READ_ALWAYS, + objectFile::WRITE_ALWAYS, + false + ), + pStruct(), + word("NO_NAME_SHAPE") + ) + ), + id_( + control.time().emplaceObject( + objectFile( + "id", + "", + objectFile::READ_IF_PRESENT, + objectFile::WRITE_ALWAYS + ), + pStruct(), + static_cast(-1) + ) + ), + propertyId_( + control.time().emplaceObject( + objectFile( + "propertyId", + "", + objectFile::READ_NEVER, + objectFile::WRITE_NEVER + ), + pStruct(), + static_cast(0) + ) + ), + diameter_( + control.time().emplaceObject( + objectFile( + "diameter", + "", + objectFile::READ_NEVER, + objectFile::WRITE_ALWAYS + ), + pStruct(), + static_cast(0.00000000001) + ) + ), + mass_( + control.time().emplaceObject( + objectFile( + "mass", + "", + objectFile::READ_NEVER, + objectFile::WRITE_ALWAYS + ), + pStruct(), + static_cast(0.0000000001) + ) + ), + accelertion_( + control.time().emplaceObject( + objectFile( + "accelertion", + "", + objectFile::READ_IF_PRESENT, + objectFile::WRITE_ALWAYS + ), + pStruct(), + zero3 + ) + ), + contactForce_( + control.time().emplaceObject( + objectFile( + "contactForce", + "", + objectFile::READ_IF_PRESENT, + objectFile::WRITE_ALWAYS + ), + pStruct(), + zero3 + ) + ), + contactTorque_( + control.time().emplaceObject( + objectFile( + "contactTorque", + "", + objectFile::READ_IF_PRESENT, + objectFile::WRITE_ALWAYS + ), + pStruct(), + zero3 + ) + ), + idHandler_(id_) +{ + + this->subscribe(pStruct()); + +} + +pFlow::uniquePtr> +pFlow::particles::getFieldObjectList()const +{ + auto objListPtr = makeUnique>(); + objListPtr().push_back( + static_cast(&id_) ); + + objListPtr().push_back( + static_cast(&propertyId_) ); + + objListPtr().push_back( + static_cast(&diameter_) ); + + objListPtr().push_back( + static_cast(&mass_) ); + + objListPtr().push_back( + static_cast(&shapeName_) ); + + return objListPtr; +} diff --git a/src/Particles/particles/particles.H b/src/Particles/particles/particles.H new file mode 100644 index 00000000..4c13d24a --- /dev/null +++ b/src/Particles/particles/particles.H @@ -0,0 +1,291 @@ +/*------------------------------- 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 __particles_H__ +#define __particles_H__ + +#include "dynamicPointStructure.H" +#include "particleIdHandler.H" +#include "demParticles.H" +namespace pFlow +{ + +class setFieldList; + +class particles +: + public eventObserver, + public demParticles +{ +protected: + + // owner repository + Time& time_; + + const word integrationMethod_; + + // dynamic point structure for particles + dynamicPointStructure dynPointStruct_; + + // - name of shapes - this is managed by particles + wordPointField& shapeName_; + + // id of particles on host + int32PointField_HD& id_; + + // property id on device + int8PointField_D& propertyId_; + + // diameter / boundig sphere size of particles on device + realPointField_D& diameter_; + + // mass on device + realPointField_D& mass_; + + // - acceleration on device + realx3PointField_D& accelertion_; + + + realx3PointField_D& contactForce_; + + + realx3PointField_D& contactTorque_; + + + // - object handling particle id + particleIdHandler idHandler_; + + virtual uniquePtr> getFieldObjectList()const; + + void zeroForce() + { + contactForce_.fill(zero3); + } + + void zeroTorque() + { + contactTorque_.fill(zero3); + } + +public: + + // type info + TypeName("particles"); + + particles(systemControl& control, const word& integrationMethod); + + inline const auto& time()const { + return time_; + } + + inline auto& time() { + return time_; + } + + inline auto integrationMethod()const + { + return integrationMethod_; + } + + inline const auto& dynPointStruct()const + { + return dynPointStruct_; + } + + inline auto& dynPointStruct() + { + return dynPointStruct_; + } + + inline const auto& pStruct()const{ + return dynPointStruct_.pStruct(); + } + + inline auto& pStruct(){ + return dynPointStruct_.pStruct(); + } + + inline auto size()const{ + return pStruct().size(); + } + + inline auto capacity() const{ + return pStruct().capacity(); + } + + inline auto activePointsMaskD()const{ + return pStruct().activePointsMaskD(); + } + + inline auto numActive()const + { + return pStruct().numActive(); + } + + inline bool allActive()const{ + return pStruct().allActive(); + } + + inline auto activeRange()const{ + return pStruct().activeRange(); + } + + inline auto activePointsMaskH()const{ + return pStruct().activePointsMaskH(); + } + + inline const auto& pointPosition()const{ + return pStruct().pointPosition(); + } + + inline const auto& position()const + { + return pStruct().pointPosition(); + } + + inline const auto& pointVelocity()const{ + return dynPointStruct().velocity(); + } + + inline const auto& velocity()const{ + return dynPointStruct().velocity(); + } + + inline const auto& id()const{ + return id_; + } + + inline auto& id(){ + return id_; + } + + inline const auto& diameter()const{ + return diameter_; + } + + inline auto& diameter(){ + return diameter_; + } + + inline const auto& mass()const{ + return mass_; + } + + inline auto& mass() { + return mass_; + } + + inline const auto& accelertion()const{ + return accelertion_; + } + + inline auto& accelertion(){ + return accelertion_; + } + + inline + realx3PointField_D& contactForce() + { + return contactForce_; + } + + inline + const realx3PointField_D& contactForce() const + { + return contactForce_; + } + + inline + realx3PointField_D& contactTorque() + { + return contactTorque_; + } + + inline + const realx3PointField_D& contactTorque() const + { + return contactTorque_; + } + + inline const auto& propertyId()const{ + return propertyId_; + } + + inline auto& propertyId(){ + return propertyId_; + } + + inline const auto& shapeName()const{ + return shapeName_; + } + + inline auto& shapName(){ + return shapeName_; + } + + bool beforeIteration() override + { + auto domain = this->control().domain(); + + auto numMarked = dynPointStruct_.markDeleteOutOfBox(domain); + /*if(numMarked) + { + output<<"\nNumber of deleted points/particles that are out of domain box: "<zeroForce(); + this->zeroTorque(); + + return true; + } + + virtual + bool insertParticles + ( + const realx3Vector& position, + const wordVector& shapes, + const setFieldList& setField + ) = 0; + + + + virtual + realx3PointField_D& rAcceleration() = 0; + + virtual + const realx3PointField_D& rAcceleration() const = 0; + + virtual + const realVector_D& boundingSphere()const = 0; + + virtual + word shapeTypeName()const = 0; + + virtual + void boundingSphereMinMax(real & minDiam, real& maxDiam)const = 0; + + + +}; // particles + +} // pFlow + +#endif //__particles_H__ diff --git a/src/Property/CMakeLists.txt b/src/Property/CMakeLists.txt new file mode 100644 index 00000000..7c36a2c1 --- /dev/null +++ b/src/Property/CMakeLists.txt @@ -0,0 +1,11 @@ + +set(SourceFiles +property.C +) + +set(link_libs Kokkos::kokkos phasicFlow) + +pFlow_add_library_install(Property SourceFiles link_libs) + + + diff --git a/src/Property/property.C b/src/Property/property.C new file mode 100644 index 00000000..6a68adb0 --- /dev/null +++ b/src/Property/property.C @@ -0,0 +1,155 @@ +/*------------------------------- 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 "property.H" +#include "dictionary.H" + +bool pFlow::property::readDictionary +( + const dictionary& dict +) +{ + + materials_ = dict.getVal("materials"); + + densities_ = dict.getVal("densities"); + + if(materials_.size() != densities_.size() ) + { + fatalErrorInFunction<< + " number of elements in material ("< + ("property", true) + ) +{ + iFstream dictStream(file); + + if(!dict_().read(dictStream)) + { + ioErrorInFile(dictStream.name(), dictStream.lineNumber()); + fatalExit; + } + + if(!readDictionary(dict_())) + { + fatalExit; + } + +} + +pFlow::property::property +( + const dictionary& dict +) +: + dict_ + ( + makeUnique(dict) + ) +{ + + if(!readDictionary(dict_())) + { + fatalExit; + } +} + diff --git a/src/Property/property.H b/src/Property/property.H new file mode 100644 index 00000000..d78a93be --- /dev/null +++ b/src/Property/property.H @@ -0,0 +1,172 @@ +/*------------------------------- 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 __property_H__ +#define __property_H__ + +#include "Vectors.H" +#include "hashMap.H" + +namespace pFlow +{ + +class dictionary; + +class property +{ +protected: + + uniquePtr dict_ = nullptr; + + // - name of materials + wordVector materials_; + + // - density of materials + realVector densities_; + + // fast mapping from name to index + wordHashMap nameIndex_; + + // - number of properties + uint32 numMaterials_ = 0; + + bool readDictionary(const dictionary& dict); + + bool writeDictionary(dictionary& dict)const; + + bool makeNameIndex(); + +public: + + // type info + TypeNameNV("property"); + + // - emptry, for reading from file + property(){} + + property(const wordVector& materils, const realVector& densities); + + property(const fileSystem& file); + + property(const dictionary& dict); + + property(const property& ) = default; + + property(property&& ) = default; + + property& operator= (const property&) = default; + + property& operator= (property&&) = default; + + ~property() = default; + + //// - Methods + + inline const auto& dict()const + { + return dict_(); + } + + inline auto numMaterials()const + { + return numMaterials_; + } + + inline const auto& materials()const{ + return materials_; + } + + // - densities + inline const auto& densities()const{ + return densities_; + } + + inline const word& material(uint32 i)const + { + return materials_[i]; + } + + inline bool material(uint32 i, word& name)const + { + if(isecond; + return true; + } + else + { + idx = 0; + return false; + } + } + + //// - IO operatoins + + // - read from dictionary + bool read(const dictionary& dict) + { + return readDictionary(dict); + } + + // - write to dictionary + bool write(dictionary& dict)const + { + return writeDictionary(dict); + } + +}; + +} + +#endif // __property_H__ diff --git a/src/demComponent/demComponent.H b/src/demComponent/demComponent.H new file mode 100644 index 00000000..60e094dc --- /dev/null +++ b/src/demComponent/demComponent.H @@ -0,0 +1,93 @@ +/*------------------------------- 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 __demComponent_H__ +#define __demComponent_H__ + +#include "systemControl.H" + + + +namespace pFlow +{ + +class demComponent +{ +protected: + + word componentName_; + + + systemControl& control_; + + + Timers timers_; + +public: + + TypeName("demComponent"); + + demComponent(const word& name, systemControl& control) + : + componentName_(name), + control_(control), + timers_(name, &control.timers()) + {} + + virtual ~demComponent() = default; + + + const auto& control()const + { + return control_; + } + + auto& control() + { + return control_; + } + + real dt()const + { + return control_.time().dt(); + } + + auto& timers(){ + return timers_; + } + + const auto& timers() const{ + return timers_; + } + + + virtual bool beforeIteration() = 0; + + + virtual bool iterate() = 0; + + + virtual bool afterIteration() = 0; + +}; + +} + +#endif \ No newline at end of file diff --git a/src/setHelpers/finalize.H b/src/setHelpers/finalize.H new file mode 100644 index 00000000..109e04a2 --- /dev/null +++ b/src/setHelpers/finalize.H @@ -0,0 +1,30 @@ +/*------------------------------- 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 __finalize_H__ +#define __finalize_H__ + +// initilized and finalize should be placed in onc scope +} +pFlow::output<< "\nFinalizing host/device execution space ...."<( + objectFile( + pointStructureFile__, + "", + objectFile::READ_ALWAYS, + objectFile::WRITE_ALWAYS + ) + ); + + +#endif diff --git a/src/setHelpers/setProperty.H b/src/setHelpers/setProperty.H new file mode 100644 index 00000000..0116c9f3 --- /dev/null +++ b/src/setHelpers/setProperty.H @@ -0,0 +1,29 @@ +/*------------------------------- 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 __setProperty_H__ +#define __setProperty_H__ + +Report(0)<<"\nReading proprties . . . "<