forked from lijiext/lammps
Merge pull request #2333 from akohlmey/iss1109
Fix memory leaks and out-of-range memory access in USER-BOCS package
This commit is contained in:
commit
45c9478f5a
|
@ -15,23 +15,27 @@
|
|||
------------------------------------------------------------------------- */
|
||||
|
||||
#include "compute_pressure_bocs.h"
|
||||
|
||||
#include <mpi.h>
|
||||
#include <cmath>
|
||||
#include <cstring>
|
||||
#include <cstdlib>
|
||||
#include "atom.h"
|
||||
#include "update.h"
|
||||
#include "domain.h"
|
||||
#include "modify.h"
|
||||
#include "fix.h"
|
||||
#include "force.h"
|
||||
#include "pair.h"
|
||||
#include "bond.h"
|
||||
|
||||
#include "angle.h"
|
||||
#include "atom.h"
|
||||
#include "bond.h"
|
||||
#include "dihedral.h"
|
||||
#include "domain.h"
|
||||
#include "error.h"
|
||||
#include "fix.h"
|
||||
#include "fmt/format.h"
|
||||
#include "force.h"
|
||||
#include "improper.h"
|
||||
#include "kspace.h"
|
||||
#include "error.h"
|
||||
#include "memory.h"
|
||||
#include "modify.h"
|
||||
#include "pair.h"
|
||||
#include "update.h"
|
||||
|
||||
|
||||
using namespace LAMMPS_NS;
|
||||
|
@ -112,6 +116,9 @@ ComputePressureBocs::ComputePressureBocs(LAMMPS *lmp, int narg, char **arg) :
|
|||
vector = new double[size_vector];
|
||||
nvirial = 0;
|
||||
vptr = NULL;
|
||||
|
||||
splines = NULL;
|
||||
spline_length = 0;
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
@ -212,13 +219,11 @@ double ComputePressureBocs::find_index(double * grid, double value)
|
|||
|
||||
if (value >= grid[i] && value <= (grid[i] + spacing)) { return i; }
|
||||
|
||||
error->all(FLERR, fmt::format("find_index could not find value in grid for value: {}", value));
|
||||
for (int i = 0; i < gridsize; ++i)
|
||||
{
|
||||
fprintf(stderr, "grid %d: %f\n",i,grid[i]);
|
||||
}
|
||||
char * errmsg = (char *) calloc(100,sizeof(char));
|
||||
sprintf(errmsg,"Value %f does not fall within spline grid.\n",value);
|
||||
error->all(FLERR,errmsg);
|
||||
|
||||
exit(1);
|
||||
}
|
||||
|
@ -233,12 +238,12 @@ double ComputePressureBocs::get_cg_p_corr(double ** grid, int basis_type,
|
|||
int i = find_index(grid[0],vCG);
|
||||
double correction, deltax = vCG - grid[0][i];
|
||||
|
||||
if (basis_type == 1)
|
||||
if (basis_type == BASIS_LINEAR_SPLINE)
|
||||
{
|
||||
correction = grid[1][i] + (deltax) *
|
||||
( grid[1][i+1] - grid[1][i] ) / ( grid[0][i+1] - grid[0][i] );
|
||||
}
|
||||
else if (basis_type == 2)
|
||||
else if (basis_type == BASIS_CUBIC_SPLINE)
|
||||
{
|
||||
correction = grid[1][i] + (grid[2][i] * deltax) +
|
||||
(grid[3][i] * pow(deltax,2)) + (grid[4][i] * pow(deltax,3));
|
||||
|
@ -257,7 +262,7 @@ double ComputePressureBocs::get_cg_p_corr(double ** grid, int basis_type,
|
|||
void ComputePressureBocs::send_cg_info(int basis_type, int sent_N_basis,
|
||||
double *sent_phi_coeff, int sent_N_mol, double sent_vavg)
|
||||
{
|
||||
if (basis_type == 0) { p_basis_type = 0; }
|
||||
if (basis_type == BASIS_ANALYTIC) { p_basis_type = BASIS_ANALYTIC; }
|
||||
else
|
||||
{
|
||||
error->all(FLERR,"Incorrect basis type passed to ComputePressureBocs\n");
|
||||
|
@ -281,8 +286,8 @@ void ComputePressureBocs::send_cg_info(int basis_type, int sent_N_basis,
|
|||
void ComputePressureBocs::send_cg_info(int basis_type,
|
||||
double ** in_splines, int gridsize)
|
||||
{
|
||||
if (basis_type == 1) { p_basis_type = 1; }
|
||||
else if (basis_type == 2) { p_basis_type = 2; }
|
||||
if (basis_type == BASIS_LINEAR_SPLINE) { p_basis_type = BASIS_LINEAR_SPLINE; }
|
||||
else if (basis_type == BASIS_CUBIC_SPLINE) { p_basis_type = BASIS_CUBIC_SPLINE; }
|
||||
else
|
||||
{
|
||||
error->all(FLERR,"Incorrect basis type passed to ComputePressureBocs\n");
|
||||
|
@ -318,11 +323,11 @@ double ComputePressureBocs::compute_scalar()
|
|||
volume = (domain->xprd * domain->yprd * domain->zprd);
|
||||
|
||||
/* MRD NJD if block */
|
||||
if ( p_basis_type == 0 )
|
||||
if ( p_basis_type == BASIS_ANALYTIC )
|
||||
{
|
||||
correction = get_cg_p_corr(N_basis,phi_coeff,N_mol,vavg,volume);
|
||||
}
|
||||
else if ( p_basis_type == 1 || p_basis_type == 2 )
|
||||
else if ( p_basis_type == BASIS_LINEAR_SPLINE || p_basis_type == BASIS_CUBIC_SPLINE )
|
||||
{
|
||||
correction = get_cg_p_corr(splines, p_basis_type, volume);
|
||||
}
|
||||
|
|
|
@ -27,7 +27,15 @@ ComputeStyle(PRESSURE/BOCS,ComputePressureBocs)
|
|||
#include "compute.h"
|
||||
|
||||
namespace LAMMPS_NS {
|
||||
// ComputePressure -> ComputePressureBocs MRD NJD
|
||||
// Enumerate the p_basis_type magic values to improve readability:
|
||||
enum{BASIS_ANALYTIC, BASIS_LINEAR_SPLINE, BASIS_CUBIC_SPLINE};
|
||||
// Enumerate the data file column names to improve readability
|
||||
enum{VOLUME, PRESSURE_CORRECTION};
|
||||
// Declare names for the number of columns in the splines data structure to improve readability
|
||||
const int NUM_LINEAR_SPLINE_COLUMNS = 2; // linear spline columns passed to compute
|
||||
const int NUM_CUBIC_SPLINE_COLUMNS = 5; // cubic spline columns passed to compute
|
||||
|
||||
// ComputePressure -> ComputePressureBocs MRD NJD
|
||||
class ComputePressureBocs : public Compute {
|
||||
public:
|
||||
ComputePressureBocs(class LAMMPS *, int, char **);
|
||||
|
|
|
@ -15,27 +15,30 @@
|
|||
------------------------------------------------------------------------- */
|
||||
|
||||
#include "fix_bocs.h"
|
||||
|
||||
#include <cstring>
|
||||
#include <cstdlib>
|
||||
#include <cmath>
|
||||
#include <vector>
|
||||
|
||||
#include "atom.h"
|
||||
#include "citeme.h"
|
||||
#include "comm.h"
|
||||
#include "compute.h"
|
||||
#include "compute_pressure_bocs.h"
|
||||
#include "domain.h"
|
||||
#include "error.h"
|
||||
#include "fix_deform.h"
|
||||
#include "fmt/format.h"
|
||||
#include "force.h"
|
||||
#include "group.h"
|
||||
#include "comm.h"
|
||||
#include "neighbor.h"
|
||||
#include "irregular.h"
|
||||
#include "modify.h"
|
||||
#include "fix_deform.h"
|
||||
#include "compute.h"
|
||||
#include "kspace.h"
|
||||
#include "update.h"
|
||||
#include "respa.h"
|
||||
#include "domain.h"
|
||||
#include "memory.h"
|
||||
#include "error.h"
|
||||
#include "citeme.h"
|
||||
|
||||
#include "compute_pressure_bocs.h"
|
||||
#include "modify.h"
|
||||
#include "neighbor.h"
|
||||
#include "respa.h"
|
||||
#include "update.h"
|
||||
|
||||
using namespace LAMMPS_NS;
|
||||
using namespace FixConst;
|
||||
|
@ -59,8 +62,12 @@ enum{NOBIAS,BIAS};
|
|||
enum{NONE,XYZ,XY,YZ,XZ};
|
||||
enum{ISO,ANISO,TRICLINIC};
|
||||
|
||||
// NB: Keep error and warning messages less than 255 chars long.
|
||||
const int MAX_MESSAGE_LENGTH = 256;
|
||||
const int NUM_INPUT_DATA_COLUMNS = 2; // columns in the pressure correction file
|
||||
|
||||
// NB:
|
||||
// - Keep error and warning messages less than 255 chars long.
|
||||
// - Allocate your char buffer to be 1 char longer than this
|
||||
const int MAX_MESSAGE_LENGTH = 255;
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
NVT,NPH,NPT integrators for improved Nose-Hoover equations of motion
|
||||
|
@ -111,6 +118,9 @@ FixBocs::FixBocs(LAMMPS *lmp, int narg, char **arg) :
|
|||
|
||||
p_match_coeffs = NULL;
|
||||
|
||||
splines = NULL;
|
||||
spline_length = 0;
|
||||
|
||||
// turn on tilt factor scaling, whenever applicable
|
||||
|
||||
dimension = domain->dimension;
|
||||
|
@ -183,10 +193,10 @@ FixBocs::FixBocs(LAMMPS *lmp, int narg, char **arg) :
|
|||
|
||||
if ( strcmp(arg[iarg], "analytic") == 0 ) {
|
||||
if (iarg + 4 > narg) {
|
||||
error->all(FLERR,"Illegal fix bocs command. basis type analytic"
|
||||
error->all(FLERR,"Illegal fix bocs command. Basis type analytic"
|
||||
" must be followed by: avg_vol n_mol n_pmatch_coeff");
|
||||
}
|
||||
p_basis_type = 0;
|
||||
p_basis_type = BASIS_ANALYTIC;
|
||||
vavg = utils::numeric(FLERR,arg[iarg+1],false,lmp);
|
||||
N_mol = utils::inumeric(FLERR,arg[iarg+2],false,lmp);
|
||||
N_p_match = utils::inumeric(FLERR,arg[iarg+3],false,lmp);
|
||||
|
@ -200,19 +210,18 @@ FixBocs::FixBocs(LAMMPS *lmp, int narg, char **arg) :
|
|||
} else if (strcmp(arg[iarg], "linear_spline") == 0 ) {
|
||||
if (iarg+2 > narg) error->all(FLERR,"Illegal fix bocs command. "
|
||||
"Supply a file name after linear_spline.");
|
||||
p_basis_type = 1;
|
||||
p_basis_type = BASIS_LINEAR_SPLINE;
|
||||
spline_length = read_F_table( arg[iarg+1], p_basis_type );
|
||||
iarg += 2;
|
||||
} else if (strcmp(arg[iarg], "cubic_spline") == 0 ) {
|
||||
if (iarg+2 > narg) error->all(FLERR,"Illegal fix bocs command. "
|
||||
"Supply a file name after cubic_spline.");
|
||||
p_basis_type = 2;
|
||||
p_basis_type = BASIS_CUBIC_SPLINE;
|
||||
spline_length = read_F_table( arg[iarg+1], p_basis_type );
|
||||
iarg += 2;
|
||||
} else {
|
||||
char errmsg[256];
|
||||
snprintf(errmsg,256,"CG basis type %s is not recognized\nSupported "
|
||||
"basis types: analytic linear_spline cubic_spline",arg[iarg]);
|
||||
std::string errmsg = fmt::format("CG basis type {} is not recognized\nSupported "
|
||||
"basis types: analytic linear_spline cubic_spline",arg[iarg]);
|
||||
error->all(FLERR,errmsg);
|
||||
} // END NJD MRD
|
||||
} else if (strcmp(arg[iarg],"tchain") == 0) {
|
||||
|
@ -244,9 +253,8 @@ FixBocs::FixBocs(LAMMPS *lmp, int narg, char **arg) :
|
|||
if (nc_pchain < 0) error->all(FLERR,"Illegal fix bocs command");
|
||||
iarg += 2;
|
||||
} else {
|
||||
char errmsg[128];
|
||||
snprintf(errmsg,128,"Illegal fix bocs command: unrecognized keyword %s"
|
||||
,arg[iarg]);
|
||||
std::string errmsg = fmt::format("Illegal fix bocs command: unrecognized keyword {}",
|
||||
arg[iarg]);
|
||||
error->all(FLERR,errmsg);
|
||||
}
|
||||
}
|
||||
|
@ -471,6 +479,12 @@ FixBocs::~FixBocs()
|
|||
}
|
||||
}
|
||||
if (p_match_coeffs) free(p_match_coeffs);
|
||||
|
||||
// Free splines memory structure
|
||||
if (splines != NULL) {
|
||||
memory->destroy(splines);
|
||||
spline_length = 0;
|
||||
}
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
@ -537,12 +551,12 @@ void FixBocs::init()
|
|||
{
|
||||
if (pressure)
|
||||
{
|
||||
if (p_basis_type == 0)
|
||||
if (p_basis_type == BASIS_ANALYTIC)
|
||||
{
|
||||
((ComputePressureBocs *)pressure)->send_cg_info(p_basis_type,
|
||||
N_p_match, p_match_coeffs, N_mol, vavg);
|
||||
}
|
||||
else if ( p_basis_type == 1 || p_basis_type == 2 )
|
||||
else if ( p_basis_type == BASIS_LINEAR_SPLINE || p_basis_type == BASIS_CUBIC_SPLINE )
|
||||
{
|
||||
((ComputePressureBocs *)pressure)->send_cg_info(p_basis_type,
|
||||
splines, spline_length);
|
||||
|
@ -630,28 +644,42 @@ void FixBocs::init()
|
|||
// NJD MRD 2 functions
|
||||
int FixBocs::read_F_table( char *filename, int p_basis_type )
|
||||
{
|
||||
FILE *fpi;
|
||||
int N_columns = 2, n_entries = 0, i;
|
||||
float f1, f2;
|
||||
int test_sscanf;
|
||||
double **data = (double **) calloc(N_columns,sizeof(double *));
|
||||
char * line = (char *) calloc(200,sizeof(char));
|
||||
|
||||
std::string message;
|
||||
double **data;
|
||||
bool badInput = false;
|
||||
char badDataMsg[MAX_MESSAGE_LENGTH];
|
||||
fpi = fopen(filename,"r");
|
||||
if (fpi)
|
||||
{
|
||||
while (fgets(line,199,fpi)) { ++n_entries; }
|
||||
|
||||
for (i = 0; i < N_columns; ++i)
|
||||
{
|
||||
data[i] = (double *) calloc(n_entries,sizeof(double));
|
||||
int numEntries = 0;
|
||||
FILE *fpi = fopen(filename,"r");
|
||||
if (fpi) {
|
||||
// Old code read the input file twice. Now we simply
|
||||
// read all the lines from the input file into a string vector,
|
||||
// then work with the data in-memory rather than do a second pass
|
||||
// through the file.
|
||||
// NB: LAMMPS coding guidelines prefer cstdio so we are intentionally
|
||||
// foregoing reading with getline
|
||||
if (comm->me == 0) {
|
||||
error->message(FLERR, fmt::format("INFO: About to read data file: {}", filename));
|
||||
}
|
||||
|
||||
// Don't need to re-open the file to make a second pass through it
|
||||
// simply rewind to beginning
|
||||
rewind(fpi);
|
||||
// Data file lines hold two floating point numbers.
|
||||
// Line length we allocate should be long enough without being too long.
|
||||
// 128 seems safe for a line we expect to be < 30 chars.
|
||||
const int MAX_F_TABLE_LINE_LENGTH = 128;
|
||||
char line[MAX_F_TABLE_LINE_LENGTH];
|
||||
std::vector<std::string> inputLines;
|
||||
while (fgets(line, MAX_F_TABLE_LINE_LENGTH, fpi)) {
|
||||
inputLines.push_back(std::string(line));
|
||||
}
|
||||
fclose(fpi);
|
||||
|
||||
numEntries = inputLines.size();
|
||||
if (comm->me == 0) {
|
||||
error->message(FLERR, fmt::format("INFO: Read {} lines from file", numEntries));
|
||||
}
|
||||
|
||||
|
||||
// Allocate memory for the two dimensional matrix
|
||||
// that holds data from the input file.
|
||||
memory->create(data, NUM_INPUT_DATA_COLUMNS, numEntries, "data");
|
||||
|
||||
double stdVolumeInterval = 0.0;
|
||||
double currVolumeInterval = 0.0;
|
||||
|
@ -659,125 +687,157 @@ int FixBocs::read_F_table( char *filename, int p_basis_type )
|
|||
// The literature indicates getting this value right in the
|
||||
// general case can be pretty complicated. I don't think it
|
||||
// needs to be complicated here, though. At least based on the
|
||||
// sample data I've seen where the volume values are fairly
|
||||
// large.
|
||||
// sample data I've seen where the volume values are fairly large.
|
||||
const double volumeIntervalTolerance = 0.001;
|
||||
n_entries = 0;
|
||||
while( fgets(line,199,fpi)) {
|
||||
++n_entries;
|
||||
test_sscanf = sscanf(line," %f , %f ",&f1, &f2);
|
||||
int lineNum = 0; // this value is only for message
|
||||
int numBadVolumeIntervals = 0; // count these for message
|
||||
float f1, f2;
|
||||
int test_sscanf;
|
||||
for (int i = 0; i < inputLines.size(); ++i) {
|
||||
lineNum++; // count each line processed now so lineNum messages can be 1-based
|
||||
test_sscanf = sscanf(inputLines.at(i).c_str()," %f , %f ",&f1, &f2);
|
||||
if (test_sscanf == 2)
|
||||
{
|
||||
data[0][n_entries-1] = (double) f1;
|
||||
data[1][n_entries-1] = (double) f2;
|
||||
if (n_entries == 2) {
|
||||
stdVolumeInterval = data[0][n_entries-1] - data[0][n_entries-2];
|
||||
//if (comm->me == 0) {
|
||||
// error->message(FLERR, fmt::format("INFO: f1 = {}, f2 = {}", f1, f2));
|
||||
//}
|
||||
data[VOLUME][i] = (double)f1;
|
||||
data[PRESSURE_CORRECTION][i] = (double)f2;
|
||||
if (i == 1)
|
||||
{
|
||||
// second entry is used to compute the validation interval used below
|
||||
stdVolumeInterval = data[VOLUME][i] - data[VOLUME][i-1];
|
||||
//if (comm->me == 0) {
|
||||
// error->message(FLERR, fmt::format("INFO: standard volume interval computed: {}", stdVolumeInterval));
|
||||
//}
|
||||
}
|
||||
else if (n_entries > 2) {
|
||||
currVolumeInterval = data[0][n_entries-1] - data[0][n_entries-2];
|
||||
else if (i > 1)
|
||||
{
|
||||
// after second entry, all intervals are validated
|
||||
currVolumeInterval = data[VOLUME][i] - data[VOLUME][i-1];
|
||||
//if (comm->me == 0) {
|
||||
// error->message(FLERR, fmt::format("INFO: current volume interval: {}", currVolumeInterval));
|
||||
//}
|
||||
if (fabs(currVolumeInterval - stdVolumeInterval) > volumeIntervalTolerance) {
|
||||
snprintf(badDataMsg,MAX_MESSAGE_LENGTH,
|
||||
"BAD VOLUME INTERVAL: spline analysis requires uniform"
|
||||
" volume distribution, found inconsistent volume"
|
||||
" differential, line %d of file %s\n\tline: %s",
|
||||
n_entries,filename,line);
|
||||
error->message(FLERR,badDataMsg);
|
||||
if (comm->me == 0) {
|
||||
message = fmt::format("Bad volume interval. Spline analysis requires uniform"
|
||||
" volume distribution, found inconsistent volume"
|
||||
" differential, line {} of file {}\n\tline: {}",
|
||||
lineNum, filename, inputLines.at(i));
|
||||
error->warning(FLERR, message);
|
||||
}
|
||||
badInput = true;
|
||||
numBadVolumeIntervals++;
|
||||
}
|
||||
// no concluding else is intentional: i = 0, first line, no interval to validate
|
||||
}
|
||||
// no else -- first entry is simply ignored
|
||||
}
|
||||
else
|
||||
{
|
||||
snprintf(badDataMsg,MAX_MESSAGE_LENGTH,
|
||||
"BAD INPUT FORMAT: did not find 2 comma separated numeric"
|
||||
" values in line %d of file %s\n\tline: %s",
|
||||
n_entries,filename,line);
|
||||
error->message(FLERR,badDataMsg);
|
||||
if (comm->me == 0) {
|
||||
message = fmt::format("Bad input format: did not find 2 comma separated numeric"
|
||||
" values in line {} of file {}\n\tline: {}",
|
||||
lineNum, filename, inputLines.at(i));
|
||||
error->warning(FLERR, message);
|
||||
}
|
||||
badInput = true;
|
||||
}
|
||||
}
|
||||
fclose(fpi);
|
||||
}
|
||||
else {
|
||||
char errmsg[MAX_MESSAGE_LENGTH];
|
||||
snprintf(errmsg,MAX_MESSAGE_LENGTH,"Unable to open file: %s\n",filename);
|
||||
error->all(FLERR,errmsg);
|
||||
}
|
||||
|
||||
if (badInput) {
|
||||
char errmsg[MAX_MESSAGE_LENGTH];
|
||||
snprintf(errmsg,MAX_MESSAGE_LENGTH,
|
||||
"Bad volume / pressure-correction data: %s\nSee details above",filename);
|
||||
error->all(FLERR,errmsg);
|
||||
}
|
||||
|
||||
if (p_basis_type == 1)
|
||||
{
|
||||
splines = (double **) calloc(2,sizeof(double *));
|
||||
splines[0] = (double *) calloc(n_entries,sizeof(double));
|
||||
splines[1] = (double *) calloc(n_entries,sizeof(double));
|
||||
int idxa, idxb;
|
||||
for (idxa = 0; idxa < 2; ++idxa)
|
||||
{
|
||||
for (idxb = 0; idxb < n_entries; ++idxb)
|
||||
if (badInput)
|
||||
{
|
||||
splines[idxa][idxb] = data[idxa][idxb];
|
||||
numBadVolumeIntervals++;
|
||||
}
|
||||
}
|
||||
|
||||
if (numBadVolumeIntervals > 0 && comm->me == 0) {
|
||||
error->message(FLERR, fmt::format("INFO: total number bad volume intervals = {}", numBadVolumeIntervals));
|
||||
}
|
||||
}
|
||||
else if (p_basis_type == 2)
|
||||
else {
|
||||
error->all(FLERR,fmt::format("ERROR: Unable to open file: {}", filename));
|
||||
}
|
||||
|
||||
if (badInput && comm->me == 0) {
|
||||
error->warning(FLERR,fmt::format("Bad volume / pressure-correction data: {}\nSee details above", filename));
|
||||
}
|
||||
|
||||
if (p_basis_type == BASIS_LINEAR_SPLINE)
|
||||
{
|
||||
spline_length = n_entries;
|
||||
build_cubic_splines(data);
|
||||
n_entries -= 1;
|
||||
spline_length = numEntries;
|
||||
numEntries = build_linear_splines(data);
|
||||
}
|
||||
else if (p_basis_type == BASIS_CUBIC_SPLINE)
|
||||
{
|
||||
spline_length = numEntries;
|
||||
numEntries = build_cubic_splines(data);
|
||||
}
|
||||
else
|
||||
{
|
||||
char errmsg[MAX_MESSAGE_LENGTH];
|
||||
snprintf(errmsg, MAX_MESSAGE_LENGTH,
|
||||
"ERROR: invalid p_basis_type value of %d in read_F_table",
|
||||
p_basis_type);
|
||||
error->all(FLERR,errmsg);
|
||||
error->all(FLERR,fmt::format("ERROR: invalid p_basis_type value of {} in read_F_table", p_basis_type));
|
||||
}
|
||||
// cleanup
|
||||
for (i = 0; i < N_columns; ++i) {
|
||||
free(data[i]);
|
||||
}
|
||||
free(data);
|
||||
return n_entries;
|
||||
|
||||
memory->destroy(data);
|
||||
return numEntries;
|
||||
}
|
||||
|
||||
void FixBocs::build_cubic_splines( double **data )
|
||||
int FixBocs::build_linear_splines(double **data) {
|
||||
//if (comm->me == 0) {
|
||||
//error->message(FLERR, fmt::format("INFO: entering build_linear_splines, spline_length = {}", spline_length));
|
||||
//}
|
||||
splines = (double **) calloc(NUM_LINEAR_SPLINE_COLUMNS,sizeof(double *));
|
||||
splines[VOLUME] = (double *) calloc(spline_length,sizeof(double));
|
||||
splines[PRESSURE_CORRECTION] = (double *) calloc(spline_length,sizeof(double));
|
||||
|
||||
for (int i = 0; i < spline_length; ++i)
|
||||
{
|
||||
splines[VOLUME][i] = data[VOLUME][i];
|
||||
splines[PRESSURE_CORRECTION][i] = data[PRESSURE_CORRECTION][i];
|
||||
}
|
||||
|
||||
if (comm->me == 0) {
|
||||
error->message(FLERR, fmt::format("INFO: leaving build_linear_splines, spline_length = {}", spline_length));
|
||||
}
|
||||
|
||||
return spline_length;
|
||||
}
|
||||
|
||||
int FixBocs::build_cubic_splines(double **data)
|
||||
{
|
||||
double *a, *b, *d, *h, *alpha, *c, *l, *mu, *z;
|
||||
//if (comm->me == 0) {
|
||||
//error->message(FLERR, fmt::format("INFO: entering build_cubic_splines, spline_length = {}", spline_length));
|
||||
//}
|
||||
int n = spline_length;
|
||||
double alpha_i;
|
||||
a = (double *) calloc(n,sizeof(double));
|
||||
b = (double *) calloc(n+1,sizeof(double));
|
||||
d = (double *) calloc(n+1,sizeof(double));
|
||||
h = (double *) calloc(n,sizeof(double));
|
||||
alpha = (double *) calloc(n,sizeof(double));
|
||||
c = (double *) calloc(n+1,sizeof(double));
|
||||
l = (double *) calloc(n,sizeof(double));
|
||||
mu = (double *) calloc(n,sizeof(double));
|
||||
z = (double *) calloc(n,sizeof(double));
|
||||
int idx;
|
||||
double *a, *b, *d, *h, *alpha, *c, *l, *mu, *z;
|
||||
// 2020-07-17 ag:
|
||||
// valgrind says that we read/write a[n] down in the
|
||||
// for(int j=n-1; j>=0; j--) loop below
|
||||
// and I agree.
|
||||
// So the size of a must be n+1, not n as was found
|
||||
// in the original code.
|
||||
memory->create(a, n+1, "a");
|
||||
memory->create(b, n+1, "b");
|
||||
memory->create(c, n+1, "c");
|
||||
memory->create(d, n+1, "d");
|
||||
|
||||
memory->create(h, n, "h");
|
||||
memory->create(alpha, n, "alpha");
|
||||
memory->create(l, n, "l");
|
||||
memory->create(mu, n, "mu");
|
||||
memory->create(z, n, "z");
|
||||
|
||||
for (int i=0; i<n; i++)
|
||||
{
|
||||
a[i] = data[1][i];
|
||||
b[i] = 0.0;
|
||||
d[i] = 0.0;
|
||||
|
||||
if (i<(n-1))
|
||||
{
|
||||
h[i] = (data[0][i+1] - data[0][i]);
|
||||
}
|
||||
|
||||
double alpha_i;
|
||||
if (i>1 && i<(n-1))
|
||||
{
|
||||
alpha_i = (3.0 / h[i]) * ( data[1][i+1] - data[1][i]) - (3.0 / h[i-1] )
|
||||
* ( data[1][i] - data[1][i-1] );
|
||||
* ( data[1][i] - data[1][i-1] );
|
||||
alpha[i-1] = alpha_i;
|
||||
}
|
||||
}
|
||||
|
@ -795,8 +855,12 @@ void FixBocs::build_cubic_splines( double **data )
|
|||
mu[n-1] = 0.0;
|
||||
z[n-1] = 0.0;
|
||||
|
||||
c[n] = 0.0;
|
||||
// 2020-07-17 ag: We've been using an uninitialized value for a[n]
|
||||
// That seems like a bad idea. This may not be the right value
|
||||
// but its a value.
|
||||
a[n] = 0.0;
|
||||
b[n] = 0.0;
|
||||
c[n] = 0.0;
|
||||
d[n] = 0.0;
|
||||
|
||||
for(int j=n-1; j>=0; j--)
|
||||
|
@ -807,21 +871,35 @@ void FixBocs::build_cubic_splines( double **data )
|
|||
|
||||
d[j] = (c[j+1]-c[j])/(3.0 * h[j]);
|
||||
}
|
||||
splines = (double **) calloc(5,sizeof(double *));
|
||||
|
||||
for ( idx = 0; idx < 5; ++idx)
|
||||
{
|
||||
splines[idx] = (double *) calloc(n-1,sizeof(double));
|
||||
}
|
||||
idx = 0;
|
||||
for ( idx = 0; idx < n - 1; ++idx)
|
||||
int numSplines = n - 1;
|
||||
memory->create(splines, NUM_CUBIC_SPLINE_COLUMNS, numSplines, "splines");
|
||||
for (int idx = 0; idx < numSplines; ++idx)
|
||||
{
|
||||
splines[0][idx] = data[0][idx];
|
||||
splines[1][idx] = a[idx];
|
||||
splines[2][idx] = b[idx];
|
||||
splines[3][idx] = c[idx];
|
||||
splines[4][idx] = d[idx];
|
||||
splines[0][idx] = data[0][idx];
|
||||
}
|
||||
|
||||
memory->destroy(a);
|
||||
memory->destroy(b);
|
||||
memory->destroy(c);
|
||||
memory->destroy(d);
|
||||
|
||||
memory->destroy(h);
|
||||
memory->destroy(alpha);
|
||||
memory->destroy(l);
|
||||
memory->destroy(mu);
|
||||
memory->destroy(z);
|
||||
|
||||
if (comm->me == 0) {
|
||||
error->message(FLERR, fmt::format("INFO: leaving build_cubic_splines, numSplines = {}", numSplines));
|
||||
}
|
||||
|
||||
// Tell the caller how many splines we created
|
||||
return numSplines;
|
||||
}
|
||||
// END NJD MRD 2 functions
|
||||
|
||||
|
@ -1519,20 +1597,21 @@ int FixBocs::modify_param(int narg, char **arg)
|
|||
|
||||
if (p_match_flag) // NJD MRD
|
||||
{
|
||||
if ( p_basis_type == 0 )
|
||||
if ( p_basis_type == BASIS_ANALYTIC )
|
||||
{
|
||||
((ComputePressureBocs *)pressure)->send_cg_info(p_basis_type, N_p_match,
|
||||
p_match_coeffs, N_mol, vavg);
|
||||
}
|
||||
else if ( p_basis_type == 1 || p_basis_type == 2 )
|
||||
else if ( p_basis_type == BASIS_LINEAR_SPLINE || p_basis_type == BASIS_CUBIC_SPLINE )
|
||||
{
|
||||
((ComputePressureBocs *)pressure)->send_cg_info(p_basis_type, splines,
|
||||
spline_length );
|
||||
((ComputePressureBocs *)pressure)->send_cg_info(p_basis_type, splines, spline_length );
|
||||
}
|
||||
}
|
||||
|
||||
if (pressure->pressflag == 0)
|
||||
error->all(FLERR,"Fix_modify pressure ID does not compute pressure");
|
||||
{
|
||||
error->all(FLERR, "Fix_modify pressure ID does not compute pressure");
|
||||
}
|
||||
return 2;
|
||||
}
|
||||
|
||||
|
|
|
@ -26,7 +26,6 @@ FixStyle(bocs,FixBocs)
|
|||
|
||||
#include "fix.h"
|
||||
|
||||
|
||||
namespace LAMMPS_NS {
|
||||
|
||||
class FixBocs : public Fix {
|
||||
|
@ -151,7 +150,8 @@ class FixBocs : public Fix {
|
|||
void nhc_press_integrate();
|
||||
|
||||
int read_F_table(char *, int);
|
||||
void build_cubic_splines(double **);
|
||||
int build_linear_splines(double **);
|
||||
int build_cubic_splines(double **);
|
||||
|
||||
virtual void nve_x(); // may be overwritten by child classes
|
||||
virtual void nve_v();
|
||||
|
@ -180,6 +180,10 @@ Self-explanatory. Check the input script syntax and compare to the
|
|||
documentation for the command. You can use -echo screen as a
|
||||
command-line option when running LAMMPS to see the offending line.
|
||||
|
||||
E: CG basis type XXX is not recognized
|
||||
|
||||
See second line of message for supported basis types.
|
||||
|
||||
E: Target temperature for fix bocs cannot be 0.0
|
||||
|
||||
Self-explanatory.
|
||||
|
|
Loading…
Reference in New Issue