moab
|
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 }