Mesh Oriented datABase  (version 5.5.1)
An array-based unstructured mesh library
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 
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(), "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 */