From d2d613c8b4d134e692bfc0a3f918dce3e1dbe98f Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Fri, 27 Jun 2025 22:20:32 -0400 Subject: [PATCH] bugfix: incorrect application of abs() to doubles @jtclemm this looks like a real bug. Can you please check how much of an impact this change has? --- src/GRANULAR/gran_sub_mod_normal.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/GRANULAR/gran_sub_mod_normal.cpp b/src/GRANULAR/gran_sub_mod_normal.cpp index 53a9a3952c..54e8f1eaa6 100644 --- a/src/GRANULAR/gran_sub_mod_normal.cpp +++ b/src/GRANULAR/gran_sub_mod_normal.cpp @@ -851,14 +851,14 @@ double GranSubModNormalMDR::calculate_forces() for (int lv1 = 0; lv1 < MDR_MAX_IT; ++lv1) { fa_tmp = deltae1D - A * 0.5 + A * sqrt(Bsq * 0.25 - square(aAdh_tmp)) * Binv; fa = fa_tmp + sqrt(MY_2PI * aAdh_tmp * gamma * Eeffinv); - if (abs(fa) < MDR_EPSILON1) { + if (fabs(fa) < MDR_EPSILON1) { break; } dfda = -aAdh_tmp * A / (B * sqrt(-square(aAdh_tmp) + Bsq * 0.25)); dfda += gamma * SQRTHALFPI / sqrt(aAdh_tmp * gamma * Eeff); aAdh_tmp = aAdh_tmp - fa / dfda; fa2 = fa_tmp + sqrt(MY_2PI * aAdh_tmp * gamma * Eeffinv); - if (abs(fa - fa2) < MDR_EPSILON2) { + if (fabs(fa - fa2) < MDR_EPSILON2) { break; } if (lv1 == MDR_MAX_IT - 1) {