[MOAB-dev] commit/MOAB: danwu: Updated NC writer and its unit test (fixed some issues caused by dummy coordinate variables). For writing a HOMME file, the header part is done. Will continue to work on the data part for non-dimension variables.
commits-noreply at bitbucket.org
commits-noreply at bitbucket.org
Tue Apr 15 10:20:43 CDT 2014
1 new commit in MOAB:
https://bitbucket.org/fathomteam/moab/commits/4160ad611092/
Changeset: 4160ad611092
Branch: ncwriter
User: danwu
Date: 2014-04-15 17:20:23
Summary: Updated NC writer and its unit test (fixed some issues caused by dummy coordinate variables). For writing a HOMME file, the header part is done. Will continue to work on the data part for non-dimension variables.
Affected #: 4 files
diff --git a/src/io/NCWriteEuler.cpp b/src/io/NCWriteEuler.cpp
index 70d5bdd..4671ed2 100644
--- a/src/io/NCWriteEuler.cpp
+++ b/src/io/NCWriteEuler.cpp
@@ -24,6 +24,7 @@ ErrorCode NCWriteEuler::write_values(std::vector<std::string>& var_names, Entity
{
Interface*& mbImpl = _writeNC->mbImpl;
std::set<std::string>& usedCoordinates = _writeNC->usedCoordinates;
+ std::set<std::string>& dummyVarNames = _writeNC->dummyVarNames;
std::map<std::string, WriteNC::VarData>& varInfo = _writeNC->varInfo;
// Start with coordinates
@@ -31,6 +32,10 @@ ErrorCode NCWriteEuler::write_values(std::vector<std::string>& var_names, Entity
setIt != usedCoordinates.end(); ++setIt) {
std::string coordName = *setIt; // Deep copy
+ // Skip dummy coordinate variables (if any)
+ if (dummyVarNames.find(coordName) != dummyVarNames.end())
+ continue;
+
std::map<std::string, WriteNC::VarData>::iterator vit = varInfo.find(coordName);
if (vit == varInfo.end())
ERRORR(MB_FAILURE, "Can't find one coordinate variable.");
@@ -56,9 +61,7 @@ ErrorCode NCWriteEuler::write_values(std::vector<std::string>& var_names, Entity
}
// Now look at requested var_names; if they have time, we will have a list, and write one at a time
- // We may also need to gather, and transpose stuff
- // Get all entities of dimension 2 from set
- // Need to reorder stuff in the order from the file, also transpose from lev dimension
+ // Need to transpose from lev dimension
ErrorCode rval;
// For each variable tag in the indexed lists, write a time step data
@@ -74,28 +77,16 @@ ErrorCode NCWriteEuler::write_values(std::vector<std::string>& var_names, Entity
// Get entities of this variable
Range ents;
switch (variableData.entLoc) {
- case WriteNC::ENTLOCVERT:
- // Vertices
- rval = mbImpl->get_entities_by_dimension(fileSet, 0, ents);
- ERRORR(rval, "Can't get entities for vertices.");
- break;
case WriteNC::ENTLOCFACE:
// Faces
rval = mbImpl->get_entities_by_dimension(fileSet, 2, ents);
ERRORR(rval, "Can't get entities for faces.");
break;
- case WriteNC::ENTLOCNSEDGE:
- case WriteNC::ENTLOCEWEDGE:
- case WriteNC::ENTLOCEDGE:
- // Edges
- rval = mbImpl->get_entities_by_dimension(fileSet, 1, ents);
- ERRORR(rval, "Can't get entities for edges.");
- break;
default:
- break;
+ ERRORR(MB_FAILURE, "Unexpected entity location type for CAM-EUL non-set variable.");
}
- // FIXME: assume now the variable has 4 dimensions as (time, lev, lat, lon)
+ // A typical variable has 4 dimensions as (time, lev, lat, lon)
// At each timestep, we need to transpose tag format (lat, lon, lev) back
// to NC format (lev, lat, lon) for writing
size_t ni = variableData.writeCounts[3]; // lon
@@ -114,7 +105,6 @@ ErrorCode NCWriteEuler::write_values(std::vector<std::string>& var_names, Entity
assert(count == (int)ents.size());
// Now write from memory directly
- // FIXME: we need to gather for multiple processors
int success = 0;
switch (variableData.varDataType) {
case NC_DOUBLE: {
diff --git a/src/io/NCWriteHOMME.cpp b/src/io/NCWriteHOMME.cpp
index 3b500b5..8625f14 100644
--- a/src/io/NCWriteHOMME.cpp
+++ b/src/io/NCWriteHOMME.cpp
@@ -22,7 +22,50 @@ NCWriteHOMME::~NCWriteHOMME()
ErrorCode NCWriteHOMME::write_values(std::vector<std::string>& var_names, EntityHandle fileSet)
{
- return MB_NOT_IMPLEMENTED;
+ Interface*& mbImpl = _writeNC->mbImpl;
+ std::set<std::string>& usedCoordinates = _writeNC->usedCoordinates;
+ std::set<std::string>& dummyVarNames = _writeNC->dummyVarNames;
+ std::map<std::string, WriteNC::VarData>& varInfo = _writeNC->varInfo;
+
+ // Start with coordinates
+ for (std::set<std::string>::iterator setIt = usedCoordinates.begin();
+ setIt != usedCoordinates.end(); ++setIt) {
+ std::string coordName = *setIt; // Deep copy
+
+ // Skip dummy coordinate variables (e.g. ncol)
+ if (dummyVarNames.find(coordName) != dummyVarNames.end())
+ continue;
+
+ std::map<std::string, WriteNC::VarData>::iterator vit = varInfo.find(coordName);
+ if (vit == varInfo.end())
+ ERRORR(MB_FAILURE, "Can't find one coordinate variable.");
+
+ WriteNC::VarData& varCoordData = vit->second;
+
+ int success = 0;
+ switch (varCoordData.varDataType) {
+ case NC_DOUBLE:
+ success = NCFUNCAP(_vara_double)(_fileId, varCoordData.varId, &varCoordData.writeStarts[0],
+ &varCoordData.writeCounts[0], (double*)(varCoordData.memoryHogs[0]));
+ ERRORS(success, "Failed to write double data.");
+ break;
+ case NC_INT:
+ success = NCFUNCAP(_vara_int)(_fileId, varCoordData.varId, &varCoordData.writeStarts[0],
+ &varCoordData.writeCounts[0], (int*)(varCoordData.memoryHogs[0]));
+ ERRORS(success, "Failed to write int data.");
+ break;
+ default:
+ success = 1;
+ break;
+ }
+ }
+
+ // Now look at requested var_names; if they have time, we will have a list, and write one at a time
+ // Need to transpose from lev dimension
+
+ // Write data, to be implemented ...
+
+ return MB_SUCCESS;
}
} /* namespace moab */
diff --git a/src/io/WriteNC.cpp b/src/io/WriteNC.cpp
index d359c83..e05dadf 100644
--- a/src/io/WriteNC.cpp
+++ b/src/io/WriteNC.cpp
@@ -533,7 +533,7 @@ ErrorCode WriteNC::collect_variable_data(std::vector<std::string>& var_names, st
currentVarData.has_tsteps = true;
// Probably will have to look at tstep_vals to match them
- sizeVar *= dimLens[j];
+ sizeVar *= dimLens[currentVarData.varDims[j]];
usedCoordinates.insert(dimName); // Collect those used, we will need to write them to the file
dbgOut.tprintf(2, " for variable %s need dimension %s with length %d\n", varname.c_str(), dimName.c_str(), dimLens[currentVarData.varDims[j]]);
}
@@ -543,7 +543,7 @@ ErrorCode WriteNC::collect_variable_data(std::vector<std::string>& var_names, st
if (currentVarData.has_tsteps) {
int index = 0;
while (true) {
- Tag indexedTag;
+ Tag indexedTag = 0;
std::stringstream ssTagNameWithIndex;
ssTagNameWithIndex << varname << index;
rval = mbImpl->tag_get_handle(ssTagNameWithIndex.str().c_str(), indexedTag);
@@ -566,18 +566,18 @@ ErrorCode WriteNC::collect_variable_data(std::vector<std::string>& var_names, st
}
else {
// Get the tag with varname
- Tag coordtag = 0;
- rval = mbImpl->tag_get_handle(varname.c_str(), coordtag);
+ Tag tag = 0;
+ rval = mbImpl->tag_get_handle(varname.c_str(), tag);
ERRORR(rval, "Can't find one tag.");
- currentVarData.varTags.push_back(coordtag); // Really, only one for these
+ currentVarData.varTags.push_back(tag); // Really, only one for these
const void* data;
- int sizeCoordinate;
- rval = mbImpl->tag_get_by_ptr(coordtag, &fileSet, 1, &data, &sizeCoordinate);
- ERRORR(rval, "Can't get coordinate values.");
+ int size;
+ rval = mbImpl->tag_get_by_ptr(tag, &fileSet, 1, &data, &size);
+ ERRORR(rval, "Can't get tag values.");
// Find the type of tag, and use it
DataType type;
- rval = mbImpl->tag_get_data_type(coordtag, type);
+ rval = mbImpl->tag_get_data_type(tag, type);
ERRORR(rval, "Can't get tag type.");
currentVarData.varDataType = NC_DOUBLE;
@@ -599,18 +599,33 @@ ErrorCode WriteNC::collect_variable_data(std::vector<std::string>& var_names, st
ERRORR(MB_FAILURE, "Can't find one coordinate variable.");
VarData& varCoordData = vit->second;
- Tag coordtag = 0;
- rval = mbImpl->tag_get_handle(coordName.c_str(), coordtag);
+ Tag coordTag = 0;
+ rval = mbImpl->tag_get_handle(coordName.c_str(), coordTag);
ERRORR(rval, "Can't find one tag.");
- varCoordData.varTags.push_back(coordtag); // Really, only one for these
+ varCoordData.varTags.push_back(coordTag); // Really, only one for these
const void* data;
int sizeCoordinate;
- rval = mbImpl->tag_get_by_ptr(coordtag, &fileSet, 1, &data, &sizeCoordinate);
+ rval = mbImpl->tag_get_by_ptr(coordTag, &fileSet, 1, &data, &sizeCoordinate);
ERRORR(rval, "Can't get coordinate values.");
dbgOut.tprintf(2, " found coordinate tag with name %s and length %d\n", coordName.c_str(),
sizeCoordinate);
+ // Get dimension length (the only dimension of this coordinate variable, with the same name)
+ assert(1 == varCoordData.varDims.size());
+ int coordDimLen = dimLens[varCoordData.varDims[0]];
+
+ if (dummyVarNames.find(coordName) != dummyVarNames.end()) {
+ // For a dummy coordinate variable, the tag size is always 1
+ // The number of coordinates should be set to dimension length, instead of 1
+ assert(1 == sizeCoordinate);
+ sizeCoordinate = coordDimLen;
+ }
+ else {
+ // The number of coordinates should be exactly the same as dimension length
+ assert(sizeCoordinate == coordDimLen);
+ }
+
// This is the length
varCoordData.sz = sizeCoordinate;
varCoordData.writeStarts.resize(1);
@@ -620,7 +635,7 @@ ErrorCode WriteNC::collect_variable_data(std::vector<std::string>& var_names, st
// Find the type of tag, and use it
DataType type;
- rval = mbImpl->tag_get_data_type(coordtag, type);
+ rval = mbImpl->tag_get_data_type(coordTag, type);
ERRORR(rval, "Can't get tag type.");
varCoordData.varDataType = NC_DOUBLE;
@@ -639,7 +654,7 @@ ErrorCode WriteNC::initialize_file(std::vector<std::string>& var_names)
// First initialize all coordinates, then fill VarData for actual variables (and dimensions)
// Check that for used coordinates we have found the tags
for (std::set<std::string>::iterator setIt = usedCoordinates.begin();
- setIt != usedCoordinates.end(); setIt++) {
+ setIt != usedCoordinates.end(); ++setIt) {
std::string coordName = *setIt; // Deep copy
std::map<std::string, VarData>::iterator vit = varInfo.find(coordName);
@@ -667,6 +682,11 @@ ErrorCode WriteNC::initialize_file(std::vector<std::string>& var_names)
example: http://www.unidata.ucar.edu/software/netcdf/docs/netcdf-c/nc_005fdef_005fvar.html#nc_005fdef_005fvar
*/
+ // Skip dummy coordinate variables (e.g. ncol)
+ if (dummyVarNames.find(coordName) != dummyVarNames.end())
+ continue;
+
+ // Define a coordinate variable
if (NCFUNC(def_var)(fileId, coordName.c_str(), varCoordData.varDataType,
1, &(varCoordData.varDims[0]), &varCoordData.varId) != NC_NOERR)
ERRORR(MB_FAILURE, "Failed to create coordinate variable.");
diff --git a/test/io/write_nc.cpp b/test/io/write_nc.cpp
index 72256e5..5261532 100644
--- a/test/io/write_nc.cpp
+++ b/test/io/write_nc.cpp
@@ -7,9 +7,11 @@ using namespace moab;
#ifdef MESHDIR
static const char example_eul[] = STRINGIFY(MESHDIR) "/io/camEul26x48x96.t3.nc";
static const char example_homme[] = STRINGIFY(MESHDIR) "/io/homme26x3458.t.3.nc";
+static const char example_homme_mapping[] = STRINGIFY(MESHDIR) "/io/HommeMapping.nc";
#else
static const char example_eul[] = "/io/camEul26x48x96.t3.nc";
static const char example_homme[] = "/io/homme26x3458.t.3.nc";
+static const char example_homme_mapping[] = "/io/HommeMapping.nc";
#endif
#ifdef USE_MPI
@@ -18,12 +20,12 @@ static const char example_homme[] = "/io/homme26x3458.t.3.nc";
#endif
// CAM-EUL
-void test_eul_read_write_T_gw();
-void test_eul_check_T_gw_values();
+void test_eul_read_write_T();
+void test_eul_check_T();
// CAM-SE (HOMME)
void test_homme_read_write_T();
-void test_homme_check_T_values();
+void test_homme_check_T();
void get_eul_options(std::string& opts);
void get_homme_options(std::string& opts);
@@ -40,10 +42,10 @@ int main(int argc, char* argv[])
argv[0] = argv[argc - argc]; // To remove the warnings in serial mode about unused variables
#endif
- result += RUN_TEST(test_eul_read_write_T_gw);
- result += RUN_TEST(test_eul_check_T_gw_values);
- //result += RUN_TEST(test_homme_read_write_T);
- //result += RUN_TEST(test_homme_check_T_values);
+ result += RUN_TEST(test_eul_read_write_T);
+ result += RUN_TEST(test_eul_check_T);
+ result += RUN_TEST(test_homme_read_write_T);
+ //result += RUN_TEST(test_homme_check_T);
#ifdef USE_MPI
fail = MPI_Finalize();
@@ -54,7 +56,9 @@ int main(int argc, char* argv[])
return result;
}
-void test_eul_read_write_T_gw()
+// We also read and write set variable gw, which is required to create the mesh
+// In test_eul_check_T_values(), we need to load the output file with mesh
+void test_eul_read_write_T()
{
Core moab;
Interface& mb = moab;
@@ -71,15 +75,16 @@ void test_eul_read_write_T_gw()
rval = mb.load_file(example_eul, &set, opts.c_str());
CHECK_ERR(rval);
- // This test will write information about variable T and gw
- // To load the output file with mesh, variable gw is required
+ // Write variables T and gw
std::string writeopts;
writeopts = std::string(";;VARIABLE=T,gw;DEBUG_IO=2;");
- rval = mb.write_file("test_eul_T_gw.nc", 0, writeopts.c_str(), &set, 1);
+ rval = mb.write_file("test_eul_T.nc", 0, writeopts.c_str(), &set, 1);
CHECK_ERR(rval);
}
-void test_eul_check_T_gw_values()
+// Check non-set variable T on some quads
+// Also check set variable gw
+void test_eul_check_T()
{
Core moab;
Interface& mb = moab;
@@ -93,7 +98,7 @@ void test_eul_check_T_gw_values()
// Load non-set variable T, set variable gw, and the mesh
opts += std::string(";VARIABLE=T,gw");
- rval = mb.load_file("test_eul_T_gw.nc", &set, opts.c_str());
+ rval = mb.load_file("test_eul_T.nc", &set, opts.c_str());
CHECK_ERR(rval);
int procs = 1;
@@ -144,6 +149,8 @@ void test_eul_check_T_gw_values()
}
}
+// We also read and write set variables lat and lon, which are are required to create the mesh
+// In test_homme_check_T_values(), we need to load the output file with mesh
void test_homme_read_write_T()
{
Core moab;
@@ -156,19 +163,21 @@ void test_homme_read_write_T()
ErrorCode rval = mb.create_meshset(MESHSET_SET, set);
CHECK_ERR(rval);
- // Load non-set variable T and the mesh
- std::string opts = orig + std::string(";DEBUG_IO=0;VARIABLE=T");
+ // Load non-set variable T, set variable lat, set variable lon, and the mesh
+ std::string opts = orig + std::string(";DEBUG_IO=0;VARIABLE=T,lat,lon");
rval = mb.load_file(example_homme, &set, opts.c_str());
CHECK_ERR(rval);
- // This test will write information about variable T
+ // Write variables T, lat and lon
std::string writeopts;
- writeopts = std::string(";;VARIABLE=T;DEBUG_IO=0;");
+ writeopts = std::string(";;VARIABLE=T,lat,lon;DEBUG_IO=0;");
rval = mb.write_file("test_homme_T.nc", 0, writeopts.c_str(), &set, 1);
CHECK_ERR(rval);
}
-void test_homme_check_T_values()
+// Check non-set variable T on some vertices
+// Also check set variables lat and lon
+void test_homme_check_T()
{
Core moab;
Interface& mb = moab;
@@ -180,8 +189,10 @@ void test_homme_check_T_values()
ErrorCode rval = mb.create_meshset(MESHSET_SET, set);
CHECK_ERR(rval);
- // Load non-set variable T and the mesh
- opts += std::string(";VARIABLE=T");
+ // Load non-set variable T, set variable lat, set variable lon, and the mesh
+ opts += std::string(";VARIABLE=T,lat,lon");
+ opts += ";CONN=";
+ opts += example_homme_mapping;
rval = mb.load_file("test_homme_T.nc", &set, opts.c_str());
CHECK_ERR(rval);
@@ -193,6 +204,40 @@ void test_homme_check_T_values()
// Only test serial case for the time being
if (1 == procs) {
+ // Get tag lat
+ Tag lat_tag;
+ rval = mb.tag_get_handle("lat", 0, MB_TYPE_OPAQUE, lat_tag, MB_TAG_SPARSE | MB_TAG_VARLEN);
+ CHECK_ERR(rval);
+
+ // Check some values of tag lat
+ const void* var_data;
+ int var_len = 0;
+ rval = mb.tag_get_by_ptr(lat_tag, &set, 1, &var_data, &var_len);
+ CHECK_ERR(rval);
+ CHECK_EQUAL(3458, var_len);
+ double* lat_val = (double*)var_data;
+ double eps = 1e-10;
+ CHECK_REAL_EQUAL(-35.2643896827547, lat_val[0], eps);
+ CHECK_REAL_EQUAL(23.8854752772335, lat_val[1728], eps);
+ CHECK_REAL_EQUAL(29.8493120043874, lat_val[1729], eps);
+ CHECK_REAL_EQUAL(38.250274171077, lat_val[3457], eps);
+
+ // Get tag lon
+ Tag lon_tag;
+ rval = mb.tag_get_handle("lon", 0, MB_TYPE_OPAQUE, lon_tag, MB_TAG_SPARSE | MB_TAG_VARLEN);
+ CHECK_ERR(rval);
+
+ // Check some values of tag lon
+ var_len = 0;
+ rval = mb.tag_get_by_ptr(lon_tag, &set, 1, &var_data, &var_len);
+ CHECK_ERR(rval);
+ CHECK_EQUAL(3458, var_len);
+ double* lon_val = (double*)var_data;
+ CHECK_REAL_EQUAL(315, lon_val[0], eps);
+ CHECK_REAL_EQUAL(202.5, lon_val[1728], eps);
+ CHECK_REAL_EQUAL(194.359423525313, lon_val[1729], eps);
+ CHECK_REAL_EQUAL(135, lon_val[3457], eps);
+
// Get tag T0
Tag Ttag0;
rval = mb.tag_get_handle("T0", 26, MB_TYPE_DOUBLE, Ttag0);
@@ -212,7 +257,7 @@ void test_homme_check_T_values()
CHECK_EQUAL((size_t)count, verts.size());
// Check some values of tag T0 on first level
- const double eps = 0.0001;
+ eps = 0.0001;
double* data = (double*) Tbuf;
CHECK_REAL_EQUAL(233.1136, data[0 * 26], eps); // First vert
CHECK_REAL_EQUAL(236.1505, data[1728 * 26], eps); // Median vert
Repository URL: https://bitbucket.org/fathomteam/moab/
--
This is a commit notification from bitbucket.org. You are receiving
this because you have the service enabled, addressing the recipient of
this email.
More information about the moab-dev
mailing list