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

modernize, avoid static buffers, use utility functions, remove debug code

Browse Source
This commit is contained in:
Axel Kohlmeyer
2021-10-22 16:00:01 -04:00
parent ede188652b
commit c08093f768
6 changed files with 35 additions and 132 deletions
Download Patch File Download Diff File Expand all files Collapse all files

14
src/MANIFOLD/fix_nve_manifold_rattle.cpp
View File

@ -140,10 +140,8 @@ FixNVEManifoldRattle::FixNVEManifoldRattle( LAMMPS *lmp, int &narg, char **arg,
if (strcmp(arg[argi], "every") == 0) {
nevery = utils::inumeric(FLERR,arg[argi+1],false,lmp);
next_output = update->ntimestep + nevery;
if (comm->me == 0) {
fprintf(screen,"Outputting every %d steps, next is %d\n",
nevery, next_output);
}
if (comm->me == 0)
utils::logmesg(lmp,"Outputting every {} steps, next is {}\n",nevery, next_output);
argi += 2;
} else if (error_on_unknown_keyword) {
error->all(FLERR,"Error parsing arg \"{}\".\n",arg[argi]);
@ -211,11 +209,9 @@ void FixNVEManifoldRattle::print_stats( const char *header )
double inv_tdiff = 1.0/( static_cast<double>(ntimestep) - stats.last_out );
stats.last_out = ntimestep;
fprintf(screen, "%s stats for time step " BIGINT_FORMAT " on %d atoms:\n",
header, ntimestep, stats.natoms);
fprintf(screen, " iters/atom: x = %f, v = %f, dofs removed %d",
x_iters * inv_tdiff, v_iters * inv_tdiff, stats.dofs_removed);
fprintf(screen,"\n");
utils::logmesg(lmp, "{} stats for time step {} on {} atoms:\n", header, ntimestep, stats.natoms);
utils::logmesg(lmp, " iters/atom: x = {}, v = {}, dofs removed = {}\n",
x_iters * inv_tdiff, v_iters * inv_tdiff, stats.dofs_removed);
}
stats.x_iters_per_atom = 0;

28
src/MANIFOLD/fix_nvt_manifold_rattle.cpp
View File

@ -32,21 +32,21 @@
------------------------------------------------------------------------- */
#include "fix_nvt_manifold_rattle.h"
#include <cstring>
#include <cmath>
#include "atom.h"
#include "force.h"
#include "update.h"
#include "error.h"
#include "group.h"
#include "citeme.h"
#include "modify.h"
#include "compute.h"
#include "error.h"
#include "force.h"
#include "group.h"
#include "modify.h"
#include "update.h"
#include "manifold.h"
#include <cmath>
#include <cstring>
using namespace LAMMPS_NS;
using namespace FixConst;
@ -115,9 +115,7 @@ FixNVTManifoldRattle::FixNVTManifoldRattle(LAMMPS *lmp, int narg, char **arg,
mtchain = utils::inumeric(FLERR, arg[argi+1],false,lmp);
argi += 2;
} else if (error_on_unknown_keyword) {
char msg[2048];
sprintf(msg,"Error parsing arg \"%s\".\n", arg[argi]);
error->all(FLERR, msg);
error->all(FLERR, "Error parsing arg \"{}\".\n", arg[argi]);
} else {
argi += 1;
}
@ -271,12 +269,8 @@ void FixNVTManifoldRattle::nhc_temp_integrate()
factor_eta = exp(-dthalf*eta_dot[0]);
if (factor_eta == 0) {
char msg[2048];
sprintf(msg, "WTF, factor_eta is 0! dthalf = %f, eta_dot[0] = %f",
dthalf, eta_dot[0]);
error->all(FLERR,msg);
}
if (factor_eta == 0)
error->all(FLERR, "factor_eta is 0! dthalf = {}, eta_dot[0] = {}", dthalf, eta_dot[0]);
nh_v_temp();

34
src/MANIFOLD/manifold_gaussian_bump.cpp
View File

@ -67,9 +67,7 @@ public:
double get_t_from_x( double xx ) const
{
if (xx < x0 || xx > x1) {
char msg[2048];
sprintf(msg, "x ( %g ) out of bounds [%g, %g]", xx, x0, x1 );
err->one(FLERR, msg);
err->one(FLERR,"x ( {} ) out of bounds [{}, {}]", xx, x0, x1 );
}
// Newton iterate to get right t.
@ -271,8 +269,6 @@ void manifold_gaussian_bump::post_param_init()
make_lut();
// test_lut();
}
@ -360,31 +356,3 @@ void manifold_gaussian_bump::lut_get_z_and_zp( double rr, double &zz,
zz = zleft * fmin + zright * frac;
zzp = zpleft * fmin + zpright * frac;
}
void manifold_gaussian_bump::test_lut()
{
double x[3], nn[3];
if (comm->me != 0) return;
FILE *fp = fopen( "test_lut_gaussian.dat", "w" );
double dx = 0.1;
for (double xx = 0; xx < 20; xx += dx) {
x[0] = xx;
x[1] = 0.0;
x[2] = 0.0;
double gg = g( x );
n( x, nn );
double taper_z;
if (xx <= rc1) {
taper_z = gaussian_bump(xx);
} else if (xx < rc2) {
taper_z = lut_get_z( xx );
} else {
taper_z = 0.0;
}
fprintf( fp, "%g %g %g %g %g %g %g\n", xx, gaussian_bump(xx), taper_z,
gg, nn[0], nn[1], nn[2] );
}
fclose(fp);
}

2
src/MANIFOLD/manifold_gaussian_bump.h
View File

@ -86,8 +86,6 @@ namespace user_manifold {
double lut_get_zp(double rr) const;
void lut_get_z_and_zp(double rr, double &zz, double &zzp) const;
void test_lut();
double taper(double);
double taper_der(double);
};

88
src/MANIFOLD/manifold_thylakoid.cpp
View File

@ -34,8 +34,6 @@ manifold_thylakoid::manifold_thylakoid( LAMMPS *lmp, int /*narg*/, char ** /*arg
// fix should call post_param_init();
}
manifold_thylakoid::~manifold_thylakoid()
{
for (std::size_t i = 0; i < parts.size(); ++i) {
@ -43,8 +41,6 @@ manifold_thylakoid::~manifold_thylakoid()
}
}
void manifold_thylakoid::post_param_init()
{
// Set coefficients:
@ -59,52 +55,21 @@ void manifold_thylakoid::post_param_init()
LB = params[1];
lB = params[2];
if (comm->me == 0) {
fprintf(screen,"My params are now: lT = %f, LT = %f, pad = %f, "
"wB = %f, LB = %f, lB = %f\n", lT, LT, pad, wB, LB, lB );
fprintf(screen,"Calling init_domains() from post_param_init().\n");
}
init_domains();
checkup();
}
void manifold_thylakoid::checkup()
{
if (comm->me == 0) {
fprintf(screen,"This is checkup of thylakoid %p\n", this);
fprintf(screen,"I have %ld parts. They are:\n", parts.size());
for (int i = 0; i < (int)parts.size(); ++i) {
fprintf(screen, "[%f, %f] x [%f, %f] x [%f, %f]\n",
parts[i]->xlo, parts[i]->xhi,
parts[i]->ylo, parts[i]->yhi,
parts[i]->zlo, parts[i]->zhi );
}
fprintf(screen,"My params are:\n");
for (int i = 0; i < NPARAMS; ++i) {
fprintf(screen,"%f\n", params[i]);
}
}
}
double manifold_thylakoid::g( const double *x )
{
int err = 0;
std::size_t idx;
thyla_part *p = get_thyla_part(x,&err,&idx);
if (err) {
char msg[2048];
sprintf(msg,"Error getting thyla_part for x = (%f, %f, %f)",x[0],x[1],x[2]);
error->one(FLERR,msg);
}
if (err) error->one(FLERR,"Error getting thyla_part for x = ({}, {}, {})",x[0],x[1],x[2]);
double con_val = p->g(x);
if (std::isfinite(con_val)) {
return con_val;
} else {
char msg[2048];
sprintf(msg,"Error, thyla_part of type %d returned %f as constraint val!",
p->type, con_val);
error->one(FLERR,msg);
error->one(FLERR,"Error, thyla_part of type {} returned {} as constraint val!", p->type, con_val);
return 0;
}
}
@ -114,20 +79,14 @@ void manifold_thylakoid::n( const double *x, double *n )
int err = 0;
std::size_t idx;
thyla_part *p = get_thyla_part(x,&err,&idx);
if (err) {
char msg[2048];
sprintf(msg,"Error getting thyla_part for x = (%f, %f, %f)",x[0],x[1],x[2]);
error->one(FLERR,msg);
}
if (err)
error->one(FLERR,"Error getting thyla_part for x = ({}, {}, {})",x[0],x[1],x[2]);
p->n(x,n);
if (std::isfinite(n[0]) && std::isfinite(n[1]) && std::isfinite(n[2])) {
return;
} else {
char msg[2048];
sprintf(msg,"Error, thyla_part of type %d returned (%f,%f,%f) as gradient!",
p->type, n[0], n[1], n[2]);
error->one(FLERR,msg);
}
} else
error->one(FLERR,"thyla_part of type {} returned ({},{},{}) as gradient!",
p->type, n[0], n[1], n[2]);
}
thyla_part *manifold_thylakoid::get_thyla_part( const double *x, int * /*err_flag*/, std::size_t *idx )
@ -140,9 +99,7 @@ thyla_part *manifold_thylakoid::get_thyla_part( const double *x, int * /*err_fla
return p;
}
}
char msg[2048];
sprintf(msg,"Could not find thyla_part for x = (%f,%f,%f)", x[0],x[1],x[2]);
error->one(FLERR,msg);
error->one(FLERR,"Could not find thyla_part for x = ({},{},{})", x[0],x[1],x[2]);
return nullptr;
}
@ -153,12 +110,9 @@ thyla_part *manifold_thylakoid::get_thyla_part( const double *x, int * /*err_fla
void manifold_thylakoid::init_domains()
{
if (wB + 2*lB > LT) {
char msg[2048];
sprintf(msg,"LT = %f not large enough to accommodate bridge with "
"wB = %f and lB = %f! %f > %f\n", LT, wB, lB, wB + 2*lB, LT);
error->one(FLERR,msg);
}
if (wB + 2*lB > LT)
error->one(FLERR,"LT = {} not large enough to accommodate bridge with "
"wB = {} and lB = {}! {} > {}\n", LT, wB, lB, wB + 2*lB, LT);
// Determine some constant coordinates:
x0 = -( 0.5*LB + lB + lT + LT + lT + pad);
@ -375,12 +329,9 @@ void manifold_thylakoid::init_domains()
parts.push_back(p);
// Check if this plane lines up with bl:
if (fabs(plr.pt[0] - bl.pt[0] + lB) > 1e-8) {
char msg[2048];
sprintf(msg,"Origins of plane left right and bridge left misaligned! %f != %f!\n",
plr.pt[0], bl.pt[0] - lB );
error->one(FLERR,msg);
}
if (fabs(plr.pt[0] - bl.pt[0] + lB) > 1e-8)
error->one(FLERR,"Origins of plane left right and bridge left misaligned! {} != {}!\n",
plr.pt[0], bl.pt[0] - lB );
// Now, for the right stack, you can mirror the other...
// To mirror them you need to invert lo[0] and hi[0] and flip their sign.
@ -428,12 +379,9 @@ void manifold_thylakoid::init_domains()
set_domain(p, prr.lo, prr.hi);
parts.push_back(p);
if (fabs(prr.pt[0] - br.pt[0] - lB) > 1e-8) {
char msg[2048];
sprintf(msg,"Origins of plane left right and bridge left misaligned! %f != %f!\n",
prr.pt[0], br.pt[0] + lB);
error->one(FLERR,msg);
}
if (fabs(prr.pt[0] - br.pt[0] - lB) > 1e-8)
error->one(FLERR,"Origins of plane left right and bridge left misaligned! {} != {}!\n",
prr.pt[0], br.pt[0] + lB);
}

1
src/MANIFOLD/manifold_thylakoid.h
View File

@ -37,7 +37,6 @@ namespace user_manifold {
virtual int nparams() { return NPARAMS; }
virtual void post_param_init();
virtual void checkup(); // Some diagnostics...
private:
void init_domains();