From 614b2f732e18fb1e446e79393f185cd659c87d1a Mon Sep 17 00:00:00 2001 From: HRN Date: Sun, 12 May 2024 19:10:04 +0330 Subject: [PATCH] develop branch update for changes in MPI branch for data transfer stage mortong indexing is added (messaging is not complete) --- .../contactLists/unsortedContactList.hpp | 5 +- .../ContactSearch/ContactSearch.hpp | 20 ++- .../boundaryContactSearch.hpp | 5 - .../contactSearch/contactSearch.hpp | 9 ++ .../methods/cellBased/NBS/NBS.hpp | 10 ++ .../methods/particleWallContactSearchs.hpp | 15 +++ src/Interaction/interaction/interaction.cpp | 2 +- .../boundarySphereInteraction.cpp | 17 ++- .../boundarySphereInteraction.hpp | 38 ++++-- .../periodicBoundarySphereInteraction.cpp | 12 +- .../periodicBoundarySphereInteraction.hpp | 3 +- .../sphereInteraction/sphereInteraction.cpp | 107 +++++++++++++--- .../sphereInteraction/sphereInteraction.hpp | 3 + .../sphereParticles/sphereParticles.cpp | 7 +- .../sphereParticles/sphereParticles.hpp | 2 +- .../dynamicPointStructure.cpp | 11 +- .../dynamicPointStructure.hpp | 2 + src/Particles/particles/particles.cpp | 11 +- src/Particles/particles/particles.hpp | 2 + src/phasicFlow/CMakeLists.txt | 2 + .../containers/VectorHD/VectorSingle.cpp | 88 ++++++++++++- .../containers/VectorHD/VectorSingle.hpp | 8 +- .../indexContainer/indexContainer.hpp | 20 ++- .../boundaryField/boundaryField.hpp | 18 ++- .../generalBoundary/generalBoundary.hpp | 4 - .../internalField/internalField.cpp | 32 ++++- .../internalField/internalField.hpp | 2 + src/phasicFlow/eventManagement/message.hpp | 12 +- .../repository/Time/baseTimeControl.cpp | 10 ++ .../repository/Time/baseTimeControl.hpp | 7 + src/phasicFlow/streams/dataIO/dataIO.cpp | 4 +- .../boundaries/boundaryBase/boundaryBase.cpp | 5 +- .../boundaries/boundaryBase/boundaryBase.hpp | 39 +++++- .../boundaryBase/scatteredFieldAccess.hpp | 10 ++ .../boundaries/boundaryExit/boundaryExit.cpp | 3 + .../boundaries/boundaryList.cpp | 36 +++--- .../boundaries/boundaryList.hpp | 14 +- .../boundaries/boundaryNone/boundaryNone.hpp | 8 +- .../boundaryPeriodic/boundaryPeriodic.cpp | 15 ++- .../boundaryPeriodic/boundaryPeriodic.hpp | 3 +- src/phasicFlow/structuredData/cells/cells.hpp | 107 +++++----------- .../domain/regularSimulationDomain.cpp | 13 +- .../internalPoints/internalPoints.cpp | 121 +++++++++++++++++- .../internalPoints/internalPoints.hpp | 4 + .../internalPoints/pointFlag.hpp | 27 +++- .../internalPoints/pointFlagKernels.hpp | 45 +++++++ .../{ => pointSorting}/mortonIndexing.cpp | 60 ++++----- .../{ => pointSorting}/mortonIndexing.hpp | 12 +- .../pointSorting/pointSorting.cpp | 52 ++++++++ .../pointSorting/pointSorting.hpp | 68 ++++++++++ .../pointStructure/pointStructure.cpp | 64 ++++++++- .../pointStructure/pointStructure.hpp | 11 +- 52 files changed, 937 insertions(+), 268 deletions(-) rename src/phasicFlow/structuredData/pointStructure/pointStructure/{ => pointSorting}/mortonIndexing.cpp (59%) rename src/phasicFlow/structuredData/pointStructure/pointStructure/{ => pointSorting}/mortonIndexing.hpp (94%) create mode 100644 src/phasicFlow/structuredData/pointStructure/pointStructure/pointSorting/pointSorting.cpp create mode 100644 src/phasicFlow/structuredData/pointStructure/pointStructure/pointSorting/pointSorting.hpp diff --git a/src/Interaction/contactLists/unsortedContactList.hpp b/src/Interaction/contactLists/unsortedContactList.hpp index 5922c4dc..8a61b2ce 100644 --- a/src/Interaction/contactLists/unsortedContactList.hpp +++ b/src/Interaction/contactLists/unsortedContactList.hpp @@ -95,7 +95,8 @@ public: // swap conainer and values swapViews(values0_, values_); swapViews(container0_, this->container_); - return UnsortedPairs::beforeBroadSearch(); + UnsortedPairs::beforeBroadSearch(); + return true; } bool afterBroadSearch() @@ -110,7 +111,7 @@ public: rpFillPairs(0,this->capacity()), *this); Kokkos::fence(); - + return true; } diff --git a/src/Interaction/contactSearch/ContactSearch/ContactSearch.hpp b/src/Interaction/contactSearch/ContactSearch/ContactSearch.hpp index d80813ca..f2dda414 100644 --- a/src/Interaction/contactSearch/ContactSearch/ContactSearch.hpp +++ b/src/Interaction/contactSearch/ContactSearch/ContactSearch.hpp @@ -160,7 +160,8 @@ public: csPairContainerType& pwPairs, bool force = false)override { - Particles().boundingSphere().updateBoundaries(DataDirection::SlaveToMaster); + if(i==0u) + Particles().boundingSphere().updateBoundaries(DataDirection::SlaveToMaster); return csBoundaries_[i].broadSearch( iter, t, @@ -183,6 +184,23 @@ public: { return ppwContactSearch_().performedSearch(); } + + + uint32 updateInterval()const override + { + return ppwContactSearch_().updateInterval(); + } + + real sizeRatio()const override + { + return ppwContactSearch_().sizeRatio(); + } + + + real cellExtent()const override + { + return ppwContactSearch_().cellExtent(); + } }; diff --git a/src/Interaction/contactSearch/boundaries/boundaryContactSearch/boundaryContactSearch.hpp b/src/Interaction/contactSearch/boundaries/boundaryContactSearch/boundaryContactSearch.hpp index 78c61bb8..55ecb47e 100644 --- a/src/Interaction/contactSearch/boundaries/boundaryContactSearch/boundaryContactSearch.hpp +++ b/src/Interaction/contactSearch/boundaries/boundaryContactSearch/boundaryContactSearch.hpp @@ -72,11 +72,6 @@ public: return contactSearch_; } - void fill(const std::any &val) override - { - return; - } - bool hearChanges( real t, real dt, diff --git a/src/Interaction/contactSearch/contactSearch/contactSearch.hpp b/src/Interaction/contactSearch/contactSearch/contactSearch.hpp index 235403c4..865ec291 100644 --- a/src/Interaction/contactSearch/contactSearch/contactSearch.hpp +++ b/src/Interaction/contactSearch/contactSearch/contactSearch.hpp @@ -140,6 +140,15 @@ public: virtual bool performedBroadSearch()const = 0; + virtual + uint32 updateInterval()const = 0; + + virtual + real sizeRatio()const = 0; + + virtual + real cellExtent()const = 0; + static uniquePtr create( const dictionary& dict, diff --git a/src/Interaction/contactSearch/methods/cellBased/NBS/NBS.hpp b/src/Interaction/contactSearch/methods/cellBased/NBS/NBS.hpp index 11d2e108..9cd5ab6f 100644 --- a/src/Interaction/contactSearch/methods/cellBased/NBS/NBS.hpp +++ b/src/Interaction/contactSearch/methods/cellBased/NBS/NBS.hpp @@ -134,6 +134,16 @@ public: { return 1; } + + real sizeRatio()const + { + return sizeRatio_; + } + + real cellExtent()const + { + return cellExtent_; + } auto getCellIterator([[maybe_unused]] uint32 lvl)const { diff --git a/src/Interaction/contactSearch/methods/particleWallContactSearchs.hpp b/src/Interaction/contactSearch/methods/particleWallContactSearchs.hpp index d72f9404..1bf5dab3 100644 --- a/src/Interaction/contactSearch/methods/particleWallContactSearchs.hpp +++ b/src/Interaction/contactSearch/methods/particleWallContactSearchs.hpp @@ -113,6 +113,21 @@ public: return false; } + uint32 updateInterval()const + { + return updateInterval_; + } + + real sizeRatio()const + { + return getMethod().sizeRatio(); + } + + real cellExtent()const + { + return getMethod().cellExtent(); + } + }; diff --git a/src/Interaction/interaction/interaction.cpp b/src/Interaction/interaction/interaction.cpp index c6099405..03abbb4b 100644 --- a/src/Interaction/interaction/interaction.cpp +++ b/src/Interaction/interaction/interaction.cpp @@ -72,7 +72,7 @@ pFlow::uniquePtr pFlow::interaction::create gSettings::sleepMiliSeconds( - 1000*(pFlowProcessors().localSize()-pFlowProcessors().localRank()-1)); + 100*(pFlowProcessors().localSize()-pFlowProcessors().localRank()-1)); pOutput.space(2)<<"Creating interaction "< +void pFlow::boundarySphereInteraction::allocatePPPairs() +{ + ppPairs_.reset(nullptr); + ppPairs_ = makeUnique(1); +} + +template +void pFlow::boundarySphereInteraction::allocatePWPairs() +{ + pwPairs_.reset(nullptr); + pwPairs_ = makeUnique(1); +} + template pFlow::boundarySphereInteraction::boundarySphereInteraction( @@ -28,8 +43,6 @@ pFlow::boundarySphereInteraction::boundarySphereInteraction( geometryMotion_(geomMotion), sphParticles_(sphPrtcls) { - ppPairs_ = makeUnique(1); - pwPairs_ = makeUnique(1); } template diff --git a/src/Interaction/sphereInteraction/boundaries/boundarySphereInteraction/boundarySphereInteraction.hpp b/src/Interaction/sphereInteraction/boundaries/boundarySphereInteraction/boundarySphereInteraction.hpp index 56058b6d..708fb61e 100644 --- a/src/Interaction/sphereInteraction/boundaries/boundarySphereInteraction/boundarySphereInteraction.hpp +++ b/src/Interaction/sphereInteraction/boundaries/boundarySphereInteraction/boundarySphereInteraction.hpp @@ -22,7 +22,7 @@ Licence: #include "virtualConstructor.hpp" #include "generalBoundary.hpp" -#include "unsortedContactList.hpp" +#include "sortedContactList.hpp" #include "sphereParticles.hpp" namespace pFlow @@ -51,7 +51,7 @@ public: using IndexType = uint32; using ContactListType = - unsortedContactList; + sortedContactList; private: @@ -60,9 +60,15 @@ private: /// const reference to sphere particles const sphereParticles& sphParticles_; - uniquePtr ppPairs_; + uniquePtr ppPairs_ = nullptr; - uniquePtr pwPairs_; + uniquePtr pwPairs_ = nullptr; + +protected: + + void allocatePPPairs(); + + void allocatePWPairs(); public: @@ -124,15 +130,30 @@ public: return pwPairs_(); } + bool ppPairsAllocated()const + { + if( ppPairs_)return true; + return false; + } + + bool pwPairsAllocated()const + { + if( pwPairs_)return true; + return false; + } + virtual bool sphereSphereInteraction( real dt, - const ContactForceModel& cfModel) + const ContactForceModel& cfModel, + uint32 step) { // for default boundary, no thing to be done - return true; + return false; } + + bool hearChanges ( real t, @@ -149,11 +170,6 @@ public: return true; } - void fill(const std::any& val)override - { - notImplementedFunction; - } - static uniquePtr create( const boundaryBase& boundary, diff --git a/src/Interaction/sphereInteraction/boundaries/periodicBoundarySphereInteraction/periodicBoundarySphereInteraction.cpp b/src/Interaction/sphereInteraction/boundaries/periodicBoundarySphereInteraction/periodicBoundarySphereInteraction.cpp index 38a4336b..3200c4e7 100644 --- a/src/Interaction/sphereInteraction/boundaries/periodicBoundarySphereInteraction/periodicBoundarySphereInteraction.cpp +++ b/src/Interaction/sphereInteraction/boundaries/periodicBoundarySphereInteraction/periodicBoundarySphereInteraction.cpp @@ -30,7 +30,10 @@ pFlow::periodicBoundarySphereInteraction::periodicBoundarySphereIntera { if(boundary.thisBoundaryIndex()%2==1) { - masterInteraction_ = true; + masterInteraction_ = true; + this->allocatePPPairs(); + this->allocatePWPairs(); + } else { @@ -42,10 +45,11 @@ template bool pFlow::periodicBoundarySphereInteraction::sphereSphereInteraction ( real dt, - const ContactForceModel &cfModel + const ContactForceModel &cfModel, + uint32 step ) { - if(!masterInteraction_) return true; + if(!masterInteraction_) return false; pFlow::periodicBoundarySIKernels::sphereSphereInteraction( dt, @@ -61,5 +65,5 @@ bool pFlow::periodicBoundarySphereInteraction::sphereSphereInteraction this->sphParticles().contactForce().deviceViewAll(), this->sphParticles().contactTorque().deviceViewAll()); - return true; + return false; } \ No newline at end of file diff --git a/src/Interaction/sphereInteraction/boundaries/periodicBoundarySphereInteraction/periodicBoundarySphereInteraction.hpp b/src/Interaction/sphereInteraction/boundaries/periodicBoundarySphereInteraction/periodicBoundarySphereInteraction.hpp index e2887787..5f1e7a45 100644 --- a/src/Interaction/sphereInteraction/boundaries/periodicBoundarySphereInteraction/periodicBoundarySphereInteraction.hpp +++ b/src/Interaction/sphereInteraction/boundaries/periodicBoundarySphereInteraction/periodicBoundarySphereInteraction.hpp @@ -83,7 +83,8 @@ public: bool sphereSphereInteraction( real dt, - const ContactForceModel& cfModel)override; + const ContactForceModel& cfModel, + uint32 step)override; }; diff --git a/src/Interaction/sphereInteraction/sphereInteraction/sphereInteraction.cpp b/src/Interaction/sphereInteraction/sphereInteraction/sphereInteraction.cpp index b4546521..2aa10ef8 100644 --- a/src/Interaction/sphereInteraction/sphereInteraction/sphereInteraction.cpp +++ b/src/Interaction/sphereInteraction/sphereInteraction/sphereInteraction.cpp @@ -137,7 +137,9 @@ pFlow::sphereInteraction::sphereInteraction ppInteractionTimer_("sphere-sphere interaction", &this->timers()), pwInteractionTimer_("sphere-wall interaction", &this->timers()), contactListMangementTimer_("contact-list management", &this->timers()), - boundaryInteractionTimer_("interaction for boundary", &this->timers()) + boundaryContactSearchTimer_("contact search for boundary", &this->timers()), + boundaryInteractionTimer_("interaction for boundary", &this->timers()), + contactListBoundaryTimer_("contact-list management for boundary", &this->timers()) { if(!createSphereInteraction()) @@ -179,11 +181,14 @@ bool pFlow::sphereInteraction::iterate() contactListMangementTimer_.pause(); } + contactListBoundaryTimer_.start(); for(uint32 i=0; i<6u; i++) { - boundaryInteraction_[i].ppPairs().beforeBroadSearch(); - boundaryInteraction_[i].pwPairs().beforeBroadSearch(); + auto& BI = boundaryInteraction_[i]; + if(BI.ppPairsAllocated()) BI.ppPairs().beforeBroadSearch(); + if(BI.pwPairsAllocated()) BI.pwPairs().beforeBroadSearch(); } + contactListBoundaryTimer_.pause(); if( sphParticles_.numActive()<=0)return true; @@ -200,21 +205,27 @@ bool pFlow::sphereInteraction::iterate() fatalExit; } + boundaryContactSearchTimer_.start(); for(uint32 i=0; i<6u; i++) { - if( !contactSearch_().boundaryBroadSearch( - i, - iter, - t, - dt, - boundaryInteraction_[i].ppPairs(), - boundaryInteraction_[i].pwPairs())) + auto& BI =boundaryInteraction_[i]; + if(BI.ppPairsAllocated()) { - fatalErrorInFunction<< - "failed to perform broadSearch for boundary index "<::iterate() contactListMangementTimer_.end(); } + contactListBoundaryTimer_.resume(); for(uint32 i=0; i<6u; i++) { - boundaryInteraction_[i].ppPairs().afterBroadSearch(); - boundaryInteraction_[i].pwPairs().afterBroadSearch(); + auto& BI = boundaryInteraction_[i]; + if(BI.ppPairsAllocated()) BI.ppPairs().afterBroadSearch(); + if(BI.pwPairsAllocated()) BI.pwPairs().afterBroadSearch(); } + contactListBoundaryTimer_.end(); + + { + boundaryInteractionTimer_.start(); + std::array requireStep{ + boundaryInteraction_[0].isBoundaryMaster(), + boundaryInteraction_[1].isBoundaryMaster(), + boundaryInteraction_[2].isBoundaryMaster(), + boundaryInteraction_[3].isBoundaryMaster(), + boundaryInteraction_[4].isBoundaryMaster(), + boundaryInteraction_[5].isBoundaryMaster()}; + int step = 1; + const auto& cfModel = this->forceModel_(); + while( std::any_of(requireStep.begin(), requireStep.end(), [](bool v) { return v==true; })) + { + for(uint32 i=0; i<6u; i++) + { + if(step==1u || requireStep[i] ) + { + requireStep[i] = boundaryInteraction_[i].sphereSphereInteraction( + dt, + this->forceModel_(), + step + ); + } + } + step++; + } + boundaryInteractionTimer_.pause(); + } ppInteractionTimer_.start(); sphereSphereInteraction(); ppInteractionTimer_.end(); @@ -239,14 +282,36 @@ bool pFlow::sphereInteraction::iterate() sphereWallInteraction(); pwInteractionTimer_.end(); - boundaryInteractionTimer_.start(); - for(uint32 i=0; i<6u; i++) { - boundaryInteraction_[i].sphereSphereInteraction( - dt, - this->forceModel_()); + boundaryInteractionTimer_.resume(); + std::array requireStep{ + !boundaryInteraction_[0].isBoundaryMaster(), + !boundaryInteraction_[1].isBoundaryMaster(), + !boundaryInteraction_[2].isBoundaryMaster(), + !boundaryInteraction_[3].isBoundaryMaster(), + !boundaryInteraction_[4].isBoundaryMaster(), + !boundaryInteraction_[5].isBoundaryMaster()}; + + int step = 2; + const auto& cfModel = this->forceModel_(); + while(std::any_of(requireStep.begin(), requireStep.end(), [](bool v) { return v==true; })) + { + for(uint32 i=0; i<6u; i++) + { + if(requireStep[i]) + { + requireStep[i] = boundaryInteraction_[i].sphereSphereInteraction( + dt, + this->forceModel_(), + step + ); + } + } + step++; } boundaryInteractionTimer_.end(); + } + return true; } diff --git a/src/Interaction/sphereInteraction/sphereInteraction/sphereInteraction.hpp b/src/Interaction/sphereInteraction/sphereInteraction/sphereInteraction.hpp index 4ccb9836..74cb2442 100644 --- a/src/Interaction/sphereInteraction/sphereInteraction/sphereInteraction.hpp +++ b/src/Interaction/sphereInteraction/sphereInteraction/sphereInteraction.hpp @@ -96,9 +96,12 @@ private: /// timer for managing contact lists (only inernal points) Timer contactListMangementTimer_; + Timer boundaryContactSearchTimer_; /// timer for boundary interaction time Timer boundaryInteractionTimer_; + Timer contactListBoundaryTimer_; + bool createSphereInteraction(); diff --git a/src/Particles/SphereParticles/sphereParticles/sphereParticles.cpp b/src/Particles/SphereParticles/sphereParticles/sphereParticles.cpp index ea29aee8..df7b98b8 100644 --- a/src/Particles/SphereParticles/sphereParticles/sphereParticles.cpp +++ b/src/Particles/SphereParticles/sphereParticles/sphereParticles.cpp @@ -389,7 +389,9 @@ pFlow::sphereParticles::sphereParticles( intPredictTimer_( "Integration-predict", &this->timers() ), intCorrectTimer_( - "Integration-correct", &this->timers() ) + "Integration-correct", &this->timers() ), + fieldUpdateTimer_( + "fieldUpdate", &this->timers() ) { auto intMethod = control.settingsDict().getVal("integrationMethod"); @@ -514,13 +516,14 @@ bool pFlow::sphereParticles::beforeIteration() rVelIntegration_().predict(dt,rVelocity_, rAcceleration_); intPredictTimer_.end(); + fieldUpdateTimer_.start(); propertyId_.updateBoundariesSlaveToMasterIfRequested(); diameter_.updateBoundariesSlaveToMasterIfRequested(); mass_.updateBoundariesSlaveToMasterIfRequested(); I_.updateBoundariesSlaveToMasterIfRequested(); rVelocity_.updateBoundariesSlaveToMasterIfRequested(); rAcceleration_.updateBoundariesSlaveToMasterIfRequested(); - + fieldUpdateTimer_.end(); return true; } diff --git a/src/Particles/SphereParticles/sphereParticles/sphereParticles.hpp b/src/Particles/SphereParticles/sphereParticles/sphereParticles.hpp index fc0aa82c..0f78e498 100644 --- a/src/Particles/SphereParticles/sphereParticles/sphereParticles.hpp +++ b/src/Particles/SphereParticles/sphereParticles/sphereParticles.hpp @@ -79,7 +79,7 @@ private: /// timer for integration computations (correction step) Timer intCorrectTimer_; - + Timer fieldUpdateTimer_; private: bool initializeParticles(); diff --git a/src/Particles/dynamicPointStructure/dynamicPointStructure.cpp b/src/Particles/dynamicPointStructure/dynamicPointStructure.cpp index c57475ec..6ad16556 100644 --- a/src/Particles/dynamicPointStructure/dynamicPointStructure.cpp +++ b/src/Particles/dynamicPointStructure/dynamicPointStructure.cpp @@ -38,6 +38,7 @@ pFlow::dynamicPointStructure::dynamicPointStructure *this, zero3 ), + velocityUpdateTimer_("velocity boundary update", &timers()), integrationMethod_ ( control.settingsDict().getVal("integrationMethod") @@ -81,11 +82,11 @@ pFlow::dynamicPointStructure::dynamicPointStructure bool pFlow::dynamicPointStructure::beforeIteration() { - return pointStructure::beforeIteration(); - /*real dt = this->dt(); - - auto& acc = time().lookupObject("acceleration"); - return predict(dt, acc);*/ + if(!pointStructure::beforeIteration())return false; + velocityUpdateTimer_.start(); + velocity_.updateBoundariesSlaveToMasterIfRequested(); + velocityUpdateTimer_.end(); + return true; } bool pFlow::dynamicPointStructure::iterate() diff --git a/src/Particles/dynamicPointStructure/dynamicPointStructure.hpp b/src/Particles/dynamicPointStructure/dynamicPointStructure.hpp index 0e6a9fcd..f82f1dca 100644 --- a/src/Particles/dynamicPointStructure/dynamicPointStructure.hpp +++ b/src/Particles/dynamicPointStructure/dynamicPointStructure.hpp @@ -44,6 +44,8 @@ private: uniquePtr integrationVel_ = nullptr; + Timer velocityUpdateTimer_; + /// @brief integration method for velocity and position word integrationMethod_; diff --git a/src/Particles/particles/particles.cpp b/src/Particles/particles/particles.cpp index 187f45da..a8fc970a 100644 --- a/src/Particles/particles/particles.cpp +++ b/src/Particles/particles/particles.cpp @@ -74,7 +74,8 @@ pFlow::particles::particles(systemControl& control) dynPointStruct_, zero3 ), - idHandler_(particleIdHandler::create(dynPointStruct_)) + idHandler_(particleIdHandler::create(dynPointStruct_)), + baseFieldBoundaryUpdateTimer_("baseFieldBoundaryUpdate",&timers()) { this->addToSubscriber(dynPointStruct_); @@ -84,18 +85,18 @@ pFlow::particles::particles(systemControl& control) bool pFlow::particles::beforeIteration() { - zeroForce(); - zeroTorque(); - if( !dynPointStruct_.beforeIteration()) { return false; } + zeroForce(); + zeroTorque(); + baseFieldBoundaryUpdateTimer_.start(); shapeIndex_.updateBoundariesSlaveToMasterIfRequested(); accelertion_.updateBoundariesSlaveToMasterIfRequested(); idHandler_().updateBoundariesSlaveToMasterIfRequested(); - + baseFieldBoundaryUpdateTimer_.end(); return true; } diff --git a/src/Particles/particles/particles.hpp b/src/Particles/particles/particles.hpp index 1f3c7ab1..12e3da8c 100644 --- a/src/Particles/particles/particles.hpp +++ b/src/Particles/particles/particles.hpp @@ -54,6 +54,8 @@ private: /// handling new ids for new particles uniquePtr idHandler_ = nullptr; + Timer baseFieldBoundaryUpdateTimer_; + /// messages for this objects static inline const message defaultMessage_{ message::DEFAULT }; diff --git a/src/phasicFlow/CMakeLists.txt b/src/phasicFlow/CMakeLists.txt index de526e20..408fc042 100644 --- a/src/phasicFlow/CMakeLists.txt +++ b/src/phasicFlow/CMakeLists.txt @@ -94,6 +94,8 @@ structuredData/pointStructure/selectors/selectorStridedRange/selectorStridedRang structuredData/pointStructure/selectors/selectorRandomPoints/selectorRandomPoints.cpp structuredData/pointStructure/selectors/selectBox/selectBox.cpp structuredData/pointStructure/selectors/selectorGeometric/selectorGeometrics.cpp +structuredData/pointStructure/pointStructure/pointSorting/mortonIndexing.cpp +structuredData/pointStructure/pointStructure/pointSorting/pointSorting.cpp triSurface/subSurface.cpp triSurface/triSurface.cpp diff --git a/src/phasicFlow/containers/VectorHD/VectorSingle.cpp b/src/phasicFlow/containers/VectorHD/VectorSingle.cpp index 24248594..dd093270 100644 --- a/src/phasicFlow/containers/VectorHD/VectorSingle.cpp +++ b/src/phasicFlow/containers/VectorHD/VectorSingle.cpp @@ -226,6 +226,32 @@ pFlow::VectorSingle::VectorSingle } } +template +pFlow::VectorSingle::VectorSingle +( + const word& name, + const ViewType1D& src +) +: + VectorSingle(name, src.size(), src.size(), RESERVE()) +{ + if constexpr(isTriviallyCopyable_) + { + copy(deviceView(), src); + } + else if constexpr( isHostAccessible_) + { + for(auto i=0u; i pFlow::VectorSingle& pFlow::VectorSingle::operator = (const VectorSingle& rhs) @@ -367,6 +393,7 @@ template INLINE_FUNCTION_H void pFlow::VectorSingle::reserve(uint32 cap) { + if(cap == capacity() ) return; changeCapacity(cap); } @@ -571,6 +598,40 @@ void pFlow::VectorSingle::assign } } + +template +INLINE_FUNCTION_H +void pFlow::VectorSingle::append(const ViewType1D& appVec) +{ + uint32 appSize = appVec.size(); + if(appSize == 0) return; + + uint32 oldS = size(); + uint32 newSize = oldS + appSize; + + changeSize(newSize); + + auto appendView = Kokkos::subview( + view_, + Kokkos::make_pair(oldS, newSize)); + + if constexpr( isTriviallyCopyable_) + { + copy(appendView, appVec); + } + else if constexpr( isHostAccessible_) + { + for(auto i=0u; i INLINE_FUNCTION_H void pFlow::VectorSingle::append @@ -813,7 +874,7 @@ bool pFlow::VectorSingle::insertSetElement template INLINE_FUNCTION_H -bool pFlow::VectorSingle::reorderItems(uint32IndexContainer indices) +bool pFlow::VectorSingle::reorderItems(const uint32IndexContainer& indices) { if(indices.size() == 0) { @@ -831,10 +892,8 @@ bool pFlow::VectorSingle::reorderItems(uint32IndexContainer indic return false; } - uint32 newSize = indices.size(); + uint32 newSize = indices.size(); - setSize(newSize); - viewType sortedView(this->name(), newSize); using policy = Kokkos::RangePolicy< execution_space,Kokkos::IndexType>; @@ -875,7 +934,24 @@ bool pFlow::VectorSingle::reorderItems(uint32IndexContainer indic Kokkos::fence(); } - copy(deviceView(), sortedView); - + setSize(newSize); + + + if constexpr( isTriviallyCopyable_ ) + { + copy(deviceView(), sortedView); + } + else if constexpr( isHostAccessible_) + { + for(auto i=0u; i& src); /// Copy assignment (perform deep copy from rhs to *this) VectorSingle& operator = (const VectorSingle& rhs) ; @@ -307,6 +310,9 @@ public: } } + INLINE_FUNCTION_H + void append(const ViewType1D& appVec); + INLINE_FUNCTION_H void append(const std::vector& appVec); @@ -331,7 +337,7 @@ public: const ViewType1D vals); INLINE_FUNCTION_H - bool reorderItems(uint32IndexContainer indices); + bool reorderItems(const uint32IndexContainer& indices); /// @brief push a new element at the end (host call only) /// resize if necessary and works on host accessible vector. diff --git a/src/phasicFlow/containers/indexContainer/indexContainer.hpp b/src/phasicFlow/containers/indexContainer/indexContainer.hpp index 72118609..4486cd2c 100644 --- a/src/phasicFlow/containers/indexContainer/indexContainer.hpp +++ b/src/phasicFlow/containers/indexContainer/indexContainer.hpp @@ -143,6 +143,16 @@ public: indexContainer(ind.data(), ind.size()) {} + indexContainer(const DeviceViewType& ind): + size_(ind.size()), + views_("indexContainer", size_) + { + copy(views_.h_view, ind); + copy(views_.d_view, ind); + min_ = pFlow::min(views_.d_view, 0, size_); + max_ = pFlow::max(views_.d_view, 0, size_); + } + /// Copy indexContainer(const indexContainer&) = default; @@ -240,13 +250,13 @@ public: void syncViews() { bool findMinMax = false; - if(views_.template need_sync()) + if(views_.template need_sync()) { Kokkos::deep_copy(views_.d_view, views_.h_view); views_.clear_sync_state(); findMinMax = true; } - else if(views_.template need_sync()) + else if(views_.template need_sync()) { Kokkos::deep_copy(views_.h_view, views_.d_view); views_.clear_sync_state(); @@ -260,6 +270,12 @@ public: } } + void syncViews(uint32 newSize) + { + size_ = newSize; + syncViews(); + } + }; diff --git a/src/phasicFlow/containers/pointField/boundaryField/boundaryField/boundaryField.hpp b/src/phasicFlow/containers/pointField/boundaryField/boundaryField/boundaryField.hpp index d30b3f18..9570f86b 100644 --- a/src/phasicFlow/containers/pointField/boundaryField/boundaryField/boundaryField.hpp +++ b/src/phasicFlow/containers/pointField/boundaryField/boundaryField/boundaryField.hpp @@ -121,7 +121,18 @@ public: return true; } - + inline + InternalFieldType& internal() + { + return internal_; + } + + inline + const InternalFieldType& internal()const + { + return internal_; + } + FieldAccessType thisField()const { return FieldAccessType( @@ -144,11 +155,6 @@ public: virtual const ProcVectorType& neighborProcField()const; - void fill(const std::any& val)override - { - return; - } - virtual void fill(const T& val) { diff --git a/src/phasicFlow/containers/pointField/boundaryField/generalBoundary/generalBoundary.hpp b/src/phasicFlow/containers/pointField/boundaryField/generalBoundary/generalBoundary.hpp index c55a7f71..6a1f6808 100644 --- a/src/phasicFlow/containers/pointField/boundaryField/generalBoundary/generalBoundary.hpp +++ b/src/phasicFlow/containers/pointField/boundaryField/generalBoundary/generalBoundary.hpp @@ -147,10 +147,6 @@ public: } const Time& time()const; - - virtual - void fill(const std::any& val)=0; - }; diff --git a/src/phasicFlow/containers/pointField/internalField/internalField.cpp b/src/phasicFlow/containers/pointField/internalField/internalField.cpp index c59faf3b..4b7c01e9 100644 --- a/src/phasicFlow/containers/pointField/internalField/internalField.cpp +++ b/src/phasicFlow/containers/pointField/internalField/internalField.cpp @@ -52,9 +52,29 @@ bool pFlow::internalField::insert(const anyList& varList) } return true; - } +template +bool pFlow::internalField::rearrange(const anyList& varList) +{ + const word eventName = message::eventName(message::ITEM_REARRANGE); + + const auto& indices = varList.getObject( + eventName); + + field_.reserve( internalPoints_.capacity()); + field_.resize(internalPoints_.size()); + if(!field_.reorderItems(indices)) + { + fatalErrorInFunction<< + "cannot reorder items in field "<< name()< pFlow::internalField::internalField ( @@ -166,15 +186,15 @@ bool pFlow::internalField:: hearChanges { // do nothing } - if(msg.equivalentTo(message::ITEM_REARRANGE)) - { - notImplementedFunction; - return false; - } if(msg.equivalentTo(message::ITEM_INSERT)) { return insert(varList); } + if(msg.equivalentTo(message::ITEM_REARRANGE)) + { + return rearrange(varList); + + } return true; } diff --git a/src/phasicFlow/containers/pointField/internalField/internalField.hpp b/src/phasicFlow/containers/pointField/internalField/internalField.hpp index 423f1ef7..ab26c663 100644 --- a/src/phasicFlow/containers/pointField/internalField/internalField.hpp +++ b/src/phasicFlow/containers/pointField/internalField/internalField.hpp @@ -73,6 +73,8 @@ protected: bool insert(const anyList& varList); + bool rearrange(const anyList& varList); + public: internalField( diff --git a/src/phasicFlow/eventManagement/message.hpp b/src/phasicFlow/eventManagement/message.hpp index 3ece7411..0593221a 100644 --- a/src/phasicFlow/eventManagement/message.hpp +++ b/src/phasicFlow/eventManagement/message.hpp @@ -48,14 +48,16 @@ public: BNDR_RESET = 10, // boundary indices reset entirely BNDR_DELETE = 11, // boundary indices deleted BNDR_APPEND = 12, // - BNDR_PROCTRANS1 = 13, // transfer of data between processors step 1 - BNDR_PROCTRANS2 = 14 // transfer of data between processors step 2 + BNDR_PROCTRANSFER_SEND = 13, // transfer of data between processors step 1 + BNDR_PROCTRANSFER_RECIEVE = 14, // transfer of data between processors step 2 + BNDR_PROCTRANSFER_WAITFILL = 15, // wait for data and fill the field + BNDR_PROC_SIZE_CHANGED = 16 }; private: - static constexpr size_t numberOfEvents_ = 15; + static constexpr size_t numberOfEvents_ = 17; std::bitset events_{0x0000}; @@ -76,7 +78,9 @@ private: "deletedIndices", "appendedIndices", "transferredIndices", - "numberRecieved" + "numToRecieve", + "insertedIndices", + "size" }; public: diff --git a/src/phasicFlow/repository/Time/baseTimeControl.cpp b/src/phasicFlow/repository/Time/baseTimeControl.cpp index 8f4f2ef2..78d04225 100644 --- a/src/phasicFlow/repository/Time/baseTimeControl.cpp +++ b/src/phasicFlow/repository/Time/baseTimeControl.cpp @@ -65,6 +65,16 @@ pFlow::baseTimeControl::baseTimeControl } +pFlow::baseTimeControl::baseTimeControl(int32 start, int32 end, int32 stride, const word &intervalPrefix) +: + isTimeStep_(true), + iRange_(start, end, max(stride,1)), + intervalPrefix_( + intervalPrefix.size()==0uL? word("interval"): intervalPrefix+"Interval" + ) +{ +} + bool pFlow::baseTimeControl::timeEvent(uint32 iter, real t, real dt) const { if(isTimeStep_) diff --git a/src/phasicFlow/repository/Time/baseTimeControl.hpp b/src/phasicFlow/repository/Time/baseTimeControl.hpp index 192973fc..5a7634d5 100644 --- a/src/phasicFlow/repository/Time/baseTimeControl.hpp +++ b/src/phasicFlow/repository/Time/baseTimeControl.hpp @@ -46,6 +46,13 @@ public: real defStartTime = 0.0 ); + baseTimeControl( + int32 start, + int32 end, + int32 stride, + const word& intervalPrefix = "" + ); + inline bool isTimeStep() const { return isTimeStep_; diff --git a/src/phasicFlow/streams/dataIO/dataIO.cpp b/src/phasicFlow/streams/dataIO/dataIO.cpp index 39fbd9e0..851ec4a5 100644 --- a/src/phasicFlow/streams/dataIO/dataIO.cpp +++ b/src/phasicFlow/streams/dataIO/dataIO.cpp @@ -22,7 +22,7 @@ template bool pFlow::writeDataAsciiBinary(iOstream& os, span data) { - if( os.isBinary() ) + if( std::is_trivially_copyable_v && os.isBinary() ) { // first write the number of data uint64 numData = data.size(); @@ -83,7 +83,7 @@ bool pFlow::readDataAsciiBinary ) { - if(is.isBinary()) + if(std::is_trivially_copyable_v && is.isBinary()) { data.clear(); // read length of data diff --git a/src/phasicFlow/structuredData/boundaries/boundaryBase/boundaryBase.cpp b/src/phasicFlow/structuredData/boundaries/boundaryBase/boundaryBase.cpp index f8307478..906d5fd9 100644 --- a/src/phasicFlow/structuredData/boundaries/boundaryBase/boundaryBase.cpp +++ b/src/phasicFlow/structuredData/boundaries/boundaryBase/boundaryBase.cpp @@ -42,6 +42,7 @@ void pFlow::boundaryBase::setNewIndices unSyncLists(); } + bool pFlow::boundaryBase::appendNewIndices ( const uint32Vector_D& newIndices @@ -144,7 +145,7 @@ bool pFlow::boundaryBase::setRemoveKeepIndices return true; } -bool pFlow::boundaryBase::transferPoints +bool pFlow::boundaryBase::transferPointsToMirror ( uint32 numTransfer, const uint32Vector_D& transferMask, @@ -203,6 +204,8 @@ pFlow::boundaryBase::boundaryBase indexList_(groupNames("indexList", dict.name())), indexListHost_(groupNames("hostIndexList",dict.name())), neighborLength_(dict.getVal("neighborLength")), + updateInetrval_(dict.getVal("updateInterval")), + boundaryExtntionLengthRatio_(dict.getVal("boundaryExtntionLengthRatio")), internal_(internal), boundaries_(bndrs), thisBoundaryIndex_(thisIndex), diff --git a/src/phasicFlow/structuredData/boundaries/boundaryBase/boundaryBase.hpp b/src/phasicFlow/structuredData/boundaries/boundaryBase/boundaryBase.hpp index 9b63c9df..978f674b 100644 --- a/src/phasicFlow/structuredData/boundaries/boundaryBase/boundaryBase.hpp +++ b/src/phasicFlow/structuredData/boundaries/boundaryBase/boundaryBase.hpp @@ -68,6 +68,12 @@ private: /// The length defined for creating neighbor list real neighborLength_; + /// time steps between successive update of boundary lists + uint32 updateInetrval_; + + /// the extra boundary extension beyound actual limits of boundary + real boundaryExtntionLengthRatio_; + /// a reference to internal points internalPoints& internal_; @@ -106,7 +112,7 @@ protected: const uint32Vector_D& removeIndices, const uint32Vector_D& keepIndices); - bool transferPoints( + bool transferPointsToMirror( uint32 numTransfer, const uint32Vector_D& transferMask, uint32 transferBoundaryIndex, @@ -126,6 +132,12 @@ protected: } } + /// Is this iter the right time for updating bounday list + bool boundaryListUpdate(uint32 iter)const + { + return iter%updateInetrval_ == 0u || iter == 0u; + } + /// Update this boundary data in two steps (1 and 2). /// This is called after calling beforeIteration for /// all boundaries, so any particle addition, deletion, @@ -133,7 +145,7 @@ protected: /// This two-step update help to have a flexible mechanism /// for data transfer, mostly for MPI related jobs. virtual - bool updataBoundary(int step) + bool updataBoundaryData(int step) { return true; } @@ -145,11 +157,13 @@ protected: /// to complete. false: if the operation is complete and no need for /// additional step in operation. virtual - bool transferData(int step) + bool transferData(uint32 iter, int step) { return false; } + + public: TypeInfo("boundaryBase"); @@ -189,20 +203,21 @@ public: /// The length from boundary plane into the domain /// where beyond that distance internal points exist. /// By conventions is it always equal to neighborLength_ + inline real neighborLengthIntoInternal()const { return neighborLength_; } /// The distance length from boundary plane - /// where neighbor particles exist in that distance. + /// where neighbor particles still exist in that distance. /// This length may be modified in each boundary type /// as required. In this case the boundaryExtensionLength /// method should also be modified accordingly. virtual real neighborLength()const { - return neighborLength_; + return (1+boundaryExtntionLengthRatio_)*neighborLength_; } /// The extention length (in vector form) for the boundary @@ -213,7 +228,7 @@ public: virtual realx3 boundaryExtensionLength()const { - return {0,0,0}; + return -boundaryExtntionLengthRatio_*neighborLength_ * boundaryPlane_.normal(); } inline @@ -351,6 +366,18 @@ public: virtual const realx3Vector_D& neighborProcPoints()const; + virtual + uint32 numToTransfer()const + { + return 0u; + } + + virtual + uint32 numToRecieve()const + { + return 0u; + } + /// - static create static uniquePtr create diff --git a/src/phasicFlow/structuredData/boundaries/boundaryBase/scatteredFieldAccess.hpp b/src/phasicFlow/structuredData/boundaries/boundaryBase/scatteredFieldAccess.hpp index f26438db..e051c803 100644 --- a/src/phasicFlow/structuredData/boundaries/boundaryBase/scatteredFieldAccess.hpp +++ b/src/phasicFlow/structuredData/boundaries/boundaryBase/scatteredFieldAccess.hpp @@ -114,6 +114,16 @@ public: return fieldVals_; } + auto& indices() + { + return indices_; + } + + const auto& indices()const + { + return indices_; + } + uint32 index(uint32 i)const { return indices_[i]; diff --git a/src/phasicFlow/structuredData/boundaries/boundaryExit/boundaryExit.cpp b/src/phasicFlow/structuredData/boundaries/boundaryExit/boundaryExit.cpp index 3b753877..9937433c 100644 --- a/src/phasicFlow/structuredData/boundaries/boundaryExit/boundaryExit.cpp +++ b/src/phasicFlow/structuredData/boundaries/boundaryExit/boundaryExit.cpp @@ -47,6 +47,9 @@ bool pFlow::boundaryExit::beforeIteration real dt ) { + + if( !boundaryListUpdate(iterNum))return true; + // nothing have to be done if(empty()) { diff --git a/src/phasicFlow/structuredData/boundaries/boundaryList.cpp b/src/phasicFlow/structuredData/boundaries/boundaryList.cpp index 6cccb194..473e523d 100644 --- a/src/phasicFlow/structuredData/boundaries/boundaryList.cpp +++ b/src/phasicFlow/structuredData/boundaries/boundaryList.cpp @@ -91,18 +91,18 @@ pFlow::boundaryList::updateNeighborLists() pFlow::boundaryList::boundaryList(pointStructure& pStruct) : ListPtr(pStruct.simDomain().sizeOfBoundaries()), pStruct_(pStruct), - timeControl_( - pStruct.simDomain().subDict("boundaries"), - "update", - pStruct.currentTime() - ) + neighborListUpdateInterval_( + pStruct.simDomain().subDict("boundaries").getVal( + "neighborListUpdateInterval" + ) + ) { } bool -pFlow::boundaryList::updateNeighborLists(uint32 iter, real t, real dt) +pFlow::boundaryList::updateNeighborLists(uint32 iter, bool force) { - if (timeControl_.timeEvent(iter, t, dt)) + if(iter%neighborListUpdateInterval_==0u || iter == 0u || force) { return updateNeighborLists(); } @@ -155,15 +155,18 @@ pFlow::boundaryList::internalDomainBox() const } bool -pFlow::boundaryList::beforeIteration(uint32 iter, real t, real dt) +pFlow::boundaryList::beforeIteration(uint32 iter, real t, real dt, bool force) { // it is time to update lists - if (timeControl_.timeEvent(iter, t, dt) && !updateNeighborLists()) + if(!updateNeighborLists(iter , force) ) { fatalErrorInFunction; return false; } + // this forces performing updating the list on each boundary + if(force) iter = 0; + for (auto bdry : *this) { if (!bdry->beforeIteration(iter, t, dt)) @@ -176,12 +179,12 @@ pFlow::boundaryList::beforeIteration(uint32 iter, real t, real dt) for (auto bdry : *this) { - bdry->updataBoundary(1); + bdry->updataBoundaryData(1); } for (auto bdry : *this) { - bdry->updataBoundary(2); + bdry->updataBoundaryData(2); } return true; @@ -215,20 +218,11 @@ pFlow::boundaryList::afterIteration(uint32 iter, real t, real dt) { if(requireStep[i]) { - requireStep[i] = this->operator[](i).transferData(step); + requireStep[i] = this->operator[](i).transferData(iter, step); } } step++; } - /*for (auto& bdry : *this) - { - if (!bdry->afterIteration(iter, t, dt)) - { - fatalErrorInFunction << "Error in afterIteration in boundary " - << bdry->name() << endl; - return false; - } - }*/ return true; } diff --git a/src/phasicFlow/structuredData/boundaries/boundaryList.hpp b/src/phasicFlow/structuredData/boundaries/boundaryList.hpp index 5eba6399..8d4c8710 100644 --- a/src/phasicFlow/structuredData/boundaries/boundaryList.hpp +++ b/src/phasicFlow/structuredData/boundaries/boundaryList.hpp @@ -17,15 +17,12 @@ Licence: implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. -----------------------------------------------------------------------------*/ - - #ifndef __boundaryList_hpp__ #define __boundaryList_hpp__ #include "domain.hpp" #include "boundaryBase.hpp" #include "ListPtr.hpp" -#include "baseTimeControl.hpp" namespace pFlow @@ -42,7 +39,7 @@ private: //// - data members pointStructure& pStruct_; - baseTimeControl timeControl_; + uint32 neighborListUpdateInterval_; domain extendedDomain_; @@ -72,7 +69,7 @@ public: /// @brief update neighbor list of boundaries based on /// the time intervals - bool updateNeighborLists(uint32 iter, real t, real dt); + bool updateNeighborLists(uint32 iter, bool force = false); bool createBoundaries(); @@ -94,11 +91,6 @@ public: return ListPtr::operator[](i); } - inline - const baseTimeControl& timeControl()const - { - return timeControl_; - } inline const auto& extendedDomain()const @@ -114,7 +106,7 @@ public: box internalDomainBox()const; - bool beforeIteration(uint32 iter, real t, real dt); + bool beforeIteration(uint32 iter, real t, real dt, bool force = false); bool iterate(uint32 iter, real t, real dt); diff --git a/src/phasicFlow/structuredData/boundaries/boundaryNone/boundaryNone.hpp b/src/phasicFlow/structuredData/boundaries/boundaryNone/boundaryNone.hpp index 372a1c28..8e103d39 100644 --- a/src/phasicFlow/structuredData/boundaries/boundaryNone/boundaryNone.hpp +++ b/src/phasicFlow/structuredData/boundaries/boundaryNone/boundaryNone.hpp @@ -43,7 +43,7 @@ public: boundaryList &bndrs, uint32 thisIndex); - ~boundaryNone() override= default; + ~boundaryNone() final= default; add_vCtor ( @@ -52,11 +52,11 @@ public: dictionary ); - bool beforeIteration(uint32 iterNum, real t, real dt) override; + bool beforeIteration(uint32 iterNum, real t, real dt) final; - bool iterate(uint32 iterNum, real t, real dt) override; + bool iterate(uint32 iterNum, real t, real dt) final; - bool afterIteration(uint32 iterNum, real t, real dt) override; + bool afterIteration(uint32 iterNum, real t, real dt) final; }; diff --git a/src/phasicFlow/structuredData/boundaries/boundaryPeriodic/boundaryPeriodic.cpp b/src/phasicFlow/structuredData/boundaries/boundaryPeriodic/boundaryPeriodic.cpp index 30fdc039..cbdecf71 100644 --- a/src/phasicFlow/structuredData/boundaries/boundaryPeriodic/boundaryPeriodic.cpp +++ b/src/phasicFlow/structuredData/boundaries/boundaryPeriodic/boundaryPeriodic.cpp @@ -33,17 +33,20 @@ pFlow::boundaryPeriodic::boundaryPeriodic ) : boundaryBase(dict, bplane, internal, bndrs, thisIndex), - mirrorBoundaryIndex_(dict.getVal("mirrorBoundaryIndex")) -{} + mirrorBoundaryIndex_(dict.getVal("mirrorBoundaryIndex")), + extensionLength_(dict.getVal("boundaryExtntionLengthRatio")) +{ + extensionLength_ = max(extensionLength_, static_cast(0.1)); +} pFlow::real pFlow::boundaryPeriodic::neighborLength() const { - return (1+extensionLength_)*boundaryBase::neighborLength(); + return (1+extensionLength_)*neighborLengthIntoInternal(); } pFlow::realx3 pFlow::boundaryPeriodic::boundaryExtensionLength() const { - return -extensionLength_*neighborLength()*boundaryBase::boundaryPlane().normal(); + return -extensionLength_*neighborLengthIntoInternal()*boundaryBase::boundaryPlane().normal(); } @@ -57,6 +60,8 @@ bool pFlow::boundaryPeriodic::beforeIteration( { return true; } + + if( !boundaryListUpdate(iterNum))return true; uint32 s = size(); uint32Vector_D transferFlags("transferFlags",s+1, s+1, RESERVE()); @@ -92,7 +97,7 @@ bool pFlow::boundaryPeriodic::beforeIteration( // to obtain the transfer vector realx3 transferVec = displacementVectroToMirror(); - return transferPoints + return transferPointsToMirror ( numTransfered, transferFlags, diff --git a/src/phasicFlow/structuredData/boundaries/boundaryPeriodic/boundaryPeriodic.hpp b/src/phasicFlow/structuredData/boundaries/boundaryPeriodic/boundaryPeriodic.hpp index b4746a5e..88983fa3 100644 --- a/src/phasicFlow/structuredData/boundaries/boundaryPeriodic/boundaryPeriodic.hpp +++ b/src/phasicFlow/structuredData/boundaries/boundaryPeriodic/boundaryPeriodic.hpp @@ -35,8 +35,7 @@ private: uint32 mirrorBoundaryIndex_; - static - inline const real extensionLength_ = 0.1; + real extensionLength_ = 0.1; public: diff --git a/src/phasicFlow/structuredData/cells/cells.hpp b/src/phasicFlow/structuredData/cells/cells.hpp index c77b1a69..c4b50813 100644 --- a/src/phasicFlow/structuredData/cells/cells.hpp +++ b/src/phasicFlow/structuredData/cells/cells.hpp @@ -28,58 +28,40 @@ Licence: namespace pFlow { -template class cells { -public: - - using CellType = triple; - -protected: +private: // - domain - box domain_{realx3(0.0), realx3(1.0)}; + box domainBox_{realx3(0.0), realx3(1.0)}; // - cell size - realx3 cellSize_{1,1,1}; - - CellType numCells_{1,1,1}; + real celldx_{1}; + int32x3 numCells_{1,1,1}; // - protected methods INLINE_FUNCTION_H void calculate() { - numCells_ = (domain_.maxPoint()-domain_.minPoint())/cellSize_ + realx3(1.0); - numCells_ = max( numCells_ , CellType(static_cast(1)) ); + numCells_ = (domainBox_.maxPoint()-domainBox_.minPoint())/celldx_ + realx3(1.0); + numCells_ = max( numCells_ , int32x3(1) ); } public: INLINE_FUNCTION_HD - cells() - {} + cells() = default; INLINE_FUNCTION_H cells(const box& domain, real cellSize) : - domain_(domain), - cellSize_(cellSize) + domainBox_(domain), + celldx_(cellSize) { calculate(); } - - INLINE_FUNCTION_H - cells(const box& domain, int32 nx, int32 ny, int32 nz) - : - domain_(domain), - cellSize_( - (domain_.maxPoint() - domain_.minPoint())/realx3(nx, ny, nz) - ), - numCells_(nx, ny, nz) - {} - INLINE_FUNCTION_HD cells(const cells&) = default; @@ -100,43 +82,36 @@ public: INLINE_FUNCTION_H void setCellSize(real cellSize) { - cellSize_ = cellSize; - calculate(); - } - - INLINE_FUNCTION_H - void setCellSize(realx3 cellSize) - { - cellSize_ = cellSize; + celldx_ = cellSize; calculate(); } INLINE_FUNCTION_HD - realx3 cellSize()const + real cellSize()const { - return cellSize_; + return celldx_; } INLINE_FUNCTION_HD - const CellType& numCells()const + const int32x3& numCells()const { return numCells_; } INLINE_FUNCTION_HD - indexType nx()const + int32 nx()const { return numCells_.x(); } INLINE_FUNCTION_HD - indexType ny()const + int32 ny()const { return numCells_.y(); } INLINE_FUNCTION_HD - indexType nz()const + int32 nz()const { return numCells_.z(); } @@ -149,22 +124,21 @@ public: static_cast(numCells_.z()); } - const auto& domain()const + const auto& domainBox()const { - return domain_; + return domainBox_; } INLINE_FUNCTION_HD - CellType pointIndex(const realx3& p)const + int32x3 pointIndex(const realx3& p)const { - return CellType( (p - domain_.minPoint())/cellSize_ ); + return int32x3( (p - domainBox_.minPoint())/celldx_ ); } INLINE_FUNCTION_HD - bool pointIndexInDomain(const realx3 p, CellType& index)const + bool pointIndexInDomain(const realx3 p, int32x3& index)const { - if( !domain_.isInside(p) ) return false; - + if(!inDomain(p))return false; index = this->pointIndex(p); return true; } @@ -172,11 +146,11 @@ public: INLINE_FUNCTION_HD bool inDomain(const realx3& p)const { - return domain_.isInside(p); + return domainBox_.isInside(p); } INLINE_FUNCTION_HD - bool isInRange(const CellType& cell)const + bool inCellRange(const int32x3& cell)const { if(cell.x()<0)return false; if(cell.y()<0)return false; @@ -188,7 +162,7 @@ public: } INLINE_FUNCTION_HD - bool isInRange(indexType i, indexType j, indexType k)const + bool inCellRange(int32 i, int32 j, int32 k)const { if(i<0)return false; if(j<0)return false; @@ -199,22 +173,7 @@ public: return true; } - INLINE_FUNCTION_HD - void extendBox( - const CellType& p1, - const CellType& p2, - const CellType& p3, - indexType extent, - CellType& minP, - CellType& maxP)const - { - minP = min( min( p1, p2), p3)-extent; - maxP = max( max( p1, p2), p3)+extent; - - minP = bound(minP); - maxP = bound(maxP); - } - + INLINE_FUNCTION_HD void extendBox( const realx3& p1, @@ -224,17 +183,17 @@ public: realx3& minP, realx3& maxP)const { - minP = min(min(p1,p2),p3) - extent*cellSize_ ; - maxP = max(max(p1,p2),p3) + extent*cellSize_ ; + minP = min(min(p1,p2),p3) - extent*celldx_ ; + maxP = max(max(p1,p2),p3) + extent*celldx_ ; minP = bound(minP); maxP = bound(maxP); } INLINE_FUNCTION_HD - CellType bound(CellType p)const + int32x3 bound(int32x3 p)const { - return CellType( + return int32x3( min( numCells_.x()-1, max(0,p.x())), min( numCells_.y()-1, max(0,p.y())), min( numCells_.z()-1, max(0,p.z())) @@ -245,9 +204,9 @@ public: realx3 bound(realx3 p)const { return realx3( - min( domain_.maxPoint().x(), max(domain_.minPoint().x(),p.x())), - min( domain_.maxPoint().y(), max(domain_.minPoint().y(),p.y())), - min( domain_.maxPoint().z(), max(domain_.minPoint().z(),p.z())) + min( domainBox_.maxPoint().x(), max(domainBox_.minPoint().x(),p.x())), + min( domainBox_.maxPoint().y(), max(domainBox_.minPoint().y(),p.y())), + min( domainBox_.maxPoint().z(), max(domainBox_.minPoint().z(),p.z())) ); } }; diff --git a/src/phasicFlow/structuredData/domain/regularSimulationDomain.cpp b/src/phasicFlow/structuredData/domain/regularSimulationDomain.cpp index e94d111e..a84325a9 100644 --- a/src/phasicFlow/structuredData/domain/regularSimulationDomain.cpp +++ b/src/phasicFlow/structuredData/domain/regularSimulationDomain.cpp @@ -31,7 +31,10 @@ bool pFlow::regularSimulationDomain::createBoundaryDicts() this->addDict("regularBoundaries", boundaries); auto& rbBoundaries = this->subDict("regularBoundaries"); - real neighborLength = boundaries.getVal("neighborLength"); + auto neighborLength = boundaries.getVal("neighborLength"); + auto boundaryExtntionLengthRatio = + boundaries.getValOrSet("boundaryExtntionLengthRatio", 0.1); + auto updateIntercal = boundaries.getValOrSet("updateInterval", 1u); for(uint32 i=0; i(points.size()); + + auto aRange = pFlagsD_.activeRange(); + uint32 emptySpots = pFlagsD_.capacity() - pFlagsD_.numActive(); + + if(emptySpots!= 0) emptySpots--; + + if( numNew > emptySpots ) + { + // increase the capacity to hold new points + aRange = pFlagsD_.activeRange(); + uint32 newCap = pFlagsD_.changeCapacity(numNew); + unSyncFlag(); + varList.emplaceBack( + msg.addAndName(message::CAP_CHANGED), + newCap); + } + + + // first check if it is possible to add to the beggining of list + if(numNew <= aRange.start()) + { + varList.emplaceBack( + msg.addAndName(message::ITEM_INSERT), + 0u, numNew); + }// check if it is possible to add to the end of the list + else if( numNew <= pFlagsD_.capacity() - aRange.end() ) + { + varList.emplaceBack( + msg.addAndName(message::ITEM_INSERT), + aRange.end(), aRange.end()+numNew); + }// we should fill the scattered empty spots + else + { + pOutput<<"numNew to be inserted "<< numNew <( + msg.addAndName(message::ITEM_INSERT), + newIndices + ); + } + + const auto& indices = varList.getObject( + message::eventName(message::ITEM_INSERT) + ); + + auto nAdded = pFlagsD_.addInternalPoints(indices.deviceView()); + unSyncFlag(); + + if(nAdded != numNew ) + { + fatalErrorInFunction; + return false; + } + + pointPosition_.reserve( pFlagsD_.capacity() ); + + if(!pointPosition_.insertSetElement(indices, points.deviceView())) + { + fatalErrorInFunction<< + "Error in inserting new positions into pointPosition"<< + " internalPoints field"<(Flag::DELETED)); + } + else + { + fill(flags_, end, cap, static_cast(Flag::DELETED)); + } + + activeRange_ = {start, end}; + numActive_ = activeRange_.numElements(); + isAllActive_ = true; + fill(flags_, activeRange_, static_cast(Flag::INTERNAL)); + + nLeft_ = nRight_ = 0; + nBottom_ = nTop_ = 0; + nRear_ = nFront_ = 0; + } + public: - friend class internalPoints; pointFlag() = default; @@ -283,6 +305,9 @@ public: } ViewType1D getActivePoints(); + + + ViewType1D getEmptyPoints(uint32 numToGet)const; /// @brief Loop over the active points and mark those out of the box diff --git a/src/phasicFlow/structuredData/pointStructure/internalPoints/pointFlagKernels.hpp b/src/phasicFlow/structuredData/pointStructure/internalPoints/pointFlagKernels.hpp index 91764cf6..de8940ba 100644 --- a/src/phasicFlow/structuredData/pointStructure/internalPoints/pointFlagKernels.hpp +++ b/src/phasicFlow/structuredData/pointStructure/internalPoints/pointFlagKernels.hpp @@ -50,6 +50,51 @@ pFlow::ViewType1D::memo return aPoints; } +template +pFlow::ViewType1D::memory_space> + pFlow::pointFlag::getEmptyPoints(uint32 numToGet)const +{ + + uint32 cap = capacity(); + using rpAPoints = Kokkos::RangePolicy>; + + ViewType1D indices("indices", cap+1); + + ViewType1D emptyPoints("emptyPoints", numToGet); + + + Kokkos::parallel_for( + "getEmptyPoints", + rPolicy(0, cap), + LAMBDA_HD(uint32 i){ + indices(i) = flags_[i] == DELETED; + }); + Kokkos::fence(); + + exclusiveScan(indices, 0u, cap, indices, 0u); + uint32 numFound = 0; + Kokkos::parallel_reduce( + "fillEmptyPoints", + rpAPoints(0, cap-1), + LAMBDA_HD(uint32 i, uint32& numFoundUpdate){ + if(indices(i)!= indices(i+1) && indices(i)< numToGet) + { + emptyPoints(indices(i)) = i; + numFoundUpdate++; + } + }, + numFound + ); + + if(numFound < numToGet) + { + return Kokkos::subview(emptyPoints, Kokkos::make_pair(0u, numFound)); + } + return emptyPoints; + +} + template pFlow::uint32 pFlow::pointFlag::markOutOfBoxDelete diff --git a/src/phasicFlow/structuredData/pointStructure/pointStructure/mortonIndexing.cpp b/src/phasicFlow/structuredData/pointStructure/pointStructure/pointSorting/mortonIndexing.cpp similarity index 59% rename from src/phasicFlow/structuredData/pointStructure/pointStructure/mortonIndexing.cpp rename to src/phasicFlow/structuredData/pointStructure/pointStructure/pointSorting/mortonIndexing.cpp index 63969747..2405bca4 100644 --- a/src/phasicFlow/structuredData/pointStructure/pointStructure/mortonIndexing.cpp +++ b/src/phasicFlow/structuredData/pointStructure/pointStructure/pointSorting/mortonIndexing.cpp @@ -17,66 +17,60 @@ Licence: implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. -----------------------------------------------------------------------------*/ - #include "mortonIndexing.hpp" #include "cells.hpp" -#include "streams.hpp" -bool pFlow::getSortedIndex( +pFlow::uint32IndexContainer pFlow::getSortedIndices( box boundingBox, - real dx, - range activeRange, - ViewType1D pos, - ViewType1D flag, - int32IndexContainer& sortedIndex) + real dx, + const ViewType1D& pos, + const pFlagTypeDevice& flag +) { + if(flag.numActive() == 0u)return uint32IndexContainer(); // obtain the morton code of the particles - cells allCells( boundingBox, dx); - int32IndexContainer index(activeRange.first, activeRange.second); + cells allCells( boundingBox, dx); + auto aRange = flag.activeRange(); - ViewType1D mortonCode("mortonCode", activeRange.second); + uint32IndexContainer sortedIndex(aRange.start(), aRange.end()); - output<<"before first kernel"< mortonCode("mortonCode", aRange.end()); - using rpMorton = - Kokkos::RangePolicy>; - int32 numActive = 0; - Kokkos::parallel_reduce + + Kokkos::parallel_for ( "mortonIndexing::getIndex::morton", - rpMorton(activeRange.first, activeRange.second), - LAMBDA_HD(int32 i, int32& sumToUpdate){ - if( flag[i] == 1 ) // active point + deviceRPolicyStatic(aRange.start(), aRange.end()), + LAMBDA_HD(uint32 i){ + if( flag.isActive(i)) // active point { auto cellInd = allCells.pointIndex(pos[i]); mortonCode[i] = xyzToMortonCode64(cellInd.x(), cellInd.y(), cellInd.z()); - sumToUpdate++; }else { mortonCode[i] = xyzToMortonCode64 ( - static_cast(-1), - static_cast(-1), - static_cast(-1) + largestPosInt32, + largestPosInt32, + largestPosInt32 ); } - }, - numActive + } ); + Kokkos::fence(); + permuteSort( mortonCode, - activeRange.first, - activeRange.second, - index.deviceView(), + aRange.start(), + aRange.end(), + sortedIndex.deviceView(), 0 ); - index.modifyOnDevice(); - index.setSize(numActive); - index.syncViews(); - sortedIndex = index; + sortedIndex.modifyOnDevice(); + sortedIndex.syncViews(flag.numActive()); - return true; + return sortedIndex; } \ No newline at end of file diff --git a/src/phasicFlow/structuredData/pointStructure/pointStructure/mortonIndexing.hpp b/src/phasicFlow/structuredData/pointStructure/pointStructure/pointSorting/mortonIndexing.hpp similarity index 94% rename from src/phasicFlow/structuredData/pointStructure/pointStructure/mortonIndexing.hpp rename to src/phasicFlow/structuredData/pointStructure/pointStructure/pointSorting/mortonIndexing.hpp index 9ccdc83e..0d8c58d7 100644 --- a/src/phasicFlow/structuredData/pointStructure/pointStructure/mortonIndexing.hpp +++ b/src/phasicFlow/structuredData/pointStructure/pointStructure/pointSorting/mortonIndexing.hpp @@ -17,24 +17,22 @@ Licence: implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. -----------------------------------------------------------------------------*/ - #ifndef __mortonIndexing_hpp__ #define __mortonIndexing_hpp__ #include "types.hpp" #include "box.hpp" #include "indexContainer.hpp" +#include "pointFlag.hpp" namespace pFlow { -bool getSortedIndex( +uint32IndexContainer getSortedIndices( box boundingBox, - real dx, - range activeRange, - ViewType1D pos, - ViewType1D flag, - int32IndexContainer& sortedIndex); + real dx, + const ViewType1D& pos, + const pFlagTypeDevice& flag); INLINE_FUNCTION_HD diff --git a/src/phasicFlow/structuredData/pointStructure/pointStructure/pointSorting/pointSorting.cpp b/src/phasicFlow/structuredData/pointStructure/pointStructure/pointSorting/pointSorting.cpp new file mode 100644 index 00000000..c8dda774 --- /dev/null +++ b/src/phasicFlow/structuredData/pointStructure/pointStructure/pointSorting/pointSorting.cpp @@ -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 "pointSorting.hpp" +#include "mortonIndexing.hpp" + +pFlow::pointSorting::pointSorting(const dictionary & dict) +: + performSorting_(dict.getValOrSet("active", Logical(false))), + timeControl_( + performSorting_()? + baseTimeControl(dict, "sorting"): + baseTimeControl(0,1,1, "sorting") + ), + dx_( + performSorting_()? + dict.getVal("dx"): + 1.0 + ) +{ + if( performSorting_() ) + REPORT(1)<<"Point sorting is "< &pos, + const pFlagTypeDevice &flag +) const +{ + return pFlow::getSortedIndices(boundingBox, dx_, pos, flag); +} diff --git a/src/phasicFlow/structuredData/pointStructure/pointStructure/pointSorting/pointSorting.hpp b/src/phasicFlow/structuredData/pointStructure/pointStructure/pointSorting/pointSorting.hpp new file mode 100644 index 00000000..d8d9f195 --- /dev/null +++ b/src/phasicFlow/structuredData/pointStructure/pointStructure/pointSorting/pointSorting.hpp @@ -0,0 +1,68 @@ +/*------------------------------- 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 __pointSorting_hpp__ +#define __pointSorting_hpp__ + +#include "baseTimeControl.hpp" +#include "indexContainer.hpp" +#include "pointFlag.hpp" + + +namespace pFlow +{ + +class box; + + +class pointSorting +{ +private: + + Logical performSorting_; + + baseTimeControl timeControl_; + + real dx_; + +public: + + explicit pointSorting(const dictionary& dict); + + bool performSorting()const + { + return performSorting_(); + } + + bool sortTime(uint32 iter, real t, real dt)const + { + return performSorting_() && timeControl_.timeEvent(iter, t, dt); + } + + uint32IndexContainer getSortedIndices( + const box& boundingBox, + const ViewType1D& pos, + const pFlagTypeDevice& flag + )const; + +}; + +} + +#endif \ No newline at end of file diff --git a/src/phasicFlow/structuredData/pointStructure/pointStructure/pointStructure.cpp b/src/phasicFlow/structuredData/pointStructure/pointStructure/pointStructure.cpp index 9a20395c..bc0e7976 100644 --- a/src/phasicFlow/structuredData/pointStructure/pointStructure/pointStructure.cpp +++ b/src/phasicFlow/structuredData/pointStructure/pointStructure/pointStructure.cpp @@ -22,6 +22,7 @@ Licence: #include "pointStructure.hpp" #include "systemControl.hpp" #include "vocabs.hpp" +#include "anyList.hpp" bool pFlow::pointStructure::setupPointStructure(const realx3Vector& points) { @@ -73,6 +74,7 @@ bool pFlow::pointStructure::setupPointStructure(const PointsTypeHost &points) } boundaries_.createBoundaries(); + boundaries_.updateNeighborLists(0u, true); return true; } @@ -111,10 +113,13 @@ pFlow::pointStructure::pointStructure ( simulationDomain::create(control) ), + pointSorting_(simulationDomain_->subDictOrCreate("pointSorting")), boundaries_ ( *this - ) + ), + boundaryUpdateTimer_("boundaryUpdate", &timers()), + boundaryDataTransferTimer_("boundaryDataTransferTimer", &timers()) { REPORT(1)<< "Reading "<< Green_Text("pointStructure")<< " from "<subDictOrCreate("pointSorting")), boundaries_ ( *this @@ -165,12 +171,58 @@ pFlow::pointStructure::pointStructure( bool pFlow::pointStructure::beforeIteration() { - if( !boundaries_.beforeIteration(currentIter(), currentTime(), dt()) ) + uint32 iter = currentIter(); + real t = currentTime(); + real deltat = dt(); + + if(pointSorting_.sortTime(iter, t, deltat)) { - fatalErrorInFunction<< - "Unable to perform beforeIteration for boundaries"< simulationDomain_ = nullptr; + pointSorting pointSorting_; + boundaryList boundaries_; - + + Timer boundaryUpdateTimer_; + + Timer boundaryDataTransferTimer_; + + + bool setupPointStructure(const realx3Vector& points); bool setupPointStructure(const PointsTypeHost& points);