From 2454691ebf5580e5273d58ae4dd277e1e18ed624 Mon Sep 17 00:00:00 2001 From: jmpearl Date: Tue, 8 Jul 2025 07:30:50 -0700 Subject: [PATCH 1/2] LSDEM: consistency improvements with Zhang 2018. mu_r and mu_t were off by factor of 2, other edits are aesthetic --- src/DEM/LinearSpringDEM.cc | 36 ++++++++++++++++++------------------ 1 file changed, 18 insertions(+), 18 deletions(-) diff --git a/src/DEM/LinearSpringDEM.cc b/src/DEM/LinearSpringDEM.cc index 906a6e23b..edaf4821a 100644 --- a/src/DEM/LinearSpringDEM.cc +++ b/src/DEM/LinearSpringDEM.cc @@ -453,15 +453,15 @@ evaluateDerivatives(const typename Dimension::Scalar /*time*/, // spring constants const auto kn = this->normalSpringConstant(); // normal const auto ks = this->tangentialSpringConstant(); // sliding - const auto kt = 0.50 * ks * shapeFactor2; - const auto kr = 0.25 * kn * shapeFactor2; + const auto kt = 2.0 * ks * shapeFactor2; + const auto kr = kn * shapeFactor2; const auto invKs = 1.0/max(ks,tiny); const auto invKt = 1.0/max(kt,tiny); const auto invKr = 1.0/max(kr,tiny); - const auto normalDampingTerms = 2.0*kn/(1.0+mNormalBeta*mNormalBeta); - const auto tangentialDampingTerms = 2.0*ks/(1.0+mTangentialBeta*mTangentialBeta); + const auto normalDampingTerms = 4.0*kn/(1.0+mNormalBeta*mNormalBeta); + const auto tangentialDampingTerms = 4.0*ks/(1.0+mTangentialBeta*mTangentialBeta); // The connectivity. const auto& nodeLists = dataBase.DEMNodeListPtrs(); @@ -636,14 +636,14 @@ evaluateDerivatives(const typename Dimension::Scalar /*time*/, const auto lj = rijMag-li; // effective quantities (mass, reduced radius) respectively, mij->mi for mi=mj - const auto mij = 2.0*(mi*mj)/(mi+mj); - const auto lij = 2.0*(li*lj)/(li+lj); + const auto mij = mi*mj/(mi+mj); + const auto lij = li*lj/(li+lj); // damping constants -- ct and cr derived quantities ala Zhang 2017 const auto Cn = std::sqrt(mij*normalDampingTerms); const auto Cs = std::sqrt(mij*tangentialDampingTerms); - const auto Ct = 0.50 * Cs * shapeFactor2; - const auto Cr = 0.25 * Cn * shapeFactor2; + const auto Ct = 2.0 * Cs * shapeFactor2; + const auto Cr = Cn * shapeFactor2; // our velocities const Vector vroti = -DEMDimension::cross(omegai,rhatij); @@ -658,9 +658,9 @@ evaluateDerivatives(const typename Dimension::Scalar /*time*/, // normal forces //------------------------------------------------------------ - const Vector fn = (kn*delta - Cn*vn)*rhatij; // normal spring - const Vector fc = Cc*shapeFactor2*lij*lij*rhatij; // normal cohesion - const Scalar fnMag = fn.magnitude(); // magnitude of normal spring force + const Vector fn = (kn*delta - Cn*vn)*rhatij; // normal spring + const Vector fc = 4.0*Cc*shapeFactor2*lij*lij*rhatij; // normal cohesion + const Scalar fnMag = fn.magnitude(); // magnitude of normal spring force // sliding //------------------------------------------------------------ @@ -771,14 +771,14 @@ evaluateDerivatives(const typename Dimension::Scalar /*time*/, const auto li = rib.magnitude(); // effective quantities - const auto mib = 2*mi; - const auto lib = 2*li; + const auto mib = mi; + const auto lib = li; // damping coefficients const auto Cn = std::sqrt(mib*normalDampingTerms); const auto Cs = std::sqrt(mib*tangentialDampingTerms); - const auto Ct = 0.50 * Cs * shapeFactor2; - const auto Cr = 0.25 * Cn * shapeFactor2; + const auto Ct = 2.0 * Cs * shapeFactor2; + const auto Cr = Cn * shapeFactor2; // velocities const Vector vroti = -DEMDimension::cross(omegai,rhatib); @@ -792,9 +792,9 @@ evaluateDerivatives(const typename Dimension::Scalar /*time*/, // normal forces //------------------------------------------------------------ - const Vector fn = (kn*delta - Cn*vn)*rhatib; // normal spring - const Vector fc = Cc*shapeFactor2*lib*lib*rhatib; // normal cohesion - const Scalar fnMag = fn.magnitude(); // magnitude of normal spring force + const Vector fn = (kn*delta - Cn*vn)*rhatib; // normal spring + const Vector fc = 4.0*Cc*shapeFactor2*lib*lib*rhatib; // normal cohesion + const Scalar fnMag = fn.magnitude(); // magnitude of normal spring force // sliding //------------------------------------------------------------ From e51cb9706276486240bc7c55598df8a12060fe68 Mon Sep 17 00:00:00 2001 From: jmpearl Date: Tue, 16 Sep 2025 07:12:40 -0700 Subject: [PATCH 2/2] updating DEM bugfix in realease readme --- RELEASE_NOTES.md | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/RELEASE_NOTES.md b/RELEASE_NOTES.md index 900c4faeb..327d89bdf 100644 --- a/RELEASE_NOTES.md +++ b/RELEASE_NOTES.md @@ -21,7 +21,10 @@ Notable changes include: * Stack allocation of tensor data; Static casting for CRTP implementation. * GeomTensor & GeomSymmetricTensor have been refactored for use on the GPU. * New Logging utility for runtime debug messages. - + + * Bug fixes + * corrected rolling and torsional coefficient in DEM which were 2x the expected value + * Build changes / improvements: * Changed `int` to `size_t` for Field and FieldList. * A python virtual environment is installed in the spheral build dir, removing the