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
NCWriteHOMME.cpp
Go to the documentation of this file.
1 /* 2  * NCWriteHOMME.cpp 3  * 4  * Created on: April 9, 2014 5  */ 6  7 #include "NCWriteHOMME.hpp" 8 #include "MBTagConventions.hpp" 9  10 namespace moab 11 { 12  13 NCWriteHOMME::~NCWriteHOMME() 14 { 15  // TODO Auto-generated destructor stub 16 } 17  18 ErrorCode NCWriteHOMME::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 32  { 33  MB_SET_ERR( MB_FAILURE, "Couldn't find 'time' dimension" ); 34  } 35  nTimeSteps = dimLens[tDim]; 36  37  // Get number of levels 38  if( ( vecIt = std::find( dimNames.begin(), dimNames.end(), "lev" ) ) != dimNames.end() ) 39  levDim = vecIt - dimNames.begin(); 40  else 41  { 42  MB_SET_ERR( MB_FAILURE, "Couldn't find 'lev' dimension" ); 43  } 44  nLevels = dimLens[levDim]; 45  46  // Get local vertices 47  rval = mbImpl->get_entities_by_dimension( _fileSet, 0, localVertsOwned );MB_CHK_SET_ERR( rval, "Trouble getting local vertices in current file set" ); 48  assert( !localVertsOwned.empty() ); 49  50 #ifdef MOAB_HAVE_MPI 51  bool& isParallel = _writeNC->isParallel; 52  if( isParallel ) 53  { 54  ParallelComm*& myPcomm = _writeNC->myPcomm; 55  int rank = myPcomm->proc_config().proc_rank(); 56  int procs = myPcomm->proc_config().proc_size(); 57  if( procs > 1 ) 58  { 59 #ifndef NDEBUG 60  unsigned int num_local_verts = localVertsOwned.size(); 61 #endif 62  rval = myPcomm->filter_pstatus( localVertsOwned, PSTATUS_NOT_OWNED, PSTATUS_NOT );MB_CHK_SET_ERR( rval, "Trouble getting owned vertices in current set" ); 63  64  // Assume that PARALLEL_RESOLVE_SHARED_ENTS option is set 65  // Verify that not all local vertices are owned by the last processor 66  if( procs - 1 == rank ) 67  assert( "PARALLEL_RESOLVE_SHARED_ENTS option is set" && localVertsOwned.size() < num_local_verts ); 68  } 69  } 70 #endif 71  72  std::vector< int > gids( localVertsOwned.size() ); 73  rval = mbImpl->tag_get_data( mGlobalIdTag, localVertsOwned, &gids[0] );MB_CHK_SET_ERR( rval, "Trouble getting global IDs on local vertices" ); 74  75  // Get localGidVertsOwned 76  std::copy( gids.rbegin(), gids.rend(), range_inserter( localGidVertsOwned ) ); 77  78  return MB_SUCCESS; 79 } 80  81 ErrorCode NCWriteHOMME::collect_variable_data( std::vector< std::string >& var_names, std::vector< int >& tstep_nums ) 82 { 83  NCWriteHelper::collect_variable_data( var_names, tstep_nums ); 84  85  std::map< std::string, WriteNC::VarData >& varInfo = _writeNC->varInfo; 86  87  for( size_t i = 0; i < var_names.size(); i++ ) 88  { 89  std::string varname = var_names[i]; 90  std::map< std::string, WriteNC::VarData >::iterator vit = varInfo.find( varname ); 91  if( vit == varInfo.end() ) MB_SET_ERR( MB_FAILURE, "Can't find variable " << varname ); 92  ; 93  94  WriteNC::VarData& currentVarData = vit->second; 95 #ifndef NDEBUG 96  std::vector< int >& varDims = currentVarData.varDims; 97 #endif 98  99  // Skip set variables, which were already processed in 100  // NCWriteHelper::collect_variable_data() 101  if( WriteNC::ENTLOCSET == currentVarData.entLoc ) continue; 102  103  // Set up writeStarts and writeCounts (maximum number of dimensions is 3) 104  currentVarData.writeStarts.resize( 3 ); 105  currentVarData.writeCounts.resize( 3 ); 106  unsigned int dim_idx = 0; 107  108  // First: time 109  if( currentVarData.has_tsteps ) 110  { 111  // Non-set variables with timesteps 112  // 3 dimensions like (time, lev, ncol) 113  // 2 dimensions like (time, ncol) 114  assert( 3 == varDims.size() || 2 == varDims.size() ); 115  116  // Time should be the first dimension 117  assert( tDim == varDims[0] ); 118  119  currentVarData.writeStarts[dim_idx] = 0; // This value is timestep dependent, will be set later 120  currentVarData.writeCounts[dim_idx] = 1; 121  dim_idx++; 122  } 123  else 124  { 125  // Non-set variables without timesteps 126  // 2 dimensions like (lev, ncol) 127  // 1 dimension like (ncol) 128  assert( 2 == varDims.size() || 1 == varDims.size() ); 129  } 130  131  // Next: lev 132  if( currentVarData.numLev > 0 ) 133  { 134  // Non-set variables with levels 135  // 3 dimensions like (time, lev, ncol) 136  // 2 dimensions like (lev, ncol) 137  assert( 3 == varDims.size() || 2 == varDims.size() ); 138  139  currentVarData.writeStarts[dim_idx] = 0; 140  currentVarData.writeCounts[dim_idx] = currentVarData.numLev; 141  dim_idx++; 142  } 143  else 144  { 145  // Non-set variables without levels 146  // 2 dimensions like (time, ncol) 147  // 1 dimension like (ncol) 148  assert( 2 == varDims.size() || 1 == varDims.size() ); 149  } 150  151  // Finally: ncol 152  switch( currentVarData.entLoc ) 153  { 154  case WriteNC::ENTLOCVERT: 155  // Vertices 156  // Start from the first localGidVerts 157  // Actually, this will be reset later for writing 158  currentVarData.writeStarts[dim_idx] = localGidVertsOwned[0] - 1; 159  currentVarData.writeCounts[dim_idx] = localGidVertsOwned.size(); 160  break; 161  default: 162  MB_SET_ERR( MB_FAILURE, "Unexpected entity location type for variable " << varname ); 163  } 164  dim_idx++; 165  166  // Get variable size 167  currentVarData.sz = 1; 168  for( std::size_t idx = 0; idx < dim_idx; idx++ ) 169  currentVarData.sz *= currentVarData.writeCounts[idx]; 170  } 171  172  return MB_SUCCESS; 173 } 174  175 ErrorCode NCWriteHOMME::write_nonset_variables( std::vector< WriteNC::VarData >& vdatas, 176  std::vector< int >& tstep_nums ) 177 { 178  Interface*& mbImpl = _writeNC->mbImpl; 179  180  int success; 181  int num_local_verts_owned = localVertsOwned.size(); 182  183  // For each indexed variable tag, write a time step data 184  for( unsigned int i = 0; i < vdatas.size(); i++ ) 185  { 186  WriteNC::VarData& variableData = vdatas[i]; 187  188  // Assume this variable is on vertices for the time being 189  switch( variableData.entLoc ) 190  { 191  case WriteNC::ENTLOCVERT: 192  // Vertices 193  break; 194  default: 195  MB_SET_ERR( MB_FAILURE, "Unexpected entity location type for variable " << variableData.varName ); 196  } 197  198  unsigned int num_timesteps; 199  unsigned int ncol_idx = 0; 200  if( variableData.has_tsteps ) 201  { 202  // Non-set variables with timesteps 203  // 3 dimensions like (time, lev, ncol) 204  // 2 dimensions like (time, ncol) 205  num_timesteps = tstep_nums.size(); 206  ncol_idx++; 207  } 208  else 209  { 210  // Non-set variables without timesteps 211  // 2 dimensions like (lev, ncol) 212  // 1 dimension like (ncol) 213  num_timesteps = 1; 214  } 215  216  unsigned int num_lev; 217  if( variableData.numLev > 0 ) 218  { 219  // Non-set variables with levels 220  // 3 dimensions like (time, lev, ncol) 221  // 2 dimensions like (lev, ncol) 222  num_lev = variableData.numLev; 223  ncol_idx++; 224  } 225  else 226  { 227  // Non-set variables without levels 228  // 2 dimensions like (time, ncol) 229  // 1 dimension like (ncol) 230  num_lev = 1; 231  } 232  233  // At each timestep, we need to transpose tag format (ncol, lev) back 234  // to NC format (lev, ncol) for writing 235  for( unsigned int t = 0; t < num_timesteps; t++ ) 236  { 237  // We will write one time step, and count will be one; start will be different 238  // Use tag_get_data instead of tag_iterate to copy tag data, as localVertsOwned 239  // might not be contiguous. We should also transpose for level so that means 240  // deep copy for transpose 241  if( tDim == variableData.varDims[0] ) variableData.writeStarts[0] = t; // This is start for time 242  std::vector< double > tag_data( num_local_verts_owned * num_lev ); 243  ErrorCode rval = mbImpl->tag_get_data( variableData.varTags[t], localVertsOwned, &tag_data[0] );MB_CHK_SET_ERR( rval, "Trouble getting tag data on owned vertices" ); 244  245 #ifdef MOAB_HAVE_PNETCDF 246  size_t nb_writes = localGidVertsOwned.psize(); 247  std::vector< int > requests( nb_writes ), statuss( nb_writes ); 248  size_t idxReq = 0; 249 #endif 250  251  // Now transpose and write copied tag data 252  // Use nonblocking put (request aggregation) 253  switch( variableData.varDataType ) 254  { 255  case NC_DOUBLE: { 256  std::vector< double > tmpdoubledata( num_local_verts_owned * num_lev ); 257  if( num_lev > 1 ) 258  { 259  // Transpose (ncol, lev) back to (lev, ncol) 260  // Note, num_local_verts_owned is not used by jik_to_kji_stride() 261  jik_to_kji_stride( num_local_verts_owned, 1, num_lev, &tmpdoubledata[0], &tag_data[0], 262  localGidVertsOwned ); 263  } 264  265  size_t indexInDoubleArray = 0; 266  size_t ic = 0; 267  for( Range::pair_iterator pair_iter = localGidVertsOwned.pair_begin(); 268  pair_iter != localGidVertsOwned.pair_end(); ++pair_iter, ic++ ) 269  { 270  EntityHandle starth = pair_iter->first; 271  EntityHandle endh = pair_iter->second; 272  variableData.writeStarts[ncol_idx] = (NCDF_SIZE)( starth - 1 ); 273  variableData.writeCounts[ncol_idx] = (NCDF_SIZE)( endh - starth + 1 ); 274  275  // Do a partial write, in each subrange 276 #ifdef MOAB_HAVE_PNETCDF 277  // Wait outside this loop 278  success = 279  NCFUNCREQP( _vara_double )( _fileId, variableData.varId, &( variableData.writeStarts[0] ), 280  &( variableData.writeCounts[0] ), 281  &( tmpdoubledata[indexInDoubleArray] ), &requests[idxReq++] ); 282 #else 283  success = NCFUNCAP( 284  _vara_double )( _fileId, variableData.varId, &( variableData.writeStarts[0] ), 285  &( variableData.writeCounts[0] ), &( tmpdoubledata[indexInDoubleArray] ) ); 286 #endif 287  if( success ) 288  MB_SET_ERR( MB_FAILURE, 289  "Failed to write double data in a loop for variable " << variableData.varName ); 290  // We need to increment the index in double array for the 291  // next subrange 292  indexInDoubleArray += ( endh - starth + 1 ) * num_lev; 293  } 294  assert( ic == localGidVertsOwned.psize() ); 295 #ifdef MOAB_HAVE_PNETCDF 296  success = ncmpi_wait_all( _fileId, requests.size(), &requests[0], &statuss[0] ); 297  if( success ) MB_SET_ERR( MB_FAILURE, "Failed on wait_all" ); 298 #endif 299  break; 300  } 301  default: 302  MB_SET_ERR( MB_NOT_IMPLEMENTED, "Writing non-double data is not implemented yet" ); 303  } 304  } 305  } 306  307  return MB_SUCCESS; 308 } 309  310 } /* namespace moab */