ICODE-ADC/lammps
4
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

stability improvements: recompute bin width for perfect match, don't test for floating point identity

Browse Source
This commit is contained in:
Axel Kohlmeyer
2022-03-05 13:46:56 -05:00
parent 15fa6049f2
commit ece25a95ad
2 changed files with 45 additions and 21 deletions
Download Patch File Download Diff File Expand all files Collapse all files

22
src/EXTRA-COMPUTE/compute_pressure_cartesian.cpp
View File

@ -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));

44
src/EXTRA-COMPUTE/compute_pressure_spherical.cpp
View File

@ -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");