From ece25a95ad213ee49f765bf26d8a4cb6562fbacf Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Sat, 5 Mar 2022 13:46:56 -0500 Subject: [PATCH] stability improvements: recompute bin width for perfect match, don't test for floating point identity --- .../compute_pressure_cartesian.cpp | 22 +++++++--- .../compute_pressure_spherical.cpp | 44 ++++++++++++------- 2 files changed, 45 insertions(+), 21 deletions(-) diff --git a/src/EXTRA-COMPUTE/compute_pressure_cartesian.cpp b/src/EXTRA-COMPUTE/compute_pressure_cartesian.cpp index e86e9de9c3..a18c16eb12 100644 --- a/src/EXTRA-COMPUTE/compute_pressure_cartesian.cpp +++ b/src/EXTRA-COMPUTE/compute_pressure_cartesian.cpp @@ -15,6 +15,7 @@ #include "atom.h" #include "citeme.h" +#include "comm.h" #include "domain.h" #include "error.h" #include "force.h" @@ -30,6 +31,7 @@ using namespace LAMMPS_NS; +#define SMALL 1.0e-10 /*----------------------------------------------------------------------------------- Contributing author: Olav Galteland (Norwegian University of Science and Technology) olav.galteland@ntnu.no @@ -82,7 +84,12 @@ ComputePressureCartesian::ComputePressureCartesian(LAMMPS *lmp, int narg, char * bin_width2 = 0.0; nbins1 = (int) ((domain->boxhi[dir1] - domain->boxlo[dir1]) / bin_width1); nbins2 = 1; - invV = bin_width1; + // adjust bin width if not a perfect match + invV = (domain->boxhi[dir1] - domain->boxlo[dir1]) / nbins1; + if ((fabs(invV - bin_width1) > SMALL) && (comm->me == 0)) + utils::logmesg(lmp, "Adjusting first bin width for compute {} from {:.6f} to {:.6f}\n", style, + bin_width1, invV); + bin_width1 = invV; if (bin_width1 <= 0.0) error->all(FLERR, "Illegal compute pressure/cartesian command. Bin width must be > 0"); @@ -101,6 +108,12 @@ ComputePressureCartesian::ComputePressureCartesian(LAMMPS *lmp, int narg, char * bin_width2 = utils::numeric(FLERR, arg[6], false, lmp); nbins2 = (int) ((domain->boxhi[dir2] - domain->boxlo[dir2]) / bin_width2); + double tmp_binwidth = (domain->boxhi[dir2] - domain->boxlo[dir2]) / nbins2; + if ((fabs(tmp_binwidth - bin_width2) > SMALL) && (comm->me == 0)) + utils::logmesg(lmp, "Adjusting second bin width for compute {} from {:.6f} to {:.6f}\n", + style, bin_width2, tmp_binwidth); + bin_width2 = tmp_binwidth; + invV *= bin_width2; if (bin_width2 <= 0.0) @@ -392,7 +405,7 @@ void ComputePressureCartesian::compute_pressure_1d(double fpair, double xi, doub bin = bin_s; // Integrate from bin_s to bin_e with step bin_step. - while (bin != bin_limit) { + while (bin < bin_limit) { // Calculating exit and entry point (xa, xb). Checking if inside current bin. if (bin == bin_s) { @@ -450,7 +463,6 @@ void ComputePressureCartesian::compute_pressure_2d(double fpair, double xi, doub double la = 0.0, lb = 0.0, l_sum = 0.0; double rij[3] = {delx, dely, delz}; double l1 = 0.0, l2, rij1, rij2; - double SMALL = 10e-5; rij1 = rij[dir1]; rij2 = rij[dir2]; @@ -458,7 +470,7 @@ void ComputePressureCartesian::compute_pressure_2d(double fpair, double xi, doub next_bin2 = (int) floor(yi / bin_width2); // Integrating along line - while (lb != 1.0) { + while (lb < 1.0) { bin1 = next_bin1; bin2 = next_bin2; @@ -501,7 +513,7 @@ void ComputePressureCartesian::compute_pressure_2d(double fpair, double xi, doub else if (bin2 >= nbins2) bin2 = nbins2 - 1; - if (bin1 + bin2 * nbins1 >= nbins1 * nbins2) error->all(FLERR, "Bin outside"); + if (bin1 + bin2 * nbins1 > nbins1 * nbins2) error->all(FLERR, "Bin outside: lb={:.16g}", lb); tpcxx[bin1 + bin2 * nbins1] += (fpair * delx * delx * (lb - la)); tpcyy[bin1 + bin2 * nbins1] += (fpair * dely * dely * (lb - la)); diff --git a/src/EXTRA-COMPUTE/compute_pressure_spherical.cpp b/src/EXTRA-COMPUTE/compute_pressure_spherical.cpp index a261c10270..86afd8f000 100644 --- a/src/EXTRA-COMPUTE/compute_pressure_spherical.cpp +++ b/src/EXTRA-COMPUTE/compute_pressure_spherical.cpp @@ -15,6 +15,7 @@ #include "atom.h" #include "citeme.h" +#include "comm.h" #include "domain.h" #include "error.h" #include "force.h" @@ -33,6 +34,8 @@ using namespace LAMMPS_NS; using namespace MathConst; +#define SMALL 1.0e-10 + /*----------------------------------------------------------------------------------- Contributing author: Olav Galteland (Norwegian University of Science and Technology) olav.galteland@ntnu.no @@ -71,6 +74,11 @@ ComputePressureSpherical::ComputePressureSpherical(LAMMPS *lmp, int narg, char * bin_width = utils::numeric(FLERR, arg[6], false, lmp); Rmax = utils::numeric(FLERR, arg[7], false, lmp); nbins = (int) (Rmax / bin_width) + 1; + double tmp_width = Rmax / nbins; + if ((fabs(bin_width - tmp_width) > SMALL) && (comm->me == 0)) + utils::logmesg(lmp, "Adjusting bin width for compute {} from {:.6f} to {:.6f}\n", style, + bin_width, tmp_width); + bin_width = tmp_width; if (bin_width <= 0.0) error->all(FLERR, "Illegal compute pressure/spherical command. Bin width must be > 0"); @@ -245,7 +253,6 @@ void ComputePressureSpherical::compute_array() double qi[3], l1, l2, l3, l4, R1, R2, Fa, Fb, l_sum; double rij, f, ririj, sqr, la, lb, sql0, lambda0; - double SMALL = 10e-12; double rsqxy, ririjxy, sqrixy, sqlxy0, A, B, C; int end_bin; @@ -328,24 +335,29 @@ void ComputePressureSpherical::compute_array() while (lb < 1.0) { l1 = lb + SMALL; bin = (int) floor(sqrt(rsq * l1 * l1 + 2.0 * ririj * l1 + sqr) / bin_width); - R1 = (bin) *bin_width; + R1 = bin * bin_width; R2 = (bin + 1) * bin_width; - l1 = lambda0 + sqrt((R1 * R1 - sql0) / rsq); - l2 = lambda0 - sqrt((R1 * R1 - sql0) / rsq); - l3 = lambda0 + sqrt((R2 * R2 - sql0) / rsq); - l4 = lambda0 - sqrt((R2 * R2 - sql0) / rsq); - - if (l4 >= 0.0 && l4 <= 1.0 && l4 > la) - lb = l4; - else if (l3 >= 0.0 && l3 <= 1.0 && l3 > la) - lb = l3; - else if (l2 >= 0.0 && l2 <= 1.0 && l2 > la) - lb = l2; - else if (l1 >= 0.0 && l1 <= 1.0 && l1 > la) - lb = l1; - else + // we must not take the square root of a negative number or divide by zero. + if ((R1*R1 < sql0) || (R2*R2 < sql0) || (rsq <= 0.0)) { lb = 1.0; + } else { + l1 = lambda0 + sqrt((R1 * R1 - sql0) / rsq); + l2 = lambda0 - sqrt((R1 * R1 - sql0) / rsq); + l3 = lambda0 + sqrt((R2 * R2 - sql0) / rsq); + l4 = lambda0 - sqrt((R2 * R2 - sql0) / rsq); + + if (l4 >= 0.0 && l4 <= 1.0 && l4 > la) + lb = l4; + else if (l3 >= 0.0 && l3 <= 1.0 && l3 > la) + lb = l3; + else if (l2 >= 0.0 && l2 <= 1.0 && l2 > la) + lb = l2; + else if (l1 >= 0.0 && l1 <= 1.0 && l1 > la) + lb = l1; + else + lb = 1.0; + } if (bin == end_bin) lb = 1.0; if (la > lb) error->all(FLERR, "Error: la > lb\n");