MOAB: Mesh Oriented datABase
(version 5.3.1)
|
00001 /* 00002 * NCWriteHelper.cpp 00003 * 00004 * Created on: Mar 28, 2014 00005 * Author: iulian 00006 */ 00007 00008 #include "NCWriteHelper.hpp" 00009 #include "NCWriteEuler.hpp" 00010 #include "NCWriteFV.hpp" 00011 #include "NCWriteHOMME.hpp" 00012 #include "NCWriteMPAS.hpp" 00013 #include "NCWriteGCRM.hpp" 00014 00015 #include "moab/WriteUtilIface.hpp" 00016 #include "MBTagConventions.hpp" 00017 00018 #include <sstream> 00019 00020 #ifdef WIN32 00021 #ifdef size_t 00022 #undef size_t 00023 #endif 00024 #endif 00025 00026 namespace moab 00027 { 00028 00029 //! Get appropriate helper instance for WriteNC class; based on some info in the file set 00030 NCWriteHelper* NCWriteHelper::get_nc_helper( WriteNC* writeNC, int fileId, const FileOptions& opts, 00031 EntityHandle fileSet ) 00032 { 00033 std::string& grid_type = writeNC->grid_type; 00034 if( grid_type == "CAM_EUL" ) 00035 return new( std::nothrow ) NCWriteEuler( writeNC, fileId, opts, fileSet ); 00036 else if( grid_type == "CAM_FV" ) 00037 return new( std::nothrow ) NCWriteFV( writeNC, fileId, opts, fileSet ); 00038 else if( grid_type == "CAM_SE" ) 00039 return new( std::nothrow ) NCWriteHOMME( writeNC, fileId, opts, fileSet ); 00040 else if( grid_type == "MPAS" ) 00041 return new( std::nothrow ) NCWriteMPAS( writeNC, fileId, opts, fileSet ); 00042 else if( grid_type == "GCRM" ) 00043 return new( std::nothrow ) NCWriteGCRM( writeNC, fileId, opts, fileSet ); 00044 00045 // Unknown NetCDF grid 00046 return NULL; 00047 } 00048 00049 ErrorCode NCWriteHelper::collect_variable_data( std::vector< std::string >& var_names, std::vector< int >& tstep_nums ) 00050 { 00051 Interface*& mbImpl = _writeNC->mbImpl; 00052 std::vector< std::string >& dimNames = _writeNC->dimNames; 00053 std::vector< int >& dimLens = _writeNC->dimLens; 00054 std::set< std::string >& usedCoordinates = _writeNC->usedCoordinates; 00055 std::set< std::string >& dummyVarNames = _writeNC->dummyVarNames; 00056 std::map< std::string, WriteNC::VarData >& varInfo = _writeNC->varInfo; 00057 DebugOutput& dbgOut = _writeNC->dbgOut; 00058 00059 ErrorCode rval; 00060 00061 usedCoordinates.clear(); 00062 00063 if( tstep_nums.empty() && nTimeSteps > 0 ) 00064 { 00065 // No timesteps input, get them all 00066 for( int i = 0; i < nTimeSteps; i++ ) 00067 tstep_nums.push_back( i ); 00068 } 00069 00070 for( size_t i = 0; i < var_names.size(); i++ ) 00071 { 00072 std::string varname = var_names[i]; 00073 std::map< std::string, WriteNC::VarData >::iterator vit = varInfo.find( varname ); 00074 if( vit == varInfo.end() ) MB_SET_ERR( MB_FAILURE, "Can't find variable " << varname ); 00075 00076 WriteNC::VarData& currentVarData = vit->second; 00077 00078 dbgOut.tprintf( 2, " for variable %s varDims.size %d \n", varname.c_str(), 00079 (int)currentVarData.varDims.size() ); 00080 for( size_t j = 0; j < currentVarData.varDims.size(); j++ ) 00081 { 00082 std::string dimName = dimNames[currentVarData.varDims[j]]; 00083 vit = varInfo.find( dimName ); 00084 if( vit == varInfo.end() ) MB_SET_ERR( MB_FAILURE, "Can't find coordinate variable " << dimName ); 00085 00086 usedCoordinates.insert( dimName ); // Collect those used, we will need to write them to the file 00087 dbgOut.tprintf( 2, " for variable %s need dimension %s with length %d\n", varname.c_str(), 00088 dimName.c_str(), dimLens[currentVarData.varDims[j]] ); 00089 } 00090 00091 // Process coordinate variables later 00092 if( usedCoordinates.find( varname ) != usedCoordinates.end() ) continue; 00093 00094 // Default has_tsteps is false 00095 if( std::find( currentVarData.varDims.begin(), currentVarData.varDims.end(), tDim ) != 00096 currentVarData.varDims.end() ) 00097 currentVarData.has_tsteps = true; 00098 00099 // Default numLev is 0 00100 if( ( std::find( currentVarData.varDims.begin(), currentVarData.varDims.end(), levDim ) != 00101 currentVarData.varDims.end() ) ) 00102 currentVarData.numLev = nLevels; 00103 00104 // Process set variables 00105 if( WriteNC::ENTLOCSET == currentVarData.entLoc ) 00106 { 00107 if( currentVarData.has_tsteps ) 00108 { 00109 // Set variables with timesteps, e.g. xtime(Time) or xtime(Time, StrLen) 00110 // TBD 00111 MB_SET_ERR( MB_NOT_IMPLEMENTED, "Writing set variables with timesteps is not implemented yet" ); 00112 } 00113 else 00114 { 00115 // Get the tag with varname 00116 Tag tag = 0; 00117 rval = mbImpl->tag_get_handle( varname.c_str(), tag );MB_CHK_SET_ERR( rval, "Can't find tag " << varname ); 00118 currentVarData.varTags.push_back( tag ); // Really, only one for these 00119 const void* data; 00120 int size; 00121 rval = mbImpl->tag_get_by_ptr( tag, &_fileSet, 1, &data, &size );MB_CHK_SET_ERR( rval, "Can't get data of tag " << varname ); 00122 00123 // Find the type of tag, and use it 00124 DataType type; 00125 rval = mbImpl->tag_get_data_type( tag, type );MB_CHK_SET_ERR( rval, "Can't get data type of tag " << varname ); 00126 00127 currentVarData.varDataType = NC_DOUBLE; 00128 if( MB_TYPE_INTEGER == type ) currentVarData.varDataType = NC_INT; 00129 00130 assert( 0 == currentVarData.memoryHogs.size() ); // Nothing so far 00131 currentVarData.memoryHogs.push_back( (void*)data ); 00132 00133 if( currentVarData.varDims.empty() ) 00134 { 00135 // Scalar variable 00136 currentVarData.writeStarts.push_back( 0 ); 00137 currentVarData.writeCounts.push_back( 1 ); 00138 } 00139 else 00140 { 00141 for( size_t j = 0; j < currentVarData.varDims.size(); j++ ) 00142 { 00143 currentVarData.writeStarts.push_back( 0 ); 00144 currentVarData.writeCounts.push_back( dimLens[currentVarData.varDims[j]] ); 00145 } 00146 } 00147 00148 // Get variable size 00149 currentVarData.sz = 1; 00150 for( std::size_t idx = 0; idx != currentVarData.writeCounts.size(); idx++ ) 00151 currentVarData.sz *= currentVarData.writeCounts[idx]; 00152 } 00153 } // if (WriteNC::ENTLOCSET == currentVarData.entLoc) 00154 // Process non-set variables 00155 else 00156 { 00157 Tag indexedTag = 0; 00158 00159 if( currentVarData.has_tsteps ) 00160 { 00161 for( unsigned int t = 0; t < tstep_nums.size(); t++ ) 00162 { 00163 std::stringstream ssTagNameWithIndex; 00164 ssTagNameWithIndex << varname << tstep_nums[t]; 00165 rval = mbImpl->tag_get_handle( ssTagNameWithIndex.str().c_str(), indexedTag );MB_CHK_SET_ERR( rval, "Can't find tag " << ssTagNameWithIndex.str() ); 00166 dbgOut.tprintf( 2, " found indexed tag %d with name %s\n", tstep_nums[t], 00167 ssTagNameWithIndex.str().c_str() ); 00168 currentVarData.varTags.push_back( indexedTag ); 00169 } 00170 } 00171 else 00172 { 00173 // This should be a user-created non-set variable without timesteps 00174 // Treat it like having one, 0th, timestep 00175 std::stringstream ssTagNameWithIndex; 00176 ssTagNameWithIndex << varname << 0; 00177 rval = mbImpl->tag_get_handle( ssTagNameWithIndex.str().c_str(), indexedTag );MB_CHK_SET_ERR( rval, "Can't find tag " << ssTagNameWithIndex.str() << " for a user-created variable" ); 00178 dbgOut.tprintf( 2, " found indexed tag 0 with name %s\n", ssTagNameWithIndex.str().c_str() ); 00179 currentVarData.varTags.push_back( indexedTag ); 00180 } 00181 00182 // The type of the tag is fixed though 00183 DataType type; 00184 rval = mbImpl->tag_get_data_type( indexedTag, type );MB_CHK_SET_ERR( rval, "Can't get data type of tag " << varname ); 00185 00186 currentVarData.varDataType = NC_DOUBLE; 00187 if( MB_TYPE_INTEGER == type ) currentVarData.varDataType = NC_INT; 00188 } 00189 } // for (size_t i = 0; i < var_names.size(); i++) 00190 00191 // Process coordinate variables here 00192 // Check that for used coordinates we have found the tags 00193 for( std::set< std::string >::iterator setIt = usedCoordinates.begin(); setIt != usedCoordinates.end(); ++setIt ) 00194 { 00195 const std::string& coordName = *setIt; 00196 00197 std::map< std::string, WriteNC::VarData >::iterator vit = varInfo.find( coordName ); 00198 if( vit == varInfo.end() ) MB_SET_ERR( MB_FAILURE, "Can't find coordinate variable " << coordName ); 00199 00200 WriteNC::VarData& varCoordData = vit->second; 00201 Tag coordTag = 0; 00202 rval = mbImpl->tag_get_handle( coordName.c_str(), coordTag );MB_CHK_SET_ERR( rval, "Can't find tag " << coordName ); 00203 varCoordData.varTags.push_back( coordTag ); // Really, only one for these 00204 00205 const void* data; 00206 int sizeCoordinate; 00207 rval = mbImpl->tag_get_by_ptr( coordTag, &_fileSet, 1, &data, &sizeCoordinate );MB_CHK_SET_ERR( rval, "Can't get coordinate values of " << coordName ); 00208 dbgOut.tprintf( 2, " found coordinate tag with name %s and length %d\n", coordName.c_str(), sizeCoordinate ); 00209 00210 // Find the type of tag, and use it 00211 DataType type; 00212 rval = mbImpl->tag_get_data_type( coordTag, type );MB_CHK_SET_ERR( rval, "Can't get data type of tag " << coordName ); 00213 varCoordData.varDataType = NC_DOUBLE; 00214 if( MB_TYPE_INTEGER == type ) varCoordData.varDataType = NC_INT; 00215 00216 // Get dimension length (the only dimension of this coordinate variable, with the same name) 00217 assert( 1 == varCoordData.varDims.size() ); 00218 int coordDimLen = dimLens[varCoordData.varDims[0]]; 00219 00220 if( dummyVarNames.find( coordName ) != dummyVarNames.end() ) 00221 { 00222 // For a dummy coordinate variable, the tag size is always 1 00223 // The number of coordinates should be set to dimension length, instead of 1 00224 assert( 1 == sizeCoordinate ); 00225 sizeCoordinate = coordDimLen; 00226 00227 // No variable data to write 00228 data = NULL; 00229 } 00230 else 00231 { 00232 // The number of coordinates should be exactly the same as dimension length 00233 // However, if timesteps spread across files and time tag has been updated, 00234 // sizeCoordinate will be larger 00235 if( varCoordData.varDims[0] != tDim ) assert( sizeCoordinate == coordDimLen ); 00236 } 00237 00238 // For time, the actual output size and values are determined by tstep_nums 00239 if( varCoordData.varDims[0] == tDim ) 00240 { 00241 // Does not apply to dummy time tag (e.g. 'Time' tag of MPAS), when timesteps 00242 // spread across files 00243 if( NULL != data ) assert( tstep_nums.size() > 0 && tstep_nums.size() <= (size_t)sizeCoordinate ); 00244 00245 sizeCoordinate = tstep_nums.size(); 00246 00247 if( NULL != data ) 00248 { 00249 assert( NC_DOUBLE == varCoordData.varDataType ); 00250 timeStepVals.resize( sizeCoordinate ); 00251 for( unsigned int t = 0; t < tstep_nums.size(); t++ ) 00252 timeStepVals[t] = ( (double*)data )[tstep_nums[t]]; 00253 00254 data = &timeStepVals[0]; 00255 } 00256 } 00257 00258 // This is the length 00259 varCoordData.sz = sizeCoordinate; 00260 varCoordData.writeStarts.resize( 1 ); 00261 varCoordData.writeStarts[0] = 0; 00262 varCoordData.writeCounts.resize( 1 ); 00263 varCoordData.writeCounts[0] = sizeCoordinate; 00264 00265 assert( 0 == varCoordData.memoryHogs.size() ); // Nothing so far 00266 varCoordData.memoryHogs.push_back( (void*)data ); 00267 } // for (std::set<std::string>::iterator setIt ... 00268 00269 return MB_SUCCESS; 00270 } 00271 00272 ErrorCode NCWriteHelper::init_file( std::vector< std::string >& var_names, std::vector< std::string >& desired_names, 00273 bool append ) 00274 { 00275 std::vector< std::string >& dimNames = _writeNC->dimNames; 00276 std::set< std::string >& usedCoordinates = _writeNC->usedCoordinates; 00277 std::set< std::string >& dummyVarNames = _writeNC->dummyVarNames; 00278 std::map< std::string, WriteNC::VarData >& varInfo = _writeNC->varInfo; 00279 std::map< std::string, WriteNC::AttData >& globalAtts = _writeNC->globalAtts; 00280 DebugOutput& dbgOut = _writeNC->dbgOut; 00281 00282 int tDim_in_dimNames = tDim; 00283 int levDim_in_dimNames = levDim; 00284 00285 // If append mode, make sure we are in define mode; a simple open will not allow creation of new 00286 // variables 00287 if( append ) 00288 { 00289 int errcode = NCFUNC( redef )( _fileId ); 00290 if( errcode != NC_NOERR ) MB_SET_ERR( MB_FAILURE, "Can't open file in redefine mode" ); 00291 } 00292 00293 // First initialize all coordinates, then fill VarData for actual variables (and dimensions) 00294 // Check that for used coordinates we have found the tags 00295 for( std::set< std::string >::iterator setIt = usedCoordinates.begin(); setIt != usedCoordinates.end(); ++setIt ) 00296 { 00297 const std::string& coordName = *setIt; 00298 00299 std::map< std::string, WriteNC::VarData >::iterator vit = varInfo.find( coordName ); 00300 if( vit == varInfo.end() ) MB_SET_ERR( MB_FAILURE, "Can't find coordinate variable " << coordName ); 00301 00302 WriteNC::VarData& varCoordData = vit->second; 00303 varCoordData.varDims.resize( 1 ); 00304 00305 // If not append, create it for sure 00306 // If append, we might already have it, including the tag / variable with the same name 00307 /* 00308 * int ncmpi_inq_dimid(int ncid, const char *name, int *idp); 00309 */ 00310 if( append ) 00311 { 00312 int dimId; 00313 if( NCFUNC( inq_dimid )( _fileId, coordName.c_str(), &dimId ) == NC_NOERR ) 00314 { // If not found, create it later 00315 varCoordData.varDims[0] = dimId; 00316 dbgOut.tprintf( 2, " file already has coordName %s dim id is %d \n", coordName.c_str(), 00317 (int)varCoordData.varDims[0] ); 00318 00319 // Update tDim and levDim to actual dimension id 00320 if( coordName == dimNames[tDim_in_dimNames] ) 00321 tDim = varCoordData.varDims[0]; 00322 else if( coordName == dimNames[levDim_in_dimNames] ) 00323 levDim = varCoordData.varDims[0]; 00324 00325 // Skip dummy coordinate variables (e.g. ncol) 00326 if( dummyVarNames.find( coordName ) != dummyVarNames.end() ) continue; 00327 00328 // Check that the coordinate is a variable too 00329 // Inquire for a variable with the same name 00330 int varId; 00331 if( NCFUNC( inq_varid )( _fileId, coordName.c_str(), &varId ) != NC_NOERR ) 00332 MB_SET_ERR( MB_FAILURE, "We do not have a variable with the same name " << coordName ); 00333 // We should also check that this variable has one dimension, and it is dimId 00334 varCoordData.varId = varId; 00335 dbgOut.tprintf( 2, " file already has coordinate %s and varId is %d \n", coordName.c_str(), varId ); 00336 00337 continue; // Maybe more checks are needed here 00338 } 00339 } 00340 00341 /* int nc_def_dim (int ncid, const char *name, size_t len, int *dimidp); 00342 * example: status = nc_def_dim(fileId, "lat", 18L, &latid); 00343 */ 00344 00345 // Actually define a dimension 00346 if( NCFUNC( def_dim )( _fileId, coordName.c_str(), (size_t)varCoordData.sz, &varCoordData.varDims[0] ) != 00347 NC_NOERR ) 00348 MB_SET_ERR( MB_FAILURE, "Failed to generate dimension " << coordName ); 00349 dbgOut.tprintf( 2, " for coordName %s dim id is %d \n", coordName.c_str(), (int)varCoordData.varDims[0] ); 00350 00351 // Update tDim and levDim to actual dimension id 00352 if( coordName == dimNames[tDim_in_dimNames] ) 00353 tDim = varCoordData.varDims[0]; 00354 else if( coordName == dimNames[levDim_in_dimNames] ) 00355 levDim = varCoordData.varDims[0]; 00356 00357 // Create a variable with the same name, and its only dimension the one we just defined 00358 /* 00359 * int nc_def_var (int ncid, const char *name, nc_type xtype, 00360 int ndims, const int dimids[], int *varidp); 00361 example: 00362 http://www.unidata.ucar.edu/software/netcdf/docs/netcdf-c/nc_005fdef_005fvar.html#nc_005fdef_005fvar 00363 */ 00364 00365 // Skip dummy coordinate variables (e.g. ncol) 00366 if( dummyVarNames.find( coordName ) != dummyVarNames.end() ) continue; 00367 00368 // Define a coordinate variable 00369 if( NCFUNC( def_var )( _fileId, coordName.c_str(), varCoordData.varDataType, 1, &( varCoordData.varDims[0] ), 00370 &varCoordData.varId ) != NC_NOERR ) 00371 MB_SET_ERR( MB_FAILURE, "Failed to create coordinate variable " << coordName ); 00372 00373 dbgOut.tprintf( 2, " for coordName %s variable id is %d \n", coordName.c_str(), varCoordData.varId ); 00374 } 00375 00376 // Now look at requested variables, and update from the index in dimNames to the actual 00377 // dimension id 00378 for( size_t i = 0; i < var_names.size(); i++ ) 00379 { 00380 std::map< std::string, WriteNC::VarData >::iterator vit = varInfo.find( var_names[i] ); 00381 if( vit == varInfo.end() ) MB_SET_ERR( MB_FAILURE, "Can't find requested variable " << var_names[i] ); 00382 00383 // Skip coordinate variables 00384 if( usedCoordinates.find( var_names[i] ) != usedCoordinates.end() ) continue; 00385 00386 WriteNC::VarData& variableData = vit->second; 00387 00388 // The index is for dimNames; we need to find out the actual dimension id (from above) 00389 int numDims = (int)variableData.varDims.size(); 00390 for( int j = 0; j < numDims; j++ ) 00391 { 00392 std::string dimName = dimNames[variableData.varDims[j]]; 00393 std::map< std::string, WriteNC::VarData >::iterator vit2 = varInfo.find( dimName ); 00394 if( vit2 == varInfo.end() ) 00395 MB_SET_ERR( MB_FAILURE, "Can't find requested coordinate variable " << dimName ); 00396 00397 WriteNC::VarData& coordData = vit2->second; 00398 // Index in dimNames to actual dimension id 00399 variableData.varDims[j] = coordData.varDims[0]; // This one, being a coordinate, is the only one 00400 dbgOut.tprintf( 2, " dimension with index %d name %s has ID %d \n", j, dimName.c_str(), 00401 variableData.varDims[j] ); 00402 } 00403 00404 // Define the variable now: 00405 int errCode = 00406 NCFUNC( def_var )( _fileId, desired_names[i].c_str(), variableData.varDataType, 00407 (int)variableData.varDims.size(), &( variableData.varDims[0] ), &variableData.varId ); 00408 if( errCode != NC_NOERR ) MB_SET_ERR( MB_FAILURE, "Failed to create requested variable " << desired_names[i] ); 00409 00410 dbgOut.tprintf( 2, " for variable %s with desired name %s variable id is %d \n", var_names[i].c_str(), 00411 desired_names[i].c_str(), variableData.varId ); 00412 // Now define the variable, with all dimensions 00413 } 00414 00415 // Define global attributes (exactly copied from the original file for the time being) 00416 // Should we modify some of them (e.g. revision_Id) later? 00417 std::map< std::string, WriteNC::AttData >::iterator attIt; 00418 for( attIt = globalAtts.begin(); attIt != globalAtts.end(); ++attIt ) 00419 { 00420 const std::string& attName = attIt->first; 00421 WriteNC::AttData& attData = attIt->second; 00422 NCDF_SIZE& attLen = attData.attLen; 00423 nc_type& attDataType = attData.attDataType; 00424 const std::string& attValue = attData.attValue; 00425 00426 switch( attDataType ) 00427 { 00428 case NC_BYTE: 00429 case NC_CHAR: 00430 if( NC_NOERR != 00431 NCFUNC( put_att_text )( _fileId, NC_GLOBAL, attName.c_str(), attLen, attValue.c_str() ) ) 00432 MB_SET_ERR( MB_FAILURE, "Failed to define text type attribute" ); 00433 break; 00434 case NC_DOUBLE: 00435 if( NC_NOERR != NCFUNC( put_att_double )( _fileId, NC_GLOBAL, attName.c_str(), NC_DOUBLE, 1, 00436 (double*)attValue.c_str() ) ) 00437 MB_SET_ERR( MB_FAILURE, "Failed to define double type attribute" ); 00438 break; 00439 case NC_FLOAT: 00440 if( NC_NOERR != NCFUNC( put_att_float )( _fileId, NC_GLOBAL, attName.c_str(), NC_FLOAT, 1, 00441 (float*)attValue.c_str() ) ) 00442 MB_SET_ERR( MB_FAILURE, "Failed to define float type attribute" ); 00443 break; 00444 case NC_INT: 00445 if( NC_NOERR != 00446 NCFUNC( put_att_int )( _fileId, NC_GLOBAL, attName.c_str(), NC_INT, 1, (int*)attValue.c_str() ) ) 00447 MB_SET_ERR( MB_FAILURE, "Failed to define int type attribute" ); 00448 break; 00449 case NC_SHORT: 00450 if( NC_NOERR != NCFUNC( put_att_short )( _fileId, NC_GLOBAL, attName.c_str(), NC_SHORT, 1, 00451 (short*)attValue.c_str() ) ) 00452 MB_SET_ERR( MB_FAILURE, "Failed to define short type attribute" ); 00453 break; 00454 default: 00455 MB_SET_ERR( MB_FAILURE, "Unknown attribute data type" ); 00456 } 00457 } 00458 00459 // Take it out of define mode 00460 if( NC_NOERR != NCFUNC( enddef )( _fileId ) ) MB_SET_ERR( MB_FAILURE, "Failed to close define mode" ); 00461 00462 return MB_SUCCESS; 00463 } 00464 00465 ErrorCode NCWriteHelper::write_values( std::vector< std::string >& var_names, std::vector< int >& tstep_nums ) 00466 { 00467 std::set< std::string >& usedCoordinates = _writeNC->usedCoordinates; 00468 std::set< std::string >& dummyVarNames = _writeNC->dummyVarNames; 00469 std::map< std::string, WriteNC::VarData >& varInfo = _writeNC->varInfo; 00470 00471 std::vector< WriteNC::VarData > vdatas; 00472 std::vector< WriteNC::VarData > vsetdatas; 00473 00474 // For set variables, include coordinates used by requested var_names 00475 for( std::set< std::string >::iterator setIt = usedCoordinates.begin(); setIt != usedCoordinates.end(); ++setIt ) 00476 { 00477 const std::string& coordName = *setIt; 00478 00479 // Skip dummy coordinate variables (if any) 00480 if( dummyVarNames.find( coordName ) != dummyVarNames.end() ) continue; 00481 00482 std::map< std::string, WriteNC::VarData >::iterator vit = varInfo.find( coordName ); 00483 if( vit == varInfo.end() ) { MB_SET_ERR( MB_FAILURE, "Can't find coordinate variable " << coordName ); } 00484 00485 vsetdatas.push_back( vit->second ); 00486 } 00487 00488 // Collect non-set and set variables from requested var_names 00489 for( unsigned int i = 0; i < var_names.size(); i++ ) 00490 { 00491 std::map< std::string, WriteNC::VarData >::iterator vit = varInfo.find( var_names[i] ); 00492 if( vit == varInfo.end() ) { MB_SET_ERR( MB_FAILURE, "Can't find requested variable " << var_names[i] ); } 00493 00494 WriteNC::VarData& variableData = vit->second; 00495 if( WriteNC::ENTLOCSET == variableData.entLoc ) 00496 { 00497 // Used coordinates has all ready been included 00498 if( usedCoordinates.find( var_names[i] ) != usedCoordinates.end() ) continue; 00499 00500 vsetdatas.push_back( variableData ); 00501 } 00502 else 00503 vdatas.push_back( variableData ); 00504 } 00505 00506 // Assume that the data ranges do not overlap across processors 00507 // While overlapped writing might still work, we should better not take that risk 00508 write_nonset_variables( vdatas, tstep_nums ); 00509 00510 // Use independent I/O mode put, since this write is only for the root processor 00511 write_set_variables( vsetdatas, tstep_nums ); 00512 00513 return MB_SUCCESS; 00514 } 00515 00516 ErrorCode NCWriteHelper::write_set_variables( std::vector< WriteNC::VarData >& vsetdatas, 00517 std::vector< int >& /* tstep_nums */ ) 00518 { 00519 int success; 00520 00521 // CAUTION: if the NetCDF ID is from a previous call to ncmpi_create rather than ncmpi_open, 00522 // all processors need to call ncmpi_begin_indep_data(). If only the root processor does so, 00523 // ncmpi_begin_indep_data() call will be blocked forever :( 00524 #ifdef MOAB_HAVE_PNETCDF 00525 // Enter independent I/O mode 00526 success = NCFUNC( begin_indep_data )( _fileId ); 00527 if( success ) MB_SET_ERR( MB_FAILURE, "Failed to begin independent I/O mode" ); 00528 #endif 00529 00530 int rank = 0; 00531 #ifdef MOAB_HAVE_MPI 00532 bool& isParallel = _writeNC->isParallel; 00533 if( isParallel ) 00534 { 00535 ParallelComm*& myPcomm = _writeNC->myPcomm; 00536 rank = myPcomm->proc_config().proc_rank(); 00537 } 00538 #endif 00539 if( 0 == rank ) 00540 { 00541 for( unsigned int i = 0; i < vsetdatas.size(); i++ ) 00542 { 00543 WriteNC::VarData& variableData = vsetdatas[i]; 00544 00545 // Set variables with timesteps, e.g. xtime(Time) or xtime(Time, StrLen) 00546 if( variableData.has_tsteps ) 00547 { 00548 MB_SET_ERR( MB_NOT_IMPLEMENTED, "Writing set variables with timesteps is not implemented yet" ); 00549 } 00550 00551 switch( variableData.varDataType ) 00552 { 00553 case NC_DOUBLE: 00554 // Independent I/O mode put 00555 success = NCFUNCP( _vara_double )( _fileId, variableData.varId, &variableData.writeStarts[0], 00556 &variableData.writeCounts[0], 00557 (double*)( variableData.memoryHogs[0] ) ); 00558 if( success ) 00559 MB_SET_ERR( MB_FAILURE, "Failed to write double data for variable " << variableData.varName ); 00560 break; 00561 case NC_INT: 00562 // Independent I/O mode put 00563 success = 00564 NCFUNCP( _vara_int )( _fileId, variableData.varId, &variableData.writeStarts[0], 00565 &variableData.writeCounts[0], (int*)( variableData.memoryHogs[0] ) ); 00566 if( success ) 00567 MB_SET_ERR( MB_FAILURE, "Failed to write int data for variable " << variableData.varName ); 00568 break; 00569 default: 00570 MB_SET_ERR( MB_NOT_IMPLEMENTED, "Writing non-double or non-int data is not implemented yet" ); 00571 } 00572 } 00573 } 00574 00575 #ifdef MOAB_HAVE_PNETCDF 00576 // End independent I/O mode 00577 success = NCFUNC( end_indep_data )( _fileId ); 00578 if( success ) MB_SET_ERR( MB_FAILURE, "Failed to end independent I/O mode" ); 00579 #endif 00580 00581 return MB_SUCCESS; 00582 } 00583 00584 ErrorCode ScdNCWriteHelper::collect_mesh_info() 00585 { 00586 Interface*& mbImpl = _writeNC->mbImpl; 00587 std::vector< std::string >& dimNames = _writeNC->dimNames; 00588 std::vector< int >& dimLens = _writeNC->dimLens; 00589 00590 ErrorCode rval; 00591 00592 // Look for time dimension 00593 std::vector< std::string >::iterator vecIt; 00594 if( ( vecIt = std::find( dimNames.begin(), dimNames.end(), "time" ) ) != dimNames.end() ) 00595 tDim = vecIt - dimNames.begin(); 00596 else if( ( vecIt = std::find( dimNames.begin(), dimNames.end(), "t" ) ) != dimNames.end() ) 00597 tDim = vecIt - dimNames.begin(); 00598 else 00599 { 00600 MB_SET_ERR( MB_FAILURE, "Couldn't find 'time' or 't' dimension" ); 00601 } 00602 nTimeSteps = dimLens[tDim]; 00603 00604 // Get number of levels 00605 if( ( vecIt = std::find( dimNames.begin(), dimNames.end(), "lev" ) ) != dimNames.end() ) 00606 levDim = vecIt - dimNames.begin(); 00607 else if( ( vecIt = std::find( dimNames.begin(), dimNames.end(), "ilev" ) ) != dimNames.end() ) 00608 levDim = vecIt - dimNames.begin(); 00609 else 00610 { 00611 MB_SET_ERR( MB_FAILURE, "Couldn't find 'lev' or 'ilev' dimension" ); 00612 } 00613 nLevels = dimLens[levDim]; 00614 00615 // __<dim_name>_LOC_MINMAX (for slon, slat, lon and lat) 00616 Tag convTag = 0; 00617 rval = mbImpl->tag_get_handle( "__slon_LOC_MINMAX", 0, MB_TYPE_INTEGER, convTag, MB_TAG_ANY );MB_CHK_SET_ERR( rval, "Trouble getting conventional tag __slon_LOC_MINMAX" ); 00618 int val[2]; 00619 rval = mbImpl->tag_get_data( convTag, &_fileSet, 1, val );MB_CHK_SET_ERR( rval, "Trouble getting data of conventional tag __slon_LOC_MINMAX" ); 00620 lDims[0] = val[0]; 00621 lDims[3] = val[1]; 00622 00623 rval = mbImpl->tag_get_handle( "__slat_LOC_MINMAX", 0, MB_TYPE_INTEGER, convTag, MB_TAG_ANY );MB_CHK_SET_ERR( rval, "Trouble getting conventional tag __slat_LOC_MINMAX" ); 00624 rval = mbImpl->tag_get_data( convTag, &_fileSet, 1, val );MB_CHK_SET_ERR( rval, "Trouble getting data of conventional tag __slat_LOC_MINMAX" ); 00625 lDims[1] = val[0]; 00626 lDims[4] = val[1]; 00627 00628 rval = mbImpl->tag_get_handle( "__lon_LOC_MINMAX", 0, MB_TYPE_INTEGER, convTag, MB_TAG_ANY );MB_CHK_SET_ERR( rval, "Trouble getting conventional tag __lon_LOC_MINMAX" ); 00629 rval = mbImpl->tag_get_data( convTag, &_fileSet, 1, val );MB_CHK_SET_ERR( rval, "Trouble getting data of conventional tag __lon_LOC_MINMAX" ); 00630 lCDims[0] = val[0]; 00631 lCDims[3] = val[1]; 00632 00633 rval = mbImpl->tag_get_handle( "__lat_LOC_MINMAX", 0, MB_TYPE_INTEGER, convTag, MB_TAG_ANY );MB_CHK_SET_ERR( rval, "Trouble getting conventional tag __lat_LOC_MINMAX" ); 00634 rval = mbImpl->tag_get_data( convTag, &_fileSet, 1, val );MB_CHK_SET_ERR( rval, "Trouble getting data of conventional tag __lat_LOC_MINMAX" ); 00635 lCDims[1] = val[0]; 00636 lCDims[4] = val[1]; 00637 00638 // Get local faces 00639 rval = mbImpl->get_entities_by_dimension( _fileSet, 2, localCellsOwned );MB_CHK_SET_ERR( rval, "Trouble getting local faces in current file set" ); 00640 assert( !localCellsOwned.empty() ); 00641 00642 #ifdef MOAB_HAVE_MPI 00643 bool& isParallel = _writeNC->isParallel; 00644 if( isParallel ) 00645 { 00646 ParallelComm*& myPcomm = _writeNC->myPcomm; 00647 int procs = myPcomm->proc_config().proc_size(); 00648 if( procs > 1 ) 00649 { 00650 rval = myPcomm->filter_pstatus( localCellsOwned, PSTATUS_NOT_OWNED, PSTATUS_NOT );MB_CHK_SET_ERR( rval, "Trouble getting owned faces in current file set" ); 00651 } 00652 } 00653 #endif 00654 00655 return MB_SUCCESS; 00656 } 00657 00658 ErrorCode ScdNCWriteHelper::collect_variable_data( std::vector< std::string >& var_names, 00659 std::vector< int >& tstep_nums ) 00660 { 00661 NCWriteHelper::collect_variable_data( var_names, tstep_nums ); 00662 00663 std::map< std::string, WriteNC::VarData >& varInfo = _writeNC->varInfo; 00664 00665 for( size_t i = 0; i < var_names.size(); i++ ) 00666 { 00667 std::string varname = var_names[i]; 00668 std::map< std::string, WriteNC::VarData >::iterator vit = varInfo.find( varname ); 00669 if( vit == varInfo.end() ) MB_SET_ERR( MB_FAILURE, "Can't find variable " << varname ); 00670 00671 WriteNC::VarData& currentVarData = vit->second; 00672 #ifndef NDEBUG 00673 std::vector< int >& varDims = currentVarData.varDims; 00674 #endif 00675 00676 // Skip set variables, which were already processed in 00677 // NCWriteHelper::collect_variable_data() 00678 if( WriteNC::ENTLOCSET == currentVarData.entLoc ) continue; 00679 00680 // Set up writeStarts and writeCounts (maximum number of dimensions is 4) 00681 currentVarData.writeStarts.resize( 4 ); 00682 currentVarData.writeCounts.resize( 4 ); 00683 unsigned int dim_idx = 0; 00684 00685 // First: time 00686 if( currentVarData.has_tsteps ) 00687 { 00688 // Non-set variables with timesteps 00689 // 4 dimensions like (time, lev, lat, lon) 00690 // 3 dimensions like (time, lat, lon) 00691 assert( 4 == varDims.size() || 3 == varDims.size() ); 00692 00693 // Time should be the first dimension 00694 assert( tDim == varDims[0] ); 00695 00696 currentVarData.writeStarts[dim_idx] = 0; // This value is timestep dependent, will be set later 00697 currentVarData.writeCounts[dim_idx] = 1; 00698 dim_idx++; 00699 } 00700 else 00701 { 00702 // Non-set variables without timesteps 00703 // 3 dimensions like (lev, lat, lon) 00704 // 2 dimensions like (lat, lon) 00705 assert( 3 == varDims.size() || 2 == varDims.size() ); 00706 } 00707 00708 // Next: lev 00709 if( currentVarData.numLev > 0 ) 00710 { 00711 // Non-set variables with levels 00712 // 4 dimensions like (time, lev, lat, lon) 00713 // 3 dimensions like (lev, lat, lon) 00714 assert( 4 == varDims.size() || 3 == varDims.size() ); 00715 00716 currentVarData.writeStarts[dim_idx] = 0; 00717 currentVarData.writeCounts[dim_idx] = currentVarData.numLev; 00718 dim_idx++; 00719 } 00720 else 00721 { 00722 // Non-set variables without levels 00723 // 3 dimensions like (time, lat, lon) 00724 // 2 dimensions like (lat, lon) 00725 assert( 3 == varDims.size() || 2 == varDims.size() ); 00726 } 00727 00728 // Finally: lat and lon 00729 switch( currentVarData.entLoc ) 00730 { 00731 case WriteNC::ENTLOCFACE: 00732 // Faces 00733 currentVarData.writeStarts[dim_idx] = lCDims[1]; 00734 currentVarData.writeCounts[dim_idx] = lCDims[4] - lCDims[1] + 1; 00735 currentVarData.writeStarts[dim_idx + 1] = lCDims[0]; 00736 currentVarData.writeCounts[dim_idx + 1] = lCDims[3] - lCDims[0] + 1; 00737 break; 00738 default: 00739 MB_SET_ERR( MB_FAILURE, "Unexpected entity location type for variable " << varname ); 00740 } 00741 dim_idx += 2; 00742 00743 // Get variable size 00744 currentVarData.sz = 1; 00745 for( std::size_t idx = 0; idx < dim_idx; idx++ ) 00746 currentVarData.sz *= currentVarData.writeCounts[idx]; 00747 } // for (size_t i = 0; i < var_names.size(); i++) 00748 00749 return MB_SUCCESS; 00750 } 00751 00752 // Write CAM-EUL and CAM-FV non-set variables on non-shared quads (e.g. T) 00753 // We assume that there are no variables on vertices and we do not support 00754 // variables on edges (e.g. US in CAM-FV) for the time being 00755 ErrorCode ScdNCWriteHelper::write_nonset_variables( std::vector< WriteNC::VarData >& vdatas, 00756 std::vector< int >& tstep_nums ) 00757 { 00758 Interface*& mbImpl = _writeNC->mbImpl; 00759 00760 int success; 00761 00762 // For each indexed variable tag, write a time step data 00763 for( unsigned int i = 0; i < vdatas.size(); i++ ) 00764 { 00765 WriteNC::VarData& variableData = vdatas[i]; 00766 00767 // Assume this variable is on faces for the time being 00768 switch( variableData.entLoc ) 00769 { 00770 case WriteNC::ENTLOCFACE: 00771 // Faces 00772 break; 00773 default: 00774 MB_SET_ERR( MB_FAILURE, "Unexpected entity location type for variable " << variableData.varName ); 00775 } 00776 00777 unsigned int num_timesteps; 00778 unsigned int lat_idx = 0; 00779 unsigned int lon_idx = 1; 00780 if( variableData.has_tsteps ) 00781 { 00782 // Non-set variables with timesteps 00783 // 4 dimensions like (time, lev, lat, lon) 00784 // 3 dimensions like (time, lat, lon) 00785 num_timesteps = tstep_nums.size(); 00786 lat_idx++; 00787 lon_idx++; 00788 } 00789 else 00790 { 00791 // Non-set variables without timesteps 00792 // 3 dimensions like (lev, lat, lon) 00793 // 2 dimensions like (lat, lon) 00794 num_timesteps = 1; 00795 } 00796 00797 unsigned int num_lev; 00798 if( variableData.numLev > 0 ) 00799 { 00800 // Non-set variables with levels 00801 // 4 dimensions like (time, lev, lat, lon) 00802 // 3 dimensions like (lev, lat, lon) 00803 num_lev = variableData.numLev; 00804 lat_idx++; 00805 lon_idx++; 00806 } 00807 else 00808 { 00809 // Non-set variables without levels 00810 // 3 dimensions like (time, lat, lon) 00811 // 2 dimensions like (lat, lon) 00812 num_lev = 1; 00813 } 00814 00815 size_t ni = variableData.writeCounts[lon_idx]; // lon 00816 size_t nj = variableData.writeCounts[lat_idx]; // lat 00817 00818 // At each timestep, we need to transpose tag format (lat, lon, lev) back 00819 // to NC format (lev, lat, lon) for writing 00820 for( unsigned int t = 0; t < num_timesteps; t++ ) 00821 { 00822 // We will write one time step, and count will be one; start will be different 00823 // Use tag_iterate to get tag data (assume that localCellsOwned is contiguous) 00824 // We should also transpose for level so that means deep copy for transpose 00825 if( tDim == variableData.varDims[0] ) variableData.writeStarts[0] = t; // This is start for time 00826 int count; 00827 void* dataptr; 00828 ErrorCode rval = mbImpl->tag_iterate( variableData.varTags[t], localCellsOwned.begin(), 00829 localCellsOwned.end(), count, dataptr );MB_CHK_SET_ERR( rval, "Failed to iterate tag on owned faces" ); 00830 assert( count == (int)localCellsOwned.size() ); 00831 00832 // Now transpose and write tag data 00833 // Use collective I/O mode put (synchronous write) for the time being, we can try 00834 // nonblocking put (request aggregation) later 00835 switch( variableData.varDataType ) 00836 { 00837 case NC_DOUBLE: { 00838 std::vector< double > tmpdoubledata( ni * nj * num_lev ); 00839 if( num_lev > 1 ) 00840 // Transpose (lat, lon, lev) back to (lev, lat, lon) 00841 jik_to_kji( ni, nj, num_lev, &tmpdoubledata[0], (double*)( dataptr ) ); 00842 success = NCFUNCAP( _vara_double )( _fileId, variableData.varId, &variableData.writeStarts[0], 00843 &variableData.writeCounts[0], &tmpdoubledata[0] ); 00844 if( success ) 00845 MB_SET_ERR( MB_FAILURE, "Failed to write double data for variable " << variableData.varName ); 00846 break; 00847 } 00848 default: 00849 MB_SET_ERR( MB_NOT_IMPLEMENTED, "Writing non-double data is not implemented yet" ); 00850 } 00851 } 00852 } 00853 00854 return MB_SUCCESS; 00855 } 00856 00857 } /* namespace moab */