Loading [MathJax]/extensions/tex2jax.js
Mesh Oriented datABase  (version 5.5.1)
An array-based unstructured mesh library
All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
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 30 NCWriteHelper* NCWriteHelper::get_nc_helper( WriteNC* writeNC, 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  593 ErrorCode ScdNCWriteHelper::collect_mesh_info() 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 */