Mesh Oriented datABase  (version 5.6.0)
An array-based unstructured mesh library
mbtempest.cpp
Go to the documentation of this file.
1 /**
2  * @file mbtempest.cpp
3  * @brief MOAB-Tempest: A powerful mesh generation and remapping tool for climate and weather applications
4  *
5  * @section overview Overview
6  * MOAB-Tempest is a command-line tool that provides mesh generation and conservative remapping capabilities
7  * for climate and weather modeling. It combines the power of MOAB (Mesh-Oriented datABase) with the
8  * TempestRemap library to enable high-performance, parallel mesh generation and remapping operations.
9  *
10  * @section features Key Features
11  * - Generation of various spherical mesh types (Cubed-Sphere, RLL, Icosahedral)
12  * - Support for high-order discretization methods (FV, CGLL, DGLL)
13  * - Conservative remapping between different mesh types
14  * - Parallel processing support via MPI
15  * - Flexible I/O with support for multiple file formats
16  * - Built-in analytical functions for testing and validation
17  *
18  * @section algorithms Supported Algorithms
19  * - Mesh Generation:
20  * - Cubed-Sphere (CS) meshes
21  * - Regular Latitude-Longitude (RLL) meshes
22  * - Icosahedral (ICO) meshes
23  * - Overlap meshes for remapping
24  * - Remapping Methods:
25  * - Finite Volume (FV)
26  * - Continuous Galerkin (CGLL)
27  * - Discontinuous Galerkin (DGLL)
28  * - Monotonic and high-order variants
29  *
30  * @section usage Basic Usage Examples
31  * @code
32  * # Generate a Cubed-Sphere mesh with resolution 25
33  * ./mbtempest --type 0 --res 25 --file cubed_sphere_mesh.h5m
34  *
35  * # Generate a RLL mesh with resolution 90x180 (lon x lat)
36  * ./mbtempest --type 1 --res 90 --file rll_mesh.h5m
37  *
38  * # Generate an Icosahedral mesh with resolution 25 (dual mesh)
39  * ./mbtempest --type 2 --res 25 --dual --file icosahedral_dual_mesh.h5m
40  *
41  * # Compute overlap between two meshes
42  * ./mbtempest --type 5 --load mesh1.h5m --load mesh2.h5m intx intersection_mesh.h5m
43  *
44  * # Generate a remapping weights file between two meshes: FV to FV (default)
45  * ./mbtempest --type 5 --load source_mesh.h5m --load target_mesh.h5m --file weights.nc
46  *
47  * # Generate remapping weights file between two meshes: making it explicit (SE to FV)
48  * ./mbtempest --type 5 --load source_mesh.h5m --load target_mesh.h5m \
49  * --order 4 --method cgll --global_id GLOBAL_DOFS \
50  * --order 1 --method fv --global_id GLOBAL_ID \
51  * --file weights_se_to_fv.nc
52  * @endcode
53  *
54  * @section options Command Line Options
55  * Run './mbtempest --help' for a complete list of available options.
56  *
57  * @section notes Notes
58  * - For parallel execution, use MPI launcher (e.g., mpirun, mpiexec)
59  * - Output formats: .h5m (MOAB), .nc (NetCDF), .exo (ExodusII)
60  * - Requires MOAB and TempestRemap libraries
61  *
62  * @author MOAB Development Team
63  * @date Created: 2023
64  */
65 
66 // standard C++ includes
67 #include <iostream>
68 #include <iomanip>
69 #include <cstdlib>
70 #include <vector>
71 #include <string>
72 #include <memory>
73 #include <sstream>
74 #include <cassert>
75 
76 // MOAB includes
77 #include "moab/Core.hpp"
78 #include "moab/IntxMesh/IntxUtils.hpp"
79 #include "moab/Remapping/TempestRemapper.hpp"
80 #include "moab/Remapping/TempestOnlineMap.hpp"
81 #include "moab/ProgOptions.hpp"
82 #include "moab/CpuTimer.hpp"
83 #include "DebugOutput.hpp"
84 
85 #ifdef MOAB_HAVE_MPI
86 // MPI includes
87 #include "moab_mpi.h"
88 #include "moab/ParallelComm.hpp"
89 #include "MBParallelConventions.h"
90 #endif
91 
92 /**
93  * @brief Context class for MOAB-TempestRemap tool configuration and state management
94  */
95 class ToolContext
96 {
97  public:
98  // Core components
99  moab::Core* const mbcore; ///< MOAB Core instance for mesh operations
100 #ifdef MOAB_HAVE_MPI
101  moab::ParallelComm* const pcomm; ///< Parallel communicator (nullptr in serial)
102 #endif
103  const int proc_id; ///< MPI process rank (0 for serial)
104  const int n_procs; ///< Total number of MPI processes (1 for serial)
105  moab::DebugOutput outputFormatter; ///< Formatter for debug output
106 
107  // Mesh and remapping configuration
108  moab::TempestRemapper::TempestMeshType meshType{ moab::TempestRemapper::DEFAULT }; ///< Type of mesh to generate
109  std::vector< std::string > inFilenames; ///< Input filenames for source and target meshes
110  std::vector< int > disc_orders; ///< Discretization orders for source and target
111  std::vector< std::string > disc_methods; ///< Discretization methods (fv, cgll, dgll) for source and target
112  std::vector< std::string > doftag_names; ///< Degree of freedom tag names for source and target
113  std::string outFilename{ "outputFile.nc" }; ///< Output filename for remapping results
114  std::string intxFilename; ///< Intersection mesh filename (optional)
115  std::string baselineFile; ///< Baseline file for verification (optional)
116  std::string variableToVerify; ///< Variable name for verification (optional)
117  std::string fvMethod{ "none" }; ///< Finite volume method specification
118 
119  // Remapping options
120  GenerateOfflineMapAlgorithmOptions mapOptions; ///< Configuration for offline map generation
121  moab::TempestOnlineMap::CAASType cassType{
122  moab::TempestOnlineMap::CAAS_NONE }; ///< Conservative and accurate advection scheme type
123  int ensureMonotonicity{ 0 }; ///< Monotonicity enforcement level (0=none, 1=basic, 2=full, 3=strict)
124  bool rrmGrids{ false }; ///< Flag to use RRM (Regional Refinement Meshes)
125  bool kdtreeSearch{ true }; ///< Enable KD-tree for spatial searches
126  bool fCheck{ false }; ///< Enable additional checking during remapping
127  bool fVolumetric{ false }; ///< Enable volumetric (3D) remapping
128  bool useGnomonicProjection{ false }; ///< Use gnomonic projection for certain operations
129  bool print_diagnostics{ false }; ///< Print detailed diagnostic information
130  bool skip_intersection{ false }; ///< Skip intersection computation (for debugging)
131  double boxeps{ 1e-7 }; ///< Epsilon for bounding box checks
132  double epsrel{ ReferenceTolerance }; ///< Relative tolerance for convergence
133 
134  // Mesh operations control
135  bool skip_io{ false }; ///< Skip file I/O operations (for testing)
136  bool computeDual{ false }; ///< Compute dual mesh
137  bool computeWeights{ false }; ///< Compute interpolation weights
138  bool verifyConservation{ false }; ///< Verify conservation properties
139  bool verifyWeights{ false }; ///< Verify interpolation weights
140  bool enforceConvexity{ false }; ///< Enforce convexity in mesh elements
141 
142  // Performance and debugging
143  std::unique_ptr< moab::CpuTimer > timer; ///< Timer for performance measurement
144  double timer_ops{ 0.0 }; ///< Operation timer value
145  std::string opName; ///< Name of current operation being timed
146  int nlayers{ 0 }; ///< Number of ghost layers for parallel operations
147  int blockSize{ 5 }; ///< Block size for vectorized operations
148 
149  // Mesh data
150  std::vector< Mesh* > meshes; ///< Collection of TempestRemap meshes
151  std::vector< moab::EntityHandle > meshsets; ///< MOAB entity sets for meshes
152 
153  /**
154  * @brief Construct a new ToolContext object with MPI support
155  * @param icore MOAB Core instance (must not be null)
156  * @param p_pcomm Parallel communicator (must not be null in MPI mode)
157  * @throw std::invalid_argument if icore is null or p_pcomm is null in MPI mode
158  */
159 #ifdef MOAB_HAVE_MPI
160  ToolContext( moab::Core* icore, moab::ParallelComm* p_pcomm )
161  : mbcore( icore ), pcomm( p_pcomm ), proc_id( p_pcomm ? p_pcomm->rank() : 0 ),
162  n_procs( p_pcomm ? p_pcomm->size() : 1 ), outputFormatter( std::cout, p_pcomm ? p_pcomm->rank() : 0, 0 )
163  {
164  if( !icore ) throw std::invalid_argument( "MOAB Core instance cannot be null" );
165  if( !p_pcomm ) throw std::invalid_argument( "ParallelComm cannot be null in MPI mode" );
166 #else
167  /**
168  * @brief Construct a new ToolContext object (serial version)
169  * @param icore MOAB Core instance (must not be null)
170  * @throw std::invalid_argument if icore is null
171  */
172  explicit ToolContext( moab::Core* icore )
173  : mbcore( icore ), proc_id( 0 ), n_procs( 1 ), outputFormatter( std::cout, 0, 0 )
174  {
175 #endif
176  // Initialize default values
177  inFilenames.reserve( 2 );
178  doftag_names = { "GLOBAL_ID", "GLOBAL_ID" };
179  disc_orders = { 1, 1 };
180  disc_methods = { "fv", "fv" };
181 
182  // Initialize timer and output formatter
183  timer = std::make_unique< moab::CpuTimer >();
184  outputFormatter.set_prefix( "[MBTempest]: " );
185 
186  // Set default map options
187  mapOptions.fNoConservation = false;
188  mapOptions.fMonotone = false;
189  mapOptions.fNoCorrectAreas = false;
190  mapOptions.fNoCheck = false;
191  mapOptions.nPin = 1;
192  mapOptions.nPout = 1;
193  }
194 
195  // Rule of Five - Delete copy/move operations as mbcore is const
196  ~ToolContext() = default;
197  ToolContext( const ToolContext& ) = delete;
198  ToolContext& operator=( const ToolContext& ) = delete;
199  ToolContext( ToolContext&& ) = delete;
200  ToolContext& operator=( ToolContext&& ) = delete;
201 
202  /**
203  * @brief Start timing an operation
204  * @param operation Name of the operation being timed
205  */
206  void timer_push( const std::string& operation )
207  {
208  timer_ops = timer->time_since_birth();
209  opName = operation;
210  }
211 
212  /**
213  * @brief Stop timing and log the operation duration
214  */
215  void timer_pop()
216  {
217  double locElapsed = timer->time_since_birth() - timer_ops;
218  double avgElapsed = locElapsed;
219  double maxElapsed = locElapsed;
220 
221 #ifdef MOAB_HAVE_MPI
222  MPI_Reduce( &locElapsed, &maxElapsed, 1, MPI_DOUBLE, MPI_MAX, 0, pcomm->comm() );
223  MPI_Reduce( &locElapsed, &avgElapsed, 1, MPI_DOUBLE, MPI_SUM, 0, pcomm->comm() );
224  avgElapsed /= n_procs;
225 #endif
226 
227  if( proc_id == 0 )
228  {
229  std::cout << "[LOG] Time taken to " << opName << ": max = " << maxElapsed << ", avg = " << avgElapsed
230  << "\n";
231  }
232  opName.clear();
233  }
234 
235  /**
236  * @brief Parse command line arguments
237  * @param argc Argument count
238  * @param argv Argument values
239  * @return moab::ErrorCode indicating success or failure
240  * @throw std::invalid_argument for invalid command line arguments
241  */
242  moab::ErrorCode ParseCLOptions( int argc, char** argv )
243  {
244  // Initialize variables for command line options
245  int imeshType = 0;
246  std::string expectedFName = "output.exo";
247  std::string expectedMethod = "fv";
248  std::string expectedFVMethod = "none";
249  std::string expectedDofTagName = "GLOBAL_ID";
250  int expectedOrder = 1;
251  int useCAAS = 0;
252  int nlayer_input = -1; // -1 means not set by user
253  bool version_info = false;
254 
255  // Print command line for debugging
256  if( proc_id == 0 )
257  {
258  std::cout << "Command line options provided to mbtempest:\n ";
259  for( int i = 0; i < argc; ++i )
260  {
261  std::cout << argv[i] << " ";
262  }
263  std::cout << "\n" << std::endl;
264  }
265 
266  // Create options object with description
267  ProgOptions opts( "mbtempest - A mesh generation and remapping tool" );
268 
269  // Mesh generation options
270  opts.addOpt< int >( "type,t",
271  "Type of mesh (default=CS; Choose from [CS=0, RLL=1, ICO=2, OVERLAP_FILES=3, "
272  "OVERLAP_MEMORY=4, OVERLAP_MOAB=5])",
273  &imeshType );
274 
275  opts.addOpt< int >( "res,r", "Resolution of the mesh (default=5)", &blockSize );
276 
277  opts.addOpt< void >( "dual,d", "Output the dual of the mesh (relevant only for ICO mesh type)", &computeDual );
278 
279  opts.addOpt< std::string >( "file,f", "Output computed mesh or remapping weights to specified filename",
280  &outFilename );
281 
282  // Input/Output options
283  opts.addOpt< std::string >(
284  "load,l", "Input mesh filenames for source and target meshes. (relevant only when computing weights)",
285  &expectedFName );
286 
287  opts.addOpt< void >( "advfront,a",
288  "Use the advancing front intersection instead of the Kd-tree based algorithm to compute "
289  "mesh intersections.",
290  &kdtreeSearch );
291 
292  opts.addOpt< std::string >( "intx,i", "Output TempestRemap intersection mesh filename", &intxFilename );
293 
294  opts.addOpt< void >(
295  "weights,w",
296  "Compute and output the weights using the overlap mesh (generally relevant only for OVERLAP mesh)",
297  &computeWeights );
298 
299  // Discretization options
300  opts.addOpt< void >(
301  "verbose,v", "Print verbose diagnostic messages during intersection and map computation (default=false)",
302  &print_diagnostics );
303 
304  opts.addOpt< std::string >( "method,m", "Discretization method for the source and target solution fields",
305  &expectedMethod );
306 
307  opts.addOpt< int >( "order,o", "Discretization orders for the source and target solution fields",
308  &expectedOrder );
309 
310  opts.addOpt< std::string >( "global_id,g",
311  "Tag name that contains the global DoF IDs for source and target solution fields",
312  &expectedDofTagName );
313 
314  // Advanced options
315  opts.addOpt< std::string >( "fvmethod",
316  "Sub-type method for FV-FV projections (invdist, delaunay, bilin, intbilin, "
317  "intbilingb, none. Default: none)",
318  &expectedFVMethod );
319 
320  opts.addOpt< void >(
321  "noconserve", "Do not apply conservation to the resultant weights (relevant only when computing weights)",
322  &mapOptions.fNoConservation );
323 
324  opts.addOpt< void >(
325  "volumetric", "Apply a volumetric projection to compute the weights (relevant only when computing weights)",
326  &fVolumetric );
327 
328  opts.addOpt< void >( "skip_intersection", "Skip mesh intersection computation.", &skip_intersection );
329 
330  opts.addOpt< void >( "skip_output", "For performance studies, skip all I/O operations.", &skip_io );
331 
332  opts.addOpt< void >( "gnomonic", "Use Gnomonic plane projections to compute coverage mesh.",
333  &useGnomonicProjection );
334 
335  opts.addOpt< void >( "enforce_convexity", "Check convexity of input meshes to compute mesh intersections",
336  &enforceConvexity );
337 
338  opts.addOpt< void >( "nobubble", "Do not use bubble on interior of spectral element nodes",
339  &mapOptions.fNoBubble );
340 
341  opts.addOpt< void >(
342  "sparseconstraints",
343  "Use sparse solver for constraints when we have high-valence (typical with high-res RLL mesh)",
344  &mapOptions.fSparseConstraints );
345 
346  opts.addOpt< void >(
347  "rrmgrids",
348  "At least one of the meshes is a regionally refined grid (relevant to accelerate intersection computation)",
349  &rrmGrids );
350 
351  opts.addOpt< void >( "checkmap", "Check the generated map for conservation and consistency", &fCheck );
352 
353  opts.addOpt< void >( "verify",
354  "Verify the accuracy of the maps by projecting analytical functions from source to target "
355  "grid by applying the maps",
356  &verifyWeights );
357 
358  opts.addOpt< std::string >( "var",
359  "Tag name of the variable to use in the verification study (error metrics for user "
360  "defined variables may not be available)",
361  &variableToVerify );
362 
363  opts.addOpt< int >( "monotonicity", "Ensure monotonicity in the weight generation. Options=[0,1,2,3]",
364  &ensureMonotonicity );
365 
366  opts.addOpt< int >( "ghost",
367  "Number of ghost layers in coverage mesh (overrides automatic selection: 0 for FV order 1, "
368  "p+1 for FV order p>1)",
369  &nlayer_input );
370 
371  opts.addOpt< double >( "boxeps", "The tolerance for boxes (default=1e-7)", &boxeps );
372 
373  opts.addOpt< int >( "limiter", "Apply nonlinear filter after linear map application", &useCAAS );
374 
375  opts.addOpt< std::string >( "baseline", "Output baseline file", &baselineFile );
376 
377  opts.addOpt< void >( "manual", "Show documentation about usage with examples" );
378 
379  opts.addOpt< void >( "version", "Show version information", &version_info );
380 
381  // Parse command line
382  opts.parseCommandLine( argc, argv );
383 
384  // Handle call for detailed information
385  if( opts.numOptSet( "manual" ) > 0 )
386  {
387  if( this->proc_id == 0 )
388  {
389  this->printHelp( argv[0] );
390  }
391  exit( 0 );
392  }
393 
394  if( version_info )
395  {
396  if( this->proc_id == 0 )
397  {
398  std::cout << "mbtempest is part of the MOAB library version " << std::string( MOAB_PACKAGE_VERSION )
399  << "\n";
400  }
401  exit( 0 );
402  }
403 
404  // Process mesh type
405  switch( imeshType )
406  {
407  case 0:
408  this->meshType = moab::TempestRemapper::CS;
409  break;
410  case 1:
411  this->meshType = moab::TempestRemapper::RLL;
412  break;
413  case 2:
414  this->meshType = moab::TempestRemapper::ICO;
415  break;
416  case 3:
417  this->meshType = moab::TempestRemapper::OVERLAP_FILES;
418  break;
419  case 4:
420  this->meshType = moab::TempestRemapper::OVERLAP_MEMORY;
421  break;
422  case 5:
423  this->meshType = moab::TempestRemapper::OVERLAP_MOAB;
424  break;
425  default:
426  this->meshType = moab::TempestRemapper::DEFAULT;
427  break;
428  }
429 
430  // Process CAAS type
431  switch( useCAAS )
432  {
433  case 1:
434  this->cassType = moab::TempestOnlineMap::CAAS_GLOBAL;
435  break;
436  case 2:
437  this->cassType = moab::TempestOnlineMap::CAAS_LOCAL;
438  break;
439  case 3:
440  this->cassType = moab::TempestOnlineMap::CAAS_LOCAL_ADJACENT;
441  break;
442  case 4:
443  this->cassType = moab::TempestOnlineMap::CAAS_QLT;
444  break;
445  default:
446  this->cassType = moab::TempestOnlineMap::CAAS_NONE;
447  break;
448  }
449 
450  // Process input files if provided
451  if( !expectedFName.empty() )
452  {
453  this->inFilenames = { expectedFName };
454  }
455 
456  // Process discretization options through processMeshOptions to handle both single and multiple values
457  // Set initial defaults that can be overridden by processMeshOptions
458  this->fvMethod = expectedFVMethod;
459  this->disc_orders = { expectedOrder, expectedOrder };
460  this->disc_methods = { expectedMethod, expectedMethod };
461  this->doftag_names = { expectedDofTagName, expectedDofTagName };
462 
463  // Let processMeshOptions handle all the discretization option processing
464  this->processMeshOptions( opts );
465 
466  // Now use the processed values for map configuration
467  this->mapOptions.nPin = this->disc_orders[0];
468  this->mapOptions.nPout = this->disc_orders[1];
469  this->mapOptions.fSourceConcave = false;
470  this->mapOptions.fTargetConcave = false;
471  this->mapOptions.strMethod = "";
472 
473  // Configure map options with the processed values - this handles all remaining setup
474  this->configureMapOptions( nlayer_input );
475 
476  // Print runtime parameters
477  this->printRuntimeParameters();
478 
479  return moab::MB_SUCCESS;
480  }
481 
482  /**
483  * @brief Get the appropriate MOAB read options based on file extension and parallel configuration
484  *
485  * @param ctx Tool context containing parallel information
486  * @param filename Input filename to determine read options
487  * @return std::string MOAB read options string
488  */
489  std::string get_file_read_options( const std::string& filename )
490  {
491  // For serial execution, return default options
492  if( n_procs <= 1 )
493  {
494  return "";
495  }
496 
497  // Extract file extension
498  const size_t last_dot = filename.find_last_of( "." );
499  if( last_dot == std::string::npos )
500  {
501  return ""; // No extension found
502  }
503 
504  const std::string extension = filename.substr( last_dot + 1 );
505 
506  // Handle H5M files
507  if( extension == "h5m" )
508  {
509  return "PARALLEL=READ_PART;PARTITION=PARALLEL_PARTITION;PARALLEL_RESOLVE_SHARED_ENTS;";
510  }
511 
512  // Handle NetCDF files
513  if( extension == "nc" )
514  {
515  // Default NetCDF options
516 #ifdef MOAB_HAVE_ZOLTAN
517  std::string netcdf_options = "PARALLEL=READ_PART;PARTITION_METHOD=RCBZOLTAN;";
518 #else
519  std::string netcdf_options = "PARALLEL=READ_PART;PARTITION_METHOD=TRIVIAL;";
520 #endif
521  // Only rank 0 needs to determine the NetCDF file type
522  if( proc_id == 0 )
523  {
524  NcFile ncFile( filename.c_str(), NcFile::ReadOnly );
525  if( !ncFile.is_valid() )
526  {
527  // Handle invalid file
528  return netcdf_options;
529  }
530 
531  // Check for different NetCDF formats
532  int format_flags = 0;
533  for( int i = 0; i < ncFile.num_dims(); i++ )
534  {
535  const std::string dim_name = ncFile.get_dim( i )->name();
536 
537  if( dim_name == "grid_size" || dim_name == "grid_corners" || dim_name == "grid_rank" )
538  {
539  format_flags |= 1; // SCRIP format
540  }
541  else if( dim_name == "nodeCount" || dim_name == "elementCount" || dim_name == "maxNodePElement" )
542  {
543  format_flags |= 2; // ESMF format
544  }
545  else if( dim_name == "nCells" || dim_name == "nEdges" || dim_name == "nVertices" ||
546  dim_name == "vertexDegree" )
547  {
548  format_flags |= 4; // MPAS format
549  }
550  }
551 
552  // Apply format-specific options
553  if( format_flags & 2 )
554  { // ESMF format
555  netcdf_options += "PARALLEL_RESOLVE_SHARED_ENTS;VARIABLE=;";
556  }
557  else if( format_flags & 1 )
558  { // SCRIP format
559  netcdf_options += ""; // no extra options necessary for now
560  }
561  else if( format_flags & 4 )
562  { // MPAS format
563  netcdf_options += "PARALLEL_RESOLVE_SHARED_ENTS;NO_EDGES;NO_MIXED_ELEMENTS;VARIABLE=;";
564  }
565  }
566 
567  // Broadcast the options to all processes
568 #ifdef MOAB_HAVE_MPI
569  int line_size = netcdf_options.size();
570  MPI_Bcast( &line_size, 1, MPI_INT, 0, MPI_COMM_WORLD );
571  if( proc_id != 0 )
572  {
573  netcdf_options.resize( line_size );
574  }
575  MPI_Bcast( const_cast< char* >( netcdf_options.data() ), line_size, MPI_CHAR, 0, MPI_COMM_WORLD );
576 #endif
577 
578  return netcdf_options;
579  }
580 
581  // Default options for other file types
582  return "PARALLEL=BCAST_DELETE;PARTITION=TRIVIAL;PARALLEL_RESOLVE_SHARED_ENTS;";
583  }
584 
585  private:
586  /**
587  * @brief Print detailed help message with usage examples
588  * @param progName Program name
589  */
590  void printHelp( const char* progName ) const
591  {
592  if( this->proc_id != 0 ) return;
593 
594  std::cout << "MOAB-Tempest: A mesh generation and remapping tool\n"
595  << "==================================================\n\n"
596  << "Usage: " << progName << " [OPTIONS]\n\n"
597  << "Mesh Generation Options:\n"
598  << " -t, --type TYPE Type of mesh to generate (required for mesh generation):\n"
599  << " 0 = Cubed-Sphere (CS)\n"
600  << " 1 = Regular Latitude-Longitude (RLL)\n"
601  << " 2 = Icosahedral (ICO)\n"
602  << " 3 = TempestRemap overlap (thin interface))\n"
603  << " 4 = MOAB with TempestRemap overlap in memory\n"
604  << " 5 = Parallel handling of Overlap meshes with MOAB (recommended)\n\n"
605  << " -r, --res N Resolution (number of elements on edge, default: 10)\n"
606  << " -f, --file FILE Output filename (default: output.h5m)\n\n"
607  << "Discretization Options:\n"
608  << " -m, --method METHOD Discretization method (default: fv):\n"
609  << " fv = Finite Volume\n"
610  << " cgll = Continuous Galerkin with Legendre-Gauss-Lobatto\n"
611  << " dgll = Discontinuous Galerkin with Legendre-Gauss-Lobatto\n\n"
612  << " -o, --order N Discretization order (default: 1, range: 1-4)\n\n"
613  << "Remapping Options:\n"
614  << " --mono N Monotonicity constraints (default: 0):\n"
615  << " 0 = No monotonicity\n"
616  << " 1 = Basic monotonicity\n"
617  << " 2 = Full monotonicity with bounds\n"
618  << " 3 = Strict monotonicity\n\n"
619  << " --limiter TYPE Nonlinear limiting (optional):\n"
620  << " none = No limiting (default)\n"
621  << " global = Global CAAS limiting\n"
622  << " local = Localized CAAS limiting\n"
623  << " qlt = Quasi-Local Tree-based limiting\n\n"
624  << "Input/Output Options:\n"
625  << " -l, --load FILE Load input mesh file (use twice for source and target)\n"
626  << " -i, --global_id TAG Global ID tag name (default: GLOBAL_ID)\n"
627  << " --diagnostics Print diagnostic information\n\n"
628  << "Miscellaneous Options:\n"
629  << " --manual Show this help message and exit\n"
630  << " --version Show version information\n\n"
631  << "Examples:\n"
632  << " # Generate a cubed-sphere mesh with resolution 25\n"
633  << " " << progName << " --type 0 --res 25 -f cs_mesh.h5m\n\n"
634  << " # Generate a latitude-longitude mesh with resolution 180\n"
635  << " " << progName << " --type 1 --res 180 -f rll_mesh.h5m\n\n"
636  << " # Create a map between two meshes with order 4\n"
637  << " " << progName << " --type 5 --load source_mesh.h5m --load target_mesh.h5m \\\n"
638  << " --method cgll --order 4 --global_id GLOBAL_DOFS \\\n"
639  << " --method fv --order 1 --limiter 1 --file map.nc\n";
640  }
641 
642  /**
643  * @brief Get mesh type as string
644  * @return String representation of mesh type
645  */
646  std::string getMeshTypeName() const
647  {
648  switch( this->meshType )
649  {
650  case moab::TempestRemapper::CS:
651  return "Cubed-Sphere";
652  case moab::TempestRemapper::RLL:
653  return "Latitude-Longitude";
654  case moab::TempestRemapper::ICO:
655  return "Icosahedral";
656  case moab::TempestRemapper::OVERLAP_FILES:
657  return "Overlap (files)";
658  case moab::TempestRemapper::OVERLAP_MEMORY:
659  return "Overlap (memory)";
660  case moab::TempestRemapper::OVERLAP_MOAB:
661  return "Overlap (MOAB)";
662  default:
663  return "Unknown";
664  }
665  }
666 
667  /**
668  * @brief Process mesh options from command line
669  * @param opts Program options
670  * @param expectedFVMethod Expected finite volume method
671  * @param nlayer_input Number of ghost layers
672  */
673  void processMeshOptions( ProgOptions& opts )
674  {
675  if( this->meshType <= moab::TempestRemapper::ICO ) return;
676 
677  // Process input files
678  std::vector< std::string > inputFiles;
679  opts.getOptAllArgs( "load,l", inputFiles );
680  if( !inputFiles.empty() )
681  {
682  this->inFilenames = inputFiles;
683  if( this->inFilenames.size() != 2 )
684  {
685  throw std::runtime_error( "Exactly two input filenames must be provided with -l/--load" );
686  }
687  }
688 
689  // Process discretization orders
690  std::vector< int > orders;
691  opts.getOptAllArgs( "order,o", orders );
692  if( !orders.empty() )
693  {
694  this->disc_orders = orders;
695  if( this->disc_orders.size() == 1 )
696  {
697  this->disc_orders.push_back( this->disc_orders[0] );
698  }
699  else if( this->disc_orders.size() != 2 )
700  {
701  throw std::runtime_error( "Must specify 1 or 2 values for order (source [target])" );
702  }
703 
704  for( const auto& order : this->disc_orders )
705  {
706  if( order < 1 || order > 4 )
707  {
708  throw std::runtime_error( "Discretization order must be between 1 and 4" );
709  }
710  }
711  }
712 
713  // Process discretization methods
714  std::vector< std::string > methods;
715  opts.getOptAllArgs( "method,m", methods );
716  if( !methods.empty() )
717  {
718  this->disc_methods = methods;
719  if( this->disc_methods.size() == 1 )
720  {
721  // Use same method for both source and target
722  this->disc_methods.push_back( this->disc_methods[0] );
723  }
724  else if( this->disc_methods.size() != 2 )
725  {
726  throw std::runtime_error( "Must specify 1 or 2 values for method (source [target])" );
727  }
728 
729  // Validate method values
730  for( const auto& method : this->disc_methods )
731  {
732  if( method != "fv" && method != "cgll" && method != "dgll" && method != "pcloud" )
733  {
734  throw std::runtime_error( "Invalid method '" + method + "'. Must be one of: fv, cgll, dgll" );
735  }
736  }
737  }
738 
739  // Process DOF tag names
740  std::vector< std::string > tags;
741  opts.getOptAllArgs( "global_id,i", tags );
742  if( !tags.empty() )
743  {
744  this->doftag_names = tags;
745  if( this->doftag_names.size() == 1 )
746  {
747  // Use same tag name for both source and target
748  this->doftag_names.push_back( this->doftag_names[0] );
749  }
750  else if( this->doftag_names.size() != 2 )
751  {
752  throw std::runtime_error( "Must specify 1 or 2 values for DOF tag names (source [target])" );
753  }
754  }
755 
756  // Process output filename if specified
757  std::string outFile;
758  if( opts.getOpt( "file,f", &outFile ) )
759  {
760  this->outFilename = outFile;
761  }
762  // Note: configureMapOptions is now called from ParseCLOptions after processMeshOptions completes
763  }
764 
765  /**
766  * @brief Print all runtime parameters in a formatted way
767  */
768  void printRuntimeParameters() const
769  {
770  if( this->proc_id != 0 ) return;
771 
772  constexpr int width = 60;
773 
774  std::cout << std::string( width, '=' ) << "\n";
775  std::cout << " MOAB-TempestRemap Runtime Configuration " << "\n";
776  std::cout << std::string( width, '=' );
777 
778  // Input files
779  if( this->meshType == moab::TempestRemapper::OVERLAP_MOAB )
780  {
781  if( !this->inFilenames.empty() )
782  {
783  std::cout << "\n\nInput Files:";
784  std::cout << "\n Source mesh: " << this->inFilenames[0];
785  std::cout << "\n Target mesh: " << this->inFilenames[1];
786  }
787 
788  std::cout << "\n\nOutput Files:";
789  if( !skip_intersection )
790  std::cout << "\n Intersection mesh: "
791  << ( this->computeWeights ? this->intxFilename : this->outFilename );
792  if( computeWeights ) std::cout << "\n Remap weights: " << this->outFilename;
793  }
794 
795  // Mesh configuration
796  std::cout << "\n\nMesh Configuration:";
797  std::cout << "\n Mesh type: " << this->getMeshTypeName();
798  if( this->meshType <= moab::TempestRemapper::ICO )
799  std::cout << "\n Resolution: " << this->blockSize;
800  if( this->meshType == moab::TempestRemapper::ICO && computeDual )
801  std::cout << "\n Compute dual: " << ( this->computeDual ? "Yes" : "No" );
802 
803  if( computeWeights )
804  {
805  std::cout << "\n Gnomonic projection: " << ( this->useGnomonicProjection ? "Yes" : "No" );
806  std::cout << "\n Intersection algorithm: " << ( this->kdtreeSearch ? "KdTree search" : "Advancing front" );
807 
808  // Discretization settings
809  std::cout << "\n\nDiscretization:";
810  std::cout << "\n Source: " << this->disc_methods[0] << " (order " << this->disc_orders[0]
811  << ")";
812  std::cout << "\n Target: " << this->disc_methods[1] << " (order " << this->disc_orders[1]
813  << ")";
814 
815  // Remapping options
816  std::cout << "\n\nRemapping Options:";
817  std::cout << "\n Method: "
818  << ( this->mapOptions.strMethod.empty() ? "Default" : this->mapOptions.strMethod );
819  std::cout << "\n Monotonicity: " << ( this->ensureMonotonicity ? "Yes" : "No" );
820  std::cout << "\n Volumetric: " << ( this->fVolumetric ? "Yes" : "No" );
821  std::cout << "\n Check consistency: " << ( this->fCheck ? "Yes" : "No" );
822  std::cout << "\n Skip intersection: " << ( this->skip_intersection ? "Yes" : "No" );
823  }
824 
825  // Parallel configuration
826  std::cout << "\n\nParallel Configuration:";
827  std::cout << "\n MPI Processes: " << this->n_procs;
828  if( this->meshType > moab::TempestRemapper::ICO ) std::cout << "\n Number of Ghost Layers: " << this->nlayers;
829 
830  std::cout << "\n\n" << std::string( width, '=' ) << "\n\n";
831  }
832 
833  /**
834  * @brief Configure map options based on command line parameters
835  * @param nlayer_input Number of ghost layers
836  */
837  void configureMapOptions( int nlayer_input )
838  {
839  // Set polynomial orders with bounds checking
840  this->mapOptions.nPin = ( this->disc_orders.empty() ) ? 1 : this->disc_orders[0];
841  this->mapOptions.nPout = ( this->disc_orders.size() > 1 ) ? this->disc_orders[1] : this->mapOptions.nPin;
842 
843  // Initialize flags
844  this->mapOptions.fSourceConcave = false;
845  this->mapOptions.fTargetConcave = false;
846  this->mapOptions.strMethod.clear();
847 
848  // Configure finite volume method if specified
849  if( this->fvMethod != "none" )
850  {
851  this->mapOptions.strMethod = this->fvMethod + ";";
852  this->mapOptions.fNoConservation = true;
853  }
854 
855  // Configure monotonicity with validation
856  this->ensureMonotonicity = std::max( 0, std::min( 3, this->ensureMonotonicity ) ); // Clamp to 0-3
857  switch( this->ensureMonotonicity )
858  {
859  case 0:
860  this->mapOptions.fMonotone = false;
861  break;
862  case 3:
863  this->mapOptions.strMethod += "mono3;";
864  this->mapOptions.fMonotone = true;
865  break;
866  case 2:
867  this->mapOptions.strMethod += "mono2;";
868  this->mapOptions.fMonotone = true;
869  break;
870  case 1:
871  default:
872  this->mapOptions.fMonotone = true;
873  break;
874  }
875 
876  // Set other options
877  this->mapOptions.fNoCorrectAreas = false;
878  this->mapOptions.fNoCheck = !this->fCheck;
879 
880  // Add volumetric flag if needed
881  if( this->fVolumetric )
882  {
883  this->mapOptions.strMethod += "volumetric;";
884  }
885 
886  // Set number of ghost layers based on method and order.
887  // FV order 1 needs 0 ghost layers; FV order p > 1 needs p+1 ghost layers.
888  if( this->fvMethod == "delaunay" || this->fvMethod == "bilin" )
889  {
890  this->skip_intersection = true;
891  this->nlayers = 3; // conservative
892  }
893  else
894  {
895  // order 1: no ghost layers
896  // order p: p+1 layers (again, being conservative)
897  this->nlayers = ( this->mapOptions.nPin > 1 ) ? this->mapOptions.nPin + 1 : 0;
898  }
899 
900  // User-supplied value always overrides the internal default (even 0 is valid).
901  if( nlayer_input >= 0 )
902  {
903  this->nlayers = nlayer_input;
904  }
905 
906  // Configure output
907  this->mapOptions.strOutputMapFile = this->outFilename;
908  this->mapOptions.strOutputFormat = "Netcdf4";
909  }
910 };
911 
912 // Forward declare some methods
913 static moab::ErrorCode CreateTempestMesh( ToolContext&, moab::TempestRemapper& remapper, Mesh* );
914 static inline constexpr double sample_constant( double dLon, double dLat ) noexcept;
915 static inline double sample_slow_harmonic( double dLon, double dLat ) noexcept;
916 static inline double sample_fast_harmonic( double dLon, double dLat ) noexcept;
917 static inline double sample_stationary_vortex( double dLon, double dLat ) noexcept;
918 
919 /////////////////////////////////////////////////////////////
920 
921 //#define MOAB_DBG
922 int main( int argc, char* argv[] )
923 {
924  try
925  {
926  NcError error( NcError::verbose_nonfatal );
927  std::stringstream sstr;
928  std::string historyStr;
929 
930  int proc_id = 0, nprocs = 1;
931 #ifdef MOAB_HAVE_MPI
932  MPI_Init( &argc, &argv );
933  MPI_Comm_rank( MPI_COMM_WORLD, &proc_id );
934  MPI_Comm_size( MPI_COMM_WORLD, &nprocs );
935 #endif
936 
937  moab::Core* mbCore = new( std::nothrow ) moab::Core;
938 
939  if( nullptr == mbCore )
940  {
941  return 1;
942  }
943 
944  // Build the history string
945  for( int ia = 0; ia < argc; ++ia )
946  historyStr += std::string( argv[ia] ) + " ";
947 
948  ToolContext* runCtx;
949 #ifdef MOAB_HAVE_MPI
950  moab::ParallelComm* pcomm = new moab::ParallelComm( mbCore, MPI_COMM_WORLD, 0 );
951 
952  runCtx = new ToolContext( mbCore, pcomm );
953  const char* writeOptions = ( nprocs > 1 ? "PARALLEL=WRITE_PART" : "" );
954 #else
955  runCtx = new ToolContext( mbCore );
956  const char* writeOptions = "";
957 #endif
958  runCtx->ParseCLOptions( argc, argv );
959 
960  const double radius_src = 1.0 /*2.0*acos(-1.0)*/;
961  const double radius_dest = 1.0 /*2.0*acos(-1.0)*/;
962 
963  moab::DebugOutput& outputFormatter = runCtx->outputFormatter;
964 
965 #ifdef MOAB_HAVE_MPI
966  moab::TempestRemapper remapper( mbCore, pcomm );
967 #else
968  moab::TempestRemapper remapper( mbCore );
969 #endif
970  remapper.meshValidate = true;
971  remapper.constructEdgeMap = true;
972  remapper.initialize();
973 
974  // Default area_method = lHuiller; Options: Girard, lHuiller, GaussQuadrature (if TR is available)
975  moab::IntxAreaUtils areaAdaptor( moab::IntxAreaUtils::lHuiller );
976 
977  Mesh* tempest_mesh = new Mesh();
978  MB_CHK_SET_ERR( CreateTempestMesh( *runCtx, remapper, tempest_mesh ), "Failed to create tempest mesh" );
979 
980  if( runCtx->meshType == moab::TempestRemapper::OVERLAP_MEMORY )
981  {
982  // Compute intersections with MOAB
983  // For the overlap method, choose between: "fuzzy", "exact" or "mixed"
984  assert( runCtx->meshes.size() == 3 );
985 
986 #ifdef MOAB_HAVE_MPI
987  MB_CHK_SET_ERR( pcomm->check_all_shared_handles(), "Failed to check all shared handles" );
988 #endif
989 
990  // Load the meshes and validate
991  MB_CHK_SET_ERR( remapper.ConvertTempestMesh( moab::Remapper::SourceMesh ), "Failed to convert source mesh" );
992  MB_CHK_SET_ERR( remapper.ConvertTempestMesh( moab::Remapper::TargetMesh ), "Failed to convert target mesh" );
993  MB_CHK_SET_ERR( remapper.ConvertTempestMesh( moab::Remapper::OverlapMesh ), "Failed to convert overlap mesh" );
994  if( !runCtx->skip_io )
995  {
996  MB_CHK_SET_ERR( mbCore->write_mesh( "tempest_intersection.h5m", &runCtx->meshsets[2], 1 ),
997  "Failed to write TempestRemap intersection mesh in MOAB format" );
998  }
999 
1000  // print verbosely about the problem setting
1001  size_t velist[6], gvelist[6];
1002  {
1003  moab::Range rintxverts, rintxelems;
1004  MB_CHK_SET_ERR( mbCore->get_entities_by_dimension( runCtx->meshsets[0], 0, rintxverts ),
1005  "Failed to get vertices" );
1006  MB_CHK_SET_ERR( mbCore->get_entities_by_dimension( runCtx->meshsets[0], 2, rintxelems ),
1007  "Failed to get elements" );
1008  velist[0] = rintxverts.size();
1009  velist[1] = rintxelems.size();
1010 
1011  moab::Range bintxverts, bintxelems;
1012  MB_CHK_SET_ERR( mbCore->get_entities_by_dimension( runCtx->meshsets[1], 0, bintxverts ),
1013  "Failed to get vertices" );
1014  MB_CHK_SET_ERR( mbCore->get_entities_by_dimension( runCtx->meshsets[1], 2, bintxelems ),
1015  "Failed to get elements" );
1016  velist[2] = bintxverts.size();
1017  velist[3] = bintxelems.size();
1018  }
1019 
1020  moab::EntityHandle intxset; // == remapper.GetMeshSet(moab::Remapper::OverlapMesh);
1021 
1022  // Compute intersections with MOAB
1023  {
1024  // Create the intersection on the sphere object
1025  runCtx->timer_push( "setup the intersector" );
1026 
1027  moab::Intx2MeshOnSphere* mbintx = new moab::Intx2MeshOnSphere( mbCore );
1028  mbintx->set_error_tolerance( runCtx->epsrel );
1029  mbintx->set_box_error( runCtx->boxeps );
1030  mbintx->set_radius_source_mesh( radius_src );
1031  mbintx->set_radius_destination_mesh( radius_dest );
1032 #ifdef MOAB_HAVE_MPI
1033  mbintx->set_parallel_comm( pcomm );
1034 #endif
1035  MB_CHK_SET_ERR( mbintx->FindMaxEdges( runCtx->meshsets[0], runCtx->meshsets[1] ),
1036  "Failed to find max edges" );
1037 
1038 #ifdef MOAB_HAVE_MPI
1039  moab::Range local_verts;
1040  MB_CHK_SET_ERR( mbintx->build_processor_euler_boxes( runCtx->meshsets[1], local_verts ),
1041  "Failed to build processor euler boxes" );
1042 
1043  runCtx->timer_pop();
1044 
1045  moab::EntityHandle covering_set;
1046  runCtx->timer_push( "communicate the mesh" );
1047  // we compute just intersection here, no need for extra ghost layers anyway
1048  // ghost layers are needed in coverage for bilinear map, which does not actually need intersection
1049  // this will be fixed in the future, bilinear map needs just coverage, not intersection
1050  // so I am not passing the ghost layer here, even though there is an option in runCtx for a ghost layer
1051  // NOTE: This is a communication-heavy kernel if mesh is distributed very differently
1052  MB_CHK_SET_ERR( mbintx->construct_covering_set( runCtx->meshsets[0], covering_set ),
1053  "Failed to construct covering set" );
1054  runCtx->timer_pop();
1055 
1056  // print verbosely about the problem setting
1057  {
1058  moab::Range cintxverts, cintxelems;
1059  MB_CHK_SET_ERR( mbCore->get_entities_by_dimension( covering_set, 0, cintxverts ),
1060  "Failed to get vertices" );
1061  MB_CHK_SET_ERR( mbCore->get_entities_by_dimension( covering_set, 2, cintxelems ),
1062  "Failed to get elements" );
1063  velist[4] = cintxverts.size();
1064  velist[5] = cintxelems.size();
1065  }
1066 
1067  MPI_Reduce( velist, gvelist, 6, MPI_UINT64_T, MPI_SUM, 0, MPI_COMM_WORLD );
1068 
1069 #else
1070  moab::EntityHandle covering_set = runCtx->meshsets[0];
1071  for( int i = 0; i < 6; i++ )
1072  gvelist[i] = velist[i];
1073 #endif
1074 
1075  if( !proc_id )
1076  {
1077  outputFormatter.printf( 0, "The source set contains %lu vertices and %lu elements \n", gvelist[0],
1078  gvelist[0] );
1079  outputFormatter.printf( 0, "The covering set contains %lu vertices and %lu elements \n", gvelist[2],
1080  gvelist[2] );
1081  outputFormatter.printf( 0, "The target set contains %lu vertices and %lu elements \n", gvelist[1],
1082  gvelist[1] );
1083  }
1084 
1085  // Now let's invoke the MOAB intersection algorithm in parallel with a
1086  // source and target mesh set representing two different decompositions
1087  runCtx->timer_push( "compute intersections with MOAB" );
1088  MB_CHK_SET_ERR( mbCore->create_meshset( moab::MESHSET_SET, intxset ), "Can't create new set" );
1089  MB_CHK_SET_ERR( mbintx->intersect_meshes( covering_set, runCtx->meshsets[1], intxset ),
1090  "Can't compute the intersection of meshes on the sphere" );
1091  runCtx->timer_pop();
1092 
1093  // free the memory
1094  delete mbintx;
1095  }
1096 
1097  {
1098  moab::Range intxelems, intxverts;
1099  MB_CHK_SET_ERR( mbCore->get_entities_by_dimension( intxset, 2, intxelems ), "Failed to get elements" );
1100  MB_CHK_SET_ERR( mbCore->get_entities_by_dimension( intxset, 0, intxverts, true ),
1101  "Failed to get vertices" );
1102  outputFormatter.printf( 0, "The intersection set contains %lu elements and %lu vertices \n",
1103  intxelems.size(), intxverts.size() );
1104 
1105  double initial_sarea =
1106  areaAdaptor.area_on_sphere( mbCore, runCtx->meshsets[0],
1107  radius_src ); // use the target to compute the initial area
1108  double initial_tarea =
1109  areaAdaptor.area_on_sphere( mbCore, runCtx->meshsets[1],
1110  radius_dest ); // use the target to compute the initial area
1111  double intx_area = areaAdaptor.area_on_sphere( mbCore, intxset, radius_src );
1112 
1113  outputFormatter.printf( 0, "mesh areas: source = %12.10f, target = %12.10f, intersection = %12.10f \n",
1114  initial_sarea, initial_tarea, intx_area );
1115  outputFormatter.printf( 0, "relative error w.r.t source = %12.10e, target = %12.10e \n",
1116  fabs( intx_area - initial_sarea ) / initial_sarea,
1117  fabs( intx_area - initial_tarea ) / initial_tarea );
1118  }
1119 
1120  // Write out our computed intersection file
1121  if( !runCtx->skip_io )
1122  {
1123  MB_CHK_SET_ERR( mbCore->write_mesh( "moab_intersection.h5m", &intxset, 1 ),
1124  "Failed to write the intersection" );
1125  }
1126 
1127  if( runCtx->computeWeights )
1128  {
1129  runCtx->timer_push( "compute weights with the Tempest meshes" );
1130  // Call to generate an offline map with the tempest meshes
1131  OfflineMap weightMap;
1132  if( GenerateOfflineMapWithMeshes( *runCtx->meshes[0], *runCtx->meshes[1], *runCtx->meshes[2],
1133  runCtx->disc_methods[0], // std::string strInputType
1134  runCtx->disc_methods[1], // std::string strOutputType,
1135  runCtx->mapOptions, weightMap ) != 0 )
1136  throw std::runtime_error( "Could not generate offline map with TempestRemap" );
1137  runCtx->timer_pop();
1138 
1139  std::map< std::string, std::string > mapAttributes;
1140  if( !runCtx->skip_io ) weightMap.Write( "outWeights.nc", mapAttributes );
1141  }
1142  }
1143  else if( runCtx->meshType == moab::TempestRemapper::OVERLAP_MOAB )
1144  {
1145  // Usage: mpiexec -n 2 tools/mbtempest -t 5 -l mycs_2.h5m -l myico_2.h5m -f myoverlap_2.h5m
1146 #ifdef MOAB_HAVE_MPI
1147  MB_CHK_SET_ERR( pcomm->check_all_shared_handles(), "Checking shared handles failed." );
1148 #endif
1149 
1150  // print verbosely about the problem setting
1151  size_t velist[4] = { 0, 0, 0, 0 }, gvelist[4] = { 0, 0, 0, 0 };
1152  {
1153  moab::Range srcverts, srcelems;
1154  MB_CHK_SET_ERR( mbCore->get_entities_by_dimension( runCtx->meshsets[0], 0, srcverts ),
1155  "Failed to get vertices" );
1156  MB_CHK_SET_ERR( mbCore->get_entities_by_dimension( runCtx->meshsets[0], 2, srcelems ),
1157  "Failed to get elements" );
1158  MB_CHK_SET_ERR( moab::IntxUtils::fix_degenerate_quads( mbCore, runCtx->meshsets[0] ),
1159  "Failed to fix degenerate quads" );
1160  if( runCtx->enforceConvexity )
1161  {
1162  MB_CHK_SET_ERR( moab::IntxUtils::enforce_convexity( mbCore, runCtx->meshsets[0], proc_id ),
1163  "Failed to enforce convexity" );
1164  }
1165  MB_CHK_SET_ERR( areaAdaptor.positive_orientation( mbCore, runCtx->meshsets[0], radius_src ),
1166  "Failed to enforce positive orientation" );
1167  velist[0] = srcverts.size();
1168  velist[1] = srcelems.size();
1169 
1170  moab::Range tgtverts, tgtelems;
1171  MB_CHK_SET_ERR( mbCore->get_entities_by_dimension( runCtx->meshsets[1], 0, tgtverts ),
1172  "Failed to get vertices" );
1173  MB_CHK_SET_ERR( mbCore->get_entities_by_dimension( runCtx->meshsets[1], 2, tgtelems ),
1174  "Failed to get elements" );
1175  MB_CHK_SET_ERR( moab::IntxUtils::fix_degenerate_quads( mbCore, runCtx->meshsets[1] ),
1176  "Failed to fix degenerate quads" );
1177  if( runCtx->enforceConvexity )
1178  {
1179  MB_CHK_SET_ERR( moab::IntxUtils::enforce_convexity( mbCore, runCtx->meshsets[1], proc_id ),
1180  "Failed to enforce convexity" );
1181  }
1182  MB_CHK_SET_ERR( areaAdaptor.positive_orientation( mbCore, runCtx->meshsets[1], radius_dest ),
1183  "Failed to enforce positive orientation" );
1184  velist[2] = tgtverts.size();
1185  velist[3] = tgtelems.size();
1186  }
1187  //MB_CHK_SET_ERR( mbCore->write_file( "source_mesh.h5m", nullptr, writeOptions, &runCtx->meshsets[0], 1 ), "Could not write source mesh" );
1188  //MB_CHK_SET_ERR( mbCore->write_file( "target_mesh.h5m", nullptr, writeOptions, &runCtx->meshsets[1], 1 ), "Could not write target mesh" );
1189 
1190  // if( runCtx->nlayers && nprocs > 1 )
1191  // {
1192  // remapper.ResetMeshSet( moab::Remapper::SourceMesh, runCtx->meshsets[3] );
1193  // runCtx->meshes[0] = remapper.GetMesh( moab::Remapper::SourceMesh ); // ?
1194  // }
1195 
1196  // First compute the covering set such that the target elements are fully covered by the
1197  // local source grid
1198  runCtx->timer_push( "construct covering set for intersection" );
1199  // if ghosting, do not use gnomonic projection
1200  if( runCtx->nlayers > 0 ) runCtx->useGnomonicProjection = false;
1201  MB_CHK_SET_ERR( remapper.ConstructCoveringSet( runCtx->epsrel, 1.0, 1.0, runCtx->boxeps, runCtx->rrmGrids,
1202  runCtx->useGnomonicProjection, runCtx->nlayers ),
1203  "Failed to construct covering set" );
1204  runCtx->timer_pop();
1205 
1206 #ifdef MOAB_HAVE_MPI
1207  MPI_Reduce( velist, gvelist, 4, MPI_UINT64_T, MPI_SUM, 0, MPI_COMM_WORLD );
1208 #else
1209  for( int i = 0; i < 4; i++ )
1210  gvelist[i] = velist[i];
1211 #endif
1212  if( !proc_id && runCtx->print_diagnostics )
1213  {
1214  outputFormatter.printf( 0, "The source set contains %lu vertices and %lu elements \n", gvelist[0],
1215  gvelist[1] );
1216  outputFormatter.printf( 0, "The target set contains %lu vertices and %lu elements \n", gvelist[2],
1217  gvelist[3] );
1218  }
1219 
1220  if( runCtx->skip_intersection )
1221  {
1222  if( !proc_id ) outputFormatter.printf( 0, "Skipping mesh intersection computation.\n" );
1223  }
1224  else
1225  {
1226  // Compute intersections with MOAB with either the Kd-tree or the advancing front algorithm
1227  runCtx->timer_push( "setup and compute mesh intersections" );
1228  MB_CHK_SET_ERR( remapper.ComputeOverlapMesh( runCtx->kdtreeSearch, false ),
1229  "Failed to compute mesh intersections" );
1230  runCtx->timer_pop();
1231  }
1232 
1233  // print some diagnostic checks to see if the overlap grid resolved the input meshes
1234  // correctly
1235  // Compute ghost overlap elements once; reused for both area diagnostics and intx file write
1236  moab::Range ghostOverlapElems;
1237 #ifdef MOAB_HAVE_MPI
1238  if( nprocs > 1 && !runCtx->skip_intersection )
1239  MB_CHK_SET_ERR( remapper.GetOverlapAugmentedEntities( ghostOverlapElems ),
1240  "Failed to get ghost overlap entities" );
1241 #endif
1242 
1243  double dTotalOverlapArea = 0.0;
1244  if( runCtx->print_diagnostics && !runCtx->skip_intersection )
1245  {
1246  // Areas for source, target, overlap meshes
1247  double local_areas[3] = { 0, 0, 0 },
1248  global_areas[3] = { 0, 0, 0 };
1249 
1250  // Helper: compute area of a meshset excluding cells with GRID_IMASK==0.
1251  // Both source and target SCRIP grids may have a land/sea mask; the intersection
1252  // only covers unmasked cells, so comparing full-mesh areas gives a misleading error.
1253  auto area_unmasked = [&]( moab::EntityHandle meshset, double radius ) -> double {
1254  moab::Tag imaskTag = 0;
1255  mbCore->tag_get_handle( "GRID_IMASK", imaskTag );
1256  if( !imaskTag ) return areaAdaptor.area_on_sphere( mbCore, meshset, radius );
1257  moab::Range cells;
1258  mbCore->get_entities_by_dimension( meshset, 2, cells );
1259  std::vector< int > masks( cells.size(), 1 );
1260  mbCore->tag_get_data( imaskTag, cells, masks.data() );
1261  moab::Range maskedCells;
1262  size_t idx = 0;
1263  for( auto it = cells.begin(); it != cells.end(); ++it, ++idx )
1264  if( !masks[idx] ) maskedCells.insert( *it );
1265  moab::Range unmasked = moab::subtract( cells, maskedCells );
1266  moab::EntityHandle tmpSet;
1267  mbCore->create_meshset( moab::MESHSET_SET, tmpSet );
1268  mbCore->add_entities( tmpSet, unmasked );
1269  double area = areaAdaptor.area_on_sphere( mbCore, tmpSet, radius );
1270  mbCore->delete_entities( &tmpSet, 1 );
1271  return area;
1272  };
1273 
1274  local_areas[0] = area_unmasked( runCtx->meshsets[0], radius_src );
1275  local_areas[1] = area_unmasked( runCtx->meshsets[1], radius_dest );
1276  // Exclude ghost overlap elements from area sum to avoid double-counting after MPI_Allreduce
1277  {
1278  moab::Range ownedOverlapElems;
1279  MB_CHK_SET_ERR( mbCore->get_entities_by_dimension( runCtx->meshsets[2], 2, ownedOverlapElems ),
1280  "Failed to get overlap elements" );
1281  ownedOverlapElems = moab::subtract( ownedOverlapElems, ghostOverlapElems );
1282  moab::EntityHandle ownedOverlapSet;
1283  MB_CHK_SET_ERR( mbCore->create_meshset( moab::MESHSET_SET, ownedOverlapSet ),
1284  "Can't create owned overlap meshset" );
1285  MB_CHK_SET_ERR( mbCore->add_entities( ownedOverlapSet, ownedOverlapElems ),
1286  "Can't add owned overlap elements" );
1287  local_areas[2] = areaAdaptor.area_on_sphere( mbCore, ownedOverlapSet, radius_src );
1288  MB_CHK_SET_ERR( mbCore->delete_entities( &ownedOverlapSet, 1 ), "Can't delete temp meshset" );
1289  }
1290 
1291 #ifdef MOAB_HAVE_MPI
1292  MPI_Allreduce( &local_areas[0], &global_areas[0], 3, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD );
1293 #else
1294  global_areas[0] = local_areas[0];
1295  global_areas[1] = local_areas[1];
1296  global_areas[2] = local_areas[2];
1297 #endif
1298  if( !proc_id )
1299  {
1300  outputFormatter.printf( 0,
1301  "initial area: source mesh = %12.14f, target mesh = "
1302  "%12.14f, overlap mesh = %12.14f\n",
1303  global_areas[0], global_areas[1], global_areas[2] );
1304  outputFormatter.printf( 0, "relative error w.r.t source = %12.14e, and target = %12.14e\n",
1305  fabs( global_areas[0] - global_areas[2] ) / global_areas[0],
1306  fabs( global_areas[1] - global_areas[2] ) / global_areas[1] );
1307  }
1308  dTotalOverlapArea = global_areas[2];
1309  }
1310 
1311  if( runCtx->intxFilename.size() && !runCtx->skip_intersection )
1312  {
1313  moab::EntityHandle writableOverlapSet;
1314  MB_CHK_SET_ERR( mbCore->create_meshset( moab::MESHSET_SET, writableOverlapSet ), "Can't create new set" );
1315  moab::EntityHandle meshOverlapSet = remapper.GetMeshSet( moab::Remapper::OverlapMesh );
1316  moab::Range ovEnts;
1317  MB_CHK_SET_ERR( mbCore->get_entities_by_dimension( meshOverlapSet, 2, ovEnts ), "Can't create new set" );
1318  MB_CHK_SET_ERR( mbCore->get_entities_by_dimension( meshOverlapSet, 0, ovEnts ), "Can't create new set" );
1319 
1320 #ifdef MOAB_HAVE_MPI
1321  // Exclude ghost overlap elements from the write: each ghost element is owned by another
1322  // rank and will be written from there. Including ghosts here causes duplicate entity
1323  // handles in the parallel HDF5 output and deadlocks the collective write.
1324  if( nprocs > 1 )
1325  {
1326  ovEnts = moab::subtract( ovEnts, ghostOverlapElems );
1327 #ifdef MOAB_DBG
1328  if( !runCtx->skip_io )
1329  {
1330  std::stringstream filename;
1331  filename << "aug_overlap" << runCtx->pcomm->rank() << ".h5m";
1332  MB_CHK_SET_ERR( mbCore->write_file( filename.str().c_str(), 0, 0, &meshOverlapSet, 1 ),
1333  "Failed to write the overlap set" );
1334  }
1335 #endif
1336  }
1337 #endif
1338  MB_CHK_SET_ERR( mbCore->add_entities( writableOverlapSet, ovEnts ), "adding local intx cells failed" );
1339 
1340 #ifdef MOAB_HAVE_MPI
1341 #ifdef MOAB_DBG
1342  if( nprocs > 1 && !runCtx->skip_io )
1343  {
1344  std::stringstream filename;
1345  filename << "writable_intx_" << runCtx->pcomm->rank() << ".h5m";
1346  MB_CHK_SET_ERR( mbCore->write_file( filename.str().c_str(), 0, 0, &writableOverlapSet, 1 ),
1347  "Failed to write the writable overlap set" );
1348  }
1349 #endif
1350 #endif
1351 
1352  size_t lastindex = runCtx->intxFilename.find_last_of( "." );
1353  sstr.str( "" );
1354  sstr << runCtx->intxFilename.substr( 0, lastindex ) << ".h5m";
1355  if( !runCtx->proc_id )
1356  std::cout << "Writing out the MOAB intersection mesh file to " << sstr.str() << std::endl;
1357 
1358  // Write out our computed intersection file
1359  if( !runCtx->skip_io )
1360  {
1361  MB_CHK_SET_ERR( mbCore->write_file( sstr.str().c_str(), nullptr, writeOptions, &writableOverlapSet, 1 ),
1362  "Failed to write the writable overlap set" );
1363  }
1364  }
1365 
1366  if( runCtx->computeWeights )
1367  {
1368  runCtx->meshes[2] = remapper.GetMesh( moab::Remapper::OverlapMesh );
1369  if( !runCtx->proc_id ) std::cout << std::endl;
1370 
1371  runCtx->timer_push( "setup computation of weights" );
1372  // Call to generate the remapping weights with the tempest meshes
1373  moab::TempestOnlineMap* weightMap = new moab::TempestOnlineMap( &remapper );
1374  runCtx->timer_pop();
1375 
1376  runCtx->timer_push( "compute weights with TempestRemap" );
1377  MB_CHK_SET_ERR( weightMap->GenerateRemappingWeights(
1378  runCtx->disc_methods[0], // std::string strInputType
1379  runCtx->disc_methods[1], // std::string strOutputType,
1380  runCtx->mapOptions, // const GenerateOfflineMapAlgorithmOptions& options
1381  runCtx->doftag_names[0], // const std::string& source_tag_name
1382  runCtx->doftag_names[1] // const std::string& target_tag_name
1383  ),
1384  "Failed to generate remapping weights" );
1385  runCtx->timer_pop();
1386 
1387  weightMap->PrintMapStatistics();
1388 
1389  // Invoke the CheckMap routine on the TempestRemap serial interface directly, if running
1390  // on a single process
1391  if( runCtx->fCheck )
1392  {
1393  const double dNormalTolerance = 1.0E-8;
1394  const double dStrictTolerance = 1.0E-12;
1395  weightMap->CheckMap( runCtx->fCheck, runCtx->fCheck, runCtx->fCheck && ( runCtx->ensureMonotonicity ),
1396  dNormalTolerance, dStrictTolerance, dTotalOverlapArea );
1397  }
1398 
1399  if( runCtx->outFilename.size() && !runCtx->skip_io )
1400  {
1401  std::map< std::string, std::string > attrMap;
1402  attrMap["MOABversion"] = std::string( MOAB_PACKAGE_VERSION );
1403  attrMap["Title"] = "MOAB-TempestRemap (mbtempest) Offline Regridding Weight Generator";
1404  attrMap["normalization"] = "ovarea";
1405  attrMap["remap_options"] = runCtx->mapOptions.strMethod;
1406  attrMap["domain_a"] = runCtx->inFilenames[0];
1407  attrMap["domain_b"] = runCtx->inFilenames[1];
1408  if( runCtx->intxFilename.size() ) attrMap["domain_aUb"] = runCtx->intxFilename;
1409  attrMap["map_aPb"] = runCtx->outFilename;
1410  attrMap["methodorder_a"] = runCtx->disc_methods[0] + ":" + std::to_string( runCtx->disc_orders[0] ) +
1411  ":" + std::string( runCtx->doftag_names[0] );
1412  attrMap["concave_a"] = runCtx->mapOptions.fSourceConcave ? "true" : "false";
1413  attrMap["methodorder_b"] = runCtx->disc_methods[1] + ":" + std::to_string( runCtx->disc_orders[1] ) +
1414  ":" + std::string( runCtx->doftag_names[1] );
1415  attrMap["concave_b"] = runCtx->mapOptions.fTargetConcave ? "true" : "false";
1416  attrMap["bubble"] = runCtx->mapOptions.fNoBubble ? "false" : "true";
1417  attrMap["history"] = historyStr;
1418 
1419  // Write the map file to disk in parallel using either HDF5 or SCRIP interface
1420  // in extra case; maybe need a better solution, just create it with the right meshset
1421  // from the beginning;
1422  MB_CHK_SET_ERR( weightMap->WriteParallelMap( runCtx->outFilename.c_str(), attrMap ),
1423  "Failed writing the parallel map to disk" );
1424  }
1425 
1426  if( runCtx->verifyWeights )
1427  {
1428  // Let us pick a sampling test function for solution evaluation
1429  // SH, SV, FH, C, USERVAR
1430  bool userVariable = false;
1431  moab::TempestOnlineMap::sample_function testFunction;
1432  if( !runCtx->variableToVerify.compare( "SH" ) )
1433  testFunction = &sample_slow_harmonic;
1434  else if( !runCtx->variableToVerify.compare( "FH" ) )
1435  testFunction = &sample_fast_harmonic;
1436  else if( !runCtx->variableToVerify.compare( "SV" ) )
1437  testFunction = &sample_stationary_vortex;
1438  else if( !runCtx->variableToVerify.compare( "C" ) )
1439  testFunction = &sample_constant;
1440  else
1441  {
1442  userVariable = runCtx->variableToVerify.size() ? true : false;
1443  testFunction = runCtx->variableToVerify.size() ? nullptr : sample_stationary_vortex;
1444  }
1445 
1446  moab::Tag srcAnalyticalFunction;
1447  moab::Tag tgtAnalyticalFunction;
1448  moab::Tag tgtProjectedFunction;
1449  if( testFunction )
1450  {
1451  runCtx->timer_push( "describe a solution on source grid" );
1452  // MB_CHK_SET_ERR( mbCore->tag_get_handle( runCtx->variableToVerify.c_str(), srcAnalyticalFunction ),
1453  // "Failed to get analytical solution on source grid" );
1454  MB_CHK_SET_ERR( weightMap->DefineAnalyticalSolution( srcAnalyticalFunction,
1455  "AnalyticalSolnSrcExact",
1456  moab::Remapper::SourceMesh, testFunction ),
1457  "Failed to define analytical solution on source grid" );
1458  runCtx->timer_pop();
1459 
1460  // runCtx->timer_push( "exchange solution on source grid" );
1461  // moab::Range& srccovEnts = remapper.GetMeshEntities( moab::Remapper::CoveringMesh );
1462  // MB_CHK_SET_ERR( pcomm->exchange_tags( srcAnalyticalFunction, srccovEnts ),
1463  // "Failed to exchange analytical solution on source grid" );
1464  // runCtx->timer_pop();
1465 
1466  runCtx->timer_push( "describe a solution on target grid" );
1467  MB_CHK_SET_ERR( weightMap->DefineAnalyticalSolution(
1468  tgtAnalyticalFunction, "AnalyticalSolnTgtExact", moab::Remapper::TargetMesh,
1469  testFunction, &tgtProjectedFunction, "ProjectedSolnTgt" ),
1470  "Failed to define analytical solution on target grid" );
1471  runCtx->timer_pop();
1472  }
1473  else
1474  {
1475  MB_CHK_SET_ERR( mbCore->tag_get_handle( runCtx->variableToVerify.c_str(), srcAnalyticalFunction ),
1476  "Failed to get analytical solution on source grid" );
1477  MB_CHK_SET_ERR( mbCore->tag_get_handle( "ProjectedSolnTgt", 1, moab::MB_TYPE_DOUBLE,
1478  tgtProjectedFunction,
1479  moab::MB_TAG_DENSE | moab::MB_TAG_CREAT ),
1480  "Failed to get projected solution on target grid" );
1481  }
1482 
1483  // if( !runCtx->skip_io )
1484  {
1485  MB_CHK_SET_ERR( mbCore->write_file( "srcWithSolnTag.h5m", nullptr, writeOptions,
1486  &runCtx->meshsets[0], 1 ),
1487  "Failed to write the source mesh with solution tag" );
1488  }
1489 
1490  runCtx->timer_push( "compute solution projection on target grid" );
1491  MB_CHK_SET_ERR( weightMap->ApplyWeights( srcAnalyticalFunction, tgtProjectedFunction, false,
1492  runCtx->cassType ),
1493  "Failed to apply weights" );
1494  runCtx->timer_pop();
1495 
1496  // if( !runCtx->skip_io )
1497  {
1498  MB_CHK_SET_ERR( mbCore->write_file( "tgtWithSolnTag2.h5m", nullptr, writeOptions,
1499  &runCtx->meshsets[1], 1 ),
1500  "Failed to write the target mesh with projected solution tag" );
1501  }
1502 
1503  if( nprocs == 1 && runCtx->baselineFile.size() )
1504  {
1505  // save the field from tgtWithSolnTag2 in a text file, and global ids for cells
1506  moab::Range tgtEntities;
1507  if( runCtx->disc_methods[1] == "pcloud" )
1508  {
1509  MB_CHK_SET_ERR( mbCore->get_entities_by_dimension( runCtx->meshsets[1], 0, tgtEntities ),
1510  "Failed to get entities by dimension" );
1511  }
1512  else
1513  {
1514  MB_CHK_SET_ERR( mbCore->get_entities_by_dimension( runCtx->meshsets[1], 2, tgtEntities ),
1515  "Failed to get entities by dimension" );
1516  }
1517  std::vector< int > globIds( tgtEntities.size() );
1518  std::vector< double > vals( tgtEntities.size() );
1519  moab::Tag projTag;
1520  MB_CHK_SET_ERR( mbCore->tag_get_handle( "ProjectedSolnTgt", projTag ),
1521  "Failed to get projected solution tag" );
1522  moab::Tag gid = mbCore->globalId_tag();
1523  MB_CHK_SET_ERR( mbCore->tag_get_data( gid, tgtEntities, &globIds[0] ), "Failed to get global ids" );
1524  MB_CHK_SET_ERR( mbCore->tag_get_data( projTag, tgtEntities, &vals[0] ),
1525  "Failed to get projected solution" );
1526  std::fstream fs;
1527  fs.open( runCtx->baselineFile.c_str(), std::fstream::out );
1528  fs << std::setprecision( 15 ); // maximum precision for doubles
1529  for( size_t i = 0; i < tgtEntities.size(); i++ )
1530  fs << globIds[i] << " " << vals[i] << "\n";
1531  fs.close();
1532  // for good measure, save the source file too, with the tag AnalyticalSolnSrcExact
1533  // it will be used later to test, along with a target file
1534  if( !runCtx->skip_io )
1535  {
1536  MB_CHK_SET_ERR( mbCore->write_file( "srcWithSolnTag.h5m", nullptr, writeOptions,
1537  &runCtx->meshsets[0], 1 ),
1538  "Failed to write the source mesh with solution tag" );
1539  }
1540  }
1541 
1542  // compute error metrics if it is a known analytical functional
1543  if( !userVariable )
1544  {
1545  runCtx->timer_push( "compute error metrics against analytical solution on target grid" );
1546  std::map< std::string, double > errMetrics;
1547  MB_CHK_SET_ERR( weightMap->ComputeMetrics( moab::Remapper::TargetMesh, tgtAnalyticalFunction,
1548  tgtProjectedFunction, errMetrics, true ),
1549  "Failed to compute error metrics" );
1550  runCtx->timer_pop();
1551  }
1552  }
1553 
1554  delete weightMap;
1555  }
1556  }
1557 
1558  // Clean up
1559  remapper.clear();
1560  delete runCtx;
1561  delete mbCore;
1562 
1563 #ifdef MOAB_HAVE_MPI
1564  MPI_Finalize();
1565 #endif
1566  return 0;
1567  }
1568  catch( const std::exception& e )
1569  {
1570  std::cerr << "[mbtempest] Fatal error: " << e.what() << std::endl;
1571 #ifdef MOAB_HAVE_MPI
1572  MPI_Abort( MPI_COMM_WORLD, 1 );
1573 #endif
1574  return 1;
1575  }
1576  catch( ... )
1577  {
1578  std::cerr << "[mbtempest] Fatal: unknown exception caught" << std::endl;
1579 #ifdef MOAB_HAVE_MPI
1580  MPI_Abort( MPI_COMM_WORLD, 1 );
1581 #endif
1582  return 1;
1583  }
1584 }
1585 
1586 ///////////////////////////////////////////////////////////////////////////////
1587 
1588 // Helper functions for each mesh type
1589 namespace
1590 {
1591 
1592 #define TR_CHK_SET_ERR( err, msg ) \
1593  if( err ) \
1594  { \
1595  std::cout << "MOAB-TempestRemap Failure. ErrorCode (" << ( err ) << ") "; \
1596  MB_CHK_SET_ERR( moab::MB_FAILURE, msg ); \
1597  }
1598 
1599 moab::ErrorCode handleOverlapMemory( ToolContext& ctx, moab::TempestRemapper& remapper, Mesh* tempest_mesh )
1600 {
1601  using namespace moab;
1602 
1603  // resize the meshsets and meshes vectors
1604  ctx.meshsets.resize( 3 );
1605  ctx.meshes.resize( 3 );
1606 
1607  ctx.meshsets[0] = remapper.GetMeshSet( Remapper::SourceMesh );
1608  ctx.meshsets[1] = remapper.GetMeshSet( Remapper::TargetMesh );
1609  ctx.meshsets[2] = remapper.GetMeshSet( Remapper::OverlapMesh );
1610 
1611  // Load and process source mesh
1612  MB_CHK_SET_ERR( remapper.LoadMesh( Remapper::SourceMesh, ctx.inFilenames[0], TempestRemapper::DEFAULT ),
1613  "Failed to load MOAB Source mesh" );
1614 
1615  // Load and process target mesh
1616  MB_CHK_SET_ERR( remapper.LoadMesh( Remapper::TargetMesh, ctx.inFilenames[1], TempestRemapper::DEFAULT ),
1617  "Failed to load MOAB Target mesh" );
1618 
1619  // Generate overlap mesh
1620  TR_CHK_SET_ERR( GenerateOverlapWithMeshes( *ctx.meshes[0], *ctx.meshes[1], *tempest_mesh, "", "NetCDF4", "exact",
1621  false ),
1622  "Failed to generate TempestRemap OverlapMesh" );
1623 
1624  remapper.SetMesh( Remapper::OverlapMesh, tempest_mesh );
1625  ctx.meshes[2] = remapper.GetMesh( Remapper::OverlapMesh );
1626 
1627  return moab::MB_SUCCESS;
1628 }
1629 
1630 moab::ErrorCode handleOverlapMOAB( ToolContext& ctx, moab::TempestRemapper& remapper )
1631 {
1632  using namespace moab;
1633 
1634  // resize the meshsets and meshes vectors
1635  ctx.meshsets.resize( 3 );
1636  ctx.meshes.resize( 3 );
1637 
1638  ctx.meshsets[0] = remapper.GetMeshSet( Remapper::SourceMesh );
1639  ctx.meshsets[1] = remapper.GetMeshSet( Remapper::TargetMesh );
1640  ctx.meshsets[2] = remapper.GetMeshSet( Remapper::OverlapMesh );
1641 
1642  constexpr double radius_src = 1.0;
1643  constexpr double radius_dest = 1.0;
1644 
1645  // Load and process target mesh
1646  {
1647  std::vector< int > metadata;
1648  std::string additional_read_opts_tgt = ctx.get_file_read_options( ctx.inFilenames[1] );
1649  if( ctx.n_procs > 1 && ctx.disc_methods[1].compare( "fv" ) != 0 ) // target discretization is cgll or dgll
1650  {
1651  // auto pcomm = new ParallelComm( ctx.mbcore, MPI_COMM_WORLD );
1652  // add one ghost layer to the target mesh
1653  // additional_read_opts_tgt = additional_read_opts_tgt + "PARALLEL_GHOSTS=3.0.2;PARALLEL_THIN_GHOST_LAYER;SKIP_AUGMENT_WITH_GHOSTS;PRINT_PARALLEL;";
1654  // additional_read_opts_tgt = additional_read_opts_tgt + "PARALLEL_COMM=1;";
1655  // additional_read_opts_tgt = additional_read_opts_tgt + "PARALLEL_GHOSTS=3.0.1;";
1656  // additional_read_opts_tgt = additional_read_opts_tgt + "PARALLEL_COMM=" + std::to_string(ctx.pcomm->get_id()) + ";";
1657  }
1658 
1659  MB_CHK_SET_ERR( remapper.LoadNativeMesh( ctx.inFilenames[1], ctx.meshsets[1], metadata,
1660  additional_read_opts_tgt.c_str() ),
1661  "Failed to load MOAB Target mesh" );
1662 
1663 #ifdef MOAB_HAVE_MPI
1664  if( ctx.n_procs > 1 && ctx.disc_methods[1].compare( "fv" ) != 0 &&
1665  false ) // target discretization is cgll or dgll
1666  {
1667  Range beforeGhost, afterGhost;
1668  ctx.mbcore->get_entities_by_dimension( ctx.meshsets[1], 2, beforeGhost );
1669 
1670  ctx.pcomm->set_debug_verbosity( 5 );
1671  MB_CHK_SET_ERR( ctx.pcomm->exchange_ghost_cells( 2, 0, 1, 0, true, true, &ctx.meshsets[1] ),
1672  "Failed to exchange ghost cells for MOAB Target mesh" );
1673  ctx.pcomm->set_debug_verbosity( 0 );
1674 
1675  ctx.mbcore->get_entities_by_dimension( ctx.meshsets[1], 2, afterGhost );
1676  std::cout << ctx.proc_id << ": N(before) = " << beforeGhost.size() << ", N(after) = " << afterGhost.size()
1677  << std::endl;
1678 
1679  std::vector< Tag > taglist;
1680  taglist.push_back( ctx.mbcore->globalId_tag() );
1681  Tag gdofTag;
1682  MB_CHK_SET_ERR( ctx.mbcore->tag_get_handle( "GLOBAL_DOFS", gdofTag ),
1683  "Failed to get global dofs tag for MOAB Target mesh" );
1684  taglist.push_back( gdofTag );
1685  MB_CHK_SET_ERR( ctx.pcomm->exchange_tags( taglist, taglist, afterGhost ),
1686  "Failed to exchange global dofs for MOAB Target mesh" );
1687  // std::set< unsigned int > commprocs;
1688  // MB_CHK_SET_ERR( ctx.pcomm->get_comm_procs( commprocs ),
1689  // "Failed to get commprocs for MOAB Target mesh" );
1690  // if (ctx.proc_id == 0)
1691  // {
1692  // std::cout << ctx.proc_id << ": commprocs = [";
1693  // for( auto p : commprocs ) std::cout << p << ", ";
1694  // std::cout << "]\n";
1695 
1696  // std::cout << ctx.proc_id << ": N(after) = " << afterGhost.size() << std::endl;
1697  // for (auto eh: afterGhost)
1698  // {
1699  // std::cout << ctx.mbcore->type_from_handle(eh) << ": " << eh << std::endl;
1700  // }
1701  // }
1702  }
1703 #endif
1704 
1705  if( !metadata.empty() )
1706  {
1707  remapper.SetMeshType( Remapper::TargetMesh, metadata );
1708  }
1709 
1710  MB_CHK_SET_ERR( IntxUtils::ScaleToRadius( ctx.mbcore, ctx.meshsets[1], radius_dest ),
1711  "Failed to preprocess MOAB Target mesh" );
1712  }
1713 
1714  // Load and process source mesh
1715  {
1716  std::vector< int > metadata;
1717  auto additional_read_opts_src = ctx.get_file_read_options( ctx.inFilenames[0] );
1718 #ifdef MOAB_HAVE_MPI
1719  if( ctx.n_procs > 1 )
1720  {
1721  // auto pcomm = new ParallelComm( ctx.mbcore, MPI_COMM_WORLD );
1722  additional_read_opts_src =
1723  additional_read_opts_src + "PARALLEL_COMM=" + std::to_string( ctx.pcomm->get_id() ) + ";";
1724  }
1725 #endif
1726  MB_CHK_SET_ERR( remapper.LoadNativeMesh( ctx.inFilenames[0], ctx.meshsets[0], metadata,
1727  additional_read_opts_src.c_str() ),
1728  "Failed to load MOAB Source mesh" );
1729 
1730  if( !metadata.empty() )
1731  {
1732  remapper.SetMeshType( Remapper::SourceMesh, metadata );
1733  }
1734 
1735  MB_CHK_SET_ERR( IntxUtils::ScaleToRadius( ctx.mbcore, ctx.meshsets[0], radius_src ),
1736  "Failed to preprocess MOAB Source mesh" );
1737  }
1738 
1739  if( ctx.computeWeights )
1740  {
1741  // Convert MOAB to TempestRemap meshes
1742  MB_CHK_SET_ERR( remapper.ConvertMeshToTempest( Remapper::SourceMesh ),
1743  "Failed to convert MOAB Source mesh to TempestRemap mesh" );
1744  ctx.meshes[0] = remapper.GetMesh( Remapper::SourceMesh );
1745 
1746  MB_CHK_SET_ERR( remapper.ConvertMeshToTempest( Remapper::TargetMesh ),
1747  "Failed to convert MOAB Target mesh to TempestRemap mesh" );
1748  ctx.meshes[1] = remapper.GetMesh( Remapper::TargetMesh );
1749  }
1750 
1751  return moab::MB_SUCCESS;
1752 }
1753 
1754 moab::ErrorCode handleOverlapFiles( ToolContext& ctx, Mesh* tempest_mesh )
1755 {
1756  ctx.timer_push( "create Tempest OverlapMesh" );
1757  TR_CHK_SET_ERR( GenerateOverlapMesh( ctx.inFilenames[0], ctx.inFilenames[1], *tempest_mesh, ctx.outFilename,
1758  "NetCDF4", "exact", true ),
1759  "Failed to create Tempest OverlapMesh" );
1760  ctx.timer_pop();
1761 
1762  // Add the overlap mesh to the list of meshes
1763  ctx.meshes.push_back( tempest_mesh );
1764  return moab::MB_SUCCESS;
1765 }
1766 
1767 /**
1768  * @brief Convert a generated TempestRemap mesh to MOAB format and write as h5m file.
1769  *
1770  * When the output filename has a .h5m extension, the mesh is converted from TempestRemap
1771  * format to MOAB format in memory and written as a native MOAB HDF5 file. This allows
1772  * the generated meshes to be loaded by mbtempest type 5 (OVERLAP_MOAB) workflows.
1773  * The TempestRemap format file is still written (with .g extension) for compatibility.
1774  */
1775 moab::ErrorCode convertAndWriteMOABMesh( ToolContext& ctx, moab::TempestRemapper& remapper, Mesh* tempest_mesh )
1776 {
1777  // Check if output filename has .h5m extension
1778  const std::string& outFile = ctx.outFilename;
1779  const size_t dot = outFile.find_last_of( "." );
1780  if( dot == std::string::npos ) return moab::MB_SUCCESS;
1781 
1782  const std::string ext = outFile.substr( dot + 1 );
1783  if( ext != "h5m" ) return moab::MB_SUCCESS;
1784 
1785  // Register the TempestRemap mesh with the remapper as SourceMesh
1786  remapper.SetMesh( moab::Remapper::SourceMesh, tempest_mesh, false );
1787 
1788  // Convert TempestRemap mesh to MOAB format
1789  ctx.timer_push( "convert TempestRemap mesh to MOAB format" );
1790  MB_CHK_SET_ERR( remapper.ConvertTempestMesh( moab::Remapper::SourceMesh ),
1791  "Failed to convert TempestRemap mesh to MOAB format" );
1792  ctx.timer_pop();
1793 
1794  // Fix degenerate quads: RLL meshes from TempestRemap have polar cells stored as
1795  // 4-node quads with duplicate vertices. Convert these to proper triangles so the
1796  // intersection algorithm can handle them correctly.
1797  moab::EntityHandle meshSet = remapper.GetMeshSet( moab::Remapper::SourceMesh );
1798  MB_CHK_SET_ERR( moab::IntxUtils::fix_degenerate_quads( ctx.mbcore, meshSet ),
1799  "Failed to fix degenerate quads in converted mesh" );
1800  ctx.timer_push( "write MOAB mesh to h5m file" );
1801  MB_CHK_SET_ERR( ctx.mbcore->write_file( outFile.c_str(), nullptr, nullptr, &meshSet, 1 ),
1802  "Failed to write MOAB mesh to h5m file" );
1803  ctx.timer_pop();
1804 
1805  if( !ctx.proc_id )
1806  ctx.outputFormatter.printf( 0, "Wrote MOAB mesh to %s\n", outFile.c_str() );
1807 
1808  return moab::MB_SUCCESS;
1809 }
1810 
1811 moab::ErrorCode handleICOMesh( ToolContext& ctx, moab::TempestRemapper& remapper, Mesh* tempest_mesh )
1812 {
1813  std::string trFilename = ctx.outFilename;
1814  const size_t dot = trFilename.find_last_of( "." );
1815  if( dot != std::string::npos && trFilename.substr( dot + 1 ) == "h5m" )
1816  trFilename = trFilename.substr( 0, dot ) + ".g";
1817 
1818  ctx.timer_push( "generate ICO mesh with TempestRemap" );
1819  TR_CHK_SET_ERR( GenerateICOMesh( *tempest_mesh, ctx.blockSize, ctx.computeDual, trFilename, "NetCDF4" ),
1820  "Failed to generate ICO mesh with TempestRemap" );
1821  ctx.timer_pop();
1822 
1823  // Add the ICO mesh to the list of meshes
1824  ctx.meshes.push_back( tempest_mesh );
1825 
1826  MB_CHK_SET_ERR( convertAndWriteMOABMesh( ctx, remapper, tempest_mesh ),
1827  "Failed to convert and write MOAB mesh" );
1828 
1829  return moab::MB_SUCCESS;
1830 }
1831 
1832 moab::ErrorCode handleRLLMesh( ToolContext& ctx, moab::TempestRemapper& remapper, Mesh* tempest_mesh )
1833 {
1834  std::string trFilename = ctx.outFilename;
1835  const size_t dot = trFilename.find_last_of( "." );
1836  if( dot != std::string::npos && trFilename.substr( dot + 1 ) == "h5m" )
1837  trFilename = trFilename.substr( 0, dot ) + ".g";
1838 
1839  ctx.timer_push( "generate RLL mesh with TempestRemap" );
1840  TR_CHK_SET_ERR( GenerateRLLMesh( *tempest_mesh, // Mesh& meshOut,
1841  ctx.blockSize * 2, ctx.blockSize, // int nLongitudes, int nLatitudes,
1842  0.0, 360.0, // double dLonBegin, double dLonEnd,
1843  -90.0, 90.0, // double dLatBegin, double dLatEnd,
1844  false, false, false, // bool fGlobalCap, bool fFlipLatLon, bool fForceGlobal,
1845  "" /*ctx.inFilename*/,
1846  "", // std::string strInputFile, std::string strInputFileLonName
1847  "", // std::string strInputFileLatName
1848  trFilename, // std::string strOutputFile
1849  "NetCDF4", // std::string strOutputFormat
1850  true // bool fVerbose
1851  ),
1852  "Failed to generate RLL mesh with TempestRemap" );
1853  ctx.timer_pop();
1854 
1855  // Add the RLL mesh to the list of meshes
1856  ctx.meshes.push_back( tempest_mesh );
1857 
1858  MB_CHK_SET_ERR( convertAndWriteMOABMesh( ctx, remapper, tempest_mesh ),
1859  "Failed to convert and write MOAB mesh" );
1860 
1861  return moab::MB_SUCCESS;
1862 }
1863 
1864 moab::ErrorCode handleCSMesh( ToolContext& ctx, moab::TempestRemapper& remapper, Mesh* tempest_mesh )
1865 {
1866  // Generate the TempestRemap Exodus mesh (always written as .g for TempestRemap compatibility)
1867  std::string trFilename = ctx.outFilename;
1868  const size_t dot = trFilename.find_last_of( "." );
1869  if( dot != std::string::npos && trFilename.substr( dot + 1 ) == "h5m" )
1870  trFilename = trFilename.substr( 0, dot ) + ".g";
1871 
1872  ctx.timer_push( "generate CS mesh with TempestRemap" );
1873  TR_CHK_SET_ERR( GenerateCSMesh( *tempest_mesh, ctx.blockSize, trFilename, "NetCDF4" ),
1874  "Failed to generate CS mesh with TempestRemap" );
1875  ctx.timer_pop();
1876 
1877  // Add the CS mesh to the list of meshes
1878  ctx.meshes.push_back( tempest_mesh );
1879 
1880  // Convert and write as MOAB h5m if requested
1881  MB_CHK_SET_ERR( convertAndWriteMOABMesh( ctx, remapper, tempest_mesh ),
1882  "Failed to convert and write MOAB mesh" );
1883 
1884  return moab::MB_SUCCESS;
1885 }
1886 
1887 } // namespace
1888 
1889 /**
1890  * @brief Creates a TempestRemap mesh based on the provided context and mesh type
1891  *
1892  * @param ctx Tool context containing configuration and state
1893  * @param remapper TempestRemap instance for mesh operations
1894  * @param tempest_mesh Output parameter for the created mesh
1895  * @return moab::ErrorCode Status of the operation
1896  */
1897 static moab::ErrorCode CreateTempestMesh( ToolContext& ctx, moab::TempestRemapper& remapper, Mesh* tempest_mesh )
1898 {
1899  using namespace moab;
1900  using RemapperType = moab::TempestRemapper;
1901 
1902  auto& outputFormatter = ctx.outputFormatter;
1903 
1904  try
1905  {
1906  switch( ctx.meshType )
1907  {
1908  case RemapperType::OVERLAP_FILES:
1909  if( !ctx.proc_id ) outputFormatter.printf( 0, "Creating TempestRemap overlap mesh ...\n" );
1910  return handleOverlapFiles( ctx, tempest_mesh );
1911 
1912  case RemapperType::OVERLAP_MEMORY:
1913  if( !ctx.proc_id )
1914  outputFormatter.printf( 0, "Convert MOAB overlap files to TempestRemap format in-memory ...\n" );
1915  return handleOverlapMemory( ctx, remapper, tempest_mesh );
1916 
1917  case RemapperType::OVERLAP_MOAB:
1918  if( !ctx.proc_id )
1919  outputFormatter.printf( 0, "Convert MOAB meshes to TempestRemap format in-memory ...\n" );
1920  return handleOverlapMOAB( ctx, remapper );
1921 
1922  case RemapperType::ICO:
1923  if( !ctx.proc_id ) outputFormatter.printf( 0, "Creating TempestRemap ICO mesh ...\n" );
1924  return handleICOMesh( ctx, remapper, tempest_mesh );
1925 
1926  case RemapperType::RLL:
1927  if( !ctx.proc_id ) outputFormatter.printf( 0, "Creating TempestRemap RLL mesh ...\n" );
1928  return handleRLLMesh( ctx, remapper, tempest_mesh );
1929 
1930  default: // Default to CS mesh
1931  if( !ctx.proc_id ) outputFormatter.printf( 0, "Creating TempestRemap CS mesh ...\n" );
1932  return handleCSMesh( ctx, remapper, tempest_mesh );
1933  }
1934  }
1935  catch( const std::exception& e )
1936  {
1937  std::cerr << "Error in CreateTempestMesh: " << e.what() << "\n";
1938  return MB_FAILURE;
1939  }
1940 }
1941 
1942 #undef MOAB_DBG
1943 
1944 ///////////////////////////////////////////////
1945 /**
1946  * @brief Sample functions for testing remapping operations
1947  *
1948  * These functions provide analytical test cases with different spatial patterns
1949  * for verifying remapping accuracy and performance.
1950  */
1951 
1952 // Constants for sample functions
1953 namespace
1954 {
1955 // Constants for sample_stationary_vortex
1956 constexpr double VORTEX_LON0 = 0.0;
1957 constexpr double VORTEX_LAT0 = 0.6;
1958 constexpr double VORTEX_R0 = 3.0;
1959 constexpr double VORTEX_D = 5.0;
1960 constexpr double VORTEX_T = 6.0;
1961 
1962 } // namespace
1963 
1964 /**
1965  * @brief Constant sample function
1966  *
1967  * @return double Always returns 1.0
1968  */
1969 static inline constexpr double sample_constant( double /*dLon*/, double /*dLat*/ ) noexcept
1970 {
1971  return 1.0;
1972 }
1973 
1974 /**
1975  * @brief Sample function with slow harmonic variation
1976  *
1977  * @param dLon Longitude in radians
1978  * @param dLat Latitude in radians
1979  * @return double Function value at (dLon, dLat)
1980  */
1981 static inline double sample_slow_harmonic( double dLon, double dLat ) noexcept
1982 {
1983  const double cosLat = std::cos( dLat );
1984  return 2.0 + cosLat * cosLat * std::cos( 2.0 * dLon );
1985 }
1986 
1987 /**
1988  * @brief Sample function with fast harmonic variation
1989  *
1990  * @param dLon Longitude in radians
1991  * @param dLat Latitude in radians
1992  * @return double Function value at (dLon, dLat)
1993  */
1994 static inline double sample_fast_harmonic( double dLon, double dLat ) noexcept
1995 {
1996  const double sin2Lat = std::sin( 2.0 * dLat );
1997  return 2.0 +
1998  sin2Lat * sin2Lat * sin2Lat * sin2Lat * sin2Lat * sin2Lat * sin2Lat * sin2Lat * std::cos( 16.0 * dLon );
1999 }
2000 
2001 /**
2002  * @brief Sample function representing a stationary vortex
2003  *
2004  * @param dLon Longitude in radians
2005  * @param dLat Latitude in radians
2006  * @return double Function value at (dLon, dLat)
2007  */
2008 static inline double sample_stationary_vortex( double dLon, double dLat ) noexcept
2009 {
2010  // Find the rotated longitude and latitude of a point on a sphere
2011  // with pole at (dLonC, dLatC)
2012  const double dSinC = std::sin( VORTEX_LAT0 );
2013  const double dCosC = std::cos( VORTEX_LAT0 );
2014  const double dSinT = std::sin( dLat );
2015  const double dCosT = std::cos( dLat );
2016 
2017  const double dTrm = dCosT * std::cos( dLon - VORTEX_LON0 );
2018  const double dX = dSinC * dTrm - dCosC * dSinT;
2019  const double dY = dCosT * std::sin( dLon - VORTEX_LON0 );
2020  const double dZ = dSinC * dSinT + dCosC * dTrm;
2021 
2022  // Calculate new longitude and latitude in rotated coordinate system
2023  double dNewLon = std::atan2( dY, dX );
2024  if( dNewLon < 0.0 )
2025  {
2026  dNewLon += 2.0 * M_PI;
2027  }
2028  const double dNewLat = std::asin( dZ );
2029 
2030  // Calculate vortex profile
2031  const double dRho = VORTEX_R0 * std::cos( dNewLat );
2032  const double dVt = 3.0 * std::sqrt( 3.0 ) / 2.0 / std::cosh( dRho ) / std::cosh( dRho ) * std::tanh( dRho );
2033 
2034  // Calculate angular velocity (avoid division by zero)
2035  const double dOmega = ( dRho == 0.0 ) ? 0.0 : ( dVt / dRho );
2036 
2037  // Return the final vortex profile
2038  return ( 1.0 - std::tanh( dRho / VORTEX_D * std::sin( dNewLon - dOmega * VORTEX_T ) ) );
2039 }
2040 
2041 ///////////////////////////////////////////////