ExchangeHalos.cpp

Go to the documentation of this file.

/** @example ExchangeHalos Driver
 * \brief Example program that shows the use case for performing tag data exchange
 * between parallel processors in order to sync data on shared entities.
 *
 * <b>This example </b>:
 *    -# Initialize MPI and instantiates MOAB
 *    -# Gets user options: Input mesh file name, vector tag length, ghost layer size etc
 *    -# Create the root and partition sets
 *    -# Instantiate ParallelComm and read the mesh file in parallel using appropriate options
 *    -# Create the required number of ghost layers as requested by the user (default = 3)
 *    -# Get 2D MPAS polygonal entities in the mesh and filter to get only the "owned" entities
 *    -# Create two tags: scalar_variable (single data/cell) and vector_variable (multiple data/cell)
 *    -# Set tag data using analytical functions for both scalar and vector fields on owned entities
 *    -# Exchange shared entity information and tags between processors
 *      -# If debugging is turned on, store mesh file and tag on root process (will not contain data on shared entities)
 *      -# Perform exchange of scalar tag data on shared entities
 *      -# Perform exchange of vector tag data on shared entities
 *      -# If debugging is turned on, store mesh file and tag on root process (will now contain data on *all* entities)
 *    -#  Destroy the MOAB instance and finalize MPI
 *
 * <b>To run: </b>
 *      mpiexec -n np ./ExchangeHalos --input <mpas_mesh_file> --nghosts <ghostlayers> --vtaglength <vector component
 * size> \
 *                    --nexchanges <number of exchange runs>
 * <b>Example:</b>
 *      mpiexec -n 4 ./ExchangeHalos --input $MOAB_DIR/MeshFiles/unittest/io/mpasx1.642.t.2.nc --nghosts 3 --vtaglength 100
 *
 * NOTE: --debug option can be added to write out extra files in h5m format to visualize some output (written from root
 * task only)
 *
 */
 
// MOAB includes
#include "moab/Core.hpp"
#include "moab/CpuTimer.hpp"
#include "moab/ProgOptions.hpp"
 
#ifndef MOAB_HAVE_MPI
#error "Please build MOAB with MPI..."
#endif
 
#include "moab/ParallelComm.hpp"
#include "MBParallelConventions.h"
 
// C++ includes
#include <iostream>
#include <string>
 
using namespace moab;
using namespace std;
#define dbgprint( MSG )                                           \
    do                                                            \
    {                                                             \
        if( context.proc_id == 0 ) std::cout << MSG << std::endl; \
    } while( false )
 
#define runchk( CODE, MSG )         \
    do                              \
    {                               \
        moab::ErrorCode err = CODE; \
        MB_CHK_SET_ERR( err, MSG ); \
    } while( false )
 
#define runchk_cont( CODE, MSG )                               \
    do                                                         \
    {                                                          \
        moab::ErrorCode err = CODE;                            \
        MB_CHK_ERR_CONT( err );                                \
        if( err ) std::cout << "Error:: " << MSG << std::endl; \
    } while( false )
 
// Run time context structure
 
/// @brief The RunttimeContext is an example specific class to store
/// the run specific input data, MOAB datastructures used during the run
/// and provides other utility functions to profile operations etc
struct RuntimeContext
{
  public:
    int dimension{ 2 };           /// dimension of the problem
    std::string input_filename;   /// input file name (nc format)
    std::string output_filename;  /// output file name (h5m format)
    int ghost_layers{ 3 };        /// number of ghost layers
    std::string scalar_tagname;   /// scalar tag name
    std::string vector_tagname;   /// vector tag name
    int vector_length{ 3 };       /// length of the vector tag components
    int num_max_exchange{ 10 };   /// total number of exchange iterations
    bool debug_output{ false };   /// write debug output information?
    int proc_id{ 1 };             /// process identifier
    int num_procs{ 1 };           /// total number of processes
    double last_counter{ 0.0 };   /// last time counter between push/pop timer
 
    // MOAB objects
    moab::Core moab_interface;
    moab::ParallelComm* parallel_communicator;
    moab::EntityHandle fileset{ 0 }, partnset{ 0 };
 
    /// @brief Constructor: allocate MOAB interface and communicator, and initialize
    /// other data members with some default values
    explicit RuntimeContext( MPI_Comm comm = MPI_COMM_WORLD );
 
    /// @brief Destructor: deallocate MOAB interface and communicator
    ~RuntimeContext();
 
    /// @brief Parse the runtime command line options
    /// @param argc - number of command line arguments
    /// @param argv - command line arguments as string list
    void ParseCLOptions( int argc, char* argv[] );
 
    /// @brief Measure and start the timer to profile a task
    /// @param operation String name of the task being measured
    inline void timer_push( const std::string& operation );
 
    /// @brief Stop the timer and store the elapsed duration
    /// @param nruns Optional argument used to average the measured time
    void timer_pop( const int nruns = 1 );
 
    /// @brief Return the last elapsed time
    /// @return last_counter from timer_pop was called
    inline double last_elapsed() const;
 
    /// @brief Load a MOAB supported file (h5m or nc format) from disk
    ///        representing an MPAS mesh
    /// @param load_ghosts Optional boolean to specify whether to load ghosts
    ///                    when reading the file (only relevant for h5m)
    /// @return Error code if any (else MB_SUCCESS)
    moab::ErrorCode load_file( bool load_ghosts = false );
 
    /// @brief Create scalar and vector tags in the MOAB mesh instance
    /// @param tagScalar Tag reference to the scalar field
    /// @param tagVector Tag reference to the vector field
    /// @param entities Entities on which both the scalar and vector fields are defined
    /// @return Error code if any (else MB_SUCCESS)
    moab::ErrorCode create_sv_tags( moab::Tag& tagScalar, moab::Tag& tagVector, moab::Range& entities );
 
    /// @brief Evaluate some closed-form Spherical Harmonic functions with an optional multiplier term
    /// @param lon Longitude in lat-lon space
    /// @param lat Latitude  in lat-lon space
    /// @param type Function type
    /// @param multiplier Optional multiplier to scale value (default=1.0)
    /// @return value of the evaluated function
    inline double evaluate_function( double lon, double lat, int type = 1, double multiplier = 1.0 ) const
    {
        switch( type )
        {
            case 1:
                return ( 2.0 + std::pow( sin( 2.0 * lat ), 16.0 ) * cos( 16.0 * lon ) ) * multiplier;
            default:
                return ( 2.0 + cos( lon ) * cos( lon ) * cos( 2.0 * lat ) ) * multiplier;
        }
    }
 
  private:
    /// @brief Compute the centroids of elements in 2D lat/lon space
    /// @param entities Entities to compute centroids
    /// @return Vector of centroids (as lat/lon)
    std::vector< double > compute_centroids( const moab::Range& entities ) const;
 
    moab::CpuTimer mTimer;
    double mTimerOps{ 0.0 };
    std::string mOpName;
};
 
//
// Start of main test program
//
int main( int argc, char** argv )
{
    // Initialize MPI first
    MPI_Init( &argc, &argv );
 
    {
        // Create our context for this example run
        RuntimeContext context;
        dbgprint( "********** Exchange halos example **********\n" );
 
        // Get the input options
        context.ParseCLOptions( argc, argv );
 
        /////////////////////////////////////////////////////////////////////////
        // Print out the input parameters in use
        dbgprint( " -- Input Parameters -- " );
        dbgprint( "    Number of Processes  = " << context.num_procs );
        dbgprint( "    Input mesh           = " << context.input_filename );
        dbgprint( "    Ghost Layers         = " << context.ghost_layers );
        dbgprint( "    Scalar Tag name      = " << context.scalar_tagname );
        dbgprint( "    Vector Tag name      = " << context.vector_tagname );
        dbgprint( "    Vector Tag length    = " << context.vector_length << endl );
        /////////////////////////////////////////////////////////////////////////
 
        // Timer storage for all phases
        double elapsed_times[4];
 
        // Read the input file specified by user, in parallel, using appropriate options
        // Supports reading partitioned h5m files and MPAS nc files directly with online Zoltan partitioning
        context.timer_push( "Read input file" );
        {
            // Load the file from disk with given options
            runchk( context.load_file( false ), "MOAB::load_file failed for filename: " << context.input_filename );
        }
        context.timer_pop();
        elapsed_times[0] = context.last_elapsed();
 
        // Let the actual measurements begin...
        dbgprint( "\n- Starting execution -\n" );
 
        // We need to set up the ghost layers requested by the user. First correct for thin layers and then
        // call `exchange_ghost_cells` to prepare the mesh for use with halo regions
        context.timer_push( "Setup ghost layers" );
        {
            // Loop over the number of ghost layers needed and ask MOAB for layers 1 at a time
            for( int ighost = 0; ighost < context.ghost_layers; ++ighost )
            {
                // Exchange ghost cells
                int ghost_dimension  = context.dimension;
                int bridge_dimension = context.dimension - 1;
                // Let us now get all ghost layers from adjacent parts
                runchk( context.parallel_communicator->exchange_ghost_cells(
                            ghost_dimension, bridge_dimension, ( ighost + 1 ), 0, true /* store_remote_handles */,
                            true /* wait_all */, &context.fileset ),
                        "Exchange ghost cells failed" );  // true to store remote handles
 
                // Ensure that all processes understand about multi-shared vertices and entities
                // in case some adjacent parts are only m layers thick (where m < context.ghost_layers)
                if( ighost < context.ghost_layers - 1 )
                    runchk( context.parallel_communicator->correct_thin_ghost_layers(),
                            "Thin layer correction failed" );
            }
        }
        context.timer_pop();
        elapsed_times[1] = context.last_elapsed();
 
        // Get the 2D MPAS elements and filter it so that we have only owned elements
        Range dimEnts;
        {
            // Get all entities of dimension = dim
            runchk( context.moab_interface.get_entities_by_dimension( context.fileset, context.dimension, dimEnts ),
                    "Getting 2D entities failed" );
            // Get only owned entities! The ghosted/shared entities will get their data when we exchange
            // So let us filter entities based on the status: NOT x NOT_OWNED = OWNED status :-)
            runchk( context.parallel_communicator->filter_pstatus( dimEnts, PSTATUS_NOT_OWNED, PSTATUS_NOT ),
                    "Filtering pstatus failed" );
 
            // Aggregate the total number of elements in the mesh
            auto numEntities     = dimEnts.size();
            int numTotalEntities = 0;
            MPI_Reduce( &numEntities, &numTotalEntities, 1, MPI_INT, MPI_SUM, 0,
                        context.parallel_communicator->proc_config().proc_comm() );
 
            // We expect the total number of elements to be constant, immaterial of number of processes.
            // If not, we have a bug!
            dbgprint( "Total number of " << context.dimension << "D elements in the mesh = " << numTotalEntities );
        }
 
        Tag tagScalar = nullptr;
        Tag tagVector = nullptr;
        // Create two tag handles: scalar_variable and vector_variable
        // Set these tags with appropriate closed form functional data
        // based on element centroid information
        runchk( context.create_sv_tags( tagScalar, tagVector, dimEnts ), "Unable to create scalar and vector tags" );
 
        // let us write out the local mesh before tag_exchange is called
        // we expect to see data only on the owned entities - and ghosted entities should have default values
        if( context.debug_output && ( context.proc_id == 0 ) )  // only on root process, for debugging
        {
            dbgprint( "> Writing to file *before* ghost exchange " );
            runchk( context.moab_interface.write_file( "exchangeHalos_output_rank0_pre.h5m", "H5M", "" ),
                    "Writing to disk failed" );
        }
 
        // Perform exchange of tag data between neighboring tasks
        dbgprint( "> Exchanging tags between processors " );
        context.timer_push( "Exchange scalar tag data" );
        for( auto irun = 0; irun < context.num_max_exchange; ++irun )
        {
            // Exchange scalar tags between processors
            runchk( context.parallel_communicator->exchange_tags( tagScalar, dimEnts ),
                    "Exchanging scalar tag between processors failed" );
        }
        context.timer_pop( context.num_max_exchange );
        elapsed_times[2] = context.last_elapsed();
 
        context.timer_push( "Exchange vector tag data" );
        for( auto irun = 0; irun < context.num_max_exchange; ++irun )
        {
            // Exchange vector tags between processors
            runchk( context.parallel_communicator->exchange_tags( tagVector, dimEnts ),
                    "Exchanging vector tag between processors failed" );
        }
        context.timer_pop( context.num_max_exchange );
        elapsed_times[3] = context.last_elapsed();
 
        // let us write out the local mesh after tag_exchange is called
        // we expect to see real data on both owned and ghost entities in halo regions (non-default values)
        if( context.debug_output && ( context.proc_id == 0 ) )  // only on root process, for debugging
        {
            dbgprint( "> Writing to file *after* ghost exchange " );
            runchk( context.moab_interface.write_file( "exchangeHalos_output_rank0_post.h5m", "H5M", "" ),
                    "Writing to disk failed" );
        }
 
        // Write out the final mesh with the tag data and mesh -- just for verification
        if( context.debug_output )
        {
            dbgprint( "> Writing out the final mesh and data in MOAB h5m format. File = " << context.output_filename );
            string write_options = ( context.num_procs > 1 ? "PARALLEL=WRITE_PART;DEBUG_IO=0;" : "" );
            // Write out to output file to visualize reduction/exchange of tag data
            runchk( context.moab_interface.write_file( context.output_filename.c_str(), "H5M", write_options.c_str() ),
                    "File write failed" );
        }
 
        // Consolidated timing results: the data is listed as follows
        // [ntasks,  nghosts,  load_mesh(I/O),  exchange_ghost_cells(setup), exchange_tags(scalar),
        // exchange_tags(vector)]
        dbgprint( "\n> Consolidated: [" << context.num_procs << ", " << context.ghost_layers << ", " << elapsed_times[0]
                                        << ", " << elapsed_times[1] << ", " << elapsed_times[2] << ", "
                                        << elapsed_times[3] << "]," );
 
        // execution finished
        dbgprint( "\n********** ExchangeHalos Example DONE! **********" );
    }
    // Done, cleanup
    MPI_Finalize();
 
    return 0;
}
 
/// Implementation details ///
 
///
RuntimeContext::RuntimeContext( MPI_Comm comm )
    : input_filename( std::string( MESH_DIR ) + std::string( "/io/mpasx1.642.t.2.nc" ) ),
      output_filename( "exchangeHalos_output.h5m" ), scalar_tagname( "scalar_variable" ),
      vector_tagname( "vector_variable" )
{
    // Create sets for the mesh and partition.  Then pass these to the load_file functions to populate the mesh.
    runchk_cont( moab_interface.create_meshset( moab::MESHSET_SET, fileset ), "Creating root set failed" );
    runchk_cont( moab_interface.create_meshset( moab::MESHSET_SET, partnset ), "Creating partition set failed" );
 
    // Create the parallel communicator object with the partition handle associated with MOAB
    parallel_communicator = moab::ParallelComm::get_pcomm( &moab_interface, partnset, &comm );
 
    proc_id   = parallel_communicator->rank();
    num_procs = parallel_communicator->size();
}
 
RuntimeContext::~RuntimeContext()
{
    delete parallel_communicator;
}
 
void RuntimeContext::ParseCLOptions( int argc, char* argv[] )
{
    ProgOptions opts;
    // Input mesh
    opts.addOpt< std::string >( "input", "Input mesh filename to load in parallel. Default=data/default_mesh_holes.h5m",
                                &input_filename );
    // Output mesh
    opts.addOpt< void >( "debug", "Should we write output file? Default=false", &debug_output );
    opts.addOpt< std::string >( "output",
                                "Output mesh filename for verification (use --debug). Default=exchangeHalos_output.h5m",
                                &output_filename );
    // Dimension of the input mesh
    // Vector tag length
    opts.addOpt< int >( "vtaglength", "Size of vector components per each entity. Default=3", &vector_length );
    // Number of halo (ghost) regions
    opts.addOpt< int >( "nghosts", "Number of ghost layers (halos) to exchange. Default=3", &ghost_layers );
    // Number of times to perform the halo exchange for timing
    opts.addOpt< int >( "nexchanges", "Number of ghost-halo exchange iterations to perform. Default=10",
                        &num_max_exchange );
 
    opts.parseCommandLine( argc, argv );
}
 
void RuntimeContext::timer_push( const std::string& operation )
{
    mTimerOps = mTimer.time_since_birth();
    mOpName   = operation;
}
 
void RuntimeContext::timer_pop( const int nruns )
{
    double locElapsed = mTimer.time_since_birth() - mTimerOps;
    double avgElapsed = 0;
    double maxElapsed = 0;
    MPI_Reduce( &locElapsed, &maxElapsed, 1, MPI_DOUBLE, MPI_MAX, 0, parallel_communicator->comm() );
    MPI_Reduce( &locElapsed, &avgElapsed, 1, MPI_DOUBLE, MPI_SUM, 0, parallel_communicator->comm() );
    if( proc_id == 0 )
    {
        avgElapsed /= num_procs;
        if( nruns > 1 )
            std::cout << "[LOG] Time taken to " << mOpName.c_str() << ", averaged over " << nruns
                      << " runs : max = " << maxElapsed / nruns << ", avg = " << avgElapsed / nruns << "\n";
        else
            std::cout << "[LOG] Time taken to " << mOpName.c_str() << " : max = " << maxElapsed
                      << ", avg = " << avgElapsed << "\n";
 
        last_counter = maxElapsed / nruns;
    }
    mOpName.clear();
}
 
double RuntimeContext::last_elapsed() const
{
    return last_counter;
}
 
moab::ErrorCode RuntimeContext::create_sv_tags( moab::Tag& tagScalar, moab::Tag& tagVector, moab::Range& entities )
{
    // Get element (centroid) coordinates so that we can evaluate some arbitrary data
    std::vector< double > entCoords = compute_centroids( entities );  // [entities * [lon, lat]]
 
    if( proc_id == 0 ) std::cout << "> Getting scalar tag handle " << scalar_tagname << "..." << std::endl;
    double defSTagValue = -1.0;
    bool createdTScalar = false;
    // Get or create the scalar tag: default name = "scalar_variable"
    // Type: double, Components: 1, Layout: Dense (all entities potentially), Default: -1.0
    runchk( moab_interface.tag_get_handle( scalar_tagname.c_str(), 1, moab::MB_TYPE_DOUBLE, tagScalar,
                                           moab::MB_TAG_CREAT | moab::MB_TAG_DENSE, &defSTagValue, &createdTScalar ),
            "Retrieving scalar tag handle failed" );
 
    // we expect to create a new tag -- fail if Tag already exists since we do not want to overwrite data
    assert( createdTScalar );
    // set the data for scalar tag with an analytical Spherical Harmonic function
    {
        std::vector< double > tagValues( entities.size(), -1.0 );
        std::generate( tagValues.begin(), tagValues.end(), [=, &entCoords]() {
            static int index = 0;
            const int offset = index++ * 2;
            return evaluate_function( entCoords[offset], entCoords[offset + 1] );
        } );
        // Set local scalar tag data for exchange
        runchk( moab_interface.tag_set_data( tagScalar, entities, tagValues.data() ),
                "Setting scalar tag data failed" );
    }
 
    if( proc_id == 0 ) std::cout << "> Getting vector tag handle " << vector_tagname << "..." << std::endl;
    std::vector< double > defVTagValue( vector_length, -1.0 );
    bool createdTVector = false;
    // Get or create the scalar tag: default name = "vector_variable"
    // Type: double, Components: vector_length, Layout: Dense (all entities potentially), Default: [-1.0,..]
    runchk( moab_interface.tag_get_handle( vector_tagname.c_str(), vector_length, moab::MB_TYPE_DOUBLE, tagVector,
                                           moab::MB_TAG_CREAT | moab::MB_TAG_DENSE, defVTagValue.data(),
                                           &createdTVector ),
            "Retrieving vector tag handle failed" );
 
    // we expect to create a new tag -- fail if Tag already exists since we do not want to overwrite data
    assert( createdTVector );
    // set the data for vector tag with an analytical Spherical Harmonic function
    // with an optional scaling for each component; just to make it look different :-)
    {
        const int veclength = vector_length;
        std::vector< double > tagValues( entities.size() * veclength, -1.0 );
        std::generate( tagValues.begin(), tagValues.end(), [=, &entCoords]() {
            static int index = 0;
            const int offset = ( index++ / veclength ) * 2;
            return this->evaluate_function( entCoords[offset], entCoords[offset + 1], 2, ( index % veclength + 1.0 ) );
        } );
        // Set local tag data for exchange
        runchk( moab_interface.tag_set_data( tagVector, entities, tagValues.data() ),
                "Setting vector tag data failed" );
    }
 
    return moab::MB_SUCCESS;
}
 
moab::ErrorCode RuntimeContext::load_file( bool load_ghosts )
{
    /// Parallel Read options:
    ///   PARALLEL = type {READ_PART} : Read on all tasks
    ///   PARTITION_METHOD = RCBZOLTAN : Use Zoltan partitioner to compute an online partition and redistribute on the
    ///   fly PARTITION = PARALLEL_PARTITION : Partition as you read based on part information stored in h5m file
    ///   PARALLEL_RESOLVE_SHARED_ENTS : Communicate to all processors to get the shared adjacencies
    ///   consistently in parallel
    ///   PARALLEL_GHOSTS : a.b.c
    ///                   : a = 2 - highest dimension of entities (2D in this case)
    ///                   : b = 1 - dimension of entities to calculate adjacencies (vertex=0, edges=1)
    ///                   : c = 3 - number of ghost layers needed (3 in this case)
    std::string read_options   = "DEBUG_IO=0;";
    std::string::size_type idx = input_filename.rfind( '.' );
    if( num_procs > 1 && idx != std::string::npos )
    {
        std::string extension = input_filename.substr( idx + 1 );
        if( !extension.compare( "nc" ) )
            // PARTITION_METHOD= [RCBZOLTAN, TRIVIAL]
            read_options += "PARALLEL=READ_PART;PARTITION_METHOD=RCBZOLTAN;"
                            "PARALLEL_RESOLVE_SHARED_ENTS;NO_EDGES;NO_MIXED_ELEMENTS;VARIABLE=;";
        else if( !extension.compare( "h5m" ) )
            read_options +=
                "PARALLEL=READ_PART;PARTITION=PARALLEL_PARTITION;"
                "PARALLEL_RESOLVE_SHARED_ENTS;" +
                ( load_ghosts ? "PARALLEL_THIN_GHOST_LAYER;PARALLEL_GHOSTS=2.1." + std::to_string( ghost_layers ) + ";"
                              : "" );
        else
        {
            std::cout << "Error unsupported file type (only h5m and nc) for this example: " << input_filename
                      << std::endl;
            return moab::MB_UNSUPPORTED_OPERATION;
        }
    }
 
    // Load the file from disk with given read options in parallel and associate all entities to fileset
    return moab_interface.load_file( input_filename.c_str(), &fileset, read_options.c_str() );
}
 
std::vector< double > RuntimeContext::compute_centroids( const moab::Range& entities ) const
{
    double node[3];
    std::vector< double > eCentroids( entities.size() * 2 );  // [lon, lat]
    size_t offset = 0;
    for( auto entity : entities )
    {
        // Get the element coordinates (centroid) on the real mesh
        runchk_cont( moab_interface.get_coords( &entity, 1, node ), "Getting entity coordinates failed" );
 
        // scale by magnitude so that element is on unit sphere
        double magnitude = std::sqrt( node[0] * node[0] + node[1] * node[1] + node[2] * node[2] );
        node[0] /= magnitude;
        node[1] /= magnitude;
        node[2] /= magnitude;
 
        // compute the spherical transformation onto unit sphere
        eCentroids[offset] = atan2( node[1], node[0] );
        if( eCentroids[offset] < 0.0 ) eCentroids[offset] += 2.0 * M_PI;
        eCentroids[offset + 1] = asin( node[2] );
 
        offset += 2;  // increment the offset
    }
    // return centroid list for elements
    return eCentroids;
}