MOAB: Mesh Oriented datABase  (version 5.3.0)
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, 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 */
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines