moab
cslam_par_test.cpp
Go to the documentation of this file.
00001 /*
00002  * cslam_par_test.cpp
00003  *  test to trigger intersection on a sphere in parallel
00004  *  it will start from an eulerian mesh (mesh + velocity from Mark Taylor)
00005  *   file: VELO00.h5m; mesh + velo at 850 milibar, in an unstructured mesh refined from
00006  *   a cube sphere grid
00007  *  the mesh is read in parallel (euler mesh);
00008  *
00009  *   lagrangian mesh is obtained using
00010  *   pos (t-dt) = pos(t) -Velo(t)*dt
00011  *   then, project on the sphere with a given radius
00012  *
00013  *  Created on: Apr 22, 2013
00014  */
00015 
00016 #include <iostream>
00017 #include <sstream>
00018 #include <time.h>
00019 #include <stdlib.h>
00020 #include <stdio.h>
00021 #include <string.h>
00022 #include "moab/Core.hpp"
00023 #include "moab/Interface.hpp"
00024 #include "Intx2MeshOnSphere.hpp"
00025 #include <math.h>
00026 #include "TestUtil.hpp"
00027 #include "moab/ParallelComm.hpp"
00028 #include "moab/ProgOptions.hpp"
00029 #include "MBParallelConventions.h"
00030 #include "moab/ReadUtilIface.hpp"
00031 #include "MBTagConventions.hpp"
00032 
00033 #include "CslamUtils.hpp"
00034 
00035 // for M_PI
00036 #include <math.h>
00037 
00038 #ifdef MESHDIR
00039 std::string TestDir( STRINGIFY(MESHDIR) );
00040 #else
00041 std::string TestDir(".");
00042 #endif
00043 
00044 using namespace moab;
00045 // some input data
00046 double EPS1=0.2; // this is for box error
00047 std::string input_mesh_file("VELO00_16p.h5m"); // input file, partitioned correctly
00048 double Radius = 1.0; // change to radius
00049 double deltaT = 1.e-6;
00050 void test_intx_in_parallel_elem_based();
00051 
00052 int main(int argc, char **argv)
00053 {
00054   MPI_Init(&argc, &argv);
00055   EPS1 = 0.000002;
00056   int result = 0;
00057 
00058   if (argc>1)
00059   {
00060     int index=1;
00061     while (index<argc)
00062     {
00063       if (!strcmp( argv[index], "-eps")) // this is for box error
00064       {
00065         EPS1=atof(argv[++index]);
00066       }
00067       if (!strcmp( argv[index], "-input"))
00068       {
00069         input_mesh_file=argv[++index];
00070       }
00071       if (!strcmp( argv[index], "-radius"))
00072       {
00073         Radius=atof(argv[++index]);
00074       }
00075       if (!strcmp( argv[index], "-deltaT"))
00076       {
00077         deltaT=atof(argv[++index]);
00078       }
00079       index++;
00080     }
00081   }
00082   std::cout << " run: -input " << input_mesh_file << "  -eps " << EPS1 <<
00083       " -radius " << Radius << " -deltaT " << deltaT << "\n";
00084 
00085   result += RUN_TEST(test_intx_in_parallel_elem_based);
00086 
00087   MPI_Finalize();
00088   return result;
00089 }
00090 // will save the LOC tag on the euler nodes
00091 ErrorCode  compute_lagrange_mesh_on_sphere(Interface * mb, EntityHandle euler_set)
00092 {
00093   ErrorCode rval = MB_SUCCESS;
00094 
00095   /*
00096    * get all quads first, then vertices, then move them on the surface of the sphere
00097    *  radius is 1, usually
00098    *  pos (t-dt) = pos(t) -Velo(t)*dt; this will be lagrange mesh, on each processor
00099    */
00100   Range quads;
00101   rval = mb->get_entities_by_type(euler_set, MBQUAD, quads);
00102   CHECK_ERR(rval);
00103 
00104   Range connecVerts;
00105   rval = mb->get_connectivity(quads, connecVerts);
00106 
00107   // the LOC tag, should be provided by the user?
00108   Tag tagh = 0;
00109   std::string tag_name("DP");
00110   rval = mb->tag_get_handle(tag_name.c_str(), 3, MB_TYPE_DOUBLE, tagh, MB_TAG_DENSE | MB_TAG_CREAT);
00111   CHECK_ERR(rval);
00112   void *data; // pointer to the DP in memory, for each vertex
00113   int count;
00114 
00115   rval = mb->tag_iterate(tagh, connecVerts.begin(), connecVerts.end(), count, data);
00116   CHECK_ERR(rval);
00117   // here we are checking contiguity
00118   assert(count == (int) connecVerts.size());
00119   double * ptr_DP=(double*)data;
00120   // get the coordinates of the old mesh, and move it around using velocity tag
00121 
00122   Tag tagv = 0;
00123   std::string velo_tag_name("VELO");
00124   rval = mb->tag_get_handle(velo_tag_name.c_str(), 3, MB_TYPE_DOUBLE, tagv, MB_TAG_DENSE);
00125   CHECK_ERR(rval);
00126 
00127   /*void *datavelo; // pointer to the VELO in memory, for each vertex
00128 
00129   rval = mb->tag_iterate(tagv, connecVerts.begin(), connecVerts.end(), count, datavelo);
00130   CHECK_ERR(rval);*/
00131   // here we are checking contiguity
00132   assert(count == (int) connecVerts.size());
00133 // now put the vertices in the right place....
00134   //int vix=0; // vertex index in new array
00135 
00136   for (Range::iterator vit=connecVerts.begin();vit!=connecVerts.end(); vit++ )
00137   {
00138     EntityHandle oldV=*vit;
00139     CartVect posi;
00140     rval = mb->get_coords(&oldV, 1, &(posi[0]) );
00141     CHECK_ERR(rval);
00142     CartVect velo;
00143     rval = mb->tag_get_data(tagv, &oldV, 1, (void *) &(velo[0]) );
00144     CHECK_ERR(rval);
00145     // do some mumbo jumbo, as in python script
00146     CartVect newPos = posi - deltaT*velo;
00147     double len1= newPos.length();
00148     newPos = Radius*newPos/len1;
00149 
00150     ptr_DP[0]=newPos[0];
00151     ptr_DP[1]=newPos[1];
00152     ptr_DP[2]=newPos[2];
00153     ptr_DP+=3; // increment to the next node
00154   }
00155 
00156   return rval;
00157 }
00158 
00159 void test_intx_in_parallel_elem_based()
00160 {
00161   std::string opts = std::string("PARALLEL=READ_PART;PARTITION=PARALLEL_PARTITION")+
00162         std::string(";PARALLEL_RESOLVE_SHARED_ENTS");
00163   Core moab;
00164   Interface & mb = moab;
00165   EntityHandle euler_set;
00166   ErrorCode rval;
00167   rval = mb.create_meshset(MESHSET_SET, euler_set);
00168   CHECK_ERR(rval);
00169   std::string example(TestDir + "/" +  input_mesh_file);
00170 
00171   rval = mb.load_file(example.c_str(), &euler_set, opts.c_str());
00172 
00173   ParallelComm* pcomm = ParallelComm::get_pcomm(&mb, 0);
00174   CHECK_ERR(rval);
00175 
00176   rval = pcomm->check_all_shared_handles();
00177   CHECK_ERR(rval);
00178 
00179   // everybody will get a DP tag, including the non owned entities; so exchange tags is not required for LOC (here)
00180   rval = compute_lagrange_mesh_on_sphere(&mb, euler_set);
00181   CHECK_ERR(rval);
00182 
00183   int rank = pcomm->proc_config().proc_rank();
00184 
00185   std::stringstream ste;
00186   ste<<"initial" << rank<<".vtk";
00187   mb.write_file(ste.str().c_str(), 0, 0, &euler_set, 1);
00188 
00189   Intx2MeshOnSphere worker(&mb);
00190 
00191   worker.SetRadius(Radius);
00192   worker.set_box_error(EPS1);//
00193   //worker.SetEntityType(MBQUAD);
00194 
00195   worker.SetErrorTolerance(Radius*1.e-8);
00196   std::cout << "error tolerance epsilon_1="<< Radius*1.e-8 << "\n";
00197   //  worker.locate_departure_points(euler_set);
00198 
00199   // we need to make sure the covering set is bigger than the euler mesh
00200   EntityHandle covering_lagr_set;
00201   rval = mb.create_meshset(MESHSET_SET, covering_lagr_set);
00202   CHECK_ERR(rval);
00203 
00204   rval = worker.create_departure_mesh_2nd_alg(euler_set, covering_lagr_set);
00205   CHECK_ERR(rval);
00206 
00207   std::stringstream ss;
00208   ss<<"partial" << rank<<".vtk";
00209   mb.write_file(ss.str().c_str(), 0, 0, &covering_lagr_set, 1);
00210   EntityHandle outputSet;
00211   rval = mb.create_meshset(MESHSET_SET, outputSet);
00212   CHECK_ERR(rval);
00213   rval = worker.intersect_meshes(covering_lagr_set, euler_set, outputSet);
00214   CHECK_ERR(rval);
00215 
00216   //std::string opts_write("PARALLEL=WRITE_PART");
00217   //rval = mb.write_file("manuf.h5m", 0, opts_write.c_str(), &outputSet, 1);
00218   std::string opts_write("");
00219   std::stringstream outf;
00220   outf<<"intersect" << rank<<".h5m";
00221   rval = mb.write_file(outf.str().c_str(), 0, 0, &outputSet, 1);
00222   double intx_area = area_on_sphere(&mb, outputSet, Radius);
00223   double arrival_area = area_on_sphere(&mb, euler_set, Radius) ;
00224   std::cout<< "On rank : " << rank << " arrival area: " << arrival_area<<
00225       "  intersection area:" << intx_area << " rel error: " << fabs((intx_area-arrival_area)/arrival_area) << "\n";
00226   CHECK_ERR(rval);
00227 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines