131 Commits

Author SHA1 Message Date
7eb707f950 Merge pull request #118 from PhasicFlow/develop
Develop
2024-10-19 12:27:41 +03:30
HRN
d3ccf354b7 code recovery regular part 2024-10-18 23:13:20 +03:30
HRN
173d3c4917 first commit after code loss - develop branch 2024-09-21 13:37:03 +03:30
8ee1fe63e6 Merge pull request #114 from ramin1728/main
binary system of particles tutorial is updated
2024-05-28 13:33:29 +03:30
ca1bc218d2 Update README.md 2024-05-28 13:26:34 +03:30
c8a28f82c9 layeredSiloFilling updated 2024-05-28 13:24:34 +04:30
b591d8fda1 binarySystemOfParticles updated 2024-05-28 13:10:19 +04:30
2f9db61983 The Tutorial is updated for version 1.x. 2024-05-28 12:47:31 +04:30
159ec5f1ff Tutorial is updated for version 1.x. 2024-05-28 12:25:10 +04:30
afa7d63c84 Tutorial is updated for version 1.x. 2024-05-28 12:23:05 +04:30
7b2c7c6199 Tutorial is updated for 1.x. 2024-05-28 12:19:37 +04:30
7b6271a165 Tutorial is updated for version 1.x. 2024-05-28 12:15:45 +04:30
526059707b Tutorials are updated. 2024-05-28 12:09:52 +04:30
180bc29df4 Tutorials are updated. 2024-05-28 12:08:25 +04:30
569646c6ac Tutorials are updated 2024-05-28 12:03:10 +04:30
3bbad4487b Merge pull request #113 from PhasicFlow/develop
Develop
2024-05-24 22:17:41 +03:30
HRN
6a66f1edfd bug fix for binary file read when dealing with multiple Fields reading from a single file 2024-05-24 22:02:46 +03:30
HRN
32a2c20613 for write to file binary MPI 2024-05-24 20:50:53 +03:30
HRN
d1b6a04292 makeing beforeIteration step-wise for boundaryBase to implement non-bloking communications 2024-05-24 00:10:49 +03:30
4f9f518b8d Merge pull request #110 from ramin1728/main
binarySystem tutorial
2024-05-22 10:00:24 +03:30
a05cd454a6 Merge pull request #112 from PhasicFlow/develop
Develop
2024-05-22 10:00:02 +03:30
d5ff1d7906 modifications for MPI+CUDA run 2024-05-22 09:56:06 +03:30
5e56671c05 layeredSiloFilling tutorial is updated for phasicFlow version 1.x. 2024-05-21 19:43:58 +04:30
afa790b04d binarySystemOfParticles 2024-05-21 16:59:15 +04:30
74833ce6a6 This toturial has been modified for version 1.x of phasicFlow. 2024-05-21 13:30:16 +04:30
888343c655 Tutorial binarySystemOfParticles are modified based on version 1.x 2024-05-19 20:50:00 +04:30
4e8b921514 bug remove for GPU run after CPU MPI parallelization
- specialization of VectorSingle for word
- dummyFile creation to solve write to file in MPI mode
2024-05-18 18:40:25 +03:30
8eba161d62 Merge pull request #109 from PhasicFlow/develop
develop branch update for changes in MPI branch for data transfer stage
2024-05-12 19:12:40 +03:30
HRN
614b2f732e develop branch update for changes in MPI branch for data transfer stage
mortong indexing is added (messaging is not complete)
2024-05-12 19:10:04 +03:30
99bb58bda7 Merge pull request #108 from PhasicFlow/develop
boundaryProcessor for dataTransfer
2024-05-05 23:10:08 +03:30
HRN
0e54e260e6 boundaryProcessor for dataTransfer
- in boundaryList in afterIteration, data transfer occurs to transfer data between prcocessors
- transferData method is added to boundaryBase
- Some modification to format the output messages in terminal
2024-05-05 23:05:23 +03:30
6a739e79b8 Merge pull request #107 from PhasicFlow/develop
Modification in streams and sphereInteraciton
2024-05-05 23:04:33 +03:30
HRN
97049760c3 Modification in streams and sphereInteraciton
1. processorOutput stream is modified to imporve visualization of output
2. some modifications for sphereInteraction after solving MPI memory leak.
2024-04-30 00:32:18 +03:30
9a69f510a5 Merge pull request #106 from PhasicFlow/develop
modification to the base code for MPI-developement up to processor bo…
2024-04-30 00:25:19 +03:30
e7215159af modification to the base code for MPI-developement up to processor boundaries 2024-04-27 09:03:19 -07:00
04fb9c8c36 Merge pull request #105 from PhasicFlow/develop
dataIO modified for virtual constructor
2024-04-24 18:06:06 +03:30
5c65b92008 dataIO modified for virtual constructor
- In new version dataIO uses virtual constructor to create the object for regular or MPI data transfer
- The instantiation is required for new data types
2024-04-20 06:11:57 -07:00
19caa17642 Merge pull request #104 from PhasicFlow/develop
zAxis is upgraded and array2D is added for support for 2D array opera…
2024-04-19 19:28:27 +03:30
e4c387a801 git ignore and CMAKE 2024-04-18 10:06:07 -07:00
ce57a18c4f Merge pull request #103 from PhasicFlow/develop
Particle insertion upgrade and pointStructure selectors
2024-04-18 00:37:03 +03:30
c432880689 zAxis is upgraded and array2D is added for support for 2D array operations
- Cylinder geometric region is corrected for minPoint and maxPoint
2024-04-17 14:04:34 -07:00
79b987c711 Particle insertion
- particleInsertion is now transfered to time folders to track particle insertion progress
- collision check is added for new inserted particles and old particles in the simulaiton.
- minor changes to dataEntry and triple
2024-04-15 13:19:03 -07:00
edcf2907e6 Merge pull request #102 from PhasicFlow/develop
pointStruture selectors refactored
2024-04-15 23:48:38 +03:30
d6798b5097 pointStructure selector
- selectors for pointStructure are updated and geometric selector is added. With this, any geometric region can also be used as a selector region for particles.
- Also, a pointField selector is added as a free funtion that can be used anywhere that is needed.
2024-04-15 13:14:29 -07:00
07bb9e9003 pointStruture selectors refactored
- geometricRegion was added as selector(for now, box, sphere, and cylinder)
- base class upgraded and new method is added to allow accessing point positions in a vector
- selection is modified to account for inactive points when selection is done based on position
- file interface for particlesPhasicFlow changed.

- TODO
- Tutorials should be updated in new version.
2024-04-13 13:13:53 -07:00
e4cc10fbd7 Merge pull request #101 from PhasicFlow/develop
Collision check and particlePosition -> random
2024-04-13 23:43:05 +03:30
e395c379cb Collision check and particlePosition -> random
- A new class is added for simple collision check
- position particles in utility is upgraded
- morton sorting is not active yet for particlesPhasicFlow
2024-04-13 07:07:36 -07:00
71f2649659 Merge pull request #100 from PhasicFlow/develop
Particle insertion is added with anyList
2024-04-13 12:13:17 +03:30
89d7e1f0ba Particle insertion (domain check)
- Domain check is added to prevent points that are outside of internalDomainBox
2024-04-13 01:41:11 -07:00
9c2a9a81b0 Particle insertion is added with anyList
- collision check is not active yet.
- variable velocity is not active yet.
- events and messages are not active yet.
2024-04-12 22:39:09 -07:00
70c1afceec Merge pull request #99 from PhasicFlow/develop
modifications to typeInfo.hpp and type system
2024-04-12 00:23:19 +03:30
97f0ddf82e modifications to typeInfo.hpp and type system
- ClassInfo is added to support for non-polymorphic typename
- removal template specialization for realx3 and etc.
- addition of new TypeInfo for triple and quadruple
- change in Logial TypeInfo to remove recursive include in some files
2024-04-11 13:48:37 -07:00
fb151f323c Merge pull request #97 from PhasicFlow/develop
Reflective boundary condition is added and tested.
2024-04-12 00:13:42 +03:30
ff8968e595 format clang-format 2024-04-08 12:33:08 -07:00
a9290e911b Reflective boundary condition is added and tested.
It requires messaging integration when changing velocity.
2024-04-05 05:31:09 -07:00
9336498e38 Merge pull request #96 from PhasicFlow/develop
Develop - Boundary conditions are created and tested
2024-04-04 23:07:16 +03:30
821dde9b1c modifications for adding boundary condtions 2024-04-04 12:33:09 -07:00
818106a34a Boundary conditions are created and tested.
- exit and periodic, follows the previous commit.
2024-04-04 12:31:19 -07:00
65b3f5a8a0 phasicFlow core changes for boundary condition addition.
- Exit boundry condition, checked.
- periodic boundary condition added, checked.

to be done:
- messaging events
- reflective
- wall-particle interaction for periodic bc.
2024-04-04 10:36:11 -07:00
89c7be0a20 Merge pull request #95 from PhasicFlow/develop
pFlowToVTK enhanced
2024-03-30 19:08:56 +03:30
ef0e752929 pFlowToVTK finalized and tested
- geometry conversion now can handle separate vtk files.
2024-03-30 08:35:16 -07:00
815b134e1e pFlowToVTK enhanced
- pointField conversion is completed, tested (finalized).
- geometry conversion is not complete yet.
2024-03-30 01:42:51 -07:00
cab49ef036 Merge pull request #94 from PhasicFlow/develop
pFlowToVTK and enhanced wall-particle search
2024-03-30 00:24:04 +03:30
f288f812fd pFlowToVTK is added
- this utility is not complete yet.
- geometry should be ajusted to be converted separately
2024-03-29 13:50:02 -07:00
19d1ad343d changes to cellsWallLevel0 to narrow down the contact list for particle-wall 2024-03-29 13:45:29 -07:00
7be70caa82 Merge pull request #93 from PhasicFlow/develop
Periodic boundary condition implementation
2024-03-26 20:39:43 +03:30
596ebb0262 Periodic boundary condition implementation
- boundaryBase class, some methods has been added, index list on host added, boundaryList friended
- generalBoundary class is created for all boundary types
- searchBoundary class is created for contact search -> (regular, none),
- pointStructure can be accessed through contactSearch

to be done:
- A new class for periodic boundary condtion should be added. This class should handel pp and pw contact searchs.
- Messages should be rearranged for boundaries.
- Events and updates should be synced, so that the same contact list can be used for internal particles and boundary particles.
2024-03-26 09:05:02 -07:00
af75d6a09c Merge pull request #92 from PhasicFlow/develop
Merge (first patch of version 1.x)
2024-03-25 01:03:20 +03:30
c7e5b30959 no changes were made, just debug 2024-03-24 14:26:41 -07:00
a900ff32bb phase1 of periodict boundary condition (serial)
upto now, particles are transfered from one boundary to another (mirro). this has been tested using iterateSphereParticles utility.
2024-03-24 14:25:03 -07:00
593cb9e6ed corrections for methods assign and append
The size change and capacity check was not correct.
A full test should be done on all the methods in this class
2024-03-24 14:22:48 -07:00
38f1c8ec59 change order of files in CMaleLists.txt
this is done for better compilation order and speeding up the compilation
2024-03-24 14:21:29 -07:00
3af1420ff8 Correction for mapperNBS 2024-03-24 03:57:48 -07:00
4c36418ec8 solver sphereGranFlow without particles insertion 2024-03-24 03:57:16 -07:00
3ada371198 iterateGeometry is finalized for version (1.x) 2024-03-24 03:01:08 -07:00
90fc1c6357 particlesPhasicFlow finalized for version (1.x) 2024-03-24 02:59:51 -07:00
93b33deffb finalized gemometryPhasicFlow for version (1.x) 2024-03-24 02:56:37 -07:00
8beaec448a modifications for version (1.0) 2024-03-24 02:49:47 -07:00
6b3a4f017b structured data has been modified for new changes in version (1.0) 2024-03-24 02:02:48 -07:00
22405db699 Trisurface modified with new data structure
This commit is based on changes to geometry class of phasicFlow (version 1.0)
2024-03-24 02:01:02 -07:00
57850119ba containers are modified for errors
This commit is before the first release of phasicFlow from develop to main (version 1.0)
2024-03-24 01:56:49 -07:00
0aa7f2dba9 Integration AB2 2024-03-24 01:56:04 -07:00
d21b7276e1 Interaction folder is finalized with adjustable box and tested for rotating drum 2024-03-24 01:46:19 -07:00
be56d8ee2e MotionModel folder completed 2024-03-22 09:42:23 -07:00
f27fbdd82c message event update for points and boundaries 2024-03-22 09:01:07 -07:00
c3821d03c0 minor correction for localProcessors, localRank and localSize 2024-03-22 08:58:01 -07:00
49af1119f9 particles and sphereParticles classes were tested for exit boundary 2024-03-22 08:26:14 -07:00
9facc73055 Update README.md 2024-03-20 11:48:25 +03:30
acfaacfe4a Geometry folder is finalized 2024-03-20 00:50:58 -07:00
df8bea213d triSurface folder content is finalized 2024-03-19 22:29:31 -07:00
c13231f4f0 Correction for non-trivially-copyable data type on host memory 2024-03-19 07:33:56 -07:00
c0cf200fa8 iterateGeometry and geometryPhasicFlow refactor 2024-02-05 21:28:30 -08:00
5b4a524afe Motion models integrated into geometryMotion, stationaryWall is added 2024-02-05 21:27:24 -08:00
9dfe98eea2 triSurfaceField refactored 2024-02-05 01:29:45 -08:00
80b61d4d73 MotionModel CRTP, rotatingAxis and vibrating 2024-02-04 23:58:57 -08:00
fd039f234f geometryPhasicFlow-pass1 and triSurface and multiTriSurface tested 2024-02-03 11:49:07 -08:00
e10ee2b6b5 periodic boundary condition is added to pointStructure 2024-02-01 12:32:15 -08:00
40d74c846f output message has been removed 2024-01-30 22:51:42 -08:00
5b9b865b08 solver iterateSphereParticles is added for testing (without particle insertion) 2024-01-29 07:58:19 -08:00
182e156786 sphereParticles tested on CPU, iteration, write to file, particle deletion 2024-01-29 07:57:19 -08:00
c0ee29e39c boundaryExit beforeIteration, not tested 2024-01-28 14:43:10 -08:00
0df384f546 sphereParticles tested 2024-01-26 01:10:10 -08:00
206df8924e utility particlesPhasicFlow updated. TODO: upgrade it further more to mandates supplying shapeName in particlesDict 2024-01-25 03:12:01 -08:00
ab6308bb0a code refactor upto pointFields and insertion etc. 2024-01-25 03:10:44 -08:00
8dc8009311 properry class refactored 2024-01-25 03:09:08 -08:00
20be76aed0 Particle base class and sphereParticle class have been updated 2024-01-25 03:07:59 -08:00
9c86fe8f31 refactor up to particles.hpp 2024-01-21 13:26:23 -08:00
46bf08fa91 typeInfo.hpp updated 2024-01-20 12:33:13 -08:00
6bc180275e pointFieild 2024-01-20 06:05:05 -08:00
0096856c98 Field class and pointStructure modified and tested 2024-01-19 22:31:08 -08:00
90a1ac8f16 refactor upto pointField read 2024-01-18 04:51:06 -08:00
f5d8daa608 MPI-parallelization upto IO file 2024-01-13 09:54:23 +03:30
280f53a230 before checking parallelIO for dictionary 2023-12-25 13:59:24 +03:30
f1baff5a59 pointStructure with boundaries, construction tested 2023-12-17 15:27:05 +03:30
143a586730 new-message-observer-subscriber and anyList for holding a list of variables of any type 2023-11-08 20:35:18 +03:30
d65aa02cbb Parallel IO before round 2 2023-10-28 12:09:40 +03:30
66fc880fcd math.hpp modified to remove ambiguis calls in cuda build 2023-10-19 12:18:29 +03:30
d2cd132b62 parallel read and write is added-we still need to consider integrating it with fields 2023-10-17 20:30:20 +03:30
6e5fe608c6 Kokkos folder updated 2023-10-03 11:15:09 +03:30
34e3f5587b correction gitignore 2023-10-01 20:20:18 +03:30
a14ce9f2bf corrections 2023-10-01 20:18:10 +03:30
6248903236 gitignor updated 2023-10-01 20:10:48 +03:30
5a94fc5688 dictionary is updated 2023-10-01 20:04:46 +03:30
59c148afa6 MPI part 2023-10-01 18:47:24 +03:30
8eec685c75 bashrc for Zoltan 2023-10-01 01:22:40 +03:30
100a211645 streams updated, processorsStream, masterStream 2023-09-30 20:18:17 +03:30
c562815bbe Update README.md 2023-09-28 23:30:50 +03:30
fa8044ad23 Update README.md 2023-09-28 23:29:51 +03:30
bb7dc17c67 processors is added and error part modified 2023-09-28 23:14:12 +03:30
503ee6be5e Types folder updated 2023-09-28 07:44:39 +03:30
1696340ada updated cmake for version-1.0 and including MPI 2023-09-27 11:28:39 +03:30
520 changed files with 36899 additions and 19064 deletions

6
.gitignore vendored
View File

@ -1,6 +1,12 @@
# files
.clang-format
.vscode
.dependencygraph
# Prerequisites # Prerequisites
*.d *.d
# Compiled Object files # Compiled Object files
*.slo *.slo
*.lo *.lo

View File

@ -1,74 +1,73 @@
cmake_minimum_required(VERSION 3.16 FATAL_ERROR) cmake_minimum_required(VERSION 3.16 FATAL_ERROR)
# set the project name and version # set the project name and version
project(phasicFlow VERSION 0.1 ) project(phasicFlow VERSION 1.0 )
set(CMAKE_CXX_STANDARD 17 CACHE STRING "" FORCE) set(CMAKE_CXX_STANDARD 17 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 Release CACHE STRING "build type" FORCE) set(CMAKE_BUILD_TYPE Debug CACHE STRING "build type" FORCE)
set(BUILD_SHARED_LIBS ON CACHE BOOL "Build using shared libraries" FORCE)
mark_as_advanced(FORCE var BUILD_SHARED_LIBS)
message(STATUS ${CMAKE_INSTALL_PREFIX}) message(STATUS ${CMAKE_INSTALL_PREFIX})
mark_as_advanced(FORCE var Kokkos_ENABLE_CUDA_LAMBDA) include(cmake/globals.cmake)
mark_as_advanced(FORCE var Kokkos_ENABLE_OPENMP)
mark_as_advanced(FORCE var Kokkos_ENABLE_SERIAL)
mark_as_advanced(FORCE var Kokkos_ENABLE_CUDA_LAMBDA)
mark_as_advanced(FORCE var BUILD_SHARED_LIBS)
option(USE_STD_PARALLEL_ALG "Use TTB std parallel algorithms" ON) #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} ${phasicFlow_BINARY_DIR}/kokkos)
Kokkos_cmake_settings()
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)
option(pFlow_Build_Cuda "Build phasicFlow and backends for Cuda execution" OFF) option(pFlow_Build_Cuda "Build phasicFlow and backends for Cuda execution" OFF)
option(pFlow_Build_Double "Build phasicFlow with double precision 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)
set(BUILD_SHARED_LIBS ON CACHE BOOL "Build using shared libraries" FORCE)
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)
set(Kokkos_ENABLE_OPENMP OFF CACHE BOOL "OpenMP execution" FORCE) set(Kokkos_ENABLE_OPENMP OFF CACHE BOOL "OpenMP execution" FORCE)
set(Kokkos_ENABLE_CUDA OFF CACHE BOOL "Cuda execution" FORCE) set(Kokkos_ENABLE_CUDA 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_ENABLE_CUDA_CONSTEXPR OFF CACHE BOOL "Enable constexpr on cuda code" FORCE)
elseif(pFlow_Build_OpenMP ) elseif(pFlow_Build_OpenMP )
set(Kokkos_ENABLE_SERIAL ON CACHE BOOL "Serial execution" FORCE) set(Kokkos_ENABLE_SERIAL ON CACHE BOOL "Serial execution" FORCE)
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_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)
set(Kokkos_ENABLE_OPENMP OFF CACHE BOOL "OpenMP execution" FORCE) set(Kokkos_ENABLE_OPENMP OFF CACHE BOOL "OpenMP execution" FORCE)
set(Kokkos_ENABLE_CUDA ON CACHE BOOL "Cuda execution" FORCE) set(Kokkos_ENABLE_CUDA ON CACHE BOOL "Cuda execution" FORCE)
set(Kokkos_ENABLE_CUDA_LAMBDA ON CACHE BOOL "Cuda execution" FORCE) set(Kokkos_ENABLE_CUDA_LAMBDA ON CACHE BOOL "Cuda execution" FORCE)
set(Kokkos_ENABLE_CUDA_CONSTEXPR ON CACHE BOOL "Enable constexpr on cuda code" FORCE)
endif() endif()
if(pFlow_Build_MPI)
find_package(MPI REQUIRED)
include(cmake/globals.cmake) endif()
message(STATUS "Valid file extensions are ${validFiles}")
include(cmake/makeLibraryGlobals.cmake) include(cmake/makeLibraryGlobals.cmake)
include(cmake/makeExecutableGlobals.cmake) include(cmake/makeExecutableGlobals.cmake)
configure_file(phasicFlowConfig.H.in phasicFlowConfig.H) configure_file(phasicFlowConfig.H.in phasicFlowConfig.H)
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}")
#main subdirectories of the code
set(Kokkos_Source_DIR)
if(DEFINED ENV{Kokkos_DIR})
set(Kokkos_Source_DIR $ENV{Kokkos_DIR})
# add_subdirectory($ENV{Kokkos_DIR} ${phasicFlow_BINARY_DIR}/kokkos)
# message(STATUS "Kokkos directory is $ENV{Kokkos_DIR}")
else()
# add_subdirectory($ENV{HOME}/Kokkos/kokkos ${phasicFlow_BINARY_DIR}/kokkos)
set(Kokkos_Source_DIR $ENV{HOME}/Kokkos/kokkos)
endif()
message(STATUS "Kokkos source directory is ${Kokkos_Source_DIR}")
add_subdirectory(${Kokkos_Source_DIR} ${phasicFlow_BINARY_DIR}/kokkos)
add_subdirectory(src) add_subdirectory(src)
@ -76,13 +75,11 @@ add_subdirectory(solvers)
add_subdirectory(utilities) add_subdirectory(utilities)
add_subdirectory(DEMSystems) #add_subdirectory(DEMSystems)
#add_subdirectory(testIO) add_subdirectory(testIO)
install(FILES "${PROJECT_BINARY_DIR}/phasicFlowConfig.H" install(FILES "${PROJECT_BINARY_DIR}/phasicFlowConfig.H"
DESTINATION include DESTINATION include)
)
include(InstallRequiredSystemLibraries) include(InstallRequiredSystemLibraries)
set(CPACK_RESOURCE_FILE_LICENSE "${CMAKE_CURRENT_SOURCE_DIR}/LICENSE") set(CPACK_RESOURCE_FILE_LICENSE "${CMAKE_CURRENT_SOURCE_DIR}/LICENSE")

View File

@ -5,8 +5,16 @@
**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 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.
**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.
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? ## How to build?
You can build PhasicFlow for CPU and GPU executions. [Here is a complete step-by-step procedure](https://github.com/PhasicFlow/phasicFlow/wiki/How-to-Build-PhasicFlow). 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 ## Online code documentation
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/) 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/)

View File

@ -1,5 +1,5 @@
export pFlow_PROJECT_VERSION=v0.1 export pFlow_PROJECT_VERSION=v-1.0
export pFlow_PROJECT=phasicFlow export pFlow_PROJECT=phasicFlow
@ -20,6 +20,8 @@ export pFlow_SRC_DIR="$pFlow_PROJECT_DIR/src"
export Kokkos_DIR="$kokkosDir" export Kokkos_DIR="$kokkosDir"
export Zoltan_DIR="$projectDir/Zoltan"
# Cleanup variables (done as final statement for a clean exit code) # Cleanup variables (done as final statement for a clean exit code)
unset projectDir unset projectDir
@ -27,5 +29,7 @@ export PATH="$pFlow_BIN_DIR:$PATH"
export LD_LIBRARY_PATH="$pFlow_LIB_DIR:$LD_LIBRARY_PATH" export LD_LIBRARY_PATH="$pFlow_LIB_DIR:$LD_LIBRARY_PATH"
export LD_LIBRARY_PATH="$Zoltan_DIR/lib:$LD_LIBRARY_PATH"
#------------------------------------------------------------------------------ #------------------------------------------------------------------------------

View File

@ -2,3 +2,52 @@
set(validFiles) set(validFiles)
list(APPEND validFiles *.C *.cpp *.cxx *.c *.cu *.H *.hpp *.hxx *.h *.cuh) list(APPEND validFiles *.C *.cpp *.cxx *.c *.cu *.H *.hpp *.hxx *.h *.cuh)
macro(Kokkos_cmake_settings)
#mark_as_advanced(FORCE var Kokkos_ENABLE_CUDA_LAMBDA)
mark_as_advanced(FORCE var Kokkos_CXX_STANDARD)
#mark_as_advanced(FORCE var Kokkos_ENABLE_CUDA_CONSTEXPR)
mark_as_advanced(FORCE var Kokkos_ENABLE_OPENMP)
mark_as_advanced(FORCE var Kokkos_ENABLE_SERIAL)
mark_as_advanced(FORCE var Kokkos_ENABLE_CUDA)
mark_as_advanced(FORCE var Kokkos_ENABLE_HIP)
mark_as_advanced(FORCE var Kokkos_ENABLE_AGGRESSIVE_VECTORIZATION)
mark_as_advanced(FORCE var Kokkos_ENABLE_BENCHMARKS)
mark_as_advanced(FORCE var Kokkos_ENABLE_COMPILE_AS_CMAKE_LANGUAGE)
mark_as_advanced(FORCE var Kokkos_ENABLE_CUDA_LDG_INTRINSIC)
mark_as_advanced(FORCE var Kokkos_ENABLE_CUDA_RELOCATABLE_DEVICE_CODE)
mark_as_advanced(FORCE var Kokkos_ENABLE_CUDA_UVM)
mark_as_advanced(FORCE var Kokkos_ENABLE_DEPRECATED_CODE_3)
mark_as_advanced(FORCE var Kokkos_ENABLE_DEPRECATED_CODE_4)
mark_as_advanced(FORCE var Kokkos_ENABLE_DEPRECATION_WARNINGS)
mark_as_advanced(FORCE var Kokkos_ENABLE_DESUL_ATOMICS_EXTERNAL)
mark_as_advanced(FORCE var Kokkos_ENABLE_EXAMPLES)
mark_as_advanced(FORCE var Kokkos_ENABLE_HEADER_SELF_CONTAINMENT_TESTS)
mark_as_advanced(FORCE var Kokkos_ENABLE_HIP_MULTIPLE_KERNEL_INSTANTIATIONS)
mark_as_advanced(FORCE var Kokkos_ENABLE_HIP_RELOCATABLE_DEVICE_CODE)
mark_as_advanced(FORCE var Kokkos_ENABLE_HPX)
mark_as_advanced(FORCE var Kokkos_ENABLE_HWLOC)
mark_as_advanced(FORCE var Kokkos_ENABLE_IMPL_CUDA_MALLOC_ASYNC)
mark_as_advanced(FORCE var Kokkos_ENABLE_IMPL_HPX_ASYNC_DISPATCH)
mark_as_advanced(FORCE var Kokkos_ENABLE_LARGE_MEM_TESTS)
mark_as_advanced(FORCE var Kokkos_ENABLE_DEBUG_DUALVIEW_MODIFY_CHECK)
mark_as_advanced(FORCE var Kokkos_ENABLE_LIBQUADMATH)
mark_as_advanced(FORCE var Kokkos_ENABLE_LIBRT)
mark_as_advanced(FORCE var Kokkos_ENABLE_MEMKIND)
mark_as_advanced(FORCE var Kokkos_ENABLE_ONEDPL)
mark_as_advanced(FORCE var Kokkos_ENABLE_OPENACC)
mark_as_advanced(FORCE var Kokkos_ENABLE_OPENMPTARGET)
mark_as_advanced(FORCE var Kokkos_ENABLE_ROCM)
mark_as_advanced(FORCE var Kokkos_ENABLE_SYCL)
mark_as_advanced(FORCE var Kokkos_ENABLE_TESTS)
mark_as_advanced(FORCE var Kokkos_ENABLE_THREADS)
mark_as_advanced(FORCE var Kokkos_ENABLE_TUNING)
mark_as_advanced(FORCE var Kokkos_ENABLE_UNSUPPORTED_ARCHS)
mark_as_advanced(FORCE var Kokkos_HPX_DIR)
mark_as_advanced(FORCE var Kokkos_HWLOC_DIR)
mark_as_advanced(FORCE var Kokkos_MEMKIND_DIR)
mark_as_advanced(FORCE var Kokkos_ROCM_DIR)
mark_as_advanced(FORCE var Kokkos_THREADS_DIR)
endmacro()

View File

@ -31,7 +31,7 @@ target_include_directories(${target_name}
message(STATUS "\nCreating make file for executable ${target_name}") message(STATUS "\nCreating make file for executable ${target_name}")
message(STATUS " ${target_name} link libraries are: ${${target_link_libs}}") message(STATUS " ${target_name} link libraries are: ${${target_link_libs}}")
message(STATUS " ${target_name} source files are: ${source_files}") message(STATUS " ${target_name} source files are: ${${source_files}}")
message(STATUS " ${target_name} include dirs are: ${includeDirs}\n") message(STATUS " ${target_name} include dirs are: ${includeDirs}\n")

View File

@ -42,7 +42,7 @@ target_include_directories(${target_name}
message(STATUS "\nCreating make file for library ${target_name}") message(STATUS "\nCreating make file for library ${target_name}")
message(STATUS " ${target_name} link libraries are: ${${target_link_libs}}") message(STATUS " ${target_name} link libraries are: ${${target_link_libs}}")
message(STATUS " ${target_name} source files are: ${source_files}") message(STATUS " ${target_name} source files are: ${source_files}")
message(STATUS " ${target_name} include dirs are: ${includeDirs}\n") message(STATUS " ${target_name} include dirs are: ${includeDirs}\n \n")
install(TARGETS ${target_name} DESTINATION lib) install(TARGETS ${target_name} DESTINATION lib)
install(FILES ${includeFiles} DESTINATION include/${target_name}) install(FILES ${includeFiles} DESTINATION include/${target_name})

View File

@ -5,5 +5,6 @@
#cmakedefine pFlow_Build_Serial #cmakedefine pFlow_Build_Serial
#cmakedefine pFlow_Build_OpenMP #cmakedefine pFlow_Build_OpenMP
#cmakedefine pFlow_Build_Cuda #cmakedefine pFlow_Build_Cuda
#cmakedefine USE_STD_PARALLEL_ALG #cmakedefine pFlow_STD_Parallel_Alg
#cmakedefine pFlow_Build_Double #cmakedefine pFlow_Build_Double
#cmakedefine pFlow_Build_MPI

View File

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

View File

@ -29,17 +29,16 @@ Licence:
* folder. * folder.
*/ */
#include "vocabs.hpp"
#include "systemControl.hpp" #include "systemControl.hpp"
#include "geometryMotion.hpp" #include "geometry.hpp"
#include "commandLine.hpp" #include "commandLine.hpp"
#include "readControlDict.hpp" //#include "readControlDict.hpp"
using pFlow::commandLine;
int main( int argc, char* argv[] ) int main( int argc, char* argv[] )
{ {
commandLine cmds( pFlow::commandLine cmds(
"iterateGeometry", "iterateGeometry",
"Performs simulation without particles, only geometry is solved"); "Performs simulation without particles, only geometry is solved");
@ -54,6 +53,8 @@ commandLine cmds(
// this should be palced in each main // this should be palced in each main
pFlow::processors::initProcessors(argc, argv);
pFlow::initialize_pFlowProcessors();
#include "initialize_Control.hpp" #include "initialize_Control.hpp"
#include "setProperty.hpp" #include "setProperty.hpp"
@ -68,6 +69,7 @@ commandLine cmds(
// this should be palced in each main // this should be palced in each main
#include "finalize.hpp" #include "finalize.hpp"
pFlow::processors::finalizeProcessors();
} }

View File

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

View File

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

View File

@ -0,0 +1,89 @@
/*------------------------------- phasicFlow ---------------------------------
O C enter of
O O E ngineering and
O O M ultiscale modeling of
OOOOOOO F luid flow
------------------------------------------------------------------------------
Copyright (C): www.cemf.ir
email: hamid.r.norouzi AT gmail.com
------------------------------------------------------------------------------
Licence:
This file is part of phasicFlow code. It is a free software for simulating
granular and multiphase flows. You can redistribute it and/or modify it under
the terms of GNU General Public License v3 or any other later versions.
phasicFlow is distributed to help others in their research in the field of
granular and multiphase flows, but WITHOUT ANY WARRANTY; without even the
implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
-----------------------------------------------------------------------------*/
/**
* \file sphereGranFlow.cpp
* \brief sphereGranFlow solver
*
* This solver simulate granular flow of cohesion-less, spherical particles.
* Particle insertion can be activated in the solver to gradually insert
* particles into the simulation. Geometry can be defined generally with
* built-in features of the code or through ASCII stl files or a combination
* of both. For more information refer to [tutorials/sphereGranFlow/]
* (https://github.com/PhasicFlow/phasicFlow/tree/main/tutorials/sphereGranFlow) folder.
*/
#include "vocabs.hpp"
#include "property.hpp"
#include "sphereParticles.hpp"
#include "systemControl.hpp"
#include "commandLine.hpp"
/**
* DEM solver for simulating granular flow of cohesion-less particles.
*
* In the root case directory just simply enter the following command to
* run the simulation. For command line options use flag -h.
*/
int main( int argc, char* argv[])
{
pFlow::commandLine cmds(
"sphereGranFlow",
"DEM solver for non-cohesive spherical particles with particle insertion "
"mechanism and moving geometry");
bool isCoupling = false;
if(!cmds.parse(argc, argv)) return 0;
// this should be palced in each main
pFlow::processors::initProcessors(argc, argv);
pFlow::initialize_pFlowProcessors();
#include "initialize_Control.hpp"
#include "setProperty.hpp"
#include "createDEMComponents.hpp"
REPORT(0)<<"\nStart of time loop . . .\n"<<END_REPORT;
do
{
sphParticles.beforeIteration();
sphParticles.iterate();
sphParticles.afterIteration();
}while(Control++);
REPORT(0)<<"\nEnd of time loop.\n"<<END_REPORT;
// this should be palced in each main
#include "finalize.hpp"
pFlow::processors::finalizeProcessors();
}

View File

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

View File

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

View File

@ -7,9 +7,9 @@ add_subdirectory(Property)
add_subdirectory(Particles) add_subdirectory(Particles)
add_subdirectory(Geometry)
add_subdirectory(Interaction) add_subdirectory(Interaction)
add_subdirectory(MotionModel) add_subdirectory(MotionModel)
add_subdirectory(Geometry)

View File

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

View File

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

View File

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

View File

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

View File

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

View File

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

View File

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

View File

@ -0,0 +1,63 @@
#ifndef __AB2Kernels_hpp__
#define __AB2Kernels_hpp__
#include "KokkosTypes.hpp"
#include "types.hpp"
namespace pFlow::AB2Kernels
{
inline
bool intAllActive(
const word& name,
real dt,
rangeU32 activeRng,
const deviceViewType1D<realx3>& y,
const deviceViewType1D<realx3>& dy,
const deviceViewType1D<realx3>& dy1
)
{
Kokkos::parallel_for(
name,
deviceRPolicyStatic (activeRng.start(), activeRng.end()),
LAMBDA_HD(uint32 i){
y[i] += dt*(static_cast<real>(1.5) * dy[i] - static_cast<real>(0.5) * dy1[i]);
dy1[i] = dy[i];
});
Kokkos::fence();
return true;
}
inline
bool intScattered
(
const word& name,
real dt,
const pFlagTypeDevice& activePoints,
const deviceViewType1D<realx3>& y,
const deviceViewType1D<realx3>& dy,
const deviceViewType1D<realx3>& dy1
)
{
Kokkos::parallel_for(
name,
deviceRPolicyStatic (activePoints.activeRange().start(), activePoints.activeRange().end()),
LAMBDA_HD(uint32 i){
if( activePoints(i))
{
y[i] += dt*(static_cast<real>(1.5) * dy[i] - static_cast<real>(0.5) * dy1[i]);
dy1[i] = dy[i];
}
});
Kokkos::fence();
return true;
}
}
#endif

View File

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

View File

@ -17,13 +17,13 @@ Licence:
implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
-----------------------------------------------------------------------------*/ -----------------------------------------------------------------------------*/
#ifndef __AdamsBashforth2_hpp__ #ifndef __AdamsBashforth2_hpp__
#define __AdamsBashforth2_hpp__ #define __AdamsBashforth2_hpp__
#include "integration.hpp" #include "integration.hpp"
#include "pointFields.hpp" #include "pointFields.hpp"
#include "boundaryIntegrationList.hpp"
namespace pFlow namespace pFlow
{ {
@ -36,41 +36,41 @@ namespace pFlow
*/ */
class AdamsBashforth2 class AdamsBashforth2
: :
public integration public integration,
public realx3PointField_D
{ {
protected: private:
/// dy at t-dt boundaryIntegrationList boundaryList_;
realx3PointField_D& dy1_;
/// Range policy for integration kernel (alias) friend class processorAB2BoundaryIntegration;
using rpIntegration = Kokkos::RangePolicy<
DefaultExecutionSpace, const auto& dy1()const
Kokkos::Schedule<Kokkos::Static>, {
Kokkos::IndexType<int32> return static_cast<const realx3PointField_D&>(*this);
>; }
auto& dy1()
{
return static_cast<realx3PointField_D&>(*this);
}
public: public:
/// Type info /// Class info
TypeInfo("AdamsBashforth2"); ClassInfo("AdamsBashforth2");
// - Constructors // - Constructors
/// Construct from components /// Construct from components
AdamsBashforth2( AdamsBashforth2(
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<AdamsBashforth2>(*this);
}
/// Destructor /// Destructor
virtual ~AdamsBashforth2()=default; ~AdamsBashforth2()final = default;
/// Add this to the virtual constructor table /// Add this to the virtual constructor table
add_vCtor( add_vCtor(
@ -81,70 +81,55 @@ public:
// - Methods // - Methods
void updateBoundariesSlaveToMasterIfRequested()override;
/// return integration method
word method()const override
{
return "AdamsBashforth2";
}
bool predict( bool predict(
real UNUSED(dt), real UNUSED(dt),
realx3Vector_D& UNUSED(y), realx3PointField_D& UNUSED(y),
realx3Vector_D& UNUSED(dy)) override; realx3PointField_D& UNUSED(dy)) final;
bool predict(
real dt,
realx3Field_D& y,
realx3PointField_D& dy) final;
bool correct( bool correct(
real dt, real dt,
realx3Vector_D& y, realx3PointField_D& y,
realx3Vector_D& dy) override; realx3PointField_D& dy) final;
bool correctPStruct(
real dt,
pointStructure& pStruct,
realx3PointField_D& vel) final;
/*bool hearChanges
(
real t,
real dt,
uint32 iter,
const message& msg,
const anyList& varList
) override;*/
bool setInitialVals( bool setInitialVals(
const int32IndexContainer& newIndices, const int32IndexContainer& newIndices,
const realx3Vector& y) override; const realx3Vector& y) final;
bool needSetInitialVals()const override bool needSetInitialVals()const final
{ {
return false; return false;
} }
/// Integrate on all points in the active range
bool intAll(
real dt,
realx3Vector_D& y,
realx3Vector_D& dy,
range activeRng);
/// Integrate on active points in the active range
template<typename activeFunctor>
bool intRange(
real dt,
realx3Vector_D& y,
realx3Vector_D& dy,
activeFunctor activeP );
}; };
template<typename activeFunctor>
bool pFlow::AdamsBashforth2::intRange(
real dt,
realx3Vector_D& y,
realx3Vector_D& dy,
activeFunctor activeP )
{
auto d_dy = dy.deviceVectorAll();
auto d_y = y.deviceVectorAll();
auto d_dy1= dy1_.deviceVectorAll();
auto activeRng = activeP.activeRange();
Kokkos::parallel_for(
"AdamsBashforth2::correct",
rpIntegration (activeRng.first, activeRng.second),
LAMBDA_HD(int32 i){
if( activeP(i))
{
d_y[i] += dt*(static_cast<real>(3.0 / 2.0) * d_dy[i] - static_cast<real>(1.0 / 2.0) * d_dy1[i]);
d_dy1[i] = d_dy[i];
}
});
Kokkos::fence();
return true;
}
} // pFlow } // pFlow

View File

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

View File

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

View File

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

View File

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

View File

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

View File

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

View File

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

View File

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

View File

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

View File

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

View File

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

View File

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

View File

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

View File

@ -0,0 +1,55 @@
#include "boundaryIntegration.hpp"
#include "pointStructure.hpp"
pFlow::boundaryIntegration::boundaryIntegration(
const boundaryBase &boundary,
const pointStructure &pStruct,
const word &method,
integration& intgrtn
)
:
generalBoundary(boundary, pStruct, "", method),
integration_(intgrtn)
{}
pFlow::uniquePtr<pFlow::boundaryIntegration> pFlow::boundaryIntegration::create(
const boundaryBase &boundary,
const pointStructure &pStruct,
const word &method,
integration& intgrtn
)
{
word bType = angleBracketsNames2(
"boundaryIntegration",
boundary.type(),
method);
word altBType{"boundaryIntegration<none>"};
if( boundaryBasevCtorSelector_.search(bType) )
{
pOutput.space(4)<<"Creating boundary "<< Green_Text(bType)<<
" for "<<boundary.name()<<endl;
return boundaryBasevCtorSelector_[bType](boundary, pStruct, method, intgrtn);
}
else if(boundaryBasevCtorSelector_.search(altBType))
{
pOutput.space(4)<<"Creating boundary "<< Green_Text(altBType)<<
" for "<<boundary.name()<<endl;
return boundaryBasevCtorSelector_[altBType](boundary, pStruct, method, intgrtn);
}
else
{
printKeys(
fatalError << "Ctor Selector "<< bType<<
" and "<< altBType << " do not exist. \n"
<<"Avaiable ones are: \n",
boundaryBasevCtorSelector_
);
fatalExit;
}
return nullptr;
}

View File

@ -0,0 +1,91 @@
#ifndef __boundaryIntegration_hpp__
#define __boundaryIntegration_hpp__
#include "generalBoundary.hpp"
#include "virtualConstructor.hpp"
#include "pointFields.hpp"
namespace pFlow
{
class integration;
class boundaryIntegration
:
public generalBoundary
{
private:
const integration& integration_;
public:
TypeInfo("boundaryIntegration<none>");
boundaryIntegration(
const boundaryBase& boundary,
const pointStructure& pStruct,
const word& method,
integration& intgrtn);
~boundaryIntegration()override = default;
create_vCtor(
boundaryIntegration,
boundaryBase,
(
const boundaryBase& boundary,
const pointStructure& pStruct,
const word& method,
integration& intgrtn
),
(boundary, pStruct, method, intgrtn)
);
add_vCtor(
boundaryIntegration,
boundaryIntegration,
boundaryBase
);
bool hearChanges(
real t,
real dt,
uint32 iter,
const message &msg,
const anyList &varList) override
{
return true;
}
const integration& Integration()const
{
return integration_;
}
virtual
bool correct(real dt, const realx3PointField_D& y, const realx3PointField_D& dy)
{
return true;
}
virtual
bool correctPStruct(real dt, const realx3PointField_D& vel)
{
return true;
}
static
uniquePtr<boundaryIntegration> create(
const boundaryBase& boundary,
const pointStructure& pStruct,
const word& method,
integration& intgrtn);
};
}
#endif

View File

@ -0,0 +1,59 @@
#include "boundaryIntegrationList.hpp"
#include "integration.hpp"
pFlow::boundaryIntegrationList::boundaryIntegrationList(
const pointStructure &pStruct,
const word &method,
integration &intgrtn
)
:
ListPtr<boundaryIntegration>(6),
boundaries_(pStruct.boundaries())
{
for(uint32 i=0; i<6; i++)
{
this->set(
i,
boundaryIntegration::create(
boundaries_[i],
pStruct,
method,
intgrtn
)
);
}
}
bool pFlow::boundaryIntegrationList::correct(real dt, realx3PointField_D &y, realx3PointField_D &dy)
{
for(auto& bndry:*this)
{
if(!bndry->correct(dt, y, dy))
{
fatalErrorInFunction<<"Error in correcting boundary "<<
bndry->boundaryName()<<endl;
return false;
}
}
return true;
}
bool pFlow::boundaryIntegrationList::correctPStruct(
real dt,
pointStructure &pStruct,
const realx3PointField_D &vel)
{
for(auto& bndry:*this)
{
if(!bndry->correctPStruct(dt, vel))
{
fatalErrorInFunction<<"Error in correcting boundary "<<
bndry->boundaryName()<<" in pointStructure."<<endl;
return false;
}
}
return true;
}

View File

@ -0,0 +1,51 @@
#ifndef __boundaryIntegrationList_hpp__
#define __boundaryIntegrationList_hpp__
#include "boundaryList.hpp"
#include "ListPtr.hpp"
#include "boundaryIntegration.hpp"
namespace pFlow
{
class integration;
class boundaryIntegrationList
:
public ListPtr<boundaryIntegration>
{
private:
const boundaryList& boundaries_;
public:
boundaryIntegrationList(
const pointStructure& pStruct,
const word& method,
integration& intgrtn
);
~boundaryIntegrationList()=default;
bool correct(
real dt,
realx3PointField_D& y,
realx3PointField_D& dy);
bool correctPStruct(
real dt,
pointStructure& pStruct,
const realx3PointField_D& vel);
};
}
#endif //__boundaryIntegrationList_hpp__

View File

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

View File

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

View File

@ -1,11 +1,37 @@
set(SourceFiles set(SourceFiles
contactSearch/methods/cellBased/NBS/mapperNBS.cpp
contactSearch/methods/cellBased/NBS/mapperNBSKernels.cpp
contactSearch/methods/cellBased/NBS/NBSLevel0.cpp
contactSearch/methods/cellBased/NBS/NBS.cpp
contactSearch/methods/cellBased/NBS/cellsWallLevel0.cpp
contactSearch/boundaries/boundaryContactSearch/boundaryContactSearch.cpp
contactSearch/boundaries/periodicBoundaryContactSearch/ppwBndryContactSearchKernels.cpp
contactSearch/boundaries/periodicBoundaryContactSearch/ppwBndryContactSearch.cpp
contactSearch/boundaries/periodicBoundaryContactSearch/wallBoundaryContactSearch.cpp
contactSearch/boundaries/periodicBoundaryContactSearch/periodicBoundaryContactSearch.cpp
contactSearch/boundaries/boundaryContactSearchList.cpp
contactSearch/contactSearch/contactSearch.cpp contactSearch/contactSearch/contactSearch.cpp
contactSearch/ContactSearch/ContactSearchs.cpp contactSearch/ContactSearch/ContactSearchs.cpp
interaction/interaction.cpp interaction/interaction.cpp
sphereInteraction/sphereInteractions.cpp sphereInteraction/sphereInteractionsLinearModels.cpp
sphereInteraction/sphereInteractionsNonLinearModels.cpp
sphereInteraction/sphereInteractionsNonLinearModModels.cpp
) )
if(pFlow_Build_MPI)
list(APPEND SourceFiles
contactSearch/boundaries/processorBoundaryContactSearch/processorBoundaryContactSearch.cpp
sphereInteraction/boundaries/processorBoundarySphereInteraction/processorBoundarySphereInteractions.cpp
contactSearch/boundaries/twoPartContactSearch/twoPartContactSearchKernels.cpp
contactSearch/boundaries/twoPartContactSearch/twoPartContactSearch.cpp
)
endif()
set(link_libs Kokkos::kokkos phasicFlow Property Particles Geometry) set(link_libs Kokkos::kokkos phasicFlow Property Particles Geometry)
pFlow_add_library_install(Interaction SourceFiles link_libs) pFlow_add_library_install(Interaction SourceFiles link_libs)

View File

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

View File

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

View File

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

View File

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

View File

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

View File

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

View File

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

View File

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

View File

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

View File

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

View File

@ -0,0 +1,77 @@
/*------------------------------- phasicFlow ---------------------------------
O C enter of
O O E ngineering and
O O M ultiscale modeling of
OOOOOOO F luid flow
------------------------------------------------------------------------------
Copyright (C): www.cemf.ir
email: hamid.r.norouzi AT gmail.com
------------------------------------------------------------------------------
Licence:
This file is part of phasicFlow code. It is a free software for simulating
granular and multiphase flows. You can redistribute it and/or modify it under
the terms of GNU General Public License v3 or any other later versions.
phasicFlow is distributed to help others in their research in the field of
granular and multiphase flows, but WITHOUT ANY WARRANTY; without even the
implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
-----------------------------------------------------------------------------*/
#include "boundaryContactSearch.hpp"
#include "contactSearch.hpp"
pFlow::boundaryContactSearch::boundaryContactSearch(
const dictionary &dict,
const boundaryBase &boundary,
const contactSearch &cSearch)
: generalBoundary(
boundary,
cSearch.pStruct(),
"",
""),
contactSearch_(cSearch),
updateInterval_(dict.getVal<uint32>("updateInterval"))
{
}
pFlow::uniquePtr<pFlow::boundaryContactSearch>
pFlow::boundaryContactSearch::create(
const dictionary &dict,
const boundaryBase &boundary,
const contactSearch &cSearch)
{
word bType = angleBracketsNames2(
"boundaryContactSearch",
pFlowProcessors().localRunTypeName(),
boundary.type());
word altBType{"boundaryContactSearch<none>"};
if( boundaryBasevCtorSelector_.search(bType) )
{
pOutput.space(4)<<"Creating contact search boundary "<< Green_Text(bType)<<
" for "<<boundary.name()<<endl;
return boundaryBasevCtorSelector_[bType](dict, boundary, cSearch);
}
else if(boundaryBasevCtorSelector_.search(altBType))
{
pOutput.space(4)<<"Creating contact search boundary "<< Green_Text(altBType)<<
" for "<<boundary.name()<<endl;
return boundaryBasevCtorSelector_[altBType](dict, boundary, cSearch);
}
else
{
printKeys
(
fatalError << "Ctor Selector "<< bType<<
" and "<< altBType << " do not exist. \n"
<<"Avaiable ones are: \n"
,
boundaryBasevCtorSelector_
);
fatalExit;
}
return nullptr;
}

View File

@ -0,0 +1,109 @@
/*------------------------------- phasicFlow ---------------------------------
O C enter of
O O E ngineering and
O O M ultiscale modeling of
OOOOOOO F luid flow
------------------------------------------------------------------------------
Copyright (C): www.cemf.ir
email: hamid.r.norouzi AT gmail.com
------------------------------------------------------------------------------
Licence:
This file is part of phasicFlow code. It is a free software for simulating
granular and multiphase flows. You can redistribute it and/or modify it under
the terms of GNU General Public License v3 or any other later versions.
phasicFlow is distributed to help others in their research in the field of
granular and multiphase flows, but WITHOUT ANY WARRANTY; without even the
implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
-----------------------------------------------------------------------------*/
#ifndef __boundaryContactSearch_hpp__
#define __boundaryContactSearch_hpp__
#include "generalBoundary.hpp"
#include "contactSearchGlobals.hpp"
#include "virtualConstructor.hpp"
namespace pFlow
{
class contactSearch;
class boundaryContactSearch
: public generalBoundary
{
private:
const contactSearch &contactSearch_;
/// @brief update interval in terms of iteration numebr
uint32 updateInterval_;
/// @brief last iteration number which contact search has been performed
uint32 lastUpdated_ = 0;
/// @brief performed search?
bool performedSearch_ = false;
public:
// type info
TypeInfo("boundaryContactSearch<none>");
boundaryContactSearch(
const dictionary &dict,
const boundaryBase &boundary,
const contactSearch &cSearch);
create_vCtor(
boundaryContactSearch,
boundaryBase,
(
const dictionary &dict,
const boundaryBase &boundary,
const contactSearch &cSearch),
(dict, boundary, cSearch));
add_vCtor(
boundaryContactSearch,
boundaryContactSearch,
boundaryBase);
const contactSearch &cSearch() const
{
return contactSearch_;
}
bool hearChanges(
real t,
real dt,
uint32 iter,
const message &msg,
const anyList &varList) override
{
if (msg.equivalentTo(message::BNDR_RESET))
{
// do nothing
}
return true;
}
virtual bool broadSearch(
uint32 iter,
real t,
real dt,
csPairContainerType &ppPairs,
csPairContainerType &pwPairs,
bool force = false)
{
return true;
}
static uniquePtr<boundaryContactSearch> create(
const dictionary &dict,
const boundaryBase &boundary,
const contactSearch &cSearch);
};
}
#endif //__boundaryContactSearch_hpp__

View File

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

View File

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

View File

@ -0,0 +1,131 @@
/*------------------------------- phasicFlow ---------------------------------
O C enter of
O O E ngineering and
O O M ultiscale modeling of
OOOOOOO F luid flow
------------------------------------------------------------------------------
Copyright (C): www.cemf.ir
email: hamid.r.norouzi AT gmail.com
------------------------------------------------------------------------------
Licence:
This file is part of phasicFlow code. It is a free software for simulating
granular and multiphase flows. You can redistribute it and/or modify it under
the terms of GNU General Public License v3 or any other later versions.
phasicFlow is distributed to help others in their research in the field of
granular and multiphase flows, but WITHOUT ANY WARRANTY; without even the
implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
-----------------------------------------------------------------------------*/
#include "periodicBoundaryContactSearch.hpp"
#include "contactSearch.hpp"
#include "particles.hpp"
#include "pointStructure.hpp"
#include "geometry.hpp"
void pFlow::periodicBoundaryContactSearch::setSearchBox()
{
auto db = pStruct().thisDomain().domainBox();
auto n1 = boundary().mirrorBoundary().boundaryPlane().normal();
auto l1 = boundary().mirrorBoundary().neighborLength();
auto n2 = boundary().boundaryPlane().normal();
auto l2 = boundary().neighborLength();
realx3 minP = db.minPoint() + (db.maxPoint()-db.minPoint())* n1+(n2*l2);
realx3 maxP = db.maxPoint() + (n1*l1);
searchBox_={minP, maxP};
}
pFlow::periodicBoundaryContactSearch::periodicBoundaryContactSearch(
const dictionary &dict,
const boundaryBase &boundary,
const contactSearch &cSearch)
:
boundaryContactSearch(dict, boundary, cSearch),
transferVec_(boundary.mirrorBoundary().displacementVectroToMirror()),
thisIndex_(thisBoundaryIndex()),
mirrorIndex_(mirrorBoundaryindex()),
diameter_(cSearch.Particles().boundingSphere())
{
if(thisIndex_%2==1)
{
masterSearch_ = true;
setSearchBox();
real minD;
real maxD;
cSearch.Particles().boundingSphereMinMax(minD, maxD);
ppContactSearch_ = makeUnique<ppwBndryContactSearch>(
searchBox_,
maxD);
const auto& geom = cSearch.Geometry();
pwContactSearch_ = makeUnique<wallBoundaryContactSearch>(
0.5,
geom.numPoints(),
geom.size(),
geom.points().deviceViewAll(),
geom.vertices().deviceViewAll(),
geom.normals().deviceViewAll());
}
else
{
masterSearch_ = false;
searchBox_={{0,0,0},{0,0,0}};
}
}
bool pFlow::periodicBoundaryContactSearch::broadSearch
(
uint32 iter,
real t,
real dt,
csPairContainerType &ppPairs,
csPairContainerType &pwPairs,
bool force
)
{
if(masterSearch_)
{
auto thisP = boundary().thisPoints();
auto thisDiams = diameter_.BoundaryField(thisIndex_).thisField();
auto mirrorP = mirrorBoundary().thisPoints();
auto mirrorDiams = diameter_.BoundaryField(mirrorIndex_).thisField();
ppContactSearch_().broadSearchPP(
ppPairs,
thisP,
thisDiams,
mirrorP,
mirrorDiams,
transferVec_);
/*pwContactSearch_().broadSearch(
pwPairs,
ppContactSearch_().searchCells(),
thisP,
thisDiams,
mirrorP,
mirrorDiams,
transferVec_,
ppContactSearch_().sizeRatio());*/
//output<<t<<" boundary pp size "<< ppPairs.size()<<endl;
//output<<t<<" boundary pw size "<< pwPairs.size()<<endl;
return true;
}else
{
return true;
}
}

View File

@ -0,0 +1,81 @@
/*------------------------------- phasicFlow ---------------------------------
O C enter of
O O E ngineering and
O O M ultiscale modeling of
OOOOOOO F luid flow
------------------------------------------------------------------------------
Copyright (C): www.cemf.ir
email: hamid.r.norouzi AT gmail.com
------------------------------------------------------------------------------
Licence:
This file is part of phasicFlow code. It is a free software for simulating
granular and multiphase flows. You can redistribute it and/or modify it under
the terms of GNU General Public License v3 or any other later versions.
phasicFlow is distributed to help others in their research in the field of
granular and multiphase flows, but WITHOUT ANY WARRANTY; without even the
implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
-----------------------------------------------------------------------------*/
#ifndef __periodicBoundaryContactSearch_hpp__
#define __periodicBoundaryContactSearch_hpp__
#include "boundaryContactSearch.hpp"
#include "box.hpp"
#include "ppwBndryContactSearch.hpp"
#include "pointFields.hpp"
#include "wallBoundaryContactSearch.hpp"
namespace pFlow
{
class periodicBoundaryContactSearch
: public boundaryContactSearch
{
private:
box searchBox_;
realx3 transferVec_;
uint32 thisIndex_;
uint32 mirrorIndex_;
uniquePtr<ppwBndryContactSearch> ppContactSearch_ = nullptr;
uniquePtr<wallBoundaryContactSearch> pwContactSearch_ = nullptr;
const realPointField_D &diameter_;
bool masterSearch_ = false;
void setSearchBox();
public:
TypeInfo("boundaryContactSearch<regular,periodic>")
periodicBoundaryContactSearch(
const dictionary &dict,
const boundaryBase &boundary,
const contactSearch &cSearch);
~periodicBoundaryContactSearch() override = default;
add_vCtor(
boundaryContactSearch,
periodicBoundaryContactSearch,
boundaryBase);
bool broadSearch(
uint32 iter,
real t,
real dt,
csPairContainerType &ppPairs,
csPairContainerType &pwPairs,
bool force = false) override;
};
}
#endif //__periodicBoundaryContactSearch_hpp__

View File

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

View File

@ -0,0 +1,86 @@
/*------------------------------- phasicFlow ---------------------------------
O C enter of
O O E ngineering and
O O M ultiscale modeling of
OOOOOOO F luid flow
------------------------------------------------------------------------------
Copyright (C): www.cemf.ir
email: hamid.r.norouzi AT gmail.com
------------------------------------------------------------------------------
Licence:
This file is part of phasicFlow code. It is a free software for simulating
granular and multiphase flows. You can redistribute it and/or modify it under
the terms of GNU General Public License v3 or any other later versions.
phasicFlow is distributed to help others in their research in the field of
granular and multiphase flows, but WITHOUT ANY WARRANTY; without even the
implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
-----------------------------------------------------------------------------*/
#ifndef __ppwBndryContactSearch_hpp__
#define __ppwBndryContactSearch_hpp__
#include "contactSearchGlobals.hpp"
#include "scatteredFieldAccess.hpp"
#include "cells.hpp"
namespace pFlow
{
class ppwBndryContactSearch
{
public:
using HeadType = deviceViewType3D<uint32>;
using NextType = deviceViewType1D<uint32>;
private:
cells searchCells_;
HeadType head_{ "periodic::head", 1, 1, 1 };
NextType next_{ "periodic::next", 1 };
real sizeRatio_ = 1.0;
uint32 nextCapacity_ = 0;
void checkAllocateNext(uint32 n);
void nullifyHead();
void nullifyNext(uint32 n);
void buildList(
const deviceScatteredFieldAccess<realx3> &points);
public:
ppwBndryContactSearch(
const box &domain,
real cellSize,
real sizeRatio = 1.0);
bool broadSearchPP(
csPairContainerType &ppPairs,
const deviceScatteredFieldAccess<realx3> &points,
const deviceScatteredFieldAccess<real> &diams,
const deviceScatteredFieldAccess<realx3> &mirrorPoints,
const deviceScatteredFieldAccess<real> &mirrorDiams,
const realx3 &transferVec);
const auto& searchCells()const
{
return searchCells_;
}
real sizeRatio()const
{
return sizeRatio_;
}
};
}
#endif //__ppwBndryContactSearch_hpp__

View File

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

View File

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

View File

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

View File

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

View File

@ -19,54 +19,124 @@ Licence:
-----------------------------------------------------------------------------*/ -----------------------------------------------------------------------------*/
#include "contactSearch.hpp" #include "contactSearch.hpp"
#include "streams.hpp"
#include "particles.hpp"
pFlow::contactSearch::contactSearch( pFlow::contactSearch::contactSearch(
const dictionary& dict, const dictionary& dict,
const box& domain, const box& extDomain,
const particles& prtcl, const particles& prtcl,
const geometry& geom, const geometry& geom,
Timers& timers) Timers& timers)
: :
interactionBase(prtcl, geom), extendedDomainBox_(extDomain),
domain_(domain), particles_(prtcl),
dict_(dict), geometry_(geom),
sphereSphereTimer_("particle-particle contact search", &timers), bTimer_("Boundary particles contact search", &timers),
sphereWallTimer_("particle-wall contact search", &timers) ppTimer_("Internal particles contact search", &timers),
dict_(dict)
{ {
} }
const pFlow::pointStructure &pFlow::contactSearch::pStruct() const
{
return particles_.pStruct();
}
pFlow::uniquePtr<pFlow::contactSearch> pFlow::contactSearch::create( bool pFlow::contactSearch::broadSearch
const dictionary& dict, (
const box& domain, const timeInfo &ti,
const particles& prtcl, csPairContainerType &ppPairs,
const geometry& geom, csPairContainerType &pwPairs,
Timers& timers) bool force
)
{ {
if(enterBroadSearch(ti, force))
{
ppTimer_.start();
if( !BroadSearch(
ti,
ppPairs,
pwPairs,
force ) )
{
fatalErrorInFunction;
performedSearch_ = false;
return false;
}
ppTimer_.end();
performedSearch_ = true;
}
else
{
performedSearch_ = false;
}
return true;
}
bool pFlow::contactSearch::boundaryBroadSearch
(
uint32 bndryIndex,
const timeInfo &ti,
csPairContainerType &ppPairs,
csPairContainerType &pwPairs,
bool force
)
{
if(enterBroadSearchBoundary(ti, force))
{
bTimer_.start();
for(uint32 i=0u; i<6u; i++)
{
if(!BoundaryBroadSearch(
i,
ti,
ppPairs,
pwPairs,
force))
{
performedSearchBoundary_ = false;
return false;
}
}
bTimer_.end();
performedSearchBoundary_ = true;
}
else
{
performedSearchBoundary_ = false;
}
return true;
}
pFlow::uniquePtr<pFlow::contactSearch> pFlow::contactSearch::create(
const dictionary &dict,
const box &extDomain,
const particles &prtcl,
const geometry &geom,
Timers &timers)
{
word baseMethName = dict.getVal<word>("method"); word baseMethName = dict.getVal<word>("method");
word wallMethod = dict.getVal<word>("wallMapping");
auto model = angleBracketsNames2("ContactSearch", baseMethName, wallMethod);
REPORT(1)<<"Selecting contact search model . . ."<<endREPORT;
auto model = angleBracketsNames("ContactSearch", baseMethName);
pOutput.space(2)<<"Selecting contact search model "<<Green_Text(model)<<endl;
if( dictionaryvCtorSelector_.search(model)) if( dictionaryvCtorSelector_.search(model))
{ {
auto objPtr = dictionaryvCtorSelector_[model] (dict, domain, prtcl, geom, timers); auto objPtr = dictionaryvCtorSelector_[model] (dict, extDomain, prtcl, geom, timers);
REPORT(2)<<"Model "<< greenText(model)<<" is created."<<endREPORT;
return objPtr; return objPtr;
} }
else else
{ {
printKeys printKeys
( (
fatalError << "Ctor Selector "<< model << " dose not exist. \n" fatalError << "Ctor Selector "<< model << " does not exist. \n"
<<"Avaiable ones are: \n\n" <<"Avaiable ones are: \n\n"
, ,
dictionaryvCtorSelector_ dictionaryvCtorSelector_

View File

@ -23,42 +23,68 @@ Licence:
#define __contactSearch_hpp__ #define __contactSearch_hpp__
#include "interactionBase.hpp" #include "contactSearchGlobals.hpp"
#include "unsortedPairs.hpp"
#include "box.hpp"
#include "dictionary.hpp" #include "dictionary.hpp"
#include "virtualConstructor.hpp"
#include "timeInfo.hpp"
#include "Timer.hpp"
namespace pFlow namespace pFlow
{ {
// - forward
class box;
class particles;
class geometry;
class pointStructure;
class contactSearch class contactSearch
:
public interactionBase
{ {
public: private:
using IdType = typename interactionBase::IdType;
using IndexType = typename interactionBase::IndexType; const box& extendedDomainBox_;
using ExecutionSpace = typename interactionBase::ExecutionSpace; /// @brief update interval in terms of iteration numebr
uint32 updateInterval_= 1;
using PairContainerType = unsortedPairs<ExecutionSpace, IdType>; /// @brief last iteration number which contact search has been performed
uint32 lastUpdated_ = 0;
protected: /// @brief performed search?
bool performedSearch_ = false;
const box& domain_; /// @brief performed search in boundaries
bool performedSearchBoundary_ = false;
/// const ref to particles
const particles& particles_;
/// const ref to geometry
const geometry& geometry_;
Timer bTimer_;
Timer ppTimer_;
dictionary dict_; dictionary dict_;
Timer sphereSphereTimer_; virtual
bool BroadSearch(
const timeInfo& ti,
csPairContainerType& ppPairs,
csPairContainerType& pwPairs,
bool force
)=0;
Timer sphereWallTimer_; virtual
bool BoundaryBroadSearch(
auto& dict() uint32 bndryIndex,
{ const timeInfo& ti,
return dict_; csPairContainerType& ppPairs,
} csPairContainerType& pwPairs,
bool force = false
)=0;
public: public:
@ -66,7 +92,7 @@ public:
contactSearch( contactSearch(
const dictionary& dict, const dictionary& dict,
const box& domain, const box& extDomain,
const particles& prtcl, const particles& prtcl,
const geometry& geom, const geometry& geom,
Timers& timers); Timers& timers);
@ -88,40 +114,103 @@ public:
(dict, domain, prtcl, geom, timers) (dict, domain, prtcl, geom, timers)
); );
const auto& domain()const inline
bool performedSearch()const
{ {
return domain_; return performedSearch_;
} }
const auto& dict()const inline
bool performedSearchBoundary()const
{
return performedSearchBoundary_;
}
inline
bool performSearch(uint32 iter, bool force = false)const
{
if((iter-lastUpdated_) % updateInterval_ == 0 || iter == 0 || force )
{
return true;
}
return false;
}
bool enterBroadSearch(const timeInfo& ti, bool force = false)const
{
return performSearch(ti.iter(), force);
}
virtual
bool enterBroadSearchBoundary(const timeInfo& ti, bool force=false)const = 0;
inline
uint32 updateInterval()const
{
return updateInterval_;
}
inline
const dictionary& dict()const
{ {
return dict_; return dict_;
} }
inline
const box& extendedDomainBox()const
{
return extendedDomainBox_;
}
inline
const particles& Particles()const
{
return particles_;
}
const pointStructure& pStruct()const;
inline
const geometry& Geometry()const
{
return geometry_;
}
inline
Timer& ppTimer()
{
return ppTimer_;
}
inline
Timer& bTimer()
{
return bTimer_;
}
virtual
bool broadSearch( bool broadSearch(
PairContainerType& ppPairs, const timeInfo& ti,
PairContainerType& pwPairs, csPairContainerType& ppPairs,
bool force = false) = 0; csPairContainerType& pwPairs,
bool force = false);
bool boundaryBroadSearch(
uint32 bndryIndex,
const timeInfo& ti,
csPairContainerType& ppPairs,
csPairContainerType& pwPairs,
bool force = false);
virtual virtual
bool ppEnterBroadSearch()const = 0; real sizeRatio()const = 0;
virtual virtual
bool pwEnterBroadSearch()const = 0; real cellExtent()const = 0;
virtual
bool ppPerformedBroadSearch()const = 0;
virtual
bool pwPerformedBroadSearch()const = 0;
static static
uniquePtr<contactSearch> create( uniquePtr<contactSearch> create(
const dictionary& dict, const dictionary& dict,
const box& domain, const box& extDomain,
const particles& prtcl, const particles& prtcl,
const geometry& geom, const geometry& geom,
Timers& timers); Timers& timers);

View File

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

View File

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

View File

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

View File

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

View File

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

View File

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

View File

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

View File

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

View File

@ -0,0 +1,69 @@
/*------------------------------- phasicFlow ---------------------------------
O C enter of
O O E ngineering and
O O M ultiscale modeling of
OOOOOOO F luid flow
------------------------------------------------------------------------------
Copyright (C): www.cemf.ir
email: hamid.r.norouzi AT gmail.com
------------------------------------------------------------------------------
Licence:
This file is part of phasicFlow code. It is a free software for simulating
granular and multiphase flows. You can redistribute it and/or modify it under
the terms of GNU General Public License v3 or any other later versions.
phasicFlow is distributed to help others in their research in the field of
granular and multiphase flows, but WITHOUT ANY WARRANTY; without even the
implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
-----------------------------------------------------------------------------*/
#include "NBS.hpp"
pFlow::NBS::NBS
(
const dictionary& dict,
const box& domainBox,
real minBSSize,
real maxBSSize,
const deviceViewType1D<realx3> &position,
const pFlagTypeDevice &flags,
const deviceViewType1D<real> &diam,
uint32 nWallPoints,
uint32 nWallElements,
const ViewType1D<realx3,memory_space>& wallPoints,
const ViewType1D<uint32x3,memory_space>& wallVertices,
const ViewType1D<realx3,memory_space>& wallNormals
)
:
particleWallContactSearchs<NBS>(
dict,
domainBox,
minBSSize,
maxBSSize,
position,
flags,
diam),
sizeRatio_(max(dict.getVal<real>("sizeRatio"), 1.0)),
cellExtent_(max(dict.getVal<real>("cellExtent"), 0.5)),
adjustableBox_(dict.getVal<Logical>("adjustableBox")),
NBSLevel0_
(
this->domainBox_,
maxBSSize,
sizeRatio_,
position,
flags,
adjustableBox_()
),
cellsWallLevel0_
(
cellExtent_,
nWallPoints,
nWallElements,
wallPoints,
wallVertices,
wallNormals
)
{
}

View File

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

View File

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

View File

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

View File

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

View File

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

View File

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

View File

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

View File

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

View File

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

View File

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

View File

@ -0,0 +1,49 @@
/*------------------------------- phasicFlow ---------------------------------
O C enter of
O O E ngineering and
O O M ultiscale modeling of
OOOOOOO F luid flow
------------------------------------------------------------------------------
Copyright (C): www.cemf.ir
email: hamid.r.norouzi AT gmail.com
------------------------------------------------------------------------------
Licence:
This file is part of phasicFlow code. It is a free software for simulating
granular and multiphase flows. You can redistribute it and/or modify it under
the terms of GNU General Public License v3 or any other later versions.
phasicFlow is distributed to help others in their research in the field of
granular and multiphase flows, but WITHOUT ANY WARRANTY; without even the
implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
-----------------------------------------------------------------------------*/
#include "phasicFlowKokkos.hpp"
#include "cells.hpp"
#include "pointFlag.hpp"
namespace pFlow::mapperNBSKernels
{
void findPointExtends(
const deviceViewType1D<realx3>& points,
const pFlagTypeDevice& flags,
realx3& minPoint,
realx3& maxPoint);
bool buildListsReduce(
const cells& searchCell,
const deviceViewType3D<uint32>& head,
const deviceViewType1D<uint32>& next,
const deviceViewType1D<realx3>& points,
const pFlagTypeDevice& flags);
bool buildLists(
const cells& searchCell,
const deviceViewType3D<uint32>& head,
const deviceViewType1D<uint32>& next,
const deviceViewType1D<realx3>& points,
const pFlagTypeDevice& flags);
}

View File

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

View File

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

View File

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

View File

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

View File

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

View File

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

View File

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

View File

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

View File

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

View File

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

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