diff --git a/src/Interaction/sphereInteraction/sphereInteraction/sphereInteraction.cpp b/src/Interaction/sphereInteraction/sphereInteraction/sphereInteraction.cpp index c757c03be..7c433c15c 100644 --- a/src/Interaction/sphereInteraction/sphereInteraction/sphereInteraction.cpp +++ b/src/Interaction/sphereInteraction/sphereInteraction/sphereInteraction.cpp @@ -31,7 +31,8 @@ bool pFlow::sphereInteraction::createSphereInteraction() rhoD.deviceView(), modelDict ); - + const_cast(sphParticles_).initializeForceChain(modelDict); + uint32 nPrtcl = sphParticles_.size(); contactSearch_ = contactSearch::create( @@ -54,6 +55,11 @@ bool pFlow::sphereInteraction::sphereSphereInteraction() { auto lastItem = ppContactList_().loopCount(); + if (sphParticles_.isForceChainActive()) + { + const_cast(sphParticles_.getForceChain()).resetPairCounter(); + } + // create the kernel functor pFlow::sphereInteractionKernels::ppInteractionFunctor ppInteraction( @@ -66,7 +72,9 @@ bool pFlow::sphereInteraction::sphereSphereInteraction() sphParticles_.velocity().deviceViewAll(), sphParticles_.rVelocity().deviceViewAll(), sphParticles_.contactForce().deviceViewAll(), - sphParticles_.contactTorque().deviceViewAll() + sphParticles_.contactTorque().deviceViewAll(), + sphParticles_.isForceChainActive() ? + const_cast(&sphParticles_.getForceChain()) : nullptr ); Kokkos::parallel_for( @@ -365,4 +373,4 @@ bool pFlow::sphereInteraction::hearChanges fatalErrorInFunction<<"Event "<< msg.eventNames()<< " is not handled in sphereInteraction"< diam_; @@ -60,7 +61,8 @@ struct ppInteractionFunctor deviceViewType1D lVel, deviceViewType1D rVel, deviceViewType1D cForce, - deviceViewType1D cTorque ) + deviceViewType1D cTorque, + forceChain* fChain = nullptr) : dt_(dt), forceModel_(forceModel), @@ -71,7 +73,8 @@ struct ppInteractionFunctor lVel_(lVel), rVel_(rVel), cForce_(cForce), // this is converted to an atomic vector - cTorque_(cTorque) // this is converted to an atomic vector + cTorque_(cTorque),// this is converted to an atomic vector + forceChainPtr_(fChain) {} INLINE_FUNCTION_HD @@ -116,7 +119,12 @@ struct ppInteractionFunctor history, FCn, FCt ); - + + if (forceChainPtr_) + { + const_cast(forceChainPtr_)->addInteraction(i, j, FCn, xi, xj); + } + forceModel_.rollingFriction( dt_, i, j, propId_i, propId_j, @@ -324,4 +332,4 @@ struct pwInteractionFunctor } -#endif //__sphereInteractionKernels_hpp__ \ No newline at end of file +#endif //__sphereInteractionKernels_hpp__ diff --git a/src/Particles/forceChain/forceChain.cpp b/src/Particles/forceChain/forceChain.cpp new file mode 100644 index 000000000..6daa581ca --- /dev/null +++ b/src/Particles/forceChain/forceChain.cpp @@ -0,0 +1,203 @@ +#include "forceChain.hpp" + + +pFlow::forceChain::forceChain( + systemControl& control, + dynamicPointStructure& dynPointStruct +) +: +forceChainActive_(false), +fieldsInitialized_(false), +forceChainFCn_( + objectFile("dummyFCn", "forceChain", objectFile::READ_NEVER, objectFile::WRITE_NEVER), + dynPointStruct, + zero3 + ), + forceChainPos_( + objectFile("dummyDist", "forceChain", objectFile::READ_NEVER, objectFile::WRITE_NEVER), + dynPointStruct, + zero3 + ), + forceChainPairs_( + objectFile("dummyPairs", "forceChain", objectFile::READ_NEVER, objectFile::WRITE_NEVER), + dynPointStruct, + zero3 + ), + pairCounter_( + objectFile("dummyCounter", "forceChain", objectFile::READ_NEVER, objectFile::WRITE_NEVER), + dynPointStruct, + 0 + ) + + + + +{ + +} + + +void pFlow::forceChain::zeroFCn() +{ + if (!forceChainActive_ || !fieldsInitialized_) return; + forceChainFCn_.fill(zero3); +} + +void pFlow::forceChain::zeroDist() +{ + if (!forceChainActive_ || !fieldsInitialized_) return; + forceChainPos_.fill(zero3); +} + +void pFlow::forceChain::zeroPairs() +{ + if (!forceChainActive_ || !fieldsInitialized_) return; + forceChainPairs_.fill(zero3); +} + +void pFlow::forceChain::zeroAll() +{ + if (!forceChainActive_ || !fieldsInitialized_) return; + zeroFCn(); + zeroDist(); + zeroPairs(); + +} + +void pFlow::forceChain::resetPairCounter() +{ + if (!forceChainActive_ || !fieldsInitialized_) return; + Kokkos::deep_copy(pairCounter_.deviceViewAll(), 0u); +} + +bool pFlow::forceChain::initializeFromDict(const dictionary& modelDict) +{ + word forceChainVal = "No"; + forceChainVal = modelDict.getValOrSet("forceChain", forceChainVal); + + if (forceChainVal == "Yes" || forceChainVal == "yes" || + forceChainVal == "True"|| forceChainVal == "true" ) + { + forceChainActive_ = true; + REPORT(1) << "ForceChain is " << Yellow_Text("active") + << " in this simulation." << END_REPORT; + } + else + { + forceChainActive_ = false; + REPORT(1) << "ForceChain is " << Yellow_Text("inactive") + << " in this simulation." << END_REPORT; + } + + return true; +} + + +void pFlow::forceChain::activateForceChain(systemControl& control, dynamicPointStructure& dynPointStruct) +{ + if (!forceChainActive_) + return; + + if (!fieldsInitialized_) + + + + { // Create the fields only when activating + new (&forceChainFCn_) realx3PointField_D( + objectFile( + "forceChainFCn", + "forceChain", + objectFile::READ_IF_PRESENT, + objectFile::WRITE_ALWAYS + ), + dynPointStruct, + zero3 + ); + + new (&forceChainPos_) realx3PointField_D( + objectFile( + "pos", + "forceChain", + objectFile::READ_IF_PRESENT, + objectFile::WRITE_ALWAYS + ), + dynPointStruct, + zero3 + ); + + new (&forceChainPairs_) realx3PointField_D( + objectFile( + "pairs", + "forceChain", + objectFile::READ_IF_PRESENT, + objectFile::WRITE_ALWAYS + ), + dynPointStruct, + zero3 + ); + + new (&pairCounter_) uint32PointField_D( + objectFile( + "pairCount", + "forceChain", + objectFile::READ_NEVER, + objectFile::WRITE_NEVER + ), + dynPointStruct, + 0 + ); + fieldsInitialized_ = true; + REPORT(1) << "ForceChain fields activated and initialized." << END_REPORT; + } +} + +void pFlow::forceChain::addInteraction(uint32 i, uint32 j, const realx3& FCn, const realx3& xi, const realx3& xj) +{ + if (!forceChainActive_ || !fieldsInitialized_) return; + + // Get device views + auto fcnView = forceChainFCn_.deviceViewAll(); + auto distView = forceChainPos_.deviceViewAll(); + auto pairsView = forceChainPairs_.deviceViewAll(); + auto counterView = pairCounter_.deviceViewAll(); + + // 1. Add forces (atomic operations) + Kokkos::atomic_add(&fcnView(i).x_, FCn.x_); + Kokkos::atomic_add(&fcnView(i).y_, FCn.y_); + Kokkos::atomic_add(&fcnView(i).z_, FCn.z_); + + Kokkos::atomic_add(&fcnView(j).x_, -FCn.x_); + Kokkos::atomic_add(&fcnView(j).y_, -FCn.y_); + Kokkos::atomic_add(&fcnView(j).z_, -FCn.z_); + + // 2. Add distances (atomic operations) + Kokkos::atomic_add(&distView(i), xi); + Kokkos::atomic_add(&distView(j), xj); + + // 3. Add pair (only once per unique pair, i < j) + if (i < j) + { + real fcn_mag = length(FCn); + + // Get a unique index using atomic counter + uint32_t pairIndex = Kokkos::atomic_fetch_add(&counterView(0), uint32_t(1)); + + // Bounds check before writing + const uint32_t maxPairs = static_cast(pairsView.extent(0)); + if (pairIndex < maxPairs) + { + // Single write of the tuple (i,j,fcn_mag) + pairsView(pairIndex) = realx3(real(i), real(j), fcn_mag); + } + } +} +bool pFlow::forceChain::isActive() const +{ + return forceChainActive_; +} + + + + + + diff --git a/src/Particles/forceChain/forceChain.hpp b/src/Particles/forceChain/forceChain.hpp new file mode 100644 index 000000000..7fa487efc --- /dev/null +++ b/src/Particles/forceChain/forceChain.hpp @@ -0,0 +1,117 @@ +#ifndef __forceChain_hpp__ +#define __forceChain_hpp__ + +#include "dynamicPointStructure.hpp" +#include "systemControl.hpp" + +namespace pFlow +{ + +class forceChain +{ +private: + + bool forceChainActive_ = false; + + bool fieldsInitialized_ = false; + + /// forceChain force field + realx3PointField_D forceChainFCn_; + + /// Distance Colliding Particles + realx3PointField_D forceChainPos_; + + /// Pairs of Colliding Particles + realx3PointField_D forceChainPairs_; + + /// pair Counter + uint32PointField_D pairCounter_; + +public: + + forceChain( + systemControl& control, + dynamicPointStructure& dynPointStruct + ); + + // Destructor + ~forceChain() = default; + + + // Member functions + void zeroFCn(); + void zeroDist(); + void zeroPairs(); + void zeroAll(); + + // Initialize force chain from dictionary + bool initializeFromDict(const dictionary& modelDict); + + // Reset pair counter (used before sphere-sphere interaction) + void resetPairCounter(); + + // Activate writing force chain + void activateForceChain(systemControl& control,dynamicPointStructure& dynPointStruct); + + void addInteraction(uint32 i, uint32 j, const realx3& FCn, const realx3& xi, const realx3& xj); + + // Getters + inline auto& forceChainFCn() + { + return forceChainFCn_; + } + + inline const auto& forceChainFCn() const + { + return forceChainFCn_; + } + + inline auto& forceChainPos() + { + return forceChainPos_; + } + + inline const auto& forceChainPos() const + { + return forceChainPos_; + } + + inline auto& forceChainPairs() + { + return forceChainPairs_; + } + + inline const auto& forceChainPairs() const + { + return forceChainPairs_; + } + + inline auto& pairCounter() + { + return pairCounter_; + } + + inline const auto& pairCounter() const + { + return pairCounter_; + } + + bool forceChainActive() const + { + return forceChainActive_; + } + + bool fieldsInitialized() const + { + return fieldsInitialized_; + } + + bool isActive() const; + //// + + +}; // forceChain + +} // namespace pFlow + +#endif //__forceChain_hpp__ diff --git a/src/Particles/particles/particles.cpp b/src/Particles/particles/particles.cpp index aefc03bac..d76d19c4d 100644 --- a/src/Particles/particles/particles.cpp +++ b/src/Particles/particles/particles.cpp @@ -54,6 +54,7 @@ pFlow::particles::particles(systemControl& control, const shape& shapes) dynPointStruct_, zero3 ), + forceChain_(control,dynPointStruct_), idHandler_(particleIdHandler::create(dynPointStruct_)), baseFieldBoundaryUpdateTimer_("baseFieldBoundaryUpdate",&timers()) { @@ -62,6 +63,15 @@ pFlow::particles::particles(systemControl& control, const shape& shapes) //idHandler_().initialIdCheck(); } +bool pFlow::particles::initializeForceChain(const dictionary& modelDict) +{ + forceChain_.initializeFromDict(modelDict); + + forceChain_.activateForceChain(control(),dynPointStruct_); + + return true; +} + pFlow::particles::~particles() { // invalidates / unsobscribe from subscriber before its actual destruction diff --git a/src/Particles/particles/particles.hpp b/src/Particles/particles/particles.hpp index 7779d7b66..5a18753ac 100644 --- a/src/Particles/particles/particles.hpp +++ b/src/Particles/particles/particles.hpp @@ -51,6 +51,9 @@ class particles /// contact torque field realx3PointField_D contactTorque_; + /// forceChain + forceChain forceChain_; + /// handling new ids for new particles uniquePtr idHandler_ = nullptr; @@ -100,6 +103,16 @@ class particles ~particles() override; + const forceChain& getForceChain() const + { + return forceChain_; + } + + bool isForceChainActive() const + { + return forceChain_.isActive(); + } + inline const auto& dynPointStruct() const { return dynPointStruct_; @@ -203,6 +216,8 @@ class particles return idHandler_().maxId(); } + bool initializeForceChain(const dictionary& modelDict); + bool beforeIteration() override; bool iterate() override;