125 Commits

Author SHA1 Message Date
5fee39cdd8 Update codingStyle.md
doxygen documentation is added
2025-03-24 17:21:50 +03:30
8fe63cca53 Create codingStyle.md 2025-03-24 14:45:05 +03:30
97e6592524 Merge pull request #192 from wanqing0421/main
fix the cuda build error with multiRotatingAxis
2025-03-20 01:03:30 +03:30
b7d47a65ad fix the cuda build error with multiRotatingAxis 2025-03-20 00:55:46 +08:00
3441b03167 Merge pull request #191 from PhasicFlow/multiRotatingAxis
Multi rotating axis
2025-03-19 18:14:09 +03:30
340f3a551a Multirotating axis motion for version 1.0 2025-03-19 18:10:50 +03:30
be7b2eb632 multiRotaingAxis debug 2025-03-16 22:08:07 +08:00
71057e9575 fill the multiRotatingAxis 2025-03-16 15:15:49 +08:00
797334af87 adapt the multiRotatingAxisMotion to v-1.0 2025-03-16 00:36:38 +08:00
5eef26a6ed Bug fix for memory leak on CPU
- the call for tbb is disabled.
- parallel sort of Kokkos is also very slow.
- for now, std::sort is used for sort, which is more performant than both avaible options. This would be changed anytime any possible solution is found.
2025-03-15 17:15:32 +03:30
892f5395bc Bug fix for observed in particles, getNth and naming for contact lists 2025-03-14 18:33:02 +03:30
b65be8881c end of file corrections 2025-03-13 23:43:52 +03:30
23783b7db9 toteBlender for v-1.0 2025-03-11 17:20:42 +03:30
389e42ee1f Merge branch 'main' of github.com:PhasicFlow/phasicFlow 2025-03-09 21:14:02 +03:30
7f7e06ae7d bug fix for compiling on Ubuntu-24.04LTS 2025-03-09 21:10:41 +03:30
b007426d5d Merge pull request #184 from wanqing0421/main
update benchmarks file and romove the unused comments in rotatingDrumMedium
2025-03-06 21:40:39 +03:30
3ff4ad1469 Merge branch 'PhasicFlow:main' into main 2025-03-07 01:20:31 +08:00
5a5a3c6daf delete the unused comments 2025-03-07 02:53:49 +08:00
bff34bbb9e Merge pull request #185 from ramin1728/RotatingDrumWithBaffles
Rotating drum with baffles
2025-03-06 18:40:59 +03:30
84197bf237 Merge branch 'PhasicFlow:main' into main 2025-03-06 18:39:55 +08:00
169dd73963 Update README.md 2025-03-06 01:06:29 +03:30
baa99c61c0 bug fix for compiling on ubuntu-24 2025-03-06 00:14:14 +03:30
08d0d62d37 ReadMe.md is corrected. 2025-03-05 23:34:37 +03:30
8dcd578a22 RotatingDrumWithBaffles is Updated. 2025-03-05 20:06:25 +03:30
eb62adc87d update benchmarks file and romove the unused comments in rotatingDrumMedium 2025-03-05 22:51:07 +08:00
5394dce7aa Merge pull request #182 from wanqing0421/main
correction for particle insertion region
2025-03-05 17:39:21 +03:30
8038a76699 Correction for settingsDict 2025-03-05 19:59:18 +08:00
c408b60f87 correction for shapes file 2025-03-05 19:56:57 +08:00
ab21acdca5 update rotatingDrumMedium tutorial 2025-03-05 08:37:52 +08:00
07b54c4077 RotatingDrumWithBaffles is Updated. 2025-03-04 22:12:54 +03:30
75f679a234 correction for particle insertion region 2025-03-04 13:57:30 +08:00
c32789f34d Update ReadMe.md 2025-03-02 16:27:00 +03:30
09c303a61e Update interaction 2025-03-02 16:21:51 +03:30
0820e003fc Update particleInsertion 2025-03-02 16:20:56 +03:30
c4c6c2fc45 Update settingsDict 2025-03-02 16:04:20 +03:30
d5dd7af06e Update particlesDict 2025-03-02 16:02:04 +03:30
b03d4825ff Update ReadMe.md 2025-03-02 15:50:04 +03:30
67df8ad206 Merge pull request #181 from ramin1728/RotaryAirLockValve
Rotary air lock valve
2025-03-02 15:47:23 +03:30
2df8133c2d RotaryAirLockValve is Updated. 2025-03-02 12:56:59 +03:30
HRN
dc0504d2fa Merge branch 'main' of github.com:PhasicFlow/phasicFlow into main 2025-02-28 11:37:11 +03:30
HRN
27dfdfa599 stationary motion now does not require the dictionary 2025-02-28 11:36:53 +03:30
8b9a9acd0c RotaryAirLockValve is Updated. 2025-02-27 20:53:20 +03:30
c87c9716ef Update README.md 2025-02-27 20:28:21 +03:30
0ed5b2337c Merge pull request #178 from wanqing0421/main - rotartyDrumSmall
update rotatingDrumSmall tutorial
2025-02-27 20:20:44 +03:30
282d9733fc correctiont for rotatingDrumSmall 2025-02-28 00:43:03 +08:00
cfd188587c revise the readme and domainDict 2025-02-27 23:18:01 +08:00
e8e1081345 update rotatingDrumSmall tutorial 2025-02-26 23:31:52 +08:00
HRN
099e85cfb1 Vector now only accepts one type of allocator, the default allocator 2025-02-26 12:19:36 +03:30
HRN
1cbeb1c963 drum-PeriodictBoundary tutorial added 2025-02-25 22:37:19 +03:30
HRN
a33ec7d8e0 corrections for readMe.md v-1.0 2025-02-24 15:09:38 +03:30
HRN
05ecf37eee box now checks for min and max point consistency 2025-02-24 14:40:29 +03:30
HRN
b44c4de3f6 reading particle position from file for partilclesPhasicFlow 2025-02-24 13:55:56 +03:30
05b256ba39 Merge pull request #170 from wanqing0421/main
New rapid filling implement
2025-02-24 12:34:05 +03:30
25b2e37d93 Merge pull request #175 from ramin1728/banarySystemOfParticles
binarySystemOfParticles is Updated.
2025-02-24 12:33:33 +03:30
a2561f0f12 Merge branch 'PhasicFlow:main' into main 2025-02-23 20:40:36 +08:00
89896c0d69 binarySystemOfParticles is Updated. 2025-02-23 15:59:37 +03:30
HRN
fd45625ce6 cmake_policy 2025-02-21 22:46:31 +03:30
HRN
3e0161a20f version control for cmake_policy 0169 2025-02-21 22:42:11 +03:30
4552d90eac Merge branch 'PhasicFlow:main' into main 2025-02-21 21:58:50 +08:00
HRN
98c8116fd3 rotaryAirlock settingsDict 2025-02-20 21:06:13 +03:30
HRN
fa1211acf8 some minor correction for homogenization silo simulation 2025-02-20 18:48:41 +03:30
HRN
12059faefc corrections for tutorial of screw conveyor 2025-02-20 17:33:55 +03:30
HRN
3954fcf4ab A new screwConveyor tutorial 2025-02-20 17:27:36 +03:30
ef1fa1ddaf Merge branch 'PhasicFlow:main' into main 2025-02-18 19:36:21 +08:00
HRN
354daab7c5 now accepts both kokkos on Home folder and automatic download 2025-02-17 18:39:33 +03:30
HRN
ed4fe6f2f5 Donwloading kokkos and installing tbb is now automatic 2025-02-17 01:13:02 +03:30
d5b9ca4c43 remove rapid filling tutorial 2025-02-16 13:08:09 +08:00
9ccc487a51 Merge branch 'PhasicFlow:main' into main 2025-02-16 12:38:31 +08:00
ae251598a4 update rapid filling 2025-02-16 12:31:11 +08:00
HRN
2a8146c43f add operator << for Set 2025-02-15 22:03:41 +03:30
HRN
fd6b3ebc60 correction for layeredSiloFilling 2025-02-15 22:02:16 +03:30
8e13c377eb Merge pull request #169 from ramin1728/main
layeredSiloFilling is Updated.
2025-02-15 19:57:39 +03:30
5e272cfa1b RotaryAirLockValve is Updated. 2025-02-14 23:30:27 +03:30
252725863f layeredSiloFilling is Updated. 2025-02-14 23:10:46 +03:30
bd4e566dc2 layeredSiloFilling 2025-02-14 23:07:13 +03:30
0532c441a8 Merge pull request #167 from PhasicFlow/develop
bug correction for the time when empty is used
2025-02-14 22:55:22 +03:30
3c1d4d57ad Merge pull request #168 from PhasicFlow/siloHemogenization
New tutorial on hemogenization silo is added
2025-02-14 22:55:09 +03:30
HRN
ff2f1d41e8 New tutorial on hemogenization silo is added 2025-02-14 22:51:46 +03:30
HRN
774afd5f37 bug correction for the time when empty is used 2025-02-14 22:50:28 +03:30
191801b344 Merge pull request #165 from ramin1728/main
binarySystemOfParticles is Updated.
2025-02-14 20:42:14 +03:30
545de300ae Merge pull request #166 from PhasicFlow/develop
edits
2025-02-14 20:40:44 +03:30
HRN
9b3c4f83b9 edits 2025-02-14 20:39:37 +03:30
b315d12357 conveyorBelt is Updated. 2025-02-11 23:35:58 +03:30
7e7184f1c5 binarySystemOfParticles is Updated. 2025-02-11 23:18:29 +03:30
29d922e3c5 A simple rapid filling demo 2025-02-10 23:12:42 +08:00
3aa6be6676 A simple rapid filling 2025-02-10 23:11:33 +08:00
9f489d07cc Merge pull request #163 from PhasicFlow/develop
global damping is activated for velocity and rVelocity in both sphere…
2025-02-10 01:11:30 +03:30
HRN
8466e02d81 global damping is activated for velocity and rVelocity in both sphere and grain solvers 2025-02-10 01:10:13 +03:30
3f1fa4ae90 Merge pull request #162 from wanqing0421/main
bug fix during the build process of cuda mode
2025-02-08 22:39:06 +03:30
f0f185983c bug fix during the build process of cuda mode 2025-02-08 23:47:21 +08:00
cb6d567dab Merge pull request #161 from PhasicFlow/develop
AB5 is added and bug is resolved when particles are being inserted
2025-02-08 18:07:39 +03:30
HRN
db9b1e62e4 AB5 is added and bug is resolved when particles are being inserted 2025-02-08 18:06:30 +03:30
3aff0f1fc6 Merge pull request #160 from PhasicFlow/develop
bug resolve, chekcForCollision is set to true for always, adjustable …
2025-02-07 23:16:48 +03:30
HRN
b9ab015eb1 bug resolve, chekcForCollision is set to true for always, adjustable search box is set to false for always, old hearChanges has been removed 2025-02-07 23:12:53 +03:30
59fbaa91ab Merge pull request #159 from PhasicFlow/develop
All messages are revisited for internal points and boundaries
2025-02-06 11:06:35 +03:30
f84484881c Merge pull request #158 from PhasicFlow/develop
pFlowToVTK now manages when Ctrl+C is used by user
2025-02-06 11:05:21 +03:30
HRN
02e0b72082 All messages are revisited for internal points and boundaries 2025-02-06 11:04:30 +03:30
HRN
edb02ecfc7 pFlowToVTK now manages when Ctrl+C is used by user 2025-02-06 10:51:13 +03:30
1540321a31 Merge pull request #157 from PhasicFlow/develop
the need to provide neighborLength in domain dictionary is lifted. No…
2025-02-03 23:51:31 +03:30
HRN
63bd9c9993 the need to provide neighborLength in domain dictionary is lifted. Now it is optional 2025-02-03 23:49:11 +03:30
f4f5f29e3c Merge pull request #156 from PhasicFlow/develop
Develop
2025-02-03 19:17:08 +03:30
HRN
fac5576df1 periodict boundary condition ->DEBUG and some minor changes boundaries 2025-02-03 19:16:14 +03:30
HRN
f5ba30b901 autoCompelte for time folders and field names 2025-02-03 19:15:08 +03:30
HRN
0f9a19d6ee reduction in boundary class size 2025-02-01 23:33:19 +03:30
d909301f32 Merge pull request #155 from PhasicFlow/develop
Develop
2025-02-01 22:18:23 +03:30
HRN
3b88b6156d Merge branch 'develop' of github.com:PhasicFlow/phasicFlow into develop 2025-02-01 22:15:34 +03:30
HRN
64c041a753 boundaryListPtr is created and other classes were changed accordingly 2025-02-01 22:14:41 +03:30
0410eb6864 Merge pull request #154 from PhasicFlow/main
merge from main
2025-01-31 13:45:24 +03:30
2d7f7ec17f Merge pull request #153 from PhasicFlow/develop
changes in contactSearch and timeControl and timeInfo
2025-01-31 13:44:07 +03:30
HRN
af2572331d some cleanup 2025-01-31 01:06:16 +03:30
HRN
f98a696fc8 changes in contactSearch and timeControl and timeInfo 2025-01-27 14:26:20 +03:30
2168456584 correctoins for running on cuda 2025-01-26 12:19:25 +03:30
6e6cabbefa Merge pull request #152 from PhasicFlow/develop
timeInfo is updated
2025-01-25 19:58:27 +03:30
HRN
42b024e1ed globalDamping is deactivated for future developement 2025-01-25 19:56:01 +03:30
HRN
debb8fd037 bug fix for stridedRange 2025-01-25 19:48:36 +03:30
HRN
0fc9eea561 timeInfo and timeValue
- timeInfo updated
- timeValue type is added = double
- AB2 bug fixed
- updateInterval for broadSearch is activated
2025-01-25 19:18:11 +03:30
HRN
1ccc321c1d Merge branch 'develop' of github.com:PhasicFlow/phasicFlow into develop 2025-01-24 21:15:51 +03:30
HRN
f4b15bc50a timeControl 2025-01-24 21:15:16 +03:30
164eedb27c Merge pull request #151 from PhasicFlow/main
update from main branch
2025-01-24 21:12:53 +03:30
2ec3dfdc7e global damping is added to the code 2025-01-20 21:02:50 +03:30
cb1faf04f8 bug fixes for build with float in cuda 2025-01-20 15:43:56 +03:30
a32786eb8a particlePhasicFlow bug-fix when flag --set-fields-only is used and no shpaeName is set 2025-01-20 15:39:17 +03:30
967bb753aa adding non-linear contact force model for grains 2025-01-20 15:37:48 +03:30
c202f9eaae AB3, AB4 added, and AB2 modified 2025-01-20 14:55:12 +03:30
HRN
8823dd3c94 file headers 2025-01-10 12:27:35 +03:30
320 changed files with 96747 additions and 24960 deletions

View File

@ -3,30 +3,17 @@ cmake_minimum_required(VERSION 3.16 FATAL_ERROR)
# set the project name and version # set the project name and version
project(phasicFlow VERSION 1.0 ) project(phasicFlow VERSION 1.0 )
set(CMAKE_CXX_STANDARD 17 CACHE STRING "" FORCE) set(CMAKE_CXX_STANDARD 20 CACHE STRING "" FORCE)
set(CMAKE_CXX_STANDARD_REQUIRED True) set(CMAKE_CXX_STANDARD_REQUIRED True)
set(CMAKE_INSTALL_PREFIX ${phasicFlow_SOURCE_DIR} CACHE PATH "Install path of phasicFlow" FORCE) set(CMAKE_INSTALL_PREFIX ${phasicFlow_SOURCE_DIR} CACHE PATH "Install path of phasicFlow" FORCE)
set(CMAKE_BUILD_TYPE Debug CACHE STRING "build type" FORCE) set(CMAKE_BUILD_TYPE Release CACHE STRING "build type")
set(BUILD_SHARED_LIBS ON CACHE BOOL "Build using shared libraries" FORCE) set(BUILD_SHARED_LIBS ON CACHE BOOL "Build using shared libraries" FORCE)
mark_as_advanced(FORCE var BUILD_SHARED_LIBS) mark_as_advanced(FORCE var BUILD_SHARED_LIBS)
message(STATUS ${CMAKE_INSTALL_PREFIX}) message(STATUS "Install prefix is:" ${CMAKE_INSTALL_PREFIX})
include(cmake/globals.cmake) include(cmake/globals.cmake)
#Kokkos directory to be included
set(Kokkos_Source_DIR)
if(DEFINED ENV{Kokkos_DIR})
set(Kokkos_Source_DIR $ENV{Kokkos_DIR})
else()
set(Kokkos_Source_DIR $ENV{HOME}/Kokkos/kokkos)
endif()
message(STATUS "Kokkos source directory is ${Kokkos_Source_DIR}")
add_subdirectory(${Kokkos_Source_DIR} ./kokkos)
Kokkos_cmake_settings()
option(pFlow_STD_Parallel_Alg "Use TTB std parallel algorithms" ON) option(pFlow_STD_Parallel_Alg "Use TTB std parallel algorithms" ON)
option(pFlow_Build_Serial "Build phasicFlow and backends for serial execution" OFF) option(pFlow_Build_Serial "Build phasicFlow and backends for serial execution" OFF)
option(pFlow_Build_OpenMP "Build phasicFlow and backends for OpenMP execution" OFF) option(pFlow_Build_OpenMP "Build phasicFlow and backends for OpenMP execution" OFF)
@ -34,6 +21,8 @@ option(pFlow_Build_Cuda "Build phasicFlow and backends for Cuda execution" OFF
option(pFlow_Build_Double "Build phasicFlow with double precision floating-oint variables" ON) option(pFlow_Build_Double "Build phasicFlow with double precision floating-oint variables" ON)
option(pFlow_Build_MPI "Build for MPI parallelization. This will enable multi-gpu run, CPU run on clusters (distributed memory machine). Use this combination Cuda+MPI, OpenMP + MPI or Serial+MPI " OFF) option(pFlow_Build_MPI "Build for MPI parallelization. This will enable multi-gpu run, CPU run on clusters (distributed memory machine). Use this combination Cuda+MPI, OpenMP + MPI or Serial+MPI " OFF)
#for installing the required packages
include(cmake/preReq.cmake)
if(pFlow_Build_Serial) if(pFlow_Build_Serial)
set(Kokkos_ENABLE_SERIAL ON CACHE BOOL "Serial execution" FORCE) set(Kokkos_ENABLE_SERIAL ON CACHE BOOL "Serial execution" FORCE)
@ -46,7 +35,8 @@ elseif(pFlow_Build_OpenMP )
set(Kokkos_ENABLE_OPENMP ON CACHE BOOL "OpenMP 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 OFF CACHE BOOL "Cuda execution" FORCE)
set(Kokkos_ENABLE_CUDA_LAMBDA OFF CACHE BOOL "Cuda execution" FORCE) set(Kokkos_ENABLE_CUDA_LAMBDA OFF CACHE BOOL "Cuda execution" FORCE)
set(Kokkos_DEFAULT_HOST_PARALLEL_EXECUTION_SPACE SERIAL CACHE STRING "" FORCE) set(Kokkos_DEFAULT_HOST_PARALLEL_EXECUTION_SPACE Serial CACHE STRING "" FORCE)
set(Kokkos_DEFAULT_DEVICE_PARALLEL_EXECUTION_SPACE OpenMP CACHE STRING "" FORCE)
set(Kokkos_ENABLE_CUDA_CONSTEXPR OFF CACHE BOOL "Enable constexpr on cuda code" FORCE) set(Kokkos_ENABLE_CUDA_CONSTEXPR OFF CACHE BOOL "Enable constexpr on cuda code" FORCE)
elseif(pFlow_Build_Cuda) elseif(pFlow_Build_Cuda)
set(Kokkos_ENABLE_SERIAL ON CACHE BOOL "Serial execution" FORCE) set(Kokkos_ENABLE_SERIAL ON CACHE BOOL "Serial execution" FORCE)
@ -65,6 +55,7 @@ include(cmake/makeExecutableGlobals.cmake)
configure_file(phasicFlowConfig.H.in phasicFlowConfig.H) configure_file(phasicFlowConfig.H.in phasicFlowConfig.H)
set(CMAKE_EXPORT_COMPILE_COMMANDS ON) set(CMAKE_EXPORT_COMPILE_COMMANDS ON)
#add a global include directory #add a global include directory
include_directories(src/setHelpers src/demComponent "${PROJECT_BINARY_DIR}") include_directories(src/setHelpers src/demComponent "${PROJECT_BINARY_DIR}")

View File

@ -78,9 +78,14 @@ pFlow::grainDEMSystem::grainDEMSystem(
REPORT(0)<< "\nCreating surface geometry for grainDEMSystem . . . "<<END_REPORT; REPORT(0)<< "\nCreating surface geometry for grainDEMSystem . . . "<<END_REPORT;
geometry_ = geometry::create(Control(), Property()); geometry_ = geometry::create(Control(), Property());
REPORT(0)<<"Reading shape dictionary ..."<<END_REPORT;
grains_ = makeUnique<grainShape>(
pFlow::shapeFile__,
&Control().caseSetup(),
Property() );
REPORT(0)<<"\nReading grain particles . . ."<<END_REPORT; REPORT(0)<<"\nReading grain particles . . ."<<END_REPORT;
particles_ = makeUnique<grainFluidParticles>(Control(), Property()); particles_ = makeUnique<grainFluidParticles>(Control(), grains_());
insertion_ = makeUnique<grainInsertion>( insertion_ = makeUnique<grainInsertion>(

View File

@ -46,6 +46,8 @@ protected:
uniquePtr<geometry> geometry_ = nullptr; uniquePtr<geometry> geometry_ = nullptr;
uniquePtr<grainShape> grains_ = nullptr;
uniquePtr<grainFluidParticles> particles_ = nullptr; uniquePtr<grainFluidParticles> particles_ = nullptr;
uniquePtr<grainInsertion> insertion_ = nullptr; uniquePtr<grainInsertion> insertion_ = nullptr;

View File

@ -35,8 +35,8 @@ void pFlow::grainFluidParticles::checkHostMemory()
pFlow::grainFluidParticles::grainFluidParticles( pFlow::grainFluidParticles::grainFluidParticles(
systemControl &control, systemControl &control,
const property &prop) const grainShape& grains)
: grainParticles(control, prop), : grainParticles(control, grains),
fluidForce_( fluidForce_(
objectFile( objectFile(
"fluidForce", "fluidForce",
@ -67,7 +67,7 @@ bool pFlow::grainFluidParticles::beforeIteration()
bool pFlow::grainFluidParticles::iterate() bool pFlow::grainFluidParticles::iterate()
{ {
const auto ti = this->TimeInfo();
accelerationTimer().start(); accelerationTimer().start();
pFlow::grainFluidParticlesKernels::acceleration( pFlow::grainFluidParticlesKernels::acceleration(
control().g(), control().g(),
@ -78,16 +78,16 @@ bool pFlow::grainFluidParticles::iterate()
contactTorque().deviceViewAll(), contactTorque().deviceViewAll(),
fluidTorque_.deviceViewAll(), fluidTorque_.deviceViewAll(),
pStruct().activePointsMaskDevice(), pStruct().activePointsMaskDevice(),
accelertion().deviceViewAll(), acceleration().deviceViewAll(),
rAcceleration().deviceViewAll() rAcceleration().deviceViewAll()
); );
accelerationTimer().end(); accelerationTimer().end();
intCorrectTimer().start(); intCorrectTimer().start();
dynPointStruct().correct(this->dt(), accelertion()); dynPointStruct().correct(ti.dt());
rVelIntegration().correct(this->dt(), rVelocity(), rAcceleration()); rVelIntegration().correct(ti.dt(), rVelocity(), rAcceleration());
intCorrectTimer().end(); intCorrectTimer().end();

View File

@ -59,7 +59,7 @@ protected:
public: public:
/// construct from systemControl and property /// construct from systemControl and property
grainFluidParticles(systemControl &control, const property& prop); grainFluidParticles(systemControl &control, const grainShape& grains);
~grainFluidParticles() override = default; ~grainFluidParticles() override = default;

View File

@ -75,12 +75,18 @@ pFlow::sphereDEMSystem::sphereDEMSystem(
propertyFile__, propertyFile__,
Control().caseSetup().path()); Control().caseSetup().path());
REPORT(0)<< "\nCreating surface geometry for sphereDEMSystem . . . "<<END_REPORT; REPORT(0)<< "\nCreating surface geometry for sphereDEMSystem . . . "<<END_REPORT;
geometry_ = geometry::create(Control(), Property()); geometry_ = geometry::create(Control(), Property());
REPORT(0)<<"Reading shapes dictionary..."<<END_REPORT;
spheres_ = makeUnique<sphereShape>(
pFlow::shapeFile__,
&Control().caseSetup(),
Property());
REPORT(0)<<"\nReading sphere particles . . ."<<END_REPORT; REPORT(0)<<"\nReading sphere particles . . ."<<END_REPORT;
particles_ = makeUnique<sphereFluidParticles>(Control(), Property()); particles_ = makeUnique<sphereFluidParticles>(Control(), spheres_());
insertion_ = makeUnique<sphereInsertion>( insertion_ = makeUnique<sphereInsertion>(

View File

@ -46,6 +46,8 @@ protected:
uniquePtr<geometry> geometry_ = nullptr; uniquePtr<geometry> geometry_ = nullptr;
uniquePtr<sphereShape> spheres_ = nullptr;
uniquePtr<sphereFluidParticles> particles_ = nullptr; uniquePtr<sphereFluidParticles> particles_ = nullptr;
uniquePtr<sphereInsertion> insertion_ = nullptr; uniquePtr<sphereInsertion> insertion_ = nullptr;

View File

@ -32,8 +32,8 @@ void pFlow::sphereFluidParticles::checkHostMemory()
pFlow::sphereFluidParticles::sphereFluidParticles( pFlow::sphereFluidParticles::sphereFluidParticles(
systemControl &control, systemControl &control,
const property &prop) const sphereShape& shpShape)
: sphereParticles(control, prop), : sphereParticles(control, shpShape),
fluidForce_( fluidForce_(
objectFile( objectFile(
"fluidForce", "fluidForce",
@ -65,7 +65,7 @@ bool pFlow::sphereFluidParticles::beforeIteration()
bool pFlow::sphereFluidParticles::iterate() bool pFlow::sphereFluidParticles::iterate()
{ {
const auto ti = this->TimeInfo();
accelerationTimer().start(); accelerationTimer().start();
pFlow::sphereFluidParticlesKernels::acceleration( pFlow::sphereFluidParticlesKernels::acceleration(
control().g(), control().g(),
@ -76,16 +76,16 @@ bool pFlow::sphereFluidParticles::iterate()
contactTorque().deviceViewAll(), contactTorque().deviceViewAll(),
fluidTorque_.deviceViewAll(), fluidTorque_.deviceViewAll(),
pStruct().activePointsMaskDevice(), pStruct().activePointsMaskDevice(),
accelertion().deviceViewAll(), acceleration().deviceViewAll(),
rAcceleration().deviceViewAll() rAcceleration().deviceViewAll()
); );
accelerationTimer().end(); accelerationTimer().end();
intCorrectTimer().start(); intCorrectTimer().start();
dynPointStruct().correct(this->dt(), accelertion()); dynPointStruct().correct(ti.dt());
rVelIntegration().correct(this->dt(), rVelocity(), rAcceleration()); rVelIntegration().correct(ti.dt(), rVelocity(), rAcceleration());
intCorrectTimer().end(); intCorrectTimer().end();

View File

@ -66,7 +66,7 @@ protected:
public: public:
/// construct from systemControl and property /// construct from systemControl and property
sphereFluidParticles(systemControl &control, const property& prop); sphereFluidParticles(systemControl &control, const sphereShape& shpShape);
/// before iteration step /// before iteration step
bool beforeIteration() override; bool beforeIteration() override;

View File

@ -1,48 +1,66 @@
<div align ="center"> <div align="center">
<img src="doc/phasicFlow_logo_github.png" style="width: 400px;"> <img src="doc/phasicFlow_logo_github.png" style="width: 400px;" alt="PhasicFlow Logo">
</div> </div>
## **PhasicFlow: High-Performance Discrete Element Method Simulations**
**PhasicFlow** is a parallel C++ code for performing DEM simulations. It can run on shared-memory multi-core computational units such as multi-core CPUs or GPUs (for now it works on CUDA-enabled GPUs). The parallelization method mainly relies on loop-level parallelization on a shared-memory computational unit. You can build and run PhasicFlow in serial mode on regular PCs, in parallel mode for multi-core CPUs, or build it for a GPU device to off-load computations to a GPU. In its current statues you can simulate millions of particles (up to 80M particles tested) on a single desktop computer. You can see the [performance tests of PhasicFlow](https://github.com/PhasicFlow/phasicFlow/wiki/Performance-of-phasicFlow) in the wiki page. PhasicFlow is a robust, open-source C++ framework designed for the efficient simulation of granular materials using the Discrete Element Method (DEM). Leveraging parallel computing paradigms, PhasicFlow is capable of executing simulations on shared-memory multi-core architectures, including CPUs and NVIDIA GPUs (CUDA-enabled). The core parallelization strategy focuses on loop-level parallelism, enabling significant performance gains on modern hardware. Users can seamlessly transition between serial execution on standard PCs, parallel execution on multi-core CPUs (OpenMP), and accelerated simulations on GPUs. Currently, PhasicFlow supports simulations involving up to 80 million particles on a single desktop workstation. Detailed performance benchmarks are available on the [PhasicFlow Wiki](https://github.com/PhasicFlow/phasicFlow/wiki/Performance-of-phasicFlow).
**MPI** parallelization with dynamic load balancing is under development. With this level of parallelization, PhasicFlow can leverage the computational power of **multi-gpu** workstations or clusters with distributed memory CPUs. **Scalable Parallelism: MPI Integration**
In summary PhasicFlow can have 6 execution modes:
1. Serial on a single CPU core,
2. Parallel on a multi-core computer/node (using OpenMP),
3. Parallel on an nvidia-GPU (using Cuda),
4. Parallel on distributed memory workstation (Using MPI)
5. Parallel on distributed memory workstations with multi-core nodes (using MPI+OpenMP)
6. Parallel on workstations with multiple GPUs (using MPI+Cuda).
## How to build?
You can build PhasicFlow for CPU and GPU executions. The latest release of PhasicFlow is v-0.1. [Here is a complete step-by-step procedure for building phasicFlow-v-0.1.](https://github.com/PhasicFlow/phasicFlow/wiki/How-to-Build-PhasicFlow).
## Online code documentation Ongoing development includes the integration of MPI-based parallelization with dynamic load balancing. This enhancement will extend PhasicFlow's capabilities to distributed memory environments, such as multi-GPU workstations and high-performance computing clusters. Upon completion, PhasicFlow will offer six distinct execution modes:
You can find a full documentation of the code, its features, and other related materials on [online documentation of the code](https://phasicflow.github.io/phasicFlow/)
## How to use PhasicFlow? 1. **Serial Execution:** Single-core CPU.
You can navigate into [tutorials folder](./tutorials) in the phasicFlow folder to see some simulation case setups. If you need more detailed discription, visit our [wiki page tutorials](https://github.com/PhasicFlow/phasicFlow/wiki/Tutorials). 2. **Shared-Memory Parallelism:** Multi-core CPU (OpenMP).
3. **GPU Acceleration:** NVIDIA GPU (CUDA).
4. **Distributed-Memory Parallelism:** MPI.
5. **Hybrid Parallelism:** MPI + OpenMP.
6. **Multi-GPU Parallelism:** MPI + CUDA.
## [PhasicFlowPlus](https://github.com/PhasicFlow/PhasicFlowPlus) # **Build and Installation**
PhasicFlowPlus is and extension to PhasicFlow for simulating particle-fluid systems using resolved and unresolved CFD-DEM. [See the repository of this package.](https://github.com/PhasicFlow/PhasicFlowPlus)
PhasicFlow can be compiled for both CPU and GPU execution.
* **Current Development (v-1.0):** Comprehensive build instructions are available [here](https://github.com/PhasicFlow/phasicFlow/wiki/How-to-build-PhasicFlow%E2%80%90v%E2%80%901.0).
* **Latest Release (v-0.1):** Detailed build instructions are available [here](https://github.com/PhasicFlow/phasicFlow/wiki/How-to-Build-PhasicFlow).
# **Comprehensive Documentation**
In-depth documentation, including code structure, features, and usage guidelines, is accessible via the [online documentation portal](https://phasicflow.github.io/phasicFlow/).
## **Tutorials and Examples**
Practical examples and simulation setups are provided in the [tutorials directory](./tutorials). For detailed explanations and step-by-step guides, please refer to the [tutorial section on the PhasicFlow Wiki](https://github.com/PhasicFlow/phasicFlow/wiki/Tutorials).
# **PhasicFlowPlus: Coupled CFD-DEM Simulations**
PhasicFlowPlus is an extension of PhasicFlow that facilitates the simulation of particle-fluid systems using resolved and unresolved CFD-DEM methods. The repository for PhasicFlowPlus can be found [here](https://github.com/PhasicFlow/PhasicFlowPlus).
## Supporting packages # How to cite PhasicFlow?
* [Kokkos](https://github.com/kokkos/kokkos) from National Technology & Engineering Solutions of Sandia, LLC (NTESS)
* [CLI11 1.8](https://github.com/CLIUtils/CLI11) from University of Cincinnati.
## How to cite PhasicFlow
If you are using PhasicFlow in your research or industrial work, cite the following [article](https://www.sciencedirect.com/science/article/pii/S0010465523001662): If you are using PhasicFlow in your research or industrial work, cite the following [article](https://www.sciencedirect.com/science/article/pii/S0010465523001662):
``` ```
@article{NOROUZI2023108821, @article
title = {PhasicFlow: A parallel, multi-architecture open-source code for DEM simulations}, {
journal = {Computer Physics Communications}, NOROUZI2023108821,
volume = {291}, title = {PhasicFlow: A parallel, multi-architecture open-source code for DEM simulations},
pages = {108821}, journal = {Computer Physics Communications},
year = {2023}, volume = {291},
issn = {0010-4655}, pages = {108821},
doi = {https://doi.org/10.1016/j.cpc.2023.108821}, year = {2023},
url = {https://www.sciencedirect.com/science/article/pii/S0010465523001662}, issn = {0010-4655},
author = {H.R. Norouzi}, doi = {https://doi.org/10.1016/j.cpc.2023.108821},
keywords = {Discrete element method, Parallel computing, CUDA, GPU, OpenMP, Granular flow} url = {https://www.sciencedirect.com/science/article/pii/S0010465523001662},
author = {H.R. Norouzi},
keywords = {Discrete element method, Parallel computing, CUDA, GPU, OpenMP, Granular flow}
} }
``` ```
# **Dependencies**
PhasicFlow relies on the following external libraries:
* **Kokkos:** A performance portability ecosystem developed by National Technology & Engineering Solutions of Sandia, LLC (NTESS). ([https://github.com/kokkos/kokkos](https://github.com/kokkos/kokkos))
* **CLI11 1.8:** A command-line interface parser developed by the University of Cincinnati. ([https://github.com/CLIUtils/CLI11](https://github.com/CLIUtils/CLI11))

View File

@ -3,57 +3,58 @@
| copyright: www.cemf.ir | | copyright: www.cemf.ir |
\* ------------------------------------------------------------------------- */ \* ------------------------------------------------------------------------- */
objectName interaction; objectName interaction;
objectType dicrionary; objectType dictionary;
fileFormat ASCII;
/*---------------------------------------------------------------------------*/
materials (glassMat wallMat); // a list of materials names materials (glassMat wallMat); // a list of materials names
densities (2500.0 2500); // density of materials [kg/m3] densities (2500.0 2500); // density of materials [kg/m3]
contactListType sortedContactList; contactListType sortedContactList;
model
{
contactForceModel nonLinearLimited;
rollingFrictionModel normal;
Yeff (1.0e6 1.0e6 // Young modulus [Pa]
1.0e6);
Geff (0.8e6 0.8e6 // Shear modulus [Pa]
0.8e6);
nu (0.25 0.25 // Poisson's ratio [-]
0.25);
en (0.97 0.85 // coefficient of normal restitution
1.00);
et (1.0 1.0 // coefficient of tangential restitution
1.0);
mu (0.65 0.65 // dynamic friction
0.65);
mur (0.1 0.1 // rolling friction
0.1);
}
contactSearch contactSearch
{ {
method NBS; method NBS;
wallMapping cellMapping;
NBSInfo updateInterval 10;
{
updateFrequency 10; // each 20 timesteps, update neighbor list
sizeRatio 1.05; // bounding box size to particle diameter (max)
}
cellMappingInfo sizeRatio 1.1;
{
updateFrequency 10; // each 20 timesteps, update neighbor list
cellExtent 0.6; // bounding box for particle-wall search (> 0.5)
}
cellExtent 0.55;
adjustableBox Yes;
} }
model
{
contactForceModel nonLinearLimited;
rollingFrictionModel normal;
/*
Property (glassMat-glassMat glassMat-wallMat
wallMat-wallMat);
*/
Yeff (1.0e6 1.0e6
1.0e6); // Young modulus [Pa]
Geff (0.8e6 0.8e6
0.8e6); // Shear modulus [Pa]
nu (0.25 0.25
0.25); // Poisson's ratio [-]
en (0.97 0.85
1.00); // coefficient of normal restitution
mu (0.65 0.65
0.65); // dynamic friction
mur (0.1 0.1
0.1); // rolling friction
}

View File

@ -3,12 +3,10 @@
| copyright: www.cemf.ir | | copyright: www.cemf.ir |
\* ------------------------------------------------------------------------- */ \* ------------------------------------------------------------------------- */
objectName particleInsertion; objectName particleInsertion;
objectType dicrionary; objectType dicrionary;
fileFormat ASCII;
active no; // is insertion active? active No; // is insertion active?
collisionCheck No; // not implemented for yes

View File

@ -2,12 +2,14 @@
| phasicFlow File | | phasicFlow File |
| copyright: www.cemf.ir | | copyright: www.cemf.ir |
\* ------------------------------------------------------------------------- */ \* ------------------------------------------------------------------------- */
objectName particleInsertion;
objectType dicrionary; objectName shapes;
fileFormat ASCII; objectType dictionary;
fileFormat ASCII;
/*---------------------------------------------------------------------------*/ /*---------------------------------------------------------------------------*/
active No; // is insertion active -> Yes or No
collisionCheck No; // is checked -> Yes or No names (glassBead); // names of shapes
diameters (0.003); // diameter of shapes
materials (glassMat); // material names for shapes

View File

@ -1,11 +0,0 @@
/* -------------------------------*- C++ -*--------------------------------- *\
| phasicFlow File |
| copyright: www.cemf.ir |
\* ------------------------------------------------------------------------- */
objectName sphereDict;
objectType sphereShape;
names (glassBead); // names of shapes
diameters (0.003); // diameter of shapes
materials (glassMat); // material names for shapes

View File

@ -0,0 +1,46 @@
/* -------------------------------*- C++ -*--------------------------------- *\
| phasicFlow File |
| copyright: www.cemf.ir |
\* ------------------------------------------------------------------------- */
objectName domainDict;
objectType dictionary;
fileFormat ASCII;
/*---------------------------------------------------------------------------*/
globalBox // Simulation domain: every particles that goes outside this domain will be deleted
{
min (-0.2 -0.2 -0.0);
max ( 0.2 0.2 1.6);
}
boundaries
{
left
{
type exit; // other options: periodic, reflective
}
right
{
type exit; // other options: periodic, reflective
}
bottom
{
type exit; // other options: periodic, reflective
}
top
{
type exit; // other options: periodic, reflective
}
rear
{
type exit; // other options: periodic, reflective
}
front
{
type exit; // other options: periodic, reflective
}
}

View File

@ -3,61 +3,84 @@
| copyright: www.cemf.ir | | copyright: www.cemf.ir |
\* ------------------------------------------------------------------------- */ \* ------------------------------------------------------------------------- */
objectName geometryDict; objectName geometryDict;
objectType dictionary; objectType dictionary;
fileFormat ASCII;
motionModel rotatingAxisMotion; motionModel rotatingAxis; // motion model: rotating object around an axis
rotatingAxisInfo // information for rotatingAxisMotion motion model
{
rotAxis
{
p1 (0.0 0.0 0.0); // first point for the axis of rotation
p2 (0.0 0.0 1.0); // second point for the axis of rotation
omega 1.256; // rotation speed (rad/s) => 12 rpm
}
}
surfaces surfaces
{ {
cylinder
{
type cylinderWall; // type of the wall
cylinder p1 (0.0 0.0 0.0); // begin point of cylinder axis
{
type cylinderWall;
p1 (0.0 0.0 0.0);
p2 (0.0 0.0 1.6);
radius1 0.2;
radius2 0.2;
resolution 24;
material wallMat;
motion rotAxis;
}
p2 (0.0 0.0 1.6); // end point of cylinder axis
wall1 radius1 0.2; // radius at p1
{
type planeWall;
p1 (-0.2 -0.2 0.0);
p2 ( 0.2 -0.2 0.0);
p3 ( 0.2 0.2 0.0);
p4 (-0.2 0.2 0.0);
material wallMat;
motion rotAxis;
}
/* radius2 0.2; // radius at p2
This is a plane wall at the front end of cylinder
*/
wall2
{
type planeWall;
p1 (-0.2 -0.2 1.6);
p2 ( 0.2 -0.2 1.6);
p3 ( 0.2 0.2 1.6);
p4 (-0.2 0.2 1.6);
material wallMat;
motion rotAxis;
}
} resolution 24; // number of divisions
// information for rotatingAxisMotion motion model material wallMat; // material name of this wall
rotatingAxisMotionInfo
{ motion rotAxis; // motion component name
rotAxis }
{
p1 (0.0 0.0 0.0); /*
p2 (0.0 0.0 1.0); This is a plane wall at the rear end of cylinder
omega 1.256; // rotation speed (rad/s) => 12 rpm */
}
wall1
{
type planeWall; // type of the wall
p1 (-0.2 -0.2 0.0); // first point of the wall
p2 ( 0.2 -0.2 0.0); // second point
p3 ( 0.2 0.2 0.0); // third point
p4 (-0.2 0.2 0.0); // fourth point
material wallMat; // material name of the wall
motion rotAxis; // motion component name
}
/*
This is a plane wall at the front end of cylinder
*/
wall2
{
type planeWall; // type of the wall
p1 (-0.2 -0.2 1.6); // first point of the wall
p2 ( 0.2 -0.2 1.6); // second point
p3 ( 0.2 0.2 1.6); // third point
p4 (-0.2 0.2 1.6); // fourth point
material wallMat; // material name of the wall
motion rotAxis; // motion component name
}
} }

View File

@ -3,42 +3,45 @@
| copyright: www.cemf.ir | | copyright: www.cemf.ir |
\* ------------------------------------------------------------------------- */ \* ------------------------------------------------------------------------- */
objectName particlesDict; objectName particlesDict;
objectType dictionary; objectType dictionary;
fileFormat ASCII;
setFields setFields
{ {
defaultValue
{
velocity realx3 (0 0 0); // linear velocity (m/s)
acceleration realx3 (0 0 0); // linear acceleration (m/s2)
rotVelocity realx3 (0 0 0); // rotational velocity (rad/s)
shapeName word glassBead; // name of the particle shape
}
defaultValue selectors
{ {}
velocity realx3 (0 0 0); // linear velocity (m/s)
acceleration realx3 (0 0 0); // linear acceleration (m/s2)
rotVelocity realx3 (0 0 0); // rotational velocity (rad/s)
shapeName word glassBead; // name of the particle shape
}
selectors
{}
} }
positionParticles positionParticles
{ {
method positionOrdered; method ordered;
maxNumberOfParticles 4000001; orderedInfo
mortonSorting Yes; {
diameter 0.003; // minimum space between centers of particles
cylinder // box for positioning particles numPoints 4000000; // number of particles in the simulation
{
p1 ( 0.0 0.0 0.01); // lower corner point of the box
p2 ( 0.0 0.0 1.59); // upper corner point of the box
radius 0.195;
}
positionOrderedInfo axisOrder (z x y); // axis order for filling the space with particles
{ }
diameter 0.003; // minimum space between centers of particles
numPoints 4000000; // number of particles in the simulation regionType cylinder; // other options: box and sphere
axisOrder (z x y); // axis order for filling the space with particles
} cylinderInfo // cylinder for positioning particles
{
p1 (0.0 0.0 0.01); // lower corner point of the box
p2 (0.0 0.0 1.59); // upper corner point of the box
radius 0.195; // radius of cylinder
}
} }

View File

@ -2,31 +2,33 @@
| phasicFlow File | | phasicFlow File |
| copyright: www.cemf.ir | | copyright: www.cemf.ir |
\* ------------------------------------------------------------------------- */ \* ------------------------------------------------------------------------- */
objectName settingsDict; objectName settingsDict;
objectType dictionary;; objectType dictionary;
fileFormat ASCII;
/*---------------------------------------------------------------------------*/
run rotatingDrum_4MParticles;
run rotatingDrum_1; dt 0.00001; // time step for integration (s)
dt 0.00001; // time step for integration (s) startTime 0; // start time for simulation
startTime 0; // start time for simulation endTime 10; // end time for simulation
endTime 10; // end time for simulation saveInterval 0.2; // time interval for saving the simulation
saveInterval 0.2; // time interval for saving the simulation timePrecision 5; // maximum number of digits for time folder
timePrecision 5; // maximum number of digits for time folder g (0 -9.8 0); // gravity vector (m/s2)
g (0 -9.8 0); // gravity vector (m/s2) includeObjects (diameter); // save necessary (i.e., required) data on disk
domain // exclude unnecessary data from saving on disk
{ excludeObjects (rVelocity.dy1 pStructPosition.dy1 pStructVelocity.dy1);
min (-0.2 -0.2 -0.0);
max ( 0.2 0.2 1.6);
}
integrationMethod AdamsBashforth3; // integration method integrationMethod AdamsBashforth2; // integration method
timersReport Yes; writeFormat binary; // data writting format (ascii or binary)
timersReportInterval 0.01; timersReport Yes;
timersReportInterval 0.05;

View File

@ -1,17 +1,55 @@
PF_cFlags="--description --help --version" PF_cFlags="--description --help --version"
AllTimeFolders=
__getAllTime(){
files=( $(ls) )
deleteFiles=(settings caseSetup cleanThisCase VTK runThisCase stl postprocess postProcess)
declare -A delk
for del in "${deleteFiles[@]}" ; do delk[$del]=1 ; done
# Tag items to remove, based on
for k in "${!files[@]}" ; do
[ "${delk[${files[$k]}]-}" ] && unset 'files[k]'
done
# Compaction
COMPREPLY=("${files[@]}")
AllTimeFolders="${files[@]}"
}
__getFields(){
__getAllTime
local -A unique_files=()
for dir in $AllTimeFolders; do
# Check if the directory exists
if [ ! -d "$dir" ]; then
continue # Skip to the next directory
fi
files_in_dir=$(find "$dir" -maxdepth 1 -type f -printf '%f\n' | sort -u)
# Add filenames to the associative array (duplicates will be overwritten)
while IFS= read -r filename; do
unique_files["$filename"]=1 # Just the key is important, value can be anything
done <<< "$files_in_dir"
done
COMPREPLY=("${!unique_files[@]}")
AllTimeFolders=
}
_pFlowToVTK(){ _pFlowToVTK(){
if [ "$3" == "--time" ]; then if [ "$3" == "--time" ] ; then
COMPREPLY=( $(ls) ) __getAllTime
elif [ "$3" == "--fields" ]; then
__getFields
else else
COMPREPLY=( $(compgen -W "$PF_cFlags --binary --no-geometry --no-particles --out-folder --time --separate-surfaces --fields" -- "$2") ) COMPREPLY=( $(compgen -W "$PF_cFlags --binary --no-geometry --no-particles --out-folder --time --separate-surfaces --fields" -- "$2") )
fi fi
} }
complete -F _pFlowToVTK pFlowToVTK complete -F _pFlowToVTK pFlowToVTK
_postprocessPhasicFlow(){ _postprocessPhasicFlow(){
if [ "$3" == "--time" ]; then if [ "$3" == "--time" ]; then
COMPREPLY=( $(ls) ) __getAllTime
else else
COMPREPLY=( $(compgen -W "$PF_cFlags --out-folder --time --zeroFolder" -- "$2") ) COMPREPLY=( $(compgen -W "$PF_cFlags --out-folder --time --zeroFolder" -- "$2") )
fi fi

54
cmake/preReq.cmake Normal file
View File

@ -0,0 +1,54 @@
if(pFlow_STD_Parallel_Alg)
# Check if libtbb-dev is installed
execute_process(
COMMAND dpkg -s libtbb-dev
RESULT_VARIABLE TBB_IS_INSTALLED
OUTPUT_QUIET
ERROR_QUIET)
if(NOT TBB_IS_INSTALLED EQUAL 0)
message(STATUS "libtbb-dev not found. Installing libtbb-dev...")
execute_process(
COMMAND sudo apt-get update
COMMAND sudo apt-get install -y libtbb-dev
RESULT_VARIABLE TBB_INSTALL_RESULT)
if(NOT TBB_INSTALL_RESULT EQUAL 0)
message(FATAL_ERROR "Failed to install libtbb-dev")
endif()
else()
message(STATUS "libtbb-dev is already installed.")
endif()
endif()
# Kokkos folder creation
set(Kokkos_Source_DIR $ENV{HOME}/Kokkos/kokkos)
if(NOT EXISTS "${Kokkos_Source_DIR}/CMakeLists.txt")
# Check CMake version and set policy CMP0169 if CMake version is 3.30 or higher
if(${CMAKE_VERSION} VERSION_GREATER_EQUAL "3.30")
cmake_policy(SET CMP0169 OLD)
endif()
include(FetchContent)
FetchContent_Declare(
kokkos
GIT_REPOSITORY https://github.com/kokkos/kokkos.git
GIT_TAG 4.4.01
)
FetchContent_GetProperties(kokkos)
if(NOT kokkos_POPULATED)
message(STATUS "Kokkos source directory not found. Downloading Kokkos version 4.4.1 ...")
FetchContent_Populate(kokkos)
set(Kokkos_Source_DIR ${kokkos_SOURCE_DIR})
endif()
endif()
message(STATUS "Kokkos source directory is ${Kokkos_Source_DIR}")
add_subdirectory(${Kokkos_Source_DIR} ./kokkos)
#Kokkos_cmake_settings()

231
doc/mdDocs/codingStyle.md Normal file
View File

@ -0,0 +1,231 @@
# PhasicFlow Coding Style Guidelines
This document outlines the coding style guidelines for the PhasicFlow codebase.
Adhering to these guidelines ensures consistency, readability, and maintainability of the project.
## 1. FormattingIndentation:
* Use 4 spaces for every logical level and block.
* Line Spacing: Leave two empty lines between sections (e.g., between functions in a .cpp file, between class members).
## 2. Naming ConventionsGeneral Naming:
* All names should start with lowercase letters, except for special names (e.g., Ergun, Hertz).
* Macro names start with Upper case or all the letters are in UPPER case.
* Compound Names: For compound names, the first word starts with a lowercase letter, and subsequent words start with an uppercase letter (e.g., boundaryBase, motionModel).
## 3. File Structure
* Header Files: Use the .hpp extension for header files.
* Source Files: Use the .cpp extension for source files.
* Header and Source File Headers: All header and source files must include a standardized header that describes the project's intention and licensing information.
* File Naming: Header and source file names should correspond to the class they contain. Aim for one class per file.
* Inline Functions: Place inline functions in a separate classNameI.hpp file to avoid cluttering the main header file.
## 4. Class DesignClass Member Order:
* Private members and methods
* Private static members and methods
* Public methods
* Public static methods
* Enumerations and Nested Classes: Declare enumerations and nested classes before all class members and methods.
* Special Functions: Each class must explicitly define all special functions:Constructor, Copy constructor and assignment operator, Move constructor and assignment operator
* Destructor: Each class must have an explicit destructor declaration:`~className() = default`; or `~className();`
* Interface classes or classes with virtual methods must have a virtual destructor.
* Virtual Method Overrides: When implementing a `virtual` method from a base class in a derived class, use the `override` keyword. The same applies to derived class destructors.
## 5. NamespacesOfficial Namespace:
The official namespace for the codebase is pFlow. All entities should be defined within this namespace.
### Example File Structure
```
src/
├── componentName1/
│ ├── componentName1.hpp
│ ├── componentName1.cpp
│ ├── componentName1I.hpp
│ └── ...
└── componentName2/
├── componentName2.hpp
├── componentName2.cpp
└── ...
```
### Example Class Structure
```C++
namespace pFlow
{
class MyClass
{
public:
enum class MyEnum
{
Value1,
Value2
};
class NestedClass
{
// ...
};
private:
int privateMember_;
void privateMethod();
static int privateStaticMember_;
static void privateStaticMethod();
public:
MyClass();
MyClass(const MyClass& other);
MyClass(MyClass&& other);
MyClass& operator=(const MyClass& other);
MyClass& operator=(MyClass&& other);
~MyClass();
void publicMethod();
static void publicStaticMethod();
};
// assuming base class has virtual methods
class DerivedClass
:
public BaseClass
{
public:
...
~DerivedClass() override;
void virtualMethod() override;
};
} // namespace pFlow
```
## 6. Doxygen Documentation
### 6.1. Ruls
provide the documentations in the header files only. In rare cases where you are generating documentations for executables, the doxygen documentation can be supplied in the .cpp file.
* **General Doxygen Style:**
* Use `///` for short documentations for methods and members.
* Use `/** */` for classes and main methods which play a significant role in the class or code.
* Place Doxygen comments *before* the declaration of the entity being documented (e.g., class, function, variable).
* Use `@param` to document function parameters, `@return` for return values, `@brief` for a short description, and `@details` for a more in-depth explanation.
* Use Markdown syntax within Doxygen comments for formatting.
* **File Headers:** Each file should contain a Doxygen comment at the top, including:
* `@file` : The name of the file.
* `@brief`: A brief description of the file's purpose.
* `@author`: The author(s) of the file.
* `@date` : The date of creation or last modification.
* **Class Documentation:**
* Use `/** */` for class documentation.
* Provide a `@brief` description of the class.
* Use `@tparam` to document template parameters.
* Document the purpose of the class, its invariants, and how it should be used.
* **Function/Method Documentation:**
* Use `///` for short documentations.
* Use `/** */` for main methods which play a significant role.
* Provide a `@brief` description of the function.
* Use `@param` to describe each parameter, including its purpose and whether it is an input, output, or input/output parameter.
* Use `@return` to describe the return value, including its meaning and possible values.
* Use `@pre` to document any preconditions that must be met before calling the function.
* Use `@post` to document any postconditions that will be true after the function returns.
* Use `@throws` to document any exceptions that the function may throw.
* Use `@details` for a more detailed explanation of the function's behavior, algorithms, or any other relevant information.
* **Variable Documentation:**
* Use `///<` for single-line documentation of variables.
### 6.2. Doxygen Documentation Examples
* **Class example**
```cpp
/**
* @brief Represents a particle in the simulation.
* @details This class stores the position, velocity, and other physical
* properties of a particle.
*/
class Particle
{
private:
Point position_; ///< The current position of the particle.
Vector velocity_; ///< The current velocity of the particle.
double mass_; ///< The mass of the particle.
public:
/// Constructs a particle with default values.
Particle();
/**
* @brief Updates the position of the particle based on its velocity
* and the given time step.
* @param deltaTime The time elapsed since the last update, in seconds.
*/
void updatePosition(const timeInfo& ti );
/// Gets the current position of the particle.
Point getPosition() const;
};
```
* **Function Example**
```cpp
/**
* @brief Calculates the distance between two points.
* @param p1 The first point.
* @param p2 The second point.
* @return The distance between the two points.
*/
double calculateDistance(const Point& p1, const Point& p2)
{
// Implementation
return 0.0;
}
/// Returns the velocity of the particle.
Vector getVelocity() const
{
return velocity_;
}
```

View File

@ -17,10 +17,17 @@ Licence:
implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
-----------------------------------------------------------------------------*/ -----------------------------------------------------------------------------*/
REPORT(0)<<"Reading shape dictionary ..."<<END_REPORT;
pFlow::grainShape grains
(
pFlow::shapeFile__,
&Control.caseSetup(),
proprties
);
// //
REPORT(0)<<"\nReading sphere particles . . ."<<END_REPORT; REPORT(0)<<"\nReading sphere particles . . ."<<END_REPORT;
pFlow::grainParticles grnParticles(Control, proprties); pFlow::grainParticles grnParticles(Control, grains);
// //
REPORT(0)<<"\nCreating particle insertion object . . ."<<END_REPORT; REPORT(0)<<"\nCreating particle insertion object . . ."<<END_REPORT;

View File

@ -18,8 +18,16 @@ Licence:
-----------------------------------------------------------------------------*/ -----------------------------------------------------------------------------*/
REPORT(0)<<"Reading shape dictionary ..."<<END_REPORT;
pFlow::sphereShape spheres
(
pFlow::shapeFile__,
&Control.caseSetup(),
proprties
);
// //
REPORT(0)<<"\nReading sphere particles . . ."<<END_REPORT; REPORT(0)<<"\nReading sphere particles . . ."<<END_REPORT;
pFlow::sphereParticles sphParticles(Control, proprties); pFlow::sphereParticles sphParticles(Control, spheres);
WARNING<<"Particle insertion has not been set yet!"<<END_WARNING; WARNING<<"Particle insertion has not been set yet!"<<END_WARNING;

View File

@ -17,10 +17,18 @@ Licence:
implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
-----------------------------------------------------------------------------*/ -----------------------------------------------------------------------------*/
REPORT(0)<<"Reading shapes dictionary..."<<END_REPORT;
pFlow::sphereShape spheres
(
pFlow::shapeFile__,
&Control.caseSetup(),
proprties
);
// //
REPORT(0)<<"\nReading sphere particles . . ."<<END_REPORT; REPORT(0)<<"\nReading sphere particles . . ."<<END_REPORT;
pFlow::sphereParticles sphParticles(Control, proprties);
pFlow::sphereParticles sphParticles(Control, spheres);
// //
REPORT(0)<<"\nCreating particle insertion object . . ."<<END_REPORT; REPORT(0)<<"\nCreating particle insertion object . . ."<<END_REPORT;

View File

@ -74,7 +74,7 @@ pFlow::initialize_pFlowProcessors();
do do
{ {
//Ping;
if(! sphInsertion.insertParticles( if(! sphInsertion.insertParticles(
Control.time().currentIter(), Control.time().currentIter(),
Control.time().currentTime(), Control.time().currentTime(),
@ -90,21 +90,25 @@ pFlow::initialize_pFlowProcessors();
// set force to zero, predict, particle deletion and etc. // set force to zero, predict, particle deletion and etc.
sphParticles.beforeIteration(); sphParticles.beforeIteration();
//Ping;
sphInteraction.beforeIteration(); sphInteraction.beforeIteration();
sphInteraction.iterate(); sphInteraction.iterate();
surfGeometry.iterate(); surfGeometry.iterate();
//Ping;
sphParticles.iterate(); sphParticles.iterate();
//Ping;
sphInteraction.afterIteration(); sphInteraction.afterIteration();
//Ping;
surfGeometry.afterIteration(); surfGeometry.afterIteration();
//Ping;
sphParticles.afterIteration(); sphParticles.afterIteration();
//Ping;
}while(Control++); }while(Control++);

View File

@ -21,7 +21,6 @@ Licence:
template<typename MotionModel> template<typename MotionModel>
bool pFlow::geometryMotion<MotionModel>::findMotionIndex() bool pFlow::geometryMotion<MotionModel>::findMotionIndex()
{ {
if(motionComponentName().size() != numSurfaces() ) if(motionComponentName().size() != numSurfaces() )
{ {
fatalErrorInFunction<< fatalErrorInFunction<<

View File

@ -28,4 +28,4 @@ template class pFlow::geometryMotion<pFlow::stationaryWall>;
template class pFlow::geometryMotion<pFlow::conveyorBeltMotion>; template class pFlow::geometryMotion<pFlow::conveyorBeltMotion>;
//template class pFlow::geometryMotion<pFlow::multiRotatingAxisMotion>; template class pFlow::geometryMotion<pFlow::multiRotatingAxisMotion>;

View File

@ -25,7 +25,7 @@ Licence:
#include "stationaryWall.hpp" #include "stationaryWall.hpp"
#include "rotatingAxisMotion.hpp" #include "rotatingAxisMotion.hpp"
#include "conveyorBeltMotion.hpp" #include "conveyorBeltMotion.hpp"
//#include "multiRotatingAxisMotion.hpp" #include "multiRotatingAxisMotion.hpp"
#include "vibratingMotion.hpp" #include "vibratingMotion.hpp"
@ -40,10 +40,7 @@ using stationaryGeometry = geometryMotion<stationaryWall>;
using conveyorBeltMotionGeometry = geometryMotion<conveyorBeltMotion>; using conveyorBeltMotionGeometry = geometryMotion<conveyorBeltMotion>;
//typedef geometryMotion<multiRotatingAxisMotion> multiRotationAxisMotionGeometry; using multiRotationAxisMotionGeometry = geometryMotion<multiRotatingAxisMotion>;
} }

View File

@ -37,7 +37,8 @@ bool intAllActive(
real dt, real dt,
realx3Field_D& y, realx3Field_D& y,
realx3PointField_D& dy, realx3PointField_D& dy,
realx3PointField_D& dy1) realx3PointField_D& dy1,
real damping = 1.0)
{ {
auto d_dy = dy.deviceView(); auto d_dy = dy.deviceView();
@ -49,7 +50,7 @@ bool intAllActive(
"AdamsBashforth2::correct", "AdamsBashforth2::correct",
rpIntegration (activeRng.start(), activeRng.end()), rpIntegration (activeRng.start(), activeRng.end()),
LAMBDA_HD(uint32 i){ LAMBDA_HD(uint32 i){
d_y[i] += dt*(static_cast<real>(1.5) * d_dy[i] - static_cast<real>(0.5) * d_dy1[i]); d_y[i] += damping * dt*(static_cast<real>(1.5) * d_dy[i] - static_cast<real>(0.5) * d_dy1[i]);
d_dy1[i] = d_dy[i]; d_dy1[i] = d_dy[i];
}); });
Kokkos::fence(); Kokkos::fence();
@ -62,7 +63,8 @@ bool intScattered
real dt, real dt,
realx3Field_D& y, realx3Field_D& y,
realx3PointField_D& dy, realx3PointField_D& dy,
realx3PointField_D& dy1 realx3PointField_D& dy1,
real damping = 1.0
) )
{ {
@ -78,7 +80,7 @@ bool intScattered
LAMBDA_HD(uint32 i){ LAMBDA_HD(uint32 i){
if( activeP(i)) if( activeP(i))
{ {
d_y[i] += dt*(static_cast<real>(1.5) * d_dy[i] - static_cast<real>(0.5) * d_dy1[i]); d_y[i] += damping * dt*(static_cast<real>(1.5) * d_dy[i] - static_cast<real>(0.5) * d_dy1[i]);
d_dy1[i] = d_dy[i]; d_dy1[i] = d_dy[i];
} }
}); });
@ -112,8 +114,11 @@ pFlow::AdamsBashforth2::AdamsBashforth2
zero3, zero3,
zero3 zero3
), ),
initialValField_(initialValField),
boundaryList_(pStruct, method, *this) boundaryList_(pStruct, method, *this)
{} {
realx3PointField_D::addEvent(message::ITEMS_INSERT);
}
void pFlow::AdamsBashforth2::updateBoundariesSlaveToMasterIfRequested() void pFlow::AdamsBashforth2::updateBoundariesSlaveToMasterIfRequested()
{ {
@ -142,18 +147,19 @@ bool pFlow::AdamsBashforth2::correct
( (
real dt, real dt,
realx3PointField_D& y, realx3PointField_D& y,
realx3PointField_D& dy realx3PointField_D& dy,
real damping
) )
{ {
auto& dy1l = dy1(); auto& dy1l = dy1();
bool success = false; bool success = false;
if(dy1l.isAllActive()) if(dy1l.isAllActive())
{ {
success = intAllActive(dt, y.field(), dy, dy1l); success = intAllActive(dt, y.field(), dy, dy1(), damping);
} }
else else
{ {
success = intScattered(dt, y.field(), dy, dy1l); success = intScattered(dt, y.field(), dy, dy1(), damping);
} }
success = success && boundaryList_.correct(dt, y, dy); success = success && boundaryList_.correct(dt, y, dy);
@ -171,11 +177,11 @@ bool pFlow::AdamsBashforth2::correctPStruct(
bool success = false; bool success = false;
if(dy1l.isAllActive()) if(dy1l.isAllActive())
{ {
success = intAllActive(dt, pStruct.pointPosition(), vel, dy1l); success = intAllActive(dt, pStruct.pointPosition(), vel, dy1());
} }
else else
{ {
success = intScattered(dt, pStruct.pointPosition(), vel, dy1l); success = intScattered(dt, pStruct.pointPosition(), vel, dy1());
} }
success = success && boundaryList_.correctPStruct(dt, pStruct, vel); success = success && boundaryList_.correctPStruct(dt, pStruct, vel);
@ -183,11 +189,21 @@ bool pFlow::AdamsBashforth2::correctPStruct(
return success; return success;
} }
/*bool pFlow::AdamsBashforth2::hearChanges
bool pFlow::AdamsBashforth2::setInitialVals( (
const int32IndexContainer& newIndices, const timeInfo &ti,
const realx3Vector& y) const message &msg,
const anyList &varList
)
{ {
return true; if(msg.equivalentTo(message::ITEMS_INSERT))
} {
return insertValues(varList, initialValField_.deviceViewAll(), dy1());
}
else
{
return realx3PointField_D::hearChanges(ti, msg, varList);
}
}*/

View File

@ -41,10 +41,14 @@ class AdamsBashforth2
{ {
private: private:
const realx3Field_D& initialValField_;
boundaryIntegrationList boundaryList_; boundaryIntegrationList boundaryList_;
friend class processorAB2BoundaryIntegration; friend class processorAB2BoundaryIntegration;
protected:
const auto& dy1()const const auto& dy1()const
{ {
return static_cast<const realx3PointField_D&>(*this); return static_cast<const realx3PointField_D&>(*this);
@ -55,6 +59,16 @@ private:
return static_cast<realx3PointField_D&>(*this); return static_cast<realx3PointField_D&>(*this);
} }
auto& initialValField()
{
return initialValField_;
}
boundaryIntegrationList& boundaryList()
{
return boundaryList_;
}
public: public:
/// Class info /// Class info
@ -70,7 +84,7 @@ public:
const realx3Field_D& initialValField); const realx3Field_D& initialValField);
/// Destructor /// Destructor
~AdamsBashforth2()final = default; ~AdamsBashforth2()override = default;
/// Add this to the virtual constructor table /// Add this to the virtual constructor table
add_vCtor( add_vCtor(
@ -102,32 +116,21 @@ public:
bool correct( bool correct(
real dt, real dt,
realx3PointField_D& y, realx3PointField_D& y,
realx3PointField_D& dy) final; realx3PointField_D& dy,
real damping = 1.0) override;
bool correctPStruct( bool correctPStruct(
real dt, real dt,
pointStructure& pStruct, pointStructure& pStruct,
realx3PointField_D& vel) final; realx3PointField_D& vel) override;
/*bool hearChanges /*bool hearChanges
( (
real t, const timeInfo& ti,
real dt,
uint32 iter,
const message& msg, const message& msg,
const anyList& varList const anyList& varList
) override;*/ ) override;*/
bool setInitialVals(
const int32IndexContainer& newIndices,
const realx3Vector& y) final;
bool needSetInitialVals()const final
{
return false;
}
}; };

View File

@ -20,90 +20,184 @@ Licence:
#include "AdamsBashforth3.hpp" #include "AdamsBashforth3.hpp"
//const real AB3_coef[] = { 23.0 / 12.0, 16.0 / 12.0, 5.0 / 12.0 }; namespace pFlow
pFlow::AdamsBashforth3::AdamsBashforth3
(
const word& baseName,
repository& owner,
const pointStructure& pStruct,
const word& method
)
:
integration(baseName, owner, pStruct, method),
history_(
owner.emplaceObject<pointField<VectorSingle,AB3History>>(
objectFile(
groupNames(baseName,"AB3History"),
"",
objectFile::READ_IF_PRESENT,
objectFile::WRITE_ALWAYS),
pStruct,
AB3History({zero3,zero3})))
{ {
} /// Range policy for integration kernel (alias)
using rpIntegration = Kokkos::RangePolicy<
DefaultExecutionSpace,
Kokkos::Schedule<Kokkos::Static>,
Kokkos::IndexType<uint32>
>;
bool pFlow::AdamsBashforth3::predict bool intAllActive(
(
real UNUSED(dt),
realx3Vector_D& UNUSED(y),
realx3Vector_D& UNUSED(dy)
)
{
return true;
}
bool pFlow::AdamsBashforth3::correct
(
real dt, real dt,
realx3Vector_D& y, realx3Field_D& y,
realx3Vector_D& dy realx3PointField_D& dy,
) realx3PointField_D& dy1,
realx3PointField_D& dy2,
real damping = 1.0)
{ {
if(this->pStruct().allActive()) auto d_dy = dy.deviceView();
{ auto d_y = y.deviceView();
return intAll(dt, y, dy, this->pStruct().activeRange()); auto d_dy1= dy1.deviceView();
} auto d_dy2= dy2.deviceView();
else auto activeRng = dy1.activeRange();
{
return intRange(dt, y, dy, this->pStruct().activePointsMaskD());
}
return true;
}
bool pFlow::AdamsBashforth3::setInitialVals(
const int32IndexContainer& newIndices,
const realx3Vector& y)
{
return true;
}
bool pFlow::AdamsBashforth3::intAll(
real dt,
realx3Vector_D& y,
realx3Vector_D& dy,
range activeRng)
{
auto d_dy = dy.deviceViewAll();
auto d_y = y.deviceViewAll();
auto d_history = history_.deviceViewAll();
Kokkos::parallel_for( Kokkos::parallel_for(
"AdamsBashforth3::correct", "AdamsBashforth3::correct",
rpIntegration (activeRng.first, activeRng.second), rpIntegration (activeRng.start(), activeRng.end()),
LAMBDA_HD(int32 i){ LAMBDA_HD(uint32 i){
auto ldy = d_dy[i]; d_y[i] += damping* dt*( static_cast<real>(23.0 / 12.0) * d_dy[i]
d_y[i] += dt*( static_cast<real>(23.0 / 12.0) * ldy - static_cast<real>(16.0 / 12.0) * d_dy1[i]
- static_cast<real>(16.0 / 12.0) * d_history[i].dy1_ + static_cast<real>(5.0 / 12.0) * d_dy2[i]) ;
+ static_cast<real>(5.0 / 12.0) * d_history[i].dy2_);
d_history[i] = {ldy ,d_history[i].dy1_}; d_dy2[i] = d_dy1[i];
d_dy1[i] = d_dy[i];
}); });
Kokkos::fence(); Kokkos::fence();
return true; return true;
} }
bool intScattered
(
real dt,
realx3Field_D& y,
realx3PointField_D& dy,
realx3PointField_D& dy1,
realx3PointField_D& dy2,
real damping = 1.0
)
{
auto d_dy = dy.deviceView();
auto d_y = y.deviceView();
auto d_dy1 = dy1.deviceView();
auto d_dy2 = dy2.deviceView();
auto activeRng = dy1.activeRange();
const auto& activeP = dy1.activePointsMaskDevice();
Kokkos::parallel_for(
"AdamsBashforth3::correct",
rpIntegration (activeRng.start(), activeRng.end()),
LAMBDA_HD(uint32 i){
if( activeP(i))
{
d_y[i] += damping * dt*( static_cast<real>(23.0 / 12.0) * d_dy[i]
- static_cast<real>(16.0 / 12.0) * d_dy1[i]
+ static_cast<real>(5.0 / 12.0) * d_dy2[i]);
d_dy2[i] = d_dy1[i];
d_dy1[i] = d_dy[i];
}
});
Kokkos::fence();
return true;
}
}
//const real AB3_coef[] = { 23.0 / 12.0, 16.0 / 12.0, 5.0 / 12.0 };
pFlow::AdamsBashforth3::AdamsBashforth3
(
const word& baseName,
pointStructure& pStruct,
const word& method,
const realx3Field_D& initialValField
)
:
AdamsBashforth2(baseName, pStruct, method, initialValField),
dy2_
(
objectFile
(
groupNames(baseName,"dy2"),
pStruct.time().integrationFolder(),
objectFile::READ_IF_PRESENT,
objectFile::WRITE_ALWAYS
),
pStruct,
zero3,
zero3
)
{
}
void pFlow::AdamsBashforth3::updateBoundariesSlaveToMasterIfRequested()
{
AdamsBashforth2::updateBoundariesSlaveToMasterIfRequested();
dy2_.updateBoundariesSlaveToMasterIfRequested();
}
bool pFlow::AdamsBashforth3::correct
(
real dt,
realx3PointField_D& y,
realx3PointField_D& dy,
real damping
)
{
bool success = false;
if(y.isAllActive())
{
success = intAllActive(dt, y.field(), dy, dy1(), dy2(), damping);
}
else
{
success = intScattered(dt, y.field(), dy, dy1(), dy2(), damping);
}
success = success && boundaryList().correct(dt, y, dy);
return success;
}
bool pFlow::AdamsBashforth3::correctPStruct(
real dt,
pointStructure &pStruct,
realx3PointField_D &vel)
{
bool success = false;
if(dy2().isAllActive())
{
success = intAllActive(dt, pStruct.pointPosition(), vel, dy1(), dy2());
}
else
{
success = intScattered(dt, pStruct.pointPosition(), vel, dy1(), dy2());
}
success = success && boundaryList().correctPStruct(dt, pStruct, vel);
return success;
}
/*bool pFlow::AdamsBashforth3::hearChanges
(
const timeInfo &ti,
const message &msg,
const anyList &varList
)
{
if(msg.equivalentTo(message::ITEMS_INSERT))
{
return insertValues(varList, initialValField().deviceViewAll(), dy1()) &&
insertValues(varList, initialValField().deviceViewAll(), dy2());
}
else
{
return realx3PointField_D::hearChanges(ti, msg, varList);
}
}*/

View File

@ -22,48 +22,14 @@ Licence:
#define __AdamsBashforth3_hpp__ #define __AdamsBashforth3_hpp__
#include "integration.hpp" #include "AdamsBashforth2.hpp"
#include "pointFields.hpp" #include "pointFields.hpp"
#include "boundaryIntegrationList.hpp"
namespace pFlow namespace pFlow
{ {
struct AB3History
{
TypeInfoNV("AB3History");
realx3 dy1_={0,0,0};
realx3 dy2_={0,0,0};
};
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;
}
/** /**
* Third order Adams-Bashforth integration method for solving ODE * Third order Adams-Bashforth integration method for solving ODE
@ -72,19 +38,26 @@ iOstream& operator<<(iOstream& str, const AB3History& ab3)
*/ */
class AdamsBashforth3 class AdamsBashforth3
: :
public integration public AdamsBashforth2
{ {
private:
friend class processorAB3BoundaryIntegration;
realx3PointField_D dy2_;
protected: protected:
/// Integration history const auto& dy2()const
pointField<VectorSingle,AB3History>& history_; {
return dy2_;
}
auto& dy2()
{
return dy2_;
}
/// Range policy for integration kernel
using rpIntegration = Kokkos::RangePolicy<
DefaultExecutionSpace,
Kokkos::Schedule<Kokkos::Static>,
Kokkos::IndexType<int32>
>;
public: public:
@ -96,17 +69,13 @@ public:
/// Construct from components /// Construct from components
AdamsBashforth3( AdamsBashforth3(
const word& baseName, const word& baseName,
repository& owner, pointStructure& pStruct,
const pointStructure& pStruct, const word& method,
const word& method); const realx3Field_D& initialValField);
uniquePtr<integration> clone()const override
{
return makeUnique<AdamsBashforth3>(*this);
}
/// Destructor /// Destructor
virtual ~AdamsBashforth3()=default; ~AdamsBashforth3() override =default;
/// Add this to the virtual constructor table /// Add this to the virtual constructor table
add_vCtor( add_vCtor(
@ -117,43 +86,39 @@ public:
// - Methods // - Methods
bool predict( void updateBoundariesSlaveToMasterIfRequested()override;
real UNUSED(dt),
realx3Vector_D & UNUSED(y),
realx3Vector_D& UNUSED(dy)) override;
bool correct(real dt, /// return integration method
realx3Vector_D & y, word method()const override
realx3Vector_D& dy) override;
bool setInitialVals(
const int32IndexContainer& newIndices,
const realx3Vector& y) override;
bool needSetInitialVals()const override
{ {
return false; return "AdamsBashforth3";
} }
/// Integrate on all points in the active range
bool intAll(
real dt,
realx3Vector_D& y,
realx3Vector_D& dy,
range activeRng);
/// Integrate on active points in the active range
template<typename activeFunctor> bool correct(
bool intRange(
real dt, real dt,
realx3Vector_D& y, realx3PointField_D& y,
realx3Vector_D& dy, realx3PointField_D& dy,
activeFunctor activeP ); real damping = 1.0) override;
bool correctPStruct(
real dt,
pointStructure& pStruct,
realx3PointField_D& vel) override;
/*bool hearChanges
(
const timeInfo& ti,
const message& msg,
const anyList& varList
) override;*/
}; };
template<typename activeFunctor> /*template<typename activeFunctor>
bool pFlow::AdamsBashforth3::intRange( bool pFlow::AdamsBashforth3::intRange(
real dt, real dt,
realx3Vector_D& y, realx3Vector_D& y,
@ -181,7 +146,7 @@ bool pFlow::AdamsBashforth3::intRange(
Kokkos::fence(); Kokkos::fence();
return true; return true;
} }*/
} // pFlow } // pFlow

View File

@ -20,96 +20,186 @@ Licence:
#include "AdamsBashforth4.hpp" #include "AdamsBashforth4.hpp"
namespace pFlow
pFlow::AdamsBashforth4::AdamsBashforth4
(
const word& baseName,
repository& owner,
const pointStructure& pStruct,
const word& method
)
:
integration(baseName, owner, pStruct, method),
history_(
owner.emplaceObject<pointField<VectorSingle,AB4History>>(
objectFile(
groupNames(baseName,"AB4History"),
"",
objectFile::READ_IF_PRESENT,
objectFile::WRITE_ALWAYS),
pStruct,
AB4History({zero3,zero3, zero3})))
{ {
} /// Range policy for integration kernel (alias)
using rpIntegration = Kokkos::RangePolicy<
DefaultExecutionSpace,
Kokkos::Schedule<Kokkos::Static>,
Kokkos::IndexType<uint32>
>;
bool pFlow::AdamsBashforth4::predict bool intAllActive(
(
real UNUSED(dt),
realx3Vector_D& UNUSED(y),
realx3Vector_D& UNUSED(dy)
)
{
return true;
}
bool pFlow::AdamsBashforth4::correct
(
real dt, real dt,
realx3Vector_D& y, realx3Field_D& y,
realx3Vector_D& dy realx3PointField_D& dy,
) realx3PointField_D& dy1,
realx3PointField_D& dy2,
realx3PointField_D& dy3,
real damping = 1.0)
{ {
if(this->pStruct().allActive()) auto d_dy = dy.deviceView();
{ auto d_y = y.deviceView();
return intAll(dt, y, dy, this->pStruct().activeRange()); auto d_dy1= dy1.deviceView();
} auto d_dy2= dy2.deviceView();
else auto d_dy3= dy3.deviceView();
{ auto activeRng = dy1.activeRange();
return intRange(dt, y, dy, this->pStruct().activePointsMaskD());
}
return true;
}
bool pFlow::AdamsBashforth4::setInitialVals(
const int32IndexContainer& newIndices,
const realx3Vector& y)
{
return true;
}
bool pFlow::AdamsBashforth4::intAll(
real dt,
realx3Vector_D& y,
realx3Vector_D& dy,
range activeRng)
{
auto d_dy = dy.deviceViewAll();
auto d_y = y.deviceViewAll();
auto d_history = history_.deviceViewAll();
Kokkos::parallel_for( Kokkos::parallel_for(
"AdamsBashforth4::correct", "AdamsBashforth4::correct",
rpIntegration (activeRng.first, activeRng.second), rpIntegration (activeRng.start(), activeRng.end()),
LAMBDA_HD(int32 i){ LAMBDA_HD(uint32 i){
d_y[i] += dt*( d_y[i] += damping * dt*(
static_cast<real>(55.0 / 24.0) * d_dy[i] static_cast<real>(55.0 / 24.0) * d_dy[i]
- static_cast<real>(59.0 / 24.0) * d_history[i].dy1_ - static_cast<real>(59.0 / 24.0) * d_dy1[i]
+ static_cast<real>(37.0 / 24.0) * d_history[i].dy2_ + static_cast<real>(37.0 / 24.0) * d_dy2[i]
- static_cast<real>( 9.0 / 24.0) * d_history[i].dy3_ - static_cast<real>( 9.0 / 24.0) * d_dy3[i]);
);
d_history[i].dy3_ = d_history[i].dy2_;
d_history[i].dy2_ = d_history[i].dy1_;
d_history[i].dy1_ = d_dy[i];
d_dy3[i] = d_dy2[i];
d_dy2[i] = d_dy1[i];
d_dy1[i] = d_dy[i];
}); });
Kokkos::fence(); Kokkos::fence();
return true; return true;
} }
bool intScattered
(
real dt,
realx3Field_D& y,
realx3PointField_D& dy,
realx3PointField_D& dy1,
realx3PointField_D& dy2,
realx3PointField_D& dy3,
real damping = 1.0
)
{
auto d_dy = dy.deviceView();
auto d_y = y.deviceView();
auto d_dy1 = dy1.deviceView();
auto d_dy2 = dy2.deviceView();
auto d_dy3 = dy3.deviceView();
auto activeRng = dy1.activeRange();
const auto& activeP = dy1.activePointsMaskDevice();
Kokkos::parallel_for(
"AdamsBashforth4::correct",
rpIntegration (activeRng.start(), activeRng.end()),
LAMBDA_HD(uint32 i){
if( activeP(i))
{
d_y[i] += damping* dt*(
static_cast<real>(55.0 / 24.0) * d_dy[i]
- static_cast<real>(59.0 / 24.0) * d_dy1[i]
+ static_cast<real>(37.0 / 24.0) * d_dy2[i]
- static_cast<real>( 9.0 / 24.0) * d_dy3[i] );
d_dy3[i] = d_dy2[i];
d_dy2[i] = d_dy1[i];
d_dy1[i] = d_dy[i];
}
});
Kokkos::fence();
return true;
}
}
pFlow::AdamsBashforth4::AdamsBashforth4
(
const word& baseName,
pointStructure& pStruct,
const word& method,
const realx3Field_D& initialValField
)
:
AdamsBashforth3(baseName, pStruct, method, initialValField),
dy3_
(
objectFile
(
groupNames(baseName,"dy3"),
pStruct.time().integrationFolder(),
objectFile::READ_IF_PRESENT,
objectFile::WRITE_ALWAYS
),
pStruct,
zero3,
zero3
)
{}
void pFlow::AdamsBashforth4::updateBoundariesSlaveToMasterIfRequested()
{
AdamsBashforth3::updateBoundariesSlaveToMasterIfRequested();
dy3_.updateBoundariesSlaveToMasterIfRequested();
}
bool pFlow::AdamsBashforth4::correct
(
real dt,
realx3PointField_D& y,
realx3PointField_D& dy,
real damping
)
{
bool success = false;
if(y.isAllActive())
{
success = intAllActive(dt, y.field(), dy, dy1(), dy2(), dy3(), damping);
}
else
{
success = intScattered(dt, y.field(), dy, dy1(), dy2(), dy3(), damping);
}
success = success && boundaryList().correct(dt, y, dy);
return success;
}
bool pFlow::AdamsBashforth4::correctPStruct(
real dt,
pointStructure &pStruct,
realx3PointField_D &vel)
{
bool success = false;
if(dy2().isAllActive())
{
success = intAllActive(dt, pStruct.pointPosition(), vel, dy1(), dy2(), dy3());
}
else
{
success = intScattered(dt, pStruct.pointPosition(), vel, dy1(), dy2(), dy3());
}
success = success && boundaryList().correctPStruct(dt, pStruct, vel);
return success;
}
/*bool pFlow::AdamsBashforth4::hearChanges
(
const timeInfo &ti,
const message &msg,
const anyList &varList
)
{
if(msg.equivalentTo(message::ITEMS_INSERT))
{
return insertValues(varList, initialValField().deviceViewAll(), dy1()) &&
insertValues(varList, initialValField().deviceViewAll(), dy2()) &&
insertValues(varList, initialValField().deviceViewAll(), dy3());
}
else
{
return realx3PointField_D::hearChanges(ti, msg, varList);
}
}*/

View File

@ -22,53 +22,14 @@ Licence:
#define __AdamsBashforth4_hpp__ #define __AdamsBashforth4_hpp__
#include "integration.hpp" #include "AdamsBashforth3.hpp"
#include "pointFields.hpp" #include "pointFields.hpp"
#include "boundaryIntegrationList.hpp"
namespace pFlow namespace pFlow
{ {
struct AB4History
{
TypeInfoNV("AB4History");
realx3 dy1_={0,0,0};
realx3 dy2_={0,0,0};
realx3 dy3_={0,0,0};
};
INLINE_FUNCTION
iIstream& operator>>(iIstream& str, AB4History& ab4)
{
str.readBegin("AB4History");
str >> ab4.dy1_;
str >> ab4.dy2_;
str >> ab4.dy3_;
str.readEnd("AB4History");
str.check(FUNCTION_NAME);
return str;
}
INLINE_FUNCTION
iOstream& operator<<(iOstream& str, const AB4History& ab4)
{
str << token::BEGIN_LIST << ab4.dy1_
<< token::SPACE << ab4.dy2_
<< token::SPACE << ab4.dy3_
<< token::END_LIST;
str.check(FUNCTION_NAME);
return str;
}
/** /**
* Fourth order Adams-Bashforth integration method for solving ODE * Fourth order Adams-Bashforth integration method for solving ODE
* *
@ -76,19 +37,25 @@ iOstream& operator<<(iOstream& str, const AB4History& ab4)
*/ */
class AdamsBashforth4 class AdamsBashforth4
: :
public integration public AdamsBashforth3
{ {
private:
friend class processorAB4BoundaryIntegration;
realx3PointField_D dy3_;
protected: protected:
/// Integration history const auto& dy3()const
pointField<VectorSingle,AB4History>& history_; {
return dy3_;
}
/// Range policy for integration kernel auto& dy3()
using rpIntegration = Kokkos::RangePolicy< {
DefaultExecutionSpace, return dy3_;
Kokkos::Schedule<Kokkos::Static>, }
Kokkos::IndexType<int32>
>;
public: public:
@ -100,17 +67,14 @@ public:
/// Construct from components /// Construct from components
AdamsBashforth4( AdamsBashforth4(
const word& baseName, const word& baseName,
repository& owner, pointStructure& pStruct,
const pointStructure& pStruct, const word& method,
const word& method); const realx3Field_D& initialValField);
uniquePtr<integration> clone()const override
{
return makeUnique<AdamsBashforth4>(*this);
}
/// Destructor /// Destructor
virtual ~AdamsBashforth4()=default; ~AdamsBashforth4() override =default;
/// Add a this to the virtual constructor table /// Add a this to the virtual constructor table
add_vCtor( add_vCtor(
@ -121,77 +85,36 @@ public:
// - Methods // - Methods
bool predict( void updateBoundariesSlaveToMasterIfRequested()override;
real UNUSED(dt),
realx3Vector_D & UNUSED(y), /// return integration method
realx3Vector_D& UNUSED(dy)) override; word method()const override
{
return "AdamsBashforth4";
}
bool correct( bool correct(
real dt, real dt,
realx3Vector_D & y, realx3PointField_D& y,
realx3Vector_D& dy) override; realx3PointField_D& dy,
real damping = 1.0) override;
bool setInitialVals( bool correctPStruct(
const int32IndexContainer& newIndices,
const realx3Vector& y) override;
bool needSetInitialVals()const override
{
return false;
}
/// Integrate on all points in the active range
bool intAll(
real dt, real dt,
realx3Vector_D& y, pointStructure& pStruct,
realx3Vector_D& dy, realx3PointField_D& vel) override;
range activeRng);
/// Integrate on active points in the active range /*bool hearChanges
template<typename activeFunctor> (
bool intRange( const timeInfo& ti,
real dt, const message& msg,
realx3Vector_D& y, const anyList& varList
realx3Vector_D& dy, ) override;*/
activeFunctor activeP );
}; };
template<typename activeFunctor>
bool pFlow::AdamsBashforth4::intRange(
real dt,
realx3Vector_D& y,
realx3Vector_D& dy,
activeFunctor activeP )
{
auto d_dy = dy.deviceViewAll();
auto d_y = y.deviceViewAll();
auto d_history = history_.deviceViewAll();
auto activeRng = activeP.activeRange();
Kokkos::parallel_for(
"AdamsBashforth4::correct",
rpIntegration (activeRng.first, activeRng.second),
LAMBDA_HD(int32 i){
if( activeP(i))
{
d_y[i] += dt*(
static_cast<real>(55.0 / 24.0) * d_dy[i]
- static_cast<real>(59.0 / 24.0) * d_history[i].dy1_
+ static_cast<real>(37.0 / 24.0) * d_history[i].dy2_
- static_cast<real>( 9.0 / 24.0) * d_history[i].dy3_
);
d_history[i].dy3_ = d_history[i].dy2_;
d_history[i].dy2_ = d_history[i].dy1_;
d_history[i].dy1_ = d_dy[i];
}
});
Kokkos::fence();
return true;
}
} // pFlow } // pFlow

View File

@ -20,93 +20,178 @@ Licence:
#include "AdamsBashforth5.hpp" #include "AdamsBashforth5.hpp"
namespace pFlow
pFlow::AdamsBashforth5::AdamsBashforth5
(
const word& baseName,
repository& owner,
const pointStructure& pStruct,
const word& method
)
:
integration(baseName, owner, pStruct, method),
history_(
owner.emplaceObject<pointField<VectorSingle,AB5History>>(
objectFile(
groupNames(baseName,"AB5History"),
"",
objectFile::READ_IF_PRESENT,
objectFile::WRITE_ALWAYS),
pStruct,
AB5History({zero3,zero3, zero3})))
{ {
} /// Range policy for integration kernel (alias)
using rpIntegration = Kokkos::RangePolicy<
DefaultExecutionSpace,
Kokkos::Schedule<Kokkos::Static>,
Kokkos::IndexType<uint32>
>;
bool pFlow::AdamsBashforth5::predict bool intAllActive(
(
real UNUSED(dt),
realx3Vector_D& UNUSED(y),
realx3Vector_D& UNUSED(dy)
)
{
return true;
}
bool pFlow::AdamsBashforth5::correct
(
real dt, real dt,
realx3Vector_D& y, realx3Field_D& y,
realx3Vector_D& dy realx3PointField_D& dy,
) realx3PointField_D& dy1,
realx3PointField_D& dy2,
realx3PointField_D& dy3,
realx3PointField_D& dy4,
real damping = 1.0)
{ {
if(this->pStruct().allActive()) auto d_dy = dy.deviceView();
{ auto d_y = y.deviceView();
return intAll(dt, y, dy, this->pStruct().activeRange()); auto d_dy1= dy1.deviceView();
} auto d_dy2= dy2.deviceView();
else auto d_dy3= dy3.deviceView();
{ auto d_dy4= dy4.deviceView();
return intRange(dt, y, dy, this->pStruct().activePointsMaskD()); auto activeRng = dy1.activeRange();
}
return true;
}
bool pFlow::AdamsBashforth5::setInitialVals(
const int32IndexContainer& newIndices,
const realx3Vector& y)
{
return true;
}
bool pFlow::AdamsBashforth5::intAll(
real dt,
realx3Vector_D& y,
realx3Vector_D& dy,
range activeRng)
{
auto d_dy = dy.deviceViewAll();
auto d_y = y.deviceViewAll();
auto d_history = history_.deviceViewAll();
Kokkos::parallel_for( Kokkos::parallel_for(
"AdamsBashforth5::correct", "AdamsBashforth5::correct",
rpIntegration (activeRng.first, activeRng.second), rpIntegration (activeRng.start(), activeRng.end()),
LAMBDA_HD(int32 i){ LAMBDA_HD(uint32 i){
d_y[i] += dt*( d_y[i] += damping* dt*(
static_cast<real>(1901.0 / 720.0) * d_dy[i] static_cast<real>(1901.0 / 720.0) * d_dy[i]
- static_cast<real>(2774.0 / 720.0) * d_history[i].dy1_ - static_cast<real>(2774.0 / 720.0) * d_dy1[i]
+ static_cast<real>(2616.0 / 720.0) * d_history[i].dy2_ + static_cast<real>(2616.0 / 720.0) * d_dy2[i]
- static_cast<real>(1274.0 / 720.0) * d_history[i].dy3_ - static_cast<real>(1274.0 / 720.0) * d_dy3[i]
+ static_cast<real>( 251.0 / 720.0) * d_history[i].dy4_ + static_cast<real>( 251.0 / 720.0) * d_dy4[i]);
);
d_history[i] = {d_dy[i] ,d_history[i].dy1_, d_history[i].dy2_, d_history[i].dy3_}; d_dy4[i] = d_dy3[i];
d_dy3[i] = d_dy2[i];
d_dy2[i] = d_dy1[i];
d_dy1[i] = d_dy[i];
}); });
Kokkos::fence(); Kokkos::fence();
return true; return true;
} }
bool intScattered
(
real dt,
realx3Field_D& y,
realx3PointField_D& dy,
realx3PointField_D& dy1,
realx3PointField_D& dy2,
realx3PointField_D& dy3,
realx3PointField_D& dy4,
real damping = 1.0
)
{
auto d_dy = dy.deviceView();
auto d_y = y.deviceView();
auto d_dy1 = dy1.deviceView();
auto d_dy2 = dy2.deviceView();
auto d_dy3 = dy3.deviceView();
auto d_dy4 = dy4.deviceView();
auto activeRng = dy1.activeRange();
const auto& activeP = dy1.activePointsMaskDevice();
Kokkos::parallel_for(
"AdamsBashforth5::correct",
rpIntegration (activeRng.start(), activeRng.end()),
LAMBDA_HD(uint32 i){
if( activeP(i))
{
d_y[i] += damping* dt*(
static_cast<real>(1901.0 / 720.0) * d_dy[i]
- static_cast<real>(2774.0 / 720.0) * d_dy1[i]
+ static_cast<real>(2616.0 / 720.0) * d_dy2[i]
- static_cast<real>(1274.0 / 720.0) * d_dy3[i]
+ static_cast<real>( 251.0 / 720.0) * d_dy4[i]);
d_dy4[i] = d_dy3[i];
d_dy3[i] = d_dy2[i];
d_dy2[i] = d_dy1[i];
d_dy1[i] = d_dy[i];
}
});
Kokkos::fence();
return true;
}
}
pFlow::AdamsBashforth5::AdamsBashforth5
(
const word &baseName,
pointStructure &pStruct,
const word &method,
const realx3Field_D &initialValField
)
:
AdamsBashforth4(baseName, pStruct, method, initialValField),
dy4_
(
objectFile
(
groupNames(baseName,"dy4"),
pStruct.time().integrationFolder(),
objectFile::READ_IF_PRESENT,
objectFile::WRITE_ALWAYS
),
pStruct,
zero3,
zero3
)
{}
void pFlow::AdamsBashforth5::updateBoundariesSlaveToMasterIfRequested()
{
AdamsBashforth4::updateBoundariesSlaveToMasterIfRequested();
dy4_.updateBoundariesSlaveToMasterIfRequested();
}
bool pFlow::AdamsBashforth5::correct
(
real dt,
realx3PointField_D &y,
realx3PointField_D &dy,
real damping
)
{
bool success = false;
if(y.isAllActive())
{
success = intAllActive(dt, y.field(), dy, dy1(), dy2(), dy3(), dy4(), damping);
}
else
{
success = intScattered(dt, y.field(), dy, dy1(), dy2(), dy3(), dy4(), damping);
}
success = success && boundaryList().correct(dt, y, dy);
return success;
}
bool pFlow::AdamsBashforth5::correctPStruct
(
real dt,
pointStructure &pStruct,
realx3PointField_D &vel
)
{
bool success = false;
if(dy2().isAllActive())
{
success = intAllActive(dt, pStruct.pointPosition(), vel, dy1(), dy2(), dy3(), dy4());
}
else
{
success = intScattered(dt, pStruct.pointPosition(), vel, dy1(), dy2(), dy3(), dy4());
}
success = success && boundaryList().correctPStruct(dt, pStruct, vel);
return success;
}

View File

@ -22,54 +22,12 @@ Licence:
#define __AdamsBashforth5_hpp__ #define __AdamsBashforth5_hpp__
#include "integration.hpp" #include "AdamsBashforth4.hpp"
#include "pointFields.hpp" #include "pointFields.hpp"
namespace pFlow namespace pFlow
{ {
struct AB5History
{
TypeInfoNV("AB5History");
realx3 dy1_={0,0,0};
realx3 dy2_={0,0,0};
realx3 dy3_={0,0,0};
realx3 dy4_={0,0,0};
};
INLINE_FUNCTION
iIstream& operator>>(iIstream& str, AB5History& ab5)
{
str.readBegin("AB5History");
str >> ab5.dy1_;
str >> ab5.dy2_;
str >> ab5.dy3_;
str >> ab5.dy4_;
str.readEnd("AB5History");
str.check(FUNCTION_NAME);
return str;
}
INLINE_FUNCTION
iOstream& operator<<(iOstream& str, const AB5History& ab5)
{
str << token::BEGIN_LIST << ab5.dy1_
<< token::SPACE << ab5.dy2_
<< token::SPACE << ab5.dy3_
<< token::SPACE << ab5.dy4_
<< token::END_LIST;
str.check(FUNCTION_NAME);
return str;
}
/** /**
* Fifth order Adams-Bashforth integration method for solving ODE * Fifth order Adams-Bashforth integration method for solving ODE
@ -78,19 +36,25 @@ iOstream& operator<<(iOstream& str, const AB5History& ab5)
*/ */
class AdamsBashforth5 class AdamsBashforth5
: :
public integration public AdamsBashforth4
{ {
private:
friend class processorAB4BoundaryIntegration;
realx3PointField_D dy4_;
protected: protected:
/// Integration history const auto& dy4()const
pointField<VectorSingle,AB5History>& history_; {
return dy4_;
}
/// Range policy for integration kernel auto& dy4()
using rpIntegration = Kokkos::RangePolicy< {
DefaultExecutionSpace, return dy4_;
Kokkos::Schedule<Kokkos::Static>, }
Kokkos::IndexType<int32>
>;
public: public:
@ -99,20 +63,19 @@ public:
// - Constructors // - Constructors
/// Construct from components
AdamsBashforth5( AdamsBashforth5(
const word& baseName, const word& baseName,
repository& owner, pointStructure& pStruct,
const pointStructure& pStruct, const word& method,
const word& method); const realx3Field_D& initialValField);
uniquePtr<integration> clone()const override
{
return makeUnique<AdamsBashforth5>(*this);
}
virtual ~AdamsBashforth5()=default;
/// Add this to the virtual constructor table /// Destructor
~AdamsBashforth5() override =default;
/// Add a this to the virtual constructor table
add_vCtor( add_vCtor(
integration, integration,
AdamsBashforth5, AdamsBashforth5,
@ -121,44 +84,29 @@ public:
// - Methods // - Methods
bool predict( void updateBoundariesSlaveToMasterIfRequested()override;
real UNUSED(dt),
realx3Vector_D & UNUSED(y), /// return integration method
realx3Vector_D& UNUSED(dy)) override; word method()const override
{
return "AdamsBashforth5";
}
bool correct( bool correct(
real dt, real dt,
realx3Vector_D & y, realx3PointField_D& y,
realx3Vector_D& dy) override; realx3PointField_D& dy,
real damping = 1.0) override;
bool setInitialVals( bool correctPStruct(
const int32IndexContainer& newIndices,
const realx3Vector& y) override;
bool needSetInitialVals()const override
{
return false;
}
/// Integrate on all points in the active range
bool intAll(
real dt, real dt,
realx3Vector_D& y, pointStructure& pStruct,
realx3Vector_D& dy, realx3PointField_D& vel) override;
range activeRng);
/// Integrate on active points in the active range
template<typename activeFunctor>
bool intRange(
real dt,
realx3Vector_D& y,
realx3Vector_D& dy,
activeFunctor activeP );
}; };
template<typename activeFunctor> /*template<typename activeFunctor>
bool pFlow::AdamsBashforth5::intRange( bool pFlow::AdamsBashforth5::intRange(
real dt, real dt,
realx3Vector_D& y, realx3Vector_D& y,
@ -190,7 +138,7 @@ bool pFlow::AdamsBashforth5::intRange(
Kokkos::fence(); Kokkos::fence();
return true; return true;
} }*/
} // pFlow } // pFlow

View File

@ -4,9 +4,10 @@ integration/integration.cpp
boundaries/boundaryIntegration.cpp boundaries/boundaryIntegration.cpp
boundaries/boundaryIntegrationList.cpp boundaries/boundaryIntegrationList.cpp
AdamsBashforth2/AdamsBashforth2.cpp AdamsBashforth2/AdamsBashforth2.cpp
#AdamsBashforth5/AdamsBashforth5.cpp AdamsBashforth3/AdamsBashforth3.cpp
#AdamsBashforth4/AdamsBashforth4.cpp AdamsBashforth4/AdamsBashforth4.cpp
#AdamsBashforth3/AdamsBashforth3.cpp AdamsBashforth5/AdamsBashforth5.cpp
#AdamsMoulton3/AdamsMoulton3.cpp #AdamsMoulton3/AdamsMoulton3.cpp
#AdamsMoulton4/AdamsMoulton4.cpp #AdamsMoulton4/AdamsMoulton4.cpp
#AdamsMoulton5/AdamsMoulton5.cpp #AdamsMoulton5/AdamsMoulton5.cpp

View File

@ -51,9 +51,7 @@ public:
); );
bool hearChanges( bool hearChanges(
real t, const timeInfo& ti,
real dt,
uint32 iter,
const message &msg, const message &msg,
const anyList &varList) override const anyList &varList) override
{ {
@ -77,6 +75,11 @@ public:
return true; return true;
} }
bool isActive()const override
{
return false;
}
static static
uniquePtr<boundaryIntegration> create( uniquePtr<boundaryIntegration> create(
const boundaryBase& boundary, const boundaryBase& boundary,

View File

@ -7,11 +7,11 @@ pFlow::boundaryIntegrationList::boundaryIntegrationList(
integration &intgrtn integration &intgrtn
) )
: :
ListPtr<boundaryIntegration>(6), boundaryListPtr<boundaryIntegration>(),
boundaries_(pStruct.boundaries()) boundaries_(pStruct.boundaries())
{ {
for(uint32 i=0; i<6; i++) ForAllBoundariesPtr(i, this)
{ {
this->set( this->set(
i, i,
@ -19,25 +19,23 @@ pFlow::boundaryIntegrationList::boundaryIntegrationList(
boundaries_[i], boundaries_[i],
pStruct, pStruct,
method, method,
intgrtn intgrtn
) )
); );
} }
} }
bool pFlow::boundaryIntegrationList::correct(real dt, realx3PointField_D &y, realx3PointField_D &dy) bool pFlow::boundaryIntegrationList::correct(real dt, realx3PointField_D &y, realx3PointField_D &dy)
{ {
for(auto& bndry:*this) ForAllActiveBoundariesPtr(i,this)
{ {
if(!bndry->correct(dt, y, dy)) if(!boundaryPtr(i)->correct(dt, y, dy))
{ {
fatalErrorInFunction<<"Error in correcting boundary "<< fatalErrorInFunction<<"Error in correcting boundary "<<
bndry->boundaryName()<<endl; boundaryPtr(i)->boundaryName()<<endl;
return false; return false;
} }
}
}
return true; return true;
} }
@ -46,14 +44,15 @@ bool pFlow::boundaryIntegrationList::correctPStruct(
pointStructure &pStruct, pointStructure &pStruct,
const realx3PointField_D &vel) const realx3PointField_D &vel)
{ {
for(auto& bndry:*this) ForAllActiveBoundariesPtr(i,this)
{ {
if(!bndry->correctPStruct(dt, vel)) if(!boundaryPtr(i)->correctPStruct(dt, vel))
{ {
fatalErrorInFunction<<"Error in correcting boundary "<< fatalErrorInFunction<<"Error in correcting boundary "<<
bndry->boundaryName()<<" in pointStructure."<<endl; boundaryPtr(i)->boundaryName()<<" in pointStructure."<<endl;
return false; return false;
} }
} }
return true; return true;
} }

View File

@ -4,7 +4,7 @@
#include "boundaryList.hpp" #include "boundaryList.hpp"
#include "ListPtr.hpp" #include "boundaryListPtr.hpp"
#include "boundaryIntegration.hpp" #include "boundaryIntegration.hpp"
@ -15,7 +15,7 @@ class integration;
class boundaryIntegrationList class boundaryIntegrationList
: :
public ListPtr<boundaryIntegration> public boundaryListPtr<boundaryIntegration>
{ {
private: private:

View File

@ -22,17 +22,39 @@ Licence:
#include "pointStructure.hpp" #include "pointStructure.hpp"
#include "repository.hpp" #include "repository.hpp"
pFlow::integration::integration bool pFlow::integration::insertValues
( (
const word& baseName, const anyList& varList,
pointStructure& pStruct, const deviceViewType1D<realx3>& srcVals,
const word&, realx3PointField_D& dstFeild
const realx3Field_D&
) )
: {
owner_(*pStruct.owner()), const word eventName = message::eventName(message::ITEMS_INSERT);
pStruct_(pStruct),
baseName_(baseName) if(!varList.checkObjectType<uint32IndexContainer>(eventName))
{
fatalErrorInFunction<<"Insertion failed in integration for "<< baseName_<<
", variable "<< eventName <<
" does not exist or the type is incorrect"<<endl;
return false;
}
const auto& indices = varList.getObject<uint32IndexContainer>(
eventName);
dstFeild.field().insertSetElement(indices, srcVals);
return true;
}
pFlow::integration::integration(
const word &baseName,
pointStructure &pStruct,
const word &,
const realx3Field_D &)
: owner_(*pStruct.owner()),
pStruct_(pStruct),
baseName_(baseName)
{} {}

View File

@ -63,6 +63,13 @@ private:
/// The base name for integration /// The base name for integration
const word baseName_; const word baseName_;
protected:
bool insertValues(
const anyList& varList,
const deviceViewType1D<realx3>& srcVals,
realx3PointField_D& dstFeild);
public: public:
/// Type info /// Type info
@ -146,22 +153,11 @@ public:
/// Correction/main integration step /// Correction/main integration step
virtual virtual
bool correct(real dt, realx3PointField_D& y, realx3PointField_D& dy) = 0; bool correct(real dt, realx3PointField_D& y, realx3PointField_D& dy, real damping = 1.0) = 0;
virtual virtual
bool correctPStruct(real dt, pointStructure& pStruct, realx3PointField_D& vel) = 0; bool correctPStruct(real dt, pointStructure& pStruct, realx3PointField_D& vel) = 0;
/// Set the initial values for new indices
virtual
bool setInitialVals(
const int32IndexContainer& newIndices,
const realx3Vector& y) = 0;
/// Check if the method requires any set initial vals
virtual
bool needSetInitialVals()const = 0;
/// Create the polymorphic object based on inputs /// Create the polymorphic object based on inputs
static static
uniquePtr<integration> create( uniquePtr<integration> create(

View File

@ -25,12 +25,6 @@ Licence:
#include "symArrays.hpp" #include "symArrays.hpp"
namespace pFlow::cfModels namespace pFlow::cfModels
{ {
@ -159,7 +153,7 @@ protected:
{ {
etha_t[i] = -2.0*log( et[i])*sqrt(kt[i]*2/7) / etha_t[i] = -2.0*log( et[i])*sqrt(kt[i]*2/7) /
sqrt(log(pow(et[i],2.0))+ pow(Pi,2.0)); sqrt(log(pow(et[i],2))+ pow(Pi,2));
} }
Vector<linearProperties> prop("prop", nElem); Vector<linearProperties> prop("prop", nElem);
@ -285,8 +279,8 @@ public:
history.overlap_t_ += Vt*dt; history.overlap_t_ += Vt*dt;
real mi = 3*Pi/4*pow(Ri,3.0)*rho_[propId_i]; real mi = 3*Pi/4*pow(Ri,3)*rho_[propId_i];
real mj = 3*Pi/4*pow(Rj,3.0)*rho_[propId_j]; real mj = 3*Pi/4*pow(Rj,3)*rho_[propId_j];
real sqrt_meff = sqrt((mi*mj)/(mi+mj)); real sqrt_meff = sqrt((mi*mj)/(mi+mj));
@ -299,29 +293,24 @@ public:
else if (addDissipationModel_==3) else if (addDissipationModel_==3)
{ {
auto pie =3.14; auto pie =3.14;
prop.en_ = exp((pow(f_,1.5)*log(prop.en_)*sqrt( (1-((pow(log(prop.en_),2))/(pow(log(prop.en_),2)+pow(pie,2))))/(1-(pow(f_,3)*(pow(log(prop.en_),2))/(pow(log(prop.en_),2)+pow(pie,2)))) ) )); prop.en_ = exp((pow(f_,static_cast<real>(1.5))*log(prop.en_)*sqrt( (1-((pow(log(prop.en_),2))/(pow(log(prop.en_),2)+pow(pie,2))))/(1-(pow(f_,3)*(pow(log(prop.en_),2))/(pow(log(prop.en_),2)+pow(pie,2)))) ) ));
} }
real ethan_ = -2.0*log(prop.en_)*sqrt(prop.kn_)/ real ethan_ = -2.0*log(prop.en_)*sqrt(prop.kn_)/
sqrt(pow(log(prop.en_),2.0)+ pow(Pi,2.0)); sqrt(pow(log(prop.en_),2)+ pow(Pi,2));
//REPORT(0)<<"\n en n is : "<<END_REPORT; //REPORT(0)<<"\n en n is : "<<END_REPORT;
//REPORT(0)<< prop.en_ <<END_REPORT; //REPORT(0)<< prop.en_ <<END_REPORT;
FCn = ( -pow(f_,3.0)*prop.kn_ * ovrlp_n - sqrt_meff * pow(f_,1.5) * ethan_ * vrn)*Nij; FCn = ( -pow(f_,3)*prop.kn_ * ovrlp_n - sqrt_meff * pow(f_,static_cast<real>(1.5)) * ethan_ * vrn)*Nij;
FCt = ( -pow(f_,3.0)*prop.kt_ * history.overlap_t_ - sqrt_meff * pow(f_,1.5) * prop.ethat_*Vt); FCt = ( -pow(f_,3)*prop.kt_ * history.overlap_t_ - sqrt_meff * pow(f_,static_cast<real>(1.5)) * prop.ethat_*Vt);
real ft = length(FCt); real ft = length(FCt);
real ft_fric = prop.mu_ * length(FCn); real ft_fric = prop.mu_ * length(FCn);
if(ft > ft_fric) if(ft > ft_fric)
{ {
if( length(history.overlap_t_) >static_cast<real>(0.0)) if( length(history.overlap_t_) >zero)
{ {
if constexpr (limited) if constexpr (limited)
{ {

View File

@ -0,0 +1,307 @@
/*------------------------------- 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 __cGNonLinearCF_hpp__
#define __cGNonLinearCF_hpp__
#include "types.hpp"
namespace pFlow::cfModels
{
template<bool limited=true>
class cGNonLinear
{
public:
struct contactForceStorage
{
realx3 overlap_t_ = 0.0;
};
struct cGNonLinearProperties
{
real Yeff_ = 1000000.0;
real Geff_ = 8000000.0;
real en_ = 1.0;
real mu_ = 0.00001;
INLINE_FUNCTION_HD
cGNonLinearProperties(){}
INLINE_FUNCTION_HD
cGNonLinearProperties(real Yeff, real Geff, real en, real mu ):
Yeff_(Yeff), Geff_(Geff), en_(en), mu_(mu)
{}
INLINE_FUNCTION_HD
cGNonLinearProperties(const cGNonLinearProperties&)=default;
INLINE_FUNCTION_HD
cGNonLinearProperties& operator=(const cGNonLinearProperties&)=default;
INLINE_FUNCTION_HD
~cGNonLinearProperties() = default;
};
protected:
using cGNonLinearArrayType = symArray<cGNonLinearProperties>;
int32 numMaterial_ = 0;
ViewType1D<real> rho_;
int32 addDissipationModel_;
cGNonLinearArrayType nonlinearProperties_;
bool readNonLinearDictionary(const dictionary& dict)
{
auto Yeff = dict.getVal<realVector>("Yeff");
auto Geff = dict.getVal<realVector>("Geff");
auto nu = dict.getVal<realVector>("nu");
auto en = dict.getVal<realVector>("en");
auto mu = dict.getVal<realVector>("mu");
auto nElem = Yeff.size();
if(nElem != nu.size())
{
fatalErrorInFunction<<
"sizes of Yeff("<<nElem<<") and nu("<<nu.size()<<") do not match.\n";
return false;
}
if(nElem != en.size())
{
fatalErrorInFunction<<
"sizes of Yeff("<<nElem<<") and en("<<en.size()<<") do not match.\n";
return false;
}
if(nElem != mu.size())
{
fatalErrorInFunction<<
"sizes of Yeff("<<nElem<<") and mu("<<mu.size()<<") do not match.\n";
return false;
}
// check if size of vector matchs a symetric array
uint32 nMat;
if( !cGNonLinearArrayType::getN(nElem, nMat) )
{
fatalErrorInFunction<<
"sizes of properties do not match a symetric array with size ("<<
numMaterial_<<"x"<<numMaterial_<<").\n";
return false;
}
else if( numMaterial_ != nMat)
{
fatalErrorInFunction<<
"size mismatch for porperties. \n";
return false;
}
Vector<cGNonLinearProperties> prop("prop",nElem);
ForAll(i,Yeff)
{
prop[i] = {Yeff[i], Geff[i], en[i], mu[i]};
}
nonlinearProperties_.assign(prop);
auto adm = dict.getVal<word>("additionalDissipationModel");
if(adm == "none")
{
addDissipationModel_ = 1;
}
else if(adm == "LU")
{
addDissipationModel_ = 2;
}
else if (adm == "GB")
{
addDissipationModel_ = 3;
}
else
{
addDissipationModel_ = 1;
}
return true;
}
static const char* modelName()
{
if constexpr (limited)
{
return "cGNonLinearLimited";
}
else
{
return "cGNonLinearNonLimited";
}
return "";
}
public:
TypeInfoNV(modelName());
INLINE_FUNCTION_HD
cGNonLinear(){}
cGNonLinear(
int32 nMaterial,
const ViewType1D<real>& rho,
const dictionary& dict)
:
numMaterial_(nMaterial),
rho_("rho",nMaterial),
nonlinearProperties_("cGNonLinearProperties",nMaterial)
{
Kokkos::deep_copy(rho_,rho);
if(!readNonLinearDictionary(dict))
{
fatalExit;
}
}
INLINE_FUNCTION_HD
cGNonLinear(const cGNonLinear&) = default;
INLINE_FUNCTION_HD
cGNonLinear(cGNonLinear&&) = default;
INLINE_FUNCTION_HD
cGNonLinear& operator=(const cGNonLinear&) = default;
INLINE_FUNCTION_HD
cGNonLinear& operator=(cGNonLinear&&) = default;
INLINE_FUNCTION_HD
~cGNonLinear()=default;
INLINE_FUNCTION_HD
int32 numMaterial()const
{
return numMaterial_;
}
//// - Methods
INLINE_FUNCTION_HD
void contactForce
(
const real dt,
const uint32 i,
const uint32 j,
const uint32 propId_i,
const uint32 propId_j,
const real Ri,
const real Rj,
const real cGFi,
const real cGFj,
const real ovrlp_n,
const realx3& Vr,
const realx3& Nij,
contactForceStorage& history,
realx3& FCn,
realx3& FCt
)const
{
const auto prop = nonlinearProperties_(propId_i,propId_j);
const real f = 2/( 1/cGFi + 1/cGFj );
real vrn = dot(Vr, Nij);
realx3 Vt = Vr - vrn*Nij;
history.overlap_t_ += Vt*dt;
real mi = 3*Pi/4*pow(Ri,static_cast<real>(3))*rho_[propId_i];
real mj = 3*Pi/4*pow(Rj,static_cast<real>(3))*rho_[propId_j];
real Reff = 1/(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);
real en = prop.en_;
if (addDissipationModel_==2)
{
en = sqrt(1+((pow(prop.en_,2)-1)*f));
}
else if (addDissipationModel_==3)
{
en = exp((pow(f,static_cast<real>(1.5))*log(prop.en_)*sqrt( (1-((pow(log(prop.en_),2))/(pow(log(prop.en_),2)+pow(Pi,2))))/(1-(pow(f,3)*(pow(log(prop.en_),2))/(pow(log(prop.en_),2)+pow(Pi,2)))) ) ));
}
real Kn = static_cast<real>(4.0/3.0) * prop.Yeff_ * sqrt(Reff*ovrlp_n);
real ethan_ = -2.0*log(en)*sqrt(Kn)/
sqrt(pow(log(en),2)+ pow(Pi,2));
FCn = ( - Kn*ovrlp_n -
sqrt_meff_K_hertz*ethan_*pow(ovrlp_n,static_cast<real>(0.25))*vrn)*Nij;
FCt = (f*f)*(- static_cast<real>(8.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<real>(8.0) * prop.Geff_ * sqrt(Reff*ovrlp_n);
FCt *= (ft_fric/ft);
history.overlap_t_ = - FCt/(f*f*kt);
}
else
{
FCt = (FCt/ft)*ft_fric;
}
}
else
{
FCt = 0.0;
}
}
}
};
} //pFlow::CFModels
#endif

View File

@ -26,11 +26,6 @@ Licence:
namespace pFlow::cfModels namespace pFlow::cfModels
{ {
@ -49,16 +44,15 @@ public:
{ {
real kn_ = 1000.0; real kn_ = 1000.0;
real kt_ = 800.0; real kt_ = 800.0;
real en_ = 0.0; real en_ = 1.0;
real ethat_ = 0.0; real mu_ = 0.00001;
real mu_ = 0.00001;
INLINE_FUNCTION_HD INLINE_FUNCTION_HD
linearProperties(){} linearProperties(){}
INLINE_FUNCTION_HD INLINE_FUNCTION_HD
linearProperties(real kn, real kt, real en, real etha_t, real mu ): linearProperties(real kn, real kt, real en, real mu ):
kn_(kn), kt_(kt), en_(en),ethat_(etha_t), mu_(mu) kn_(kn), kt_(kt), en_(en), mu_(mu)
{} {}
INLINE_FUNCTION_HD INLINE_FUNCTION_HD
@ -93,7 +87,6 @@ protected:
auto kn = dict.getVal<realVector>("kn"); auto kn = dict.getVal<realVector>("kn");
auto kt = dict.getVal<realVector>("kt"); auto kt = dict.getVal<realVector>("kt");
auto en = dict.getVal<realVector>("en"); auto en = dict.getVal<realVector>("en");
auto et = dict.getVal<realVector>("et");
auto mu = dict.getVal<realVector>("mu"); auto mu = dict.getVal<realVector>("mu");
@ -114,12 +107,6 @@ protected:
return false; return false;
} }
if(nElem != et.size())
{
fatalErrorInFunction<<
"sizes of kn("<<nElem<<") and et("<<et.size()<<") do not match.\n";
return false;
}
if(nElem != mu.size()) if(nElem != mu.size())
{ {
@ -148,30 +135,16 @@ protected:
} }
realVector etha_t("etha_t", nElem);
ForAll(i , kn)
{
etha_t[i] = -2.0*log( et[i])*sqrt(kt[i]*2/7) /
sqrt(log(pow(et[i],2.0))+ pow(Pi,2.0));
}
Vector<linearProperties> prop("prop", nElem); Vector<linearProperties> prop("prop", nElem);
ForAll(i,kn) ForAll(i,kn)
{ {
prop[i] = {kn[i], kt[i], en[i], etha_t[i], mu[i] }; prop[i] = {kn[i], kt[i], en[i], mu[i] };
} }
linearProperties_.assign(prop); linearProperties_.assign(prop);
auto adm = dict.getVal<word>("additionalDissipationModel");
auto adm = dict.getVal<word>("additionalDissipationModel");
if(adm == "none") if(adm == "none")
{ {
addDissipationModel_ = 1; addDissipationModel_ = 1;
@ -282,27 +255,40 @@ public:
history.overlap_t_ += Vt*dt; history.overlap_t_ += Vt*dt;
real mi = 3*Pi/4*pow(Ri,3.0)*rho_[propId_i]; real mi = 3*Pi/4*pow(Ri,3)*rho_[propId_i];
real mj = 3*Pi/4*pow(Rj,3.0)*rho_[propId_j]; real mj = 3*Pi/4*pow(Rj,3)*rho_[propId_j];
real sqrt_meff = sqrt((mi*mj)/(mi+mj)); real sqrt_meff = sqrt((mi*mj)/(mi+mj));
real en = prop.en_;
if (addDissipationModel_==2) if (addDissipationModel_==2)
{ {
prop.en_ = sqrt(1+((pow(prop.en_,2)-1)*f_)); en = sqrt(1+((pow(prop.en_,2)-1)*f_));
} }
else if (addDissipationModel_==3) else if (addDissipationModel_==3)
{ {
auto pie =3.14;
prop.en_ = exp((pow(f_,1.5)*log(prop.en_)*sqrt( (1-((pow(log(prop.en_),2))/(pow(log(prop.en_),2)+pow(pie,2))))/(1-(pow(f_,3)*(pow(log(prop.en_),2))/(pow(log(prop.en_),2)+pow(pie,2)))) ) )); en = exp(
(
pow(f_,static_cast<real>(1.5))*
log(prop.en_)*
sqrt(
(1-((pow(log(prop.en_),static_cast<real>(2.0))
)/
(
pow(log(prop.en_),static_cast<real>(2.0))+
pow(Pi,static_cast<real>(2.0)))))/
(1-(pow(f_,3)*(pow(log(prop.en_),2))/
(pow(log(prop.en_),static_cast<real>(2.0))+
pow(Pi,static_cast<real>(2.0))))) ) ));
} }
real ethan_ = -2.0*log(prop.en_)*sqrt(prop.kn_)/
sqrt(pow(log(prop.en_),2.0)+ pow(Pi,2.0)); real ethan_ = -2.0*log(en)*sqrt(prop.kn_)/
sqrt(pow(log(en),2)+ pow(Pi,2));
FCn = ( -pow(f_,1.0)*prop.kn_ * ovrlp_n - sqrt_meff * pow(f_,0.5) * ethan_ * vrn)*Nij; FCn = ( -f_*prop.kn_ * ovrlp_n - sqrt_meff * pow(f_,half) * ethan_ * vrn)*Nij;
FCt = ( -pow(f_,1.0)*prop.kt_ * history.overlap_t_ - sqrt_meff * pow(f_,0.5) * prop.ethat_*Vt); FCt = ( -f_*prop.kt_ * history.overlap_t_);

View File

@ -139,10 +139,10 @@ protected:
ForAll(i , kn) ForAll(i , kn)
{ {
etha_n[i] = -2.0*log(en[i])*sqrt(kn[i])/ etha_n[i] = -2.0*log(en[i])*sqrt(kn[i])/
sqrt(pow(log(en[i]),2.0)+ pow(Pi,2.0)); sqrt(pow(log(en[i]),static_cast<real>(2.0))+ pow(Pi,static_cast<real>(2.0)));
etha_t[i] = -2.0*log( et[i]*sqrt(kt[i]) )/ etha_t[i] = -2.0*log( et[i]*sqrt(kt[i]) )/
sqrt(pow(log(et[i]),2.0)+ pow(Pi,2.0)); sqrt(pow(log(et[i]),static_cast<real>(2.0))+ pow(Pi,static_cast<real>(2.0)));
} }
Vector<linearProperties> prop("prop", nElem); Vector<linearProperties> prop("prop", nElem);
@ -243,8 +243,8 @@ public:
history.overlap_t_ += Vt*dt; history.overlap_t_ += Vt*dt;
real mi = 3*Pi/4*pow(Ri,3.0)*rho_[propId_i]; real mi = 3*Pi/4*pow(Ri,static_cast<real>(3.0))*rho_[propId_i];
real mj = 3*Pi/4*pow(Rj,3.0)*rho_[propId_j]; real mj = 3*Pi/4*pow(Rj,static_cast<real>(3.0))*rho_[propId_j];
real sqrt_meff = sqrt((mi*mj)/(mi+mj)); real sqrt_meff = sqrt((mi*mj)/(mi+mj));

View File

@ -131,7 +131,7 @@ protected:
// we take out sqrt(meff*K_hertz) here and then consider this term // we take out sqrt(meff*K_hertz) here and then consider this term
// when calculating damping part. // when calculating damping part.
etha_n[i] = -2.2664*log(en[i])/ etha_n[i] = -2.2664*log(en[i])/
sqrt(pow(log(en[i]),2.0)+ pow(Pi,2.0)); sqrt(pow(log(en[i]),static_cast<real>(2.0))+ pow(Pi,static_cast<real>(2.0)));
// no damping for tangential part // no damping for tangential part
@ -255,7 +255,7 @@ public:
// apply friction // apply friction
if(ft > ft_fric) if(ft > ft_fric)
{ {
if( length(history.overlap_t_) >0.0) if( length(history.overlap_t_) >zero)
{ {
if constexpr (limited) if constexpr (limited)
{ {

View File

@ -23,6 +23,7 @@ Licence:
#include "cGAbsoluteLinearCF.hpp" #include "cGAbsoluteLinearCF.hpp"
#include "cGRelativeLinearCF.hpp" #include "cGRelativeLinearCF.hpp"
#include "cGNonLinearCF.hpp"
#include "grainRolling.hpp" #include "grainRolling.hpp"
@ -38,6 +39,10 @@ using nonLimitedCGAbsoluteLinearGrainRolling = grainRolling<cGAbsoluteLinear<fal
using limitedCGRelativeLinearGrainRolling = grainRolling<cGRelativeLinear<true>>; using limitedCGRelativeLinearGrainRolling = grainRolling<cGRelativeLinear<true>>;
using nonLimitedCGRelativeLinearGrainRolling = grainRolling<cGRelativeLinear<false>>; using nonLimitedCGRelativeLinearGrainRolling = grainRolling<cGRelativeLinear<false>>;
using limitedCGNonLinearGrainRolling = grainRolling<cGNonLinear<true>>;
using nonLimitedCGNonLinearGrainRolling = grainRolling<cGNonLinear<false>>;
} }

View File

@ -96,10 +96,10 @@ public:
realx3 w_hat = wi-wj; realx3 w_hat = wi-wj;
real w_hat_mag = length(w_hat); real w_hat_mag = length(w_hat);
if( !equal(w_hat_mag,0.0) ) if( !equal(w_hat_mag,static_cast<real>(0.0)) )
w_hat /= w_hat_mag; w_hat /= w_hat_mag;
else else
w_hat = 0.0; w_hat = static_cast<real>(0.0);
auto Reff = (Ri*Rj)/(Ri+Rj); auto Reff = (Ri*Rj)/(Ri+Rj);

View File

@ -94,10 +94,10 @@ public:
realx3 w_hat = wi-wj; realx3 w_hat = wi-wj;
real w_hat_mag = length(w_hat); real w_hat_mag = length(w_hat);
if( !equal(w_hat_mag,0.0) ) if( !equal(w_hat_mag,static_cast<real>(0.0)) )
w_hat /= w_hat_mag; w_hat /= w_hat_mag;
else else
w_hat = 0.0; w_hat = static_cast<real>(0.0);
auto Reff = (Ri*Rj)/(Ri+Rj); auto Reff = (Ri*Rj)/(Ri+Rj);

View File

@ -80,13 +80,17 @@ public:
TypeInfoNV("sortedContactList"); TypeInfoNV("sortedContactList");
sortedContactList(uint32 initialSize =1)
explicit sortedContactList(uint32 initialSize =1)
: :
SortedPairs(initialSize), sortedContactList("sortedContactList", initialSize)
values_("values", SortedPairs::capacity()), {}
sortedPairs0_("sortedPairs0", SortedPairs::capacity()),
values0_("values0", SortedPairs::capacity()) sortedContactList(const word& name, uint32 initialSize =1)
:
SortedPairs(name, initialSize),
values_(groupNames(name, "values"), SortedPairs::capacity()),
sortedPairs0_(groupNames(name, "sortedPairs0"), SortedPairs::capacity()),
values0_(groupNames(name, "values0"), SortedPairs::capacity())
{} {}
bool beforeBroadSearch() bool beforeBroadSearch()

View File

@ -110,11 +110,11 @@ public:
// constructors // constructors
explicit sortedPairs(uint32 initialSize =1) explicit sortedPairs(const word& name, uint32 initialSize =1)
: :
UnsortedPairs(initialSize), UnsortedPairs(initialSize),
flags_("flags_",UnsortedPairs::capacity()+1), flags_( groupNames(name, "flags_"), UnsortedPairs::capacity()+1),
sortedPairs_("sortedPairs_",UnsortedPairs::capacity()) sortedPairs_(groupNames(name, "sortedPairs_"), UnsortedPairs::capacity())
{} {}
@ -193,7 +193,7 @@ public:
if( capacity+1 > flags_.size() ) if( capacity+1 > flags_.size() )
{ {
reallocNoInit(flags_, capacity+1); reallocInit(flags_, capacity+1);
} }
// fill the flags // fill the flags
@ -219,7 +219,7 @@ public:
{ {
// get more space to prevent reallocations in next iterations // get more space to prevent reallocations in next iterations
uint32 len = size_*1.1+1; uint32 len = size_*1.1+1;
reallocNoInit(sortedPairs_, len); reallocInit(sortedPairs_, len);
} }
Kokkos::parallel_for( Kokkos::parallel_for(
@ -232,6 +232,7 @@ public:
sort(sortedPairs_, 0, size_ ); sort(sortedPairs_, 0, size_ );
} }
INLINE_FUNCTION_HD INLINE_FUNCTION_HD

View File

@ -82,11 +82,16 @@ public:
TypeInfoNV("unsortedContactList"); TypeInfoNV("unsortedContactList");
explicit unsortedContactList(uint32 capacity=1) explicit unsortedContactList(uint32 capacity=1)
:
unsortedContactList("unsortedContactList", capacity)
{}
unsortedContactList(const word& name, uint32 capacity=1)
: :
UnsortedPairs(capacity), UnsortedPairs(capacity),
values_("values", UnsortedPairs::capacity()), values_(groupNames(name, "values"), UnsortedPairs::capacity()),
container0_(capacity), container0_(capacity),
values0_("values0",container0_.capacity()) values0_(groupNames(name, "values0"),container0_.capacity())
{} {}

View File

@ -194,7 +194,9 @@ public:
{ {
uint newCap = container_.capacity()+len; uint newCap = container_.capacity()+len;
this->clear(); this->clear();
//output<<"----------------before "<<capacity()<< " " << size()<<endl;
container_.rehash(newCap); container_.rehash(newCap);
//output<<"----------------after "<<capacity()<< " " << size()<<endl;
} }
INLINE_FUNCTION_H INLINE_FUNCTION_H

View File

@ -60,9 +60,9 @@ private:
bool force = false) override bool force = false) override
{ {
const auto& position = Particles().pointPosition().deviceViewAll(); auto position = Particles().pointPosition().deviceViewAll();
const auto& flags = Particles().dynPointStruct().activePointsMaskDevice(); auto flags = Particles().dynPointStruct().activePointsMaskDevice();
const auto& diam = Particles().boundingSphere().deviceViewAll(); auto diam = Particles().boundingSphere().deviceViewAll();
return ppwContactSearch_().broadSearch( return ppwContactSearch_().broadSearch(
ti.iter(), ti.iter(),
@ -134,7 +134,7 @@ public:
ppwContactSearch_ = ppwContactSearch_ =
makeUnique<SearchMethodType> makeUnique<SearchMethodType>
( (
dict(), csDict,
this->extendedDomainBox(), this->extendedDomainBox(),
minD, minD,
maxD, maxD,
@ -154,10 +154,8 @@ public:
ContactSearch, ContactSearch,
dictionary); dictionary);
bool enterBroadSearchBoundary(const timeInfo& ti, bool force=false)const override bool enterBroadSearchBoundary(const timeInfo& ti, bool force=false)const override;
{
return enterBroadSearch(ti, force) || csBoundaries_.boundariesUpdated();
}
real sizeRatio()const override real sizeRatio()const override
{ {
@ -171,7 +169,12 @@ public:
}; };
template <class searchMethod>
inline bool ContactSearch<searchMethod>::enterBroadSearchBoundary(const timeInfo &ti, bool force) const
{
return performSearch(ti.iter(),force) || csBoundaries_.boundariesUpdated();
} }
}
#endif //__ContactSearch_hpp__ #endif //__ContactSearch_hpp__

View File

@ -73,9 +73,7 @@ public:
} }
bool hearChanges( bool hearChanges(
real t, const timeInfo& ti,
real dt,
uint32 iter,
const message &msg, const message &msg,
const anyList &varList) override const anyList &varList) override
{ {
@ -83,8 +81,10 @@ public:
if (msg.equivalentTo(message::BNDR_RESET)) if (msg.equivalentTo(message::BNDR_RESET))
{ {
// do nothing // do nothing
return true;
} }
return true; fatalErrorInFunction;
return false;
} }
virtual bool broadSearch( virtual bool broadSearch(
@ -98,6 +98,11 @@ public:
return true; return true;
} }
bool isActive()const override
{
return false;
}
static uniquePtr<boundaryContactSearch> create( static uniquePtr<boundaryContactSearch> create(
const dictionary &dict, const dictionary &dict,
const boundaryBase &boundary, const boundaryBase &boundary,

View File

@ -1,3 +1,23 @@
/*------------------------------- phasicFlow ---------------------------------
O C enter of
O O E ngineering and
O O M ultiscale modeling of
OOOOOOO F luid flow
------------------------------------------------------------------------------
Copyright (C): www.cemf.ir
email: hamid.r.norouzi AT gmail.com
------------------------------------------------------------------------------
Licence:
This file is part of phasicFlow code. It is a free software for simulating
granular and multiphase flows. You can redistribute it and/or modify it under
the terms of GNU General Public License v3 or any other later versions.
phasicFlow is distributed to help others in their research in the field of
granular and multiphase flows, but WITHOUT ANY WARRANTY; without even the
implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
-----------------------------------------------------------------------------*/
#include "boundaryContactSearchList.hpp" #include "boundaryContactSearchList.hpp"
#include "boundaryList.hpp" #include "boundaryList.hpp"
@ -5,7 +25,7 @@ void pFlow::boundaryContactSearchList::setList(
const dictionary &dict, const dictionary &dict,
const contactSearch &cSearch) const contactSearch &cSearch)
{ {
for(auto i=0; i<boundaries_.size(); i++) ForAllBoundariesPtr(i, this)
{ {
this->set this->set
( (
@ -20,7 +40,7 @@ pFlow::boundaryContactSearchList::boundaryContactSearchList(
const boundaryList& bndrs, const boundaryList& bndrs,
const contactSearch &cSearch) const contactSearch &cSearch)
: :
ListPtr(bndrs.size()), boundaryListPtr(),
boundaries_(bndrs) boundaries_(bndrs)
{ {
setList(dict, cSearch); setList(dict, cSearch);

View File

@ -1,4 +1,24 @@
#include "ListPtr.hpp" /*------------------------------- 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 "boundaryListPtr.hpp"
#include "boundaryContactSearch.hpp" #include "boundaryContactSearch.hpp"
namespace pFlow namespace pFlow
@ -9,7 +29,7 @@ class contactSearch;
class boundaryContactSearchList class boundaryContactSearchList
: :
public ListPtr<boundaryContactSearch> public boundaryListPtr<boundaryContactSearch>
{ {
private: private:

View File

@ -74,6 +74,11 @@ public:
csPairContainerType &ppPairs, csPairContainerType &ppPairs,
csPairContainerType &pwPairs, csPairContainerType &pwPairs,
bool force = false) override; bool force = false) override;
bool isActive()const override
{
return true;
}
}; };
} }

View File

@ -11,7 +11,7 @@ pFlow::wallBoundaryContactSearch::wallBoundaryContactSearch
const ViewType1D<realx3, memory_space> &normals const ViewType1D<realx3, memory_space> &normals
) )
: :
cellExtent_( max(cellExtent, 0.5 ) ), cellExtent_( max(cellExtent, half ) ),
numElements_(numElements), numElements_(numElements),
numPoints_(numPoints), numPoints_(numPoints),
vertices_(vertices), vertices_(vertices),

View File

@ -31,11 +31,11 @@ pFlow::contactSearch::contactSearch(
Timers& timers) Timers& timers)
: :
extendedDomainBox_(extDomain), extendedDomainBox_(extDomain),
updateInterval_(dict.getValMax<uint32>("updateInterval", 1u)),
particles_(prtcl), particles_(prtcl),
geometry_(geom), geometry_(geom),
bTimer_("Boundary particles contact search", &timers), bTimer_("Boundary particles contact search", &timers),
ppTimer_("Internal particles contact search", &timers), ppTimer_("Internal particles contact search", &timers)
dict_(dict)
{ {
} }
@ -69,6 +69,7 @@ bool pFlow::contactSearch::broadSearch
} }
ppTimer_.end(); ppTimer_.end();
performedSearch_ = true; performedSearch_ = true;
lastUpdated_ = ti.currentIter();
} }
else else
{ {
@ -92,6 +93,7 @@ bool pFlow::contactSearch::boundaryBroadSearch
bTimer_.start(); bTimer_.start();
for(uint32 i=0u; i<6u; i++) for(uint32 i=0u; i<6u; i++)
{ {
//output<<" boundarySearch "<< i <<" for iter "<< ti.iter()<<endl;
if(!BoundaryBroadSearch( if(!BoundaryBroadSearch(
i, i,
ti, ti,

View File

@ -67,8 +67,6 @@ private:
Timer ppTimer_; Timer ppTimer_;
dictionary dict_;
virtual virtual
bool BroadSearch( bool BroadSearch(
const timeInfo& ti, const timeInfo& ti,
@ -150,12 +148,6 @@ public:
return updateInterval_; return updateInterval_;
} }
inline
const dictionary& dict()const
{
return dict_;
}
inline inline
const box& extendedDomainBox()const const box& extendedDomainBox()const
{ {

View File

@ -44,9 +44,9 @@ pFlow::NBS::NBS
position, position,
flags, flags,
diam), diam),
sizeRatio_(max(dict.getVal<real>("sizeRatio"), 1.0)), sizeRatio_(max(dict.getVal<real>("sizeRatio"), one)),
cellExtent_(max(dict.getVal<real>("cellExtent"), 0.5)), cellExtent_(max(dict.getVal<real>("cellExtent"), half)),
adjustableBox_(dict.getVal<Logical>("adjustableBox")), adjustableBox_(false),//adjustableBox_(dict.getVal<Logical>("adjustableBox")),
NBSLevel0_ NBSLevel0_
( (
this->domainBox_, this->domainBox_,

View File

@ -31,7 +31,7 @@ pFlow::cellsWallLevel0::cellsWallLevel0
const ViewType1D<realx3, memory_space>& normals const ViewType1D<realx3, memory_space>& normals
) )
: :
cellExtent_( max(cellExtent, 0.5 ) ), cellExtent_( max(cellExtent, half ) ),
numElements_(numElements), numElements_(numElements),
numPoints_(numPoints), numPoints_(numPoints),
vertices_(vertices), vertices_(vertices),

View File

@ -23,7 +23,7 @@ Licence:
#include "streams.hpp" #include "streams.hpp"
pFlow::uint32 pFlow::mapperNBS::checkInterval_ = 1000; pFlow::uint32 pFlow::mapperNBS::checkInterval_ = 1000;
pFlow::real pFlow::mapperNBS::enlargementFactor_ = 1.1; pFlow::real pFlow::mapperNBS::enlargementFactor_ = 1.5;
bool pFlow::mapperNBS::setSearchBox bool pFlow::mapperNBS::setSearchBox
( (

View File

@ -21,7 +21,7 @@ Licence:
-----------------------------------------------------------------------------*/ -----------------------------------------------------------------------------*/
#include "mapperNBSKernels.hpp" #include "mapperNBSKernels.hpp"
#include "streams.hpp"
void pFlow::mapperNBSKernels::findPointExtends void pFlow::mapperNBSKernels::findPointExtends
( (
@ -153,7 +153,10 @@ bool pFlow::mapperNBSKernels::buildLists
const pFlagTypeDevice &flags const pFlagTypeDevice &flags
) )
{ {
auto aRange = flags.activeRange(); auto aRange = flags.activeRange();
auto pp = points;
if(flags.isAllActive() ) if(flags.isAllActive() )
{ {
Kokkos::parallel_for Kokkos::parallel_for
@ -162,7 +165,7 @@ bool pFlow::mapperNBSKernels::buildLists
deviceRPolicyStatic(aRange.start(), aRange.end()), deviceRPolicyStatic(aRange.start(), aRange.end()),
LAMBDA_HD(uint32 i) LAMBDA_HD(uint32 i)
{ {
auto ind = searchCell.pointIndex(points[i]); auto ind = searchCell.pointIndex(pp[i]);
uint32 old = Kokkos::atomic_exchange(&head(ind.x(), ind.y(), ind.z()), i); uint32 old = Kokkos::atomic_exchange(&head(ind.x(), ind.y(), ind.z()), i);
next[i] = old; next[i] = old;
} }

View File

@ -152,22 +152,24 @@ public:
return false; return false;
} }
bool hearChanges bool hearChanges
( (
real t, const timeInfo& ti,
real dt,
uint32 iter,
const message& msg, const message& msg,
const anyList& varList const anyList& varList
) override ) override
{ {
if(msg.equivalentTo(message::BNDR_RESET))
{
return true;
}
fatalErrorInFunction<<"Event "<< msg.eventNames() << " is not handled !"<<endl;
return false;
}
pOutput<<"Function (hearChanges in boundarySphereInteractions)is not implmented Message "<< bool isActive()const override
msg <<endl<<" name "<< this->boundaryName() <<" type "<< this->type()<<endl;; {
//notImplementedFunction; return false;
return true;
} }
static static

View File

@ -7,11 +7,11 @@ pFlow::boundaryGrainInteractionList<CFModel, gMModel>::boundaryGrainInteractionL
const gMModel &geomMotion const gMModel &geomMotion
) )
: :
ListPtr<boundaryGrainInteraction<CFModel,gMModel>>(6), boundaryListPtr<boundaryGrainInteraction<CFModel,gMModel>>(),
boundaries_(grnPrtcls.pStruct().boundaries()) boundaries_(grnPrtcls.pStruct().boundaries())
{ {
//gSettings::sleepMiliSeconds(1000*pFlowProcessors().localRank()); //gSettings::sleepMiliSeconds(1000*pFlowProcessors().localRank());
for(uint32 i=0; i<6; i++) ForAllBoundariesPtr(i, this)
{ {
this->set( this->set(
i, i,

View File

@ -3,7 +3,7 @@
#include "boundaryList.hpp" #include "boundaryList.hpp"
#include "ListPtr.hpp" #include "boundaryListPtr.hpp"
#include "boundaryGrainInteraction.hpp" #include "boundaryGrainInteraction.hpp"
@ -14,7 +14,7 @@ namespace pFlow
template<typename contactForceModel,typename geometryMotionModel> template<typename contactForceModel,typename geometryMotionModel>
class boundaryGrainInteractionList class boundaryGrainInteractionList
: :
public ListPtr<boundaryGrainInteraction<contactForceModel,geometryMotionModel>> public boundaryListPtr<boundaryGrainInteraction<contactForceModel,geometryMotionModel>>
{ {
private: private:

View File

@ -78,14 +78,16 @@ public:
~periodicBoundaryGrainInteraction()override = default; ~periodicBoundaryGrainInteraction()override = default;
bool grainGrainInteraction( bool grainGrainInteraction(
real dt, real dt,
const ContactForceModel& cfModel, const ContactForceModel& cfModel,
uint32 step)override; uint32 step)override;
bool isActive()const override
{
return true;
}
}; };
} }

View File

@ -41,9 +41,9 @@ bool pFlow::grainInteraction<cFM,gMM, cLT>::createGrainInteraction()
geometryMotion_, geometryMotion_,
timers()); timers());
ppContactList_ = makeUnique<ContactListType>(nPrtcl+1); ppContactList_ = makeUnique<ContactListType>("Grain-Grain",nPrtcl+1);
pwContactList_ = makeUnique<ContactListType>(nPrtcl/5+1); pwContactList_ = makeUnique<ContactListType>("Grain-wall",nPrtcl/5+1);
return true; return true;
} }
@ -344,16 +344,12 @@ bool pFlow::grainInteraction<cFM,gMM, cLT>::afterIteration()
template<typename cFM,typename gMM,template <class, class, class> class cLT> template<typename cFM,typename gMM,template <class, class, class> class cLT>
bool pFlow::grainInteraction<cFM,gMM, cLT>::hearChanges bool pFlow::grainInteraction<cFM,gMM, cLT>::hearChanges
( (
real t, const timeInfo& ti,
real dt,
uint32 iter,
const message& msg, const message& msg,
const anyList& varList const anyList& varList
) )
{ {
if(msg.equivalentTo(message::ITEM_REARRANGE)) fatalErrorInFunction<<"Event "<< msg.eventNames()<<
{ " is not handled in grainInteraction"<<endl;
notImplementedFunction; return false;
}
return true;
} }

View File

@ -152,9 +152,7 @@ public:
/// Check for changes in the point structures. (overriden from observer) /// Check for changes in the point structures. (overriden from observer)
bool hearChanges( bool hearChanges(
real t, const timeInfo& ti,
real dt,
uint32 iter,
const message& msg, const message& msg,
const anyList& varList)override; const anyList& varList)override;

View File

@ -41,3 +41,9 @@ Licence:
createBoundaryGrainInteraction(ForceModel, GeomModel) createBoundaryGrainInteraction(ForceModel, GeomModel)
createInteraction(pFlow::cfModels::limitedCGNonLinearGrainRolling, pFlow::rotationAxisMotionGeometry);
createInteraction(pFlow::cfModels::nonLimitedCGNonLinearGrainRolling, pFlow::rotationAxisMotionGeometry);
createInteraction(pFlow::cfModels::limitedCGNonLinearGrainRolling, pFlow::stationaryGeometry);
createInteraction(pFlow::cfModels::nonLimitedCGNonLinearGrainRolling, pFlow::stationaryGeometry);

View File

@ -51,7 +51,7 @@ private:
const geometry& geometry_; const geometry& geometry_;
static inline static inline
const message msg_ = message::ITEM_REARRANGE; const message msg_ = message::ITEMS_REARRANGE;
public: public:

View File

@ -148,7 +148,7 @@ public:
const ContactForceModel& cfModel, const ContactForceModel& cfModel,
uint32 step) uint32 step)
{ {
// for default boundary, no thing to be done // for default boundary, nothing to be done
return false; return false;
} }
@ -156,20 +156,26 @@ public:
bool hearChanges bool hearChanges
( (
real t, const timeInfo& ti,
real dt,
uint32 iter,
const message& msg, const message& msg,
const anyList& varList const anyList& varList
) override ) override
{ {
if(msg.equivalentTo(message::BNDR_RESET))
{
// do nothing
return true;
}
pOutput<<"Function (hearChanges in boundarySphereInteractions)is not implmented Message "<< pOutput<<"Function (hearChanges in boundarySphereInteractions)is not implmented Message "<<
msg <<endl<<" name "<< this->boundaryName() <<" type "<< this->type()<<endl;; msg <<endl<<" name "<< this->boundaryName() <<" type "<< this->type()<<endl;
//notImplementedFunction;
return true; return true;
} }
bool isActive()const override
{
return false;
}
static static
uniquePtr<BoundarySphereInteractionType> create( uniquePtr<BoundarySphereInteractionType> create(
const boundaryBase& boundary, const boundaryBase& boundary,

View File

@ -7,11 +7,11 @@ pFlow::boundarySphereInteractionList<CFModel, gMModel>::boundarySphereInteractio
const gMModel &geomMotion const gMModel &geomMotion
) )
: :
ListPtr<boundarySphereInteraction<CFModel,gMModel>>(6), boundaryListPtr<boundarySphereInteraction<CFModel,gMModel>>(),
boundaries_(sphPrtcls.pStruct().boundaries()) boundaries_(sphPrtcls.pStruct().boundaries())
{ {
//gSettings::sleepMiliSeconds(1000*pFlowProcessors().localRank()); output<<boundaries_.size()<<endl;
for(uint32 i=0; i<6; i++) ForAllBoundariesPtr(i, this)
{ {
this->set( this->set(
i, i,

View File

@ -3,7 +3,7 @@
#include "boundaryList.hpp" #include "boundaryList.hpp"
#include "ListPtr.hpp" #include "boundaryListPtr.hpp"
#include "boundarySphereInteraction.hpp" #include "boundarySphereInteraction.hpp"
@ -14,7 +14,7 @@ namespace pFlow
template<typename contactForceModel,typename geometryMotionModel> template<typename contactForceModel,typename geometryMotionModel>
class boundarySphereInteractionList class boundarySphereInteractionList
: :
public ListPtr<boundarySphereInteraction<contactForceModel,geometryMotionModel>> public boundaryListPtr<boundarySphereInteraction<contactForceModel,geometryMotionModel>>
{ {
private: private:

View File

@ -78,14 +78,16 @@ public:
~periodicBoundarySphereInteraction()override = default; ~periodicBoundarySphereInteraction()override = default;
bool sphereSphereInteraction( bool sphereSphereInteraction(
real dt, real dt,
const ContactForceModel& cfModel, const ContactForceModel& cfModel,
uint32 step)override; uint32 step)override;
bool isActive()const override
{
return true;
}
}; };
} }

View File

@ -41,9 +41,9 @@ bool pFlow::sphereInteraction<cFM,gMM, cLT>::createSphereInteraction()
geometryMotion_, geometryMotion_,
timers()); timers());
ppContactList_ = makeUnique<ContactListType>(nPrtcl+1); ppContactList_ = makeUnique<ContactListType>("sphere-sphere",nPrtcl+1);
pwContactList_ = makeUnique<ContactListType>(nPrtcl/5+1); pwContactList_ = makeUnique<ContactListType>("sphere-wall",nPrtcl/5+1);
return true; return true;
} }
@ -193,7 +193,9 @@ bool pFlow::sphereInteraction<cFM,gMM, cLT>::iterate()
{ {
contactListMangementBoundaryTimer_.start(); contactListMangementBoundaryTimer_.start();
ComputationTimer().start(); ComputationTimer().start();
for(uint32 i=0; i<6u; i++)
ForAllActiveBoundaries(i, boundaryInteraction_)
//for(size_t i=0; i<6; i++)
{ {
if(activeBoundaries_[i]) if(activeBoundaries_[i])
{ {
@ -202,6 +204,7 @@ bool pFlow::sphereInteraction<cFM,gMM, cLT>::iterate()
BI.pwPairs().beforeBroadSearch(); BI.pwPairs().beforeBroadSearch();
} }
} }
ComputationTimer().end(); ComputationTimer().end();
contactListMangementBoundaryTimer_.pause(); contactListMangementBoundaryTimer_.pause();
} }
@ -219,7 +222,8 @@ bool pFlow::sphereInteraction<cFM,gMM, cLT>::iterate()
fatalExit; fatalExit;
} }
for(uint32 i=0; i<6u; i++) ForAllActiveBoundaries(i, boundaryInteraction_)
//for(size_t i=0; i<6; i++)
{ {
if(activeBoundaries_[i]) if(activeBoundaries_[i])
{ {
@ -253,7 +257,9 @@ bool pFlow::sphereInteraction<cFM,gMM, cLT>::iterate()
{ {
contactListMangementBoundaryTimer_.resume(); contactListMangementBoundaryTimer_.resume();
ComputationTimer().start(); ComputationTimer().start();
for(uint32 i=0; i<6u; i++)
ForAllActiveBoundaries(i, boundaryInteraction_)
//for(size_t i=0; i<6; i++)
{ {
if(activeBoundaries_[i]) if(activeBoundaries_[i])
{ {
@ -262,6 +268,7 @@ bool pFlow::sphereInteraction<cFM,gMM, cLT>::iterate()
BI.pwPairs().afterBroadSearch(); BI.pwPairs().afterBroadSearch();
} }
} }
ComputationTimer().end(); ComputationTimer().end();
contactListMangementBoundaryTimer_.end(); contactListMangementBoundaryTimer_.end();
} }
@ -274,7 +281,8 @@ bool pFlow::sphereInteraction<cFM,gMM, cLT>::iterate()
ComputationTimer().start(); ComputationTimer().start();
while(requireStep.anyElement(true) && step <= 10) while(requireStep.anyElement(true) && step <= 10)
{ {
for(uint32 i=0; i<6u; i++) ForAllBoundaries(i, boundaryInteraction_)
//for(size_t i=0; i<6; i++)
{ {
if(requireStep[i] ) if(requireStep[i] )
{ {
@ -313,7 +321,8 @@ bool pFlow::sphereInteraction<cFM,gMM, cLT>::iterate()
const auto& cfModel = this->forceModel_(); const auto& cfModel = this->forceModel_();
while( requireStep.anyElement(true) && step < 20 ) while( requireStep.anyElement(true) && step < 20 )
{ {
for(uint32 i=0; i<6u; i++) ForAllBoundaries(i, boundaryInteraction_)
//for(size_t i=0; i<6; i++)
{ {
if(requireStep[i]) if(requireStep[i])
{ {
@ -342,16 +351,18 @@ bool pFlow::sphereInteraction<cFM,gMM, cLT>::afterIteration()
template<typename cFM,typename gMM,template <class, class, class> class cLT> template<typename cFM,typename gMM,template <class, class, class> class cLT>
bool pFlow::sphereInteraction<cFM,gMM, cLT>::hearChanges bool pFlow::sphereInteraction<cFM,gMM, cLT>::hearChanges
( (
real t, const timeInfo& ti,
real dt,
uint32 iter,
const message& msg, const message& msg,
const anyList& varList const anyList& varList
) )
{ {
if(msg.equivalentTo(message::ITEM_REARRANGE)) if(msg.equivalentTo(message::ITEMS_REARRANGE))
{ {
notImplementedFunction; notImplementedFunction;
return false;
} }
return true;
fatalErrorInFunction<<"Event "<< msg.eventNames()<<
" is not handled in sphereInteraction"<<endl;
return false;
} }

View File

@ -152,9 +152,7 @@ public:
/// Check for changes in the point structures. (overriden from observer) /// Check for changes in the point structures. (overriden from observer)
bool hearChanges( bool hearChanges(
real t, const timeInfo& ti,
real dt,
uint32 iter,
const message& msg, const message& msg,
const anyList& varList)override; const anyList& varList)override;

View File

@ -58,5 +58,5 @@ createInteraction(pFlow::cfModels::limitedNonLinearModNormalRolling, pFlow::conv
createInteraction(pFlow::cfModels::nonLimitedNonLinearModNormalRolling,pFlow::conveyorBeltMotionGeometry); createInteraction(pFlow::cfModels::nonLimitedNonLinearModNormalRolling,pFlow::conveyorBeltMotionGeometry);
// multiRotationAxisMotionGeometry // multiRotationAxisMotionGeometry
//createInteraction(pFlow::cfModels::limitedNonLinearModNormalRolling, pFlow::multiRotationAxisMotionGeometry); createInteraction(pFlow::cfModels::limitedNonLinearModNormalRolling, pFlow::multiRotationAxisMotionGeometry);
//createInteraction(pFlow::cfModels::nonLimitedNonLinearModNormalRolling,pFlow::multiRotationAxisMotionGeometry); createInteraction(pFlow::cfModels::nonLimitedNonLinearModNormalRolling,pFlow::multiRotationAxisMotionGeometry);

View File

@ -14,8 +14,8 @@ entities/stationary/stationary.cpp
conveyorBeltMotion/conveyorBeltMotion.cpp conveyorBeltMotion/conveyorBeltMotion.cpp
entities/conveyorBelt/conveyorBelt.cpp entities/conveyorBelt/conveyorBelt.cpp
#entities/multiRotatingAxis/multiRotatingAxis.cpp entities/multiRotatingAxis/multiRotatingAxis.cpp
#multiRotatingAxisMotion/multiRotatingAxisMotion.cpp multiRotatingAxisMotion/multiRotatingAxisMotion.cpp
) )

View File

@ -32,7 +32,6 @@ bool pFlow::MotionModel<Model, Component>::impl_nameToIndex(const word& name, ui
indx = static_cast<uint32>(i); indx = static_cast<uint32>(i);
return true; return true;
} }
} }
template<typename Model, typename Component> template<typename Model, typename Component>

View File

@ -22,31 +22,32 @@ Licence:
#include "multiRotatingAxisMotion.hpp" #include "multiRotatingAxisMotion.hpp"
#include "dictionary.hpp" #include "dictionary.hpp"
/// Construct from dictionary
FUNCTION_H FUNCTION_H
pFlow::multiRotatingAxis::multiRotatingAxis pFlow::multiRotatingAxis::multiRotatingAxis(const dictionary& dict)
( :
multiRotatingAxisMotion* axisMotion rotatingAxis(dict)
) {}
{
//axisList_ = axisMotion->getAxisListPtr();
}
FUNCTION_H FUNCTION_H
pFlow::multiRotatingAxis::multiRotatingAxis pFlow::multiRotatingAxis::multiRotatingAxis
( (
multiRotatingAxisMotion* axisMotion, multiRotatingAxis* axisListPtr,
const wordList& componentsNames,
const dictionary& dict const dictionary& dict
) )
:
rotatingAxis(dict),
axisList_(axisListPtr)
{ {
if(!read(axisMotion, dict)) if(!read(dict, componentsNames))
{ {
fatalErrorInFunction<< fatalErrorInFunction<<
" error in reading rotatingAxis from dictionary "<< dict.globalName()<<endl; " error in reading rotatingAxis from dictionary "<< dict.globalName()<<endl;
fatalExit; fatalExit;
} }
//axisList_ = axisMotion->getAxisListPtr();
} }
@ -54,22 +55,29 @@ pFlow::multiRotatingAxis::multiRotatingAxis
FUNCTION_H FUNCTION_H
bool pFlow::multiRotatingAxis::read bool pFlow::multiRotatingAxis::read
( (
multiRotatingAxisMotion* axisMotion, const dictionary& dict,
const dictionary& dict const wordList& componentNames
) )
{ {
if(!rotatingAxis::read(dict))return false;
word rotAxis = dict.getValOrSet<word>("rotationAxis", "none"); word rotAxis = dict.getValOrSet<word>("rotationAxis", "none");
if(rotAxis == "none") if(rotAxis == "none")
{ {
parentAxisIndex_ = -1; parentAxisIndex_ = static_cast<uint32>(-1);
} }
else else
{ {
parentAxisIndex_ = axisMotion-> nameToIndex(rotAxis); if( auto i = componentNames.findi(rotAxis); i != -1 )
{
parentAxisIndex_ = i;
}
else
{
fatalErrorInFunction<<"crotationAxis "<< rotAxis<<" in dictionary "<<
dict.globalName()<<" is not found in list of axis names "<< componentNames<<endl;
return false;
}
} }
return true; return true;
@ -78,8 +86,8 @@ bool pFlow::multiRotatingAxis::read
FUNCTION_H FUNCTION_H
bool pFlow::multiRotatingAxis::write bool pFlow::multiRotatingAxis::write
( (
const multiRotatingAxisMotion* axisMotion, dictionary& dict,
dictionary& dict const wordList& componentNames
) const ) const
{ {
if( !rotatingAxis::write(dict) ) return false; if( !rotatingAxis::write(dict) ) return false;
@ -90,10 +98,8 @@ bool pFlow::multiRotatingAxis::write
} }
else else
{ {
auto rotAxis = axisMotion->indexToName(parentAxisIndex_); dict.add("rotationAxis", componentNames[parentAxisIndex_]);
dict.add("rotationAxis", rotAxis);
} }
return true; return true;
} }

View File

@ -24,7 +24,7 @@ Licence:
#include "rotatingAxis.hpp" #include "rotatingAxis.hpp"
#include "KokkosTypes.hpp" #include "KokkosTypes.hpp"
#include "List.hpp"
namespace pFlow namespace pFlow
{ {
@ -79,26 +79,31 @@ class multiRotatingAxis
protected: protected:
/// This is device pointer to all axes /// This is device pointer to all axes
multiRotatingAxis* axisList_; multiRotatingAxis* axisList_ = nullptr;
/// Index of parent axis /// Index of parent axis
int32 parentAxisIndex_ = -1; uint32 parentAxisIndex_ = static_cast<uint32>(-1);
public: public:
TypeInfoNV("multiRotatingAxis");
// - Constructors // - Constructors
/// Empty Constructor /// Empty Constructor
INLINE_FUNCTION_HD FUNCTION_HD
multiRotatingAxis(){} multiRotatingAxis() = default;
/// Empty with list of axes /// Construct from dictionary
FUNCTION_H FUNCTION_H
multiRotatingAxis(multiRotatingAxisMotion* axisMotion); explicit multiRotatingAxis(const dictionary& dict);
/// Construct from dictionary and list of axes /// Construct from dictionary and list of axes
FUNCTION_H FUNCTION_H
multiRotatingAxis(multiRotatingAxisMotion* axisMotion, const dictionary& dict); multiRotatingAxis(
multiRotatingAxis* axisListPtr,
const wordList& componentsNames,
const dictionary& dict);
/// Copy constructor /// Copy constructor
FUNCTION_HD FUNCTION_HD
@ -123,11 +128,11 @@ public:
while(parIndex != -1) while(parIndex != -1)
{ {
auto& ax = axisList_[parIndex]; auto& ax = axisList_[parIndex];
parentVel += ax.linTangentialVelocityPoint(p); parentVel += ax.linVelocityPoint(p);
parIndex = ax.parentAxisIndex(); parIndex = ax.parentAxisIndex();
} }
return parentVel + rotatingAxis::linTangentialVelocityPoint(p); return parentVel + rotatingAxis::linVelocityPoint(p);
} }
/// Translate point p for dt seconds based on the axis information /// Translate point p for dt seconds based on the axis information
@ -143,7 +148,7 @@ public:
} }
auto parIndex = parentAxisIndex_; auto parIndex = parentAxisIndex_;
while(parIndex != -1) while(parIndex != static_cast<uint32>(-1))
{ {
auto& ax = axisList_[parIndex]; auto& ax = axisList_[parIndex];
newP = pFlow::rotate(newP, ax, dt); newP = pFlow::rotate(newP, ax, dt);
@ -157,12 +162,12 @@ public:
INLINE_FUNCTION_HD INLINE_FUNCTION_HD
bool hasParent()const bool hasParent()const
{ {
return parentAxisIndex_ > -1; return parentAxisIndex_ != static_cast<uint32>(-1);
} }
/// Return the index of parent axis /// Return the index of parent axis
INLINE_FUNCTION_HD INLINE_FUNCTION_HD
int32 parentAxisIndex()const uint32 parentAxisIndex()const
{ {
return parentAxisIndex_; return parentAxisIndex_;
} }
@ -182,6 +187,7 @@ public:
* It is assumed that the axis with deepest level (with more parrents) is * It is assumed that the axis with deepest level (with more parrents) is
* moved first and then the axis with lower levels. * moved first and then the axis with lower levels.
*/ */
INLINE_FUNCTION_HD
void move(real dt) void move(real dt)
{ {
@ -201,11 +207,12 @@ public:
/// Read from dictionary /// Read from dictionary
FUNCTION_H FUNCTION_H
bool read(multiRotatingAxisMotion* axisMotion, const dictionary& dict); bool read(const dictionary& dict, const wordList& componentNames);
/// Write to dictionary /// Write to dictionary
FUNCTION_H FUNCTION_H
bool write(const multiRotatingAxisMotion* axisMotion, dictionary& dict) const; bool write(dictionary& dict, const wordList& componentNames) const;
}; };

View File

@ -47,7 +47,7 @@ pFlow::rotatingAxis::rotatingAxis
timeInterval(), timeInterval(),
line(p1, p2), line(p1, p2),
omega_(omega), omega_(omega),
rotating_(!equal(omega,0.0)) rotating_(!equal(omega,static_cast<real>(0.0)))
{ {
} }
@ -58,7 +58,7 @@ pFlow::real pFlow::rotatingAxis::setOmega(real omega)
{ {
auto tmp = omega_; auto tmp = omega_;
omega_ = omega; omega_ = omega;
rotating_ = !equal(omega, 0.0); rotating_ = !equal(omega, static_cast<real>(0.0));
return tmp; return tmp;
} }
@ -72,7 +72,7 @@ bool pFlow::rotatingAxis::read
if(!timeInterval::read(dict))return false; if(!timeInterval::read(dict))return false;
if(!line::read(dict)) return false; if(!line::read(dict)) return false;
real omega = dict.getValOrSet("omega", 0.0); real omega = dict.getValOrSet("omega", static_cast<real>(0.0));
setOmega(omega); setOmega(omega);
return true; return true;

View File

@ -19,40 +19,63 @@ Licence:
-----------------------------------------------------------------------------*/ -----------------------------------------------------------------------------*/
#include "multiRotatingAxisMotion.hpp" #include "multiRotatingAxisMotion.hpp"
#include "dictionary.hpp"
#include "vocabs.hpp"
void pFlow::multiRotatingAxisMotion::impl_setTime
bool pFlow::multiRotatingAxisMotion::readDictionary
( (
const dictionary& dict uint32 iter,
) real t,
real dt
)const
{ {
auto motion = motionComponents_.deviceViewAll();
Kokkos::parallel_for(
"multiRotatingAxisMotion::impl_setTime",
deviceRPolicyStatic(0, numComponents_),
LAMBDA_D(uint32 i){
motion[i].setTime(t);
});
Kokkos::fence();
}
auto motionModel = dict.getVal<word>("motionModel"); bool pFlow::multiRotatingAxisMotion::impl_move(uint32 iter, real t , real dt ) const
{
auto motion = motionComponents_.deviceViewAll();
Kokkos::parallel_for(
"multiRotatingAxisMotion::impl_move",
deviceRPolicyStatic(0, numComponents_),
LAMBDA_D(uint32 i){
motion[i].move(dt);
});
Kokkos::fence();
return true;
}
if(motionModel != "multiRotatingAxisMotion") bool pFlow::multiRotatingAxisMotion::impl_readDictionary(const dictionary &dict)
{
auto modelName = dict.getVal<word>("motionModel");
if(modelName != getTypeName<ModelComponent>())
{ {
fatalErrorInFunction<< fatalErrorInFunction<<
" motionModel should be multiRotatingAxisMotion, but found " " motionModel should be "<< Yellow_Text(getTypeName<ModelComponent>())<<
<< motionModel <<endl; ", but found "<< Yellow_Text(modelName)<<endl;
return false; return false;
} }
auto& motionInfo = dict.subDict("multiRotatingAxisMotionInfo"); auto& motionInfo = dict.subDict(modelName+"Info");
auto axisNames = motionInfo.dictionaryKeywords(); auto compNames = motionInfo.dictionaryKeywords();
wordList rotationAxis;
// first check if wordList rotationAxisNames;
for(auto& aName: axisNames) // in the first round read all dictionaries
for(auto& compName: compNames)
{ {
auto& axDict = motionInfo.subDict(aName); auto& axDict = motionInfo.subDict(compName);
if(auto axPtr = makeUnique<rotatingAxis>(axDict); axPtr) if(auto axPtr = makeUnique<rotatingAxis>(axDict); axPtr)
{ {
rotationAxis.push_back( rotationAxisNames.push_back(
axDict.getValOrSet<word>("rotationAxis", "none")); axDict.getValOrSet<word>("rotationAxis", "none"));
} }
else else
@ -63,26 +86,26 @@ bool pFlow::multiRotatingAxisMotion::readDictionary
} }
} }
if( !axisNames.search("none") ) if( !compNames.search("none") )
{ {
axisNames.push_back("none"); compNames.push_back("none");
rotationAxis.push_back("none"); rotationAxisNames.push_back("none");
} }
using intPair = std::pair<int32, int32>; using intPair = std::pair<int32, int32>;
std::vector<intPair> numRotAxis; std::vector<intPair> numRotAxis;
for(size_t i=0; i< axisNames.size(); i++) for(size_t i=0; i< compNames.size(); i++)
{ {
word rotAxis = rotationAxis[i]; word rotAxis = rotationAxisNames[i];
int32 n=0; int32 n=0;
while(rotAxis != "none") while(rotAxis != "none")
{ {
n++; n++;
if(int32 iAxis = axisNames.findi(rotAxis) ; iAxis != -1) if(int32 iAxis = compNames.findi(rotAxis) ; iAxis != -1)
{ {
rotAxis = rotationAxis[iAxis]; rotAxis = rotationAxisNames[iAxis];
}else }else
{ {
fatalErrorInFunction<< fatalErrorInFunction<<
@ -98,60 +121,73 @@ bool pFlow::multiRotatingAxisMotion::readDictionary
auto compareFunc = [](const intPair& a, const intPair& b) auto compareFunc = [](const intPair& a, const intPair& b)
{ return a.first > b.first; }; { return a.first > b.first; };
algorithms::STD::sort(numRotAxis.data(), numRotAxis.size(), compareFunc); std::sort(numRotAxis.begin(), numRotAxis.end(), compareFunc);
Vector<int> sortedIndex;
sortedIndex_.clear(); componentNames_.clear();
axisName_.clear();
output<<compNames<<endl;
for(auto ax:numRotAxis) for(auto ax:numRotAxis)
{ {
axisName_.push_back(axisNames[ax.second]); componentNames_.push_back(compNames[ax.second]);
sortedIndex_.push_back(ax.second); sortedIndex.push_back(ax.second);
} }
numAxis_ = axisName_.size(); numComponents_ = componentNames_.size();
axis_.clear(); motionComponents_.reserve(numComponents_);
axis_.reserve(numAxis_); sortedIndex_.assign(sortedIndex);
Vector<ModelComponent> components("Read::modelComponent",
compNames.size()+1,
0,
RESERVE());
// create the actual axis vector for(auto& compName: componentNames_)
for(auto& aName: axisName_)
{ {
if(aName != "none")
if(compName != "none")
{ {
auto& axDict = motionInfo.subDict(aName); auto& compDict = motionInfo.subDict(compName);
axis_.push_back( components.emplace_back(
multiRotatingAxis(this, axDict)); motionComponents_.data(),
componentNames_,
compDict);
} }
else else
{ {
axis_.push_back( components.emplace_back(impl_noneComponent());
multiRotatingAxis(this));
} }
} }
return true; motionComponents_.assign(components);
return true;
} }
bool pFlow::multiRotatingAxisMotion::writeDictionary
bool pFlow::multiRotatingAxisMotion::impl_writeDictionary
( (
dictionary& dict dictionary& dict
)const )const
{ {
dict.add("motionModel", "multiRotatingAxisMotion"); word modelName = "multiRotatingAxis";
auto& motionInfo = dict.subDictOrCreate("multiRotatingAxisMotionInfo"); dict.add("motionModel", modelName );
ForAll(i, axis_) auto modelDictName = modelName+"Info";
auto& motionInfo = dict.subDictOrCreate(modelDictName);
auto hostComponents = motionComponents_.hostView();
ForAll(i, motionComponents_)
{ {
auto& axDict = motionInfo.subDictOrCreate(axisName_[i]); auto& axDict = motionInfo.subDictOrCreate(componentNames_[i]);
if( !axis_.hostVectorAll()[i].write(this,axDict)) if( !hostComponents[i].write(axDict, componentNames_))
{ {
fatalErrorInFunction<< fatalErrorInFunction<<
" error in writing axis "<< axisName_[i] << " to dicrionary " " error in writing axis "<< componentNames_[i] << " to dicrionary "
<< motionInfo.globalName()<<endl; << motionInfo.globalName()<<endl;
return false; return false;
} }
@ -160,79 +196,52 @@ bool pFlow::multiRotatingAxisMotion::writeDictionary
return true; return true;
} }
pFlow::multiRotatingAxisMotion::multiRotatingAxisMotion() pFlow::multiRotatingAxisMotion::multiRotatingAxisMotion(
{} const objectFile &objf,
repository *owner)
pFlow::multiRotatingAxisMotion::multiRotatingAxisMotion : fileDictionary(objf, owner)
(
const dictionary& dict
)
{ {
if(! readDictionary(dict) )
if(! getModel().impl_readDictionary(*this) )
{ {
fatalErrorInFunction;
fatalExit; fatalExit;
} }
} }
FUNCTION_H pFlow::multiRotatingAxisMotion::multiRotatingAxisMotion
bool pFlow::multiRotatingAxisMotion::move(real t, real dt)
{
// every thing is done on host
for(int32 i=0; i<numAxis_; i++)
{
auto& ax = axis_[sortedIndex_[i]];
ax.setTime(t);
ax.setAxisList(getAxisListPtrHost());
ax.move(dt);
}
// transfer to device
axis_.modifyOnHost();
axis_.syncViews();
return true;
}
bool pFlow::multiRotatingAxisMotion::read
( (
iIstream& is const objectFile &objf,
const dictionary &dict,
repository *owner
) )
:
fileDictionary(objf, dict, owner)
{ {
// create an empty file dictionary if(!getModel().impl_readDictionary(*this) )
dictionary motionInfo(motionModelFile__, true);
// read dictionary from stream
if( !motionInfo.read(is) )
{ {
ioErrorInFile(is.name(), is.lineNumber()) << fatalErrorInFunction;
" error in reading dictionray " << motionModelFile__ <<" from file. \n"; fatalExit;
return false;
} }
if( !readDictionary(motionInfo) ) return false;
return true;
} }
bool pFlow::multiRotatingAxisMotion::write bool pFlow::multiRotatingAxisMotion::write
( (
iOstream& os iOstream &os,
)const const IOPattern &iop
) const
{ {
// create an empty file dictionary // a global dictionary
dictionary motionInfo(motionModelFile__, true); dictionary newDict(fileDictionary::dictionary::name(), true);
if( iop.thisProcWriteData() )
if( !writeDictionary(motionInfo))
{ {
return false; if( !getModel().impl_writeDictionary(newDict) ||
!newDict.write(os))
{
fatalErrorInFunction<<
" error in writing to dictionary "<< newDict.globalName()<<endl;
return false;
}
} }
return true;
if( !motionInfo.write(os) )
{
ioErrorInFile( os.name(), os.lineNumber() )<<
" error in writing dictionray to file. \n";
return false;
}
return true;
} }

View File

@ -22,18 +22,16 @@ Licence:
#define __multiRotatingAxisMotion_hpp__ #define __multiRotatingAxisMotion_hpp__
#include "types.hpp" #include "MotionModel.hpp"
#include "typeInfo.hpp"
#include "VectorDual.hpp"
#include "List.hpp"
#include "multiRotatingAxis.hpp" #include "multiRotatingAxis.hpp"
#include "fileDictionary.hpp"
namespace pFlow namespace pFlow
{ {
// forward // forward
class dictionary; // class dictionary;
/** /**
* Rotating axis motion model for walls * Rotating axis motion model for walls
@ -63,200 +61,55 @@ multiRotatingAxisMotionInfo
* *
*/ */
class multiRotatingAxisMotion class multiRotatingAxisMotion
:
public fileDictionary,
public MotionModel<multiRotatingAxisMotion, multiRotatingAxis>
{ {
public:
/** Motion model class to be passed to computational units/kernels for
* transfing points and returning velocities at various positions
*/
class Model
{
protected:
deviceViewType1D<multiRotatingAxis> axis_;
int32 numAxis_=0;
public:
INLINE_FUNCTION_HD
Model(deviceViewType1D<multiRotatingAxis> 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].pointTangentialVel(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 axis_[n].transferPoint(p, dt);
}
INLINE_FUNCTION_HD int32 numComponents()const
{
return numAxis_;
}
};
protected: protected:
using axisVector_HD = VectorDual<multiRotatingAxis>; VectorSingle<int32> sortedIndex_;
/// Vector of multiRotaingAxis axes friend MotionModel<multiRotatingAxisMotion, multiRotatingAxis>;
axisVector_HD axis_;
/// Sorted index based on number of parrents each axis ha /// is the geometry attached to this component moving
VectorDual<int32> sortedIndex_; bool impl_isMoving()const
{
return true;
}
/// List of axes names
wordList axisName_;
/// Number of axes
label numAxis_= 0;
/// Read from a dictionary /// Read from dictionary
bool readDictionary(const dictionary& dict); bool impl_readDictionary(const dictionary& dict);
/// Write to a dictionary bool impl_writeDictionary(dictionary& dict)const;
bool writeDictionary(dictionary& dict)const;
public: public:
/// Type info TypeInfo("multiRotatingAxisMotion");
TypeInfoNV("multiRotatingAxisMotion");
// - Constructor multiRotatingAxisMotion(const objectFile& objf, repository* owner);
/// Empty constructor multiRotatingAxisMotion(
FUNCTION_H const objectFile& objf,
multiRotatingAxisMotion(); const dictionary& dict,
repository* owner);
/// Construct with dictionary using fileDictionary::write;
FUNCTION_H
multiRotatingAxisMotion(const dictionary& dict);
/// Copy constructor bool write(iOstream& os, const IOPattern& iop)const override;
FUNCTION_H
multiRotatingAxisMotion(const multiRotatingAxisMotion&) = default;
/// No Move static
multiRotatingAxisMotion(multiRotatingAxisMotion&&) = delete; multiRotatingAxis noneComponent()
{
/// Copy assignment return multiRotatingAxis();
FUNCTION_H }
multiRotatingAxisMotion& operator=(const multiRotatingAxisMotion&) = default;
/// No move assignment
multiRotatingAxisMotion& operator=(multiRotatingAxisMotion&&) = delete;
/// Destructor
FUNCTION_H
~multiRotatingAxisMotion() = default;
// - Methods
/// Retrun motion model at time t
Model getModel(real t)
{
for(int32 i= 0; i<numAxis_; i++ )
{
axis_[i].setTime(t);
axis_[i].setAxisList(getAxisListPtrDevice());
}
axis_.modifyOnHost();
axis_.syncViews();
return Model(axis_.deviceVector(), numAxis_);
}
/// Pointer to axis list on host side
INLINE_FUNCTION_H
multiRotatingAxis* getAxisListPtrHost()
{
return axis_.hostVectorAll().data();
}
/// Pointer to axis list on device
INLINE_FUNCTION_H
multiRotatingAxis* getAxisListPtrDevice()
{
return axis_.deviceVectorAll().data();
}
/// Name of motion component to index
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;
}
}
/// Index of motion component to component name
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 <<endl<<
" size of axes_ is "<<numAxis_<<endl;
fatalExit;
return "";
}
}
/// Is moving
INLINE_FUNCTION_HD
bool isMoving()const
{
return true;
}
/// Move points
FUNCTION_H
bool move(real t, real dt);
// - IO operation
/// Read from input stream is
FUNCTION_H
bool read(iIstream& is);
/// Write to output stream os
FUNCTION_H
bool write(iOstream& os)const;
// TODO: make this method protected
void impl_setTime(uint32 iter, real t, real dt)const;
/// move the component itself
bool impl_move(uint32 iter, real t, real dt)const;
}; };
} // pFlow } // pFlow

View File

@ -29,7 +29,7 @@ pFlow::stationaryWall::stationaryWall
: :
fileDictionary(objf, owner) fileDictionary(objf, owner)
{ {
const auto& dummy = this->subDictOrCreate("stationaryInfo");
if(!impl_readDictionary(*this) ) if(!impl_readDictionary(*this) )
{ {
fatalErrorInFunction; fatalErrorInFunction;
@ -46,6 +46,8 @@ pFlow::stationaryWall::stationaryWall
: :
fileDictionary(objf, dict, owner) fileDictionary(objf, dict, owner)
{ {
const auto& dummy = this->subDictOrCreate("stationaryInfo");
if(!impl_readDictionary(*this) ) if(!impl_readDictionary(*this) )
{ {
fatalErrorInFunction; fatalErrorInFunction;

View File

@ -6,6 +6,9 @@ particles/shape/shape.cpp
particles/particles.cpp particles/particles.cpp
particles/particleIdHandler/particleIdHandler.cpp particles/particleIdHandler/particleIdHandler.cpp
particles/regularParticleIdHandler/regularParticleIdHandler.cpp particles/regularParticleIdHandler/regularParticleIdHandler.cpp
globalDamping/globalDamping.cpp
SphereParticles/sphereShape/sphereShape.cpp SphereParticles/sphereShape/sphereShape.cpp
SphereParticles/sphereParticles/sphereParticles.cpp SphereParticles/sphereParticles/sphereParticles.cpp
SphereParticles/sphereParticles/sphereParticlesKernels.cpp SphereParticles/sphereParticles/sphereParticlesKernels.cpp

View File

@ -50,9 +50,7 @@ public:
const grainParticles& Particles()const; const grainParticles& Particles()const;
bool hearChanges( bool hearChanges(
real t, const timeInfo& ti,
real dt,
uint32 iter,
const message &msg, const message &msg,
const anyList &varList) override const anyList &varList) override
{ {
@ -65,6 +63,11 @@ public:
return true; return true;
} }
bool isActive()const override
{
return false;
}
static static
uniquePtr<boundaryGrainParticles> create( uniquePtr<boundaryGrainParticles> create(
const boundaryBase &boundary, const boundaryBase &boundary,

View File

@ -5,10 +5,10 @@ pFlow::boundaryGrainParticlesList::boundaryGrainParticlesList(
grainParticles &prtcls grainParticles &prtcls
) )
: :
ListPtr(bndrs.size()), boundaryListPtr(),
boundaries_(bndrs) boundaries_(bndrs)
{ {
for(auto i=0; i<boundaries_.size(); i++) ForAllBoundariesPtr(i, this)
{ {
this->set this->set
( (
@ -16,4 +16,5 @@ pFlow::boundaryGrainParticlesList::boundaryGrainParticlesList(
boundaryGrainParticles::create(boundaries_[i], prtcls) boundaryGrainParticles::create(boundaries_[i], prtcls)
); );
} }
} }

View File

@ -3,7 +3,7 @@
#ifndef __boundaryGrainParticlesList_hpp__ #ifndef __boundaryGrainParticlesList_hpp__
#define __boundaryGrainParticlesList_hpp__ #define __boundaryGrainParticlesList_hpp__
#include "ListPtr.hpp" #include "boundaryListPtr.hpp"
#include "boundaryList.hpp" #include "boundaryList.hpp"
#include "boundaryGrainParticles.hpp" #include "boundaryGrainParticles.hpp"
@ -12,7 +12,7 @@ namespace pFlow
class boundaryGrainParticlesList class boundaryGrainParticlesList
: :
public ListPtr<boundaryGrainParticles> public boundaryListPtr<boundaryGrainParticles>
{ {
private: private:

View File

@ -142,16 +142,11 @@ pFlow::grainParticles::getParticlesInfoFromShape(
pFlow::grainParticles::grainParticles( pFlow::grainParticles::grainParticles(
systemControl &control, systemControl &control,
const property& prop const grainShape& gShape
) )
: :
particles(control), particles(control, gShape),
grains_ grains_(gShape),
(
shapeFile__,
&control.caseSetup(),
prop
),
propertyId_ propertyId_
( (
objectFile objectFile
@ -253,7 +248,7 @@ pFlow::grainParticles::grainParticles(
"rVelocity", "rVelocity",
dynPointStruct(), dynPointStruct(),
intMethod, intMethod,
rVelocity_.field() rAcceleration_.field()
); );
if( !rVelIntegration_ ) if( !rVelIntegration_ )
@ -280,7 +275,7 @@ bool pFlow::grainParticles::beforeIteration()
particles::beforeIteration(); particles::beforeIteration();
intPredictTimer_.start(); intPredictTimer_.start();
auto dt = this->dt(); auto dt = this->dt();
dynPointStruct().predict(dt, accelertion()); dynPointStruct().predict(dt);
rVelIntegration_().predict(dt,rVelocity_, rAcceleration_); rVelIntegration_().predict(dt,rVelocity_, rAcceleration_);
intPredictTimer_.end(); intPredictTimer_.end();
@ -301,8 +296,8 @@ bool pFlow::grainParticles::beforeIteration()
bool pFlow::grainParticles::iterate() bool pFlow::grainParticles::iterate()
{ {
timeInfo ti = TimeInfo(); const timeInfo ti = TimeInfo();
realx3 g = control().g(); const realx3 g = control().g();
particles::iterate(); particles::iterate();
accelerationTimer_.start(); accelerationTimer_.start();
@ -313,25 +308,28 @@ bool pFlow::grainParticles::iterate()
I().deviceViewAll(), I().deviceViewAll(),
contactTorque().deviceViewAll(), contactTorque().deviceViewAll(),
dynPointStruct().activePointsMaskDevice(), dynPointStruct().activePointsMaskDevice(),
accelertion().deviceViewAll(), acceleration().deviceViewAll(),
rAcceleration().deviceViewAll() rAcceleration().deviceViewAll()
); );
for(auto& bndry:boundaryGrainParticles_)
ForAllActiveBoundaries(i,boundaryGrainParticles_)
{ {
bndry->acceleration(ti, g); boundaryGrainParticles_[i].acceleration(ti, g);
} }
accelerationTimer_.end(); accelerationTimer_.end();
intCorrectTimer_.start(); intCorrectTimer_.start();
if(!dynPointStruct().correct(dt(), accelertion())) if(!dynPointStruct().correct(ti.dt()))
{ {
return false; return false;
} }
real damping = dynPointStruct().dampingFactor(ti);
if(!rVelIntegration_().correct( if(!rVelIntegration_().correct(
dt(), ti.dt(),
rVelocity_, rVelocity_,
rAcceleration_)) rAcceleration_,
damping))
{ {
return false; return false;
} }

View File

@ -48,7 +48,7 @@ public:
private: private:
/// reference to shapes /// reference to shapes
ShapeType grains_; const ShapeType& grains_;
/// property id on device /// property id on device
uint32PointField_D propertyId_; uint32PointField_D propertyId_;
@ -121,7 +121,7 @@ protected:
public: public:
/// construct from systemControl and property /// construct from systemControl and property
grainParticles(systemControl& control, const property& prop); grainParticles(systemControl& control, const grainShape& gShape);
~grainParticles() override = default; ~grainParticles() override = default;
@ -172,9 +172,7 @@ public:
} }
bool hearChanges( bool hearChanges(
real t, const timeInfo& ti,
real dt,
uint32 iter,
const message& msg, const message& msg,
const anyList& varList const anyList& varList
) override ) override

View File

@ -81,7 +81,7 @@ pFlow::insertion::readInsertionDict()
if (active_) if (active_)
{ {
checkForCollision_ = getVal<Logical>("checkForCollision"); //checkForCollision_ = getVal<Logical>("checkForCollision");
REPORT(1) << "Particle insertion mechanism is " << Yellow_Text("active") REPORT(1) << "Particle insertion mechanism is " << Yellow_Text("active")
<< " in the simulation." << END_REPORT; << " in the simulation." << END_REPORT;

View File

@ -42,7 +42,7 @@ private:
Logical active_ = "No"; Logical active_ = "No";
/// Check for collision? It is not active now /// Check for collision? It is not active now
Logical checkForCollision_ = "No"; Logical checkForCollision_ = "Yes";
/// if increase velocity in case particles are failed /// if increase velocity in case particles are failed
/// to be inserted due to collision /// to be inserted due to collision

View File

@ -24,6 +24,7 @@ Licence:
#include "particles.hpp" #include "particles.hpp"
#include "twoPartEntry.hpp" #include "twoPartEntry.hpp"
#include "types.hpp" #include "types.hpp"
#include "shape.hpp"
namespace pFlow namespace pFlow
{ {

View File

@ -173,7 +173,7 @@ public:
inline bool insertionTime(uint32 iter, real t, real dt) const inline bool insertionTime(uint32 iter, real t, real dt) const
{ {
return tControl_.timeEvent(iter, t, dt); return tControl_.eventTime(iter, t, dt);
} }
uint32 numberToBeInserted(uint32 iter, real t, real dt); uint32 numberToBeInserted(uint32 iter, real t, real dt);

View File

@ -50,9 +50,7 @@ public:
const sphereParticles& Particles()const; const sphereParticles& Particles()const;
bool hearChanges( bool hearChanges(
real t, const timeInfo& ti,
real dt,
uint32 iter,
const message &msg, const message &msg,
const anyList &varList) override const anyList &varList) override
{ {
@ -65,6 +63,11 @@ public:
return true; return true;
} }
bool isActive()const override
{
return false;
}
static static
uniquePtr<boundarySphereParticles> create( uniquePtr<boundarySphereParticles> create(
const boundaryBase &boundary, const boundaryBase &boundary,

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