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

MAINT: Turned 'global' options into a 'thermo yes'/'thermo no' option that enables dumping of thermo data to the netcdf file.

Browse Source
This commit is contained in:
Lars Pastewka
2017-06-08 22:43:10 +02:00
parent 36c8b26fef
commit be476c9e1d
2 changed files with 97 additions and 251 deletions
Download Patch File Download Diff File Expand all files Collapse all files

324
src/USER-NETCDF/dump_netcdf.cpp
View File

@ -55,11 +55,13 @@
#include "universe.h"
#include "variable.h"
#include "force.h"
#include "output.h"
#include "thermo.h"
using namespace LAMMPS_NS;
using namespace MathConst;
enum{INT,DOUBLE}; // same as in dump_custom.cpp
enum{INT,FLOAT,BIGINT}; // same as in thermo.cpp
const char NC_FRAME_STR[] = "frame";
const char NC_SPATIAL_STR[] = "spatial";
@ -208,15 +210,15 @@ DumpNetCDF::DumpNetCDF(LAMMPS *lmp, int narg, char **arg) :
perat[inc].field[idim] = i;
}
n_perframe = 0;
perframe = NULL;
n_buffer = 0;
int_buffer = NULL;
double_buffer = NULL;
double_precision = false;
thermo = false;
thermovar = NULL;
framei = 0;
}
@ -227,8 +229,7 @@ DumpNetCDF::~DumpNetCDF()
closefile();
delete [] perat;
if (n_perframe > 0)
delete [] perframe;
if (thermovar) delete [] thermovar;
if (int_buffer) memory->sfree(int_buffer);
if (double_buffer) memory->sfree(double_buffer);
@ -238,6 +239,11 @@ DumpNetCDF::~DumpNetCDF()
void DumpNetCDF::openfile()
{
if (thermo && !singlefile_opened) {
if (thermovar) delete [] thermovar;
thermovar = new int[output->thermo->nfield];
}
// now the computes and fixes have been initialized, so we can query
// for the size of vector quantities
for (int i = 0; i < n_perat; i++) {
@ -289,30 +295,30 @@ void DumpNetCDF::openfile()
// dimensions
NCERRX( nc_inq_dimid(ncid, NC_FRAME_STR, &frame_dim), NC_FRAME_STR );
NCERRX( nc_inq_dimid(ncid, NC_SPATIAL_STR, &spatial_dim),
NC_SPATIAL_STR );
NC_SPATIAL_STR );
NCERRX( nc_inq_dimid(ncid, NC_VOIGT_STR, &Voigt_dim), NC_VOIGT_STR );
NCERRX( nc_inq_dimid(ncid, NC_ATOM_STR, &atom_dim), NC_ATOM_STR );
NCERRX( nc_inq_dimid(ncid, NC_CELL_SPATIAL_STR, &cell_spatial_dim),
NC_CELL_SPATIAL_STR );
NC_CELL_SPATIAL_STR );
NCERRX( nc_inq_dimid(ncid, NC_CELL_ANGULAR_STR, &cell_angular_dim),
NC_CELL_ANGULAR_STR );
NC_CELL_ANGULAR_STR );
NCERRX( nc_inq_dimid(ncid, NC_LABEL_STR, &label_dim), NC_LABEL_STR );
// default variables
NCERRX( nc_inq_varid(ncid, NC_SPATIAL_STR, &spatial_var),
NC_SPATIAL_STR );
NC_SPATIAL_STR );
NCERRX( nc_inq_varid(ncid, NC_CELL_SPATIAL_STR, &cell_spatial_var),
NC_CELL_SPATIAL_STR);
NC_CELL_SPATIAL_STR);
NCERRX( nc_inq_varid(ncid, NC_CELL_ANGULAR_STR, &cell_angular_var),
NC_CELL_ANGULAR_STR);
NC_CELL_ANGULAR_STR);
NCERRX( nc_inq_varid(ncid, NC_TIME_STR, &time_var), NC_TIME_STR );
NCERRX( nc_inq_varid(ncid, NC_CELL_ORIGIN_STR, &cell_origin_var),
NC_CELL_ORIGIN_STR );
NC_CELL_ORIGIN_STR );
NCERRX( nc_inq_varid(ncid, NC_CELL_LENGTHS_STR, &cell_lengths_var),
NC_CELL_LENGTHS_STR);
NC_CELL_LENGTHS_STR);
NCERRX( nc_inq_varid(ncid, NC_CELL_ANGLES_STR, &cell_angles_var),
NC_CELL_ANGLES_STR);
NC_CELL_ANGLES_STR);
// variables specified in the input file
for (int i = 0; i < n_perat; i++) {
@ -334,9 +340,12 @@ void DumpNetCDF::openfile()
}
// perframe variables
for (int i = 0; i < n_perframe; i++) {
NCERRX( nc_inq_varid(ncid, perframe[i].name, &perframe[i].var),
perframe[i].name );
if (thermo) {
Thermo *th = output->thermo;
for (int i = 0; i < th->nfield; i++) {
NCERRX( nc_inq_varid(ncid, th->keyword[i], &thermovar[i]),
th->keyword[i] );
}
}
size_t nframes;
@ -355,48 +364,48 @@ void DumpNetCDF::openfile()
singlefile_opened = 1;
NCERRX( nc_create(filename, NC_64BIT_DATA, &ncid),
filename );
filename );
// dimensions
NCERRX( nc_def_dim(ncid, NC_FRAME_STR, NC_UNLIMITED, &frame_dim),
NC_FRAME_STR );
NC_FRAME_STR );
NCERRX( nc_def_dim(ncid, NC_SPATIAL_STR, 3, &spatial_dim),
NC_SPATIAL_STR );
NC_SPATIAL_STR );
NCERRX( nc_def_dim(ncid, NC_VOIGT_STR, 6, &Voigt_dim),
NC_VOIGT_STR );
NC_VOIGT_STR );
NCERRX( nc_def_dim(ncid, NC_ATOM_STR, ntotalgr, &atom_dim),
NC_ATOM_STR );
NC_ATOM_STR );
NCERRX( nc_def_dim(ncid, NC_CELL_SPATIAL_STR, 3, &cell_spatial_dim),
NC_CELL_SPATIAL_STR );
NC_CELL_SPATIAL_STR );
NCERRX( nc_def_dim(ncid, NC_CELL_ANGULAR_STR, 3, &cell_angular_dim),
NC_CELL_ANGULAR_STR );
NC_CELL_ANGULAR_STR );
NCERRX( nc_def_dim(ncid, NC_LABEL_STR, 10, &label_dim),
NC_LABEL_STR );
NC_LABEL_STR );
// default variables
dims[0] = spatial_dim;
NCERRX( nc_def_var(ncid, NC_SPATIAL_STR, NC_CHAR, 1, dims, &spatial_var),
NC_SPATIAL_STR );
NC_SPATIAL_STR );
NCERRX( nc_def_var(ncid, NC_CELL_SPATIAL_STR, NC_CHAR, 1, dims,
&cell_spatial_var), NC_CELL_SPATIAL_STR );
&cell_spatial_var), NC_CELL_SPATIAL_STR );
dims[0] = spatial_dim;
dims[1] = label_dim;
NCERRX( nc_def_var(ncid, NC_CELL_ANGULAR_STR, NC_CHAR, 2, dims,
&cell_angular_var), NC_CELL_ANGULAR_STR );
&cell_angular_var), NC_CELL_ANGULAR_STR );
dims[0] = frame_dim;
NCERRX( nc_def_var(ncid, NC_TIME_STR, NC_DOUBLE, 1, dims, &time_var),
NC_TIME_STR);
NC_TIME_STR);
dims[0] = frame_dim;
dims[1] = cell_spatial_dim;
NCERRX( nc_def_var(ncid, NC_CELL_ORIGIN_STR, NC_DOUBLE, 2, dims,
&cell_origin_var), NC_CELL_ORIGIN_STR );
&cell_origin_var), NC_CELL_ORIGIN_STR );
NCERRX( nc_def_var(ncid, NC_CELL_LENGTHS_STR, NC_DOUBLE, 2, dims,
&cell_lengths_var), NC_CELL_LENGTHS_STR );
&cell_lengths_var), NC_CELL_LENGTHS_STR );
dims[0] = frame_dim;
dims[1] = cell_angular_dim;
NCERRX( nc_def_var(ncid, NC_CELL_ANGLES_STR, NC_DOUBLE, 2, dims,
&cell_angles_var), NC_CELL_ANGLES_STR );
&cell_angles_var), NC_CELL_ANGLES_STR );
// variables specified in the input file
dims[0] = frame_dim;
@ -423,17 +432,17 @@ void DumpNetCDF::openfile()
// this is a tensor in Voigt notation
dims[2] = Voigt_dim;
NCERRX( nc_def_var(ncid, perat[i].name, xtype, 2, dims+1,
&perat[i].var), perat[i].name );
&perat[i].var), perat[i].name );
}
else if (perat[i].dims == 3) {
// this is a vector, we need to store x-, y- and z-coordinates
dims[2] = spatial_dim;
NCERRX( nc_def_var(ncid, perat[i].name, xtype, 2, dims+1,
&perat[i].var), perat[i].name );
&perat[i].var), perat[i].name );
}
else if (perat[i].dims == 1) {
NCERRX( nc_def_var(ncid, perat[i].name, xtype, 1, dims+1,
&perat[i].var), perat[i].name );
&perat[i].var), perat[i].name );
}
else {
char errstr[1024];
@ -448,17 +457,17 @@ void DumpNetCDF::openfile()
// this is a tensor in Voigt notation
dims[2] = Voigt_dim;
NCERRX( nc_def_var(ncid, perat[i].name, xtype, 3, dims,
&perat[i].var), perat[i].name );
&perat[i].var), perat[i].name );
}
else if (perat[i].dims == 3) {
// this is a vector, we need to store x-, y- and z-coordinates
dims[2] = spatial_dim;
NCERRX( nc_def_var(ncid, perat[i].name, xtype, 3, dims,
&perat[i].var), perat[i].name );
&perat[i].var), perat[i].name );
}
else if (perat[i].dims == 1) {
NCERRX( nc_def_var(ncid, perat[i].name, xtype, 2, dims,
&perat[i].var), perat[i].name );
&perat[i].var), perat[i].name );
}
else {
char errstr[1024];
@ -471,14 +480,21 @@ void DumpNetCDF::openfile()
}
// perframe variables
for (int i = 0; i < n_perframe; i++) {
if (perframe[i].type == THIS_IS_A_BIGINT) {
NCERRX( nc_def_var(ncid, perframe[i].name, NC_LONG, 1, dims,
&perframe[i].var), perframe[i].name );
}
else {
NCERRX( nc_def_var(ncid, perframe[i].name, NC_DOUBLE, 1, dims,
&perframe[i].var), perframe[i].name );
if (thermo) {
Thermo *th = output->thermo;
for (int i = 0; i < th->nfield; i++) {
if (th->vtype[i] == FLOAT) {
NCERRX( nc_def_var(ncid, th->keyword[i], NC_DOUBLE, 1, dims,
&thermovar[i]), th->keyword[i] );
}
else if (th->vtype[i] == INT) {
NCERRX( nc_def_var(ncid, th->keyword[i], NC_INT, 1, dims,
&thermovar[i]), th->keyword[i] );
}
else if (th->vtype[i] == BIGINT) {
NCERRX( nc_def_var(ncid, th->keyword[i], NC_LONG, 1, dims,
&thermovar[i]), th->keyword[i] );
}
}
}
@ -622,46 +638,30 @@ void DumpNetCDF::write()
start[0] = framei-1;
start[1] = 0;
for (int i = 0; i < n_perframe; i++) {
if (perframe[i].type == THIS_IS_A_BIGINT) {
bigint data;
(this->*perframe[i].compute)((void*) &data);
if (filewriter)
if (thermo) {
Thermo *th = output->thermo;
for (int i = 0; i < th->nfield; i++) {
th->call_vfunc(i);
if (filewriter) {
if (th->vtype[i] == FLOAT) {
NCERRX( nc_put_var1_double(ncid, thermovar[i], start,
&th->dvalue),
th->keyword[i] );
}
else if (th->vtype[i] == INT) {
NCERRX( nc_put_var1_int(ncid, thermovar[i], start, &th->ivalue),
th->keyword[i] );
}
else if (th->vtype[i] == BIGINT) {
#if defined(LAMMPS_SMALLBIG) || defined(LAMMPS_BIGBIG)
NCERR( nc_put_var1_long(ncid, perframe[i].var, start, &data) );
NCERRX( nc_put_var1_long(ncid, thermovar[i], start, &th->bivalue),
th->keyword[i] );
#else
NCERR( nc_put_var1_int(ncid, perframe[i].var, start, &data) );
NCERRX( nc_put_var1_int(ncid, thermovar[i], start, &th->bivalue),
th->keyword[i] );
#endif
}
else {
double data;
int j = perframe[i].index;
int idim = perframe[i].dim;
if (perframe[i].type == THIS_IS_A_COMPUTE) {
if (idim >= 0) {
modify->compute[j]->compute_vector();
data = modify->compute[j]->vector[idim];
}
else
data = modify->compute[j]->compute_scalar();
}
else if (perframe[i].type == THIS_IS_A_FIX) {
if (idim >= 0) {
data = modify->fix[j]->compute_vector(idim);
}
else
data = modify->fix[j]->compute_scalar();
}
else if (perframe[i].type == THIS_IS_A_VARIABLE) {
j = input->variable->find(perframe[i].id);
data = input->variable->compute_equal(j);
}
if (filewriter)
NCERR( nc_put_var1_double(ncid, perframe[i].var, start, &data) );
}
}
@ -908,126 +908,19 @@ int DumpNetCDF::modify_param(int narg, char **arg)
iarg++;
return 2;
}
else if (strcmp(arg[iarg],"global") == 0) {
// "perframe" quantities, i.e. not per-atom stuff
else if (strcmp(arg[iarg],"thermo") == 0) {
iarg++;
n_perframe = narg-iarg;
perframe = new nc_perframe_t[n_perframe];
for (int i = 0; iarg < narg; iarg++, i++) {
int n;
char *suffix=NULL;
if (!strcmp(arg[iarg],"step")) {
perframe[i].type = THIS_IS_A_BIGINT;
perframe[i].compute = &DumpNetCDF::compute_step;
strcpy(perframe[i].name, arg[iarg]);
}
else if (!strcmp(arg[iarg],"elapsed")) {
perframe[i].type = THIS_IS_A_BIGINT;
perframe[i].compute = &DumpNetCDF::compute_elapsed;
strcpy(perframe[i].name, arg[iarg]);
}
else if (!strcmp(arg[iarg],"elaplong")) {
perframe[i].type = THIS_IS_A_BIGINT;
perframe[i].compute = &DumpNetCDF::compute_elapsed_long;
strcpy(perframe[i].name, arg[iarg]);
}
else {
n = strlen(arg[iarg]);
if (n > 2) {
suffix = new char[n-1];
strcpy(suffix, arg[iarg]+2);
}
else {
char errstr[1024];
sprintf(errstr, "perframe quantity '%s' must thermo quantity or "
"compute, fix or variable", arg[iarg]);
error->all(FLERR,errstr);
}
if (!strncmp(arg[iarg], "c_", 2)) {
int idim = -1;
char *ptr = strchr(suffix, '[');
if (ptr) {
if (suffix[strlen(suffix)-1] != ']')
error->all(FLERR,"Missing ']' in dump modify command");
*ptr = '\0';
idim = ptr[1] - '1';
}
n = modify->find_compute(suffix);
if (n < 0)
error->all(FLERR,"Could not find dump modify compute ID");
if (modify->compute[n]->peratom_flag != 0)
error->all(FLERR,"Dump modify compute ID computes per-atom info");
if (idim >= 0 && modify->compute[n]->vector_flag == 0)
error->all(FLERR,"Dump modify compute ID does not compute vector");
if (idim < 0 && modify->compute[n]->scalar_flag == 0)
error->all(FLERR,"Dump modify compute ID does not compute scalar");
perframe[i].type = THIS_IS_A_COMPUTE;
perframe[i].dim = idim;
perframe[i].index = n;
strcpy(perframe[i].name, arg[iarg]);
}
else if (!strncmp(arg[iarg], "f_", 2)) {
int idim = -1;
char *ptr = strchr(suffix, '[');
if (ptr) {
if (suffix[strlen(suffix)-1] != ']')
error->all(FLERR,"Missing ']' in dump modify command");
*ptr = '\0';
idim = ptr[1] - '1';
}
n = modify->find_fix(suffix);
if (n < 0)
error->all(FLERR,"Could not find dump modify fix ID");
if (modify->fix[n]->peratom_flag != 0)
error->all(FLERR,"Dump modify fix ID computes per-atom info");
if (idim >= 0 && modify->fix[n]->vector_flag == 0)
error->all(FLERR,"Dump modify fix ID does not compute vector");
if (idim < 0 && modify->fix[n]->scalar_flag == 0)
error->all(FLERR,"Dump modify fix ID does not compute vector");
perframe[i].type = THIS_IS_A_FIX;
perframe[i].dim = idim;
perframe[i].index = n;
strcpy(perframe[i].name, arg[iarg]);
}
else if (!strncmp(arg[iarg], "v_", 2)) {
n = input->variable->find(suffix);
if (n < 0)
error->all(FLERR,"Could not find dump modify variable ID");
if (!input->variable->equalstyle(n))
error->all(FLERR,"Dump modify variable must be of style equal");
perframe[i].type = THIS_IS_A_VARIABLE;
perframe[i].dim = 1;
perframe[i].index = n;
strcpy(perframe[i].name, arg[iarg]);
strcpy(perframe[i].id, suffix);
}
else {
char errstr[1024];
sprintf(errstr, "perframe quantity '%s' must be compute, fix or "
"variable", arg[iarg]);
error->all(FLERR,errstr);
}
delete [] suffix;
}
if (iarg >= narg)
error->all(FLERR,"expected 'yes' or 'no' after 'thermo' keyword.");
if (strcmp(arg[iarg],"yes") == 0) {
thermo = true;
}
return narg;
else if (strcmp(arg[iarg],"no") == 0) {
thermo = false;
}
else error->all(FLERR,"expected 'yes' or 'no' after 'thermo' keyword.");
iarg++;
return 2;
} else return 0;
}
@ -1101,41 +994,14 @@ void DumpNetCDF::ncerr(int err, const char *descr, int line)
char errstr[1024];
if (descr) {
sprintf(errstr, "NetCDF failed with error '%s' (while accessing '%s') "
" in line %i of %s.", nc_strerror(err), descr, line, __FILE__);
" in line %i of %s.", nc_strerror(err), descr, line, __FILE__);
}
else {
sprintf(errstr, "NetCDF failed with error '%s' in line %i of %s.",
nc_strerror(err), line, __FILE__);
nc_strerror(err), line, __FILE__);
}
error->one(FLERR,errstr);
}
}
/* ----------------------------------------------------------------------
one method for every keyword thermo can output
called by compute() or evaluate_keyword()
compute will have already been called
set ivalue/dvalue/bivalue if value is int/double/bigint
customize a new keyword by adding a method
------------------------------------------------------------------------- */
void DumpNetCDF::compute_step(void *r)
{
*((bigint *) r) = update->ntimestep;
}
/* ---------------------------------------------------------------------- */
void DumpNetCDF::compute_elapsed(void *r)
{
*((bigint *) r) = update->ntimestep - update->firststep;
}
/* ---------------------------------------------------------------------- */
void DumpNetCDF::compute_elapsed_long(void *r)
{
*((bigint *) r) = update->ntimestep - update->beginstep;
}
#endif /* defined(LMP_HAS_NETCDF) */