mbtempest.cpp File Reference

MOAB-Tempest: A powerful mesh generation and remapping tool for climate and weather applications. More...

#include <iostream>
#include <iomanip>
#include <cstdlib>
#include <vector>
#include <string>
#include <memory>
#include <sstream>
#include <cassert>
#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"

Include dependency graph for mbtempest.cpp:

Go to the source code of this file.

Classes
class	ToolContext
	Context class for MOAB-TempestRemap tool configuration and state management. More...

Namespaces
	anonymous_namespace{mbtempest.cpp}
	Sample functions for testing remapping operations.

Macros
#define	TR_CHK_SET_ERR(err, msg)

Functions
static moab::ErrorCode	CreateTempestMesh (ToolContext &ctx, moab::TempestRemapper &remapper, Mesh *tempest_mesh)
	Creates a TempestRemap mesh based on the provided context and mesh type. More...
static constexpr double	sample_constant (double, double) noexcept
	Constant sample function. More...
static double	sample_slow_harmonic (double dLon, double dLat) noexcept
	Sample function with slow harmonic variation. More...
static double	sample_fast_harmonic (double dLon, double dLat) noexcept
	Sample function with fast harmonic variation. More...
static double	sample_stationary_vortex (double dLon, double dLat) noexcept
	Sample function representing a stationary vortex. More...
int	main (int argc, char *argv[])
moab::ErrorCode	anonymous_namespace{mbtempest.cpp}::handleOverlapMOAB (ToolContext &ctx, moab::TempestRemapper &remapper)
moab::ErrorCode	anonymous_namespace{mbtempest.cpp}::convertAndWriteMOABMesh (ToolContext &ctx, moab::TempestRemapper &remapper, Mesh *tempest_mesh)
	Convert a generated TempestRemap mesh to MOAB format and write as h5m file. More...
moab::ErrorCode	anonymous_namespace{mbtempest.cpp}::handleICOMesh (ToolContext &ctx, moab::TempestRemapper &remapper, Mesh *tempest_mesh)
moab::ErrorCode	anonymous_namespace{mbtempest.cpp}::handleRLLMesh (ToolContext &ctx, moab::TempestRemapper &remapper, Mesh *tempest_mesh)
moab::ErrorCode	anonymous_namespace{mbtempest.cpp}::handleCSMesh (ToolContext &ctx, moab::TempestRemapper &remapper, Mesh *tempest_mesh)

Variables
constexpr double	anonymous_namespace{mbtempest.cpp}::VORTEX_LON0 = 0.0
constexpr double	anonymous_namespace{mbtempest.cpp}::VORTEX_LAT0 = 0.6
constexpr double	anonymous_namespace{mbtempest.cpp}::VORTEX_R0 = 3.0
constexpr double	anonymous_namespace{mbtempest.cpp}::VORTEX_D = 5.0
constexpr double	anonymous_namespace{mbtempest.cpp}::VORTEX_T = 6.0

Detailed Description

MOAB-Tempest: A powerful mesh generation and remapping tool for climate and weather applications.

Overview

MOAB-Tempest is a command-line tool that provides mesh generation and conservative remapping capabilities for climate and weather modeling. It combines the power of MOAB (Mesh-Oriented datABase) with the TempestRemap library to enable high-performance, parallel mesh generation and remapping operations.

Key Features

• Generation of various spherical mesh types (Cubed-Sphere, RLL, Icosahedral)
• Support for high-order discretization methods (FV, CGLL, DGLL)
• Conservative remapping between different mesh types
• Parallel processing support via MPI
• Flexible I/O with support for multiple file formats
• Built-in analytical functions for testing and validation

Supported Algorithms

• Mesh Generation:
  • Cubed-Sphere (CS) meshes
  • Regular Latitude-Longitude (RLL) meshes
  • Icosahedral (ICO) meshes
  • Overlap meshes for remapping
• Remapping Methods:
  • Finite Volume (FV)
  • Continuous Galerkin (CGLL)
  • Discontinuous Galerkin (DGLL)
  • Monotonic and high-order variants

Basic Usage Examples

# Generate a Cubed-Sphere mesh with resolution 25
./mbtempest --type 0 --res 25 --file cubed_sphere_mesh.h5m
 
# Generate a RLL mesh with resolution 90x180 (lon x lat)
./mbtempest --type 1 --res 90 --file rll_mesh.h5m
 
# Generate an Icosahedral mesh with resolution 25 (dual mesh)
./mbtempest --type 2 --res 25 --dual --file icosahedral_dual_mesh.h5m
 
# Compute overlap between two meshes
./mbtempest --type 5 --load mesh1.h5m --load mesh2.h5m intx intersection_mesh.h5m
 
# Generate a remapping weights file between two meshes: FV to FV (default)
./mbtempest --type 5 --load source_mesh.h5m --load target_mesh.h5m --file weights.nc
 
# Generate remapping weights file between two meshes: making it explicit (SE to FV)
./mbtempest --type 5 --load source_mesh.h5m --load target_mesh.h5m \
            --order 4 --method cgll --global_id GLOBAL_DOFS \
            --order 1 --method fv --global_id GLOBAL_ID \
            --file weights_se_to_fv.nc

Command Line Options

Run './mbtempest –help' for a complete list of available options.

Notes

• For parallel execution, use MPI launcher (e.g., mpirun, mpiexec)
• Output formats: .h5m (MOAB), .nc (NetCDF), .exo (ExodusII)
• Requires MOAB and TempestRemap libraries

Author

MOAB Development Team

Date

Created: 2023

Definition in file mbtempest.cpp.

Macro Definition Documentation

TR_CHK_SET_ERR

#define TR_CHK_SET_ERR	(		err,
			msg
	)

Value:

    if( err )                                                                     \
    {                                                                             \
        std::cout << "MOAB-TempestRemap Failure. ErrorCode (" << ( err ) << ") "; \
        MB_CHK_SET_ERR( moab::MB_FAILURE, msg );                                  \
    }

Definition at line 1675 of file mbtempest.cpp.

Function Documentation

CreateTempestMesh()

static moab::ErrorCode CreateTempestMesh ( ToolContext &  ctx, moab::TempestRemapper &  remapper, Mesh *  tempest_mesh  )

static

Creates a TempestRemap mesh based on the provided context and mesh type.

Parameters

ctx	Tool context containing configuration and state
remapper	TempestRemap instance for mesh operations
tempest_mesh	Output parameter for the created mesh

Returns

moab::ErrorCode Status of the operation

Definition at line 1989 of file mbtempest.cpp.

{
    using namespace moab;
    using RemapperType = moab::TempestRemapper;
 
    auto& outputFormatter = ctx.outputFormatter;
 
    try
    {
        switch( ctx.meshType )
        {
            case RemapperType::OVERLAP_FILES:
#ifdef MOAB_HAVE_NETCDF
                if( !ctx.proc_id ) outputFormatter.printf( 0, "Creating TempestRemap overlap mesh ...\n" );
                return handleOverlapFiles( ctx, tempest_mesh );
#else
                MB_CHK_SET_ERR( moab::MB_FAILURE,
                                "OVERLAP_FILES mode requires NetCDF (TempestRemap file-based overlap "
                                "generation); build with NetCDF or use OVERLAP_MOAB mode instead" );
#endif
 
            case RemapperType::OVERLAP_MEMORY:
#ifdef MOAB_HAVE_NETCDF
                if( !ctx.proc_id )
                    outputFormatter.printf( 0, "Convert MOAB overlap files to TempestRemap format in-memory ...\n" );
                return handleOverlapMemory( ctx, remapper, tempest_mesh );
#else
                MB_CHK_SET_ERR( moab::MB_FAILURE,
                                "OVERLAP_MEMORY mode requires NetCDF (TempestRemap-native mesh loading); "
                                "build with NetCDF or use OVERLAP_MOAB mode instead" );
#endif
 
            case RemapperType::OVERLAP_MOAB:
                if( !ctx.proc_id )
                    outputFormatter.printf( 0, "Convert MOAB meshes to TempestRemap format in-memory ...\n" );
                return handleOverlapMOAB( ctx, remapper );
 
            case RemapperType::ICO:
                if( !ctx.proc_id ) outputFormatter.printf( 0, "Creating TempestRemap ICO mesh ...\n" );
                return handleICOMesh( ctx, remapper, tempest_mesh );
 
            case RemapperType::RLL:
                if( !ctx.proc_id ) outputFormatter.printf( 0, "Creating TempestRemap RLL mesh ...\n" );
                return handleRLLMesh( ctx, remapper, tempest_mesh );
 
            default:  // Default to CS mesh
                if( !ctx.proc_id ) outputFormatter.printf( 0, "Creating TempestRemap CS mesh ...\n" );
                return handleCSMesh( ctx, remapper, tempest_mesh );
        }
    }
    catch( const std::exception& e )
    {
        std::cerr << "Error in CreateTempestMesh: " << e.what() << "\n";
        return MB_FAILURE;
    }
}

References anonymous_namespace{mbtempest.cpp}::handleCSMesh(), anonymous_namespace{mbtempest.cpp}::handleICOMesh(), anonymous_namespace{mbtempest.cpp}::handleOverlapMOAB(), anonymous_namespace{mbtempest.cpp}::handleRLLMesh(), MB_CHK_SET_ERR, ToolContext::meshType, ToolContext::outputFormatter, moab::DebugOutput::printf(), and ToolContext::proc_id.

Referenced by main().

main()

int main	(	int	argc,
		char *	argv[]
	)

Definition at line 970 of file mbtempest.cpp.

{
    try
    {
#ifdef MOAB_HAVE_NETCDF
    NcError error( NcError::verbose_nonfatal );
#endif
    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( nullptr == 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();
    // The command-line choice supersedes the remapper default, so that intersection,
    // coverage construction and orientation fixups all use one formula.
    remapper.SetAreaMethod( runCtx->areaMethod );
    // --rrmgrids: the two domains need not coincide, so cells may be partially covered.
    remapper.SetRegionalMesh( runCtx->rrmGrids );
 
    // Area formula selected by --area_method (default: Van Oosterom-Strackee)
    moab::IntxAreaUtils areaAdaptor( runCtx->areaMethod );
 
    Mesh* tempest_mesh = new Mesh();
    MB_CHK_SET_ERR( CreateTempestMesh( *runCtx, remapper, tempest_mesh ), "Failed to create tempest mesh" );
 
    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
        MB_CHK_SET_ERR( pcomm->check_all_shared_handles(), "Failed to check all shared handles" );
#endif
 
        // Load the meshes and validate
        MB_CHK_SET_ERR( remapper.ConvertTempestMesh( moab::Remapper::SourceMesh ), "Failed to convert source mesh" );
        MB_CHK_SET_ERR( remapper.ConvertTempestMesh( moab::Remapper::TargetMesh ), "Failed to convert target mesh" );
        MB_CHK_SET_ERR( remapper.ConvertTempestMesh( moab::Remapper::OverlapMesh ), "Failed to convert overlap mesh" );
        if( !runCtx->skip_io )
        {
            MB_CHK_SET_ERR( mbCore->write_mesh( "tempest_intersection.h5m", &runCtx->meshsets[2], 1 ),
                            "Failed to write TempestRemap intersection mesh in MOAB format" );
        }
 
        // print verbosely about the problem setting
        size_t velist[6], gvelist[6];
        {
            moab::Range rintxverts, rintxelems;
            MB_CHK_SET_ERR( mbCore->get_entities_by_dimension( runCtx->meshsets[0], 0, rintxverts ),
                            "Failed to get vertices" );
            MB_CHK_SET_ERR( mbCore->get_entities_by_dimension( runCtx->meshsets[0], 2, rintxelems ),
                            "Failed to get elements" );
            velist[0] = rintxverts.size();
            velist[1] = rintxelems.size();
 
            moab::Range bintxverts, bintxelems;
            MB_CHK_SET_ERR( mbCore->get_entities_by_dimension( runCtx->meshsets[1], 0, bintxverts ),
                            "Failed to get vertices" );
            MB_CHK_SET_ERR( mbCore->get_entities_by_dimension( runCtx->meshsets[1], 2, bintxelems ),
                            "Failed to get elements" );
            velist[2] = bintxverts.size();
            velist[3] = 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, runCtx->areaMethod );
            mbintx->set_error_tolerance( runCtx->epsrel );
            mbintx->set_box_error( runCtx->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
            MB_CHK_SET_ERR( mbintx->FindMaxEdges( runCtx->meshsets[0], runCtx->meshsets[1] ),
                            "Failed to find max edges" );
 
#ifdef MOAB_HAVE_MPI
            moab::Range local_verts;
            MB_CHK_SET_ERR( mbintx->build_processor_euler_boxes( runCtx->meshsets[1], local_verts ),
                            "Failed to build processor euler boxes" );
 
            runCtx->timer_pop();
 
            moab::EntityHandle covering_set;
            runCtx->timer_push( "communicate the mesh" );
            // we compute just intersection here, no need for extra ghost layers anyway
            // ghost layers are needed in coverage for bilinear map, which does not actually need intersection
            // this will be fixed in the future, bilinear map needs just coverage, not intersection
            // so I am not passing the ghost layer here, even though there is an option in runCtx for a ghost layer
            // NOTE: This is a communication-heavy kernel if mesh is distributed very differently
            MB_CHK_SET_ERR( mbintx->construct_covering_set( runCtx->meshsets[0], covering_set ),
                            "Failed to construct covering set" );
            runCtx->timer_pop();
 
            // print verbosely about the problem setting
            {
                moab::Range cintxverts, cintxelems;
                MB_CHK_SET_ERR( mbCore->get_entities_by_dimension( covering_set, 0, cintxverts ),
                                "Failed to get vertices" );
                MB_CHK_SET_ERR( mbCore->get_entities_by_dimension( covering_set, 2, cintxelems ),
                                "Failed to get elements" );
                velist[4] = cintxverts.size();
                velist[5] = cintxelems.size();
            }
 
            MPI_Reduce( velist, gvelist, 6, MPI_UINT64_T, MPI_SUM, 0, MPI_COMM_WORLD );
 
#else
            moab::EntityHandle covering_set = runCtx->meshsets[0];
            for( int i = 0; i < 6; i++ )
                gvelist[i] = velist[i];
#endif
 
            if( !proc_id )
            {
                outputFormatter.printf( 0, "The source set contains %lu vertices and %lu elements \n", gvelist[0],
                                        gvelist[0] );
                outputFormatter.printf( 0, "The covering set contains %lu vertices and %lu elements \n", gvelist[2],
                                        gvelist[2] );
                outputFormatter.printf( 0, "The target set contains %lu vertices and %lu elements \n", gvelist[1],
                                        gvelist[1] );
            }
 
            // 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" );
            MB_CHK_SET_ERR( mbCore->create_meshset( moab::MESHSET_SET, intxset ), "Can't create new set" );
            MB_CHK_SET_ERR( mbintx->intersect_meshes( covering_set, runCtx->meshsets[1], intxset ),
                            "Can't compute the intersection of meshes on the sphere" );
            runCtx->timer_pop();
 
            // free the memory
            delete mbintx;
        }
 
        {
            moab::Range intxelems, intxverts;
            MB_CHK_SET_ERR( mbCore->get_entities_by_dimension( intxset, 2, intxelems ), "Failed to get elements" );
            MB_CHK_SET_ERR( mbCore->get_entities_by_dimension( intxset, 0, intxverts, true ),
                            "Failed to get vertices" );
            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
        if( !runCtx->skip_io )
        {
            MB_CHK_SET_ERR( mbCore->write_mesh( "moab_intersection.h5m", &intxset, 1 ),
                            "Failed to write the intersection" );
        }
 
        if( runCtx->computeWeights )
        {
            runCtx->timer_push( "compute weights with the Tempest meshes" );
            // Call to generate an offline map with the tempest meshes
            OfflineMap weightMap;
            if( 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 ) != 0 )
                throw std::runtime_error( "Could not generate offline map with TempestRemap" );
            runCtx->timer_pop();
 
            std::map< std::string, std::string > mapAttributes;
#ifdef MOAB_HAVE_NETCDF
            if( !runCtx->skip_io ) weightMap.Write( "outWeights.nc", mapAttributes );
#else
            (void)mapAttributes;
            MB_CHK_SET_ERR( moab::MB_FAILURE,
                            "Writing a TempestRemap OfflineMap (OVERLAP_MEMORY path) requires NetCDF; use the "
                            "MOAB (OVERLAP_MOAB) workflow, which writes weights via PnetCDF/HDF5 instead" );
#endif
        }
    }
    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
        MB_CHK_SET_ERR( pcomm->check_all_shared_handles(), "Checking shared handles failed." );
#endif
 
        // print verbosely about the problem setting
        size_t velist[4] = { 0, 0, 0, 0 }, gvelist[4] = { 0, 0, 0, 0 };
        {
            moab::Range srcverts, srcelems;
            MB_CHK_SET_ERR( mbCore->get_entities_by_dimension( runCtx->meshsets[0], 0, srcverts ),
                            "Failed to get vertices" );
            MB_CHK_SET_ERR( mbCore->get_entities_by_dimension( runCtx->meshsets[0], 2, srcelems ),
                            "Failed to get elements" );
            MB_CHK_SET_ERR( moab::IntxUtils::fix_degenerate_quads( mbCore, runCtx->meshsets[0] ),
                            "Failed to fix degenerate quads" );
            if( runCtx->enforceConvexity )
            {
                MB_CHK_SET_ERR( moab::IntxUtils::enforce_convexity( mbCore, runCtx->meshsets[0], proc_id ),
                                "Failed to enforce convexity" );
            }
            MB_CHK_SET_ERR( areaAdaptor.positive_orientation( mbCore, runCtx->meshsets[0], radius_src ),
                            "Failed to enforce positive orientation" );
            velist[0] = srcverts.size();
            velist[1] = srcelems.size();
 
            moab::Range tgtverts, tgtelems;
            MB_CHK_SET_ERR( mbCore->get_entities_by_dimension( runCtx->meshsets[1], 0, tgtverts ),
                            "Failed to get vertices" );
            MB_CHK_SET_ERR( mbCore->get_entities_by_dimension( runCtx->meshsets[1], 2, tgtelems ),
                            "Failed to get elements" );
            MB_CHK_SET_ERR( moab::IntxUtils::fix_degenerate_quads( mbCore, runCtx->meshsets[1] ),
                            "Failed to fix degenerate quads" );
            if( runCtx->enforceConvexity )
            {
                MB_CHK_SET_ERR( moab::IntxUtils::enforce_convexity( mbCore, runCtx->meshsets[1], proc_id ),
                                "Failed to enforce convexity" );
            }
            MB_CHK_SET_ERR( areaAdaptor.positive_orientation( mbCore, runCtx->meshsets[1], radius_dest ),
                            "Failed to enforce positive orientation" );
            velist[2] = tgtverts.size();
            velist[3] = tgtelems.size();
        }
        //MB_CHK_SET_ERR( mbCore->write_file( "source_mesh.h5m", nullptr, writeOptions, &runCtx->meshsets[0], 1 ), "Could not write source mesh" );
        //MB_CHK_SET_ERR( mbCore->write_file( "target_mesh.h5m", nullptr, writeOptions, &runCtx->meshsets[1], 1 ), "Could not write target mesh" );
 
        // if( runCtx->nlayers && nprocs > 1 )
        // {
        //     remapper.ResetMeshSet( moab::Remapper::SourceMesh, runCtx->meshsets[3] );
        //     runCtx->meshes[0] = remapper.GetMesh( moab::Remapper::SourceMesh );  //  ?
        // }
 
        // First compute the covering set such that the target elements are fully covered by the
        // local source grid
        runCtx->timer_push( "construct covering set for intersection" );
        // Ghosting and the gnomonic-projection coverage path are mutually
        // exclusive.  Ghost layers are now enabled by default for the
        // bilinear-family methods, so warn rather than silently dropping an
        // explicit --gnomonic request.
        if( runCtx->nlayers > 0 && runCtx->useGnomonicProjection )
        {
            if( !proc_id )
                std::cout << "  [WARNING] --gnomonic is ignored when ghost layers are used (nlayers = "
                          << runCtx->nlayers << "); pass '--ghost 0' to force the gnomonic coverage path.\n";
            runCtx->useGnomonicProjection = false;
        }
        MB_CHK_SET_ERR( remapper.ConstructCoveringSet( runCtx->epsrel, 1.0, 1.0, runCtx->boxeps, runCtx->rrmGrids,
                                                       runCtx->useGnomonicProjection, runCtx->nlayers ),
                        "Failed to construct covering set" );
        runCtx->timer_pop();
 
#ifdef MOAB_HAVE_MPI
        MPI_Reduce( velist, gvelist, 4, MPI_UINT64_T, MPI_SUM, 0, MPI_COMM_WORLD );
#else
        for( int i = 0; i < 4; i++ )
            gvelist[i] = velist[i];
#endif
        if( !proc_id && runCtx->print_diagnostics )
        {
            outputFormatter.printf( 0, "The source set contains %lu vertices and %lu elements \n", gvelist[0],
                                    gvelist[1] );
            outputFormatter.printf( 0, "The target set contains %lu vertices and %lu elements \n", gvelist[2],
                                    gvelist[3] );
        }
 
        if( runCtx->skip_intersection )
        {
            if( !proc_id ) outputFormatter.printf( 0, "Skipping mesh intersection computation.\n" );
        }
        else
        {
            // Compute intersections with MOAB with either the Kd-tree or the advancing front algorithm
            runCtx->timer_push( "setup and compute mesh intersections" );
            MB_CHK_SET_ERR( remapper.ComputeOverlapMesh( runCtx->kdtreeSearch, false ),
                            "Failed to compute mesh intersections" );
            runCtx->timer_pop();
        }
 
        // print some diagnostic checks to see if the overlap grid resolved the input meshes
        // correctly
        // Compute ghost overlap elements once; reused for both area diagnostics and intx file write
        moab::Range ghostOverlapElems;
#ifdef MOAB_HAVE_MPI
        if( nprocs > 1 && !runCtx->skip_intersection )
            MB_CHK_SET_ERR( remapper.GetOverlapAugmentedEntities( ghostOverlapElems ),
                            "Failed to get ghost overlap entities" );
#endif
 
        double dTotalOverlapArea = 0.0;
        if( runCtx->print_diagnostics && !runCtx->skip_intersection )
        {
             // Areas for source, target, overlap meshes
            double local_areas[3]  = { 0, 0, 0 },
                   global_areas[3] = { 0, 0, 0 };
 
            // Helper: compute area of a meshset excluding cells with GRID_IMASK==0.
            // Both source and target SCRIP grids may have a land/sea mask; the intersection
            // only covers unmasked cells, so comparing full-mesh areas gives a misleading error.
            auto area_unmasked = [&]( moab::EntityHandle meshset, double radius ) -> double {
                moab::Tag imaskTag = 0;
                mbCore->tag_get_handle( "GRID_IMASK", imaskTag );
                if( !imaskTag ) return areaAdaptor.area_on_sphere( mbCore, meshset, radius );
                moab::Range cells;
                mbCore->get_entities_by_dimension( meshset, 2, cells );
                std::vector< int > masks( cells.size(), 1 );
                mbCore->tag_get_data( imaskTag, cells, masks.data() );
                moab::Range maskedCells;
                size_t idx = 0;
                for( auto it = cells.begin(); it != cells.end(); ++it, ++idx )
                    if( !masks[idx] ) maskedCells.insert( *it );
                moab::Range unmasked = moab::subtract( cells, maskedCells );
                moab::EntityHandle tmpSet;
                mbCore->create_meshset( moab::MESHSET_SET, tmpSet );
                mbCore->add_entities( tmpSet, unmasked );
                double area = areaAdaptor.area_on_sphere( mbCore, tmpSet, radius );
                mbCore->delete_entities( &tmpSet, 1 );
                return area;
            };
 
            local_areas[0] = area_unmasked( runCtx->meshsets[0], radius_src );
            local_areas[1] = area_unmasked( runCtx->meshsets[1], radius_dest );
            // Exclude ghost overlap elements from area sum to avoid double-counting after MPI_Allreduce
            {
                moab::Range ownedOverlapElems;
                MB_CHK_SET_ERR( mbCore->get_entities_by_dimension( runCtx->meshsets[2], 2, ownedOverlapElems ),
                                "Failed to get overlap elements" );
                ownedOverlapElems = moab::subtract( ownedOverlapElems, ghostOverlapElems );
                moab::EntityHandle ownedOverlapSet;
                MB_CHK_SET_ERR( mbCore->create_meshset( moab::MESHSET_SET, ownedOverlapSet ),
                                "Can't create owned overlap meshset" );
                MB_CHK_SET_ERR( mbCore->add_entities( ownedOverlapSet, ownedOverlapElems ),
                                "Can't add owned overlap elements" );
                local_areas[2] = areaAdaptor.area_on_sphere( mbCore, ownedOverlapSet, radius_src );
                MB_CHK_SET_ERR( mbCore->delete_entities( &ownedOverlapSet, 1 ), "Can't delete temp meshset" );
            }
 
#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, "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->rrmGrids )
                {
                    // For coincident domains the shortfall below is roundoff and the
                    // relative errors above already say so.  For a regional mesh it is
                    // the physically meaningful quantity: how much of each mesh the
                    // other one actually covers.
                    outputFormatter.printf( 0,
                                            "regional mesh: overlap covers %6.2f%% of source and %6.2f%% of "
                                            "target area\n",
                                            100.0 * global_areas[2] / global_areas[0],
                                            100.0 * global_areas[2] / global_areas[1] );
                }
            }
            dTotalOverlapArea = global_areas[2];
        }
 
        if( runCtx->intxFilename.size() && !runCtx->skip_intersection )
        {
            moab::EntityHandle writableOverlapSet;
            MB_CHK_SET_ERR( mbCore->create_meshset( moab::MESHSET_SET, writableOverlapSet ), "Can't create new set" );
            moab::EntityHandle meshOverlapSet = remapper.GetMeshSet( moab::Remapper::OverlapMesh );
            moab::Range ovEnts;
            MB_CHK_SET_ERR( mbCore->get_entities_by_dimension( meshOverlapSet, 2, ovEnts ), "Can't create new set" );
            MB_CHK_SET_ERR( mbCore->get_entities_by_dimension( meshOverlapSet, 0, ovEnts ), "Can't create new set" );
 
#ifdef MOAB_HAVE_MPI
            // Exclude ghost overlap elements from the write: each ghost element is owned by another
            // rank and will be written from there. Including ghosts here causes duplicate entity
            // handles in the parallel HDF5 output and deadlocks the collective write.
            if( nprocs > 1 )
            {
                ovEnts = moab::subtract( ovEnts, ghostOverlapElems );
#ifdef MOAB_DBG
                if( !runCtx->skip_io )
                {
                    std::stringstream filename;
                    filename << "aug_overlap" << runCtx->pcomm->rank() << ".h5m";
                    MB_CHK_SET_ERR( mbCore->write_file( filename.str().c_str(), 0, 0, &meshOverlapSet, 1 ),
                                    "Failed to write the overlap set" );
                }
#endif
            }
#endif
            MB_CHK_SET_ERR( mbCore->add_entities( writableOverlapSet, ovEnts ), "adding local intx cells failed" );
 
#ifdef MOAB_HAVE_MPI
#ifdef MOAB_DBG
            if( nprocs > 1 && !runCtx->skip_io )
            {
                std::stringstream filename;
                filename << "writable_intx_" << runCtx->pcomm->rank() << ".h5m";
                MB_CHK_SET_ERR( mbCore->write_file( filename.str().c_str(), 0, 0, &writableOverlapSet, 1 ),
                                "Failed to write the writable overlap set" );
            }
#endif
#endif
 
            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
            if( !runCtx->skip_io )
            {
                MB_CHK_SET_ERR( mbCore->write_file( sstr.str().c_str(), nullptr, writeOptions, &writableOverlapSet, 1 ),
                                "Failed to write the writable overlap set" );
            }
        }
 
        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" );
            MB_CHK_SET_ERR( 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
                                ),
                            "Failed to generate remapping weights" );
            runCtx->timer_pop();
 
            weightMap->PrintMapStatistics();
 
            // Invoke the CheckMap routine on the TempestRemap serial interface directly, if running
            // on a single process
            if( runCtx->fCheck )
            {
                const double dNormalTolerance = 1.0E-8;
                const double dStrictTolerance = 1.0E-12;
                weightMap->CheckMap( runCtx->fCheck, runCtx->fCheck, runCtx->fCheck && ( runCtx->ensureMonotonicity ),
                                     dNormalTolerance, dStrictTolerance, dTotalOverlapArea );
            }
 
            if( runCtx->outFilename.size() && !runCtx->skip_io )
            {
                std::map< std::string, std::string > attrMap;
                attrMap["MOABversion"]   = std::string( MOAB_PACKAGE_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];
                if( runCtx->intxFilename.size() ) 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
                // in extra case; maybe need a better solution, just create it with the right meshset
                // from the beginning;
                MB_CHK_SET_ERR( weightMap->WriteParallelMap( runCtx->outFilename.c_str(), attrMap ),
                                "Failed writing the parallel map to disk" );
            }
 
            if( runCtx->verifyWeights )
            {
                // Let us pick a sampling test function for solution evaluation
                // SH, SV, FH, C, USERVAR
                bool userVariable = false;
                moab::TempestOnlineMap::sample_function testFunction;
                if( !runCtx->variableToVerify.compare( "SH" ) )
                    testFunction = &sample_slow_harmonic;
                else if( !runCtx->variableToVerify.compare( "FH" ) )
                    testFunction = &sample_fast_harmonic;
                else if( !runCtx->variableToVerify.compare( "SV" ) )
                    testFunction = &sample_stationary_vortex;
                else if( !runCtx->variableToVerify.compare( "C" ) )
                    testFunction = &sample_constant;
                else
                {
                    userVariable = runCtx->variableToVerify.size() ? true : false;
                    testFunction = runCtx->variableToVerify.size() ? nullptr : sample_stationary_vortex;
                }
 
                moab::Tag srcAnalyticalFunction;
                moab::Tag tgtAnalyticalFunction;
                moab::Tag tgtProjectedFunction;
                if( testFunction )
                {
                    runCtx->timer_push( "describe a solution on source grid" );
                    // MB_CHK_SET_ERR( mbCore->tag_get_handle( runCtx->variableToVerify.c_str(), srcAnalyticalFunction ),
                    //                 "Failed to get analytical solution on source grid" );
                    MB_CHK_SET_ERR( weightMap->DefineAnalyticalSolution( srcAnalyticalFunction,
                                                                         "AnalyticalSolnSrcExact",
                                                                         moab::Remapper::SourceMesh, testFunction ),
                                    "Failed to define analytical solution on source grid" );
                    runCtx->timer_pop();
 
                    // runCtx->timer_push( "exchange solution on source grid" );
                    // moab::Range& srccovEnts = remapper.GetMeshEntities( moab::Remapper::CoveringMesh );
                    // MB_CHK_SET_ERR( pcomm->exchange_tags( srcAnalyticalFunction, srccovEnts ),
                    //                 "Failed to exchange analytical solution on source grid" );
                    // runCtx->timer_pop();
 
                    runCtx->timer_push( "describe a solution on target grid" );
                    MB_CHK_SET_ERR( weightMap->DefineAnalyticalSolution(
                                        tgtAnalyticalFunction, "AnalyticalSolnTgtExact", moab::Remapper::TargetMesh,
                                        testFunction, &tgtProjectedFunction, "ProjectedSolnTgt" ),
                                    "Failed to define analytical solution on target grid" );
                    runCtx->timer_pop();
                }
                else
                {
                    MB_CHK_SET_ERR( mbCore->tag_get_handle( runCtx->variableToVerify.c_str(), srcAnalyticalFunction ),
                                    "Failed to get analytical solution on source grid" );
                    MB_CHK_SET_ERR( mbCore->tag_get_handle( "ProjectedSolnTgt", 1, moab::MB_TYPE_DOUBLE,
                                                            tgtProjectedFunction,
                                                            moab::MB_TAG_DENSE | moab::MB_TAG_CREAT ),
                                    "Failed to get projected solution on target grid" );
                }
 
                // if( !runCtx->skip_io )
                {
                    MB_CHK_SET_ERR( mbCore->write_file( "srcWithSolnTag.h5m", nullptr, writeOptions,
                                                        &runCtx->meshsets[0], 1 ),
                                    "Failed to write the source mesh with solution tag" );
                }
 
                runCtx->timer_push( "compute solution projection on target grid" );
                MB_CHK_SET_ERR( weightMap->ApplyWeights( srcAnalyticalFunction, tgtProjectedFunction, false,
                                                         runCtx->cassType ),
                                "Failed to apply weights" );
                runCtx->timer_pop();
 
                // if( !runCtx->skip_io )
                {
                    MB_CHK_SET_ERR( mbCore->write_file( "tgtWithSolnTag2.h5m", nullptr, writeOptions,
                                                        &runCtx->meshsets[1], 1 ),
                                    "Failed to write the target mesh with projected solution tag" );
                }
 
                if( nprocs == 1 && runCtx->baselineFile.size() )
                {
                    // save the field from tgtWithSolnTag2 in a text file, and global ids for cells
                    moab::Range tgtEntities;
                    if( runCtx->disc_methods[1] == "pcloud" )
                    {
                        MB_CHK_SET_ERR( mbCore->get_entities_by_dimension( runCtx->meshsets[1], 0, tgtEntities ),
                                        "Failed to get entities by dimension" );
                    }
                    else
                    {
                        MB_CHK_SET_ERR( mbCore->get_entities_by_dimension( runCtx->meshsets[1], 2, tgtEntities ),
                                        "Failed to get entities by dimension" );
                    }
                    std::vector< int > globIds( tgtEntities.size() );
                    std::vector< double > vals( tgtEntities.size() );
                    moab::Tag projTag;
                    MB_CHK_SET_ERR( mbCore->tag_get_handle( "ProjectedSolnTgt", projTag ),
                                    "Failed to get projected solution tag" );
                    moab::Tag gid = mbCore->globalId_tag();
                    MB_CHK_SET_ERR( mbCore->tag_get_data( gid, tgtEntities, &globIds[0] ), "Failed to get global ids" );
                    MB_CHK_SET_ERR( mbCore->tag_get_data( projTag, tgtEntities, &vals[0] ),
                                    "Failed to get projected solution" );
                    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 < tgtEntities.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
                    if( !runCtx->skip_io )
                    {
                        MB_CHK_SET_ERR( mbCore->write_file( "srcWithSolnTag.h5m", nullptr, writeOptions,
                                                            &runCtx->meshsets[0], 1 ),
                                        "Failed to write the source mesh with solution tag" );
                    }
                }
 
                // compute error metrics if it is a known analytical functional
                if( !userVariable )
                {
                    runCtx->timer_push( "compute error metrics against analytical solution on target grid" );
                    std::map< std::string, double > errMetrics;
                    MB_CHK_SET_ERR( weightMap->ComputeMetrics( moab::Remapper::TargetMesh, tgtAnalyticalFunction,
                                                               tgtProjectedFunction, errMetrics, true ),
                                    "Failed to compute error metrics" );
                    runCtx->timer_pop();
                }
            }
 
            delete weightMap;
        }
    }
 
    // Clean up
    remapper.clear();
    delete runCtx;
    delete mbCore;
 
#ifdef MOAB_HAVE_MPI
    MPI_Finalize();
#endif
    return 0;
    }
    catch( const std::exception& e )
    {
        std::cerr << "[mbtempest] Fatal error: " << e.what() << std::endl;
#ifdef MOAB_HAVE_MPI
        MPI_Abort( MPI_COMM_WORLD, 1 );
#endif
        return 1;
    }
    catch( ... )
    {
        std::cerr << "[mbtempest] Fatal: unknown exception caught" << std::endl;
#ifdef MOAB_HAVE_MPI
        MPI_Abort( MPI_COMM_WORLD, 1 );
#endif
        return 1;
    }
}

References moab::Core::add_entities(), moab::TempestOnlineMap::ApplyWeights(), moab::IntxAreaUtils::area_on_sphere(), ToolContext::areaMethod, ToolContext::baselineFile, moab::Range::begin(), ToolContext::boxeps, ToolContext::cassType, moab::ParallelComm::check_all_shared_handles(), moab::TempestRemapper::clear(), moab::TempestOnlineMap::ComputeMetrics(), moab::TempestRemapper::ComputeOverlapMesh(), ToolContext::computeWeights, moab::TempestRemapper::ConstructCoveringSet(), moab::TempestRemapper::constructEdgeMap, moab::TempestRemapper::ConvertTempestMesh(), moab::Core::create_meshset(), CreateTempestMesh(), moab::TempestOnlineMap::DefineAnalyticalSolution(), moab::Core::delete_entities(), ToolContext::disc_methods, ToolContext::disc_orders, ToolContext::doftag_names, moab::Range::end(), moab::IntxUtils::enforce_convexity(), ToolContext::enforceConvexity, ToolContext::ensureMonotonicity, ToolContext::epsrel, error(), ToolContext::fCheck, moab::Intx2Mesh::FindMaxEdges(), moab::IntxUtils::fix_degenerate_quads(), moab::TempestOnlineMap::GenerateRemappingWeights(), moab::Core::get_entities_by_dimension(), moab::TempestRemapper::GetMesh(), moab::TempestRemapper::GetMeshSet(), moab::TempestRemapper::GetOverlapAugmentedEntities(), moab::Core::globalId_tag(), ToolContext::inFilenames, moab::TempestRemapper::initialize(), moab::Range::insert(), moab::Intx2Mesh::intersect_meshes(), ToolContext::intxFilename, ToolContext::kdtreeSearch, ToolContext::mapOptions, MB_CHK_SET_ERR, MB_TAG_CREAT, MB_TAG_DENSE, MB_TYPE_DOUBLE, ToolContext::meshes, MESHSET_SET, ToolContext::meshsets, ToolContext::meshType, moab::TempestRemapper::meshValidate, MOAB_PACKAGE_VERSION, ToolContext::nlayers, ToolContext::outFilename, ToolContext::outputFormatter, moab::TempestRemapper::OVERLAP_MEMORY, moab::TempestRemapper::OVERLAP_MOAB, moab::Remapper::OverlapMesh, ToolContext::ParseCLOptions(), moab::IntxAreaUtils::positive_orientation(), ToolContext::print_diagnostics, moab::DebugOutput::printf(), moab::TempestOnlineMap::PrintMapStatistics(), ToolContext::proc_id, radius, ToolContext::rrmGrids, sample_constant(), sample_fast_harmonic(), sample_slow_harmonic(), sample_stationary_vortex(), moab::Intx2Mesh::set_box_error(), moab::Intx2Mesh::set_error_tolerance(), moab::Intx2MeshOnSphere::set_radius_destination_mesh(), moab::Intx2MeshOnSphere::set_radius_source_mesh(), moab::TempestRemapper::SetAreaMethod(), moab::TempestRemapper::SetRegionalMesh(), moab::Range::size(), ToolContext::skip_intersection, ToolContext::skip_io, moab::Remapper::SourceMesh, moab::subtract(), moab::Core::tag_get_data(), moab::Core::tag_get_handle(), moab::Remapper::TargetMesh, ToolContext::timer_pop(), ToolContext::timer_push(), ToolContext::ToolContext(), ToolContext::useGnomonicProjection, ToolContext::variableToVerify, ToolContext::verifyWeights, moab::Core::write_file(), moab::Core::write_mesh(), and moab::TempestOnlineMap::WriteParallelMap().

sample_constant()

static constexpr double sample_constant ( double  , double    )

inlinestaticconstexprnoexcept

Constant sample function.

Returns

double Always returns 1.0

Definition at line 2073 of file mbtempest.cpp.

{
    return 1.0;
}

Referenced by main().

sample_fast_harmonic()

static double sample_fast_harmonic ( double  dLon, double  dLat  )

inlinestaticnoexcept

Sample function with fast harmonic variation.

Parameters

dLon	Longitude in radians
dLat	Latitude in radians

Returns

double Function value at (dLon, dLat)

Definition at line 2098 of file mbtempest.cpp.

{
    const double sin2Lat = std::sin( 2.0 * dLat );
    return 2.0 +
           sin2Lat * sin2Lat * sin2Lat * sin2Lat * sin2Lat * sin2Lat * sin2Lat * sin2Lat * std::cos( 16.0 * dLon );
}

Referenced by main().

sample_slow_harmonic()

static double sample_slow_harmonic ( double  dLon, double  dLat  )

inlinestaticnoexcept

Sample function with slow harmonic variation.

Parameters

dLon	Longitude in radians
dLat	Latitude in radians

Returns

double Function value at (dLon, dLat)

Definition at line 2085 of file mbtempest.cpp.

{
    const double cosLat = std::cos( dLat );
    return 2.0 + cosLat * cosLat * std::cos( 2.0 * dLon );
}

Referenced by main().

sample_stationary_vortex()

static double sample_stationary_vortex ( double  dLon, double  dLat  )

inlinestaticnoexcept

Sample function representing a stationary vortex.

Parameters

dLon	Longitude in radians
dLat	Latitude in radians

Returns

double Function value at (dLon, dLat)

Definition at line 2112 of file mbtempest.cpp.

{
    // Find the rotated longitude and latitude of a point on a sphere
    // with pole at (dLonC, dLatC)
    const double dSinC = std::sin( VORTEX_LAT0 );
    const double dCosC = std::cos( VORTEX_LAT0 );
    const double dSinT = std::sin( dLat );
    const double dCosT = std::cos( dLat );
 
    const double dTrm = dCosT * std::cos( dLon - VORTEX_LON0 );
    const double dX   = dSinC * dTrm - dCosC * dSinT;
    const double dY   = dCosT * std::sin( dLon - VORTEX_LON0 );
    const double dZ   = dSinC * dSinT + dCosC * dTrm;
 
    // Calculate new longitude and latitude in rotated coordinate system
    double dNewLon = std::atan2( dY, dX );
    if( dNewLon < 0.0 )
    {
        dNewLon += 2.0 * M_PI;
    }
    const double dNewLat = std::asin( dZ );
 
    // Calculate vortex profile
    const double dRho = VORTEX_R0 * std::cos( dNewLat );
    const double dVt  = 3.0 * std::sqrt( 3.0 ) / 2.0 / std::cosh( dRho ) / std::cosh( dRho ) * std::tanh( dRho );
 
    // Calculate angular velocity (avoid division by zero)
    const double dOmega = ( dRho == 0.0 ) ? 0.0 : ( dVt / dRho );
 
    // Return the final vortex profile
    return ( 1.0 - std::tanh( dRho / VORTEX_D * std::sin( dNewLon - dOmega * VORTEX_T ) ) );
}

References anonymous_namespace{mbtempest.cpp}::VORTEX_D, anonymous_namespace{mbtempest.cpp}::VORTEX_LAT0, anonymous_namespace{mbtempest.cpp}::VORTEX_LON0, anonymous_namespace{mbtempest.cpp}::VORTEX_R0, and anonymous_namespace{mbtempest.cpp}::VORTEX_T.

Referenced by main().