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
NCWriteGCRM.cpp
Go to the documentation of this file.
1 /* 2  * NCWriteGCRM.cpp 3  * 4  * Created on: April 9, 2014 5  */ 6  7 #include "NCWriteGCRM.hpp" 8 #include "MBTagConventions.hpp" 9  10 namespace moab 11 { 12  13 NCWriteGCRM::~NCWriteGCRM() 14 { 15  // TODO Auto-generated destructor stub 16 } 17  18 ErrorCode NCWriteGCRM::collect_mesh_info() 19 { 20  Interface*& mbImpl = _writeNC->mbImpl; 21  std::vector< std::string >& dimNames = _writeNC->dimNames; 22  std::vector< int >& dimLens = _writeNC->dimLens; 23  Tag& mGlobalIdTag = _writeNC->mGlobalIdTag; 24  25  ErrorCode rval; 26  27  // Look for time dimension 28  std::vector< std::string >::iterator vecIt; 29  if( ( vecIt = std::find( dimNames.begin(), dimNames.end(), "Time" ) ) != dimNames.end() ) 30  tDim = vecIt - dimNames.begin(); 31  else if( ( vecIt = std::find( dimNames.begin(), dimNames.end(), "time" ) ) != dimNames.end() ) 32  tDim = vecIt - dimNames.begin(); 33  else 34  { 35  MB_SET_ERR( MB_FAILURE, "Couldn't find 'Time' or 'time' dimension" ); 36  } 37  nTimeSteps = dimLens[tDim]; 38  39  // Get number of levels 40  if( ( vecIt = std::find( dimNames.begin(), dimNames.end(), "layers" ) ) != dimNames.end() ) 41  levDim = vecIt - dimNames.begin(); 42  else 43  { 44  MB_SET_ERR( MB_FAILURE, "Couldn't find 'layers' dimension" ); 45  } 46  nLevels = dimLens[levDim]; 47  48  // Get local vertices 49  rval = mbImpl->get_entities_by_dimension( _fileSet, 0, localVertsOwned );MB_CHK_SET_ERR( rval, "Trouble getting local vertices in current file set" ); 50  assert( !localVertsOwned.empty() ); 51  52  // Get local edges 53  rval = mbImpl->get_entities_by_dimension( _fileSet, 1, localEdgesOwned );MB_CHK_SET_ERR( rval, "Trouble getting local edges in current file set" ); 54  // There are no edges if NO_EDGES read option is set 55  56  // Get local cells 57  rval = mbImpl->get_entities_by_dimension( _fileSet, 2, localCellsOwned );MB_CHK_SET_ERR( rval, "Trouble getting local cells in current file set" ); 58  assert( !localCellsOwned.empty() ); 59  60 #ifdef MOAB_HAVE_MPI 61  bool& isParallel = _writeNC->isParallel; 62  if( isParallel ) 63  { 64  ParallelComm*& myPcomm = _writeNC->myPcomm; 65  int rank = myPcomm->proc_config().proc_rank(); 66  int procs = myPcomm->proc_config().proc_size(); 67  if( procs > 1 ) 68  { 69 #ifndef NDEBUG 70  unsigned int num_local_verts = localVertsOwned.size(); 71 #endif 72  rval = myPcomm->filter_pstatus( localVertsOwned, PSTATUS_NOT_OWNED, PSTATUS_NOT );MB_CHK_SET_ERR( rval, "Trouble getting owned vertices in current file set" ); 73  74  // Assume that PARALLEL_RESOLVE_SHARED_ENTS option is set 75  // Verify that not all local vertices are owned by the last processor 76  if( procs - 1 == rank ) 77  assert( "PARALLEL_RESOLVE_SHARED_ENTS option is set" && localVertsOwned.size() < num_local_verts ); 78  79  if( !localEdgesOwned.empty() ) 80  { 81  rval = myPcomm->filter_pstatus( localEdgesOwned, PSTATUS_NOT_OWNED, PSTATUS_NOT );MB_CHK_SET_ERR( rval, "Trouble getting owned edges in current file set" ); 82  } 83  84  rval = myPcomm->filter_pstatus( localCellsOwned, PSTATUS_NOT_OWNED, PSTATUS_NOT );MB_CHK_SET_ERR( rval, "Trouble getting owned cells in current file set" ); 85  } 86  } 87 #endif 88  89  std::vector< int > gids( localVertsOwned.size() ); 90  rval = mbImpl->tag_get_data( mGlobalIdTag, localVertsOwned, &gids[0] );MB_CHK_SET_ERR( rval, "Trouble getting global IDs on local vertices" ); 91  92  // Get localGidVertsOwned 93  std::copy( gids.rbegin(), gids.rend(), range_inserter( localGidVertsOwned ) ); 94  95  if( !localEdgesOwned.empty() ) 96  { 97  gids.resize( localEdgesOwned.size() ); 98  rval = mbImpl->tag_get_data( mGlobalIdTag, localEdgesOwned, &gids[0] );MB_CHK_SET_ERR( rval, "Trouble getting global IDs on local edges" ); 99  100  // Get localGidEdgesOwned 101  std::copy( gids.rbegin(), gids.rend(), range_inserter( localGidEdgesOwned ) ); 102  } 103  104  gids.resize( localCellsOwned.size() ); 105  rval = mbImpl->tag_get_data( mGlobalIdTag, localCellsOwned, &gids[0] );MB_CHK_SET_ERR( rval, "Trouble getting global IDs on local cells" ); 106  107  // Get localGidCellsOwned 108  std::copy( gids.rbegin(), gids.rend(), range_inserter( localGidCellsOwned ) ); 109  110  return MB_SUCCESS; 111 } 112  113 ErrorCode NCWriteGCRM::collect_variable_data( std::vector< std::string >& var_names, std::vector< int >& tstep_nums ) 114 { 115  NCWriteHelper::collect_variable_data( var_names, tstep_nums ); 116  117  std::vector< std::string >& dimNames = _writeNC->dimNames; 118  std::vector< int >& dimLens = _writeNC->dimLens; 119  120  // Dimension indices for other optional levels 121  std::vector< unsigned int > opt_lev_dims; 122  123  unsigned int lev_idx; 124  std::vector< std::string >::iterator vecIt; 125  126  // Get index of interface levels 127  if( ( vecIt = std::find( dimNames.begin(), dimNames.end(), "interfaces" ) ) != dimNames.end() ) 128  { 129  lev_idx = vecIt - dimNames.begin(); 130  opt_lev_dims.push_back( lev_idx ); 131  } 132  133  std::map< std::string, WriteNC::VarData >& varInfo = _writeNC->varInfo; 134  135  for( size_t i = 0; i < var_names.size(); i++ ) 136  { 137  std::string varname = var_names[i]; 138  std::map< std::string, WriteNC::VarData >::iterator vit = varInfo.find( varname ); 139  if( vit == varInfo.end() ) MB_SET_ERR( MB_FAILURE, "Can't find variable " << varname ); 140  141  WriteNC::VarData& currentVarData = vit->second; 142  std::vector< int >& varDims = currentVarData.varDims; 143  144  // Skip edge variables, if there are no edges 145  if( localEdgesOwned.empty() && currentVarData.entLoc == WriteNC::ENTLOCEDGE ) continue; 146  147  // If layers dimension is not found, try other optional levels such as interfaces 148  if( std::find( varDims.begin(), varDims.end(), levDim ) == varDims.end() ) 149  { 150  for( unsigned int j = 0; j < opt_lev_dims.size(); j++ ) 151  { 152  if( std::find( varDims.begin(), varDims.end(), opt_lev_dims[j] ) != varDims.end() ) 153  { 154  currentVarData.numLev = dimLens[opt_lev_dims[j]]; 155  break; 156  } 157  } 158  } 159  160  // Skip set variables, which were already processed in 161  // NCWriteHelper::collect_variable_data() 162  if( WriteNC::ENTLOCSET == currentVarData.entLoc ) continue; 163  164  // Set up writeStarts and writeCounts (maximum number of dimensions is 3) 165  currentVarData.writeStarts.resize( 3 ); 166  currentVarData.writeCounts.resize( 3 ); 167  unsigned int dim_idx = 0; 168  169  // First: time 170  if( currentVarData.has_tsteps ) 171  { 172  // Non-set variables with timesteps 173  // 3 dimensions like (time, cells, layers) 174  // 2 dimensions like (time, cells) 175  assert( 3 == varDims.size() || 2 == varDims.size() ); 176  177  // Time should be the first dimension 178  assert( tDim == varDims[0] ); 179  180  currentVarData.writeStarts[dim_idx] = 0; // This value is timestep dependent, will be set later 181  currentVarData.writeCounts[dim_idx] = 1; 182  dim_idx++; 183  } 184  else 185  { 186  // Non-set variables without timesteps 187  // 2 dimensions like (cells, layers) 188  // 1 dimension like (cells) 189  assert( 2 == varDims.size() || 1 == varDims.size() ); 190  } 191  192  // Next: corners / cells / edges 193  switch( currentVarData.entLoc ) 194  { 195  case WriteNC::ENTLOCVERT: 196  // Vertices 197  // Start from the first localGidVerts 198  // Actually, this will be reset later for writing 199  currentVarData.writeStarts[dim_idx] = localGidVertsOwned[0] - 1; 200  currentVarData.writeCounts[dim_idx] = localGidVertsOwned.size(); 201  break; 202  case WriteNC::ENTLOCFACE: 203  // Faces 204  // Start from the first localGidCells 205  // Actually, this will be reset later for writing 206  currentVarData.writeStarts[dim_idx] = localGidCellsOwned[0] - 1; 207  currentVarData.writeCounts[dim_idx] = localGidCellsOwned.size(); 208  break; 209  case WriteNC::ENTLOCEDGE: 210  // Edges 211  // Start from the first localGidEdges 212  // Actually, this will be reset later for writing 213  currentVarData.writeStarts[dim_idx] = localGidEdgesOwned[0] - 1; 214  currentVarData.writeCounts[dim_idx] = localGidEdgesOwned.size(); 215  break; 216  default: 217  MB_SET_ERR( MB_FAILURE, "Unexpected entity location type for variable " << varname ); 218  } 219  dim_idx++; 220  221  // Finally: layers or other optional levels, it is possible that there is no 222  // level dimension (numLev is 0) for this non-set variable 223  if( currentVarData.numLev > 0 ) 224  { 225  // Non-set variables with levels 226  // 3 dimensions like (time, cells, layers) 227  // 2 dimensions like (cells, layers) 228  assert( 3 == varDims.size() || 2 == varDims.size() ); 229  230  currentVarData.writeStarts[dim_idx] = 0; 231  currentVarData.writeCounts[dim_idx] = currentVarData.numLev; 232  dim_idx++; 233  } 234  else 235  { 236  // Non-set variables without levels 237  // 2 dimensions like (time, cells) 238  // 1 dimension like (cells) 239  assert( 2 == varDims.size() || 1 == varDims.size() ); 240  } 241  242  // Get variable size 243  currentVarData.sz = 1; 244  for( std::size_t idx = 0; idx < dim_idx; idx++ ) 245  currentVarData.sz *= currentVarData.writeCounts[idx]; 246  } // for (size_t i = 0; i < var_names.size(); i++) 247  248  return MB_SUCCESS; 249 } 250  251 ErrorCode NCWriteGCRM::write_nonset_variables( std::vector< WriteNC::VarData >& vdatas, std::vector< int >& tstep_nums ) 252 { 253  Interface*& mbImpl = _writeNC->mbImpl; 254  255  int success; 256  257  // For each indexed variable tag, write a time step data 258  for( unsigned int i = 0; i < vdatas.size(); i++ ) 259  { 260  WriteNC::VarData& variableData = vdatas[i]; 261  262  // Skip edge variables, if there are no edges 263  if( localEdgesOwned.empty() && variableData.entLoc == WriteNC::ENTLOCEDGE ) continue; 264  265  // Get local owned entities of this variable 266  Range* pLocalEntsOwned = NULL; 267  Range* pLocalGidEntsOwned = NULL; 268  switch( variableData.entLoc ) 269  { 270  case WriteNC::ENTLOCVERT: 271  // Vertices 272  pLocalEntsOwned = &localVertsOwned; 273  pLocalGidEntsOwned = &localGidVertsOwned; 274  break; 275  case WriteNC::ENTLOCEDGE: 276  // Edges 277  pLocalEntsOwned = &localEdgesOwned; 278  pLocalGidEntsOwned = &localGidEdgesOwned; 279  break; 280  case WriteNC::ENTLOCFACE: 281  // Cells 282  pLocalEntsOwned = &localCellsOwned; 283  pLocalGidEntsOwned = &localGidCellsOwned; 284  break; 285  default: 286  MB_SET_ERR( MB_FAILURE, "Unexpected entity location type for variable " << variableData.varName ); 287  } 288  289  unsigned int num_timesteps; 290  unsigned int ents_idx = 0; 291  if( variableData.has_tsteps ) 292  { 293  // Non-set variables with timesteps 294  // 3 dimensions like (time, cells, layers) 295  // 2 dimensions like (time, cells) 296  num_timesteps = tstep_nums.size(); 297  ents_idx++; 298  } 299  else 300  { 301  // Non-set variables without timesteps 302  // 2 dimensions like (cells, layers) 303  // 1 dimension like (cells) 304  num_timesteps = 1; 305  } 306  307  unsigned int num_lev; 308  if( variableData.numLev > 0 ) 309  { 310  // Non-set variables with levels 311  // 3 dimensions like (time, cells, layers) 312  // 2 dimensions like (cells, layers) 313  num_lev = variableData.numLev; 314  } 315  else 316  { 317  // Non-set variables without levels 318  // 2 dimensions like (time, cells) 319  // 1 dimension like (cells) 320  num_lev = 1; 321  } 322  323  for( unsigned int t = 0; t < num_timesteps; t++ ) 324  { 325  // We will write one time step, and count will be one; start will be different 326  // Use tag_get_data instead of tag_iterate to copy tag data, as localEntsOwned 327  // might not be contiguous. 328  if( tDim == variableData.varDims[0] ) variableData.writeStarts[0] = t; // This is start for time 329  std::vector< double > tag_data( pLocalEntsOwned->size() * num_lev ); 330  ErrorCode rval = mbImpl->tag_get_data( variableData.varTags[t], *pLocalEntsOwned, &tag_data[0] );MB_CHK_SET_ERR( rval, "Trouble getting tag data on owned entities" ); 331  332 #ifdef MOAB_HAVE_PNETCDF 333  size_t nb_writes = pLocalGidEntsOwned->psize(); 334  std::vector< int > requests( nb_writes ), statuss( nb_writes ); 335  size_t idxReq = 0; 336 #endif 337  338  // Now write copied tag data 339  // Use nonblocking put (request aggregation) 340  switch( variableData.varDataType ) 341  { 342  case NC_DOUBLE: { 343  size_t indexInDoubleArray = 0; 344  size_t ic = 0; 345  for( Range::pair_iterator pair_iter = pLocalGidEntsOwned->pair_begin(); 346  pair_iter != pLocalGidEntsOwned->pair_end(); ++pair_iter, ic++ ) 347  { 348  EntityHandle starth = pair_iter->first; 349  EntityHandle endh = pair_iter->second; 350  variableData.writeStarts[ents_idx] = (NCDF_SIZE)( starth - 1 ); 351  variableData.writeCounts[ents_idx] = (NCDF_SIZE)( endh - starth + 1 ); 352  353  // Do a partial write, in each subrange 354 #ifdef MOAB_HAVE_PNETCDF 355  // Wait outside this loop 356  success = 357  NCFUNCREQP( _vara_double )( _fileId, variableData.varId, &( variableData.writeStarts[0] ), 358  &( variableData.writeCounts[0] ), 359  &( tag_data[indexInDoubleArray] ), &requests[idxReq++] ); 360 #else 361  success = NCFUNCAP( 362  _vara_double )( _fileId, variableData.varId, &( variableData.writeStarts[0] ), 363  &( variableData.writeCounts[0] ), &( tag_data[indexInDoubleArray] ) ); 364 #endif 365  if( success ) 366  MB_SET_ERR( MB_FAILURE, 367  "Failed to write double data in a loop for variable " << variableData.varName ); 368  // We need to increment the index in double array for the 369  // next subrange 370  indexInDoubleArray += ( endh - starth + 1 ) * num_lev; 371  } 372  assert( ic == pLocalGidEntsOwned->psize() ); 373 #ifdef MOAB_HAVE_PNETCDF 374  success = ncmpi_wait_all( _fileId, requests.size(), &requests[0], &statuss[0] ); 375  if( success ) MB_SET_ERR( MB_FAILURE, "Failed on wait_all" ); 376 #endif 377  break; 378  } 379  default: 380  MB_SET_ERR( MB_NOT_IMPLEMENTED, "Writing non-double data is not implemented yet" ); 381  } 382  } 383  } 384  385  return MB_SUCCESS; 386 } 387  388 } /* namespace moab */