Mesh Oriented datABase  (version 5.4.1)
Array-based unstructured mesh datastructure
NCWriteHOMME.cpp
Go to the documentation of this file.
00001 /*
00002  * NCWriteHOMME.cpp
00003  *
00004  *  Created on: April 9, 2014
00005  */
00006 
00007 #include "NCWriteHOMME.hpp"
00008 #include "MBTagConventions.hpp"
00009 
00010 namespace moab
00011 {
00012 
00013 NCWriteHOMME::~NCWriteHOMME()
00014 {
00015     // TODO Auto-generated destructor stub
00016 }
00017 
00018 ErrorCode NCWriteHOMME::collect_mesh_info()
00019 {
00020     Interface*& mbImpl                   = _writeNC->mbImpl;
00021     std::vector< std::string >& dimNames = _writeNC->dimNames;
00022     std::vector< int >& dimLens          = _writeNC->dimLens;
00023     Tag& mGlobalIdTag                    = _writeNC->mGlobalIdTag;
00024 
00025     ErrorCode rval;
00026 
00027     // Look for time dimension
00028     std::vector< std::string >::iterator vecIt;
00029     if( ( vecIt = std::find( dimNames.begin(), dimNames.end(), "time" ) ) != dimNames.end() )
00030         tDim = vecIt - dimNames.begin();
00031     else
00032     {
00033         MB_SET_ERR( MB_FAILURE, "Couldn't find 'time' dimension" );
00034     }
00035     nTimeSteps = dimLens[tDim];
00036 
00037     // Get number of levels
00038     if( ( vecIt = std::find( dimNames.begin(), dimNames.end(), "lev" ) ) != dimNames.end() )
00039         levDim = vecIt - dimNames.begin();
00040     else
00041     {
00042         MB_SET_ERR( MB_FAILURE, "Couldn't find 'lev' dimension" );
00043     }
00044     nLevels = dimLens[levDim];
00045 
00046     // Get local vertices
00047     rval = mbImpl->get_entities_by_dimension( _fileSet, 0, localVertsOwned );MB_CHK_SET_ERR( rval, "Trouble getting local vertices in current file set" );
00048     assert( !localVertsOwned.empty() );
00049 
00050 #ifdef MOAB_HAVE_MPI
00051     bool& isParallel = _writeNC->isParallel;
00052     if( isParallel )
00053     {
00054         ParallelComm*& myPcomm = _writeNC->myPcomm;
00055         int rank               = myPcomm->proc_config().proc_rank();
00056         int procs              = myPcomm->proc_config().proc_size();
00057         if( procs > 1 )
00058         {
00059 #ifndef NDEBUG
00060             unsigned int num_local_verts = localVertsOwned.size();
00061 #endif
00062             rval = myPcomm->filter_pstatus( localVertsOwned, PSTATUS_NOT_OWNED, PSTATUS_NOT );MB_CHK_SET_ERR( rval, "Trouble getting owned vertices in current set" );
00063 
00064             // Assume that PARALLEL_RESOLVE_SHARED_ENTS option is set
00065             // Verify that not all local vertices are owned by the last processor
00066             if( procs - 1 == rank )
00067                 assert( "PARALLEL_RESOLVE_SHARED_ENTS option is set" && localVertsOwned.size() < num_local_verts );
00068         }
00069     }
00070 #endif
00071 
00072     std::vector< int > gids( localVertsOwned.size() );
00073     rval = mbImpl->tag_get_data( mGlobalIdTag, localVertsOwned, &gids[0] );MB_CHK_SET_ERR( rval, "Trouble getting global IDs on local vertices" );
00074 
00075     // Get localGidVertsOwned
00076     std::copy( gids.rbegin(), gids.rend(), range_inserter( localGidVertsOwned ) );
00077 
00078     return MB_SUCCESS;
00079 }
00080 
00081 ErrorCode NCWriteHOMME::collect_variable_data( std::vector< std::string >& var_names, std::vector< int >& tstep_nums )
00082 {
00083     NCWriteHelper::collect_variable_data( var_names, tstep_nums );
00084 
00085     std::map< std::string, WriteNC::VarData >& varInfo = _writeNC->varInfo;
00086 
00087     for( size_t i = 0; i < var_names.size(); i++ )
00088     {
00089         std::string varname                                     = var_names[i];
00090         std::map< std::string, WriteNC::VarData >::iterator vit = varInfo.find( varname );
00091         if( vit == varInfo.end() ) MB_SET_ERR( MB_FAILURE, "Can't find variable " << varname );
00092         ;
00093 
00094         WriteNC::VarData& currentVarData = vit->second;
00095 #ifndef NDEBUG
00096         std::vector< int >& varDims = currentVarData.varDims;
00097 #endif
00098 
00099         // Skip set variables, which were already processed in
00100         // NCWriteHelper::collect_variable_data()
00101         if( WriteNC::ENTLOCSET == currentVarData.entLoc ) continue;
00102 
00103         // Set up writeStarts and writeCounts (maximum number of dimensions is 3)
00104         currentVarData.writeStarts.resize( 3 );
00105         currentVarData.writeCounts.resize( 3 );
00106         unsigned int dim_idx = 0;
00107 
00108         // First: time
00109         if( currentVarData.has_tsteps )
00110         {
00111             // Non-set variables with timesteps
00112             // 3 dimensions like (time, lev, ncol)
00113             // 2 dimensions like (time, ncol)
00114             assert( 3 == varDims.size() || 2 == varDims.size() );
00115 
00116             // Time should be the first dimension
00117             assert( tDim == varDims[0] );
00118 
00119             currentVarData.writeStarts[dim_idx] = 0;  // This value is timestep dependent, will be set later
00120             currentVarData.writeCounts[dim_idx] = 1;
00121             dim_idx++;
00122         }
00123         else
00124         {
00125             // Non-set variables without timesteps
00126             // 2 dimensions like (lev, ncol)
00127             // 1 dimension like (ncol)
00128             assert( 2 == varDims.size() || 1 == varDims.size() );
00129         }
00130 
00131         // Next: lev
00132         if( currentVarData.numLev > 0 )
00133         {
00134             // Non-set variables with levels
00135             // 3 dimensions like (time, lev, ncol)
00136             // 2 dimensions like (lev, ncol)
00137             assert( 3 == varDims.size() || 2 == varDims.size() );
00138 
00139             currentVarData.writeStarts[dim_idx] = 0;
00140             currentVarData.writeCounts[dim_idx] = currentVarData.numLev;
00141             dim_idx++;
00142         }
00143         else
00144         {
00145             // Non-set variables without levels
00146             // 2 dimensions like (time, ncol)
00147             // 1 dimension like (ncol)
00148             assert( 2 == varDims.size() || 1 == varDims.size() );
00149         }
00150 
00151         // Finally: ncol
00152         switch( currentVarData.entLoc )
00153         {
00154             case WriteNC::ENTLOCVERT:
00155                 // Vertices
00156                 // Start from the first localGidVerts
00157                 // Actually, this will be reset later for writing
00158                 currentVarData.writeStarts[dim_idx] = localGidVertsOwned[0] - 1;
00159                 currentVarData.writeCounts[dim_idx] = localGidVertsOwned.size();
00160                 break;
00161             default:
00162                 MB_SET_ERR( MB_FAILURE, "Unexpected entity location type for variable " << varname );
00163         }
00164         dim_idx++;
00165 
00166         // Get variable size
00167         currentVarData.sz = 1;
00168         for( std::size_t idx = 0; idx < dim_idx; idx++ )
00169             currentVarData.sz *= currentVarData.writeCounts[idx];
00170     }
00171 
00172     return MB_SUCCESS;
00173 }
00174 
00175 ErrorCode NCWriteHOMME::write_nonset_variables( std::vector< WriteNC::VarData >& vdatas,
00176                                                 std::vector< int >& tstep_nums )
00177 {
00178     Interface*& mbImpl = _writeNC->mbImpl;
00179 
00180     int success;
00181     int num_local_verts_owned = localVertsOwned.size();
00182 
00183     // For each indexed variable tag, write a time step data
00184     for( unsigned int i = 0; i < vdatas.size(); i++ )
00185     {
00186         WriteNC::VarData& variableData = vdatas[i];
00187 
00188         // Assume this variable is on vertices for the time being
00189         switch( variableData.entLoc )
00190         {
00191             case WriteNC::ENTLOCVERT:
00192                 // Vertices
00193                 break;
00194             default:
00195                 MB_SET_ERR( MB_FAILURE, "Unexpected entity location type for variable " << variableData.varName );
00196         }
00197 
00198         unsigned int num_timesteps;
00199         unsigned int ncol_idx = 0;
00200         if( variableData.has_tsteps )
00201         {
00202             // Non-set variables with timesteps
00203             // 3 dimensions like (time, lev, ncol)
00204             // 2 dimensions like (time, ncol)
00205             num_timesteps = tstep_nums.size();
00206             ncol_idx++;
00207         }
00208         else
00209         {
00210             // Non-set variables without timesteps
00211             // 2 dimensions like (lev, ncol)
00212             // 1 dimension like (ncol)
00213             num_timesteps = 1;
00214         }
00215 
00216         unsigned int num_lev;
00217         if( variableData.numLev > 0 )
00218         {
00219             // Non-set variables with levels
00220             // 3 dimensions like (time, lev, ncol)
00221             // 2 dimensions like (lev, ncol)
00222             num_lev = variableData.numLev;
00223             ncol_idx++;
00224         }
00225         else
00226         {
00227             // Non-set variables without levels
00228             // 2 dimensions like (time, ncol)
00229             // 1 dimension like (ncol)
00230             num_lev = 1;
00231         }
00232 
00233         // At each timestep, we need to transpose tag format (ncol, lev) back
00234         // to NC format (lev, ncol) for writing
00235         for( unsigned int t = 0; t < num_timesteps; t++ )
00236         {
00237             // We will write one time step, and count will be one; start will be different
00238             // Use tag_get_data instead of tag_iterate to copy tag data, as localVertsOwned
00239             // might not be contiguous. We should also transpose for level so that means
00240             // deep copy for transpose
00241             if( tDim == variableData.varDims[0] ) variableData.writeStarts[0] = t;  // This is start for time
00242             std::vector< double > tag_data( num_local_verts_owned * num_lev );
00243             ErrorCode rval = mbImpl->tag_get_data( variableData.varTags[t], localVertsOwned, &tag_data[0] );MB_CHK_SET_ERR( rval, "Trouble getting tag data on owned vertices" );
00244 
00245 #ifdef MOAB_HAVE_PNETCDF
00246             size_t nb_writes = localGidVertsOwned.psize();
00247             std::vector< int > requests( nb_writes ), statuss( nb_writes );
00248             size_t idxReq = 0;
00249 #endif
00250 
00251             // Now transpose and write copied tag data
00252             // Use nonblocking put (request aggregation)
00253             switch( variableData.varDataType )
00254             {
00255                 case NC_DOUBLE: {
00256                     std::vector< double > tmpdoubledata( num_local_verts_owned * num_lev );
00257                     if( num_lev > 1 )
00258                     {
00259                         // Transpose (ncol, lev) back to (lev, ncol)
00260                         // Note, num_local_verts_owned is not used by jik_to_kji_stride()
00261                         jik_to_kji_stride( num_local_verts_owned, 1, num_lev, &tmpdoubledata[0], &tag_data[0],
00262                                            localGidVertsOwned );
00263                     }
00264 
00265                     size_t indexInDoubleArray = 0;
00266                     size_t ic                 = 0;
00267                     for( Range::pair_iterator pair_iter = localGidVertsOwned.pair_begin();
00268                          pair_iter != localGidVertsOwned.pair_end(); ++pair_iter, ic++ )
00269                     {
00270                         EntityHandle starth                = pair_iter->first;
00271                         EntityHandle endh                  = pair_iter->second;
00272                         variableData.writeStarts[ncol_idx] = (NCDF_SIZE)( starth - 1 );
00273                         variableData.writeCounts[ncol_idx] = (NCDF_SIZE)( endh - starth + 1 );
00274 
00275                         // Do a partial write, in each subrange
00276 #ifdef MOAB_HAVE_PNETCDF
00277                         // Wait outside this loop
00278                         success =
00279                             NCFUNCREQP( _vara_double )( _fileId, variableData.varId, &( variableData.writeStarts[0] ),
00280                                                         &( variableData.writeCounts[0] ),
00281                                                         &( tmpdoubledata[indexInDoubleArray] ), &requests[idxReq++] );
00282 #else
00283                         success = NCFUNCAP(
00284                             _vara_double )( _fileId, variableData.varId, &( variableData.writeStarts[0] ),
00285                                             &( variableData.writeCounts[0] ), &( tmpdoubledata[indexInDoubleArray] ) );
00286 #endif
00287                         if( success )
00288                             MB_SET_ERR( MB_FAILURE,
00289                                         "Failed to write double data in a loop for variable " << variableData.varName );
00290                         // We need to increment the index in double array for the
00291                         // next subrange
00292                         indexInDoubleArray += ( endh - starth + 1 ) * num_lev;
00293                     }
00294                     assert( ic == localGidVertsOwned.psize() );
00295 #ifdef MOAB_HAVE_PNETCDF
00296                     success = ncmpi_wait_all( _fileId, requests.size(), &requests[0], &statuss[0] );
00297                     if( success ) MB_SET_ERR( MB_FAILURE, "Failed on wait_all" );
00298 #endif
00299                     break;
00300                 }
00301                 default:
00302                     MB_SET_ERR( MB_NOT_IMPLEMENTED, "Writing non-double data is not implemented yet" );
00303             }
00304         }
00305     }
00306 
00307     return MB_SUCCESS;
00308 }
00309 
00310 } /* namespace moab */
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines