Mesh Oriented datABase  (version 5.5.1)
An array-based unstructured mesh library
mbtempest.cpp
Go to the documentation of this file.
1 /*
2  * Usage: MOAB-Tempest tool
3  *
4  * Generate a Cubed-Sphere mesh: ./mbtempest -t 0 -res 25 -f cubed_sphere_mesh.exo
5  * Generate a RLL mesh: ./mbtempest -t 1 -res 25 -f rll_mesh.exo
6  * Generate a Icosahedral-Sphere mesh: ./mbtempest -t 2 -res 25 <-dual> -f icosa_mesh.exo
7  *
8  * Now you can compute the intersections between the meshes too!
9  *
10  * Generate the overlap mesh: ./mbtempest -t 3 -l cubed_sphere_mesh.exo -l rll_mesh.exo -f
11  * overlap_mesh.exo
12  *
13  */
14 
15 // standard C++ includes
16 #include <iostream>
17 #include <iomanip>
18 #include <cstdlib>
19 #include <vector>
20 #include <string>
21 #include <sstream>
22 #include <cassert>
23 
24 // MOAB includes
25 #include "moab/Core.hpp"
26 #include "moab/IntxMesh/IntxUtils.hpp"
27 #include "moab/Remapping/TempestRemapper.hpp"
28 #include "moab/Remapping/TempestOnlineMap.hpp"
29 #include "moab/ProgOptions.hpp"
30 #include "moab/CpuTimer.hpp"
31 #include "DebugOutput.hpp"
32 
33 //#ifndef MOAB_HAVE_MPI
34 // #error mbtempest tool requires MPI configuration
35 //#endif
36 
37 #ifdef MOAB_HAVE_MPI
38 // MPI includes
39 #include "moab_mpi.h"
40 #include "moab/ParallelComm.hpp"
41 #include "MBParallelConventions.h"
42 #endif
43 
44 struct ToolContext
45 {
46  moab::Core* mbcore;
47 #ifdef MOAB_HAVE_MPI
48  moab::ParallelComm* pcomm;
49 #endif
50  const int proc_id, n_procs; // MPI process id and number of processes
51  moab::DebugOutput outputFormatter; // Output formatter
52  int blockSize; // Resolution of the mesh
53  int nlayers; // Number of ghost layers in the source mesh
54  std::vector< std::string > inFilenames; // Input mesh filenames
55  std::vector< Mesh* > meshes; // Mesh references in TempestRemap format
56  std::vector< moab::EntityHandle > meshsets; // Meshsets in MOAB format
57  std::vector< int > disc_orders; // Discretization orders
58  std::vector< std::string > disc_methods; // Discretization method names (fv, cgll, dgll, etc.)
59  std::vector< std::string > doftag_names; // DoF tag names (GLOBAL_ID, etc.)
60  std::string fvMethod; // FV method name (invdist, delaunay, bilin, intbilin, intbilingb, none)
61  std::string outFilename; // Output filename
62  std::string intxFilename; // Output intersection mesh filename
63  std::string baselineFile; // Output baseline file (for verification)
64  std::string variableToVerify; // Variable name for verification
65  moab::TempestRemapper::TempestMeshType
66  meshType; // Mesh type (CS, RLL, ICO, OVERLAP_FILES, OVERLAP_MEMORY, OVERLAP_MOAB)
67  bool skip_io; // Skip I/O operations; strictly to measure the performance of the intersection/map/SpMV algorithms
68  bool computeDual; // Compute the dual of the mesh
69  bool computeWeights; // Compute the projection weights
70  bool verifyWeights; // Verify the projection weights
71  bool enforceConvexity; // Enforce convexity of the input meshes
72  int ensureMonotonicity; // Ensure monotonicity in the weight generation process
73  bool rrmGrids; // Flag specifying that we are dealing with regionally refined grids (possibly nonoverlapping)
74  bool kdtreeSearch; // Use Kd-tree based search for computing mesh intersections instead of advancing front
75  bool fCheck; // Check the generated map for conservation and consistency
76  bool fVolumetric; // Apply a volumetric projection to compute the weights
77  bool useGnomonicProjection; // Use Gnomonic plane projections to compute coverage mesh
78  moab::TempestOnlineMap::CAASType cassType; // CAAS filter type
79  GenerateOfflineMapAlgorithmOptions mapOptions; // Offline map options
80  bool print_diagnostics; // Print diagnostics
81  double boxeps; // Box error tolerance
82  double epsrel; // Relative error tolerance
83 
84 #ifdef MOAB_HAVE_MPI
85  ToolContext( moab::Core* icore, moab::ParallelComm* p_pcomm )
86  : mbcore( icore ), pcomm( p_pcomm ), proc_id( pcomm->rank() ), n_procs( pcomm->size() ),
87  outputFormatter( std::cout, pcomm->rank(), 0 ),
88 #else
89  ToolContext( moab::Core* icore )
90  : mbcore( icore ), proc_id( 0 ), n_procs( 1 ), outputFormatter( std::cout, 0, 0 ),
91 #endif
92  blockSize( 5 ), nlayers( 0 ), fvMethod( "none" ), outFilename( "outputFile.nc" ), intxFilename( "" ),
93  baselineFile( "" ), variableToVerify( "" ), meshType( moab::TempestRemapper::DEFAULT ), skip_io( false ),
94  computeDual( false ), computeWeights( false ), verifyWeights( false ), enforceConvexity( false ),
95  ensureMonotonicity( 0 ), rrmGrids( false ), kdtreeSearch( true ), fCheck( false ), fVolumetric( false ),
96  useGnomonicProjection( false ), cassType( moab::TempestOnlineMap::CAAS_NONE ), print_diagnostics( false ),
97  boxeps( 1e-7 ), // Box error tolerance default value
98  epsrel( ReferenceTolerance ) // ReferenceTolerance is defined in Defines.h in TempestRemap
99  {
100  inFilenames.resize( 2 );
101  doftag_names.resize( 2 );
102  timer = new moab::CpuTimer();
103 
104  outputFormatter.set_prefix( "[MBTempest]: " );
105  }
106 
107  ~ToolContext()
108  {
109  // for (unsigned i=0; i < meshes.size(); ++i) delete meshes[i];
110  meshes.clear();
111  inFilenames.clear();
112  disc_orders.clear();
113  disc_methods.clear();
114  doftag_names.clear();
115  outFilename.clear();
116  intxFilename.clear();
117  baselineFile.clear();
118  variableToVerify.clear();
119  meshsets.clear();
120  delete timer;
121  }
122 
123  void timer_push( std::string operation )
124  {
125  timer_ops = timer->time_since_birth();
126  opName = operation;
127  }
128 
129  void timer_pop()
130  {
131  double locElapsed = timer->time_since_birth() - timer_ops, avgElapsed = 0, maxElapsed = 0;
132 #ifdef MOAB_HAVE_MPI
133  MPI_Reduce( &locElapsed, &maxElapsed, 1, MPI_DOUBLE, MPI_MAX, 0, pcomm->comm() );
134  MPI_Reduce( &locElapsed, &avgElapsed, 1, MPI_DOUBLE, MPI_SUM, 0, pcomm->comm() );
135 #else
136  maxElapsed = locElapsed;
137  avgElapsed = locElapsed;
138 #endif
139  if( !proc_id )
140  {
141  avgElapsed /= n_procs;
142  std::cout << "[LOG] Time taken to " << opName.c_str() << ": max = " << maxElapsed
143  << ", avg = " << avgElapsed << "\n";
144  }
145  // std::cout << "\n[LOG" << proc_id << "] Time taken to " << opName << " = " <<
146  // timer->time_since_birth() - timer_ops << std::endl;
147  opName.clear();
148  }
149 
150  void ParseCLOptions( int argc, char* argv[] )
151  {
152  ProgOptions opts;
153  int imeshType = 0;
154  std::string expectedFName = "output.exo";
155  std::string expectedMethod = "fv";
156  std::string expectedFVMethod = "none";
157  std::string expectedDofTagName = "GLOBAL_ID";
158  int expectedOrder = 1;
159  int useCAAS = 0;
160  int nlayer_input = 0;
161 
162  if( !proc_id )
163  {
164  std::cout << "Command line options provided to mbtempest:\n ";
165  for( int iarg = 0; iarg < argc; ++iarg )
166  std::cout << argv[iarg] << " ";
167  std::cout << std::endl << std::endl;
168  }
169 
170  opts.addOpt< int >( "type,t",
171  "Type of mesh (default=CS; Choose from [CS=0, RLL=1, ICO=2, OVERLAP_FILES=3, "
172  "OVERLAP_MEMORY=4, OVERLAP_MOAB=5])",
173  &imeshType );
174 
175  opts.addOpt< int >( "res,r", "Resolution of the mesh (default=5)", &blockSize );
176 
177  opts.addOpt< void >( "dual,d", "Output the dual of the mesh (relevant only for ICO mesh type)", &computeDual );
178 
179  opts.addOpt< std::string >( "file,f", "Output computed mesh or remapping weights to specified filename",
180  &outFilename );
181  opts.addOpt< std::string >(
182  "load,l", "Input mesh filenames for source and target meshes. (relevant only when computing weights)",
183  &expectedFName );
184 
185  opts.addOpt< void >( "advfront,a",
186  "Use the advancing front intersection instead of the Kd-tree based algorithm "
187  "to compute mesh intersections. (relevant only when computing weights)" );
188 
189  opts.addOpt< std::string >( "intx,i", "Output TempestRemap intersection mesh filename", &intxFilename );
190 
191  opts.addOpt< void >( "weights,w",
192  "Compute and output the weights using the overlap mesh (generally "
193  "relevant only for OVERLAP mesh)",
194  &computeWeights );
195 
196  opts.addOpt< void >( "verbose,v",
197  "Print verbose diagnostic messages during intersection and map computation (default=false)",
198  &print_diagnostics );
199 
200  opts.addOpt< std::string >( "method,m", "Discretization method for the source and target solution fields",
201  &expectedMethod );
202 
203  opts.addOpt< int >( "order,o", "Discretization orders for the source and target solution fields",
204  &expectedOrder );
205 
206  opts.addOpt< std::string >( "global_id,g",
207  "Tag name that contains the global DoF IDs for source and target solution fields",
208  &expectedDofTagName );
209 
210  opts.addOpt< std::string >( "fvmethod",
211  "Sub-type method for FV-FV projections (invdist, delaunay, bilin, "
212  "intbilin, intbilingb, none. Default: none)",
213  &expectedFVMethod );
214 
215  opts.addOpt< void >( "noconserve",
216  "Do not apply conservation to the resultant weights (relevant only "
217  "when computing weights)",
218  &mapOptions.fNoConservation );
219 
220  opts.addOpt< void >( "volumetric",
221  "Apply a volumetric projection to compute the weights (relevant only "
222  "when computing weights)",
223  &fVolumetric );
224 
225  opts.addOpt< void >( "skip_output", "For performance studies, skip all I/O operations.", &skip_io );
226 
227  opts.addOpt< void >( "gnomonic", "Use Gnomonic plane projections to compute coverage mesh.",
228  &useGnomonicProjection );
229 
230  opts.addOpt< void >( "enforce_convexity", "check convexity of input meshes to compute mesh intersections",
231  &enforceConvexity );
232 
233  opts.addOpt< void >( "nobubble", "do not use bubble on interior of spectral element nodes",
234  &mapOptions.fNoBubble );
235 
236  opts.addOpt< void >(
237  "sparseconstraints",
238  "Use sparse solver for constraints when we have high-valence (typical with high-res RLL mesh)",
239  &mapOptions.fSparseConstraints );
240 
241  opts.addOpt< void >( "rrmgrids",
242  "At least one of the meshes is a regionally refined grid (relevant to "
243  "accelerate intersection computation)",
244  &rrmGrids );
245 
246  opts.addOpt< void >( "checkmap", "Check the generated map for conservation and consistency", &fCheck );
247 
248  opts.addOpt< void >( "verify",
249  "Verify the accuracy of the maps by projecting analytical functions "
250  "from source to target "
251  "grid by applying the maps",
252  &verifyWeights );
253 
254  opts.addOpt< std::string >( "var",
255  "Tag name of the variable to use in the verification study "
256  "(error metrics for user defined variables may not be available)",
257  &variableToVerify );
258 
259  opts.addOpt< int >( "monotonicity", "Ensure monotonicity in the weight generation. Options=[0,1,2,3]",
260  &ensureMonotonicity );
261 
262  opts.addOpt< int >( "ghost", "Number of ghost layers in coverage mesh (default=1)", &nlayer_input );
263 
264  opts.addOpt< double >( "boxeps", "The tolerance for boxes (default=1e-7)", &boxeps );
265 
266  opts.addOpt< int >( "caas", "apply CAAS nonlinear filter after linear map application", &useCAAS );
267 
268  opts.addOpt< std::string >( "baseline", "Output baseline file", &baselineFile );
269 
270  opts.parseCommandLine( argc, argv );
271 
272  // By default - use Kd-tree based search; if user asks for advancing front, disable Kd-tree
273  // algorithm
274  kdtreeSearch = opts.numOptSet( "advfront,a" ) == 0;
275 
276  switch( imeshType )
277  {
278  case 0:
279  meshType = moab::TempestRemapper::CS;
280  break;
281 
282  case 1:
283  meshType = moab::TempestRemapper::RLL;
284  break;
285 
286  case 2:
287  meshType = moab::TempestRemapper::ICO;
288  break;
289 
290  case 3:
291  meshType = moab::TempestRemapper::OVERLAP_FILES;
292  break;
293 
294  case 4:
295  meshType = moab::TempestRemapper::OVERLAP_MEMORY;
296  break;
297 
298  case 5:
299  meshType = moab::TempestRemapper::OVERLAP_MOAB;
300  break;
301 
302  default:
303  meshType = moab::TempestRemapper::DEFAULT;
304  break;
305  }
306 
307  // decipher whether we want to use CAAS filter when applying the projection
308  switch( useCAAS )
309  {
310  case moab::TempestOnlineMap::CAAS_GLOBAL:
311  cassType = moab::TempestOnlineMap::CAAS_GLOBAL;
312  break;
313  case moab::TempestOnlineMap::CAAS_LOCAL:
314  cassType = moab::TempestOnlineMap::CAAS_LOCAL;
315  break;
316  case moab::TempestOnlineMap::CAAS_LOCAL_ADJACENT:
317  cassType = moab::TempestOnlineMap::CAAS_LOCAL_ADJACENT;
318  break;
319  case moab::TempestOnlineMap::CAAS_QLT:
320  cassType = moab::TempestOnlineMap::CAAS_QLT;
321  break;
322  default:
323  cassType = moab::TempestOnlineMap::CAAS_NONE;
324  break;
325  }
326 
327  if( meshType > moab::TempestRemapper::ICO ) // compute overlap mesh and maps possibly
328  {
329  opts.getOptAllArgs( "load,l", inFilenames );
330  opts.getOptAllArgs( "order,o", disc_orders );
331  opts.getOptAllArgs( "method,m", disc_methods );
332  opts.getOptAllArgs( "global_id,i", doftag_names );
333 
334  assert( inFilenames.size() == 2 );
335  assert( ensureMonotonicity >= 0 && ensureMonotonicity <= 3 );
336 
337  // get discretization order parameters
338  if( disc_orders.size() == 0 ) disc_orders.resize( 2, 1 );
339  if( disc_orders.size() == 1 ) disc_orders.push_back( disc_orders[0] );
340 
341  // get discretization method parameters
342  if( disc_methods.size() == 0 ) disc_methods.resize( 2, "fv" );
343  if( disc_methods.size() == 1 ) disc_methods.push_back( disc_methods[0] );
344 
345  // get DoF tagname parameters
346  if( doftag_names.size() == 0 ) doftag_names.resize( 2, "GLOBAL_ID" );
347  if( doftag_names.size() == 1 ) doftag_names.push_back( doftag_names[0] );
348 
349  assert( disc_orders.size() == 2 );
350  assert( disc_methods.size() == 2 );
351  assert( doftag_names.size() == 2 );
352 
353  // for computing maps and overlaps, set discretization orders
354  mapOptions.nPin = disc_orders[0];
355  mapOptions.nPout = disc_orders[1];
356  mapOptions.fSourceConcave = false;
357  mapOptions.fTargetConcave = false;
358 
359  mapOptions.strMethod = "";
360 
361  if( expectedFVMethod != "none" )
362  {
363  mapOptions.strMethod += expectedFVMethod + ";";
364  fvMethod = expectedFVMethod;
365 
366  // These FV projection methods are non-conservative; specify it explicitly
367  mapOptions.fNoConservation = true;
368  }
369  switch( ensureMonotonicity )
370  {
371  case 0:
372  mapOptions.fMonotone = false;
373  break;
374  case 3:
375  mapOptions.strMethod += "mono3;";
376  break;
377  case 2:
378  mapOptions.strMethod += "mono2;";
379  break;
380  default:
381  mapOptions.fMonotone = true;
382  }
383  mapOptions.fNoCorrectAreas = false;
384  mapOptions.fNoCheck = !fCheck;
385 
386  //assert( fVolumetric && fInverseDistanceMap == false ); // both options cannot be active
387  if( fVolumetric ) mapOptions.strMethod += "volumetric;";
388 
389  // For global meshes, this default should work out of the box.
390  if( !fvMethod.compare( "bilin" ) ) nlayers = 3;
391  else
392  nlayers = ( mapOptions.nPin > 1 ? mapOptions.nPin + 1 : 0 );
393  if ( nlayer_input )
394  nlayers = std::max( nlayer_input, nlayers );
395  }
396 
397  // clear temporary string name
398  expectedFName.clear();
399 
400  mapOptions.strOutputMapFile = outFilename;
401  mapOptions.strOutputFormat = "Netcdf4";
402  }
403 
404  private:
405  moab::CpuTimer* timer;
406  double timer_ops;
407  std::string opName;
408 };
409 
410 // Forward declare some methods
411 static moab::ErrorCode CreateTempestMesh( ToolContext&, moab::TempestRemapper& remapper, Mesh* );
412 inline double sample_slow_harmonic( double dLon, double dLat );
413 inline double sample_fast_harmonic( double dLon, double dLat );
414 inline double sample_constant( double dLon, double dLat );
415 inline double sample_stationary_vortex( double dLon, double dLat );
416 
417 std::string get_file_read_options( ToolContext& ctx, std::string filename )
418 {
419  // if running in serial, blank option may suffice in general
420  std::string opts = "";
421  if( ctx.n_procs > 1 )
422  {
423  size_t lastindex = filename.find_last_of( "." );
424  std::string extension = filename.substr( lastindex + 1, filename.size() );
425  if( extension == "h5m" ) return "PARALLEL=READ_PART;PARTITION=PARALLEL_PARTITION;PARALLEL_RESOLVE_SHARED_ENTS;";
426  // "PARALLEL_GHOSTS=2.0.3;PARALLEL_THIN_GHOST_LAYER;";
427  else if( extension == "nc" )
428  {
429  // set default set of options
430  opts = "PARALLEL=READ_PART;PARTITION_METHOD=RCBZOLTAN;";
431 
432  if( ctx.proc_id == 0 )
433  {
434  {
435  NcFile ncFile( filename.c_str(), NcFile::ReadOnly );
436  if( !ncFile.is_valid() )
437  _EXCEPTION1( "Unable to open grid file \"%s\" for reading", filename.c_str() );
438 
439  // Check for dimension names "grid_size", "grid_rank" and "grid_corners"
440  int iSCRIPFormat = 0, iMPASFormat = 0, iESMFFormat = 0;
441  for( int i = 0; i < ncFile.num_dims(); i++ )
442  {
443  std::string strDimName = ncFile.get_dim( i )->name();
444  if( strDimName == "grid_size" || strDimName == "grid_corners" || strDimName == "grid_rank" )
445  iSCRIPFormat++;
446  if( strDimName == "nodeCount" || strDimName == "elementCount" ||
447  strDimName == "maxNodePElement" )
448  iESMFFormat++;
449  if( strDimName == "nCells" || strDimName == "nEdges" || strDimName == "nVertices" ||
450  strDimName == "vertexDegree" )
451  iMPASFormat++;
452  }
453  ncFile.close();
454 
455  if( iESMFFormat == 3 ) // Input from a NetCDF ESMF file
456  {
457  opts = "PARALLEL=READ_PART;PARTITION_METHOD=RCBZOLTAN;PARALLEL_RESOLVE_SHARED_ENTS;VARIABLE=;";
458  }
459  if( iSCRIPFormat == 3 ) // Input from a NetCDF SCRIP file
460  {
461  opts = "PARALLEL=READ_PART;PARTITION_METHOD=RCBZOLTAN;";
462  }
463  if( iMPASFormat == 4 ) // Input from a NetCDF SCRIP file
464  {
465  opts = "PARALLEL=READ_PART;PARTITION_METHOD=RCBZOLTAN;"
466  "PARALLEL_RESOLVE_SHARED_ENTS;NO_EDGES;NO_MIXED_ELEMENTS;VARIABLE=;";
467  }
468  }
469  }
470 
471 #ifdef MOAB_HAVE_MPI
472  int line_size = opts.size();
473  MPI_Bcast( &line_size, 1, MPI_INT, 0, MPI_COMM_WORLD );
474  if( ctx.proc_id != 0 ) opts.resize( line_size );
475  MPI_Bcast( const_cast< char* >( opts.data() ), line_size, MPI_CHAR, 0, MPI_COMM_WORLD );
476 #endif
477  }
478  else
479  return "PARALLEL=BCAST_DELETE;PARTITION=TRIVIAL;PARALLEL_RESOLVE_SHARED_ENTS;";
480  }
481  return opts;
482 }
483 //#define MOAB_DBG
484 int main( int argc, char* argv[] )
485 {
486  moab::ErrorCode rval;
487  NcError error( NcError::verbose_nonfatal );
488  std::stringstream sstr;
489  std::string historyStr;
490 
491  int proc_id = 0, nprocs = 1;
492 #ifdef MOAB_HAVE_MPI
493  MPI_Init( &argc, &argv );
494  MPI_Comm_rank( MPI_COMM_WORLD, &proc_id );
495  MPI_Comm_size( MPI_COMM_WORLD, &nprocs );
496 #endif
497 
498  moab::Core* mbCore = new( std::nothrow ) moab::Core;
499 
500  if( nullptr == mbCore )
501  {
502  return 1;
503  }
504 
505  // Build the history string
506  for( int ia = 0; ia < argc; ++ia )
507  historyStr += std::string( argv[ia] ) + " ";
508 
509  ToolContext* runCtx;
510 #ifdef MOAB_HAVE_MPI
511  moab::ParallelComm* pcomm = new moab::ParallelComm( mbCore, MPI_COMM_WORLD, 0 );
512 
513  runCtx = new ToolContext( mbCore, pcomm );
514  const char* writeOptions = ( nprocs > 1 ? "PARALLEL=WRITE_PART" : "" );
515 #else
516  runCtx = new ToolContext( mbCore );
517  const char* writeOptions = "";
518 #endif
519  runCtx->ParseCLOptions( argc, argv );
520 
521  const double radius_src = 1.0 /*2.0*acos(-1.0)*/;
522  const double radius_dest = 1.0 /*2.0*acos(-1.0)*/;
523 
524  moab::DebugOutput& outputFormatter = runCtx->outputFormatter;
525 
526 #ifdef MOAB_HAVE_MPI
527  moab::TempestRemapper remapper( mbCore, pcomm );
528 #else
529  moab::TempestRemapper remapper( mbCore );
530 #endif
531  remapper.meshValidate = true;
532  remapper.constructEdgeMap = true;
533  remapper.initialize();
534 
535  // Default area_method = lHuiller; Options: Girard, lHuiller, GaussQuadrature (if TR is available)
536  moab::IntxAreaUtils areaAdaptor( moab::IntxAreaUtils::GaussQuadrature );
537 
538  Mesh* tempest_mesh = new Mesh();
539  rval = CreateTempestMesh( *runCtx, remapper, tempest_mesh );MB_CHK_ERR( rval );
540 
541  if( runCtx->meshType == moab::TempestRemapper::OVERLAP_MEMORY )
542  {
543  // Compute intersections with MOAB
544  // For the overlap method, choose between: "fuzzy", "exact" or "mixed"
545  assert( runCtx->meshes.size() == 3 );
546 
547 #ifdef MOAB_HAVE_MPI
548  rval = pcomm->check_all_shared_handles();MB_CHK_ERR( rval );
549 #endif
550 
551  // Load the meshes and validate
552  rval = remapper.ConvertTempestMesh( moab::Remapper::SourceMesh );MB_CHK_ERR( rval );
553  rval = remapper.ConvertTempestMesh( moab::Remapper::TargetMesh );MB_CHK_ERR( rval );
554  rval = remapper.ConvertTempestMesh( moab::Remapper::OverlapMesh );MB_CHK_ERR( rval );
555  if( !runCtx->skip_io )
556  {
557  rval = mbCore->write_mesh( "tempest_intersection.h5m", &runCtx->meshsets[2], 1 );MB_CHK_ERR( rval );
558  }
559 
560  // print verbosely about the problem setting
561  size_t velist[6], gvelist[6];
562  {
563  moab::Range rintxverts, rintxelems;
564  rval = mbCore->get_entities_by_dimension( runCtx->meshsets[0], 0, rintxverts );MB_CHK_ERR( rval );
565  rval = mbCore->get_entities_by_dimension( runCtx->meshsets[0], 2, rintxelems );MB_CHK_ERR( rval );
566  velist[0] = rintxverts.size();
567  velist[1] = rintxelems.size();
568 
569  moab::Range bintxverts, bintxelems;
570  rval = mbCore->get_entities_by_dimension( runCtx->meshsets[1], 0, bintxverts );MB_CHK_ERR( rval );
571  rval = mbCore->get_entities_by_dimension( runCtx->meshsets[1], 2, bintxelems );MB_CHK_ERR( rval );
572  velist[2] = bintxverts.size();
573  velist[3] = bintxelems.size();
574  }
575 
576  moab::EntityHandle intxset; // == remapper.GetMeshSet(moab::Remapper::OverlapMesh);
577 
578  // Compute intersections with MOAB
579  {
580  // Create the intersection on the sphere object
581  runCtx->timer_push( "setup the intersector" );
582 
583  moab::Intx2MeshOnSphere* mbintx = new moab::Intx2MeshOnSphere( mbCore );
584  mbintx->set_error_tolerance( runCtx->epsrel );
585  mbintx->set_box_error( runCtx->boxeps );
586  mbintx->set_radius_source_mesh( radius_src );
587  mbintx->set_radius_destination_mesh( radius_dest );
588 #ifdef MOAB_HAVE_MPI
589  mbintx->set_parallel_comm( pcomm );
590 #endif
591  rval = mbintx->FindMaxEdges( runCtx->meshsets[0], runCtx->meshsets[1] );MB_CHK_ERR( rval );
592 
593 #ifdef MOAB_HAVE_MPI
594  moab::Range local_verts;
595  rval = mbintx->build_processor_euler_boxes( runCtx->meshsets[1], local_verts );MB_CHK_ERR( rval );
596 
597  runCtx->timer_pop();
598 
599  moab::EntityHandle covering_set;
600  runCtx->timer_push( "communicate the mesh" );
601  rval = mbintx->construct_covering_set( runCtx->meshsets[0], covering_set );MB_CHK_ERR( rval ); // lots of communication if mesh is distributed very differently
602  runCtx->timer_pop();
603 
604  // print verbosely about the problem setting
605  {
606  moab::Range cintxverts, cintxelems;
607  rval = mbCore->get_entities_by_dimension( covering_set, 0, cintxverts );MB_CHK_ERR( rval );
608  rval = mbCore->get_entities_by_dimension( covering_set, 2, cintxelems );MB_CHK_ERR( rval );
609  velist[4] = cintxverts.size();
610  velist[5] = cintxelems.size();
611  }
612 
613  MPI_Reduce( velist, gvelist, 6, MPI_UINT64_T, MPI_SUM, 0, MPI_COMM_WORLD );
614 
615 #else
616  moab::EntityHandle covering_set = runCtx->meshsets[0];
617  for( int i = 0; i < 6; i++ )
618  gvelist[i] = velist[i];
619 #endif
620 
621  if( !proc_id )
622  {
623  outputFormatter.printf( 0, "The source set contains %lu vertices and %lu elements \n", gvelist[0],
624  gvelist[0] );
625  outputFormatter.printf( 0, "The covering set contains %lu vertices and %lu elements \n", gvelist[2],
626  gvelist[2] );
627  outputFormatter.printf( 0, "The target set contains %lu vertices and %lu elements \n", gvelist[1],
628  gvelist[1] );
629  }
630 
631  // Now let's invoke the MOAB intersection algorithm in parallel with a
632  // source and target mesh set representing two different decompositions
633  runCtx->timer_push( "compute intersections with MOAB" );
634  rval = mbCore->create_meshset( moab::MESHSET_SET, intxset );MB_CHK_SET_ERR( rval, "Can't create new set" );
635  rval = mbintx->intersect_meshes( covering_set, runCtx->meshsets[1], intxset );MB_CHK_SET_ERR( rval, "Can't compute the intersection of meshes on the sphere" );
636  runCtx->timer_pop();
637 
638  // free the memory
639  delete mbintx;
640  }
641 
642  {
643  moab::Range intxelems, intxverts;
644  rval = mbCore->get_entities_by_dimension( intxset, 2, intxelems );MB_CHK_ERR( rval );
645  rval = mbCore->get_entities_by_dimension( intxset, 0, intxverts, true );MB_CHK_ERR( rval );
646  outputFormatter.printf( 0, "The intersection set contains %lu elements and %lu vertices \n",
647  intxelems.size(), intxverts.size() );
648 
649  moab::IntxAreaUtils areaAdaptorHuiller( moab::IntxAreaUtils::lHuiller ); // lHuiller, GaussQuadrature
650  double initial_sarea =
651  areaAdaptorHuiller.area_on_sphere( mbCore, runCtx->meshsets[0],
652  radius_src ); // use the target to compute the initial area
653  double initial_tarea =
654  areaAdaptorHuiller.area_on_sphere( mbCore, runCtx->meshsets[1],
655  radius_dest ); // use the target to compute the initial area
656  double intx_area = areaAdaptorHuiller.area_on_sphere( mbCore, intxset, radius_src );
657 
658  outputFormatter.printf( 0, "mesh areas: source = %12.10f, target = %12.10f, intersection = %12.10f \n",
659  initial_sarea, initial_tarea, intx_area );
660  outputFormatter.printf( 0, "relative error w.r.t source = %12.10e, target = %12.10e \n",
661  fabs( intx_area - initial_sarea ) / initial_sarea,
662  fabs( intx_area - initial_tarea ) / initial_tarea );
663  }
664 
665  // Write out our computed intersection file
666  if( !runCtx->skip_io )
667  {
668  rval = mbCore->write_mesh( "moab_intersection.h5m", &intxset, 1 );MB_CHK_ERR( rval );
669  }
670 
671  if( runCtx->computeWeights )
672  {
673  runCtx->timer_push( "compute weights with the Tempest meshes" );
674  // Call to generate an offline map with the tempest meshes
675  OfflineMap weightMap;
676  int err = GenerateOfflineMapWithMeshes( *runCtx->meshes[0], *runCtx->meshes[1], *runCtx->meshes[2],
677  runCtx->disc_methods[0], // std::string strInputType
678  runCtx->disc_methods[1], // std::string strOutputType,
679  runCtx->mapOptions, weightMap );
680  runCtx->timer_pop();
681 
682  std::map< std::string, std::string > mapAttributes;
683  if( err )
684  {
685  rval = moab::MB_FAILURE;
686  }
687  else
688  {
689  if( !runCtx->skip_io ) weightMap.Write( "outWeights.nc", mapAttributes );
690  }
691  }
692  }
693  else if( runCtx->meshType == moab::TempestRemapper::OVERLAP_MOAB )
694  {
695  // Usage: mpiexec -n 2 tools/mbtempest -t 5 -l mycs_2.h5m -l myico_2.h5m -f myoverlap_2.h5m
696 #ifdef MOAB_HAVE_MPI
697  rval = pcomm->check_all_shared_handles();MB_CHK_ERR( rval );
698 #endif
699 
700  // print verbosely about the problem setting
701  size_t velist[4] = {}, gvelist[4] = {};
702  {
703  moab::Range srcverts, srcelems;
704  rval = mbCore->get_entities_by_dimension( runCtx->meshsets[0], 0, srcverts );MB_CHK_ERR( rval );
705  rval = mbCore->get_entities_by_dimension( runCtx->meshsets[0], 2, srcelems );MB_CHK_ERR( rval );
706  rval = moab::IntxUtils::fix_degenerate_quads( mbCore, runCtx->meshsets[0] );MB_CHK_ERR( rval );
707  if( runCtx->enforceConvexity )
708  {
709  rval = moab::IntxUtils::enforce_convexity( mbCore, runCtx->meshsets[0], proc_id );MB_CHK_ERR( rval );
710  }
711  rval = areaAdaptor.positive_orientation( mbCore, runCtx->meshsets[0], radius_src );MB_CHK_ERR( rval );
712  // if( !proc_id )
713  // outputFormatter.printf( 0, "The source set contains %lu vertices and %lu elements \n", srcverts.size(),
714  // srcelems.size() );
715  velist[0] = srcverts.size();
716  velist[1] = srcelems.size();
717 
718  moab::Range tgtverts, tgtelems;
719  rval = mbCore->get_entities_by_dimension( runCtx->meshsets[1], 0, tgtverts );MB_CHK_ERR( rval );
720  rval = mbCore->get_entities_by_dimension( runCtx->meshsets[1], 2, tgtelems );MB_CHK_ERR( rval );
721  rval = moab::IntxUtils::fix_degenerate_quads( mbCore, runCtx->meshsets[1] );MB_CHK_ERR( rval );
722  if( runCtx->enforceConvexity )
723  {
724  rval = moab::IntxUtils::enforce_convexity( mbCore, runCtx->meshsets[1], proc_id );MB_CHK_ERR( rval );
725  }
726  rval = areaAdaptor.positive_orientation( mbCore, runCtx->meshsets[1], radius_dest );MB_CHK_ERR( rval );
727  // if( !proc_id )
728  // outputFormatter.printf( 0, "The target set contains %lu vertices and %lu elements \n", tgtverts.size(),
729  // tgtelems.size() );
730  velist[2] = tgtverts.size();
731  velist[3] = tgtelems.size();
732  }
733  //rval = mbCore->write_file( "source_mesh.h5m", nullptr, writeOptions, &runCtx->meshsets[0], 1 );MB_CHK_ERR( rval );
734  //rval = mbCore->write_file( "target_mesh.h5m", nullptr, writeOptions, &runCtx->meshsets[1], 1 );MB_CHK_ERR( rval );
735 
736  // if( runCtx->nlayers && nprocs > 1 )
737  // {
738  // remapper.ResetMeshSet( moab::Remapper::SourceMesh, runCtx->meshsets[3] );
739  // runCtx->meshes[0] = remapper.GetMesh( moab::Remapper::SourceMesh ); // ?
740  // }
741 
742  // First compute the covering set such that the target elements are fully covered by the
743  // local source grid
744  runCtx->timer_push( "construct covering set for intersection" );
745  rval = remapper.ConstructCoveringSet( runCtx->epsrel, 1.0, 1.0, runCtx->boxeps, runCtx->rrmGrids,
746  runCtx->useGnomonicProjection, runCtx->nlayers );MB_CHK_ERR( rval );
747  runCtx->timer_pop();
748 
749 #ifdef MOAB_HAVE_MPI
750  MPI_Reduce( velist, gvelist, 4, MPI_UINT64_T, MPI_SUM, 0, MPI_COMM_WORLD );
751 #else
752  for( int i = 0; i < 4; i++ )
753  gvelist[i] = velist[i];
754 #endif
755 
756  if( !proc_id && runCtx->print_diagnostics )
757  {
758  outputFormatter.printf( 0, "The source set contains %lu vertices and %lu elements \n", gvelist[0],
759  gvelist[1] );
760  outputFormatter.printf( 0, "The target set contains %lu vertices and %lu elements \n", gvelist[2],
761  gvelist[3] );
762  }
763 
764  // Compute intersections with MOAB with either the Kd-tree or the advancing front algorithm
765  runCtx->timer_push( "setup and compute mesh intersections" );
766  rval = remapper.ComputeOverlapMesh( runCtx->kdtreeSearch, false, runCtx->nlayers );MB_CHK_ERR( rval );
767  runCtx->timer_pop();
768 
769  // print some diagnostic checks to see if the overlap grid resolved the input meshes
770  // correctly
771  double dTotalOverlapArea = 0.0;
772  if( runCtx->print_diagnostics )
773  {
774  moab::IntxAreaUtils areaAdaptorHuiller(
775  moab::IntxAreaUtils::GaussQuadrature ); // lHuiller, GaussQuadrature
776  double local_areas[3],
777  global_areas[3]; // Array for Initial area, and through Method 1 and Method 2
778  // local_areas[0] = area_on_sphere_lHuiller ( mbCore, runCtx->meshsets[1], radius_src );
779  local_areas[0] = areaAdaptorHuiller.area_on_sphere( mbCore, runCtx->meshsets[0], radius_src );
780  local_areas[1] = areaAdaptorHuiller.area_on_sphere( mbCore, runCtx->meshsets[1], radius_dest );
781  local_areas[2] = areaAdaptorHuiller.area_on_sphere( mbCore, runCtx->meshsets[2], radius_src );
782 
783 #ifdef MOAB_HAVE_MPI
784  MPI_Allreduce( &local_areas[0], &global_areas[0], 3, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD );
785 #else
786  global_areas[0] = local_areas[0];
787  global_areas[1] = local_areas[1];
788  global_areas[2] = local_areas[2];
789 #endif
790  if( !proc_id )
791  {
792  outputFormatter.printf( 0,
793  "initial area: source mesh = %12.14f, target mesh = "
794  "%12.14f, overlap mesh = %12.14f\n",
795  global_areas[0], global_areas[1], global_areas[2] );
796  // outputFormatter.printf ( 0, " area with l'Huiller: %12.14f with Girard:
797  // %12.14f\n", global_areas[2], global_areas[3] ); outputFormatter.printf ( 0, "
798  // relative difference areas = %12.10e\n", fabs ( global_areas[2] - global_areas[3]
799  // ) / global_areas[2] );
800  outputFormatter.printf( 0, "relative error w.r.t source = %12.14e, and target = %12.14e\n",
801  fabs( global_areas[0] - global_areas[2] ) / global_areas[0],
802  fabs( global_areas[1] - global_areas[2] ) / global_areas[1] );
803  }
804  dTotalOverlapArea = global_areas[2];
805  }
806 
807  if( runCtx->intxFilename.size() )
808  {
809  moab::EntityHandle writableOverlapSet;
810  rval = mbCore->create_meshset( moab::MESHSET_SET, writableOverlapSet );MB_CHK_SET_ERR( rval, "Can't create new set" );
811  moab::EntityHandle meshOverlapSet = remapper.GetMeshSet( moab::Remapper::OverlapMesh );
812  moab::Range ovEnts;
813  rval = mbCore->get_entities_by_dimension( meshOverlapSet, 2, ovEnts );MB_CHK_SET_ERR( rval, "Can't create new set" );
814  rval = mbCore->get_entities_by_dimension( meshOverlapSet, 0, ovEnts );MB_CHK_SET_ERR( rval, "Can't create new set" );
815 
816 #ifdef MOAB_HAVE_MPI
817  // Do not remove ghosted entities if we still haven't computed weights
818  // Remove ghosted entities from overlap set before writing the new mesh set to file
819  if( nprocs > 1 )
820  {
821  moab::Range ghostedEnts;
822  rval = remapper.GetOverlapAugmentedEntities( ghostedEnts );MB_CHK_ERR( rval );
823  ovEnts = moab::subtract( ovEnts, ghostedEnts );
824 #ifdef MOAB_DBG
825  if( !runCtx->skip_io )
826  {
827  std::stringstream filename;
828  filename << "aug_overlap" << runCtx->pcomm->rank() << ".h5m";
829  rval = runCtx->mbcore->write_file( filename.str().c_str(), 0, 0, &meshOverlapSet, 1 );MB_CHK_ERR( rval );
830  }
831 #endif
832  }
833 #endif
834  rval = mbCore->add_entities( writableOverlapSet, ovEnts );MB_CHK_SET_ERR( rval, "adding local intx cells failed" );
835 
836 #ifdef MOAB_HAVE_MPI
837 #ifdef MOAB_DBG
838  if( nprocs > 1 && !runCtx->skip_io )
839  {
840  std::stringstream filename;
841  filename << "writable_intx_" << runCtx->pcomm->rank() << ".h5m";
842  rval = runCtx->mbcore->write_file( filename.str().c_str(), 0, 0, &writableOverlapSet, 1 );MB_CHK_ERR( rval );
843  }
844 #endif
845 #endif
846 
847  size_t lastindex = runCtx->intxFilename.find_last_of( "." );
848  sstr.str( "" );
849  sstr << runCtx->intxFilename.substr( 0, lastindex ) << ".h5m";
850  if( !runCtx->proc_id )
851  std::cout << "Writing out the MOAB intersection mesh file to " << sstr.str() << std::endl;
852 
853  // Write out our computed intersection file
854  if( !runCtx->skip_io )
855  {
856  rval = mbCore->write_file( sstr.str().c_str(), nullptr, writeOptions, &writableOverlapSet, 1 );MB_CHK_ERR( rval );
857  }
858  }
859 
860  if( runCtx->computeWeights )
861  {
862  runCtx->meshes[2] = remapper.GetMesh( moab::Remapper::OverlapMesh );
863  if( !runCtx->proc_id ) std::cout << std::endl;
864 
865  runCtx->timer_push( "setup computation of weights" );
866  // Call to generate the remapping weights with the tempest meshes
867  moab::TempestOnlineMap* weightMap = new moab::TempestOnlineMap( &remapper );
868  runCtx->timer_pop();
869 
870  runCtx->timer_push( "compute weights with TempestRemap" );
871  rval = weightMap->GenerateRemappingWeights(
872  runCtx->disc_methods[0], // std::string strInputType
873  runCtx->disc_methods[1], // std::string strOutputType,
874  runCtx->mapOptions, // const GenerateOfflineMapAlgorithmOptions& options
875  runCtx->doftag_names[0], // const std::string& source_tag_name
876  runCtx->doftag_names[1] // const std::string& target_tag_name
877  );MB_CHK_ERR( rval );
878  runCtx->timer_pop();
879 
880  weightMap->PrintMapStatistics();
881 
882  // Invoke the CheckMap routine on the TempestRemap serial interface directly, if running
883  // on a single process
884  if( runCtx->fCheck )
885  {
886  const double dNormalTolerance = 1.0E-8;
887  const double dStrictTolerance = 1.0E-12;
888  weightMap->CheckMap( runCtx->fCheck, runCtx->fCheck, runCtx->fCheck && ( runCtx->ensureMonotonicity ),
889  dNormalTolerance, dStrictTolerance, dTotalOverlapArea );
890  }
891 
892  if( runCtx->outFilename.size() )
893  {
894  std::map< std::string, std::string > attrMap;
895  attrMap["MOABversion"] = std::string( MOAB_VERSION );
896  attrMap["Title"] = "MOAB-TempestRemap (mbtempest) Offline Regridding Weight Generator";
897  attrMap["normalization"] = "ovarea";
898  attrMap["remap_options"] = runCtx->mapOptions.strMethod;
899  attrMap["domain_a"] = runCtx->inFilenames[0];
900  attrMap["domain_b"] = runCtx->inFilenames[1];
901  if( runCtx->intxFilename.size() ) attrMap["domain_aUb"] = runCtx->intxFilename;
902  attrMap["map_aPb"] = runCtx->outFilename;
903  attrMap["methodorder_a"] = runCtx->disc_methods[0] + ":" + std::to_string( runCtx->disc_orders[0] ) +
904  ":" + std::string( runCtx->doftag_names[0] );
905  attrMap["concave_a"] = runCtx->mapOptions.fSourceConcave ? "true" : "false";
906  attrMap["methodorder_b"] = runCtx->disc_methods[1] + ":" + std::to_string( runCtx->disc_orders[1] ) +
907  ":" + std::string( runCtx->doftag_names[1] );
908  attrMap["concave_b"] = runCtx->mapOptions.fTargetConcave ? "true" : "false";
909  attrMap["bubble"] = runCtx->mapOptions.fNoBubble ? "false" : "true";
910  attrMap["history"] = historyStr;
911 
912  // Write the map file to disk in parallel using either HDF5 or SCRIP interface
913  // in extra case; maybe need a better solution, just create it with the right meshset
914  // from the beginning;
915  if( !runCtx->skip_io )
916  {
917  rval = weightMap->WriteParallelMap( runCtx->outFilename.c_str(), attrMap );MB_CHK_ERR( rval );
918  }
919  }
920 
921  if( runCtx->verifyWeights )
922  {
923  // Let us pick a sampling test function for solution evaluation
924  // SH, SV, FH, USERVAR
925  bool userVariable = false;
926  moab::TempestOnlineMap::sample_function testFunction;
927  if( !runCtx->variableToVerify.compare( "SH" ) )
928  testFunction = &sample_slow_harmonic;
929  else if( !runCtx->variableToVerify.compare( "FH" ) )
930  testFunction = &sample_fast_harmonic;
931  else if( !runCtx->variableToVerify.compare( "SV" ) )
932  testFunction = &sample_stationary_vortex;
933  else
934  {
935  userVariable = runCtx->variableToVerify.size() ? true : false;
936  testFunction = runCtx->variableToVerify.size() ? nullptr : sample_stationary_vortex;
937  }
938 
939  moab::Tag srcAnalyticalFunction;
940  moab::Tag tgtAnalyticalFunction;
941  moab::Tag tgtProjectedFunction;
942  if( testFunction )
943  {
944  runCtx->timer_push( "describe a solution on source grid" );
945  rval = weightMap->DefineAnalyticalSolution( srcAnalyticalFunction, "AnalyticalSolnSrcExact",
946  moab::Remapper::SourceMesh, testFunction );MB_CHK_ERR( rval );
947  runCtx->timer_pop();
948 
949  runCtx->timer_push( "describe a solution on target grid" );
950 
951  rval = weightMap->DefineAnalyticalSolution( tgtAnalyticalFunction, "AnalyticalSolnTgtExact",
952  moab::Remapper::TargetMesh, testFunction,
953  &tgtProjectedFunction, "ProjectedSolnTgt" );MB_CHK_ERR( rval );
954  // rval = mbCore->write_file ( "tgtWithSolnTag.h5m", nullptr, writeOptions,
955  // &runCtx->meshsets[1], 1 ); MB_CHK_ERR ( rval );
956  runCtx->timer_pop();
957  }
958  else
959  {
960  rval = mbCore->tag_get_handle( runCtx->variableToVerify.c_str(), srcAnalyticalFunction );MB_CHK_ERR( rval );
961 
962  rval = mbCore->tag_get_handle( "ProjectedSolnTgt", 1, moab::MB_TYPE_DOUBLE, tgtProjectedFunction,
963  moab::MB_TAG_DENSE | moab::MB_TAG_CREAT );MB_CHK_ERR( rval );
964  }
965 
966  if( !runCtx->skip_io )
967  {
968  rval = mbCore->write_file( "srcWithSolnTag.h5m", nullptr, writeOptions, &runCtx->meshsets[0], 1 );MB_CHK_ERR( rval );
969  }
970 
971  runCtx->timer_push( "compute solution projection on target grid" );
972  rval = weightMap->ApplyWeights( srcAnalyticalFunction, tgtProjectedFunction, false, runCtx->cassType );MB_CHK_ERR( rval );
973  runCtx->timer_pop();
974 
975  if( !runCtx->skip_io )
976  {
977  rval = mbCore->write_file( "tgtWithSolnTag2.h5m", nullptr, writeOptions, &runCtx->meshsets[1], 1 );MB_CHK_ERR( rval );
978  }
979 
980  if( nprocs == 1 && runCtx->baselineFile.size() )
981  {
982  // save the field from tgtWithSolnTag2 in a text file, and global ids for cells
983  moab::Range tgtCells;
984  rval = mbCore->get_entities_by_dimension( runCtx->meshsets[1], 2, tgtCells );MB_CHK_ERR( rval );
985  std::vector< int > globIds;
986  globIds.resize( tgtCells.size() );
987  std::vector< double > vals;
988  vals.resize( tgtCells.size() );
989  moab::Tag projTag;
990  rval = mbCore->tag_get_handle( "ProjectedSolnTgt", projTag );MB_CHK_ERR( rval );
991  moab::Tag gid = mbCore->globalId_tag();
992  rval = mbCore->tag_get_data( gid, tgtCells, &globIds[0] );MB_CHK_ERR( rval );
993  rval = mbCore->tag_get_data( projTag, tgtCells, &vals[0] );MB_CHK_ERR( rval );
994  std::fstream fs;
995  fs.open( runCtx->baselineFile.c_str(), std::fstream::out );
996  fs << std::setprecision( 15 ); // maximum precision for doubles
997  for( size_t i = 0; i < tgtCells.size(); i++ )
998  fs << globIds[i] << " " << vals[i] << "\n";
999  fs.close();
1000  // for good measure, save the source file too, with the tag AnalyticalSolnSrcExact
1001  // it will be used later to test, along with a target file
1002  if( !runCtx->skip_io )
1003  {
1004  rval = mbCore->write_file( "srcWithSolnTag.h5m", nullptr, writeOptions, &runCtx->meshsets[0], 1 );MB_CHK_ERR( rval );
1005  }
1006  }
1007 
1008  // compute error metrics if it is a known analytical functional
1009  if( !userVariable )
1010  {
1011  runCtx->timer_push( "compute error metrics against analytical solution on target grid" );
1012  std::map< std::string, double > errMetrics;
1013  rval = weightMap->ComputeMetrics( moab::Remapper::TargetMesh, tgtAnalyticalFunction,
1014  tgtProjectedFunction, errMetrics, true );MB_CHK_ERR( rval );
1015  runCtx->timer_pop();
1016  }
1017  }
1018 
1019  delete weightMap;
1020  }
1021  }
1022 
1023  // Clean up
1024  remapper.clear();
1025  delete runCtx;
1026  delete mbCore;
1027 
1028 #ifdef MOAB_HAVE_MPI
1029  MPI_Finalize();
1030 #endif
1031  exit( 0 );
1032 }
1033 
1034 static moab::ErrorCode CreateTempestMesh( ToolContext& ctx, moab::TempestRemapper& remapper, Mesh* tempest_mesh )
1035 {
1036  moab::ErrorCode rval = moab::MB_SUCCESS;
1037  int err;
1038  moab::DebugOutput& outputFormatter = ctx.outputFormatter;
1039 
1040  if( !ctx.proc_id )
1041  {
1042  outputFormatter.printf( 0, "Creating TempestRemap Mesh object ...\n" );
1043  }
1044 
1045  if( ctx.meshType == moab::TempestRemapper::OVERLAP_FILES )
1046  {
1047  ctx.timer_push( "create Tempest OverlapMesh" );
1048  // For the overlap method, choose between: "fuzzy", "exact" or "mixed"
1049  err = GenerateOverlapMesh( ctx.inFilenames[0], ctx.inFilenames[1], *tempest_mesh, ctx.outFilename, "NetCDF4",
1050  "exact", true );
1051  if( err )
1052  rval = moab::MB_FAILURE;
1053  else
1054  ctx.meshes.push_back( tempest_mesh );
1055  ctx.timer_pop();
1056  }
1057  else if( ctx.meshType == moab::TempestRemapper::OVERLAP_MEMORY )
1058  {
1059  // Load the meshes and validate
1060  ctx.meshsets.resize( 3 );
1061  ctx.meshes.resize( 3 );
1062  ctx.meshsets[0] = remapper.GetMeshSet( moab::Remapper::SourceMesh );
1063  ctx.meshsets[1] = remapper.GetMeshSet( moab::Remapper::TargetMesh );
1064  ctx.meshsets[2] = remapper.GetMeshSet( moab::Remapper::OverlapMesh );
1065 
1066  ctx.timer_push( "load MOAB Source mesh" );
1067  // First the source
1068  rval = remapper.LoadMesh( moab::Remapper::SourceMesh, ctx.inFilenames[0], moab::TempestRemapper::DEFAULT );MB_CHK_ERR( rval );
1069  ctx.meshes[0] = remapper.GetMesh( moab::Remapper::SourceMesh );
1070  ctx.timer_pop();
1071 
1072  ctx.timer_push( "load MOAB Target mesh" );
1073  // Next the target
1074  rval = remapper.LoadMesh( moab::Remapper::TargetMesh, ctx.inFilenames[1], moab::TempestRemapper::DEFAULT );MB_CHK_ERR( rval );
1075  ctx.meshes[1] = remapper.GetMesh( moab::Remapper::TargetMesh );
1076  ctx.timer_pop();
1077 
1078  ctx.timer_push( "generate TempestRemap OverlapMesh" );
1079  // Now let us construct the overlap mesh, by calling TempestRemap interface directly
1080  // For the overlap method, choose between: "fuzzy", "exact" or "mixed"
1081  err = GenerateOverlapWithMeshes( *ctx.meshes[0], *ctx.meshes[1], *tempest_mesh, "" /*ctx.outFilename*/,
1082  "NetCDF4", "exact", false );
1083  ctx.timer_pop();
1084  if( err )
1085  rval = moab::MB_FAILURE;
1086  else
1087  {
1088  remapper.SetMesh( moab::Remapper::OverlapMesh, tempest_mesh );
1089  ctx.meshes[2] = remapper.GetMesh( moab::Remapper::OverlapMesh );
1090  // ctx.meshes.push_back(*tempest_mesh);
1091  }
1092  }
1093  else if( ctx.meshType == moab::TempestRemapper::OVERLAP_MOAB )
1094  {
1095  ctx.meshsets.resize( 3 );
1096  ctx.meshes.resize( 3 );
1097  ctx.meshsets[0] = remapper.GetMeshSet( moab::Remapper::SourceMesh );
1098  ctx.meshsets[1] = remapper.GetMeshSet( moab::Remapper::TargetMesh );
1099  ctx.meshsets[2] = remapper.GetMeshSet( moab::Remapper::OverlapMesh );
1100 
1101  const double radius_src = 1.0 /*2.0*acos(-1.0)*/;
1102  const double radius_dest = 1.0 /*2.0*acos(-1.0)*/;
1103 
1104  std::vector< int > smetadata, tmetadata;
1105  //const char* additional_read_opts = ( ctx.n_procs > 1 ? "NO_SET_CONTAINING_PARENTS;" : "" );
1106  std::string additional_read_opts_src = get_file_read_options( ctx, ctx.inFilenames[0] );
1107 
1108  ctx.timer_push( "load MOAB Source mesh" );
1109  // Load the source mesh and validate
1110  rval =
1111  remapper.LoadNativeMesh( ctx.inFilenames[0], ctx.meshsets[0], smetadata, additional_read_opts_src.c_str() );MB_CHK_ERR( rval );
1112  if( smetadata.size() ) remapper.SetMeshType( moab::Remapper::SourceMesh, smetadata );
1113  ctx.timer_pop();
1114 
1115  ctx.timer_push( "preprocess MOAB Source mesh" );
1116  // Rescale the radius of both to compute the intersection
1117  rval = moab::IntxUtils::ScaleToRadius( ctx.mbcore, ctx.meshsets[0], radius_src );MB_CHK_ERR( rval );
1118  // if order >=2, ghost at least this many layers (order -1) it may be too much
1119 #ifdef MOAB_HAVE_MPI
1120  if( ctx.nlayers && ctx.n_procs > 1 )
1121  {
1122  // get order -1 ghost layers; actually it should be decided by the mesh
1123  // if the mesh has holes, it could be more
1124 
1125  moab::EntityHandle set_with_ghosts;
1126  rval = remapper.GhostLayers( ctx.meshsets[0], ctx.nlayers, set_with_ghosts );MB_CHK_ERR( rval );
1127  remapper.SetMeshSet( moab::Remapper::SourceMeshWithGhosts, set_with_ghosts );
1128 #ifdef MOAB_DBG
1129  if( !runCtx->skip_io )
1130  {
1131  // write the new source sets, after layers were decided, should see the ghosts now
1132  std::stringstream filename;
1133  filename << "expand_source" << ctx.pcomm->rank() << ".h5m";
1134  rval = ctx.mbcore->write_file( filename.str().c_str(), 0, 0, &( ctx.meshsets[0] ), 1 );MB_CHK_ERR( rval );
1135  }
1136 #endif
1137  }
1138 #endif
1139  ctx.timer_pop();
1140 
1141  // Load the target mesh and validate
1142  std::string addititional_read_opts_tgt = get_file_read_options( ctx, ctx.inFilenames[1] );
1143  // addititional_read_opts_tgt += "PARALLEL_GHOSTS=2.0.3;PARALLEL_THIN_GHOST_LAYER;";
1144 
1145  ctx.timer_push( "load MOAB Target mesh" );
1146  rval = remapper.LoadNativeMesh( ctx.inFilenames[1], ctx.meshsets[1], tmetadata,
1147  addititional_read_opts_tgt.c_str() );MB_CHK_ERR( rval );
1148  if( tmetadata.size() ) remapper.SetMeshType( moab::Remapper::TargetMesh, tmetadata );
1149  ctx.timer_pop();
1150 
1151  ctx.timer_push( "preprocess MOAB Target mesh" );
1152  rval = moab::IntxUtils::ScaleToRadius( ctx.mbcore, ctx.meshsets[1], radius_dest );MB_CHK_ERR( rval );
1153  ctx.timer_pop();
1154 
1155  if( ctx.computeWeights )
1156  {
1157  ctx.timer_push( "convert MOAB meshes to TempestRemap meshes in memory" );
1158  // convert MOAB representation to TempestRemap's Mesh for source
1159  rval = remapper.ConvertMeshToTempest( moab::Remapper::SourceMesh );MB_CHK_ERR( rval );
1160  ctx.meshes[0] = remapper.GetMesh( moab::Remapper::SourceMesh );
1161 
1162  // convert MOAB representation to TempestRemap's Mesh for target
1163  rval = remapper.ConvertMeshToTempest( moab::Remapper::TargetMesh );MB_CHK_ERR( rval );
1164  ctx.meshes[1] = remapper.GetMesh( moab::Remapper::TargetMesh );
1165  ctx.timer_pop();
1166  }
1167  }
1168  else if( ctx.meshType == moab::TempestRemapper::ICO )
1169  {
1170  ctx.timer_push( "generate ICO mesh with TempestRemap" );
1171  err = GenerateICOMesh( *tempest_mesh, ctx.blockSize, ctx.computeDual, ctx.outFilename, "NetCDF4" );
1172  ctx.timer_pop();
1173 
1174  if( err )
1175  rval = moab::MB_FAILURE;
1176  else
1177  ctx.meshes.push_back( tempest_mesh );
1178  }
1179  else if( ctx.meshType == moab::TempestRemapper::RLL )
1180  {
1181  ctx.timer_push( "generate RLL mesh with TempestRemap" );
1182  err = GenerateRLLMesh( *tempest_mesh, // Mesh& meshOut,
1183  ctx.blockSize * 2, ctx.blockSize, // int nLongitudes, int nLatitudes,
1184  0.0, 360.0, // double dLonBegin, double dLonEnd,
1185  -90.0, 90.0, // double dLatBegin, double dLatEnd,
1186  false, false, false, // bool fGlobalCap, bool fFlipLatLon, bool fForceGlobal,
1187  "" /*ctx.inFilename*/, "",
1188  "", // std::string strInputFile, std::string strInputFileLonName, std::string
1189  // strInputFileLatName,
1190  ctx.outFilename, "NetCDF4", // std::string strOutputFile, std::string strOutputFormat
1191  true // bool fVerbose
1192  );
1193  ctx.timer_pop();
1194 
1195  if( err )
1196  rval = moab::MB_FAILURE;
1197  else
1198  ctx.meshes.push_back( tempest_mesh );
1199  }
1200  else // default
1201  {
1202  ctx.timer_push( "generate CS mesh with TempestRemap" );
1203  err = GenerateCSMesh( *tempest_mesh, ctx.blockSize, ctx.outFilename, "NetCDF4" );
1204  ctx.timer_pop();
1205  if( err )
1206  rval = moab::MB_FAILURE;
1207  else
1208  ctx.meshes.push_back( tempest_mesh );
1209  }
1210 
1211  if( ctx.meshType != moab::TempestRemapper::OVERLAP_MOAB && !tempest_mesh )
1212  {
1213  std::cout << "Tempest Mesh is not a complete object; Quitting...";
1214  exit( -1 );
1215  }
1216 
1217  return rval;
1218 }
1219 
1220 #undef MOAB_DBG
1221 ///////////////////////////////////////////////
1222 // Test functions
1223 
1224 double sample_slow_harmonic( double dLon, double dLat )
1225 {
1226  return ( 2.0 + cos( dLat ) * cos( dLat ) * cos( 2.0 * dLon ) );
1227 }
1228 
1229 double sample_fast_harmonic( double dLon, double dLat )
1230 {
1231  return ( 2.0 + pow( sin( 2.0 * dLat ), 16.0 ) * cos( 16.0 * dLon ) );
1232  // return (2.0 + pow(cos(2.0 * dLat), 16.0) * cos(16.0 * dLon));
1233 }
1234 
1235 double sample_constant( double /*dLon*/, double /*dLat*/ )
1236 {
1237  return 1.0;
1238 }
1239 
1240 double sample_stationary_vortex( double dLon, double dLat )
1241 {
1242  const double dLon0 = 0.0;
1243  const double dLat0 = 0.6;
1244  const double dR0 = 3.0;
1245  const double dD = 5.0;
1246  const double dT = 6.0;
1247 
1248  /// Find the rotated longitude and latitude of a point on a sphere
1249  /// with pole at (dLonC, dLatC).
1250  {
1251  double dSinC = sin( dLat0 );
1252  double dCosC = cos( dLat0 );
1253  double dCosT = cos( dLat );
1254  double dSinT = sin( dLat );
1255 
1256  double dTrm = dCosT * cos( dLon - dLon0 );
1257  double dX = dSinC * dTrm - dCosC * dSinT;
1258  double dY = dCosT * sin( dLon - dLon0 );
1259  double dZ = dSinC * dSinT + dCosC * dTrm;
1260 
1261  dLon = atan2( dY, dX );
1262  if( dLon < 0.0 )
1263  {
1264  dLon += 2.0 * M_PI;
1265  }
1266  dLat = asin( dZ );
1267  }
1268 
1269  double dRho = dR0 * cos( dLat );
1270  double dVt = 3.0 * sqrt( 3.0 ) / 2.0 / cosh( dRho ) / cosh( dRho ) * tanh( dRho );
1271 
1272  double dOmega;
1273  if( dRho == 0.0 )
1274  {
1275  dOmega = 0.0;
1276  }
1277  else
1278  {
1279  dOmega = dVt / dRho;
1280  }
1281 
1282  return ( 1.0 - tanh( dRho / dD * sin( dLon - dOmega * dT ) ) );
1283 }
1284 
1285 ///////////////////////////////////////////////