Merge pull request #269 from LLNL/bugfix/demSolidBCs

DEM Solid Boundary Bug Fix

RELEASE_NOTES.md

@@ -23,7 +23,10 @@ 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:
     * Improved the target export functionality.
 
@@ -32,6 +35,10 @@ Notable changes include:
     * Updating header lists for including Spheral modules in external projects.
     * Adding effective viscous pressure back to FSISPH.
     * 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 
+      * fixed solid boundary restartability for moving bcs
 
 Version v2024.01.00 -- Release date 2024-01-19
 ==============================================

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

src/DEM/DEM.py

@@ -18,6 +18,7 @@ def LinearSpringDEM(dataBase,
                     cohesiveTensileStrength = 0.0,
                     shapeFactor = 0.0,
                     stepsPerCollision = 25,
+                    enableFastTimeStepping = True,
                     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 = True,
         xmin = (-1e100, -1e100, -1e100),
         xmax = ( 1e100,  1e100,  1e100)):
     return LinearSpringDEM(dataBase,
@@ -97,6 +101,7 @@ def DEM(dataBase,
                            cohesiveTensileStrength,
                            shapeFactor,
                            stepsPerCollision,
+                           enableFastTimeStepping,
                            xmin,
                            xmax)
 

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 <iostream>
 #include <stdexcept>
 #include <limits.h>
@@ -81,6 +79,7 @@ DEMBase(const DataBase<Dimension>& dataBase,
         const Vector& xmax):
   Physics<Dimension>(),
   mDataBase(dataBase),
+  mNewSolidBoundaryIndex(0),
   mSolidBoundaries(),
   mCycle(0),
   mContactRemovalFrequency((int)stepsPerCollision),
@@ -316,9 +315,12 @@ registerState(DataBase<Dimension>& dataBase,
   auto rollingDisplacementPolicy = make_policy<ReplaceAndIncrementPairFieldList<Dimension,std::vector<Vector>>>();
   auto torsionalDisplacementPolicy = make_policy<ReplaceAndIncrementPairFieldList<Dimension,std::vector<Scalar>>>();
 
+  // solid boundary conditions w/ properties that need to be integrated
   auto boundaryPolicy = make_policy<DEMBoundaryPolicy<Dimension>>(mSolidBoundaries);
-
-  state.enroll(DEMFieldNames::solidBoundaries,boundaryPolicy);
+  state.enroll(DEMFieldNames::solidBoundaryPolicy,boundaryPolicy);
+
+  for (auto ibc = 0u; ibc < this->numSolidBoundaries(); ++ibc){
+    mSolidBoundaries[ibc]->registerState(dataBase,state);}
 
   state.enroll(mTimeStepMask);
   state.enroll(mass);
@@ -410,6 +412,21 @@ preStepInitialize(const DataBase<Dimension>& dataBase,
   TIME_END("DEMpreStepInitialize");
 }
 
+//------------------------------------------------------------------------------
+// This method is called once at the beginning of a timestep, after all state registration.
+//------------------------------------------------------------------------------
+// template<typename Dimension>
+// void
+// DEMBase<Dimension>::
+// finalize(const DataBase<Dimension>& dataBase, 
+//                   State<Dimension>& state,
+//                   StateDerivatives<Dimension>& derivatives) {
+//   TIME_BEGIN("DEMfinalize");
+
+
+//   TIME_END("DEMfinalize");
+// }
+
 //------------------------------------------------------------------------------
 // Call before deriv evaluation
 //------------------------------------------------------------------------------
@@ -421,8 +438,9 @@ initialize(const Scalar  time,
            const DataBase<Dimension>& dataBase,
                  State<Dimension>& state,
                  StateDerivatives<Dimension>& derivs){
+TIME_BEGIN("DEMinitialize");
 
-
+TIME_END("DEMinitialize");
 }
 
 
@@ -576,6 +594,8 @@ void
 DEMBase<Dimension>::
 initializeOverlap(const DataBase<Dimension>& dataBase, const int startingCompositeParticleIndex){
 
+  TIME_BEGIN("DEMinitializeOverlap");
+
   const auto& connectivityMap = dataBase.connectivityMap();
   const auto& pairs = connectivityMap.nodePairList();
   const auto  numPairs = pairs.size();
@@ -616,6 +636,7 @@ initializeOverlap(const DataBase<Dimension>& dataBase, const int startingComposi
       mEquilibriumOverlap(storeNodeList,storeNode)[storeContact] = delta0;
     }
   }
+  TIME_END("DEMinitializeOverlap");
 }
 //------------------------------------------------------------------------------
 // Redistribution methods -- before we redistribute, we are going to make sure
@@ -628,6 +649,9 @@ template<typename Dimension>
 void
 DEMBase<Dimension>::
 initializeBeforeRedistribution(){
+
+  TIME_BEGIN("DEMinitializeBeforeRedistribution");
+
   this->prepNeighborIndicesForRedistribution();
 
   this->prepPairFieldListForRedistribution(mShearDisplacement);
@@ -641,6 +665,9 @@ initializeBeforeRedistribution(){
   this->prepPairFieldListForRedistribution(mNewRollingDisplacement);
   this->prepPairFieldListForRedistribution(mDDtTorsionalDisplacement);
   this->prepPairFieldListForRedistribution(mNewTorsionalDisplacement);
+
+  TIME_END("DEMinitializeBeforeRedistribution");
+
 }
 
 template<typename Dimension>
@@ -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<typename Dimension>
 template<typename Value>
@@ -752,6 +786,8 @@ void
 DEMBase<Dimension>::
 updateContactMap(const DataBase<Dimension>& dataBase){
 
+  TIME_BEGIN("DEMupdateContactMap");
+
   const auto& uniqueIndex = dataBase.DEMUniqueIndex();
   const auto& radius = dataBase.DEMParticleRadius();
   const auto& position = dataBase.DEMPosition();
@@ -846,7 +882,7 @@ updateContactMap(const DataBase<Dimension>& 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);
@@ -861,7 +897,7 @@ updateContactMap(const DataBase<Dimension>& 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
@@ -870,6 +906,7 @@ updateContactMap(const DataBase<Dimension>& dataBase){
       }   // loop nodes
     }     // loop nodelists
   }       // loop solid boundaries
+  TIME_END("DEMupdateContactMap");
 }         // method
 
 //------------------------------------------------------------------------------
@@ -879,7 +916,7 @@ template<typename Dimension>
 void
 DEMBase<Dimension>::
 identifyInactiveContacts(const DataBase<Dimension>& dataBase){
-
+  TIME_BEGIN("DEMidentifyInactiveContacts");
   const auto bufferDistance = dataBase.maxNeighborSearchBuffer();
   const auto  numNodeLists = dataBase.numNodeLists();
   const auto& nodeListPtrs = dataBase.DEMNodeListPtrs();
@@ -934,5 +971,6 @@ identifyInactiveContacts(const DataBase<Dimension>& dataBase){
 
   } // loop particle bc contacts
 }   // omp parallel region
+TIME_END("DEMidentifyInactiveContacts");
 }   // class
 }   // namespace

src/DEM/DEMBase.hh

@@ -87,6 +87,14 @@ public:
                            const DataBase<Dimension>& dataBase,
                            const State<Dimension>& state,
                            StateDerivatives<Dimension>& derivs) const override;
+
+  // hook after the intergrator step
+  // virtual
+  // void finalize(const Scalar time,
+  //               const Scalar dt,
+  //               const DataBase<Dimension>& dataBase,
+  //               const State<Dimension>& state,
+  //                     StateDerivatives<Dimension>& derivs) const override;
 
   // Apply boundary conditions to the physics specific fields.
   virtual
@@ -192,12 +200,13 @@ public:
                              const Vector vrotj) const;
 
   // Solid Bounderies 
+  int newSolidBoundaryIndex() const;
   void appendSolidBoundary(SolidBoundaryBase<Dimension>& boundary);
   void clearSolidBoundaries();
+  void removeSolidBoundary(const SolidBoundaryBase<Dimension>& boundary);
   bool haveSolidBoundary(const SolidBoundaryBase<Dimension>& boundary) const;
   unsigned int numSolidBoundaries() const;
   const std::vector<SolidBoundaryBase<Dimension>*>& solidBoundaryConditions() const;
-  int getSolidBoundaryUniqueIndex(const int x) const;
 
   // counts
   unsigned int numParticleParticleContacts() const;
@@ -216,13 +225,14 @@ protected:
 
   const DataBase<Dimension>& mDataBase;
 
+  int mNewSolidBoundaryIndex;
   std::vector<SolidBoundaryBase<Dimension>*> 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 +260,7 @@ protected:
   FieldList<Dimension,std::vector<Scalar>> mDDtTorsionalDisplacement;  // derivative to evolve frictional spring displacement
   FieldList<Dimension,std::vector<Scalar>> mNewTorsionalDisplacement;  // handles rotation of frictional spring and reset on slip
 
+  // map to storage location from connectivityMap to pairwise fieldlists
   std::vector<ContactIndex> mContactStorageIndices;
 
   // The restart registration.

src/DEM/DEMBaseInline.hh

@@ -321,8 +321,6 @@ rollingMoment(const Dim<3>::Vector rhatij,
   return DEMDimension<Dim<3>>::cross((vroti + vrotj),rhatij).unitVector();
 }
 
-
-
 //------------------------------------------------------------------------------
 // Add a Boundary condition to the end of the current boundary list.
 //------------------------------------------------------------------------------
@@ -331,7 +329,9 @@ inline
 void
 DEMBase<Dimension>::
 appendSolidBoundary(SolidBoundaryBase<Dimension>& boundary) {
-    mSolidBoundaries.push_back(&boundary);
+  mNewSolidBoundaryIndex -= 1;
+  boundary.uniqueIndex(mNewSolidBoundaryIndex);
+  mSolidBoundaries.push_back(&boundary);
 }
 
 //------------------------------------------------------------------------------
@@ -345,6 +345,15 @@ clearSolidBoundaries() {
   mSolidBoundaries = std::vector<SolidBoundaryBase<Dimension>*>();
 }
 
+
+template<typename Dimension>
+inline
+int
+DEMBase<Dimension>::
+newSolidBoundaryIndex() const {
+  return mNewSolidBoundaryIndex;
+}
+
 //------------------------------------------------------------------------------
 // Test if the given Boundary condition is listed in the physics package.
 //------------------------------------------------------------------------------
@@ -356,6 +365,16 @@ haveSolidBoundary(const SolidBoundaryBase<Dimension>& boundary) const {
   return std::count(mSolidBoundaries.begin(), mSolidBoundaries.end(), &boundary) > 0;
 }
 
+template<typename Dimension>
+inline
+void
+DEMBase<Dimension>::
+removeSolidBoundary(const SolidBoundaryBase<Dimension>& boundary) {
+  const auto bcPtr = std::find(mSolidBoundaries.begin(),mSolidBoundaries.end(),&boundary);
+  if (bcPtr != mSolidBoundaries.end()) mSolidBoundaries.erase(bcPtr);
+}
+
+
 template<typename Dimension>
 inline
 unsigned int
@@ -401,12 +420,4 @@ numParticleBoundaryContacts() const {
   return (mContactStorageIndices.size()-this->numParticleParticleContacts());
 }
 
-template<typename Dimension>
-inline
-int
-DEMBase<Dimension>::
-getSolidBoundaryUniqueIndex(const int x) const {
-  return -x-1;
-}
-
 }

src/DEM/DEMFieldNames.cc

@@ -18,4 +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::solidBoundaries = "solid boundaries";
+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";

src/DEM/DEMFieldNames.hh

@@ -23,7 +23,9 @@ 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;
 };
 
 }