mbtempest.cpp

Go to the documentation of this file.

/*
 * Usage: MOAB-Tempest tool
 *
 * Generate a Cubed-Sphere mesh: ./mbtempest -t 0 -res 25 -f cubed_sphere_mesh.exo
 * Generate a RLL mesh: ./mbtempest -t 1 -res 25 -f rll_mesh.exo
 * Generate a Icosahedral-Sphere mesh: ./mbtempest -t 2 -res 25 <-dual> -f icosa_mesh.exo
 *
 * Now you can compute the intersections between the meshes too!
 *
 * Generate the overlap mesh: ./mbtempest -t 3 -l cubed_sphere_mesh.exo -l rll_mesh.exo -f
 * overlap_mesh.exo
 *
 */
 
// standard C++ includes
#include <iostream>
#include <iomanip>
#include <cstdlib>
#include <vector>
#include <string>
#include <sstream>
#include <cassert>
 
// MOAB includes
#include "moab/Core.hpp"
#include "moab/IntxMesh/IntxUtils.hpp"
#include "moab/Remapping/TempestRemapper.hpp"
#include "moab/Remapping/TempestOnlineMap.hpp"
#include "moab/ProgOptions.hpp"
#include "moab/CpuTimer.hpp"
#include "DebugOutput.hpp"
 
//#ifndef MOAB_HAVE_MPI
// #error mbtempest tool requires MPI configuration
//#endif
 
#ifdef MOAB_HAVE_MPI
// MPI includes
#include "moab_mpi.h"
#include "moab/ParallelComm.hpp"
#include "MBParallelConventions.h"
#endif
 
struct ToolContext
{
 moab::Core* mbcore;
#ifdef MOAB_HAVE_MPI
 moab::ParallelComm* pcomm;
#endif
 const int proc_id, n_procs;
 moab::DebugOutput outputFormatter;
 int blockSize;
 std::vector< std::string > inFilenames;
 std::vector< Mesh* > meshes;
 std::vector< moab::EntityHandle > meshsets;
 std::vector< int > disc_orders;
 std::vector< std::string > disc_methods;
 std::vector< std::string > doftag_names;
 std::string fvMethod;
 std::string outFilename;
 std::string intxFilename;
 std::string baselineFile;
 moab::TempestRemapper::TempestMeshType meshType;
 bool computeDual;
 bool computeWeights;
 bool verifyWeights;
 bool enforceConvexity;
 int ensureMonotonicity;
 bool rrmGrids;
 bool kdtreeSearch;
 bool fCheck;
 bool fVolumetric;
 bool useGnomonicProjection;
 GenerateOfflineMapAlgorithmOptions mapOptions;
 
#ifdef MOAB_HAVE_MPI
 ToolContext( moab::Core* icore, moab::ParallelComm* p_pcomm )
 : mbcore( icore ), pcomm( p_pcomm ), proc_id( pcomm->rank() ), n_procs( pcomm->size() ),
 outputFormatter( std::cout, pcomm->rank(), 0 ),
#else
 ToolContext( moab::Core* icore )
 : mbcore( icore ), proc_id( 0 ), n_procs( 1 ), outputFormatter( std::cout, 0, 0 ),
#endif
 blockSize( 5 ), fvMethod( "none" ), outFilename( "outputFile.nc" ), intxFilename( "intxFile.h5m" ),
 baselineFile( "" ), meshType( moab::TempestRemapper::DEFAULT ), computeDual( false ), computeWeights( false ),
 verifyWeights( false ), enforceConvexity( false ), ensureMonotonicity( 0 ), rrmGrids( false ),
 kdtreeSearch( true ), fCheck( n_procs > 1 ? false : true ), fVolumetric( false ),
 useGnomonicProjection( false )
 {
 inFilenames.resize( 2 );
 doftag_names.resize( 2 );
 timer = new moab::CpuTimer();
 
 outputFormatter.set_prefix( "[MBTempest]: " );
 }
 
 ~ToolContext()
 {
 // for (unsigned i=0; i < meshes.size(); ++i) delete meshes[i];
 meshes.clear();
 inFilenames.clear();
 disc_orders.clear();
 disc_methods.clear();
 doftag_names.clear();
 outFilename.clear();
 intxFilename.clear();
 baselineFile.clear();
 meshsets.clear();
 delete timer;
 }
 
 void timer_push( std::string operation )
 {
 timer_ops = timer->time_since_birth();
 opName = operation;
 }
 
 void timer_pop()
 {
 double locElapsed = timer->time_since_birth() - timer_ops, avgElapsed = 0, maxElapsed = 0;
#ifdef MOAB_HAVE_MPI
 MPI_Reduce( &locElapsed, &maxElapsed, 1, MPI_DOUBLE, MPI_MAX, 0, pcomm->comm() );
 MPI_Reduce( &locElapsed, &avgElapsed, 1, MPI_DOUBLE, MPI_SUM, 0, pcomm->comm() );
#else
 maxElapsed = locElapsed;
 avgElapsed = locElapsed;
#endif
 if( !proc_id )
 {
 avgElapsed /= n_procs;
 std::cout << "[LOG] Time taken to " << opName.c_str() << ": max = " << maxElapsed
 << ", avg = " << avgElapsed << "\n";
 }
 // std::cout << "\n[LOG" << proc_id << "] Time taken to " << opName << " = " <<
 // timer->time_since_birth() - timer_ops << std::endl;
 opName.clear();
 }
 
 void ParseCLOptions( int argc, char* argv[] )
 {
 ProgOptions opts;
 int imeshType = 0;
 std::string expectedFName = "output.exo";
 std::string expectedMethod = "fv";
 std::string expectedFVMethod = "none";
 std::string expectedDofTagName = "GLOBAL_ID";
 int expectedOrder = 1;
 
 if( !proc_id )
 {
 std::cout << "Command line options provided to mbtempest:\n ";
 for( int iarg = 0; iarg < argc; ++iarg )
 std::cout << argv[iarg] << " ";
 std::cout << std::endl << std::endl;
 }
 
 opts.addOpt< int >( "type,t",
 "Type of mesh (default=CS; Choose from [CS=0, RLL=1, ICO=2, OVERLAP_FILES=3, "
 "OVERLAP_MEMORY=4, OVERLAP_MOAB=5])",
 &imeshType );
 
 opts.addOpt< int >( "res,r", "Resolution of the mesh (default=5)", &blockSize );
 
 opts.addOpt< void >( "dual,d", "Output the dual of the mesh (relevant only for ICO mesh type)", &computeDual );
 
 opts.addOpt< std::string >( "file,f", "Output computed mesh or remapping weights to specified filename",
 &outFilename );
 opts.addOpt< std::string >(
 "load,l", "Input mesh filenames for source and target meshes. (relevant only when computing weights)",
 &expectedFName );
 
 opts.addOpt< void >( "advfront,a",
 "Use the advancing front intersection instead of the Kd-tree based algorithm "
 "to compute mesh intersections. (relevant only when computing weights)" );
 
 opts.addOpt< std::string >( "intx,i", "Output TempestRemap intersection mesh filename", &intxFilename );
 
 opts.addOpt< void >( "weights,w",
 "Compute and output the weights using the overlap mesh (generally "
 "relevant only for OVERLAP mesh)",
 &computeWeights );
 
 opts.addOpt< std::string >( "method,m", "Discretization method for the source and target solution fields",
 &expectedMethod );
 
 opts.addOpt< int >( "order,o", "Discretization orders for the source and target solution fields",
 &expectedOrder );
 
 opts.addOpt< std::string >( "global_id,g",
 "Tag name that contains the global DoF IDs for source and target solution fields",
 &expectedDofTagName );
 
 opts.addOpt< void >( "noconserve",
 "Do not apply conservation to the resultant weights (relevant only "
 "when computing weights)",
 &mapOptions.fNoConservation );
 
 opts.addOpt< void >( "volumetric",
 "Apply a volumetric projection to compute the weights (relevant only "
 "when computing weights)",
 &fVolumetric );
 
 opts.addOpt< int >( "monotonicity", "Ensure monotonicity in the weight generation. Options=[0,1,2,3]",
 &ensureMonotonicity );
 
 opts.addOpt< void >( "gnomonic", "Use Gnomonic plane projections to compute coverage mesh.",
 &useGnomonicProjection );
 
 opts.addOpt< std::string >( "fvmethod",
 "Sub-type method for FV-FV projections (invdist, delaunay, bilin, "
 "intbilin, intbilingb, none. Default: none)",
 &expectedFVMethod );
 
 opts.addOpt< void >( "enforce_convexity", "check convexity of input meshes to compute mesh intersections",
 &enforceConvexity );
 
 opts.addOpt< void >( "nobubble", "do not use bubble on interior of spectral element nodes",
 &mapOptions.fNoBubble );
 
 opts.addOpt< void >(
 "sparseconstraints",
 "Use sparse solver for constraints when we have high-valence (typical with high-res RLL mesh)",
 &mapOptions.fSparseConstraints );
 
 opts.addOpt< void >( "rrmgrids",
 "At least one of the meshes is a regionally refined grid (relevant to "
 "accelerate intersection computation)",
 &rrmGrids );
 
 opts.addOpt< void >( "checkmap", "Check the generated map for conservation and consistency", &fCheck );
 
 opts.addOpt< void >( "verify",
 "Verify the accuracy of the maps by projecting analytical functions "
 "from source to target "
 "grid by applying the maps",
 &verifyWeights );
 
 opts.addOpt< std::string >( "baseline", "Output baseline file", &baselineFile );
 
 opts.parseCommandLine( argc, argv );
 
 // By default - use Kd-tree based search; if user asks for advancing front, disable Kd-tree
 // algorithm
 kdtreeSearch = opts.numOptSet( "advfront,a" ) == 0;
 
 switch( imeshType )
 {
 case 0:
 meshType = moab::TempestRemapper::CS;
 break;
 
 case 1:
 meshType = moab::TempestRemapper::RLL;
 break;
 
 case 2:
 meshType = moab::TempestRemapper::ICO;
 break;
 
 case 3:
 meshType = moab::TempestRemapper::OVERLAP_FILES;
 break;
 
 case 4:
 meshType = moab::TempestRemapper::OVERLAP_MEMORY;
 break;
 
 case 5:
 meshType = moab::TempestRemapper::OVERLAP_MOAB;
 break;
 
 default:
 meshType = moab::TempestRemapper::DEFAULT;
 break;
 }
 
 if( meshType > moab::TempestRemapper::ICO ) // compute overlap mesh and maps possibly
 {
 opts.getOptAllArgs( "load,l", inFilenames );
 opts.getOptAllArgs( "order,o", disc_orders );
 opts.getOptAllArgs( "method,m", disc_methods );
 opts.getOptAllArgs( "global_id,i", doftag_names );
 
 assert( inFilenames.size() == 2 );
 assert( disc_orders.size() == 2 );
 assert( disc_methods.size() == 2 );
 assert( ensureMonotonicity >= 0 && ensureMonotonicity <= 3 );
 
 // get discretization order parameters
 if( disc_orders.size() == 0 ) disc_orders.resize( 2, 1 );
 if( disc_orders.size() == 1 ) disc_orders.push_back( 1 );
 
 // get discretization method parameters
 if( disc_methods.size() == 0 ) disc_methods.resize( 2, "fv" );
 if( disc_methods.size() == 1 ) disc_methods.push_back( "fv" );
 
 // get DoF tagname parameters
 if( doftag_names.size() == 0 ) doftag_names.resize( 2, "GLOBAL_ID" );
 if( doftag_names.size() == 1 ) doftag_names.push_back( "GLOBAL_ID" );
 
 // for computing maps and overlaps, set discretization orders
 mapOptions.nPin = disc_orders[0];
 mapOptions.nPout = disc_orders[1];
 mapOptions.fSourceConcave = false;
 mapOptions.fTargetConcave = false;
 
 mapOptions.strMethod = "";
 
 if( expectedFVMethod != "none" )
 {
 mapOptions.strMethod += expectedFVMethod + ";";
 fvMethod = expectedFVMethod;
 
 // These FV projection methods are non-conservative; specify it explicitly
 mapOptions.fNoConservation = true;
 }
 switch( ensureMonotonicity )
 {
 case 0:
 mapOptions.fMonotone = false;
 break;
 case 3:
 mapOptions.strMethod += "mono3;";
 break;
 case 2:
 mapOptions.strMethod += "mono2;";
 break;
 default:
 mapOptions.fMonotone = true;
 }
 mapOptions.fNoCorrectAreas = false;
 mapOptions.fNoCheck = !fCheck;
 
 //assert( fVolumetric && fInverseDistanceMap == false ); // both options cannot be active
 if( fVolumetric ) mapOptions.strMethod += "volumetric;";
 }
 
 // clear temporary string name
 expectedFName.clear();
 
 mapOptions.strOutputMapFile = outFilename;
 mapOptions.strOutputFormat = "Netcdf4";
 }
 
 private:
 moab::CpuTimer* timer;
 double timer_ops;
 std::string opName;
};
 
// Forward declare some methods
static moab::ErrorCode CreateTempestMesh( ToolContext&, moab::TempestRemapper& remapper, Mesh* );
inline double sample_slow_harmonic( double dLon, double dLat );
inline double sample_fast_harmonic( double dLon, double dLat );
inline double sample_constant( double dLon, double dLat );
inline double sample_stationary_vortex( double dLon, double dLat );
 
std::string get_file_read_options( ToolContext& ctx, std::string filename )
{
 // if running in serial, blank option may suffice in general
 std::string opts = "";
 if( ctx.n_procs > 1 )
 {
 size_t lastindex = filename.find_last_of( "." );
 std::string extension = filename.substr( lastindex + 1, filename.size() );
 if( extension == "h5m" )
 return "PARALLEL=READ_PART;PARTITION=PARALLEL_PARTITION;PARALLEL_RESOLVE_SHARED_ENTS;";
 else if( extension == "nc" )
 {
 // set default set of options
 opts = "PARALLEL=READ_PART;PARTITION_METHOD=RCBZOLTAN;";
 
 if( ctx.proc_id == 0 )
 {
 {
 NcFile ncFile( filename.c_str(), NcFile::ReadOnly );
 if( !ncFile.is_valid() )
 _EXCEPTION1( "Unable to open grid file \"%s\" for reading", filename.c_str() );
 
 // Check for dimension names "grid_size", "grid_rank" and "grid_corners"
 int iSCRIPFormat = 0, iMPASFormat = 0, iESMFFormat = 0;
 for( int i = 0; i < ncFile.num_dims(); i++ )
 {
 std::string strDimName = ncFile.get_dim( i )->name();
 if( strDimName == "grid_size" || strDimName == "grid_corners" || strDimName == "grid_rank" )
 iSCRIPFormat++;
 if( strDimName == "nodeCount" || strDimName == "elementCount" ||
 strDimName == "maxNodePElement" )
 iESMFFormat++;
 if( strDimName == "nCells" || strDimName == "nEdges" || strDimName == "nVertices" ||
 strDimName == "vertexDegree" )
 iMPASFormat++;
 }
 ncFile.close();
 
 if( iESMFFormat == 3 ) // Input from a NetCDF ESMF file
 {
 opts = "PARALLEL=READ_PART;PARTITION_METHOD=RCBZOLTAN;PARALLEL_RESOLVE_SHARED_ENTS;VARIABLE=;";
 }
 if( iSCRIPFormat == 3 ) // Input from a NetCDF SCRIP file
 {
 opts = "PARALLEL=READ_PART;PARTITION_METHOD=RCBZOLTAN;";
 }
 if( iMPASFormat == 4 ) // Input from a NetCDF SCRIP file
 {
 opts = "PARALLEL=READ_PART;PARTITION_METHOD=RCBZOLTAN;"
 "PARALLEL_RESOLVE_SHARED_ENTS;NO_EDGES;NO_MIXED_ELEMENTS;VARIABLE=;";
 }
 }
 }
 
#ifdef MOAB_HAVE_MPI
 int line_size = opts.size();
 MPI_Bcast( &line_size, 1, MPI_INT, 0, MPI_COMM_WORLD );
 if( ctx.proc_id != 0 ) opts.resize( line_size );
 MPI_Bcast( const_cast< char* >( opts.data() ), line_size, MPI_CHAR, 0, MPI_COMM_WORLD );
#endif
 }
 else
 return "PARALLEL=BCAST_DELETE;PARTITION=TRIVIAL;PARALLEL_RESOLVE_SHARED_ENTS;";
 }
 return opts;
}
 
int main( int argc, char* argv[] )
{
 moab::ErrorCode rval;
 NcError error( NcError::verbose_nonfatal );
 std::stringstream sstr;
 std::string historyStr;
 
 int proc_id = 0, nprocs = 1;
#ifdef MOAB_HAVE_MPI
 MPI_Init( &argc, &argv );
 MPI_Comm_rank( MPI_COMM_WORLD, &proc_id );
 MPI_Comm_size( MPI_COMM_WORLD, &nprocs );
#endif
 
 moab::Core* mbCore = new( std::nothrow ) moab::Core;
 
 if( NULL == mbCore )
 {
 return 1;
 }
 
 // Build the history string
 for( int ia = 0; ia < argc; ++ia )
 historyStr += std::string( argv[ia] ) + " ";
 
 ToolContext* runCtx;
#ifdef MOAB_HAVE_MPI
 moab::ParallelComm* pcomm = new moab::ParallelComm( mbCore, MPI_COMM_WORLD, 0 );
 
 runCtx = new ToolContext( mbCore, pcomm );
 const char* writeOptions = ( nprocs > 1 ? "PARALLEL=WRITE_PART" : "" );
#else
 runCtx = new ToolContext( mbCore );
 const char* writeOptions = "";
#endif
 runCtx->ParseCLOptions( argc, argv );
 
 const double radius_src = 1.0 /*2.0*acos(-1.0)*/;
 const double radius_dest = 1.0 /*2.0*acos(-1.0)*/;
 
 moab::DebugOutput& outputFormatter = runCtx->outputFormatter;
 
#ifdef MOAB_HAVE_MPI
 moab::TempestRemapper remapper( mbCore, pcomm );
#else
 moab::TempestRemapper remapper( mbCore );
#endif
 remapper.meshValidate = true;
 remapper.constructEdgeMap = true;
 remapper.initialize();
 
 // Default area_method = lHuiller; Options: Girard, GaussQuadrature (if TR is available)
#ifdef MOAB_HAVE_TEMPESTREMAP
 moab::IntxAreaUtils areaAdaptor( moab::IntxAreaUtils::GaussQuadrature );
#else
 moab::IntxAreaUtils areaAdaptor( moab::IntxAreaUtils::lHuiller );
#endif
 
 Mesh* tempest_mesh = new Mesh();
 runCtx->timer_push( "create Tempest mesh" );
 rval = CreateTempestMesh( *runCtx, remapper, tempest_mesh );MB_CHK_ERR( rval );
 runCtx->timer_pop();
 
 const double epsrel = ReferenceTolerance; // ReferenceTolerance is defined in Defines.h in tempestremap;
 // Defines.h is included in SparseMatrix.h
 // SparseMatrix.h is included in OfflineMap.h
 // OfflineMap.h is included in TempestOnlineMap.hpp
 // TempestOnlineMap.hpp is included in this file, and is part of MOAB
 // Some constant parameters
 
 const double boxeps = 1e-1;
 
 if( runCtx->meshType == moab::TempestRemapper::OVERLAP_MEMORY )
 {
 // Compute intersections with MOAB
 // For the overlap method, choose between: "fuzzy", "exact" or "mixed"
 assert( runCtx->meshes.size() == 3 );
 
#ifdef MOAB_HAVE_MPI
 rval = pcomm->check_all_shared_handles();MB_CHK_ERR( rval );
#endif
 
 // Load the meshes and validate
 rval = remapper.ConvertTempestMesh( moab::Remapper::SourceMesh );MB_CHK_ERR( rval );
 rval = remapper.ConvertTempestMesh( moab::Remapper::TargetMesh );MB_CHK_ERR( rval );
 rval = remapper.ConvertTempestMesh( moab::Remapper::OverlapMesh );MB_CHK_ERR( rval );
 rval = mbCore->write_mesh( "tempest_intersection.h5m", &runCtx->meshsets[2], 1 );MB_CHK_ERR( rval );
 
 // print verbosely about the problem setting
 {
 moab::Range rintxverts, rintxelems;
 rval = mbCore->get_entities_by_dimension( runCtx->meshsets[0], 0, rintxverts );MB_CHK_ERR( rval );
 rval = mbCore->get_entities_by_dimension( runCtx->meshsets[0], 2, rintxelems );MB_CHK_ERR( rval );
 outputFormatter.printf( 0, "The source set contains %lu vertices and %lu elements \n", rintxverts.size(),
 rintxelems.size() );
 
 moab::Range bintxverts, bintxelems;
 rval = mbCore->get_entities_by_dimension( runCtx->meshsets[1], 0, bintxverts );MB_CHK_ERR( rval );
 rval = mbCore->get_entities_by_dimension( runCtx->meshsets[1], 2, bintxelems );MB_CHK_ERR( rval );
 outputFormatter.printf( 0, "The target set contains %lu vertices and %lu elements \n", bintxverts.size(),
 bintxelems.size() );
 }
 
 moab::EntityHandle intxset; // == remapper.GetMeshSet(moab::Remapper::OverlapMesh);
 
 // Compute intersections with MOAB
 {
 // Create the intersection on the sphere object
 runCtx->timer_push( "setup the intersector" );
 
 moab::Intx2MeshOnSphere* mbintx = new moab::Intx2MeshOnSphere( mbCore );
 mbintx->set_error_tolerance( epsrel );
 mbintx->set_box_error( boxeps );
 mbintx->set_radius_source_mesh( radius_src );
 mbintx->set_radius_destination_mesh( radius_dest );
#ifdef MOAB_HAVE_MPI
 mbintx->set_parallel_comm( pcomm );
#endif
 rval = mbintx->FindMaxEdges( runCtx->meshsets[0], runCtx->meshsets[1] );MB_CHK_ERR( rval );
 
#ifdef MOAB_HAVE_MPI
 moab::Range local_verts;
 rval = mbintx->build_processor_euler_boxes( runCtx->meshsets[1], local_verts );MB_CHK_ERR( rval );
 
 runCtx->timer_pop();
 
 moab::EntityHandle covering_set;
 runCtx->timer_push( "communicate the mesh" );
 rval = mbintx->construct_covering_set( runCtx->meshsets[0], covering_set );MB_CHK_ERR( rval ); // lots of communication if mesh is distributed very differently
 runCtx->timer_pop();
#else
 moab::EntityHandle covering_set = runCtx->meshsets[0];
#endif
 // Now let's invoke the MOAB intersection algorithm in parallel with a
 // source and target mesh set representing two different decompositions
 runCtx->timer_push( "compute intersections with MOAB" );
 rval = mbCore->create_meshset( moab::MESHSET_SET, intxset );MB_CHK_SET_ERR( rval, "Can't create new set" );
 rval = mbintx->intersect_meshes( covering_set, runCtx->meshsets[1], intxset );MB_CHK_SET_ERR( rval, "Can't compute the intersection of meshes on the sphere" );
 runCtx->timer_pop();
 
 // free the memory
 delete mbintx;
 }
 
 {
 moab::Range intxelems, intxverts;
 rval = mbCore->get_entities_by_dimension( intxset, 2, intxelems );MB_CHK_ERR( rval );
 rval = mbCore->get_entities_by_dimension( intxset, 0, intxverts, true );MB_CHK_ERR( rval );
 outputFormatter.printf( 0, "The intersection set contains %lu elements and %lu vertices \n",
 intxelems.size(), intxverts.size() );
 
 double initial_sarea =
 areaAdaptor.area_on_sphere( mbCore, runCtx->meshsets[0],
 radius_src ); // use the target to compute the initial area
 double initial_tarea =
 areaAdaptor.area_on_sphere( mbCore, runCtx->meshsets[1],
 radius_dest ); // use the target to compute the initial area
 double intx_area = areaAdaptor.area_on_sphere( mbCore, intxset, radius_src );
 
 outputFormatter.printf( 0, "mesh areas: source = %12.10f, target = %12.10f, intersection = %12.10f \n",
 initial_sarea, initial_tarea, intx_area );
 outputFormatter.printf( 0, "relative error w.r.t source = %12.10e, target = %12.10e \n",
 fabs( intx_area - initial_sarea ) / initial_sarea,
 fabs( intx_area - initial_tarea ) / initial_tarea );
 }
 
 // Write out our computed intersection file
 rval = mbCore->write_mesh( "moab_intersection.h5m", &intxset, 1 );MB_CHK_ERR( rval );
 
 if( runCtx->computeWeights )
 {
 runCtx->timer_push( "compute weights with the Tempest meshes" );
 // Call to generate an offline map with the tempest meshes
 OfflineMap weightMap;
 int err = GenerateOfflineMapWithMeshes( *runCtx->meshes[0], *runCtx->meshes[1], *runCtx->meshes[2],
 runCtx->disc_methods[0], // std::string strInputType
 runCtx->disc_methods[1], // std::string strOutputType,
 runCtx->mapOptions, weightMap );
 runCtx->timer_pop();
 
 std::map< std::string, std::string > mapAttributes;
 if( err )
 {
 rval = moab::MB_FAILURE;
 }
 else
 {
 weightMap.Write( "outWeights.nc", mapAttributes );
 }
 }
 }
 else if( runCtx->meshType == moab::TempestRemapper::OVERLAP_MOAB )
 {
 // Usage: mpiexec -n 2 tools/mbtempest -t 5 -l mycs_2.h5m -l myico_2.h5m -f myoverlap_2.h5m
#ifdef MOAB_HAVE_MPI
 rval = pcomm->check_all_shared_handles();MB_CHK_ERR( rval );
#endif
 // print verbosely about the problem setting
 {
 moab::Range srcverts, srcelems;
 rval = mbCore->get_entities_by_dimension( runCtx->meshsets[0], 0, srcverts );MB_CHK_ERR( rval );
 rval = mbCore->get_entities_by_dimension( runCtx->meshsets[0], 2, srcelems );MB_CHK_ERR( rval );
 rval = moab::IntxUtils::fix_degenerate_quads( mbCore, runCtx->meshsets[0] );MB_CHK_ERR( rval );
 if( runCtx->enforceConvexity )
 {
 rval = moab::IntxUtils::enforce_convexity( mbCore, runCtx->meshsets[0], proc_id );MB_CHK_ERR( rval );
 }
 rval = areaAdaptor.positive_orientation( mbCore, runCtx->meshsets[0], radius_src );MB_CHK_ERR( rval );
 if( !proc_id )
 outputFormatter.printf( 0, "The source set contains %lu vertices and %lu elements \n", srcverts.size(),
 srcelems.size() );
 
 moab::Range tgtverts, tgtelems;
 rval = mbCore->get_entities_by_dimension( runCtx->meshsets[1], 0, tgtverts );MB_CHK_ERR( rval );
 rval = mbCore->get_entities_by_dimension( runCtx->meshsets[1], 2, tgtelems );MB_CHK_ERR( rval );
 rval = moab::IntxUtils::fix_degenerate_quads( mbCore, runCtx->meshsets[1] );MB_CHK_ERR( rval );
 if( runCtx->enforceConvexity )
 {
 rval = moab::IntxUtils::enforce_convexity( mbCore, runCtx->meshsets[1], proc_id );MB_CHK_ERR( rval );
 }
 rval = areaAdaptor.positive_orientation( mbCore, runCtx->meshsets[1], radius_dest );MB_CHK_ERR( rval );
 if( !proc_id )
 outputFormatter.printf( 0, "The target set contains %lu vertices and %lu elements \n", tgtverts.size(),
 tgtelems.size() );
 }
 //rval = mbCore->write_file( "source_mesh.h5m", NULL, writeOptions, &runCtx->meshsets[0], 1 );MB_CHK_ERR( rval );
 //rval = mbCore->write_file( "target_mesh.h5m", NULL, writeOptions, &runCtx->meshsets[1], 1 );MB_CHK_ERR( rval );
 
 // First compute the covering set such that the target elements are fully covered by the
 // lcoal source grid
 runCtx->timer_push( "construct covering set for intersection" );
 rval =
 remapper.ConstructCoveringSet( epsrel, 1.0, 1.0, boxeps, runCtx->rrmGrids, runCtx->useGnomonicProjection );MB_CHK_ERR( rval );
 runCtx->timer_pop();
 
 // Compute intersections with MOAB with either the Kd-tree or the advancing front algorithm
 runCtx->timer_push( "setup and compute mesh intersections" );
 rval = remapper.ComputeOverlapMesh( runCtx->kdtreeSearch, false );MB_CHK_ERR( rval );
 runCtx->timer_pop();
 
 // print some diagnostic checks to see if the overlap grid resolved the input meshes
 // correctly
 {
 double local_areas[3],
 global_areas[3]; // Array for Initial area, and through Method 1 and Method 2
 // local_areas[0] = area_on_sphere_lHuiller ( mbCore, runCtx->meshsets[1], radius_src );
 local_areas[0] = areaAdaptor.area_on_sphere( mbCore, runCtx->meshsets[0], radius_src );
 local_areas[1] = areaAdaptor.area_on_sphere( mbCore, runCtx->meshsets[1], radius_dest );
 local_areas[2] = areaAdaptor.area_on_sphere( mbCore, runCtx->meshsets[2], radius_src );
 
#ifdef MOAB_HAVE_MPI
 MPI_Allreduce( &local_areas[0], &global_areas[0], 3, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD );
#else
 global_areas[0] = local_areas[0];
 global_areas[1] = local_areas[1];
 global_areas[2] = local_areas[2];
#endif
 if( !proc_id )
 {
 outputFormatter.printf( 0,
 "initial area: source mesh = %12.14f, target mesh = "
 "%12.14f, overlap mesh = %12.14f\n",
 global_areas[0], global_areas[1], global_areas[2] );
 // outputFormatter.printf ( 0, " area with l'Huiller: %12.14f with Girard:
 // %12.14f\n", global_areas[2], global_areas[3] ); outputFormatter.printf ( 0, "
 // relative difference areas = %12.10e\n", fabs ( global_areas[2] - global_areas[3]
 // ) / global_areas[2] );
 outputFormatter.printf( 0, "relative error w.r.t source = %12.14e, and target = %12.14e\n",
 fabs( global_areas[0] - global_areas[2] ) / global_areas[0],
 fabs( global_areas[1] - global_areas[2] ) / global_areas[1] );
 }
 }
 
 if( runCtx->intxFilename.size() )
 {
 moab::EntityHandle writableOverlapSet;
 rval = mbCore->create_meshset( moab::MESHSET_SET, writableOverlapSet );MB_CHK_SET_ERR( rval, "Can't create new set" );
 moab::EntityHandle meshOverlapSet = remapper.GetMeshSet( moab::Remapper::OverlapMesh );
 moab::Range ovEnts;
 rval = mbCore->get_entities_by_dimension( meshOverlapSet, 2, ovEnts );MB_CHK_SET_ERR( rval, "Can't create new set" );
 rval = mbCore->get_entities_by_dimension( meshOverlapSet, 0, ovEnts );MB_CHK_SET_ERR( rval, "Can't create new set" );
 
#ifdef MOAB_HAVE_MPI
 // Do not remove ghosted entities if we still haven't computed weights
 // Remove ghosted entities from overlap set before writing the new mesh set to file
 if( nprocs > 1 )
 {
 moab::Range ghostedEnts;
 rval = remapper.GetOverlapAugmentedEntities( ghostedEnts );MB_CHK_ERR( rval );
 ovEnts = moab::subtract( ovEnts, ghostedEnts );
 }
#endif
 rval = mbCore->add_entities( writableOverlapSet, ovEnts );MB_CHK_SET_ERR( rval, "Deleting ghosted entities failed" );
 
 size_t lastindex = runCtx->intxFilename.find_last_of( "." );
 sstr.str( "" );
 sstr << runCtx->intxFilename.substr( 0, lastindex ) << ".h5m";
 if( !runCtx->proc_id )
 std::cout << "Writing out the MOAB intersection mesh file to " << sstr.str() << std::endl;
 
 // Write out our computed intersection file
 rval = mbCore->write_file( sstr.str().c_str(), NULL, writeOptions, &writableOverlapSet, 1 );MB_CHK_ERR( rval );
 }
 
 if( runCtx->computeWeights )
 {
 runCtx->meshes[2] = remapper.GetMesh( moab::Remapper::OverlapMesh );
 if( !runCtx->proc_id ) std::cout << std::endl;
 
 runCtx->timer_push( "setup computation of weights" );
 // Call to generate the remapping weights with the tempest meshes
 moab::TempestOnlineMap* weightMap = new moab::TempestOnlineMap( &remapper );
 runCtx->timer_pop();
 
 runCtx->timer_push( "compute weights with TempestRemap" );
 rval = weightMap->GenerateRemappingWeights(
 runCtx->disc_methods[0], // std::string strInputType
 runCtx->disc_methods[1], // std::string strOutputType,
 runCtx->mapOptions, // const GenerateOfflineMapAlgorithmOptions& options
 runCtx->doftag_names[0], // const std::string& source_tag_name
 runCtx->doftag_names[1] // const std::string& target_tag_name
 );MB_CHK_ERR( rval );
 runCtx->timer_pop();
 
 // Invoke the CheckMap routine on the TempestRemap serial interface directly, if running
 // on a single process
 if( nprocs == 1 )
 {
 const double dNormalTolerance = 1.0E-8;
 const double dStrictTolerance = 1.0E-12;
 weightMap->CheckMap( runCtx->fCheck, runCtx->fCheck, runCtx->fCheck && ( runCtx->ensureMonotonicity ),
 dNormalTolerance, dStrictTolerance );
 }
 
 if( runCtx->outFilename.size() )
 {
 std::map< std::string, std::string > attrMap;
 attrMap["MOABversion"] = std::string( MOAB_VERSION );
 attrMap["Title"] = "MOAB-TempestRemap (mbtempest) Offline Regridding Weight Generator";
 attrMap["normalization"] = "ovarea";
 attrMap["remap_options"] = runCtx->mapOptions.strMethod;
 attrMap["domain_a"] = runCtx->inFilenames[0];
 attrMap["domain_b"] = runCtx->inFilenames[1];
 attrMap["domain_aUb"] = runCtx->intxFilename;
 attrMap["map_aPb"] = runCtx->outFilename;
 attrMap["methodorder_a"] = runCtx->disc_methods[0] + ":" + std::to_string( runCtx->disc_orders[0] ) +
 ":" + std::string( runCtx->doftag_names[0] );
 attrMap["concave_a"] = runCtx->mapOptions.fSourceConcave ? "true" : "false";
 attrMap["methodorder_b"] = runCtx->disc_methods[1] + ":" + std::to_string( runCtx->disc_orders[1] ) +
 ":" + std::string( runCtx->doftag_names[1] );
 attrMap["concave_b"] = runCtx->mapOptions.fTargetConcave ? "true" : "false";
 attrMap["bubble"] = runCtx->mapOptions.fNoBubble ? "false" : "true";
 attrMap["history"] = historyStr;
 
 // Write the map file to disk in parallel using either HDF5 or SCRIP interface
 rval = weightMap->WriteParallelMap( runCtx->outFilename.c_str(), attrMap );MB_CHK_ERR( rval );
 }
 
 if( runCtx->verifyWeights )
 {
 // Let us pick a sampling test function for solution evaluation
 moab::TempestOnlineMap::sample_function testFunction =
 &sample_stationary_vortex; // &sample_slow_harmonic;
 
 runCtx->timer_push( "describe a solution on source grid" );
 moab::Tag srcAnalyticalFunction;
 rval = weightMap->DefineAnalyticalSolution( srcAnalyticalFunction, "AnalyticalSolnSrcExact",
 moab::Remapper::SourceMesh, testFunction );MB_CHK_ERR( rval );
 runCtx->timer_pop();
 // rval = mbCore->write_file ( "srcWithSolnTag.h5m", NULL, writeOptions,
 // &runCtx->meshsets[0], 1 ); MB_CHK_ERR ( rval );
 
 runCtx->timer_push( "describe a solution on target grid" );
 moab::Tag tgtAnalyticalFunction;
 moab::Tag tgtProjectedFunction;
 rval = weightMap->DefineAnalyticalSolution( tgtAnalyticalFunction, "AnalyticalSolnTgtExact",
 moab::Remapper::TargetMesh, testFunction,
 &tgtProjectedFunction, "ProjectedSolnTgt" );MB_CHK_ERR( rval );
 // rval = mbCore->write_file ( "tgtWithSolnTag.h5m", NULL, writeOptions,
 // &runCtx->meshsets[1], 1 ); MB_CHK_ERR ( rval );
 runCtx->timer_pop();
 
 runCtx->timer_push( "compute solution projection on target grid" );
 rval = weightMap->ApplyWeights( srcAnalyticalFunction, tgtProjectedFunction );MB_CHK_ERR( rval );
 runCtx->timer_pop();
 rval = mbCore->write_file( "tgtWithSolnTag2.h5m", NULL, writeOptions, &runCtx->meshsets[1], 1 );MB_CHK_ERR( rval );
 
 if( nprocs == 1 && runCtx->baselineFile.size() )
 {
 // save the field from tgtWithSolnTag2 in a text file, and global ids for cells
 moab::Range tgtCells;
 rval = mbCore->get_entities_by_dimension( runCtx->meshsets[1], 2, tgtCells );MB_CHK_ERR( rval );
 std::vector< int > globIds;
 globIds.resize( tgtCells.size() );
 std::vector< double > vals;
 vals.resize( tgtCells.size() );
 moab::Tag projTag;
 rval = mbCore->tag_get_handle( "ProjectedSolnTgt", projTag );MB_CHK_ERR( rval );
 moab::Tag gid = mbCore->globalId_tag();
 rval = mbCore->tag_get_data( gid, tgtCells, &globIds[0] );MB_CHK_ERR( rval );
 rval = mbCore->tag_get_data( projTag, tgtCells, &vals[0] );MB_CHK_ERR( rval );
 std::fstream fs;
 fs.open( runCtx->baselineFile.c_str(), std::fstream::out );
 fs << std::setprecision( 15 ); // maximum precision for doubles
 for( size_t i = 0; i < tgtCells.size(); i++ )
 fs << globIds[i] << " " << vals[i] << "\n";
 fs.close();
 // for good measure, save the source file too, with the tag AnalyticalSolnSrcExact
 // it will be used later to test, along with a target file
 rval = mbCore->write_file( "srcWithSolnTag.h5m", NULL, writeOptions, &runCtx->meshsets[0], 1 );MB_CHK_ERR( rval );
 }
 runCtx->timer_push( "compute error metrics against analytical solution on target grid" );
 std::map< std::string, double > errMetrics;
 rval = weightMap->ComputeMetrics( moab::Remapper::TargetMesh, tgtAnalyticalFunction,
 tgtProjectedFunction, errMetrics, true );MB_CHK_ERR( rval );
 runCtx->timer_pop();
 }
 
 delete weightMap;
 }
 }
 
 // Clean up
 remapper.clear();
 delete runCtx;
 delete mbCore;
 
#ifdef MOAB_HAVE_MPI
 MPI_Finalize();
#endif
 exit( 0 );
}
 
static moab::ErrorCode CreateTempestMesh( ToolContext& ctx, moab::TempestRemapper& remapper, Mesh* tempest_mesh )
{
 moab::ErrorCode rval = moab::MB_SUCCESS;
 int err;
 moab::DebugOutput& outputFormatter = ctx.outputFormatter;
 
 if( !ctx.proc_id )
 {
 outputFormatter.printf( 0, "Creating TempestRemap Mesh object ...\n" );
 }
 
 if( ctx.meshType == moab::TempestRemapper::OVERLAP_FILES )
 {
 // For the overlap method, choose between: "fuzzy", "exact" or "mixed"
 err = GenerateOverlapMesh( ctx.inFilenames[0], ctx.inFilenames[1], *tempest_mesh, ctx.outFilename, "NetCDF4",
 "exact", true );
 
 if( err )
 {
 rval = moab::MB_FAILURE;
 }
 else
 {
 ctx.meshes.push_back( tempest_mesh );
 }
 }
 else if( ctx.meshType == moab::TempestRemapper::OVERLAP_MEMORY )
 {
 // Load the meshes and validate
 ctx.meshsets.resize( 3 );
 ctx.meshes.resize( 3 );
 ctx.meshsets[0] = remapper.GetMeshSet( moab::Remapper::SourceMesh );
 ctx.meshsets[1] = remapper.GetMeshSet( moab::Remapper::TargetMesh );
 ctx.meshsets[2] = remapper.GetMeshSet( moab::Remapper::OverlapMesh );
 
 // First the source
 rval = remapper.LoadMesh( moab::Remapper::SourceMesh, ctx.inFilenames[0], moab::TempestRemapper::DEFAULT );MB_CHK_ERR( rval );
 ctx.meshes[0] = remapper.GetMesh( moab::Remapper::SourceMesh );
 
 // Next the target
 rval = remapper.LoadMesh( moab::Remapper::TargetMesh, ctx.inFilenames[1], moab::TempestRemapper::DEFAULT );MB_CHK_ERR( rval );
 ctx.meshes[1] = remapper.GetMesh( moab::Remapper::TargetMesh );
 
 // Now let us construct the overlap mesh, by calling TempestRemap interface directly
 // For the overlap method, choose between: "fuzzy", "exact" or "mixed"
 err = GenerateOverlapWithMeshes( *ctx.meshes[0], *ctx.meshes[1], *tempest_mesh, "" /*ctx.outFilename*/,
 "NetCDF4", "exact", false );
 
 if( err )
 {
 rval = moab::MB_FAILURE;
 }
 else
 {
 remapper.SetMesh( moab::Remapper::OverlapMesh, tempest_mesh );
 ctx.meshes[2] = remapper.GetMesh( moab::Remapper::OverlapMesh );
 // ctx.meshes.push_back(*tempest_mesh);
 }
 }
 else if( ctx.meshType == moab::TempestRemapper::OVERLAP_MOAB )
 {
 ctx.meshsets.resize( 3 );
 ctx.meshes.resize( 3 );
 ctx.meshsets[0] = remapper.GetMeshSet( moab::Remapper::SourceMesh );
 ctx.meshsets[1] = remapper.GetMeshSet( moab::Remapper::TargetMesh );
 ctx.meshsets[2] = remapper.GetMeshSet( moab::Remapper::OverlapMesh );
 
 const double radius_src = 1.0 /*2.0*acos(-1.0)*/;
 const double radius_dest = 1.0 /*2.0*acos(-1.0)*/;
 
 //const char* additional_read_opts = ( ctx.n_procs > 1 ? "NO_SET_CONTAINING_PARENTS;" : "" );
 std::string additional_read_opts_src = get_file_read_options( ctx, ctx.inFilenames[0] );
 std::vector< int > smetadata, tmetadata;
 // Load the source mesh and validate
 rval =
 remapper.LoadNativeMesh( ctx.inFilenames[0], ctx.meshsets[0], smetadata, additional_read_opts_src.c_str() );MB_CHK_ERR( rval );
 if( smetadata.size() )
 {
 remapper.SetMeshType( moab::Remapper::SourceMesh, smetadata );
 }
 // Rescale the radius of both to compute the intersection
 rval = moab::IntxUtils::ScaleToRadius( ctx.mbcore, ctx.meshsets[0], radius_src );MB_CHK_ERR( rval );
 rval = remapper.ConvertMeshToTempest( moab::Remapper::SourceMesh );MB_CHK_ERR( rval );
 ctx.meshes[0] = remapper.GetMesh( moab::Remapper::SourceMesh );
 
 // Load the target mesh and validate
 std::string addititional_read_opts_tgt = get_file_read_options( ctx, ctx.inFilenames[1] );
 rval = remapper.LoadNativeMesh( ctx.inFilenames[1], ctx.meshsets[1], tmetadata,
 addititional_read_opts_tgt.c_str() );MB_CHK_ERR( rval );
 if( tmetadata.size() )
 {
 remapper.SetMeshType( moab::Remapper::TargetMesh, tmetadata );
 }
 rval = moab::IntxUtils::ScaleToRadius( ctx.mbcore, ctx.meshsets[1], radius_dest );MB_CHK_ERR( rval );
 rval = remapper.ConvertMeshToTempest( moab::Remapper::TargetMesh );MB_CHK_ERR( rval );
 ctx.meshes[1] = remapper.GetMesh( moab::Remapper::TargetMesh );
 
 // ctx.meshes[0]->Write( "SourceMeshMBTR.g" );
 // ctx.meshes[1]->Write( "TargetMeshMBTR.g" );
 }
 else if( ctx.meshType == moab::TempestRemapper::ICO )
 {
 err = GenerateICOMesh( *tempest_mesh, ctx.blockSize, ctx.computeDual, ctx.outFilename, "NetCDF4" );
 
 if( err )
 {
 rval = moab::MB_FAILURE;
 }
 else
 {
 ctx.meshes.push_back( tempest_mesh );
 }
 }
 else if( ctx.meshType == moab::TempestRemapper::RLL )
 {
 err = GenerateRLLMesh( *tempest_mesh, // Mesh& meshOut,
 ctx.blockSize * 2, ctx.blockSize, // int nLongitudes, int nLatitudes,
 0.0, 360.0, // double dLonBegin, double dLonEnd,
 -90.0, 90.0, // double dLatBegin, double dLatEnd,
 false, false, false, // bool fGlobalCap, bool fFlipLatLon, bool fForceGlobal,
 "" /*ctx.inFilename*/, "",
 "", // std::string strInputFile, std::string strInputFileLonName, std::string
 // strInputFileLatName,
 ctx.outFilename, "NetCDF4", // std::string strOutputFile, std::string strOutputFormat
 true // bool fVerbose
 );
 
 if( err )
 {
 rval = moab::MB_FAILURE;
 }
 else
 {
 ctx.meshes.push_back( tempest_mesh );
 }
 }
 else // default
 {
 err = GenerateCSMesh( *tempest_mesh, ctx.blockSize, ctx.outFilename, "NetCDF4" );
 
 if( err )
 {
 rval = moab::MB_FAILURE;
 }
 else
 {
 ctx.meshes.push_back( tempest_mesh );
 }
 }
 
 if( ctx.meshType != moab::TempestRemapper::OVERLAP_MOAB && !tempest_mesh )
 {
 std::cout << "Tempest Mesh is not a complete object; Quitting...";
 exit( -1 );
 }
 
 return rval;
}
 
///////////////////////////////////////////////
// Test functions
 
double sample_slow_harmonic( double dLon, double dLat )
{
 return ( 2.0 + cos( dLat ) * cos( dLat ) * cos( 2.0 * dLon ) );
}
 
double sample_fast_harmonic( double dLon, double dLat )
{
 return ( 2.0 + pow( sin( 2.0 * dLat ), 16.0 ) * cos( 16.0 * dLon ) );
 // return (2.0 + pow(cos(2.0 * dLat), 16.0) * cos(16.0 * dLon));
}
 
double sample_constant( double /*dLon*/, double /*dLat*/ )
{
 return 1.0;
}
 
double sample_stationary_vortex( double dLon, double dLat )
{
 const double dLon0 = 0.0;
 const double dLat0 = 0.6;
 const double dR0 = 3.0;
 const double dD = 5.0;
 const double dT = 6.0;
 
 /// Find the rotated longitude and latitude of a point on a sphere
 /// with pole at (dLonC, dLatC).
 {
 double dSinC = sin( dLat0 );
 double dCosC = cos( dLat0 );
 double dCosT = cos( dLat );
 double dSinT = sin( dLat );
 
 double dTrm = dCosT * cos( dLon - dLon0 );
 double dX = dSinC * dTrm - dCosC * dSinT;
 double dY = dCosT * sin( dLon - dLon0 );
 double dZ = dSinC * dSinT + dCosC * dTrm;
 
 dLon = atan2( dY, dX );
 if( dLon < 0.0 )
 {
 dLon += 2.0 * M_PI;
 }
 dLat = asin( dZ );
 }
 
 double dRho = dR0 * cos( dLat );
 double dVt = 3.0 * sqrt( 3.0 ) / 2.0 / cosh( dRho ) / cosh( dRho ) * tanh( dRho );
 
 double dOmega;
 if( dRho == 0.0 )
 {
 dOmega = 0.0;
 }
 else
 {
 dOmega = dVt / dRho;
 }
 
 return ( 1.0 - tanh( dRho / dD * sin( dLon - dOmega * dT ) ) );
}
 
///////////////////////////////////////////////