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