Mesh Oriented datABase  (version 5.5.1)
An array-based unstructured mesh library
Classes | Functions
mbtempest.cpp File Reference
#include <iostream>
#include <iomanip>
#include <cstdlib>
#include <vector>
#include <string>
#include <sstream>
#include <cassert>
#include "moab/Core.hpp"
#include "moab/IntxMesh/IntxUtils.hpp"
#include "moab/Remapping/TempestRemapper.hpp"
#include "moab/Remapping/TempestOnlineMap.hpp"
#include "moab/ProgOptions.hpp"
#include "moab/CpuTimer.hpp"
#include "DebugOutput.hpp"
+ Include dependency graph for mbtempest.cpp:

Go to the source code of this file.

Classes

struct  ToolContext
 

Functions

static moab::ErrorCode CreateTempestMesh (ToolContext &, moab::TempestRemapper &remapper, Mesh *)
 
double sample_slow_harmonic (double dLon, double dLat)
 
double sample_fast_harmonic (double dLon, double dLat)
 
double sample_constant (double dLon, double dLat)
 
double sample_stationary_vortex (double dLon, double dLat)
 
std::string get_file_read_options (ToolContext &ctx, std::string filename)
 
int main (int argc, char *argv[])
 

Function Documentation

◆ CreateTempestMesh()

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

Definition at line 1034 of file mbtempest.cpp.

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 }

References ToolContext::blockSize, ToolContext::computeDual, ToolContext::computeWeights, moab::TempestRemapper::ConvertMeshToTempest(), moab::TempestRemapper::DEFAULT, ErrorCode, get_file_read_options(), moab::TempestRemapper::GetMesh(), moab::TempestRemapper::GetMeshSet(), moab::TempestRemapper::ICO, ToolContext::inFilenames, moab::TempestRemapper::LoadMesh(), moab::Remapper::LoadNativeMesh(), MB_CHK_ERR, MB_SUCCESS, ToolContext::mbcore, ToolContext::meshes, ToolContext::meshsets, ToolContext::meshType, ToolContext::n_procs, ToolContext::nlayers, ToolContext::outFilename, ToolContext::outputFormatter, moab::TempestRemapper::OVERLAP_FILES, moab::TempestRemapper::OVERLAP_MEMORY, moab::TempestRemapper::OVERLAP_MOAB, moab::Remapper::OverlapMesh, moab::DebugOutput::printf(), ToolContext::proc_id, moab::TempestRemapper::RLL, moab::IntxUtils::ScaleToRadius(), moab::TempestRemapper::SetMesh(), moab::TempestRemapper::SetMeshSet(), moab::TempestRemapper::SetMeshType(), moab::Remapper::SourceMesh, moab::Remapper::SourceMeshWithGhosts, moab::Remapper::TargetMesh, ToolContext::timer_pop(), ToolContext::timer_push(), and moab::Core::write_file().

Referenced by main().

◆ get_file_read_options()

std::string get_file_read_options ( ToolContext &  ctx,
std::string  filename 
)

Definition at line 417 of file mbtempest.cpp.

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 }

References ToolContext::n_procs, ncFile, and ToolContext::proc_id.

Referenced by CreateTempestMesh().

◆ main()

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

Definition at line 484 of file mbtempest.cpp.

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 }

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

◆ sample_constant()

double sample_constant ( double  dLon,
double  dLat 
)
inline

Definition at line 1235 of file mbtempest.cpp.

1236 {
1237  return 1.0;
1238 }

◆ sample_fast_harmonic()

double sample_fast_harmonic ( double  dLon,
double  dLat 
)
inline

Definition at line 1229 of file mbtempest.cpp.

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 }

Referenced by main().

◆ sample_slow_harmonic()

double sample_slow_harmonic ( double  dLon,
double  dLat 
)
inline

Definition at line 1224 of file mbtempest.cpp.

1225 {
1226  return ( 2.0 + cos( dLat ) * cos( dLat ) * cos( 2.0 * dLon ) );
1227 }

Referenced by main().

◆ sample_stationary_vortex()

double sample_stationary_vortex ( double  dLon,
double  dLat 
)
inline

Find the rotated longitude and latitude of a point on a sphere with pole at (dLonC, dLatC).

Definition at line 1240 of file mbtempest.cpp.

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 }

Referenced by main().