1
6
7 #include "NCWriteGCRM.hpp"
8 #include "MBTagConventions.hpp"
9
10 namespace moab
11 {
12
13 NCWriteGCRM::~NCWriteGCRM()
14 {
15
16 }
17
18 ErrorCode NCWriteGCRM::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
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
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
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
53 rval = mbImpl->get_entities_by_dimension( _fileSet, 1, localEdgesOwned );MB_CHK_SET_ERR( rval, "Trouble getting local edges in current file set" );
54
55
56
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
75
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
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
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
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
121 std::vector< unsigned int > opt_lev_dims;
122
123 unsigned int lev_idx;
124 std::vector< std::string >::iterator vecIt;
125
126
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
145 if( localEdgesOwned.empty() && currentVarData.entLoc == WriteNC::ENTLOCEDGE ) continue;
146
147
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
161
162 if( WriteNC::ENTLOCSET == currentVarData.entLoc ) continue;
163
164
165 currentVarData.writeStarts.resize( 3 );
166 currentVarData.writeCounts.resize( 3 );
167 unsigned int dim_idx = 0;
168
169
170 if( currentVarData.has_tsteps )
171 {
172
173
174
175 assert( 3 == varDims.size() || 2 == varDims.size() );
176
177
178 assert( tDim == varDims[0] );
179
180 currentVarData.writeStarts[dim_idx] = 0;
181 currentVarData.writeCounts[dim_idx] = 1;
182 dim_idx++;
183 }
184 else
185 {
186
187
188
189 assert( 2 == varDims.size() || 1 == varDims.size() );
190 }
191
192
193 switch( currentVarData.entLoc )
194 {
195 case WriteNC::ENTLOCVERT:
196
197
198
199 currentVarData.writeStarts[dim_idx] = localGidVertsOwned[0] - 1;
200 currentVarData.writeCounts[dim_idx] = localGidVertsOwned.size();
201 break;
202 case WriteNC::ENTLOCFACE:
203
204
205
206 currentVarData.writeStarts[dim_idx] = localGidCellsOwned[0] - 1;
207 currentVarData.writeCounts[dim_idx] = localGidCellsOwned.size();
208 break;
209 case WriteNC::ENTLOCEDGE:
210
211
212
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
222
223 if( currentVarData.numLev > 0 )
224 {
225
226
227
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
237
238
239 assert( 2 == varDims.size() || 1 == varDims.size() );
240 }
241
242
243 currentVarData.sz = 1;
244 for( std::size_t idx = 0; idx < dim_idx; idx++ )
245 currentVarData.sz *= currentVarData.writeCounts[idx];
246 }
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
258 for( unsigned int i = 0; i < vdatas.size(); i++ )
259 {
260 WriteNC::VarData& variableData = vdatas[i];
261
262
263 if( localEdgesOwned.empty() && variableData.entLoc == WriteNC::ENTLOCEDGE ) continue;
264
265
266 Range* pLocalEntsOwned = NULL;
267 Range* pLocalGidEntsOwned = NULL;
268 switch( variableData.entLoc )
269 {
270 case WriteNC::ENTLOCVERT:
271
272 pLocalEntsOwned = &localVertsOwned;
273 pLocalGidEntsOwned = &localGidVertsOwned;
274 break;
275 case WriteNC::ENTLOCEDGE:
276
277 pLocalEntsOwned = &localEdgesOwned;
278 pLocalGidEntsOwned = &localGidEdgesOwned;
279 break;
280 case WriteNC::ENTLOCFACE:
281
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
294
295
296 num_timesteps = tstep_nums.size();
297 ents_idx++;
298 }
299 else
300 {
301
302
303
304 num_timesteps = 1;
305 }
306
307 unsigned int num_lev;
308 if( variableData.numLev > 0 )
309 {
310
311
312
313 num_lev = variableData.numLev;
314 }
315 else
316 {
317
318
319
320 num_lev = 1;
321 }
322
323 for( unsigned int t = 0; t < num_timesteps; t++ )
324 {
325
326
327
328 if( tDim == variableData.varDims[0] ) variableData.writeStarts[0] = t;
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
339
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
354 #ifdef MOAB_HAVE_PNETCDF
355
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
369
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 }