MOAB: Mesh Oriented datABase  (version 5.2.1)
mbtempest.cpp
Go to the documentation of this file.
00001 /*
00002  * Usage: MOAB-Tempest tool
00003  *
00004  * Generate a Cubed-Sphere mesh: ./mbtempest -t 0 -res 25 -f cubed_sphere_mesh.exo
00005  * Generate a RLL mesh: ./mbtempest -t 1 -res 25 -f rll_mesh.exo
00006  * Generate a Icosahedral-Sphere mesh: ./mbtempest -t 2 -res 25 <-dual> -f icosa_mesh.exo
00007  *
00008  * Now you can compute the intersections between the meshes too!
00009  *
00010  * Generate the overlap mesh: ./mbtempest -t 3 -l cubed_sphere_mesh.exo -l rll_mesh.exo -f
00011  * overlap_mesh.exo
00012  *
00013  */
00014 
00015 #include <iostream>
00016 #include <cstdlib>
00017 #include <vector>
00018 #include <string>
00019 #include <sstream>
00020 #include <cassert>
00021 
00022 #include "moab/Core.hpp"
00023 #include "moab/IntxMesh/IntxUtils.hpp"
00024 #include "moab/Remapping/TempestRemapper.hpp"
00025 #include "moab/Remapping/TempestOnlineMap.hpp"
00026 #include "moab/ProgOptions.hpp"
00027 #include "moab/CpuTimer.hpp"
00028 #include "DebugOutput.hpp"
00029 
00030 //#ifndef MOAB_HAVE_MPI
00031 //    #error mbtempest tool requires MPI configuration
00032 //#endif
00033 
00034 #ifdef MOAB_HAVE_MPI
00035 // MPI includes
00036 #include "moab_mpi.h"
00037 #include "moab/ParallelComm.hpp"
00038 #include "MBParallelConventions.h"
00039 #endif
00040 
00041 struct ToolContext
00042 {
00043     moab::Interface* mbcore;
00044 #ifdef MOAB_HAVE_MPI
00045     moab::ParallelComm* pcomm;
00046 #endif
00047     const int proc_id, n_procs;
00048     moab::DebugOutput outputFormatter;
00049     int blockSize;
00050     std::vector< std::string > inFilenames;
00051     std::vector< Mesh* > meshes;
00052     std::vector< moab::EntityHandle > meshsets;
00053     std::vector< int > disc_orders;
00054     std::vector< std::string > disc_methods;
00055     std::vector< std::string > doftag_names;
00056     std::string outFilename;
00057     std::string intxFilename;
00058     moab::TempestRemapper::TempestMeshType meshType;
00059     bool computeDual;
00060     bool computeWeights;
00061     bool verifyWeights;
00062     bool enforceConvexity;
00063     int ensureMonotonicity;
00064     bool fNoConservation;
00065     bool fVolumetric;
00066     bool rrmGrids;
00067     bool kdtreeSearch;
00068     bool fNoBubble, fInputConcave, fOutputConcave, fNoCheck;
00069 
00070 #ifdef MOAB_HAVE_MPI
00071     ToolContext( moab::Interface* icore, moab::ParallelComm* p_pcomm )
00072         : mbcore( icore ), pcomm( p_pcomm ), proc_id( pcomm->rank() ), n_procs( pcomm->size() ),
00073           outputFormatter( std::cout, pcomm->rank(), 0 ),
00074 #else
00075     ToolContext ( moab::Interface* icore ) :
00076             mbcore(icore),
00077             proc_id ( 0 ), n_procs ( 1 ),
00078             outputFormatter ( std::cout, 0, 0 ),
00079 #endif
00080           blockSize( 5 ), outFilename( "output.exo" ), intxFilename( "" ), meshType( moab::TempestRemapper::DEFAULT ),
00081           computeDual( false ), computeWeights( false ), verifyWeights( false ), enforceConvexity(false), ensureMonotonicity( 0 ),
00082           fNoConservation( false ), fVolumetric( false ), rrmGrids( false ), kdtreeSearch( true ), fNoBubble( false ),
00083           fInputConcave( false ), fOutputConcave( false ), fNoCheck( false )
00084     {
00085         inFilenames.resize( 2 );
00086         doftag_names.resize( 2 );
00087         timer = new moab::CpuTimer();
00088 
00089         outputFormatter.set_prefix( "[MBTempest]: " );
00090     }
00091 
00092     ~ToolContext()
00093     {
00094         // for (unsigned i=0; i < meshes.size(); ++i) delete meshes[i];
00095         meshes.clear();
00096         inFilenames.clear();
00097         disc_orders.clear();
00098         disc_methods.clear();
00099         doftag_names.clear();
00100         outFilename.clear();
00101         intxFilename.clear();
00102         meshsets.clear();
00103         delete timer;
00104     }
00105 
00106     void timer_push( std::string operation )
00107     {
00108         timer_ops = timer->time_since_birth();
00109         opName    = operation;
00110     }
00111 
00112     void timer_pop()
00113     {
00114         double locElapsed = timer->time_since_birth() - timer_ops, avgElapsed = 0, maxElapsed = 0;
00115 #ifdef MOAB_HAVE_MPI
00116         MPI_Reduce( &locElapsed, &maxElapsed, 1, MPI_DOUBLE, MPI_MAX, 0, pcomm->comm() );
00117         MPI_Reduce( &locElapsed, &avgElapsed, 1, MPI_DOUBLE, MPI_SUM, 0, pcomm->comm() );
00118 #else
00119         maxElapsed = locElapsed;
00120         avgElapsed = locElapsed;
00121 #endif
00122         if( !proc_id )
00123         {
00124             avgElapsed /= n_procs;
00125             std::cout << "[LOG] Time taken to " << opName.c_str() << ": max = " << maxElapsed
00126                       << ", avg = " << avgElapsed << "\n";
00127         }
00128         // std::cout << "\n[LOG" << proc_id << "] Time taken to " << opName << " = " <<
00129         // timer->time_since_birth() - timer_ops << std::endl;
00130         opName.clear();
00131     }
00132 
00133     void ParseCLOptions( int argc, char* argv[] )
00134     {
00135         ProgOptions opts;
00136         int imeshType                  = 0;
00137         std::string expectedFName      = "output.exo";
00138         std::string expectedMethod     = "fv";
00139         std::string expectedDofTagName = "GLOBAL_ID";
00140         int expectedOrder              = 1;
00141 
00142         if( !proc_id )
00143         {
00144             std::cout << "Command line options provided to mbtempest:\n  ";
00145             for( int iarg = 0; iarg < argc; ++iarg )
00146                 std::cout << argv[iarg] << " ";
00147             std::cout << std::endl << std::endl;
00148         }
00149 
00150         opts.addOpt< int >( "type,t",
00151                             "Type of mesh (default=CS; Choose from [CS=0, RLL=1, ICO=2, OVERLAP_FILES=3, "
00152                             "OVERLAP_MEMORY=4, OVERLAP_MOAB=5])",
00153                             &imeshType );
00154         opts.addOpt< int >( "res,r", "Resolution of the mesh (default=5)", &blockSize );
00155         opts.addOpt< void >( "dual,d", "Output the dual of the mesh (generally relevant only for ICO mesh)",
00156                              &computeDual );
00157         opts.addOpt< void >( "weights,w",
00158                              "Compute and output the weights using the overlap mesh (generally "
00159                              "relevant only for OVERLAP mesh)",
00160                              &computeWeights );
00161         opts.addOpt< void >( "noconserve,c",
00162                              "Do not apply conservation to the resultant weights (relevant only "
00163                              "when computing weights)",
00164                              &fNoConservation );
00165         opts.addOpt< void >( "volumetric,v",
00166                              "Apply a volumetric projection to compute the weights (relevant only "
00167                              "when computing weights)",
00168                              &fVolumetric );
00169         opts.addOpt< void >( "rrmgrids",
00170                              "At least one of the meshes is a regionally refined grid (relevant to "
00171                              "accelerate intersection computation)",
00172                              &rrmGrids );
00173         opts.addOpt< void >( "nocheck", "Do not check the generated map for conservation and consistency", &fNoCheck );
00174         opts.addOpt< void >( "advfront,a",
00175                              "Use the advancing front intersection instead of the Kd-tree based algorithm "
00176                              "to compute mesh intersections" );
00177         opts.addOpt< void >( "verify",
00178                              "Verify the accuracy of the maps by projecting analytical functions "
00179                              "from source to target "
00180                              "grid by applying the maps",
00181                              &verifyWeights );
00182         opts.addOpt< int >( "monotonic,n", "Ensure monotonicity in the weight generation", &ensureMonotonicity );
00183         opts.addOpt< std::string >( "load,l", "Input mesh filenames (a source and target mesh)", &expectedFName );
00184         opts.addOpt< int >( "order,o", "Discretization orders for the source and target solution fields",
00185                             &expectedOrder );
00186         opts.addOpt< std::string >( "method,m", "Discretization method for the source and target solution fields",
00187                                     &expectedMethod );
00188         opts.addOpt< std::string >( "global_id,g",
00189                                     "Tag name that contains the global DoF IDs for source and target solution fields",
00190                                     &expectedDofTagName );
00191         opts.addOpt< std::string >( "file,f", "Output remapping weights filename", &outFilename );
00192         opts.addOpt< std::string >( "intx,i", "Output TempestRemap intersection mesh filename", &intxFilename );
00193 
00194         opts.addOpt< void >( "enforce_convexity",
00195                                 "check convexity of input meshes to compute mesh intersections" , &enforceConvexity);
00196 
00197         opts.parseCommandLine( argc, argv );
00198 
00199         // By default - use Kd-tree based search; if user asks for advancing front, disable Kd-tree
00200         // algorithm
00201         kdtreeSearch = opts.numOptSet( "advfront,a" ) == 0;
00202 
00203         switch( imeshType )
00204         {
00205             case 0:
00206                 meshType = moab::TempestRemapper::CS;
00207                 break;
00208 
00209             case 1:
00210                 meshType = moab::TempestRemapper::RLL;
00211                 break;
00212 
00213             case 2:
00214                 meshType = moab::TempestRemapper::ICO;
00215                 break;
00216 
00217             case 3:
00218                 meshType = moab::TempestRemapper::OVERLAP_FILES;
00219                 break;
00220 
00221             case 4:
00222                 meshType = moab::TempestRemapper::OVERLAP_MEMORY;
00223                 break;
00224 
00225             case 5:
00226                 meshType = moab::TempestRemapper::OVERLAP_MOAB;
00227                 break;
00228 
00229             default:
00230                 meshType = moab::TempestRemapper::DEFAULT;
00231                 break;
00232         }
00233 
00234         if( meshType > moab::TempestRemapper::ICO )
00235         {
00236             opts.getOptAllArgs( "load,l", inFilenames );
00237             opts.getOptAllArgs( "order,o", disc_orders );
00238             opts.getOptAllArgs( "method,m", disc_methods );
00239             opts.getOptAllArgs( "global_id,i", doftag_names );
00240 
00241             if( disc_orders.size() == 0 ) { disc_orders.resize( 2, 1 ); }
00242 
00243             if( disc_orders.size() == 1 ) { disc_orders.push_back( 1 ); }
00244 
00245             if( disc_methods.size() == 0 ) { disc_methods.resize( 2, "fv" ); }
00246 
00247             if( disc_methods.size() == 1 ) { disc_methods.push_back( "fv" ); }
00248 
00249             if( doftag_names.size() == 0 ) { doftag_names.resize( 2, "GLOBAL_ID" ); }
00250 
00251             if( doftag_names.size() == 1 ) { doftag_names.push_back( "GLOBAL_ID" ); }
00252 
00253             assert( inFilenames.size() == 2 );
00254             assert( disc_orders.size() == 2 );
00255             assert( disc_methods.size() == 2 );
00256         }
00257 
00258         expectedFName.clear();
00259     }
00260 
00261   private:
00262     moab::CpuTimer* timer;
00263     double timer_ops;
00264     std::string opName;
00265 };
00266 
00267 // Forward declare some methods
00268 static moab::ErrorCode CreateTempestMesh( ToolContext&, moab::TempestRemapper& remapper, Mesh* );
00269 inline double sample_slow_harmonic( double dLon, double dLat );
00270 inline double sample_fast_harmonic( double dLon, double dLat );
00271 inline double sample_constant( double dLon, double dLat );
00272 inline double sample_stationary_vortex( double dLon, double dLat );
00273 
00274 int main( int argc, char* argv[] )
00275 {
00276     moab::ErrorCode rval;
00277     NcError error( NcError::verbose_nonfatal );
00278     std::stringstream sstr;
00279 
00280     int proc_id = 0, nprocs = 1;
00281 #ifdef MOAB_HAVE_MPI
00282     MPI_Init( &argc, &argv );
00283     MPI_Comm_rank( MPI_COMM_WORLD, &proc_id );
00284     MPI_Comm_size( MPI_COMM_WORLD, &nprocs );
00285 #endif
00286 
00287     moab::Interface* mbCore = new( std::nothrow ) moab::Core;
00288 
00289     if( NULL == mbCore ) { return 1; }
00290 
00291     ToolContext* runCtx;
00292 #ifdef MOAB_HAVE_MPI
00293     moab::ParallelComm* pcomm = new moab::ParallelComm( mbCore, MPI_COMM_WORLD, 0 );
00294 
00295     runCtx                   = new ToolContext( mbCore, pcomm );
00296     const char* writeOptions = ( nprocs > 1 ? "PARALLEL=WRITE_PART" : "" );
00297 #else
00298     runCtx                   = new ToolContext( mbCore );
00299     const char* writeOptions = "";
00300 #endif
00301     runCtx->ParseCLOptions( argc, argv );
00302 
00303     const double radius_src  = 1.0 /*2.0*acos(-1.0)*/;
00304     const double radius_dest = 1.0 /*2.0*acos(-1.0)*/;
00305 
00306     moab::DebugOutput& outputFormatter = runCtx->outputFormatter;
00307 
00308 #ifdef MOAB_HAVE_MPI
00309     moab::TempestRemapper remapper( mbCore, pcomm );
00310 #else
00311     moab::TempestRemapper remapper( mbCore );
00312 #endif
00313     remapper.meshValidate     = true;
00314     remapper.constructEdgeMap = false;
00315     remapper.initialize();
00316 
00317     // Default area_method = lHuiller; Options: Girard, GaussQuadrature (if TR is available)
00318 #ifdef MOAB_HAVE_TEMPESTREMAP
00319     moab::IntxAreaUtils areaAdaptor( moab::IntxAreaUtils::GaussQuadrature );
00320 #else
00321     moab::IntxAreaUtils areaAdaptor( moab::IntxAreaUtils::lHuiller );
00322 #endif
00323 
00324     Mesh* tempest_mesh = new Mesh();
00325     runCtx->timer_push( "create Tempest mesh" );
00326     rval = CreateTempestMesh( *runCtx, remapper, tempest_mesh );MB_CHK_ERR( rval );
00327     runCtx->timer_pop();
00328 
00329     // Some constant parameters
00330     const double epsrel = 1.e-15;
00331     const double boxeps = 1e-6;
00332 
00333     if( runCtx->meshType == moab::TempestRemapper::OVERLAP_MEMORY )
00334     {
00335         // Compute intersections with MOAB
00336         // For the overlap method, choose between: "fuzzy", "exact" or "mixed"
00337         assert( runCtx->meshes.size() == 3 );
00338 
00339 #ifdef MOAB_HAVE_MPI
00340         rval = pcomm->check_all_shared_handles();MB_CHK_ERR( rval );
00341 #endif
00342 
00343         // Load the meshes and validate
00344         rval = remapper.ConvertTempestMesh( moab::Remapper::SourceMesh );MB_CHK_ERR( rval );
00345         rval = remapper.ConvertTempestMesh( moab::Remapper::TargetMesh );MB_CHK_ERR( rval );
00346         remapper.SetMeshType( moab::Remapper::OverlapMesh, moab::TempestRemapper::OVERLAP_FILES );
00347         rval = remapper.ConvertTempestMesh( moab::Remapper::OverlapMesh );MB_CHK_ERR( rval );
00348         rval = mbCore->write_mesh( "tempest_intersection.h5m", &runCtx->meshsets[2], 1 );MB_CHK_ERR( rval );
00349 
00350         // print verbosely about the problem setting
00351         {
00352             moab::Range rintxverts, rintxelems;
00353             rval = mbCore->get_entities_by_dimension( runCtx->meshsets[0], 0, rintxverts );MB_CHK_ERR( rval );
00354             rval = mbCore->get_entities_by_dimension( runCtx->meshsets[0], 2, rintxelems );MB_CHK_ERR( rval );
00355             outputFormatter.printf( 0, "The source set contains %lu vertices and %lu elements \n", rintxverts.size(),
00356                                     rintxelems.size() );
00357 
00358             moab::Range bintxverts, bintxelems;
00359             rval = mbCore->get_entities_by_dimension( runCtx->meshsets[1], 0, bintxverts );MB_CHK_ERR( rval );
00360             rval = mbCore->get_entities_by_dimension( runCtx->meshsets[1], 2, bintxelems );MB_CHK_ERR( rval );
00361             outputFormatter.printf( 0, "The target set contains %lu vertices and %lu elements \n", bintxverts.size(),
00362                                     bintxelems.size() );
00363         }
00364 
00365         moab::EntityHandle intxset;  // == remapper.GetMeshSet(moab::Remapper::OverlapMesh);
00366 
00367         // Compute intersections with MOAB
00368         {
00369             // Create the intersection on the sphere object
00370             runCtx->timer_push( "setup the intersector" );
00371 
00372             moab::Intx2MeshOnSphere* mbintx = new moab::Intx2MeshOnSphere( mbCore );
00373             mbintx->set_error_tolerance( epsrel );
00374             mbintx->set_box_error( boxeps );
00375             mbintx->set_radius_source_mesh( radius_src );
00376             mbintx->set_radius_destination_mesh( radius_dest );
00377 #ifdef MOAB_HAVE_MPI
00378             mbintx->set_parallel_comm( pcomm );
00379 #endif
00380             rval = mbintx->FindMaxEdges( runCtx->meshsets[0], runCtx->meshsets[1] );MB_CHK_ERR( rval );
00381 
00382 #ifdef MOAB_HAVE_MPI
00383             moab::Range local_verts;
00384             rval = mbintx->build_processor_euler_boxes( runCtx->meshsets[1], local_verts );MB_CHK_ERR( rval );
00385 
00386             runCtx->timer_pop();
00387 
00388             moab::EntityHandle covering_set;
00389             runCtx->timer_push( "communicate the mesh" );
00390             rval = mbintx->construct_covering_set( runCtx->meshsets[0], covering_set );MB_CHK_ERR( rval );  // lots of communication if mesh is distributed very differently
00391             runCtx->timer_pop();
00392 #else
00393             moab::EntityHandle covering_set = runCtx->meshsets[0];
00394 #endif
00395             // Now let's invoke the MOAB intersection algorithm in parallel with a
00396             // source and target mesh set representing two different decompositions
00397             runCtx->timer_push( "compute intersections with MOAB" );
00398             rval = mbCore->create_meshset( moab::MESHSET_SET, intxset );MB_CHK_SET_ERR( rval, "Can't create new set" );
00399             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" );
00400             runCtx->timer_pop();
00401 
00402             // free the memory
00403             delete mbintx;
00404         }
00405 
00406         {
00407             moab::Range intxelems, intxverts;
00408             rval = mbCore->get_entities_by_dimension( intxset, 2, intxelems );MB_CHK_ERR( rval );
00409             rval = mbCore->get_entities_by_dimension( intxset, 0, intxverts, true );MB_CHK_ERR( rval );
00410             outputFormatter.printf( 0, "The intersection set contains %lu elements and %lu vertices \n",
00411                                     intxelems.size(), intxverts.size() );
00412 
00413             double initial_sarea =
00414                 areaAdaptor.area_on_sphere( mbCore, runCtx->meshsets[0],
00415                                             radius_src );  // use the target to compute the initial area
00416             double initial_tarea =
00417                 areaAdaptor.area_on_sphere( mbCore, runCtx->meshsets[1],
00418                                             radius_dest );  // use the target to compute the initial area
00419             double intx_area = areaAdaptor.area_on_sphere( mbCore, intxset, radius_src );
00420 
00421             outputFormatter.printf( 0, "mesh areas: source = %12.10f, target = %12.10f, intersection = %12.10f \n",
00422                                     initial_sarea, initial_tarea, intx_area );
00423             outputFormatter.printf( 0, "relative error w.r.t source = %12.10e, target = %12.10e \n",
00424                                     fabs( intx_area - initial_sarea ) / initial_sarea,
00425                                     fabs( intx_area - initial_tarea ) / initial_tarea );
00426         }
00427 
00428         // Write out our computed intersection file
00429         rval = mbCore->write_mesh( "moab_intersection.h5m", &intxset, 1 );MB_CHK_ERR( rval );
00430 
00431         if( runCtx->computeWeights )
00432         {
00433             runCtx->timer_push( "compute weights with the Tempest meshes" );
00434             // Call to generate an offline map with the tempest meshes
00435             OfflineMap weightMap;
00436             int err = GenerateOfflineMapWithMeshes(
00437                 weightMap, *runCtx->meshes[0], *runCtx->meshes[1], *runCtx->meshes[2], "",
00438                 "",  // std::string strInputMeta, std::string strOutputMeta,
00439                 runCtx->disc_methods[0],
00440                 runCtx->disc_methods[1],                         // std::string strInputType, std::string strOutputType,
00441                 runCtx->disc_orders[0], runCtx->disc_orders[1],  // int nPin=4, int nPout=4,
00442                 runCtx->fNoBubble, true,
00443                 runCtx->ensureMonotonicity  // bool fNoBubble = false, bool fCorrectAreas = false,
00444                                             // int fMonotoneTypeID = 0
00445             );
00446             runCtx->timer_pop();
00447 
00448             std::map< std::string, std::string > mapAttributes;
00449             if( err ) { rval = moab::MB_FAILURE; }
00450             else
00451             {
00452                 weightMap.Write( "outWeights.nc", mapAttributes );
00453             }
00454         }
00455     }
00456     else if( runCtx->meshType == moab::TempestRemapper::OVERLAP_MOAB )
00457     {
00458         // Usage: mpiexec -n 2 tools/mbtempest -t 5 -l mycs_2.h5m -l myico_2.h5m -f myoverlap_2.h5m
00459 #ifdef MOAB_HAVE_MPI
00460         rval = pcomm->check_all_shared_handles();MB_CHK_ERR( rval );
00461 #endif
00462         // print verbosely about the problem setting
00463         {
00464             moab::Range rintxverts, rintxelems;
00465             rval = mbCore->get_entities_by_dimension( runCtx->meshsets[0], 0, rintxverts );MB_CHK_ERR( rval );
00466             rval = mbCore->get_entities_by_dimension( runCtx->meshsets[0], 2, rintxelems );MB_CHK_ERR( rval );
00467             rval = moab::IntxUtils::fix_degenerate_quads( mbCore, runCtx->meshsets[0] );MB_CHK_ERR( rval );
00468             if (runCtx->enforceConvexity)
00469             {
00470                 rval = moab::IntxUtils::enforce_convexity( mbCore, runCtx->meshsets[0], proc_id );MB_CHK_ERR( rval );
00471             }
00472             rval = areaAdaptor.positive_orientation( mbCore, runCtx->meshsets[0], radius_src );MB_CHK_ERR( rval );
00473             if( !proc_id )
00474                 outputFormatter.printf( 0, "The source set contains %lu vertices and %lu elements \n",
00475                                         rintxverts.size(), rintxelems.size() );
00476 
00477             moab::Range bintxverts, bintxelems;
00478             rval = mbCore->get_entities_by_dimension( runCtx->meshsets[1], 0, bintxverts );MB_CHK_ERR( rval );
00479             rval = mbCore->get_entities_by_dimension( runCtx->meshsets[1], 2, bintxelems );MB_CHK_ERR( rval );
00480             rval = moab::IntxUtils::fix_degenerate_quads( mbCore, runCtx->meshsets[1] );MB_CHK_ERR( rval );
00481             if (runCtx->enforceConvexity)
00482             {
00483                 rval = moab::IntxUtils::enforce_convexity( mbCore, runCtx->meshsets[1], proc_id );MB_CHK_ERR( rval );
00484             }
00485             rval = areaAdaptor.positive_orientation( mbCore, runCtx->meshsets[1], radius_dest );MB_CHK_ERR( rval );
00486             if( !proc_id )
00487                 outputFormatter.printf( 0, "The target set contains %lu vertices and %lu elements \n",
00488                                         bintxverts.size(), bintxelems.size() );
00489         }
00490 
00491         // First compute the covering set such that the target elements are fully covered by the
00492         // lcoal source grid
00493         runCtx->timer_push( "construct covering set for intersection" );
00494         rval = remapper.ConstructCoveringSet( epsrel, 1.0, 1.0, boxeps, runCtx->rrmGrids );MB_CHK_ERR( rval );
00495         runCtx->timer_pop();
00496 
00497         // Compute intersections with MOAB with either the Kd-tree or the advancing front algorithm
00498         runCtx->timer_push( "setup and compute mesh intersections" );
00499         rval = remapper.ComputeOverlapMesh( runCtx->kdtreeSearch, false );MB_CHK_ERR( rval );
00500         runCtx->timer_pop();
00501 
00502         // print some diagnostic checks to see if the overlap grid resolved the input meshes
00503         // correctly
00504         {
00505             double local_areas[3],
00506                 global_areas[3];  // Array for Initial area, and through Method 1 and Method 2
00507             // local_areas[0] = area_on_sphere_lHuiller ( mbCore, runCtx->meshsets[1], radius_src );
00508             local_areas[0] = areaAdaptor.area_on_sphere( mbCore, runCtx->meshsets[0], radius_src );
00509             local_areas[1] = areaAdaptor.area_on_sphere( mbCore, runCtx->meshsets[1], radius_dest );
00510             local_areas[2] = areaAdaptor.area_on_sphere( mbCore, runCtx->meshsets[2], radius_src );
00511 
00512 #ifdef MOAB_HAVE_MPI
00513             MPI_Allreduce( &local_areas[0], &global_areas[0], 3, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD );
00514 #else
00515             global_areas[0]                 = local_areas[0];
00516             global_areas[1]                 = local_areas[1];
00517             global_areas[2]                 = local_areas[2];
00518 #endif
00519             if( !proc_id )
00520             {
00521                 outputFormatter.printf( 0,
00522                                         "initial area: source mesh = %12.14f, target mesh = "
00523                                         "%12.14f, overlap mesh = %12.14f\n",
00524                                         global_areas[0], global_areas[1], global_areas[2] );
00525                 // outputFormatter.printf ( 0, " area with l'Huiller: %12.14f with Girard:
00526                 // %12.14f\n", global_areas[2], global_areas[3] ); outputFormatter.printf ( 0, "
00527                 // relative difference areas = %12.10e\n", fabs ( global_areas[2] - global_areas[3]
00528                 // ) / global_areas[2] );
00529                 outputFormatter.printf( 0, "relative error w.r.t source = %12.14e, and target = %12.14e\n",
00530                                         fabs( global_areas[0] - global_areas[2] ) / global_areas[0],
00531                                         fabs( global_areas[1] - global_areas[2] ) / global_areas[1] );
00532             }
00533         }
00534 
00535         if( runCtx->intxFilename.size() )
00536         {
00537             moab::EntityHandle writableOverlapSet;
00538             rval = mbCore->create_meshset( moab::MESHSET_SET, writableOverlapSet );MB_CHK_SET_ERR( rval, "Can't create new set" );
00539             moab::EntityHandle meshOverlapSet = remapper.GetMeshSet( moab::Remapper::OverlapMesh );
00540             moab::Range ovEnts;
00541             rval = mbCore->get_entities_by_dimension( meshOverlapSet, 2, ovEnts );MB_CHK_SET_ERR( rval, "Can't create new set" );
00542             rval = mbCore->get_entities_by_dimension( meshOverlapSet, 0, ovEnts );MB_CHK_SET_ERR( rval, "Can't create new set" );
00543 
00544 #ifdef MOAB_HAVE_MPI
00545             // Do not remove ghosted entities if we still haven't computed weights
00546             // Remove ghosted entities from overlap set before writing the new mesh set to file
00547             if( nprocs > 1 )
00548             {
00549                 moab::Range ghostedEnts;
00550                 rval = remapper.GetOverlapAugmentedEntities( ghostedEnts );MB_CHK_ERR( rval );
00551                 ovEnts = moab::subtract( ovEnts, ghostedEnts );
00552             }
00553 #endif
00554             rval = mbCore->add_entities( writableOverlapSet, ovEnts );MB_CHK_SET_ERR( rval, "Deleting ghosted entities failed" );
00555 
00556             size_t lastindex = runCtx->intxFilename.find_last_of( "." );
00557             sstr.str( "" );
00558             sstr << runCtx->intxFilename.substr( 0, lastindex ) << ".h5m";
00559             if( !runCtx->proc_id )
00560                 std::cout << "Writing out the MOAB intersection mesh file to " << sstr.str() << std::endl;
00561 
00562             // Write out our computed intersection file
00563             rval = mbCore->write_file( sstr.str().c_str(), NULL, writeOptions, &writableOverlapSet, 1 );MB_CHK_ERR( rval );
00564         }
00565 
00566         if( runCtx->computeWeights )
00567         {
00568             runCtx->meshes[2] = remapper.GetMesh( moab::Remapper::OverlapMesh );
00569             if( !runCtx->proc_id ) std::cout << std::endl;
00570 
00571             runCtx->timer_push( "setup computation of weights" );
00572             // Call to generate the remapping weights with the tempest meshes
00573             moab::TempestOnlineMap* weightMap = new moab::TempestOnlineMap( &remapper );
00574             runCtx->timer_pop();
00575 
00576             runCtx->timer_push( "compute weights with TempestRemap" );
00577 
00578             rval = weightMap->GenerateRemappingWeights(
00579                 runCtx->disc_methods[0],
00580                 runCtx->disc_methods[1],                         // std::string strInputType, std::string strOutputType,
00581                 runCtx->disc_orders[0], runCtx->disc_orders[1],  // int nPin=4, int nPout=4,
00582                 runCtx->fNoBubble,
00583                 runCtx->ensureMonotonicity,  // bool fNoBubble=true, int fMonotoneTypeID=0,
00584                 runCtx->fVolumetric, runCtx->fNoConservation,
00585                 runCtx->fNoCheck,  // bool fVolumetric=false, bool fNoConservation=false, bool
00586                                    // fNoCheck=false,
00587                 runCtx->doftag_names[0],
00588                 runCtx->doftag_names[1],  // std::string source_tag_name, std::string
00589                                           // target_tag_name,
00590                 "",                       //"",                                                       // std::string
00591                                           // strVariables="",
00592                 "", "",                   // std::string strInputData="", std::string strOutputData="",
00593                 "", true,                 // std::string strNColName="", bool fOutputDouble=true,
00594                 "", false,
00595                 0.0,  // std::string strPreserveVariables="", bool fPreserveAll=false, double
00596                       // dFillValueOverride=0.0,
00597                 runCtx->fInputConcave,
00598                 runCtx->fOutputConcave  // bool fInputConcave = false, bool fOutputConcave = false
00599             );MB_CHK_ERR( rval );
00600             runCtx->timer_pop();
00601 
00602             // Invoke the CheckMap routine on the TempestRemap serial interface directly, if running
00603             // on a single process
00604             if( nprocs == 1 )
00605             {
00606                 const double dNormalTolerance = 1.0E-8;
00607                 const double dStrictTolerance = 1.0E-12;
00608                 weightMap->CheckMap( !runCtx->fNoCheck, !runCtx->fNoCheck,
00609                                      !runCtx->fNoCheck && ( runCtx->ensureMonotonicity ), dNormalTolerance,
00610                                      dStrictTolerance );
00611             }
00612 
00613             if( runCtx->outFilename.size() )
00614             {
00615                                 // Write the map file to disk in parallel using either HDF5 or SCRIP interface
00616                                 rval = weightMap->WriteParallelMap( runCtx->outFilename.c_str() );MB_CHK_ERR( rval );
00617 
00618                 // Write out the metadata information for the map file
00619                 if( proc_id == 0 )
00620                 {
00621                     size_t lastindex      = runCtx->outFilename.find_last_of( "." );
00622                     sstr.str( "" );
00623                     sstr << runCtx->outFilename.substr( 0, lastindex ) << ".meta";
00624 
00625                     std::ofstream metafile( sstr.str() );
00626                     metafile << "Generator = MOAB-TempestRemap (mbtempest) Offline Regridding "
00627                                 "Weight Generator"
00628                              << std::endl;
00629                     metafile << "domain_a = " << runCtx->inFilenames[0] << std::endl;
00630                     metafile << "domain_b = " << runCtx->inFilenames[1] << std::endl;
00631                     metafile << "domain_aNb = "
00632                              << ( runCtx->intxFilename.size() ? runCtx->intxFilename : "outOverlap.h5m" ) << std::endl;
00633                     metafile << "map_aPb = " << runCtx->outFilename << std::endl;
00634                     metafile << "type_src = " << runCtx->disc_methods[0] << std::endl;
00635                     metafile << "np_src = " << runCtx->disc_orders[0] << std::endl;
00636                     metafile << "concave_src = " << ( runCtx->fInputConcave ? "true" : "false" ) << std::endl;
00637                     metafile << "type_dst = " << runCtx->disc_methods[1] << std::endl;
00638                     metafile << "np_dst = " << runCtx->disc_orders[1] << std::endl;
00639                     metafile << "concave_dst = " << ( runCtx->fOutputConcave ? "true" : "false" ) << std::endl;
00640                     metafile << "mono_type = " << runCtx->ensureMonotonicity << std::endl;
00641                     metafile << "bubble = " << ( runCtx->fNoBubble ? "false" : "true" ) << std::endl;
00642                     metafile << "version = "
00643                              << "MOAB v5.1.0+" << std::endl;
00644                     metafile.close();
00645                 }
00646             }
00647 
00648             if( runCtx->verifyWeights )
00649             {
00650                 // Let us pick a sampling test function for solution evaluation
00651                 moab::TempestOnlineMap::sample_function testFunction =
00652                     &sample_stationary_vortex;  // &sample_slow_harmonic;
00653 
00654                 runCtx->timer_push( "describe a solution on source grid" );
00655                 moab::Tag srcAnalyticalFunction;
00656                 rval = weightMap->DefineAnalyticalSolution( srcAnalyticalFunction, "AnalyticalSolnSrcExact",
00657                                                             moab::Remapper::SourceMesh, testFunction );MB_CHK_ERR( rval );
00658                 runCtx->timer_pop();
00659                 // rval = mbCore->write_file ( "srcWithSolnTag.h5m", NULL, writeOptions,
00660                 // &runCtx->meshsets[0], 1 ); MB_CHK_ERR ( rval );
00661 
00662                 runCtx->timer_push( "describe a solution on target grid" );
00663                 moab::Tag tgtAnalyticalFunction;
00664                 moab::Tag tgtProjectedFunction;
00665                 rval = weightMap->DefineAnalyticalSolution( tgtAnalyticalFunction, "AnalyticalSolnTgtExact",
00666                                                             moab::Remapper::TargetMesh, testFunction,
00667                                                             &tgtProjectedFunction, "ProjectedSolnTgt" );MB_CHK_ERR( rval );
00668                 // rval = mbCore->write_file ( "tgtWithSolnTag.h5m", NULL, writeOptions,
00669                 // &runCtx->meshsets[1], 1 ); MB_CHK_ERR ( rval );
00670                 runCtx->timer_pop();
00671 
00672                 runCtx->timer_push( "compute solution projection on target grid" );
00673                 rval = weightMap->ApplyWeights( srcAnalyticalFunction, tgtProjectedFunction );MB_CHK_ERR( rval );
00674                 runCtx->timer_pop();
00675                 rval = mbCore->write_file( "tgtWithSolnTag2.h5m", NULL, writeOptions, &runCtx->meshsets[1], 1 );MB_CHK_ERR( rval );
00676 
00677                 runCtx->timer_push( "compute error metrics against analytical solution on target grid" );
00678                 std::map< std::string, double > errMetrics;
00679                 rval = weightMap->ComputeMetrics( moab::Remapper::TargetMesh, tgtAnalyticalFunction,
00680                                                   tgtProjectedFunction, errMetrics, true );MB_CHK_ERR( rval );
00681                 runCtx->timer_pop();
00682             }
00683 
00684             delete weightMap;
00685         }
00686     }
00687 
00688     // Clean up
00689     remapper.clear();
00690     delete runCtx;
00691     delete mbCore;
00692 
00693 #ifdef MOAB_HAVE_MPI
00694     MPI_Finalize();
00695 #endif
00696     exit( 0 );
00697 }
00698 
00699 static moab::ErrorCode CreateTempestMesh( ToolContext& ctx, moab::TempestRemapper& remapper, Mesh* tempest_mesh )
00700 {
00701     moab::ErrorCode rval = moab::MB_SUCCESS;
00702     int err;
00703     moab::DebugOutput& outputFormatter = ctx.outputFormatter;
00704 
00705     if( !ctx.proc_id ) { outputFormatter.printf( 0, "Creating TempestRemap Mesh object ...\n" ); }
00706 
00707     if( ctx.meshType == moab::TempestRemapper::OVERLAP_FILES )
00708     {
00709         // For the overlap method, choose between: "fuzzy", "exact" or "mixed"
00710         err = GenerateOverlapMesh( ctx.inFilenames[0], ctx.inFilenames[1], *tempest_mesh, ctx.outFilename, "NetCDF4",
00711                                    "exact", true );
00712 
00713         if( err ) { rval = moab::MB_FAILURE; }
00714         else
00715         {
00716             ctx.meshes.push_back( tempest_mesh );
00717         }
00718     }
00719     else if( ctx.meshType == moab::TempestRemapper::OVERLAP_MEMORY )
00720     {
00721         // Load the meshes and validate
00722         ctx.meshsets.resize( 3 );
00723         ctx.meshes.resize( 3 );
00724         ctx.meshsets[0] = remapper.GetMeshSet( moab::Remapper::SourceMesh );
00725         ctx.meshsets[1] = remapper.GetMeshSet( moab::Remapper::TargetMesh );
00726         ctx.meshsets[2] = remapper.GetMeshSet( moab::Remapper::OverlapMesh );
00727 
00728         // First the source
00729         rval = remapper.LoadMesh( moab::Remapper::SourceMesh, ctx.inFilenames[0], moab::TempestRemapper::DEFAULT );MB_CHK_ERR( rval );
00730         ctx.meshes[0] = remapper.GetMesh( moab::Remapper::SourceMesh );
00731 
00732         // Next the target
00733         rval = remapper.LoadMesh( moab::Remapper::TargetMesh, ctx.inFilenames[1], moab::TempestRemapper::DEFAULT );MB_CHK_ERR( rval );
00734         ctx.meshes[1] = remapper.GetMesh( moab::Remapper::TargetMesh );
00735 
00736         // Now let us construct the overlap mesh, by calling TempestRemap interface directly
00737         // For the overlap method, choose between: "fuzzy", "exact" or "mixed"
00738         err = GenerateOverlapWithMeshes( *ctx.meshes[0], *ctx.meshes[1], *tempest_mesh, "" /*ctx.outFilename*/,
00739                                          "NetCDF4", "exact", false );
00740 
00741         if( err ) { rval = moab::MB_FAILURE; }
00742         else
00743         {
00744             remapper.SetMesh( moab::Remapper::OverlapMesh, tempest_mesh );
00745             ctx.meshes[2] = remapper.GetMesh( moab::Remapper::OverlapMesh );
00746             // ctx.meshes.push_back(*tempest_mesh);
00747         }
00748     }
00749     else if( ctx.meshType == moab::TempestRemapper::OVERLAP_MOAB )
00750     {
00751         ctx.meshsets.resize( 3 );
00752         ctx.meshes.resize( 3 );
00753         ctx.meshsets[0] = remapper.GetMeshSet( moab::Remapper::SourceMesh );
00754         ctx.meshsets[1] = remapper.GetMeshSet( moab::Remapper::TargetMesh );
00755         ctx.meshsets[2] = remapper.GetMeshSet( moab::Remapper::OverlapMesh );
00756 
00757         const double radius_src  = 1.0 /*2.0*acos(-1.0)*/;
00758         const double radius_dest = 1.0 /*2.0*acos(-1.0)*/;
00759 
00760         const char* additional_read_opts = ( ctx.n_procs > 1 ? "NO_SET_CONTAINING_PARENTS;" : "" );
00761         // Load the source mesh and validate
00762         rval = remapper.LoadNativeMesh( ctx.inFilenames[0], ctx.meshsets[0], additional_read_opts );MB_CHK_ERR( rval );
00763         // Rescale the radius of both to compute the intersection
00764         rval = moab::IntxUtils::ScaleToRadius( ctx.mbcore, ctx.meshsets[0], radius_src );MB_CHK_ERR( rval );
00765         rval = remapper.ConvertMeshToTempest( moab::Remapper::SourceMesh );MB_CHK_ERR( rval );
00766         ctx.meshes[0] = remapper.GetMesh( moab::Remapper::SourceMesh );
00767 
00768         // Load the target mesh and validate
00769         rval = remapper.LoadNativeMesh( ctx.inFilenames[1], ctx.meshsets[1], additional_read_opts );MB_CHK_ERR( rval );
00770         rval = moab::IntxUtils::ScaleToRadius( ctx.mbcore, ctx.meshsets[1], radius_dest );MB_CHK_ERR( rval );
00771         rval = remapper.ConvertMeshToTempest( moab::Remapper::TargetMesh );MB_CHK_ERR( rval );
00772         ctx.meshes[1] = remapper.GetMesh( moab::Remapper::TargetMesh );
00773     }
00774     else if( ctx.meshType == moab::TempestRemapper::ICO )
00775     {
00776         err = GenerateICOMesh( *tempest_mesh, ctx.blockSize, ctx.computeDual, ctx.outFilename, "NetCDF4" );
00777 
00778         if( err ) { rval = moab::MB_FAILURE; }
00779         else
00780         {
00781             ctx.meshes.push_back( tempest_mesh );
00782         }
00783     }
00784     else if( ctx.meshType == moab::TempestRemapper::RLL )
00785     {
00786         err = GenerateRLLMesh( *tempest_mesh,                     // Mesh& meshOut,
00787                                ctx.blockSize * 2, ctx.blockSize,  // int nLongitudes, int nLatitudes,
00788                                0.0, 360.0,                        // double dLonBegin, double dLonEnd,
00789                                -90.0, 90.0,                       // double dLatBegin, double dLatEnd,
00790                                false, false, false,  // bool fGlobalCap, bool fFlipLatLon, bool fForceGlobal,
00791                                "" /*ctx.inFilename*/, "",
00792                                "",  // std::string strInputFile, std::string strInputFileLonName, std::string
00793                                     // strInputFileLatName,
00794                                ctx.outFilename, "NetCDF4",  // std::string strOutputFile, std::string strOutputFormat
00795                                true                         // bool fVerbose
00796         );
00797 
00798         if( err ) { rval = moab::MB_FAILURE; }
00799         else
00800         {
00801             ctx.meshes.push_back( tempest_mesh );
00802         }
00803     }
00804     else  // default
00805     {
00806         err = GenerateCSMesh( *tempest_mesh, ctx.blockSize, ctx.outFilename, "NetCDF4" );
00807 
00808         if( err ) { rval = moab::MB_FAILURE; }
00809         else
00810         {
00811             ctx.meshes.push_back( tempest_mesh );
00812         }
00813     }
00814 
00815     if( ctx.meshType != moab::TempestRemapper::OVERLAP_MOAB && !tempest_mesh )
00816     {
00817         std::cout << "Tempest Mesh is not a complete object; Quitting...";
00818         exit( -1 );
00819     }
00820 
00821     return rval;
00822 }
00823 
00824 ///////////////////////////////////////////////
00825 //         Test functions
00826 
00827 double sample_slow_harmonic( double dLon, double dLat )
00828 {
00829     return ( 2.0 + cos( dLat ) * cos( dLat ) * cos( 2.0 * dLon ) );
00830 }
00831 
00832 double sample_fast_harmonic( double dLon, double dLat )
00833 {
00834     return ( 2.0 + pow( sin( 2.0 * dLat ), 16.0 ) * cos( 16.0 * dLon ) );
00835     // return (2.0 + pow(cos(2.0 * dLat), 16.0) * cos(16.0 * dLon));
00836 }
00837 
00838 double sample_constant( double /*dLon*/, double /*dLat*/ )
00839 {
00840     return 1.0;
00841 }
00842 
00843 double sample_stationary_vortex( double dLon, double dLat )
00844 {
00845     const double dLon0 = 0.0;
00846     const double dLat0 = 0.6;
00847     const double dR0   = 3.0;
00848     const double dD    = 5.0;
00849     const double dT    = 6.0;
00850 
00851     ///     Find the rotated longitude and latitude of a point on a sphere
00852     ///     with pole at (dLonC, dLatC).
00853     {
00854         double dSinC = sin( dLat0 );
00855         double dCosC = cos( dLat0 );
00856         double dCosT = cos( dLat );
00857         double dSinT = sin( dLat );
00858 
00859         double dTrm = dCosT * cos( dLon - dLon0 );
00860         double dX   = dSinC * dTrm - dCosC * dSinT;
00861         double dY   = dCosT * sin( dLon - dLon0 );
00862         double dZ   = dSinC * dSinT + dCosC * dTrm;
00863 
00864         dLon = atan2( dY, dX );
00865         if( dLon < 0.0 ) { dLon += 2.0 * M_PI; }
00866         dLat = asin( dZ );
00867     }
00868 
00869     double dRho = dR0 * cos( dLat );
00870     double dVt  = 3.0 * sqrt( 3.0 ) / 2.0 / cosh( dRho ) / cosh( dRho ) * tanh( dRho );
00871 
00872     double dOmega;
00873     if( dRho == 0.0 ) { dOmega = 0.0; }
00874     else
00875     {
00876         dOmega = dVt / dRho;
00877     }
00878 
00879     return ( 1.0 - tanh( dRho / dD * sin( dLon - dOmega * dT ) ) );
00880 }
00881 
00882 ///////////////////////////////////////////////
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines