From 0df384f54622c0e112da3288eaaa06763fffaac3 Mon Sep 17 00:00:00 2001 From: Hamidreza Norouzi Date: Fri, 26 Jan 2024 01:10:10 -0800 Subject: [PATCH] sphereParticles tested --- src/Particles/CMakeLists.txt | 1 + .../sphereParticles/sphereParticles.cpp | 132 ++++++++++++------ .../sphereParticles/sphereParticles.hpp | 30 ++++ .../sphereParticlesKernels.cpp | 100 +++++++++++++ .../sphereParticlesKernels.hpp | 54 +++---- .../sphereShape/sphereShape.cpp | 9 +- .../sphereShape/sphereShape.hpp | 4 +- .../dynamicPointStructure.cpp | 31 ---- src/Particles/particles/baseShapeNames.cpp | 7 +- src/Particles/particles/baseShapeNames.hpp | 2 +- src/Particles/particles/particles.cpp | 93 ------------ src/Particles/particles/particles.hpp | 58 ++------ src/Particles/particles/shape.cpp | 15 +- src/Particles/particles/shape.hpp | 6 +- .../containers/Vector/stdVectorHelper.hpp | 2 +- .../streams/dataIO/dataIOTemplate.cpp | 6 +- 16 files changed, 277 insertions(+), 273 deletions(-) create mode 100644 src/Particles/SphereParticles/sphereParticles/sphereParticlesKernels.cpp diff --git a/src/Particles/CMakeLists.txt b/src/Particles/CMakeLists.txt index 1b01ef94..b5c78045 100644 --- a/src/Particles/CMakeLists.txt +++ b/src/Particles/CMakeLists.txt @@ -7,6 +7,7 @@ particles/particles.cpp #particles/particleIdHandler.cpp SphereParticles/sphereShape/sphereShape.cpp SphereParticles/sphereParticles/sphereParticles.cpp +SphereParticles/sphereParticles/sphereParticlesKernels.cpp #Insertion/shapeMixture/shapeMixture.cpp #Insertion/insertion/insertion.cpp #Insertion/Insertion/Insertions.cpp diff --git a/src/Particles/SphereParticles/sphereParticles/sphereParticles.cpp b/src/Particles/SphereParticles/sphereParticles/sphereParticles.cpp index 04f4f19c..fcb7f757 100644 --- a/src/Particles/SphereParticles/sphereParticles/sphereParticles.cpp +++ b/src/Particles/SphereParticles/sphereParticles/sphereParticles.cpp @@ -21,9 +21,10 @@ Licence: #include "sphereParticles.hpp" #include "systemControl.hpp" #include "vocabs.hpp" -//#include "setFieldList.hpp" -//#include "sphereParticlesKernels.hpp" +#include "sphereParticlesKernels.hpp" + +//#include "setFieldList.hpp" /*pFlow::uniquePtr> pFlow::sphereParticles::getFieldObjectList()const { @@ -205,16 +206,36 @@ bool pFlow::sphereParticles::initInertia() exeSpace, Kokkos::IndexType>; - auto aPointsMask = dynPointStruct().activePointsMaskDevice(); - auto aRange = aPointsMask.activeRange(); + auto [minIndex, maxIndex] = minMax(shapeIndex().internal()); - auto field_I = I_.fieldDevice(); - auto field_shapeIndex = shapeIndex().fieldDevice(); - - const auto& shp = getShapes(); - realVector_D I ("I", shp.Inertia()); - auto d_I = I.deviceVector(); + if( !spheres_.indexValid(maxIndex) ) + { + fatalErrorInFunction<< + "the maximum value of shapeIndex is "<< maxIndex << + " which is not valid."<timers() ) { - + auto intMethod = control.settingsDict().getVal("integrationMethod"); REPORT(1)<<"Creating integration method "<needSetInitialVals()) - { - - auto [minInd, maxInd] = pStruct().activeRange(); - int32IndexContainer indexHD(minInd, maxInd); - - auto n = indexHD.size(); - auto index = indexHD.indicesHost(); - - realx3Vector rvel(n,RESERVE()); - const auto hrVel = rVelocity_.hostVector(); - - for(auto i=0; isetInitialVals(indexHD, rvel); - - }*/ - if(!initInertia()) { + fatalErrorInFunction; fatalExit; } @@ -437,18 +473,16 @@ bool pFlow::sphereParticles::iterate() particles::iterate(); accelerationTimer_.start(); - WARNING<<"pFlow::sphereParticlesKernels::acceleration"<, + Kokkos::IndexType>; + +void pFlow::sphereParticlesKernels::addMassDiamInertiaProp +( + deviceViewType1D shapeIndex, + deviceViewType1D mass, + deviceViewType1D diameter, + deviceViewType1D I, + deviceViewType1D propertyId, + pFlagTypeDevice incld, + deviceViewType1D src_mass, + deviceViewType1D src_diameter, + deviceViewType1D src_I, + deviceViewType1D src_propertyId +) +{ + auto aRange = incld.activeRange(); + + Kokkos::parallel_for( + "particles::initInertia", + policy(aRange.start(), aRange.end()), + LAMBDA_HD(uint32 i) + { + if(incld(i)) + { + uint32 index = shapeIndex[i]; + I[i] = src_I[index]; + diameter[i] = src_diameter[index]; + mass[i] = src_mass[index]; + propertyId[i] = src_propertyId[index]; + } + }); + +} + +void pFlow::sphereParticlesKernels::acceleration +( + realx3 g, + deviceViewType1D mass, + deviceViewType1D force, + deviceViewType1D I, + deviceViewType1D torque, + pFlagTypeDevice incld, + deviceViewType1D lAcc, + deviceViewType1D rAcc +) +{ + + auto activeRange = incld.activeRange(); + if(incld.isAllActive()) + { + Kokkos::parallel_for( + "pFlow::sphereParticlesKernels::acceleration", + policy(activeRange.start(), activeRange.end()), + LAMBDA_HD(uint32 i){ + lAcc[i] = force[i]/mass[i] + g; + rAcc[i] = torque[i]/I[i]; + }); + } + else + { + Kokkos::parallel_for( + "pFlow::sphereParticlesKernels::acceleration", + policy(activeRange.start(), activeRange.end()), + LAMBDA_HD(uint32 i){ + if(incld(i)) + { + lAcc[i] = force[i]/mass[i] + g; + rAcc[i] = torque[i]/I[i]; + } + }); + + } + + Kokkos::fence(); +} diff --git a/src/Particles/SphereParticles/sphereParticles/sphereParticlesKernels.hpp b/src/Particles/SphereParticles/sphereParticles/sphereParticlesKernels.hpp index 2edeff2b..f2ebbe22 100644 --- a/src/Particles/SphereParticles/sphereParticles/sphereParticlesKernels.hpp +++ b/src/Particles/SphereParticles/sphereParticles/sphereParticlesKernels.hpp @@ -21,56 +21,36 @@ Licence: #ifndef __sphereParticlesKernels_hpp__ #define __sphereParticlesKernels_hpp__ +#include "types.hpp" +#include "pointFlag.hpp" + namespace pFlow::sphereParticlesKernels { -using rpAcceleration = Kokkos::RangePolicy< - DefaultExecutionSpace, - Kokkos::Schedule, - Kokkos::IndexType - >; +void addMassDiamInertiaProp( + deviceViewType1D shapeIndex, + deviceViewType1D mass, + deviceViewType1D diameter, + deviceViewType1D I, + deviceViewType1D propertyId, + pFlagTypeDevice incld, + deviceViewType1D src_mass, + deviceViewType1D src_diameter, + deviceViewType1D src_I, + deviceViewType1D src_propertyId +); -template void acceleration( realx3 g, deviceViewType1D mass, deviceViewType1D force, deviceViewType1D I, deviceViewType1D torque, - IncludeFunctionType incld, + pFlagTypeDevice incld, deviceViewType1D lAcc, deviceViewType1D rAcc - ) -{ +); - auto activeRange = incld.activeRange(); - if(incld.allActive()) - { - Kokkos::parallel_for( - "pFlow::sphereParticlesKernels::acceleration", - rpAcceleration(activeRange.first, activeRange.second), - LAMBDA_HD(int32 i){ - lAcc[i] = force[i]/mass[i] + g; - rAcc[i] = torque[i]/I[i]; - }); - } - else - { - Kokkos::parallel_for( - "pFlow::sphereParticlesKernels::acceleration", - rpAcceleration(activeRange.first, activeRange.second), - LAMBDA_HD(int32 i){ - if(incld(i)) - { - lAcc[i] = force[i]/mass[i] + g; - rAcc[i] = torque[i]/I[i]; - } - }); - - } - - Kokkos::fence(); -} } diff --git a/src/Particles/SphereParticles/sphereShape/sphereShape.cpp b/src/Particles/SphereParticles/sphereShape/sphereShape.cpp index f07942ce..75ec21e3 100644 --- a/src/Particles/SphereParticles/sphereShape/sphereShape.cpp +++ b/src/Particles/SphereParticles/sphereShape/sphereShape.cpp @@ -21,9 +21,9 @@ Licence: #include "sphereShape.hpp" -bool pFlow::sphereShape::readDictionary() +bool pFlow::sphereShape::readFromDictionary3() { - + diameters_ = getVal("diameters"); if(diameters_.size() != numShapes() ) @@ -61,6 +61,11 @@ pFlow::sphereShape::sphereShape shape(fileName, owner, prop) { + if(!readFromDictionary3()) + { + fatalExit; + fatalErrorInFunction; + } } pFlow::real pFlow::sphereShape::maxBoundingSphere() const diff --git a/src/Particles/SphereParticles/sphereShape/sphereShape.hpp b/src/Particles/SphereParticles/sphereShape/sphereShape.hpp index 6c03bfd4..b59625ce 100644 --- a/src/Particles/SphereParticles/sphereShape/sphereShape.hpp +++ b/src/Particles/SphereParticles/sphereShape/sphereShape.hpp @@ -26,8 +26,6 @@ Licence: namespace pFlow { - - class sphereShape : public shape @@ -37,7 +35,7 @@ private: // - diameter of spheres realVector diameters_; - bool readDictionary(); + bool readFromDictionary3(); protected: diff --git a/src/Particles/dynamicPointStructure/dynamicPointStructure.cpp b/src/Particles/dynamicPointStructure/dynamicPointStructure.cpp index d99766e2..122d2392 100644 --- a/src/Particles/dynamicPointStructure/dynamicPointStructure.cpp +++ b/src/Particles/dynamicPointStructure/dynamicPointStructure.cpp @@ -76,37 +76,6 @@ pFlow::dynamicPointStructure::dynamicPointStructure fatalExit; } - - WARNING << "Initialization of integrationPos_ and integrationVel_ should be doen!"<needSetInitialVals()) return; - - - - auto [minInd, maxInd] = pStruct().activeRange(); - int32IndexContainer indexHD(minInd, maxInd); - - auto n = indexHD.size(); - auto index = indexHD.indicesHost(); - - realx3Vector pos(n,RESERVE()); - realx3Vector vel(n,RESERVE()); - const auto hVel = velocity().hostVector(); - const auto hPos = pStruct().pointPosition().hostVector(); - - for(auto i=0; isetInitialVals(indexHD, pos); - - REPORT(2)<< "Initializing the required vectors for velocity integratoin\n "<setInitialVals(indexHD, vel);*/ - } bool pFlow::dynamicPointStructure::predict diff --git a/src/Particles/particles/baseShapeNames.cpp b/src/Particles/particles/baseShapeNames.cpp index c8200f0c..f20efc06 100644 --- a/src/Particles/particles/baseShapeNames.cpp +++ b/src/Particles/particles/baseShapeNames.cpp @@ -40,7 +40,7 @@ bool pFlow::baseShapeNames::createHashNames() } -bool pFlow::baseShapeNames::readFromDictionary() +bool pFlow::baseShapeNames::readFromDictionary1() { shapeNames_ = getVal("names"); @@ -69,18 +69,21 @@ pFlow::baseShapeNames::baseShapeNames ) { - if( !readFromDictionary() ) + if( !readFromDictionary1() ) { + fatalErrorInFunction; fatalExit; } if( !createHashNames()) { + fatalErrorInFunction; fatalExit; } } + bool pFlow::baseShapeNames::writeToDict(dictionary &dict)const { diff --git a/src/Particles/particles/baseShapeNames.hpp b/src/Particles/particles/baseShapeNames.hpp index 9c427c32..8f4ca44e 100644 --- a/src/Particles/particles/baseShapeNames.hpp +++ b/src/Particles/particles/baseShapeNames.hpp @@ -46,7 +46,7 @@ private: bool createHashNames(); - bool readFromDictionary(); + bool readFromDictionary1(); protected: diff --git a/src/Particles/particles/particles.cpp b/src/Particles/particles/particles.cpp index ff32aa4f..cbda5c61 100644 --- a/src/Particles/particles/particles.cpp +++ b/src/Particles/particles/particles.cpp @@ -21,60 +21,6 @@ Licence: #include "particles.hpp" -bool pFlow::particles::initMassDiameter()const -{ - - auto [minIndex, maxIndex] = minMax(shapeIndex_.internal()); - - const auto& shp = getShapes(); - - if( !shp.indexValid(maxIndex) ) - { - fatalErrorInFunction<< - "the maximum value of shape index is "<< maxIndex << - " which is not valid."<>; - - auto field_diameter = diameter_.fieldDevice(); - auto field_mass = mass_.fieldDevice(); - auto field_propId = propertyId_.fieldDevice(); - auto field_shapeIndex = shapeIndex_.fieldDevice(); - - auto d_d = d.deviceVector(); - auto d_mass = mass.deviceVector(); - auto d_propId = propId.deviceVector(); - - Kokkos::parallel_for( - "particles::initMassDiameter", - policy(aRange.start(), aRange.end()), - LAMBDA_HD(uint32 i) - { - if(aPointsMask(i)) - { - uint32 index = field_shapeIndex[i]; - field_diameter[i] = d_d[index]; - field_mass[i] = d_mass[index]; - field_propId[i] = d_propId[index]; - } - }); - Kokkos::fence(); - - return true; -} pFlow::particles::particles ( @@ -96,18 +42,6 @@ pFlow::particles::particles dynPointStruct_, static_cast(-1) ), - propertyId_ - ( - objectFile - ( - "propertyId", - "", - objectFile::READ_ALWAYS, - objectFile::WRITE_NEVER - ), - dynPointStruct_, - static_cast(0) - ), shapeIndex_ ( objectFile @@ -120,26 +54,6 @@ pFlow::particles::particles dynPointStruct_, 0 ), - diameter_ - ( - objectFile( - "diameter", - "", - objectFile::READ_NEVER, - objectFile::WRITE_NEVER), - dynPointStruct_, - 0.00000000001 - ), - mass_ - ( - objectFile( - "mass", - "", - objectFile::READ_NEVER, - objectFile::WRITE_NEVER), - dynPointStruct_, - 0.0000000001 - ), accelertion_ ( objectFile( @@ -167,14 +81,7 @@ pFlow::particles::particles dynPointStruct_, zero3) { - this->addToSubscriber(dynPointStruct_); - - if(!initMassDiameter()) - { - fatalExit; - } - } bool pFlow::particles::beforeIteration() diff --git a/src/Particles/particles/particles.hpp b/src/Particles/particles/particles.hpp index 5f9e7ea8..6b847e63 100644 --- a/src/Particles/particles/particles.hpp +++ b/src/Particles/particles/particles.hpp @@ -43,18 +43,9 @@ private: /// id of particles on host uint32PointField_D id_; - /// property id on device - uint8PointField_D propertyId_; - - /// property id on device + uint32PointField_D shapeIndex_; - /// diameter / boundig sphere size of particles on device - realPointField_D diameter_; - - /// mass of particles field - realPointField_D mass_; - /// acceleration on device realx3PointField_D accelertion_; @@ -73,8 +64,7 @@ private: message::ITEM_DELETE ); - bool initMassDiameter()const; - + void zeroForce() { contactForce_.fill(zero3); @@ -152,29 +142,6 @@ public: return dynPointStruct_.velocity(); } - inline - const auto& diameter()const - { - return diameter_; - } - - inline - auto& diameter() - { - return diameter_; - } - - inline - const auto& mass()const - { - return mass_; - } - - inline auto& mass() - { - return mass_; - } - inline const auto& accelertion()const { @@ -211,18 +178,6 @@ public: return contactTorque_; } - inline - const auto& propertyId()const - { - return propertyId_; - } - - inline - auto& propertyId() - { - return propertyId_; - } - bool beforeIteration() override; bool iterate()override; @@ -237,7 +192,14 @@ public: const setFieldList& setField ) = 0;*/ + virtual + const uint32PointField_D& propertyId()const = 0; + virtual + const realPointField_D& diameter()const = 0; + + virtual + const realPointField_D& mass()const = 0; virtual realx3PointField_D& rAcceleration() = 0; @@ -255,7 +217,7 @@ public: const shape& getShapes()const = 0; virtual - void boundingSphereMinMax(real & minDiam, real& maxDiam)const; + void boundingSphereMinMax(real & minDiam, real& maxDiam)const = 0; diff --git a/src/Particles/particles/shape.cpp b/src/Particles/particles/shape.cpp index b73fe9c0..bcf66e74 100644 --- a/src/Particles/particles/shape.cpp +++ b/src/Particles/particles/shape.cpp @@ -7,24 +7,23 @@ bool pFlow::shape::findPropertyIds() ForAll( i, materialNames_) { - if(uint32 propId; property_.nameToIndex(materialNames_[i], propId) ) { - shapePropertyIds_[i] = static_cast(propId); + shapePropertyIds_[i] = propId; } else { fatalErrorInFunction<<"Material name "<< materialNames_[i]<< - "is not valid in dictionary "<("materials"); if(materialNames_.size() != numShapes() ) @@ -49,13 +48,15 @@ pFlow::shape::shape property_(prop) { - if( !readFromDictionary() ) + if( !readFromDictionary2() ) { + fatalErrorInFunction; fatalExit; } if(!findPropertyIds()) { + fatalErrorInFunction; fatalExit; } diff --git a/src/Particles/particles/shape.hpp b/src/Particles/particles/shape.hpp index a85a85ae..1ce5d8ce 100644 --- a/src/Particles/particles/shape.hpp +++ b/src/Particles/particles/shape.hpp @@ -37,7 +37,7 @@ private: const property& property_; /// list of property ids of shapes (index) - uint8Vector shapePropertyIds_; + uint32Vector shapePropertyIds_; /// list of material names wordVector materialNames_; @@ -45,7 +45,7 @@ private: bool findPropertyIds(); - bool readFromDictionary(); + bool readFromDictionary2(); protected: @@ -96,7 +96,7 @@ public: } inline - bool propIdValid(int8 propId)const + bool propIdValid(uint32 propId)const { return static_cast(propId) < property_.numMaterials(); } diff --git a/src/phasicFlow/containers/Vector/stdVectorHelper.hpp b/src/phasicFlow/containers/Vector/stdVectorHelper.hpp index 60a71873..37a85a6f 100644 --- a/src/phasicFlow/containers/Vector/stdVectorHelper.hpp +++ b/src/phasicFlow/containers/Vector/stdVectorHelper.hpp @@ -111,7 +111,7 @@ bool readStdVector std::vector& vec ) { - + return readDataAsciiBinary(is, vec); } diff --git a/src/phasicFlow/streams/dataIO/dataIOTemplate.cpp b/src/phasicFlow/streams/dataIO/dataIOTemplate.cpp index 1682d58e..74efd812 100644 --- a/src/phasicFlow/streams/dataIO/dataIOTemplate.cpp +++ b/src/phasicFlow/streams/dataIO/dataIOTemplate.cpp @@ -130,12 +130,15 @@ bool pFlow::readDataAscii token firstToken(is); size_t len = 0; + bool lenFound = false; if( firstToken.isInt64()) { len = firstToken.int64Token(); + lenFound = true; vec.reserve(len); firstToken = token(is); } + T val{}; if( firstToken.isPunctuation() ) // start of vector @@ -164,6 +167,7 @@ bool pFlow::readDataAscii is >> lastToken; is.fatalCheck(FUNCTION_NAME); + } } else @@ -174,7 +178,7 @@ bool pFlow::readDataAscii return false; } - if(len>0&& len != vec.size()) + if(lenFound && len>0 && len != vec.size()) { warningInFunction<<"vector length specified "<< len << " is different from number of elements "<< vec.size()<