235 Commits

Author SHA1 Message Date
HRN
774afd5f37 bug correction for the time when empty is used 2025-02-14 22:50:28 +03:30
191801b344 Merge pull request #165 from ramin1728/main
binarySystemOfParticles is Updated.
2025-02-14 20:42:14 +03:30
545de300ae Merge pull request #166 from PhasicFlow/develop
edits
2025-02-14 20:40:44 +03:30
HRN
9b3c4f83b9 edits 2025-02-14 20:39:37 +03:30
b315d12357 conveyorBelt is Updated. 2025-02-11 23:35:58 +03:30
7e7184f1c5 binarySystemOfParticles is Updated. 2025-02-11 23:18:29 +03:30
9f489d07cc Merge pull request #163 from PhasicFlow/develop
global damping is activated for velocity and rVelocity in both sphere…
2025-02-10 01:11:30 +03:30
HRN
8466e02d81 global damping is activated for velocity and rVelocity in both sphere and grain solvers 2025-02-10 01:10:13 +03:30
3f1fa4ae90 Merge pull request #162 from wanqing0421/main
bug fix during the build process of cuda mode
2025-02-08 22:39:06 +03:30
f0f185983c bug fix during the build process of cuda mode 2025-02-08 23:47:21 +08:00
cb6d567dab Merge pull request #161 from PhasicFlow/develop
AB5 is added and bug is resolved when particles are being inserted
2025-02-08 18:07:39 +03:30
HRN
db9b1e62e4 AB5 is added and bug is resolved when particles are being inserted 2025-02-08 18:06:30 +03:30
3aff0f1fc6 Merge pull request #160 from PhasicFlow/develop
bug resolve, chekcForCollision is set to true for always, adjustable …
2025-02-07 23:16:48 +03:30
HRN
b9ab015eb1 bug resolve, chekcForCollision is set to true for always, adjustable search box is set to false for always, old hearChanges has been removed 2025-02-07 23:12:53 +03:30
59fbaa91ab Merge pull request #159 from PhasicFlow/develop
All messages are revisited for internal points and boundaries
2025-02-06 11:06:35 +03:30
f84484881c Merge pull request #158 from PhasicFlow/develop
pFlowToVTK now manages when Ctrl+C is used by user
2025-02-06 11:05:21 +03:30
HRN
02e0b72082 All messages are revisited for internal points and boundaries 2025-02-06 11:04:30 +03:30
HRN
edb02ecfc7 pFlowToVTK now manages when Ctrl+C is used by user 2025-02-06 10:51:13 +03:30
1540321a31 Merge pull request #157 from PhasicFlow/develop
the need to provide neighborLength in domain dictionary is lifted. No…
2025-02-03 23:51:31 +03:30
HRN
63bd9c9993 the need to provide neighborLength in domain dictionary is lifted. Now it is optional 2025-02-03 23:49:11 +03:30
f4f5f29e3c Merge pull request #156 from PhasicFlow/develop
Develop
2025-02-03 19:17:08 +03:30
HRN
fac5576df1 periodict boundary condition ->DEBUG and some minor changes boundaries 2025-02-03 19:16:14 +03:30
HRN
f5ba30b901 autoCompelte for time folders and field names 2025-02-03 19:15:08 +03:30
HRN
0f9a19d6ee reduction in boundary class size 2025-02-01 23:33:19 +03:30
d909301f32 Merge pull request #155 from PhasicFlow/develop
Develop
2025-02-01 22:18:23 +03:30
HRN
3b88b6156d Merge branch 'develop' of github.com:PhasicFlow/phasicFlow into develop 2025-02-01 22:15:34 +03:30
HRN
64c041a753 boundaryListPtr is created and other classes were changed accordingly 2025-02-01 22:14:41 +03:30
0410eb6864 Merge pull request #154 from PhasicFlow/main
merge from main
2025-01-31 13:45:24 +03:30
2d7f7ec17f Merge pull request #153 from PhasicFlow/develop
changes in contactSearch and timeControl and timeInfo
2025-01-31 13:44:07 +03:30
HRN
af2572331d some cleanup 2025-01-31 01:06:16 +03:30
HRN
f98a696fc8 changes in contactSearch and timeControl and timeInfo 2025-01-27 14:26:20 +03:30
2168456584 correctoins for running on cuda 2025-01-26 12:19:25 +03:30
6e6cabbefa Merge pull request #152 from PhasicFlow/develop
timeInfo is updated
2025-01-25 19:58:27 +03:30
HRN
42b024e1ed globalDamping is deactivated for future developement 2025-01-25 19:56:01 +03:30
HRN
debb8fd037 bug fix for stridedRange 2025-01-25 19:48:36 +03:30
HRN
0fc9eea561 timeInfo and timeValue
- timeInfo updated
- timeValue type is added = double
- AB2 bug fixed
- updateInterval for broadSearch is activated
2025-01-25 19:18:11 +03:30
HRN
1ccc321c1d Merge branch 'develop' of github.com:PhasicFlow/phasicFlow into develop 2025-01-24 21:15:51 +03:30
HRN
f4b15bc50a timeControl 2025-01-24 21:15:16 +03:30
164eedb27c Merge pull request #151 from PhasicFlow/main
update from main branch
2025-01-24 21:12:53 +03:30
2ec3dfdc7e global damping is added to the code 2025-01-20 21:02:50 +03:30
cb1faf04f8 bug fixes for build with float in cuda 2025-01-20 15:43:56 +03:30
a32786eb8a particlePhasicFlow bug-fix when flag --set-fields-only is used and no shpaeName is set 2025-01-20 15:39:17 +03:30
967bb753aa adding non-linear contact force model for grains 2025-01-20 15:37:48 +03:30
c202f9eaae AB3, AB4 added, and AB2 modified 2025-01-20 14:55:12 +03:30
HRN
8823dd3c94 file headers 2025-01-10 12:27:35 +03:30
3d6fa28221 Merge pull request #150 from PhasicFlow/develop
Auto-complete for command line
2025-01-09 21:03:00 +03:30
HRN
2155771044 autoComplete file 2025-01-09 21:02:01 +03:30
HRN
bfa3759fc9 Auto complete for command line is added to phasicFlow commands 2025-01-09 20:58:33 +03:30
bab1984b8e Merge pull request #149 from PhasicFlow/develop
Binary conversion from pointFiled files to vtk format
2025-01-09 19:20:25 +03:30
HRN
66c3cda39e Binary conversion from pointFiled files to vtk format
- all vtk files are stored in the /particles/ folder
- terminal output is modified
- time series file are added. files with extntions .vtk.series can be loaded into paraview
2025-01-09 19:16:18 +03:30
75975a43de Merge pull request #148 from Alireza77h/main
bug fixes for cuda run in plus side
2025-01-08 11:39:03 +03:30
f74ec1d1ac bug fixes for cuda run in plus side 2025-01-08 08:42:45 +03:30
f169e3fc89 Merge pull request #147 from PhasicFlow/develop
bug fix for memory error in PhasicFlowPlus-fluid interaction
2025-01-04 15:49:27 +03:30
HRN
d39e1ec27b bug fix for memory error in PhasicFlowPlus-fluid interaction 2025-01-04 15:46:30 +03:30
dbb3b5c3a2 Merge pull request #143 from ramin1728/V-blender
updated V1.0 V-blender
2024-12-28 17:45:11 +03:30
c02027f5d2 updated V1.0 V-blender 2024-12-28 14:58:28 +03:30
31ebc5963c updated V1.0 toteblender 2024-12-28 14:09:09 +03:30
0d5f8a98a2 updated V1.0 screwConveyor 2024-12-28 13:46:09 +03:30
018770a1f1 updated V1.0 RotaryDrumWithBaffle 2024-12-28 13:31:17 +03:30
59145425b1 Merge pull request #139 from PhasicFlow/develop
grainDEMSystem folder is added
2024-12-28 13:25:06 +03:30
HRN
ff6c3454c9 grainDEMSystem folder is added 2024-12-28 13:22:12 +03:30
a63ed8b0fa Merge pull request #138 from ramin1728/rotaryairlock
RotaryAirLockValve
2024-12-28 13:20:17 +03:30
278afc8ab6 updated V1.0 RotaryAirLock 2024-12-28 13:09:06 +03:30
6179f6ea3b t 2024-12-28 13:00:42 +03:30
e4eb7d74be Merge pull request #137 from PhasicFlow/develop
Some changes to be compatible with phasicFlowPlus-v1.0
2024-12-27 19:29:31 +03:30
HRN
1f242d8911 Some changes to be compatible with phasicFlowPlus-v1.0 2024-12-27 19:26:09 +03:30
46e862c5e5 Merge pull request #136 from PhasicFlow/develop
sphereDEMSystem is updated for version 1.0
2024-12-26 22:59:03 +03:30
HRN
16f8ab4572 sphereDEMSystem is updated for version 1.0 2024-12-26 19:07:39 +03:30
809d4ca617 Merge pull request #135 from PhasicFlow/develop
DEMSystem updated for version 1.0
2024-12-25 18:28:05 +03:30
HRN
bc22012ecd DEMSystem updated for version 1.0 2024-12-25 18:26:53 +03:30
51d85a0b79 Merge pull request #133 from wanqing0421/main
Add the conveyor belt model and tutorial
2024-12-25 09:39:07 +03:30
bf6ecfb3ff Merge pull request #134 from ramin1728/rotatingDrumMedium
updated V1.0 rotatingDrumMedium
2024-12-25 09:34:07 +03:30
73e4295a25 Undo changes in bashrc and change the velocity_ to tangentVelocity_ 2024-12-25 08:38:09 +08:00
a99278447e update conveyor belt model tutorial 2024-12-24 22:52:18 +08:00
dae2af5354 updated V1.0 rotatingDrumMedium 2024-12-24 18:14:19 +03:30
12b3413306 Merge pull request #131 from ramin1728/rotatingDrumSmall
updated V1.0 rotatingDrumSmall
2024-12-24 18:08:32 +03:30
ebff41619e updated V1.0 rotatingDrumSmall 2024-12-24 18:03:38 +03:30
50750c2574 update conveyor belt model 2024-12-24 17:11:49 +08:00
89a47d1a15 update conveyor belt model 2024-12-24 16:59:42 +08:00
ef2a13dd7e add conveyor belt model 2024-12-24 16:39:40 +08:00
5318610c1f add converyor belt model 2024-12-24 15:37:19 +08:00
acdad47823 Merge pull request #126 from ramin1728/main
binarySystemOfParticles updated V1.0
2024-12-03 16:26:04 +03:30
1008ea8c9a layeredSiloFilling is updated. 2024-12-03 15:49:59 +03:30
93c146391c Tutorial is Updated 2024-12-03 13:31:17 +03:30
5db98b9488 updated V1.0 2024-12-03 12:19:36 +03:30
c5ed2ad1e9 Merge pull request #125 from PhasicFlow/develop
Develop
2024-11-24 19:02:35 +03:30
a6a2d120e7 Merge pull request #123 from ramin1728/main
binarySystemOfParticles tutorial Updated for v1
2024-11-24 19:01:49 +03:30
7f0cb9cdb5 Merge pull request #124 from PhasicFlow/main
merge from main to develop to update everyting
2024-11-24 19:00:53 +03:30
24fe136bb2 Add files via upload 2024-11-24 16:57:44 +03:30
c21d63dacf Add files via upload 2024-11-24 16:13:50 +03:30
a247243d45 Add files via upload 2024-11-24 16:12:56 +03:30
c8132d5067 Add files via upload 2024-11-24 15:58:26 +03:30
51e43f089e Add files via upload 2024-11-24 15:54:02 +03:30
a606e48e66 Merge pull request #122 from PhasicFlow/develop
bug fix for postprocessPhasicFlow
2024-11-22 21:45:47 +03:30
HRN
3e9ac75c93 Merge branch 'develop' of github.com:PhasicFlow/phasicFlow into develop 2024-11-22 21:41:31 +03:30
HRN
5df18532b4 bug fix for postProcess 2024-11-22 21:33:46 +03:30
a113a350d2 Merge pull request #120 from Alireza77h/main
course graining tutorial added
2024-11-18 20:43:42 +03:30
64cda707a9 Merge pull request #121 from PhasicFlow/develop
Utility postProcess is modified to be used in Version 1.0
2024-11-18 20:40:55 +03:30
fdb50d5af8 Merge branch 'main' into develop 2024-11-18 20:40:36 +03:30
HRN
75fba2710e Utility postProcess is modified to be used in Version 1.0
- counting dictionary is added to postProcess
- Code need to be modified to cleaned (fields in the repository should be removed onces postProcess is passed the next time step)
2024-11-18 20:27:44 +03:30
5a00361df6 cours graining tutorial added 2024-11-17 14:59:58 +03:30
192ff67bff Update bashrc
version is added to the end of file name
2024-11-17 14:47:49 +03:30
66539ae97a Merge pull request #119 from Alireza77h/main
course graining added
2024-11-17 14:42:06 +03:30
2bd95b933f course graining added 2024-11-17 10:19:40 +03:30
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
673 changed files with 48667 additions and 21798 deletions

6
.gitignore vendored
View File

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

View File

@ -1,75 +1,73 @@
cmake_minimum_required(VERSION 3.16 FATAL_ERROR)
# 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_REQUIRED True)
set(CMAKE_INSTALL_PREFIX ${phasicFlow_SOURCE_DIR} CACHE PATH "Install path of phasicFlow" FORCE)
set(CMAKE_BUILD_TYPE Release CACHE STRING "build type" FORCE)
set(BUILD_SHARED_LIBS ON CACHE BOOL "Build using shared libraries" FORCE)
mark_as_advanced(FORCE var BUILD_SHARED_LIBS)
message(STATUS ${CMAKE_INSTALL_PREFIX})
mark_as_advanced(FORCE var Kokkos_ENABLE_CUDA_LAMBDA)
mark_as_advanced(FORCE var Kokkos_ENABLE_OPENMP)
mark_as_advanced(FORCE var Kokkos_ENABLE_SERIAL)
mark_as_advanced(FORCE var Kokkos_ENABLE_CUDA_LAMBDA)
mark_as_advanced(FORCE var BUILD_SHARED_LIBS)
include(cmake/globals.cmake)
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} ./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_OpenMP "Build phasicFlow and backends for OpenMP 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)
set(Kokkos_ENABLE_SERIAL ON CACHE BOOL "Serial execution" FORCE)
set(Kokkos_ENABLE_OPENMP OFF CACHE BOOL "OpenMP execution" FORCE)
set(Kokkos_ENABLE_CUDA OFF CACHE BOOL "Cuda execution" FORCE)
set(Kokkos_ENABLE_CUDA_LAMBDA OFF CACHE BOOL "Cuda execution" FORCE)
set(Kokkos_ENABLE_CUDA_CONSTEXPR OFF CACHE BOOL "Enable constexpr on cuda code" FORCE)
elseif(pFlow_Build_OpenMP )
set(Kokkos_ENABLE_SERIAL ON CACHE BOOL "Serial execution" FORCE)
set(Kokkos_ENABLE_OPENMP ON CACHE BOOL "OpenMP execution" FORCE)
set(Kokkos_ENABLE_CUDA OFF CACHE BOOL "Cuda execution" FORCE)
set(Kokkos_ENABLE_CUDA_LAMBDA OFF CACHE BOOL "Cuda execution" FORCE)
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)
set(Kokkos_ENABLE_SERIAL ON CACHE BOOL "Serial execution" FORCE)
set(Kokkos_ENABLE_OPENMP OFF CACHE BOOL "OpenMP execution" FORCE)
set(Kokkos_ENABLE_CUDA ON CACHE BOOL "Cuda execution" FORCE)
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()
include(cmake/globals.cmake)
message(STATUS "Valid file extensions are ${validFiles}")
if(pFlow_Build_MPI)
find_package(MPI REQUIRED)
endif()
include(cmake/makeLibraryGlobals.cmake)
include(cmake/makeExecutableGlobals.cmake)
configure_file(phasicFlowConfig.H.in phasicFlowConfig.H)
set(CMAKE_EXPORT_COMPILE_COMMANDS ON)
#add a global include directory
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(solvers)
@ -77,12 +75,9 @@ add_subdirectory(solvers)
add_subdirectory(utilities)
add_subdirectory(DEMSystems)
#add_subdirectory(testIO)
install(FILES "${PROJECT_BINARY_DIR}/phasicFlowConfig.H"
DESTINATION include
)
DESTINATION include)
include(InstallRequiredSystemLibraries)
set(CPACK_RESOURCE_FILE_LICENSE "${CMAKE_CURRENT_SOURCE_DIR}/LICENSE")

View File

@ -1,9 +1,11 @@
set(SourceFiles
DEMSystem/DEMSystem.cpp
sphereDEMSystem/sphereDEMSystem.cpp
sphereDEMSystem/sphereFluidParticles.cpp
domainDistribute/domainDistribute.cpp
DEMSystem/DEMSystem.cpp
sphereDEMSystem/sphereFluidParticles.cpp
sphereDEMSystem/sphereDEMSystem.cpp
grainDEMSystem/grainFluidParticles.cpp
grainDEMSystem/grainDEMSystem.cpp
)
set(link_libs Kokkos::kokkos phasicFlow Particles Geometry Property Interaction Interaction Utilities)

View File

@ -33,14 +33,14 @@ pFlow::DEMSystem::DEMSystem(
ControlDict_()
{
REPORT(0)<<"Initializing host/device execution spaces . . . \n";
REPORT(1)<<"Host execution space is "<< greenText(DefaultHostExecutionSpace::name())<<endREPORT;
REPORT(1)<<"Device execution space is "<<greenText(DefaultExecutionSpace::name())<<endREPORT;
REPORT(0)<<"Initializing host/device execution spaces . . . \n";
REPORT(1)<<"Host execution space is "<< Green_Text(DefaultHostExecutionSpace::name())<<END_REPORT;
REPORT(1)<<"Device execution space is "<<Green_Text(DefaultExecutionSpace::name())<<END_REPORT;
// initialize Kokkos
Kokkos::initialize( argc, argv );
REPORT(0)<<"\nCreating Control repository . . ."<<endREPORT;
REPORT(0)<<"\nCreating Control repository . . ."<<END_REPORT;
Control_ = makeUnique<systemControl>(
ControlDict_.startTime(),
ControlDict_.endTime(),
@ -87,4 +87,3 @@ pFlow::uniquePtr<pFlow::DEMSystem>
return nullptr;
}

View File

@ -25,6 +25,7 @@ Licence:
#include "types.hpp"
#include "span.hpp"
#include "box.hpp"
#include "virtualConstructor.hpp"
#include "uniquePtr.hpp"
#include "systemControl.hpp"
@ -60,6 +61,7 @@ public:
DEMSystem(const DEMSystem&)=delete;
/// @brief no assignment
DEMSystem& operator = (const DEMSystem&)=delete;
create_vCtor(
@ -111,19 +113,34 @@ public:
int32 numParInDomain(int32 di)const = 0;
virtual
std::vector<int32> numParInDomain()const = 0;
std::vector<int32> numParInDomains()const = 0;
virtual
span<const int32> parIndexInDomain(int32 di)const = 0;
span<const int32> parIndexInDomain(int32 domIndx)const = 0;
virtual
span<real> parDiameter() = 0;
span<real> diameter() = 0;
virtual
span<realx3> parVelocity() = 0;
span<real> courseGrainFactor() = 0;
virtual
span<realx3> parPosition() = 0;
span<realx3> acceleration()=0;
virtual
span<realx3> velocity() = 0;
virtual
span<realx3> position() = 0;
virtual
span<realx3> rAcceleration()=0;
virtual
span<realx3> rVelocity() = 0;
virtual
span<realx3> rPosition() = 0;
virtual
span<realx3> parFluidForce() = 0;
@ -153,7 +170,6 @@ public:
bool iterate(real upToTime) = 0;
static
uniquePtr<DEMSystem>
create(
@ -162,8 +178,6 @@ public:
int argc,
char* argv[]);
};

View File

@ -33,7 +33,7 @@ void pFlow::domainDistribute::clcDomains(const std::vector<box>& domains)
pFlow::domainDistribute::domainDistribute(
const Vector<box>& domains,
const std::vector<box>& domains,
real maxBoundingBox)
:
numDomains_(domains.size()),
@ -47,10 +47,10 @@ maxBoundingBoxSize_(maxBoundingBox)
}
bool pFlow::domainDistribute::locateParticles(
ViewType1D<realx3,HostSpace> points, includeMask mask)
ViewType1D<realx3,HostSpace> points, const pFlagTypeHost& mask)
{
range activeRange = mask.activeRange();
const rangeU32 activeRange = mask.activeRange();
for(int32 di=0; di<numDomains_; di++)
@ -59,7 +59,7 @@ bool pFlow::domainDistribute::locateParticles(
}
for(int32 i=activeRange.first; i<activeRange.second; i++)
for(int32 i=activeRange.start(); i<activeRange.end(); i++)
{
if(mask(i))
{

View File

@ -43,19 +43,16 @@ protected:
int32Vector_H numParInDomain_;
real maxBoundingBoxSize_;
real domainExtension_ = 1.0;
using includeMask = typename pointStructure::activePointsHost;
void clcDomains(const std::vector<box>& domains);
public:
domainDistribute(
const Vector<box>& domains,
const std::vector<box>& domains,
real maxBoundingBox);
~domainDistribute()=default;
@ -78,7 +75,7 @@ public:
{
return
span<const int32>(
particlesInDomains_[di].hostVectorAll().data(),
particlesInDomains_[di].hostViewAll().data(),
numParInDomain_[di]
);
}
@ -91,7 +88,7 @@ public:
//template<typename includeMask>
bool locateParticles(
ViewType1D<realx3,HostSpace> points, includeMask mask);
ViewType1D<realx3,HostSpace> points, const pFlagTypeHost& mask);
};

View File

@ -0,0 +1,273 @@
/*------------------------------- 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 "grainDEMSystem.hpp"
#include "vocabs.hpp"
bool pFlow::grainDEMSystem::loop()
{
do
{
//
if(! insertion_().insertParticles(
Control().time().currentIter(),
Control().time().currentTime(),
Control().time().dt() ) )
{
fatalError<<
"particle insertion failed in grainDEMSystem.\n";
return false;
}
geometry_->beforeIteration();
interaction_->beforeIteration();
particles_->beforeIteration();
interaction_->iterate();
particles_->iterate();
geometry_->iterate();
particles_->afterIteration();
geometry_->afterIteration();
}while(Control()++);
return true;
}
pFlow::grainDEMSystem::grainDEMSystem(
word demSystemName,
const std::vector<box>& domains,
int argc,
char* argv[],
bool requireRVel)
:
DEMSystem(demSystemName, domains, argc, argv),
requireRVel_(requireRVel)
{
REPORT(0)<<"\nReading proprties . . . "<<END_REPORT;
property_ = makeUnique<property>(
propertyFile__,
Control().caseSetup().path());
REPORT(0)<< "\nCreating surface geometry for grainDEMSystem . . . "<<END_REPORT;
geometry_ = geometry::create(Control(), Property());
REPORT(0)<<"Reading shape dictionary ..."<<END_REPORT;
grains_ = makeUnique<grainShape>(
pFlow::shapeFile__,
&Control().caseSetup(),
Property() );
REPORT(0)<<"\nReading grain particles . . ."<<END_REPORT;
particles_ = makeUnique<grainFluidParticles>(Control(), grains_());
insertion_ = makeUnique<grainInsertion>(
particles_(),
particles_().grains());
REPORT(0)<<"\nCreating interaction model for grain-grain contact and grain-wall contact . . ."<<END_REPORT;
interaction_ = interaction::create(
Control(),
Particles(),
Geometry());
auto maxD = maxBounndingSphereSize();
particleDistribution_ = makeUnique<domainDistribute>(domains, maxD);
}
pFlow::grainDEMSystem::~grainDEMSystem()
{
}
bool pFlow::grainDEMSystem::updateParticleDistribution(
real extentFraction,
const std::vector<box> domains)
{
if(!particleDistribution_->changeDomainsSize(
extentFraction,
maxBounndingSphereSize(),
domains))
{
fatalErrorInFunction<<
"cannot change the domain size"<<endl;
return false;
}
if(!particleDistribution_->locateParticles(
positionHost_,
particles_->pStruct().activePointsMaskHost()))
{
fatalErrorInFunction<<
"error in locating particles among sub-domains"<<endl;
return false;
}
return true;
}
pFlow::int32
pFlow::grainDEMSystem::numParInDomain(int32 di)const
{
return particleDistribution_().numParInDomain(di);
}
std::vector<pFlow::int32>
pFlow::grainDEMSystem::numParInDomains()const
{
const auto& distribute = particleDistribution_();
int32 numDomains = distribute.numDomains();
std::vector<int32> nums(numDomains);
for(int32 i=0; i<numDomains; i++)
{
nums[i] = distribute.numParInDomain(i);
}
return nums;
}
pFlow::span<const pFlow::int32>
pFlow::grainDEMSystem::parIndexInDomain(int32 di)const
{
return particleDistribution_->particlesInDomain(di);
}
pFlow::span<pFlow::real> pFlow::grainDEMSystem::diameter()
{
return span<real>(diameterHost_.data(), diameterHost_.size());
}
pFlow::span<pFlow::real> pFlow::grainDEMSystem::courseGrainFactor()
{
return span<real>(particles_->courseGrainFactorHost().data(), particles_->courseGrainFactorHost().size());
}
pFlow::span<pFlow::realx3> pFlow::grainDEMSystem::acceleration()
{
return span<realx3>(nullptr, 0);
}
pFlow::span<pFlow::realx3> pFlow::grainDEMSystem::velocity()
{
return span<realx3>(velocityHost_.data(), velocityHost_.size());
}
pFlow::span<pFlow::realx3> pFlow::grainDEMSystem::position()
{
return span<realx3>(positionHost_.data(), positionHost_.size());
}
pFlow::span<pFlow::realx3> pFlow::grainDEMSystem::rAcceleration()
{
return span<realx3>(nullptr, 0);
}
pFlow::span<pFlow::realx3> pFlow::grainDEMSystem::rVelocity()
{
return span<realx3>(rVelocityHost_.data(), rVelocityHost_.size());
}
pFlow::span<pFlow::realx3> pFlow::grainDEMSystem::rPosition()
{
return span<realx3>(nullptr, 0);
}
pFlow::span<pFlow::realx3> pFlow::grainDEMSystem::parFluidForce()
{
auto& hVec = particles_->fluidForceHost();
return span<realx3>(hVec.data(), hVec.size());
}
pFlow::span<pFlow::realx3> pFlow::grainDEMSystem::parFluidTorque()
{
auto& hVec = particles_->fluidTorqueHost();
return span<realx3>(hVec.data(), hVec.size());
}
bool pFlow::grainDEMSystem::sendFluidForceToDEM()
{
particles_->fluidForceHostUpdatedSync();
return true;
}
bool pFlow::grainDEMSystem::sendFluidTorqueToDEM()
{
particles_->fluidTorqueHostUpdatedSync();
return true;
}
bool pFlow::grainDEMSystem::beforeIteration()
{
velocityHost_ = std::as_const(particles_()).velocity().hostView();
positionHost_ = std::as_const(particles_()).pointPosition().hostView();
diameterHost_ = particles_->diameter().hostView();
if(requireRVel_)
rVelocityHost_ = std::as_const(particles_()).rVelocity().hostView();
return true;
}
bool pFlow::grainDEMSystem::iterate(
real upToTime,
real timeToWrite,
word timeName)
{
Control().time().setStopAt(upToTime);
Control().time().setOutputToFile(timeToWrite, timeName);
return loop();
return true;
}
bool pFlow::grainDEMSystem::iterate(real upToTime)
{
Control().time().setStopAt(upToTime);
return loop();
return true;
}
pFlow::real
pFlow::grainDEMSystem::maxBounndingSphereSize()const
{
real minD, maxD;
particles_->boundingSphereMinMax(minD, maxD);
return maxD;
}

View File

@ -0,0 +1,165 @@
/*------------------------------- 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 __grainDEMSystem_hpp__
#define __grainDEMSystem_hpp__
#include "DEMSystem.hpp"
#include "property.hpp"
#include "uniquePtr.hpp"
#include "geometry.hpp"
#include "grainFluidParticles.hpp"
#include "interaction.hpp"
#include "Insertions.hpp"
#include "domainDistribute.hpp"
namespace pFlow
{
class grainDEMSystem
:
public DEMSystem
{
protected:
// protected members
uniquePtr<property> property_ = nullptr;
uniquePtr<geometry> geometry_ = nullptr;
uniquePtr<grainShape> grains_ = nullptr;
uniquePtr<grainFluidParticles> particles_ = nullptr;
uniquePtr<grainInsertion> insertion_ = nullptr;
uniquePtr<interaction> interaction_ = nullptr;
uniquePtr<domainDistribute> particleDistribution_=nullptr;
// to be used for CPU communications
ViewType1D<realx3, HostSpace> velocityHost_;
ViewType1D<realx3, HostSpace> positionHost_;
ViewType1D<real, HostSpace> diameterHost_;
bool requireRVel_ = false;
ViewType1D<realx3, HostSpace> rVelocityHost_;
// protected member functions
auto& Property()
{
return property_();
}
auto& Geometry()
{
return geometry_();
}
auto& Particles()
{
return particles_();
}
auto& Interaction()
{
return interaction_();
}
bool loop();
public:
TypeInfo("grainDEMSystem");
grainDEMSystem(
word demSystemName,
const std::vector<box>& domains,
int argc,
char* argv[],
bool requireRVel=false);
virtual ~grainDEMSystem();
grainDEMSystem(const grainDEMSystem&)=delete;
grainDEMSystem& operator = (const grainDEMSystem&)=delete;
add_vCtor(
DEMSystem,
grainDEMSystem,
word);
bool updateParticleDistribution(real extentFraction, const std::vector<box> domains) override;
int32 numParInDomain(int32 di)const override;
std::vector<int32> numParInDomains()const override;
span<const int32> parIndexInDomain(int32 di)const override;
span<real> diameter() override;
span<real> courseGrainFactor() override;
span<realx3> acceleration() override;
span<realx3> velocity() override;
span<realx3> position() override;
span<realx3> rAcceleration() override;
span<realx3> rVelocity() override;
span<realx3> rPosition() override;
span<realx3> parFluidForce() override;
span<realx3> parFluidTorque() override;
bool sendFluidForceToDEM() override;
bool sendFluidTorqueToDEM() override;
bool beforeIteration() override;
bool iterate(
real upToTime,
real timeToWrite,
word timeName) override;
bool iterate(real upToTime) override;
real maxBounndingSphereSize()const override;
};
} // pFlow
#endif

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.
-----------------------------------------------------------------------------*/
#include "grainFluidParticles.hpp"
#include "grainFluidParticlesKernels.hpp"
void pFlow::grainFluidParticles::checkHostMemory()
{
if(fluidForce_.size()!=fluidForceHost_.size())
{
resizeNoInit(fluidForceHost_, fluidForce_.size());
resizeNoInit(fluidTorqueHost_, fluidTorque_.size());
}
// copy the data (if required) from device to host
courseGrainFactorHost_ = coarseGrainFactor().hostView();
}
pFlow::grainFluidParticles::grainFluidParticles(
systemControl &control,
const grainShape& grains)
: grainParticles(control, grains),
fluidForce_(
objectFile(
"fluidForce",
"",
objectFile::READ_IF_PRESENT,
objectFile::WRITE_ALWAYS),
dynPointStruct(),
zero3),
fluidTorque_(
objectFile(
"fluidTorque",
"",
objectFile::READ_IF_PRESENT,
objectFile::WRITE_NEVER),
dynPointStruct(),
zero3)
{
checkHostMemory();
}
bool pFlow::grainFluidParticles::beforeIteration()
{
grainParticles::beforeIteration();
checkHostMemory();
return true;
}
bool pFlow::grainFluidParticles::iterate()
{
const auto ti = this->TimeInfo();
accelerationTimer().start();
pFlow::grainFluidParticlesKernels::acceleration(
control().g(),
mass().deviceViewAll(),
contactForce().deviceViewAll(),
fluidForce_.deviceViewAll(),
I().deviceViewAll(),
contactTorque().deviceViewAll(),
fluidTorque_.deviceViewAll(),
pStruct().activePointsMaskDevice(),
acceleration().deviceViewAll(),
rAcceleration().deviceViewAll()
);
accelerationTimer().end();
intCorrectTimer().start();
dynPointStruct().correct(ti.dt());
rVelIntegration().correct(ti.dt(), rVelocity(), rAcceleration());
intCorrectTimer().end();
return true;
}
void pFlow::grainFluidParticles::fluidForceHostUpdatedSync()
{
copy(fluidForce_.deviceView(), fluidForceHost_);
return;
}
void pFlow::grainFluidParticles::fluidTorqueHostUpdatedSync()
{
copy(fluidTorque_.deviceView(), fluidTorqueHost_);
return;
}

View File

@ -0,0 +1,105 @@
/*------------------------------- 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.
-----------------------------------------------------------------------------*/
/*!
@class pFlow::grainFluidParticles
@brief Class for managing grain particles with fluid interactions
This is a top-level class that contains the essential components for
defining grain prticles with fluid interactions in a DEM simulation.
*/
#ifndef __grainFluidParticles_hpp__
#define __grainFluidParticles_hpp__
#include "grainParticles.hpp"
namespace pFlow
{
class grainFluidParticles
:
public grainParticles
{
protected:
/// pointField of rotational acceleration of particles on device
realx3PointField_D fluidForce_;
hostViewType1D<realx3> fluidForceHost_;
realx3PointField_D fluidTorque_;
hostViewType1D<realx3> fluidTorqueHost_;
hostViewType1D<real> courseGrainFactorHost_;
void checkHostMemory();
public:
/// construct from systemControl and property
grainFluidParticles(systemControl &control, const grainShape& grains);
~grainFluidParticles() override = default;
/// before iteration step
bool beforeIteration() override;
/// iterate particles
bool iterate() override;
auto& fluidForce()
{
return fluidForce_;
}
auto& fluidTorque()
{
return fluidTorque_;
}
auto& fluidForceHost()
{
return fluidForceHost_;
}
auto& fluidTorqueHost()
{
return fluidTorqueHost_;
}
auto& courseGrainFactorHost()
{
return courseGrainFactorHost_;
}
void fluidForceHostUpdatedSync();
void fluidTorqueHostUpdatedSync();
}; //grainFluidParticles
} // pFlow
#endif //__sphereFluidParticles_hpp__

View File

@ -0,0 +1,79 @@
/*------------------------------- 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 __grainFluidParticlesKernels_hpp__
#define __grainFluidParticlesKernels_hpp__
namespace pFlow::grainFluidParticlesKernels
{
using rpAcceleration = Kokkos::RangePolicy<
DefaultExecutionSpace,
Kokkos::Schedule<Kokkos::Static>,
Kokkos::IndexType<int32>
>;
template<typename IncludeFunctionType>
void acceleration(
realx3 g,
deviceViewType1D<real> mass,
deviceViewType1D<realx3> force,
deviceViewType1D<realx3> fluidForce,
deviceViewType1D<real> I,
deviceViewType1D<realx3> torque,
deviceViewType1D<realx3> fluidTorque,
IncludeFunctionType incld,
deviceViewType1D<realx3> lAcc,
deviceViewType1D<realx3> rAcc
)
{
auto activeRange = incld.activeRange();
if(incld.isAllActive())
{
Kokkos::parallel_for(
"pFlow::grainParticlesKernels::acceleration",
rpAcceleration(activeRange.first, activeRange.second),
LAMBDA_HD(int32 i){
lAcc[i] = (force[i]+fluidForce[i])/mass[i] + g;
rAcc[i] = (torque[i]+fluidTorque[i])/I[i];
});
}
else
{
Kokkos::parallel_for(
"pFlow::grainParticlesKernels::acceleration",
rpAcceleration(activeRange.first, activeRange.second),
LAMBDA_HD(int32 i){
if(incld(i))
{
lAcc[i] = (force[i]+fluidForce[i])/mass[i] + g;
rAcc[i] = (torque[i]+fluidTorque[i])/I[i];
}
});
}
Kokkos::fence();
}
}
#endif

View File

@ -19,6 +19,7 @@ Licence:
-----------------------------------------------------------------------------*/
#include "sphereDEMSystem.hpp"
#include "vocabs.hpp"
bool pFlow::sphereDEMSystem::loop()
{
@ -26,17 +27,16 @@ bool pFlow::sphereDEMSystem::loop()
do
{
//
if(! insertion_().insertParticles(
Control().time().currentIter(),
Control().time().currentTime(),
Control().time().dt() ) )
{
fatalError<<
"particle insertion failed in sphereDFlow solver.\n";
"particle insertion failed in sphereDEMSystem.\n";
return false;
}
geometry_->beforeIteration();
interaction_->beforeIteration();
@ -63,29 +63,37 @@ pFlow::sphereDEMSystem::sphereDEMSystem(
word demSystemName,
const std::vector<box>& domains,
int argc,
char* argv[])
char* argv[],
bool requireRVel)
:
DEMSystem(demSystemName, domains, argc, argv)
DEMSystem(demSystemName, domains, argc, argv),
requireRVel_(requireRVel)
{
REPORT(0)<<"\nReading proprties . . . "<<endREPORT;
REPORT(0)<<"\nReading proprties . . . "<<END_REPORT;
property_ = makeUnique<property>(
Control().caseSetup().path()+propertyFile__);
propertyFile__,
Control().caseSetup().path());
REPORT(0)<< "\nCreating surface geometry for sphereDEMSystem . . . "<<endREPORT;
REPORT(0)<< "\nCreating surface geometry for sphereDEMSystem . . . "<<END_REPORT;
geometry_ = geometry::create(Control(), Property());
REPORT(0)<<"Reading shapes dictionary..."<<END_REPORT;
spheres_ = makeUnique<sphereShape>(
pFlow::shapeFile__,
&Control().caseSetup(),
Property());
REPORT(0)<<"\nReading sphere particles . . ."<<endREPORT;
particles_ = makeUnique<sphereFluidParticles>(Control(), Property());
REPORT(0)<<"\nReading sphere particles . . ."<<END_REPORT;
particles_ = makeUnique<sphereFluidParticles>(Control(), spheres_());
insertion_ = makeUnique<sphereInsertion>(
Control().caseSetup().path()+insertionFile__,
particles_(),
particles_().shapes());
particles_().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;
interaction_ = interaction::create(
Control(),
Particles(),
@ -98,6 +106,7 @@ pFlow::sphereDEMSystem::sphereDEMSystem(
}
pFlow::sphereDEMSystem::~sphereDEMSystem()
{
@ -119,8 +128,8 @@ bool pFlow::sphereDEMSystem::updateParticleDistribution(
}
if(!particleDistribution_->locateParticles(
parPosition_,
particles_->pStruct().activePointsMaskH()))
positionHost_,
particles_->pStruct().activePointsMaskHost()))
{
fatalErrorInFunction<<
"error in locating particles among sub-domains"<<endl;
@ -137,7 +146,7 @@ pFlow::int32
}
std::vector<pFlow::int32>
pFlow::sphereDEMSystem::numParInDomain()const
pFlow::sphereDEMSystem::numParInDomains()const
{
const auto& distribute = particleDistribution_();
int32 numDomains = distribute.numDomains();
@ -156,31 +165,56 @@ pFlow::sphereDEMSystem::parIndexInDomain(int32 di)const
return particleDistribution_->particlesInDomain(di);
}
pFlow::span<pFlow::real> pFlow::sphereDEMSystem::parDiameter()
pFlow::span<pFlow::real> pFlow::sphereDEMSystem::diameter()
{
return span<real>(parDiameter_.data(), parDiameter_.size());
return span<real>(diameterHost_.data(), diameterHost_.size());
}
pFlow::span<pFlow::realx3> pFlow::sphereDEMSystem::parVelocity()
pFlow::span<pFlow::real> pFlow::sphereDEMSystem::courseGrainFactor()
{
return span<realx3>(parVelocity_.data(), parVelocity_.size());
return span<real>(nullptr, 0);
}
pFlow::span<pFlow::realx3> pFlow::sphereDEMSystem::parPosition()
pFlow::span<pFlow::realx3> pFlow::sphereDEMSystem::acceleration()
{
return span<realx3>(parPosition_.data(), parPosition_.size());
return span<realx3>(nullptr, 0);
}
pFlow::span<pFlow::realx3> pFlow::sphereDEMSystem::velocity()
{
return span<realx3>(velocityHost_.data(), velocityHost_.size());
}
pFlow::span<pFlow::realx3> pFlow::sphereDEMSystem::position()
{
return span<realx3>(positionHost_.data(), positionHost_.size());
}
pFlow::span<pFlow::realx3> pFlow::sphereDEMSystem::rAcceleration()
{
return span<realx3>(nullptr, 0);
}
pFlow::span<pFlow::realx3> pFlow::sphereDEMSystem::rVelocity()
{
return span<realx3>(rVelocityHost_.data(), rVelocityHost_.size());
}
pFlow::span<pFlow::realx3> pFlow::sphereDEMSystem::rPosition()
{
return span<realx3>(nullptr, 0);
}
pFlow::span<pFlow::realx3> pFlow::sphereDEMSystem::parFluidForce()
{
auto& hVec = particles_->fluidForceHostAll();
auto& hVec = particles_->fluidForceHost();
return span<realx3>(hVec.data(), hVec.size());
}
pFlow::span<pFlow::realx3> pFlow::sphereDEMSystem::parFluidTorque()
{
auto& hVec = particles_->fluidTorqueHostAll();
auto& hVec = particles_->fluidTorqueHost();
return span<realx3>(hVec.data(), hVec.size());
}
@ -198,9 +232,14 @@ bool pFlow::sphereDEMSystem::sendFluidTorqueToDEM()
bool pFlow::sphereDEMSystem::beforeIteration()
{
parVelocity_ = particles_->dynPointStruct().velocityHostAll();
parPosition_ = particles_->dynPointStruct().pointPositionHostAll();
parDiameter_ = particles_->diameter().hostVectorAll();
velocityHost_ = std::as_const(particles_()).velocity().hostView();
positionHost_ = std::as_const(particles_()).pointPosition().hostView();
diameterHost_ = particles_->diameter().hostView();
if(requireRVel_)
rVelocityHost_ = std::as_const(particles_()).rVelocity().hostView();
return true;
}

View File

@ -42,24 +42,31 @@ protected:
// protected members
uniquePtr<property> property_ = nullptr;
uniquePtr<property> property_ = nullptr;
uniquePtr<geometry> geometry_ = nullptr;
uniquePtr<geometry> geometry_ = nullptr;
uniquePtr<sphereFluidParticles> particles_ = nullptr;
uniquePtr<sphereShape> spheres_ = nullptr;
uniquePtr<sphereInsertion> insertion_ = nullptr;
uniquePtr<sphereFluidParticles> particles_ = nullptr;
uniquePtr<interaction> interaction_ = nullptr;
uniquePtr<sphereInsertion> insertion_ = nullptr;
uniquePtr<domainDistribute> particleDistribution_=nullptr;
uniquePtr<interaction> interaction_ = nullptr;
uniquePtr<domainDistribute> particleDistribution_=nullptr;
// to be used for CPU communications
ViewType1D<realx3, HostSpace> parVelocity_;
ViewType1D<realx3, HostSpace> velocityHost_;
ViewType1D<realx3, HostSpace> parPosition_;
ViewType1D<realx3, HostSpace> positionHost_;
ViewType1D<real, HostSpace> diameterHost_;
bool requireRVel_ = false;
ViewType1D<realx3, HostSpace> rVelocityHost_;
ViewType1D<real, HostSpace> parDiameter_;
// protected member functions
auto& Property()
@ -92,7 +99,8 @@ public:
word demSystemName,
const std::vector<box>& domains,
int argc,
char* argv[]);
char* argv[],
bool requireRVel=false);
virtual ~sphereDEMSystem();
@ -110,15 +118,25 @@ public:
int32 numParInDomain(int32 di)const override;
std::vector<int32> numParInDomain()const override;
std::vector<int32> numParInDomains()const override;
span<const int32> parIndexInDomain(int32 di)const override;
span<real> parDiameter() override;
span<real> diameter() override;
span<realx3> parVelocity() override;
span<real> courseGrainFactor() override;
span<realx3> parPosition() override;
span<realx3> acceleration() override;
span<realx3> velocity() override;
span<realx3> position() override;
span<realx3> rAcceleration() override;
span<realx3> rVelocity() override;
span<realx3> rPosition() override;
span<realx3> parFluidForce() override;

View File

@ -21,85 +21,85 @@ Licence:
#include "sphereFluidParticles.hpp"
#include "sphereFluidParticlesKernels.hpp"
void pFlow::sphereFluidParticles::checkHostMemory()
{
if(fluidForce_.size()!=fluidForceHost_.size())
{
resizeNoInit(fluidForceHost_, fluidForce_.size());
resizeNoInit(fluidTorqueHost_, fluidTorque_.size());
}
}
pFlow::sphereFluidParticles::sphereFluidParticles(
systemControl &control,
const property& prop
)
:
sphereParticles(control, prop),
fluidForce_(
this->time().emplaceObject<realx3PointField_HD>(
objectFile(
"fluidForce",
"",
objectFile::READ_IF_PRESENT,
objectFile::WRITE_ALWAYS
),
pStruct(),
zero3
)
),
fluidTorque_(
this->time().emplaceObject<realx3PointField_HD>(
objectFile(
"fluidTorque",
"",
objectFile::READ_IF_PRESENT,
objectFile::WRITE_ALWAYS
),
pStruct(),
zero3
)
)
{}
systemControl &control,
const sphereShape& shpShape)
: sphereParticles(control, shpShape),
fluidForce_(
objectFile(
"fluidForce",
"",
objectFile::READ_IF_PRESENT,
objectFile::WRITE_ALWAYS),
dynPointStruct(),
zero3),
fluidTorque_(
objectFile(
"fluidTorque",
"",
objectFile::READ_IF_PRESENT,
objectFile::WRITE_NEVER),
dynPointStruct(),
zero3)
{
checkHostMemory();
}
bool pFlow::sphereFluidParticles::beforeIteration()
{
sphereParticles::beforeIteration();
checkHostMemory();
return true;
}
bool pFlow::sphereFluidParticles::iterate()
{
accelerationTimer_.start();
const auto ti = this->TimeInfo();
accelerationTimer().start();
pFlow::sphereFluidParticlesKernels::acceleration(
control().g(),
mass().deviceVectorAll(),
contactForce().deviceVectorAll(),
fluidForce().deviceVectorAll(),
I().deviceVectorAll(),
contactTorque().deviceVectorAll(),
fluidTorque().deviceVectorAll(),
pStruct().activePointsMaskD(),
accelertion().deviceVectorAll(),
rAcceleration().deviceVectorAll()
mass().deviceViewAll(),
contactForce().deviceViewAll(),
fluidForce_.deviceViewAll(),
I().deviceViewAll(),
contactTorque().deviceViewAll(),
fluidTorque_.deviceViewAll(),
pStruct().activePointsMaskDevice(),
acceleration().deviceViewAll(),
rAcceleration().deviceViewAll()
);
accelerationTimer_.end();
accelerationTimer().end();
intCorrectTimer_.start();
intCorrectTimer().start();
dynPointStruct_.correct(this->dt(), accelertion_);
dynPointStruct().correct(ti.dt());
rVelIntegration_().correct(this->dt(), rVelocity_, rAcceleration_);
rVelIntegration().correct(ti.dt(), rVelocity(), rAcceleration());
intCorrectTimer_.end();
intCorrectTimer().end();
return true;
}
void pFlow::sphereFluidParticles::fluidForceHostUpdatedSync()
{
fluidForce_.modifyOnHost();
fluidForce_.syncViews();
copy(fluidForce_.deviceView(), fluidForceHost_);
return;
}
void pFlow::sphereFluidParticles::fluidTorqueHostUpdatedSync()
{
fluidTorque_.modifyOnHost();
fluidTorque_.syncViews();
copy(fluidTorque_.deviceView(), fluidTorqueHost_);
return;
}

View File

@ -43,12 +43,17 @@ class sphereFluidParticles
protected:
/// pointField of rotational acceleration of particles on device
realx3PointField_HD& fluidForce_;
realx3PointField_D fluidForce_;
realx3PointField_HD& fluidTorque_;
hostViewType1D<realx3> fluidForceHost_;
realx3PointField_D fluidTorque_;
void zeroFluidForce_H()
hostViewType1D<realx3> fluidTorqueHost_;
void checkHostMemory();
/*void zeroFluidForce_H()
{
fluidForce_.fillHost(zero3);
}
@ -56,12 +61,12 @@ protected:
void zeroFluidTorque_H()
{
fluidTorque_.fillHost(zero3);
}
}*/
public:
/// construct from systemControl and property
sphereFluidParticles(systemControl &control, const property& prop);
sphereFluidParticles(systemControl &control, const sphereShape& shpShape);
/// before iteration step
bool beforeIteration() override;
@ -81,17 +86,16 @@ public:
}
auto& fluidForceHostAll()
auto& fluidForceHost()
{
return fluidForce_.hostVectorAll();
return fluidForceHost_;
}
auto& fluidTorqueHostAll()
auto& fluidTorqueHost()
{
return fluidTorque_.hostVectorAll();
return fluidTorqueHost_;
}
void fluidForceHostUpdatedSync();
void fluidTorqueHostUpdatedSync();

View File

@ -46,7 +46,7 @@ void acceleration(
{
auto activeRange = incld.activeRange();
if(incld.allActive())
if(incld.isAllActive())
{
Kokkos::parallel_for(
"pFlow::sphereParticlesKernels::acceleration",

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.
**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?
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
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/)

71
cmake/autoComplete Normal file
View File

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

View File

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

View File

@ -2,3 +2,52 @@
set(validFiles)
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,8 +31,8 @@ target_include_directories(${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} source files are: ${source_files}")
message(STATUS " ${target_name} include dirs are: ${includeDirs}\n")
#message(STATUS " ${target_name} source files are: ${${source_files}}")
#message(STATUS " ${target_name} include dirs are: ${includeDirs}\n")
install(TARGETS ${target_name} DESTINATION bin)

View File

@ -41,8 +41,8 @@ target_include_directories(${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} source files are: ${source_files}")
message(STATUS " ${target_name} include dirs are: ${includeDirs}\n")
#message(STATUS " ${target_name} source files are: ${source_files}")
#message(STATUS " ${target_name} include dirs are: ${includeDirs}\n \n")
install(TARGETS ${target_name} DESTINATION lib)
install(FILES ${includeFiles} DESTINATION include/${target_name})

View File

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

View File

@ -1,7 +1,9 @@
#add_subdirectory(iterateSphereParticles)
add_subdirectory(iterateSphereParticles)
add_subdirectory(iterateGeometry)
add_subdirectory(sphereGranFlow)
add_subdirectory(grainGranFlow)

View File

@ -0,0 +1,7 @@
set(source_files
grainGranFlow.cpp
)
set(link_lib Kokkos::kokkos phasicFlow Particles Geometry Property Interaction Interaction Utilities)
pFlow_make_executable_install(grainGranFlow source_files link_lib)

View File

@ -0,0 +1,57 @@
/*------------------------------- 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.
-----------------------------------------------------------------------------*/
REPORT(0)<<"Reading shape dictionary ..."<<END_REPORT;
pFlow::grainShape grains
(
pFlow::shapeFile__,
&Control.caseSetup(),
proprties
);
//
REPORT(0)<<"\nReading sphere particles . . ."<<END_REPORT;
pFlow::grainParticles grnParticles(Control, grains);
//
REPORT(0)<<"\nCreating particle insertion object . . ."<<END_REPORT;
/*auto& sphInsertion =
Control.caseSetup().emplaceObject<sphereInsertion>(
objectFile(
insertionFile__,
"",
objectFile::READ_ALWAYS,
objectFile::WRITE_ALWAYS
),
sphParticles,
sphParticles.shapes()
);*/
auto grnInsertion = pFlow::grainInsertion(
grnParticles,
grnParticles.grains());
REPORT(0)<<"\nCreating interaction model for sphere-sphere contact and sphere-wall contact . . ."<<END_REPORT;
auto interactionPtr = pFlow::interaction::create(
Control,
grnParticles,
surfGeometry
);
auto& grnInteraction = interactionPtr();

View File

@ -0,0 +1,123 @@
/*------------------------------- phasicFlow ---------------------------------
O C enter of
O O E ngineering and
O O M ultiscale modeling of
OOOOOOO F luid flow
------------------------------------------------------------------------------
Copyright (C): www.cemf.ir
email: hamid.r.norouzi AT gmail.com
------------------------------------------------------------------------------
Licence:
This file is part of phasicFlow code. It is a free software for simulating
granular and multiphase flows. You can redistribute it and/or modify it under
the terms of GNU General Public License v3 or any other later versions.
phasicFlow is distributed to help others in their research in the field of
granular and multiphase flows, but WITHOUT ANY WARRANTY; without even the
implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
-----------------------------------------------------------------------------*/
/**
* \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 "phasicFlowKokkos.hpp"
#include "systemControl.hpp"
#include "commandLine.hpp"
#include "property.hpp"
#include "geometry.hpp"
#include "grainParticles.hpp"
#include "interaction.hpp"
#include "Insertions.hpp"
//#include "readControlDict.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(
"grainGranFlow",
"DEM solver for non-cohesive spherical grain 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 "setSurfaceGeometry.hpp"
#include "createDEMComponents.hpp"
REPORT(0)<<"\nStart of time loop . . .\n"<<END_REPORT;
do
{
if(! grnInsertion.insertParticles(
Control.time().currentIter(),
Control.time().currentTime(),
Control.time().dt() ) )
{
fatalError<<
"particle insertion failed in grain DFlow solver.\n";
return 1;
}
// set force to zero
surfGeometry.beforeIteration();
// set force to zero, predict, particle deletion and etc.
grnParticles.beforeIteration();
grnInteraction.beforeIteration();
grnInteraction.iterate();
surfGeometry.iterate();
grnParticles.iterate();
grnInteraction.afterIteration();
surfGeometry.afterIteration();
grnParticles.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

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

@ -0,0 +1,33 @@
/*------------------------------- 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.
-----------------------------------------------------------------------------*/
REPORT(0)<<"Reading shape dictionary ..."<<END_REPORT;
pFlow::sphereShape spheres
(
pFlow::shapeFile__,
&Control.caseSetup(),
proprties
);
//
REPORT(0)<<"\nReading sphere particles . . ."<<END_REPORT;
pFlow::sphereParticles sphParticles(Control, spheres);
WARNING<<"Particle insertion has not been set yet!"<<END_WARNING;

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

@ -17,13 +17,21 @@ Licence:
implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
-----------------------------------------------------------------------------*/
REPORT(0)<<"Reading shapes dictionary..."<<END_REPORT;
pFlow::sphereShape spheres
(
pFlow::shapeFile__,
&Control.caseSetup(),
proprties
);
//
REPORT(0)<<"\nReading sphere particles . . ."<<END_REPORT;
pFlow::sphereParticles sphParticles(Control, spheres);
//
REPORT(0)<<"\nReading sphere particles . . ."<<endREPORT;
sphereParticles sphParticles(Control, proprties);
//
REPORT(0)<<"\nCreating particle insertion object . . ."<<endREPORT;
REPORT(0)<<"\nCreating particle insertion object . . ."<<END_REPORT;
/*auto& sphInsertion =
Control.caseSetup().emplaceObject<sphereInsertion>(
objectFile(
@ -36,13 +44,12 @@ REPORT(0)<<"\nCreating particle insertion object . . ."<<endREPORT;
sphParticles.shapes()
);*/
auto sphInsertion = sphereInsertion(
Control.caseSetup().path()+insertionFile__,
auto sphInsertion = pFlow::sphereInsertion(
sphParticles,
sphParticles.shapes());
sphParticles.spheres());
REPORT(0)<<"\nCreating interaction model for sphere-sphere contact and sphere-wall contact . . ."<<endREPORT;
auto interactionPtr = interaction::create(
REPORT(0)<<"\nCreating interaction model for sphere-sphere contact and sphere-wall contact . . ."<<END_REPORT;
auto interactionPtr = pFlow::interaction::create(
Control,
sphParticles,
surfGeometry

View File

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

View File

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

View File

@ -19,32 +19,40 @@ Licence:
-----------------------------------------------------------------------------*/
#include "geometry.hpp"
#include "systemControl.hpp"
#include "vocabs.hpp"
bool pFlow::geometry::findPropertyId()
bool pFlow::geometry::createPropertyId()
{
int8Vector propId(0, surface().capacity(),RESERVE());
propId.clear();
uint32 pId;
ForAll(matI, materialName_)
if(materialName_.size() != numSurfaces() )
{
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<<
"material name for the geometry is invalid: "<< materialName_[matI]<<endl;
return false;
fatalErrorInFunction<<
"Property/material name is invalid "<<materialName_[i]<<
". 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);
@ -53,83 +61,363 @@ bool pFlow::geometry::findPropertyId()
}
void pFlow::geometry::zeroForce()
{
contactForceWall_.fill(zero3);
}
pFlow::geometry::geometry
(
systemControl& control,
const property& prop
)
:
demGeometry(control),
multiTriSurface
(
objectFile
(
triSurfaceFile__,
"",
objectFile::READ_ALWAYS,
objectFile::WRITE_ALWAYS
),
&control.geometry()
),
demComponent
(
"geometry",
control
),
wallProperty_(prop),
geometryRepository_(control.geometry()),
triSurface_(
control.geometry().emplaceObject<multiTriSurface>(
objectFile(
"triSurface",
"",
objectFile::READ_ALWAYS,
objectFile::WRITE_ALWAYS
)
)
propertyId_
(
objectFile
(
"propertyId",
"",
objectFile::READ_NEVER,
objectFile::WRITE_NEVER
),
motionComponentName_(
control.geometry().emplaceObject<wordField>(
objectFile(
"motionComponentName",
"",
objectFile::READ_ALWAYS,
objectFile::WRITE_ALWAYS
),
"motionNamesList"
)
*this,
0u
),
contactForceWall_
(
objectFile
(
"contactForcWall",
"",
objectFile::READ_IF_PRESENT,
objectFile::WRITE_NEVER
),
materialName_(
control.geometry().emplaceObject<wordField>(
objectFile(
"materialName",
"",
objectFile::READ_ALWAYS,
objectFile::WRITE_ALWAYS
),
"materialNamesList"
)
*this,
zero3
),
normalStressWall_
(
objectFile
(
"normalStressWall",
"",
objectFile::READ_IF_PRESENT,
objectFile::WRITE_NEVER
),
propertyId_(
control.geometry().emplaceObject<int8TriSurfaceField_D>(
objectFile(
"propertyId",
"",
objectFile::READ_NEVER,
objectFile::WRITE_NEVER),
surface(),
0 ) ),
contactForceWall_(
control.geometry().emplaceObject<realx3TriSurfaceField_D>(
objectFile(
"contactForceWall",
"",
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) )
*this,
zero3
),
shearStressWall_
(
objectFile
(
"shearStressWall",
"",
objectFile::READ_IF_PRESENT,
objectFile::WRITE_NEVER
),
*this,
zero3
)
{
if(!findPropertyId())
readWholeObject_ = false;
if( !IOobject::readObject() )
{
fatalErrorInFunction<<
"Error in reading from file "<<IOobject::path()<<endl;
fatalExit;
}
readWholeObject_ = true;
if( this->numSurfaces() != motionComponentName_.size() )
{
fatalErrorInFunction<<
"Number of surfaces is not equal to number of motion component names"<<endl;
fatalExit;
}
if(!createPropertyId())
{
fatalExit;
}
}
pFlow::geometry::geometry
(
systemControl &control,
const property &prop,
multiTriSurface &surf,
const wordVector &motionCompName,
const wordVector &materialName,
const dictionary& motionDict
)
:
multiTriSurface
(
objectFile
(
triSurfaceFile__,
"",
objectFile::READ_NEVER,
objectFile::WRITE_ALWAYS
),
&control.geometry(),
surf
),
demComponent
(
"geometry",
control
),
wallProperty_
(
prop
),
propertyId_
(
objectFile
(
"propertyId",
"",
objectFile::READ_NEVER,
objectFile::WRITE_NEVER
),
*this,
0u
),
contactForceWall_
(
objectFile
(
"contactForcWall",
"",
objectFile::READ_IF_PRESENT,
objectFile::WRITE_NEVER
),
*this,
zero3
),
normalStressWall_
(
objectFile
(
"normalStressWall",
"",
objectFile::READ_IF_PRESENT,
objectFile::WRITE_NEVER
),
*this,
zero3
),
shearStressWall_
(
objectFile
(
"shearStressWall",
"",
objectFile::READ_IF_PRESENT,
objectFile::WRITE_NEVER
),
*this,
zero3
)
{
motionComponentName_.assign(motionCompName);
materialName_.assign(materialName);
if( this->numSurfaces() != motionComponentName_.size() )
{
fatalErrorInFunction<<
"Number of surfaces ("<< 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,
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()
{
@ -296,48 +539,4 @@ bool pFlow::geometry::afterIteration()
Kokkos::fence();
return true;
}
pFlow::uniquePtr<pFlow::geometry>
pFlow::geometry::create(
systemControl& control,
const property& prop,
const dictionary& dict,
const multiTriSurface& triSurface,
const wordVector& motionCompName,
const wordVector& propName)
{
word model = dict.getVal<word>("motionModel");
auto geomModel = angleBracketsNames("geometry", model);
REPORT(1)<< "Selecting geometry model . . ."<<endREPORT;
if( dictionaryvCtorSelector_.search(geomModel) )
{
auto objPtr = dictionaryvCtorSelector_[geomModel]
(
control,
prop,
dict,
triSurface,
motionCompName,
propName
);
REPORT(2)<<"Model "<< greenText(geomModel)<<" is created.\n"<<endREPORT;
return objPtr;
}
else
{
printKeys
(
fatalError << "Ctor Selector "<< yellowText(geomModel) << " dose not exist. \n"
<<"Avaiable ones are: \n\n"
,
dictionaryvCtorSelector_
);
fatalExit;
}
return nullptr;
}
}*/

View File

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

View File

@ -21,41 +21,116 @@ Licence:
template<typename MotionModel>
bool pFlow::geometryMotion<MotionModel>::findMotionIndex()
{
motionIndex_.clear();
triMotionIndex_.reserve( this->surface().capacity() );
triMotionIndex_.clear();
ForAll( surfI, motionComponentName_)
if(motionComponentName().size() != numSurfaces() )
{
auto mName = motionComponentName_[surfI];
auto mInd = motionModel_.nameToIndex(mName);
motionIndex_.push_back(mInd);
// fill motionIndex for triangles of the surface
int32 surfSize = this->surface().surfNumTriangles(surfI);
for(int32 i=0; i<surfSize; i++)
fatalErrorInFunction<<
"size of motion component names in the triSurface is not"<<
" equal to size of number of sub-surfaces"<<endl;
return false;
}
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());
pointMotionIndex_.clear();
ForAll(surfI, motionIndex_)
{
auto nP = triSurface_.surfNumPoints(surfI);
for(int32 i=0; i<nP; i++)
{
pointMotionIndex_.push_back(motionIndex_[surfI]);
}
}
pointMotionIndex_.syncViews();
surfMotionIndex_.assign(surfMotionIndex);
triMotionIndex_.assign(triMotionIndex);
pointMotionIndex_.assign(pointMotionIndex);
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>
pFlow::geometryMotion<MotionModel>::geometryMotion
@ -65,144 +140,75 @@ pFlow::geometryMotion<MotionModel>::geometryMotion
)
:
geometry(control, prop),
motionModel_(
this->owner().template emplaceObject<MotionModel>(
objectFile(
motionModelFile__,
"",
objectFile::READ_ALWAYS,
objectFile::WRITE_ALWAYS
)
)
motionModel_
(
objectFile
(
motionModelFile__,
"",
objectFile::READ_ALWAYS,
objectFile::WRITE_ALWAYS
),
owner()
),
moveGeomTimer_("move geometry", &this->timers())
{
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() )
if(!findMotionIndex())
{
moveGeomTimer_.start();
moveGeometry();
moveGeomTimer_.end();
fatalExit;
}
return true;
}
template<typename MotionModel>
bool pFlow::geometryMotion<MotionModel>::afterIteration()
template <typename MotionModelType>
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();
return true;
if(!findMotionIndex())
{
fatalExit;
}
}
template<typename MotionModel>
bool pFlow::geometryMotion<MotionModel>::moveGeometry()
{
real dt = this->dt();
real t = this->currentTime();
auto pointMIndex= pointMotionIndex_.deviceVector();
auto mModel = motionModel_.getModel(t);
realx3* points = triSurface_.pointsData_D();
auto numPoints = triSurface_.numPoints();
Kokkos::parallel_for(
"geometryMotion<MotionModel>::movePoints",
numPoints,
LAMBDA_HD(int32 i){
auto newPos = mModel.transferPoint(pointMIndex[i], points[i], dt);
points[i] = newPos;
});
Kokkos::fence();
// move the motion components
motionModel_.move(t,dt);
// end of calculations
moveGeomTimer_.end();
return true;
}
template<typename MotionModel>
bool pFlow::geometryMotion<MotionModel>::iterate()
{
if( motionModel_.isMoving() )
{
moveGeomTimer_.start();
moveGeometry();
this->calculateNormals();
moveGeomTimer_.end();
}
return true;
}

View File

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

View File

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

View File

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

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,87 +19,191 @@ Licence:
-----------------------------------------------------------------------------*/
#include "AdamsBashforth2.hpp"
#include "pointStructure.hpp"
#include "Time.hpp"
#include "vocabs.hpp"
//const real AB2_coef[] = { 3.0 / 2.0, 1.0 / 2.0};
pFlow::AdamsBashforth2::AdamsBashforth2
(
const word& baseName,
repository& owner,
const pointStructure& pStruct,
const word& method
)
:
integration(baseName, owner, pStruct, method),
dy1_(
owner.emplaceObject<realx3PointField_D>(
objectFile(
groupNames(baseName,"dy1"),
"",
objectFile::READ_IF_PRESENT,
objectFile::WRITE_ALWAYS),
pStruct,
zero3))
namespace pFlow
{
}
/// Range policy for integration kernel (alias)
using rpIntegration = Kokkos::RangePolicy<
DefaultExecutionSpace,
Kokkos::Schedule<Kokkos::Static>,
Kokkos::IndexType<uint32>
>;
bool pFlow::AdamsBashforth2::predict
(
real UNUSED(dt),
realx3Vector_D& UNUSED(y),
realx3Vector_D& UNUSED(dy)
)
{
return true;
}
bool pFlow::AdamsBashforth2::correct
(
bool intAllActive(
real dt,
realx3Vector_D& y,
realx3Vector_D& dy
)
{
if(this->pStruct().allActive())
{
return intAll(dt, y, dy, this->pStruct().activeRange());
}
else
{
return intRange(dt, y, dy, this->pStruct().activePointsMaskD());
}
return true;
}
bool pFlow::AdamsBashforth2::setInitialVals(
const int32IndexContainer& newIndices,
const realx3Vector& y)
{
return true;
}
bool pFlow::AdamsBashforth2::intAll(
real dt,
realx3Vector_D& y,
realx3Vector_D& dy,
range activeRng)
realx3Field_D& y,
realx3PointField_D& dy,
realx3PointField_D& dy1,
real damping = 1.0)
{
auto d_dy = dy.deviceVectorAll();
auto d_y = y.deviceVectorAll();
auto d_dy1= dy1_.deviceVectorAll();
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.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]);
rpIntegration (activeRng.start(), activeRng.end()),
LAMBDA_HD(uint32 i){
d_y[i] += damping * dt*(static_cast<real>(1.5) * d_dy[i] - static_cast<real>(0.5) * d_dy1[i]);
d_dy1[i] = d_dy[i];
});
Kokkos::fence();
return true;
}
bool intScattered
(
real dt,
realx3Field_D& y,
realx3PointField_D& dy,
realx3PointField_D& dy1,
real damping = 1.0
)
{
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] += damping * 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
(
const word& baseName,
pointStructure& pStruct,
const word& method,
const realx3Field_D& initialValField
)
:
integration(baseName, pStruct, method, initialValField),
realx3PointField_D
(
objectFile
(
groupNames(baseName,"dy1"),
pStruct.time().integrationFolder(),
objectFile::READ_IF_PRESENT,
objectFile::WRITE_ALWAYS
),
pStruct,
zero3,
zero3
),
initialValField_(initialValField),
boundaryList_(pStruct, method, *this)
{
realx3PointField_D::addEvent(message::ITEMS_INSERT);
}
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
(
real dt,
realx3Field_D &y,
realx3PointField_D &dy
)
{
return true;
}
bool pFlow::AdamsBashforth2::correct
(
real dt,
realx3PointField_D& y,
realx3PointField_D& dy,
real damping
)
{
auto& dy1l = dy1();
bool success = false;
if(dy1l.isAllActive())
{
success = intAllActive(dt, y.field(), dy, dy1(), damping);
}
else
{
success = intScattered(dt, y.field(), dy, dy1(), damping);
}
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, dy1());
}
else
{
success = intScattered(dt, pStruct.pointPosition(), vel, dy1());
}
success = success && boundaryList_.correctPStruct(dt, pStruct, vel);
return success;
}
/*bool pFlow::AdamsBashforth2::hearChanges
(
const timeInfo &ti,
const message &msg,
const anyList &varList
)
{
if(msg.equivalentTo(message::ITEMS_INSERT))
{
return insertValues(varList, initialValField_.deviceViewAll(), dy1());
}
else
{
return realx3PointField_D::hearChanges(ti, msg, varList);
}
}*/

View File

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

View File

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

View File

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

View File

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

View File

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

View File

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

View File

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

View File

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

View File

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

View File

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

View File

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

View File

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

View File

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

View File

@ -1,15 +1,23 @@
list(APPEND SourceFiles
integration/integration.cpp
AdamsBashforth5/AdamsBashforth5.cpp
AdamsBashforth4/AdamsBashforth4.cpp
AdamsBashforth3/AdamsBashforth3.cpp
boundaries/boundaryIntegration.cpp
boundaries/boundaryIntegrationList.cpp
AdamsBashforth2/AdamsBashforth2.cpp
AdamsMoulton3/AdamsMoulton3.cpp
AdamsMoulton4/AdamsMoulton4.cpp
AdamsMoulton5/AdamsMoulton5.cpp
AdamsBashforth3/AdamsBashforth3.cpp
AdamsBashforth4/AdamsBashforth4.cpp
AdamsBashforth5/AdamsBashforth5.cpp
#AdamsMoulton3/AdamsMoulton3.cpp
#AdamsMoulton4/AdamsMoulton4.cpp
#AdamsMoulton5/AdamsMoulton5.cpp
)
if(pFlow_Build_MPI)
list(APPEND SourceFiles
AdamsBashforth2/processorAB2BoundaryIntegration.cpp)
endif()
set(link_libs Kokkos::kokkos phasicFlow)
pFlow_add_library_install(Integration SourceFiles link_libs)

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,94 @@
#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(
const timeInfo& ti,
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;
}
bool isActive()const override
{
return false;
}
static
uniquePtr<boundaryIntegration> create(
const boundaryBase& boundary,
const pointStructure& pStruct,
const word& method,
integration& intgrtn);
};
}
#endif

View File

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

View File

@ -0,0 +1,51 @@
#ifndef __boundaryIntegrationList_hpp__
#define __boundaryIntegrationList_hpp__
#include "boundaryList.hpp"
#include "boundaryListPtr.hpp"
#include "boundaryIntegration.hpp"
namespace pFlow
{
class integration;
class boundaryIntegrationList
:
public boundaryListPtr<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,57 @@ Licence:
-----------------------------------------------------------------------------*/
#include "integration.hpp"
#include "pointStructure.hpp"
#include "repository.hpp"
pFlow::integration::integration
bool pFlow::integration::insertValues
(
const word& baseName,
repository& owner,
const pointStructure& pStruct,
const word& method
const anyList& varList,
const deviceViewType1D<realx3>& srcVals,
realx3PointField_D& dstFeild
)
:
owner_(owner),
baseName_(baseName),
pStruct_(pStruct)
{
CONSUME_PARAM(method);
const word eventName = message::eventName(message::ITEMS_INSERT);
if(!varList.checkObjectType<uint32IndexContainer>(eventName))
{
fatalErrorInFunction<<"Insertion failed in integration for "<< baseName_<<
", variable "<< eventName <<
" does not exist or the type is incorrect"<<endl;
return false;
}
const auto& indices = varList.getObject<uint32IndexContainer>(
eventName);
dstFeild.field().insertSetElement(indices, srcVals);
return true;
}
pFlow::integration::integration(
const word &baseName,
pointStructure &pStruct,
const word &,
const realx3Field_D &)
: owner_(*pStruct.owner()),
pStruct_(pStruct),
baseName_(baseName)
{}
pFlow::uniquePtr<pFlow::integration>
pFlow::integration::create(
pFlow::integration::create
(
const word& baseName,
repository& owner,
const pointStructure& pStruct,
const word& method)
pointStructure& pStruct,
const word& method,
const realx3Field_D& initialValField
)
{
if( wordvCtorSelector_.search(method) )
{
return wordvCtorSelector_[method] (baseName, owner, pStruct, method);
return wordvCtorSelector_[method] (baseName, pStruct, method, initialValField);
}
else
{

View File

@ -23,14 +23,16 @@ Licence:
#include "virtualConstructor.hpp"
#include "Vectors.hpp"
#include "pointStructure.hpp"
#include "repository.hpp"
#include "pointFields.hpp"
namespace pFlow
{
// - Forward
class repository;
class pointStructure;
/**
* Base class for integrating the first order ODE (IVP)
*
@ -48,18 +50,25 @@ namespace pFlow
*/
class integration
{
protected:
private:
// - Protected data members
/// The owner repository that all fields are storred in
repository& owner_;
/// A reference to pointStructure
const pointStructure& pStruct_;
/// The base name for integration
const word baseName_;
/// A reference to pointStructure
const pointStructure& pStruct_;
protected:
bool insertValues(
const anyList& varList,
const deviceViewType1D<realx3>& srcVals,
realx3PointField_D& dstFeild);
public:
@ -72,9 +81,9 @@ public:
/// Construct from components
integration(
const word& baseName,
repository& owner,
const pointStructure& pStruct,
const word& method);
pointStructure& pStruct,
const word& method,
const realx3Field_D& initialValField);
/// Copy constructor
integration(const integration&) = default;
@ -88,22 +97,22 @@ public:
/// Move assignment
integration& operator = (integration&&) = default;
/// Polymorphic copy/cloning
virtual
uniquePtr<integration> clone()const=0;
/// Destructor
virtual ~integration()=default;
/// Add a virtual constructor
create_vCtor(
create_vCtor
(
integration,
word,
(const word& baseName,
repository& owner,
const pointStructure& pStruct,
const word& method),
(baseName, owner, pStruct, method) );
(
const word& baseName,
pointStructure& pStruct,
const word& method,
const realx3Field_D& initialValField
),
(baseName, pStruct, method, initialValField)
);
// - Methods
@ -129,32 +138,33 @@ public:
return owner_;
}
virtual
void updateBoundariesSlaveToMasterIfRequested() = 0;
/// return integration method
virtual
word method()const = 0 ;
/// Prediction step in integration
virtual
bool predict(real dt, realx3Vector_D& y, realx3Vector_D& dy) = 0;
bool predict(real dt, realx3PointField_D& y, realx3PointField_D& dy) = 0;
virtual
bool predict(real dt, realx3Field_D& y, realx3PointField_D& dy) = 0;
/// Correction/main integration step
virtual
bool correct(real dt, realx3Vector_D& y, realx3Vector_D& dy) = 0;
bool correct(real dt, realx3PointField_D& y, realx3PointField_D& dy, real damping = 1.0) = 0;
/// Set the initial values for new indices
virtual
bool setInitialVals(
const int32IndexContainer& newIndices,
const realx3Vector& y) = 0;
/// Check if the method requires any set initial vals
virtual
bool needSetInitialVals()const = 0;
bool correctPStruct(real dt, pointStructure& pStruct, realx3PointField_D& vel) = 0;
/// Create the polymorphic object based on inputs
static
uniquePtr<integration> create(
const word& baseName,
repository& owner,
const pointStructure& pStruct,
const word& method);
pointStructure& pStruct,
const word& method,
const realx3Field_D& initialValField);
}; // integration

View File

@ -1,11 +1,43 @@
set(SourceFiles
contactSearch/methods/cellBased/NBS/mapperNBS.cpp
contactSearch/methods/cellBased/NBS/mapperNBSKernels.cpp
contactSearch/methods/cellBased/NBS/NBSLevel0.cpp
contactSearch/methods/cellBased/NBS/NBS.cpp
contactSearch/methods/cellBased/NBS/cellsWallLevel0.cpp
contactSearch/boundaries/boundaryContactSearch/boundaryContactSearch.cpp
contactSearch/boundaries/periodicBoundaryContactSearch/ppwBndryContactSearchKernels.cpp
contactSearch/boundaries/periodicBoundaryContactSearch/ppwBndryContactSearch.cpp
contactSearch/boundaries/periodicBoundaryContactSearch/wallBoundaryContactSearch.cpp
contactSearch/boundaries/periodicBoundaryContactSearch/periodicBoundaryContactSearch.cpp
contactSearch/boundaries/boundaryContactSearchList.cpp
contactSearch/contactSearch/contactSearch.cpp
contactSearch/ContactSearch/ContactSearchs.cpp
interaction/interaction.cpp
sphereInteraction/sphereInteractions.cpp
sphereInteraction/sphereInteractionsLinearModels.cpp
sphereInteraction/sphereInteractionsNonLinearModels.cpp
sphereInteraction/sphereInteractionsNonLinearModModels.cpp
grainInteraction/grainInteractionsLinearModels.cpp
grainInteraction/grainInteractionsNonLinearModels.cpp
grainInteraction/grainInteractionsNonLinearModModels.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)
pFlow_add_library_install(Interaction SourceFiles link_libs)

View File

@ -0,0 +1,339 @@
/*------------------------------- 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 __cGAbsoluteLinearCF_hpp__
#define __cGAbsoluteLinearCF_hpp__
#include "types.hpp"
#include "symArrays.hpp"
namespace pFlow::cfModels
{
template<bool limited=true>
class cGAbsoluteLinear
{
public:
struct contactForceStorage
{
realx3 overlap_t_ = 0.0;
};
struct linearProperties
{
real kn_ = 1000.0;
real kt_ = 800.0;
real en_ = 0.0;
real ethat_ = 0.0;
real mu_ = 0.00001;
INLINE_FUNCTION_HD
linearProperties(){}
INLINE_FUNCTION_HD
linearProperties(real kn, real kt, real en, real etha_t, real mu ):
kn_(kn), kt_(kt), en_(en),ethat_(etha_t), mu_(mu)
{}
INLINE_FUNCTION_HD
linearProperties(const linearProperties&)=default;
INLINE_FUNCTION_HD
linearProperties& operator=(const linearProperties&)=default;
INLINE_FUNCTION_HD
~linearProperties() = default;
};
protected:
using LinearArrayType = symArray<linearProperties>;
int32 numMaterial_ = 0;
ViewType1D<real> rho_;
LinearArrayType linearProperties_;
int32 addDissipationModel_;
bool readLinearDictionary(const dictionary& dict)
{
auto kn = dict.getVal<realVector>("kn");
auto kt = dict.getVal<realVector>("kt");
auto en = dict.getVal<realVector>("en");
auto et = dict.getVal<realVector>("et");
auto mu = dict.getVal<realVector>("mu");
auto nElem = kn.size();
if(nElem != kt.size())
{
fatalErrorInFunction<<
"sizes of kn("<<nElem<<") and kt("<<kt.size()<<") do not match.\n";
return false;
}
if(nElem != kt.size())
{
fatalErrorInFunction<<
"sizes of kn("<<nElem<<") and kt("<<kt.size()<<") do not match.\n";
return false;
}
if(nElem != en.size())
{
fatalErrorInFunction<<
"sizes of kn("<<nElem<<") and en("<<en.size()<<") do not match.\n";
return false;
}
if(nElem != et.size())
{
fatalErrorInFunction<<
"sizes of kn("<<nElem<<") and et("<<et.size()<<") do not match.\n";
return false;
}
if(nElem != mu.size())
{
fatalErrorInFunction<<
"sizes of kn("<<nElem<<") and mu("<<mu.size()<<") do not match.\n";
return false;
}
// check if size of vector matchs a symetric array
uint32 nMat;
if( !LinearArrayType::getN(nElem, nMat) )
{
fatalErrorInFunction<<
"sizes of properties do not match a symetric array with size ("<<
numMaterial_<<"x"<<numMaterial_<<").\n";
return false;
}
else if( numMaterial_ != nMat)
{
fatalErrorInFunction<<
"size mismatch for porperties. \n"<<
"you supplied "<< numMaterial_<<" items in materials list and "<<
nMat << " for other properties.\n";
return false;
}
realVector etha_t("etha_t", nElem);
ForAll(i , kn)
{
etha_t[i] = -2.0*log( et[i])*sqrt(kt[i]*2/7) /
sqrt(log(pow(et[i],2))+ pow(Pi,2));
}
Vector<linearProperties> prop("prop", nElem);
ForAll(i,kn)
{
prop[i] = {kn[i], kt[i], en[i], etha_t[i], mu[i] };
}
linearProperties_.assign(prop);
auto adm = dict.getVal<word>("additionalDissipationModel");
if(adm == "none")
{
addDissipationModel_ = 1;
}
else if(adm == "LU")
{
addDissipationModel_ = 2;
}
else if (adm == "GB")
{
addDissipationModel_ = 3;
}
else
{
addDissipationModel_ = 1;
}
return true;
}
static const char* modelName()
{
if constexpr (limited)
{
return "cGAbsoluteLinearLimited";
}
else
{
return "cGAbsoluteLinearNonLimited";
}
return "";
}
public:
TypeInfoNV(modelName());
INLINE_FUNCTION_HD
cGAbsoluteLinear(){}
cGAbsoluteLinear(int32 nMaterial, const ViewType1D<real>& rho, const dictionary& dict)
:
numMaterial_(nMaterial),
rho_("rho",nMaterial),
linearProperties_("linearProperties",nMaterial)
{
Kokkos::deep_copy(rho_,rho);
if(!readLinearDictionary(dict))
{
fatalExit;
}
}
INLINE_FUNCTION_HD
cGAbsoluteLinear(const cGAbsoluteLinear&) = default;
INLINE_FUNCTION_HD
cGAbsoluteLinear(cGAbsoluteLinear&&) = default;
INLINE_FUNCTION_HD
cGAbsoluteLinear& operator=(const cGAbsoluteLinear&) = default;
INLINE_FUNCTION_HD
cGAbsoluteLinear& operator=(cGAbsoluteLinear&&) = default;
INLINE_FUNCTION_HD
~cGAbsoluteLinear()=default;
INLINE_FUNCTION_HD
int32 numMaterial()const
{
return numMaterial_;
}
//// - Methods
INLINE_FUNCTION_HD
void contactForce
(
const real dt,
const uint32 i,
const uint32 j,
const uint32 propId_i,
const uint32 propId_j,
const real Ri,
const real Rj,
const real cGFi,
const real cGFj,
const real ovrlp_n,
const realx3& Vr,
const realx3& Nij,
contactForceStorage& history,
realx3& FCn,
realx3& FCt
)const
{
auto prop = linearProperties_(propId_i,propId_j);
real f_ = ( cGFi + cGFj )/2 ;
real vrn = dot(Vr, Nij);
realx3 Vt = Vr - vrn*Nij;
history.overlap_t_ += Vt*dt;
real mi = 3*Pi/4*pow(Ri,3)*rho_[propId_i];
real mj = 3*Pi/4*pow(Rj,3)*rho_[propId_j];
real sqrt_meff = sqrt((mi*mj)/(mi+mj));
// disipation model
if (addDissipationModel_==2)
{
prop.en_ = sqrt(1+((pow(prop.en_,2)-1)*f_));
}
else if (addDissipationModel_==3)
{
auto pie =3.14;
prop.en_ = exp((pow(f_,static_cast<real>(1.5))*log(prop.en_)*sqrt( (1-((pow(log(prop.en_),2))/(pow(log(prop.en_),2)+pow(pie,2))))/(1-(pow(f_,3)*(pow(log(prop.en_),2))/(pow(log(prop.en_),2)+pow(pie,2)))) ) ));
}
real ethan_ = -2.0*log(prop.en_)*sqrt(prop.kn_)/
sqrt(pow(log(prop.en_),2)+ pow(Pi,2));
//REPORT(0)<<"\n en n is : "<<END_REPORT;
//REPORT(0)<< prop.en_ <<END_REPORT;
FCn = ( -pow(f_,3)*prop.kn_ * ovrlp_n - sqrt_meff * pow(f_,static_cast<real>(1.5)) * ethan_ * vrn)*Nij;
FCt = ( -pow(f_,3)*prop.kt_ * history.overlap_t_ - sqrt_meff * pow(f_,static_cast<real>(1.5)) * prop.ethat_*Vt);
real ft = length(FCt);
real ft_fric = prop.mu_ * length(FCn);
if(ft > ft_fric)
{
if( length(history.overlap_t_) >zero)
{
if constexpr (limited)
{
FCt *= (ft_fric/ft);
history.overlap_t_ = - (FCt/prop.kt_);
}
else
{
FCt = (FCt/ft)*ft_fric;
}
//cout<<"friction is applied here \n";
}
else
{
FCt = 0.0;
}
}
}
};
} //pFlow::cfModels
#endif

View File

@ -0,0 +1,307 @@
/*------------------------------- phasicFlow ---------------------------------
O C enter of
O O E ngineering and
O O M ultiscale modeling of
OOOOOOO F luid flow
------------------------------------------------------------------------------
Copyright (C): www.cemf.ir
email: hamid.r.norouzi AT gmail.com
------------------------------------------------------------------------------
Licence:
This file is part of phasicFlow code. It is a free software for simulating
granular and multiphase flows. You can redistribute it and/or modify it under
the terms of GNU General Public License v3 or any other later versions.
phasicFlow is distributed to help others in their research in the field of
granular and multiphase flows, but WITHOUT ANY WARRANTY; without even the
implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
-----------------------------------------------------------------------------*/
#ifndef __cGNonLinearCF_hpp__
#define __cGNonLinearCF_hpp__
#include "types.hpp"
namespace pFlow::cfModels
{
template<bool limited=true>
class cGNonLinear
{
public:
struct contactForceStorage
{
realx3 overlap_t_ = 0.0;
};
struct cGNonLinearProperties
{
real Yeff_ = 1000000.0;
real Geff_ = 8000000.0;
real en_ = 1.0;
real mu_ = 0.00001;
INLINE_FUNCTION_HD
cGNonLinearProperties(){}
INLINE_FUNCTION_HD
cGNonLinearProperties(real Yeff, real Geff, real en, real mu ):
Yeff_(Yeff), Geff_(Geff), en_(en), mu_(mu)
{}
INLINE_FUNCTION_HD
cGNonLinearProperties(const cGNonLinearProperties&)=default;
INLINE_FUNCTION_HD
cGNonLinearProperties& operator=(const cGNonLinearProperties&)=default;
INLINE_FUNCTION_HD
~cGNonLinearProperties() = default;
};
protected:
using cGNonLinearArrayType = symArray<cGNonLinearProperties>;
int32 numMaterial_ = 0;
ViewType1D<real> rho_;
int32 addDissipationModel_;
cGNonLinearArrayType nonlinearProperties_;
bool readNonLinearDictionary(const dictionary& dict)
{
auto Yeff = dict.getVal<realVector>("Yeff");
auto Geff = dict.getVal<realVector>("Geff");
auto nu = dict.getVal<realVector>("nu");
auto en = dict.getVal<realVector>("en");
auto mu = dict.getVal<realVector>("mu");
auto nElem = Yeff.size();
if(nElem != nu.size())
{
fatalErrorInFunction<<
"sizes of Yeff("<<nElem<<") and nu("<<nu.size()<<") do not match.\n";
return false;
}
if(nElem != en.size())
{
fatalErrorInFunction<<
"sizes of Yeff("<<nElem<<") and en("<<en.size()<<") do not match.\n";
return false;
}
if(nElem != mu.size())
{
fatalErrorInFunction<<
"sizes of Yeff("<<nElem<<") and mu("<<mu.size()<<") do not match.\n";
return false;
}
// check if size of vector matchs a symetric array
uint32 nMat;
if( !cGNonLinearArrayType::getN(nElem, nMat) )
{
fatalErrorInFunction<<
"sizes of properties do not match a symetric array with size ("<<
numMaterial_<<"x"<<numMaterial_<<").\n";
return false;
}
else if( numMaterial_ != nMat)
{
fatalErrorInFunction<<
"size mismatch for porperties. \n";
return false;
}
Vector<cGNonLinearProperties> prop("prop",nElem);
ForAll(i,Yeff)
{
prop[i] = {Yeff[i], Geff[i], en[i], mu[i]};
}
nonlinearProperties_.assign(prop);
auto adm = dict.getVal<word>("additionalDissipationModel");
if(adm == "none")
{
addDissipationModel_ = 1;
}
else if(adm == "LU")
{
addDissipationModel_ = 2;
}
else if (adm == "GB")
{
addDissipationModel_ = 3;
}
else
{
addDissipationModel_ = 1;
}
return true;
}
static const char* modelName()
{
if constexpr (limited)
{
return "cGNonLinearLimited";
}
else
{
return "cGNonLinearNonLimited";
}
return "";
}
public:
TypeInfoNV(modelName());
INLINE_FUNCTION_HD
cGNonLinear(){}
cGNonLinear(
int32 nMaterial,
const ViewType1D<real>& rho,
const dictionary& dict)
:
numMaterial_(nMaterial),
rho_("rho",nMaterial),
nonlinearProperties_("cGNonLinearProperties",nMaterial)
{
Kokkos::deep_copy(rho_,rho);
if(!readNonLinearDictionary(dict))
{
fatalExit;
}
}
INLINE_FUNCTION_HD
cGNonLinear(const cGNonLinear&) = default;
INLINE_FUNCTION_HD
cGNonLinear(cGNonLinear&&) = default;
INLINE_FUNCTION_HD
cGNonLinear& operator=(const cGNonLinear&) = default;
INLINE_FUNCTION_HD
cGNonLinear& operator=(cGNonLinear&&) = default;
INLINE_FUNCTION_HD
~cGNonLinear()=default;
INLINE_FUNCTION_HD
int32 numMaterial()const
{
return numMaterial_;
}
//// - Methods
INLINE_FUNCTION_HD
void contactForce
(
const real dt,
const uint32 i,
const uint32 j,
const uint32 propId_i,
const uint32 propId_j,
const real Ri,
const real Rj,
const real cGFi,
const real cGFj,
const real ovrlp_n,
const realx3& Vr,
const realx3& Nij,
contactForceStorage& history,
realx3& FCn,
realx3& FCt
)const
{
const auto prop = nonlinearProperties_(propId_i,propId_j);
const real f = 2/( 1/cGFi + 1/cGFj );
real vrn = dot(Vr, Nij);
realx3 Vt = Vr - vrn*Nij;
history.overlap_t_ += Vt*dt;
real mi = 3*Pi/4*pow(Ri,static_cast<real>(3))*rho_[propId_i];
real mj = 3*Pi/4*pow(Rj,static_cast<real>(3))*rho_[propId_j];
real Reff = 1/(1/Ri + 1/Rj);
real K_hertz = 4.0/3.0*prop.Yeff_*sqrt(Reff);
real sqrt_meff_K_hertz = sqrt((mi*mj)/(mi+mj) * K_hertz);
real en = prop.en_;
if (addDissipationModel_==2)
{
en = sqrt(1+((pow(prop.en_,2)-1)*f));
}
else if (addDissipationModel_==3)
{
en = exp((pow(f,static_cast<real>(1.5))*log(prop.en_)*sqrt( (1-((pow(log(prop.en_),2))/(pow(log(prop.en_),2)+pow(Pi,2))))/(1-(pow(f,3)*(pow(log(prop.en_),2))/(pow(log(prop.en_),2)+pow(Pi,2)))) ) ));
}
real Kn = static_cast<real>(4.0/3.0) * prop.Yeff_ * sqrt(Reff*ovrlp_n);
real ethan_ = -2.0*log(en)*sqrt(Kn)/
sqrt(pow(log(en),2)+ pow(Pi,2));
FCn = ( - Kn*ovrlp_n -
sqrt_meff_K_hertz*ethan_*pow(ovrlp_n,static_cast<real>(0.25))*vrn)*Nij;
FCt = (f*f)*(- static_cast<real>(8.0) * prop.Geff_ * sqrt(Reff*ovrlp_n) ) * history.overlap_t_;
real ft = length(FCt);
real ft_fric = prop.mu_ * length(FCn);
// apply friction
if(ft > ft_fric)
{
if( length(history.overlap_t_) >0.0)
{
if constexpr (limited)
{
real kt = static_cast<real>(8.0) * prop.Geff_ * sqrt(Reff*ovrlp_n);
FCt *= (ft_fric/ft);
history.overlap_t_ = - FCt/(f*f*kt);
}
else
{
FCt = (FCt/ft)*ft_fric;
}
}
else
{
FCt = 0.0;
}
}
}
};
} //pFlow::CFModels
#endif

View File

@ -0,0 +1,326 @@
/*------------------------------- 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 __cGRelativeLinearCF_hpp__
#define __cGRelativeLinearCF_hpp__
#include "types.hpp"
#include "symArrays.hpp"
namespace pFlow::cfModels
{
template<bool limited=true>
class cGRelativeLinear
{
public:
struct contactForceStorage
{
realx3 overlap_t_ = 0.0;
};
struct linearProperties
{
real kn_ = 1000.0;
real kt_ = 800.0;
real en_ = 1.0;
real mu_ = 0.00001;
INLINE_FUNCTION_HD
linearProperties(){}
INLINE_FUNCTION_HD
linearProperties(real kn, real kt, real en, real mu ):
kn_(kn), kt_(kt), en_(en), mu_(mu)
{}
INLINE_FUNCTION_HD
linearProperties(const linearProperties&)=default;
INLINE_FUNCTION_HD
linearProperties& operator=(const linearProperties&)=default;
INLINE_FUNCTION_HD
~linearProperties() = default;
};
protected:
using LinearArrayType = symArray<linearProperties>;
int32 numMaterial_ = 0;
ViewType1D<real> rho_;
LinearArrayType linearProperties_;
int32 addDissipationModel_;
bool readLinearDictionary(const dictionary& dict)
{
auto kn = dict.getVal<realVector>("kn");
auto kt = dict.getVal<realVector>("kt");
auto en = dict.getVal<realVector>("en");
auto mu = dict.getVal<realVector>("mu");
auto nElem = kn.size();
if(nElem != kt.size())
{
fatalErrorInFunction<<
"sizes of kn("<<nElem<<") and kt("<<kt.size()<<") do not match.\n";
return false;
}
if(nElem != en.size())
{
fatalErrorInFunction<<
"sizes of kn("<<nElem<<") and en("<<en.size()<<") do not match.\n";
return false;
}
if(nElem != mu.size())
{
fatalErrorInFunction<<
"sizes of kn("<<nElem<<") and mu("<<mu.size()<<") do not match.\n";
return false;
}
// check if size of vector matchs a symetric array
uint32 nMat;
if( !LinearArrayType::getN(nElem, nMat) )
{
fatalErrorInFunction<<
"sizes of properties do not match a symetric array with size ("<<
numMaterial_<<"x"<<numMaterial_<<").\n";
return false;
}
else if( numMaterial_ != nMat)
{
fatalErrorInFunction<<
"size mismatch for porperties. \n"<<
"you supplied "<< numMaterial_<<" items in materials list and "<<
nMat << " for other properties.\n";
return false;
}
Vector<linearProperties> prop("prop", nElem);
ForAll(i,kn)
{
prop[i] = {kn[i], kt[i], en[i], mu[i] };
}
linearProperties_.assign(prop);
auto adm = dict.getVal<word>("additionalDissipationModel");
if(adm == "none")
{
addDissipationModel_ = 1;
}
else if(adm == "LU")
{
addDissipationModel_ = 2;
}
else if (adm == "GB")
{
addDissipationModel_ = 3;
}
else
{
addDissipationModel_ = 1;
}
return true;
}
static const char* modelName()
{
if constexpr (limited)
{
return "cGRelativeLinearLimited";
}
else
{
return "cGRelativeLinearNonLimited";
}
return "";
}
public:
TypeInfoNV(modelName());
INLINE_FUNCTION_HD
cGRelativeLinear(){}
cGRelativeLinear(int32 nMaterial, const ViewType1D<real>& rho, const dictionary& dict)
:
numMaterial_(nMaterial),
rho_("rho",nMaterial),
linearProperties_("linearProperties",nMaterial)
{
Kokkos::deep_copy(rho_,rho);
if(!readLinearDictionary(dict))
{
fatalExit;
}
}
INLINE_FUNCTION_HD
cGRelativeLinear(const cGRelativeLinear&) = default;
INLINE_FUNCTION_HD
cGRelativeLinear(cGRelativeLinear&&) = default;
INLINE_FUNCTION_HD
cGRelativeLinear& operator=(const cGRelativeLinear&) = default;
INLINE_FUNCTION_HD
cGRelativeLinear& operator=(cGRelativeLinear&&) = default;
INLINE_FUNCTION_HD
~cGRelativeLinear()=default;
INLINE_FUNCTION_HD
int32 numMaterial()const
{
return numMaterial_;
}
//// - Methods
INLINE_FUNCTION_HD
void contactForce
(
const real dt,
const uint32 i,
const uint32 j,
const uint32 propId_i,
const uint32 propId_j,
const real Ri,
const real Rj,
const real cGFi,
const real cGFj,
const real ovrlp_n,
const realx3& Vr,
const realx3& Nij,
contactForceStorage& history,
realx3& FCn,
realx3& FCt
)const
{
auto prop = linearProperties_(propId_i,propId_j);
real f_ = ( cGFi + cGFj )/2 ;
real vrn = dot(Vr, Nij);
realx3 Vt = Vr - vrn*Nij;
history.overlap_t_ += Vt*dt;
real mi = 3*Pi/4*pow(Ri,3)*rho_[propId_i];
real mj = 3*Pi/4*pow(Rj,3)*rho_[propId_j];
real sqrt_meff = sqrt((mi*mj)/(mi+mj));
real en = prop.en_;
if (addDissipationModel_==2)
{
en = sqrt(1+((pow(prop.en_,2)-1)*f_));
}
else if (addDissipationModel_==3)
{
en = exp(
(
pow(f_,static_cast<real>(1.5))*
log(prop.en_)*
sqrt(
(1-((pow(log(prop.en_),static_cast<real>(2.0))
)/
(
pow(log(prop.en_),static_cast<real>(2.0))+
pow(Pi,static_cast<real>(2.0)))))/
(1-(pow(f_,3)*(pow(log(prop.en_),2))/
(pow(log(prop.en_),static_cast<real>(2.0))+
pow(Pi,static_cast<real>(2.0))))) ) ));
}
real ethan_ = -2.0*log(en)*sqrt(prop.kn_)/
sqrt(pow(log(en),2)+ pow(Pi,2));
FCn = ( -f_*prop.kn_ * ovrlp_n - sqrt_meff * pow(f_,half) * ethan_ * vrn)*Nij;
FCt = ( -f_*prop.kt_ * history.overlap_t_);
real ft = length(FCt);
real ft_fric = prop.mu_ * length(FCn);
if(ft > ft_fric)
{
if( length(history.overlap_t_) >static_cast<real>(0.0))
{
if constexpr (limited)
{
FCt *= (ft_fric/ft);
history.overlap_t_ = - (FCt/prop.kt_);
}
else
{
FCt = (FCt/ft)*ft_fric;
}
//cout<<"friction is applied here \n";
}
else
{
FCt = 0.0;
}
}
}
};
} //pFlow::cfModels
#endif

View File

@ -133,18 +133,19 @@ protected:
return false;
}
realVector etha_n(nElem);
realVector etha_t(nElem);
realVector etha_n("etha_n", nElem);
realVector etha_t("etha_t", nElem);
ForAll(i , kn)
{
etha_n[i] = -2.0*log(en[i])*sqrt(kn[i])/
sqrt(pow(log(en[i]),2.0)+ pow(Pi,2.0));
sqrt(pow(log(en[i]),static_cast<real>(2.0))+ pow(Pi,static_cast<real>(2.0)));
etha_t[i] = -2.0*log( et[i]*sqrt(kt[i]) )/
sqrt(pow(log(et[i]),2.0)+ pow(Pi,2.0));
sqrt(pow(log(et[i]),static_cast<real>(2.0))+ pow(Pi,static_cast<real>(2.0)));
}
Vector<linearProperties> prop(nElem);
Vector<linearProperties> prop("prop", nElem);
ForAll(i,kn)
{
prop[i] = {kn[i], kt[i], etha_n[i], etha_t[i], mu[i]};
@ -219,10 +220,10 @@ public:
void contactForce
(
const real dt,
const int32 i,
const int32 j,
const int32 propId_i,
const int32 propId_j,
const uint32 i,
const uint32 j,
const uint32 propId_i,
const uint32 propId_j,
const real Ri,
const real Rj,
const real ovrlp_n,
@ -242,8 +243,8 @@ public:
history.overlap_t_ += Vt*dt;
real mi = 3*Pi/4*pow(Ri,3.0)*rho_[propId_i];
real mj = 3*Pi/4*pow(Rj,3.0)*rho_[propId_j];
real mi = 3*Pi/4*pow(Ri,static_cast<real>(3.0))*rho_[propId_i];
real mj = 3*Pi/4*pow(Rj,static_cast<real>(3.0))*rho_[propId_j];
real sqrt_meff = sqrt((mi*mj)/(mi+mj));

View File

@ -121,7 +121,7 @@ protected:
}
realVector etha_n(nElem);
realVector etha_n("etha_n",nElem);
ForAll(i , en)
{
@ -131,13 +131,13 @@ protected:
// we take out sqrt(meff*K_hertz) here and then consider this term
// when calculating damping part.
etha_n[i] = -2.2664*log(en[i])/
sqrt(pow(log(en[i]),2.0)+ pow(Pi,2.0));
sqrt(pow(log(en[i]),static_cast<real>(2.0))+ pow(Pi,static_cast<real>(2.0)));
// no damping for tangential part
}
Vector<nonLinearProperties> prop(nElem);
Vector<nonLinearProperties> prop("prop",nElem);
ForAll(i,Yeff)
{
prop[i] = {Yeff[i], Geff[i], etha_n[i], mu[i]};
@ -214,10 +214,10 @@ public:
void contactForce
(
const real dt,
const int32 i,
const int32 j,
const int32 propId_i,
const int32 propId_j,
const uint32 i,
const uint32 j,
const uint32 propId_i,
const uint32 propId_j,
const real Ri,
const real Rj,
const real ovrlp_n,
@ -255,7 +255,7 @@ public:
// apply friction
if(ft > ft_fric)
{
if( length(history.overlap_t_) >0.0)
if( length(history.overlap_t_) >zero)
{
if constexpr (limited)
{

View File

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

View File

@ -0,0 +1,50 @@
/*------------------------------- 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 __grainContactForceModels_hpp__
#define __grainContactForceModels_hpp__
#include "cGAbsoluteLinearCF.hpp"
#include "cGRelativeLinearCF.hpp"
#include "cGNonLinearCF.hpp"
#include "grainRolling.hpp"
namespace pFlow::cfModels
{
using limitedCGAbsoluteLinearGrainRolling = grainRolling<cGAbsoluteLinear<true>>;
using nonLimitedCGAbsoluteLinearGrainRolling = grainRolling<cGAbsoluteLinear<false>>;
using limitedCGRelativeLinearGrainRolling = grainRolling<cGRelativeLinear<true>>;
using nonLimitedCGRelativeLinearGrainRolling = grainRolling<cGRelativeLinear<false>>;
using limitedCGNonLinearGrainRolling = grainRolling<cGNonLinear<true>>;
using nonLimitedCGNonLinearGrainRolling = grainRolling<cGNonLinear<false>>;
}
#endif //__grainContactForceModels_hpp__

View File

@ -0,0 +1,119 @@
/*------------------------------- 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 __grainRolling_hpp__
#define __grainRolling_hpp__
namespace pFlow::cfModels
{
template<typename contactForceModel>
class grainRolling
:
public contactForceModel
{
public:
using contactForceStorage =
typename contactForceModel::contactForceStorage;
realSymArray_D mur_;
bool readGrainDict(const dictionary& dict)
{
auto mur = dict.getVal<realVector>("mur");
uint32 nMat;
if(!realSymArray_D::getN(mur.size(),nMat) || nMat != this->numMaterial())
{
fatalErrorInFunction<<
"wrong number of values supplied in mur. \n";
return false;
}
mur_.assign(mur);
return true;
}
public:
TypeInfoNV(word("normal<"+contactForceModel::TYPENAME()+">"));
grainRolling(int32 nMaterial, const ViewType1D<real>& rho, const dictionary& dict)
:
contactForceModel(nMaterial, rho, dict),
mur_("mur", nMaterial)
{
if(!readGrainDict(dict))
{
fatalExit;
}
}
INLINE_FUNCTION_HD
void rollingFriction
(
const real dt,
const uint32 i,
const uint32 j,
const uint32 propId_i,
const uint32 propId_j,
const real Ri,
const real Rj,
const real cGFi,
const real cGFj,
const realx3& wi,
const realx3& wj,
const realx3& Nij,
const realx3& FCn,
realx3& Mri,
realx3& Mrj
)const
{
realx3 w_hat = wi-wj;
real w_hat_mag = length(w_hat);
if( !equal(w_hat_mag,static_cast<real>(0.0)) )
w_hat /= w_hat_mag;
else
w_hat = static_cast<real>(0.0);
auto Reff = (Ri*Rj)/(Ri+Rj);
Mri = ( -mur_(propId_i,propId_j) *length(FCn) * Reff ) * w_hat ;
//removing the normal part
// Mri = Mri - ( (Mri .dot. nij)*nij )
Mrj = -Mri;
}
};
}
#endif

View File

@ -76,10 +76,10 @@ public:
void rollingFriction
(
const real dt,
const int32 i,
const int32 j,
const int32 propId_i,
const int32 propId_j,
const uint32 i,
const uint32 j,
const uint32 propId_i,
const uint32 propId_j,
const real Ri,
const real Rj,
const realx3& wi,
@ -94,10 +94,10 @@ public:
realx3 w_hat = wi-wj;
real w_hat_mag = length(w_hat);
if( !equal(w_hat_mag,0.0) )
if( !equal(w_hat_mag,static_cast<real>(0.0)) )
w_hat /= w_hat_mag;
else
w_hat = 0.0;
w_hat = static_cast<real>(0.0);
auto Reff = (Ri*Rj)/(Ri+Rj);

View File

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

View File

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

View File

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

View File

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

View File

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

View File

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

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

View File

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

View File

@ -0,0 +1,63 @@
/*------------------------------- phasicFlow ---------------------------------
O C enter of
O O E ngineering and
O O M ultiscale modeling of
OOOOOOO F luid flow
------------------------------------------------------------------------------
Copyright (C): www.cemf.ir
email: hamid.r.norouzi AT gmail.com
------------------------------------------------------------------------------
Licence:
This file is part of phasicFlow code. It is a free software for simulating
granular and multiphase flows. You can redistribute it and/or modify it under
the terms of GNU General Public License v3 or any other later versions.
phasicFlow is distributed to help others in their research in the field of
granular and multiphase flows, but WITHOUT ANY WARRANTY; without even the
implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
-----------------------------------------------------------------------------*/
#include "boundaryListPtr.hpp"
#include "boundaryContactSearch.hpp"
namespace pFlow
{
class boundaryList;
class contactSearch;
class boundaryContactSearchList
:
public boundaryListPtr<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,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 __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;
bool isActive()const override
{
return true;
}
};
}
#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, half ) ),
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,126 @@ Licence:
-----------------------------------------------------------------------------*/
#include "contactSearch.hpp"
#include "streams.hpp"
#include "particles.hpp"
pFlow::contactSearch::contactSearch(
const dictionary& dict,
const box& domain,
const box& extDomain,
const particles& prtcl,
const geometry& geom,
Timers& timers)
:
interactionBase(prtcl, geom),
domain_(domain),
dict_(dict),
sphereSphereTimer_("particle-particle contact search", &timers),
sphereWallTimer_("particle-wall contact search", &timers)
extendedDomainBox_(extDomain),
updateInterval_(dict.getValMax<uint32>("updateInterval", 1u)),
particles_(prtcl),
geometry_(geom),
bTimer_("Boundary particles contact search", &timers),
ppTimer_("Internal particles contact search", &timers)
{
}
const pFlow::pointStructure &pFlow::contactSearch::pStruct() const
{
return particles_.pStruct();
}
pFlow::uniquePtr<pFlow::contactSearch> pFlow::contactSearch::create(
const dictionary& dict,
const box& domain,
const particles& prtcl,
const geometry& geom,
Timers& timers)
bool pFlow::contactSearch::broadSearch
(
const timeInfo &ti,
csPairContainerType &ppPairs,
csPairContainerType &pwPairs,
bool force
)
{
if(enterBroadSearch(ti, force))
{
ppTimer_.start();
if( !BroadSearch(
ti,
ppPairs,
pwPairs,
force ) )
{
fatalErrorInFunction;
performedSearch_ = false;
return false;
}
ppTimer_.end();
performedSearch_ = true;
lastUpdated_ = ti.currentIter();
}
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++)
{
//output<<" boundarySearch "<< i <<" for iter "<< ti.iter()<<endl;
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 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))
{
auto objPtr = dictionaryvCtorSelector_[model] (dict, domain, prtcl, geom, timers);
REPORT(2)<<"Model "<< greenText(model)<<" is created."<<endREPORT;
auto objPtr = dictionaryvCtorSelector_[model] (dict, extDomain, prtcl, geom, timers);
return objPtr;
}
else
{
printKeys
(
fatalError << "Ctor Selector "<< model << " dose not exist. \n"
fatalError << "Ctor Selector "<< model << " does not exist. \n"
<<"Avaiable ones are: \n\n"
,
dictionaryvCtorSelector_

View File

@ -23,42 +23,66 @@ Licence:
#define __contactSearch_hpp__
#include "interactionBase.hpp"
#include "unsortedPairs.hpp"
#include "box.hpp"
#include "contactSearchGlobals.hpp"
#include "dictionary.hpp"
#include "virtualConstructor.hpp"
#include "timeInfo.hpp"
#include "Timer.hpp"
namespace pFlow
{
// - forward
class box;
class particles;
class geometry;
class pointStructure;
class contactSearch
:
public interactionBase
{
public:
using IdType = typename interactionBase::IdType;
private:
using IndexType = typename interactionBase::IndexType;
const box& extendedDomainBox_;
using ExecutionSpace = typename interactionBase::ExecutionSpace;
/// @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;
dictionary dict_;
/// const ref to particles
const particles& particles_;
Timer sphereSphereTimer_;
/// const ref to geometry
const geometry& geometry_;
Timer sphereWallTimer_;
Timer bTimer_;
auto& dict()
{
return dict_;
}
Timer ppTimer_;
virtual
bool BroadSearch(
const timeInfo& ti,
csPairContainerType& ppPairs,
csPairContainerType& pwPairs,
bool force
)=0;
virtual
bool BoundaryBroadSearch(
uint32 bndryIndex,
const timeInfo& ti,
csPairContainerType& ppPairs,
csPairContainerType& pwPairs,
bool force = false
)=0;
public:
@ -66,7 +90,7 @@ public:
contactSearch(
const dictionary& dict,
const box& domain,
const box& extDomain,
const particles& prtcl,
const geometry& geom,
Timers& timers);
@ -88,40 +112,97 @@ public:
(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 dict_;
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 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_;
}
bool broadSearch(
PairContainerType& ppPairs,
PairContainerType& pwPairs,
bool force = false) = 0;
const timeInfo& ti,
csPairContainerType& ppPairs,
csPairContainerType& pwPairs,
bool force = false);
bool boundaryBroadSearch(
uint32 bndryIndex,
const timeInfo& ti,
csPairContainerType& ppPairs,
csPairContainerType& pwPairs,
bool force = false);
virtual
bool ppEnterBroadSearch()const = 0;
real sizeRatio()const = 0;
virtual
bool pwEnterBroadSearch()const = 0;
virtual
bool ppPerformedBroadSearch()const = 0;
virtual
bool pwPerformedBroadSearch()const = 0;
real cellExtent()const = 0;
static
uniquePtr<contactSearch> create(
const dictionary& dict,
const box& domain,
const box& extDomain,
const particles& prtcl,
const geometry& geom,
Timers& timers);

View File

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

View File

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

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

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