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
NCWriteMPAS.cpp
Go to the documentation of this file.
1 /* 2  * NCWriteMPAS.cpp 3  * 4  * Created on: April 9, 2014 5  */ 6  7 #include "NCWriteMPAS.hpp" 8 #include "MBTagConventions.hpp" 9  10 namespace moab 11 { 12  13 NCWriteMPAS::~NCWriteMPAS() 14 { 15  // TODO Auto-generated destructor stub 16 } 17  18 ErrorCode NCWriteMPAS::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(), "nVertLevels" ) ) != dimNames.end() ) 41  levDim = vecIt - dimNames.begin(); 42  else 43  { 44  MB_SET_ERR( MB_FAILURE, "Couldn't find 'nVertLevels' 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 NCWriteMPAS::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 vertex levels P1 127  if( ( vecIt = std::find( dimNames.begin(), dimNames.end(), "nVertLevelsP1" ) ) != dimNames.end() ) 128  { 129  lev_idx = vecIt - dimNames.begin(); 130  opt_lev_dims.push_back( lev_idx ); 131  } 132  133  // Get index of vertex levels P2 134  if( ( vecIt = std::find( dimNames.begin(), dimNames.end(), "nVertLevelsP2" ) ) != dimNames.end() ) 135  { 136  lev_idx = vecIt - dimNames.begin(); 137  opt_lev_dims.push_back( lev_idx ); 138  } 139  140  // Get index of soil levels 141  if( ( vecIt = std::find( dimNames.begin(), dimNames.end(), "nSoilLevels" ) ) != dimNames.end() ) 142  { 143  lev_idx = vecIt - dimNames.begin(); 144  opt_lev_dims.push_back( lev_idx ); 145  } 146  147  std::map< std::string, WriteNC::VarData >& varInfo = _writeNC->varInfo; 148  149  for( size_t i = 0; i < var_names.size(); i++ ) 150  { 151  std::string varname = var_names[i]; 152  std::map< std::string, WriteNC::VarData >::iterator vit = varInfo.find( varname ); 153  if( vit == varInfo.end() ) MB_SET_ERR( MB_FAILURE, "Can't find variable " << varname ); 154  155  WriteNC::VarData& currentVarData = vit->second; 156  std::vector< int >& varDims = currentVarData.varDims; 157  158  // Skip edge variables, if there are no edges 159  if( localEdgesOwned.empty() && currentVarData.entLoc == WriteNC::ENTLOCEDGE ) continue; 160  161  // If nVertLevels dimension is not found, try other optional levels such as nVertLevelsP1 162  if( std::find( varDims.begin(), varDims.end(), levDim ) == varDims.end() ) 163  { 164  for( unsigned int j = 0; j < opt_lev_dims.size(); j++ ) 165  { 166  if( std::find( varDims.begin(), varDims.end(), opt_lev_dims[j] ) != varDims.end() ) 167  { 168  currentVarData.numLev = dimLens[opt_lev_dims[j]]; 169  break; 170  } 171  } 172  } 173  174  // Skip set variables, which were already processed in 175  // NCWriteHelper::collect_variable_data() 176  if( WriteNC::ENTLOCSET == currentVarData.entLoc ) continue; 177  178  // Set up writeStarts and writeCounts (maximum number of dimensions is 3) 179  currentVarData.writeStarts.resize( 3 ); 180  currentVarData.writeCounts.resize( 3 ); 181  unsigned int dim_idx = 0; 182  183  // First: Time 184  if( currentVarData.has_tsteps ) 185  { 186  // Non-set variables with timesteps 187  // 3 dimensions like (Time, nCells, nVertLevels) 188  // 2 dimensions like (Time, nCells) 189  assert( 3 == varDims.size() || 2 == varDims.size() ); 190  191  // Time should be the first dimension 192  assert( tDim == varDims[0] ); 193  194  currentVarData.writeStarts[dim_idx] = 0; // This value is timestep dependent, will be set later 195  currentVarData.writeCounts[dim_idx] = 1; 196  dim_idx++; 197  } 198  else 199  { 200  // Non-set variables without timesteps 201  // 2 dimensions like (nCells, nVertLevels) 202  // 1 dimension like (nCells) 203  assert( 2 == varDims.size() || 1 == varDims.size() ); 204  } 205  206  // Next: nVertices / nCells / nEdges 207  switch( currentVarData.entLoc ) 208  { 209  case WriteNC::ENTLOCVERT: 210  // Vertices 211  // Start from the first localGidVerts 212  // Actually, this will be reset later for writing 213  currentVarData.writeStarts[dim_idx] = localGidVertsOwned[0] - 1; 214  currentVarData.writeCounts[dim_idx] = localGidVertsOwned.size(); 215  break; 216  case WriteNC::ENTLOCFACE: 217  // Faces 218  // Start from the first localGidCells 219  // Actually, this will be reset later for writing 220  currentVarData.writeStarts[dim_idx] = localGidCellsOwned[0] - 1; 221  currentVarData.writeCounts[dim_idx] = localGidCellsOwned.size(); 222  break; 223  case WriteNC::ENTLOCEDGE: 224  // Edges 225  // Start from the first localGidEdges 226  // Actually, this will be reset later for writing 227  currentVarData.writeStarts[dim_idx] = localGidEdgesOwned[0] - 1; 228  currentVarData.writeCounts[dim_idx] = localGidEdgesOwned.size(); 229  break; 230  default: 231  MB_SET_ERR( MB_FAILURE, "Unexpected entity location type for variable " << varname ); 232  } 233  dim_idx++; 234  235  // Finally: nVertLevels or other optional levels, it is possible that there is no 236  // level dimension (numLev is 0) for this non-set variable, e.g. (Time, nCells) 237  if( currentVarData.numLev > 0 ) 238  { 239  // Non-set variables with levels 240  // 3 dimensions like (Time, nCells, nVertLevels) 241  // 2 dimensions like (nCells, nVertLevels) 242  assert( 3 == varDims.size() || 2 == varDims.size() ); 243  244  currentVarData.writeStarts[dim_idx] = 0; 245  currentVarData.writeCounts[dim_idx] = currentVarData.numLev; 246  dim_idx++; 247  } 248  else 249  { 250  // Non-set variables without levels 251  // 2 dimensions like (Time, nCells) 252  // 1 dimension like (nCells) 253  assert( 2 == varDims.size() || 1 == varDims.size() ); 254  } 255  256  // Get variable size 257  currentVarData.sz = 1; 258  for( std::size_t idx = 0; idx < dim_idx; idx++ ) 259  currentVarData.sz *= currentVarData.writeCounts[idx]; 260  } 261  262  return MB_SUCCESS; 263 } 264  265 ErrorCode NCWriteMPAS::write_nonset_variables( std::vector< WriteNC::VarData >& vdatas, std::vector< int >& tstep_nums ) 266 { 267  Interface*& mbImpl = _writeNC->mbImpl; 268  269  int success; 270  271  // For each variable tag in the indexed lists, write a time step data 272  for( unsigned int i = 0; i < vdatas.size(); i++ ) 273  { 274  WriteNC::VarData& variableData = vdatas[i]; 275  276  // Skip edge variables, if there are no edges 277  if( localEdgesOwned.empty() && variableData.entLoc == WriteNC::ENTLOCEDGE ) continue; 278  279  // Get local owned entities of this variable 280  Range* pLocalEntsOwned = NULL; 281  Range* pLocalGidEntsOwned = NULL; 282  switch( variableData.entLoc ) 283  { 284  case WriteNC::ENTLOCVERT: 285  // Vertices 286  pLocalEntsOwned = &localVertsOwned; 287  pLocalGidEntsOwned = &localGidVertsOwned; 288  break; 289  case WriteNC::ENTLOCEDGE: 290  // Edges 291  pLocalEntsOwned = &localEdgesOwned; 292  pLocalGidEntsOwned = &localGidEdgesOwned; 293  break; 294  case WriteNC::ENTLOCFACE: 295  // Cells 296  pLocalEntsOwned = &localCellsOwned; 297  pLocalGidEntsOwned = &localGidCellsOwned; 298  break; 299  default: 300  MB_SET_ERR( MB_FAILURE, "Unexpected entity location type for variable " << variableData.varName ); 301  } 302  303  unsigned int num_timesteps; 304  unsigned int ents_idx = 0; 305  if( variableData.has_tsteps ) 306  { 307  // Non-set variables with timesteps 308  // 3 dimensions like (Time, nCells, nVertLevels) 309  // 2 dimensions like (Time, nCells) 310  num_timesteps = tstep_nums.size(); 311  ents_idx++; 312  } 313  else 314  { 315  // Non-set variables without timesteps 316  // 2 dimensions like (nCells, nVertLevels) 317  // 1 dimension like (nCells) 318  num_timesteps = 1; 319  } 320  321  unsigned int num_lev; 322  if( variableData.numLev > 0 ) 323  { 324  // Non-set variables with levels 325  // 3 dimensions like (Time, nCells, nVertLevels) 326  // 2 dimensions like (nCells, nVertLevels) 327  num_lev = variableData.numLev; 328  } 329  else 330  { 331  // Non-set variables without levels 332  // 2 dimensions like (Time, nCells) 333  // 1 dimension like (nCells) 334  num_lev = 1; 335  } 336  337  for( unsigned int t = 0; t < num_timesteps; t++ ) 338  { 339  // We will write one time step, and count will be one; start will be different 340  // Use tag_get_data instead of tag_iterate to copy tag data, as localEntsOwned 341  // might not be contiguous. 342  if( tDim == variableData.varDims[0] ) variableData.writeStarts[0] = t; // This is start for time 343  std::vector< double > tag_data( pLocalEntsOwned->size() * num_lev ); 344  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" ); 345  346 #ifdef MOAB_HAVE_PNETCDF 347  size_t nb_writes = pLocalGidEntsOwned->psize(); 348  std::vector< int > requests( nb_writes ), statuss( nb_writes ); 349  size_t idxReq = 0; 350 #endif 351  352  // Now write copied tag data 353  // Use nonblocking put (request aggregation) 354  switch( variableData.varDataType ) 355  { 356  case NC_DOUBLE: { 357  size_t indexInDoubleArray = 0; 358  size_t ic = 0; 359  for( Range::pair_iterator pair_iter = pLocalGidEntsOwned->pair_begin(); 360  pair_iter != pLocalGidEntsOwned->pair_end(); ++pair_iter, ic++ ) 361  { 362  EntityHandle starth = pair_iter->first; 363  EntityHandle endh = pair_iter->second; 364  variableData.writeStarts[ents_idx] = (NCDF_SIZE)( starth - 1 ); 365  variableData.writeCounts[ents_idx] = (NCDF_SIZE)( endh - starth + 1 ); 366  367  // Do a partial write, in each subrange 368 #ifdef MOAB_HAVE_PNETCDF 369  // Wait outside this loop 370  success = 371  NCFUNCREQP( _vara_double )( _fileId, variableData.varId, &( variableData.writeStarts[0] ), 372  &( variableData.writeCounts[0] ), 373  &( tag_data[indexInDoubleArray] ), &requests[idxReq++] ); 374 #else 375  success = NCFUNCAP( 376  _vara_double )( _fileId, variableData.varId, &( variableData.writeStarts[0] ), 377  &( variableData.writeCounts[0] ), &( tag_data[indexInDoubleArray] ) ); 378 #endif 379  if( success ) 380  MB_SET_ERR( MB_FAILURE, 381  "Failed to write double data in a loop for variable " << variableData.varName ); 382  // We need to increment the index in double array for the 383  // next subrange 384  indexInDoubleArray += ( endh - starth + 1 ) * num_lev; 385  } 386  assert( ic == pLocalGidEntsOwned->psize() ); 387 #ifdef MOAB_HAVE_PNETCDF 388  success = ncmpi_wait_all( _fileId, requests.size(), &requests[0], &statuss[0] ); 389  if( success ) MB_SET_ERR( MB_FAILURE, "Failed on wait_all" ); 390 #endif 391  break; 392  } 393  default: 394  MB_SET_ERR( MB_NOT_IMPLEMENTED, "Writing non-double data is not implemented yet" ); 395  } 396  } 397  } 398  399  return MB_SUCCESS; 400 } 401  402 } /* namespace moab */