From 83307cee5ae63f00f6f2af8d8336bde5bd7a779d Mon Sep 17 00:00:00 2001 From: jmpearl Date: Sun, 11 Feb 2024 12:40:10 -0800 Subject: [PATCH 01/38] adding Spheral::Vector an an any class, we might want to switch to a std::variant and visitor format for the state contents --- src/DataBase/StateBase.cc | 27 +++++++++++++++++++++++++-- 1 file changed, 25 insertions(+), 2 deletions(-) diff --git a/src/DataBase/StateBase.cc b/src/DataBase/StateBase.cc index c11aa57dd..42111c6f0 100644 --- a/src/DataBase/StateBase.cc +++ b/src/DataBase/StateBase.cc @@ -140,7 +140,16 @@ operator==(const StateBase& rhs) const { result = false; } } catch (const boost::bad_any_cast&) { - std::cerr << "StateBase::operator== WARNING: unable to compare values for " << lhsItr->first << "\n"; + try { + auto lhsPtr = boost::any_cast(lhsItr->second); + auto rhsPtr = boost::any_cast(rhsItr->second); + if (*lhsPtr != *rhsPtr) { + cerr << "Vector for " << lhsItr->first << " don't match." << endl; + result = false; + } + } catch (const boost::bad_any_cast&) { + std::cerr << "StateBase::operator== WARNING: unable to compare values for " << lhsItr->first << "\n"; + } } } } @@ -391,12 +400,18 @@ assign(const StateBase& rhs) { const auto rhsptr = boost::any_cast*>(anyrhs); *lhsptr = *rhsptr; } catch(const boost::bad_any_cast&) { + try { + std::cout << "key= " << itr -> first << std::endl; + auto lhsptr = boost::any_cast(anylhs); + const auto rhsptr = boost::any_cast(anyrhs); + *lhsptr = *rhsptr; + } catch(const boost::bad_any_cast&) { // We'll assume other things don't need to be assigned... // VERIFY2(false, "StateBase::assign ERROR: unknown type for key " << itr->first << "\n"); + } } } } - // Copy the connectivity (by reference). This thing is too // big to carry around separate copies! if (rhs.mConnectivityMapPtr != NULL) { @@ -446,8 +461,16 @@ copyState() { itr->second = clone.get(); } catch (const boost::bad_any_cast&) { + try { + auto ptr = boost::any_cast(anythingPtr); + auto clone = std::shared_ptr(new Vector(*ptr)); + mCache.push_back(clone); + itr->second = clone.get(); + + } catch (const boost::bad_any_cast&) { // We'll assume other things don't need to be copied... // VERIFY2(false, "StateBase::copyState ERROR: unrecognized type for " << itr->first << "\n"); + } } } } From 65b29fafb7f9cc17990eccbc88022b5e2ccc1de7 Mon Sep 17 00:00:00 2001 From: jmpearl Date: Sun, 11 Feb 2024 12:45:16 -0800 Subject: [PATCH 02/38] first crack at fixing the moving bc issue --- src/DEM/DEMFieldNames.cc | 3 +- src/DEM/DEMFieldNames.hh | 1 + .../CircularPlaneSolidBoundary.cc | 8 ++++ .../CircularPlaneSolidBoundary.hh | 4 ++ .../SolidBoundary/CylinderSolidBoundary.cc | 9 ++++- .../SolidBoundary/CylinderSolidBoundary.hh | 4 ++ .../InfinitePlaneSolidBoundary.cc | 10 +++++ .../InfinitePlaneSolidBoundary.hh | 6 ++- .../RectangularPlaneSolidBoundary.cc | 7 ++++ .../RectangularPlaneSolidBoundary.hh | 4 ++ src/DEM/SolidBoundary/SolidBoundaryBase.hh | 6 +++ src/DEM/SolidBoundary/SphereSolidBoundary.cc | 8 ++++ src/DEM/SolidBoundary/SphereSolidBoundary.hh | 5 +++ src/PYB11/DEM/SolidBoundaries.py | 40 +++++++++++++++++++ .../DEM/SolidBoundaryBaseAbstractMethods.py | 7 ++++ 15 files changed, 119 insertions(+), 3 deletions(-) diff --git a/src/DEM/DEMFieldNames.cc b/src/DEM/DEMFieldNames.cc index ce4603288..ae4b05361 100644 --- a/src/DEM/DEMFieldNames.cc +++ b/src/DEM/DEMFieldNames.cc @@ -18,4 +18,5 @@ const std::string Spheral::DEMFieldNames::shearDisplacement = "shear displacemen const std::string Spheral::DEMFieldNames::rollingDisplacement = "rolling displacement"; const std::string Spheral::DEMFieldNames::torsionalDisplacement = "torsional displacement"; const std::string Spheral::DEMFieldNames::equilibriumOverlap = "equilibrium overlap"; -const std::string Spheral::DEMFieldNames::solidBoundaries = "solid boundaries"; \ No newline at end of file +const std::string Spheral::DEMFieldNames::solidBoundaries = "solid boundaries"; +const std::string Spheral::DEMFieldNames::solidBoundaryPolicy = "solid boundary policy"; \ No newline at end of file diff --git a/src/DEM/DEMFieldNames.hh b/src/DEM/DEMFieldNames.hh index e3109e33c..4acdafb26 100644 --- a/src/DEM/DEMFieldNames.hh +++ b/src/DEM/DEMFieldNames.hh @@ -24,6 +24,7 @@ struct DEMFieldNames { static const std::string torsionalDisplacement; static const std::string equilibriumOverlap; static const std::string solidBoundaries; + static const std::string solidBoundaryPolicy; }; } diff --git a/src/DEM/SolidBoundary/CircularPlaneSolidBoundary.cc b/src/DEM/SolidBoundary/CircularPlaneSolidBoundary.cc index ec2d559dc..2b565c195 100644 --- a/src/DEM/SolidBoundary/CircularPlaneSolidBoundary.cc +++ b/src/DEM/SolidBoundary/CircularPlaneSolidBoundary.cc @@ -49,6 +49,14 @@ velocity(const Vector& position) const { return mVelocity; } +template +void +CircularPlaneSolidBoundary:: +registerState(DataBase& dataBase, + State& state, + const std::string& boundaryKey) { +} + template void CircularPlaneSolidBoundary:: diff --git a/src/DEM/SolidBoundary/CircularPlaneSolidBoundary.hh b/src/DEM/SolidBoundary/CircularPlaneSolidBoundary.hh index f1e2c0937..703bd52d7 100644 --- a/src/DEM/SolidBoundary/CircularPlaneSolidBoundary.hh +++ b/src/DEM/SolidBoundary/CircularPlaneSolidBoundary.hh @@ -34,6 +34,10 @@ public: virtual Vector distance(const Vector& position) const override; virtual Vector velocity(const Vector& position) const override; + virtual void registerState(DataBase& dataBase, + State& state, + const std::string& boundaryKey) override; + virtual void update(const double multiplier, const double time, const double dt) override; diff --git a/src/DEM/SolidBoundary/CylinderSolidBoundary.cc b/src/DEM/SolidBoundary/CylinderSolidBoundary.cc index 082258876..20ac59838 100644 --- a/src/DEM/SolidBoundary/CylinderSolidBoundary.cc +++ b/src/DEM/SolidBoundary/CylinderSolidBoundary.cc @@ -31,7 +31,6 @@ CylinderSolidBoundary:: ~CylinderSolidBoundary(){ } - template typename Dimension::Vector CylinderSolidBoundary:: @@ -51,6 +50,14 @@ velocity(const Vector& position) const { return mVelocity; } +template +void +CylinderSolidBoundary:: +registerState(DataBase& dataBase, + State& state, + const std::string& boundaryKey) { +} + template void CylinderSolidBoundary:: diff --git a/src/DEM/SolidBoundary/CylinderSolidBoundary.hh b/src/DEM/SolidBoundary/CylinderSolidBoundary.hh index da95e2cfc..7acd1eed1 100644 --- a/src/DEM/SolidBoundary/CylinderSolidBoundary.hh +++ b/src/DEM/SolidBoundary/CylinderSolidBoundary.hh @@ -35,6 +35,10 @@ public: virtual Vector distance(const Vector& position) const override; virtual Vector velocity(const Vector& position) const override; + virtual void registerState(DataBase& dataBase, + State& state, + const std::string& boundaryKey) override; + virtual void update(const double multiplier, const double time, const double dt) override; diff --git a/src/DEM/SolidBoundary/InfinitePlaneSolidBoundary.cc b/src/DEM/SolidBoundary/InfinitePlaneSolidBoundary.cc index a15cb7a44..71e8a2890 100644 --- a/src/DEM/SolidBoundary/InfinitePlaneSolidBoundary.cc +++ b/src/DEM/SolidBoundary/InfinitePlaneSolidBoundary.cc @@ -40,6 +40,16 @@ velocity(const Vector& position) const { return mVelocity; } +template +void +InfinitePlaneSolidBoundary:: +registerState(DataBase& dataBase, + State& state, + const std::string& boundaryKey) { + const auto pointKey = boundaryKey +"_point"; + state.enrollAny(pointKey,mPoint); +} + template void InfinitePlaneSolidBoundary:: diff --git a/src/DEM/SolidBoundary/InfinitePlaneSolidBoundary.hh b/src/DEM/SolidBoundary/InfinitePlaneSolidBoundary.hh index b80539818..782f6b65c 100644 --- a/src/DEM/SolidBoundary/InfinitePlaneSolidBoundary.hh +++ b/src/DEM/SolidBoundary/InfinitePlaneSolidBoundary.hh @@ -31,7 +31,11 @@ public: virtual Vector distance(const Vector& position) const override; virtual Vector velocity(const Vector& position) const override; - + + virtual void registerState(DataBase& dataBase, + State& state, + const std::string& boundaryKey) override; + virtual void update(const double multiplier, const double time, const double dt) override; diff --git a/src/DEM/SolidBoundary/RectangularPlaneSolidBoundary.cc b/src/DEM/SolidBoundary/RectangularPlaneSolidBoundary.cc index 5290b957e..9ff3257a3 100644 --- a/src/DEM/SolidBoundary/RectangularPlaneSolidBoundary.cc +++ b/src/DEM/SolidBoundary/RectangularPlaneSolidBoundary.cc @@ -44,6 +44,13 @@ velocity(const Vector& position) const { return mVelocity; } +template +void +RectangularPlaneSolidBoundary:: +registerState(DataBase& dataBase, + State& state, + const std::string& boundaryKey) { +} template void RectangularPlaneSolidBoundary:: diff --git a/src/DEM/SolidBoundary/RectangularPlaneSolidBoundary.hh b/src/DEM/SolidBoundary/RectangularPlaneSolidBoundary.hh index 3442a6791..60d0d406a 100644 --- a/src/DEM/SolidBoundary/RectangularPlaneSolidBoundary.hh +++ b/src/DEM/SolidBoundary/RectangularPlaneSolidBoundary.hh @@ -35,6 +35,10 @@ public: virtual Vector distance(const Vector& position) const override; virtual Vector velocity(const Vector& position) const override; + virtual void registerState(DataBase& dataBase, + State& state, + const std::string& boundaryKey) override; + virtual void update(const double multiplier, const double time, const double dt) override; diff --git a/src/DEM/SolidBoundary/SolidBoundaryBase.hh b/src/DEM/SolidBoundary/SolidBoundaryBase.hh index 2bd843c4f..0a65e6fe4 100644 --- a/src/DEM/SolidBoundary/SolidBoundaryBase.hh +++ b/src/DEM/SolidBoundary/SolidBoundaryBase.hh @@ -17,6 +17,8 @@ #ifndef __Spheral_SolidBoundaryBase_hh__ #define __Spheral_SolidBoundaryBase_hh__ +#include + namespace Spheral { template class State; @@ -39,6 +41,10 @@ public: virtual Vector distance(const Vector& position) const = 0; virtual Vector velocity(const Vector& position) const = 0; + virtual void registerState(DataBase& dataBase, + State& state, + const std::string& boundaryKey) = 0; + virtual void update(const double multiplier, const double t, const double dt) = 0; diff --git a/src/DEM/SolidBoundary/SphereSolidBoundary.cc b/src/DEM/SolidBoundary/SphereSolidBoundary.cc index b0781cff4..850e97d45 100644 --- a/src/DEM/SolidBoundary/SphereSolidBoundary.cc +++ b/src/DEM/SolidBoundary/SphereSolidBoundary.cc @@ -71,6 +71,14 @@ velocity(const Vector& position) const { return mVelocity; } +template +void +SphereSolidBoundary:: +registerState(DataBase& dataBase, + State& state, + const std::string& boundaryKey) { +} + template void SphereSolidBoundary:: diff --git a/src/DEM/SolidBoundary/SphereSolidBoundary.hh b/src/DEM/SolidBoundary/SphereSolidBoundary.hh index 0aa78f63b..5db055556 100644 --- a/src/DEM/SolidBoundary/SphereSolidBoundary.hh +++ b/src/DEM/SolidBoundary/SphereSolidBoundary.hh @@ -35,6 +35,11 @@ public: virtual Vector distance(const Vector& position) const override; virtual Vector velocity(const Vector& position) const override; + virtual void registerState(DataBase& dataBase, + State& state, + const std::string& boundaryKey) override; + + virtual void update(const double multiplier, const double time, const double dt) override; diff --git a/src/PYB11/DEM/SolidBoundaries.py b/src/PYB11/DEM/SolidBoundaries.py index 8923e5d28..31509ebdd 100644 --- a/src/PYB11/DEM/SolidBoundaries.py +++ b/src/PYB11/DEM/SolidBoundaries.py @@ -33,6 +33,14 @@ def pyinit(point = "const Vector&", normal = "const Vector&"): "solid planar boundary" + @PYB11virtual + def registerState(self, + dataBase = "DataBase<%(Dimension)s>&", + state = "State<%(Dimension)s>&", + boundaryKey = "const std::string&"): + "Register the state solid bc expects to use and evolve." + return "void" + @PYB11virtual def update(self, multiplier = "const double", @@ -77,6 +85,14 @@ def pyinit(point = "const Vector&", basis = "const Tensor&"): "solid planar boundary" + @PYB11virtual + def registerState(self, + dataBase = "DataBase<%(Dimension)s>&", + state = "State<%(Dimension)s>&", + boundaryKey = "const std::string&"): + "Register the state solid bc expects to use and evolve." + return "void" + @PYB11virtual def update(self, multiplier = "const double", @@ -121,6 +137,14 @@ def pyinit(point = "const Vector&", extent = "const Scalar"): "solid planar boundary" + @PYB11virtual + def registerState(self, + dataBase = "DataBase<%(Dimension)s>&", + state = "State<%(Dimension)s>&", + boundaryKey = "const std::string&"): + "Register the state solid bc expects to use and evolve." + return "void" + @PYB11virtual def update(self, multiplier = "const double", @@ -166,6 +190,14 @@ def pyinit(point = "const Vector&", length = "const Scalar"): "solid planar boundary" + @PYB11virtual + def registerState(self, + dataBase = "DataBase<%(Dimension)s>&", + state = "State<%(Dimension)s>&", + boundaryKey = "const std::string&"): + "Register the state solid bc expects to use and evolve." + return "void" + @PYB11virtual def update(self, multiplier = "const double", @@ -213,6 +245,14 @@ def pyinit(center = "const Vector&", clipAxis = "const Vector&"): "solid planar boundary" + @PYB11virtual + def registerState(self, + dataBase = "DataBase<%(Dimension)s>&", + state = "State<%(Dimension)s>&", + boundaryKey = "const std::string&"): + "Register the state solid bc expects to use and evolve." + return "void" + @PYB11virtual def update(self, multiplier = "const double", diff --git a/src/PYB11/DEM/SolidBoundaryBaseAbstractMethods.py b/src/PYB11/DEM/SolidBoundaryBaseAbstractMethods.py index d370930de..777c5d0fc 100644 --- a/src/PYB11/DEM/SolidBoundaryBaseAbstractMethods.py +++ b/src/PYB11/DEM/SolidBoundaryBaseAbstractMethods.py @@ -18,6 +18,13 @@ def distance(self, "distance vector to bc." return "Vector" + def registerState(self, + dataBase = "DataBase<%(Dimension)s>&", + state = "State<%(Dimension)s>&", + boundaryKey = "const std::string&"): + "Register the state solid bc expects to use and evolve." + return "void" + def update(self, multiplier = "const double", t = "const double", From 7eeaa7b4b75562152b3bc456e74472a8a67cd730 Mon Sep 17 00:00:00 2001 From: jmpearl Date: Sun, 11 Feb 2024 12:45:47 -0800 Subject: [PATCH 03/38] letting the solid BCs enroll there own stuff --- src/DEM/DEMBase.cc | 40 +++++++++++++++++++++++++++++++++++++--- src/DEM/DEMBase.hh | 15 ++++++++++++--- 2 files changed, 49 insertions(+), 6 deletions(-) diff --git a/src/DEM/DEMBase.cc b/src/DEM/DEMBase.cc index 6057c68f1..fc684fe5b 100644 --- a/src/DEM/DEMBase.cc +++ b/src/DEM/DEMBase.cc @@ -316,9 +316,21 @@ registerState(DataBase& dataBase, auto rollingDisplacementPolicy = make_policy>>(); auto torsionalDisplacementPolicy = make_policy>>(); + // solid boundary conditions w/ properties that need to be integrated auto boundaryPolicy = make_policy>(mSolidBoundaries); + state.enroll(DEMFieldNames::solidBoundaryPolicy,boundaryPolicy); + const auto& solidBoundaries = this->solidBoundaryConditions(); + const auto numSolidBoundaries = this->numSolidBoundaries(); + for (auto ibc = 0u; ibc < numSolidBoundaries; ++ibc){ + const auto name = "SolidBoundary_" + std::to_string(this->getSolidBoundaryUniqueIndex(ibc)); + solidBoundaries[ibc]->registerState(dataBase,state,name); + } - state.enroll(DEMFieldNames::solidBoundaries,boundaryPolicy); + // for (ConstSolidBoundaryIterator solidboundItr = mSolidBoundaries.begin(); + // solidboundItr != mSolidBoundaries.end(); + // ++solidboundItr){ + // (*solidboundItr)->registerState(dataBase,state); + // } state.enroll(mTimeStepMask); state.enroll(mass); @@ -410,6 +422,21 @@ preStepInitialize(const DataBase& dataBase, TIME_END("DEMpreStepInitialize"); } +//------------------------------------------------------------------------------ +// This method is called once at the beginning of a timestep, after all state registration. +//------------------------------------------------------------------------------ +// template +// void +// DEMBase:: +// finalize(const DataBase& dataBase, +// State& state, +// StateDerivatives& derivatives) { +// TIME_BEGIN("DEMfinalize"); + + +// TIME_END("DEMfinalize"); +// } + //------------------------------------------------------------------------------ // Call before deriv evaluation //------------------------------------------------------------------------------ @@ -650,7 +677,14 @@ finalizeAfterRedistribution() { } //------------------------------------------------------------------------------ -// redistribution sub function +// redistribution sub function -- we store the pairwise fields in one of the +// nodes. Prior to redistribution, we don't know where the new domain +// boundaries will be so we want to store the pairwise data in both nodes. +// pairs that span a domain boundary already do this so we just need to make +// sure all the internal contacts data is stored in both nodes. We do this +// by looking through the contacts for each node. If the pair node is an +// internal node than it needs to add the current (storage node) as a contact. +// The NeighborIndices needs special treatement so it gets its own method //------------------------------------------------------------------------------ template template @@ -861,7 +895,7 @@ updateContactMap(const DataBase& dataBase){ // now add our contact #pragma omp critical { - mContactStorageIndices.push_back(ContactIndex(nodeListi, // storage nodelist index + mContactStorageIndices.push_back(ContactIndex(nodeListi, // storage nodelist index i, // storage node index storageContactIndex, // storage contact index ibc)); // bc index diff --git a/src/DEM/DEMBase.hh b/src/DEM/DEMBase.hh index 5b158f0d3..1604943d9 100644 --- a/src/DEM/DEMBase.hh +++ b/src/DEM/DEMBase.hh @@ -87,6 +87,14 @@ public: const DataBase& dataBase, const State& state, StateDerivatives& derivs) const override; + + // hook after the intergrator step + // virtual + // void finalize(const Scalar time, + // const Scalar dt, + // const DataBase& dataBase, + // const State& state, + // StateDerivatives& derivs) const override; // Apply boundary conditions to the physics specific fields. virtual @@ -218,11 +226,11 @@ protected: std::vector*> mSolidBoundaries; - int mCycle; - int mContactRemovalFrequency; + int mCycle; // current cycle + int mContactRemovalFrequency; // how often do we clear out old contacts // number of steps per collision time-scale - Scalar mStepsPerCollision; + Scalar mStepsPerCollision; // Optional bounding box for generating the mesh. Vector mxmin, mxmax; @@ -250,6 +258,7 @@ protected: FieldList> mDDtTorsionalDisplacement; // derivative to evolve frictional spring displacement FieldList> mNewTorsionalDisplacement; // handles rotation of frictional spring and reset on slip + // map to storage location from connectivityMap to pairwise fieldlists std::vector mContactStorageIndices; // The restart registration. From f92099acd5996c6068ea37e9445affa80095a071 Mon Sep 17 00:00:00 2001 From: jmpearl Date: Sun, 11 Feb 2024 13:55:23 -0800 Subject: [PATCH 04/38] updating DEM timers --- src/DEM/DEMBase.cc | 18 ++++++++++++++++-- src/DEM/LinearSpringDEM.cc | 13 ++++++++++++- 2 files changed, 28 insertions(+), 3 deletions(-) diff --git a/src/DEM/DEMBase.cc b/src/DEM/DEMBase.cc index fc684fe5b..9d67219ac 100644 --- a/src/DEM/DEMBase.cc +++ b/src/DEM/DEMBase.cc @@ -448,8 +448,9 @@ initialize(const Scalar time, const DataBase& dataBase, State& state, StateDerivatives& derivs){ +TIME_BEGIN("DEMinitialize"); - +TIME_END("DEMinitialize"); } @@ -603,6 +604,8 @@ void DEMBase:: initializeOverlap(const DataBase& dataBase, const int startingCompositeParticleIndex){ + TIME_BEGIN("DEMinitializeOverlap"); + const auto& connectivityMap = dataBase.connectivityMap(); const auto& pairs = connectivityMap.nodePairList(); const auto numPairs = pairs.size(); @@ -643,6 +646,7 @@ initializeOverlap(const DataBase& dataBase, const int startingComposi mEquilibriumOverlap(storeNodeList,storeNode)[storeContact] = delta0; } } + TIME_END("DEMinitializeOverlap"); } //------------------------------------------------------------------------------ // Redistribution methods -- before we redistribute, we are going to make sure @@ -655,6 +659,9 @@ template void DEMBase:: initializeBeforeRedistribution(){ + + TIME_BEGIN("DEMinitializeBeforeRedistribution"); + this->prepNeighborIndicesForRedistribution(); this->prepPairFieldListForRedistribution(mShearDisplacement); @@ -668,6 +675,9 @@ initializeBeforeRedistribution(){ this->prepPairFieldListForRedistribution(mNewRollingDisplacement); this->prepPairFieldListForRedistribution(mDDtTorsionalDisplacement); this->prepPairFieldListForRedistribution(mNewTorsionalDisplacement); + + TIME_END("DEMinitializeBeforeRedistribution"); + } template @@ -786,6 +796,8 @@ void DEMBase:: updateContactMap(const DataBase& dataBase){ + TIME_BEGIN("DEMupdateContactMap"); + const auto& uniqueIndex = dataBase.DEMUniqueIndex(); const auto& radius = dataBase.DEMParticleRadius(); const auto& position = dataBase.DEMPosition(); @@ -904,6 +916,7 @@ updateContactMap(const DataBase& dataBase){ } // loop nodes } // loop nodelists } // loop solid boundaries + TIME_END("DEMupdateContactMap"); } // method //------------------------------------------------------------------------------ @@ -913,7 +926,7 @@ template void DEMBase:: identifyInactiveContacts(const DataBase& dataBase){ - + TIME_BEGIN("DEMidentifyInactiveContacts"); const auto bufferDistance = dataBase.maxNeighborSearchBuffer(); const auto numNodeLists = dataBase.numNodeLists(); const auto& nodeListPtrs = dataBase.DEMNodeListPtrs(); @@ -968,5 +981,6 @@ identifyInactiveContacts(const DataBase& dataBase){ } // loop particle bc contacts } // omp parallel region +TIME_END("DEMidentifyInactiveContacts"); } // class } // namespace diff --git a/src/DEM/LinearSpringDEM.cc b/src/DEM/LinearSpringDEM.cc index 05cef1c15..bed1bc730 100644 --- a/src/DEM/LinearSpringDEM.cc +++ b/src/DEM/LinearSpringDEM.cc @@ -33,6 +33,8 @@ #include "DEM/ContactStorageLocation.hh" #include "DEM/SolidBoundary/SolidBoundaryBase.hh" +#include "Utilities/Timer.hh" + #ifdef _OPENMP #include "omp.h" #endif @@ -115,6 +117,8 @@ dt(const DataBase& dataBase, const StateDerivatives& /*derivs*/, const typename Dimension::Scalar /*currentTime*/) const{ + TIME_BEGIN("LinearSpringDEMdt"); + // Get some useful fluid variables from the DataBase. const auto mass = state.fields(HydroFieldNames::mass, 0.0); const auto position = state.fields(HydroFieldNames::position, Vector::zero); @@ -268,6 +272,7 @@ dt(const DataBase& dataBase, } } + TIME_END("LinearSpringDEMdt"); return result; } @@ -279,8 +284,10 @@ template void LinearSpringDEM:: initializeProblemStartup(DataBase& dataBase){ + TIME_BEGIN("LinearSpringDEMinitializeProblemStartup"); DEMBase::initializeProblemStartup(dataBase); this->setMomentOfInertia(); + TIME_END("LinearSpringDEMinitializeProblemStartup"); } //------------------------------------------------------------------------------ @@ -291,9 +298,11 @@ void LinearSpringDEM:: registerState(DataBase& dataBase, State& state) { + TIME_BEGIN("LinearSpringDEMregisterState"); DEMBase::registerState(dataBase,state); dataBase.resizeDEMFieldList(mMomentOfInertia, 0.0, DEMFieldNames::momentOfInertia, false); state.enroll(mMomentOfInertia); + TIME_END("LinearSpringDEMregisterState"); } //------------------------------------------------------------------------------ // evaluate the derivatives @@ -306,7 +315,7 @@ evaluateDerivatives(const typename Dimension::Scalar /*time*/, const DataBase& dataBase, const State& state, StateDerivatives& derivatives) const{ - + TIME_BEGIN("LinearSpringDEMevaluateDerivatives"); this->resizeDerivativePairFieldLists(derivatives); // A few useful constants we'll use in the following loop. @@ -796,6 +805,8 @@ evaluateDerivatives(const typename Dimension::Scalar /*time*/, DxDt(nodeListi,i) = veli; } // loop nodes } // loop nodelists + + TIME_END("LinearSpringDEMevaluateDerivatives"); } // method From 79b9bb2b1bba26657bab221f6801783a2dee93db Mon Sep 17 00:00:00 2001 From: jmpearl Date: Sun, 11 Feb 2024 13:55:51 -0800 Subject: [PATCH 05/38] adding more solid bc fields to state --- src/DEM/SolidBoundary/InfinitePlaneSolidBoundary.cc | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/src/DEM/SolidBoundary/InfinitePlaneSolidBoundary.cc b/src/DEM/SolidBoundary/InfinitePlaneSolidBoundary.cc index 71e8a2890..00edba84f 100644 --- a/src/DEM/SolidBoundary/InfinitePlaneSolidBoundary.cc +++ b/src/DEM/SolidBoundary/InfinitePlaneSolidBoundary.cc @@ -47,7 +47,11 @@ registerState(DataBase& dataBase, State& state, const std::string& boundaryKey) { const auto pointKey = boundaryKey +"_point"; + const auto velocityKey = boundaryKey +"_velocity"; + const auto normalKey = boundaryKey +"_normal"; state.enrollAny(pointKey,mPoint); + state.enrollAny(velocityKey,mVelocity); + state.enrollAny(normalKey,mNormal); } template From 4db2b6740cd97f6e0654f160b73efb1182e5af9c Mon Sep 17 00:00:00 2001 From: jmpearl Date: Fri, 1 Mar 2024 09:33:54 -0800 Subject: [PATCH 06/38] updating the DEM SB update approach --- src/DEM/SolidBoundary/CircularPlaneSolidBoundary.cc | 6 ++++++ src/DEM/SolidBoundary/RectangularPlaneSolidBoundary.cc | 4 ++++ src/DEM/SolidBoundary/SolidBoundaryBase.hh | 9 ++++++++- 3 files changed, 18 insertions(+), 1 deletion(-) diff --git a/src/DEM/SolidBoundary/CircularPlaneSolidBoundary.cc b/src/DEM/SolidBoundary/CircularPlaneSolidBoundary.cc index 2b565c195..fdd6f0aae 100644 --- a/src/DEM/SolidBoundary/CircularPlaneSolidBoundary.cc +++ b/src/DEM/SolidBoundary/CircularPlaneSolidBoundary.cc @@ -55,6 +55,12 @@ CircularPlaneSolidBoundary:: registerState(DataBase& dataBase, State& state, const std::string& boundaryKey) { + const auto pointKey = boundaryKey +"_point"; + const auto velocityKey = boundaryKey +"_velocity"; + const auto normalKey = boundaryKey +"_normal"; + state.enrollAny(pointKey,mPoint); + state.enrollAny(pointKey,mVelocity); + state.enrollAny(pointKey,mNormal); } template diff --git a/src/DEM/SolidBoundary/RectangularPlaneSolidBoundary.cc b/src/DEM/SolidBoundary/RectangularPlaneSolidBoundary.cc index 9ff3257a3..d438348fa 100644 --- a/src/DEM/SolidBoundary/RectangularPlaneSolidBoundary.cc +++ b/src/DEM/SolidBoundary/RectangularPlaneSolidBoundary.cc @@ -50,6 +50,10 @@ RectangularPlaneSolidBoundary:: registerState(DataBase& dataBase, State& state, const std::string& boundaryKey) { + const auto pointKey = boundaryKey +"_point"; + const auto velocityKey = boundaryKey +"_velocity"; + state.enrollAny(pointKey,mPoint); + state.enrollAny(velocityKey,mVelocity); } template void diff --git a/src/DEM/SolidBoundary/SolidBoundaryBase.hh b/src/DEM/SolidBoundary/SolidBoundaryBase.hh index 0a65e6fe4..43ab437e5 100644 --- a/src/DEM/SolidBoundary/SolidBoundaryBase.hh +++ b/src/DEM/SolidBoundary/SolidBoundaryBase.hh @@ -10,7 +10,14 @@ // is implemented in the DEMBase class. The boundaries // themselves simply define the surface and how it evolves // in time. -// +//----------------------------------------------------------------------------// +// ToDo +// -- add complete registration methods +// -- set the unique index of the bc so no std::str input req. for regState +// -- make restartable +// -- stateBase add scalar and tensor (talk to Mike about different pattern) +// -- tests +//----------------------------------------------------------------------------- // J.M. Pearl 2023 //----------------------------------------------------------------------------// From 74318c8846b9c1f772902dd00dca1bb047aea480 Mon Sep 17 00:00:00 2001 From: jmpearl Date: Sat, 2 Mar 2024 10:26:45 -0800 Subject: [PATCH 07/38] adding unique id as a solid bc property --- src/DEM/DEMBase.cc | 3 +-- src/DEM/SolidBoundary/SolidBoundaryBase.hh | 11 +++++++++- .../SolidBoundary/SolidBoundaryBaseInline.hh | 20 +++++++++++++++++++ 3 files changed, 31 insertions(+), 3 deletions(-) create mode 100644 src/DEM/SolidBoundary/SolidBoundaryBaseInline.hh diff --git a/src/DEM/DEMBase.cc b/src/DEM/DEMBase.cc index 9d67219ac..77f165bd5 100644 --- a/src/DEM/DEMBase.cc +++ b/src/DEM/DEMBase.cc @@ -18,7 +18,6 @@ #include "DataBase/ReplaceState.hh" #include "DataBase/DataBase.hh" - #include "Field/FieldList.hh" #include "Field/NodeIterators.hh" @@ -47,7 +46,6 @@ #include "omp.h" #endif - #include #include #include @@ -319,6 +317,7 @@ registerState(DataBase& dataBase, // solid boundary conditions w/ properties that need to be integrated auto boundaryPolicy = make_policy>(mSolidBoundaries); state.enroll(DEMFieldNames::solidBoundaryPolicy,boundaryPolicy); + const auto& solidBoundaries = this->solidBoundaryConditions(); const auto numSolidBoundaries = this->numSolidBoundaries(); for (auto ibc = 0u; ibc < numSolidBoundaries; ++ibc){ diff --git a/src/DEM/SolidBoundary/SolidBoundaryBase.hh b/src/DEM/SolidBoundary/SolidBoundaryBase.hh index 43ab437e5..f02737e62 100644 --- a/src/DEM/SolidBoundary/SolidBoundaryBase.hh +++ b/src/DEM/SolidBoundary/SolidBoundaryBase.hh @@ -55,8 +55,17 @@ public: virtual void update(const double multiplier, const double t, const double dt) = 0; -}; + + void uniqueIndex(const int uId); + int uniqueIndex() const; + +private: +//--------------------------- Public Interface ---------------------------// +int mUniqueIndex; +}; } +#include "SolidBoundaryBaseInline.hh" + #endif diff --git a/src/DEM/SolidBoundary/SolidBoundaryBaseInline.hh b/src/DEM/SolidBoundary/SolidBoundaryBaseInline.hh new file mode 100644 index 000000000..edd8dc21e --- /dev/null +++ b/src/DEM/SolidBoundary/SolidBoundaryBaseInline.hh @@ -0,0 +1,20 @@ +namespace Spheral { +//------------------------------------------------------------------------------ +// set/get unique index for the solid bc +//------------------------------------------------------------------------------ +template +inline +void +SolidBoundaryBase:: +uniqueIndex(const int uId){ + mUniqueIndex = uId; +} +template +inline +int +SolidBoundaryBase:: +uniqueIndex() const { + return mUniqueIndex; +}; +} + From 2283c71642c578b5c42c6fdc86a28b51ceba714d Mon Sep 17 00:00:00 2001 From: jmpearl Date: Sat, 2 Mar 2024 12:06:29 -0800 Subject: [PATCH 08/38] Improving DEM solidBC unique indices to accomdate mid simulation removal. This improvement addresses a bug found by Sean Wiggins where particles were erronously getting stuck on solid boundaries if one was deleted. The previous unique index scheme caused pairwise integrate properties to inadvertently get transfered to another solid BC if one was deleted. This commit also gives the user the option to delete a solid bc by passing that bc to the DEM package --- src/DEM/DEMBase.cc | 5 +++-- src/DEM/DEMBase.hh | 4 +++- src/DEM/DEMBaseInline.hh | 40 +++++++++++++++++++++++++++++----------- src/PYB11/DEM/DEMBase.py | 9 ++++----- 4 files changed, 39 insertions(+), 19 deletions(-) diff --git a/src/DEM/DEMBase.cc b/src/DEM/DEMBase.cc index 77f165bd5..c09250635 100644 --- a/src/DEM/DEMBase.cc +++ b/src/DEM/DEMBase.cc @@ -79,6 +79,7 @@ DEMBase(const DataBase& dataBase, const Vector& xmax): Physics(), mDataBase(dataBase), + mNewSolidBoundaryIndex(0), mSolidBoundaries(), mCycle(0), mContactRemovalFrequency((int)stepsPerCollision), @@ -321,7 +322,7 @@ registerState(DataBase& dataBase, const auto& solidBoundaries = this->solidBoundaryConditions(); const auto numSolidBoundaries = this->numSolidBoundaries(); for (auto ibc = 0u; ibc < numSolidBoundaries; ++ibc){ - const auto name = "SolidBoundary_" + std::to_string(this->getSolidBoundaryUniqueIndex(ibc)); + const auto name = "SolidBoundary_" + std::to_string(solidBoundaries[ibc]->uniqueIndex()); solidBoundaries[ibc]->registerState(dataBase,state,name); } @@ -891,7 +892,7 @@ updateContactMap(const DataBase& dataBase){ if (disBc.magnitude() < Ri*(1+bufferDistance)){ // create a unique index for the boundary condition - const auto uId_bc = this->getSolidBoundaryUniqueIndex(ibc); + const auto uId_bc = solidBoundaryi->uniqueIndex(); // check to see if it already exists const auto neighborContacts = mNeighborIndices(nodeListi,i); diff --git a/src/DEM/DEMBase.hh b/src/DEM/DEMBase.hh index 1604943d9..2053eb918 100644 --- a/src/DEM/DEMBase.hh +++ b/src/DEM/DEMBase.hh @@ -200,12 +200,13 @@ public: const Vector vrotj) const; // Solid Bounderies + int newSolidBoundaryIndex() const; void appendSolidBoundary(SolidBoundaryBase& boundary); void clearSolidBoundaries(); + void removeSolidBoundary(const SolidBoundaryBase& boundary); bool haveSolidBoundary(const SolidBoundaryBase& boundary) const; unsigned int numSolidBoundaries() const; const std::vector*>& solidBoundaryConditions() const; - int getSolidBoundaryUniqueIndex(const int x) const; // counts unsigned int numParticleParticleContacts() const; @@ -224,6 +225,7 @@ protected: const DataBase& mDataBase; + int mNewSolidBoundaryIndex; std::vector*> mSolidBoundaries; int mCycle; // current cycle diff --git a/src/DEM/DEMBaseInline.hh b/src/DEM/DEMBaseInline.hh index dffe855e2..0a7f3b94d 100644 --- a/src/DEM/DEMBaseInline.hh +++ b/src/DEM/DEMBaseInline.hh @@ -321,8 +321,6 @@ rollingMoment(const Dim<3>::Vector rhatij, return DEMDimension>::cross((vroti + vrotj),rhatij).unitVector(); } - - //------------------------------------------------------------------------------ // Add a Boundary condition to the end of the current boundary list. //------------------------------------------------------------------------------ @@ -331,7 +329,9 @@ inline void DEMBase:: appendSolidBoundary(SolidBoundaryBase& boundary) { - mSolidBoundaries.push_back(&boundary); + mNewSolidBoundaryIndex -= 1; + boundary.uniqueIndex(mNewSolidBoundaryIndex); + mSolidBoundaries.push_back(&boundary); } //------------------------------------------------------------------------------ @@ -345,6 +345,15 @@ clearSolidBoundaries() { mSolidBoundaries = std::vector*>(); } + +template +inline +int +DEMBase:: +newSolidBoundaryIndex() const { + return mNewSolidBoundaryIndex; +} + //------------------------------------------------------------------------------ // Test if the given Boundary condition is listed in the physics package. //------------------------------------------------------------------------------ @@ -356,6 +365,23 @@ haveSolidBoundary(const SolidBoundaryBase& boundary) const { return std::count(mSolidBoundaries.begin(), mSolidBoundaries.end(), &boundary) > 0; } +template +inline +void +DEMBase:: +removeSolidBoundary(const SolidBoundaryBase& boundary) { + + const auto bcPtr = std::find(mSolidBoundaries.begin(),mSolidBoundaries.end(),&boundary); + if (bcPtr != mSolidBoundaries.end()){ + //const int bcIndex = std::distance(mSolidBoundaries.begin(),bcPtr); + mSolidBoundaries.erase(bcPtr); + //std::cout << "removing solid bc " << bcIndex << std::endl; + } else { + std:: cout << "THEY AINT HOME"; + } +} + + template inline unsigned int @@ -401,12 +427,4 @@ numParticleBoundaryContacts() const { return (mContactStorageIndices.size()-this->numParticleParticleContacts()); } -template -inline -int -DEMBase:: -getSolidBoundaryUniqueIndex(const int x) const { - return -x-1; -} - } diff --git a/src/PYB11/DEM/DEMBase.py b/src/PYB11/DEM/DEMBase.py index 0627717eb..0969720fa 100644 --- a/src/PYB11/DEM/DEMBase.py +++ b/src/PYB11/DEM/DEMBase.py @@ -113,11 +113,10 @@ def numSolidBoundaries(self): def haveSolidBoundary(boundary = "const SolidBoundaryBase<%(Dimension)s>&"): "is this boundary being tracked?" return "bool" - - @PYB11const - def getSolidBoundaryUniqueIndex(x="const int"): - "Unique index for neighborIndices pairFieldList (returns -x-1)" - return "int" + + def removeSolidBoundary(boundary = "const SolidBoundaryBase<%(Dimension)s>&"): + "remove the specified solid boundary" + return "void" #........................................................................... # Properties From e0bce89889f4d0b58589f5582a5599af65996841 Mon Sep 17 00:00:00 2001 From: jmpearl Date: Sat, 2 Mar 2024 12:14:11 -0800 Subject: [PATCH 09/38] adding to / modifying python exposure --- src/PYB11/DEM/DEMBase.py | 10 ++++++---- 1 file changed, 6 insertions(+), 4 deletions(-) diff --git a/src/PYB11/DEM/DEMBase.py b/src/PYB11/DEM/DEMBase.py index 0969720fa..3f2620e59 100644 --- a/src/PYB11/DEM/DEMBase.py +++ b/src/PYB11/DEM/DEMBase.py @@ -104,10 +104,10 @@ def clearSolidBoundaries(self): "remove all solid boundaries from the dem package" return "void" - @PYB11const - def numSolidBoundaries(self): - "return the number of solid boundaries being tracked" - return "unsigned int" + # @PYB11const + # def numSolidBoundaries(self): + # "return the number of solid boundaries being tracked" + # return "unsigned int" @PYB11const def haveSolidBoundary(boundary = "const SolidBoundaryBase<%(Dimension)s>&"): @@ -144,6 +144,8 @@ def removeSolidBoundary(boundary = "const SolidBoundaryBase<%(Dimension)s>&"): DDtShearDisplacement = PYB11property("const FieldList<%(Dimension)s, vector>&","DDtShearDisplacement", returnpolicy="reference_internal") isActiveContact = PYB11property("const FieldList<%(Dimension)s, vector>&","isActiveContact", returnpolicy="reference_internal") + newSolidBoundaryIndex = PYB11property("int", "newSolidBoundaryIndex", doc="index of the most recent solid bc added to the package") + numSolidBoundaries = PYB11property("unsigned int", "numSolidBoundaries", doc="number of solid boundaries") numContacts = PYB11property("unsigned int", "numContacts", doc="Total number of contacts") numParticleParticleContacts = PYB11property("unsigned int", "numParticleParticleContacts", doc="Number of interactions with other dem particles") numParticleBoundaryContacts = PYB11property("unsigned int", "numParticleBoundaryContacts", doc="Number interactions with solid boundaries") From 79150e9269868144872253723006bcb4e213b41a Mon Sep 17 00:00:00 2001 From: jmpearl Date: Sat, 2 Mar 2024 12:37:39 -0800 Subject: [PATCH 10/38] moving the solid BC name into the solid bc --- src/DEM/DEMBase.cc | 10 +++------- .../SolidBoundary/CircularPlaneSolidBoundary.cc | 4 ++-- .../SolidBoundary/CircularPlaneSolidBoundary.hh | 3 +-- src/DEM/SolidBoundary/CylinderSolidBoundary.cc | 3 +-- src/DEM/SolidBoundary/CylinderSolidBoundary.hh | 3 +-- .../SolidBoundary/InfinitePlaneSolidBoundary.cc | 4 ++-- .../SolidBoundary/InfinitePlaneSolidBoundary.hh | 3 +-- .../RectangularPlaneSolidBoundary.cc | 4 ++-- .../RectangularPlaneSolidBoundary.hh | 3 +-- src/DEM/SolidBoundary/SolidBoundaryBase.hh | 3 +-- src/DEM/SolidBoundary/SphereSolidBoundary.cc | 3 +-- src/DEM/SolidBoundary/SphereSolidBoundary.hh | 3 +-- src/PYB11/DEM/SolidBoundaries.py | 15 +++++---------- src/PYB11/DEM/SolidBoundaryBaseAbstractMethods.py | 3 +-- 14 files changed, 23 insertions(+), 41 deletions(-) diff --git a/src/DEM/DEMBase.cc b/src/DEM/DEMBase.cc index c09250635..604ebc8d8 100644 --- a/src/DEM/DEMBase.cc +++ b/src/DEM/DEMBase.cc @@ -319,13 +319,9 @@ registerState(DataBase& dataBase, auto boundaryPolicy = make_policy>(mSolidBoundaries); state.enroll(DEMFieldNames::solidBoundaryPolicy,boundaryPolicy); - const auto& solidBoundaries = this->solidBoundaryConditions(); - const auto numSolidBoundaries = this->numSolidBoundaries(); - for (auto ibc = 0u; ibc < numSolidBoundaries; ++ibc){ - const auto name = "SolidBoundary_" + std::to_string(solidBoundaries[ibc]->uniqueIndex()); - solidBoundaries[ibc]->registerState(dataBase,state,name); - } - + for (auto ibc = 0u; ibc < this->numSolidBoundaries(); ++ibc){ + mSolidBoundaries[ibc]->registerState(dataBase,state);} + //const auto name = "SolidBoundary_" + std::to_string(solidBoundaries[ibc]->uniqueIndex()); // for (ConstSolidBoundaryIterator solidboundItr = mSolidBoundaries.begin(); // solidboundItr != mSolidBoundaries.end(); // ++solidboundItr){ diff --git a/src/DEM/SolidBoundary/CircularPlaneSolidBoundary.cc b/src/DEM/SolidBoundary/CircularPlaneSolidBoundary.cc index fdd6f0aae..de6b6002a 100644 --- a/src/DEM/SolidBoundary/CircularPlaneSolidBoundary.cc +++ b/src/DEM/SolidBoundary/CircularPlaneSolidBoundary.cc @@ -53,8 +53,8 @@ template void CircularPlaneSolidBoundary:: registerState(DataBase& dataBase, - State& state, - const std::string& boundaryKey) { + State& state) { + const auto boundaryKey = "CircularPlaneSolidBoundary_" + std::to_string(std::abs(this->uniqueIndex())); const auto pointKey = boundaryKey +"_point"; const auto velocityKey = boundaryKey +"_velocity"; const auto normalKey = boundaryKey +"_normal"; diff --git a/src/DEM/SolidBoundary/CircularPlaneSolidBoundary.hh b/src/DEM/SolidBoundary/CircularPlaneSolidBoundary.hh index 703bd52d7..aaa75f9d1 100644 --- a/src/DEM/SolidBoundary/CircularPlaneSolidBoundary.hh +++ b/src/DEM/SolidBoundary/CircularPlaneSolidBoundary.hh @@ -35,8 +35,7 @@ public: virtual Vector velocity(const Vector& position) const override; virtual void registerState(DataBase& dataBase, - State& state, - const std::string& boundaryKey) override; + State& state) override; virtual void update(const double multiplier, const double time, diff --git a/src/DEM/SolidBoundary/CylinderSolidBoundary.cc b/src/DEM/SolidBoundary/CylinderSolidBoundary.cc index 20ac59838..80c5bd30d 100644 --- a/src/DEM/SolidBoundary/CylinderSolidBoundary.cc +++ b/src/DEM/SolidBoundary/CylinderSolidBoundary.cc @@ -54,8 +54,7 @@ template void CylinderSolidBoundary:: registerState(DataBase& dataBase, - State& state, - const std::string& boundaryKey) { + State& state) { } template diff --git a/src/DEM/SolidBoundary/CylinderSolidBoundary.hh b/src/DEM/SolidBoundary/CylinderSolidBoundary.hh index 7acd1eed1..e4cbe30b0 100644 --- a/src/DEM/SolidBoundary/CylinderSolidBoundary.hh +++ b/src/DEM/SolidBoundary/CylinderSolidBoundary.hh @@ -36,8 +36,7 @@ public: virtual Vector velocity(const Vector& position) const override; virtual void registerState(DataBase& dataBase, - State& state, - const std::string& boundaryKey) override; + State& state) override; virtual void update(const double multiplier, const double time, diff --git a/src/DEM/SolidBoundary/InfinitePlaneSolidBoundary.cc b/src/DEM/SolidBoundary/InfinitePlaneSolidBoundary.cc index 00edba84f..a755eca0c 100644 --- a/src/DEM/SolidBoundary/InfinitePlaneSolidBoundary.cc +++ b/src/DEM/SolidBoundary/InfinitePlaneSolidBoundary.cc @@ -44,8 +44,8 @@ template void InfinitePlaneSolidBoundary:: registerState(DataBase& dataBase, - State& state, - const std::string& boundaryKey) { + State& state) { + const auto boundaryKey = "InfinitePlaneSolidBoundary_" + std::to_string(std::abs(this->uniqueIndex())); const auto pointKey = boundaryKey +"_point"; const auto velocityKey = boundaryKey +"_velocity"; const auto normalKey = boundaryKey +"_normal"; diff --git a/src/DEM/SolidBoundary/InfinitePlaneSolidBoundary.hh b/src/DEM/SolidBoundary/InfinitePlaneSolidBoundary.hh index 782f6b65c..eff5365f0 100644 --- a/src/DEM/SolidBoundary/InfinitePlaneSolidBoundary.hh +++ b/src/DEM/SolidBoundary/InfinitePlaneSolidBoundary.hh @@ -33,8 +33,7 @@ public: virtual Vector velocity(const Vector& position) const override; virtual void registerState(DataBase& dataBase, - State& state, - const std::string& boundaryKey) override; + State& state) override; virtual void update(const double multiplier, const double time, diff --git a/src/DEM/SolidBoundary/RectangularPlaneSolidBoundary.cc b/src/DEM/SolidBoundary/RectangularPlaneSolidBoundary.cc index d438348fa..57d0abc2a 100644 --- a/src/DEM/SolidBoundary/RectangularPlaneSolidBoundary.cc +++ b/src/DEM/SolidBoundary/RectangularPlaneSolidBoundary.cc @@ -48,8 +48,8 @@ template void RectangularPlaneSolidBoundary:: registerState(DataBase& dataBase, - State& state, - const std::string& boundaryKey) { + State& state) { + const auto boundaryKey = "RectangularPlaneSolidBoundary_" + std::to_string(std::abs(this->uniqueIndex())); const auto pointKey = boundaryKey +"_point"; const auto velocityKey = boundaryKey +"_velocity"; state.enrollAny(pointKey,mPoint); diff --git a/src/DEM/SolidBoundary/RectangularPlaneSolidBoundary.hh b/src/DEM/SolidBoundary/RectangularPlaneSolidBoundary.hh index 60d0d406a..bef29204a 100644 --- a/src/DEM/SolidBoundary/RectangularPlaneSolidBoundary.hh +++ b/src/DEM/SolidBoundary/RectangularPlaneSolidBoundary.hh @@ -36,8 +36,7 @@ public: virtual Vector velocity(const Vector& position) const override; virtual void registerState(DataBase& dataBase, - State& state, - const std::string& boundaryKey) override; + State& state) override; virtual void update(const double multiplier, const double time, diff --git a/src/DEM/SolidBoundary/SolidBoundaryBase.hh b/src/DEM/SolidBoundary/SolidBoundaryBase.hh index f02737e62..6960fd5cd 100644 --- a/src/DEM/SolidBoundary/SolidBoundaryBase.hh +++ b/src/DEM/SolidBoundary/SolidBoundaryBase.hh @@ -49,8 +49,7 @@ public: virtual Vector velocity(const Vector& position) const = 0; virtual void registerState(DataBase& dataBase, - State& state, - const std::string& boundaryKey) = 0; + State& state) = 0; virtual void update(const double multiplier, const double t, diff --git a/src/DEM/SolidBoundary/SphereSolidBoundary.cc b/src/DEM/SolidBoundary/SphereSolidBoundary.cc index 850e97d45..df79479eb 100644 --- a/src/DEM/SolidBoundary/SphereSolidBoundary.cc +++ b/src/DEM/SolidBoundary/SphereSolidBoundary.cc @@ -75,8 +75,7 @@ template void SphereSolidBoundary:: registerState(DataBase& dataBase, - State& state, - const std::string& boundaryKey) { + State& state) { } template diff --git a/src/DEM/SolidBoundary/SphereSolidBoundary.hh b/src/DEM/SolidBoundary/SphereSolidBoundary.hh index 5db055556..3d260aee7 100644 --- a/src/DEM/SolidBoundary/SphereSolidBoundary.hh +++ b/src/DEM/SolidBoundary/SphereSolidBoundary.hh @@ -36,8 +36,7 @@ public: virtual Vector velocity(const Vector& position) const override; virtual void registerState(DataBase& dataBase, - State& state, - const std::string& boundaryKey) override; + State& state) override; virtual void update(const double multiplier, diff --git a/src/PYB11/DEM/SolidBoundaries.py b/src/PYB11/DEM/SolidBoundaries.py index 31509ebdd..621c6a09a 100644 --- a/src/PYB11/DEM/SolidBoundaries.py +++ b/src/PYB11/DEM/SolidBoundaries.py @@ -36,8 +36,7 @@ def pyinit(point = "const Vector&", @PYB11virtual def registerState(self, dataBase = "DataBase<%(Dimension)s>&", - state = "State<%(Dimension)s>&", - boundaryKey = "const std::string&"): + state = "State<%(Dimension)s>&"): "Register the state solid bc expects to use and evolve." return "void" @@ -88,8 +87,7 @@ def pyinit(point = "const Vector&", @PYB11virtual def registerState(self, dataBase = "DataBase<%(Dimension)s>&", - state = "State<%(Dimension)s>&", - boundaryKey = "const std::string&"): + state = "State<%(Dimension)s>&"): "Register the state solid bc expects to use and evolve." return "void" @@ -140,8 +138,7 @@ def pyinit(point = "const Vector&", @PYB11virtual def registerState(self, dataBase = "DataBase<%(Dimension)s>&", - state = "State<%(Dimension)s>&", - boundaryKey = "const std::string&"): + state = "State<%(Dimension)s>&"): "Register the state solid bc expects to use and evolve." return "void" @@ -193,8 +190,7 @@ def pyinit(point = "const Vector&", @PYB11virtual def registerState(self, dataBase = "DataBase<%(Dimension)s>&", - state = "State<%(Dimension)s>&", - boundaryKey = "const std::string&"): + state = "State<%(Dimension)s>&"): "Register the state solid bc expects to use and evolve." return "void" @@ -248,8 +244,7 @@ def pyinit(center = "const Vector&", @PYB11virtual def registerState(self, dataBase = "DataBase<%(Dimension)s>&", - state = "State<%(Dimension)s>&", - boundaryKey = "const std::string&"): + state = "State<%(Dimension)s>&"): "Register the state solid bc expects to use and evolve." return "void" diff --git a/src/PYB11/DEM/SolidBoundaryBaseAbstractMethods.py b/src/PYB11/DEM/SolidBoundaryBaseAbstractMethods.py index 777c5d0fc..e7278fe1a 100644 --- a/src/PYB11/DEM/SolidBoundaryBaseAbstractMethods.py +++ b/src/PYB11/DEM/SolidBoundaryBaseAbstractMethods.py @@ -20,8 +20,7 @@ def distance(self, def registerState(self, dataBase = "DataBase<%(Dimension)s>&", - state = "State<%(Dimension)s>&", - boundaryKey = "const std::string&"): + state = "State<%(Dimension)s>&"): "Register the state solid bc expects to use and evolve." return "void" From ffda8daf95d255f178ec5ccc4237afce285ef0f8 Mon Sep 17 00:00:00 2001 From: jmpearl Date: Sat, 2 Mar 2024 12:39:29 -0800 Subject: [PATCH 11/38] DEM cleaning out old commented out sections --- src/DEM/DEMBase.cc | 6 ------ src/PYB11/DEM/DEMBase.py | 8 +------- 2 files changed, 1 insertion(+), 13 deletions(-) diff --git a/src/DEM/DEMBase.cc b/src/DEM/DEMBase.cc index 604ebc8d8..03867b6e2 100644 --- a/src/DEM/DEMBase.cc +++ b/src/DEM/DEMBase.cc @@ -321,12 +321,6 @@ registerState(DataBase& dataBase, for (auto ibc = 0u; ibc < this->numSolidBoundaries(); ++ibc){ mSolidBoundaries[ibc]->registerState(dataBase,state);} - //const auto name = "SolidBoundary_" + std::to_string(solidBoundaries[ibc]->uniqueIndex()); - // for (ConstSolidBoundaryIterator solidboundItr = mSolidBoundaries.begin(); - // solidboundItr != mSolidBoundaries.end(); - // ++solidboundItr){ - // (*solidboundItr)->registerState(dataBase,state); - // } state.enroll(mTimeStepMask); state.enroll(mass); diff --git a/src/PYB11/DEM/DEMBase.py b/src/PYB11/DEM/DEMBase.py index 3f2620e59..3cc312c69 100644 --- a/src/PYB11/DEM/DEMBase.py +++ b/src/PYB11/DEM/DEMBase.py @@ -58,8 +58,7 @@ def initialize(time = "const Scalar", state = "State<%(Dimension)s>&", derivs = "StateDerivatives<%(Dimension)s>&"): "Initialize the DEM before we start a derivative evaluation." - return "void" - + return "void" @PYB11virtual @PYB11const @@ -104,11 +103,6 @@ def clearSolidBoundaries(self): "remove all solid boundaries from the dem package" return "void" - # @PYB11const - # def numSolidBoundaries(self): - # "return the number of solid boundaries being tracked" - # return "unsigned int" - @PYB11const def haveSolidBoundary(boundary = "const SolidBoundaryBase<%(Dimension)s>&"): "is this boundary being tracked?" From f163ecd44c82bae83e45def82c65fac2c70a1e97 Mon Sep 17 00:00:00 2001 From: jmpearl Date: Sat, 2 Mar 2024 12:40:43 -0800 Subject: [PATCH 12/38] DEM more cleaning --- src/DEM/DEMBaseInline.hh | 9 +-------- 1 file changed, 1 insertion(+), 8 deletions(-) diff --git a/src/DEM/DEMBaseInline.hh b/src/DEM/DEMBaseInline.hh index 0a7f3b94d..f716956dc 100644 --- a/src/DEM/DEMBaseInline.hh +++ b/src/DEM/DEMBaseInline.hh @@ -370,15 +370,8 @@ inline void DEMBase:: removeSolidBoundary(const SolidBoundaryBase& boundary) { - const auto bcPtr = std::find(mSolidBoundaries.begin(),mSolidBoundaries.end(),&boundary); - if (bcPtr != mSolidBoundaries.end()){ - //const int bcIndex = std::distance(mSolidBoundaries.begin(),bcPtr); - mSolidBoundaries.erase(bcPtr); - //std::cout << "removing solid bc " << bcIndex << std::endl; - } else { - std:: cout << "THEY AINT HOME"; - } + if (bcPtr != mSolidBoundaries.end()) mSolidBoundaries.erase(bcPtr); } From c50e00b8526c573b20aed63aeef4cf9966fefc7a Mon Sep 17 00:00:00 2001 From: jmpearl Date: Sun, 3 Mar 2024 07:49:39 -0800 Subject: [PATCH 13/38] DEM cleaning out some print diagnostics --- src/DataBase/StateBase.cc | 1 - 1 file changed, 1 deletion(-) diff --git a/src/DataBase/StateBase.cc b/src/DataBase/StateBase.cc index 42111c6f0..cae66eab2 100644 --- a/src/DataBase/StateBase.cc +++ b/src/DataBase/StateBase.cc @@ -401,7 +401,6 @@ assign(const StateBase& rhs) { *lhsptr = *rhsptr; } catch(const boost::bad_any_cast&) { try { - std::cout << "key= " << itr -> first << std::endl; auto lhsptr = boost::any_cast(anylhs); const auto rhsptr = boost::any_cast(anyrhs); *lhsptr = *rhsptr; From bd522f9a51b9b484df78377b6ebe7fbc3ac46eb4 Mon Sep 17 00:00:00 2001 From: jmpearl Date: Sun, 3 Mar 2024 07:50:58 -0800 Subject: [PATCH 14/38] DEM: adding field to track max overlap --- src/DEM/DEMFieldNames.cc | 1 + src/DEM/DEMFieldNames.hh | 1 + src/DEM/LinearSpringDEM.cc | 47 ++++++++++++++++++-- src/DEM/LinearSpringDEM.hh | 7 +++ src/DEM/LinearSpringDEMInline.hh | 16 +++++++ src/DataBase/CMakeLists.txt | 1 + src/DataBase/MaxReplaceState.hh | 51 ++++++++++++++++++++++ src/DataBase/MaxReplaceStateInline.hh | 62 +++++++++++++++++++++++++++ 8 files changed, 183 insertions(+), 3 deletions(-) create mode 100644 src/DataBase/MaxReplaceState.hh create mode 100644 src/DataBase/MaxReplaceStateInline.hh diff --git a/src/DEM/DEMFieldNames.cc b/src/DEM/DEMFieldNames.cc index ae4b05361..b0420b608 100644 --- a/src/DEM/DEMFieldNames.cc +++ b/src/DEM/DEMFieldNames.cc @@ -18,5 +18,6 @@ const std::string Spheral::DEMFieldNames::shearDisplacement = "shear displacemen const std::string Spheral::DEMFieldNames::rollingDisplacement = "rolling displacement"; const std::string Spheral::DEMFieldNames::torsionalDisplacement = "torsional displacement"; const std::string Spheral::DEMFieldNames::equilibriumOverlap = "equilibrium overlap"; +const std::string Spheral::DEMFieldNames::maximumOverlap = "maximum overlap"; const std::string Spheral::DEMFieldNames::solidBoundaries = "solid boundaries"; const std::string Spheral::DEMFieldNames::solidBoundaryPolicy = "solid boundary policy"; \ No newline at end of file diff --git a/src/DEM/DEMFieldNames.hh b/src/DEM/DEMFieldNames.hh index 4acdafb26..3358a9925 100644 --- a/src/DEM/DEMFieldNames.hh +++ b/src/DEM/DEMFieldNames.hh @@ -23,6 +23,7 @@ struct DEMFieldNames { static const std::string rollingDisplacement; static const std::string torsionalDisplacement; static const std::string equilibriumOverlap; + static const std::string maximumOverlap; static const std::string solidBoundaries; static const std::string solidBoundaryPolicy; }; diff --git a/src/DEM/LinearSpringDEM.cc b/src/DEM/LinearSpringDEM.cc index bed1bc730..c591e1197 100644 --- a/src/DEM/LinearSpringDEM.cc +++ b/src/DEM/LinearSpringDEM.cc @@ -18,7 +18,7 @@ #include "DataBase/StateDerivatives.hh" #include "DataBase/DataBase.hh" #include "DataBase/IncrementState.hh" -#include "DataBase/ReplaceState.hh" +#include "DataBase/MaxReplaceState.hh" #include "Field/FieldList.hh" #include "Neighbor/ConnectivityMap.hh" @@ -95,8 +95,12 @@ LinearSpringDEM(const DataBase& dataBase, mShapeFactor(shapeFactor), mNormalBeta(M_PI/std::log(std::max(normalRestitutionCoefficient,1.0e-3))), mTangentialBeta(M_PI/std::log(std::max(tangentialRestitutionCoefficient,1.0e-3))), - mMomentOfInertia(FieldStorageType::CopyFields) { + mMomentOfInertia(FieldStorageType::CopyFields), + mMaximumOverlap(FieldStorageType::CopyFields), + mNewMaximumOverlap(FieldStorageType::CopyFields) { mMomentOfInertia = dataBase.newDEMFieldList(0.0, DEMFieldNames::momentOfInertia); + mMaximumOverlap = dataBase.newDEMFieldList(0.0, DEMFieldNames::maximumOverlap); + mNewMaximumOverlap = dataBase.newDEMFieldList(0.0,MaxReplaceState::prefix() + DEMFieldNames::maximumOverlap); } //------------------------------------------------------------------------------ @@ -299,11 +303,39 @@ LinearSpringDEM:: registerState(DataBase& dataBase, State& state) { TIME_BEGIN("LinearSpringDEMregisterState"); + DEMBase::registerState(dataBase,state); + dataBase.resizeDEMFieldList(mMomentOfInertia, 0.0, DEMFieldNames::momentOfInertia, false); + dataBase.resizeDEMFieldList(mMaximumOverlap, 0.0, DEMFieldNames::maximumOverlap, false); + + auto maxOverlapPolicy = make_policy>(); + state.enroll(mMomentOfInertia); + state.enroll(mMaximumOverlap, maxOverlapPolicy); + TIME_END("LinearSpringDEMregisterState"); } + +//------------------------------------------------------------------------------ +// Register the state we need/are going to evolve. +//------------------------------------------------------------------------------ +template +void +LinearSpringDEM:: +registerDerivatives(DataBase& dataBase, + StateDerivatives& derivs) { + TIME_BEGIN("LinearSpringDEMregisterDerivs"); + + DEMBase::registerDerivatives(dataBase,derivs); + + dataBase.resizeDEMFieldList(mNewMaximumOverlap, 0.0, DEMFieldNames::maximumOverlap, false); + + derivs.enroll(mNewMaximumOverlap); + + TIME_END("LinearSpringDEMregisterDerivs"); +} + //------------------------------------------------------------------------------ // evaluate the derivatives //------------------------------------------------------------------------------ @@ -387,10 +419,12 @@ evaluateDerivatives(const typename Dimension::Scalar /*time*/, auto DxDt = derivatives.fields(IncrementState::prefix() + HydroFieldNames::position, Vector::zero); auto DvDt = derivatives.fields(HydroFieldNames::hydroAcceleration, Vector::zero); auto DomegaDt = derivatives.fields(IncrementState::prefix() + DEMFieldNames::angularVelocity, DEMDimension::zero); + auto newMaximumOverlap = derivatives.fields(MaxReplaceState::prefix() + DEMFieldNames::maximumOverlap, 0.0); CHECK(DxDt.size() == numNodeLists); CHECK(DvDt.size() == numNodeLists); CHECK(DomegaDt.size() == numNodeLists); + CHECK(newMaximumOverlap.size() == numNodeLists); // Get the deriv pairFieldLists auto DDtShearDisplacement = derivatives.fields(ReplaceAndIncrementPairFieldList>::incrementPrefix() + DEMFieldNames::shearDisplacement, std::vector()); @@ -421,6 +455,7 @@ evaluateDerivatives(const typename Dimension::Scalar /*time*/, typename SpheralThreads::FieldListStack threadStack; auto DvDt_thread = DvDt.threadCopy(threadStack); auto DomegaDt_thread = DomegaDt.threadCopy(threadStack); + auto newMaxOverlap_thread = newMaximumOverlap.threadCopy(threadStack, ThreadReduction::MAX); //------------------------------------ // particle-particle contacts @@ -458,7 +493,7 @@ evaluateDerivatives(const typename Dimension::Scalar /*time*/, // if so do the things if (delta0 > 0.0){ - + // get remaining state for node i const auto cIdi = compositeIndex(nodeListi,i); const auto uIdi = uniqueIndices(nodeListi,i); @@ -482,10 +517,12 @@ evaluateDerivatives(const typename Dimension::Scalar /*time*/, // Get the derivs from node i auto& DvDti = DvDt_thread(nodeListi, i); auto& DomegaDti = DomegaDt_thread(nodeListi, i); + auto& maxOverlapi = newMaxOverlap_thread(nodeListi,i); // Get the derivs from node j auto& DvDtj = DvDt_thread(nodeListj, j); auto& DomegaDtj = DomegaDt_thread(nodeListj, j); + auto& maxOverlapj = newMaxOverlap_thread(nodeListj,j); // storage sign, this makes pairwise values i-j independent const int storageSign = (uIdi <= uIdj ? 1 : -1); @@ -608,6 +645,10 @@ evaluateDerivatives(const typename Dimension::Scalar /*time*/, DomegaDti += (Msliding*li - (Mtorsion + Mrolling) * lij)/Ii; DomegaDtj += (Msliding*lj + (Mtorsion + Mrolling) * lij)/Ij; + // update max overlaps + maxOverlapi = max(maxOverlapi,delta); + maxOverlapj = max(maxOverlapj,delta); + // for spring updates newShearDisplacement(nodeListi,i)[contacti] = storageSign*newDeltaSlidij; DDtShearDisplacement(nodeListi,i)[contacti] = storageSign*vs; diff --git a/src/DEM/LinearSpringDEM.hh b/src/DEM/LinearSpringDEM.hh index b932dc7d5..a2bf7fecb 100644 --- a/src/DEM/LinearSpringDEM.hh +++ b/src/DEM/LinearSpringDEM.hh @@ -65,6 +65,9 @@ public: virtual void registerState(DataBase& dataBase, State& state) override; + virtual void registerDerivatives(DataBase& dataBase, + StateDerivatives& derivs) override; + virtual void evaluateDerivatives(const Scalar time, const Scalar dt, const DataBase& dataBase, @@ -123,6 +126,8 @@ public: // get methods for class FieldLists const FieldList& momentOfInertia() const; + const FieldList& maximumOverlap() const; + const FieldList& newMaximumOverlap() const; //**************************************************************************** @@ -149,6 +154,8 @@ private: // field Lists FieldList mMomentOfInertia; + FieldList mMaximumOverlap; + FieldList mNewMaximumOverlap; // FieldList mOptimalSpringConstant; // No default constructor, copying, or assignment. diff --git a/src/DEM/LinearSpringDEMInline.hh b/src/DEM/LinearSpringDEMInline.hh index f87e570b5..4451a3d41 100644 --- a/src/DEM/LinearSpringDEMInline.hh +++ b/src/DEM/LinearSpringDEMInline.hh @@ -257,4 +257,20 @@ momentOfInertia() const { return mMomentOfInertia; } +template +inline +const FieldList& +LinearSpringDEM:: +maximumOverlap() const { + return mMaximumOverlap; +} + +template +inline +const FieldList& +LinearSpringDEM:: +newMaximumOverlap() const { + return mNewMaximumOverlap; +} + } diff --git a/src/DataBase/CMakeLists.txt b/src/DataBase/CMakeLists.txt index f7cf8171d..e6dc6cee5 100644 --- a/src/DataBase/CMakeLists.txt +++ b/src/DataBase/CMakeLists.txt @@ -25,6 +25,7 @@ set(DataBase_headers IncrementStateInline.hh IncrementBoundedState.hh IncrementBoundedStateInline.hh + MaxReplaceState.hh PureReplaceState.hh PureReplaceStateInline.hh ReplaceState.hh diff --git a/src/DataBase/MaxReplaceState.hh b/src/DataBase/MaxReplaceState.hh new file mode 100644 index 000000000..099a946a9 --- /dev/null +++ b/src/DataBase/MaxReplaceState.hh @@ -0,0 +1,51 @@ +//---------------------------------Spheral++----------------------------------// +// MaxReplaceState -- Replaces the state with the max of the state/deriv value +// +// J.M. Pearl 2024 +//----------------------------------------------------------------------------// +#ifndef __Spheral_MaxReplaceState_hh__ +#define __Spheral_MaxReplaceState_hh__ + +#include "FieldUpdatePolicy.hh" +#include "Utilities/DBC.hh" + +namespace Spheral { + +// Forward declarations. +template class StateDerivatives; + +template +class MaxReplaceState: public FieldUpdatePolicy { +public: + //--------------------------- Public Interface ---------------------------// + // Useful typedefs + using KeyType = typename FieldUpdatePolicy::KeyType; + + // Constructors, destructor. + MaxReplaceState(std::initializer_list depends = {}); + virtual ~MaxReplaceState() {} + + // Overload the methods describing how to update Fields. + virtual void update(const KeyType& key, + State& state, + StateDerivatives& derivs, + const double multiplier, + const double t, + const double dt) override; + + // Equivalence. + virtual bool operator==(const UpdatePolicyBase& rhs) const override; + + static const std::string prefix() { return "new "; } + +private: + //--------------------------- Private Interface ---------------------------// + MaxReplaceState(const MaxReplaceState& rhs); + MaxReplaceState& operator=(const MaxReplaceState& rhs); +}; + +} + +#include "MaxReplaceStateInline.hh" + +#endif diff --git a/src/DataBase/MaxReplaceStateInline.hh b/src/DataBase/MaxReplaceStateInline.hh new file mode 100644 index 000000000..a69e0c813 --- /dev/null +++ b/src/DataBase/MaxReplaceStateInline.hh @@ -0,0 +1,62 @@ +//---------------------------------Spheral++----------------------------------// +// MaxReplaceState -- Replaces the state with the max of the state/deriv value +// +// J.M. Pearl 2024 +//----------------------------------------------------------------------------// +#include "State.hh" +#include "StateDerivatives.hh" +#include "Field/Field.hh" +#include "Utilities/DBC.hh" + +namespace Spheral { + +//------------------------------------------------------------------------------ +// Constructors. +//------------------------------------------------------------------------------ +template +inline +MaxReplaceState:: +MaxReplaceState(std::initializer_list depends): + FieldUpdatePolicy(depends) { +} + +//------------------------------------------------------------------------------ +// Update the field. +//------------------------------------------------------------------------------ +template +inline +void +MaxReplaceState:: +update(const KeyType& key, + State& state, + StateDerivatives& derivs, + const double /*multiplier*/, + const double /*t*/, + const double /*dt*/) { + + // Find the matching replacement field from the StateDerivatives. + const auto replaceKey = prefix() + key; + auto& f = state.field(key, ValueType()); + const auto& df = derivs.field(replaceKey, ValueType()); + + // Loop over the internal values of the field. + const auto n = f.nodeList().numInternalNodes(); +#pragma omp parallel for + for (auto i = 0u; i < n; ++i) { + f(i) = std::max(df(i),f(i)); + } +} + +//------------------------------------------------------------------------------ +// Equivalence operator. +//------------------------------------------------------------------------------ +template +inline +bool +MaxReplaceState:: +operator==(const UpdatePolicyBase& rhs) const { + return dynamic_cast*>(&rhs) != nullptr; +} + +} + From 8edbf2bbc850a2ff373ec74928d739af75561e59 Mon Sep 17 00:00:00 2001 From: jmpearl Date: Sun, 3 Mar 2024 21:12:19 -0800 Subject: [PATCH 15/38] LinearSpringDEM replacing redundent code with inline functions --- src/DEM/LinearSpringDEM.cc | 142 ++++++------------------------- src/DEM/LinearSpringDEM.hh | 37 ++++++-- src/DEM/LinearSpringDEMInline.hh | 63 ++++++++++++++ 3 files changed, 122 insertions(+), 120 deletions(-) diff --git a/src/DEM/LinearSpringDEM.cc b/src/DEM/LinearSpringDEM.cc index c591e1197..479e201d8 100644 --- a/src/DEM/LinearSpringDEM.cc +++ b/src/DEM/LinearSpringDEM.cc @@ -357,8 +357,8 @@ evaluateDerivatives(const typename Dimension::Scalar /*time*/, const auto muD = this->dynamicFrictionCoefficient(); const auto muS = this->staticFrictionCoefficient(); - const auto muT = this->torsionalFrictionCoefficient(); - const auto muR = this->rollingFrictionCoefficient(); + const auto muT = this->torsionalFrictionCoefficient() * shapeFactor * muS; + const auto muR = this->rollingFrictionCoefficient() * shapeFactor; const auto Cc = this->cohesiveTensileStrength(); @@ -535,6 +535,7 @@ evaluateDerivatives(const typename Dimension::Scalar /*time*/, // boolean checks const auto isBondedParticle = (cIdi == cIdj); + const auto allowSliding = (!isBondedParticle); // effective delta const auto delta = delta0 - overlapij; @@ -562,8 +563,8 @@ evaluateDerivatives(const typename Dimension::Scalar /*time*/, const Vector vij = vi-vj + li*vroti - lj*vrotj; - const Scalar vn = vij.dot(rhatij); // normal velocity - const Vector vs = vij - vn*rhatij; // sliding velocity + const Scalar vn = vij.dot(rhatij); // normal velocity + const Vector vs = vij - vn*rhatij; // sliding velocity const Vector vr = -li*vroti - lj*vrotj; // rolling velocity const Scalar vt = -lij*DEMDimension::dot(omegai-omegaj,rhatij); // torsion velocity @@ -575,73 +576,32 @@ evaluateDerivatives(const typename Dimension::Scalar /*time*/, // sliding //------------------------------------------------------------ - // project onto new tangential plane -- maintain magnitude - Vector newDeltaSlidij = (deltaSlidij - rhatij.dot(deltaSlidij)*rhatij).unitVector()*deltaSlidij.magnitude(); - - // spring dashpot - const Vector ft0spring = - ks*newDeltaSlidij; - const Vector ft0damp = - Cs*vs; - Vector ft = ft0spring + ft0damp; - - // static friction limit - if (!isBondedParticle and (ft.magnitude() > muS*fnMag)){ - - const Scalar ftDynamic = muD*fnMag; - ft = ftDynamic*ft.unitVector(); - - newDeltaSlidij = ( ft0damp.magnitude() > ftDynamic ? - Vector::zero : - -(ft-ft0damp)*invKs ); - - } + Vector newDeltaSlidij, fs; + this->slidingSpringDamper(ks,Cs,muS,muD,deltaSlidij,vs,fnMag,invKs,rhatij,allowSliding, + newDeltaSlidij,fs); // outputs // torsion //------------------------------------------------------------ - // since we use a scalar no need to modify here - auto newDeltaTorsij = deltaTorsij; - - // spring dashpot - const Scalar Mt0spring = - kt*newDeltaTorsij; - const Scalar Mt0damp = - Ct*vt; - Scalar MtorsionMag = (Mt0spring + Mt0damp); - const Scalar MtStatic = muT*shapeFactor*muS*fnMag; - - // limit to static - if (!isBondedParticle and (std::abs(MtorsionMag) > MtStatic)){ - MtorsionMag = (MtorsionMag > 0.0 ? 1.0 : -1.0)*MtStatic; - newDeltaTorsij = (std::abs(Mt0damp) > MtStatic ? 0.0 : -(MtorsionMag-Mt0damp)*invKt); - } + Scalar newDeltaTorsij, ft; + this->slidingSpringDamper(kt,Ct,muT,muT,deltaTorsij,vt,fnMag,invKt,allowSliding, + newDeltaTorsij,ft); // output // rolling //------------------------------------------------------------ - // project onto new tangential plane -- maintain magnitude - Vector newDeltaRollij = (deltaRollij - rhatij.dot(deltaRollij)*rhatij).unitVector()*deltaRollij.magnitude(); - - // spring dashpot - const Vector Mr0spring = - kr*newDeltaRollij; - const Vector Mr0damp = - Cr*vr; - Vector effectiveRollingForce = (Mr0spring + Mr0damp); - const Scalar MrStatic = muR*shapeFactor*fnMag; - - // limit to static - if (!isBondedParticle and (effectiveRollingForce.magnitude() > MrStatic)){ - effectiveRollingForce = MrStatic*effectiveRollingForce.unitVector(); - newDeltaRollij = (Mr0damp.magnitude() > MrStatic ? - Vector::zero : - -(effectiveRollingForce-Mr0damp)*invKr); - } - + Vector newDeltaRollij, fr; + this->slidingSpringDamper(kr,Cr,muR,muR,deltaRollij,vr,fnMag,invKr,rhatij,allowSliding, + newDeltaRollij,fr); // outputs // accelerations //------------------------------------------------------------ // Rectilinear Acceleration - const Vector fij = fn - fc + ft; + const Vector fij = fn - fc + fs; DvDti += fij/mi; DvDtj -= fij/mj; // angular acceleration const auto Msliding = -DEMDimension::cross(rhatij,fij); - const auto Mrolling = -DEMDimension::cross(rhatij,effectiveRollingForce); - const auto Mtorsion = MtorsionMag * this->torsionMoment(rhatij,omegai,omegaj); // rename torsionDirection + const auto Mrolling = -DEMDimension::cross(rhatij,fr); + const auto Mtorsion = ft * this->torsionMoment(rhatij,omegai,omegaj); // rename torsionDirection DomegaDti += (Msliding*li - (Mtorsion + Mrolling) * lij)/Ii; DomegaDtj += (Msliding*lj + (Mtorsion + Mrolling) * lij)/Ij; @@ -725,12 +685,7 @@ evaluateDerivatives(const typename Dimension::Scalar /*time*/, const auto mib = 2*mi; const auto lib = 2*li; - // damping constants -- ct and cr derived quantities ala Zhang 2017 - //----------------------------To-Do----------------------------------- - // oh man okay, solid bc's were overdamping. The 0.5 factor got the - // coeff of restitution right (i guess i missed a factor of 2 somewhere) - // need to go back and do the math later. - //-------------------------------------------------------------------- + // damping coefficients const auto Cn = std::sqrt(mib*normalDampingTerms); const auto Cs = std::sqrt(mib*tangentialDampingTerms); const auto Ct = 0.50 * Cs * shapeFactor2; @@ -754,62 +709,21 @@ evaluateDerivatives(const typename Dimension::Scalar /*time*/, // sliding //------------------------------------------------------------ - // project onto new tangential plane -- maintain magnitude - Vector newDeltaSlidib = (deltaSlidib - rhatib.dot(deltaSlidib)*rhatib).unitVector()*deltaSlidib.magnitude(); - - // spring dashpot - const Vector ft0spring = - ks*newDeltaSlidib; - const Vector ft0damp = - Cs*vs; - Vector ft = ft0spring + ft0damp; - - // static friction limit - if (ft.magnitude() > muS*fnMag){ - - const Scalar ftDynamic = muD*fnMag; - ft = ftDynamic*ft.unitVector(); - - newDeltaSlidib = ( ft0damp.magnitude() > ftDynamic ? - Vector::zero : - -(ft-ft0damp)*invKs ); - - } + Vector newDeltaSlidib, ft; + this->slidingSpringDamper(ks,Cs,muS,muD,deltaSlidib,vs,fnMag,invKs,rhatib,true, + newDeltaSlidib,ft); // outputs // torsion //------------------------------------------------------------ - // since we use a scalar no need to modify here - auto newDeltaTorsib = deltaTorsib; - - // spring dashpot - const Scalar Mt0spring = - kt*newDeltaTorsib; - const Scalar Mt0damp = - Ct*vt; - Scalar MtorsionMag = (Mt0spring + Mt0damp); - const Scalar MtStatic = muT*shapeFactor*muS*fnMag; - - // limit to static - if (std::abs(MtorsionMag) > MtStatic){ - MtorsionMag = (MtorsionMag > 0.0 ? 1.0 : -1.0)*MtStatic; - newDeltaTorsib = (std::abs(Mt0damp) > MtStatic ? 0.0 : -(MtorsionMag-Mt0damp)*invKt); - } + Scalar newDeltaTorsib, MtorsionMag; + this->slidingSpringDamper(kt,Ct,muT,muT,deltaTorsib,vt,fnMag,invKt,true, + newDeltaTorsib, MtorsionMag); // output // rolling //------------------------------------------------------------ - // project onto new tangential plane -- maintain magnitude - Vector newDeltaRollib = (deltaRollib - rhatib.dot(deltaRollib)*rhatib).unitVector()*deltaRollib.magnitude(); - - // spring dashpot - const Vector Mr0spring = - kr*newDeltaRollib; - const Vector Mr0damp = - Cr*vr; - Vector effectiveRollingForce = (Mr0spring + Mr0damp); - const Scalar MrStatic = muR*shapeFactor*fnMag; - - // limit to static - if (effectiveRollingForce.magnitude() > MrStatic){ - effectiveRollingForce = MrStatic*effectiveRollingForce.unitVector(); - newDeltaRollib = (Mr0damp.magnitude() > MrStatic ? - Vector::zero : - -(effectiveRollingForce-Mr0damp)*invKr); - } - + Vector newDeltaRollib, fr; + this->slidingSpringDamper(kr,Cr,muR,muR,deltaRollib,vr,fnMag,invKr,rhatib,true, + newDeltaRollib, fr); // outputs // accelerations //------------------------------------------------------------ // Rectilinear Acceleration @@ -818,7 +732,7 @@ evaluateDerivatives(const typename Dimension::Scalar /*time*/, // angular acceleration const auto Msliding = -DEMDimension::cross(rhatib,fib); - const auto Mrolling = -DEMDimension::cross(rhatib,effectiveRollingForce); + const auto Mrolling = -DEMDimension::cross(rhatib,fr); const auto Mtorsion = MtorsionMag * this->torsionMoment(rhatib,omegai,0*omegai); // rename torsionDirection DomegaDti += (Msliding*li - (Mtorsion + Mrolling) * lib)/Ii; diff --git a/src/DEM/LinearSpringDEM.hh b/src/DEM/LinearSpringDEM.hh index a2bf7fecb..290b8fa5d 100644 --- a/src/DEM/LinearSpringDEM.hh +++ b/src/DEM/LinearSpringDEM.hh @@ -73,12 +73,37 @@ public: const DataBase& dataBase, const State& state, StateDerivatives& derivs) const override; - virtual - void applyGhostBoundaries(State& state, - StateDerivatives& derivs) override; - virtual - void enforceBoundaries(State& state, - StateDerivatives& derivs) override; + + virtual void applyGhostBoundaries(State& state, + StateDerivatives& derivs) override; + + virtual void enforceBoundaries(State& state, + StateDerivatives& derivs) override; + + void slidingSpringDamper(const Scalar k, + const Scalar C, + const Scalar mus, + const Scalar mud, + const Vector& x, + const Vector& DxDt, + const Scalar fnMag, + const Scalar invK, + const Vector& rhatij, + const bool allowSiding, + Vector& xNew, + Vector& force) const; + + void slidingSpringDamper(const Scalar k, + const Scalar C, + const Scalar mus, + const Scalar mud, + const Scalar x, + const Scalar DxDt, + const Scalar fnMag, + const Scalar invK, + const bool allowSiding, + Scalar& xNew, + Scalar& force) const; // set/gets Scalar normalSpringConstant() const; diff --git a/src/DEM/LinearSpringDEMInline.hh b/src/DEM/LinearSpringDEMInline.hh index 4451a3d41..6668bfe99 100644 --- a/src/DEM/LinearSpringDEMInline.hh +++ b/src/DEM/LinearSpringDEMInline.hh @@ -246,6 +246,69 @@ momentOfInertia(const Dim<3>::Scalar m, const Dim<3>::Scalar R) const { } +//------------------------------------------------------------------------------ +// friction functions +//------------------------------------------------------------------------------ +template +inline +void +LinearSpringDEM:: +slidingSpringDamper(const typename Dimension::Scalar k, + const typename Dimension::Scalar C, + const typename Dimension::Scalar mus, + const typename Dimension::Scalar mud, + const typename Dimension::Vector& x, + const typename Dimension::Vector& DxDt, + const typename Dimension::Scalar fnMag, + const typename Dimension::Scalar invK, + const typename Dimension::Vector& rhatij, + const bool allowSliding, + typename Dimension::Vector& xNew, + typename Dimension::Vector& forceTotal) const{ + + xNew = (x - rhatij.dot(x)*rhatij).unitVector()*x.magnitude(); + + const Vector forceSpring = - k * xNew; + const Vector forceDamper = - C * DxDt; + + forceTotal = forceSpring + forceDamper; + + if (allowSliding and (forceTotal.magnitude() > mus * fnMag)){ + forceTotal = mud*fnMag*forceTotal.unitVector(); + xNew = (forceDamper.magnitude() > mud*fnMag ? + Vector::zero : + -(forceTotal-forceDamper)*invK ); + } +} + +template +inline +void +LinearSpringDEM:: +slidingSpringDamper(const typename Dimension::Scalar k, + const typename Dimension::Scalar C, + const typename Dimension::Scalar mus, + const typename Dimension::Scalar mud, + const typename Dimension::Scalar x, + const typename Dimension::Scalar DxDt, + const typename Dimension::Scalar fnMag, + const typename Dimension::Scalar invK, + const bool allowSliding, + typename Dimension::Scalar& xNew, + typename Dimension::Scalar& forceTotal) const{ + xNew = x; + + const Scalar forceSpring = - k * xNew; + const Scalar forceDamper = - C * DxDt; + + forceTotal = forceSpring + forceDamper; + + if (allowSliding and (std::abs(forceTotal) > mus * fnMag)){ + forceTotal = (forceTotal > 0.0 ? 1.0 : -1.0) * mud * fnMag; + xNew = (std::abs(forceDamper) > mud * fnMag ? 0.0 : -(forceTotal-forceDamper)*invK); + } +} + //------------------------------------------------------------------------------ // FieldList //------------------------------------------------------------------------------ From c9093660f7316b01a6b3105d03800e431c110820 Mon Sep 17 00:00:00 2001 From: jmpearl Date: Sun, 3 Mar 2024 21:17:22 -0800 Subject: [PATCH 16/38] DEM adding ats tests for rolling and torsion on surfaces --- .../rollingParticleTest.py | 302 +++++++++++++++++ .../torsionalParticleTest.py | 313 ++++++++++++++++++ 2 files changed, 615 insertions(+) create mode 100644 tests/functional/DEM/LinearSpringDEM/SolidBoundaryCondition/rollingParticleTest.py create mode 100644 tests/functional/DEM/LinearSpringDEM/SolidBoundaryCondition/torsionalParticleTest.py diff --git a/tests/functional/DEM/LinearSpringDEM/SolidBoundaryCondition/rollingParticleTest.py b/tests/functional/DEM/LinearSpringDEM/SolidBoundaryCondition/rollingParticleTest.py new file mode 100644 index 000000000..8a58078e1 --- /dev/null +++ b/tests/functional/DEM/LinearSpringDEM/SolidBoundaryCondition/rollingParticleTest.py @@ -0,0 +1,302 @@ +#ATS:DEM3dImpact = test( SELF, "--clearDirectories True --checkConservation True --goalTime 1.0", label="DEM impacting squares -- 3-D (parallel)", np=8) + +import os, sys, shutil, mpi, random +from math import * +from Spheral3d import * +from SpheralTestUtilities import * +from findLastRestart import * +from GenerateNodeDistribution3d import * +from GenerateDEMfromSPHGenerator import GenerateDEMfromSPHGenerator3d + +import numpy as np + +sys.path.insert(0, '..') +from DEMConservationTracker import TrackConservation3d as TrackConservation + +if mpi.procs > 1: + from PeanoHilbertDistributeNodes import distributeNodes3d +else: + from DistributeNodes import distributeNodes3d + +title("DEM 3d Drop Test with Particle Generation") +# this tests the ability to generate particle on the fly +# during the course of a simulation using a periodic +# work function. It also tests the solid boundary condition + +#------------------------------------------------------------------------------- +# Generic problem parameters +#------------------------------------------------------------------------------- +commandLine(numParticlePerLength = 3, # number of particles on a side of the box + normalSpringConstant=1.0, # spring constant for LDS model + normalRestitutionCoefficient=0.55, # restitution coefficient to get damping const + tangentialSpringConstant=0.2857, # spring constant for LDS model + tangentialRestitutionCoefficient=0.55, # restitution coefficient to get damping const + dynamicFriction = 1.0, # static friction coefficient sliding + staticFriction = 1.0, # dynamic friction coefficient sliding + rollingFriction = 1.05, # static friction coefficient for rolling + torsionalFriction = 1.3, # static friction coefficient for torsion + cohesiveTensileStrength = 0.0, # units of pressure + shapeFactor = 0.1, # in [0,1] shape factor from Zhang 2018, 0 - no torsion or rolling + + particleRadius = 0.10, # particle radius + particleDensity = 2.60, + particleVelocity = 0.1, + + neighborSearchBuffer = 0.1, # multiplicative buffer to radius for neighbor search algo + nPerh = 1.01, + + # integration + IntegratorConstructor = VerletIntegrator, + stepsPerCollision = 50, # replaces CFL for DEM + goalTime = 25.0, + dt = 1e-8, + dtMin = 1.0e-8, + dtMax = 0.1, + dtGrowth = 2.0, + steps = None, + maxSteps = None, + statsStep = 10, + domainIndependent = False, + rigorousBoundaries = False, + dtverbose = False, + + # output control + vizCycle = None, + vizTime = 1.0, + clearDirectories = False, + restoreCycle = None, + restartStep = 10000, + redistributeStep = 100000000000000, + dataDir = "dumps-DEM-impactingSquares-3d", + + # ats + checkRestart = False, + checkConservation = False, # turn on error checking for momentum conservation + conservationErrorThreshold = 2e-14, # relative error for momentum conservation + ) + +#------------------------------------------------------------------------------- +# file things +#------------------------------------------------------------------------------- +testName = "DEM-ImpactingSquares-3d" +restartDir = os.path.join(dataDir, "restarts") +vizDir = os.path.join(dataDir, "visit") +restartBaseName = os.path.join(restartDir, testName) +vizBaseName = testName + +if vizCycle is None and vizTime is None: + vizBaseName=None + +#------------------------------------------------------------------------------- +# Check if the necessary output directories exist. If not, create them. +#------------------------------------------------------------------------------- +if mpi.rank == 0: + if clearDirectories and os.path.exists(dataDir): + shutil.rmtree(dataDir) + if not os.path.exists(restartDir): + os.makedirs(restartDir) + if not os.path.exists(vizDir): + os.makedirs(vizDir) +mpi.barrier() + +#------------------------------------------------------------------------------- +# If we're restarting, find the set of most recent restart files. +#------------------------------------------------------------------------------- +if restoreCycle is None: + restoreCycle = findLastRestart(restartBaseName) + + +#------------------------------------------------------------------------------- +# This doesn't really matter kernel filler for neighbor algo +#------------------------------------------------------------------------------- +WT = TableKernel(WendlandC2Kernel(), 1000) + +#------------------------------------------------------------------------------- +# Make the NodeList. +#------------------------------------------------------------------------------- +units = CGuS() +nodes1 = makeDEMNodeList("nodeList1", + neighborSearchBuffer = neighborSearchBuffer, + kernelExtent = WT.kernelExtent) +nodeSet = [nodes1] +for nodes in nodeSet: + output("nodes.name") + output("nodes.hmin") + output("nodes.hmax") + output("nodes.hminratio") + output("nodes.nodesPerSmoothingScale") + +#------------------------------------------------------------------------------- +# Set the node properties. +#------------------------------------------------------------------------------- +if restoreCycle is None: + generator0 = GenerateNodeDistribution3d(2, 2, 1, + rho = 1.0, + distributionType = "lattice", + xmin = (0, 0.0, 1.0), + xmax = ( 1.0,1.0, 2.0), + nNodePerh = nPerh) + + # # really simple bar shaped particle + def DEMParticleGenerator(xi,yi,zi,Hi,mi,Ri): + xout = [xi] + yout = [yi] + zout = [particleRadius] + mout = [particleDensity * 4.0/3.0*np.pi*particleRadius**3] + Rout = [particleRadius] + return xout,yout,zout,mout,Rout + + generator1 = GenerateDEMfromSPHGenerator3d(WT, + generator0, + DEMParticleGenerator=DEMParticleGenerator) + + distributeNodes3d((nodes1, generator1)) + +#------------------------------------------------------------------------------- +# Construct a DataBase to hold our node list +#------------------------------------------------------------------------------- +db = DataBase() +output("db") +for nodes in nodeSet: + db.appendNodeList(nodes) +output("db.numNodeLists") +output("db.numDEMNodeLists") +output("db.numFluidNodeLists") + + +#------------------------------------------------------------------------------- +# PhysicsPackage : DEM +#------------------------------------------------------------------------------- +dem = DEM(db, + normalSpringConstant = normalSpringConstant, + normalRestitutionCoefficient = normalRestitutionCoefficient, + tangentialSpringConstant = tangentialSpringConstant, + tangentialRestitutionCoefficient = tangentialRestitutionCoefficient, + dynamicFrictionCoefficient = dynamicFriction, + staticFrictionCoefficient = staticFriction, + rollingFrictionCoefficient = rollingFriction, + torsionalFrictionCoefficient = torsionalFriction, + cohesiveTensileStrength = cohesiveTensileStrength, + shapeFactor = shapeFactor, + stepsPerCollision = stepsPerCollision) + +packages = [dem] + +solidWall = InfinitePlaneSolidBoundary(Vector(0.0, 0.0, 0.0), Vector( 0.0, 0.0, 1.0)) +#solidWall2 = CylinderSolidBoundary(Vector(0.0, 0.0, -10.0),Vector(0.0, 0.0,1.0),5.0,15.0) + +dem.appendSolidBoundary(solidWall) +#dem.appendSolidBoundary(solidWall2) +# #------------------------------------------------------------------------------- +# # PhysicsPackage : gravity +# #------------------------------------------------------------------------------- +gravity = ConstantAcceleration(a0 = Vector(0.0, 0.0, -1.00), + nodeList = nodes1) +packages += [gravity] + +# #------------------------------------------------------------------------------- +# # Create boundary conditions. +# #------------------------------------------------------------------------------- +# plane1 = Plane(Vector(0.0, 0.0, 0.0), Vector( 0.0, 0.0, 1.0)) +# bc1 = ReflectingBoundary(plane1) +# bcSet = [bc1] + +# for p in packages: +# for bc in bcSet: +# p.appendBoundary(bc) + +#------------------------------------------------------------------------------- +# Fields and Variables +#------------------------------------------------------------------------------- +numNodeLists = db.numNodeLists +nodeLists = db.nodeLists() + +position = db.DEMPosition +mass = db.DEMMass +velocity = db.DEMVelocity +H = db.DEMHfield +radius = db.DEMParticleRadius +compositeParticleIndex = db.DEMCompositeParticleIndex + +uniqueIndex = db.DEMUniqueIndex +omega = dem.omega + +# pure rolling +if mpi.rank == 0 : + for i in range(nodes.numInternalNodes): + if i > 0: + velocity[0][i]=Vector(particleVelocity,0,0) + omega[0][i]=Vector(0,particleVelocity/particleRadius,0) + +#------------------------------------------------------------------------------- +# Initial Conditions +#------------------------------------------------------------------------------- +#------------------------------------------------------------------------------- +# Construct a time integrator, and add the physics packages. +#------------------------------------------------------------------------------- + +integrator = IntegratorConstructor(db) +for p in packages: + integrator.appendPhysicsPackage(p) +integrator.lastDt = dt +integrator.dtMin = dtMin +integrator.dtMax = dtMax +integrator.dtGrowth = dtGrowth +integrator.domainDecompositionIndependent = domainIndependent +integrator.verbose = dtverbose +integrator.rigorousBoundaries = rigorousBoundaries + +integrator.cullGhostNodes = False + +output("integrator") +output("integrator.havePhysicsPackage(dem)") +output("integrator.lastDt") +output("integrator.dtMin") +output("integrator.dtMax") +output("integrator.dtGrowth") +output("integrator.domainDecompositionIndependent") +output("integrator.rigorousBoundaries") +output("integrator.verbose") + +#------------------------------------------------------------------------------- +# Periodic Work Function: Track conseravation +#------------------------------------------------------------------------------- + +conservation = TrackConservation(db, + dem, + verbose=True) + +periodicWork = [(conservation.periodicWorkFunction,1)] + + +#------------------------------------------------------------------------------- +# Make the problem controller. +#------------------------------------------------------------------------------- +from SpheralPointmeshSiloDump import dumpPhysicsState +control = SpheralController(integrator, + iterateInitialH = False, + initializeDerivatives = True, + statsStep = statsStep, + restartStep = restartStep, + redistributeStep=redistributeStep, + restartBaseName = restartBaseName, + restoreCycle = restoreCycle, + vizBaseName = vizBaseName, + vizMethod = dumpPhysicsState, + vizGhosts=True, + vizDir = vizDir, + vizStep = vizCycle, + vizTime = vizTime, + periodicWork=periodicWork) +output("control") + +#control.redistribute = PeanoHilbertOrderRedistributeNodes(db.maxKernelExtent,workBalance=False) +#------------------------------------------------------------------------------- +# Advance to the end time. +#------------------------------------------------------------------------------- + +if not steps is None: + control.step(steps) +else: + control.advance(goalTime, maxSteps) + diff --git a/tests/functional/DEM/LinearSpringDEM/SolidBoundaryCondition/torsionalParticleTest.py b/tests/functional/DEM/LinearSpringDEM/SolidBoundaryCondition/torsionalParticleTest.py new file mode 100644 index 000000000..1704c98f1 --- /dev/null +++ b/tests/functional/DEM/LinearSpringDEM/SolidBoundaryCondition/torsionalParticleTest.py @@ -0,0 +1,313 @@ +#ATS:DEM3dImpact = test( SELF, "--clearDirectories True --checkConservation True --goalTime 1.0", label="DEM impacting squares -- 3-D (parallel)", np=8) + +import os, sys, shutil, mpi, random +from math import * +from Spheral3d import * +from SpheralTestUtilities import * +from findLastRestart import * +from GenerateNodeDistribution3d import * +from GenerateDEMfromSPHGenerator import GenerateDEMfromSPHGenerator3d + +import numpy as np + +sys.path.insert(0, '..') +from DEMConservationTracker import TrackConservation3d as TrackConservation + +if mpi.procs > 1: + from PeanoHilbertDistributeNodes import distributeNodes3d +else: + from DistributeNodes import distributeNodes3d + +title("DEM 3d Drop Test with Particle Generation") +# this tests the ability to generate particle on the fly +# during the course of a simulation using a periodic +# work function. It also tests the solid boundary condition + +#------------------------------------------------------------------------------- +# Generic problem parameters +#------------------------------------------------------------------------------- +commandLine(numParticlePerLength = 3, # number of particles on a side of the box + normalSpringConstant=100.0, # spring constant for LDS model + normalRestitutionCoefficient=0.55, # restitution coefficient to get damping const + tangentialSpringConstant=20.857, # spring constant for LDS model + tangentialRestitutionCoefficient=1.00, # restitution coefficient to get damping const + dynamicFriction = 1.0, # static friction coefficient sliding + staticFriction = 1.0, # dynamic friction coefficient sliding + rollingFriction = 1.05, # static friction coefficient for rolling + torsionalFriction = 1.3, # static friction coefficient for torsion + cohesiveTensileStrength = 0.0, # units of pressure + shapeFactor = 0.1, # in [0,1] shape factor from Zhang 2018, 0 - no torsion or rolling + + particleRadius = 1.0, # particle radius + particleDensity = 2.60, + particleVelocity = 0.1, + + neighborSearchBuffer = 0.1, # multiplicative buffer to radius for neighbor search algo + nPerh = 1.01, + + + useSolidBoundary = True, + + # integration + IntegratorConstructor = VerletIntegrator, + stepsPerCollision = 100, # replaces CFL for DEM + goalTime = 25.0, + dt = 1e-8, + dtMin = 1.0e-8, + dtMax = 0.1, + dtGrowth = 2.0, + steps = None, + maxSteps = None, + statsStep = 10, + domainIndependent = False, + rigorousBoundaries = False, + dtverbose = False, + + # output control + vizCycle = None, + vizTime = 1.0, + clearDirectories = False, + restoreCycle = None, + restartStep = 10000, + redistributeStep = 100000000000000, + dataDir = "dumps-DEM-impactingSquares-3d", + + # ats + checkRestart = False, + checkConservation = False, # turn on error checking for momentum conservation + conservationErrorThreshold = 2e-14, # relative error for momentum conservation + ) + +#------------------------------------------------------------------------------- +# file things +#------------------------------------------------------------------------------- +testName = "DEM-ImpactingSquares-3d" +restartDir = os.path.join(dataDir, "restarts") +vizDir = os.path.join(dataDir, "visit") +restartBaseName = os.path.join(restartDir, testName) +vizBaseName = testName + +if vizCycle is None and vizTime is None: + vizBaseName=None + +#------------------------------------------------------------------------------- +# Check if the necessary output directories exist. If not, create them. +#------------------------------------------------------------------------------- +if mpi.rank == 0: + if clearDirectories and os.path.exists(dataDir): + shutil.rmtree(dataDir) + if not os.path.exists(restartDir): + os.makedirs(restartDir) + if not os.path.exists(vizDir): + os.makedirs(vizDir) +mpi.barrier() + +#------------------------------------------------------------------------------- +# If we're restarting, find the set of most recent restart files. +#------------------------------------------------------------------------------- +if restoreCycle is None: + restoreCycle = findLastRestart(restartBaseName) + + +#------------------------------------------------------------------------------- +# This doesn't really matter kernel filler for neighbor algo +#------------------------------------------------------------------------------- +WT = TableKernel(WendlandC2Kernel(), 1000) + +#------------------------------------------------------------------------------- +# Make the NodeList. +#------------------------------------------------------------------------------- +units = CGuS() +nodes1 = makeDEMNodeList("nodeList1", + neighborSearchBuffer = neighborSearchBuffer, + kernelExtent = WT.kernelExtent) +nodeSet = [nodes1] +for nodes in nodeSet: + output("nodes.name") + output("nodes.hmin") + output("nodes.hmax") + output("nodes.hminratio") + output("nodes.nodesPerSmoothingScale") + +#------------------------------------------------------------------------------- +# Set the node properties. +#------------------------------------------------------------------------------- +if restoreCycle is None: + generator0 = GenerateNodeDistribution3d(2, 2, 1, + rho = 1.0, + distributionType = "lattice", + xmin = (0, 0.0, 0), + xmax = ( particleRadius*5,particleRadius*5, particleRadius*5), + nNodePerh = nPerh) + + # # really simple bar shaped particle + def DEMParticleGenerator(xi,yi,zi,Hi,mi,Ri): + xout = [xi] + yout = [yi] + zout = [particleRadius] + mout = [particleDensity * 4.0/3.0*np.pi*particleRadius**3] + Rout = [particleRadius] + return xout,yout,zout,mout,Rout + + generator1 = GenerateDEMfromSPHGenerator3d(WT, + generator0, + DEMParticleGenerator=DEMParticleGenerator) + + distributeNodes3d((nodes1, generator1)) + +#------------------------------------------------------------------------------- +# Construct a DataBase to hold our node list +#------------------------------------------------------------------------------- +db = DataBase() +output("db") +for nodes in nodeSet: + db.appendNodeList(nodes) +output("db.numNodeLists") +output("db.numDEMNodeLists") +output("db.numFluidNodeLists") + + +#------------------------------------------------------------------------------- +# PhysicsPackage : DEM +#------------------------------------------------------------------------------- +dem = DEM(db, + normalSpringConstant = normalSpringConstant, + normalRestitutionCoefficient = normalRestitutionCoefficient, + tangentialSpringConstant = tangentialSpringConstant, + tangentialRestitutionCoefficient = tangentialRestitutionCoefficient, + dynamicFrictionCoefficient = dynamicFriction, + staticFrictionCoefficient = staticFriction, + rollingFrictionCoefficient = rollingFriction, + torsionalFrictionCoefficient = torsionalFriction, + cohesiveTensileStrength = cohesiveTensileStrength, + shapeFactor = shapeFactor, + stepsPerCollision = stepsPerCollision) + + +packages = [dem] + + +# #------------------------------------------------------------------------------- +# # PhysicsPackage : gravity +# #------------------------------------------------------------------------------- +gravity = ConstantAcceleration(a0 = Vector(0.0, 0.0, -1.00), + nodeList = nodes1) +#packages += [gravity] + +if useSolidBoundary: + solidWall = InfinitePlaneSolidBoundary(Vector(0.0, 0.0, 0.0), Vector( 0.0, 0.0, 1.0)) + dem.appendSolidBoundary(solidWall) +else: + bcs = [ReflectingBoundary(Plane(Vector(0.0, 0.0, 0.0), Vector( 0.0, 0.0, 1.0)))] + for package in packages: + for bc in bcs: + package.appendBoundary(bc) +# #------------------------------------------------------------------------------- +# # Create boundary conditions. +# #------------------------------------------------------------------------------- +# plane1 = Plane(Vector(0.0, 0.0, 0.0), Vector( 0.0, 0.0, 1.0)) +# bc1 = ReflectingBoundary(plane1) +# bcSet = [bc1] + +# for p in packages: +# for bc in bcSet: +# p.appendBoundary(bc) + +#------------------------------------------------------------------------------- +# Fields and Variables +#------------------------------------------------------------------------------- +numNodeLists = db.numNodeLists +nodeLists = db.nodeLists() + +position = db.DEMPosition +mass = db.DEMMass +velocity = db.DEMVelocity +H = db.DEMHfield +radius = db.DEMParticleRadius +compositeParticleIndex = db.DEMCompositeParticleIndex + +uniqueIndex = db.DEMUniqueIndex +omega = dem.omega + + +#------------------------------------------------------------------------------- +# Initial Conditions +#------------------------------------------------------------------------------- +#------------------------------------------------------------------------------- +# Construct a time integrator, and add the physics packages. +#------------------------------------------------------------------------------- + +integrator = IntegratorConstructor(db) +for p in packages: + integrator.appendPhysicsPackage(p) +integrator.lastDt = dt +integrator.dtMin = dtMin +integrator.dtMax = dtMax +integrator.dtGrowth = dtGrowth +integrator.domainDecompositionIndependent = domainIndependent +integrator.verbose = dtverbose +integrator.rigorousBoundaries = rigorousBoundaries + +integrator.cullGhostNodes = False + +output("integrator") +output("integrator.havePhysicsPackage(dem)") +output("integrator.lastDt") +output("integrator.dtMin") +output("integrator.dtMax") +output("integrator.dtGrowth") +output("integrator.domainDecompositionIndependent") +output("integrator.rigorousBoundaries") +output("integrator.verbose") + +#------------------------------------------------------------------------------- +# Periodic Work Function: Track conseravation +#------------------------------------------------------------------------------- + +conservation = TrackConservation(db, + dem, + verbose=True) + +periodicWork = [(conservation.periodicWorkFunction,1)] + +for i in range(nodes.numInternalNodes): + if i > 0: + velocity[0][i]=Vector(0,0,-0.1) + #omega[0][i]=Vector(0,0,particleVelocity/particleRadius) + +#------------------------------------------------------------------------------- +# Make the problem controller. +#------------------------------------------------------------------------------- +from SpheralPointmeshSiloDump import dumpPhysicsState +control = SpheralController(integrator, + iterateInitialH = False, + initializeDerivatives = True, + statsStep = statsStep, + restartStep = restartStep, + redistributeStep=redistributeStep, + restartBaseName = restartBaseName, + restoreCycle = restoreCycle, + vizBaseName = vizBaseName, + vizMethod = dumpPhysicsState, + vizGhosts=True, + vizDir = vizDir, + vizStep = vizCycle, + vizTime = vizTime, + periodicWork=periodicWork) +output("control") + +#control.redistribute = PeanoHilbertOrderRedistributeNodes(db.maxKernelExtent,workBalance=False) +#------------------------------------------------------------------------------- +# Advance to the end time. +#------------------------------------------------------------------------------- + +if not steps is None: + control.step(steps) +else: + control.advance(10.0, maxSteps) + # for i in range(nodes.numInternalNodes): + # if i > 0: + # velocity[0][i]=Vector(0,0,0) + # omega[0][i]=Vector(0,0,particleVelocity/particleRadius) + control.advance(10.0+goalTime, maxSteps) + From 717fdd943c3d25796e1fc6688b63991ae86eeef7 Mon Sep 17 00:00:00 2001 From: jmpearl Date: Sun, 3 Mar 2024 21:18:19 -0800 Subject: [PATCH 17/38] DEM 3d impact squares - making it reproducable --- .../DEM/LinearSpringDEM/ImpactingSquares/impactingSquares-3d.py | 2 ++ 1 file changed, 2 insertions(+) diff --git a/tests/functional/DEM/LinearSpringDEM/ImpactingSquares/impactingSquares-3d.py b/tests/functional/DEM/LinearSpringDEM/ImpactingSquares/impactingSquares-3d.py index 732f513b8..49640ab04 100644 --- a/tests/functional/DEM/LinearSpringDEM/ImpactingSquares/impactingSquares-3d.py +++ b/tests/functional/DEM/LinearSpringDEM/ImpactingSquares/impactingSquares-3d.py @@ -16,6 +16,8 @@ else: from DistributeNodes import distributeNodes3d +random.seed(0) + title("DEM 3d Impacting Squares") # This tests the conservation properties of the DEM package when # distribution across multiple processors From 423d72ae500cc61aa12cf1f36733836d9f43efd6 Mon Sep 17 00:00:00 2001 From: jmpearl Date: Sun, 3 Mar 2024 21:35:08 -0800 Subject: [PATCH 18/38] DEM adding Scalars to allowed state fields --- src/DataBase/StateBase.cc | 21 ++++++++++++++++++--- 1 file changed, 18 insertions(+), 3 deletions(-) diff --git a/src/DataBase/StateBase.cc b/src/DataBase/StateBase.cc index cae66eab2..4ae5e7c99 100644 --- a/src/DataBase/StateBase.cc +++ b/src/DataBase/StateBase.cc @@ -148,7 +148,16 @@ operator==(const StateBase& rhs) const { result = false; } } catch (const boost::bad_any_cast&) { - std::cerr << "StateBase::operator== WARNING: unable to compare values for " << lhsItr->first << "\n"; + try { + auto lhsPtr = boost::any_cast(lhsItr->second); + auto rhsPtr = boost::any_cast(rhsItr->second); + if (*lhsPtr != *rhsPtr) { + cerr << "Scalar for " << lhsItr->first << " don't match." << endl; + result = false; + } + } catch (const boost::bad_any_cast&) { + std::cerr << "StateBase::operator== WARNING: unable to compare values for " << lhsItr->first << "\n"; + } } } } @@ -405,8 +414,14 @@ assign(const StateBase& rhs) { const auto rhsptr = boost::any_cast(anyrhs); *lhsptr = *rhsptr; } catch(const boost::bad_any_cast&) { - // We'll assume other things don't need to be assigned... - // VERIFY2(false, "StateBase::assign ERROR: unknown type for key " << itr->first << "\n"); + try { + auto lhsptr = boost::any_cast(anylhs); + const auto rhsptr = boost::any_cast(anyrhs); + *lhsptr = *rhsptr; + } catch(const boost::bad_any_cast&) { + // We'll assume other things don't need to be assigned... + // VERIFY2(false, "StateBase::assign ERROR: unknown type for key " << itr->first << "\n"); + } } } } From 7a4d59a71b1fbfb099c0fe48da277e8d99ee9adb Mon Sep 17 00:00:00 2001 From: jmpearl Date: Mon, 4 Mar 2024 20:22:59 -0800 Subject: [PATCH 19/38] DEM misplaced semicolon caused debug build failure --- src/DEM/SolidBoundary/SolidBoundaryBaseInline.hh | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/DEM/SolidBoundary/SolidBoundaryBaseInline.hh b/src/DEM/SolidBoundary/SolidBoundaryBaseInline.hh index edd8dc21e..a784b94ac 100644 --- a/src/DEM/SolidBoundary/SolidBoundaryBaseInline.hh +++ b/src/DEM/SolidBoundary/SolidBoundaryBaseInline.hh @@ -15,6 +15,6 @@ int SolidBoundaryBase:: uniqueIndex() const { return mUniqueIndex; -}; +} } From 821cc1526ef2ee04d43721b4d2e1383e2b3a42ad Mon Sep 17 00:00:00 2001 From: jmpearl Date: Thu, 7 Mar 2024 06:18:19 -0800 Subject: [PATCH 20/38] DEMs max overlap needed to be restartable" --- src/DEM/LinearSpringDEM.cc | 8 ++++++++ 1 file changed, 8 insertions(+) diff --git a/src/DEM/LinearSpringDEM.cc b/src/DEM/LinearSpringDEM.cc index 479e201d8..b98b5c61b 100644 --- a/src/DEM/LinearSpringDEM.cc +++ b/src/DEM/LinearSpringDEM.cc @@ -259,6 +259,10 @@ dt(const DataBase& dataBase, // Ensure no point moves further than the buffer distance in one timestep + //-------------------------------------------------------------------------------------- + // NOTE: it would be nice if this wasn't based on the absolute velocity for cases + // where we have a blob of dem particles moving at elevated speeds + //-------------------------------------------------------------------------------------- const auto numNodeLists = position.size(); for (auto k = 0u; k < numNodeLists; ++k) { const auto n = position[k]->size(); @@ -841,6 +845,8 @@ LinearSpringDEM:: dumpState(FileIO& file, const string& pathName) const { DEMBase::dumpState(file,pathName); file.write(mMomentOfInertia, pathName + "/momentOfInertia"); + file.write(mMaximumOverlap, pathName + "/maximumOverlap"); + file.write(mNewMaximumOverlap, pathName + "/newMaximumOverlap"); } //------------------------------------------------------------------------------ @@ -852,6 +858,8 @@ LinearSpringDEM:: restoreState(const FileIO& file, const string& pathName) { DEMBase::restoreState(file,pathName); file.read(mMomentOfInertia, pathName + "/momentOfInertia"); + file.read(mMaximumOverlap, pathName + "/maximumOverlap"); + file.read(mNewMaximumOverlap, pathName + "/newMaximumOverlap"); } } // namespace From c4eb76d8d61725b9ea56dcc0f0bc77bbdc239a8a Mon Sep 17 00:00:00 2001 From: jmpearl Date: Mon, 18 Mar 2024 12:24:51 -0700 Subject: [PATCH 21/38] making fast time stepping optional to allows integrator's updateBCFreq option --- src/DEM/DEM.py | 5 ++++ src/DEM/LinearSpringDEM.cc | 51 +++++++++++++++++++++++++++----- src/DEM/LinearSpringDEM.hh | 18 +++++++++++ src/DEM/LinearSpringDEMInline.hh | 34 +++++++++++++++++++++ src/PYB11/DEM/LinearSpringDEM.py | 5 +++- 5 files changed, 104 insertions(+), 9 deletions(-) diff --git a/src/DEM/DEM.py b/src/DEM/DEM.py index f24ae5f9b..fd70c08e9 100644 --- a/src/DEM/DEM.py +++ b/src/DEM/DEM.py @@ -18,6 +18,7 @@ def LinearSpringDEM(dataBase, cohesiveTensileStrength = 0.0, shapeFactor = 0.0, stepsPerCollision = 25, + enableFastTimeStepping = False, xmin = (-1e100, -1e100, -1e100), xmax = ( 1e100, 1e100, 1e100)): @@ -32,6 +33,7 @@ def LinearSpringDEM(dataBase, assert dynamicFrictionCoefficient >= 0, "dynamicFrictionCoefficient must be positive" assert rollingFrictionCoefficient >= 0, "rollingFrictionCoefficient must be positive" assert torsionalFrictionCoefficient >= 0, "torsionalFrictionCoefficient must be positive" + assert isinstance(enableFastTimeStepping,bool) #if (stepsPerCollision < 10) print("WARNING: stepsPerCollision is very low, reccomended is 25-50") @@ -60,6 +62,7 @@ def LinearSpringDEM(dataBase, "cohesiveTensileStrength" : cohesiveTensileStrength, "shapeFactor" : shapeFactor, "stepsPerCollision" : stepsPerCollision, + "enableFastTimeStepping" : enableFastTimeStepping, "xmin" : eval("Vector%id(%g, %g, %g)" % xmin), "xmax" : eval("Vector%id(%g, %g, %g)" % xmax)} @@ -83,6 +86,7 @@ def DEM(dataBase, cohesiveTensileStrength=0.0, shapeFactor=0.0, stepsPerCollision = 25, + enableFastTimeStepping = False, xmin = (-1e100, -1e100, -1e100), xmax = ( 1e100, 1e100, 1e100)): return LinearSpringDEM(dataBase, @@ -97,6 +101,7 @@ def DEM(dataBase, cohesiveTensileStrength, shapeFactor, stepsPerCollision, + enableFastTimeStepping, xmin, xmax) diff --git a/src/DEM/LinearSpringDEM.cc b/src/DEM/LinearSpringDEM.cc index b98b5c61b..bf54fabb2 100644 --- a/src/DEM/LinearSpringDEM.cc +++ b/src/DEM/LinearSpringDEM.cc @@ -80,9 +80,11 @@ LinearSpringDEM(const DataBase& dataBase, const Scalar cohesiveTensileStrength, const Scalar shapeFactor, const Scalar stepsPerCollision, + const bool enableFastTimeStepping, const Vector& xmin, const Vector& xmax): DEMBase(dataBase,stepsPerCollision,xmin,xmax), + mEnableFastTimeStepping(enableFastTimeStepping), mNormalSpringConstant(normalSpringConstant), mNormalRestitutionCoefficient(normalRestitutionCoefficient), mTangentialSpringConstant(tangentialSpringConstant), @@ -95,12 +97,17 @@ LinearSpringDEM(const DataBase& dataBase, mShapeFactor(shapeFactor), mNormalBeta(M_PI/std::log(std::max(normalRestitutionCoefficient,1.0e-3))), mTangentialBeta(M_PI/std::log(std::max(tangentialRestitutionCoefficient,1.0e-3))), + mCollisionDuration(0.0), mMomentOfInertia(FieldStorageType::CopyFields), mMaximumOverlap(FieldStorageType::CopyFields), mNewMaximumOverlap(FieldStorageType::CopyFields) { mMomentOfInertia = dataBase.newDEMFieldList(0.0, DEMFieldNames::momentOfInertia); mMaximumOverlap = dataBase.newDEMFieldList(0.0, DEMFieldNames::maximumOverlap); mNewMaximumOverlap = dataBase.newDEMFieldList(0.0,MaxReplaceState::prefix() + DEMFieldNames::maximumOverlap); + + const auto mass = dataBase.DEMMass(); + const auto minMass = mass.min(); + mCollisionDuration = M_PI*std::sqrt(0.5*minMass/mNormalSpringConstant * (1.0 + 1.0/(mNormalBeta*mNormalBeta))); } //------------------------------------------------------------------------------ @@ -118,20 +125,50 @@ typename LinearSpringDEM::TimeStepType LinearSpringDEM:: dt(const DataBase& dataBase, const State& state, - const StateDerivatives& /*derivs*/, - const typename Dimension::Scalar /*currentTime*/) const{ + const StateDerivatives& derivs, + const typename Dimension::Scalar currentTime) const { TIME_BEGIN("LinearSpringDEMdt"); + auto dtMin = std::numeric_limits::max(); + TimeStepType result(dtMin, "DEM error, this message should not get to the end"); + + if (this->enableFastTimeStepping()){ + result = this->variableTimeStep(dataBase, + state, + derivs, + currentTime); + }else{ + result = this->fixedTimeStep(); + } + + TIME_END("LinearSpringDEMdt"); + return result; +} + +//------------------------------------------------------------------------------ +// set our timestep (for now its a constant single value) +//------------------------------------------------------------------------------ +template +typename LinearSpringDEM::TimeStepType +LinearSpringDEM:: +fixedTimeStep() const { + return make_pair(mCollisionDuration/this->stepsPerCollision(),("fixed-dt Linear Spring DEM vote for time step")); +} + +template +typename LinearSpringDEM::TimeStepType +LinearSpringDEM:: +variableTimeStep(const DataBase& dataBase, + const State& state, + const StateDerivatives& /*derivs*/, + const typename Dimension::Scalar /*currentTime*/) const { + // Get some useful fluid variables from the DataBase. const auto mass = state.fields(HydroFieldNames::mass, 0.0); const auto position = state.fields(HydroFieldNames::position, Vector::zero); const auto velocity = state.fields(HydroFieldNames::velocity, Vector::zero); const auto r = state.fields(DEMFieldNames::particleRadius, 0.0); - //const auto& connectivityMap = dataBase.connectivityMap(this->requireGhostConnectivity(), - // this->requireOverlapConnectivity(), - // this->requireIntersectionConnectivity()); - const auto& contacts = this->contactStorageIndices(); const unsigned int numP2PContacts = this->numParticleParticleContacts(); @@ -280,11 +317,9 @@ dt(const DataBase& dataBase, } } - TIME_END("LinearSpringDEMdt"); return result; } - //------------------------------------------------------------------------------ // method that fires once on startup //------------------------------------------------------------------------------ diff --git a/src/DEM/LinearSpringDEM.hh b/src/DEM/LinearSpringDEM.hh index 290b8fa5d..6ac98bbcd 100644 --- a/src/DEM/LinearSpringDEM.hh +++ b/src/DEM/LinearSpringDEM.hh @@ -49,6 +49,7 @@ public: const Scalar cohesiveTensileStrength, const Scalar shapeFactor, const Scalar stepsPerCollision, + const bool enableFastTimeStepping, const Vector& xmin, const Vector& xmax); @@ -80,6 +81,15 @@ public: virtual void enforceBoundaries(State& state, StateDerivatives& derivs) override; + // sub-methods for dt + TimeStepType variableTimeStep(const DataBase& dataBase, + const State& state, + const StateDerivatives& derivs, + const Scalar time) const; + + TimeStepType fixedTimeStep() const; + + // generalized spring damper functions (inlined) void slidingSpringDamper(const Scalar k, const Scalar C, const Scalar mus, @@ -106,6 +116,9 @@ public: Scalar& force) const; // set/gets + bool enableFastTimeStepping() const; + void enableFastTimeStepping(bool x); + Scalar normalSpringConstant() const; void normalSpringConstant(Scalar x); @@ -142,6 +155,9 @@ public: Scalar tangentialBeta() const; void tangentialBeta(Scalar x); + Scalar collisionDuration() const; + void collisionDuration(Scalar x); + // set moment of inertia on start up void setMomentOfInertia(); @@ -163,6 +179,7 @@ public: //**************************************************************************** private: //--------------------------- Private Interface ---------------------------// + Scalar mEnableFastTimeStepping; Scalar mNormalSpringConstant; Scalar mNormalRestitutionCoefficient; Scalar mTangentialSpringConstant; @@ -176,6 +193,7 @@ private: Scalar mNormalBeta; Scalar mTangentialBeta; + Scalar mCollisionDuration; // field Lists FieldList mMomentOfInertia; diff --git a/src/DEM/LinearSpringDEMInline.hh b/src/DEM/LinearSpringDEMInline.hh index 6668bfe99..90bb53884 100644 --- a/src/DEM/LinearSpringDEMInline.hh +++ b/src/DEM/LinearSpringDEMInline.hh @@ -1,5 +1,25 @@ namespace Spheral { + +//------------------------------------------------------------------------------ +// set/get to activate/deactivate fast timestepping +//------------------------------------------------------------------------------ +template +inline +bool +LinearSpringDEM:: +enableFastTimeStepping() const { + return mEnableFastTimeStepping; +} +template +inline +void +LinearSpringDEM:: +enableFastTimeStepping(bool x) { + mEnableFastTimeStepping = x; +} + + //------------------------------------------------------------------------------ // set/get our spring constant //------------------------------------------------------------------------------ @@ -217,6 +237,20 @@ tangentialBeta(typename Dimension::Scalar x) { mTangentialBeta = x; } +template +inline +typename Dimension::Scalar +LinearSpringDEM:: +collisionDuration() const { + return mCollisionDuration; +} +template +inline +void +LinearSpringDEM:: +collisionDuration(typename Dimension::Scalar x) { + mCollisionDuration = x; +} //------------------------------------------------------------------------------ // moment of interia specializations diff --git a/src/PYB11/DEM/LinearSpringDEM.py b/src/PYB11/DEM/LinearSpringDEM.py index 7500c626c..1eb288cb6 100644 --- a/src/PYB11/DEM/LinearSpringDEM.py +++ b/src/PYB11/DEM/LinearSpringDEM.py @@ -26,6 +26,7 @@ def pyinit(dataBase = "const DataBase<%(Dimension)s>&", cohesiveTensileStrength = "const Scalar", shapeFactor = "const Scalar", stepsPerCollision = "const Scalar", + enableFastTimeStepping = "const bool", xmin = "const Vector&", xmax = "const Vector&"): "DEMBase constructor" @@ -68,7 +69,8 @@ def momentOfInertia(massi = "const Scalar", def setMomentOfInertia(self): return "void" - + + enableFastTimeStepping = PYB11property("bool", "enableFastTimeStepping", "enableFastTimeStepping", doc="activate fast time stepping") normalSpringConstant = PYB11property("Scalar", "normalSpringConstant", "normalSpringConstant", doc="normal spring constant") normalRestitutionCoefficient = PYB11property("Scalar", "normalRestitutionCoefficient", "normalRestitutionCoefficient", doc="normal restitution coefficient") tangentialSpringConstant = PYB11property("Scalar", "tangentialSpringConstant", "tangentialSpringConstant", doc="tangential spring constant") @@ -83,5 +85,6 @@ def setMomentOfInertia(self): shapeFactor = PYB11property("Scalar", "shapeFactor", "shapeFactor", doc="shape factor - simple approach to non-spherical particles") normalBeta = PYB11property("Scalar", "normalBeta", "normalBeta", doc="a damping parameter") tangentialBeta = PYB11property("Scalar", "tangentialBeta", "tangentialBeta", doc="a damping parameter") + collisionDuration = PYB11property("Scalar", "collisionDuration", "collisionDuration", doc="duration of a contact") momentOfInertia = PYB11property("const FieldList<%(Dimension)s, Scalar>&","momentOfInertia", returnpolicy="reference_internal") \ No newline at end of file From c7094c54a7f247062c2bbfe5ad53a2d6da5bbb29 Mon Sep 17 00:00:00 2001 From: jmpearl Date: Mon, 18 Mar 2024 15:02:19 -0700 Subject: [PATCH 22/38] some DEM tests assume fast time stepping is on --- src/DEM/DEM.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/DEM/DEM.py b/src/DEM/DEM.py index fd70c08e9..98ba8ac91 100644 --- a/src/DEM/DEM.py +++ b/src/DEM/DEM.py @@ -18,7 +18,7 @@ def LinearSpringDEM(dataBase, cohesiveTensileStrength = 0.0, shapeFactor = 0.0, stepsPerCollision = 25, - enableFastTimeStepping = False, + enableFastTimeStepping = True, xmin = (-1e100, -1e100, -1e100), xmax = ( 1e100, 1e100, 1e100)): @@ -86,7 +86,7 @@ def DEM(dataBase, cohesiveTensileStrength=0.0, shapeFactor=0.0, stepsPerCollision = 25, - enableFastTimeStepping = False, + enableFastTimeStepping = True, xmin = (-1e100, -1e100, -1e100), xmax = ( 1e100, 1e100, 1e100)): return LinearSpringDEM(dataBase, From 525f9e4ec9f20ddd05e826e4ed09211e02d68e06 Mon Sep 17 00:00:00 2001 From: jmpearl Date: Tue, 19 Mar 2024 10:32:33 -0700 Subject: [PATCH 23/38] DEM update ATS and tangentialDampTerm Coeffs --- src/DEM/LinearSpringDEM.cc | 2 +- tests/functional/DEM/DEMTests.ats | 1 + ...onalParticleTest.py => pureTorision-3d.py} | 155 +++++++++++------- 3 files changed, 102 insertions(+), 56 deletions(-) rename tests/functional/DEM/LinearSpringDEM/SolidBoundaryCondition/{torsionalParticleTest.py => pureTorision-3d.py} (70%) diff --git a/src/DEM/LinearSpringDEM.cc b/src/DEM/LinearSpringDEM.cc index bf54fabb2..d59adebe0 100644 --- a/src/DEM/LinearSpringDEM.cc +++ b/src/DEM/LinearSpringDEM.cc @@ -412,7 +412,7 @@ evaluateDerivatives(const typename Dimension::Scalar /*time*/, const auto invKr = 1.0/max(kr,tiny); const auto normalDampingTerms = 2.0*kn/(1.0+mNormalBeta*mNormalBeta); - const auto tangentialDampingTerms = 4.0/5.0*ks/(1.0+mTangentialBeta*mTangentialBeta); + const auto tangentialDampingTerms = 2.0*ks/(1.0+mTangentialBeta*mTangentialBeta); // The connectivity. const auto& nodeLists = dataBase.DEMNodeListPtrs(); diff --git a/tests/functional/DEM/DEMTests.ats b/tests/functional/DEM/DEMTests.ats index d2204d41e..be2784348 100644 --- a/tests/functional/DEM/DEMTests.ats +++ b/tests/functional/DEM/DEMTests.ats @@ -13,3 +13,4 @@ source("LinearSpringDEM/ImpactingSquares/impactingSquares-3d.py") # test restitution coefficient and friction against solid boundaries source("LinearSpringDEM/SolidBoundaryCondition/singleParticleBoundaryCollision-3d.py") source("LinearSpringDEM/SolidBoundaryCondition/finiteSolidPlanes-3d.py") +source("LinearSpringDEM/SolidBoundaryCondition/pureTorsion-3d.py") diff --git a/tests/functional/DEM/LinearSpringDEM/SolidBoundaryCondition/torsionalParticleTest.py b/tests/functional/DEM/LinearSpringDEM/SolidBoundaryCondition/pureTorision-3d.py similarity index 70% rename from tests/functional/DEM/LinearSpringDEM/SolidBoundaryCondition/torsionalParticleTest.py rename to tests/functional/DEM/LinearSpringDEM/SolidBoundaryCondition/pureTorision-3d.py index 1704c98f1..9fffdd123 100644 --- a/tests/functional/DEM/LinearSpringDEM/SolidBoundaryCondition/torsionalParticleTest.py +++ b/tests/functional/DEM/LinearSpringDEM/SolidBoundaryCondition/pureTorision-3d.py @@ -1,4 +1,4 @@ -#ATS:DEM3dImpact = test( SELF, "--clearDirectories True --checkConservation True --goalTime 1.0", label="DEM impacting squares -- 3-D (parallel)", np=8) +#ATS:DEM3dPureTorsion = test( SELF, "--clearDirectories True --checkRestitutionCoefficient True --checkNaturalFrequency True", label="DEM puretorsion test -- 3-D (serial)") import os, sys, shutil, mpi, random from math import * @@ -18,19 +18,16 @@ else: from DistributeNodes import distributeNodes3d -title("DEM 3d Drop Test with Particle Generation") -# this tests the ability to generate particle on the fly -# during the course of a simulation using a periodic -# work function. It also tests the solid boundary condition - +title("DEM 3d Pure Torsion Test") +# This tests the natural freq and restitution coefficient for Zhang et.al. formulation #------------------------------------------------------------------------------- # Generic problem parameters #------------------------------------------------------------------------------- commandLine(numParticlePerLength = 3, # number of particles on a side of the box - normalSpringConstant=100.0, # spring constant for LDS model + normalSpringConstant=100.0, # spring constant for LDS model normalRestitutionCoefficient=0.55, # restitution coefficient to get damping const - tangentialSpringConstant=20.857, # spring constant for LDS model - tangentialRestitutionCoefficient=1.00, # restitution coefficient to get damping const + tangentialSpringConstant=20.857, # spring constant for LDS model + tangentialRestitutionCoefficient=0.55, # restitution coefficient to get damping const dynamicFriction = 1.0, # static friction coefficient sliding staticFriction = 1.0, # dynamic friction coefficient sliding rollingFriction = 1.05, # static friction coefficient for rolling @@ -38,11 +35,11 @@ cohesiveTensileStrength = 0.0, # units of pressure shapeFactor = 0.1, # in [0,1] shape factor from Zhang 2018, 0 - no torsion or rolling - particleRadius = 1.0, # particle radius + particleRadius = 1.0, # particle radius particleDensity = 2.60, particleVelocity = 0.1, - neighborSearchBuffer = 0.1, # multiplicative buffer to radius for neighbor search algo + neighborSearchBuffer = 0.1, # multiplicative buffer to radius for neighbor search algo nPerh = 1.01, @@ -50,7 +47,7 @@ # integration IntegratorConstructor = VerletIntegrator, - stepsPerCollision = 100, # replaces CFL for DEM + stepsPerCollision = 50, goalTime = 25.0, dt = 1e-8, dtMin = 1.0e-8, @@ -73,15 +70,16 @@ dataDir = "dumps-DEM-impactingSquares-3d", # ats - checkRestart = False, - checkConservation = False, # turn on error checking for momentum conservation - conservationErrorThreshold = 2e-14, # relative error for momentum conservation + checkRestitutionCoefficient = False, + threshRestitutionCoefficient = 0.05, + checkNaturalFrequency = False, + threshNaturalFrequency = 0.1, ) #------------------------------------------------------------------------------- # file things #------------------------------------------------------------------------------- -testName = "DEM-ImpactingSquares-3d" +testName = "DEM-PureTorsionTest-3d" restartDir = os.path.join(dataDir, "restarts") vizDir = os.path.join(dataDir, "visit") restartBaseName = os.path.join(restartDir, testName) @@ -140,7 +138,7 @@ xmax = ( particleRadius*5,particleRadius*5, particleRadius*5), nNodePerh = nPerh) - # # really simple bar shaped particle + # transforms particle properties from generator def DEMParticleGenerator(xi,yi,zi,Hi,mi,Ri): xout = [xi] yout = [yi] @@ -187,31 +185,27 @@ def DEMParticleGenerator(xi,yi,zi,Hi,mi,Ri): packages = [dem] -# #------------------------------------------------------------------------------- -# # PhysicsPackage : gravity -# #------------------------------------------------------------------------------- +#------------------------------------------------------------------------------- +# PhysicsPackage : gravity +#------------------------------------------------------------------------------- gravity = ConstantAcceleration(a0 = Vector(0.0, 0.0, -1.00), nodeList = nodes1) -#packages += [gravity] - +packages.append(gravity) +#------------------------------------------------------------------------------- +# Create boundary conditions. +#------------------------------------------------------------------------------- +# implement boundary condition using the DEM packages solid wall feature if useSolidBoundary: + solidWall = InfinitePlaneSolidBoundary(Vector(0.0, 0.0, 0.0), Vector( 0.0, 0.0, 1.0)) dem.appendSolidBoundary(solidWall) + +# implement boundary condition using Spheral's ghost particle reflection else: bcs = [ReflectingBoundary(Plane(Vector(0.0, 0.0, 0.0), Vector( 0.0, 0.0, 1.0)))] for package in packages: for bc in bcs: package.appendBoundary(bc) -# #------------------------------------------------------------------------------- -# # Create boundary conditions. -# #------------------------------------------------------------------------------- -# plane1 = Plane(Vector(0.0, 0.0, 0.0), Vector( 0.0, 0.0, 1.0)) -# bc1 = ReflectingBoundary(plane1) -# bcSet = [bc1] - -# for p in packages: -# for bc in bcSet: -# p.appendBoundary(bc) #------------------------------------------------------------------------------- # Fields and Variables @@ -229,10 +223,6 @@ def DEMParticleGenerator(xi,yi,zi,Hi,mi,Ri): uniqueIndex = db.DEMUniqueIndex omega = dem.omega - -#------------------------------------------------------------------------------- -# Initial Conditions -#------------------------------------------------------------------------------- #------------------------------------------------------------------------------- # Construct a time integrator, and add the physics packages. #------------------------------------------------------------------------------- @@ -261,19 +251,21 @@ def DEMParticleGenerator(xi,yi,zi,Hi,mi,Ri): output("integrator.verbose") #------------------------------------------------------------------------------- -# Periodic Work Function: Track conseravation +# Periodic Work Function: track resitution coefficient #------------------------------------------------------------------------------- - -conservation = TrackConservation(db, - dem, - verbose=True) - -periodicWork = [(conservation.periodicWorkFunction,1)] - -for i in range(nodes.numInternalNodes): - if i > 0: - velocity[0][i]=Vector(0,0,-0.1) - #omega[0][i]=Vector(0,0,particleVelocity/particleRadius) +class OmegaTracker: + def __init__(self): + self.maxOmega = 0.0 + self.period = 0.0 + self.omegan = 0.0 + def periodicWorkFunction(self,cycle,time,dt): + omegai = omega[0][1][2] + if omegai < self.maxOmega: + self.maxOmega = omegai + self.period = time - 15 + self.omegan = pi / self.period +omegaTracker = OmegaTracker() +periodicWork=[(omegaTracker.periodicWorkFunction,1)] #------------------------------------------------------------------------------- # Make the problem controller. @@ -296,7 +288,6 @@ def DEMParticleGenerator(xi,yi,zi,Hi,mi,Ri): periodicWork=periodicWork) output("control") -#control.redistribute = PeanoHilbertOrderRedistributeNodes(db.maxKernelExtent,workBalance=False) #------------------------------------------------------------------------------- # Advance to the end time. #------------------------------------------------------------------------------- @@ -304,10 +295,64 @@ def DEMParticleGenerator(xi,yi,zi,Hi,mi,Ri): if not steps is None: control.step(steps) else: - control.advance(10.0, maxSteps) - # for i in range(nodes.numInternalNodes): - # if i > 0: - # velocity[0][i]=Vector(0,0,0) - # omega[0][i]=Vector(0,0,particleVelocity/particleRadius) - control.advance(10.0+goalTime, maxSteps) + # settle dem particles on solid bc in grav field + control.advance(15.0, maxSteps) + + # give 3 particles torsional rotation + for i in range(nodes.numInternalNodes): + if i > 0: + velocity[0][i]=Vector(0,0,0) + omega[0][i]=Vector(0,0,particleVelocity/particleRadius) + + # run to the goal time + control.advance(15.0+goalTime, maxSteps) + +if checkRestitutionCoefficient or checkNaturalFrequency: + + # get the sliding damping constant + beta = pi/log(min(max(tangentialRestitutionCoefficient,1e-4),0.9999)) + nu = 2*tangentialSpringConstant/(1+beta*beta) + mi = particleDensity * 4.0/3.0*np.pi*particleRadius**3 + C = sqrt(2*mi*nu) + + # get the natural frequency and effective resitution coefficient of torsion + omegan = sqrt(5*tangentialSpringConstant*shapeFactor*shapeFactor/mi) + alphan = 5 * C*shapeFactor*shapeFactor/ (2 * mi) + analyticRestitution = exp(-alphan*pi/sqrt(omegan**2-alphan**2)) + + print("") + print("==============================================================") + + if checkNaturalFrequency: + print("") + print(" Checking Torsional Natural Frequency ") + print("") + print(" analytic natural freq : %g" % omegan) + print(" numerical natural freq : %g" % omegaTracker.omegan) + + relativeErrorNaturalFrequency = (abs(omegaTracker.omegan-omegan)/omegan) + + print(" relative error : %g" % relativeErrorNaturalFrequency) + + if relativeErrorNaturalFrequency > threshNaturalFrequency: + raise ValueError(" natural frequency is not within error bounds ") + + if checkRestitutionCoefficient: + + numericalRestitutionCoefficient = (-omegaTracker.maxOmega/particleVelocity*particleRadius) + + print("") + print(" Checking Torsional Restitution Coefficient ") + print("") + print(" analytic restitution coefficient : %g" % analyticRestitution) + print(" numerical restitution coefficient : %g" % numericalRestitutionCoefficient) + + relativeErrorRestitution = (abs(numericalRestitutionCoefficient-analyticRestitution)/analyticRestitution) + + print(" relative error : %g" % relativeErrorRestitution) + + if relativeErrorRestitution > threshRestitutionCoefficient: + raise ValueError(" restitution coefficient is not within error bounds ") + + print("==============================================================") \ No newline at end of file From 913600c10cd84c2faff79ac2c3e4fcef7c81907f Mon Sep 17 00:00:00 2001 From: jmpearl Date: Tue, 19 Mar 2024 10:42:49 -0700 Subject: [PATCH 24/38] DEM update release notes --- RELEASE_NOTES.md | 10 ++++++++-- 1 file changed, 8 insertions(+), 2 deletions(-) diff --git a/RELEASE_NOTES.md b/RELEASE_NOTES.md index 45896fb40..0c22251bf 100644 --- a/RELEASE_NOTES.md +++ b/RELEASE_NOTES.md @@ -9,13 +9,19 @@ Notable changes include: * Adding an optional second-stage problem start-up hook to the Physics package interface: Physics::initializeProblemStartupDependencies. The idea is to keep basic sizing of arrays and such in the first stage (Physics::initializeProblemStartup), while this new hook is used for updating any initial Physics state (and therefore provides a State and StateDerivatives object). - + * DEM + * new field list to track max particle overlap + * user can optional turn off fast time stepping + * Build changes / improvements: * - + * Bug Fixes / improvements: * Fixed bug with ConstantBoundary in the presence of porosity with the new porosity models introduced in v2024.01.00. * Initial volumes for damage models were incorrectly not taking into account pore space when computing failure statistics for seeding flaws. Fixed. + * DEM + * fixed bug in solid boundary unique indices that causes particle sticking + * fixed bug in solid boundary update policies Version v2024.01.00 -- Release date 2024-01-19 ============================================== From a31b02fdedb468e88ac8a07ea2804e6430566204 Mon Sep 17 00:00:00 2001 From: jmpearl Date: Tue, 19 Mar 2024 14:58:31 -0700 Subject: [PATCH 25/38] DEM bugfix misnamed ats --- .../{pureTorision-3d.py => pureTorsion-3d.py} | 10 ++++------ 1 file changed, 4 insertions(+), 6 deletions(-) rename tests/functional/DEM/LinearSpringDEM/SolidBoundaryCondition/{pureTorision-3d.py => pureTorsion-3d.py} (97%) diff --git a/tests/functional/DEM/LinearSpringDEM/SolidBoundaryCondition/pureTorision-3d.py b/tests/functional/DEM/LinearSpringDEM/SolidBoundaryCondition/pureTorsion-3d.py similarity index 97% rename from tests/functional/DEM/LinearSpringDEM/SolidBoundaryCondition/pureTorision-3d.py rename to tests/functional/DEM/LinearSpringDEM/SolidBoundaryCondition/pureTorsion-3d.py index 9fffdd123..0302d1c63 100644 --- a/tests/functional/DEM/LinearSpringDEM/SolidBoundaryCondition/pureTorision-3d.py +++ b/tests/functional/DEM/LinearSpringDEM/SolidBoundaryCondition/pureTorsion-3d.py @@ -1,6 +1,6 @@ -#ATS:DEM3dPureTorsion = test( SELF, "--clearDirectories True --checkRestitutionCoefficient True --checkNaturalFrequency True", label="DEM puretorsion test -- 3-D (serial)") +#ATS:DEM3dPT = test( SELF, "--clearDirectories True --checkRestitutionCoefficient True --checkNaturalFrequency True", label="DEM pure torsion test -- 3-D (serial)") -import os, sys, shutil, mpi, random +import os, sys, shutil, mpi from math import * from Spheral3d import * from SpheralTestUtilities import * @@ -10,15 +10,13 @@ import numpy as np -sys.path.insert(0, '..') -from DEMConservationTracker import TrackConservation3d as TrackConservation - if mpi.procs > 1: from PeanoHilbertDistributeNodes import distributeNodes3d else: from DistributeNodes import distributeNodes3d title("DEM 3d Pure Torsion Test") + # This tests the natural freq and restitution coefficient for Zhang et.al. formulation #------------------------------------------------------------------------------- # Generic problem parameters @@ -67,7 +65,7 @@ restoreCycle = None, restartStep = 10000, redistributeStep = 100000000000000, - dataDir = "dumps-DEM-impactingSquares-3d", + dataDir = "dumps-DEM-PureTorsionTest-3d", # ats checkRestitutionCoefficient = False, From 41c8178c5eb4573cb013ab6e868e19d66f27f4b6 Mon Sep 17 00:00:00 2001 From: jmpearl Date: Wed, 20 Mar 2024 08:20:51 -0700 Subject: [PATCH 26/38] DEM ats bugfix for particle-wall and more checks added to particle-particle --- .../singleParticleBoundaryCollision-3d.py | 6 ++--- .../twoParticleCollision-3d.py | 22 +++++++++++++++++-- 2 files changed, 23 insertions(+), 5 deletions(-) diff --git a/tests/functional/DEM/LinearSpringDEM/SolidBoundaryCondition/singleParticleBoundaryCollision-3d.py b/tests/functional/DEM/LinearSpringDEM/SolidBoundaryCondition/singleParticleBoundaryCollision-3d.py index f02deb3cb..c2e04c9ff 100644 --- a/tests/functional/DEM/LinearSpringDEM/SolidBoundaryCondition/singleParticleBoundaryCollision-3d.py +++ b/tests/functional/DEM/LinearSpringDEM/SolidBoundaryCondition/singleParticleBoundaryCollision-3d.py @@ -301,16 +301,16 @@ if boolCheckSlidingFrictionX: if omega[0][0].magnitude() > omega0: raise ValueError("particles are rotating faster post-collision") - if omega[0][0].y > omegaThreshold or omega[0][0].z > omegaThreshold: + if abs(omega[0][0].x) > omegaThreshold or abs(omega[0][0].z) > omegaThreshold: raise ValueError("erroneous spin-up in perpendicular direction") if boolCheckSlidingFrictionY: if omega[0][0].magnitude() > omega0: raise ValueError("particles are rotating faster post-collision") - if omega[0][0].x > omegaThreshold or omega[0][0].z > omegaThreshold: + if abs(omega[0][0].y) > omegaThreshold or abs(omega[0][0].z) > omegaThreshold: raise ValueError("erroneous spin-up in perpendicular direction") if boolCheckTorsionalFriction: if omega[0][0].magnitude() > omega0: raise ValueError("particles are rotating faster post-collision") - if omega[0][0].x > omegaThreshold or omega[0][0].y > omegaThreshold: + if abs(omega[0][0].x) > omegaThreshold or abs(omega[0][0].y) > omegaThreshold: raise ValueError("erroneous spin-up in perpendicular direction") diff --git a/tests/functional/DEM/LinearSpringDEM/TwoParticleCollision/twoParticleCollision-3d.py b/tests/functional/DEM/LinearSpringDEM/TwoParticleCollision/twoParticleCollision-3d.py index 9030f7422..33fbcb2f0 100644 --- a/tests/functional/DEM/LinearSpringDEM/TwoParticleCollision/twoParticleCollision-3d.py +++ b/tests/functional/DEM/LinearSpringDEM/TwoParticleCollision/twoParticleCollision-3d.py @@ -342,15 +342,33 @@ if conservation.deltaRotationalMomentumZ() > conservationErrorThreshold: raise ValueError("rotational momentum -z conservation error, %g, exceeds bounds" % conservation.deltaRotationalMomentumZ()) -if boolCheckSlidingFriction or boolCheckRollingFriction or boolCheckTorsionalFriction: +if boolCheckSlidingFriction or boolCheckRollingFriction: # check for non-physical behavior #------------------------------------------------------------- if omega[0][0].magnitude()+omega[0][1].magnitude() > 2*omega0: raise ValueError("particles are rotating faster post-collision") - + if abs(omega[0][0].x) > omegaThreshold or abs(omega[0][0].y) > omegaThreshold: + raise ValueError("erroneous spin-up of particle 0 in perpendicular direction") + if abs(omega[0][1].x) > omegaThreshold or abs(omega[0][1].y) > omegaThreshold: + raise ValueError("erroneous spin-up of particle 1 in perpendicular direction") + +if boolCheckTorsionalFriction: +# check for non-physical behavior +#------------------------------------------------------------- + if omega[0][0].magnitude()+omega[0][1].magnitude() > 2*omega0: + raise ValueError("particles are rotating faster post-collision") + if abs(omega[0][0].z) > omegaThreshold or abs(omega[0][0].y) > omegaThreshold: + raise ValueError("erroneous spin-up of particle 0 in perpendicular direction") + if abs(omega[0][1].z) > omegaThreshold or abs(omega[0][1].y) > omegaThreshold: + raise ValueError("erroneous spin-up of particle 1 in perpendicular direction") + if boolCheckTorsionalObjectivity: # to satify objectivity omega (along axis) should not change when equal #------------------------------------------------------------- omegaError = (2*omega0 - omega[0][0][0] - omega[0][1][0]) / (2*omega0) if omegaError > torsionalObjectivityThreshold: raise ValueError("torsional objectivity failure with relative angular velocity error, %g, exceeds bounds" % omegaError) + if abs(omega[0][0].z) > omegaThreshold or abs(omega[0][0].y) > omegaThreshold: + raise ValueError("erroneous spin-up of particle 0 in perpendicular direction") + if abs(omega[0][1].z) > omegaThreshold or abs(omega[0][1].y) > omegaThreshold: + raise ValueError("erroneous spin-up of particle 1 in perpendicular direction") \ No newline at end of file From 3c4a248b42e0e8b3f6622d3f575b7b22dface650 Mon Sep 17 00:00:00 2001 From: jmpearl Date: Wed, 20 Mar 2024 08:24:48 -0700 Subject: [PATCH 27/38] DEM ats particle-particle bugfix --- tests/functional/DEM/DEMTests.ats | 4 ++-- .../TwoParticleCollision/twoParticleCollision-3d.py | 1 + 2 files changed, 3 insertions(+), 2 deletions(-) diff --git a/tests/functional/DEM/DEMTests.ats b/tests/functional/DEM/DEMTests.ats index be2784348..4d269c4a6 100644 --- a/tests/functional/DEM/DEMTests.ats +++ b/tests/functional/DEM/DEMTests.ats @@ -10,7 +10,7 @@ source("LinearSpringDEM/TwoParticleCollision/fiveParticleCollision-2d.py") source("LinearSpringDEM/ImpactingSquares/impactingSquares-2d.py") source("LinearSpringDEM/ImpactingSquares/impactingSquares-3d.py") -# test restitution coefficient and friction against solid boundaries +# test restitution coefficient, natural frequency, and friction against solid boundaries source("LinearSpringDEM/SolidBoundaryCondition/singleParticleBoundaryCollision-3d.py") source("LinearSpringDEM/SolidBoundaryCondition/finiteSolidPlanes-3d.py") -source("LinearSpringDEM/SolidBoundaryCondition/pureTorsion-3d.py") +source("LinearSpringDEM/SolidBoundaryCondition/pureTorsion-3d.py") \ No newline at end of file diff --git a/tests/functional/DEM/LinearSpringDEM/TwoParticleCollision/twoParticleCollision-3d.py b/tests/functional/DEM/LinearSpringDEM/TwoParticleCollision/twoParticleCollision-3d.py index 33fbcb2f0..27192d7dd 100644 --- a/tests/functional/DEM/LinearSpringDEM/TwoParticleCollision/twoParticleCollision-3d.py +++ b/tests/functional/DEM/LinearSpringDEM/TwoParticleCollision/twoParticleCollision-3d.py @@ -77,6 +77,7 @@ checkConservation = False, # turn on error checking for momentum conservation restitutionErrorThreshold = 0.02, # relative error actual restitution vs nominal conservationErrorThreshold = 1e-15, # relative error for momentum conservation + omegaThreshold = 1e-14, # tolerance for erroneous spin up in inactive directions torsionalObjectivityThreshold = 1e-10 # relative error bounds on torsion objectivity test ) From fd350bd4c11dd73208be22e144cff8a50fb0379e Mon Sep 17 00:00:00 2001 From: jmpearl Date: Wed, 3 Apr 2024 06:34:55 -0700 Subject: [PATCH 28/38] adding a test for simple constant moving bc --- tests/functional/DEM/DEMTests.ats | 5 +- ...ingleParticleMovingBCRestitutionTest-3d.py | 265 ++++++++++++++++++ 2 files changed, 269 insertions(+), 1 deletion(-) create mode 100644 tests/functional/DEM/LinearSpringDEM/MovingSolidBoundaries/singleParticleMovingBCRestitutionTest-3d.py diff --git a/tests/functional/DEM/DEMTests.ats b/tests/functional/DEM/DEMTests.ats index 4d269c4a6..a69ed5a01 100644 --- a/tests/functional/DEM/DEMTests.ats +++ b/tests/functional/DEM/DEMTests.ats @@ -13,4 +13,7 @@ source("LinearSpringDEM/ImpactingSquares/impactingSquares-3d.py") # test restitution coefficient, natural frequency, and friction against solid boundaries source("LinearSpringDEM/SolidBoundaryCondition/singleParticleBoundaryCollision-3d.py") source("LinearSpringDEM/SolidBoundaryCondition/finiteSolidPlanes-3d.py") -source("LinearSpringDEM/SolidBoundaryCondition/pureTorsion-3d.py") \ No newline at end of file +source("LinearSpringDEM/SolidBoundaryCondition/pureTorsion-3d.py") + +# test moving solid boundaries +source("LinearSpringDEM/MovingSolidBoundaries/singleParticleBoundaryCollision-3d.py") \ No newline at end of file diff --git a/tests/functional/DEM/LinearSpringDEM/MovingSolidBoundaries/singleParticleMovingBCRestitutionTest-3d.py b/tests/functional/DEM/LinearSpringDEM/MovingSolidBoundaries/singleParticleMovingBCRestitutionTest-3d.py new file mode 100644 index 000000000..22a4e0f39 --- /dev/null +++ b/tests/functional/DEM/LinearSpringDEM/MovingSolidBoundaries/singleParticleMovingBCRestitutionTest-3d.py @@ -0,0 +1,265 @@ +#ATS:DEM3dMSPBC1 = test( SELF, "--clearDirectories True --boolCheckRestitutionCoefficient True --normalRestitutionCoefficient 1.0 --steps 100", label="DEM perfectly elastic collision with solid boundary -- 3-D (serial)") +#ATS:DEM3dMSPBC2 = test( SELF, "--clearDirectories True --boolCheckRestitutionCoefficient True --normalRestitutionCoefficient 0.5 --steps 100", label="DEM inelastic collision with solid boundary -- 3-D (serial)") + +import os, shutil, mpi +from math import * +from Spheral3d import * +from SpheralTestUtilities import * +from findLastRestart import * +from GenerateNodeDistribution3d import * +from GenerateDEMfromSPHGenerator import GenerateDEMfromSPHGenerator3d + +if mpi.procs > 1: + from PeanoHilbertDistributeNodes import distributeNodes3d +else: + from DistributeNodes import distributeNodes3d + +title("DEM Moving Boundary Restitution Coefficient Test") + +#------------------------------------------------------------------------------- +# Generic problem parameters +#------------------------------------------------------------------------------- +commandLine(vImpact = 1.0, # impact velocity (negative z direction) + vbc = 0.2, # velocity of bc (positive z direction) + + h0 = 1.00, # initial height above the solid bc plane + radius = 0.95, # particle radius + normalSpringConstant=10000.0, # spring constant for LDS model + normalRestitutionCoefficient=1.00, # restitution coefficient to get damping const + tangentialSpringConstant=2857.0, # spring constant for LDS model + tangentialRestitutionCoefficient=0.55, # restitution coefficient to get damping const + dynamicFriction = 1.0, # static friction coefficient sliding + staticFriction = 1.0, # dynamic friction coefficient sliding + rollingFriction = 1.05, # static friction coefficient for rolling + torsionalFriction = 1.3, # static friction coefficient for torsion + cohesiveTensileStrength = 0.0, # units of pressure + shapeFactor = 0.5, # in [0,1] shape factor from Zhang 2018, 0 - no torsion or rolling + + neighborSearchBuffer = 0.1, # multiplicative buffer to radius for neighbor search algo + + # integration + IntegratorConstructor = VerletIntegrator, # Verlet one integrator to garenteee conservation + stepsPerCollision = 50, # replaces CFL for DEM + goalTime = 3.0, + dt = 1e-8, + dtMin = 1.0e-8, + dtMax = 0.1, + dtGrowth = 2.0, + steps = None, + maxSteps = None, + statsStep = 10, + domainIndependent = False, + rigorousBoundaries = False, + dtverbose = False, + + # output control + vizCycle = None, + vizTime = 0.1, + clearDirectories = False, + restoreCycle = None, + restartStep = 1000, + redistributeStep = 500, + dataDir = "dumps-DEM-particle-boundary-3d", + + # ats parameters + boolCheckRestitutionCoefficient=False, # turn on error checking for restitution coefficient + boolCheckSlidingFrictionX=False, # checks sliding friction reduces relative rotation + boolCheckSlidingFrictionY=False, # checks rolling friction reduces relative rotation + boolCheckTorsionalFriction=False, # checks torsional friction reduces relative rotation + restitutionErrorThreshold = 0.02, # relative error actual restitution vs nominal + omegaThreshold = 1e-14, # theshold for perpendicular components that should stay zero + ) + +#------------------------------------------------------------------------------- +# check for bad inputs +#------------------------------------------------------------------------------- +assert mpi.procs == 1 +assert h0 > radius +assert shapeFactor <= 1.0 and shapeFactor >= 0.0 +assert dynamicFriction >= 0.0 +assert staticFriction >= 0.0 +assert torsionalFriction >= 0.0 +assert rollingFriction >= 0.0 +assert cohesiveTensileStrength >= 0.0 + +#------------------------------------------------------------------------------- +# file things +#------------------------------------------------------------------------------- +testName = "DEM-SingleParticleMovingBoundaryCollision-3d" + +dataDir = os.path.join(dataDir, + "restitutionCoefficient=%s" % normalRestitutionCoefficient) + +restartDir = os.path.join(dataDir, "restarts") +vizDir = os.path.join(dataDir, "visit") +restartBaseName = os.path.join(restartDir, testName) +vizBaseName = testName + +if vizCycle is None and vizTime is None: + vizBaseName=None + +#------------------------------------------------------------------------------- +# Check if the necessary output directories exist. If not, create them. +#------------------------------------------------------------------------------- +if mpi.rank == 0: + if clearDirectories and os.path.exists(dataDir): + shutil.rmtree(dataDir) + if not os.path.exists(restartDir): + os.makedirs(restartDir) + if not os.path.exists(vizDir): + os.makedirs(vizDir) +mpi.barrier() + +#------------------------------------------------------------------------------- +# If we're restarting, find the set of most recent restart files. +#------------------------------------------------------------------------------- +if restoreCycle is None: + restoreCycle = findLastRestart(restartBaseName) + +#------------------------------------------------------------------------------- +# This doesn't really matter kernel filler for neighbor algo +#------------------------------------------------------------------------------- +WT = TableKernel(WendlandC2Kernel(), 1000) + +#------------------------------------------------------------------------------- +# Make the NodeList. +#------------------------------------------------------------------------------- +units = CGuS() +nodes1 = makeDEMNodeList("nodeList1", + hmin = 1.0e-30, + hmax = 1.0e30, + hminratio = 100.0, + neighborSearchBuffer = neighborSearchBuffer, + kernelExtent = WT.kernelExtent) +nodeSet = [nodes1] +for nodes in nodeSet: + output("nodes.name") + output("nodes.hmin") + output("nodes.hmax") + output("nodes.hminratio") + output("nodes.nodesPerSmoothingScale") + +#------------------------------------------------------------------------------- +# Set the node properties. (gen 2 particles visit doesn't like just one) +#------------------------------------------------------------------------------- +generator0 = GenerateNodeDistribution3d(1, 1, 1, + rho = 1.0, + distributionType = "lattice", + xmin = (-1.0, -1.0, -1+h0), + xmax = (1.0, 1.0, 1+h0)) + +generator1 = GenerateDEMfromSPHGenerator3d(WT, + generator0) +distributeNodes3d((nodes1, generator1)) + +#------------------------------------------------------------------------------- +# Construct a DataBase to hold our node list +#------------------------------------------------------------------------------- +db = DataBase() +output("db") +for nodes in nodeSet: + db.appendNodeList(nodes) +output("db.numNodeLists") +output("db.numDEMNodeLists") +output("db.numFluidNodeLists") + + +#------------------------------------------------------------------------------- +# DEM +#------------------------------------------------------------------------------- +dem = DEM(db, + normalSpringConstant = normalSpringConstant, + normalRestitutionCoefficient = normalRestitutionCoefficient, + tangentialSpringConstant = tangentialSpringConstant, + tangentialRestitutionCoefficient = tangentialRestitutionCoefficient, + dynamicFrictionCoefficient = dynamicFriction, + staticFrictionCoefficient = staticFriction, + rollingFrictionCoefficient = rollingFriction, + torsionalFrictionCoefficient = torsionalFriction, + cohesiveTensileStrength =cohesiveTensileStrength, + shapeFactor = shapeFactor, + stepsPerCollision = stepsPerCollision) + +packages = [dem] + +# set the solid bcs and add to dem package +solidWall = InfinitePlaneSolidBoundary(Vector(0.0, 0.0, 0.0), Vector( 0.0, 0.0, 1.0)) +solidWall.velocity = Vector(0.0,0.0,vbc) + +dem.appendSolidBoundary(solidWall) + +#------------------------------------------------------------------------------- +# initial conditions +#------------------------------------------------------------------------------- +velocity = nodes1.velocity() +particleRadius = nodes1.particleRadius() + +velocity[0] = Vector(0.0,0.0,-vImpact) +particleRadius[0] = radius + +#------------------------------------------------------------------------------- +# Construct a time integrator, and add the physics packages. +#------------------------------------------------------------------------------- +integrator = IntegratorConstructor(db) +for p in packages: + integrator.appendPhysicsPackage(p) +integrator.lastDt = dt +integrator.dtMin = dtMin +integrator.dtMax = dtMax +integrator.dtGrowth = dtGrowth +integrator.domainDecompositionIndependent = domainIndependent +integrator.verbose = dtverbose +integrator.rigorousBoundaries = rigorousBoundaries + +integrator.cullGhostNodes = False + +output("integrator") +output("integrator.havePhysicsPackage(dem)") +output("integrator.lastDt") +output("integrator.dtMin") +output("integrator.dtMax") +output("integrator.dtGrowth") +output("integrator.domainDecompositionIndependent") +output("integrator.rigorousBoundaries") +output("integrator.verbose") + +#------------------------------------------------------------------------------- +# Make the problem controller. +#------------------------------------------------------------------------------- +from SpheralPointmeshSiloDump import dumpPhysicsState +control = SpheralController(integrator, WT, + iterateInitialH = False, + initializeDerivatives = True, + statsStep = statsStep, + restartStep = restartStep, + restartBaseName = restartBaseName, + restoreCycle = restoreCycle, + vizBaseName = vizBaseName, + vizMethod=dumpPhysicsState, + vizDir = vizDir, + vizStep = vizCycle, + vizTime = vizTime) +output("control") + +#------------------------------------------------------------------------------- +# Advance to the end time. +#------------------------------------------------------------------------------- +if not steps is None: + control.step(steps) +else: + control.advance(goalTime, maxSteps) + +#------------------------------------------------------------------------------- +# Great success? +#------------------------------------------------------------------------------- +if boolCheckRestitutionCoefficient: +# check our restitution coefficient is correct +#------------------------------------------------------------- + vijPostImpact = -velocity[0].z + vbc + vijPreImpact = vImpact + vbc + restitutionEff = vijPostImpact/vijPreImpact + restitutionError = abs(restitutionEff + normalRestitutionCoefficient)/normalRestitutionCoefficient + if restitutionError > restitutionErrorThreshold: + print(" final velocity = {0}".format(vijPostImpact)) + print(" initial velocity = {0}".format(vijPreImpact)) + raise ValueError(" relative restitution coefficient error, %g, exceeds bounds" % restitutionError) \ No newline at end of file From 85bbb293c0427c4c27b65e024458f77515657c60 Mon Sep 17 00:00:00 2001 From: jmpearl Date: Wed, 3 Apr 2024 07:23:38 -0700 Subject: [PATCH 29/38] DEM: updating some solid bcs to allow rectilinear motion --- src/DEM/SolidBoundary/CylinderSolidBoundary.cc | 7 +++++++ src/DEM/SolidBoundary/SphereSolidBoundary.cc | 13 ++++++++++++- 2 files changed, 19 insertions(+), 1 deletion(-) diff --git a/src/DEM/SolidBoundary/CylinderSolidBoundary.cc b/src/DEM/SolidBoundary/CylinderSolidBoundary.cc index 80c5bd30d..0a3f84520 100644 --- a/src/DEM/SolidBoundary/CylinderSolidBoundary.cc +++ b/src/DEM/SolidBoundary/CylinderSolidBoundary.cc @@ -55,6 +55,13 @@ void CylinderSolidBoundary:: registerState(DataBase& dataBase, State& state) { + const auto boundaryKey = "CylinderSolidBoundary_" + std::to_string(std::abs(this->uniqueIndex())); + const auto pointKey = boundaryKey +"_point"; + const auto velocityKey = boundaryKey +"_velocity"; + //const auto normalKey = boundaryKey +"_normal"; + state.enrollAny(pointKey,mPoint); + state.enrollAny(pointKey,mVelocity); + //state.enrollAny(pointKey,mNormal); } template diff --git a/src/DEM/SolidBoundary/SphereSolidBoundary.cc b/src/DEM/SolidBoundary/SphereSolidBoundary.cc index df79479eb..25f8f68ec 100644 --- a/src/DEM/SolidBoundary/SphereSolidBoundary.cc +++ b/src/DEM/SolidBoundary/SphereSolidBoundary.cc @@ -75,7 +75,17 @@ template void SphereSolidBoundary:: registerState(DataBase& dataBase, - State& state) { + State& state) { + + const auto boundaryKey = "SphereSolidBoundary_" + std::to_string(std::abs(this->uniqueIndex())); + const auto pointKey = boundaryKey +"_point"; + const auto clipPointKey = boundaryKey +"_clipPoint"; + const auto velocityKey = boundaryKey +"_velocity"; + + state.enrollAny(pointKey,mCenter); + state.enrollAny(clipPointKey,mClipPoint); + state.enrollAny(pointKey,mVelocity); + } template @@ -83,6 +93,7 @@ void SphereSolidBoundary:: update(const double multiplier, const double t, const double dt) { mCenter += multiplier*mVelocity; + mClipPoint += multiplier*mVelocity; } From bd43bd07cb006a7622726af4da79f3a5a6deb9e9 Mon Sep 17 00:00:00 2001 From: jmpearl Date: Wed, 3 Apr 2024 07:40:22 -0700 Subject: [PATCH 30/38] DEM: test rename --- tests/functional/DEM/DEMTests.ats | 2 +- ...onTest-3d.py => singleParticleMovingBoundaryCollision-3d.py} | 0 2 files changed, 1 insertion(+), 1 deletion(-) rename tests/functional/DEM/LinearSpringDEM/MovingSolidBoundaries/{singleParticleMovingBCRestitutionTest-3d.py => singleParticleMovingBoundaryCollision-3d.py} (100%) diff --git a/tests/functional/DEM/DEMTests.ats b/tests/functional/DEM/DEMTests.ats index a69ed5a01..bcb3ff032 100644 --- a/tests/functional/DEM/DEMTests.ats +++ b/tests/functional/DEM/DEMTests.ats @@ -16,4 +16,4 @@ source("LinearSpringDEM/SolidBoundaryCondition/finiteSolidPlanes-3d.py") source("LinearSpringDEM/SolidBoundaryCondition/pureTorsion-3d.py") # test moving solid boundaries -source("LinearSpringDEM/MovingSolidBoundaries/singleParticleBoundaryCollision-3d.py") \ No newline at end of file +source("LinearSpringDEM/MovingSolidBoundaries/singleParticleMovingBoundaryCollision-3d.py") \ No newline at end of file diff --git a/tests/functional/DEM/LinearSpringDEM/MovingSolidBoundaries/singleParticleMovingBCRestitutionTest-3d.py b/tests/functional/DEM/LinearSpringDEM/MovingSolidBoundaries/singleParticleMovingBoundaryCollision-3d.py similarity index 100% rename from tests/functional/DEM/LinearSpringDEM/MovingSolidBoundaries/singleParticleMovingBCRestitutionTest-3d.py rename to tests/functional/DEM/LinearSpringDEM/MovingSolidBoundaries/singleParticleMovingBoundaryCollision-3d.py From f845f16ac8a45ce764f6aa2d435e33cf0cdb506f Mon Sep 17 00:00:00 2001 From: jmpearl Date: Mon, 8 Apr 2024 10:08:07 -0700 Subject: [PATCH 31/38] DEM: solid bc issue reproducer --- .../MovingSolidBoundaries/rotatingDrum-2d.py | 293 ++++++++++++++++++ 1 file changed, 293 insertions(+) create mode 100644 tests/functional/DEM/LinearSpringDEM/MovingSolidBoundaries/rotatingDrum-2d.py diff --git a/tests/functional/DEM/LinearSpringDEM/MovingSolidBoundaries/rotatingDrum-2d.py b/tests/functional/DEM/LinearSpringDEM/MovingSolidBoundaries/rotatingDrum-2d.py new file mode 100644 index 000000000..020c8456b --- /dev/null +++ b/tests/functional/DEM/LinearSpringDEM/MovingSolidBoundaries/rotatingDrum-2d.py @@ -0,0 +1,293 @@ +import os, shutil, mpi +from math import * +from Spheral2d import * +from SpheralTestUtilities import * +from findLastRestart import * +from GenerateNodeDistribution2d import * +from GenerateDEMfromSPHGenerator import GenerateDEMfromSPHGenerator2d + +if mpi.procs > 1: + from PeanoHilbertDistributeNodes import distributeNodes2d +else: + from DistributeNodes import distributeNodes2d + +title("DEM Rotating Drum") + +class RotatingDrumBoundary(SphereSolidBoundary2d): + def __init__(self, + center, + radius, + omega): + SphereSolidBoundary2d.__init__(self, + center = center, + radius = radius, + clipPoint = center + Vector(0.0,radius*1.1), + clipAxis = Vector(0,1)) + return + +print("THERE") +#------------------------------------------------------------------------------- +# Generic problem parameters +#------------------------------------------------------------------------------- +commandLine(omegaDrum = 0.1, # angular velocity of drum + radiusDrum = 10.0, # radius of the drum + yThresh = 0.2, # level to initial fill to + yClip = 0.0, # level to clip at after settling + g0 = 1.0, # grav acceleration + + radius = 0.5, # particle radius + normalSpringConstant=10000.0, # spring constant for LDS model + normalRestitutionCoefficient=1.00, # restitution coefficient to get damping const + tangentialSpringConstant=2857.0, # spring constant for LDS model + tangentialRestitutionCoefficient=0.55, # restitution coefficient to get damping const + dynamicFriction = 1.0, # static friction coefficient sliding + staticFriction = 1.0, # dynamic friction coefficient sliding + rollingFriction = 1.05, # static friction coefficient for rolling + torsionalFriction = 1.3, # static friction coefficient for torsion + cohesiveTensileStrength = 0.0, # units of pressure + shapeFactor = 0.5, # in [0,1] shape factor from Zhang 2018, 0 - no torsion or rolling + + neighborSearchBuffer = 0.1, # multiplicative buffer to radius for neighbor search algo + + # integration + IntegratorConstructor = VerletIntegrator, # Verlet one integrator to garenteee conservation + stepsPerCollision = 50, # replaces CFL for DEM + goalTime = 3.0, + dt = 1e-8, + dtMin = 1.0e-8, + dtMax = 0.1, + dtGrowth = 2.0, + steps = None, + maxSteps = None, + statsStep = 10, + domainIndependent = False, + rigorousBoundaries = False, + dtverbose = False, + + # output control + vizCycle = None, + vizTime = 0.1, + clearDirectories = False, + restoreCycle = None, + restartStep = 1000, + redistributeStep = 500, + dataDir = "dumps-DEM-rotating-drum-2d", + + # ats parameters + boolCheckRestitutionCoefficient=False, # turn on error checking for restitution coefficient + boolCheckSlidingFrictionX=False, # checks sliding friction reduces relative rotation + boolCheckSlidingFrictionY=False, # checks rolling friction reduces relative rotation + boolCheckTorsionalFriction=False, # checks torsional friction reduces relative rotation + restitutionErrorThreshold = 0.02, # relative error actual restitution vs nominal + omegaThreshold = 1e-14, # theshold for perpendicular components that should stay zero + ) + +#------------------------------------------------------------------------------- +# check for bad inputs +#------------------------------------------------------------------------------- +assert mpi.procs == 1 +assert radius < radiusDrum/2.0 +assert shapeFactor <= 1.0 and shapeFactor >= 0.0 +assert dynamicFriction >= 0.0 +assert staticFriction >= 0.0 +assert torsionalFriction >= 0.0 +assert rollingFriction >= 0.0 +assert cohesiveTensileStrength >= 0.0 + +#------------------------------------------------------------------------------- +# file things +#------------------------------------------------------------------------------- +testName = "DEM-RotatingDrum-2d" + +dataDir = os.path.join(dataDir, + "restitutionCoefficient=%s" % normalRestitutionCoefficient) + +restartDir = os.path.join(dataDir, "restarts") +vizDir = os.path.join(dataDir, "visit") +restartBaseName = os.path.join(restartDir, testName) +vizBaseName = testName + +if vizCycle is None and vizTime is None: + vizBaseName=None + +#------------------------------------------------------------------------------- +# Check if the necessary output directories exist. If not, create them. +#------------------------------------------------------------------------------- +if mpi.rank == 0: + if clearDirectories and os.path.exists(dataDir): + shutil.rmtree(dataDir) + if not os.path.exists(restartDir): + os.makedirs(restartDir) + if not os.path.exists(vizDir): + os.makedirs(vizDir) +mpi.barrier() + +#------------------------------------------------------------------------------- +# If we're restarting, find the set of most recent restart files. +#------------------------------------------------------------------------------- +if restoreCycle is None: + restoreCycle = findLastRestart(restartBaseName) + +#------------------------------------------------------------------------------- +# This doesn't really matter kernel filler for neighbor algo +#------------------------------------------------------------------------------- +WT = TableKernel(WendlandC2Kernel(), 1000) + +#------------------------------------------------------------------------------- +# Make the NodeList. +#------------------------------------------------------------------------------- +units = CGuS() +nodes1 = makeDEMNodeList("nodeList1", + hmin = 1.0e-30, + hmax = 1.0e30, + hminratio = 100.0, + neighborSearchBuffer = neighborSearchBuffer, + kernelExtent = WT.kernelExtent) +nodeSet = [nodes1] +for nodes in nodeSet: + output("nodes.name") + output("nodes.hmin") + output("nodes.hmax") + output("nodes.hminratio") + output("nodes.nodesPerSmoothingScale") + +#------------------------------------------------------------------------------- +# Set the node properties. (gen 2 particles visit doesn't like just one) +#------------------------------------------------------------------------------- +def rejecter(x,y,m,H): + xnew,ynew,mnew,Hnew = [],[],[],[] + for i in range(len(x)): + if y[i] < yThresh and sqrt(x[i]**2+y[i]**2) < radiusDrum - 1.10*radius: + xnew.append(x[i]) + ynew.append(y[i]) + mnew.append(m[i]) + Hnew.append(H[i]) + + return xnew,ynew,mnew,Hnew + +nx = int(2*radiusDrum / (radius * 1.02)) +ny = nx + +generator0 = GenerateNodeDistribution2d(nx, ny, + rho = 1.0, + distributionType = "lattice", + xmin = (-radiusDrum, -radiusDrum), + xmax = ( radiusDrum, radiusDrum), + rejecter=rejecter) + +generator1 = GenerateDEMfromSPHGenerator2d(WT, + generator0) +distributeNodes2d((nodes1, generator1)) + +#------------------------------------------------------------------------------- +# Construct a DataBase to hold our node list +#------------------------------------------------------------------------------- +db = DataBase() +output("db") +for nodes in nodeSet: + db.appendNodeList(nodes) +output("db.numNodeLists") +output("db.numDEMNodeLists") +output("db.numFluidNodeLists") + + +#------------------------------------------------------------------------------- +# Physics Package: DEM +#------------------------------------------------------------------------------- +dem = DEM(db, + normalSpringConstant = normalSpringConstant, + normalRestitutionCoefficient = normalRestitutionCoefficient, + tangentialSpringConstant = tangentialSpringConstant, + tangentialRestitutionCoefficient = tangentialRestitutionCoefficient, + dynamicFrictionCoefficient = dynamicFriction, + staticFrictionCoefficient = staticFriction, + rollingFrictionCoefficient = rollingFriction, + torsionalFrictionCoefficient = torsionalFriction, + cohesiveTensileStrength =cohesiveTensileStrength, + shapeFactor = shapeFactor, + stepsPerCollision = stepsPerCollision) + +packages = [dem] + +# set yup a circle (clip outside circle to deactivate clipping (i know not the best)) +solidWall = RotatingDrumBoundary(center = Vector(0.0, 0.0), + radius = radiusDrum, + omega = 0.0) + +# solidWall = SphereSolidBoundary(center = Vector(0.0, 0.0), +# radius = radiusDrum, +# clipPoint = Vector(0.0,radiusDrum*1.1), +# clipAxis = Vector(0,1)) + +print("here") +dem.appendSolidBoundary(solidWall) +print("HERE") +#------------------------------------------------------------------------------- +# Gravity: DEM +#------------------------------------------------------------------------------- +gravity = ConstantAcceleration(a0 = Vector(0.0,-g0), + nodeList = nodes1) +packages += [gravity] + +#------------------------------------------------------------------------------- +# initial conditions +#------------------------------------------------------------------------------- +# velocity = nodes1.velocity() +# particleRadius = nodes1.particleRadius() + +# velocity[0] = Vector(0.0,0.0,-vImpact) +# particleRadius[0] = radius + +#------------------------------------------------------------------------------- +# Construct a time integrator, and add the physics packages. +#------------------------------------------------------------------------------- +integrator = IntegratorConstructor(db) +for p in packages: + integrator.appendPhysicsPackage(p) +integrator.lastDt = dt +integrator.dtMin = dtMin +integrator.dtMax = dtMax +integrator.dtGrowth = dtGrowth +integrator.domainDecompositionIndependent = domainIndependent +integrator.verbose = dtverbose +integrator.rigorousBoundaries = rigorousBoundaries + +integrator.cullGhostNodes = False + +output("integrator") +output("integrator.havePhysicsPackage(dem)") +output("integrator.lastDt") +output("integrator.dtMin") +output("integrator.dtMax") +output("integrator.dtGrowth") +output("integrator.domainDecompositionIndependent") +output("integrator.rigorousBoundaries") +output("integrator.verbose") + +#------------------------------------------------------------------------------- +# Make the problem controller. +#------------------------------------------------------------------------------- +print("POOPSIE") +from SpheralPointmeshSiloDump import dumpPhysicsState +control = SpheralController(integrator, WT, + iterateInitialH = False, + initializeDerivatives = True, + statsStep = statsStep, + restartStep = restartStep, + restartBaseName = restartBaseName, + restoreCycle = restoreCycle, + vizBaseName = vizBaseName, + vizMethod=dumpPhysicsState, + vizDir = vizDir, + vizStep = vizCycle, + vizTime = vizTime) +output("control") + +#------------------------------------------------------------------------------- +# Advance to the end time. +#------------------------------------------------------------------------------- +if not steps is None: + control.step(steps) +else: + control.advance(goalTime, maxSteps) + From b5e2a768bd135d2ea4e299e7900ca04933264c36 Mon Sep 17 00:00:00 2001 From: jmpearl Date: Mon, 8 Apr 2024 10:09:09 -0700 Subject: [PATCH 32/38] updating some solid bc stuff --- src/DEM/SolidBoundary/SolidBoundaryBase.hh | 2 +- .../SolidBoundary/SolidBoundaryBaseInline.hh | 2 +- src/PYB11/DEM/SolidBoundaries.py | 27 ++++++++++++++----- 3 files changed, 22 insertions(+), 9 deletions(-) diff --git a/src/DEM/SolidBoundary/SolidBoundaryBase.hh b/src/DEM/SolidBoundary/SolidBoundaryBase.hh index 6960fd5cd..849c7668d 100644 --- a/src/DEM/SolidBoundary/SolidBoundaryBase.hh +++ b/src/DEM/SolidBoundary/SolidBoundaryBase.hh @@ -55,7 +55,7 @@ public: const double t, const double dt) = 0; - void uniqueIndex(const int uId); + void uniqueIndex(int uId); int uniqueIndex() const; private: diff --git a/src/DEM/SolidBoundary/SolidBoundaryBaseInline.hh b/src/DEM/SolidBoundary/SolidBoundaryBaseInline.hh index a784b94ac..8a4e4a12c 100644 --- a/src/DEM/SolidBoundary/SolidBoundaryBaseInline.hh +++ b/src/DEM/SolidBoundary/SolidBoundaryBaseInline.hh @@ -6,7 +6,7 @@ template inline void SolidBoundaryBase:: -uniqueIndex(const int uId){ +uniqueIndex(int uId){ mUniqueIndex = uId; } template diff --git a/src/PYB11/DEM/SolidBoundaries.py b/src/PYB11/DEM/SolidBoundaries.py index 621c6a09a..19f0da6fe 100644 --- a/src/PYB11/DEM/SolidBoundaries.py +++ b/src/PYB11/DEM/SolidBoundaries.py @@ -12,9 +12,11 @@ class SolidBoundaryBase: typedef typename %(Dimension)s::Scalar Scalar; typedef typename %(Dimension)s::Vector Vector; """ - def pyinit(): + def pyinit(self): "constructor for base class DEM solid boundary conditions" + uniqueIndex = PYB11property("int", "uniqueIndex", "uniqueIndex", doc="unique index for solid boundary") + PYB11inject(SolidBoundaryBaseAbstractMethods, SolidBoundaryBase, pure_virtual=True) #------------------------------------------------------------------------------- @@ -29,7 +31,8 @@ class InfinitePlaneSolidBoundary(SolidBoundaryBase): typedef typename %(Dimension)s::Vector Vector; """ - def pyinit(point = "const Vector&", + def pyinit(self, + point = "const Vector&", normal = "const Vector&"): "solid planar boundary" @@ -79,7 +82,8 @@ class RectangularPlaneSolidBoundary(SolidBoundaryBase): typedef typename %(Dimension)s::Tensor Tensor; """ - def pyinit(point = "const Vector&", + def pyinit(self, + point = "const Vector&", extent = "const Vector&", basis = "const Tensor&"): "solid planar boundary" @@ -130,7 +134,8 @@ class CircularPlaneSolidBoundary(SolidBoundaryBase): typedef typename %(Dimension)s::Vector Vector; """ - def pyinit(point = "const Vector&", + def pyinit(self, + point = "const Vector&", normal = "const Vector&", extent = "const Scalar"): "solid planar boundary" @@ -181,7 +186,8 @@ class CylinderSolidBoundary(SolidBoundaryBase): typedef typename %(Dimension)s::Vector Vector; """ - def pyinit(point = "const Vector&", + def pyinit(self, + point = "const Vector&", axis = "const Vector&", radius = "const Scalar", length = "const Scalar"): @@ -235,7 +241,8 @@ class SphereSolidBoundary(SolidBoundaryBase): typedef typename %(Dimension)s::Vector Vector; """ - def pyinit(center = "const Vector&", + def pyinit(self, + center = "const Vector&", radius = "const Scalar", clipPoint = "const Vector&", clipAxis = "const Vector&"): @@ -274,4 +281,10 @@ def distance(self, center = PYB11property("const Vector&", "center", "center", returnpolicy="reference_internal", doc="center of sphere") radius = PYB11property("Scalar", "radius", "radius", doc="radius of sphere") clipPoint = PYB11property("const Vector&", "clipPoint", "clipPoint", returnpolicy="reference_internal", doc="point on clip plane") - clipAxis = PYB11property("const Vector&", "clipAxis", "clipAxis", returnpolicy="reference_internal", doc="normal in clip plane") \ No newline at end of file + clipAxis = PYB11property("const Vector&", "clipAxis", "clipAxis", returnpolicy="reference_internal", doc="normal in clip plane") + +PYB11inject(SolidBoundaryBaseAbstractMethods, SphereSolidBoundary, virtual=True, pure_virtual=False) +PYB11inject(SolidBoundaryBaseAbstractMethods, InfinitePlaneSolidBoundary, virtual=True, pure_virtual=False) +PYB11inject(SolidBoundaryBaseAbstractMethods, RectangularPlaneSolidBoundary, virtual=True, pure_virtual=False) +PYB11inject(SolidBoundaryBaseAbstractMethods, CircularPlaneSolidBoundary, virtual=True, pure_virtual=False) +PYB11inject(SolidBoundaryBaseAbstractMethods, CylinderSolidBoundary, virtual=True, pure_virtual=False) From 891fff874a238928ca91e04dbd266e2699a94642 Mon Sep 17 00:00:00 2001 From: jmpearl Date: Mon, 15 Apr 2024 06:35:08 -0700 Subject: [PATCH 33/38] merge resolution --- RELEASE_NOTES.md | 5 ----- 1 file changed, 5 deletions(-) diff --git a/RELEASE_NOTES.md b/RELEASE_NOTES.md index 8c368fdfa..9ba54c5d5 100644 --- a/RELEASE_NOTES.md +++ b/RELEASE_NOTES.md @@ -28,13 +28,8 @@ Notable changes include: * user can optional turn off fast time stepping * Build changes / improvements: -<<<<<<< HEAD - * - -======= * Improved the target export functionality. ->>>>>>> develop * Bug Fixes / improvements: * Fixed bug with ConstantBoundary in the presence of porosity with the new porosity models introduced in v2024.01.00. * Updating header lists for including Spheral modules in external projects. From 557887ea62e7a07106c441e37bc178dc43ff3a8b Mon Sep 17 00:00:00 2001 From: jmpearl Date: Mon, 15 Apr 2024 11:02:25 -0700 Subject: [PATCH 34/38] DEM fixing a name collision --- src/DEM/LinearSpringDEM.cc | 4 ++-- .../CircularPlaneSolidBoundary.cc | 2 +- .../CircularPlaneSolidBoundary.hh | 2 +- .../SolidBoundary/CylinderSolidBoundary.cc | 2 +- .../SolidBoundary/CylinderSolidBoundary.hh | 2 +- .../InfinitePlaneSolidBoundary.cc | 2 +- .../InfinitePlaneSolidBoundary.hh | 2 +- .../RectangularPlaneSolidBoundary.cc | 2 +- .../RectangularPlaneSolidBoundary.hh | 2 +- src/DEM/SolidBoundary/SolidBoundaryBase.hh | 2 +- src/DEM/SolidBoundary/SphereSolidBoundary.cc | 2 +- src/DEM/SolidBoundary/SphereSolidBoundary.hh | 2 +- src/PYB11/DEM/SolidBoundaries.py | 20 +++++++++---------- .../DEM/SolidBoundaryBaseAbstractMethods.py | 2 +- 14 files changed, 24 insertions(+), 24 deletions(-) diff --git a/src/DEM/LinearSpringDEM.cc b/src/DEM/LinearSpringDEM.cc index d59adebe0..56e3e2f92 100644 --- a/src/DEM/LinearSpringDEM.cc +++ b/src/DEM/LinearSpringDEM.cc @@ -260,7 +260,7 @@ variableTimeStep(const DataBase& dataBase, //solid boundary and distance vector to particle i const auto& solidBoundary = solidBoundaries[bci]; const auto rib = solidBoundary->distance(ri); - const auto vb = solidBoundary->velocity(ri); + const auto vb = solidBoundary->localVelocity(ri); // Compare closing speed to separation const auto vib = vi-vb; @@ -712,7 +712,7 @@ evaluateDerivatives(const typename Dimension::Scalar /*time*/, const auto deltaTorsib = torsionalDisplacement(nodeListi,i)[contacti]; // velocity of boundary @ ri - const auto vb = solidBoundary->velocity(ri); + const auto vb = solidBoundary->localVelocity(ri); // line of action for the contact const auto rhatib = rib.unitVector(); diff --git a/src/DEM/SolidBoundary/CircularPlaneSolidBoundary.cc b/src/DEM/SolidBoundary/CircularPlaneSolidBoundary.cc index de6b6002a..41441195f 100644 --- a/src/DEM/SolidBoundary/CircularPlaneSolidBoundary.cc +++ b/src/DEM/SolidBoundary/CircularPlaneSolidBoundary.cc @@ -45,7 +45,7 @@ distance(const Vector& position) const { template typename Dimension::Vector CircularPlaneSolidBoundary:: -velocity(const Vector& position) const { +localVelocity(const Vector& position) const { return mVelocity; } diff --git a/src/DEM/SolidBoundary/CircularPlaneSolidBoundary.hh b/src/DEM/SolidBoundary/CircularPlaneSolidBoundary.hh index aaa75f9d1..15ee17b02 100644 --- a/src/DEM/SolidBoundary/CircularPlaneSolidBoundary.hh +++ b/src/DEM/SolidBoundary/CircularPlaneSolidBoundary.hh @@ -32,7 +32,7 @@ public: ~CircularPlaneSolidBoundary(); virtual Vector distance(const Vector& position) const override; - virtual Vector velocity(const Vector& position) const override; + virtual Vector localVelocity(const Vector& position) const override; virtual void registerState(DataBase& dataBase, State& state) override; diff --git a/src/DEM/SolidBoundary/CylinderSolidBoundary.cc b/src/DEM/SolidBoundary/CylinderSolidBoundary.cc index 0a3f84520..fffbae75f 100644 --- a/src/DEM/SolidBoundary/CylinderSolidBoundary.cc +++ b/src/DEM/SolidBoundary/CylinderSolidBoundary.cc @@ -46,7 +46,7 @@ distance(const Vector& position) const { template typename Dimension::Vector CylinderSolidBoundary:: -velocity(const Vector& position) const { +localVelocity(const Vector& position) const { return mVelocity; } diff --git a/src/DEM/SolidBoundary/CylinderSolidBoundary.hh b/src/DEM/SolidBoundary/CylinderSolidBoundary.hh index e4cbe30b0..b53e046a2 100644 --- a/src/DEM/SolidBoundary/CylinderSolidBoundary.hh +++ b/src/DEM/SolidBoundary/CylinderSolidBoundary.hh @@ -33,7 +33,7 @@ public: ~CylinderSolidBoundary(); virtual Vector distance(const Vector& position) const override; - virtual Vector velocity(const Vector& position) const override; + virtual Vector localVelocity(const Vector& position) const override; virtual void registerState(DataBase& dataBase, State& state) override; diff --git a/src/DEM/SolidBoundary/InfinitePlaneSolidBoundary.cc b/src/DEM/SolidBoundary/InfinitePlaneSolidBoundary.cc index a755eca0c..743f1e564 100644 --- a/src/DEM/SolidBoundary/InfinitePlaneSolidBoundary.cc +++ b/src/DEM/SolidBoundary/InfinitePlaneSolidBoundary.cc @@ -36,7 +36,7 @@ distance(const Vector& position) const { template typename Dimension::Vector InfinitePlaneSolidBoundary:: -velocity(const Vector& position) const { +localVelocity(const Vector& position) const { return mVelocity; } diff --git a/src/DEM/SolidBoundary/InfinitePlaneSolidBoundary.hh b/src/DEM/SolidBoundary/InfinitePlaneSolidBoundary.hh index eff5365f0..ea7c84a92 100644 --- a/src/DEM/SolidBoundary/InfinitePlaneSolidBoundary.hh +++ b/src/DEM/SolidBoundary/InfinitePlaneSolidBoundary.hh @@ -30,7 +30,7 @@ public: ~InfinitePlaneSolidBoundary(); virtual Vector distance(const Vector& position) const override; - virtual Vector velocity(const Vector& position) const override; + virtual Vector localVelocity(const Vector& position) const override; virtual void registerState(DataBase& dataBase, State& state) override; diff --git a/src/DEM/SolidBoundary/RectangularPlaneSolidBoundary.cc b/src/DEM/SolidBoundary/RectangularPlaneSolidBoundary.cc index 57d0abc2a..9b04932de 100644 --- a/src/DEM/SolidBoundary/RectangularPlaneSolidBoundary.cc +++ b/src/DEM/SolidBoundary/RectangularPlaneSolidBoundary.cc @@ -40,7 +40,7 @@ distance(const Vector& position) const { template typename Dimension::Vector RectangularPlaneSolidBoundary:: -velocity(const Vector& position) const { +localVelocity(const Vector& position) const { return mVelocity; } diff --git a/src/DEM/SolidBoundary/RectangularPlaneSolidBoundary.hh b/src/DEM/SolidBoundary/RectangularPlaneSolidBoundary.hh index bef29204a..82d1ee5d3 100644 --- a/src/DEM/SolidBoundary/RectangularPlaneSolidBoundary.hh +++ b/src/DEM/SolidBoundary/RectangularPlaneSolidBoundary.hh @@ -33,7 +33,7 @@ public: ~RectangularPlaneSolidBoundary(); virtual Vector distance(const Vector& position) const override; - virtual Vector velocity(const Vector& position) const override; + virtual Vector localVelocity(const Vector& position) const override; virtual void registerState(DataBase& dataBase, State& state) override; diff --git a/src/DEM/SolidBoundary/SolidBoundaryBase.hh b/src/DEM/SolidBoundary/SolidBoundaryBase.hh index 849c7668d..9b350bd8c 100644 --- a/src/DEM/SolidBoundary/SolidBoundaryBase.hh +++ b/src/DEM/SolidBoundary/SolidBoundaryBase.hh @@ -46,7 +46,7 @@ public: virtual ~SolidBoundaryBase(); virtual Vector distance(const Vector& position) const = 0; - virtual Vector velocity(const Vector& position) const = 0; + virtual Vector localVelocity(const Vector& position) const = 0; virtual void registerState(DataBase& dataBase, State& state) = 0; diff --git a/src/DEM/SolidBoundary/SphereSolidBoundary.cc b/src/DEM/SolidBoundary/SphereSolidBoundary.cc index 25f8f68ec..91b6be843 100644 --- a/src/DEM/SolidBoundary/SphereSolidBoundary.cc +++ b/src/DEM/SolidBoundary/SphereSolidBoundary.cc @@ -67,7 +67,7 @@ distance(const Vector& position) const { template typename Dimension::Vector SphereSolidBoundary:: -velocity(const Vector& position) const { +localVelocity(const Vector& position) const { return mVelocity; } diff --git a/src/DEM/SolidBoundary/SphereSolidBoundary.hh b/src/DEM/SolidBoundary/SphereSolidBoundary.hh index 3d260aee7..c957053f9 100644 --- a/src/DEM/SolidBoundary/SphereSolidBoundary.hh +++ b/src/DEM/SolidBoundary/SphereSolidBoundary.hh @@ -33,7 +33,7 @@ public: ~SphereSolidBoundary(); virtual Vector distance(const Vector& position) const override; - virtual Vector velocity(const Vector& position) const override; + virtual Vector localVelocity(const Vector& position) const override; virtual void registerState(DataBase& dataBase, State& state) override; diff --git a/src/PYB11/DEM/SolidBoundaries.py b/src/PYB11/DEM/SolidBoundaries.py index 19f0da6fe..50511b524 100644 --- a/src/PYB11/DEM/SolidBoundaries.py +++ b/src/PYB11/DEM/SolidBoundaries.py @@ -53,7 +53,7 @@ def update(self, @PYB11virtual @PYB11const - def velocity(self, + def localVelocity(self, position = "const Vector&"): "velocity of bc." return "Vector" @@ -105,7 +105,7 @@ def update(self, @PYB11virtual @PYB11const - def velocity(self, + def localVelocity(self, position = "const Vector&"): "velocity of bc." return "Vector" @@ -157,7 +157,7 @@ def update(self, @PYB11virtual @PYB11const - def velocity(self, + def localVelocity(self, position = "const Vector&"): "velocity of bc." return "Vector" @@ -210,7 +210,7 @@ def update(self, @PYB11virtual @PYB11const - def velocity(self, + def localVelocity(self, position = "const Vector&"): "velocity of bc." return "Vector" @@ -265,7 +265,7 @@ def update(self, @PYB11virtual @PYB11const - def velocity(self, + def localVelocity(self, position = "const Vector&"): "velocity of bc." return "Vector" @@ -283,8 +283,8 @@ def distance(self, clipPoint = PYB11property("const Vector&", "clipPoint", "clipPoint", returnpolicy="reference_internal", doc="point on clip plane") clipAxis = PYB11property("const Vector&", "clipAxis", "clipAxis", returnpolicy="reference_internal", doc="normal in clip plane") -PYB11inject(SolidBoundaryBaseAbstractMethods, SphereSolidBoundary, virtual=True, pure_virtual=False) -PYB11inject(SolidBoundaryBaseAbstractMethods, InfinitePlaneSolidBoundary, virtual=True, pure_virtual=False) -PYB11inject(SolidBoundaryBaseAbstractMethods, RectangularPlaneSolidBoundary, virtual=True, pure_virtual=False) -PYB11inject(SolidBoundaryBaseAbstractMethods, CircularPlaneSolidBoundary, virtual=True, pure_virtual=False) -PYB11inject(SolidBoundaryBaseAbstractMethods, CylinderSolidBoundary, virtual=True, pure_virtual=False) +#PYB11inject(SolidBoundaryBaseAbstractMethods, SphereSolidBoundary, virtual=True, pure_virtual=False) +#PYB11inject(SolidBoundaryBaseAbstractMethods, InfinitePlaneSolidBoundary, virtual=True, pure_virtual=False) +#PYB11inject(SolidBoundaryBaseAbstractMethods, RectangularPlaneSolidBoundary, virtual=True, pure_virtual=False) +#PYB11inject(SolidBoundaryBaseAbstractMethods, CircularPlaneSolidBoundary, virtual=True, pure_virtual=False) +#PYB11inject(SolidBoundaryBaseAbstractMethods, CylinderSolidBoundary, virtual=True, pure_virtual=False) diff --git a/src/PYB11/DEM/SolidBoundaryBaseAbstractMethods.py b/src/PYB11/DEM/SolidBoundaryBaseAbstractMethods.py index e7278fe1a..584c314fb 100644 --- a/src/PYB11/DEM/SolidBoundaryBaseAbstractMethods.py +++ b/src/PYB11/DEM/SolidBoundaryBaseAbstractMethods.py @@ -7,7 +7,7 @@ class SolidBoundaryBaseAbstractMethods: @PYB11const - def velocity(self, + def localVelocity(self, position = "const Vector&"): "velocity of bc." return "Vector" From ba25094e19d2caf85fb5c50436e1a46c5076cb70 Mon Sep 17 00:00:00 2001 From: jmpearl Date: Mon, 15 Apr 2024 12:49:44 -0700 Subject: [PATCH 35/38] DEM BC updates RotatingDrum working --- src/DEM/CMakeLists.txt | 2 + .../ClippedSphereSolidBoundary.cc | 109 ++++++++++++++++++ .../ClippedSphereSolidBoundary.hh | 84 ++++++++++++++ .../ClippedSphereSolidBoundaryInline.hh | 87 ++++++++++++++ .../ClippedSphereSolidBoundaryInst.cc.py | 11 ++ src/DEM/SolidBoundary/SphereSolidBoundary.cc | 49 ++------ src/DEM/SolidBoundary/SphereSolidBoundary.hh | 22 ++-- .../SphereSolidBoundaryInline.hh | 38 ++---- src/PYB11/DEM/DEM_PYB11.py | 2 + src/PYB11/DEM/SolidBoundaries.py | 56 ++++++++- .../MovingSolidBoundaries/rotatingDrum-2d.py | 54 ++++----- 11 files changed, 394 insertions(+), 120 deletions(-) create mode 100644 src/DEM/SolidBoundary/ClippedSphereSolidBoundary.cc create mode 100644 src/DEM/SolidBoundary/ClippedSphereSolidBoundary.hh create mode 100644 src/DEM/SolidBoundary/ClippedSphereSolidBoundaryInline.hh create mode 100644 src/DEM/SolidBoundary/ClippedSphereSolidBoundaryInst.cc.py diff --git a/src/DEM/CMakeLists.txt b/src/DEM/CMakeLists.txt index 28606c290..acfb721eb 100644 --- a/src/DEM/CMakeLists.txt +++ b/src/DEM/CMakeLists.txt @@ -9,6 +9,7 @@ set(DEM_inst SolidBoundary/RectangularPlaneSolidBoundary SolidBoundary/CylinderSolidBoundary SolidBoundary/SphereSolidBoundary + SolidBoundary/ClippedSphereSolidBoundary ReplacePairFieldList IncrementPairFieldList ReplaceAndIncrementPairFieldList @@ -29,6 +30,7 @@ set(DEM_headers SolidBoundary/RectangularPlaneSolidBoundary.hh SolidBoundary/CylinderSolidBoundary.hh SolidBoundary/SphereSolidBoundary.hh + SolidBoundary/ClippedSphereSolidBoundary.hh ReplacePairFieldList.hh ReplaceAndIncrementPairFieldList.hh IncrementPairFieldList.hh diff --git a/src/DEM/SolidBoundary/ClippedSphereSolidBoundary.cc b/src/DEM/SolidBoundary/ClippedSphereSolidBoundary.cc new file mode 100644 index 000000000..06232761a --- /dev/null +++ b/src/DEM/SolidBoundary/ClippedSphereSolidBoundary.cc @@ -0,0 +1,109 @@ +//---------------------------------Spheral++----------------------------------// +// ClippedSphereSolidBoundary -- N-dimensional spherical solid boundary for DEM. +// +// J.M. Pearl 2023 +//----------------------------------------------------------------------------// + +#include "DataBase/DataBase.hh" +#include "DataBase/State.hh" +#include "DataBase/StateDerivatives.hh" + +#include "DEM/SolidBoundary/ClippedSphereSolidBoundary.hh" + +#include + +namespace Spheral { + +template +ClippedSphereSolidBoundary:: +ClippedSphereSolidBoundary(const Vector& center, + const Scalar radius, + const Vector& clipPoint, + const Vector& clipAxis): + SolidBoundaryBase(), + mCenter(center), + mRadius(radius), + mClipPoint(clipPoint), + mClipAxis(clipAxis), + mClipIntersectionRadius(0.0), + mVelocity(Vector::zero){ + this->setClipIntersectionRadius(); + mClipAxis = mClipAxis.unitVector(); +} + +template +ClippedSphereSolidBoundary:: +~ClippedSphereSolidBoundary(){ +} + + +template +typename Dimension::Vector +ClippedSphereSolidBoundary:: +distance(const Vector& position) const { + + // contacting sphere + const auto contactPoint = (position - mCenter).unitVector()*mRadius + mCenter; + Vector dist = position - contactPoint; + + const auto planeSignedDistance = (contactPoint - mClipPoint).dot(mClipAxis); + + // if contant pt above clip plane check for edge contact + if (planeSignedDistance > 0.0){ + + // break into perp and in-plane components + const auto q = (position-mClipPoint); + const auto qnMag = q.dot(mClipAxis); + const auto qn = qnMag * mClipAxis; + const auto qr = q - qn; + + // if outside circle enforce planar solid bc + dist = min(qr.magnitude() - mClipIntersectionRadius,0.0)*qr.unitVector() + qn; + + } + return dist; +} + +template +typename Dimension::Vector +ClippedSphereSolidBoundary:: +localVelocity(const Vector& position) const { + return mVelocity; +} + +template +void +ClippedSphereSolidBoundary:: +registerState(DataBase& dataBase, + State& state) { + + const auto boundaryKey = "ClippedSphereSolidBoundary_" + std::to_string(std::abs(this->uniqueIndex())); + const auto pointKey = boundaryKey +"_point"; + const auto clipPointKey = boundaryKey +"_clipPoint"; + const auto velocityKey = boundaryKey +"_velocity"; + + state.enrollAny(pointKey,mCenter); + state.enrollAny(clipPointKey,mClipPoint); + state.enrollAny(pointKey,mVelocity); + +} + +template +void +ClippedSphereSolidBoundary:: +update(const double multiplier, const double t, const double dt) { + mCenter += multiplier*mVelocity; + mClipPoint += multiplier*mVelocity; +} + + +template +void +ClippedSphereSolidBoundary:: +setClipIntersectionRadius() { + const auto rcMag = (mClipPoint - mCenter).dot(mClipAxis); + mClipIntersectionRadius = (rcMag < mRadius ? std::sqrt(mRadius*mRadius-rcMag*rcMag) : 0.0); + mClipPoint = rcMag * mClipAxis + mCenter; +} + +} \ No newline at end of file diff --git a/src/DEM/SolidBoundary/ClippedSphereSolidBoundary.hh b/src/DEM/SolidBoundary/ClippedSphereSolidBoundary.hh new file mode 100644 index 000000000..325afdaa9 --- /dev/null +++ b/src/DEM/SolidBoundary/ClippedSphereSolidBoundary.hh @@ -0,0 +1,84 @@ +//---------------------------------Spheral++----------------------------------// +// ClippedSphereSolidBoundary -- cylinder with finite length solid boundary for DEM +// +// J.M. Pearl 2023 +//----------------------------------------------------------------------------// + +#ifndef __Spheral_ClippedSphereSolidBoundary_hh__ +#define __Spheral_ClippedSphereSolidBoundary_hh__ + +#include "DEM/SolidBoundary/SolidBoundaryBase.hh" + +namespace Spheral { + +template class State; +template class StateDerivatives; +template class DataBase; + +template +class ClippedSphereSolidBoundary : public SolidBoundaryBase { + + typedef typename Dimension::Scalar Scalar; + typedef typename Dimension::Vector Vector; + typedef typename Dimension::Tensor Tensor; + +public: + //--------------------------- Public Interface ---------------------------// + + ClippedSphereSolidBoundary(const Vector& center, + const Scalar radius, + const Vector& clipPoint, + const Vector& clipAxis); + + ~ClippedSphereSolidBoundary(); + + virtual Vector distance(const Vector& position) const override; + virtual Vector localVelocity(const Vector& position) const override; + + virtual void registerState(DataBase& dataBase, + State& state) override; + + + virtual void update(const double multiplier, + const double time, + const double dt) override; + + const Vector& center() const; + void center(const Vector& value); + + Scalar radius() const; + void radius(Scalar value); + + const Vector& clipPoint() const; + void clipPoint(const Vector& value); + + const Vector& clipAxis() const; + void clipAxis(const Vector& value); + + const Vector& velocity() const; + void velocity(const Vector& value); + + void setClipIntersectionRadius(); +protected: + //-------------------------- Protected Interface --------------------------// + Vector mCenter; + Scalar mRadius; + Vector mClipPoint; + Vector mClipAxis; + Scalar mClipIntersectionRadius; + + Vector mVelocity; + +private: + //--------------------------- Private Interface ---------------------------// + // No default constructor, copying, or assignment. + ClippedSphereSolidBoundary(); + ClippedSphereSolidBoundary(const ClippedSphereSolidBoundary&); + ClippedSphereSolidBoundary& operator=(const ClippedSphereSolidBoundary&); +}; + +} + +#include "ClippedSphereSolidBoundaryInline.hh" + +#endif diff --git a/src/DEM/SolidBoundary/ClippedSphereSolidBoundaryInline.hh b/src/DEM/SolidBoundary/ClippedSphereSolidBoundaryInline.hh new file mode 100644 index 000000000..220b6bc8b --- /dev/null +++ b/src/DEM/SolidBoundary/ClippedSphereSolidBoundaryInline.hh @@ -0,0 +1,87 @@ +namespace Spheral { + +template +inline +const typename Dimension::Vector& +ClippedSphereSolidBoundary:: +center() const { + return mCenter; +} + +template +inline +void +ClippedSphereSolidBoundary:: +center(const typename Dimension::Vector& value) { + mCenter=value; +} + +template +inline +typename Dimension::Scalar +ClippedSphereSolidBoundary:: +radius() const { + return mRadius; +} + +template +inline +void +ClippedSphereSolidBoundary:: +radius(typename Dimension::Scalar value) { + mRadius=value; +} + + +template +inline +const typename Dimension::Vector& +ClippedSphereSolidBoundary:: +clipPoint() const { + return mClipPoint; +} + +template +inline +void +ClippedSphereSolidBoundary:: +clipPoint(const typename Dimension::Vector& value) { + mClipPoint=value; + this->setClipIntersectionRadius(); +} + + +template +inline +const typename Dimension::Vector& +ClippedSphereSolidBoundary:: +clipAxis() const { + return mClipAxis; +} + +template +inline +void +ClippedSphereSolidBoundary:: +clipAxis(const typename Dimension::Vector& value) { + mClipAxis=value.unitVector(); + this->setClipIntersectionRadius(); +} + +template +inline +const typename Dimension::Vector& +ClippedSphereSolidBoundary:: +velocity() const { + return mVelocity; +} + +template +inline +void +ClippedSphereSolidBoundary:: +velocity(const typename Dimension::Vector& value) { + mVelocity=value; +} + +} \ No newline at end of file diff --git a/src/DEM/SolidBoundary/ClippedSphereSolidBoundaryInst.cc.py b/src/DEM/SolidBoundary/ClippedSphereSolidBoundaryInst.cc.py new file mode 100644 index 000000000..c76e1cedc --- /dev/null +++ b/src/DEM/SolidBoundary/ClippedSphereSolidBoundaryInst.cc.py @@ -0,0 +1,11 @@ +text = """ +//------------------------------------------------------------------------------ +// Explict instantiation. +//------------------------------------------------------------------------------ +#include "DEM/SolidBoundary/ClippedSphereSolidBoundary.cc" +#include "Geometry/Dimension.hh" + +namespace Spheral { + template class ClippedSphereSolidBoundary< Dim< %(ndim)s > >; +} +""" diff --git a/src/DEM/SolidBoundary/SphereSolidBoundary.cc b/src/DEM/SolidBoundary/SphereSolidBoundary.cc index 91b6be843..2c7b3a94c 100644 --- a/src/DEM/SolidBoundary/SphereSolidBoundary.cc +++ b/src/DEM/SolidBoundary/SphereSolidBoundary.cc @@ -18,17 +18,12 @@ template SphereSolidBoundary:: SphereSolidBoundary(const Vector& center, const Scalar radius, - const Vector& clipPoint, - const Vector& clipAxis): + const RotationType& angularVelocity): SolidBoundaryBase(), mCenter(center), mRadius(radius), - mClipPoint(clipPoint), - mClipAxis(clipAxis), - mClipIntersectionRadius(0.0), - mVelocity(Vector::zero){ - this->setClipIntersectionRadius(); - mClipAxis = mClipAxis.unitVector(); + mVelocity(Vector::zero), + mAngularVelocity(angularVelocity){ } template @@ -41,34 +36,16 @@ template typename Dimension::Vector SphereSolidBoundary:: distance(const Vector& position) const { - - // contacting sphere - const auto contactPoint = (position - mCenter).unitVector()*mRadius + mCenter; - Vector dist = position - contactPoint; - - const auto planeSignedDistance = (contactPoint - mClipPoint).dot(mClipAxis); - - // if contant pt above clip plane check for edge contact - if (planeSignedDistance > 0.0){ - - // break into perp and in-plane components - const auto q = (position-mClipPoint); - const auto qnMag = q.dot(mClipAxis); - const auto qn = qnMag * mClipAxis; - const auto qr = q - qn; - - // if outside circle enforce planar solid bc - dist = min(qr.magnitude() - mClipIntersectionRadius,0.0)*qr.unitVector() + qn; - - } - return dist; + const auto p = position - mCenter; + return p - p.unitVector()*mRadius; } template typename Dimension::Vector SphereSolidBoundary:: localVelocity(const Vector& position) const { - return mVelocity; + const auto rVector = (position - mCenter).unitVector()*mRadius; + return mVelocity + DEMDimension::cross(mAngularVelocity,rVector); } template @@ -79,11 +56,9 @@ registerState(DataBase& dataBase, const auto boundaryKey = "SphereSolidBoundary_" + std::to_string(std::abs(this->uniqueIndex())); const auto pointKey = boundaryKey +"_point"; - const auto clipPointKey = boundaryKey +"_clipPoint"; const auto velocityKey = boundaryKey +"_velocity"; state.enrollAny(pointKey,mCenter); - state.enrollAny(clipPointKey,mClipPoint); state.enrollAny(pointKey,mVelocity); } @@ -93,17 +68,7 @@ void SphereSolidBoundary:: update(const double multiplier, const double t, const double dt) { mCenter += multiplier*mVelocity; - mClipPoint += multiplier*mVelocity; } -template -void -SphereSolidBoundary:: -setClipIntersectionRadius() { - const auto rcMag = (mClipPoint - mCenter).dot(mClipAxis); - mClipIntersectionRadius = (rcMag < mRadius ? std::sqrt(mRadius*mRadius-rcMag*rcMag) : 0.0); - mClipPoint = rcMag * mClipAxis + mCenter; -} - } \ No newline at end of file diff --git a/src/DEM/SolidBoundary/SphereSolidBoundary.hh b/src/DEM/SolidBoundary/SphereSolidBoundary.hh index c957053f9..8e44b7e79 100644 --- a/src/DEM/SolidBoundary/SphereSolidBoundary.hh +++ b/src/DEM/SolidBoundary/SphereSolidBoundary.hh @@ -7,6 +7,7 @@ #ifndef __Spheral_SphereSolidBoundary_hh__ #define __Spheral_SphereSolidBoundary_hh__ +#include "DEM/DEMDimension.hh" #include "DEM/SolidBoundary/SolidBoundaryBase.hh" namespace Spheral { @@ -21,14 +22,12 @@ class SphereSolidBoundary : public SolidBoundaryBase { typedef typename Dimension::Scalar Scalar; typedef typename Dimension::Vector Vector; typedef typename Dimension::Tensor Tensor; - + typedef typename DEMDimension::AngularVector RotationType; public: //--------------------------- Public Interface ---------------------------// - SphereSolidBoundary(const Vector& center, const Scalar radius, - const Vector& clipPoint, - const Vector& clipAxis); + const RotationType& angularVelocity); ~SphereSolidBoundary(); @@ -49,26 +48,21 @@ public: Scalar radius() const; void radius(Scalar value); - const Vector& clipPoint() const; - void clipPoint(const Vector& value); - - const Vector& clipAxis() const; - void clipAxis(const Vector& value); - const Vector& velocity() const; void velocity(const Vector& value); - void setClipIntersectionRadius(); + const RotationType& angularVelocity() const; + void angularVelocity(const RotationType& value); + protected: //-------------------------- Protected Interface --------------------------// Vector mCenter; Scalar mRadius; - Vector mClipPoint; - Vector mClipAxis; - Scalar mClipIntersectionRadius; Vector mVelocity; + RotationType mAngularVelocity; + private: //--------------------------- Private Interface ---------------------------// // No default constructor, copying, or assignment. diff --git a/src/DEM/SolidBoundary/SphereSolidBoundaryInline.hh b/src/DEM/SolidBoundary/SphereSolidBoundaryInline.hh index 58b99a7ef..f8d4b6014 100644 --- a/src/DEM/SolidBoundary/SphereSolidBoundaryInline.hh +++ b/src/DEM/SolidBoundary/SphereSolidBoundaryInline.hh @@ -32,56 +32,36 @@ radius(typename Dimension::Scalar value) { mRadius=value; } - template inline const typename Dimension::Vector& SphereSolidBoundary:: -clipPoint() const { - return mClipPoint; -} - -template -inline -void -SphereSolidBoundary:: -clipPoint(const typename Dimension::Vector& value) { - mClipPoint=value; - this->setClipIntersectionRadius(); -} - - -template -inline -const typename Dimension::Vector& -SphereSolidBoundary:: -clipAxis() const { - return mClipAxis; +velocity() const { + return mVelocity; } template inline void SphereSolidBoundary:: -clipAxis(const typename Dimension::Vector& value) { - mClipAxis=value.unitVector(); - this->setClipIntersectionRadius(); +velocity(const typename Dimension::Vector& value) { + mVelocity=value; } template inline -const typename Dimension::Vector& +const typename DEMDimension::AngularVector& SphereSolidBoundary:: -velocity() const { - return mVelocity; +angularVelocity() const { + return mAngularVelocity; } template inline void SphereSolidBoundary:: -velocity(const typename Dimension::Vector& value) { - mVelocity=value; +angularVelocity(const typename DEMDimension::AngularVector& value) { + mAngularVelocity=value; } } \ No newline at end of file diff --git a/src/PYB11/DEM/DEM_PYB11.py b/src/PYB11/DEM/DEM_PYB11.py index 7983d462f..fe7923454 100644 --- a/src/PYB11/DEM/DEM_PYB11.py +++ b/src/PYB11/DEM/DEM_PYB11.py @@ -26,6 +26,7 @@ '"DEM/SolidBoundary/CircularPlaneSolidBoundary.hh"', '"DEM/SolidBoundary/CylinderSolidBoundary.hh"', '"DEM/SolidBoundary/SphereSolidBoundary.hh"', + '"DEM/SolidBoundary/ClippedSphereSolidBoundary.hh"', '"FileIO/FileIO.hh"'] #------------------------------------------------------------------------------- @@ -46,6 +47,7 @@ CircularPlaneSolidBoundary%(ndim)id = PYB11TemplateClass(CircularPlaneSolidBoundary, template_parameters="%(Dimension)s") CylinderSolidBoundary%(ndim)id = PYB11TemplateClass(CylinderSolidBoundary, template_parameters="%(Dimension)s") SphereSolidBoundary%(ndim)id = PYB11TemplateClass(SphereSolidBoundary, template_parameters="%(Dimension)s") +ClippedSphereSolidBoundary%(ndim)id = PYB11TemplateClass(ClippedSphereSolidBoundary, template_parameters="%(Dimension)s") ''' % {"ndim" : ndim, "Dimension" : "Dim<" + str(ndim) + ">"}) diff --git a/src/PYB11/DEM/SolidBoundaries.py b/src/PYB11/DEM/SolidBoundaries.py index 50511b524..65ae21fe8 100644 --- a/src/PYB11/DEM/SolidBoundaries.py +++ b/src/PYB11/DEM/SolidBoundaries.py @@ -230,12 +230,65 @@ def distance(self, #------------------------------------------------------------------------------- -# Cylinder solid boundary. In 2d this would be two planes. +# Sphere solid boundary. In 2d this would be a circle. #------------------------------------------------------------------------------- @PYB11template("Dimension") @PYB11module("SpheralDEM") class SphereSolidBoundary(SolidBoundaryBase): + PYB11typedefs = """ + typedef typename %(Dimension)s::Scalar Scalar; + typedef typename %(Dimension)s::Vector Vector; + typedef typename DEMDimension<%(Dimension)s>::AngularVector RotationType; + """ + + def pyinit(self, + center = "const Vector&", + radius = "const Scalar", + angularVelocity = "const RotationType&"): + "solid planar boundary" + + @PYB11virtual + def registerState(self, + dataBase = "DataBase<%(Dimension)s>&", + state = "State<%(Dimension)s>&"): + "Register the state solid bc expects to use and evolve." + return "void" + + @PYB11virtual + def update(self, + multiplier = "const double", + t = "const double", + dt = "const double",): + "distance vector to bc." + return "void" + + @PYB11virtual + @PYB11const + def localVelocity(self, + position = "const Vector&"): + "velocity of bc." + return "Vector" + + @PYB11virtual + @PYB11const + def distance(self, + position = "const Vector&"): + "distance vector to bc." + return "Vector" + + velocity = PYB11property("const Vector&", "velocity", "velocity", returnpolicy="reference_internal", doc="velocity of plane") + center = PYB11property("const Vector&", "center", "center", returnpolicy="reference_internal", doc="center of sphere") + radius = PYB11property("Scalar", "radius", "radius", doc="radius of sphere") + angularVelocity = PYB11property("const RotationType&", "angularVelocity", "angularVelocity", doc="rotation about center point") + +#------------------------------------------------------------------------------- +# Sphere solid boundary intersected with an infinite plane. +#------------------------------------------------------------------------------- +@PYB11template("Dimension") +@PYB11module("SpheralDEM") +class ClippedSphereSolidBoundary(SolidBoundaryBase): + PYB11typedefs = """ typedef typename %(Dimension)s::Scalar Scalar; typedef typename %(Dimension)s::Vector Vector; @@ -283,6 +336,7 @@ def distance(self, clipPoint = PYB11property("const Vector&", "clipPoint", "clipPoint", returnpolicy="reference_internal", doc="point on clip plane") clipAxis = PYB11property("const Vector&", "clipAxis", "clipAxis", returnpolicy="reference_internal", doc="normal in clip plane") + #PYB11inject(SolidBoundaryBaseAbstractMethods, SphereSolidBoundary, virtual=True, pure_virtual=False) #PYB11inject(SolidBoundaryBaseAbstractMethods, InfinitePlaneSolidBoundary, virtual=True, pure_virtual=False) #PYB11inject(SolidBoundaryBaseAbstractMethods, RectangularPlaneSolidBoundary, virtual=True, pure_virtual=False) diff --git a/tests/functional/DEM/LinearSpringDEM/MovingSolidBoundaries/rotatingDrum-2d.py b/tests/functional/DEM/LinearSpringDEM/MovingSolidBoundaries/rotatingDrum-2d.py index 020c8456b..f6df34ad5 100644 --- a/tests/functional/DEM/LinearSpringDEM/MovingSolidBoundaries/rotatingDrum-2d.py +++ b/tests/functional/DEM/LinearSpringDEM/MovingSolidBoundaries/rotatingDrum-2d.py @@ -11,33 +11,18 @@ else: from DistributeNodes import distributeNodes2d -title("DEM Rotating Drum") - -class RotatingDrumBoundary(SphereSolidBoundary2d): - def __init__(self, - center, - radius, - omega): - SphereSolidBoundary2d.__init__(self, - center = center, - radius = radius, - clipPoint = center + Vector(0.0,radius*1.1), - clipAxis = Vector(0,1)) - return - -print("THERE") #------------------------------------------------------------------------------- # Generic problem parameters #------------------------------------------------------------------------------- -commandLine(omegaDrum = 0.1, # angular velocity of drum +commandLine(omegaDrum = 0.25, # angular velocity of drum radiusDrum = 10.0, # radius of the drum yThresh = 0.2, # level to initial fill to yClip = 0.0, # level to clip at after settling - g0 = 1.0, # grav acceleration + g0 = 5.0, # grav acceleration radius = 0.5, # particle radius normalSpringConstant=10000.0, # spring constant for LDS model - normalRestitutionCoefficient=1.00, # restitution coefficient to get damping const + normalRestitutionCoefficient=0.55, # restitution coefficient to get damping const tangentialSpringConstant=2857.0, # spring constant for LDS model tangentialRestitutionCoefficient=0.55, # restitution coefficient to get damping const dynamicFriction = 1.0, # static friction coefficient sliding @@ -45,14 +30,16 @@ def __init__(self, rollingFriction = 1.05, # static friction coefficient for rolling torsionalFriction = 1.3, # static friction coefficient for torsion cohesiveTensileStrength = 0.0, # units of pressure - shapeFactor = 0.5, # in [0,1] shape factor from Zhang 2018, 0 - no torsion or rolling + shapeFactor = 0.1, # in [0,1] shape factor from Zhang 2018, 0 - no torsion or rolling neighborSearchBuffer = 0.1, # multiplicative buffer to radius for neighbor search algo # integration IntegratorConstructor = VerletIntegrator, # Verlet one integrator to garenteee conservation stepsPerCollision = 50, # replaces CFL for DEM - goalTime = 3.0, + updateBoundaryFrequency = 10, # CAREFUL: make sure fast time stepping is off for DEM + settleTime = 5.0, # time simulated before we start spinning + goalTime = 15.0, # duration of spin cycle dt = 1e-8, dtMin = 1.0e-8, dtMax = 0.1, @@ -85,7 +72,6 @@ def __init__(self, #------------------------------------------------------------------------------- # check for bad inputs #------------------------------------------------------------------------------- -assert mpi.procs == 1 assert radius < radiusDrum/2.0 assert shapeFactor <= 1.0 and shapeFactor >= 0.0 assert dynamicFriction >= 0.0 @@ -205,23 +191,17 @@ def rejecter(x,y,m,H): torsionalFrictionCoefficient = torsionalFriction, cohesiveTensileStrength =cohesiveTensileStrength, shapeFactor = shapeFactor, - stepsPerCollision = stepsPerCollision) + stepsPerCollision = stepsPerCollision, + enableFastTimeStepping = False) packages = [dem] -# set yup a circle (clip outside circle to deactivate clipping (i know not the best)) -solidWall = RotatingDrumBoundary(center = Vector(0.0, 0.0), +solidWall = SphereSolidBoundary(center = Vector(0.0, 0.0), radius = radiusDrum, - omega = 0.0) - -# solidWall = SphereSolidBoundary(center = Vector(0.0, 0.0), -# radius = radiusDrum, -# clipPoint = Vector(0.0,radiusDrum*1.1), -# clipAxis = Vector(0,1)) + angularVelocity = 0.0) -print("here") dem.appendSolidBoundary(solidWall) -print("HERE") + #------------------------------------------------------------------------------- # Gravity: DEM #------------------------------------------------------------------------------- @@ -251,7 +231,7 @@ def rejecter(x,y,m,H): integrator.domainDecompositionIndependent = domainIndependent integrator.verbose = dtverbose integrator.rigorousBoundaries = rigorousBoundaries - +integrator.updateBoundaryFrequency = 10 integrator.cullGhostNodes = False output("integrator") @@ -289,5 +269,11 @@ def rejecter(x,y,m,H): if not steps is None: control.step(steps) else: - control.advance(goalTime, maxSteps) + control.advance(settleTime, maxSteps) +solidWall.angularVelocity = omegaDrum + +if not steps is None: + control.step(steps) +else: + control.advance(settleTime+goalTime, maxSteps) From ab6f92891165dc6438fba08bebe5cc3bc145d254 Mon Sep 17 00:00:00 2001 From: jmpearl Date: Sun, 12 May 2024 08:09:46 -0700 Subject: [PATCH 36/38] whoopsie --- .../DEM/LinearSpringDEM/MovingSolidBoundaries/rotatingDrum-2d.py | 1 - 1 file changed, 1 deletion(-) diff --git a/tests/functional/DEM/LinearSpringDEM/MovingSolidBoundaries/rotatingDrum-2d.py b/tests/functional/DEM/LinearSpringDEM/MovingSolidBoundaries/rotatingDrum-2d.py index f6df34ad5..3503c6352 100644 --- a/tests/functional/DEM/LinearSpringDEM/MovingSolidBoundaries/rotatingDrum-2d.py +++ b/tests/functional/DEM/LinearSpringDEM/MovingSolidBoundaries/rotatingDrum-2d.py @@ -247,7 +247,6 @@ def rejecter(x,y,m,H): #------------------------------------------------------------------------------- # Make the problem controller. #------------------------------------------------------------------------------- -print("POOPSIE") from SpheralPointmeshSiloDump import dumpPhysicsState control = SpheralController(integrator, WT, iterateInitialH = False, From c517fad8a75ffef721cb6f6c42940cad2effc043 Mon Sep 17 00:00:00 2001 From: jmpearl Date: Mon, 13 May 2024 03:21:43 -0700 Subject: [PATCH 37/38] DEM restart for solid bcs --- .../CircularPlaneSolidBoundary.cc | 29 ++++++++++++++++ .../CircularPlaneSolidBoundary.hh | 4 +++ .../ClippedSphereSolidBoundary.cc | 33 +++++++++++++++++++ .../ClippedSphereSolidBoundary.hh | 4 +++ .../SolidBoundary/CylinderSolidBoundary.cc | 31 +++++++++++++++++ .../SolidBoundary/CylinderSolidBoundary.hh | 4 +++ .../InfinitePlaneSolidBoundary.cc | 26 +++++++++++++++ .../InfinitePlaneSolidBoundary.hh | 4 +++ .../RectangularPlaneSolidBoundary.cc | 29 ++++++++++++++++ .../RectangularPlaneSolidBoundary.hh | 4 +++ src/DEM/SolidBoundary/SolidBoundaryBase.cc | 5 ++- src/DEM/SolidBoundary/SolidBoundaryBase.hh | 10 ++++++ src/DEM/SolidBoundary/SphereSolidBoundary.cc | 27 +++++++++++++++ src/DEM/SolidBoundary/SphereSolidBoundary.hh | 4 +++ 14 files changed, 213 insertions(+), 1 deletion(-) diff --git a/src/DEM/SolidBoundary/CircularPlaneSolidBoundary.cc b/src/DEM/SolidBoundary/CircularPlaneSolidBoundary.cc index 41441195f..efda0cf5e 100644 --- a/src/DEM/SolidBoundary/CircularPlaneSolidBoundary.cc +++ b/src/DEM/SolidBoundary/CircularPlaneSolidBoundary.cc @@ -5,12 +5,17 @@ // J.M. Pearl 2023 //----------------------------------------------------------------------------// +#include "FileIO/FileIO.hh" + #include "DataBase/DataBase.hh" #include "DataBase/State.hh" #include "DataBase/StateDerivatives.hh" #include "DEM/SolidBoundary/CircularPlaneSolidBoundary.hh" +#include +using std::string; + namespace Spheral { template @@ -71,4 +76,28 @@ update(const double multiplier, const double t, const double dt) { } +//------------------------------------------------------------------------------ +// Restart +//------------------------------------------------------------------------------ +template +void +CircularPlaneSolidBoundary:: +dumpState(FileIO& file, const string& pathName) const { + file.write(mPoint, pathName + "/point"); + file.write(mNormal, pathName + "/normal"); + file.write(mExtent, pathName + "/extent"); + file.write(mVelocity, pathName + "/velocity"); +} + + +template +void +CircularPlaneSolidBoundary:: +restoreState(const FileIO& file, const string& pathName) { + file.read(mPoint, pathName + "/point"); + file.read(mNormal, pathName + "/normal"); + file.read(mExtent, pathName + "/extent"); + file.read(mVelocity, pathName + "/velocity"); +} + } \ No newline at end of file diff --git a/src/DEM/SolidBoundary/CircularPlaneSolidBoundary.hh b/src/DEM/SolidBoundary/CircularPlaneSolidBoundary.hh index 15ee17b02..173321759 100644 --- a/src/DEM/SolidBoundary/CircularPlaneSolidBoundary.hh +++ b/src/DEM/SolidBoundary/CircularPlaneSolidBoundary.hh @@ -53,6 +53,10 @@ public: const Vector& velocity() const; void velocity(const Vector& value); + virtual std::string label() const { return "CircularPlaneSolidBoundary" ; } + virtual void dumpState(FileIO& file, const std::string& pathName) const override; + virtual void restoreState(const FileIO& file, const std::string& pathName) override; + protected: //-------------------------- Protected Interface --------------------------// Vector mPoint; diff --git a/src/DEM/SolidBoundary/ClippedSphereSolidBoundary.cc b/src/DEM/SolidBoundary/ClippedSphereSolidBoundary.cc index 06232761a..0be718247 100644 --- a/src/DEM/SolidBoundary/ClippedSphereSolidBoundary.cc +++ b/src/DEM/SolidBoundary/ClippedSphereSolidBoundary.cc @@ -4,6 +4,8 @@ // J.M. Pearl 2023 //----------------------------------------------------------------------------// +#include "FileIO/FileIO.hh" + #include "DataBase/DataBase.hh" #include "DataBase/State.hh" #include "DataBase/StateDerivatives.hh" @@ -12,6 +14,9 @@ #include +#include +using std::string; + namespace Spheral { template @@ -106,4 +111,32 @@ setClipIntersectionRadius() { mClipPoint = rcMag * mClipAxis + mCenter; } +//------------------------------------------------------------------------------ +// Restart +//------------------------------------------------------------------------------ +template +void +ClippedSphereSolidBoundary:: +dumpState(FileIO& file, const string& pathName) const { + file.write(mCenter, pathName + "/center"); + file.write(mRadius, pathName + "/radius"); + file.write(mClipPoint, pathName + "/clipPoint"); + file.write(mClipAxis, pathName + "/clipAxis"); + file.write(mClipIntersectionRadius, pathName + "/clipIntersectionRadius"); + file.write(mVelocity, pathName + "/velocity"); +} + + +template +void +ClippedSphereSolidBoundary:: +restoreState(const FileIO& file, const string& pathName) { + file.read(mCenter, pathName + "/center"); + file.read(mRadius, pathName + "/radius"); + file.read(mClipPoint, pathName + "/clipPoint"); + file.read(mClipAxis, pathName + "/clipAxis"); + file.read(mClipIntersectionRadius, pathName + "/clipIntersectionRadius"); + file.read(mVelocity, pathName + "/velocity"); +} + } \ No newline at end of file diff --git a/src/DEM/SolidBoundary/ClippedSphereSolidBoundary.hh b/src/DEM/SolidBoundary/ClippedSphereSolidBoundary.hh index 325afdaa9..2e479685e 100644 --- a/src/DEM/SolidBoundary/ClippedSphereSolidBoundary.hh +++ b/src/DEM/SolidBoundary/ClippedSphereSolidBoundary.hh @@ -59,6 +59,10 @@ public: void velocity(const Vector& value); void setClipIntersectionRadius(); + + virtual std::string label() const { return "ClippedSphereSolidBoundary" ; } + virtual void dumpState(FileIO& file, const std::string& pathName) const override; + virtual void restoreState(const FileIO& file, const std::string& pathName) override; protected: //-------------------------- Protected Interface --------------------------// Vector mCenter; diff --git a/src/DEM/SolidBoundary/CylinderSolidBoundary.cc b/src/DEM/SolidBoundary/CylinderSolidBoundary.cc index fffbae75f..8280ef1c4 100644 --- a/src/DEM/SolidBoundary/CylinderSolidBoundary.cc +++ b/src/DEM/SolidBoundary/CylinderSolidBoundary.cc @@ -4,12 +4,17 @@ // J.M. Pearl 2023 //----------------------------------------------------------------------------// +#include "FileIO/FileIO.hh" + #include "DataBase/DataBase.hh" #include "DataBase/State.hh" #include "DataBase/StateDerivatives.hh" #include "DEM/SolidBoundary/CylinderSolidBoundary.hh" +#include +using std::string; + namespace Spheral { template @@ -71,5 +76,31 @@ update(const double multiplier, const double t, const double dt) { mPoint += multiplier*mVelocity; } +//------------------------------------------------------------------------------ +// Restart +//------------------------------------------------------------------------------ +template +void +CylinderSolidBoundary:: +dumpState(FileIO& file, const string& pathName) const { + file.write(mPoint, pathName + "/point"); + file.write(mAxis, pathName + "/axis"); + file.write(mRadius, pathName + "/radius"); + file.write(mLength, pathName + "/length"); + file.write(mVelocity, pathName + "/velocity"); +} + + +template +void +CylinderSolidBoundary:: +restoreState(const FileIO& file, const string& pathName) { + file.read(mPoint, pathName + "/point"); + file.read(mAxis, pathName + "/axis"); + file.read(mRadius, pathName + "/radius"); + file.read(mLength, pathName + "/length"); + file.read(mVelocity, pathName + "/velocity"); +} + } \ No newline at end of file diff --git a/src/DEM/SolidBoundary/CylinderSolidBoundary.hh b/src/DEM/SolidBoundary/CylinderSolidBoundary.hh index b53e046a2..7da17e8d8 100644 --- a/src/DEM/SolidBoundary/CylinderSolidBoundary.hh +++ b/src/DEM/SolidBoundary/CylinderSolidBoundary.hh @@ -57,6 +57,10 @@ public: const Vector& velocity() const; void velocity(const Vector& value); + virtual std::string label() const { return "CylinderSolidBoundary" ; } + virtual void dumpState(FileIO& file, const std::string& pathName) const override; + virtual void restoreState(const FileIO& file, const std::string& pathName) override; + protected: //-------------------------- Protected Interface --------------------------// Vector mPoint; diff --git a/src/DEM/SolidBoundary/InfinitePlaneSolidBoundary.cc b/src/DEM/SolidBoundary/InfinitePlaneSolidBoundary.cc index 743f1e564..eb29ba5c5 100644 --- a/src/DEM/SolidBoundary/InfinitePlaneSolidBoundary.cc +++ b/src/DEM/SolidBoundary/InfinitePlaneSolidBoundary.cc @@ -4,12 +4,17 @@ // J.M. Pearl 2023 //----------------------------------------------------------------------------// +#include "FileIO/FileIO.hh" + #include "DataBase/DataBase.hh" #include "DataBase/State.hh" #include "DataBase/StateDerivatives.hh" #include "DEM/SolidBoundary/InfinitePlaneSolidBoundary.hh" +#include +using std::string; + namespace Spheral { template @@ -61,5 +66,26 @@ update(const double multiplier, const double t, const double dt) { mPoint += multiplier*mVelocity; } +//------------------------------------------------------------------------------ +// Restart +//------------------------------------------------------------------------------ +template +void +InfinitePlaneSolidBoundary:: +dumpState(FileIO& file, const string& pathName) const { + file.write(mPoint, pathName + "/point"); + file.write(mNormal, pathName + "/normal"); + file.write(mVelocity, pathName + "/velocity"); +} + + +template +void +InfinitePlaneSolidBoundary:: +restoreState(const FileIO& file, const string& pathName) { + file.read(mPoint, pathName + "/point"); + file.read(mNormal, pathName + "/normal"); + file.read(mVelocity, pathName + "/velocity"); +} } \ No newline at end of file diff --git a/src/DEM/SolidBoundary/InfinitePlaneSolidBoundary.hh b/src/DEM/SolidBoundary/InfinitePlaneSolidBoundary.hh index ea7c84a92..2c80a5be6 100644 --- a/src/DEM/SolidBoundary/InfinitePlaneSolidBoundary.hh +++ b/src/DEM/SolidBoundary/InfinitePlaneSolidBoundary.hh @@ -48,6 +48,10 @@ public: const Vector& velocity() const; void velocity(const Vector& value); + virtual std::string label() const { return "InfinitePlaneSolidBoundary" ; } + virtual void dumpState(FileIO& file, const std::string& pathName) const override; + virtual void restoreState(const FileIO& file, const std::string& pathName) override; + protected: //-------------------------- Protected Interface --------------------------// Vector mPoint; diff --git a/src/DEM/SolidBoundary/RectangularPlaneSolidBoundary.cc b/src/DEM/SolidBoundary/RectangularPlaneSolidBoundary.cc index 9b04932de..9e38fc6e2 100644 --- a/src/DEM/SolidBoundary/RectangularPlaneSolidBoundary.cc +++ b/src/DEM/SolidBoundary/RectangularPlaneSolidBoundary.cc @@ -5,11 +5,16 @@ // J.M. Pearl 2023 //----------------------------------------------------------------------------// +#include "FileIO/FileIO.hh" + #include "DataBase/DataBase.hh" #include "DataBase/State.hh" #include "DataBase/StateDerivatives.hh" #include "DEM/SolidBoundary/RectangularPlaneSolidBoundary.hh" +#include +using std::string; + namespace Spheral { template @@ -62,5 +67,29 @@ update(const double multiplier, const double t, const double dt) { mPoint += multiplier*mVelocity; } +//------------------------------------------------------------------------------ +// Restart +//------------------------------------------------------------------------------ +template +void +RectangularPlaneSolidBoundary:: +dumpState(FileIO& file, const string& pathName) const { + file.write(mPoint, pathName + "/point"); + file.write(mBasis, pathName + "/basis"); + file.write(mExtent, pathName + "/extent"); + file.write(mVelocity, pathName + "/velocity"); +} + + +template +void +RectangularPlaneSolidBoundary:: +restoreState(const FileIO& file, const string& pathName) { + file.read(mPoint, pathName + "/point"); + file.read(mBasis, pathName + "/basis"); + file.read(mExtent, pathName + "/extent"); + file.read(mVelocity, pathName + "/velocity"); +} + } \ No newline at end of file diff --git a/src/DEM/SolidBoundary/RectangularPlaneSolidBoundary.hh b/src/DEM/SolidBoundary/RectangularPlaneSolidBoundary.hh index 82d1ee5d3..62c75e917 100644 --- a/src/DEM/SolidBoundary/RectangularPlaneSolidBoundary.hh +++ b/src/DEM/SolidBoundary/RectangularPlaneSolidBoundary.hh @@ -54,6 +54,10 @@ public: const Vector& velocity() const; void velocity(const Vector& value); + virtual std::string label() const { return "RectangularPlaneSolidBoundary" ; } + virtual void dumpState(FileIO& file, const std::string& pathName) const override; + virtual void restoreState(const FileIO& file, const std::string& pathName) override; + protected: //-------------------------- Protected Interface --------------------------// Vector mPoint; diff --git a/src/DEM/SolidBoundary/SolidBoundaryBase.cc b/src/DEM/SolidBoundary/SolidBoundaryBase.cc index 5bc9f98eb..5e919452a 100644 --- a/src/DEM/SolidBoundary/SolidBoundaryBase.cc +++ b/src/DEM/SolidBoundary/SolidBoundaryBase.cc @@ -20,7 +20,10 @@ namespace Spheral { template SolidBoundaryBase:: -SolidBoundaryBase(){} +SolidBoundaryBase(): + mUniqueIndex(-1), + mRestart(registerWithRestart(*this)){ +} template SolidBoundaryBase:: ~SolidBoundaryBase(){} diff --git a/src/DEM/SolidBoundary/SolidBoundaryBase.hh b/src/DEM/SolidBoundary/SolidBoundaryBase.hh index 9b350bd8c..ef8b89f5c 100644 --- a/src/DEM/SolidBoundary/SolidBoundaryBase.hh +++ b/src/DEM/SolidBoundary/SolidBoundaryBase.hh @@ -24,6 +24,8 @@ #ifndef __Spheral_SolidBoundaryBase_hh__ #define __Spheral_SolidBoundaryBase_hh__ +#include "DataOutput/registerWithRestart.hh" + #include namespace Spheral { @@ -31,6 +33,7 @@ namespace Spheral { template class State; template class StateDerivatives; template class DataBase; +class FileIO; template class SolidBoundaryBase { @@ -58,10 +61,17 @@ public: void uniqueIndex(int uId); int uniqueIndex() const; + // restartability will default to no-op + virtual std::string label() const { return "SolidBoundaryBase" ; } + virtual void dumpState(FileIO& file, const std::string& pathName) const {}; + virtual void restoreState(const FileIO& file, const std::string& pathName) {}; + private: //--------------------------- Public Interface ---------------------------// int mUniqueIndex; +RestartRegistrationType mRestart; + }; } diff --git a/src/DEM/SolidBoundary/SphereSolidBoundary.cc b/src/DEM/SolidBoundary/SphereSolidBoundary.cc index 2c7b3a94c..3a600a1a1 100644 --- a/src/DEM/SolidBoundary/SphereSolidBoundary.cc +++ b/src/DEM/SolidBoundary/SphereSolidBoundary.cc @@ -4,6 +4,8 @@ // J.M. Pearl 2023 //----------------------------------------------------------------------------// +#include "FileIO/FileIO.hh" + #include "DataBase/DataBase.hh" #include "DataBase/State.hh" #include "DataBase/StateDerivatives.hh" @@ -11,6 +13,8 @@ #include "DEM/SolidBoundary/SphereSolidBoundary.hh" #include +#include +using std::string; namespace Spheral { @@ -70,5 +74,28 @@ update(const double multiplier, const double t, const double dt) { mCenter += multiplier*mVelocity; } +//------------------------------------------------------------------------------ +// Restart +//------------------------------------------------------------------------------ +template +void +SphereSolidBoundary:: +dumpState(FileIO& file, const string& pathName) const { + file.write(mAngularVelocity, pathName + "/omega"); + file.write(mCenter, pathName + "/center"); + file.write(mRadius, pathName + "/radius"); + file.write(mVelocity, pathName + "/velocity"); +} + + +template +void +SphereSolidBoundary:: +restoreState(const FileIO& file, const string& pathName) { + file.read(mAngularVelocity, pathName + "/omega"); + file.read(mCenter, pathName + "/center"); + file.read(mRadius, pathName + "/radius"); + file.read(mVelocity, pathName + "/velocity"); +} } \ No newline at end of file diff --git a/src/DEM/SolidBoundary/SphereSolidBoundary.hh b/src/DEM/SolidBoundary/SphereSolidBoundary.hh index 8e44b7e79..c6a932e48 100644 --- a/src/DEM/SolidBoundary/SphereSolidBoundary.hh +++ b/src/DEM/SolidBoundary/SphereSolidBoundary.hh @@ -54,6 +54,10 @@ public: const RotationType& angularVelocity() const; void angularVelocity(const RotationType& value); + virtual std::string label() const { return "SphereSolidBoundary" ; } + virtual void dumpState(FileIO& file, const std::string& pathName) const override; + virtual void restoreState(const FileIO& file, const std::string& pathName) override; + protected: //-------------------------- Protected Interface --------------------------// Vector mCenter; From df5d638c2b08300c5609ad23017f726136befad5 Mon Sep 17 00:00:00 2001 From: jmpearl Date: Mon, 13 May 2024 15:49:48 -0700 Subject: [PATCH 38/38] updating release notes --- RELEASE_NOTES.md | 1 + 1 file changed, 1 insertion(+) diff --git a/RELEASE_NOTES.md b/RELEASE_NOTES.md index 9ba54c5d5..fe35d2250 100644 --- a/RELEASE_NOTES.md +++ b/RELEASE_NOTES.md @@ -38,6 +38,7 @@ Notable changes include: * DEM * fixed bug in solid boundary unique indices that causes particle sticking * fixed bug in solid boundary update policies + * fixed solid boundary restartability for moving bcs Version v2024.01.00 -- Release date 2024-01-19 ==============================================