|
MOAB: Mesh Oriented datABase
(version 5.4.1)
|
|
MOAB: Mesh Oriented datABase
(version 5.4.1)
|
#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) |
| int | main (int argc, char *argv[]) |
| static moab::ErrorCode CreateTempestMesh | ( | ToolContext & | ctx, |
| moab::TempestRemapper & | remapper, | ||
| Mesh * | tempest_mesh | ||
| ) | [static] |
Definition at line 786 of file mbtempest.cpp.
References ToolContext::blockSize, ToolContext::computeDual, moab::TempestRemapper::ConvertMeshToTempest(), moab::TempestRemapper::DEFAULT, ErrorCode, 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::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::SetMeshType(), moab::Remapper::SourceMesh, and moab::Remapper::TargetMesh.
Referenced by main().
{
moab::ErrorCode rval = moab::MB_SUCCESS;
int err;
moab::DebugOutput& outputFormatter = ctx.outputFormatter;
if( !ctx.proc_id )
{
outputFormatter.printf( 0, "Creating TempestRemap Mesh object ...\n" );
}
if( ctx.meshType == moab::TempestRemapper::OVERLAP_FILES )
{
// For the overlap method, choose between: "fuzzy", "exact" or "mixed"
err = GenerateOverlapMesh( ctx.inFilenames[0], ctx.inFilenames[1], *tempest_mesh, ctx.outFilename, "NetCDF4",
"exact", true );
if( err )
{
rval = moab::MB_FAILURE;
}
else
{
ctx.meshes.push_back( tempest_mesh );
}
}
else if( ctx.meshType == moab::TempestRemapper::OVERLAP_MEMORY )
{
// Load the meshes and validate
ctx.meshsets.resize( 3 );
ctx.meshes.resize( 3 );
ctx.meshsets[0] = remapper.GetMeshSet( moab::Remapper::SourceMesh );
ctx.meshsets[1] = remapper.GetMeshSet( moab::Remapper::TargetMesh );
ctx.meshsets[2] = remapper.GetMeshSet( moab::Remapper::OverlapMesh );
// First the source
rval = remapper.LoadMesh( moab::Remapper::SourceMesh, ctx.inFilenames[0], moab::TempestRemapper::DEFAULT );MB_CHK_ERR( rval );
ctx.meshes[0] = remapper.GetMesh( moab::Remapper::SourceMesh );
// Next the target
rval = remapper.LoadMesh( moab::Remapper::TargetMesh, ctx.inFilenames[1], moab::TempestRemapper::DEFAULT );MB_CHK_ERR( rval );
ctx.meshes[1] = remapper.GetMesh( moab::Remapper::TargetMesh );
// Now let us construct the overlap mesh, by calling TempestRemap interface directly
// For the overlap method, choose between: "fuzzy", "exact" or "mixed"
err = GenerateOverlapWithMeshes( *ctx.meshes[0], *ctx.meshes[1], *tempest_mesh, "" /*ctx.outFilename*/,
"NetCDF4", "exact", false );
if( err )
{
rval = moab::MB_FAILURE;
}
else
{
remapper.SetMesh( moab::Remapper::OverlapMesh, tempest_mesh );
ctx.meshes[2] = remapper.GetMesh( moab::Remapper::OverlapMesh );
// ctx.meshes.push_back(*tempest_mesh);
}
}
else if( ctx.meshType == moab::TempestRemapper::OVERLAP_MOAB )
{
ctx.meshsets.resize( 3 );
ctx.meshes.resize( 3 );
ctx.meshsets[0] = remapper.GetMeshSet( moab::Remapper::SourceMesh );
ctx.meshsets[1] = remapper.GetMeshSet( moab::Remapper::TargetMesh );
ctx.meshsets[2] = remapper.GetMeshSet( moab::Remapper::OverlapMesh );
const double radius_src = 1.0 /*2.0*acos(-1.0)*/;
const double radius_dest = 1.0 /*2.0*acos(-1.0)*/;
const char* additional_read_opts = ( ctx.n_procs > 1 ? "NO_SET_CONTAINING_PARENTS;" : "" );
std::vector< int > smetadata, tmetadata;
// Load the source mesh and validate
rval = remapper.LoadNativeMesh( ctx.inFilenames[0], ctx.meshsets[0], smetadata, additional_read_opts );MB_CHK_ERR( rval );
if( smetadata.size() )
{
// remapper.SetMeshType( moab::Remapper::SourceMesh,
// static_cast< moab::TempestRemapper::TempestMeshType >( smetadata[0] ), &smetadata
// );
remapper.SetMeshType( moab::Remapper::SourceMesh,
static_cast< moab::TempestRemapper::TempestMeshType >( smetadata[0] ) );
}
// Rescale the radius of both to compute the intersection
rval = moab::IntxUtils::ScaleToRadius( ctx.mbcore, ctx.meshsets[0], radius_src );MB_CHK_ERR( rval );
rval = remapper.ConvertMeshToTempest( moab::Remapper::SourceMesh );MB_CHK_ERR( rval );
ctx.meshes[0] = remapper.GetMesh( moab::Remapper::SourceMesh );
// Load the target mesh and validate
rval = remapper.LoadNativeMesh( ctx.inFilenames[1], ctx.meshsets[1], tmetadata, additional_read_opts );MB_CHK_ERR( rval );
if( tmetadata.size() )
{
// remapper.SetMeshType( moab::Remapper::TargetMesh,
// static_cast< moab::TempestRemapper::TempestMeshType >( tmetadata[0] ), &tmetadata
// );
remapper.SetMeshType( moab::Remapper::TargetMesh,
static_cast< moab::TempestRemapper::TempestMeshType >( tmetadata[0] ) );
}
rval = moab::IntxUtils::ScaleToRadius( ctx.mbcore, ctx.meshsets[1], radius_dest );MB_CHK_ERR( rval );
rval = remapper.ConvertMeshToTempest( moab::Remapper::TargetMesh );MB_CHK_ERR( rval );
ctx.meshes[1] = remapper.GetMesh( moab::Remapper::TargetMesh );
}
else if( ctx.meshType == moab::TempestRemapper::ICO )
{
err = GenerateICOMesh( *tempest_mesh, ctx.blockSize, ctx.computeDual, ctx.outFilename, "NetCDF4" );
if( err )
{
rval = moab::MB_FAILURE;
}
else
{
ctx.meshes.push_back( tempest_mesh );
}
}
else if( ctx.meshType == moab::TempestRemapper::RLL )
{
err = GenerateRLLMesh( *tempest_mesh, // Mesh& meshOut,
ctx.blockSize * 2, ctx.blockSize, // int nLongitudes, int nLatitudes,
0.0, 360.0, // double dLonBegin, double dLonEnd,
-90.0, 90.0, // double dLatBegin, double dLatEnd,
false, false, false, // bool fGlobalCap, bool fFlipLatLon, bool fForceGlobal,
"" /*ctx.inFilename*/, "",
"", // std::string strInputFile, std::string strInputFileLonName, std::string
// strInputFileLatName,
ctx.outFilename, "NetCDF4", // std::string strOutputFile, std::string strOutputFormat
true // bool fVerbose
);
if( err )
{
rval = moab::MB_FAILURE;
}
else
{
ctx.meshes.push_back( tempest_mesh );
}
}
else // default
{
err = GenerateCSMesh( *tempest_mesh, ctx.blockSize, ctx.outFilename, "NetCDF4" );
if( err )
{
rval = moab::MB_FAILURE;
}
else
{
ctx.meshes.push_back( tempest_mesh );
}
}
if( ctx.meshType != moab::TempestRemapper::OVERLAP_MOAB && !tempest_mesh )
{
std::cout << "Tempest Mesh is not a complete object; Quitting...";
exit( -1 );
}
return rval;
}
| int main | ( | int | argc, |
| char * | argv[] | ||
| ) |
Definition at line 347 of file mbtempest.cpp.
References moab::Interface::add_entities(), moab::TempestOnlineMap::ApplyWeights(), moab::IntxAreaUtils::area_on_sphere(), ToolContext::baselineFile, 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::Interface::create_meshset(), CreateTempestMesh(), moab::TempestOnlineMap::DefineAnalyticalSolution(), ToolContext::disc_methods, ToolContext::disc_orders, ToolContext::doftag_names, moab::IntxUtils::enforce_convexity(), ToolContext::enforceConvexity, ToolContext::ensureMonotonicity, moab::error(), ErrorCode, ToolContext::fCheck, moab::Intx2Mesh::FindMaxEdges(), moab::IntxUtils::fix_degenerate_quads(), moab::IntxAreaUtils::GaussQuadrature, moab::TempestOnlineMap::GenerateRemappingWeights(), moab::Interface::get_entities_by_dimension(), moab::TempestRemapper::GetMesh(), moab::TempestRemapper::GetMeshSet(), moab::TempestRemapper::GetOverlapAugmentedEntities(), moab::Interface::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, ToolContext::meshes, MESHSET_SET, ToolContext::meshsets, ToolContext::meshType, moab::TempestRemapper::meshValidate, MPI_COMM_WORLD, ToolContext::outFilename, ToolContext::outputFormatter, moab::TempestRemapper::OVERLAP_FILES, moab::TempestRemapper::OVERLAP_MEMORY, moab::TempestRemapper::OVERLAP_MOAB, moab::Remapper::OverlapMesh, ToolContext::ParseCLOptions(), moab::IntxAreaUtils::positive_orientation(), moab::DebugOutput::printf(), ToolContext::proc_id, ToolContext::rrmGrids, 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::TempestRemapper::SetMeshType(), moab::Range::size(), moab::Remapper::SourceMesh, moab::subtract(), moab::Interface::tag_get_data(), moab::Interface::tag_get_handle(), moab::Remapper::TargetMesh, ToolContext::timer_pop(), ToolContext::timer_push(), ToolContext::verifyWeights, moab::Interface::write_file(), moab::Interface::write_mesh(), and moab::TempestOnlineMap::WriteParallelMap().
{
moab::ErrorCode rval;
NcError error( NcError::verbose_nonfatal );
std::stringstream sstr;
int proc_id = 0, nprocs = 1;
#ifdef MOAB_HAVE_MPI
MPI_Init( &argc, &argv );
MPI_Comm_rank( MPI_COMM_WORLD, &proc_id );
MPI_Comm_size( MPI_COMM_WORLD, &nprocs );
#endif
moab::Interface* mbCore = new( std::nothrow ) moab::Core;
if( NULL == mbCore )
{
return 1;
}
ToolContext* runCtx;
#ifdef MOAB_HAVE_MPI
moab::ParallelComm* pcomm = new moab::ParallelComm( mbCore, MPI_COMM_WORLD, 0 );
runCtx = new ToolContext( mbCore, pcomm );
const char* writeOptions = ( nprocs > 1 ? "PARALLEL=WRITE_PART" : "" );
#else
runCtx = new ToolContext( mbCore );
const char* writeOptions = "";
#endif
runCtx->ParseCLOptions( argc, argv );
const double radius_src = 1.0 /*2.0*acos(-1.0)*/;
const double radius_dest = 1.0 /*2.0*acos(-1.0)*/;
moab::DebugOutput& outputFormatter = runCtx->outputFormatter;
#ifdef MOAB_HAVE_MPI
moab::TempestRemapper remapper( mbCore, pcomm );
#else
moab::TempestRemapper remapper( mbCore );
#endif
remapper.meshValidate = true;
remapper.constructEdgeMap = false;
remapper.initialize();
// Default area_method = lHuiller; Options: Girard, GaussQuadrature (if TR is available)
#ifdef MOAB_HAVE_TEMPESTREMAP
moab::IntxAreaUtils areaAdaptor( moab::IntxAreaUtils::GaussQuadrature );
#else
moab::IntxAreaUtils areaAdaptor( moab::IntxAreaUtils::lHuiller );
#endif
Mesh* tempest_mesh = new Mesh();
runCtx->timer_push( "create Tempest mesh" );
rval = CreateTempestMesh( *runCtx, remapper, tempest_mesh );MB_CHK_ERR( rval );
runCtx->timer_pop();
const double epsrel = ReferenceTolerance; // ReferenceTolerance is defined in Defines.h in tempestremap;
// Defines.h is included in SparseMatrix.h
// SparseMatrix.h is included in OfflineMap.h
// OfflineMap.h is included in TempestOnlineMap.hpp
// TempestOnlineMap.hpp is included in this file, and is part of MOAB
// Some constant parameters
const double boxeps = 1e-6;
if( runCtx->meshType == moab::TempestRemapper::OVERLAP_MEMORY )
{
// Compute intersections with MOAB
// For the overlap method, choose between: "fuzzy", "exact" or "mixed"
assert( runCtx->meshes.size() == 3 );
#ifdef MOAB_HAVE_MPI
rval = pcomm->check_all_shared_handles();MB_CHK_ERR( rval );
#endif
// Load the meshes and validate
rval = remapper.ConvertTempestMesh( moab::Remapper::SourceMesh );MB_CHK_ERR( rval );
rval = remapper.ConvertTempestMesh( moab::Remapper::TargetMesh );MB_CHK_ERR( rval );
remapper.SetMeshType( moab::Remapper::OverlapMesh, moab::TempestRemapper::OVERLAP_FILES );
rval = remapper.ConvertTempestMesh( moab::Remapper::OverlapMesh );MB_CHK_ERR( rval );
rval = mbCore->write_mesh( "tempest_intersection.h5m", &runCtx->meshsets[2], 1 );MB_CHK_ERR( rval );
// print verbosely about the problem setting
{
moab::Range rintxverts, rintxelems;
rval = mbCore->get_entities_by_dimension( runCtx->meshsets[0], 0, rintxverts );MB_CHK_ERR( rval );
rval = mbCore->get_entities_by_dimension( runCtx->meshsets[0], 2, rintxelems );MB_CHK_ERR( rval );
outputFormatter.printf( 0, "The source set contains %lu vertices and %lu elements \n", rintxverts.size(),
rintxelems.size() );
moab::Range bintxverts, bintxelems;
rval = mbCore->get_entities_by_dimension( runCtx->meshsets[1], 0, bintxverts );MB_CHK_ERR( rval );
rval = mbCore->get_entities_by_dimension( runCtx->meshsets[1], 2, bintxelems );MB_CHK_ERR( rval );
outputFormatter.printf( 0, "The target set contains %lu vertices and %lu elements \n", bintxverts.size(),
bintxelems.size() );
}
moab::EntityHandle intxset; // == remapper.GetMeshSet(moab::Remapper::OverlapMesh);
// Compute intersections with MOAB
{
// Create the intersection on the sphere object
runCtx->timer_push( "setup the intersector" );
moab::Intx2MeshOnSphere* mbintx = new moab::Intx2MeshOnSphere( mbCore );
mbintx->set_error_tolerance( epsrel );
mbintx->set_box_error( boxeps );
mbintx->set_radius_source_mesh( radius_src );
mbintx->set_radius_destination_mesh( radius_dest );
#ifdef MOAB_HAVE_MPI
mbintx->set_parallel_comm( pcomm );
#endif
rval = mbintx->FindMaxEdges( runCtx->meshsets[0], runCtx->meshsets[1] );MB_CHK_ERR( rval );
#ifdef MOAB_HAVE_MPI
moab::Range local_verts;
rval = mbintx->build_processor_euler_boxes( runCtx->meshsets[1], local_verts );MB_CHK_ERR( rval );
runCtx->timer_pop();
moab::EntityHandle covering_set;
runCtx->timer_push( "communicate the mesh" );
rval = mbintx->construct_covering_set( runCtx->meshsets[0], covering_set );MB_CHK_ERR( rval ); // lots of communication if mesh is distributed very differently
runCtx->timer_pop();
#else
moab::EntityHandle covering_set = runCtx->meshsets[0];
#endif
// Now let's invoke the MOAB intersection algorithm in parallel with a
// source and target mesh set representing two different decompositions
runCtx->timer_push( "compute intersections with MOAB" );
rval = mbCore->create_meshset( moab::MESHSET_SET, intxset );MB_CHK_SET_ERR( rval, "Can't create new set" );
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" );
runCtx->timer_pop();
// free the memory
delete mbintx;
}
{
moab::Range intxelems, intxverts;
rval = mbCore->get_entities_by_dimension( intxset, 2, intxelems );MB_CHK_ERR( rval );
rval = mbCore->get_entities_by_dimension( intxset, 0, intxverts, true );MB_CHK_ERR( rval );
outputFormatter.printf( 0, "The intersection set contains %lu elements and %lu vertices \n",
intxelems.size(), intxverts.size() );
double initial_sarea =
areaAdaptor.area_on_sphere( mbCore, runCtx->meshsets[0],
radius_src ); // use the target to compute the initial area
double initial_tarea =
areaAdaptor.area_on_sphere( mbCore, runCtx->meshsets[1],
radius_dest ); // use the target to compute the initial area
double intx_area = areaAdaptor.area_on_sphere( mbCore, intxset, radius_src );
outputFormatter.printf( 0, "mesh areas: source = %12.10f, target = %12.10f, intersection = %12.10f \n",
initial_sarea, initial_tarea, intx_area );
outputFormatter.printf( 0, "relative error w.r.t source = %12.10e, target = %12.10e \n",
fabs( intx_area - initial_sarea ) / initial_sarea,
fabs( intx_area - initial_tarea ) / initial_tarea );
}
// Write out our computed intersection file
rval = mbCore->write_mesh( "moab_intersection.h5m", &intxset, 1 );MB_CHK_ERR( rval );
if( runCtx->computeWeights )
{
runCtx->timer_push( "compute weights with the Tempest meshes" );
// Call to generate an offline map with the tempest meshes
OfflineMap weightMap;
int err = GenerateOfflineMapWithMeshes( *runCtx->meshes[0], *runCtx->meshes[1], *runCtx->meshes[2],
runCtx->disc_methods[0], // std::string strInputType
runCtx->disc_methods[1], // std::string strOutputType,
runCtx->mapOptions, weightMap );
runCtx->timer_pop();
std::map< std::string, std::string > mapAttributes;
if( err )
{
rval = moab::MB_FAILURE;
}
else
{
weightMap.Write( "outWeights.nc", mapAttributes );
}
}
}
else if( runCtx->meshType == moab::TempestRemapper::OVERLAP_MOAB )
{
// Usage: mpiexec -n 2 tools/mbtempest -t 5 -l mycs_2.h5m -l myico_2.h5m -f myoverlap_2.h5m
#ifdef MOAB_HAVE_MPI
rval = pcomm->check_all_shared_handles();MB_CHK_ERR( rval );
#endif
// print verbosely about the problem setting
{
moab::Range srcverts, srcelems;
rval = mbCore->get_entities_by_dimension( runCtx->meshsets[0], 0, srcverts );MB_CHK_ERR( rval );
rval = mbCore->get_entities_by_dimension( runCtx->meshsets[0], 2, srcelems );MB_CHK_ERR( rval );
rval = moab::IntxUtils::fix_degenerate_quads( mbCore, runCtx->meshsets[0] );MB_CHK_ERR( rval );
if( runCtx->enforceConvexity )
{
rval = moab::IntxUtils::enforce_convexity( mbCore, runCtx->meshsets[0], proc_id );MB_CHK_ERR( rval );
}
rval = areaAdaptor.positive_orientation( mbCore, runCtx->meshsets[0], radius_src );MB_CHK_ERR( rval );
if( !proc_id )
outputFormatter.printf( 0, "The source set contains %lu vertices and %lu elements \n",
srcverts.size(), srcelems.size() );
moab::Range tgtverts, tgtelems;
rval = mbCore->get_entities_by_dimension( runCtx->meshsets[1], 0, tgtverts );MB_CHK_ERR( rval );
rval = mbCore->get_entities_by_dimension( runCtx->meshsets[1], 2, tgtelems );MB_CHK_ERR( rval );
rval = moab::IntxUtils::fix_degenerate_quads( mbCore, runCtx->meshsets[1] );MB_CHK_ERR( rval );
if( runCtx->enforceConvexity )
{
rval = moab::IntxUtils::enforce_convexity( mbCore, runCtx->meshsets[1], proc_id );MB_CHK_ERR( rval );
}
rval = areaAdaptor.positive_orientation( mbCore, runCtx->meshsets[1], radius_dest );MB_CHK_ERR( rval );
if( !proc_id )
outputFormatter.printf( 0, "The target set contains %lu vertices and %lu elements \n",
tgtverts.size(), tgtelems.size() );
}
// First compute the covering set such that the target elements are fully covered by the
// lcoal source grid
runCtx->timer_push( "construct covering set for intersection" );
rval = remapper.ConstructCoveringSet( epsrel, 1.0, 1.0, boxeps, runCtx->rrmGrids );MB_CHK_ERR( rval );
runCtx->timer_pop();
// Compute intersections with MOAB with either the Kd-tree or the advancing front algorithm
runCtx->timer_push( "setup and compute mesh intersections" );
rval = remapper.ComputeOverlapMesh( runCtx->kdtreeSearch, false );MB_CHK_ERR( rval );
runCtx->timer_pop();
// print some diagnostic checks to see if the overlap grid resolved the input meshes
// correctly
{
double local_areas[3],
global_areas[3]; // Array for Initial area, and through Method 1 and Method 2
// local_areas[0] = area_on_sphere_lHuiller ( mbCore, runCtx->meshsets[1], radius_src );
local_areas[0] = areaAdaptor.area_on_sphere( mbCore, runCtx->meshsets[0], radius_src );
local_areas[1] = areaAdaptor.area_on_sphere( mbCore, runCtx->meshsets[1], radius_dest );
local_areas[2] = areaAdaptor.area_on_sphere( mbCore, runCtx->meshsets[2], radius_src );
#ifdef MOAB_HAVE_MPI
MPI_Allreduce( &local_areas[0], &global_areas[0], 3, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD );
#else
global_areas[0] = local_areas[0];
global_areas[1] = local_areas[1];
global_areas[2] = local_areas[2];
#endif
if( !proc_id )
{
outputFormatter.printf( 0,
"initial area: source mesh = %12.14f, target mesh = "
"%12.14f, overlap mesh = %12.14f\n",
global_areas[0], global_areas[1], global_areas[2] );
// outputFormatter.printf ( 0, " area with l'Huiller: %12.14f with Girard:
// %12.14f\n", global_areas[2], global_areas[3] ); outputFormatter.printf ( 0, "
// relative difference areas = %12.10e\n", fabs ( global_areas[2] - global_areas[3]
// ) / global_areas[2] );
outputFormatter.printf( 0, "relative error w.r.t source = %12.14e, and target = %12.14e\n",
fabs( global_areas[0] - global_areas[2] ) / global_areas[0],
fabs( global_areas[1] - global_areas[2] ) / global_areas[1] );
}
}
if( runCtx->intxFilename.size() )
{
moab::EntityHandle writableOverlapSet;
rval = mbCore->create_meshset( moab::MESHSET_SET, writableOverlapSet );MB_CHK_SET_ERR( rval, "Can't create new set" );
moab::EntityHandle meshOverlapSet = remapper.GetMeshSet( moab::Remapper::OverlapMesh );
moab::Range ovEnts;
rval = mbCore->get_entities_by_dimension( meshOverlapSet, 2, ovEnts );MB_CHK_SET_ERR( rval, "Can't create new set" );
rval = mbCore->get_entities_by_dimension( meshOverlapSet, 0, ovEnts );MB_CHK_SET_ERR( rval, "Can't create new set" );
#ifdef MOAB_HAVE_MPI
// Do not remove ghosted entities if we still haven't computed weights
// Remove ghosted entities from overlap set before writing the new mesh set to file
if( nprocs > 1 )
{
moab::Range ghostedEnts;
rval = remapper.GetOverlapAugmentedEntities( ghostedEnts );MB_CHK_ERR( rval );
ovEnts = moab::subtract( ovEnts, ghostedEnts );
}
#endif
rval = mbCore->add_entities( writableOverlapSet, ovEnts );MB_CHK_SET_ERR( rval, "Deleting ghosted entities failed" );
size_t lastindex = runCtx->intxFilename.find_last_of( "." );
sstr.str( "" );
sstr << runCtx->intxFilename.substr( 0, lastindex ) << ".h5m";
if( !runCtx->proc_id )
std::cout << "Writing out the MOAB intersection mesh file to " << sstr.str() << std::endl;
// Write out our computed intersection file
rval = mbCore->write_file( sstr.str().c_str(), NULL, writeOptions, &writableOverlapSet, 1 );MB_CHK_ERR( rval );
}
if( runCtx->computeWeights )
{
runCtx->meshes[2] = remapper.GetMesh( moab::Remapper::OverlapMesh );
if( !runCtx->proc_id ) std::cout << std::endl;
runCtx->timer_push( "setup computation of weights" );
// Call to generate the remapping weights with the tempest meshes
moab::TempestOnlineMap* weightMap = new moab::TempestOnlineMap( &remapper );
runCtx->timer_pop();
runCtx->timer_push( "compute weights with TempestRemap" );
rval = weightMap->GenerateRemappingWeights(
runCtx->disc_methods[0], // std::string strInputType
runCtx->disc_methods[1], // std::string strOutputType,
runCtx->mapOptions, // const GenerateOfflineMapAlgorithmOptions& options
runCtx->doftag_names[0], // const std::string& source_tag_name
runCtx->doftag_names[1] // const std::string& target_tag_name
);MB_CHK_ERR( rval );
runCtx->timer_pop();
// Invoke the CheckMap routine on the TempestRemap serial interface directly, if running
// on a single process
if( nprocs == 1 )
{
const double dNormalTolerance = 1.0E-8;
const double dStrictTolerance = 1.0E-12;
weightMap->CheckMap( runCtx->fCheck, runCtx->fCheck, runCtx->fCheck && ( runCtx->ensureMonotonicity ),
dNormalTolerance, dStrictTolerance );
}
if( runCtx->outFilename.size() )
{
// Write the map file to disk in parallel using either HDF5 or SCRIP interface
rval = weightMap->WriteParallelMap( runCtx->outFilename.c_str() );MB_CHK_ERR( rval );
// Write out the metadata information for the map file
if( proc_id == 0 )
{
size_t lastindex = runCtx->outFilename.find_last_of( "." );
sstr.str( "" );
sstr << runCtx->outFilename.substr( 0, lastindex ) << ".meta";
std::ofstream metafile( sstr.str() );
metafile << "Generator = MOAB-TempestRemap (mbtempest) Offline Regridding "
"Weight Generator"
<< std::endl;
metafile << "domain_a = " << runCtx->inFilenames[0] << std::endl;
metafile << "domain_b = " << runCtx->inFilenames[1] << std::endl;
metafile << "domain_aNb = "
<< ( runCtx->intxFilename.size() ? runCtx->intxFilename : "outOverlap.h5m" ) << std::endl;
metafile << "map_aPb = " << runCtx->outFilename << std::endl;
metafile << "type_src = " << runCtx->disc_methods[0] << std::endl;
metafile << "np_src = " << runCtx->disc_orders[0] << std::endl;
metafile << "concave_src = " << ( runCtx->mapOptions.fSourceConcave ? "true" : "false" )
<< std::endl;
metafile << "type_dst = " << runCtx->disc_methods[1] << std::endl;
metafile << "np_dst = " << runCtx->disc_orders[1] << std::endl;
metafile << "concave_dst = " << ( runCtx->mapOptions.fTargetConcave ? "true" : "false" )
<< std::endl;
metafile << "mono_type = " << runCtx->ensureMonotonicity << std::endl;
metafile << "bubble = " << ( runCtx->mapOptions.fNoBubble ? "false" : "true" ) << std::endl;
metafile << "version = "
<< "MOAB v5.1.0+" << std::endl;
metafile.close();
}
}
if( runCtx->verifyWeights )
{
// Let us pick a sampling test function for solution evaluation
moab::TempestOnlineMap::sample_function testFunction =
&sample_stationary_vortex; // &sample_slow_harmonic;
runCtx->timer_push( "describe a solution on source grid" );
moab::Tag srcAnalyticalFunction;
rval = weightMap->DefineAnalyticalSolution( srcAnalyticalFunction, "AnalyticalSolnSrcExact",
moab::Remapper::SourceMesh, testFunction );MB_CHK_ERR( rval );
runCtx->timer_pop();
// rval = mbCore->write_file ( "srcWithSolnTag.h5m", NULL, writeOptions,
// &runCtx->meshsets[0], 1 ); MB_CHK_ERR ( rval );
runCtx->timer_push( "describe a solution on target grid" );
moab::Tag tgtAnalyticalFunction;
moab::Tag tgtProjectedFunction;
rval = weightMap->DefineAnalyticalSolution( tgtAnalyticalFunction, "AnalyticalSolnTgtExact",
moab::Remapper::TargetMesh, testFunction,
&tgtProjectedFunction, "ProjectedSolnTgt" );MB_CHK_ERR( rval );
// rval = mbCore->write_file ( "tgtWithSolnTag.h5m", NULL, writeOptions,
// &runCtx->meshsets[1], 1 ); MB_CHK_ERR ( rval );
runCtx->timer_pop();
runCtx->timer_push( "compute solution projection on target grid" );
rval = weightMap->ApplyWeights( srcAnalyticalFunction, tgtProjectedFunction );MB_CHK_ERR( rval );
runCtx->timer_pop();
rval = mbCore->write_file( "tgtWithSolnTag2.h5m", NULL, writeOptions, &runCtx->meshsets[1], 1 );MB_CHK_ERR( rval );
if( nprocs == 1 && runCtx->baselineFile.size() )
{
// save the field from tgtWithSolnTag2 in a text file, and global ids for cells
moab::Range tgtCells;
rval = mbCore->get_entities_by_dimension( runCtx->meshsets[1], 2, tgtCells );MB_CHK_ERR( rval );
std::vector< int > globIds;
globIds.resize( tgtCells.size() );
std::vector< double > vals;
vals.resize( tgtCells.size() );
moab::Tag projTag;
rval = mbCore->tag_get_handle( "ProjectedSolnTgt", projTag );MB_CHK_ERR( rval );
moab::Tag gid = mbCore->globalId_tag();
rval = mbCore->tag_get_data( gid, tgtCells, &globIds[0] );MB_CHK_ERR( rval );
rval = mbCore->tag_get_data( projTag, tgtCells, &vals[0] );MB_CHK_ERR( rval );
std::fstream fs;
fs.open( runCtx->baselineFile.c_str(), std::fstream::out );
fs << std::setprecision( 15 ); // maximum precision for doubles
for( size_t i = 0; i < tgtCells.size(); i++ )
fs << globIds[i] << " " << vals[i] << "\n";
fs.close();
// for good measure, save the source file too, with the tag AnalyticalSolnSrcExact
// it will be used later to test, along with a target file
rval = mbCore->write_file( "srcWithSolnTag.h5m", NULL, writeOptions, &runCtx->meshsets[0], 1 );MB_CHK_ERR( rval );
}
runCtx->timer_push( "compute error metrics against analytical solution on target grid" );
std::map< std::string, double > errMetrics;
rval = weightMap->ComputeMetrics( moab::Remapper::TargetMesh, tgtAnalyticalFunction,
tgtProjectedFunction, errMetrics, true );MB_CHK_ERR( rval );
runCtx->timer_pop();
}
delete weightMap;
}
}
// Clean up
remapper.clear();
delete runCtx;
delete mbCore;
#ifdef MOAB_HAVE_MPI
MPI_Finalize();
#endif
exit( 0 );
}
| double sample_constant | ( | double | dLon, |
| double | dLat | ||
| ) | [inline] |
Definition at line 960 of file mbtempest.cpp.
{
return 1.0;
}
| double sample_fast_harmonic | ( | double | dLon, |
| double | dLat | ||
| ) | [inline] |
Definition at line 954 of file mbtempest.cpp.
{
return ( 2.0 + pow( sin( 2.0 * dLat ), 16.0 ) * cos( 16.0 * dLon ) );
// return (2.0 + pow(cos(2.0 * dLat), 16.0) * cos(16.0 * dLon));
}
| double sample_slow_harmonic | ( | double | dLon, |
| double | dLat | ||
| ) | [inline] |
Definition at line 949 of file mbtempest.cpp.
{
return ( 2.0 + cos( dLat ) * cos( dLat ) * cos( 2.0 * dLon ) );
}
| 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 965 of file mbtempest.cpp.
Referenced by main().
{
const double dLon0 = 0.0;
const double dLat0 = 0.6;
const double dR0 = 3.0;
const double dD = 5.0;
const double dT = 6.0;
/// Find the rotated longitude and latitude of a point on a sphere
/// with pole at (dLonC, dLatC).
{
double dSinC = sin( dLat0 );
double dCosC = cos( dLat0 );
double dCosT = cos( dLat );
double dSinT = sin( dLat );
double dTrm = dCosT * cos( dLon - dLon0 );
double dX = dSinC * dTrm - dCosC * dSinT;
double dY = dCosT * sin( dLon - dLon0 );
double dZ = dSinC * dSinT + dCosC * dTrm;
dLon = atan2( dY, dX );
if( dLon < 0.0 )
{
dLon += 2.0 * M_PI;
}
dLat = asin( dZ );
}
double dRho = dR0 * cos( dLat );
double dVt = 3.0 * sqrt( 3.0 ) / 2.0 / cosh( dRho ) / cosh( dRho ) * tanh( dRho );
double dOmega;
if( dRho == 0.0 )
{
dOmega = 0.0;
}
else
{
dOmega = dVt / dRho;
}
return ( 1.0 - tanh( dRho / dD * sin( dLon - dOmega * dT ) ) );
}
1.7.6.1