[MOAB-dev] commit/MOAB: danwu: Updated NC writer to support user-created non-set variables that might not have time or level dimension.

commits-noreply at bitbucket.org commits-noreply at bitbucket.org
Mon May 12 17:47:13 CDT 2014


1 new commit in MOAB:

https://bitbucket.org/fathomteam/moab/commits/96f41d5bc215/
Changeset:   96f41d5bc215
Branch:      master
User:        danwu
Date:        2014-05-13 00:47:00
Summary:     Updated NC writer to support user-created non-set variables that might not have time or level dimension.

Affected #:  4 files

diff --git a/src/io/NCWriteGCRM.cpp b/src/io/NCWriteGCRM.cpp
index 36160d2..b9f412c 100644
--- a/src/io/NCWriteGCRM.cpp
+++ b/src/io/NCWriteGCRM.cpp
@@ -130,69 +130,94 @@ ErrorCode NCWriteGCRM::collect_variable_data(std::vector<std::string>& var_names
       ERRORR(MB_FAILURE, "Can't find one variable.");
 
     WriteNC::VarData& currentVarData = vit->second;
+    std::vector<int>& varDims = currentVarData.varDims;
 
     // Skip edge variables, if there are no edges
     if (localEdgesOwned.empty() && currentVarData.entLoc == WriteNC::ENTLOCEDGE)
       continue;
 
-    std::vector<int>& varDims = currentVarData.varDims;
+    // Skip set variables, which were already processed in NCWriteHelper::collect_variable_data()
+    if (WriteNC::ENTLOCSET == currentVarData.entLoc)
+      continue;
 
+    // Set up writeStarts and writeCounts (maximum number of dimensions is 3)
+    currentVarData.writeStarts.resize(3);
+    currentVarData.writeCounts.resize(3);
+    unsigned int dim_idx = 0;
+
+    // First: time
     if (currentVarData.has_tsteps) {
-      // Support non-set variables with 3 dimensions like (Time, cells, layers), or
-      // 2 dimensions like (Time, cells)
+      // Non-set variables with timesteps
+      // 3 dimensions like (time, cells, layers)
+      // 2 dimensions like (time, cells)
       assert(3 == varDims.size() || 2 == varDims.size());
 
       // Time should be the first dimension
       assert(tDim == 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: nVertices / cells / nEdges
-      switch (currentVarData.entLoc) {
-        case WriteNC::ENTLOCVERT:
-          // Vertices
-          // Start from the first localGidVerts
-          // Actually, this will be reset later for writing
-          currentVarData.writeStarts[1] = localGidVertsOwned[0] - 1;
-          currentVarData.writeCounts[1] = localGidVertsOwned.size();
-          break;
-        case WriteNC::ENTLOCFACE:
-          // Faces
-          // Start from the first localGidCells
-          // Actually, this will be reset later for writing
-          currentVarData.writeStarts[1] = localGidCellsOwned[0] - 1;
-          currentVarData.writeCounts[1] = localGidCellsOwned.size();
-          break;
-        case WriteNC::ENTLOCEDGE:
-          // Edges
-          // Start from the first localGidEdges
-          // Actually, this will be reset later for writing
-          currentVarData.writeStarts[1] = localGidEdgesOwned[0] - 1;
-          currentVarData.writeCounts[1] = localGidEdgesOwned.size();
-          break;
-        default:
-          ERRORR(MB_FAILURE, "Unexpected entity location type for GCRM non-set variable.");
-      }
+      currentVarData.writeStarts[dim_idx] = 0; // This value is timestep dependent, will be set later
+      currentVarData.writeCounts[dim_idx] = 1;
+      dim_idx++;
+    }
+    else {
+      // Non-set variables without timesteps
+      // 2 dimensions like (cells, layers)
+      // 1 dimension like (cells)
+      assert(2 == varDims.size() || 1 == varDims.size());
+    }
 
-      // Finally: layers or other optional levels, it is possible that there is no
-      // level dimension (numLev is 0) for this non-set variable
-      if (currentVarData.numLev < 1)
-        currentVarData.numLev = 1;
-      currentVarData.writeStarts[2] = 0;
-      currentVarData.writeCounts[2] = currentVarData.numLev;
+    // Next: corners / cells / edges
+    switch (currentVarData.entLoc) {
+      case WriteNC::ENTLOCVERT:
+        // Vertices
+        // Start from the first localGidVerts
+        // Actually, this will be reset later for writing
+        currentVarData.writeStarts[dim_idx] = localGidVertsOwned[0] - 1;
+        currentVarData.writeCounts[dim_idx] = localGidVertsOwned.size();
+        break;
+      case WriteNC::ENTLOCFACE:
+        // Faces
+        // Start from the first localGidCells
+        // Actually, this will be reset later for writing
+        currentVarData.writeStarts[dim_idx] = localGidCellsOwned[0] - 1;
+        currentVarData.writeCounts[dim_idx] = localGidCellsOwned.size();
+        break;
+      case WriteNC::ENTLOCEDGE:
+        // Edges
+        // Start from the first localGidEdges
+        // Actually, this will be reset later for writing
+        currentVarData.writeStarts[dim_idx] = localGidEdgesOwned[0] - 1;
+        currentVarData.writeCounts[dim_idx] = localGidEdgesOwned.size();
+        break;
+      default:
+        ERRORR(MB_FAILURE, "Unexpected entity location type for GCRM non-set variable.");
+    }
+    dim_idx++;
+
+    // Finally: layers or other optional levels, it is possible that there is no
+    // level dimension (numLev is 0) for this non-set variable
+    if (currentVarData.numLev > 0) {
+      // Non-set variables with levels
+      // 3 dimensions like (time, cells, layers)
+      // 2 dimensions like (cells, layers)
+      assert(3 == varDims.size() || 2 == varDims.size());
+
+      currentVarData.writeStarts[dim_idx] = 0;
+      currentVarData.writeCounts[dim_idx] = currentVarData.numLev;
+      dim_idx++;
+    }
+    else {
+      // Non-set variables without levels
+      // 2 dimensions like (time, cells)
+      // 1 dimension like (cells)
+      assert(2 == varDims.size() || 1 == varDims.size());
     }
 
     // Get variable size
     currentVarData.sz = 1;
-    for (std::size_t idx = 0; idx != 3; idx++)
+    for (std::size_t idx = 0; idx < dim_idx; idx++)
       currentVarData.sz *= currentVarData.writeCounts[idx];
-  }
+  } // for (size_t i = 0; i < var_names.size(); i++)
 
   return MB_SUCCESS;
 }
@@ -202,8 +227,6 @@ ErrorCode NCWriteGCRM::write_nonset_variables(std::vector<WriteNC::VarData>& vda
   Interface*& mbImpl = _writeNC->mbImpl;
 
   int success;
-  Range* pLocalEntsOwned = NULL;
-  Range* pLocalGidEntsOwned = NULL;
 
   // For each indexed variable tag, write a time step data
   for (unsigned int i = 0; i < vdatas.size(); i++) {
@@ -213,10 +236,9 @@ ErrorCode NCWriteGCRM::write_nonset_variables(std::vector<WriteNC::VarData>& vda
     if (localEdgesOwned.empty() && variableData.entLoc == WriteNC::ENTLOCEDGE)
       continue;
 
-    // Time should be the first dimension
-    assert(tDim == variableData.varDims[0]);
-
     // Get local owned entities of this variable
+    Range* pLocalEntsOwned = NULL;
+    Range* pLocalGidEntsOwned = NULL;
     switch (variableData.entLoc) {
       case WriteNC::ENTLOCVERT:
         // Vertices
@@ -237,15 +259,45 @@ ErrorCode NCWriteGCRM::write_nonset_variables(std::vector<WriteNC::VarData>& vda
         ERRORR(MB_FAILURE, "Unexpected entity location type for GCRM non-set variable.");
     }
 
-    // A typical GCRM non-set variable has 3 dimensions as (time, cells, layers)
+    unsigned int num_timesteps;
+    unsigned int ents_idx = 0;
+    if (variableData.has_tsteps) {
+      // Non-set variables with timesteps
+      // 3 dimensions like (time, cells, layers)
+      // 2 dimensions like (time, cells)
+      num_timesteps = tstep_nums.size();
+      ents_idx++;
+    }
+    else {
+      // Non-set variables without timesteps
+      // 2 dimensions like (cells, layers)
+      // 1 dimension like (cells)
+      num_timesteps = 1;
+    }
+
+    unsigned int num_lev;
+    if (variableData.numLev > 0) {
+      // Non-set variables with levels
+      // 3 dimensions like (time, cells, layers)
+      // 2 dimensions like (cells, layers)
+      num_lev = variableData.numLev;
+    }
+    else {
+      // Non-set variables without levels
+      // 2 dimensions like (time, cells)
+      // 1 dimension like (cells)
+      num_lev = 1;
+    }
+
     for (unsigned int t = 0; t < tstep_nums.size(); t++) {
       // We will write one time step, and count will be one; start will be different
       // Use tag_get_data instead of tag_iterate to copy tag data, as localEntsOwned
       // might not be contiguous.
-      variableData.writeStarts[0] = t; // This is start for time
-      std::vector<double> tag_data(pLocalEntsOwned->size() * variableData.numLev);
+      if (tDim == variableData.varDims[0])
+        variableData.writeStarts[0] = t; // This is start for time
+      std::vector<double> tag_data(pLocalEntsOwned->size() * num_lev);
       ErrorCode rval = mbImpl->tag_get_data(variableData.varTags[t], *pLocalEntsOwned, &tag_data[0]);
-      ERRORR(rval, "Trouble getting tag data on owned vertices.");
+      ERRORR(rval, "Trouble getting tag data on owned entities.");
 
 #ifdef PNETCDF_FILE
       size_t nb_writes = pLocalGidEntsOwned->psize();
@@ -263,8 +315,8 @@ ErrorCode NCWriteGCRM::write_nonset_variables(std::vector<WriteNC::VarData>& vda
               pair_iter != pLocalGidEntsOwned->pair_end(); ++pair_iter, ic++) {
             EntityHandle starth = pair_iter->first;
             EntityHandle endh = pair_iter->second;
-            variableData.writeStarts[1] = (NCDF_SIZE)(starth - 1);
-            variableData.writeCounts[1] = (NCDF_SIZE)(endh - starth + 1);
+            variableData.writeStarts[ents_idx] = (NCDF_SIZE)(starth - 1);
+            variableData.writeCounts[ents_idx] = (NCDF_SIZE)(endh - starth + 1);
 
             // Do a partial write, in each subrange
 #ifdef PNETCDF_FILE
@@ -280,7 +332,7 @@ ErrorCode NCWriteGCRM::write_nonset_variables(std::vector<WriteNC::VarData>& vda
             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;
+            indexInDoubleArray += (endh - starth + 1) * num_lev;
           }
           assert(ic == pLocalGidEntsOwned->psize());
 #ifdef PNETCDF_FILE

diff --git a/src/io/NCWriteHOMME.cpp b/src/io/NCWriteHOMME.cpp
index d7b865d..7adb909 100644
--- a/src/io/NCWriteHOMME.cpp
+++ b/src/io/NCWriteHOMME.cpp
@@ -94,42 +94,73 @@ ErrorCode NCWriteHOMME::collect_variable_data(std::vector<std::string>& var_name
       ERRORR(MB_FAILURE, "Can't find one variable.");
 
     WriteNC::VarData& currentVarData = vit->second;
+    std::vector<int>& varDims = currentVarData.varDims;
+
+    // Skip set variables, which were already processed in NCWriteHelper::collect_variable_data()
+    if (WriteNC::ENTLOCSET == currentVarData.entLoc)
+      continue;
+
+    // Set up writeStarts and writeCounts (maximum number of dimensions is 3)
+    currentVarData.writeStarts.resize(3);
+    currentVarData.writeCounts.resize(3);
+    unsigned int dim_idx = 0;
+
+    // First: time
     if (currentVarData.has_tsteps) {
-      // Support non-set variables with 3 dimensions like (time, lev, ncol)
-      assert(3 == currentVarData.varDims.size());
+      // Non-set variables with timesteps
+      // 3 dimensions like (time, lev, ncol)
+      // 2 dimensions like (time, ncol)
+      assert(3 == varDims.size() || 2 == 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] = localGidVertsOwned.size();
-          break;
-        default:
-          ERRORR(MB_FAILURE, "Unexpected entity location type for HOMME non-set variable.");
-      }
+      assert(tDim == varDims[0]);
+
+      currentVarData.writeStarts[dim_idx] = 0; // This value is timestep dependent, will be set later
+      currentVarData.writeCounts[dim_idx] = 1;
+      dim_idx++;
+    }
+    else {
+      // Non-set variables without timesteps
+      // 2 dimensions like (lev, ncol)
+      // 1 dimension like (ncol)
+      assert(2 == varDims.size() || 1 == varDims.size());
+    }
+
+    // Next: lev
+    if (currentVarData.numLev > 0) {
+      // Non-set variables with levels
+      // 3 dimensions like (time, lev, ncol)
+      // 2 dimensions like (lev, ncol)
+      assert(3 == varDims.size() || 2 == varDims.size());
+
+      currentVarData.writeStarts[dim_idx] = 0;
+      currentVarData.writeCounts[dim_idx] = currentVarData.numLev;
+      dim_idx++;
+    }
+    else {
+      // Non-set variables without levels
+      // 2 dimensions like (time, ncol)
+      // 1 dimension like (ncol)
+      assert(2 == varDims.size() || 1 == varDims.size());
+    }
+
+    // Finally: ncol
+    switch (currentVarData.entLoc) {
+      case WriteNC::ENTLOCVERT:
+        // Vertices
+        // Start from the first localGidVerts
+        // Actually, this will be reset later for writing
+        currentVarData.writeStarts[dim_idx] = localGidVertsOwned[0] - 1;
+        currentVarData.writeCounts[dim_idx] = localGidVertsOwned.size();
+        break;
+      default:
+        ERRORR(MB_FAILURE, "Unexpected entity location type for HOMME non-set variable.");
     }
+    dim_idx++;
 
     // Get variable size
     currentVarData.sz = 1;
-    for (std::size_t idx = 0; idx != currentVarData.writeCounts.size(); idx++)
+    for (std::size_t idx = 0; idx < dim_idx; idx++)
       currentVarData.sz *= currentVarData.writeCounts[idx];
   }
 
@@ -147,9 +178,6 @@ ErrorCode NCWriteHOMME::write_nonset_variables(std::vector<WriteNC::VarData>& vd
   for (unsigned int i = 0; i < vdatas.size(); i++) {
     WriteNC::VarData& variableData = vdatas[i];
 
-    // 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:
@@ -159,16 +187,47 @@ ErrorCode NCWriteHOMME::write_nonset_variables(std::vector<WriteNC::VarData>& vd
         ERRORR(MB_FAILURE, "Unexpected entity location type for HOMME non-set variable.");
     }
 
-    // A typical HOMME non-set variable has 3 dimensions as (time, lev, ncol)
+    unsigned int num_timesteps;
+    unsigned int ncol_idx = 0;
+    if (variableData.has_tsteps) {
+      // Non-set variables with timesteps
+      // 3 dimensions like (time, lev, ncol)
+      // 2 dimensions like (time, ncol)
+      num_timesteps = tstep_nums.size();
+      ncol_idx++;
+    }
+    else {
+      // Non-set variables without timesteps
+      // 2 dimensions like (lev, ncol)
+      // 1 dimension like (ncol)
+      num_timesteps = 1;
+    }
+
+    unsigned int num_lev;
+    if (variableData.numLev > 0) {
+      // Non-set variables with levels
+      // 3 dimensions like (time, lev, ncol)
+      // 2 dimensions like (lev, ncol)
+      num_lev = variableData.numLev;
+      ncol_idx++;
+    }
+    else {
+      // Non-set variables without levels
+      // 2 dimensions like (time, ncol)
+      // 1 dimension like (ncol)
+      num_lev = 1;
+    }
+
     // At each timestep, we need to transpose tag format (ncol, lev) back
     // to NC format (lev, ncol) for writing
-    for (unsigned int t = 0; t < tstep_nums.size(); t++) {
+    for (unsigned int t = 0; t < num_timesteps; t++) {
       // We will write one time step, and count will be one; start will be different
       // Use tag_get_data instead of tag_iterate to copy tag data, as localVertsOwned
       // might not be contiguous. We should also transpose for level so that means
       // deep copy for transpose
-      variableData.writeStarts[0] = t; // This is start for time
-      std::vector<double> tag_data(num_local_verts_owned * variableData.numLev);
+      if (tDim == variableData.varDims[0])
+        variableData.writeStarts[0] = t; // This is start for time
+      std::vector<double> tag_data(num_local_verts_owned * num_lev);
       ErrorCode rval = mbImpl->tag_get_data(variableData.varTags[t], localVertsOwned, &tag_data[0]);
       ERRORR(rval, "Trouble getting tag data on owned vertices.");
 
@@ -182,9 +241,10 @@ ErrorCode NCWriteHOMME::write_nonset_variables(std::vector<WriteNC::VarData>& vd
       // Use nonblocking put (request aggregation)
       switch (variableData.varDataType) {
         case NC_DOUBLE: {
-          std::vector<double> tmpdoubledata(num_local_verts_owned * variableData.numLev);
-          // Transpose (ncol, lev) back to (lev, ncol)
-          jik_to_kji(num_local_verts_owned, 1, variableData.numLev, &tmpdoubledata[0], &tag_data[0]);
+          std::vector<double> tmpdoubledata(num_local_verts_owned * num_lev);
+          if (num_lev > 1)
+            // Transpose (ncol, lev) back to (lev, ncol)
+            jik_to_kji(num_local_verts_owned, 1, num_lev, &tmpdoubledata[0], &tag_data[0]);
 
           size_t indexInDoubleArray = 0;
           size_t ic = 0;
@@ -192,8 +252,8 @@ ErrorCode NCWriteHOMME::write_nonset_variables(std::vector<WriteNC::VarData>& vd
               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);
+            variableData.writeStarts[ncol_idx] = (NCDF_SIZE)(starth - 1);
+            variableData.writeCounts[ncol_idx] = (NCDF_SIZE)(endh - starth + 1);
 
             // Do a partial write, in each subrange
 #ifdef PNETCDF_FILE
@@ -209,7 +269,7 @@ ErrorCode NCWriteHOMME::write_nonset_variables(std::vector<WriteNC::VarData>& vd
             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;
+            indexInDoubleArray += (endh - starth + 1) * num_lev;
           }
           assert(ic == localGidVertsOwned.psize());
 #ifdef PNETCDF_FILE

diff --git a/src/io/NCWriteHelper.cpp b/src/io/NCWriteHelper.cpp
index b315637..635b595 100644
--- a/src/io/NCWriteHelper.cpp
+++ b/src/io/NCWriteHelper.cpp
@@ -95,68 +95,88 @@ ErrorCode NCWriteHelper::collect_variable_data(std::vector<std::string>& var_nam
     if ((std::find(currentVarData.varDims.begin(), currentVarData.varDims.end(), levDim) != currentVarData.varDims.end()))
       currentVarData.numLev = nLevels;
 
-    // Process non-coordinate variables with time steps
-    if (currentVarData.has_tsteps) {
-      // It might be a set variable, e.g. xtime(Time) or xtime(Time, StrLen)
-      for (unsigned int t = 0; t < tstep_nums.size(); t++) {
-        Tag indexedTag = 0;
-        std::stringstream ssTagNameWithIndex;
-        ssTagNameWithIndex << varname << tstep_nums[t];
-        rval = mbImpl->tag_get_handle(ssTagNameWithIndex.str().c_str(), indexedTag);
+    // Process set variables
+    if (WriteNC::ENTLOCSET == currentVarData.entLoc) {
+      if (currentVarData.has_tsteps) {
+        // Set variables with timesteps, e.g. xtime(Time) or xtime(Time, StrLen)
+        // TBD
+      }
+      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.");
-        dbgOut.tprintf(2, "    found indexed tag %d with name %s\n", tstep_nums[t], ssTagNameWithIndex.str().c_str());
-        currentVarData.varTags.push_back(indexedTag);
+        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.");
 
-        // The type of the tag is fixed though
+        // Find the type of tag, and use it
         DataType type;
-        rval = mbImpl->tag_get_data_type(indexedTag, 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 {
+          for (size_t j = 0; j < currentVarData.varDims.size(); j++) {
+            currentVarData.writeStarts.push_back(0);
+            currentVarData.writeCounts.push_back(dimLens[currentVarData.varDims[j]]);
+          }
+        }
+
+        // Get variable size
+        currentVarData.sz = 1;
+        for (std::size_t idx = 0; idx != currentVarData.writeCounts.size(); idx++)
+          currentVarData.sz *= currentVarData.writeCounts[idx];
       }
-    }
-    // Process non-coordinate variables without time steps
+    } // if (WriteNC::ENTLOCSET == currentVarData.entLoc)
+    // Process non-set variables
     else {
-      // Should be a set variable
-      assert(WriteNC::ENTLOCSET == currentVarData.entLoc);
-
-      // 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
+      Tag indexedTag = 0;
+
+      if (currentVarData.has_tsteps) {
+        for (unsigned int t = 0; t < tstep_nums.size(); t++) {
+          std::stringstream ssTagNameWithIndex;
+          ssTagNameWithIndex << varname << tstep_nums[t];
+          rval = mbImpl->tag_get_handle(ssTagNameWithIndex.str().c_str(), indexedTag);
+          ERRORR(rval, "Can't find one tag.");
+          dbgOut.tprintf(2, "    found indexed tag %d with name %s\n", tstep_nums[t], ssTagNameWithIndex.str().c_str());
+          currentVarData.varTags.push_back(indexedTag);
+        }
+      }
+      else {
+        // This should be a user-created non-set variable without timesteps
+        // Treat it like having one, 0th, timestep
+        std::stringstream ssTagNameWithIndex;
+        ssTagNameWithIndex << varname << 0;
+        rval = mbImpl->tag_get_handle(ssTagNameWithIndex.str().c_str(), indexedTag);
+        ERRORR(rval, "Can't find tag for a user-created variable.");
+        dbgOut.tprintf(2, "    found indexed tag 0 with name %s\n", ssTagNameWithIndex.str().c_str());
+        currentVarData.varTags.push_back(indexedTag);
+      }
+
+      // The type of the tag is fixed though
       DataType type;
-      rval = mbImpl->tag_get_data_type(tag, 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;
-
-      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 {
-        for (size_t j = 0; j < currentVarData.varDims.size(); j++) {
-          currentVarData.writeStarts.push_back(0);
-          currentVarData.writeCounts.push_back(dimLens[currentVarData.varDims[j]]);
-        }
-      }
     }
-  }
+  } // for (size_t i = 0; i < var_names.size(); i++)
 
   // Process coordinate variables here
   // Check that for used coordinates we have found the tags
@@ -231,7 +251,7 @@ ErrorCode NCWriteHelper::collect_variable_data(std::vector<std::string>& var_nam
 
     assert(0 == varCoordData.memoryHogs.size()); // Nothing so far
     varCoordData.memoryHogs.push_back((void*)data);
-  }
+  } // for (std::set<std::string>::iterator setIt ...
 
   return MB_SUCCESS;
 }
@@ -626,44 +646,75 @@ ErrorCode ScdNCWriteHelper::collect_variable_data(std::vector<std::string>& var_
       ERRORR(MB_FAILURE, "Can't find one variable.");
 
     WriteNC::VarData& currentVarData = vit->second;
+    std::vector<int>& varDims = currentVarData.varDims;
+
+    // Skip set variables, which were already processed in NCWriteHelper::collect_variable_data()
+    if (WriteNC::ENTLOCSET == currentVarData.entLoc)
+      continue;
+
+    // Set up writeStarts and writeCounts (maximum number of dimensions is 4)
+    currentVarData.writeStarts.resize(4);
+    currentVarData.writeCounts.resize(4);
+    unsigned int dim_idx = 0;
+
+    // First: time
     if (currentVarData.has_tsteps) {
-      // Support non-set variables with 4 dimensions like (time, lev, lat, lon)
-      assert(4 == currentVarData.varDims.size());
+      // Non-set variables with timesteps
+      // 4 dimensions like (time, lev, lat, lon)
+      // 3 dimensions like (time, lat, lon)
+      assert(4 == varDims.size() || 3 == 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 and lon
-      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, "Unexpected entity location type for structured mesh non-set variable.");
-      }
+      assert(tDim == varDims[0]);
+
+      currentVarData.writeStarts[dim_idx] = 0; // This value is timestep dependent, will be set later
+      currentVarData.writeCounts[dim_idx] = 1;
+      dim_idx++;
+    }
+    else {
+      // Non-set variables without timesteps
+      // 3 dimensions like (lev, lat, lon)
+      // 2 dimensions like (lat, lon)
+      assert(3 == varDims.size() || 2 == varDims.size());
+    }
+
+    // Next: lev
+    if (currentVarData.numLev > 0) {
+      // Non-set variables with levels
+      // 4 dimensions like (time, lev, lat, lon)
+      // 3 dimensions like (lev, lat, lon)
+      assert(4 == varDims.size() || 3 == varDims.size());
+
+      currentVarData.writeStarts[dim_idx] = 0;
+      currentVarData.writeCounts[dim_idx] = currentVarData.numLev;
+      dim_idx++;
+    }
+    else {
+      // Non-set variables without levels
+      // 3 dimensions like (time, lat, lon)
+      // 2 dimensions like (lat, lon)
+      assert(3 == varDims.size() || 2 == varDims.size());
+    }
+
+    // Finally: lat and lon
+    switch (currentVarData.entLoc) {
+      case WriteNC::ENTLOCFACE:
+        // Faces
+        currentVarData.writeStarts[dim_idx] = lCDims[1];
+        currentVarData.writeCounts[dim_idx] = lCDims[4] - lCDims[1] + 1;
+        currentVarData.writeStarts[dim_idx + 1] = lCDims[0];
+        currentVarData.writeCounts[dim_idx + 1] = lCDims[3] - lCDims[0] + 1;
+        break;
+      default:
+        ERRORR(MB_FAILURE, "Unexpected entity location type for structured mesh non-set variable.");
     }
+    dim_idx += 2;
 
     // Get variable size
     currentVarData.sz = 1;
-    for (std::size_t idx = 0; idx != currentVarData.writeCounts.size(); idx++)
+    for (std::size_t idx = 0; idx < dim_idx; idx++)
       currentVarData.sz *= currentVarData.writeCounts[idx];
-  }
+  } // for (size_t i = 0; i < var_names.size(); i++)
 
   return MB_SUCCESS;
 }
@@ -681,9 +732,6 @@ ErrorCode ScdNCWriteHelper::write_nonset_variables(std::vector<WriteNC::VarData>
   for (unsigned int i = 0; i < vdatas.size(); i++) {
     WriteNC::VarData& variableData = vdatas[i];
 
-    // Time should be the first dimension
-    assert(tDim == variableData.varDims[0]);
-
     // Assume this variable is on faces for the time being
     switch (variableData.entLoc) {
       case WriteNC::ENTLOCFACE:
@@ -693,20 +741,51 @@ ErrorCode ScdNCWriteHelper::write_nonset_variables(std::vector<WriteNC::VarData>
         ERRORR(MB_FAILURE, "Unexpected entity location type for structured mesh non-set variable.");
     }
 
-    // A typical CAM-EUL or CAM-FV non-set 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
+    unsigned int num_timesteps;
+    unsigned int lat_idx = 0;
+    unsigned int lon_idx = 1;
+    if (variableData.has_tsteps) {
+      // Non-set variables with timesteps
+      // 4 dimensions like (time, lev, lat, lon)
+      // 3 dimensions like (time, lat, lon)
+      num_timesteps = tstep_nums.size();
+      lat_idx++;
+      lon_idx++;
+    }
+    else {
+      // Non-set variables without timesteps
+      // 3 dimensions like (lev, lat, lon)
+      // 2 dimensions like (lat, lon)
+      num_timesteps = 1;
+    }
+
+    unsigned int num_lev;
+    if (variableData.numLev > 0) {
+      // Non-set variables with levels
+      // 4 dimensions like (time, lev, lat, lon)
+      // 3 dimensions like (lev, lat, lon)
+      num_lev = variableData.numLev;
+      lat_idx++;
+      lon_idx++;
+    }
+    else {
+      // Non-set variables without levels
+      // 3 dimensions like (time, lat, lon)
+      // 2 dimensions like (lat, lon)
+      num_lev = 1;
+    }
 
-    variableData.writeCounts[0] = 1; // We will write one time step
+    size_t ni = variableData.writeCounts[lon_idx]; // lon
+    size_t nj = variableData.writeCounts[lat_idx]; // lat
 
-    for (unsigned int t = 0; t < tstep_nums.size(); t++) {
+    // At each timestep, we need to transpose tag format (lat, lon, lev) back
+    // to NC format (lev, lat, lon) for writing
+    for (unsigned int t = 0; t < num_timesteps; t++) {
       // We will write one time step, and count will be one; start will be different
       // Use tag_iterate to get tag data (assume that localCellsOwned is contiguous)
       // We should also transpose for level so that means deep copy for transpose
-      variableData.writeStarts[0] = t; // This is start for time
+      if (tDim == variableData.varDims[0])
+        variableData.writeStarts[0] = t; // This is start for time
       int count;
       void* dataptr;
       ErrorCode rval = mbImpl->tag_iterate(variableData.varTags[t], localCellsOwned.begin(), localCellsOwned.end(), count, dataptr);
@@ -718,9 +797,10 @@ ErrorCode ScdNCWriteHelper::write_nonset_variables(std::vector<WriteNC::VarData>
       // nonblocking put (request aggregation) later
       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));
+          std::vector<double> tmpdoubledata(ni*nj*num_lev);
+          if (num_lev > 1)
+            // Transpose (lat, lon, lev) back to (lev, lat, lon)
+            jik_to_kji(ni, nj, num_lev, &tmpdoubledata[0], (double*)(dataptr));
           success = NCFUNCAP(_vara_double)(_fileId, variableData.varId,
                     &variableData.writeStarts[0], &variableData.writeCounts[0],
                     &tmpdoubledata[0]);

diff --git a/src/io/NCWriteMPAS.cpp b/src/io/NCWriteMPAS.cpp
index ac0f70a..6a40dad 100644
--- a/src/io/NCWriteMPAS.cpp
+++ b/src/io/NCWriteMPAS.cpp
@@ -157,13 +157,13 @@ ErrorCode NCWriteMPAS::collect_variable_data(std::vector<std::string>& var_names
       ERRORR(MB_FAILURE, "Can't find one variable.");
 
     WriteNC::VarData& currentVarData = vit->second;
+    std::vector<int>& varDims = currentVarData.varDims;
 
     // Skip edge variables, if there are no edges
     if (localEdgesOwned.empty() && currentVarData.entLoc == WriteNC::ENTLOCEDGE)
       continue;
 
     // If nVertLevels dimension is not found, try other optional levels such as nVertLevelsP1
-    std::vector<int>& varDims = currentVarData.varDims;
     if (std::find(varDims.begin(), varDims.end(), levDim) == varDims.end()) {
       for (unsigned int j = 0; j < opt_lev_dims.size(); j++) {
         if (std::find(varDims.begin(), varDims.end(), opt_lev_dims[j]) != varDims.end()) {
@@ -173,60 +173,86 @@ ErrorCode NCWriteMPAS::collect_variable_data(std::vector<std::string>& var_names
       }
     }
 
+    // Skip set variables, which were already processed in NCWriteHelper::collect_variable_data()
+    if (WriteNC::ENTLOCSET == currentVarData.entLoc)
+      continue;
+
+    // Set up writeStarts and writeCounts (maximum number of dimensions is 3)
+    currentVarData.writeStarts.resize(3);
+    currentVarData.writeCounts.resize(3);
+    unsigned int dim_idx = 0;
+
+    // First: Time
     if (currentVarData.has_tsteps) {
-      // Support non-set variables with 3 dimensions like (Time, nCells, nVertLevels), or
+      // Non-set variables with timesteps
+      // 3 dimensions like (Time, nCells, nVertLevels)
       // 2 dimensions like (Time, nCells)
       assert(3 == varDims.size() || 2 == varDims.size());
 
       // Time should be the first dimension
       assert(tDim == 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: nVertices / nCells / nEdges
-      switch (currentVarData.entLoc) {
-        case WriteNC::ENTLOCVERT:
-          // Vertices
-          // Start from the first localGidVerts
-          // Actually, this will be reset later for writing
-          currentVarData.writeStarts[1] = localGidVertsOwned[0] - 1;
-          currentVarData.writeCounts[1] = localGidVertsOwned.size();
-          break;
-        case WriteNC::ENTLOCFACE:
-          // Faces
-          // Start from the first localGidCells
-          // Actually, this will be reset later for writing
-          currentVarData.writeStarts[1] = localGidCellsOwned[0] - 1;
-          currentVarData.writeCounts[1] = localGidCellsOwned.size();
-          break;
-        case WriteNC::ENTLOCEDGE:
-          // Edges
-          // Start from the first localGidEdges
-          // Actually, this will be reset later for writing
-          currentVarData.writeStarts[1] = localGidEdgesOwned[0] - 1;
-          currentVarData.writeCounts[1] = localGidEdgesOwned.size();
-          break;
-        default:
-          ERRORR(MB_FAILURE, "Unexpected entity location type for MPAS non-set variable.");
-      }
+      currentVarData.writeStarts[dim_idx] = 0; // This value is timestep dependent, will be set later
+      currentVarData.writeCounts[dim_idx] = 1;
+      dim_idx++;
+    }
+    else {
+      // Non-set variables without timesteps
+      // 2 dimensions like (nCells, nVertLevels)
+      // 1 dimension like (nCells)
+      assert(2 == varDims.size() || 1 == varDims.size());
+    }
+
+    // Next: nVertices / nCells / nEdges
+    switch (currentVarData.entLoc) {
+      case WriteNC::ENTLOCVERT:
+        // Vertices
+        // Start from the first localGidVerts
+        // Actually, this will be reset later for writing
+        currentVarData.writeStarts[dim_idx] = localGidVertsOwned[0] - 1;
+        currentVarData.writeCounts[dim_idx] = localGidVertsOwned.size();
+        break;
+      case WriteNC::ENTLOCFACE:
+        // Faces
+        // Start from the first localGidCells
+        // Actually, this will be reset later for writing
+        currentVarData.writeStarts[dim_idx] = localGidCellsOwned[0] - 1;
+        currentVarData.writeCounts[dim_idx] = localGidCellsOwned.size();
+        break;
+      case WriteNC::ENTLOCEDGE:
+        // Edges
+        // Start from the first localGidEdges
+        // Actually, this will be reset later for writing
+        currentVarData.writeStarts[dim_idx] = localGidEdgesOwned[0] - 1;
+        currentVarData.writeCounts[dim_idx] = localGidEdgesOwned.size();
+        break;
+      default:
+        ERRORR(MB_FAILURE, "Unexpected entity location type for MPAS non-set variable.");
+    }
+    dim_idx++;
+
+    // Finally: nVertLevels or other optional levels, it is possible that there is no
+    // level dimension (numLev is 0) for this non-set variable, e.g. (Time, nCells)
+    if (currentVarData.numLev > 0) {
+      // Non-set variables with levels
+      // 3 dimensions like (Time, nCells, nVertLevels)
+      // 2 dimensions like (nCells, nVertLevels)
+      assert(3 == varDims.size() || 2 == varDims.size());
 
-      // Finally: nVertLevels or other optional levels, it is possible that there is no
-      // level dimension (numLev is 0) for this non-set variable, e.g. (Time, nCells)
-      if (currentVarData.numLev < 1)
-        currentVarData.numLev = 1;
-      currentVarData.writeStarts[2] = 0;
-      currentVarData.writeCounts[2] = currentVarData.numLev;
+      currentVarData.writeStarts[dim_idx] = 0;
+      currentVarData.writeCounts[dim_idx] = currentVarData.numLev;
+      dim_idx++;
+    }
+    else {
+      // Non-set variables without levels
+      // 2 dimensions like (Time, nCells)
+      // 1 dimension like (nCells)
+      assert(2 == varDims.size() || 1 == varDims.size());
     }
 
     // Get variable size
     currentVarData.sz = 1;
-    for (std::size_t idx = 0; idx != 3; idx++)
+    for (std::size_t idx = 0; idx < dim_idx; idx++)
       currentVarData.sz *= currentVarData.writeCounts[idx];
   }
 
@@ -238,8 +264,6 @@ ErrorCode NCWriteMPAS::write_nonset_variables(std::vector<WriteNC::VarData>& vda
   Interface*& mbImpl = _writeNC->mbImpl;
 
   int success;
-  Range* pLocalEntsOwned = NULL;
-  Range* pLocalGidEntsOwned = NULL;
 
   // For each variable tag in the indexed lists, write a time step data
   for (unsigned int i = 0; i < vdatas.size(); i++) {
@@ -249,10 +273,9 @@ ErrorCode NCWriteMPAS::write_nonset_variables(std::vector<WriteNC::VarData>& vda
     if (localEdgesOwned.empty() && variableData.entLoc == WriteNC::ENTLOCEDGE)
       continue;
 
-    // Time should be the first dimension
-    assert(tDim == variableData.varDims[0]);
-
     // Get local owned entities of this variable
+    Range* pLocalEntsOwned = NULL;
+    Range* pLocalGidEntsOwned = NULL;
     switch (variableData.entLoc) {
       case WriteNC::ENTLOCVERT:
         // Vertices
@@ -273,15 +296,45 @@ ErrorCode NCWriteMPAS::write_nonset_variables(std::vector<WriteNC::VarData>& vda
         ERRORR(MB_FAILURE, "Unexpected entity location type for MPAS non-set variable.");
     }
 
-    // A typical MPAS non-set variable has 3 dimensions as (Time, nCells, nVertLevels)
-    for (unsigned int t = 0; t < tstep_nums.size(); t++) {
+    unsigned int num_timesteps;
+    unsigned int ents_idx = 0;
+    if (variableData.has_tsteps) {
+      // Non-set variables with timesteps
+      // 3 dimensions like (Time, nCells, nVertLevels)
+      // 2 dimensions like (Time, nCells)
+      num_timesteps = tstep_nums.size();
+      ents_idx++;
+    }
+    else {
+      // Non-set variables without timesteps
+      // 2 dimensions like (nCells, nVertLevels)
+      // 1 dimension like (nCells)
+      num_timesteps = 1;
+    }
+
+    unsigned int num_lev;
+    if (variableData.numLev > 0) {
+      // Non-set variables with levels
+      // 3 dimensions like (Time, nCells, nVertLevels)
+      // 2 dimensions like (nCells, nVertLevels)
+      num_lev = variableData.numLev;
+    }
+    else {
+      // Non-set variables without levels
+      // 2 dimensions like (Time, nCells)
+      // 1 dimension like (nCells)
+      num_lev = 1;
+    }
+
+    for (unsigned int t = 0; t < num_timesteps; t++) {
       // We will write one time step, and count will be one; start will be different
       // Use tag_get_data instead of tag_iterate to copy tag data, as localEntsOwned
       // might not be contiguous.
-      variableData.writeStarts[0] = t; // This is start for time
-      std::vector<double> tag_data(pLocalEntsOwned->size() * variableData.numLev);
+      if (tDim == variableData.varDims[0])
+        variableData.writeStarts[0] = t; // This is start for time
+      std::vector<double> tag_data(pLocalEntsOwned->size() * num_lev);
       ErrorCode rval = mbImpl->tag_get_data(variableData.varTags[t], *pLocalEntsOwned, &tag_data[0]);
-      ERRORR(rval, "Trouble getting tag data on owned vertices.");
+      ERRORR(rval, "Trouble getting tag data on owned entities.");
 
 #ifdef PNETCDF_FILE
       size_t nb_writes = pLocalGidEntsOwned->psize();
@@ -299,8 +352,8 @@ ErrorCode NCWriteMPAS::write_nonset_variables(std::vector<WriteNC::VarData>& vda
               pair_iter != pLocalGidEntsOwned->pair_end(); ++pair_iter, ic++) {
             EntityHandle starth = pair_iter->first;
             EntityHandle endh = pair_iter->second;
-            variableData.writeStarts[1] = (NCDF_SIZE)(starth - 1);
-            variableData.writeCounts[1] = (NCDF_SIZE)(endh - starth + 1);
+            variableData.writeStarts[ents_idx] = (NCDF_SIZE)(starth - 1);
+            variableData.writeCounts[ents_idx] = (NCDF_SIZE)(endh - starth + 1);
 
             // Do a partial write, in each subrange
 #ifdef PNETCDF_FILE
@@ -316,7 +369,7 @@ ErrorCode NCWriteMPAS::write_nonset_variables(std::vector<WriteNC::VarData>& vda
             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;
+            indexInDoubleArray += (endh - starth + 1) * num_lev;
           }
           assert(ic == pLocalGidEntsOwned->psize());
 #ifdef PNETCDF_FILE

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