diff --git a/src/NETCDF/dump_netcdf.cpp b/src/NETCDF/dump_netcdf.cpp index 5f30c941ca..137d6368c2 100644 --- a/src/NETCDF/dump_netcdf.cpp +++ b/src/NETCDF/dump_netcdf.cpp @@ -19,6 +19,7 @@ #if defined(LMP_HAS_NETCDF) #include "dump_netcdf.h" +#include "netcdf_units.h" #include "atom.h" #include "comm.h" @@ -43,6 +44,9 @@ using namespace LAMMPS_NS; using namespace MathConst; +using NetCDFUnits::Quantity; +using NetCDFUnits::get_unit_for; +using NetCDFUnits::LMP_MAX_VAR_DIMS; static const char NC_FRAME_STR[] = "frame"; static const char NC_SPATIAL_STR[] = "spatial"; @@ -63,7 +67,6 @@ static const char NC_SCALE_FACTOR_STR[] = "scale_factor"; static constexpr int THIS_IS_A_FIX = -1; static constexpr int THIS_IS_A_COMPUTE = -2; static constexpr int THIS_IS_A_VARIABLE = -3; -static constexpr int THIS_IS_A_BIGINT = -4; /* ---------------------------------------------------------------------- */ @@ -102,6 +105,7 @@ DumpNetCDF::DumpNetCDF(LAMMPS *lmp, int narg, char **arg) : int ndims = 1; std::string mangled = earg[i]; bool constant = false; + int quantity = Quantity::UNKNOWN; // name mangling // in the AMBER specification @@ -109,26 +113,32 @@ DumpNetCDF::DumpNetCDF(LAMMPS *lmp, int narg, char **arg) : idim = mangled[0] - 'x'; ndims = 3; mangled = "coordinates"; + quantity = Quantity::DISTANCE; } else if ((mangled == "vx") || (mangled == "vy") || (mangled == "vz")) { idim = mangled[1] - 'x'; ndims = 3; mangled = "velocities"; + quantity = Quantity::VELOCITY; } else if ((mangled == "xs") || (mangled == "ys") || (mangled == "zs")) { idim = mangled[0] - 'x'; ndims = 3; mangled = "scaled_coordinates"; + // no unit for scaled coordinates } else if ((mangled == "xu") || (mangled == "yu") || (mangled == "zu")) { idim = mangled[0] - 'x'; ndims = 3; mangled = "unwrapped_coordinates"; + quantity = Quantity::DISTANCE; } else if ((mangled == "fx") || (mangled == "fy") || (mangled == "fz")) { idim = mangled[1] - 'x'; ndims = 3; mangled = "forces"; + quantity = Quantity::FORCE; } else if ((mangled == "mux") || (mangled == "muy") || (mangled == "muz")) { idim = mangled[2] - 'x'; ndims = 3; mangled = "mu"; + quantity = Quantity::DIPOLE_MOMENT; } else if (utils::strmatch(mangled, "^c_")) { std::size_t found = mangled.find('['); if (found != std::string::npos) { @@ -175,13 +185,14 @@ DumpNetCDF::DumpNetCDF(LAMMPS *lmp, int narg, char **arg) : perat[inc].constant = constant; perat[inc].ndumped = 0; perat[inc].field[idim] = i; + perat[inc].quantity = quantity; } n_buffer = 0; int_buffer = nullptr; double_buffer = nullptr; - double_precision = false; + type_nc_real = NC_FLOAT; thermo = false; thermovar = nullptr; @@ -196,7 +207,7 @@ DumpNetCDF::~DumpNetCDF() closefile(); delete[] perat; - if (thermovar) delete[] thermovar; + delete[] thermovar; if (int_buffer) memory->sfree(int_buffer); if (double_buffer) memory->sfree(double_buffer); @@ -224,7 +235,7 @@ void DumpNetCDF::openfile() } if (thermo && !singlefile_opened) { - if (thermovar) delete[] thermovar; + delete[] thermovar; thermovar = new int[output->thermo->nfield]; } @@ -290,18 +301,18 @@ void DumpNetCDF::openfile() NCERRX( nc_inq_dimid(ncid, NC_LABEL_STR, &label_dim), NC_LABEL_STR ); for (int i = 0; i < n_perat; i++) { - int dims = perat[i].dims; - if (vector_dim[dims] < 0) { + int dim = perat[i].dims; + if (vector_dim[dim] < 0) { char dimstr[1024]; - if (dims == 3) { + if (dim == 3) { strcpy(dimstr, NC_SPATIAL_STR); - } else if (dims == 6) { + } else if (dim == 6) { strcpy(dimstr, NC_VOIGT_STR); } else { - sprintf(dimstr, "vec%i", dims); + sprintf(dimstr, "vec%i", dim); } - if (dims != 1) { - NCERRX( nc_inq_dimid(ncid, dimstr, &vector_dim[dims]), dimstr ); + if (dim != 1) { + NCERRX( nc_inq_dimid(ncid, dimstr, &vector_dim[dim]), dimstr ); } } } @@ -339,9 +350,8 @@ void DumpNetCDF::openfile() if (framei != 0 && !multifile) error->all(FLERR,"at keyword requires use of 'append yes'"); - int dims[NC_MAX_VAR_DIMS]; - size_t index[NC_MAX_VAR_DIMS], count[NC_MAX_VAR_DIMS]; - double d[1]; + int dims[LMP_MAX_VAR_DIMS]; + size_t index[LMP_MAX_VAR_DIMS], count[LMP_MAX_VAR_DIMS]; if (singlefile_opened) return; singlefile_opened = 1; @@ -373,22 +383,22 @@ void DumpNetCDF::openfile() } // default variables - dims[0] = 0; + dims[0] = vector_dim[3]; NCERRX( nc_def_var(ncid, NC_SPATIAL_STR, NC_CHAR, 1, dims, &spatial_var), NC_SPATIAL_STR ); NCERRX( nc_def_var(ncid, NC_CELL_SPATIAL_STR, NC_CHAR, 1, dims, &cell_spatial_var), NC_CELL_SPATIAL_STR ); - dims[0] = 0; + dims[0] = vector_dim[3]; dims[1] = label_dim; NCERRX( nc_def_var(ncid, NC_CELL_ANGULAR_STR, NC_CHAR, 2, dims, &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); + NCERRX( nc_def_var(ncid, NC_TIME_STR, type_nc_real, 1, dims, &time_var), 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 ); - NCERRX( nc_def_var(ncid, NC_CELL_LENGTHS_STR, NC_DOUBLE, 2, dims, &cell_lengths_var), NC_CELL_LENGTHS_STR ); + NCERRX( nc_def_var(ncid, NC_CELL_ORIGIN_STR, type_nc_real, 2, dims, &cell_origin_var), NC_CELL_ORIGIN_STR ); + NCERRX( nc_def_var(ncid, NC_CELL_LENGTHS_STR, type_nc_real, 2, dims, &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 ); + NCERRX( nc_def_var(ncid, NC_CELL_ANGLES_STR, type_nc_real, 2, dims, &cell_angles_var), NC_CELL_ANGLES_STR ); // variables specified in the input file dims[0] = frame_dim; @@ -397,7 +407,6 @@ void DumpNetCDF::openfile() for (int i = 0; i < n_perat; i++) { nc_type xtype; - // Type mangling if (vtype[perat[i].field[0]] == Dump::INT) { xtype = NC_INT; @@ -406,10 +415,7 @@ void DumpNetCDF::openfile() } else if (vtype[perat[i].field[0]] == Dump::STRING) { error->all(FLERR,"Dump netcdf currently does not support dumping string properties"); } else { - if (double_precision) - xtype = NC_DOUBLE; - else - xtype = NC_FLOAT; + xtype = type_nc_real; } if (perat[i].constant) { @@ -430,6 +436,11 @@ void DumpNetCDF::openfile() NCERRX( nc_def_var(ncid, perat[i].name, xtype, 3, dims, &perat[i].var), perat[i].name ); } } + + std::string unit = get_unit_for(update->unit_style, perat[i].quantity, error); + if (!unit.empty()) { + NCERR( nc_put_att_text(ncid, perat[i].var, NC_UNITS_STR, unit.size(), unit.c_str()) ); + } } // perframe variables @@ -437,7 +448,7 @@ void DumpNetCDF::openfile() Thermo *th = output->thermo; for (int i = 0; i < th->nfield; i++) { if (th->vtype[i] == Thermo::FLOAT) { - NCERRX( nc_def_var(ncid, th->keyword[i], NC_DOUBLE, 1, dims, + NCERRX( nc_def_var(ncid, th->keyword[i], type_nc_real, 1, dims, &thermovar[i]), th->keyword[i] ); } else if (th->vtype[i] == Thermo::INT) { NCERRX( nc_def_var(ncid, th->keyword[i], NC_INT, 1, dims, @@ -461,43 +472,18 @@ void DumpNetCDF::openfile() NCERR( nc_put_att_text(ncid, NC_GLOBAL, "program", 6, "LAMMPS") ); NCERR( nc_put_att_text(ncid, NC_GLOBAL, "programVersion",strlen(lmp->version), lmp->version) ); - // units - if (!strcmp(update->unit_style, "lj")) { - NCERR( nc_put_att_text(ncid, time_var, NC_UNITS_STR, 2, "lj") ); - NCERR( nc_put_att_text(ncid, cell_origin_var, NC_UNITS_STR, 2, "lj") ); - NCERR( nc_put_att_text(ncid, cell_lengths_var, NC_UNITS_STR, 2, "lj") ); - } else if (!strcmp(update->unit_style, "real")) { - NCERR( nc_put_att_text(ncid, time_var, NC_UNITS_STR, 11, "femtosecond") ); - NCERR( nc_put_att_text(ncid, cell_origin_var, NC_UNITS_STR, 8, "Angstrom") ); - NCERR( nc_put_att_text(ncid, cell_lengths_var, NC_UNITS_STR, 8, "Angstrom") ); - } else if (!strcmp(update->unit_style, "metal")) { - NCERR( nc_put_att_text(ncid, time_var, NC_UNITS_STR, 10, "picosecond") ); - NCERR( nc_put_att_text(ncid, cell_origin_var, NC_UNITS_STR, 8, "Angstrom") ); - NCERR( nc_put_att_text(ncid, cell_lengths_var, NC_UNITS_STR, 8, "Angstrom") ); - } else if (!strcmp(update->unit_style, "si")) { - NCERR( nc_put_att_text(ncid, time_var, NC_UNITS_STR, 6, "second") ); - NCERR( nc_put_att_text(ncid, cell_origin_var, NC_UNITS_STR, 5, "meter") ); - NCERR( nc_put_att_text(ncid, cell_lengths_var, NC_UNITS_STR, 5, "meter") ); - } else if (!strcmp(update->unit_style, "cgs")) { - NCERR( nc_put_att_text(ncid, time_var, NC_UNITS_STR, 6, "second") ); - NCERR( nc_put_att_text(ncid, cell_origin_var, NC_UNITS_STR, 10, "centimeter") ); - NCERR( nc_put_att_text(ncid, cell_lengths_var, NC_UNITS_STR, 10, "centimeter") ); - } else if (!strcmp(update->unit_style, "electron")) { - NCERR( nc_put_att_text(ncid, time_var, NC_UNITS_STR, 11, "femtosecond") ); - NCERR( nc_put_att_text(ncid, cell_origin_var, NC_UNITS_STR, 4, "Bohr") ); - NCERR( nc_put_att_text(ncid, cell_lengths_var, NC_UNITS_STR, 4, "Bohr") ); - } else { - error->all(FLERR,"Unsupported unit style: {}", update->unit_style); - } + // units & scale + std::string unit = get_unit_for(update->unit_style, Quantity::TIME, error); + NCERR( nc_put_att_text(ncid, time_var, NC_UNITS_STR, unit.size(), unit.c_str()) ); - NCERR( nc_put_att_text(ncid, cell_angles_var, NC_UNITS_STR,6, "degree") ); + unit = get_unit_for(update->unit_style, Quantity::DISTANCE, error); + NCERR( nc_put_att_text(ncid, cell_origin_var, NC_UNITS_STR, unit.size(), unit.c_str()) ); + NCERR( nc_put_att_text(ncid, cell_lengths_var, NC_UNITS_STR, unit.size(), unit.c_str()) ); - d[0] = update->dt; - NCERR( nc_put_att_double(ncid, time_var, NC_SCALE_FACTOR_STR,NC_DOUBLE, 1, d) ); - d[0] = 1.0; - NCERR( nc_put_att_double(ncid, cell_origin_var, NC_SCALE_FACTOR_STR,NC_DOUBLE, 1, d) ); - d[0] = 1.0; - NCERR( nc_put_att_double(ncid, cell_lengths_var, NC_SCALE_FACTOR_STR,NC_DOUBLE, 1, d) ); + NCERR( nc_put_att_text(ncid, cell_angles_var, NC_UNITS_STR, 6, "degree") ); + + float scale[1] = {static_cast(update->dt)}; + NCERR( nc_put_att_float(ncid, time_var, NC_SCALE_FACTOR_STR, NC_FLOAT, 1, scale) ); /* * Finished with definition @@ -735,8 +721,8 @@ void DumpNetCDF::write_header(bigint n) void DumpNetCDF::write_data(int n, double *mybuf) { - size_t start[NC_MAX_VAR_DIMS], count[NC_MAX_VAR_DIMS]; - ptrdiff_t stride[NC_MAX_VAR_DIMS]; + size_t start[LMP_MAX_VAR_DIMS], count[LMP_MAX_VAR_DIMS]; + ptrdiff_t stride[LMP_MAX_VAR_DIMS]; if (!int_buffer) { n_buffer = n; @@ -872,7 +858,12 @@ int DumpNetCDF::modify_param(int narg, char **arg) if (strcmp(arg[iarg],"double") == 0) { iarg++; if (iarg >= narg) error->all(FLERR,"expected 'yes' or 'no' after 'double' keyword."); - double_precision = utils::logical(FLERR,arg[iarg],false,lmp) == 1; + + if (utils::logical(FLERR,arg[iarg],false,lmp) == 1) + type_nc_real = NC_DOUBLE; + else + type_nc_real = NC_FLOAT; + iarg++; return 2; } else if (strcmp(arg[iarg],"at") == 0) { @@ -897,10 +888,10 @@ int DumpNetCDF::modify_param(int narg, char **arg) void DumpNetCDF::ncerr(int err, const char *descr, int line) { if (err != NC_NOERR) { - if (descr) error->one(FLERR,"NetCDF failed with error '{}' (while accessing '{}') " - " in line {} of {}.", nc_strerror(err), descr, line, __FILE__); - else error->one(FLERR,"NetCDF failed with error '{}' in line {} of {}.", - nc_strerror(err), line, __FILE__); + if (descr) error->one(__FILE__, line, "NetCDF failed with error '{}' (while accessing '{}') ", + nc_strerror(err), descr); + else error->one(__FILE__, line,"NetCDF failed with error '{}' in line {} of {}.", + nc_strerror(err)); } } diff --git a/src/NETCDF/dump_netcdf.h b/src/NETCDF/dump_netcdf.h index dd9c50873e..f3a4e81d9c 100644 --- a/src/NETCDF/dump_netcdf.h +++ b/src/NETCDF/dump_netcdf.h @@ -24,15 +24,12 @@ DumpStyle(netcdf,DumpNetCDF); #else #ifndef LMP_DUMP_NETCDF_H -#define LMP_DUMP_NETCDFC_H +#define LMP_DUMP_NETCDF_H #include "dump_custom.h" namespace LAMMPS_NS { -const int NC_FIELD_NAME_MAX = 100; -const int DUMP_NC_MAX_DIMS = 100; - class DumpNetCDF : public DumpCustom { public: DumpNetCDF(class LAMMPS *, int, char **); @@ -40,12 +37,16 @@ class DumpNetCDF : public DumpCustom { virtual void write(); private: + static constexpr int NC_FIELD_NAME_MAX = 100; + static constexpr int DUMP_NC_MAX_DIMS = 100; + // per-atoms quantities (positions, velocities, etc.) struct nc_perat_t { int dims; // number of dimensions int field[DUMP_NC_MAX_DIMS]; // field indices corresponding to the dim. char name[NC_FIELD_NAME_MAX]; // field name int var; // NetCDF variable + int quantity; // type of the quantity bool constant; // is this property per file (not per frame) int ndumped; // number of enties written for this prop. @@ -62,8 +63,8 @@ class DumpNetCDF : public DumpCustom { int *thermovar; // NetCDF variables for thermo output - bool double_precision; // write everything as double precision - bool thermo; // write thermo output to netcdf file + int type_nc_real; // netcdf type to use for real variables: float or double + bool thermo; // write thermo output to netcdf file bigint n_buffer; // size of buffer bigint *int_buffer; // buffer for passing data to netcdf diff --git a/src/NETCDF/dump_netcdf_mpiio.cpp b/src/NETCDF/dump_netcdf_mpiio.cpp index 0a76203f96..a1c9d20e61 100644 --- a/src/NETCDF/dump_netcdf_mpiio.cpp +++ b/src/NETCDF/dump_netcdf_mpiio.cpp @@ -19,6 +19,7 @@ #if defined(LMP_HAS_PNETCDF) #include "dump_netcdf_mpiio.h" +#include "netcdf_units.h" #include "atom.h" #include "comm.h" @@ -43,6 +44,9 @@ using namespace LAMMPS_NS; using namespace MathConst; +using NetCDFUnits::Quantity; +using NetCDFUnits::get_unit_for; +using NetCDFUnits::LMP_MAX_VAR_DIMS; static const char NC_FRAME_STR[] = "frame"; static const char NC_SPATIAL_STR[] = "spatial"; @@ -63,7 +67,6 @@ static const char NC_SCALE_FACTOR_STR[] = "scale_factor"; static constexpr int THIS_IS_A_FIX = -1; static constexpr int THIS_IS_A_COMPUTE = -2; static constexpr int THIS_IS_A_VARIABLE = -3; -static constexpr int THIS_IS_A_BIGINT = -4; /* ---------------------------------------------------------------------- */ @@ -101,7 +104,7 @@ DumpNetCDFMPIIO::DumpNetCDFMPIIO(LAMMPS *lmp, int narg, char **arg) : int idim = 0; int ndims = 1; std::string mangled = earg[i]; - bool constant = false; + int quantity = Quantity::UNKNOWN; // name mangling // in the AMBER specification @@ -109,26 +112,32 @@ DumpNetCDFMPIIO::DumpNetCDFMPIIO(LAMMPS *lmp, int narg, char **arg) : idim = mangled[0] - 'x'; ndims = 3; mangled = "coordinates"; + quantity = Quantity::DISTANCE; } else if ((mangled == "vx") || (mangled == "vy") || (mangled == "vz")) { idim = mangled[1] - 'x'; ndims = 3; mangled = "velocities"; + quantity = Quantity::VELOCITY; } else if ((mangled == "xs") || (mangled == "ys") || (mangled == "zs")) { idim = mangled[0] - 'x'; ndims = 3; mangled = "scaled_coordinates"; + // no unit for scaled coordinates } else if ((mangled == "xu") || (mangled == "yu") || (mangled == "zu")) { idim = mangled[0] - 'x'; ndims = 3; mangled = "unwrapped_coordinates"; + quantity = Quantity::DISTANCE; } else if ((mangled == "fx") || (mangled == "fy") || (mangled == "fz")) { idim = mangled[1] - 'x'; ndims = 3; mangled = "forces"; + quantity = Quantity::FORCE; } else if ((mangled == "mux") || (mangled == "muy") || (mangled == "muz")) { idim = mangled[2] - 'x'; ndims = 3; mangled = "mu"; + quantity = Quantity::DIPOLE_MOMENT; } else if (utils::strmatch(mangled, "^c_")) { std::size_t found = mangled.find('['); if (found != std::string::npos) { @@ -173,13 +182,14 @@ DumpNetCDFMPIIO::DumpNetCDFMPIIO(LAMMPS *lmp, int narg, char **arg) : } perat[inc].field[idim] = i; + perat[inc].quantity = quantity; } n_buffer = 0; int_buffer = nullptr; double_buffer = nullptr; - double_precision = false; + type_nc_real = NC_FLOAT; thermo = false; thermovar = nullptr; @@ -194,7 +204,7 @@ DumpNetCDFMPIIO::~DumpNetCDFMPIIO() closefile(); delete[] perat; - if (thermovar) delete[] thermovar; + delete[] thermovar; if (int_buffer) memory->sfree(int_buffer); if (double_buffer) memory->sfree(double_buffer); @@ -211,8 +221,7 @@ void DumpNetCDFMPIIO::openfile() char *ptr = strchr(filestar,'*'); *ptr = '\0'; if (padflag == 0) - sprintf(filecurrent,"%s" BIGINT_FORMAT "%s", - filestar,update->ntimestep,ptr+1); + sprintf(filecurrent,"%s" BIGINT_FORMAT "%s", filestar,update->ntimestep,ptr+1); else { char bif[8],pad[16]; strcpy(bif,BIGINT_FORMAT); @@ -223,7 +232,7 @@ void DumpNetCDFMPIIO::openfile() } if (thermo && !singlefile_opened) { - if (thermovar) delete[] thermovar; + delete[] thermovar; thermovar = new int[output->thermo->nfield]; } @@ -275,9 +284,6 @@ void DumpNetCDFMPIIO::openfile() if (!platform::file_is_readable(filecurrent)) error->all(FLERR, "cannot append to non-existent file {}", filecurrent); - MPI_Offset index[NC_MAX_VAR_DIMS], count[NC_MAX_VAR_DIMS]; - double d[1]; - if (singlefile_opened) return; singlefile_opened = 1; @@ -291,18 +297,18 @@ void DumpNetCDFMPIIO::openfile() NCERRX( ncmpi_inq_dimid(ncid, NC_LABEL_STR, &label_dim), NC_LABEL_STR ); for (int i = 0; i < n_perat; i++) { - int dims = perat[i].dims; - if (vector_dim[dims] < 0) { + int dim = perat[i].dims; + if (vector_dim[dim] < 0) { char dimstr[1024]; - if (dims == 3) { + if (dim == 3) { strcpy(dimstr, NC_SPATIAL_STR); - } else if (dims == 6) { + } else if (dim == 6) { strcpy(dimstr, NC_VOIGT_STR); } else { - sprintf(dimstr, "vec%i", dims); + sprintf(dimstr, "vec%i", dim); } - if (dims != 1) { - NCERRX( ncmpi_inq_dimid(ncid, dimstr, &vector_dim[dims]), dimstr ); + if (dim != 1) { + NCERRX( ncmpi_inq_dimid(ncid, dimstr, &vector_dim[dim]), dimstr ); } } } @@ -340,9 +346,8 @@ void DumpNetCDFMPIIO::openfile() if (framei != 0 && !multifile) error->all(FLERR,"at keyword requires use of 'append yes'"); - int dims[NC_MAX_VAR_DIMS]; - MPI_Offset index[NC_MAX_VAR_DIMS], count[NC_MAX_VAR_DIMS]; - double d[1]; + int dims[LMP_MAX_VAR_DIMS]; + MPI_Offset index[LMP_MAX_VAR_DIMS], count[LMP_MAX_VAR_DIMS]; if (singlefile_opened) return; singlefile_opened = 1; @@ -356,19 +361,24 @@ void DumpNetCDFMPIIO::openfile() NCERRX( ncmpi_def_dim(ncid, NC_CELL_ANGULAR_STR, 3, &cell_angular_dim), NC_CELL_ANGULAR_STR ); NCERRX( ncmpi_def_dim(ncid, NC_LABEL_STR, 10, &label_dim), NC_LABEL_STR ); + if (vector_dim[3] < 0) + NCERRX( ncmpi_def_dim(ncid, NC_SPATIAL_STR, 3, &vector_dim[3]), NC_SPATIAL_STR ); + if (vector_dim[6] < 0) + NCERRX( ncmpi_def_dim(ncid, NC_VOIGT_STR, 6, &vector_dim[6]), NC_VOIGT_STR ); + for (int i = 0; i < n_perat; i++) { - int dims = perat[i].dims; - if (vector_dim[dims] < 0) { + int dim = perat[i].dims; + if (vector_dim[dim] < 0) { char dimstr[1024]; - if (dims == 3) { + if (dim == 3) { strcpy(dimstr, NC_SPATIAL_STR); - } else if (dims == 6) { + } else if (dim == 6) { strcpy(dimstr, NC_VOIGT_STR); } else { - sprintf(dimstr, "vec%i", dims); + sprintf(dimstr, "vec%i", dim); } - if (dims != 1) { - NCERRX( ncmpi_def_dim(ncid, dimstr, dims, &vector_dim[dims]), dimstr ); + if (dim != 1) { + NCERRX( ncmpi_def_dim(ncid, dimstr, dim, &vector_dim[dim]), dimstr ); } } } @@ -380,16 +390,15 @@ void DumpNetCDFMPIIO::openfile() dims[0] = vector_dim[3]; dims[1] = label_dim; NCERRX( ncmpi_def_var(ncid, NC_CELL_ANGULAR_STR, NC_CHAR, 2, dims, &cell_angular_var), NC_CELL_ANGULAR_STR ); - dims[0] = frame_dim; - NCERRX( ncmpi_def_var(ncid, NC_TIME_STR, NC_DOUBLE, 1, dims, &time_var), NC_TIME_STR); + NCERRX( ncmpi_def_var(ncid, NC_TIME_STR, type_nc_real, 1, dims, &time_var), NC_TIME_STR); dims[0] = frame_dim; dims[1] = cell_spatial_dim; - NCERRX( ncmpi_def_var(ncid, NC_CELL_ORIGIN_STR, NC_DOUBLE, 2, dims, &cell_origin_var), NC_CELL_ORIGIN_STR ); - NCERRX( ncmpi_def_var(ncid, NC_CELL_LENGTHS_STR, NC_DOUBLE, 2, dims, &cell_lengths_var), NC_CELL_LENGTHS_STR ); + NCERRX( ncmpi_def_var(ncid, NC_CELL_ORIGIN_STR, type_nc_real, 2, dims, &cell_origin_var), NC_CELL_ORIGIN_STR ); + NCERRX( ncmpi_def_var(ncid, NC_CELL_LENGTHS_STR, type_nc_real, 2, dims, &cell_lengths_var), NC_CELL_LENGTHS_STR ); dims[0] = frame_dim; dims[1] = cell_angular_dim; - NCERRX( ncmpi_def_var(ncid, NC_CELL_ANGLES_STR, NC_DOUBLE, 2, dims, &cell_angles_var), NC_CELL_ANGLES_STR ); + NCERRX( ncmpi_def_var(ncid, NC_CELL_ANGLES_STR, type_nc_real, 2, dims, &cell_angles_var), NC_CELL_ANGLES_STR ); // variables specified in the input file dims[0] = frame_dim; @@ -405,10 +414,7 @@ void DumpNetCDFMPIIO::openfile() } else if (vtype[perat[i].field[0]] == Dump::BIGINT) { xtype = NC_INT64; } else { - if (double_precision) - xtype = NC_DOUBLE; - else - xtype = NC_FLOAT; + xtype = type_nc_real; } if (perat[i].dims == 1) { @@ -418,6 +424,11 @@ void DumpNetCDFMPIIO::openfile() dims[2] = vector_dim[perat[i].dims]; NCERRX( ncmpi_def_var(ncid, perat[i].name, xtype, 3, dims, &perat[i].var), perat[i].name ); } + + std::string unit = get_unit_for(update->unit_style, perat[i].quantity, error); + if (!unit.empty()) { + NCERR( ncmpi_put_att_text(ncid, perat[i].var, NC_UNITS_STR, unit.size(), unit.c_str()) ); + } } // perframe variables @@ -425,7 +436,7 @@ void DumpNetCDFMPIIO::openfile() Thermo *th = output->thermo; for (int i = 0; i < th->nfield; i++) { if (th->vtype[i] == Thermo::FLOAT) { - NCERRX( ncmpi_def_var(ncid, th->keyword[i], NC_DOUBLE, 1, dims, &thermovar[i]), th->keyword[i] ); + NCERRX( ncmpi_def_var(ncid, th->keyword[i], type_nc_real, 1, dims, &thermovar[i]), th->keyword[i] ); } else if (th->vtype[i] == Thermo::INT) { NCERRX( ncmpi_def_var(ncid, th->keyword[i], NC_INT, 1, dims, &thermovar[i]), th->keyword[i] ); } else if (th->vtype[i] == Thermo::BIGINT) { @@ -445,43 +456,18 @@ void DumpNetCDFMPIIO::openfile() NCERR( ncmpi_put_att_text(ncid, NC_GLOBAL, "program", 6, "LAMMPS") ); NCERR( ncmpi_put_att_text(ncid, NC_GLOBAL, "programVersion", strlen(lmp->version), lmp->version) ); - // units - if (!strcmp(update->unit_style, "lj")) { - NCERR( ncmpi_put_att_text(ncid, time_var, NC_UNITS_STR, 2, "lj") ); - NCERR( ncmpi_put_att_text(ncid, cell_origin_var, NC_UNITS_STR, 2, "lj") ); - NCERR( ncmpi_put_att_text(ncid, cell_lengths_var, NC_UNITS_STR, 2, "lj") ); - } else if (!strcmp(update->unit_style, "real")) { - NCERR( ncmpi_put_att_text(ncid, time_var, NC_UNITS_STR, 11, "femtosecond") ); - NCERR( ncmpi_put_att_text(ncid, cell_origin_var, NC_UNITS_STR, 8, "Angstrom") ); - NCERR( ncmpi_put_att_text(ncid, cell_lengths_var, NC_UNITS_STR, 8, "Angstrom") ); - } else if (!strcmp(update->unit_style, "metal")) { - NCERR( ncmpi_put_att_text(ncid, time_var, NC_UNITS_STR, 10, "picosecond") ); - NCERR( ncmpi_put_att_text(ncid, cell_origin_var, NC_UNITS_STR, 8, "Angstrom") ); - NCERR( ncmpi_put_att_text(ncid, cell_lengths_var, NC_UNITS_STR, 8, "Angstrom") ); - } else if (!strcmp(update->unit_style, "si")) { - NCERR( ncmpi_put_att_text(ncid, time_var, NC_UNITS_STR, 6, "second") ); - NCERR( ncmpi_put_att_text(ncid, cell_origin_var, NC_UNITS_STR, 5, "meter") ); - NCERR( ncmpi_put_att_text(ncid, cell_lengths_var, NC_UNITS_STR, 5, "meter") ); - } else if (!strcmp(update->unit_style, "cgs")) { - NCERR( ncmpi_put_att_text(ncid, time_var, NC_UNITS_STR, 6, "second") ); - NCERR( ncmpi_put_att_text(ncid, cell_origin_var, NC_UNITS_STR, 10, "centimeter") ); - NCERR( ncmpi_put_att_text(ncid, cell_lengths_var, NC_UNITS_STR, 10, "centimeter") ); - } else if (!strcmp(update->unit_style, "electron")) { - NCERR( ncmpi_put_att_text(ncid, time_var, NC_UNITS_STR, 11, "femtosecond") ); - NCERR( ncmpi_put_att_text(ncid, cell_origin_var, NC_UNITS_STR, 4, "Bohr") ); - NCERR( ncmpi_put_att_text(ncid, cell_lengths_var, NC_UNITS_STR, 4, "Bohr") ); - } else { - error->all(FLERR,"Unsupported unit style: {}", update->unit_style); - } + // units & scale + std::string unit = get_unit_for(update->unit_style, Quantity::TIME, error); + NCERR( ncmpi_put_att_text(ncid, time_var, NC_UNITS_STR, unit.size(), unit.c_str()) ); + + unit = get_unit_for(update->unit_style, Quantity::DISTANCE, error); + NCERR( ncmpi_put_att_text(ncid, cell_origin_var, NC_UNITS_STR, unit.size(), unit.c_str()) ); + NCERR( ncmpi_put_att_text(ncid, cell_lengths_var, NC_UNITS_STR, unit.size(), unit.c_str()) ); NCERR( ncmpi_put_att_text(ncid, cell_angles_var, NC_UNITS_STR, 6, "degree") ); - d[0] = update->dt; - NCERR( ncmpi_put_att_double(ncid, time_var, NC_SCALE_FACTOR_STR, NC_DOUBLE, 1, d) ); - d[0] = 1.0; - NCERR( ncmpi_put_att_double(ncid, cell_origin_var, NC_SCALE_FACTOR_STR, NC_DOUBLE, 1, d) ); - d[0] = 1.0; - NCERR( ncmpi_put_att_double(ncid, cell_lengths_var, NC_SCALE_FACTOR_STR, NC_DOUBLE, 1, d) ); + float scale[1] = {static_cast(update->dt)}; + NCERR( ncmpi_put_att_float(ncid, time_var, NC_SCALE_FACTOR_STR, NC_FLOAT, 1, scale) ); /* * Finished with definition @@ -502,16 +488,13 @@ void DumpNetCDFMPIIO::openfile() index[1] = 0; count[0] = 1; count[1] = 5; - NCERR( ncmpi_put_vara_text(ncid, cell_angular_var, index, count, - "alpha") ); + NCERR( ncmpi_put_vara_text(ncid, cell_angular_var, index, count, "alpha") ); index[0] = 1; count[1] = 4; - NCERR( ncmpi_put_vara_text(ncid, cell_angular_var, index, count, - "beta") ); + NCERR( ncmpi_put_vara_text(ncid, cell_angular_var, index, count, "beta") ); index[0] = 2; count[1] = 5; - NCERR( ncmpi_put_vara_text(ncid, cell_angular_var, index, count, - "gamma") ); + NCERR( ncmpi_put_vara_text(ncid, cell_angular_var, index, count, "gamma") ); } NCERR( ncmpi_end_indep_data(ncid) ); @@ -753,8 +736,7 @@ void DumpNetCDFMPIIO::write_time_and_cell() void DumpNetCDFMPIIO::write_data(int n, double *mybuf) { - MPI_Offset start[NC_MAX_VAR_DIMS], count[NC_MAX_VAR_DIMS]; - MPI_Offset stride[NC_MAX_VAR_DIMS]; + MPI_Offset start[LMP_MAX_VAR_DIMS], count[LMP_MAX_VAR_DIMS], stride[LMP_MAX_VAR_DIMS]; if (!int_buffer) { n_buffer = std::max(1, n); @@ -867,7 +849,12 @@ int DumpNetCDFMPIIO::modify_param(int narg, char **arg) if (strcmp(arg[iarg],"double") == 0) { iarg++; if (iarg >= narg) error->all(FLERR,"expected 'yes' or 'no' after 'double' keyword."); - double_precision = utils::logical(FLERR,arg[iarg],false,lmp) == 1; + + if (utils::logical(FLERR,arg[iarg],false,lmp) == 1) + type_nc_real = NC_DOUBLE; + else + type_nc_real = NC_FLOAT; + iarg++; return 2; } else if (strcmp(arg[iarg],"at") == 0) { @@ -892,10 +879,9 @@ int DumpNetCDFMPIIO::modify_param(int narg, char **arg) void DumpNetCDFMPIIO::ncerr(int err, const char *descr, int line) { if (err != NC_NOERR) { - if (descr) error->one(FLERR,"NetCDF failed with error '{}' (while accessing '{}') " - " in line {} of {}.", ncmpi_strerror(err), descr, line, __FILE__); - else error->one(FLERR,"NetCDF failed with error '{}' in line {} of {}.", - ncmpi_strerror(err), line, __FILE__); + if (descr) error->one(__FILE__, line, "NetCDF failed with error '{}' (while accessing '{}') ", + ncmpi_strerror(err), descr); + else error->one(__FILE__, line,"NetCDF failed with error '{}'.", ncmpi_strerror(err)); } } diff --git a/src/NETCDF/dump_netcdf_mpiio.h b/src/NETCDF/dump_netcdf_mpiio.h index 56c07fc3d3..ec6cbaec04 100644 --- a/src/NETCDF/dump_netcdf_mpiio.h +++ b/src/NETCDF/dump_netcdf_mpiio.h @@ -30,9 +30,6 @@ DumpStyle(netcdf/mpiio,DumpNetCDFMPIIO); namespace LAMMPS_NS { -const int NC_MPIIO_FIELD_NAME_MAX = 100; -const int DUMP_NC_MPIIO_MAX_DIMS = 100; - class DumpNetCDFMPIIO : public DumpCustom { public: DumpNetCDFMPIIO(class LAMMPS *, int, char **); @@ -40,16 +37,18 @@ class DumpNetCDFMPIIO : public DumpCustom { virtual void write(); private: + static constexpr int NC_MPIIO_FIELD_NAME_MAX = 100; + static constexpr int DUMP_NC_MPIIO_MAX_DIMS = 100; + // per-atoms quantities (positions, velocities, etc.) struct nc_perat_t { int dims; // number of dimensions int field[DUMP_NC_MPIIO_MAX_DIMS]; // field indices corresponding to the dim. char name[NC_MPIIO_FIELD_NAME_MAX]; // field name int var; // NetCDF variable + int quantity; // type of the quantity }; - typedef void (DumpNetCDFMPIIO::*funcptr_t)(void *); - int framei; // current frame index int blocki; // current block index int ndata; // number of data blocks to expect @@ -61,8 +60,8 @@ class DumpNetCDFMPIIO : public DumpCustom { int *thermovar; // NetCDF variables for thermo output - bool double_precision; // write everything as double precision - bool thermo; // write thermo output to netcdf file + int type_nc_real; // netcdf type to use for real variables: float or double + bool thermo; // write thermo output to netcdf file bigint n_buffer; // size of buffer bigint *int_buffer; // buffer for passing data to netcdf diff --git a/src/NETCDF/netcdf_units.cpp b/src/NETCDF/netcdf_units.cpp new file mode 100644 index 0000000000..0ee0ebbde0 --- /dev/null +++ b/src/NETCDF/netcdf_units.cpp @@ -0,0 +1,145 @@ +/* ---------------------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + https://www.lammps.org/, Sandia National Laboratories + Steve Plimpton, sjplimp@sandia.gov + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +/* ---------------------------------------------------------------------- + Contributing author: Lars Pastewka (University of Freiburg) +------------------------------------------------------------------------- */ + +#if defined(LMP_HAS_NETCDF) || defined(LMP_HAS_PNETCDF) + +#include "netcdf_units.h" + +#include "error.h" + +using namespace LAMMPS_NS; + +std::string NetCDFUnits::get_unit_for(const char *unit_style, int quantity, Error *error) +{ + if (!strcmp(unit_style, "lj")) { + if (quantity == Quantity::UNKNOWN) { + return ""; + } else { + return "lj"; + } + } else if (!strcmp(unit_style, "real")) { + switch (quantity) { + case Quantity::UNKNOWN: + return ""; + case Quantity::TIME: + return "femtosecond"; + case Quantity::DISTANCE: + return "angstrom"; + case Quantity::VELOCITY: + return "angstrom/femtosecond"; + case Quantity::FORCE: + return "(Kcal/mol)/angstrom)"; + case Quantity::DIPOLE_MOMENT: + return "e * angstrom"; + } + } else if (!strcmp(unit_style, "metal")) { + switch (quantity) { + case Quantity::UNKNOWN: + return ""; + case Quantity::TIME: + return "picosecond"; + case Quantity::DISTANCE: + return "angstrom"; + case Quantity::VELOCITY: + return "angstrom/picosecond"; + case Quantity::FORCE: + return "eV/angstrom"; + case Quantity::DIPOLE_MOMENT: + return "e * angstrom"; + } + } else if (!strcmp(unit_style, "si")) { + switch (quantity) { + case Quantity::UNKNOWN: + return ""; + case Quantity::TIME: + return "second"; + case Quantity::DISTANCE: + return "meter"; + case Quantity::VELOCITY: + return "meter/second"; + case Quantity::FORCE: + return "Newton"; + case Quantity::DIPOLE_MOMENT: + return "Coulomb * meter"; + } + } else if (!strcmp(unit_style, "cgs")) { + switch (quantity) { + case Quantity::UNKNOWN: + return ""; + case Quantity::TIME: + return "second"; + case Quantity::DISTANCE: + return "centimeter"; + case Quantity::VELOCITY: + return "centimeter/second"; + case Quantity::FORCE: + return "dynes"; + case Quantity::DIPOLE_MOMENT: + return "statcoul * cm"; + } + } else if (!strcmp(unit_style, "electron")) { + switch (quantity) { + case Quantity::UNKNOWN: + return ""; + case Quantity::TIME: + return "femtoseconds"; + case Quantity::DISTANCE: + return "Bohr"; + case Quantity::VELOCITY: + return "Bohr/atomic time units"; + case Quantity::FORCE: + return "Hartree/Bohr"; + case Quantity::DIPOLE_MOMENT: + return "Debye"; + } + } else if (!strcmp(unit_style, "micro")) { + switch (quantity) { + case Quantity::UNKNOWN: + return ""; + case Quantity::TIME: + return "microseconds"; + case Quantity::DISTANCE: + return "micrometers"; + case Quantity::VELOCITY: + return "micrometers/microsecond"; + case Quantity::FORCE: + return "picogram * micrometer/microsecond^2"; + case Quantity::DIPOLE_MOMENT: + return "picocoulomb * micrometer"; + } + } else if (!strcmp(unit_style, "nano")) { + switch (quantity) { + case Quantity::UNKNOWN: + return ""; + case Quantity::TIME: + return "nanoseconds"; + case Quantity::DISTANCE: + return "nanometers"; + case Quantity::VELOCITY: + return "nanometers/nanosecond"; + case Quantity::FORCE: + return "attogram * nanometer/nanosecond^2"; + case Quantity::DIPOLE_MOMENT: + return "e * nanometer"; + } + } + + error->all(FLERR, "Unsupported unit style: {}", unit_style); + return ""; +} + +#endif diff --git a/src/NETCDF/netcdf_units.h b/src/NETCDF/netcdf_units.h new file mode 100644 index 0000000000..85f9b05888 --- /dev/null +++ b/src/NETCDF/netcdf_units.h @@ -0,0 +1,49 @@ +/* -*- c++ -*- ---------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + https://www.lammps.org/, Sandia National Laboratories + Steve Plimpton, sjplimp@sandia.gov + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +/* ---------------------------------------------------------------------- + Contributing authors: Lars Pastewka (University of Freiburg), Guillaume Fraux (EPFL) +------------------------------------------------------------------------- */ + +#ifndef LMP_NETCDF_UNITS_H +#define LMP_NETCDF_UNITS_H + +#if defined(LMP_HAS_NETCDF) || defined(LMP_HAS_PNETCDF) + +#include + +namespace LAMMPS_NS { +class Error; + +namespace NetCDFUnits { + // type of quantity for per-atom values (used to get the unit) + enum Quantity { + UNKNOWN = 0, + TIME, + DISTANCE, + VELOCITY, + FORCE, + DIPOLE_MOMENT, + }; + + // for compatibility with older NetCDF versions + static constexpr int LMP_MAX_VAR_DIMS = 1024; + + // get the name of the unit for the given `quantity` in the given LAMMPS + // `unit_style` any error will be reported through `error` + std::string get_unit_for(const char *unit_style, int quantity, Error *error); +} // namespace NetCDFUnits +} // namespace LAMMPS_NS + +#endif +#endif diff --git a/unittest/formats/CMakeLists.txt b/unittest/formats/CMakeLists.txt index a17707fdf1..be8e055adb 100644 --- a/unittest/formats/CMakeLists.txt +++ b/unittest/formats/CMakeLists.txt @@ -116,6 +116,16 @@ target_link_libraries(test_dump_local PRIVATE lammps GTest::GMock) add_test(NAME DumpLocal COMMAND test_dump_local WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR}) set_tests_properties(DumpLocal PROPERTIES ENVIRONMENT "LAMMPS_POTENTIALS=${LAMMPS_POTENTIALS_DIR}") +if(PKG_NETCDF) + find_program(NCDUMP NAMES ncdump ncdump.exe) + add_executable(test_dump_netcdf test_dump_netcdf.cpp) + target_link_libraries(test_dump_netcdf PRIVATE lammps GTest::GMock) + add_test(NAME DumpNetCDF COMMAND test_dump_netcdf WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR}) + if(NOT (NCDUMP STREQUAL "NCDUMP-NOTFOUND")) + set_tests_properties(DumpNetCDF PROPERTIES ENVIRONMENT "LAMMPS_POTENTIALS=${LAMMPS_POTENTIALS_DIR};NCDUMP_BINARY=${NCDUMP}") + endif() +endif() + if(BUILD_TOOLS) set_tests_properties(DumpAtom PROPERTIES ENVIRONMENT "BINARY2TXT_BINARY=$") set_tests_properties(DumpCustom PROPERTIES ENVIRONMENT "BINARY2TXT_BINARY=$") diff --git a/unittest/formats/test_dump_netcdf.cpp b/unittest/formats/test_dump_netcdf.cpp new file mode 100644 index 0000000000..8b4110b352 --- /dev/null +++ b/unittest/formats/test_dump_netcdf.cpp @@ -0,0 +1,412 @@ +/* ---------------------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + https://www.lammps.org/, Sandia National Laboratories + Steve Plimpton, sjplimp@sandia.gov + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +#include "../testing/core.h" +#include "../testing/systems/melt.h" +#include "../testing/utils.h" +#include "fmt/format.h" +#include "library.h" +#include "output.h" +#include "thermo.h" +#include "utils.h" +#include "version.h" +#include "gmock/gmock.h" +#include "gtest/gtest.h" + +#include +#include + +using ::testing::Eq; + +char *NCDUMP_BINARY = nullptr; +bool verbose = false; + +class DumpNetCDFTest : public MeltTest { + std::string dump_style = "netcdf"; + +public: + void set_style(const std::string &new_style) { dump_style = new_style; } + + void enable_triclinic() + { + BEGIN_HIDE_OUTPUT(); + command("change_box all triclinic"); + END_HIDE_OUTPUT(); + } + + std::string dump_filename(std::string ident) + { + return fmt::format("dump_{}_{}.nc", dump_style, ident); + } + + void generate_dump(std::string dump_file, std::string fields, std::string dump_modify_options, + int ntimesteps) + { + BEGIN_HIDE_OUTPUT(); + command(fmt::format("dump id all {} 1 {} {}", dump_style, dump_file, fields)); + + if (!dump_modify_options.empty()) { + command(fmt::format("dump_modify id {}", dump_modify_options)); + } + + command(fmt::format("run {} post no", ntimesteps)); + END_HIDE_OUTPUT(); + } + + void continue_dump(int ntimesteps) + { + BEGIN_HIDE_OUTPUT(); + command(fmt::format("run {} pre no post no", ntimesteps)); + END_HIDE_OUTPUT(); + } + + void close_dump() + { + BEGIN_HIDE_OUTPUT(); + command("undump id"); + END_HIDE_OUTPUT(); + } + + std::string convert_binary_to_text(std::string binary_file) + { + BEGIN_HIDE_OUTPUT(); + std::string cmdline = fmt::format("{0} {1} > {1}.txt", NCDUMP_BINARY, binary_file); + system(cmdline.c_str()); + END_HIDE_OUTPUT(); + return fmt::format("{}.txt", binary_file); + } +}; + +TEST_F(DumpNetCDFTest, run0_plain) +{ + if (!lammps_has_style(lmp, "dump", "netcdf")) GTEST_SKIP(); + auto dump_file = dump_filename("run0"); + auto fields = "id type proc procp1 mass x y z ix iy iz xu yu zu vx vy vz fx fy fz"; + set_style("netcdf"); + generate_dump(dump_file, fields, "", 0); + + ASSERT_FILE_EXISTS(dump_file); + if (NCDUMP_BINARY) { + auto converted_file = convert_binary_to_text(dump_file); + auto lines = read_lines(converted_file); + auto header = utils::split_words(lines[0]); + ASSERT_EQ(lines.size(), 233); + ASSERT_THAT(header[0], Eq("netcdf")); + ASSERT_THAT(header[1] + ".nc", Eq(dump_file)); + + // check dimensions section + auto section = std::find(lines.begin(), lines.end(), "dimensions:"); + for (auto line = ++section; line < lines.end(); ++line) { + auto words = utils::split_words(*line); + if ((words.size() < 1) || (words[0] == "variables:")) break; + if (words[0] == "atom") ASSERT_THAT(words[2], Eq("32")); + if (words[0] == "label") ASSERT_THAT(words[2], Eq("10")); + if (words[0] == "Voigt") ASSERT_THAT(words[2], Eq("6")); + if (words[0] == "spatial") ASSERT_THAT(words[2], Eq("3")); + } + + // check variables section + section = std::find(lines.begin(), lines.end(), "variables:"); + for (auto line = ++section; line < lines.end(); ++line) { + auto words = utils::split_words(*line); + if ((words.size() < 2) || (words[0] == "data:")) break; + if (words[0] == "time:units") ASSERT_THAT(words[2], Eq("lj")); + if (words[0] == "time:scale_factor") ASSERT_THAT(words[2], Eq("0.005f")); + if (words[0] == "cell_origin:units") ASSERT_THAT(words[2], Eq("lj")); + if (words[0] == "cell_angles:units") ASSERT_THAT(words[2], Eq("degree")); + if (words[1] == "id(frame,") ASSERT_THAT(words[2], Eq("atom)")); + if (words[1] == "type(frame,") ASSERT_THAT(words[2], Eq("atom)")); + if (words[1] == "proc(frame,") ASSERT_THAT(words[2], Eq("atom)")); + if (words[1] == "procp1(frame,") ASSERT_THAT(words[2], Eq("atom)")); + if (words[1] == "mass(frame,") ASSERT_THAT(words[2], Eq("atom)")); + if (words[1] == "ix(frame,") ASSERT_THAT(words[2], Eq("atom)")); + if (words[1] == "iy(frame,") ASSERT_THAT(words[2], Eq("atom)")); + if (words[1] == "iz(frame,") ASSERT_THAT(words[2], Eq("atom)")); + if (words[0] == ":Conventions") ASSERT_THAT(words[2], Eq("AMBER")); + if (words[0] == ":ConventionVersion") ASSERT_THAT(words[2], Eq("1.0")); + if (words[0] == ":program") ASSERT_THAT(words[2], Eq("LAMMPS")); + if (words[0] == ":programVersion") ASSERT_THAT(words[2], Eq(LAMMPS_VERSION)); + } + + // check data section + section = std::find(lines.begin(), lines.end(), "data:"); + for (auto line = ++section; line < lines.end(); ++line) { + auto words = utils::split_words(*line); + if (words.size() > 0) { + if (words[0] == "spatial") ASSERT_THAT(words[2], Eq("xyz")); + if (words[0] == "cell_spatial") ASSERT_THAT(words[2], Eq("abc")); + if (words[0] == "cell_origin") { + ++line; + words = utils::split_words(*line); + ASSERT_THAT(words[0], Eq("0,")); + ASSERT_THAT(words[1], Eq("0,")); + ASSERT_THAT(words[2], Eq("0")); + } + if (words[0] == "cell_lengths") { + ++line; + words = utils::split_words(*line); + ASSERT_THAT(words[0], Eq("3.359192,")); + ASSERT_THAT(words[1], Eq("3.359192,")); + ASSERT_THAT(words[2], Eq("3.359192")); + } + if (words[0] == "cell_angles") { + ++line; + words = utils::split_words(*line); + ASSERT_THAT(words[0], Eq("90,")); + ASSERT_THAT(words[1], Eq("90,")); + ASSERT_THAT(words[2], Eq("90")); + } + if (words[0] == "id") { + ++line; + words = utils::split_words(*line); + ASSERT_THAT(words[0], Eq("1,")); + ASSERT_THAT(words[1], Eq("2,")); + ASSERT_THAT(words[2], Eq("3,")); + ASSERT_THAT(words[3], Eq("4,")); + ASSERT_THAT(words[4], Eq("5,")); + ASSERT_THAT(words[5], Eq("6,")); + ASSERT_THAT(words[6], Eq("7,")); + ASSERT_THAT(words[7], Eq("8,")); + ASSERT_THAT(words[8], Eq("9,")); + ASSERT_THAT(words[9], Eq("10,")); + } + if (words[0] == "mass") { + ++line; + words = utils::split_words(*line); + ASSERT_THAT(words[0], Eq("1,")); + ASSERT_THAT(words[1], Eq("1,")); + ASSERT_THAT(words[2], Eq("1,")); + ASSERT_THAT(words[3], Eq("1,")); + ASSERT_THAT(words[4], Eq("1,")); + ASSERT_THAT(words[5], Eq("1,")); + ASSERT_THAT(words[6], Eq("1,")); + ASSERT_THAT(words[7], Eq("1,")); + ASSERT_THAT(words[8], Eq("1,")); + ASSERT_THAT(words[9], Eq("1,")); + } + if (words[0] == "coordinates") { + ++line; + words = utils::split_words(*line); + ASSERT_THAT(words[0], Eq("0,")); + ASSERT_THAT(words[1], Eq("0,")); + ASSERT_THAT(words[2], Eq("0,")); + ++line; + words = utils::split_words(*line); + ASSERT_THAT(words[0], Eq("0.8397981,")); + ASSERT_THAT(words[1], Eq("0.8397981,")); + ASSERT_THAT(words[2], Eq("0,")); + ++line; + words = utils::split_words(*line); + ASSERT_THAT(words[0], Eq("0.8397981,")); + ASSERT_THAT(words[1], Eq("0,")); + ASSERT_THAT(words[2], Eq("0.8397981,")); + ++line; + words = utils::split_words(*line); + ASSERT_THAT(words[0], Eq("0,")); + ASSERT_THAT(words[1], Eq("0.8397981,")); + ASSERT_THAT(words[2], Eq("0.8397981,")); + } + if (words[0] == "ix") { + ++line; + words = utils::split_words(*line); + ASSERT_THAT(words[0], Eq("0,")); + ASSERT_THAT(words[1], Eq("0,")); + ASSERT_THAT(words[2], Eq("0,")); + ASSERT_THAT(words[3], Eq("0,")); + ASSERT_THAT(words[4], Eq("0,")); + ASSERT_THAT(words[5], Eq("0,")); + ASSERT_THAT(words[6], Eq("0,")); + ASSERT_THAT(words[7], Eq("0,")); + ASSERT_THAT(words[8], Eq("0,")); + ASSERT_THAT(words[9], Eq("0,")); + } + } + } + delete_file(converted_file); + } + delete_file(dump_file); +} + +TEST_F(DumpNetCDFTest, run0_mpi) +{ + if (!lammps_has_style(lmp, "dump", "netcdf/mpiio")) GTEST_SKIP(); + auto dump_file = dump_filename("mpi0"); + auto fields = "id type proc procp1 mass x y z ix iy iz xu yu zu vx vy vz fx fy fz"; + set_style("netcdf/mpiio"); + generate_dump(dump_file, fields, "", 0); + + ASSERT_FILE_EXISTS(dump_file); + if (NCDUMP_BINARY) { + auto converted_file = convert_binary_to_text(dump_file); + auto lines = read_lines(converted_file); + auto header = utils::split_words(lines[0]); + ASSERT_EQ(lines.size(), 234); + ASSERT_THAT(header[0], Eq("netcdf")); + ASSERT_THAT(header[1] + ".nc", Eq(dump_file)); + + // check dimensions section + auto section = std::find(lines.begin(), lines.end(), "dimensions:"); + for (auto line = ++section; line < lines.end(); ++line) { + auto words = utils::split_words(*line); + if ((words.size() < 1) || (words[0] == "variables:")) break; + if (words[0] == "atom") ASSERT_THAT(words[2], Eq("32")); + if (words[0] == "label") ASSERT_THAT(words[2], Eq("10")); + if (words[0] == "Voigt") ASSERT_THAT(words[2], Eq("6")); + if (words[0] == "spatial") ASSERT_THAT(words[2], Eq("3")); + } + + // check variables section + section = std::find(lines.begin(), lines.end(), "variables:"); + for (auto line = ++section; line < lines.end(); ++line) { + auto words = utils::split_words(*line); + if ((words.size() < 2) || (words[0] == "data:")) break; + if (words[0] == "time:units") ASSERT_THAT(words[2], Eq("lj")); + if (words[0] == "time:scale_factor") ASSERT_THAT(words[2], Eq("0.005f")); + if (words[0] == "cell_origin:units") ASSERT_THAT(words[2], Eq("lj")); + if (words[0] == "cell_angles:units") ASSERT_THAT(words[2], Eq("degree")); + if (words[1] == "id(frame,") ASSERT_THAT(words[2], Eq("atom)")); + if (words[1] == "type(frame,") ASSERT_THAT(words[2], Eq("atom)")); + if (words[1] == "proc(frame,") ASSERT_THAT(words[2], Eq("atom)")); + if (words[1] == "procp1(frame,") ASSERT_THAT(words[2], Eq("atom)")); + if (words[1] == "mass(frame,") ASSERT_THAT(words[2], Eq("atom)")); + if (words[1] == "ix(frame,") ASSERT_THAT(words[2], Eq("atom)")); + if (words[1] == "iy(frame,") ASSERT_THAT(words[2], Eq("atom)")); + if (words[1] == "iz(frame,") ASSERT_THAT(words[2], Eq("atom)")); + if (words[0] == ":Conventions") ASSERT_THAT(words[2], Eq("AMBER")); + if (words[0] == ":ConventionVersion") ASSERT_THAT(words[2], Eq("1.0")); + if (words[0] == ":program") ASSERT_THAT(words[2], Eq("LAMMPS")); + if (words[0] == ":programVersion") ASSERT_THAT(words[2], Eq(LAMMPS_VERSION)); + } + + // check data section + section = std::find(lines.begin(), lines.end(), "data:"); + for (auto line = ++section; line < lines.end(); ++line) { + auto words = utils::split_words(*line); + if (words.size() > 0) { + if (words[0] == "spatial") ASSERT_THAT(words[2], Eq("xyz")); + if (words[0] == "cell_spatial") ASSERT_THAT(words[2], Eq("abc")); + if (words[0] == "cell_origin") { + ++line; + words = utils::split_words(*line); + ASSERT_THAT(words[0], Eq("0,")); + ASSERT_THAT(words[1], Eq("0,")); + ASSERT_THAT(words[2], Eq("0")); + } + if (words[0] == "cell_lengths") { + ++line; + words = utils::split_words(*line); + ASSERT_THAT(words[0], Eq("3.359192,")); + ASSERT_THAT(words[1], Eq("3.359192,")); + ASSERT_THAT(words[2], Eq("3.359192")); + } + if (words[0] == "cell_angles") { + ++line; + words = utils::split_words(*line); + ASSERT_THAT(words[0], Eq("90,")); + ASSERT_THAT(words[1], Eq("90,")); + ASSERT_THAT(words[2], Eq("90")); + } + if (words[0] == "id") { + ++line; + words = utils::split_words(*line); + ASSERT_THAT(words[0], Eq("1,")); + ASSERT_THAT(words[1], Eq("2,")); + ASSERT_THAT(words[2], Eq("3,")); + ASSERT_THAT(words[3], Eq("4,")); + ASSERT_THAT(words[4], Eq("5,")); + ASSERT_THAT(words[5], Eq("6,")); + ASSERT_THAT(words[6], Eq("7,")); + ASSERT_THAT(words[7], Eq("8,")); + ASSERT_THAT(words[8], Eq("9,")); + ASSERT_THAT(words[9], Eq("10,")); + } + if (words[0] == "mass") { + ++line; + words = utils::split_words(*line); + ASSERT_THAT(words[0], Eq("1,")); + ASSERT_THAT(words[1], Eq("1,")); + ASSERT_THAT(words[2], Eq("1,")); + ASSERT_THAT(words[3], Eq("1,")); + ASSERT_THAT(words[4], Eq("1,")); + ASSERT_THAT(words[5], Eq("1,")); + ASSERT_THAT(words[6], Eq("1,")); + ASSERT_THAT(words[7], Eq("1,")); + ASSERT_THAT(words[8], Eq("1,")); + ASSERT_THAT(words[9], Eq("1,")); + } + if (words[0] == "coordinates") { + ++line; + words = utils::split_words(*line); + ASSERT_THAT(words[0], Eq("0,")); + ASSERT_THAT(words[1], Eq("0,")); + ASSERT_THAT(words[2], Eq("0,")); + ++line; + words = utils::split_words(*line); + ASSERT_THAT(words[0], Eq("0.8397981,")); + ASSERT_THAT(words[1], Eq("0.8397981,")); + ASSERT_THAT(words[2], Eq("0,")); + ++line; + words = utils::split_words(*line); + ASSERT_THAT(words[0], Eq("0.8397981,")); + ASSERT_THAT(words[1], Eq("0,")); + ASSERT_THAT(words[2], Eq("0.8397981,")); + ++line; + words = utils::split_words(*line); + ASSERT_THAT(words[0], Eq("0,")); + ASSERT_THAT(words[1], Eq("0.8397981,")); + ASSERT_THAT(words[2], Eq("0.8397981,")); + } + if (words[0] == "ix") { + ++line; + words = utils::split_words(*line); + ASSERT_THAT(words[0], Eq("0,")); + ASSERT_THAT(words[1], Eq("0,")); + ASSERT_THAT(words[2], Eq("0,")); + ASSERT_THAT(words[3], Eq("0,")); + ASSERT_THAT(words[4], Eq("0,")); + ASSERT_THAT(words[5], Eq("0,")); + ASSERT_THAT(words[6], Eq("0,")); + ASSERT_THAT(words[7], Eq("0,")); + ASSERT_THAT(words[8], Eq("0,")); + ASSERT_THAT(words[9], Eq("0,")); + } + } + } + delete_file(converted_file); + } + delete_file(dump_file); +} + +int main(int argc, char **argv) +{ + MPI_Init(&argc, &argv); + ::testing::InitGoogleMock(&argc, argv); + + // handle arguments passed via environment variable + if (const char *var = getenv("TEST_ARGS")) { + std::vector env = utils::split_words(var); + for (auto arg : env) { + if (arg == "-v") { + verbose = true; + } + } + } + + NCDUMP_BINARY = getenv("NCDUMP_BINARY"); + + if ((argc > 1) && (strcmp(argv[1], "-v") == 0)) verbose = true; + + int rv = RUN_ALL_TESTS(); + MPI_Finalize(); + return rv; +}