Mesh Oriented datABase  (version 5.5.1)
An array-based unstructured mesh library
NCWriteHelper.cpp
Go to the documentation of this file.
1 /*
2  * NCWriteHelper.cpp
3  *
4  * Created on: Mar 28, 2014
5  * Author: iulian
6  */
7 
8 #include "NCWriteHelper.hpp"
9 #include "NCWriteEuler.hpp"
10 #include "NCWriteFV.hpp"
11 #include "NCWriteHOMME.hpp"
12 #include "NCWriteMPAS.hpp"
13 #include "NCWriteGCRM.hpp"
14 
15 #include "moab/WriteUtilIface.hpp"
16 #include "MBTagConventions.hpp"
17 
18 #include <sstream>
19 
20 #ifdef WIN32
21 #ifdef size_t
22 #undef size_t
23 #endif
24 #endif
25 
26 namespace moab
27 {
28 
29 //! Get appropriate helper instance for WriteNC class; based on some info in the file set
31  int fileId,
32  const FileOptions& opts,
33  EntityHandle fileSet )
34 {
35  std::string& grid_type = writeNC->grid_type;
36  if( grid_type == "CAM_EUL" )
37  return new( std::nothrow ) NCWriteEuler( writeNC, fileId, opts, fileSet );
38  else if( grid_type == "CAM_FV" )
39  return new( std::nothrow ) NCWriteFV( writeNC, fileId, opts, fileSet );
40  else if( grid_type == "CAM_SE" )
41  return new( std::nothrow ) NCWriteHOMME( writeNC, fileId, opts, fileSet );
42  else if( grid_type == "MPAS" )
43  return new( std::nothrow ) NCWriteMPAS( writeNC, fileId, opts, fileSet );
44  else if( grid_type == "GCRM" )
45  return new( std::nothrow ) NCWriteGCRM( writeNC, fileId, opts, fileSet );
46 
47  // Unknown NetCDF grid
48  return NULL;
49 }
50 
51 ErrorCode NCWriteHelper::collect_variable_data( std::vector< std::string >& var_names, std::vector< int >& tstep_nums )
52 {
53  Interface*& mbImpl = _writeNC->mbImpl;
54  std::vector< std::string >& dimNames = _writeNC->dimNames;
55  std::vector< int >& dimLens = _writeNC->dimLens;
56  std::set< std::string >& usedCoordinates = _writeNC->usedCoordinates;
57  std::set< std::string >& dummyVarNames = _writeNC->dummyVarNames;
58  std::map< std::string, WriteNC::VarData >& varInfo = _writeNC->varInfo;
59  DebugOutput& dbgOut = _writeNC->dbgOut;
60 
61  ErrorCode rval;
62 
63  usedCoordinates.clear();
64 
65  if( tstep_nums.empty() && nTimeSteps > 0 )
66  {
67  // No timesteps input, get them all
68  for( int i = 0; i < nTimeSteps; i++ )
69  tstep_nums.push_back( i );
70  }
71 
72  for( size_t i = 0; i < var_names.size(); i++ )
73  {
74  std::string varname = var_names[i];
75  std::map< std::string, WriteNC::VarData >::iterator vit = varInfo.find( varname );
76  if( vit == varInfo.end() ) MB_SET_ERR( MB_FAILURE, "Can't find variable " << varname );
77 
78  WriteNC::VarData& currentVarData = vit->second;
79 
80  dbgOut.tprintf( 2, " for variable %s varDims.size %d \n", varname.c_str(),
81  (int)currentVarData.varDims.size() );
82  for( size_t j = 0; j < currentVarData.varDims.size(); j++ )
83  {
84  std::string dimName = dimNames[currentVarData.varDims[j]];
85  vit = varInfo.find( dimName );
86  if( vit == varInfo.end() ) MB_SET_ERR( MB_FAILURE, "Can't find coordinate variable " << dimName );
87 
88  usedCoordinates.insert( dimName ); // Collect those used, we will need to write them to the file
89  dbgOut.tprintf( 2, " for variable %s need dimension %s with length %d\n", varname.c_str(),
90  dimName.c_str(), dimLens[currentVarData.varDims[j]] );
91  }
92 
93  // Process coordinate variables later
94  if( usedCoordinates.find( varname ) != usedCoordinates.end() ) continue;
95 
96  // Default has_tsteps is false
97  if( std::find( currentVarData.varDims.begin(), currentVarData.varDims.end(), tDim ) !=
98  currentVarData.varDims.end() )
99  currentVarData.has_tsteps = true;
100 
101  // Default numLev is 0
102  if( ( std::find( currentVarData.varDims.begin(), currentVarData.varDims.end(), levDim ) !=
103  currentVarData.varDims.end() ) )
104  currentVarData.numLev = nLevels;
105 
106  // Process set variables
107  if( WriteNC::ENTLOCSET == currentVarData.entLoc )
108  {
109  if( currentVarData.has_tsteps )
110  {
111  // Set variables with timesteps, e.g. xtime(Time) or xtime(Time, StrLen)
112  // TBD
113  MB_SET_ERR( MB_NOT_IMPLEMENTED, "Writing set variables with timesteps is not implemented yet" );
114  }
115  else
116  {
117  // Get the tag with varname
118  Tag tag = 0;
119  rval = mbImpl->tag_get_handle( varname.c_str(), tag );MB_CHK_SET_ERR( rval, "Can't find tag " << varname );
120  currentVarData.varTags.push_back( tag ); // Really, only one for these
121  const void* data;
122  int size;
123  rval = mbImpl->tag_get_by_ptr( tag, &_fileSet, 1, &data, &size );MB_CHK_SET_ERR( rval, "Can't get data of tag " << varname );
124 
125  // Find the type of tag, and use it
126  DataType type;
127  rval = mbImpl->tag_get_data_type( tag, type );MB_CHK_SET_ERR( rval, "Can't get data type of tag " << varname );
128 
129  currentVarData.varDataType = NC_DOUBLE;
130  if( MB_TYPE_INTEGER == type ) currentVarData.varDataType = NC_INT;
131 
132  assert( 0 == currentVarData.memoryHogs.size() ); // Nothing so far
133  currentVarData.memoryHogs.push_back( (void*)data );
134 
135  if( currentVarData.varDims.empty() )
136  {
137  // Scalar variable
138  currentVarData.writeStarts.push_back( 0 );
139  currentVarData.writeCounts.push_back( 1 );
140  }
141  else
142  {
143  for( size_t j = 0; j < currentVarData.varDims.size(); j++ )
144  {
145  currentVarData.writeStarts.push_back( 0 );
146  currentVarData.writeCounts.push_back( dimLens[currentVarData.varDims[j]] );
147  }
148  }
149 
150  // Get variable size
151  currentVarData.sz = 1;
152  for( std::size_t idx = 0; idx != currentVarData.writeCounts.size(); idx++ )
153  currentVarData.sz *= currentVarData.writeCounts[idx];
154  }
155  } // if (WriteNC::ENTLOCSET == currentVarData.entLoc)
156  // Process non-set variables
157  else
158  {
159  Tag indexedTag = 0;
160 
161  if( currentVarData.has_tsteps )
162  {
163  for( unsigned int t = 0; t < tstep_nums.size(); t++ )
164  {
165  std::stringstream ssTagNameWithIndex;
166  ssTagNameWithIndex << varname << tstep_nums[t];
167  rval = mbImpl->tag_get_handle( ssTagNameWithIndex.str().c_str(), indexedTag );MB_CHK_SET_ERR( rval, "Can't find tag " << ssTagNameWithIndex.str() );
168  dbgOut.tprintf( 2, " found indexed tag %d with name %s\n", tstep_nums[t],
169  ssTagNameWithIndex.str().c_str() );
170  currentVarData.varTags.push_back( indexedTag );
171  }
172  }
173  else
174  {
175  // This should be a user-created non-set variable without timesteps
176  // Treat it like having one, 0th, timestep
177  std::stringstream ssTagNameWithIndex;
178  ssTagNameWithIndex << varname << 0;
179  rval = mbImpl->tag_get_handle( ssTagNameWithIndex.str().c_str(), indexedTag );MB_CHK_SET_ERR( rval, "Can't find tag " << ssTagNameWithIndex.str() << " for a user-created variable" );
180  dbgOut.tprintf( 2, " found indexed tag 0 with name %s\n", ssTagNameWithIndex.str().c_str() );
181  currentVarData.varTags.push_back( indexedTag );
182  }
183 
184  // The type of the tag is fixed though
185  DataType type;
186  rval = mbImpl->tag_get_data_type( indexedTag, type );MB_CHK_SET_ERR( rval, "Can't get data type of tag " << varname );
187 
188  currentVarData.varDataType = NC_DOUBLE;
189  if( MB_TYPE_INTEGER == type ) currentVarData.varDataType = NC_INT;
190  }
191  } // for (size_t i = 0; i < var_names.size(); i++)
192 
193  // Process coordinate variables here
194  // Check that for used coordinates we have found the tags
195  for( std::set< std::string >::iterator setIt = usedCoordinates.begin(); setIt != usedCoordinates.end(); ++setIt )
196  {
197  const std::string& coordName = *setIt;
198 
199  std::map< std::string, WriteNC::VarData >::iterator vit = varInfo.find( coordName );
200  if( vit == varInfo.end() ) MB_SET_ERR( MB_FAILURE, "Can't find coordinate variable " << coordName );
201 
202  WriteNC::VarData& varCoordData = vit->second;
203  Tag coordTag = 0;
204  rval = mbImpl->tag_get_handle( coordName.c_str(), coordTag );MB_CHK_SET_ERR( rval, "Can't find tag " << coordName );
205  varCoordData.varTags.push_back( coordTag ); // Really, only one for these
206 
207  const void* data;
208  int sizeCoordinate;
209  rval = mbImpl->tag_get_by_ptr( coordTag, &_fileSet, 1, &data, &sizeCoordinate );MB_CHK_SET_ERR( rval, "Can't get coordinate values of " << coordName );
210  dbgOut.tprintf( 2, " found coordinate tag with name %s and length %d\n", coordName.c_str(), sizeCoordinate );
211 
212  // Find the type of tag, and use it
213  DataType type;
214  rval = mbImpl->tag_get_data_type( coordTag, type );MB_CHK_SET_ERR( rval, "Can't get data type of tag " << coordName );
215  varCoordData.varDataType = NC_DOUBLE;
216  if( MB_TYPE_INTEGER == type ) varCoordData.varDataType = NC_INT;
217 
218  // Get dimension length (the only dimension of this coordinate variable, with the same name)
219  assert( 1 == varCoordData.varDims.size() );
220  int coordDimLen = dimLens[varCoordData.varDims[0]];
221 
222  if( dummyVarNames.find( coordName ) != dummyVarNames.end() )
223  {
224  // For a dummy coordinate variable, the tag size is always 1
225  // The number of coordinates should be set to dimension length, instead of 1
226  assert( 1 == sizeCoordinate );
227  sizeCoordinate = coordDimLen;
228 
229  // No variable data to write
230  data = NULL;
231  }
232  else
233  {
234  // The number of coordinates should be exactly the same as dimension length
235  // However, if timesteps spread across files and time tag has been updated,
236  // sizeCoordinate will be larger
237  if( varCoordData.varDims[0] != tDim ) assert( sizeCoordinate == coordDimLen );
238  }
239 
240  // For time, the actual output size and values are determined by tstep_nums
241  if( varCoordData.varDims[0] == tDim )
242  {
243  // Does not apply to dummy time tag (e.g. 'Time' tag of MPAS), when timesteps
244  // spread across files
245  if( NULL != data ) assert( tstep_nums.size() > 0 && tstep_nums.size() <= (size_t)sizeCoordinate );
246 
247  sizeCoordinate = tstep_nums.size();
248 
249  if( NULL != data )
250  {
251  assert( NC_DOUBLE == varCoordData.varDataType );
252  timeStepVals.resize( sizeCoordinate );
253  for( unsigned int t = 0; t < tstep_nums.size(); t++ )
254  timeStepVals[t] = ( (double*)data )[tstep_nums[t]];
255 
256  data = &timeStepVals[0];
257  }
258  }
259 
260  // This is the length
261  varCoordData.sz = sizeCoordinate;
262  varCoordData.writeStarts.resize( 1 );
263  varCoordData.writeStarts[0] = 0;
264  varCoordData.writeCounts.resize( 1 );
265  varCoordData.writeCounts[0] = sizeCoordinate;
266 
267  assert( 0 == varCoordData.memoryHogs.size() ); // Nothing so far
268  varCoordData.memoryHogs.push_back( (void*)data );
269  } // for (std::set<std::string>::iterator setIt ...
270 
271  return MB_SUCCESS;
272 }
273 
274 ErrorCode NCWriteHelper::init_file( std::vector< std::string >& var_names,
275  std::vector< std::string >& desired_names,
276  bool append )
277 {
278  std::vector< std::string >& dimNames = _writeNC->dimNames;
279  std::set< std::string >& usedCoordinates = _writeNC->usedCoordinates;
280  std::set< std::string >& dummyVarNames = _writeNC->dummyVarNames;
281  std::map< std::string, WriteNC::VarData >& varInfo = _writeNC->varInfo;
282  std::map< std::string, WriteNC::AttData >& globalAtts = _writeNC->globalAtts;
283  DebugOutput& dbgOut = _writeNC->dbgOut;
284 
285  int tDim_in_dimNames = tDim;
286  int levDim_in_dimNames = levDim;
287 
288  // If append mode, make sure we are in define mode; a simple open will not allow creation of new
289  // variables
290  if( append )
291  {
292  int errcode = NCFUNC( redef )( _fileId );
293  if( errcode != NC_NOERR ) MB_SET_ERR( MB_FAILURE, "Can't open file in redefine mode" );
294  }
295 
296  // First initialize all coordinates, then fill VarData for actual variables (and dimensions)
297  // Check that for used coordinates we have found the tags
298  for( std::set< std::string >::iterator setIt = usedCoordinates.begin(); setIt != usedCoordinates.end(); ++setIt )
299  {
300  const std::string& coordName = *setIt;
301 
302  std::map< std::string, WriteNC::VarData >::iterator vit = varInfo.find( coordName );
303  if( vit == varInfo.end() ) MB_SET_ERR( MB_FAILURE, "Can't find coordinate variable " << coordName );
304 
305  WriteNC::VarData& varCoordData = vit->second;
306  varCoordData.varDims.resize( 1 );
307 
308  // If not append, create it for sure
309  // If append, we might already have it, including the tag / variable with the same name
310  /*
311  * int ncmpi_inq_dimid(int ncid, const char *name, int *idp);
312  */
313  if( append )
314  {
315  int dimId;
316  if( NCFUNC( inq_dimid )( _fileId, coordName.c_str(), &dimId ) == NC_NOERR )
317  { // If not found, create it later
318  varCoordData.varDims[0] = dimId;
319  dbgOut.tprintf( 2, " file already has coordName %s dim id is %d \n", coordName.c_str(),
320  (int)varCoordData.varDims[0] );
321 
322  // Update tDim and levDim to actual dimension id
323  if( coordName == dimNames[tDim_in_dimNames] )
324  tDim = varCoordData.varDims[0];
325  else if( coordName == dimNames[levDim_in_dimNames] )
326  levDim = varCoordData.varDims[0];
327 
328  // Skip dummy coordinate variables (e.g. ncol)
329  if( dummyVarNames.find( coordName ) != dummyVarNames.end() ) continue;
330 
331  // Check that the coordinate is a variable too
332  // Inquire for a variable with the same name
333  int varId;
334  if( NCFUNC( inq_varid )( _fileId, coordName.c_str(), &varId ) != NC_NOERR )
335  MB_SET_ERR( MB_FAILURE, "We do not have a variable with the same name " << coordName );
336  // We should also check that this variable has one dimension, and it is dimId
337  varCoordData.varId = varId;
338  dbgOut.tprintf( 2, " file already has coordinate %s and varId is %d \n", coordName.c_str(), varId );
339 
340  continue; // Maybe more checks are needed here
341  }
342  }
343 
344  /* int nc_def_dim (int ncid, const char *name, size_t len, int *dimidp);
345  * example: status = nc_def_dim(fileId, "lat", 18L, &latid);
346  */
347 
348  // Actually define a dimension
349  if( NCFUNC( def_dim )( _fileId, coordName.c_str(), (size_t)varCoordData.sz, &varCoordData.varDims[0] ) !=
350  NC_NOERR )
351  MB_SET_ERR( MB_FAILURE, "Failed to generate dimension " << coordName );
352  dbgOut.tprintf( 2, " for coordName %s dim id is %d \n", coordName.c_str(), (int)varCoordData.varDims[0] );
353 
354  // Update tDim and levDim to actual dimension id
355  if( coordName == dimNames[tDim_in_dimNames] )
356  tDim = varCoordData.varDims[0];
357  else if( coordName == dimNames[levDim_in_dimNames] )
358  levDim = varCoordData.varDims[0];
359 
360  // Create a variable with the same name, and its only dimension the one we just defined
361  /*
362  * int nc_def_var (int ncid, const char *name, nc_type xtype,
363  int ndims, const int dimids[], int *varidp);
364  example:
365  http://www.unidata.ucar.edu/software/netcdf/docs/netcdf-c/nc_005fdef_005fvar.html#nc_005fdef_005fvar
366  */
367 
368  // Skip dummy coordinate variables (e.g. ncol)
369  if( dummyVarNames.find( coordName ) != dummyVarNames.end() ) continue;
370 
371  // Define a coordinate variable
372  if( NCFUNC( def_var )( _fileId, coordName.c_str(), varCoordData.varDataType, 1, &( varCoordData.varDims[0] ),
373  &varCoordData.varId ) != NC_NOERR )
374  MB_SET_ERR( MB_FAILURE, "Failed to create coordinate variable " << coordName );
375 
376  dbgOut.tprintf( 2, " for coordName %s variable id is %d \n", coordName.c_str(), varCoordData.varId );
377  }
378 
379  // Now look at requested variables, and update from the index in dimNames to the actual
380  // dimension id
381  for( size_t i = 0; i < var_names.size(); i++ )
382  {
383  std::map< std::string, WriteNC::VarData >::iterator vit = varInfo.find( var_names[i] );
384  if( vit == varInfo.end() ) MB_SET_ERR( MB_FAILURE, "Can't find requested variable " << var_names[i] );
385 
386  // Skip coordinate variables
387  if( usedCoordinates.find( var_names[i] ) != usedCoordinates.end() ) continue;
388 
389  WriteNC::VarData& variableData = vit->second;
390 
391  // The index is for dimNames; we need to find out the actual dimension id (from above)
392  int numDims = (int)variableData.varDims.size();
393  for( int j = 0; j < numDims; j++ )
394  {
395  std::string dimName = dimNames[variableData.varDims[j]];
396  std::map< std::string, WriteNC::VarData >::iterator vit2 = varInfo.find( dimName );
397  if( vit2 == varInfo.end() )
398  MB_SET_ERR( MB_FAILURE, "Can't find requested coordinate variable " << dimName );
399 
400  WriteNC::VarData& coordData = vit2->second;
401  // Index in dimNames to actual dimension id
402  variableData.varDims[j] = coordData.varDims[0]; // This one, being a coordinate, is the only one
403  dbgOut.tprintf( 2, " dimension with index %d name %s has ID %d \n", j, dimName.c_str(),
404  variableData.varDims[j] );
405  }
406 
407  // Define the variable now:
408  int errCode =
409  NCFUNC( def_var )( _fileId, desired_names[i].c_str(), variableData.varDataType,
410  (int)variableData.varDims.size(), &( variableData.varDims[0] ), &variableData.varId );
411  if( errCode != NC_NOERR ) MB_SET_ERR( MB_FAILURE, "Failed to create requested variable " << desired_names[i] );
412 
413  dbgOut.tprintf( 2, " for variable %s with desired name %s variable id is %d \n", var_names[i].c_str(),
414  desired_names[i].c_str(), variableData.varId );
415  // Now define the variable, with all dimensions
416  }
417 
418  // Define global attributes (exactly copied from the original file for the time being)
419  // Should we modify some of them (e.g. revision_Id) later?
420  std::map< std::string, WriteNC::AttData >::iterator attIt;
421  for( attIt = globalAtts.begin(); attIt != globalAtts.end(); ++attIt )
422  {
423  const std::string& attName = attIt->first;
424  WriteNC::AttData& attData = attIt->second;
425  NCDF_SIZE& attLen = attData.attLen;
426  nc_type& attDataType = attData.attDataType;
427  const std::string& attValue = attData.attValue;
428 
429  switch( attDataType )
430  {
431  case NC_BYTE:
432  case NC_CHAR:
433  if( NC_NOERR !=
434  NCFUNC( put_att_text )( _fileId, NC_GLOBAL, attName.c_str(), attLen, attValue.c_str() ) )
435  MB_SET_ERR( MB_FAILURE, "Failed to define text type attribute" );
436  break;
437  case NC_DOUBLE:
438  if( NC_NOERR != NCFUNC( put_att_double )( _fileId, NC_GLOBAL, attName.c_str(), NC_DOUBLE, 1,
439  (double*)attValue.c_str() ) )
440  MB_SET_ERR( MB_FAILURE, "Failed to define double type attribute" );
441  break;
442  case NC_FLOAT:
443  if( NC_NOERR != NCFUNC( put_att_float )( _fileId, NC_GLOBAL, attName.c_str(), NC_FLOAT, 1,
444  (float*)attValue.c_str() ) )
445  MB_SET_ERR( MB_FAILURE, "Failed to define float type attribute" );
446  break;
447  case NC_INT:
448  if( NC_NOERR !=
449  NCFUNC( put_att_int )( _fileId, NC_GLOBAL, attName.c_str(), NC_INT, 1, (int*)attValue.c_str() ) )
450  MB_SET_ERR( MB_FAILURE, "Failed to define int type attribute" );
451  break;
452  case NC_SHORT:
453  if( NC_NOERR != NCFUNC( put_att_short )( _fileId, NC_GLOBAL, attName.c_str(), NC_SHORT, 1,
454  (short*)attValue.c_str() ) )
455  MB_SET_ERR( MB_FAILURE, "Failed to define short type attribute" );
456  break;
457  default:
458  MB_SET_ERR( MB_FAILURE, "Unknown attribute data type" );
459  }
460  }
461 
462  // Take it out of define mode
463  if( NC_NOERR != NCFUNC( enddef )( _fileId ) ) MB_SET_ERR( MB_FAILURE, "Failed to close define mode" );
464 
465  return MB_SUCCESS;
466 }
467 
468 ErrorCode NCWriteHelper::write_values( std::vector< std::string >& var_names, std::vector< int >& tstep_nums )
469 {
470  std::set< std::string >& usedCoordinates = _writeNC->usedCoordinates;
471  std::set< std::string >& dummyVarNames = _writeNC->dummyVarNames;
472  std::map< std::string, WriteNC::VarData >& varInfo = _writeNC->varInfo;
473 
474  std::vector< WriteNC::VarData > vdatas;
475  std::vector< WriteNC::VarData > vsetdatas;
476 
477  // For set variables, include coordinates used by requested var_names
478  for( std::set< std::string >::iterator setIt = usedCoordinates.begin(); setIt != usedCoordinates.end(); ++setIt )
479  {
480  const std::string& coordName = *setIt;
481 
482  // Skip dummy coordinate variables (if any)
483  if( dummyVarNames.find( coordName ) != dummyVarNames.end() ) continue;
484 
485  std::map< std::string, WriteNC::VarData >::iterator vit = varInfo.find( coordName );
486  if( vit == varInfo.end() )
487  {
488  MB_SET_ERR( MB_FAILURE, "Can't find coordinate variable " << coordName );
489  }
490 
491  vsetdatas.push_back( vit->second );
492  }
493 
494  // Collect non-set and set variables from requested var_names
495  for( unsigned int i = 0; i < var_names.size(); i++ )
496  {
497  std::map< std::string, WriteNC::VarData >::iterator vit = varInfo.find( var_names[i] );
498  if( vit == varInfo.end() )
499  {
500  MB_SET_ERR( MB_FAILURE, "Can't find requested variable " << var_names[i] );
501  }
502 
503  WriteNC::VarData& variableData = vit->second;
504  if( WriteNC::ENTLOCSET == variableData.entLoc )
505  {
506  // Used coordinates has all ready been included
507  if( usedCoordinates.find( var_names[i] ) != usedCoordinates.end() ) continue;
508 
509  vsetdatas.push_back( variableData );
510  }
511  else
512  vdatas.push_back( variableData );
513  }
514 
515  // Assume that the data ranges do not overlap across processors
516  // While overlapped writing might still work, we should better not take that risk
517  write_nonset_variables( vdatas, tstep_nums );
518 
519  // Use independent I/O mode put, since this write is only for the root processor
520  write_set_variables( vsetdatas, tstep_nums );
521 
522  return MB_SUCCESS;
523 }
524 
525 ErrorCode NCWriteHelper::write_set_variables( std::vector< WriteNC::VarData >& vsetdatas,
526  std::vector< int >& /* tstep_nums */ )
527 {
528  int success;
529 
530  // CAUTION: if the NetCDF ID is from a previous call to ncmpi_create rather than ncmpi_open,
531  // all processors need to call ncmpi_begin_indep_data(). If only the root processor does so,
532  // ncmpi_begin_indep_data() call will be blocked forever :(
533 #ifdef MOAB_HAVE_PNETCDF
534  // Enter independent I/O mode
535  success = NCFUNC( begin_indep_data )( _fileId );
536  if( success ) MB_SET_ERR( MB_FAILURE, "Failed to begin independent I/O mode" );
537 #endif
538 
539  int rank = 0;
540 #ifdef MOAB_HAVE_MPI
541  bool& isParallel = _writeNC->isParallel;
542  if( isParallel )
543  {
544  ParallelComm*& myPcomm = _writeNC->myPcomm;
545  rank = myPcomm->proc_config().proc_rank();
546  }
547 #endif
548  if( 0 == rank )
549  {
550  for( unsigned int i = 0; i < vsetdatas.size(); i++ )
551  {
552  WriteNC::VarData& variableData = vsetdatas[i];
553 
554  // Set variables with timesteps, e.g. xtime(Time) or xtime(Time, StrLen)
555  if( variableData.has_tsteps )
556  {
557  MB_SET_ERR( MB_NOT_IMPLEMENTED, "Writing set variables with timesteps is not implemented yet" );
558  }
559 
560  switch( variableData.varDataType )
561  {
562  case NC_DOUBLE:
563  // Independent I/O mode put
564  success = NCFUNCP( _vara_double )( _fileId, variableData.varId, &variableData.writeStarts[0],
565  &variableData.writeCounts[0],
566  (double*)( variableData.memoryHogs[0] ) );
567  if( success )
568  MB_SET_ERR( MB_FAILURE, "Failed to write double data for variable " << variableData.varName );
569  break;
570  case NC_INT:
571  // Independent I/O mode put
572  success =
573  NCFUNCP( _vara_int )( _fileId, variableData.varId, &variableData.writeStarts[0],
574  &variableData.writeCounts[0], (int*)( variableData.memoryHogs[0] ) );
575  if( success )
576  MB_SET_ERR( MB_FAILURE, "Failed to write int data for variable " << variableData.varName );
577  break;
578  default:
579  MB_SET_ERR( MB_NOT_IMPLEMENTED, "Writing non-double or non-int data is not implemented yet" );
580  }
581  }
582  }
583 
584 #ifdef MOAB_HAVE_PNETCDF
585  // End independent I/O mode
586  success = NCFUNC( end_indep_data )( _fileId );
587  if( success ) MB_SET_ERR( MB_FAILURE, "Failed to end independent I/O mode" );
588 #endif
589 
590  return MB_SUCCESS;
591 }
592 
594 {
595  Interface*& mbImpl = _writeNC->mbImpl;
596  std::vector< std::string >& dimNames = _writeNC->dimNames;
597  std::vector< int >& dimLens = _writeNC->dimLens;
598 
599  ErrorCode rval;
600 
601  // Look for time dimension
602  std::vector< std::string >::iterator vecIt;
603  if( ( vecIt = std::find( dimNames.begin(), dimNames.end(), "time" ) ) != dimNames.end() )
604  tDim = vecIt - dimNames.begin();
605  else if( ( vecIt = std::find( dimNames.begin(), dimNames.end(), "t" ) ) != dimNames.end() )
606  tDim = vecIt - dimNames.begin();
607  else
608  {
609  MB_SET_ERR( MB_FAILURE, "Couldn't find 'time' or 't' dimension" );
610  }
611  nTimeSteps = dimLens[tDim];
612 
613  // Get number of levels
614  if( ( vecIt = std::find( dimNames.begin(), dimNames.end(), "lev" ) ) != dimNames.end() )
615  levDim = vecIt - dimNames.begin();
616  else if( ( vecIt = std::find( dimNames.begin(), dimNames.end(), "ilev" ) ) != dimNames.end() )
617  levDim = vecIt - dimNames.begin();
618  else
619  {
620  MB_SET_ERR( MB_FAILURE, "Couldn't find 'lev' or 'ilev' dimension" );
621  }
622  nLevels = dimLens[levDim];
623 
624  // __<dim_name>_LOC_MINMAX (for slon, slat, lon and lat)
625  Tag convTag = 0;
626  rval = mbImpl->tag_get_handle( "__slon_LOC_MINMAX", 0, MB_TYPE_INTEGER, convTag, MB_TAG_ANY );MB_CHK_SET_ERR( rval, "Trouble getting conventional tag __slon_LOC_MINMAX" );
627  int val[2];
628  rval = mbImpl->tag_get_data( convTag, &_fileSet, 1, val );MB_CHK_SET_ERR( rval, "Trouble getting data of conventional tag __slon_LOC_MINMAX" );
629  lDims[0] = val[0];
630  lDims[3] = val[1];
631 
632  rval = mbImpl->tag_get_handle( "__slat_LOC_MINMAX", 0, MB_TYPE_INTEGER, convTag, MB_TAG_ANY );MB_CHK_SET_ERR( rval, "Trouble getting conventional tag __slat_LOC_MINMAX" );
633  rval = mbImpl->tag_get_data( convTag, &_fileSet, 1, val );MB_CHK_SET_ERR( rval, "Trouble getting data of conventional tag __slat_LOC_MINMAX" );
634  lDims[1] = val[0];
635  lDims[4] = val[1];
636 
637  rval = mbImpl->tag_get_handle( "__lon_LOC_MINMAX", 0, MB_TYPE_INTEGER, convTag, MB_TAG_ANY );MB_CHK_SET_ERR( rval, "Trouble getting conventional tag __lon_LOC_MINMAX" );
638  rval = mbImpl->tag_get_data( convTag, &_fileSet, 1, val );MB_CHK_SET_ERR( rval, "Trouble getting data of conventional tag __lon_LOC_MINMAX" );
639  lCDims[0] = val[0];
640  lCDims[3] = val[1];
641 
642  rval = mbImpl->tag_get_handle( "__lat_LOC_MINMAX", 0, MB_TYPE_INTEGER, convTag, MB_TAG_ANY );MB_CHK_SET_ERR( rval, "Trouble getting conventional tag __lat_LOC_MINMAX" );
643  rval = mbImpl->tag_get_data( convTag, &_fileSet, 1, val );MB_CHK_SET_ERR( rval, "Trouble getting data of conventional tag __lat_LOC_MINMAX" );
644  lCDims[1] = val[0];
645  lCDims[4] = val[1];
646 
647  // Get local faces
648  rval = mbImpl->get_entities_by_dimension( _fileSet, 2, localCellsOwned );MB_CHK_SET_ERR( rval, "Trouble getting local faces in current file set" );
649  assert( !localCellsOwned.empty() );
650 
651 #ifdef MOAB_HAVE_MPI
652  bool& isParallel = _writeNC->isParallel;
653  if( isParallel )
654  {
655  ParallelComm*& myPcomm = _writeNC->myPcomm;
656  int procs = myPcomm->proc_config().proc_size();
657  if( procs > 1 )
658  {
659  rval = myPcomm->filter_pstatus( localCellsOwned, PSTATUS_NOT_OWNED, PSTATUS_NOT );MB_CHK_SET_ERR( rval, "Trouble getting owned faces in current file set" );
660  }
661  }
662 #endif
663 
664  return MB_SUCCESS;
665 }
666 
667 ErrorCode ScdNCWriteHelper::collect_variable_data( std::vector< std::string >& var_names,
668  std::vector< int >& tstep_nums )
669 {
670  NCWriteHelper::collect_variable_data( var_names, tstep_nums );
671 
672  std::map< std::string, WriteNC::VarData >& varInfo = _writeNC->varInfo;
673 
674  for( size_t i = 0; i < var_names.size(); i++ )
675  {
676  std::string varname = var_names[i];
677  std::map< std::string, WriteNC::VarData >::iterator vit = varInfo.find( varname );
678  if( vit == varInfo.end() ) MB_SET_ERR( MB_FAILURE, "Can't find variable " << varname );
679 
680  WriteNC::VarData& currentVarData = vit->second;
681 #ifndef NDEBUG
682  std::vector< int >& varDims = currentVarData.varDims;
683 #endif
684 
685  // Skip set variables, which were already processed in
686  // NCWriteHelper::collect_variable_data()
687  if( WriteNC::ENTLOCSET == currentVarData.entLoc ) continue;
688 
689  // Set up writeStarts and writeCounts (maximum number of dimensions is 4)
690  currentVarData.writeStarts.resize( 4 );
691  currentVarData.writeCounts.resize( 4 );
692  unsigned int dim_idx = 0;
693 
694  // First: time
695  if( currentVarData.has_tsteps )
696  {
697  // Non-set variables with timesteps
698  // 4 dimensions like (time, lev, lat, lon)
699  // 3 dimensions like (time, lat, lon)
700  assert( 4 == varDims.size() || 3 == varDims.size() );
701 
702  // Time should be the first dimension
703  assert( tDim == varDims[0] );
704 
705  currentVarData.writeStarts[dim_idx] = 0; // This value is timestep dependent, will be set later
706  currentVarData.writeCounts[dim_idx] = 1;
707  dim_idx++;
708  }
709  else
710  {
711  // Non-set variables without timesteps
712  // 3 dimensions like (lev, lat, lon)
713  // 2 dimensions like (lat, lon)
714  assert( 3 == varDims.size() || 2 == varDims.size() );
715  }
716 
717  // Next: lev
718  if( currentVarData.numLev > 0 )
719  {
720  // Non-set variables with levels
721  // 4 dimensions like (time, lev, lat, lon)
722  // 3 dimensions like (lev, lat, lon)
723  assert( 4 == varDims.size() || 3 == varDims.size() );
724 
725  currentVarData.writeStarts[dim_idx] = 0;
726  currentVarData.writeCounts[dim_idx] = currentVarData.numLev;
727  dim_idx++;
728  }
729  else
730  {
731  // Non-set variables without levels
732  // 3 dimensions like (time, lat, lon)
733  // 2 dimensions like (lat, lon)
734  assert( 3 == varDims.size() || 2 == varDims.size() );
735  }
736 
737  // Finally: lat and lon
738  switch( currentVarData.entLoc )
739  {
740  case WriteNC::ENTLOCFACE:
741  // Faces
742  currentVarData.writeStarts[dim_idx] = lCDims[1];
743  currentVarData.writeCounts[dim_idx] = lCDims[4] - lCDims[1] + 1;
744  currentVarData.writeStarts[dim_idx + 1] = lCDims[0];
745  currentVarData.writeCounts[dim_idx + 1] = lCDims[3] - lCDims[0] + 1;
746  break;
747  default:
748  MB_SET_ERR( MB_FAILURE, "Unexpected entity location type for variable " << varname );
749  }
750  dim_idx += 2;
751 
752  // Get variable size
753  currentVarData.sz = 1;
754  for( std::size_t idx = 0; idx < dim_idx; idx++ )
755  currentVarData.sz *= currentVarData.writeCounts[idx];
756  } // for (size_t i = 0; i < var_names.size(); i++)
757 
758  return MB_SUCCESS;
759 }
760 
761 // Write CAM-EUL and CAM-FV non-set variables on non-shared quads (e.g. T)
762 // We assume that there are no variables on vertices and we do not support
763 // variables on edges (e.g. US in CAM-FV) for the time being
764 ErrorCode ScdNCWriteHelper::write_nonset_variables( std::vector< WriteNC::VarData >& vdatas,
765  std::vector< int >& tstep_nums )
766 {
767  Interface*& mbImpl = _writeNC->mbImpl;
768 
769  int success;
770 
771  // For each indexed variable tag, write a time step data
772  for( unsigned int i = 0; i < vdatas.size(); i++ )
773  {
774  WriteNC::VarData& variableData = vdatas[i];
775 
776  // Assume this variable is on faces for the time being
777  switch( variableData.entLoc )
778  {
779  case WriteNC::ENTLOCFACE:
780  // Faces
781  break;
782  default:
783  MB_SET_ERR( MB_FAILURE, "Unexpected entity location type for variable " << variableData.varName );
784  }
785 
786  unsigned int num_timesteps;
787  unsigned int lat_idx = 0;
788  unsigned int lon_idx = 1;
789  if( variableData.has_tsteps )
790  {
791  // Non-set variables with timesteps
792  // 4 dimensions like (time, lev, lat, lon)
793  // 3 dimensions like (time, lat, lon)
794  num_timesteps = tstep_nums.size();
795  lat_idx++;
796  lon_idx++;
797  }
798  else
799  {
800  // Non-set variables without timesteps
801  // 3 dimensions like (lev, lat, lon)
802  // 2 dimensions like (lat, lon)
803  num_timesteps = 1;
804  }
805 
806  unsigned int num_lev;
807  if( variableData.numLev > 0 )
808  {
809  // Non-set variables with levels
810  // 4 dimensions like (time, lev, lat, lon)
811  // 3 dimensions like (lev, lat, lon)
812  num_lev = variableData.numLev;
813  lat_idx++;
814  lon_idx++;
815  }
816  else
817  {
818  // Non-set variables without levels
819  // 3 dimensions like (time, lat, lon)
820  // 2 dimensions like (lat, lon)
821  num_lev = 1;
822  }
823 
824  size_t ni = variableData.writeCounts[lon_idx]; // lon
825  size_t nj = variableData.writeCounts[lat_idx]; // lat
826 
827  // At each timestep, we need to transpose tag format (lat, lon, lev) back
828  // to NC format (lev, lat, lon) for writing
829  for( unsigned int t = 0; t < num_timesteps; t++ )
830  {
831  // We will write one time step, and count will be one; start will be different
832  // Use tag_iterate to get tag data (assume that localCellsOwned is contiguous)
833  // We should also transpose for level so that means deep copy for transpose
834  if( tDim == variableData.varDims[0] ) variableData.writeStarts[0] = t; // This is start for time
835  int count;
836  void* dataptr;
837  ErrorCode rval = mbImpl->tag_iterate( variableData.varTags[t], localCellsOwned.begin(),
838  localCellsOwned.end(), count, dataptr );MB_CHK_SET_ERR( rval, "Failed to iterate tag on owned faces" );
839  assert( count == (int)localCellsOwned.size() );
840 
841  // Now transpose and write tag data
842  // Use collective I/O mode put (synchronous write) for the time being, we can try
843  // nonblocking put (request aggregation) later
844  switch( variableData.varDataType )
845  {
846  case NC_DOUBLE: {
847  std::vector< double > tmpdoubledata( ni * nj * num_lev );
848  if( num_lev > 1 )
849  // Transpose (lat, lon, lev) back to (lev, lat, lon)
850  jik_to_kji( ni, nj, num_lev, &tmpdoubledata[0], (double*)( dataptr ) );
851  success = NCFUNCAP( _vara_double )( _fileId, variableData.varId, &variableData.writeStarts[0],
852  &variableData.writeCounts[0], &tmpdoubledata[0] );
853  if( success )
854  MB_SET_ERR( MB_FAILURE, "Failed to write double data for variable " << variableData.varName );
855  break;
856  }
857  default:
858  MB_SET_ERR( MB_NOT_IMPLEMENTED, "Writing non-double data is not implemented yet" );
859  }
860  }
861  }
862 
863  return MB_SUCCESS;
864 }
865 
866 } /* namespace moab */