moab
diffusion.cpp
Go to the documentation of this file.
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)", &gtol);
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      &parallelWrite);
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 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines