TempestRemapper.hpp

Go to the documentation of this file.

/*
 * =====================================================================================
 *
 * Filename: TempestRemapper.hpp
 *
 * Description: Interface to the TempestRemap library to enable intersection and
 * high-order conservative remapping of climate solution from
 * arbitrary resolution of source and target grids on the sphere.
 *
 * Author: Vijay S. Mahadevan (vijaysm), [email protected]
 *
 * =====================================================================================
 */
 
#ifndef MB_TEMPESTREMAPPER_HPP
#define MB_TEMPESTREMAPPER_HPP
 
#include "moab/Remapping/Remapper.hpp"
#include "moab/IntxMesh/Intx2MeshOnSphere.hpp"
#include "moab/IntxMesh/IntxUtils.hpp"
 
// Tempest includes
#ifdef MOAB_HAVE_TEMPESTREMAP
#include "netcdfcpp.h"
#include "TempestRemapAPI.h"
#else
#error "This tool depends on TempestRemap library. Reconfigure using --with-tempestremap"
#endif
 
namespace moab
{
 
// Forward declare our friend, the mapper
class TempestOnlineMap;
 
class TempestRemapper : public Remapper
{
 public:
#ifdef MOAB_HAVE_MPI
 TempestRemapper( moab::Interface* mbInt, moab::ParallelComm* pcomm = NULL, bool offlineMode = false )
 : Remapper( mbInt, pcomm ),
#else
 TempestRemapper( moab::Interface* mbInt, bool offlineMode = false )
 : Remapper( mbInt ),
#endif
 offlineWorkflow( offlineMode ), meshValidate( false ), constructEdgeMap( false ), m_source_type( DEFAULT ),
 m_target_type( DEFAULT )
 {
 }
 
 virtual ~TempestRemapper();
 
 // Mesh type with a correspondence to Tempest/Climate formats
 enum TempestMeshType
 {
 DEFAULT = -1,
 CS = 0,
 RLL = 1,
 ICO = 2,
 ICOD = 3,
 OVERLAP_FILES = 4,
 OVERLAP_MEMORY = 5,
 OVERLAP_MOAB = 6
 };
 
 friend class TempestOnlineMap;
 
 public:
 /// <summary>
 /// Initialize the TempestRemapper object internal datastructures including the mesh sets
 /// and TempestRemap mesh references.
 /// </summary>
 virtual ErrorCode initialize( bool initialize_fsets = true );
 
 /// <summary>
 /// Deallocate and clear any memory initialized in the TempestRemapper object
 /// </summary>
 virtual ErrorCode clear();
 
 /// <summary>
 /// Generate a mesh in memory of given type (CS/RLL/ICO/MPAS(structured)) and store it
 /// under the context specified by the user.
 /// </summary>
 moab::ErrorCode GenerateMesh( Remapper::IntersectionContext ctx, TempestMeshType type );
 
 /// <summary>
 /// Load a mesh from disk of given type and store it under the context specified by the
 /// user.
 /// </summary>
 moab::ErrorCode LoadMesh( Remapper::IntersectionContext ctx, std::string inputFilename, TempestMeshType type );
 
 /// <summary>
 /// Construct a source covering mesh such that it completely encompasses the target grid in
 /// parallel. This operation is critical to ensure that the parallel advancing-front
 /// intersection algorithm can the intersection mesh only locally without any process
 /// communication.
 /// </summary>
 moab::ErrorCode ConstructCoveringSet( double tolerance = 1e-8,
 double radius_src = 1.0,
 double radius_tgt = 1.0,
 double boxeps = 0.1,
 bool regional_mesh = false,
 bool gnomonic = true,
 int order = 1 );
 
 /// <summary>
 /// Compute the intersection mesh between the source and target grids that have been
 /// instantiated in the Remapper. This function invokes the parallel advancing-front
 /// intersection algorithm internally for spherical meshes and can handle arbitrary
 /// unstructured grids (CS, RLL, ICO, MPAS) with and without holes.
 /// </summary>
 moab::ErrorCode ComputeOverlapMesh( bool kdtree_search = true, bool use_tempest = false, int nLayers = 0 );
 
 /* Converters between MOAB and Tempest representations */
 
 /// <summary>
 /// Convert the TempestRemap mesh object to a corresponding MOAB mesh representation
 /// according to the intersection context.
 /// </summary>
 moab::ErrorCode ConvertTempestMesh( Remapper::IntersectionContext ctx );
 
 /// <summary>
 /// Convert the MOAB mesh representation to a corresponding TempestRemap mesh object
 /// according to the intersection context.
 /// </summary>
 moab::ErrorCode ConvertMeshToTempest( Remapper::IntersectionContext ctx );
 
 /// <summary>
 /// Get the TempestRemap mesh object according to the intersection context.
 /// </summary>
 Mesh* GetMesh( Remapper::IntersectionContext ctx );
 
 /// <summary>
 /// Set the TempestRemap mesh object according to the intersection context.
 /// </summary>
 void SetMesh( Remapper::IntersectionContext ctx, Mesh* mesh, bool overwrite = true );
 
 void SetMeshSet( Remapper::IntersectionContext ctx /* Remapper::CoveringMesh*/,
 moab::EntityHandle mset,
 moab::Range* entities = nullptr );
 /// <summary>
 /// Get the covering mesh (TempestRemap) object.
 /// </summary>
 Mesh* GetCoveringMesh();
 
 /// <summary>
 /// Get the MOAB mesh set corresponding to the intersection context.
 /// </summary>
 moab::EntityHandle& GetMeshSet( Remapper::IntersectionContext ctx );
 
 /// <summary>
 /// Const overload. Get the MOAB mesh set corresponding to the intersection context.
 /// </summary>
 moab::EntityHandle GetMeshSet( Remapper::IntersectionContext ctx ) const;
 
 /// <summary>
 /// Get the mesh element entities corresponding to the intersection context.
 /// </summary>
 moab::Range& GetMeshEntities( Remapper::IntersectionContext ctx );
 
 /// <summary>
 /// Const overload. Get the mesh element entities corresponding to the intersection context.
 /// </summary>
 const moab::Range& GetMeshEntities( Remapper::IntersectionContext ctx ) const;
 
 /// <summary>
 /// Get the mesh vertices corresponding to the intersection context. Useful for point-cloud
 /// meshes
 /// </summary>
 moab::Range& GetMeshVertices( Remapper::IntersectionContext ctx );
 
 /// <summary>
 /// Const overload. Get the mesh vertices corresponding to the intersection context. Useful
 /// for point-cloud meshes.
 /// </summary>
 const moab::Range& GetMeshVertices( Remapper::IntersectionContext ctx ) const;
 
 /// <summary>
 /// Get access to the underlying source covering set if available. Else return the source
 /// set.
 /// </summary>
 moab::EntityHandle& GetCoveringSet();
 
 /// <summary>
 /// Set the mesh type corresponding to the intersection context
 /// </summary>
 void SetMeshType( Remapper::IntersectionContext ctx, const std::vector< int >& metadata );
 
 /// <summary>
 /// reconstruct mesh, used now only for IO; need a better solution maybe
 /// </summary>
 void ResetMeshSet( Remapper::IntersectionContext ctx, moab::EntityHandle meshSet );
 
 /// <summary>
 /// Get the mesh type corresponding to the intersection context
 /// </summary>
 TempestMeshType GetMeshType( Remapper::IntersectionContext ctx ) const;
 
 /// <summary>
 /// Get the global ID corresponding to the local entity ID according to the context (source,
 /// target, intersection)
 /// </summary>
 int GetGlobalID( Remapper::IntersectionContext ctx, int localID );
 
 /// <summary>
 /// Get the local ID corresponding to the global entity ID according to the context (source,
 /// target, intersection)
 /// </summary>
 int GetLocalID( Remapper::IntersectionContext ctx, int globalID );
 
 /// <summary>
 /// Gather the overlap mesh and asssociated source/target data and write it out to disk
 /// using the TempestRemap output interface. This information can then be used with the
 /// "GenerateOfflineMap" tool in TempestRemap as needed.
 /// </summary>
 moab::ErrorCode WriteTempestIntersectionMesh( std::string strOutputFileName,
 const bool fAllParallel,
 const bool fInputConcave,
 const bool fOutputConcave );
 
 /// <summary>
 /// Generate the necessary metadata and specifically the GLL node numbering for DoFs for a
 /// CS mesh. This negates the need for running external code like HOMME to output the
 /// numbering needed for computing maps. The functionality is used through the `mbconvert`
 /// tool to compute processor-invariant Global DoF IDs at GLL nodes.
 /// </summary>
 moab::ErrorCode GenerateCSMeshMetadata( const int ntot_elements,
 moab::Range& entities,
 moab::Range* secondary_entities,
 const std::string& dofTagName,
 int nP );
 
 /// <summary>
 /// Generate the necessary metadata for DoF node numbering in a given mesh.
 /// Currently, only the functionality to generate numbering on CS grids is supported.
 /// </summary>
 moab::ErrorCode GenerateMeshMetadata( Mesh& mesh,
 const int ntot_elements,
 moab::Range& entities,
 moab::Range* secondary_entities,
 const std::string dofTagName,
 int nP );
 
 /// <summary>
 /// Compute the local and global IDs for elements in source/target/coverage meshes.
 /// </summary>
 moab::ErrorCode ComputeGlobalLocalMaps();
 
 /// <summary>
 /// Get all the ghosted overlap entities that were accumulated to enable conservation in
 /// parallel
 /// </summary>
 moab::ErrorCode GetOverlapAugmentedEntities( moab::Range& sharedGhostEntities );
 
#ifndef MOAB_HAVE_MPI
 /// <summary>
 /// Internal method to assign vertex and element global IDs if one does not exist already
 /// </summary>
 moab::ErrorCode assign_vertex_element_IDs( Tag idtag,
 EntityHandle this_set,
 const int dimension = 2,
 const int start_id = 1 );
#endif
 // <summary>
 /// Get the masks that could have been defined
 /// </summary>
 ErrorCode GetIMasks( Remapper::IntersectionContext ctx, std::vector< int >& masks );
 
 public: // public members
 const bool offlineWorkflow; // check whether we are in an offline workflow context (mbtempest)
 bool meshValidate; // Validate the mesh after loading from file
 
 bool constructEdgeMap; // Construct the edge map within the TempestRemap datastructures
 
 static const bool verbose = true;
 
 private:
 moab::ErrorCode convert_overlap_mesh_sorted_by_source();
 
 // private methods
 moab::ErrorCode load_tempest_mesh_private( std::string inputFilename, Mesh** tempest_mesh );
 
 moab::ErrorCode convert_mesh_to_tempest_private( Mesh* mesh,
 moab::EntityHandle meshset,
 moab::Range& entities,
 moab::Range* pverts );
 
 moab::ErrorCode convert_tempest_mesh_private( TempestMeshType type,
 Mesh* mesh,
 moab::EntityHandle& meshset,
 moab::Range& entities,
 moab::Range* vertices );
 
 moab::ErrorCode augment_overlap_set();
 
 /* Source meshset, mesh and entity references */
 Mesh* m_source;
 TempestMeshType m_source_type;
 moab::Range m_source_entities;
 moab::Range m_source_vertices;
 moab::EntityHandle m_source_set;
 int max_source_edges;
 bool point_cloud_source;
 std::vector< int > m_source_metadata;
 
 /* Target meshset, mesh and entity references */
 Mesh* m_target;
 TempestMeshType m_target_type;
 moab::Range m_target_entities;
 moab::Range m_target_vertices;
 moab::EntityHandle m_target_set;
 int max_target_edges;
 bool point_cloud_target;
 std::vector< int > m_target_metadata;
 
 /* Overlap meshset, mesh and entity references */
 Mesh* m_overlap;
 TempestMeshType m_overlap_type;
 moab::Range m_overlap_entities;
 moab::EntityHandle m_overlap_set;
 std::vector< std::pair< int, int > > m_sorted_overlap_order;
 
 /* Intersection context on a sphere */
 moab::Intx2MeshOnSphere* mbintx;
 
 /* Parallel - migrated mesh that is in the local view */
 Mesh* m_covering_source;
 moab::EntityHandle m_covering_source_set;
 moab::Range m_covering_source_entities;
 moab::Range m_covering_source_vertices;
 
 /* local to glboal and global to local ID maps */
 std::map< int, int > gid_to_lid_src, gid_to_lid_covsrc, gid_to_lid_tgt;
 std::map< int, int > lid_to_gid_src, lid_to_gid_covsrc, lid_to_gid_tgt;
 
 IntxAreaUtils::AreaMethod m_area_method;
 
 // this is a copy of the original source/target set,
 // with some ghosts created for parallel cases
 moab::EntityHandle m_source_set_with_ghosts; // source+ghost
 moab::EntityHandle m_target_set_with_ghosts; // target+ghost
 
 bool rrmgrids;
 bool is_parallel, is_root;
 int rank, size;
};
 
// Inline functions
inline Mesh* TempestRemapper::GetMesh( Remapper::IntersectionContext ctx )
{
 switch( ctx )
 {
 case Remapper::SourceMesh:
 return m_source;
 case Remapper::TargetMesh:
 return m_target;
 case Remapper::OverlapMesh:
 return m_overlap;
 case Remapper::CoveringMesh:
 return m_covering_source;
 case Remapper::DEFAULT: // Remapper::SourceMeshWithGhosts does not need one
 default:
 return NULL;
 }
}
 
inline void TempestRemapper::SetMesh( Remapper::IntersectionContext ctx, Mesh* mesh, bool overwrite )
{
 switch( ctx )
 {
 case Remapper::SourceMesh:
 if( !overwrite && m_source ) return;
 if( overwrite && m_source ) delete m_source;
 m_source = mesh;
 break;
 case Remapper::TargetMesh:
 if( !overwrite && m_target ) return;
 if( overwrite && m_target ) delete m_target;
 m_target = mesh;
 break;
 case Remapper::OverlapMesh:
 if( !overwrite && m_overlap ) return;
 if( overwrite && m_overlap ) delete m_overlap;
 m_overlap = mesh;
 break;
 case Remapper::CoveringMesh:
 if( !overwrite && m_covering_source ) return;
 if( overwrite && m_covering_source ) delete m_covering_source;
 m_covering_source = mesh;
 break;
 case Remapper::DEFAULT: // Remapper::SourceMeshWithGhosts does not need one
 default:
 break;
 }
}
 
inline void TempestRemapper::ResetMeshSet( Remapper::IntersectionContext ctx, moab::EntityHandle meshSet )
{
 switch( ctx )
 {
 case Remapper::SourceMesh:
 delete m_source;
 m_source = new Mesh;
 m_source_set = meshSet;
 convert_mesh_to_tempest_private( m_source, m_source_set, m_source_entities, &m_source_vertices );
 m_source->CalculateFaceAreas( false ); // fInputConcave is false ?
 break;
 case Remapper::TargetMesh:
 // not needed yet
 break;
 case Remapper::OverlapMesh:
 // not needed yet
 break;
 case Remapper::CoveringMesh:
 // not needed yet
 break;
 case Remapper::DEFAULT:
 default:
 break;
 }
}
 
 
inline moab::EntityHandle& TempestRemapper::GetMeshSet( Remapper::IntersectionContext ctx )
{
 switch( ctx )
 {
 case Remapper::SourceMesh:
 return m_source_set;
 case Remapper::TargetMesh:
 return m_target_set;
 case Remapper::OverlapMesh:
 return m_overlap_set;
 case Remapper::CoveringMesh:
 return m_covering_source_set;
 case Remapper::SourceMeshWithGhosts:
 return m_source_set_with_ghosts;
 case Remapper::TargetMeshWithGhosts:
 return m_target_set_with_ghosts;
 case Remapper::DEFAULT:
 default:
 MB_SET_ERR_RET_VAL( "Invalid context passed to GetMeshSet", m_overlap_set );
 }
}
 
inline moab::EntityHandle TempestRemapper::GetMeshSet( Remapper::IntersectionContext ctx ) const
{
 switch( ctx )
 {
 case Remapper::SourceMesh:
 return m_source_set;
 case Remapper::TargetMesh:
 return m_target_set;
 case Remapper::OverlapMesh:
 return m_overlap_set;
 case Remapper::CoveringMesh:
 return m_covering_source_set;
 case Remapper::SourceMeshWithGhosts:
 return m_source_set_with_ghosts;
 case Remapper::TargetMeshWithGhosts:
 return m_target_set_with_ghosts;
 case Remapper::DEFAULT:
 default:
 MB_SET_ERR_RET_VAL( "Invalid context passed to GetMeshSet", m_overlap_set );
 }
}
 
inline moab::Range& TempestRemapper::GetMeshEntities( Remapper::IntersectionContext ctx )
{
 switch( ctx )
 {
 case Remapper::SourceMesh:
 return m_source_entities;
 case Remapper::TargetMesh:
 return m_target_entities;
 case Remapper::OverlapMesh:
 return m_overlap_entities;
 case Remapper::CoveringMesh:
 return m_covering_source_entities;
 case Remapper::DEFAULT:
 default:
 MB_SET_ERR_RET_VAL( "Invalid context passed to GetMeshSet", m_overlap_entities );
 }
}
 
inline const moab::Range& TempestRemapper::GetMeshEntities( Remapper::IntersectionContext ctx ) const
{
 switch( ctx )
 {
 case Remapper::SourceMesh:
 return m_source_entities;
 case Remapper::TargetMesh:
 return m_target_entities;
 case Remapper::OverlapMesh:
 return m_overlap_entities;
 case Remapper::CoveringMesh:
 return m_covering_source_entities;
 case Remapper::DEFAULT:
 default:
 MB_SET_ERR_RET_VAL( "Invalid context passed to GetMeshSet", m_overlap_entities );
 }
}
 
inline moab::Range& TempestRemapper::GetMeshVertices( Remapper::IntersectionContext ctx )
{
 switch( ctx )
 {
 case Remapper::SourceMesh:
 return m_source_vertices;
 case Remapper::TargetMesh:
 return m_target_vertices;
 case Remapper::CoveringMesh:
 return m_covering_source_vertices;
 case Remapper::DEFAULT:
 default:
 MB_SET_ERR_RET_VAL( "Invalid context passed to GetMeshSet", m_source_vertices );
 }
}
 
inline const moab::Range& TempestRemapper::GetMeshVertices( Remapper::IntersectionContext ctx ) const
{
 switch( ctx )
 {
 case Remapper::SourceMesh:
 return m_source_vertices;
 case Remapper::TargetMesh:
 return m_target_vertices;
 case Remapper::CoveringMesh:
 return m_covering_source_vertices;
 case Remapper::DEFAULT:
 default:
 MB_SET_ERR_RET_VAL( "Invalid context passed to GetMeshSet", m_source_vertices );
 }
}
 
inline void TempestRemapper::SetMeshType( Remapper::IntersectionContext ctx, const std::vector< int >& metadata )
{
 switch( ctx )
 {
 case Remapper::SourceMesh:
 m_source_type = static_cast< moab::TempestRemapper::TempestMeshType >( metadata[0] );
 if( metadata[0] == 1 ) // RLL mesh
 {
 m_source_metadata.resize( 2 );
 m_source_metadata[0] = metadata[1];
 m_source_metadata[1] = metadata[2];
 }
 else
 {
 m_source_metadata.resize( 1 );
 m_source_metadata[0] = metadata[1];
 }
 break;
 case Remapper::TargetMesh:
 m_target_type = static_cast< moab::TempestRemapper::TempestMeshType >( metadata[0] );
 if( metadata[0] == 1 ) // RLL mesh
 {
 m_target_metadata.resize( 2 );
 m_target_metadata[0] = metadata[1];
 m_target_metadata[1] = metadata[2];
 }
 else
 {
 m_target_metadata.resize( 1 );
 m_target_metadata[0] = metadata[1];
 }
 break;
 case Remapper::OverlapMesh:
 m_overlap_type = moab::TempestRemapper::OVERLAP_FILES;
 default:
 break;
 }
}
 
inline TempestRemapper::TempestMeshType TempestRemapper::GetMeshType( Remapper::IntersectionContext ctx ) const
{
 switch( ctx )
 {
 case Remapper::SourceMesh:
 return m_source_type;
 case Remapper::TargetMesh:
 return m_target_type;
 case Remapper::OverlapMesh:
 return m_overlap_type;
 case Remapper::DEFAULT: // not need yet case Remapper::SourceMeshWithGhosts
 default:
 return TempestRemapper::DEFAULT;
 }
}
 
inline Mesh* TempestRemapper::GetCoveringMesh()
{
 return m_covering_source;
}
 
inline moab::EntityHandle& TempestRemapper::GetCoveringSet()
{
 return m_covering_source_set;
}
 
inline int TempestRemapper::GetGlobalID( Remapper::IntersectionContext ctx, int localID )
{
 switch( ctx )
 {
 case Remapper::SourceMesh:
 return lid_to_gid_src[localID];
 case Remapper::TargetMesh:
 return lid_to_gid_tgt[localID];
 case Remapper::CoveringMesh:
 return lid_to_gid_covsrc[localID];
 case Remapper::OverlapMesh:
 case Remapper::SourceMeshWithGhosts:
 case Remapper::DEFAULT:
 default:
 return -1;
 }
}
 
inline int TempestRemapper::GetLocalID( Remapper::IntersectionContext ctx, int globalID )
{
 switch( ctx )
 {
 case Remapper::SourceMesh:
 return gid_to_lid_src[globalID];
 case Remapper::TargetMesh:
 return gid_to_lid_tgt[globalID];
 case Remapper::CoveringMesh:
 return gid_to_lid_covsrc[globalID];
 case Remapper::DEFAULT:
 case Remapper::OverlapMesh:
 case Remapper::SourceMeshWithGhosts:
 default:
 return -1;
 }
}
 
} // namespace moab
 
#endif // MB_TEMPESTREMAPPER_HPP