Skip to content

LinearSpringDEM bugfix #384

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Draft
wants to merge 4 commits into
base: develop
Choose a base branch
from
Draft

LinearSpringDEM bugfix #384

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
2454691
LSDEM: consistency improvements with Zhang 2018. mu_r and mu_t were o…
jmpearl Jul 8, 2025
eefb630
Merge branch 'develop' into bugfix/DEM
jmpearl Jul 15, 2025
5d0482b
Merge branch 'develop' into bugfix/DEM
jmpearl Sep 8, 2025
e51cb97
updating DEM bugfix in realease readme
jmpearl Sep 16, 2025
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
5 changes: 4 additions & 1 deletion RELEASE_NOTES.md
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
36 changes: 18 additions & 18 deletions src/DEM/LinearSpringDEM.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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();
Expand Down Expand Up @@ -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<Dimension>::cross(omegai,rhatij);
Expand All @@ -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
//------------------------------------------------------------
Expand Down Expand Up @@ -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<Dimension>::cross(omegai,rhatib);
Expand All @@ -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
//------------------------------------------------------------
Expand Down