Mesh Oriented datABase  (version 5.5.1)
An array-based unstructured mesh library
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 
14 {
15  // TODO Auto-generated destructor stub
16 }
17 
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 */