diff --git a/RELEASE_NOTES.md b/RELEASE_NOTES.md index f853ab37c..fe35d2250 100644 --- a/RELEASE_NOTES.md +++ b/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 ============================================== 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/DEM.py b/src/DEM/DEM.py index f24ae5f9b..98ba8ac91 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 = 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) diff --git a/src/DEM/DEMBase.cc b/src/DEM/DEMBase.cc index 6057c68f1..03867b6e2 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 @@ -81,6 +79,7 @@ DEMBase(const DataBase& dataBase, const Vector& xmax): Physics(), mDataBase(dataBase), + mNewSolidBoundaryIndex(0), mSolidBoundaries(), mCycle(0), mContactRemovalFrequency((int)stepsPerCollision), @@ -316,9 +315,12 @@ 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::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& 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 //------------------------------------------------------------------------------ @@ -421,8 +438,9 @@ initialize(const Scalar time, const DataBase& dataBase, State& state, StateDerivatives& derivs){ +TIME_BEGIN("DEMinitialize"); - +TIME_END("DEMinitialize"); } @@ -576,6 +594,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(); @@ -616,6 +636,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 @@ -628,6 +649,9 @@ template void DEMBase:: 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 @@ -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 @@ -752,6 +786,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(); @@ -846,7 +882,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); @@ -861,7 +897,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 @@ -870,6 +906,7 @@ updateContactMap(const DataBase& dataBase){ } // loop nodes } // loop nodelists } // loop solid boundaries + TIME_END("DEMupdateContactMap"); } // method //------------------------------------------------------------------------------ @@ -879,7 +916,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(); @@ -934,5 +971,6 @@ identifyInactiveContacts(const DataBase& dataBase){ } // loop particle bc contacts } // omp parallel region +TIME_END("DEMidentifyInactiveContacts"); } // class } // namespace diff --git a/src/DEM/DEMBase.hh b/src/DEM/DEMBase.hh index 5b158f0d3..2053eb918 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 @@ -192,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; @@ -216,13 +225,14 @@ protected: const DataBase& mDataBase; + int mNewSolidBoundaryIndex; 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 +260,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. diff --git a/src/DEM/DEMBaseInline.hh b/src/DEM/DEMBaseInline.hh index dffe855e2..f716956dc 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,16 @@ 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()) mSolidBoundaries.erase(bcPtr); +} + + template inline unsigned int @@ -401,12 +420,4 @@ numParticleBoundaryContacts() const { return (mContactStorageIndices.size()-this->numParticleParticleContacts()); } -template -inline -int -DEMBase:: -getSolidBoundaryUniqueIndex(const int x) const { - return -x-1; -} - } diff --git a/src/DEM/DEMFieldNames.cc b/src/DEM/DEMFieldNames.cc index ce4603288..b0420b608 100644 --- a/src/DEM/DEMFieldNames.cc +++ b/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"; \ No newline at end of file +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 e3109e33c..3358a9925 100644 --- a/src/DEM/DEMFieldNames.hh +++ b/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; }; } diff --git a/src/DEM/LinearSpringDEM.cc b/src/DEM/LinearSpringDEM.cc index 05cef1c15..56e3e2f92 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" @@ -33,6 +33,8 @@ #include "DEM/ContactStorageLocation.hh" #include "DEM/SolidBoundary/SolidBoundaryBase.hh" +#include "Utilities/Timer.hh" + #ifdef _OPENMP #include "omp.h" #endif @@ -78,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), @@ -93,8 +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))), - mMomentOfInertia(FieldStorageType::CopyFields) { + 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))); } //------------------------------------------------------------------------------ @@ -112,18 +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(); @@ -215,7 +260,7 @@ dt(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; @@ -251,6 +296,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(); @@ -271,7 +320,6 @@ dt(const DataBase& dataBase, return result; } - //------------------------------------------------------------------------------ // method that fires once on startup //------------------------------------------------------------------------------ @@ -279,8 +327,10 @@ template void LinearSpringDEM:: initializeProblemStartup(DataBase& dataBase){ + TIME_BEGIN("LinearSpringDEMinitializeProblemStartup"); DEMBase::initializeProblemStartup(dataBase); this->setMomentOfInertia(); + TIME_END("LinearSpringDEMinitializeProblemStartup"); } //------------------------------------------------------------------------------ @@ -291,10 +341,40 @@ void 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 //------------------------------------------------------------------------------ @@ -306,7 +386,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. @@ -316,8 +396,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(); @@ -332,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(); @@ -378,10 +458,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()); @@ -412,6 +494,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 @@ -449,7 +532,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); @@ -473,10 +556,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); @@ -489,6 +574,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; @@ -516,8 +602,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 @@ -529,76 +615,39 @@ 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; + // 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; @@ -663,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(); @@ -675,12 +724,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; @@ -704,62 +748,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 @@ -768,7 +771,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; @@ -796,6 +799,8 @@ evaluateDerivatives(const typename Dimension::Scalar /*time*/, DxDt(nodeListi,i) = veli; } // loop nodes } // loop nodelists + + TIME_END("LinearSpringDEMevaluateDerivatives"); } // method @@ -875,6 +880,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"); } //------------------------------------------------------------------------------ @@ -886,6 +893,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 diff --git a/src/DEM/LinearSpringDEM.hh b/src/DEM/LinearSpringDEM.hh index b932dc7d5..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); @@ -65,19 +66,59 @@ 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, 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; + + // 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, + 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 + bool enableFastTimeStepping() const; + void enableFastTimeStepping(bool x); + Scalar normalSpringConstant() const; void normalSpringConstant(Scalar x); @@ -114,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(); @@ -123,6 +167,8 @@ public: // get methods for class FieldLists const FieldList& momentOfInertia() const; + const FieldList& maximumOverlap() const; + const FieldList& newMaximumOverlap() const; //**************************************************************************** @@ -133,6 +179,7 @@ public: //**************************************************************************** private: //--------------------------- Private Interface ---------------------------// + Scalar mEnableFastTimeStepping; Scalar mNormalSpringConstant; Scalar mNormalRestitutionCoefficient; Scalar mTangentialSpringConstant; @@ -146,9 +193,12 @@ private: Scalar mNormalBeta; Scalar mTangentialBeta; + Scalar mCollisionDuration; // 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..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 @@ -246,6 +280,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 //------------------------------------------------------------------------------ @@ -257,4 +354,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/DEM/SolidBoundary/CircularPlaneSolidBoundary.cc b/src/DEM/SolidBoundary/CircularPlaneSolidBoundary.cc index ec2d559dc..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 @@ -45,10 +50,24 @@ distance(const Vector& position) const { template typename Dimension::Vector CircularPlaneSolidBoundary:: -velocity(const Vector& position) const { +localVelocity(const Vector& position) const { return mVelocity; } +template +void +CircularPlaneSolidBoundary:: +registerState(DataBase& dataBase, + 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"; + state.enrollAny(pointKey,mPoint); + state.enrollAny(pointKey,mVelocity); + state.enrollAny(pointKey,mNormal); +} + template void CircularPlaneSolidBoundary:: @@ -57,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 f1e2c0937..173321759 100644 --- a/src/DEM/SolidBoundary/CircularPlaneSolidBoundary.hh +++ b/src/DEM/SolidBoundary/CircularPlaneSolidBoundary.hh @@ -32,7 +32,10 @@ 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; virtual void update(const double multiplier, const double time, @@ -50,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 new file mode 100644 index 000000000..0be718247 --- /dev/null +++ b/src/DEM/SolidBoundary/ClippedSphereSolidBoundary.cc @@ -0,0 +1,142 @@ +//---------------------------------Spheral++----------------------------------// +// ClippedSphereSolidBoundary -- N-dimensional spherical solid boundary for DEM. +// +// J.M. Pearl 2023 +//----------------------------------------------------------------------------// + +#include "FileIO/FileIO.hh" + +#include "DataBase/DataBase.hh" +#include "DataBase/State.hh" +#include "DataBase/StateDerivatives.hh" + +#include "DEM/SolidBoundary/ClippedSphereSolidBoundary.hh" + +#include + +#include +using std::string; + +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; +} + +//------------------------------------------------------------------------------ +// 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 new file mode 100644 index 000000000..2e479685e --- /dev/null +++ b/src/DEM/SolidBoundary/ClippedSphereSolidBoundary.hh @@ -0,0 +1,88 @@ +//---------------------------------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(); + + 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; + 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/CylinderSolidBoundary.cc b/src/DEM/SolidBoundary/CylinderSolidBoundary.cc index 082258876..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 @@ -31,7 +36,6 @@ CylinderSolidBoundary:: ~CylinderSolidBoundary(){ } - template typename Dimension::Vector CylinderSolidBoundary:: @@ -47,10 +51,24 @@ distance(const Vector& position) const { template typename Dimension::Vector CylinderSolidBoundary:: -velocity(const Vector& position) const { +localVelocity(const Vector& position) const { return mVelocity; } +template +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 void CylinderSolidBoundary:: @@ -58,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 da95e2cfc..7da17e8d8 100644 --- a/src/DEM/SolidBoundary/CylinderSolidBoundary.hh +++ b/src/DEM/SolidBoundary/CylinderSolidBoundary.hh @@ -33,7 +33,10 @@ 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; virtual void update(const double multiplier, const double time, @@ -54,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 a15cb7a44..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 @@ -36,10 +41,24 @@ distance(const Vector& position) const { template typename Dimension::Vector InfinitePlaneSolidBoundary:: -velocity(const Vector& position) const { +localVelocity(const Vector& position) const { return mVelocity; } +template +void +InfinitePlaneSolidBoundary:: +registerState(DataBase& dataBase, + 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"; + state.enrollAny(pointKey,mPoint); + state.enrollAny(velocityKey,mVelocity); + state.enrollAny(normalKey,mNormal); +} + template void InfinitePlaneSolidBoundary:: @@ -47,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 b80539818..2c80a5be6 100644 --- a/src/DEM/SolidBoundary/InfinitePlaneSolidBoundary.hh +++ b/src/DEM/SolidBoundary/InfinitePlaneSolidBoundary.hh @@ -30,8 +30,11 @@ 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; + virtual void update(const double multiplier, const double time, const double dt) override; @@ -45,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 5290b957e..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 @@ -40,10 +45,21 @@ distance(const Vector& position) const { template typename Dimension::Vector RectangularPlaneSolidBoundary:: -velocity(const Vector& position) const { +localVelocity(const Vector& position) const { return mVelocity; } +template +void +RectangularPlaneSolidBoundary:: +registerState(DataBase& dataBase, + 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); + state.enrollAny(velocityKey,mVelocity); +} template void RectangularPlaneSolidBoundary:: @@ -51,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 3442a6791..62c75e917 100644 --- a/src/DEM/SolidBoundary/RectangularPlaneSolidBoundary.hh +++ b/src/DEM/SolidBoundary/RectangularPlaneSolidBoundary.hh @@ -33,7 +33,10 @@ 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; virtual void update(const double multiplier, const double time, @@ -51,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 2bd843c4f..ef8b89f5c 100644 --- a/src/DEM/SolidBoundary/SolidBoundaryBase.hh +++ b/src/DEM/SolidBoundary/SolidBoundaryBase.hh @@ -10,18 +10,30 @@ // 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 //----------------------------------------------------------------------------// #ifndef __Spheral_SolidBoundaryBase_hh__ #define __Spheral_SolidBoundaryBase_hh__ +#include "DataOutput/registerWithRestart.hh" + +#include + namespace Spheral { template class State; template class StateDerivatives; template class DataBase; +class FileIO; template class SolidBoundaryBase { @@ -37,13 +49,32 @@ 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; virtual void update(const double multiplier, const double t, const double dt) = 0; -}; + + 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; + +}; } +#include "SolidBoundaryBaseInline.hh" + #endif diff --git a/src/DEM/SolidBoundary/SolidBoundaryBaseInline.hh b/src/DEM/SolidBoundary/SolidBoundaryBaseInline.hh new file mode 100644 index 000000000..8a4e4a12c --- /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(int uId){ + mUniqueIndex = uId; +} +template +inline +int +SolidBoundaryBase:: +uniqueIndex() const { + return mUniqueIndex; +} +} + diff --git a/src/DEM/SolidBoundary/SphereSolidBoundary.cc b/src/DEM/SolidBoundary/SphereSolidBoundary.cc index b0781cff4..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 { @@ -18,17 +22,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,51 +40,62 @@ template typename Dimension::Vector SphereSolidBoundary:: distance(const Vector& position) const { + const auto p = position - mCenter; + return p - p.unitVector()*mRadius; +} - // contacting sphere - const auto contactPoint = (position - mCenter).unitVector()*mRadius + mCenter; - Vector dist = position - contactPoint; - - const auto planeSignedDistance = (contactPoint - mClipPoint).dot(mClipAxis); +template +typename Dimension::Vector +SphereSolidBoundary:: +localVelocity(const Vector& position) const { + const auto rVector = (position - mCenter).unitVector()*mRadius; + return mVelocity + DEMDimension::cross(mAngularVelocity,rVector); +} - // if contant pt above clip plane check for edge contact - if (planeSignedDistance > 0.0){ +template +void +SphereSolidBoundary:: +registerState(DataBase& dataBase, + State& state) { - // 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; + const auto boundaryKey = "SphereSolidBoundary_" + std::to_string(std::abs(this->uniqueIndex())); + const auto pointKey = boundaryKey +"_point"; + const auto velocityKey = boundaryKey +"_velocity"; - // if outside circle enforce planar solid bc - dist = min(qr.magnitude() - mClipIntersectionRadius,0.0)*qr.unitVector() + qn; + state.enrollAny(pointKey,mCenter); + state.enrollAny(pointKey,mVelocity); - } - return dist; } template -typename Dimension::Vector +void SphereSolidBoundary:: -velocity(const Vector& position) const { - return mVelocity; +update(const double multiplier, const double t, const double dt) { + mCenter += multiplier*mVelocity; } +//------------------------------------------------------------------------------ +// Restart +//------------------------------------------------------------------------------ template void SphereSolidBoundary:: -update(const double multiplier, const double t, const double dt) { - mCenter += multiplier*mVelocity; +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:: -setClipIntersectionRadius() { - const auto rcMag = (mClipPoint - mCenter).dot(mClipAxis); - mClipIntersectionRadius = (rcMag < mRadius ? std::sqrt(mRadius*mRadius-rcMag*rcMag) : 0.0); - mClipPoint = rcMag * mClipAxis + mCenter; +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 0aa78f63b..c6a932e48 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,19 +22,21 @@ 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(); 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; + virtual void update(const double multiplier, const double time, @@ -45,26 +48,25 @@ 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); + + 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; 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/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; +} + +} + diff --git a/src/DataBase/StateBase.cc b/src/DataBase/StateBase.cc index c11aa57dd..4ae5e7c99 100644 --- a/src/DataBase/StateBase.cc +++ b/src/DataBase/StateBase.cc @@ -140,7 +140,25 @@ 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&) { + 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"; + } + } } } } @@ -391,12 +409,23 @@ 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&) { + 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"); + } + } } } } - // Copy the connectivity (by reference). This thing is too // big to carry around separate copies! if (rhs.mConnectivityMapPtr != NULL) { @@ -446,8 +475,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"); + } } } } diff --git a/src/PYB11/DEM/DEMBase.py b/src/PYB11/DEM/DEMBase.py index 0627717eb..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,20 +103,14 @@ 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?" 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 @@ -145,6 +138,8 @@ def getSolidBoundaryUniqueIndex(x="const int"): 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") 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/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 diff --git a/src/PYB11/DEM/SolidBoundaries.py b/src/PYB11/DEM/SolidBoundaries.py index 8923e5d28..65ae21fe8 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,10 +31,18 @@ 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" + @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", @@ -43,7 +53,7 @@ def update(self, @PYB11virtual @PYB11const - def velocity(self, + def localVelocity(self, position = "const Vector&"): "velocity of bc." return "Vector" @@ -72,11 +82,19 @@ 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" + @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", @@ -87,7 +105,7 @@ def update(self, @PYB11virtual @PYB11const - def velocity(self, + def localVelocity(self, position = "const Vector&"): "velocity of bc." return "Vector" @@ -116,11 +134,19 @@ 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" + @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", @@ -131,7 +157,7 @@ def update(self, @PYB11virtual @PYB11const - def velocity(self, + def localVelocity(self, position = "const Vector&"): "velocity of bc." return "Vector" @@ -160,12 +186,20 @@ 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"): "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", @@ -176,7 +210,7 @@ def update(self, @PYB11virtual @PYB11const - def velocity(self, + def localVelocity(self, position = "const Vector&"): "velocity of bc." return "Vector" @@ -196,23 +230,84 @@ 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; """ - def pyinit(center = "const Vector&", + def pyinit(self, + center = "const Vector&", radius = "const Scalar", clipPoint = "const Vector&", clipAxis = "const Vector&"): "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", @@ -223,7 +318,7 @@ def update(self, @PYB11virtual @PYB11const - def velocity(self, + def localVelocity(self, position = "const Vector&"): "velocity of bc." return "Vector" @@ -239,4 +334,11 @@ 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) diff --git a/src/PYB11/DEM/SolidBoundaryBaseAbstractMethods.py b/src/PYB11/DEM/SolidBoundaryBaseAbstractMethods.py index d370930de..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" @@ -18,6 +18,12 @@ def distance(self, "distance vector to bc." return "Vector" + def registerState(self, + dataBase = "DataBase<%(Dimension)s>&", + state = "State<%(Dimension)s>&"): + "Register the state solid bc expects to use and evolve." + return "void" + def update(self, multiplier = "const double", t = "const double", diff --git a/tests/functional/DEM/DEMTests.ats b/tests/functional/DEM/DEMTests.ats index d2204d41e..bcb3ff032 100644 --- a/tests/functional/DEM/DEMTests.ats +++ b/tests/functional/DEM/DEMTests.ats @@ -10,6 +10,10 @@ 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") + +# test moving solid boundaries +source("LinearSpringDEM/MovingSolidBoundaries/singleParticleMovingBoundaryCollision-3d.py") \ No newline at end of file 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 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..3503c6352 --- /dev/null +++ b/tests/functional/DEM/LinearSpringDEM/MovingSolidBoundaries/rotatingDrum-2d.py @@ -0,0 +1,278 @@ +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 + +#------------------------------------------------------------------------------- +# Generic problem parameters +#------------------------------------------------------------------------------- +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 = 5.0, # grav acceleration + + radius = 0.5, # particle radius + normalSpringConstant=10000.0, # spring constant for LDS model + 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 + 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 + + 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 + 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, + 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 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, + enableFastTimeStepping = False) + +packages = [dem] + +solidWall = SphereSolidBoundary(center = Vector(0.0, 0.0), + radius = radiusDrum, + angularVelocity = 0.0) + +dem.appendSolidBoundary(solidWall) + +#------------------------------------------------------------------------------- +# 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.updateBoundaryFrequency = 10 +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(settleTime, maxSteps) + +solidWall.angularVelocity = omegaDrum + +if not steps is None: + control.step(steps) +else: + control.advance(settleTime+goalTime, maxSteps) diff --git a/tests/functional/DEM/LinearSpringDEM/MovingSolidBoundaries/singleParticleMovingBoundaryCollision-3d.py b/tests/functional/DEM/LinearSpringDEM/MovingSolidBoundaries/singleParticleMovingBoundaryCollision-3d.py new file mode 100644 index 000000000..22a4e0f39 --- /dev/null +++ b/tests/functional/DEM/LinearSpringDEM/MovingSolidBoundaries/singleParticleMovingBoundaryCollision-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 diff --git a/tests/functional/DEM/LinearSpringDEM/SolidBoundaryCondition/pureTorsion-3d.py b/tests/functional/DEM/LinearSpringDEM/SolidBoundaryCondition/pureTorsion-3d.py new file mode 100644 index 000000000..0302d1c63 --- /dev/null +++ b/tests/functional/DEM/LinearSpringDEM/SolidBoundaryCondition/pureTorsion-3d.py @@ -0,0 +1,356 @@ +#ATS:DEM3dPT = test( SELF, "--clearDirectories True --checkRestitutionCoefficient True --checkNaturalFrequency True", label="DEM pure torsion test -- 3-D (serial)") + +import os, sys, shutil, mpi +from math import * +from Spheral3d import * +from SpheralTestUtilities import * +from findLastRestart import * +from GenerateNodeDistribution3d import * +from GenerateDEMfromSPHGenerator import GenerateDEMfromSPHGenerator3d + +import numpy as np + +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 +#------------------------------------------------------------------------------- +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=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 = 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 = 50, + 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-PureTorsionTest-3d", + + # ats + checkRestitutionCoefficient = False, + threshRestitutionCoefficient = 0.05, + checkNaturalFrequency = False, + threshNaturalFrequency = 0.1, + ) + +#------------------------------------------------------------------------------- +# file things +#------------------------------------------------------------------------------- +testName = "DEM-PureTorsionTest-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) + + # transforms particle properties from generator + 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.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) + +#------------------------------------------------------------------------------- +# 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 + +#------------------------------------------------------------------------------- +# 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 resitution coefficient +#------------------------------------------------------------------------------- +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. +#------------------------------------------------------------------------------- +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") + +#------------------------------------------------------------------------------- +# Advance to the end time. +#------------------------------------------------------------------------------- + +if not steps is None: + control.step(steps) +else: + + # 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 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/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..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 ) @@ -342,15 +343,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