moab::NCHelperHOMME Class Reference

Child helper class for HOMME grid (CAM_SE) More...

#include <NCHelperHOMME.hpp>

Inheritance diagram for moab::NCHelperHOMME:

Collaboration diagram for moab::NCHelperHOMME:

Public Member Functions
	NCHelperHOMME (ReadNC *readNC, int fileId, const FileOptions &opts, EntityHandle fileSet)
Public Member Functions inherited from moab::UcdNCHelper
	UcdNCHelper (ReadNC *readNC, int fileId, const FileOptions &opts, EntityHandle fileSet)
virtual	~UcdNCHelper ()
Public Member Functions inherited from moab::NCHelper
	NCHelper (ReadNC *readNC, int fileId, const FileOptions &opts, EntityHandle fileSet)
virtual	~NCHelper ()
ErrorCode	create_conventional_tags (const std::vector< int > &tstep_nums)
	Create NC conventional tags. More...
ErrorCode	update_time_tag_vals ()
	Update time tag values if timesteps spread across files. More...

Static Public Member Functions
static bool	can_read_file (ReadNC *readNC, int fileId)
Static Public Member Functions inherited from moab::NCHelper
static ReadNC::NCFormatType	get_nc_format (ReadNC *readNC, int fileId)
	Get appropriate format to read the file. More...
static std::string	get_default_ncformat_options (ReadNC::NCFormatType format)
	Get appropriate format to read the file. More...
static NCHelper *	get_nc_helper (ReadNC *readNC, int fileId, const FileOptions &opts, EntityHandle fileSet)
	Get appropriate helper instance for ReadNC class. More...

Private Member Functions
virtual ErrorCode	init_mesh_vals ()
	Implementation of NCHelper::init_mesh_vals() More...
virtual ErrorCode	check_existing_mesh ()
	Implementation of NCHelper::check_existing_mesh() More...
virtual ErrorCode	create_mesh (Range &faces)
	Implementation of NCHelper::create_mesh() More...
virtual std::string	get_mesh_type_name ()
	Implementation of NCHelper::get_mesh_type_name() More...
virtual ErrorCode	read_ucd_variables_to_nonset_allocate (std::vector< ReadNC::VarData > &vdatas, std::vector< int > &tstep_nums)
	Implementation of UcdNCHelper::read_ucd_variables_to_nonset_allocate() More...
virtual ErrorCode	read_ucd_variables_to_nonset (std::vector< ReadNC::VarData > &vdatas, std::vector< int > &tstep_nums)
	Implementation of UcdNCHelper::read_ucd_variables_to_nonset() More...

Private Attributes
int	_spectralOrder
int	connectId
bool	isConnFile

Additional Inherited Members
Protected Member Functions inherited from moab::UcdNCHelper
template<typename T >
void	kji_to_jik_stride (size_t, size_t nj, size_t nk, void *dest, T *source, Range &localGid)
	This version takes as input the moab range, from which we actually need just the size of each sequence, for a proper transpose of the data. More...
Protected Member Functions inherited from moab::NCHelper
ErrorCode	read_variables_setup (std::vector< std::string > &var_names, std::vector< int > &tstep_nums, std::vector< ReadNC::VarData > &vdatas, std::vector< ReadNC::VarData > &vsetdatas)
	Separate set and non-set variables (common to scd mesh and ucd mesh) More...
ErrorCode	read_variables_to_set (std::vector< ReadNC::VarData > &vdatas, std::vector< int > &tstep_nums)
	Read set variables (common to scd mesh and ucd mesh) More...
ErrorCode	read_coordinate (const char *var_name, int lmin, int lmax, std::vector< double > &cvals)
ErrorCode	get_tag_to_set (ReadNC::VarData &var_data, int tstep_num, Tag &tagh)
ErrorCode	get_tag_to_nonset (ReadNC::VarData &var_data, int tstep_num, Tag &tagh, int num_lev)
ErrorCode	create_attrib_string (const std::map< std::string, ReadNC::AttData > &attMap, std::string &attString, std::vector< int > &attLen)
	Create a character string attString of attMap. with '\0' terminating each attribute name, ';' separating the data type and value, and ';' separating one name/data type/value from the next'. attLen stores the end position for each name/data type/ value. More...
ErrorCode	create_dummy_variables ()
	For a dimension that does not have a corresponding coordinate variable (e.g. ncol for HOMME), create a dummy variable with a sparse tag to store the dimension length. More...
Protected Attributes inherited from moab::UcdNCHelper
int	nCells
	Dimensions of global grid in file. More...
int	nEdges
int	nVertices
int	nLocalCells
	Dimensions of my local part of grid. More...
int	nLocalEdges
int	nLocalVertices
std::vector< double >	xVertVals
	Coordinate values for vertices. More...
std::vector< double >	yVertVals
std::vector< double >	zVertVals
int	cDim
	Dimension numbers for nCells, nEdges and nVertices. More...
int	eDim
int	vDim
Range	localGidCells
	Local global ID for cells, edges and vertices. More...
Range	localGidEdges
Range	localGidVerts
Protected Attributes inherited from moab::NCHelper
ReadNC *	_readNC
	Allow NCHelper to directly access members of ReadNC. More...
int	_fileId
	Cache some information from ReadNC. More...
const FileOptions &	_opts
EntityHandle	_fileSet
int	nTimeSteps
	Dimensions of time and level. More...
int	nLevels
std::vector< double >	tVals
	Values for time and level. More...
std::vector< double >	levVals
int	tDim
	Dimension numbers for time and level. More...
int	levDim
std::set< std::string >	ignoredVarNames
	Ignored variables. More...
std::set< std::string >	dummyVarNames
	Dummy variables. More...

Detailed Description

Child helper class for HOMME grid (CAM_SE)

Definition at line 18 of file NCHelperHOMME.hpp.

Constructor & Destructor Documentation

NCHelperHOMME()

moab::NCHelperHOMME::NCHelperHOMME	(	ReadNC *	readNC,
		int	fileId,
		const FileOptions &	opts,
		EntityHandle	fileSet
	)

Definition at line 11 of file NCHelperHOMME.cpp.

    : UcdNCHelper( readNC, fileId, opts, fileSet ), _spectralOrder( -1 ), connectId( -1 ), isConnFile( false )
{
    // Calculate spectral order
    std::map< std::string, ReadNC::AttData >::iterator attIt = readNC->globalAtts.find( "np" );
    if( attIt != readNC->globalAtts.end() )
    {
        int success = NCFUNC( get_att_int )( readNC->fileId, attIt->second.attVarId, attIt->second.attName.c_str(),
                                             &_spectralOrder );
        if( 0 == success ) _spectralOrder--;  // Spectral order is one less than np
    }
    else
    {
        // As can_read_file() returns true and there is no global attribute "np", it should be a
        // connectivity file
        isConnFile     = true;
        _spectralOrder = 3;  // Assume np is 4
    }
}

References _spectralOrder, moab::ReadNC::fileId, moab::ReadNC::globalAtts, isConnFile, and NCFUNC.

Member Function Documentation

can_read_file()

bool moab::NCHelperHOMME::can_read_file ( ReadNC *  readNC, int  fileId  )

static

Definition at line 31 of file NCHelperHOMME.cpp.

{
    // If global attribute "np" exists then it should be the HOMME grid
    if( readNC->globalAtts.find( "np" ) != readNC->globalAtts.end() )
    {
        // Make sure it is CAM grid
        std::map< std::string, ReadNC::AttData >::iterator attIt = readNC->globalAtts.find( "source" );
        if( attIt == readNC->globalAtts.end() ) return false;
        unsigned int sz = attIt->second.attLen;
        std::string att_data;
        att_data.resize( sz + 1 );
        att_data[sz] = '\000';
        int success =
            NCFUNC( get_att_text )( fileId, attIt->second.attVarId, attIt->second.attName.c_str(), &att_data[0] );
        if( success ) return false;
        if( att_data.find( "CAM" ) == std::string::npos ) return false;
 
        return true;
    }
    else
    {
        // If dimension names "ncol" AND "ncorners" AND "ncells" exist, then it should be the HOMME
        // connectivity file In this case, the mesh can still be created
        std::vector< std::string >& dimNames = readNC->dimNames;
        if( ( std::find( dimNames.begin(), dimNames.end(), std::string( "ncol" ) ) != dimNames.end() ) &&
            ( std::find( dimNames.begin(), dimNames.end(), std::string( "ncorners" ) ) != dimNames.end() ) &&
            ( std::find( dimNames.begin(), dimNames.end(), std::string( "ncells" ) ) != dimNames.end() ) )
            return true;
    }
 
    return false;
}

References moab::ReadNC::dimNames, moab::ReadNC::globalAtts, and NCFUNC.

Referenced by moab::NCHelper::get_nc_format().

check_existing_mesh()

ErrorCode moab::NCHelperHOMME::check_existing_mesh ( )

privatevirtual

Implementation of NCHelper::check_existing_mesh()

Implements moab::NCHelper.

Definition at line 222 of file NCHelperHOMME.cpp.

{
    Interface*& mbImpl = _readNC->mbImpl;
    Tag& mGlobalIdTag  = _readNC->mGlobalIdTag;
    bool& noMesh       = _readNC->noMesh;
 
    if( noMesh && localGidVerts.empty() )
    {
        // We need to populate localGidVerts range with the gids of vertices from current file set
        // localGidVerts is important in reading the variable data into the nodes
        // Also, for our purposes, localGidVerts is truly the GLOBAL_ID tag data, not other
        // file_id tags that could get passed around in other scenarios for parallel reading
 
        // Get all vertices from current file set (it is the input set in no_mesh scenario)
        Range local_verts;
        MB_CHK_SET_ERR( mbImpl->get_entities_by_dimension( _fileSet, 0, local_verts ),
                        "Trouble getting local vertices in current file set" );
 
        if( !local_verts.empty() )
        {
            std::vector< int > gids( local_verts.size() );
 
            // !IMPORTANT : this has to be the GLOBAL_ID tag
            MB_CHK_SET_ERR( mbImpl->tag_get_data( mGlobalIdTag, local_verts, &gids[0] ),
                            "Trouble getting local gid values of vertices" );
 
            // Restore localGidVerts
            std::copy( gids.rbegin(), gids.rend(), range_inserter( localGidVerts ) );
            nLocalVertices = localGidVerts.size();
        }
    }
 
    return MB_SUCCESS;
}

References moab::NCHelper::_fileSet, moab::NCHelper::_readNC, moab::Range::empty(), moab::Interface::get_entities_by_dimension(), moab::UcdNCHelper::localGidVerts, MB_CHK_SET_ERR, MB_SUCCESS, moab::ReadNC::mbImpl, moab::ReadNC::mGlobalIdTag, moab::UcdNCHelper::nLocalVertices, moab::ReadNC::noMesh, moab::Range::size(), and moab::Interface::tag_get_data().

create_mesh()

ErrorCode moab::NCHelperHOMME::create_mesh ( Range &  faces )

privatevirtual

Implementation of NCHelper::create_mesh()

Implements moab::NCHelper.

Definition at line 257 of file NCHelperHOMME.cpp.

{
    Interface*& mbImpl         = _readNC->mbImpl;
    std::string& fileName      = _readNC->fileName;
    Tag& mGlobalIdTag          = _readNC->mGlobalIdTag;
    const Tag*& mpFileIdTag    = _readNC->mpFileIdTag;
    DebugOutput& dbgOut        = _readNC->dbgOut;
    bool& spectralMesh         = _readNC->spectralMesh;
    int& gatherSetRank         = _readNC->gatherSetRank;
    int& trivialPartitionShift = _readNC->trivialPartitionShift;
 
    int rank  = 0;
    int procs = 1;
#ifdef MOAB_HAVE_MPI
    bool& isParallel = _readNC->isParallel;
    if( isParallel )
    {
        ParallelComm*& myPcomm = _readNC->myPcomm;
        rank                   = myPcomm->proc_config().proc_rank();
        procs                  = myPcomm->proc_config().proc_size();
    }
#endif
 
    int success = 0;
 
    // Need to get/read connectivity data before creating elements
    std::string conn_fname;
 
    if( isConnFile )
    {
        // Connectivity file has already been read
        connectId = _readNC->fileId;
    }
    else
    {
        // Try to open the connectivity file through CONN option, if used
        if( MB_SUCCESS != _opts.get_str_option( "CONN", conn_fname ) )
        {
            // Default convention for reading HOMME is a file HommeMapping.nc in same dir as data
            // file
            conn_fname = std::string( fileName );
            size_t idx = conn_fname.find_last_of( "/" );
            if( idx != std::string::npos )
                conn_fname = conn_fname.substr( 0, idx ).append( "/HommeMapping.nc" );
            else
                conn_fname = "HommeMapping.nc";
        }
        // Format-aware open via the runtime dispatch layer. Mirrors the
        // open path in ReadNC::load_file — probe the on-disk format on
        // rank 0, broadcast, pick the backend, then dispatch.
        int connFormat = NCFMT_UNKNOWN;
#ifdef MOAB_HAVE_MPI
        if( isParallel )
        {
            ParallelComm*& myPcomm = _readNC->myPcomm;
            int rank               = myPcomm->proc_config().proc_rank();
            if( rank == 0 ) connFormat = mbnc_detect_format( conn_fname.c_str() );
            MPI_Bcast( &connFormat, 1, MPI_INT, 0, myPcomm->proc_config().proc_comm() );
        }
        else
#endif
        {
            connFormat = mbnc_detect_format( conn_fname.c_str() );
        }
 
#ifdef MOAB_HAVE_MPI
        const int conn_mpi_size = isParallel ? _readNC->myPcomm->proc_config().proc_size() : 1;
#else
        const int conn_mpi_size = 1;
#endif
 
        const NcBackend connBackend = mbnc_choose_backend_for_read( connFormat, conn_mpi_size );
        if( connBackend == NCB_NONE )
        {
            MB_SET_ERR( MB_FAILURE, "Cannot find a compatible reader for HOMME connectivity file '"
                                        << conn_fname << "'" );
        }
 
#ifdef MOAB_HAVE_MPI
        if( connBackend == NCB_NETCDF_SERIAL || conn_mpi_size == 1 )
        {
            success = mbnc_open( conn_fname.c_str(), 0, &connectId );
        }
        else
        {
            ParallelComm*& myPcomm = _readNC->myPcomm;
            success                = mbnc_open_par( connBackend, myPcomm->proc_config().proc_comm(), MPI_INFO_NULL,
                                                    conn_fname.c_str(), 0, &connectId );
        }
#else
        success = mbnc_open( conn_fname.c_str(), 0, &connectId );
#endif
        if( success ) MB_SET_ERR( MB_FAILURE, "Failed on open" );
    }
 
    std::vector< std::string > conn_names;
    std::vector< int > conn_vals;
    MB_CHK_SET_ERR( _readNC->get_dimensions( connectId, conn_names, conn_vals ),
                    "Failed to get dimensions for connectivity" );
 
    // Read connectivity into temporary variable
    int num_fine_quads   = 0;
    int num_coarse_quads = 0;
    int start_idx        = 0;
    std::vector< std::string >::iterator vit;
    int idx = 0;
    if( ( vit = std::find( conn_names.begin(), conn_names.end(), "ncells" ) ) != conn_names.end() )
        idx = vit - conn_names.begin();
    else if( ( vit = std::find( conn_names.begin(), conn_names.end(), "ncenters" ) ) != conn_names.end() )
        idx = vit - conn_names.begin();
    else
    {
        MB_SET_ERR( MB_FAILURE, "Failed to get number of quads" );
    }
    int num_quads = conn_vals[idx];
    if( !isConnFile && num_quads != nCells )
    {
        dbgOut.tprintf( 1,
                        "Warning: number of quads from %s and cells from %s are inconsistent; "
                        "num_quads = %d, nCells = %d.\n",
                        conn_fname.c_str(), fileName.c_str(), num_quads, nCells );
    }
 
    // Get the connectivity into tmp_conn2 and permute into tmp_conn
    int cornerVarId;
    success = NCFUNC( inq_varid )( connectId, "element_corners", &cornerVarId );
    if( success ) MB_SET_ERR( MB_FAILURE, "Failed to get variable id of 'element_corners'" );
    NCDF_SIZE tmp_starts[2] = { 0, 0 };
    NCDF_SIZE tmp_counts[2] = { 4, static_cast< NCDF_SIZE >( num_quads ) };
    std::vector< int > tmp_conn( 4 * num_quads ), tmp_conn2( 4 * num_quads );
    success = NCFUNCAG( _vara_int )( connectId, cornerVarId, tmp_starts, tmp_counts, &tmp_conn2[0] );
    if( success ) MB_SET_ERR( MB_FAILURE, "Failed to get temporary connectivity" );
    if( isConnFile )
    {
        // This data/connectivity file will be closed later in ReadNC::load_file()
    }
    else
    {
        success = NCFUNC( close )( connectId );
        if( success ) MB_SET_ERR( MB_FAILURE, "Failed on close" );
    }
    // Permute the connectivity
    for( int i = 0; i < num_quads; i++ )
    {
        tmp_conn[4 * i]     = tmp_conn2[i];
        tmp_conn[4 * i + 1] = tmp_conn2[i + 1 * num_quads];
        tmp_conn[4 * i + 2] = tmp_conn2[i + 2 * num_quads];
        tmp_conn[4 * i + 3] = tmp_conn2[i + 3 * num_quads];
    }
 
    // Need to know whether we'll be creating gather mesh later, to make sure
    // we allocate enough space in one shot
    bool create_gathers = false;
    if( rank == gatherSetRank ) create_gathers = true;
 
    // Shift rank to obtain a rotated trivial partition
    int shifted_rank = rank;
    if( procs >= 2 && trivialPartitionShift > 0 ) shifted_rank = ( rank + trivialPartitionShift ) % procs;
 
    // Compute the number of local quads, accounting for coarse or fine representation
    // spectral_unit is the # fine quads per coarse quad, or spectralOrder^2
    int spectral_unit = ( spectralMesh ? _spectralOrder * _spectralOrder : 1 );
    // num_coarse_quads is the number of quads instantiated in MOAB; if !spectralMesh,
    // num_coarse_quads = num_fine_quads
    num_coarse_quads = int( std::floor( 1.0 * num_quads / ( spectral_unit * procs ) ) );
    // start_idx is the starting index in the HommeMapping connectivity list for this proc, before
    // converting to coarse quad representation
    start_idx = 4 * shifted_rank * num_coarse_quads * spectral_unit;
    // iextra = # coarse quads extra after equal split over procs
    int iextra = num_quads % ( procs * spectral_unit );
    if( shifted_rank < iextra ) num_coarse_quads++;
    start_idx += 4 * spectral_unit * std::min( shifted_rank, iextra );
    // num_fine_quads is the number of quads in the connectivity list in HommeMapping file assigned
    // to this proc
    num_fine_quads = spectral_unit * num_coarse_quads;
 
    // Now create num_coarse_quads
    EntityHandle* conn_arr;
    EntityHandle start_vertex;
    Range tmp_range;
 
    // Read connectivity into that space
    EntityHandle* sv_ptr = NULL;
    EntityHandle start_quad;
    SpectralMeshTool smt( mbImpl, _spectralOrder );
    if( !spectralMesh )
    {
        MB_CHK_SET_ERR( _readNC->readMeshIface->get_element_connect(
                            num_coarse_quads, 4, MBQUAD, 0, start_quad, conn_arr,
                            // Might have to create gather mesh later
                            ( create_gathers ? num_coarse_quads + num_quads : num_coarse_quads ) ),
                        "Failed to create local quads" );
        tmp_range.insert( start_quad, start_quad + num_coarse_quads - 1 );
        int* tmp_conn_end = ( &tmp_conn[start_idx + 4 * num_fine_quads - 1] ) + 1;
        std::copy( &tmp_conn[start_idx], tmp_conn_end, conn_arr );
        std::copy( conn_arr, conn_arr + 4 * num_fine_quads, range_inserter( localGidVerts ) );
    }
    else
    {
        MB_CHK_SET_ERR( smt.create_spectral_elems( &tmp_conn[0], num_fine_quads, 2, tmp_range, start_idx,
                                                   &localGidVerts ),
                        "Failed to create spectral elements" );
        int count, v_per_e;
        MB_CHK_SET_ERR( mbImpl->connect_iterate( tmp_range.begin(), tmp_range.end(), conn_arr, v_per_e, count ),
                        "Failed to get connectivity of spectral elements" );
        MB_CHK_SET_ERR( mbImpl->tag_iterate( smt.spectral_vertices_tag( true ), tmp_range.begin(), tmp_range.end(),
                                             count, (void*&)sv_ptr ),
                        "Failed to get fine connectivity of spectral elements" );
    }
 
    // Create vertices
    nLocalVertices = localGidVerts.size();
    std::vector< double* > arrays;
    MB_CHK_SET_ERR(
        _readNC->readMeshIface->get_node_coords( 3, nLocalVertices, 0, start_vertex, arrays,
                                                 // Might have to create gather mesh later
                                                 ( create_gathers ? nLocalVertices + nVertices : nLocalVertices ) ),
        "Failed to create local vertices" );
 
    // Set vertex coordinates
    Range::iterator rit;
    double* xptr = arrays[0];
    double* yptr = arrays[1];
    double* zptr = arrays[2];
    int i;
    for( i = 0, rit = localGidVerts.begin(); i < nLocalVertices; i++, ++rit )
    {
        assert( *rit < xVertVals.size() + 1 );
        xptr[i] = xVertVals[( *rit ) - 1];  // lon
        yptr[i] = yVertVals[( *rit ) - 1];  // lat
    }
 
    // Convert lon/lat/rad to x/y/z
    const double pideg = acos( -1.0 ) / 180.0;
    double rad         = ( isConnFile ) ? 8000.0 : 8000.0 + levVals[0];
    for( i = 0; i < nLocalVertices; i++ )
    {
        double cosphi = cos( pideg * yptr[i] );
        double zmult  = sin( pideg * yptr[i] );
        double xmult  = cosphi * cos( xptr[i] * pideg );
        double ymult  = cosphi * sin( xptr[i] * pideg );
        xptr[i]       = rad * xmult;
        yptr[i]       = rad * ymult;
        zptr[i]       = rad * zmult;
    }
 
    // Get ptr to gid memory for vertices
    Range vert_range( start_vertex, start_vertex + nLocalVertices - 1 );
    void* data;
    int count;
    MB_CHK_SET_ERR( mbImpl->tag_iterate( mGlobalIdTag, vert_range.begin(), vert_range.end(), count, data ),
                    "Failed to iterate global id tag on local vertices" );
    assert( count == nLocalVertices );
    int* gid_data = (int*)data;
    std::copy( localGidVerts.begin(), localGidVerts.end(), gid_data );
 
    // Duplicate global id data, which will be used to resolve sharing
    if( mpFileIdTag )
    {
        MB_CHK_SET_ERR( mbImpl->tag_iterate( *mpFileIdTag, vert_range.begin(), vert_range.end(), count, data ),
                        "Failed to iterate file id tag on local vertices" );
        assert( count == nLocalVertices );
        int bytes_per_tag = 4;
        MB_CHK_SET_ERR( mbImpl->tag_get_bytes( *mpFileIdTag, bytes_per_tag ),
                        "Can't get number of bytes for file id tag" );
        if( 4 == bytes_per_tag )
        {
            gid_data = (int*)data;
            std::copy( localGidVerts.begin(), localGidVerts.end(), gid_data );
        }
        else if( 8 == bytes_per_tag )
        {  // Should be a handle tag on 64 bit machine?
            long* handle_tag_data = (long*)data;
            std::copy( localGidVerts.begin(), localGidVerts.end(), handle_tag_data );
        }
    }
 
    // Create map from file ids to vertex handles, used later to set connectivity
    std::map< EntityHandle, EntityHandle > vert_handles;
    for( rit = localGidVerts.begin(), i = 0; rit != localGidVerts.end(); ++rit, i++ )
        vert_handles[*rit] = start_vertex + i;
 
    // Compute proper handles in connectivity using offset
    for( int q = 0; q < 4 * num_coarse_quads; q++ )
    {
        conn_arr[q] = vert_handles[conn_arr[q]];
        assert( conn_arr[q] );
    }
    if( spectralMesh )
    {
        int verts_per_quad = ( _spectralOrder + 1 ) * ( _spectralOrder + 1 );
        for( int q = 0; q < verts_per_quad * num_coarse_quads; q++ )
        {
            sv_ptr[q] = vert_handles[sv_ptr[q]];
            assert( sv_ptr[q] );
        }
    }
 
    // Add new vertices and quads to current file set
    faces.merge( tmp_range );
    tmp_range.insert( start_vertex, start_vertex + nLocalVertices - 1 );
    MB_CHK_SET_ERR( mbImpl->add_entities( _fileSet, tmp_range ),
                    "Failed to add new vertices and quads to current file set" );
 
    // Mark the set with the spectral order
    Tag sporder;
    MB_CHK_SET_ERR( mbImpl->tag_get_handle( "SPECTRAL_ORDER", 1, MB_TYPE_INTEGER, sporder,
                                            MB_TAG_SPARSE | MB_TAG_CREAT ),
                    "Trouble creating SPECTRAL_ORDER tag" );
    MB_CHK_SET_ERR( mbImpl->tag_set_data( sporder, &_fileSet, 1, &_spectralOrder ),
                    "Trouble setting data to SPECTRAL_ORDER tag" );
 
    if( create_gathers )
    {
        EntityHandle gather_set;
        MB_CHK_SET_ERR( _readNC->readMeshIface->create_gather_set( gather_set ), "Failed to create gather set" );
 
        // Create vertices
        arrays.clear();
        // Don't need to specify allocation number here, because we know enough verts were created
        // before
        MB_CHK_SET_ERR( _readNC->readMeshIface->get_node_coords( 3, nVertices, 0, start_vertex, arrays ),
                        "Failed to create gather set vertices" );
 
        xptr = arrays[0];
        yptr = arrays[1];
        zptr = arrays[2];
        for( i = 0; i < nVertices; i++ )
        {
            double cosphi = cos( pideg * yVertVals[i] );
            double zmult  = sin( pideg * yVertVals[i] );
            double xmult  = cosphi * cos( xVertVals[i] * pideg );
            double ymult  = cosphi * sin( xVertVals[i] * pideg );
            xptr[i]       = rad * xmult;
            yptr[i]       = rad * ymult;
            zptr[i]       = rad * zmult;
        }
 
        // Get ptr to gid memory for vertices
        Range gather_set_verts_range( start_vertex, start_vertex + nVertices - 1 );
        MB_CHK_SET_ERR( mbImpl->tag_iterate( mGlobalIdTag, gather_set_verts_range.begin(), gather_set_verts_range.end(),
                                             count, data ),
                        "Failed to iterate global id tag on gather set vertices" );
        assert( count == nVertices );
        gid_data = (int*)data;
        for( int j = 1; j <= nVertices; j++ )
            gid_data[j - 1] = j;
        // Set the file id tag too, it should be bigger something not interfering with global id
        if( mpFileIdTag )
        {
            MB_CHK_SET_ERR( mbImpl->tag_iterate( *mpFileIdTag, gather_set_verts_range.begin(),
                                                 gather_set_verts_range.end(), count, data ),
                            "Failed to iterate file id tag on gather set vertices" );
            assert( count == nVertices );
            int bytes_per_tag = 4;
            MB_CHK_SET_ERR( mbImpl->tag_get_bytes( *mpFileIdTag, bytes_per_tag ),
                            "Can't get number of bytes for file id tag" );
            if( 4 == bytes_per_tag )
            {
                gid_data = (int*)data;
                for( int j = 1; j <= nVertices; j++ )
                    gid_data[j - 1] = nVertices + j;  // Bigger than global id tag
            }
            else if( 8 == bytes_per_tag )
            {  // Should be a handle tag on 64 bit machine?
                long* handle_tag_data = (long*)data;
                for( int j = 1; j <= nVertices; j++ )
                    handle_tag_data[j - 1] = nVertices + j;  // Bigger than global id tag
            }
        }
 
        MB_CHK_SET_ERR( mbImpl->add_entities( gather_set, gather_set_verts_range ),
                        "Failed to add vertices to the gather set" );
 
        // Create quads
        Range gather_set_quads_range;
        // Don't need to specify allocation number here, because we know enough quads were created
        // before
        MB_CHK_SET_ERR( _readNC->readMeshIface->get_element_connect( num_quads, 4, MBQUAD, 0, start_quad, conn_arr ),
                        "Failed to create gather set quads" );
        gather_set_quads_range.insert( start_quad, start_quad + num_quads - 1 );
        int* tmp_conn_end = ( &tmp_conn[4 * num_quads - 1] ) + 1;
        std::copy( &tmp_conn[0], tmp_conn_end, conn_arr );
        for( i = 0; i != 4 * num_quads; i++ )
            conn_arr[i] += start_vertex - 1;  // Connectivity array is shifted by where the gather verts start
        MB_CHK_SET_ERR( mbImpl->add_entities( gather_set, gather_set_quads_range ),
                        "Failed to add quads to the gather set" );
    }
 
    return MB_SUCCESS;
}

References moab::NCHelper::_fileSet, moab::NCHelper::_opts, moab::NCHelper::_readNC, _spectralOrder, moab::Interface::add_entities(), moab::Range::begin(), moab::Interface::connect_iterate(), connectId, moab::ReadUtilIface::create_gather_set(), moab::SpectralMeshTool::create_spectral_elems(), moab::ReadNC::dbgOut, moab::Range::end(), moab::ReadNC::fileId, moab::ReadNC::fileName, moab::ReadNC::gatherSetRank, moab::ReadNC::get_dimensions(), moab::ReadUtilIface::get_element_connect(), moab::ReadUtilIface::get_node_coords(), moab::FileOptions::get_str_option(), moab::Range::insert(), isConnFile, moab::ReadNC::isParallel, moab::NCHelper::levVals, moab::UcdNCHelper::localGidVerts, MB_CHK_SET_ERR, MB_SET_ERR, MB_SUCCESS, MB_TAG_CREAT, MB_TAG_SPARSE, MB_TYPE_INTEGER, moab::ReadNC::mbImpl, moab::mbnc_choose_backend_for_read(), moab::mbnc_detect_format(), moab::mbnc_open(), MBQUAD, moab::Range::merge(), moab::ReadNC::mGlobalIdTag, moab::ReadNC::mpFileIdTag, moab::NCB_NETCDF_SERIAL, moab::NCB_NONE, NCDF_SIZE, moab::UcdNCHelper::nCells, moab::NCFMT_UNKNOWN, NCFUNC, NCFUNCAG, moab::UcdNCHelper::nLocalVertices, moab::UcdNCHelper::nVertices, moab::pideg, moab::ProcConfig::proc_comm(), moab::ParallelComm::proc_config(), moab::ProcConfig::proc_rank(), moab::ProcConfig::proc_size(), moab::ReadNC::readMeshIface, moab::Range::size(), moab::SpectralMeshTool::spectral_vertices_tag(), moab::ReadNC::spectralMesh, moab::Interface::tag_get_bytes(), moab::Interface::tag_get_handle(), moab::Interface::tag_iterate(), moab::Interface::tag_set_data(), moab::DebugOutput::tprintf(), moab::ReadNC::trivialPartitionShift, moab::UcdNCHelper::xVertVals, and moab::UcdNCHelper::yVertVals.

get_mesh_type_name()

virtual std::string moab::NCHelperHOMME::get_mesh_type_name ( )

inlineprivatevirtual

Implementation of NCHelper::get_mesh_type_name()

Implements moab::NCHelper.

Definition at line 32 of file NCHelperHOMME.hpp.

    {
        return "CAM_SE";
    }

init_mesh_vals()

ErrorCode moab::NCHelperHOMME::init_mesh_vals ( )

privatevirtual

Implementation of NCHelper::init_mesh_vals()

Implements moab::NCHelper.

Definition at line 64 of file NCHelperHOMME.cpp.

{
    std::vector< std::string >& dimNames              = _readNC->dimNames;
    std::vector< int >& dimLens                       = _readNC->dimLens;
    std::map< std::string, ReadNC::VarData >& varInfo = _readNC->varInfo;
 
    unsigned int idx;
    std::vector< std::string >::iterator vit;
 
    // Look for time dimension
    if( isConnFile )
    {
        // Connectivity file might not have time dimension
    }
    else
    {
        if( ( vit = std::find( dimNames.begin(), dimNames.end(), "time" ) ) != dimNames.end() )
            idx = vit - dimNames.begin();
        else if( ( vit = std::find( dimNames.begin(), dimNames.end(), "t" ) ) != dimNames.end() )
            idx = vit - dimNames.begin();
        else
        {
            MB_SET_ERR( MB_FAILURE, "Couldn't find 'time' or 't' dimension" );
        }
        tDim       = idx;
        nTimeSteps = dimLens[idx];
    }
 
    // Get number of vertices (labeled as number of columns)
    if( ( vit = std::find( dimNames.begin(), dimNames.end(), "ncol" ) ) != dimNames.end() )
        idx = vit - dimNames.begin();
    else
    {
        MB_SET_ERR( MB_FAILURE, "Couldn't find 'ncol' dimension" );
    }
    vDim      = idx;
    nVertices = dimLens[idx];
 
    // Set number of cells
    nCells = nVertices - 2;
 
    // Get number of levels
    if( isConnFile )
    {
        // Connectivity file might not have level dimension
    }
    else
    {
        if( ( vit = std::find( dimNames.begin(), dimNames.end(), "lev" ) ) != dimNames.end() )
            idx = vit - dimNames.begin();
        else if( ( vit = std::find( dimNames.begin(), dimNames.end(), "ilev" ) ) != dimNames.end() )
            idx = vit - dimNames.begin();
        else
        {
            MB_SET_ERR( MB_FAILURE, "Couldn't find 'lev' or 'ilev' dimension" );
        }
        levDim  = idx;
        nLevels = dimLens[idx];
    }
 
    // Store lon values in xVertVals
    std::map< std::string, ReadNC::VarData >::iterator vmit;
    if( ( vmit = varInfo.find( "lon" ) ) != varInfo.end() && ( *vmit ).second.varDims.size() == 1 )
    {
        MB_CHK_SET_ERR( read_coordinate( "lon", 0, nVertices - 1, xVertVals ), "Trouble reading 'lon' variable" );
    }
    else
    {
        MB_SET_ERR( MB_FAILURE, "Couldn't find 'lon' variable" );
    }
 
    // Store lat values in yVertVals
    if( ( vmit = varInfo.find( "lat" ) ) != varInfo.end() && ( *vmit ).second.varDims.size() == 1 )
    {
        MB_CHK_SET_ERR( read_coordinate( "lat", 0, nVertices - 1, yVertVals ), "Trouble reading 'lat' variable" );
    }
    else
    {
        MB_SET_ERR( MB_FAILURE, "Couldn't find 'lat' variable" );
    }
 
    // Store lev values in levVals
    if( isConnFile )
    {
        // Connectivity file might not have level variable
    }
    else
    {
        if( ( vmit = varInfo.find( "lev" ) ) != varInfo.end() && ( *vmit ).second.varDims.size() == 1 )
        {
            MB_CHK_SET_ERR( read_coordinate( "lev", 0, nLevels - 1, levVals ), "Trouble reading 'lev' variable" );
 
            // Decide whether down is positive
            char posval[10] = { 0 };
            int success     = NCFUNC( get_att_text )( _fileId, ( *vmit ).second.varId, "positive", posval );
            if( 0 == success && !strcmp( posval, "down" ) )
            {
                for( std::vector< double >::iterator dvit = levVals.begin(); dvit != levVals.end(); ++dvit )
                    ( *dvit ) *= -1.0;
            }
        }
        else
        {
            MB_SET_ERR( MB_FAILURE, "Couldn't find 'lev' variable" );
        }
    }
 
    // Store time coordinate values in tVals
    if( isConnFile )
    {
        // Connectivity file might not have time variable
    }
    else
    {
        if( ( vmit = varInfo.find( "time" ) ) != varInfo.end() && ( *vmit ).second.varDims.size() == 1 )
        {
            MB_CHK_SET_ERR( read_coordinate( "time", 0, nTimeSteps - 1, tVals ), "Trouble reading 'time' variable" );
        }
        else if( ( vmit = varInfo.find( "t" ) ) != varInfo.end() && ( *vmit ).second.varDims.size() == 1 )
        {
            MB_CHK_SET_ERR( read_coordinate( "t", 0, nTimeSteps - 1, tVals ), "Trouble reading 't' variable" );
        }
        else
        {
            // If expected time variable does not exist, set dummy values to tVals
            for( int t = 0; t < nTimeSteps; t++ )
                tVals.push_back( (double)t );
        }
    }
 
    // For each variable, determine the entity location type and number of levels
    std::map< std::string, ReadNC::VarData >::iterator mit;
    for( mit = varInfo.begin(); mit != varInfo.end(); ++mit )
    {
        ReadNC::VarData& vd = ( *mit ).second;
 
        // Default entLoc is ENTLOCSET
        if( std::find( vd.varDims.begin(), vd.varDims.end(), tDim ) != vd.varDims.end() )
        {
            if( std::find( vd.varDims.begin(), vd.varDims.end(), vDim ) != vd.varDims.end() )
                vd.entLoc = ReadNC::ENTLOCVERT;
        }
 
        // Default numLev is 0
        if( std::find( vd.varDims.begin(), vd.varDims.end(), levDim ) != vd.varDims.end() ) vd.numLev = nLevels;
    }
 
    // Hack: create dummy variables for dimensions (like ncol) with no corresponding coordinate
    // variables
    MB_CHK_SET_ERR( create_dummy_variables(), "Failed to create dummy variables" );
 
    return MB_SUCCESS;
}

References moab::NCHelper::_fileId, moab::NCHelper::_readNC, moab::NCHelper::create_dummy_variables(), moab::ReadNC::dimLens, moab::ReadNC::dimNames, moab::ReadNC::VarData::entLoc, moab::ReadNC::ENTLOCVERT, isConnFile, moab::NCHelper::levDim, moab::NCHelper::levVals, MB_CHK_SET_ERR, MB_SET_ERR, MB_SUCCESS, moab::UcdNCHelper::nCells, NCFUNC, moab::NCHelper::nLevels, moab::NCHelper::nTimeSteps, moab::ReadNC::VarData::numLev, moab::UcdNCHelper::nVertices, moab::NCHelper::read_coordinate(), moab::NCHelper::tDim, moab::NCHelper::tVals, moab::ReadNC::VarData::varDims, moab::ReadNC::varInfo, moab::UcdNCHelper::vDim, moab::UcdNCHelper::xVertVals, and moab::UcdNCHelper::yVertVals.

read_ucd_variables_to_nonset()

ErrorCode moab::NCHelperHOMME::read_ucd_variables_to_nonset ( std::vector< ReadNC::VarData > &  vdatas, std::vector< int > &  tstep_nums  )

privatevirtual

Implementation of UcdNCHelper::read_ucd_variables_to_nonset()

Implements moab::UcdNCHelper.

Definition at line 838 of file NCHelperHOMME.cpp.

{
    DebugOutput& dbgOut = _readNC->dbgOut;
 
    MB_CHK_SET_ERR( read_ucd_variables_to_nonset_allocate( vdatas, tstep_nums ),
                    "Trouble allocating space to read non-set variables" );
 
    // Finally, read into that space
    int success;
    for( unsigned int i = 0; i < vdatas.size(); i++ )
    {
        std::size_t sz = vdatas[i].sz;
 
        // A typical supported variable: float T(time, lev, ncol)
        // For tag values, need transpose (lev, ncol) to (ncol, lev)
        size_t ni = vdatas[i].readCounts[2];  // ncol
        size_t nj = 1;                        // Here we should just set nj to 1
        size_t nk = vdatas[i].readCounts[1];  // lev
 
        for( unsigned int t = 0; t < tstep_nums.size(); t++ )
        {
            // Tag data for this timestep
            void* data = vdatas[i].varDatas[t];
 
            // Set readStart for each timestep along time dimension
            vdatas[i].readStarts[0] = tstep_nums[t];
 
            switch( vdatas[i].varDataType )
            {
                case NC_FLOAT:
                case NC_DOUBLE: {
                    // Read float as double
                    std::vector< double > tmpdoubledata( sz );
 
                    // In the case of ucd mesh, and on multiple proc,
                    // we need to read as many times as subranges we have in the
                    // localGidVerts range;
                    // basically, we have to give a different point
                    // for data to start, for every subrange :(
                    size_t indexInDoubleArray = 0;
                    size_t ic                 = 0;
                    for( Range::pair_iterator pair_iter = localGidVerts.pair_begin();
                         pair_iter != localGidVerts.pair_end(); ++pair_iter, ic++ )
                    {
                        EntityHandle starth     = pair_iter->first;
                        EntityHandle endh       = pair_iter->second;  // Inclusive
                        vdatas[i].readStarts[2] = (NCDF_SIZE)( starth - 1 );
                        vdatas[i].readCounts[2] = (NCDF_SIZE)( endh - starth + 1 );
 
                        success = NCFUNCAG( _vara_double )( _fileId, vdatas[i].varId, &( vdatas[i].readStarts[0] ),
                                                            &( vdatas[i].readCounts[0] ),
                                                            &( tmpdoubledata[indexInDoubleArray] ) );
                        if( success )
                            MB_SET_ERR( MB_FAILURE,
                                        "Failed to read double data in a loop for variable " << vdatas[i].varName );
                        // We need to increment the index in double array for the
                        // next subrange
                        indexInDoubleArray += ( endh - starth + 1 ) * 1 * vdatas[i].numLev;
                    }
                    assert( ic == localGidVerts.psize() );
 
                    if( vdatas[i].numLev > 1 )
                        // Transpose (lev, ncol) to (ncol, lev)
                        kji_to_jik_stride( ni, nj, nk, data, &tmpdoubledata[0], localGidVerts );
                    else
                    {
                        for( std::size_t idx = 0; idx != tmpdoubledata.size(); idx++ )
                            ( (double*)data )[idx] = tmpdoubledata[idx];
                    }
 
                    break;
                }
                default:
                    MB_SET_ERR( MB_FAILURE, "Unexpected data type for variable " << vdatas[i].varName );
            }
        }
    }
 
    // Debug output, if requested
    if( 1 == dbgOut.get_verbosity() )
    {
        dbgOut.printf( 1, "Read variables: %s", vdatas.begin()->varName.c_str() );
        for( unsigned int i = 1; i < vdatas.size(); i++ )
            dbgOut.printf( 1, ", %s ", vdatas[i].varName.c_str() );
        dbgOut.tprintf( 1, "\n" );
    }
 
    return MB_SUCCESS;
}

References moab::NCHelper::_fileId, moab::NCHelper::_readNC, moab::ReadNC::dbgOut, moab::DebugOutput::get_verbosity(), moab::UcdNCHelper::kji_to_jik_stride(), moab::UcdNCHelper::localGidVerts, MB_CHK_SET_ERR, MB_SET_ERR, MB_SUCCESS, NCDF_SIZE, NCFUNCAG, moab::Range::pair_begin(), moab::Range::pair_end(), moab::DebugOutput::printf(), moab::Range::psize(), read_ucd_variables_to_nonset_allocate(), and moab::DebugOutput::tprintf().

read_ucd_variables_to_nonset_allocate()

ErrorCode moab::NCHelperHOMME::read_ucd_variables_to_nonset_allocate ( std::vector< ReadNC::VarData > &  vdatas, std::vector< int > &  tstep_nums  )

privatevirtual

Implementation of UcdNCHelper::read_ucd_variables_to_nonset_allocate()

Implements moab::UcdNCHelper.

Definition at line 649 of file NCHelperHOMME.cpp.

{
    Interface*& mbImpl          = _readNC->mbImpl;
    std::vector< int >& dimLens = _readNC->dimLens;
    DebugOutput& dbgOut         = _readNC->dbgOut;
 
    Range* range = NULL;
 
    // Get vertices
    Range verts;
    MB_CHK_SET_ERR( mbImpl->get_entities_by_dimension( _fileSet, 0, verts ),
                    "Trouble getting vertices in current file set" );
    assert( "Should only have a single vertex subrange, since they were read in one shot" && verts.psize() == 1 );
 
    for( unsigned int i = 0; i < vdatas.size(); i++ )
    {
        // Support non-set variables with 3 dimensions like (time, lev, ncol)
        assert( 3 == vdatas[i].varDims.size() );
 
        // For a non-set variable, time should be the first dimension
        assert( tDim == vdatas[i].varDims[0] );
 
        // Set up readStarts and readCounts
        vdatas[i].readStarts.resize( 3 );
        vdatas[i].readCounts.resize( 3 );
 
        // First: time
        vdatas[i].readStarts[0] = 0;  // This value is timestep dependent, will be set later
        vdatas[i].readCounts[0] = 1;
 
        // Next: lev
        vdatas[i].readStarts[1] = 0;
        vdatas[i].readCounts[1] = vdatas[i].numLev;
 
        // Finally: ncol
        switch( vdatas[i].entLoc )
        {
            case ReadNC::ENTLOCVERT:
                // Vertices
                // Start from the first localGidVerts
                // Actually, this will be reset later on in a loop
                vdatas[i].readStarts[2] = localGidVerts[0] - 1;
                vdatas[i].readCounts[2] = nLocalVertices;
                range                   = &verts;
                break;
            default:
                MB_SET_ERR( MB_FAILURE, "Unexpected entity location type for variable " << vdatas[i].varName );
        }
 
        // Get variable size
        vdatas[i].sz = 1;
        for( std::size_t idx = 0; idx != 3; idx++ )
            vdatas[i].sz *= vdatas[i].readCounts[idx];
 
        for( unsigned int t = 0; t < tstep_nums.size(); t++ )
        {
            dbgOut.tprintf( 2, "Reading variable %s, time step %d\n", vdatas[i].varName.c_str(), tstep_nums[t] );
 
            if( tstep_nums[t] >= dimLens[tDim] )
            {
                MB_SET_ERR( MB_INDEX_OUT_OF_RANGE, "Wrong value for timestep number " << tstep_nums[t] );
            }
 
            // Get the tag to read into
            if( !vdatas[i].varTags[t] )
            {
                MB_CHK_SET_ERR( get_tag_to_nonset( vdatas[i], tstep_nums[t], vdatas[i].varTags[t], vdatas[i].numLev ),
                                "Trouble getting tag for variable " << vdatas[i].varName );
            }
 
            // Get ptr to tag space
            void* data;
            int count;
            MB_CHK_SET_ERR( mbImpl->tag_iterate( vdatas[i].varTags[t], range->begin(), range->end(), count, data ),
                            "Failed to iterate tag for variable " << vdatas[i].varName );
            assert( (unsigned)count == range->size() );
            vdatas[i].varDatas[t] = data;
        }
    }
 
    return MB_SUCCESS;
}

References moab::NCHelper::_fileSet, moab::NCHelper::_readNC, moab::Range::begin(), moab::ReadNC::dbgOut, moab::ReadNC::dimLens, moab::Range::end(), moab::ReadNC::ENTLOCVERT, moab::Interface::get_entities_by_dimension(), moab::NCHelper::get_tag_to_nonset(), moab::UcdNCHelper::localGidVerts, MB_CHK_SET_ERR, MB_INDEX_OUT_OF_RANGE, MB_SET_ERR, MB_SUCCESS, moab::ReadNC::mbImpl, moab::UcdNCHelper::nLocalVertices, moab::Range::psize(), moab::Range::size(), moab::Interface::tag_iterate(), moab::NCHelper::tDim, and moab::DebugOutput::tprintf().

Referenced by read_ucd_variables_to_nonset().

Member Data Documentation

_spectralOrder

int moab::NCHelperHOMME::_spectralOrder

private

Definition at line 51 of file NCHelperHOMME.hpp.

Referenced by create_mesh(), and NCHelperHOMME().

connectId

int moab::NCHelperHOMME::connectId

private

Definition at line 52 of file NCHelperHOMME.hpp.

Referenced by create_mesh().

isConnFile

bool moab::NCHelperHOMME::isConnFile

private

Definition at line 53 of file NCHelperHOMME.hpp.

Referenced by create_mesh(), init_mesh_vals(), and NCHelperHOMME().

---

The documentation for this class was generated from the following files:

• NCHelperHOMME.hpp
• NCHelperHOMME.cpp