moab
|
00001 /* 00002 * diffusion.cpp 00003 * 00004 * Created on: Aug 12, 2013 00005 * 00006 */ 00007 00008 /* trigger a diffusion computation in serial, first; 00009 we will start from an unstructured mesh on a sphere; 00010 do case 1: give some tracers a "slotted cylinder" ; compute the initial concentrations 00011 along a slotted cylinder, and see how it varies in time; 00012 use formula from Nair & Lauritzen paper: 00013 A class of deformational flow test cases for linear transport problems 00014 on the sphere; see CSLAM Utils case1 00015 scalar fields are defined at page 4-5 in the paper 00016 00017 steps: 00018 first create an initial tracer field, and save it as field on a sphere 00019 (initial time step 0) 00020 then use velocities (non-divergent, first), to transport the tracer on 00021 the sphere, according to some flow fields 00022 */ 00023 // copy from intx_mpas test 00024 00025 #include <iostream> 00026 #include <sstream> 00027 #include <time.h> 00028 #include <stdlib.h> 00029 #include <stdio.h> 00030 #include <string.h> 00031 #include "moab/Core.hpp" 00032 #include "moab/Interface.hpp" 00033 #include "Intx2MeshOnSphere.hpp" 00034 #include <math.h> 00035 #include "moab/ProgOptions.hpp" 00036 #include "MBTagConventions.hpp" 00037 #include "TestUtil.hpp" 00038 #include "moab/ParallelComm.hpp" 00039 00040 #include "CslamUtils.hpp" 00041 00042 const char BRIEF_DESC[] = 00043 "Simulate a transport problem in a semi-Lagrangian formulation\n"; 00044 std::ostringstream LONG_DESC; 00045 00046 // non smooth scalar field 00047 // some input data 00048 double gtol = 1.e-9; // this is for geometry tolerance 00049 00050 double radius = 1.;// in m: 6371220. 00051 00052 int numSteps = 50; // number of times with velocity displayed at points 00053 double T = 5; 00054 00055 int case_number = 1; // 1, 2 (non-divergent) 3 divergent 00056 00057 moab::Tag corrTag; 00058 bool writeFiles = false; 00059 bool parallelWrite = false; 00060 bool velocity = false; 00061 int field_type = 1 ; // 1 quasi smooth, 2 - smooth, 3 non-smooth, 00062 #ifdef MESHDIR 00063 std::string TestDir( STRINGIFY(MESHDIR) ); 00064 #else 00065 std::string TestDir("."); 00066 #endif 00067 00068 // for M_PI 00069 #include <math.h> 00070 00071 #define STRINGIFY_(X) #X 00072 #define STRINGIFY(X) STRINGIFY_(X) 00073 00074 using namespace moab; 00075 ErrorCode add_field_value(Interface * mb, EntityHandle euler_set, int rank, Tag & tagTracer, Tag & tagElem, Tag & tagArea) 00076 { 00077 ErrorCode rval = MB_SUCCESS; 00078 00079 /* 00080 * get all plys first, then vertices, then move them on the surface of the sphere 00081 * radius is 1., most of the time 00082 * 00083 */ 00084 Range polygons; 00085 rval = mb->get_entities_by_dimension(euler_set, 2, polygons); 00086 if (MB_SUCCESS != rval) 00087 return rval; 00088 00089 Range connecVerts; 00090 rval = mb->get_connectivity(polygons, connecVerts); 00091 if (MB_SUCCESS != rval) 00092 return rval; 00093 00094 00095 00096 void *data; // pointer to the LOC in memory, for each vertex 00097 int count; 00098 00099 rval = mb->tag_iterate(tagTracer, connecVerts.begin(), connecVerts.end(), count, data); 00100 CHECK_ERR(rval); 00101 // here we are checking contiguity 00102 assert(count == (int) connecVerts.size()); 00103 double * ptr_DP=(double*)data; 00104 // lambda is for longitude, theta for latitude 00105 // param will be: (la1, te1), (la2, te2), b, c; hmax=1, r=1/2 00106 // nondivergent flow, page 5, case 1, (la1, te1) = (M_PI, M_PI/3) 00107 // (la2, te2) = (M_PI, -M_PI/3) 00108 // la1, te1 la2 te2 b c hmax r 00109 if (field_type==1) // quasi smooth 00110 { 00111 double params[] = { M_PI, M_PI/3, M_PI, -M_PI/3, 0.1, 0.9, 1., 0.5}; 00112 for (Range::iterator vit=connecVerts.begin();vit!=connecVerts.end(); vit++ ) 00113 { 00114 EntityHandle oldV=*vit; 00115 CartVect posi; 00116 rval = mb->get_coords(&oldV, 1, &(posi[0]) ); 00117 CHECK_ERR(rval); 00118 00119 SphereCoords sphCoord = cart_to_spherical(posi); 00120 00121 ptr_DP[0]=quasi_smooth_field(sphCoord.lon, sphCoord.lat, params);; 00122 00123 ptr_DP++; // increment to the next node 00124 } 00125 } 00126 else if (2 == field_type) // smooth 00127 { 00128 CartVect p1, p2; 00129 SphereCoords spr; 00130 spr.R = 1; 00131 spr.lat = M_PI/3; 00132 spr.lon= M_PI; 00133 p1 = spherical_to_cart(spr); 00134 spr.lat = -M_PI/3; 00135 p2 = spherical_to_cart(spr); 00136 // x1, y1, z1, x2, y2, z2, h_max, b0 00137 double params[] = { p1[0], p1[1], p1[2], p2[0], p2[1], p2[2], 1, 5.}; 00138 for (Range::iterator vit=connecVerts.begin();vit!=connecVerts.end(); vit++ ) 00139 { 00140 EntityHandle oldV=*vit; 00141 CartVect posi; 00142 rval = mb->get_coords(&oldV, 1, &(posi[0]) ); 00143 CHECK_ERR(rval); 00144 00145 SphereCoords sphCoord = cart_to_spherical(posi); 00146 00147 ptr_DP[0]=smooth_field(sphCoord.lon, sphCoord.lat, params);; 00148 00149 ptr_DP++; // increment to the next node 00150 } 00151 } 00152 else if (3 == field_type) // slotted 00153 { 00154 // la1, te1, la2, te2, b, c, r 00155 double params[] = { M_PI, M_PI/3, M_PI, -M_PI/3, 0.1, 0.9, 0.5};// no h_max 00156 for (Range::iterator vit=connecVerts.begin();vit!=connecVerts.end(); vit++ ) 00157 { 00158 EntityHandle oldV=*vit; 00159 CartVect posi; 00160 rval = mb->get_coords(&oldV, 1, &(posi[0]) ); 00161 CHECK_ERR(rval); 00162 00163 SphereCoords sphCoord = cart_to_spherical(posi); 00164 00165 ptr_DP[0]=slotted_cylinder_field(sphCoord.lon, sphCoord.lat, params);; 00166 00167 ptr_DP++; // increment to the next node 00168 } 00169 } 00170 00171 // add average value for quad/polygon (average corners) 00172 // do some averages 00173 00174 00175 Range::iterator iter = polygons.begin(); 00176 double local_mass = 0.; // this is total mass on one proc 00177 while (iter != polygons.end()) 00178 { 00179 rval = mb->tag_iterate(tagElem, iter, polygons.end(), count, data); 00180 CHECK_ERR(rval); 00181 double * ptr=(double*)data; 00182 00183 rval = mb->tag_iterate(tagArea, iter, polygons.end(), count, data); 00184 CHECK_ERR(rval); 00185 double * ptrArea=(double*)data; 00186 for (int i=0; i<count; i++, iter++, ptr++, ptrArea++) 00187 { 00188 const moab::EntityHandle * conn = NULL; 00189 int num_nodes = 0; 00190 rval = mb->get_connectivity(*iter, conn, num_nodes); 00191 CHECK_ERR(rval); 00192 if (num_nodes==0) 00193 return MB_FAILURE; 00194 std::vector<double> nodeVals(num_nodes); 00195 double average=0.; 00196 rval = mb->tag_get_data(tagTracer, conn, num_nodes, &nodeVals[0] ); 00197 CHECK_ERR(rval); 00198 for (int j=0; j<num_nodes; j++) 00199 average+=nodeVals[j]; 00200 average/=num_nodes; 00201 *ptr = average; 00202 00203 // now get area 00204 std::vector<double> coords; 00205 coords.resize(3*num_nodes); 00206 rval = mb->get_coords(conn, num_nodes, &coords[0]); 00207 CHECK_ERR(rval); 00208 *ptrArea = area_spherical_polygon_lHuiller (&coords[0], num_nodes, radius); 00209 00210 // we should have used some 00211 // total mass: 00212 local_mass += *ptrArea * average; 00213 } 00214 00215 } 00216 double total_mass=0.; 00217 int mpi_err = MPI_Reduce(&local_mass, &total_mass, 1, MPI_DOUBLE, MPI_SUM, 0, MPI_COMM_WORLD); 00218 if (MPI_SUCCESS != mpi_err) return MB_FAILURE; 00219 00220 if (rank==0) 00221 std::cout << "initial total mass:" << total_mass << "\n"; 00222 00223 // now we can delete the tags? not yet 00224 return MB_SUCCESS; 00225 } 00226 00227 ErrorCode compute_velocity_case1(Interface * mb, EntityHandle euler_set, Tag & tagh, int rank, int tStep) 00228 { 00229 ErrorCode rval = MB_SUCCESS; 00230 00231 Range polygons; 00232 rval = mb->get_entities_by_dimension(euler_set, 2, polygons); 00233 if (MB_SUCCESS != rval) 00234 return rval; 00235 00236 Range connecVerts; 00237 rval = mb->get_connectivity(polygons, connecVerts); 00238 if (MB_SUCCESS != rval) 00239 return rval; 00240 00241 void *data; // pointer to the velo in memory, for each vertex 00242 int count; 00243 00244 rval = mb->tag_iterate(tagh, connecVerts.begin(), connecVerts.end(), count, data); 00245 CHECK_ERR(rval); 00246 // here we are checking contiguity 00247 assert(count == (int) connecVerts.size()); 00248 double * ptr_velo=(double*)data; 00249 // lambda is for longitude, theta for latitude 00250 00251 for (Range::iterator vit=connecVerts.begin();vit!=connecVerts.end(); vit++ ) 00252 { 00253 EntityHandle oldV=*vit; 00254 CartVect posi; 00255 rval = mb->get_coords(&oldV, 1, &(posi[0]) ); 00256 CHECK_ERR(rval); 00257 CartVect velo ; 00258 double t = T * tStep/numSteps; // 00259 velocity_case1(posi, t, velo); 00260 00261 ptr_velo[0]= velo[0]; 00262 ptr_velo[1]= velo[1]; 00263 ptr_velo[2]= velo[2]; 00264 00265 // increment to the next node 00266 ptr_velo+=3;// to next velocity 00267 } 00268 00269 std::stringstream velos; 00270 velos<<"Tracer" << rank<<"_"<<tStep<< ".vtk"; 00271 rval = mb->write_file(velos.str().c_str(), 0, 0, &euler_set, 1, &tagh, 1); 00272 CHECK_ERR(rval); 00273 00274 00275 return MB_SUCCESS; 00276 } 00277 ErrorCode create_lagr_mesh(Interface * mb, EntityHandle euler_set, EntityHandle lagr_set) 00278 { 00279 // create the handle tag for the corresponding element / vertex 00280 00281 EntityHandle dum = 0; 00282 00283 ErrorCode rval = mb->tag_get_handle(CORRTAGNAME, 00284 1, MB_TYPE_HANDLE, corrTag, 00285 MB_TAG_DENSE|MB_TAG_CREAT, &dum); 00286 CHECK_ERR(rval); 00287 00288 // give the same global id to new verts and cells created in the lagr(departure) mesh 00289 Tag gid; 00290 rval = mb->tag_get_handle(GLOBAL_ID_TAG_NAME, 1, MB_TYPE_INTEGER, gid, MB_TAG_DENSE); 00291 CHECK_ERR(rval); 00292 00293 Range polys; 00294 rval = mb->get_entities_by_dimension(euler_set, 2, polys); 00295 CHECK_ERR(rval); 00296 00297 Range connecVerts; 00298 rval = mb->get_connectivity(polys, connecVerts); 00299 CHECK_ERR(rval); 00300 00301 std::map<EntityHandle, EntityHandle> newNodes; 00302 for (Range::iterator vit = connecVerts.begin(); vit != connecVerts.end(); vit++) 00303 { 00304 EntityHandle oldV = *vit; 00305 CartVect posi; 00306 rval = mb->get_coords(&oldV, 1, &(posi[0])); 00307 CHECK_ERR(rval); 00308 int global_id; 00309 rval = mb->tag_get_data(gid, &oldV, 1, &global_id); 00310 CHECK_ERR(rval); 00311 EntityHandle new_vert; 00312 rval = mb->create_vertex(&(posi[0]), new_vert); // duplicate the position 00313 CHECK_ERR(rval); 00314 newNodes[oldV] = new_vert; 00315 // set also the correspondent tag :) 00316 rval = mb->tag_set_data(corrTag, &oldV, 1, &new_vert); 00317 CHECK_ERR(rval); 00318 // also the other side 00319 // need to check if we really need this; the new vertex will never need the old vertex 00320 // we have the global id which is the same 00321 /*rval = mb->tag_set_data(corrTag, &new_vert, 1, &oldV); 00322 CHECK_ERR(rval);*/ 00323 // set the global id on the corresponding vertex the same as the initial vertex 00324 rval = mb->tag_set_data(gid, &new_vert, 1, &global_id); 00325 CHECK_ERR(rval); 00326 00327 } 00328 for (Range::iterator it = polys.begin(); it != polys.end(); it++) 00329 { 00330 EntityHandle q = *it; 00331 int nnodes; 00332 const EntityHandle * conn; 00333 rval = mb->get_connectivity(q, conn, nnodes); 00334 CHECK_ERR(rval); 00335 int global_id; 00336 rval = mb->tag_get_data(gid, &q, 1, &global_id); 00337 CHECK_ERR(rval); 00338 EntityType typeElem = mb->type_from_handle(q); 00339 std::vector<EntityHandle> new_conn(nnodes); 00340 for (int i = 0; i < nnodes; i++) 00341 { 00342 EntityHandle v1 = conn[i]; 00343 new_conn[i] = newNodes[v1]; 00344 } 00345 EntityHandle newElement; 00346 rval = mb->create_element(typeElem, &new_conn[0], nnodes, newElement); 00347 CHECK_ERR(rval); 00348 //set the corresponding tag; not sure we need this one, from old to new 00349 /*rval = mb->tag_set_data(corrTag, &q, 1, &newElement); 00350 CHECK_ERR(rval);*/ 00351 rval = mb->tag_set_data(corrTag, &newElement, 1, &q); 00352 CHECK_ERR(rval); 00353 // set the global id 00354 rval = mb->tag_set_data(gid, &newElement, 1, &global_id); 00355 CHECK_ERR(rval); 00356 00357 rval = mb->add_entities(lagr_set, &newElement, 1); 00358 CHECK_ERR(rval); 00359 } 00360 00361 return MB_SUCCESS; 00362 } 00363 ErrorCode compute_tracer_case1(Interface * mb, Intx2MeshOnSphere & worker, EntityHandle euler_set, 00364 EntityHandle lagr_set, EntityHandle out_set, Tag & tagElem, Tag & tagArea, int rank, 00365 int tStep, Range & connecVerts) 00366 { 00367 ErrorCode rval = MB_SUCCESS; 00368 00369 if (!corrTag) 00370 return MB_FAILURE; 00371 double t = tStep * T / numSteps; // numSteps is global; so is T 00372 double delta_t = T / numSteps; // this is global too, actually 00373 Range polys; 00374 rval = mb->get_entities_by_dimension(euler_set, 2, polys); 00375 CHECK_ERR(rval); 00376 00377 // change coordinates of lagr mesh vertices 00378 for (Range::iterator vit = connecVerts.begin(); vit != connecVerts.end(); 00379 vit++) 00380 { 00381 EntityHandle oldV = *vit; 00382 CartVect posi; 00383 rval = mb->get_coords(&oldV, 1, &(posi[0])); 00384 CHECK_ERR(rval); 00385 // cslam utils, case 1 00386 CartVect newPos; 00387 departure_point_case1(posi, t, delta_t, newPos); 00388 newPos = radius * newPos; // do we need this? the radius should be 1 00389 EntityHandle new_vert; 00390 rval = mb->tag_get_data(corrTag, &oldV, 1, &new_vert); 00391 CHECK_ERR(rval); 00392 // set the new position for the new vertex 00393 rval = mb->set_coords(&new_vert, 1, &(newPos[0])); 00394 CHECK_ERR(rval); 00395 } 00396 00397 // if in parallel, we have to move some elements to another proc, and receive other cells 00398 // from other procs 00399 // lagr and euler are preserved 00400 if (writeFiles) // so if need to write lagr files too 00401 { 00402 00403 } 00404 EntityHandle covering_set; 00405 rval = worker.create_departure_mesh_3rd_alg(lagr_set, covering_set); 00406 if (writeFiles) // so if write 00407 { 00408 std::stringstream departureMesh; 00409 departureMesh << "Departure" << rank << "_" << tStep << ".vtk"; 00410 rval = mb->write_file(departureMesh.str().c_str(), 0, 0, &lagr_set, 1); 00411 CHECK_ERR(rval); 00412 00413 std::stringstream newTracer; 00414 newTracer << "Tracer" << rank << "_" << tStep << ".vtk"; 00415 rval = mb->write_file(newTracer.str().c_str(), 0, 0, &euler_set, 1); 00416 CHECK_ERR(rval); 00417 00418 std::stringstream lagr_cover; 00419 lagr_cover << "Cover" << rank << "_" << tStep << ".vtk"; 00420 rval = mb->write_file(lagr_cover.str().c_str(), 0, 0, &covering_set, 1); 00421 CHECK_ERR(rval); 00422 00423 } 00424 // so we have now the departure at the previous time 00425 // intersect the 2 meshes (what about some checking of convexity?) for sufficient 00426 // small dt, it is not an issue; 00427 00428 // std::cout << "error tolerance epsilon_1=" << gtol << "\n"; 00429 00430 rval = worker.intersect_meshes(covering_set, euler_set, out_set); 00431 CHECK_ERR(rval); 00432 if (writeFiles) // so if write 00433 { 00434 std::stringstream intx_mesh; 00435 intx_mesh << "Intx" << rank << "_" << tStep << ".vtk"; 00436 rval = mb->write_file(intx_mesh.str().c_str(), 0, 0, &out_set, 1); 00437 } 00438 00439 // serially: lagr is the same order as euler; 00440 // we need to update now the tracer information on each element, based on 00441 // initial value and areas of each resulting polygons 00442 if (parallelWrite && tStep==1) 00443 { 00444 std::stringstream resTrace; 00445 resTrace << "Tracer" << "_" << tStep-1 << ".h5m"; 00446 rval = mb->write_file(resTrace.str().c_str(), 0, "PARALLEL=WRITE_PART", &euler_set, 1, &tagElem, 1); 00447 } 00448 rval = worker.update_tracer_data(out_set, tagElem, tagArea); 00449 CHECK_ERR(rval); 00450 00451 if (parallelWrite) 00452 { 00453 std::stringstream resTrace; 00454 resTrace << "Tracer" << "_" << tStep << ".h5m"; 00455 rval = mb->write_file(resTrace.str().c_str(), 0, "PARALLEL=WRITE_PART", &euler_set, 1, &tagElem, 1); 00456 } 00457 00458 if (writeFiles) // so if write 00459 { 00460 std::stringstream newIntx; 00461 newIntx << "newIntx" << rank << "_" << tStep << ".vtk"; 00462 rval = mb->write_file(newIntx.str().c_str(), 0, 0, &out_set, 1); 00463 CHECK_ERR(rval); 00464 } 00465 // delete now the polygons and the elements of out_set 00466 // also, all verts that are not in euler set or lagr_set 00467 Range allVerts; 00468 rval = mb->get_entities_by_dimension(0, 0, allVerts); 00469 CHECK_ERR(rval); 00470 00471 Range allElems; 00472 rval = mb->get_entities_by_dimension(0, 2, allElems); 00473 CHECK_ERR(rval); 00474 // add to polys range the lagr polys 00475 rval = mb->get_entities_by_dimension(lagr_set, 2, polys); // do not delete lagr set either, with its vertices 00476 CHECK_ERR(rval); 00477 // add to the connecVerts range all verts, from all initial polys 00478 Range vertsToStay; 00479 rval = mb->get_connectivity(polys, vertsToStay); 00480 CHECK_ERR(rval); 00481 00482 Range todeleteVerts = subtract(allVerts, vertsToStay); 00483 00484 Range todeleteElem = subtract(allElems, polys); 00485 00486 // empty the out mesh set 00487 rval = mb->clear_meshset(&out_set, 1); 00488 CHECK_ERR(rval); 00489 00490 rval = mb->delete_entities(todeleteElem); 00491 CHECK_ERR(rval); 00492 rval = mb->delete_entities(todeleteVerts); 00493 CHECK_ERR(rval); 00494 if (rank==0) 00495 std::cout << " step: " << tStep << "\n"; 00496 return rval; 00497 } 00498 int main(int argc, char **argv) 00499 { 00500 00501 MPI_Init(&argc, &argv); 00502 LONG_DESC << "This program simulates a transport problem on a sphere" 00503 " according to a benchmark from a Nair & Lauritzen paper.\n" 00504 << "It starts with a partitioned mesh on a sphere, add a tracer, and steps through.\n" << 00505 "The flow reverses after half time, and it should return to original configuration, if the integration was exact. "; 00506 ProgOptions opts(LONG_DESC.str(), BRIEF_DESC); 00507 00508 // read a homme file, partitioned in 16 so far 00509 std::string fileN= TestDir + "/HN16.h5m"; 00510 const char *filename_mesh1 = fileN.c_str(); 00511 00512 opts.addOpt<double>("gtolerance,g", 00513 "geometric absolute tolerance (used for point concidence on the sphere)", >ol); 00514 00515 std::string input_file; 00516 opts.addOpt<std::string>("input_file,i", "input mesh file, partitioned", 00517 &input_file); 00518 std::string extra_read_opts; 00519 opts.addOpt<std::string>("extra_read_options,O", "extra read options ", 00520 &extra_read_opts); 00521 //int field_type; 00522 opts.addOpt<int>("field_type,f", 00523 "field type-- 1: quasi-smooth; 2: smooth; 3: slotted cylinders (non-smooth)", &field_type); 00524 00525 opts.addOpt<int>("num_steps,n", 00526 "number of steps ", &numSteps); 00527 00528 //bool reorder = false; 00529 opts.addOpt<void>("write_debug_files,w", "write debugging files during simulation ", 00530 &writeFiles); 00531 00532 opts.addOpt<void>("write_velocity_files,v", "Reorder mesh to group entities by partition", 00533 &velocity); 00534 00535 opts.addOpt<void>("write_result_in_parallel,p", "write tracer result files", 00536 ¶llelWrite); 00537 00538 opts.parseCommandLine(argc, argv); 00539 00540 if (!input_file.empty()) 00541 filename_mesh1=input_file.c_str(); 00542 00543 // read in parallel, in the "euler_set", the initial mesh 00544 std::string optsRead = std::string("PARALLEL=READ_PART;PARTITION=PARALLEL_PARTITION")+ 00545 std::string(";PARALLEL_RESOLVE_SHARED_ENTS")+extra_read_opts; 00546 Core moab; 00547 Interface & mb = moab; 00548 EntityHandle euler_set; 00549 ErrorCode rval; 00550 rval = mb.create_meshset(MESHSET_SET, euler_set); 00551 CHECK_ERR(rval); 00552 00553 rval = mb.load_file(filename_mesh1, &euler_set, optsRead.c_str()); 00554 00555 ParallelComm* pcomm = ParallelComm::get_pcomm(&mb, 0); 00556 CHECK_ERR(rval); 00557 00558 rval = pcomm->check_all_shared_handles(); 00559 CHECK_ERR(rval); 00560 00561 int rank = pcomm->proc_config().proc_rank(); 00562 00563 if (0==rank) 00564 { 00565 std::cout << " case 1: use -gtol " << gtol << 00566 " -R " << radius << " -input " << filename_mesh1 << " -f " << field_type << 00567 " numSteps: " << numSteps << "\n"; 00568 std::cout<<" write debug results: " << (writeFiles ? "yes" : "no") << "\n"; 00569 std::cout<< " write tracer in parallel: " << ( parallelWrite ? "yes" : "no") << "\n"; 00570 std::cout <<" output velocity: " << (velocity? "yes" : "no") << "\n"; 00571 } 00572 00573 // tagTracer is the value at nodes 00574 Tag tagTracer = 0; 00575 std::string tag_name("Tracer"); 00576 rval = mb.tag_get_handle(tag_name.c_str(), 1, MB_TYPE_DOUBLE, tagTracer, MB_TAG_DENSE | MB_TAG_CREAT); 00577 CHECK_ERR(rval); 00578 00579 // tagElem is the average computed at each element, from nodal values 00580 Tag tagElem = 0; 00581 std::string tag_name2("TracerAverage"); 00582 rval = mb.tag_get_handle(tag_name2.c_str(), 1, MB_TYPE_DOUBLE, tagElem, MB_TAG_DENSE | MB_TAG_CREAT); 00583 CHECK_ERR(rval); 00584 00585 // area of the euler element is fixed, store it; it is used to recompute the averages at each 00586 // time step 00587 Tag tagArea = 0; 00588 std::string tag_name4("Area"); 00589 rval = mb.tag_get_handle(tag_name4.c_str(), 1, MB_TYPE_DOUBLE, tagArea, MB_TAG_DENSE | MB_TAG_CREAT); 00590 CHECK_ERR(rval); 00591 00592 // add a field value, quasi smooth first 00593 rval = add_field_value(&mb, euler_set, rank, tagTracer, tagElem, tagArea); 00594 CHECK_ERR(rval); 00595 00596 // iniVals are used for 1-norm error computation 00597 Range redEls; 00598 rval = mb.get_entities_by_dimension(euler_set, 2, redEls); 00599 CHECK_ERR(rval); 00600 std::vector<double> iniVals(redEls.size()); 00601 rval = mb.tag_get_data(tagElem, redEls, &iniVals[0]); 00602 CHECK_ERR(rval); 00603 00604 Tag tagh = 0; 00605 std::string tag_name3("Case1"); 00606 rval = mb.tag_get_handle(tag_name3.c_str(), 3, MB_TYPE_DOUBLE, tagh, MB_TAG_DENSE | MB_TAG_CREAT); 00607 CHECK_ERR(rval); 00608 EntityHandle out_set, lagr_set; 00609 rval = mb.create_meshset(MESHSET_SET, out_set); 00610 CHECK_ERR(rval); 00611 rval = mb.create_meshset(MESHSET_SET, lagr_set); 00612 CHECK_ERR(rval); 00613 // copy the initial mesh in the lagrangian set 00614 // initial vertices will be at the same position as euler; 00615 00616 rval = create_lagr_mesh(&mb, euler_set, lagr_set); 00617 CHECK_ERR(rval); 00618 00619 Intx2MeshOnSphere worker(&mb); 00620 worker.SetRadius(radius); 00621 00622 worker.SetErrorTolerance(gtol); 00623 00624 Range local_verts; 00625 rval = worker.build_processor_euler_boxes(euler_set, local_verts);// output also the local_verts 00626 // these stay fixed for one run 00627 // other things from intersection might need to change, like input blue set (departure set) 00628 // so we need also a method to clean memory 00629 CHECK_ERR(rval); 00630 00631 for (int i=1; i<numSteps+1; i++) 00632 { 00633 // time depends on i; t = i*T/numSteps: ( 0, T/numSteps, 2*T/numSteps, ..., T ) 00634 // this is really just to create some plots; it is not really needed to proceed 00635 // the compute_tracer_case1 method actually computes the departure point position 00636 if (velocity) 00637 { 00638 rval = compute_velocity_case1(&mb, euler_set, tagh, rank, i); 00639 CHECK_ERR(rval); 00640 } 00641 00642 // this is to actually compute concentrations at time step i, using the 00643 // current concentrations 00644 // 00645 rval = compute_tracer_case1(&mb, worker, euler_set, lagr_set, out_set, 00646 tagElem, tagArea, rank, i, local_verts); 00647 CHECK_ERR(rval); 00648 00649 } 00650 00651 //final vals and 1-norm 00652 Range::iterator iter = redEls.begin(); 00653 double norm1 = 0.; 00654 int count =0; 00655 void * data; 00656 int j=0;// index in iniVals 00657 while (iter != redEls.end()) 00658 { 00659 rval = mb.tag_iterate(tagElem, iter, redEls.end(), count, data); 00660 CHECK_ERR(rval); 00661 double * ptrTracer=(double*)data; 00662 00663 rval = mb.tag_iterate(tagArea, iter, redEls.end(), count, data); 00664 CHECK_ERR(rval); 00665 double * ptrArea=(double*)data; 00666 for (int i=0; i<count; i++, iter++, ptrTracer++, ptrArea++, j++) 00667 { 00668 //double area = *ptrArea; 00669 norm1+=fabs(*ptrTracer - iniVals[j])* (*ptrArea); 00670 } 00671 } 00672 00673 double total_norm1=0; 00674 int mpi_err = MPI_Reduce(&norm1, &total_norm1, 1, MPI_DOUBLE, MPI_SUM, 0, MPI_COMM_WORLD); 00675 if (MPI_SUCCESS != mpi_err) return 1; 00676 if (0==rank) 00677 std::cout << " numSteps:" << numSteps << " 1-norm:" << total_norm1 << "\n"; 00678 MPI_Finalize(); 00679 return 0; 00680 }