[MOAB-dev] commit/MOAB: danwu: Merged ncwriter into master
commits-noreply at bitbucket.org
commits-noreply at bitbucket.org
Mon Apr 28 17:02:18 CDT 2014
1 new commit in MOAB:
https://bitbucket.org/fathomteam/moab/commits/7e36def04100/
Changeset: 7e36def04100
Branch: master
User: danwu
Date: 2014-04-29 00:02:09
Summary: Merged ncwriter into master
Affected #: 20 files
diff --git a/.gitignore b/.gitignore
index ba9fe1e..d686ffa 100644
--- a/.gitignore
+++ b/.gitignore
@@ -155,6 +155,7 @@ test/io/*.g
test/io/gmsh_test
test/io/ideas_test
test/io/nastran_test
+test/io/*.nc
test/io/read_cgm_basic_test
test/io/read_cgm_connectivity_test
test/io/read_cgm_group_test
@@ -168,6 +169,7 @@ test/io/smf_test
test/io/stl_test
test/io/tqdcfr
test/io/vtk_test
+test/io/write_nc
test/kd_tree_test
test/kd_tree_time
test/kd_tree_tool
@@ -187,6 +189,7 @@ test/parallel/*.h5m
test/parallel/mbparallelcomm_test
test/parallel/mhdf_parallel
test/parallel/mpastrvpart
+test/parallel/*.nc
test/parallel/par_coupler_test
test/parallel/par_intx_sph
test/parallel/parallel_hdf5_test
@@ -206,6 +209,7 @@ test/parallel/structured3
test/parallel/uber_parallel_test
test/parallel/ucdtrvpart
test/parallel/*.vtk
+test/parallel/write_nc_par
test/perf/adj_time
test/perf/perf
test/perf/perftool
diff --git a/src/ReaderWriterSet.cpp b/src/ReaderWriterSet.cpp
index 5015a7f..1f254b0 100644
--- a/src/ReaderWriterSet.cpp
+++ b/src/ReaderWriterSet.cpp
@@ -48,6 +48,7 @@
#ifdef NETCDF_FILE
# include "ReadNCDF.hpp"
# include "WriteNCDF.hpp"
+# include "WriteNC.hpp"
# include "WriteSLAC.hpp"
# include "ReadNC.hpp"
# include "ReadGCRM.hpp"
@@ -105,7 +106,7 @@ ReaderWriterSet::ReaderWriterSet( Core* mdb, Error* handler )
const char* exo_sufxs[] = { "exo", "exoII", "exo2", "g", "gen", NULL };
register_factory( ReadNCDF::factory, WriteNCDF::factory, "Exodus II", exo_sufxs, "EXODUS" );
register_factory( ReadGCRM::factory, NULL, "GCRM NC", "nc", "GCRM" );
- register_factory( ReadNC::factory, NULL, "Climate NC", "nc", "NC" );
+ register_factory( ReadNC::factory, WriteNC::factory, "Climate NC", "nc", "NC" );
#endif
#ifdef CGNS_FILE
diff --git a/src/io/Makefile.am b/src/io/Makefile.am
index 09fe3d3..1269523 100644
--- a/src/io/Makefile.am
+++ b/src/io/Makefile.am
@@ -17,11 +17,17 @@ if NETCDF_FILE
WriteNCDF.cpp WriteNCDF.hpp \
WriteSLAC.cpp WriteSLAC.hpp \
ReadNC.cpp ReadNC.hpp \
+ WriteNC.cpp WriteNC.hpp \
NCHelper.cpp NCHelper.hpp \
NCHelperEuler.cpp NCHelperEuler.hpp \
NCHelperFV.cpp NCHelperFV.hpp \
NCHelperHOMME.cpp NCHelperHOMME.hpp \
NCHelperMPAS.cpp NCHelperMPAS.hpp \
+ NCWriteHelper.cpp NCWriteHelper.hpp \
+ NCWriteEuler.cpp NCWriteEuler.hpp \
+ NCWriteFV.cpp NCWriteFV.hpp \
+ NCWriteHOMME.cpp NCWriteHOMME.hpp \
+ NCWriteMPAS.cpp NCWriteMPAS.hpp \
ReadGCRM.cpp ReadGCRM.hpp
else
MOAB_NETCDF_SRCS =
@@ -30,11 +36,17 @@ endif
if PNETCDF_FILE
if !NETCDF_FILE
MOAB_NETCDF_SRCS += ReadNC.cpp ReadNC.hpp \
+ WriteNC.cpp WriteNC.hpp \
NCHelper.cpp NCHelper.hpp \
NCHelperEuler.cpp NCHelperEuler.hpp \
NCHelperFV.cpp NCHelperFV.hpp \
NCHelperHOMME.cpp NCHelperHOMME.hpp \
NCHelperMPAS.cpp NCHelperMPAS.hpp \
+ NCWriteHelper.cpp NCWriteHelper.hpp \
+ NCWriteEuler.cpp NCWriteEuler.hpp \
+ NCWriteFV.cpp NCWriteFV.hpp \
+ NCWriteHOMME.cpp NCWriteHOMME.hpp \
+ NCWriteMPAS.cpp NCWriteMPAS.hpp \
ReadGCRM.cpp ReadGCRM.hpp
endif
endif
diff --git a/src/io/NCHelper.cpp b/src/io/NCHelper.cpp
index 9ef1332..0e970d9 100644
--- a/src/io/NCHelper.cpp
+++ b/src/io/NCHelper.cpp
@@ -192,16 +192,16 @@ ErrorCode NCHelper::create_conventional_tags(const std::vector<int>& tstep_nums)
unsigned int varDimSz = varInfo[mapIter->first].varDims.size();
if (varDimSz == 0)
continue;
- varInfo[mapIter->first].varTags.resize(varDimSz, 0);
+ std::vector<Tag> varDimTags(varDimSz);
for (unsigned int i = 0; i != varDimSz; i++) {
Tag tmptag = 0;
std::string tmptagname = dimNames[varInfo[mapIter->first].varDims[i]];
mbImpl->tag_get_handle(tmptagname.c_str(), 0, MB_TYPE_OPAQUE, tmptag, MB_TAG_ANY);
- varInfo[mapIter->first].varTags[i] = tmptag;
+ varDimTags[i] = tmptag;
}
rval = mbImpl->tag_get_handle(tag_name.c_str(), varDimSz, MB_TYPE_HANDLE, varNamesDimsTag, MB_TAG_SPARSE | MB_TAG_CREAT);
ERRORR(rval, "Trouble creating __<var_name>_DIMS tag.");
- rval = mbImpl->tag_set_data(varNamesDimsTag, &_fileSet, 1, &(varInfo[mapIter->first].varTags[0]));
+ rval = mbImpl->tag_set_data(varNamesDimsTag, &_fileSet, 1, &(varDimTags[0]));
ERRORR(rval, "Trouble setting data for __<var_name>_DIMS tag.");
if (MB_SUCCESS == rval)
dbgOut.tprintf(2, "Tag created for variable %s\n", tag_name.c_str());
@@ -252,21 +252,37 @@ ErrorCode NCHelper::create_conventional_tags(const std::vector<int>& tstep_nums)
Tag varAttTag = 0;
rval = mbImpl->tag_get_handle(tag_name.c_str(), 0, MB_TYPE_OPAQUE, varAttTag, MB_TAG_CREAT | MB_TAG_SPARSE | MB_TAG_VARLEN);
ERRORR(rval, "Trouble creating __<var_name>_ATTRIBS tag.");
+
std::string varAttVal;
std::vector<int> varAttLen;
- rval = create_attrib_string(mapIter->second.varAtts, varAttVal, varAttLen);
- ERRORR(rval, "Trouble creating attribute strings.");
+ if (mapIter->second.numAtts < 1) {
+ if (dummyVarNames.find(mapIter->first) != dummyVarNames.end()) {
+ // This variable is a dummy dimension variable
+ varAttVal = "DUMMY_VAR";
+ }
+ else {
+ // This variable has no attributes
+ varAttVal = "NO_ATTRIBS";
+ }
+ }
+ else {
+ rval = create_attrib_string(mapIter->second.varAtts, varAttVal, varAttLen);
+ ERRORR(rval, "Trouble creating attribute string.");
+ }
const void* varAttPtr = varAttVal.c_str();
int varAttSz = varAttVal.size();
+ if (0 == varAttSz)
+ varAttSz = 1;
rval = mbImpl->tag_set_by_ptr(varAttTag, &_fileSet, 1, &varAttPtr, &varAttSz);
ERRORR(rval, "Trouble setting data for __<var_name>_ATTRIBS tag.");
if (MB_SUCCESS == rval)
dbgOut.tprintf(2, "Tag created for variable %s\n", tag_name.c_str());
- if (varAttLen.size() == 0)
- varAttLen.push_back(0);
+
ssTagName << "_LEN";
tag_name = ssTagName.str();
Tag varAttLenTag = 0;
+ if (0 == varAttLen.size())
+ varAttLen.push_back(0);
rval = mbImpl->tag_get_handle(tag_name.c_str(), varAttLen.size(), MB_TYPE_INTEGER, varAttLenTag, MB_TAG_SPARSE | MB_TAG_CREAT);
ERRORR(rval, "Trouble creating __<var_name>_ATTRIBS_LEN tag.");
rval = mbImpl->tag_set_data(varAttLenTag, &_fileSet, 1, &varAttLen[0]);
@@ -324,6 +340,11 @@ ErrorCode NCHelper::read_variable_setup(std::vector<std::string>& var_names, std
dummyVarNames.find(vd.varName) != dummyVarNames.end())
continue;
+ // Dimension variables were read before creating conventional tags
+ std::vector<std::string>& dimNames = _readNC->dimNames;
+ if (std::find(dimNames.begin(), dimNames.end(), vd.varName) != dimNames.end())
+ continue;
+
if (vd.entLoc == ReadNC::ENTLOCSET)
vsetdatas.push_back(vd);
else
diff --git a/src/io/NCWriteEuler.cpp b/src/io/NCWriteEuler.cpp
new file mode 100644
index 0000000..4b3f1ff
--- /dev/null
+++ b/src/io/NCWriteEuler.cpp
@@ -0,0 +1,23 @@
+/*
+ * NCWriteEuler.cpp
+ *
+ * Created on: Mar 28, 2014
+ */
+
+#include "NCWriteEuler.hpp"
+#include "moab/WriteUtilIface.hpp"
+
+#define ERRORR(rval, str) \
+ if (MB_SUCCESS != rval) { _writeNC->mWriteIface->report_error("%s", str); return rval; }
+
+#define ERRORS(err, str) \
+ if (err) { _writeNC->mWriteIface->report_error("%s", str); return MB_FAILURE; }
+
+namespace moab {
+
+NCWriteEuler::~NCWriteEuler()
+{
+ // TODO Auto-generated destructor stub
+}
+
+} /* namespace moab */
diff --git a/src/io/NCWriteEuler.hpp b/src/io/NCWriteEuler.hpp
new file mode 100644
index 0000000..2cb6486
--- /dev/null
+++ b/src/io/NCWriteEuler.hpp
@@ -0,0 +1,27 @@
+/*
+ * NCWriteEuler.hpp
+ *
+ * nc write helper for euler type data (CAM)
+ * Created on: Mar 28, 2014
+ *
+ */
+
+#ifndef NCWRITEEULER_HPP_
+#define NCWRITEEULER_HPP_
+
+#include "NCWriteHelper.hpp"
+
+namespace moab {
+
+class NCWriteEuler: public ScdNCWriteHelper
+{
+public:
+ NCWriteEuler(WriteNC* writeNC, int fileId, const FileOptions& opts, EntityHandle fileSet)
+: ScdNCWriteHelper(writeNC, fileId, opts, fileSet) {}
+
+ virtual ~NCWriteEuler();
+};
+
+} // namespace moab
+
+#endif
diff --git a/src/io/NCWriteFV.cpp b/src/io/NCWriteFV.cpp
new file mode 100644
index 0000000..8de489f
--- /dev/null
+++ b/src/io/NCWriteFV.cpp
@@ -0,0 +1,23 @@
+/*
+ * NCWriteFV.cpp
+ *
+ * Created on: April 9, 2014
+ */
+
+#include "NCWriteFV.hpp"
+#include "moab/WriteUtilIface.hpp"
+
+#define ERRORR(rval, str) \
+ if (MB_SUCCESS != rval) { _writeNC->mWriteIface->report_error("%s", str); return rval; }
+
+#define ERRORS(err, str) \
+ if (err) { _writeNC->mWriteIface->report_error("%s", str); return MB_FAILURE; }
+
+namespace moab {
+
+NCWriteFV::~NCWriteFV()
+{
+ // TODO Auto-generated destructor stub
+}
+
+} /* namespace moab */
diff --git a/src/io/NCWriteFV.hpp b/src/io/NCWriteFV.hpp
new file mode 100644
index 0000000..d3286fc
--- /dev/null
+++ b/src/io/NCWriteFV.hpp
@@ -0,0 +1,27 @@
+/*
+ * NCWriteFV.hpp
+ *
+ * nc write helper for FV type data (CAM)
+ * Created on: April 9, 2014
+ *
+ */
+
+#ifndef NCWRITEFV_HPP_
+#define NCWRITEFV_HPP_
+
+#include "NCWriteHelper.hpp"
+
+namespace moab {
+
+class NCWriteFV: public ScdNCWriteHelper
+{
+public:
+ NCWriteFV(WriteNC* writeNC, int fileId, const FileOptions& opts, EntityHandle fileSet)
+: ScdNCWriteHelper(writeNC, fileId, opts, fileSet) {}
+
+ virtual ~NCWriteFV();
+};
+
+} // namespace moab
+
+#endif // NCWRITEFV_HPP_
diff --git a/src/io/NCWriteHOMME.cpp b/src/io/NCWriteHOMME.cpp
new file mode 100644
index 0000000..9320f62
--- /dev/null
+++ b/src/io/NCWriteHOMME.cpp
@@ -0,0 +1,316 @@
+/*
+ * NCWriteHOMME.cpp
+ *
+ * Created on: April 9, 2014
+ */
+
+#include "NCWriteHOMME.hpp"
+#include "moab/WriteUtilIface.hpp"
+#include "MBTagConventions.hpp"
+
+#define ERRORR(rval, str) \
+ if (MB_SUCCESS != rval) { _writeNC->mWriteIface->report_error("%s", str); return rval; }
+
+#define ERRORS(err, str) \
+ if (err) { _writeNC->mWriteIface->report_error("%s", str); return MB_FAILURE; }
+
+namespace moab {
+
+NCWriteHOMME::~NCWriteHOMME()
+{
+ // TODO Auto-generated destructor stub
+}
+
+ErrorCode NCWriteHOMME::collect_mesh_info()
+{
+ Interface*& mbImpl = _writeNC->mbImpl;
+ std::vector<std::string>& dimNames = _writeNC->dimNames;
+ std::vector<int>& dimLens = _writeNC->dimLens;
+ Tag& mGlobalIdTag = _writeNC->mGlobalIdTag;
+
+ ErrorCode rval;
+
+ // Look for time dimension
+ std::vector<std::string>::iterator vecIt;
+ if ((vecIt = std::find(dimNames.begin(), dimNames.end(), "time")) != dimNames.end())
+ tDim = vecIt - dimNames.begin();
+ else {
+ ERRORR(MB_FAILURE, "Couldn't find 'time' dimension.");
+ }
+ nTimeSteps = dimLens[tDim];
+
+ // Get number of levels
+ if ((vecIt = std::find(dimNames.begin(), dimNames.end(), "lev")) != dimNames.end())
+ levDim = vecIt - dimNames.begin();
+ else {
+ ERRORR(MB_FAILURE, "Couldn't find 'lev' dimension.");
+ }
+ nLevels = dimLens[levDim];
+
+ Range local_verts;
+ rval = mbImpl->get_entities_by_dimension(_fileSet, 0, local_verts);
+ ERRORR(rval, "Trouble getting local vertices in current file set.");
+ assert(!local_verts.empty());
+
+#ifdef USE_MPI
+ bool& isParallel = _writeNC->isParallel;
+ if (isParallel) {
+ ParallelComm*& myPcomm = _writeNC->myPcomm;
+ int rank = myPcomm->proc_config().proc_rank();
+ int procs = myPcomm->proc_config().proc_size();
+ if (procs > 1) {
+ rval = myPcomm->filter_pstatus(local_verts, PSTATUS_NOT_OWNED, PSTATUS_NOT, -1, &localVertsOwned);
+ ERRORR(rval, "Trouble getting owned vertices in set.");
+ // Assume that PARALLEL_RESOLVE_SHARED_ENTS option is set
+ // We should avoid writing in parallel with overlapped data
+ if (rank > 0)
+ assert("PARALLEL_RESOLVE_SHARED_ENTS option is set" && localVertsOwned.size() < local_verts.size());
+ }
+ else
+ localVertsOwned = local_verts;
+ }
+ else
+ localVertsOwned = local_verts; // Not running in parallel, but still with MPI
+#else
+ localVertsOwned = local_verts;
+#endif
+
+ std::vector<int> gids(localVertsOwned.size());
+ rval = mbImpl->tag_get_data(mGlobalIdTag, localVertsOwned, &gids[0]);
+ ERRORR(rval, "Trouble getting global IDs on local vertices.");
+
+ // Restore localGidVerts
+ std::copy(gids.rbegin(), gids.rend(), range_inserter(localGidVertsOwned));
+ nLocalVerticesOwned = localGidVertsOwned.size();
+
+ return MB_SUCCESS;
+}
+
+ErrorCode NCWriteHOMME::collect_variable_data(std::vector<std::string>& var_names)
+{
+ NCWriteHelper::collect_variable_data(var_names);
+
+ std::map<std::string, WriteNC::VarData>& varInfo = _writeNC->varInfo;
+
+ for (size_t i = 0; i < var_names.size(); i++) {
+ std::string varname = var_names[i];
+ std::map<std::string, WriteNC::VarData>::iterator vit = varInfo.find(varname);
+ if (vit == varInfo.end())
+ ERRORR(MB_FAILURE, "Can't find one variable.");
+
+ WriteNC::VarData& currentVarData = vit->second;
+ if (currentVarData.has_tsteps) {
+ // Support non-set variables with 3 dimensions like (time, lev, ncol)
+ assert(3 == currentVarData.varDims.size());
+
+ // Time should be the first dimension
+ assert(tDim == currentVarData.varDims[0]);
+
+ // Set up writeStarts and writeCounts
+ currentVarData.writeStarts.resize(3);
+ currentVarData.writeCounts.resize(3);
+
+ // First: time
+ currentVarData.writeStarts[0] = 0; // This value is timestep dependent, will be set later
+ currentVarData.writeCounts[0] = 1;
+
+ // Next: lev
+ currentVarData.writeStarts[1] = 0;
+ currentVarData.writeCounts[1] = currentVarData.numLev;
+
+ // Finally: ncol
+ switch (currentVarData.entLoc) {
+ case WriteNC::ENTLOCVERT:
+ // Vertices
+ // Start from the first localGidVerts
+ // Actually, this will be reset later for writing
+ currentVarData.writeStarts[2] = localGidVertsOwned[0] - 1;
+ currentVarData.writeCounts[2] = nLocalVerticesOwned;
+ break;
+ default:
+ ERRORR(MB_FAILURE, "Unexpected entity location type for HOMME non-set variable.");
+ }
+ }
+
+ // Get variable size
+ currentVarData.sz = 1;
+ for (std::size_t idx = 0; idx != currentVarData.writeCounts.size(); idx++)
+ currentVarData.sz *= currentVarData.writeCounts[idx];
+ }
+
+ return MB_SUCCESS;
+}
+
+ErrorCode NCWriteHOMME::write_values(std::vector<std::string>& var_names)
+{
+ 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;
+
+ int success;
+
+ // 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
+ // For each variable tag in the indexed lists, write a time step data
+ // Assume the first dimension is time (need to check); if not, just write regularly
+ for (size_t i = 0; i < var_names.size(); i++) {
+ std::map<std::string, WriteNC::VarData>::iterator vit = varInfo.find(var_names[i]);
+ if (vit == varInfo.end())
+ ERRORR(MB_FAILURE, "Can't find variable requested.");
+
+ WriteNC::VarData& variableData = vit->second;
+ int numTimeSteps = (int)variableData.varTags.size();
+ if (variableData.has_tsteps) {
+ // Time should be the first dimension
+ assert(tDim == variableData.varDims[0]);
+
+ // Assume this variable is on vertices for the time being
+ switch (variableData.entLoc) {
+ case WriteNC::ENTLOCVERT:
+ // Vertices
+ break;
+ default:
+ ERRORR(MB_FAILURE, "Unexpected entity location type for HOMME non-set variable.");
+ }
+
+ // A typical variable has 3 dimensions as (time, lev, ncol)
+ // At each timestep, we need to transpose tag format (ncol, lev) back
+ // to NC format (lev, ncol) for writing
+ // FIXME: Should use tstep_nums (from writing options) later
+ for (int j = 0; j < numTimeSteps; j++) {
+ // We will write one time step, and count will be one; start will be different
+ // Use tag_get_data instead of tag_iterate to get values, as localVertsOwned
+ // might not be contiguous. We should also transpose for level so that means
+ // deep copy for transpose
+ variableData.writeStarts[0] = j; // This is time, again
+ std::vector<double> tag_data(nLocalVerticesOwned * variableData.numLev);
+ ErrorCode rval = mbImpl->tag_get_data(variableData.varTags[j], localVertsOwned, &tag_data[0]);
+ ERRORR(rval, "Trouble getting tag data on owned vertices.");
+
+#ifdef PNETCDF_FILE
+ size_t nb_writes = localGidVertsOwned.psize();
+ std::vector<int> requests(nb_writes), statuss(nb_writes);
+ size_t idxReq = 0;
+#endif
+
+ // Now write from memory directly
+ switch (variableData.varDataType) {
+ case NC_DOUBLE: {
+ std::vector<double> tmpdoubledata(nLocalVerticesOwned * variableData.numLev);
+ // Transpose (ncol, lev) back to (lev, ncol)
+ jik_to_kji(nLocalVerticesOwned, 1, variableData.numLev, &tmpdoubledata[0], &tag_data[0]);
+
+ size_t indexInDoubleArray = 0;
+ size_t ic = 0;
+ for (Range::pair_iterator pair_iter = localGidVertsOwned.pair_begin();
+ pair_iter != localGidVertsOwned.pair_end(); ++pair_iter, ic++) {
+ EntityHandle starth = pair_iter->first;
+ EntityHandle endh = pair_iter->second;
+ variableData.writeStarts[2] = (NCDF_SIZE)(starth - 1);
+ variableData.writeCounts[2] = (NCDF_SIZE)(endh - starth + 1);
+
+ // Do a partial write, in each subrange
+#ifdef PNETCDF_FILE
+ // Wait outside this loop
+ success = NCFUNCREQP(_vara_double)(_fileId, variableData.varId,
+ &(variableData.writeStarts[0]), &(variableData.writeCounts[0]),
+ &(tmpdoubledata[indexInDoubleArray]), &requests[idxReq++]);
+#else
+ success = NCFUNCAP(_vara_double)(_fileId, variableData.varId,
+ &(variableData.writeStarts[0]), &(variableData.writeCounts[0]),
+ &(tmpdoubledata[indexInDoubleArray]));
+#endif
+ ERRORS(success, "Failed to read double data in loop");
+ // We need to increment the index in double array for the
+ // next subrange
+ indexInDoubleArray += (endh - starth + 1) * variableData.numLev;
+ }
+ assert(ic == localGidVertsOwned.psize());
+#ifdef PNETCDF_FILE
+ success = ncmpi_wait_all(_fileId, requests.size(), &requests[0], &statuss[0]);
+ ERRORS(success, "Failed on wait_all.");
+#endif
+ break;
+ }
+ default:
+ ERRORR(MB_FAILURE, "Not implemented yet.");
+ }
+ }
+ } // if (variableData.has_tsteps)
+ else {
+ switch (variableData.varDataType) {
+ case NC_DOUBLE:
+ success = NCFUNCAP(_vara_double)(_fileId, variableData.varId, &variableData.writeStarts[0],
+ &variableData.writeCounts[0], (double*)(variableData.memoryHogs[0]));
+ ERRORS(success, "Failed to write double data.");
+ break;
+ default:
+ ERRORR(MB_FAILURE, "Not implemented yet.");
+ }
+ }
+ }
+
+ // Write coordinates used by requested var_names
+ // Use independent I/O mode put, since this write is only for the root processor
+ // CAUTION: if the NetCDF ID is from a previous call to ncmpi_create rather than ncmpi_open,
+ // all processors need to call ncmpi_begin_indep_data(). If only the root processor does so,
+ // ncmpi_begin_indep_data() call will be blocked forever :(
+#ifdef PNETCDF_FILE
+ // Enter independent I/O mode
+ success = NCFUNC(begin_indep_data)(_fileId);
+ ERRORS(success, "Failed to begin independent I/O mode.");
+#endif
+
+ int rank = 0;
+#ifdef USE_MPI
+ bool& isParallel = _writeNC->isParallel;
+ if (isParallel) {
+ ParallelComm*& myPcomm = _writeNC->myPcomm;
+ rank = myPcomm->proc_config().proc_rank();
+ }
+#endif
+ if (0 == rank) {
+ for (std::set<std::string>::iterator setIt = usedCoordinates.begin();
+ setIt != usedCoordinates.end(); ++setIt) {
+ const std::string& coordName = *setIt;
+
+ // 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;
+
+ switch (varCoordData.varDataType) {
+ case NC_DOUBLE:
+ // Independent I/O mode put
+ success = NCFUNCP(_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:
+ // Independent I/O mode put
+ success = NCFUNCP(_vara_int)(_fileId, varCoordData.varId, &varCoordData.writeStarts[0],
+ &varCoordData.writeCounts[0], (int*)(varCoordData.memoryHogs[0]));
+ ERRORS(success, "Failed to write int data.");
+ break;
+ default:
+ ERRORR(MB_FAILURE, "Not implemented yet.");
+ }
+ }
+ }
+
+#ifdef PNETCDF_FILE
+ // End independent I/O mode
+ success = NCFUNC(end_indep_data)(_fileId);
+ ERRORS(success, "Failed to end independent I/O mode.");
+#endif
+
+ return MB_SUCCESS;
+}
+
+} /* namespace moab */
diff --git a/src/io/NCWriteHOMME.hpp b/src/io/NCWriteHOMME.hpp
new file mode 100644
index 0000000..7bb99d9
--- /dev/null
+++ b/src/io/NCWriteHOMME.hpp
@@ -0,0 +1,37 @@
+/*
+ * NCWriteHOMME.hpp
+ *
+ * nc write helper for HOMME type data (CAM)
+ * Created on: April 9, 2014
+ *
+ */
+
+#ifndef NCWRITEHOMME_HPP_
+#define NCWRITEHOMME_HPP_
+
+#include "NCWriteHelper.hpp"
+
+namespace moab {
+
+class NCWriteHOMME: public UcdNCWriteHelper
+{
+public:
+ NCWriteHOMME(WriteNC* writeNC, int fileId, const FileOptions& opts, EntityHandle fileSet)
+: UcdNCWriteHelper(writeNC, fileId, opts, fileSet) {}
+
+ virtual ~NCWriteHOMME();
+
+private:
+ //! Implementation of NCWriteHelper::collect_mesh_info()
+ virtual ErrorCode collect_mesh_info();
+
+ //! Collect data for specified variables
+ virtual ErrorCode collect_variable_data(std::vector<std::string>& var_names);
+
+ //! Implementation of NCWriteHelper::write_values()
+ virtual ErrorCode write_values(std::vector<std::string>& var_names);
+};
+
+} // namespace moab
+
+#endif
diff --git a/src/io/NCWriteHelper.cpp b/src/io/NCWriteHelper.cpp
new file mode 100644
index 0000000..4f95ef8
--- /dev/null
+++ b/src/io/NCWriteHelper.cpp
@@ -0,0 +1,622 @@
+/*
+ * NCWriteHelper.cpp
+ *
+ * Created on: Mar 28, 2014
+ * Author: iulian
+ */
+
+#include "NCWriteHelper.hpp"
+#include "NCWriteEuler.hpp"
+#include "NCWriteFV.hpp"
+#include "NCWriteHOMME.hpp"
+#include "NCWriteMPAS.hpp"
+
+#include "moab/WriteUtilIface.hpp"
+
+#include <sstream>
+
+#define ERRORR(rval, str) \
+ if (MB_SUCCESS != rval) { _writeNC->mWriteIface->report_error("%s", str); return rval; }
+
+#define ERRORS(err, str) \
+ if (err) { _writeNC->mWriteIface->report_error("%s", str); return MB_FAILURE; }
+
+namespace moab {
+
+//! Get appropriate helper instance for WriteNC class; based on some info in the file set
+NCWriteHelper* NCWriteHelper::get_nc_helper(WriteNC* writeNC, int fileId, const FileOptions& opts, EntityHandle fileSet)
+{
+ std::string& grid_type = writeNC->grid_type;
+ if (grid_type == "CAM_EUL")
+ return new (std::nothrow) NCWriteEuler(writeNC, fileId, opts, fileSet);
+ else if (grid_type == "CAM_FV")
+ return new (std::nothrow) NCWriteFV(writeNC, fileId, opts, fileSet);
+ else if (grid_type == "CAM_SE")
+ return new (std::nothrow) NCWriteHOMME(writeNC, fileId, opts, fileSet);
+ else if (grid_type == "MPAS")
+ return new (std::nothrow) NCWriteMPAS(writeNC, fileId, opts, fileSet);
+
+ // Unknown NetCDF grid
+ return NULL;
+}
+
+ErrorCode NCWriteHelper::collect_variable_data(std::vector<std::string>& var_names)
+{
+ Interface*& mbImpl = _writeNC->mbImpl;
+ std::vector<std::string>& dimNames = _writeNC->dimNames;
+ std::vector<int>& dimLens = _writeNC->dimLens;
+ std::set<std::string>& usedCoordinates = _writeNC->usedCoordinates;
+ std::set<std::string>& dummyVarNames = _writeNC->dummyVarNames;
+ std::map<std::string, WriteNC::VarData>& varInfo = _writeNC->varInfo;
+ DebugOutput& dbgOut = _writeNC->dbgOut;
+
+ ErrorCode rval;
+
+ usedCoordinates.clear();
+
+ for (size_t i = 0; i < var_names.size(); i++) {
+ std::string varname = var_names[i];
+ std::map<std::string, WriteNC::VarData>::iterator vit = varInfo.find(varname);
+ if (vit == varInfo.end())
+ ERRORR(MB_FAILURE, "Can't find one variable.");
+
+ WriteNC::VarData& currentVarData = vit->second;
+
+ currentVarData.has_tsteps = false;
+ if ((std::find(currentVarData.varDims.begin(), currentVarData.varDims.end(), tDim) != currentVarData.varDims.end())
+ && (currentVarData.varDims.size() > 1)) // So it is not time itself
+ currentVarData.has_tsteps = true;
+
+ currentVarData.numLev = 1;
+ if ((std::find(currentVarData.varDims.begin(), currentVarData.varDims.end(), levDim) != currentVarData.varDims.end()))
+ currentVarData.numLev = nLevels;
+
+ dbgOut.tprintf(2, " for variable %s varDims.size %d \n", varname.c_str(), (int)currentVarData.varDims.size());
+ for (size_t j = 0; j < currentVarData.varDims.size(); j++) {
+ std::string dimName = dimNames[currentVarData.varDims[j]];
+ vit = varInfo.find(dimName);
+ if (vit == varInfo.end())
+ ERRORR(MB_FAILURE, "Can't find one coordinate variable.");
+
+ 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]]);
+ }
+
+ // Process coordinate variables later
+ if (usedCoordinates.find(varname) != usedCoordinates.end())
+ continue;
+
+ // Process non-set variables with time steps
+ if (currentVarData.has_tsteps) {
+ int index = 0;
+ // FIXME: Should use tstep_nums (from writing options) later
+ while (true) {
+ Tag indexedTag = 0;
+ std::stringstream ssTagNameWithIndex;
+ ssTagNameWithIndex << varname << index;
+ rval = mbImpl->tag_get_handle(ssTagNameWithIndex.str().c_str(), indexedTag);
+ if (MB_SUCCESS != rval)
+ break;
+ dbgOut.tprintf(2, " found indexed tag %d with name %s\n", index, ssTagNameWithIndex.str().c_str());
+ currentVarData.varTags.push_back(indexedTag);
+ index++; // We should get out of the loop at some point
+
+ // The type of the tag is fixed though
+ DataType type;
+ rval = mbImpl->tag_get_data_type(indexedTag, type);
+ ERRORR(rval, "Can't get tag type.");
+
+ currentVarData.varDataType = NC_DOUBLE;
+ if (MB_TYPE_INTEGER == type)
+ currentVarData.varDataType = NC_INT;
+ }
+ }
+ // Process set variables that are not coordinate variables
+ else {
+ // Get the tag with varname
+ Tag tag = 0;
+ rval = mbImpl->tag_get_handle(varname.c_str(), tag);
+ ERRORR(rval, "Can't find one tag.");
+ currentVarData.varTags.push_back(tag); // Really, only one for these
+ const void* data;
+ 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(tag, type);
+ ERRORR(rval, "Can't get tag type.");
+
+ currentVarData.varDataType = NC_DOUBLE;
+ if (MB_TYPE_INTEGER == type)
+ currentVarData.varDataType = NC_INT;
+
+ assert(0 == currentVarData.memoryHogs.size()); // Nothing so far
+ currentVarData.memoryHogs.push_back((void*)data);
+
+ if (currentVarData.varDims.empty()) {
+ // Scalar variable
+ currentVarData.writeStarts.push_back(0);
+ currentVarData.writeCounts.push_back(1);
+ }
+ else {
+ assert(1 == currentVarData.varDims.size());
+ currentVarData.writeStarts.push_back(0);
+ currentVarData.writeCounts.push_back(dimLens[currentVarData.varDims[0]]);
+ }
+ }
+ }
+
+ // Process coordinate variables
+ // Check that for used coordinates we have found the tags
+ for (std::set<std::string>::iterator setIt = usedCoordinates.begin();
+ setIt != usedCoordinates.end(); ++setIt) {
+ const std::string& coordName = *setIt;
+
+ 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;
+ 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
+
+ const void* data;
+ int 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);
+ varCoordData.writeStarts[0] = 0;
+ varCoordData.writeCounts.resize(1);
+ varCoordData.writeCounts[0] = sizeCoordinate;
+
+ // Find the type of tag, and use it
+ DataType type;
+ rval = mbImpl->tag_get_data_type(coordTag, type);
+ ERRORR(rval, "Can't get tag type.");
+
+ varCoordData.varDataType = NC_DOUBLE;
+ if (MB_TYPE_INTEGER == type)
+ varCoordData.varDataType = NC_INT;
+
+ assert(0 == varCoordData.memoryHogs.size()); // Nothing so far
+ varCoordData.memoryHogs.push_back((void*)data);
+ }
+
+ return MB_SUCCESS;
+}
+
+ErrorCode NCWriteHelper::init_file(std::vector<std::string>& var_names)
+{
+ std::vector<std::string>& dimNames = _writeNC->dimNames;
+ std::set<std::string>& usedCoordinates = _writeNC->usedCoordinates;
+ std::set<std::string>& dummyVarNames = _writeNC->dummyVarNames;
+ std::map<std::string, WriteNC::VarData>& varInfo = _writeNC->varInfo;
+ std::map<std::string, WriteNC::AttData>& globalAtts = _writeNC->globalAtts;
+ DebugOutput& dbgOut = _writeNC->dbgOut;
+
+ int tDim_in_dimNames = tDim;
+ int levDim_in_dimNames = levDim;
+
+ // 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) {
+ const std::string& coordName = *setIt;
+
+ 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;
+ varCoordData.varDims.resize(1);
+
+ /* int nc_def_dim (int ncid, const char *name, size_t len, int *dimidp);
+ * example: status = nc_def_dim(fileId, "lat", 18L, &latid);
+ */
+
+ // Actually define a dimension
+ if (NCFUNC(def_dim)(_fileId, coordName.c_str(), (size_t)varCoordData.sz,
+ &varCoordData.varDims[0]) != NC_NOERR)
+ ERRORR(MB_FAILURE, "Failed to generate dimension.");
+
+ dbgOut.tprintf(2, " for coordName %s dim id is %d \n", coordName.c_str(), (int)varCoordData.varDims[0]);
+
+ // Update tDim and levDim to actual dimension id
+ if (coordName == dimNames[tDim_in_dimNames])
+ tDim = varCoordData.varDims[0];
+ else if (coordName == dimNames[levDim_in_dimNames])
+ levDim = varCoordData.varDims[0];
+
+ // Create a variable with the same name, and its only dimension the one we just defined
+ /*
+ * int nc_def_var (int ncid, const char *name, nc_type xtype,
+ int ndims, const int dimids[], int *varidp);
+ 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.");
+
+ dbgOut.tprintf(2, " for coordName %s variable id is %d \n", coordName.c_str(), varCoordData.varId);
+ }
+
+ // Now look at requested variables, and update from the index in dimNames to the actual dimension id
+ for (size_t i = 0; i < var_names.size(); i++) {
+ std::map<std::string, WriteNC::VarData>::iterator vit = varInfo.find(var_names[i]);
+ if (vit == varInfo.end())
+ ERRORR(MB_FAILURE, "Can't find variable requested.");
+
+ // Skip coordinate variables
+ if (usedCoordinates.find(var_names[i]) != usedCoordinates.end())
+ continue;
+
+ WriteNC::VarData& variableData = vit->second;
+ int numDims = (int)variableData.varDims.size();
+ // The index is for dimNames; we need to find out the actual dimension id (from above)
+ for (int j = 0; j < numDims; j++) {
+ std::string dimName = dimNames[variableData.varDims[j]];
+ std::map<std::string, WriteNC::VarData>::iterator vit2 = varInfo.find(dimName);
+ if (vit2 == varInfo.end())
+ ERRORR(MB_FAILURE, "Can't find coordinate variable requested.");
+
+ WriteNC::VarData& coordData = vit2->second;
+ // Index in dimNames to actual dimension id
+ variableData.varDims[j] = coordData.varDims[0]; // This one, being a coordinate, is the only one
+ dbgOut.tprintf(2, " dimension with index %d name %s has ID %d \n",
+ j, dimName.c_str(), variableData.varDims[j]);
+ }
+
+ // Define the variable now:
+ if (NCFUNC(def_var)(_fileId, var_names[i].c_str(), variableData.varDataType,
+ (int)variableData.varDims.size(), &(variableData.varDims[0]),
+ &variableData.varId) != NC_NOERR)
+ ERRORR(MB_FAILURE, "Failed to create coordinate variable.");
+
+ dbgOut.tprintf(2, " for variable %s variable id is %d \n", var_names[i].c_str(), variableData.varId);
+ // Now define the variable, with all dimensions
+ }
+
+ // Define global attributes (exactly copied from the original file for the time being)
+ // Should we modify some of them (e.g. revision_Id) later?
+ std::map<std::string, WriteNC::AttData>::iterator attIt;
+ for (attIt = globalAtts.begin(); attIt != globalAtts.end(); ++attIt) {
+ const std::string& attName = attIt->first;
+ WriteNC::AttData& attData = attIt->second;
+ NCDF_SIZE& attLen = attData.attLen;
+ nc_type& attDataType = attData.attDataType;
+ const std::string& attValue = attData.attValue;
+
+ switch (attDataType) {
+ case NC_BYTE:
+ case NC_CHAR:
+ if (NC_NOERR != NCFUNC(put_att_text)(_fileId, NC_GLOBAL, attName.c_str(), attLen, attValue.c_str()))
+ ERRORR(MB_FAILURE, "Failed to define text type attribute.");
+ break;
+ case NC_DOUBLE:
+ if (NC_NOERR != NCFUNC(put_att_double)(_fileId, NC_GLOBAL, attName.c_str(), NC_DOUBLE, 1, (double*)attValue.c_str()))
+ ERRORR(MB_FAILURE, "Failed to define double type attribute.");
+ break;
+ case NC_FLOAT:
+ if (NC_NOERR != NCFUNC(put_att_float)(_fileId, NC_GLOBAL, attName.c_str(), NC_FLOAT, 1, (float*)attValue.c_str()))
+ ERRORR(MB_FAILURE, "Failed to define float type attribute.");
+ break;
+ case NC_INT:
+ if (NC_NOERR != NCFUNC(put_att_int)(_fileId, NC_GLOBAL, attName.c_str(), NC_INT, 1, (int*)attValue.c_str()))
+ ERRORR(MB_FAILURE, "Failed to define int type attribute.");
+ break;
+ case NC_SHORT:
+ if (NC_NOERR != NCFUNC(put_att_short)(_fileId, NC_GLOBAL, attName.c_str(), NC_SHORT, 1, (short*)attValue.c_str()))
+ ERRORR(MB_FAILURE, "Failed to define short type attribute.");
+ break;
+ default:
+ ERRORR(MB_FAILURE, "Unknown attribute data type.");
+ }
+ }
+
+ // Take it out of define mode
+ if (NC_NOERR != NCFUNC(enddef)(_fileId))
+ ERRORR(MB_FAILURE, "Failed to close define mode.");
+
+ return MB_SUCCESS;
+}
+
+ErrorCode ScdNCWriteHelper::collect_mesh_info()
+{
+ Interface*& mbImpl = _writeNC->mbImpl;
+ std::vector<std::string>& dimNames = _writeNC->dimNames;
+ std::vector<int>& dimLens = _writeNC->dimLens;
+
+ ErrorCode rval;
+
+ // Look for time dimension
+ std::vector<std::string>::iterator vecIt;
+ if ((vecIt = std::find(dimNames.begin(), dimNames.end(), "time")) != dimNames.end())
+ tDim = vecIt - dimNames.begin();
+ else if ((vecIt = std::find(dimNames.begin(), dimNames.end(), "t")) != dimNames.end())
+ tDim = vecIt - dimNames.begin();
+ else {
+ ERRORR(MB_FAILURE, "Couldn't find 'time' or 't' dimension.");
+ }
+ nTimeSteps = dimLens[tDim];
+
+ // Get number of levels
+ if ((vecIt = std::find(dimNames.begin(), dimNames.end(), "lev")) != dimNames.end())
+ levDim = vecIt - dimNames.begin();
+ else if ((vecIt = std::find(dimNames.begin(), dimNames.end(), "ilev")) != dimNames.end())
+ levDim = vecIt - dimNames.begin();
+ else {
+ ERRORR(MB_FAILURE, "Couldn't find 'lev' or 'ilev' dimension.");
+ }
+ nLevels = dimLens[levDim];
+
+ // __<dim_name>_LOC_MINMAX (for slon, slat, lon and lat)
+ Tag convTag = 0;
+ rval = mbImpl->tag_get_handle("__slon_LOC_MINMAX", 0, MB_TYPE_INTEGER, convTag, MB_TAG_ANY);
+ ERRORR(rval, "Trouble getting conventional tag __slon_LOC_MINMAX.");
+ int val[2];
+ rval = mbImpl->tag_get_data(convTag, &_fileSet, 1, val);
+ ERRORR(rval, "Trouble getting values for conventional tag __slon_LOC_MINMAX.");
+ lDims[0] = val[0];
+ lDims[3] = val[1];
+
+ rval = mbImpl->tag_get_handle("__slat_LOC_MINMAX", 0, MB_TYPE_INTEGER, convTag, MB_TAG_ANY);
+ ERRORR(rval, "Trouble getting conventional tag __slat_LOC_MINMAX.");
+ rval = mbImpl->tag_get_data(convTag, &_fileSet, 1, val);
+ ERRORR(rval, "Trouble getting values for conventional tag __slat_LOC_MINMAX.");
+ lDims[1] = val[0];
+ lDims[4] = val[1];
+
+ rval = mbImpl->tag_get_handle("__lon_LOC_MINMAX", 0, MB_TYPE_INTEGER, convTag, MB_TAG_ANY);
+ ERRORR(rval, "Trouble getting conventional tag __lon_LOC_MINMAX.");
+ rval = mbImpl->tag_get_data(convTag, &_fileSet, 1, val);
+ ERRORR(rval, "Trouble getting values for conventional tag __lon_LOC_MINMAX.");
+ lCDims[0] = val[0];
+ lCDims[3] = val[1];
+
+ rval = mbImpl->tag_get_handle("__lat_LOC_MINMAX", 0, MB_TYPE_INTEGER, convTag, MB_TAG_ANY);
+ ERRORR(rval, "Trouble getting conventional tag __lat_LOC_MINMAX.");
+ rval = mbImpl->tag_get_data(convTag, &_fileSet, 1, val);
+ ERRORR(rval, "Trouble getting values for conventional tag __lat_LOC_MINMAX.");
+ lCDims[1] = val[0];
+ lCDims[4] = val[1];
+
+ return MB_SUCCESS;
+}
+
+ErrorCode ScdNCWriteHelper::collect_variable_data(std::vector<std::string>& var_names)
+{
+ NCWriteHelper::collect_variable_data(var_names);
+
+ std::map<std::string, WriteNC::VarData>& varInfo = _writeNC->varInfo;
+
+ for (size_t i = 0; i < var_names.size(); i++) {
+ std::string varname = var_names[i];
+ std::map<std::string, WriteNC::VarData>::iterator vit = varInfo.find(varname);
+ if (vit == varInfo.end())
+ ERRORR(MB_FAILURE, "Can't find one variable.");
+
+ WriteNC::VarData& currentVarData = vit->second;
+ if (currentVarData.has_tsteps) {
+ // Support non-set variables with 4 dimensions like (time, lev, lat, lon)
+ assert(4 == currentVarData.varDims.size());
+
+ // Time should be the first dimension
+ assert(tDim == currentVarData.varDims[0]);
+
+ // Set up writeStarts and writeCounts
+ currentVarData.writeStarts.resize(4);
+ currentVarData.writeCounts.resize(4);
+
+ // First: time
+ currentVarData.writeStarts[0] = 0; // This value is timestep dependent, will be set later
+ currentVarData.writeCounts[0] = 1;
+
+ // Next: lev
+ currentVarData.writeStarts[1] = 0;
+ currentVarData.writeCounts[1] = currentVarData.numLev;
+
+ // Finally: lat (or slat) and lon (or slon)
+ switch (currentVarData.entLoc) {
+ case WriteNC::ENTLOCFACE:
+ // Faces
+ currentVarData.writeStarts[2] = lCDims[1];
+ currentVarData.writeCounts[2] = lCDims[4] - lCDims[1] + 1;
+ currentVarData.writeStarts[3] = lCDims[0];
+ currentVarData.writeCounts[3] = lCDims[3] - lCDims[0] + 1;
+ break;
+ default:
+ ERRORR(MB_FAILURE, "Not implemented yet.");
+ }
+ }
+
+ // Get variable size
+ currentVarData.sz = 1;
+ for (std::size_t idx = 0; idx != currentVarData.writeCounts.size(); idx++)
+ currentVarData.sz *= currentVarData.writeCounts[idx];
+ }
+
+ return MB_SUCCESS;
+}
+
+ErrorCode ScdNCWriteHelper::write_values(std::vector<std::string>& var_names)
+{
+ 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;
+
+ ErrorCode rval;
+ int success;
+
+ // Now look at requested var_names; if they have time, we will have a list, and write one at a time
+ // If not, just write regularly
+ // For parallel write, we assume that the data ranges do not overlap across processors (this is true
+ // for Euler and FV variables on quads)
+ // While overlapped writing might still work, we should better not take that risk
+ // Use collective I/O mode put (synchronous write) for the time being, we can try nonblocking put
+ // (request aggregation) later
+ for (size_t i = 0; i < var_names.size(); i++) {
+ std::map<std::string, WriteNC::VarData>::iterator vit = varInfo.find(var_names[i]);
+ if (vit == varInfo.end())
+ ERRORR(MB_FAILURE, "Can't find variable requested.");
+
+ WriteNC::VarData& variableData = vit->second;
+ int numTimeSteps = (int)variableData.varTags.size();
+ if (variableData.has_tsteps) {
+ // Time should be the first dimension
+ assert(tDim == variableData.varDims[0]);
+
+ // Get entities of this variable
+ Range ents;
+ switch (variableData.entLoc) {
+ case WriteNC::ENTLOCFACE:
+ // Faces
+ rval = mbImpl->get_entities_by_dimension(_fileSet, 2, ents);
+ ERRORR(rval, "Can't get entities for faces.");
+ break;
+ default:
+ ERRORR(MB_FAILURE, "Not implemented yet.");
+ }
+
+ // 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
+ size_t nj = variableData.writeCounts[2]; // lat
+ size_t nk = variableData.writeCounts[1]; // lev
+
+ variableData.writeCounts[0] = 1; // We will write one time step
+ // FIXME: Should use tstep_nums (from writing options) later
+ for (int j = 0; j < numTimeSteps; j++) {
+ // We will write one time step, and count will be one; start will be different
+ // We will write values directly from tag_iterate, but we should also transpose for level
+ // so that means deep copy for transpose
+ variableData.writeStarts[0] = j; // This is time, again
+ int count;
+ void* dataptr;
+ rval = mbImpl->tag_iterate(variableData.varTags[j], ents.begin(), ents.end(), count, dataptr);
+ assert(count == (int)ents.size());
+
+ // Now write from memory directly
+ switch (variableData.varDataType) {
+ case NC_DOUBLE: {
+ std::vector<double> tmpdoubledata(ni*nj*nk);
+ // Transpose (lat, lon, lev) back to (lev, lat, lon)
+ jik_to_kji(ni, nj, nk, &tmpdoubledata[0], (double*)(dataptr));
+ success = NCFUNCAP(_vara_double)(_fileId, variableData.varId,
+ &variableData.writeStarts[0], &variableData.writeCounts[0],
+ &tmpdoubledata[0]);
+ ERRORS(success, "Failed to write double data.");
+ break;
+ }
+ default:
+ ERRORR(MB_FAILURE, "Not implemented yet.");
+ }
+ }
+ }
+ else {
+ switch (variableData.varDataType) {
+ case NC_DOUBLE:
+ success = NCFUNCAP(_vara_double)(_fileId, variableData.varId, &variableData.writeStarts[0],
+ &variableData.writeCounts[0], (double*)(variableData.memoryHogs[0]));
+ ERRORS(success, "Failed to write double data.");
+ break;
+ default:
+ ERRORR(MB_FAILURE, "Not implemented yet.");
+ }
+ }
+ }
+
+ // Write coordinates used by requested var_names
+ // Use independent I/O mode put, since this write is only for the root processor
+ // CAUTION: if the NetCDF ID is from a previous call to ncmpi_create rather than ncmpi_open,
+ // all processors need to call ncmpi_begin_indep_data(). If only the root processor does so,
+ // ncmpi_begin_indep_data() call will be blocked forever :(
+#ifdef PNETCDF_FILE
+ // Enter independent I/O mode
+ success = NCFUNC(begin_indep_data)(_fileId);
+ ERRORS(success, "Failed to begin independent I/O mode.");
+#endif
+
+ int rank = 0;
+#ifdef USE_MPI
+ bool& isParallel = _writeNC->isParallel;
+ if (isParallel) {
+ ParallelComm*& myPcomm = _writeNC->myPcomm;
+ rank = myPcomm->proc_config().proc_rank();
+ }
+#endif
+ if (0 == rank) {
+ for (std::set<std::string>::iterator setIt = usedCoordinates.begin();
+ setIt != usedCoordinates.end(); ++setIt) {
+ const std::string& coordName = *setIt;
+
+ // 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.");
+
+ WriteNC::VarData& varCoordData = vit->second;
+
+ switch (varCoordData.varDataType) {
+ case NC_DOUBLE:
+ // Independent I/O mode put
+ success = NCFUNCP(_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:
+ // Independent I/O mode put
+ success = NCFUNCP(_vara_int)(_fileId, varCoordData.varId, &varCoordData.writeStarts[0],
+ &varCoordData.writeCounts[0], (int*)(varCoordData.memoryHogs[0]));
+ ERRORS(success, "Failed to write int data.");
+ break;
+ default:
+ ERRORR(MB_FAILURE, "Not implemented yet.");
+ }
+ }
+ }
+
+#ifdef PNETCDF_FILE
+ // End independent I/O mode
+ success = NCFUNC(end_indep_data)(_fileId);
+ ERRORS(success, "Failed to end independent I/O mode.");
+#endif
+
+ return MB_SUCCESS;
+}
+
+} /* namespace moab */
diff --git a/src/io/NCWriteHelper.hpp b/src/io/NCWriteHelper.hpp
new file mode 100644
index 0000000..b05fe12
--- /dev/null
+++ b/src/io/NCWriteHelper.hpp
@@ -0,0 +1,125 @@
+/*
+ * NCWriteHelper.hpp
+ *
+ * Purpose : Climate NC writer file helper; abstract, will be implemented for each type
+ *
+ * Created on: Mar 28, 2014
+ */
+
+#ifndef NCWRITEHELPER_HPP_
+#define NCWRITEHELPER_HPP_
+#include "WriteNC.hpp"
+
+namespace moab {
+
+class NCWriteHelper
+{
+public:
+ NCWriteHelper(WriteNC* writeNC, int fileId, const FileOptions& opts, EntityHandle fileSet)
+: _writeNC(writeNC), _fileId(fileId), _opts(opts), _fileSet(fileSet),
+ nTimeSteps(0), nLevels(1), tDim(-1), levDim(-1) {}
+ virtual ~NCWriteHelper() {};
+
+ //! Get appropriate helper instance for WriteNC class based on some info in the file set
+ static NCWriteHelper* get_nc_helper(WriteNC* writeNC, int fileId, const FileOptions& opts, EntityHandle fileSet);
+
+ //! Collect necessary info about local mesh
+ virtual ErrorCode collect_mesh_info() = 0;
+
+ //! Collect data for specified variables
+ virtual ErrorCode collect_variable_data(std::vector<std::string>& var_names);
+
+ //! Take the info from VarData and write first the coordinates, then the actual variables
+ virtual ErrorCode write_values(std::vector<std::string>& var_names) = 0;
+
+ //! Initialize file: this is where all defines are done
+ //! The VarData dimension ids are filled up after define
+ ErrorCode init_file(std::vector<std::string>& var_names);
+
+protected:
+ template <typename T> void jik_to_kji(size_t ni, size_t nj, size_t nk, T* dest, T* source)
+ {
+ size_t nik = ni * nk, nij = ni * nj;
+ for (std::size_t k = 0; k != nk; k++)
+ for (std::size_t j = 0; j != nj; j++)
+ for (std::size_t i = 0; i != ni; i++)
+ dest[k*nij + j*ni + i] = source[j*nik + i*nk + k];
+ }
+
+protected:
+ //! Allow NCWriteHelper to directly access members of WriteNC
+ WriteNC* _writeNC;
+
+ //! Cache some information from ReadNC
+ int _fileId;
+ const FileOptions& _opts;
+ EntityHandle _fileSet;
+
+ //! Dimensions of time and level
+ int nTimeSteps, nLevels;
+
+ //! Dimension numbers for time and level
+ int tDim, levDim;
+};
+
+//! Child helper class for scd mesh, e.g. CAM_EL or CAM_FV
+class ScdNCWriteHelper : public NCWriteHelper
+{
+public:
+ ScdNCWriteHelper(WriteNC* writeNC, int fileId, const FileOptions& opts, EntityHandle fileSet)
+: NCWriteHelper(writeNC, fileId, opts, fileSet)
+ {
+ for (unsigned int i = 0; i < 6; i++) {
+ lDims[i] = -1;
+ lCDims[i] = -1;
+ }
+ }
+ virtual ~ScdNCWriteHelper() {}
+
+private:
+ //! Implementation of NCWriteHelper::collect_mesh_info()
+ virtual ErrorCode collect_mesh_info();
+
+ //! Collect data for specified variables
+ virtual ErrorCode collect_variable_data(std::vector<std::string>& var_names);
+
+ //! Implementation of NCWriteHelper::write_values()
+ virtual ErrorCode write_values(std::vector<std::string>& var_names);
+
+protected:
+ //! Dimensions of my local part of grid
+ int lDims[6];
+
+ //! Center dimensions of my local part of grid
+ int lCDims[6];
+};
+
+//! Child helper class for ucd mesh, e.g. CAM_SE (HOMME) or MPAS
+class UcdNCWriteHelper : public NCWriteHelper
+{
+public:
+ UcdNCWriteHelper(WriteNC* writeNC, int fileId, const FileOptions& opts, EntityHandle fileSet)
+: NCWriteHelper(writeNC, fileId, opts, fileSet),
+ nLocalCellsOwned(0), nLocalEdgesOwned(0), nLocalVerticesOwned(0),
+ cDim(-1), eDim(-1), vDim(-1) {}
+ virtual ~UcdNCWriteHelper() {}
+
+protected:
+ //! Dimensions of my local owned part of grid
+ int nLocalCellsOwned;
+ int nLocalEdgesOwned;
+ int nLocalVerticesOwned;
+
+ //! Dimension numbers for nCells, nEdges and nVertices
+ int cDim, eDim, vDim;
+
+ //! Local owned cells, edges and vertices
+ Range localCellsOwned, localEdgesOwned, localVertsOwned;
+
+ //! Local global ID for owned cells, edges and vertices
+ Range localGidCellsOwned, localGidEdgesOwned, localGidVertsOwned;
+};
+
+} // namespace moab
+
+#endif
diff --git a/src/io/NCWriteMPAS.cpp b/src/io/NCWriteMPAS.cpp
new file mode 100644
index 0000000..72dc35f
--- /dev/null
+++ b/src/io/NCWriteMPAS.cpp
@@ -0,0 +1,33 @@
+/*
+ * NCWriteMPAS.cpp
+ *
+ * Created on: April 9, 2014
+ */
+
+#include "NCWriteMPAS.hpp"
+#include "moab/WriteUtilIface.hpp"
+
+#define ERRORR(rval, str) \
+ if (MB_SUCCESS != rval) { _writeNC->mWriteIface->report_error("%s", str); return rval; }
+
+#define ERRORS(err, str) \
+ if (err) { _writeNC->mWriteIface->report_error("%s", str); return MB_FAILURE; }
+
+namespace moab {
+
+NCWriteMPAS::~NCWriteMPAS()
+{
+ // TODO Auto-generated destructor stub
+}
+
+ErrorCode NCWriteMPAS::collect_mesh_info()
+{
+ return MB_NOT_IMPLEMENTED;
+}
+
+ErrorCode NCWriteMPAS::write_values(std::vector<std::string>& /* var_names */)
+{
+ return MB_NOT_IMPLEMENTED;
+}
+
+} /* namespace moab */
diff --git a/src/io/NCWriteMPAS.hpp b/src/io/NCWriteMPAS.hpp
new file mode 100644
index 0000000..702984f
--- /dev/null
+++ b/src/io/NCWriteMPAS.hpp
@@ -0,0 +1,34 @@
+/*
+ * NCWriteMPAS.hpp
+ *
+ * nc write helper for MPAS type data (CAM)
+ * Created on: April 9, 2014
+ *
+ */
+
+#ifndef NCWRITEMPAS_HPP_
+#define NCWRITEMPAS_HPP_
+
+#include "NCWriteHelper.hpp"
+
+namespace moab {
+
+class NCWriteMPAS: public UcdNCWriteHelper
+{
+public:
+ NCWriteMPAS(WriteNC* writeNC, int fileId, const FileOptions& opts, EntityHandle fileSet)
+: UcdNCWriteHelper(writeNC, fileId, opts, fileSet) {}
+
+ virtual ~NCWriteMPAS();
+
+private:
+ //! Implementation of NCWriteHelper::collect_mesh_info()
+ virtual ErrorCode collect_mesh_info();
+
+ //! Collect data for specified variables
+ virtual ErrorCode write_values(std::vector<std::string>& var_names);
+};
+
+} // namespace moab
+
+#endif
diff --git a/src/io/ReadNC.cpp b/src/io/ReadNC.cpp
index a0e02b4..0a41cfa 100644
--- a/src/io/ReadNC.cpp
+++ b/src/io/ReadNC.cpp
@@ -118,6 +118,28 @@ ErrorCode ReadNC::load_file(const char* file_name, const EntityHandle* file_set,
ERRORR(rval, "Trouble checking mesh from last read.\n");
}
+ // Create some conventional tags, e.g. __NUM_DIMS
+ // Keep a flag on the file set to prevent conventional tags from being created again on a second read
+ Tag convTagsCreated = 0;
+ int def_val = 0;
+ rval = mbImpl->tag_get_handle("__CONV_TAGS_CREATED", 1, MB_TYPE_INTEGER, convTagsCreated,
+ MB_TAG_SPARSE | MB_TAG_CREAT, &def_val);
+ ERRORR(rval, "Trouble getting _CONV_TAGS_CREATED tag.");
+ int create_conv_tags_flag = 0;
+ rval = mbImpl->tag_get_data(convTagsCreated, &tmp_set, 1, &create_conv_tags_flag);
+ if (0 == create_conv_tags_flag) {
+ // Read dimension variables to create tags like __<var_name>_DIMS
+ rval = myHelper->read_variables(dimNames, tstep_nums);
+ ERRORR(rval, "Trouble reading dimension variables.");
+
+ rval = myHelper->create_conventional_tags(tstep_nums);
+ ERRORR(rval, "Trouble creating NC conventional tags.");
+
+ create_conv_tags_flag = 1;
+ rval = mbImpl->tag_set_data(convTagsCreated, &tmp_set, 1, &create_conv_tags_flag);
+ ERRORR(rval, "Trouble setting data for _CONV_TAGS_CREATED tag.");
+ }
+
// Create mesh vertex/edge/face sequences
Range faces;
if (!noMesh) {
@@ -125,26 +147,10 @@ ErrorCode ReadNC::load_file(const char* file_name, const EntityHandle* file_set,
ERRORR(rval, "Trouble creating mesh.");
}
- // Read variables onto grid
+ // Read specified variables onto grid
if (!noVars) {
rval = myHelper->read_variables(var_names, tstep_nums);
- if (MB_FAILURE == rval)
- return rval;
- }
- else {
- // Read dimension variables by default (the dimensions that are also variables)
- std::vector<std::string> dim_var_names;
- for (unsigned int i = 0; i < dimNames.size(); i++) {
- std::map<std::string, VarData>::iterator mit = varInfo.find(dimNames[i]);
- if (mit != varInfo.end())
- dim_var_names.push_back(dimNames[i]);
- }
-
- if (!dim_var_names.empty()) {
- rval = myHelper->read_variables(dim_var_names, tstep_nums);
- if (MB_FAILURE == rval)
- return rval;
- }
+ ERRORR(rval, "Trouble reading specified variables.");
}
#ifdef USE_MPI
@@ -175,12 +181,6 @@ ErrorCode ReadNC::load_file(const char* file_name, const EntityHandle* file_set,
}
#endif
- // Create NC conventional tags when loading header info only
- if (noMesh && noVars) {
- rval = myHelper->create_conventional_tags(tstep_nums);
- ERRORR(rval, "Trouble creating NC conventional tags.");
- }
-
mbImpl->release_interface(scdi);
scdi = NULL;
This diff is so big that we needed to truncate the remainder.
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