cslam_par_test.cpp

Go to the documentation of this file.

/*
 * cslam_par_test.cpp
 * test to trigger intersection on a sphere in parallel
 * it will start from an eulerian mesh (mesh + velocity from Mark Taylor)
 * file: VELO00.h5m; mesh + velo at 850 milibar, in an unstructured mesh refined from
 * a cube sphere grid
 * the mesh is read in parallel (euler mesh);
 *
 * lagrangian mesh is obtained using
 * pos (t-dt) = pos(t) -Velo(t)*dt
 * then, project on the sphere with a given radius
 *
 * Created on: Apr 22, 2013
 */
 
#include <iostream>
#include <sstream>
#include <ctime>
#include <cstdlib>
#include <cstdio>
#include <cstring>
#include "moab/Core.hpp"
#include "moab/Interface.hpp"
#include "moab/IntxMesh/Intx2MeshOnSphere.hpp"
#include <cmath>
#include "TestUtil.hpp"
#include "moab/ParallelComm.hpp"
#include "moab/ProgOptions.hpp"
#include "MBParallelConventions.h"
#include "moab/ReadUtilIface.hpp"
#include "MBTagConventions.hpp"
 
#include "moab/IntxMesh/IntxUtils.hpp"
#include "IntxUtilsCSLAM.hpp"
 
// for M_PI
#include <cmath>
 
using namespace moab;
// some input data
double EPS1 = 0.2; // this is for box error
std::string input_mesh_file( "VELO00_16p.h5m" ); // input file, partitioned correctly
double Radius = 1.0; // change to radius
double deltaT = 1.e-6;
void test_intx_in_parallel_elem_based();
 
int main( int argc, char** argv )
{
 MPI_Init( &argc, &argv );
 EPS1 = 0.000002;
 int result = 0;
 
 if( argc > 1 )
 {
 int index = 1;
 while( index < argc )
 {
 if( !strcmp( argv[index], "-eps" ) ) // this is for box error
 {
 EPS1 = atof( argv[++index] );
 }
 if( !strcmp( argv[index], "-input" ) )
 {
 input_mesh_file = argv[++index];
 }
 if( !strcmp( argv[index], "-radius" ) )
 {
 Radius = atof( argv[++index] );
 }
 if( !strcmp( argv[index], "-deltaT" ) )
 {
 deltaT = atof( argv[++index] );
 }
 index++;
 }
 }
 std::cout << " run: -input " << input_mesh_file << " -eps " << EPS1 << " -radius " << Radius << " -deltaT "
 << deltaT << "\n";
 
 result += RUN_TEST( test_intx_in_parallel_elem_based );
 
 MPI_Finalize();
 return result;
}
// will save the LOC tag on the euler nodes
ErrorCode compute_lagrange_mesh_on_sphere( Interface* mb, EntityHandle euler_set )
{
 /*
 * get all quads first, then vertices, then move them on the surface of the sphere
 * radius is 1, usually
 * pos (t-dt) = pos(t) -Velo(t)*dt; this will be lagrange mesh, on each processor
 */
 Range quads;
 ErrorCode rval = mb->get_entities_by_type( euler_set, MBQUAD, quads );MB_CHK_ERR( rval );
 
 Range connecVerts;
 rval = mb->get_connectivity( quads, connecVerts );MB_CHK_ERR( rval );
 
 // the LOC tag, should be provided by the user?
 Tag tagh = 0;
 std::string tag_name( "DP" );
 rval = mb->tag_get_handle( tag_name.c_str(), 3, MB_TYPE_DOUBLE, tagh, MB_TAG_DENSE | MB_TAG_CREAT );MB_CHK_ERR( rval );
 void* data; // pointer to the DP in memory, for each vertex
 int count;
 
 rval = mb->tag_iterate( tagh, connecVerts.begin(), connecVerts.end(), count, data );MB_CHK_ERR( rval );
 // here we are checking contiguity
 assert( count == (int)connecVerts.size() );
 double* ptr_DP = (double*)data;
 // get the coordinates of the old mesh, and move it around using velocity tag
 
 Tag tagv = 0;
 std::string velo_tag_name( "VELO" );
 rval = mb->tag_get_handle( velo_tag_name.c_str(), 3, MB_TYPE_DOUBLE, tagv, MB_TAG_DENSE );MB_CHK_ERR( rval );
 
 /*void *datavelo; // pointer to the VELO in memory, for each vertex
 
 rval = mb->tag_iterate(tagv, connecVerts.begin(), connecVerts.end(), count, datavelo);MB_CHK_ERR(rval);*/
 // here we are checking contiguity
 assert( count == (int)connecVerts.size() );
 // now put the vertices in the right place....
 // int vix=0; // vertex index in new array
 
 for( Range::iterator vit = connecVerts.begin(); vit != connecVerts.end(); ++vit )
 {
 EntityHandle oldV = *vit;
 CartVect posi;
 rval = mb->get_coords( &oldV, 1, &( posi[0] ) );MB_CHK_ERR( rval );
 CartVect velo;
 rval = mb->tag_get_data( tagv, &oldV, 1, (void*)&( velo[0] ) );MB_CHK_ERR( rval );
 // do some mumbo jumbo, as in python script
 CartVect newPos = posi - deltaT * velo;
 double len1 = newPos.length();
 newPos = Radius * newPos / len1;
 
 ptr_DP[0] = newPos[0];
 ptr_DP[1] = newPos[1];
 ptr_DP[2] = newPos[2];
 ptr_DP += 3; // increment to the next node
 }
 
 return rval;
}
 
void test_intx_in_parallel_elem_based()
{
 std::string opts = std::string( "PARALLEL=READ_PART;PARTITION=PARALLEL_PARTITION" ) +
 std::string( ";PARALLEL_RESOLVE_SHARED_ENTS" );
 Core moab;
 Interface& mb = moab;
 EntityHandle euler_set;
 ErrorCode rval;
 rval = mb.create_meshset( MESHSET_SET, euler_set );MB_CHK_ERR_RET( rval );
 std::string example( TestDir + "unittest/" + input_mesh_file );
 
 rval = mb.load_file( example.c_str(), &euler_set, opts.c_str() );
 
 ParallelComm* pcomm = ParallelComm::get_pcomm( &mb, 0 );MB_CHK_ERR_RET( rval );
 
 rval = pcomm->check_all_shared_handles();MB_CHK_ERR_RET( rval );
 
 // everybody will get a DP tag, including the non owned entities; so exchange tags is not
 // required for LOC (here)
 rval = compute_lagrange_mesh_on_sphere( &mb, euler_set );MB_CHK_ERR_RET( rval );
 
 int rank = pcomm->proc_config().proc_rank();
 
 std::stringstream ste;
 ste << "initial" << rank << ".vtk";
 mb.write_file( ste.str().c_str(), 0, 0, &euler_set, 1 );
 
 Intx2MeshOnSphere worker( &mb );
 
 worker.set_radius_source_mesh( Radius );
 worker.set_radius_destination_mesh( Radius );
 worker.set_box_error( EPS1 ); //
 // worker.SetEntityType(MBQUAD);
 
 worker.set_error_tolerance( Radius * 1.e-8 );
 std::cout << "error tolerance epsilon_1=" << Radius * 1.e-8 << "\n";
 // worker.locate_departure_points(euler_set);
 
 rval = worker.FindMaxEdges( euler_set, euler_set ); // departure will be the same max_edges
 // we need to make sure the covering set is bigger than the euler mesh
 EntityHandle covering_lagr_set;
 rval = mb.create_meshset( MESHSET_SET, covering_lagr_set );MB_CHK_ERR_RET( rval );
 
 rval = worker.create_departure_mesh_2nd_alg( euler_set, covering_lagr_set );MB_CHK_ERR_RET( rval );
 
 std::stringstream ss;
 ss << "partial" << rank << ".vtk";
 mb.write_file( ss.str().c_str(), 0, 0, &covering_lagr_set, 1 );
 EntityHandle outputSet;
 rval = mb.create_meshset( MESHSET_SET, outputSet );MB_CHK_ERR_RET( rval );
 rval = worker.intersect_meshes( covering_lagr_set, euler_set, outputSet );MB_CHK_ERR_RET( rval );
 
 // std::string opts_write("PARALLEL=WRITE_PART");
 // rval = mb.write_file("manuf.h5m", 0, opts_write.c_str(), &outputSet, 1);
 // std::string opts_write("");
 std::stringstream outf;
 outf << "intersect" << rank << ".h5m";
 rval = mb.write_file( outf.str().c_str(), 0, 0, &outputSet, 1 );MB_CHK_ERR_RET( rval );
 
 moab::IntxAreaUtils sphAreaUtils;
 double intx_area = sphAreaUtils.area_on_sphere( &mb, outputSet, Radius );
 double arrival_area = sphAreaUtils.area_on_sphere( &mb, euler_set, Radius );
 std::cout << "On rank : " << rank << " arrival area: " << arrival_area << " intersection area:" << intx_area
 << " rel error: " << fabs( ( intx_area - arrival_area ) / arrival_area ) << "\n";
}