Mesh Oriented datABase  (version 5.5.1)
An array-based unstructured mesh library
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 
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
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],
263  }
264 
265  size_t indexInDoubleArray = 0;
266  size_t ic = 0;
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 */