MOAB: Mesh Oriented datABase  (version 5.4.1)
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 <ctime>
00028 #include <cstdlib>
00029 #include <cstdio>
00030 #include <cstring>
00031 #include <cmath>  // for M_PI
00032 
00033 #include "moab/Core.hpp"
00034 #include "moab/Interface.hpp"
00035 #include "moab/IntxMesh/Intx2MeshOnSphere.hpp"
00036 #include "moab/ProgOptions.hpp"
00037 #include "MBTagConventions.hpp"
00038 #include "moab/ParallelComm.hpp"
00039 #include "moab/IntxMesh/IntxUtils.hpp"
00040 #include "IntxUtilsCSLAM.hpp"
00041 
00042 #include "TestUtil.hpp"
00043 
00044 const char BRIEF_DESC[] = "Simulate a transport problem in a semi-Lagrangian formulation\n";
00045 std::ostringstream LONG_DESC;
00046 
00047 // non smooth scalar field
00048 // some input data
00049 double gtol   = 1.e-9;  // this is for geometry tolerance
00050 double radius = 1.;     // in m:  6371220.
00051 
00052 int numSteps       = 3;  // number of times with velocity displayed at points
00053 double T           = 1;
00054 int case_number    = 1;  // 1, 2 (non-divergent) 3 divergent
00055 bool writeFiles    = false;
00056 bool parallelWrite = false;
00057 bool velocity      = false;
00058 int field_type     = 1;  // 1 quasi smooth, 2 - smooth, 3 non-smooth,
00059 
00060 using namespace moab;
00061 ErrorCode add_field_value( Interface* mb, EntityHandle euler_set, int rank, Tag& tagTracer, Tag& tagElem, Tag& tagArea )
00062 {
00063     /*
00064      * get all plys first, then vertices, then move them on the surface of the sphere
00065      *  radius is 1., most of the time
00066      *
00067      */
00068     Range polygons;
00069     ErrorCode rval = mb->get_entities_by_dimension( euler_set, 2, polygons );
00070     if( MB_SUCCESS != rval ) return rval;
00071 
00072     Range connecVerts;
00073     rval = mb->get_connectivity( polygons, connecVerts );
00074     if( MB_SUCCESS != rval ) return rval;
00075 
00076     void* data;  // pointer to the LOC in memory, for each vertex
00077     int count;
00078 
00079     rval = mb->tag_iterate( tagTracer, connecVerts.begin(), connecVerts.end(), count, data );MB_CHK_ERR( rval );
00080     // here we are checking contiguity
00081     assert( count == (int)connecVerts.size() );
00082     double* ptr_DP = (double*)data;
00083     // lambda is for longitude, theta for latitude
00084     // param will be: (la1, te1), (la2, te2), b, c; hmax=1, r=1/2
00085     // nondivergent flow, page 5, case 1, (la1, te1) = (M_PI, M_PI/3)
00086     //                                    (la2, te2) = (M_PI, -M_PI/3)
00087     //                 la1,    te1    la2    te2     b     c  hmax  r
00088     if( field_type == 1 )  // quasi smooth
00089     {
00090         double params[] = { M_PI, M_PI / 3, M_PI, -M_PI / 3, 0.1, 0.9, 1., 0.5 };
00091         for( Range::iterator vit = connecVerts.begin(); vit != connecVerts.end(); ++vit )
00092         {
00093             EntityHandle oldV = *vit;
00094             CartVect posi;
00095             rval = mb->get_coords( &oldV, 1, &( posi[0] ) );MB_CHK_ERR( rval );
00096 
00097             moab::IntxUtils::SphereCoords sphCoord = moab::IntxUtils::cart_to_spherical( posi );
00098 
00099             ptr_DP[0] = IntxUtilsCSLAM::quasi_smooth_field( sphCoord.lon, sphCoord.lat, params );
00100             ;
00101 
00102             ptr_DP++;  // increment to the next node
00103         }
00104     }
00105     else if( 2 == field_type )  // smooth
00106     {
00107         CartVect p1, p2;
00108         moab::IntxUtils::SphereCoords spr;
00109         spr.R   = 1;
00110         spr.lat = M_PI / 3;
00111         spr.lon = M_PI;
00112         p1      = moab::IntxUtils::spherical_to_cart( spr );
00113         spr.lat = -M_PI / 3;
00114         p2      = moab::IntxUtils::spherical_to_cart( spr );
00115         //                  x1,    y1,     z1,    x2,   y2,    z2,   h_max, b0
00116         double params[] = { p1[0], p1[1], p1[2], p2[0], p2[1], p2[2], 1, 5. };
00117         for( Range::iterator vit = connecVerts.begin(); vit != connecVerts.end(); ++vit )
00118         {
00119             EntityHandle oldV = *vit;
00120             CartVect posi;
00121             rval = mb->get_coords( &oldV, 1, &( posi[0] ) );MB_CHK_ERR( rval );
00122 
00123             moab::IntxUtils::SphereCoords sphCoord = moab::IntxUtils::cart_to_spherical( posi );
00124 
00125             ptr_DP[0] = IntxUtilsCSLAM::smooth_field( sphCoord.lon, sphCoord.lat, params );
00126             ;
00127 
00128             ptr_DP++;  // increment to the next node
00129         }
00130     }
00131     else if( 3 == field_type )  // slotted
00132     {
00133         //                   la1, te1,   la2, te2,       b,   c,   r
00134         double params[] = { M_PI, M_PI / 3, M_PI, -M_PI / 3, 0.1, 0.9, 0.5 };  // no h_max
00135         for( Range::iterator vit = connecVerts.begin(); vit != connecVerts.end(); ++vit )
00136         {
00137             EntityHandle oldV = *vit;
00138             CartVect posi;
00139             rval = mb->get_coords( &oldV, 1, &( posi[0] ) );MB_CHK_ERR( rval );
00140 
00141             moab::IntxUtils::SphereCoords sphCoord = moab::IntxUtils::cart_to_spherical( posi );
00142 
00143             ptr_DP[0] = IntxUtilsCSLAM::slotted_cylinder_field( sphCoord.lon, sphCoord.lat, params );
00144             ;
00145 
00146             ptr_DP++;  // increment to the next node
00147         }
00148     }
00149 
00150     // add average value for quad/polygon (average corners)
00151     // do some averages
00152 
00153     Range::iterator iter = polygons.begin();
00154     double local_mass    = 0.;  // this is total mass on one proc
00155     while( iter != polygons.end() )
00156     {
00157         rval = mb->tag_iterate( tagElem, iter, polygons.end(), count, data );MB_CHK_ERR( rval );
00158         double* ptr = (double*)data;
00159 
00160         rval = mb->tag_iterate( tagArea, iter, polygons.end(), count, data );MB_CHK_ERR( rval );
00161         double* ptrArea = (double*)data;
00162         for( int i = 0; i < count; i++, ++iter, ptr++, ptrArea++ )
00163         {
00164             const moab::EntityHandle* conn = NULL;
00165             int num_nodes                  = 0;
00166             rval                           = mb->get_connectivity( *iter, conn, num_nodes );MB_CHK_ERR( rval );
00167             if( num_nodes == 0 ) return MB_FAILURE;
00168             std::vector< double > nodeVals( num_nodes );
00169             double average = 0.;
00170             rval           = mb->tag_get_data( tagTracer, conn, num_nodes, &nodeVals[0] );MB_CHK_ERR( rval );
00171             for( int j = 0; j < num_nodes; j++ )
00172                 average += nodeVals[j];
00173             average /= num_nodes;
00174             *ptr = average;
00175 
00176             // now get area
00177             std::vector< double > coords;
00178             coords.resize( 3 * num_nodes );
00179             rval = mb->get_coords( conn, num_nodes, &coords[0] );MB_CHK_ERR( rval );
00180 
00181             moab::IntxAreaUtils sphAreaUtils;
00182             *ptrArea = sphAreaUtils.area_spherical_polygon( &coords[0], num_nodes, radius );
00183 
00184             // we should have used some
00185             // total mass:
00186             local_mass += *ptrArea * average;
00187         }
00188     }
00189     double total_mass = 0.;
00190     int mpi_err       = MPI_Reduce( &local_mass, &total_mass, 1, MPI_DOUBLE, MPI_SUM, 0, MPI_COMM_WORLD );
00191     if( MPI_SUCCESS != mpi_err ) return MB_FAILURE;
00192 
00193     if( rank == 0 ) std::cout << "initial total mass:" << total_mass << "\n";
00194 
00195     // now we can delete the tags? not yet
00196     return MB_SUCCESS;
00197 }
00198 
00199 ErrorCode compute_velocity_case1( Interface* mb, EntityHandle euler_set, Tag& tagh, int rank, int tStep )
00200 {
00201     Range polygons;
00202     ErrorCode rval = mb->get_entities_by_dimension( euler_set, 2, polygons );
00203     if( MB_SUCCESS != rval ) return rval;
00204 
00205     Range connecVerts;
00206     rval = mb->get_connectivity( polygons, connecVerts );
00207     if( MB_SUCCESS != rval ) return rval;
00208 
00209     void* data;  // pointer to the velo in memory, for each vertex
00210     int count;
00211 
00212     rval = mb->tag_iterate( tagh, connecVerts.begin(), connecVerts.end(), count, data );MB_CHK_ERR( rval );
00213     // here we are checking contiguity
00214     assert( count == (int)connecVerts.size() );
00215     double* ptr_velo = (double*)data;
00216     // lambda is for longitude, theta for latitude
00217 
00218     for( Range::iterator vit = connecVerts.begin(); vit != connecVerts.end(); ++vit )
00219     {
00220         EntityHandle oldV = *vit;
00221         CartVect posi;
00222         rval = mb->get_coords( &oldV, 1, &( posi[0] ) );MB_CHK_ERR( rval );
00223         CartVect velo;
00224         double t = T * tStep / numSteps;  //
00225         IntxUtilsCSLAM::velocity_case1( posi, t, velo );
00226 
00227         ptr_velo[0] = velo[0];
00228         ptr_velo[1] = velo[1];
00229         ptr_velo[2] = velo[2];
00230 
00231         // increment to the next node
00232         ptr_velo += 3;  // to next velocity
00233     }
00234 
00235     std::stringstream velos;
00236     velos << "Tracer" << rank << "_" << tStep << ".vtk";
00237     rval = mb->write_file( velos.str().c_str(), 0, 0, &euler_set, 1, &tagh, 1 );MB_CHK_ERR( rval );
00238 
00239     return MB_SUCCESS;
00240 }
00241 
00242 ErrorCode compute_tracer_case1( Interface* mb,
00243                                 Intx2MeshOnSphere& worker,
00244                                 EntityHandle euler_set,
00245                                 EntityHandle lagr_set,
00246                                 EntityHandle out_set,
00247                                 Tag& tagElem,
00248                                 Tag& tagArea,
00249                                 int rank,
00250                                 int tStep,
00251                                 Range& connecVerts )
00252 {
00253     ErrorCode rval;
00254     EntityHandle dum = 0;
00255     Tag corrTag;
00256     rval = mb->tag_get_handle( CORRTAGNAME, 1, MB_TYPE_HANDLE, corrTag, MB_TAG_DENSE, &dum );MB_CHK_ERR( rval );
00257 
00258     double t       = tStep * T / numSteps;  // numSteps is global; so is T
00259     double delta_t = T / numSteps;          // this is global too, actually
00260     Range polys;
00261     rval = mb->get_entities_by_dimension( euler_set, 2, polys );MB_CHK_ERR( rval );
00262 
00263     // change coordinates of lagr mesh vertices
00264     for( Range::iterator vit = connecVerts.begin(); vit != connecVerts.end(); ++vit )
00265     {
00266         EntityHandle oldV = *vit;
00267         CartVect posi;
00268         rval = mb->get_coords( &oldV, 1, &( posi[0] ) );MB_CHK_ERR( rval );
00269         // Intx utils, case 1
00270         CartVect newPos;
00271         IntxUtilsCSLAM::departure_point_case1( posi, t, delta_t, newPos );
00272         newPos = radius * newPos;  // do we need this? the radius should be 1
00273         EntityHandle new_vert;
00274         rval = mb->tag_get_data( corrTag, &oldV, 1, &new_vert );MB_CHK_ERR( rval );
00275         // set the new position for the new vertex
00276         rval = mb->set_coords( &new_vert, 1, &( newPos[0] ) );MB_CHK_ERR( rval );
00277     }
00278 
00279     // if in parallel, we have to move some elements to another proc, and receive other cells
00280     // from other procs
00281     // lagr and euler are preserved
00282     EntityHandle covering_set;
00283     rval = worker.create_departure_mesh_3rd_alg( lagr_set, covering_set );MB_CHK_ERR( rval );
00284     if( writeFiles )  // so if write
00285     {
00286         std::stringstream departureMesh;
00287         departureMesh << "Departure" << rank << "_" << tStep << ".vtk";
00288         rval = mb->write_file( departureMesh.str().c_str(), 0, 0, &lagr_set, 1 );MB_CHK_ERR( rval );
00289 
00290         std::stringstream newTracer;
00291         newTracer << "Tracer" << rank << "_" << tStep << ".vtk";
00292         rval = mb->write_file( newTracer.str().c_str(), 0, 0, &euler_set, 1 );MB_CHK_ERR( rval );
00293 
00294         std::stringstream lagr_cover;
00295         lagr_cover << "Cover" << rank << "_" << tStep << ".vtk";
00296         rval = mb->write_file( lagr_cover.str().c_str(), 0, 0, &covering_set, 1 );MB_CHK_ERR( rval );
00297     }
00298     // so we have now the departure at the previous time
00299     // intersect the 2 meshes (what about some checking of convexity?) for sufficient
00300     // small dt, it is not an issue;
00301 
00302     rval = worker.intersect_meshes( covering_set, euler_set, out_set );MB_CHK_ERR( rval );
00303     if( writeFiles )  // so if write
00304     {
00305         std::stringstream intx_mesh;
00306         intx_mesh << "Intx" << rank << "_" << tStep << ".vtk";
00307         rval = mb->write_file( intx_mesh.str().c_str(), 0, 0, &out_set, 1 );MB_CHK_ERR( rval );
00308     }
00309 
00310     // serially: lagr is the same order as euler;
00311     // we need to update now the tracer information on each element, based on
00312     // initial value and areas of each resulting polygons
00313     if( parallelWrite && tStep == 1 )
00314     {
00315         std::stringstream resTrace;
00316         resTrace << "Tracer"
00317                  << "_" << tStep - 1 << ".h5m";
00318         rval = mb->write_file( resTrace.str().c_str(), 0, "PARALLEL=WRITE_PART", &euler_set, 1, &tagElem, 1 );MB_CHK_ERR( rval );
00319     }
00320     rval = worker.update_tracer_data( out_set, tagElem, tagArea );MB_CHK_ERR( rval );
00321 
00322     if( parallelWrite )
00323     {
00324         std::stringstream resTrace;
00325         resTrace << "Tracer"
00326                  << "_" << tStep << ".h5m";
00327         rval = mb->write_file( resTrace.str().c_str(), 0, "PARALLEL=WRITE_PART", &euler_set, 1, &tagElem, 1 );
00328     }
00329 
00330     if( writeFiles )  // so if write
00331     {
00332         std::stringstream newIntx;
00333         newIntx << "newIntx" << rank << "_" << tStep << ".vtk";
00334         rval = mb->write_file( newIntx.str().c_str(), 0, 0, &out_set, 1 );MB_CHK_ERR( rval );
00335     }
00336     // delete now the polygons and the elements of out_set
00337     // also, all verts that are not in euler set or lagr_set
00338     Range allVerts;
00339     rval = mb->get_entities_by_dimension( 0, 0, allVerts );MB_CHK_ERR( rval );
00340 
00341     Range allElems;
00342     rval = mb->get_entities_by_dimension( 0, 2, allElems );MB_CHK_ERR( rval );
00343 
00344     // add to polys range the lagr polys
00345     // do not delete lagr set either, with its vertices
00346     rval = mb->get_entities_by_dimension( lagr_set, 2, polys );MB_CHK_ERR( rval );
00347     // add to the connecVerts range all verts, from all initial polys
00348     Range vertsToStay;
00349     rval = mb->get_connectivity( polys, vertsToStay );MB_CHK_ERR( rval );
00350 
00351     Range todeleteVerts = subtract( allVerts, vertsToStay );
00352 
00353     Range todeleteElem = subtract( allElems, polys );
00354 
00355     // empty the out mesh set
00356     rval = mb->clear_meshset( &out_set, 1 );MB_CHK_ERR( rval );
00357 
00358     rval = mb->delete_entities( todeleteElem );MB_CHK_ERR( rval );
00359     rval = mb->delete_entities( todeleteVerts );MB_CHK_ERR( rval );
00360     if( rank == 0 ) std::cout << " step: " << tStep << "\n";
00361     return rval;
00362 }
00363 
00364 int main( int argc, char** argv )
00365 {
00366     MPI_Init( &argc, &argv );
00367     LONG_DESC << "This program simulates a transport problem on a sphere"
00368                  " according to a benchmark from a Nair & Lauritzen paper.\n"
00369               << "It starts with a partitioned mesh on a sphere, add a tracer, and steps through.\n"
00370               << "The flow reverses after half time, and it should return to original "
00371                  "configuration, if the integration was exact. ";
00372     ProgOptions opts( LONG_DESC.str(), BRIEF_DESC );
00373 
00374     // read a homme file, partitioned in 16 so far
00375     std::string fileN          = TestDir + "unittest/mbcslam/fine4.h5m";
00376     const char* filename_mesh1 = fileN.c_str();
00377 
00378     opts.addOpt< double >( "gtolerance,g", "geometric absolute tolerance (used for point concidence on the sphere)",
00379                            &gtol );
00380 
00381     std::string input_file;
00382     opts.addOpt< std::string >( "input_file,i", "input mesh file, partitioned", &input_file );
00383     std::string extra_read_opts;
00384     opts.addOpt< std::string >( "extra_read_options,O", "extra read options ", &extra_read_opts );
00385     // int field_type;
00386     opts.addOpt< int >( "field_type,f", "field type--  1: quasi-smooth; 2: smooth; 3: slotted cylinders (non-smooth)",
00387                         &field_type );
00388 
00389     opts.addOpt< int >( "num_steps,n", "number of  steps ", &numSteps );
00390 
00391     // bool reorder = false;
00392     opts.addOpt< void >( "write_debug_files,w", "write debugging files during simulation ", &writeFiles );
00393 
00394     opts.addOpt< void >( "write_velocity_files,v", "Reorder mesh to group entities by partition", &velocity );
00395 
00396     opts.addOpt< void >( "write_result_in_parallel,p", "write tracer result files", &parallelWrite );
00397 
00398     opts.parseCommandLine( argc, argv );
00399 
00400     if( !input_file.empty() ) filename_mesh1 = input_file.c_str();
00401 
00402     // read in parallel, in the "euler_set", the initial mesh
00403     std::string optsRead = std::string( "PARALLEL=READ_PART;PARTITION=PARALLEL_PARTITION" ) +
00404                            std::string( ";PARALLEL_RESOLVE_SHARED_ENTS" ) + extra_read_opts;
00405     Core moab;
00406     Interface& mb = moab;
00407     EntityHandle euler_set;
00408     ErrorCode rval;
00409     rval = mb.create_meshset( MESHSET_SET, euler_set );MB_CHK_ERR( rval );
00410 
00411     rval = mb.load_file( filename_mesh1, &euler_set, optsRead.c_str() );MB_CHK_ERR( rval );
00412 
00413     ParallelComm* pcomm = ParallelComm::get_pcomm( &mb, 0 );
00414 
00415     rval = pcomm->check_all_shared_handles();MB_CHK_ERR( rval );
00416 
00417     int rank = pcomm->proc_config().proc_rank();
00418 
00419     if( 0 == rank )
00420     {
00421         std::cout << " case 1: use -gtol " << gtol << " -R " << radius << " -input " << filename_mesh1 << " -f "
00422                   << field_type << " numSteps: " << numSteps << "\n";
00423         std::cout << " write debug results: " << ( writeFiles ? "yes" : "no" ) << "\n";
00424         std::cout << " write tracer in parallel: " << ( parallelWrite ? "yes" : "no" ) << "\n";
00425         std::cout << " output velocity: " << ( velocity ? "yes" : "no" ) << "\n";
00426     }
00427 
00428     // tagTracer is the value at nodes
00429     Tag tagTracer = 0;
00430     std::string tag_name( "Tracer" );
00431     rval = mb.tag_get_handle( tag_name.c_str(), 1, MB_TYPE_DOUBLE, tagTracer, MB_TAG_DENSE | MB_TAG_CREAT );MB_CHK_ERR( rval );
00432 
00433     // tagElem is the average computed at each element, from nodal values
00434     Tag tagElem = 0;
00435     std::string tag_name2( "TracerAverage" );
00436     rval = mb.tag_get_handle( tag_name2.c_str(), 1, MB_TYPE_DOUBLE, tagElem, MB_TAG_DENSE | MB_TAG_CREAT );MB_CHK_ERR( rval );
00437 
00438     // area of the euler element is fixed, store it; it is used to recompute the averages at each
00439     // time step
00440     Tag tagArea = 0;
00441     std::string tag_name4( "Area" );
00442     rval = mb.tag_get_handle( tag_name4.c_str(), 1, MB_TYPE_DOUBLE, tagArea, MB_TAG_DENSE | MB_TAG_CREAT );MB_CHK_ERR( rval );
00443 
00444     // add a field value, quasi smooth first
00445     rval = add_field_value( &mb, euler_set, rank, tagTracer, tagElem, tagArea );MB_CHK_ERR( rval );
00446 
00447     // iniVals are used for 1-norm error computation
00448     Range redEls;
00449     rval = mb.get_entities_by_dimension( euler_set, 2, redEls );MB_CHK_ERR( rval );
00450     std::vector< double > iniVals( redEls.size() );
00451     rval = mb.tag_get_data( tagElem, redEls, &iniVals[0] );MB_CHK_ERR( rval );
00452 
00453     Tag tagh = 0;
00454     std::string tag_name3( "Case1" );
00455     rval = mb.tag_get_handle( tag_name3.c_str(), 3, MB_TYPE_DOUBLE, tagh, MB_TAG_DENSE | MB_TAG_CREAT );MB_CHK_ERR( rval );
00456 
00457     EntityHandle out_set, lagr_set;
00458     rval = mb.create_meshset( MESHSET_SET, out_set );MB_CHK_ERR( rval );
00459     rval = mb.create_meshset( MESHSET_SET, lagr_set );MB_CHK_ERR( rval );
00460     // copy the initial mesh in the lagrangian set
00461     // initial vertices will be at the same position as euler;
00462 
00463     rval = IntxUtilsCSLAM::deep_copy_set( &mb, euler_set, lagr_set );MB_CHK_ERR( rval );
00464 
00465     Intx2MeshOnSphere worker( &mb );
00466     worker.set_radius_source_mesh( radius );
00467     worker.set_radius_destination_mesh( radius );
00468     worker.set_error_tolerance( gtol );
00469     worker.set_parallel_comm( pcomm );
00470 
00471     rval = worker.FindMaxEdges( lagr_set, euler_set );MB_CHK_ERR( rval );
00472     Range local_verts;
00473     // output also the local_verts
00474     rval = worker.build_processor_euler_boxes( euler_set, local_verts );MB_CHK_ERR( rval );
00475     // these stay fixed for one run
00476     // other things from intersection might need to change, like input blue set (departure set)
00477     // so we need also a method to clean memory
00478 
00479     for( int i = 1; i < numSteps + 1; i++ )
00480     {
00481         // time depends on i; t = i*T/numSteps: ( 0, T/numSteps, 2*T/numSteps, ..., T )
00482         // this is really just to create some plots; it is not really needed to proceed
00483         // the compute_tracer_case1 method actually computes the departure point position
00484         if( velocity )
00485         {
00486             rval = compute_velocity_case1( &mb, euler_set, tagh, rank, i );MB_CHK_ERR( rval );
00487         }
00488 
00489         // this is to actually compute concentrations at time step i, using the
00490         //  current concentrations
00491         rval =
00492             compute_tracer_case1( &mb, worker, euler_set, lagr_set, out_set, tagElem, tagArea, rank, i, local_verts );MB_CHK_ERR( rval );
00493     }
00494 
00495     // final vals and 1-norm
00496     Range::iterator iter = redEls.begin();
00497     double norm1         = 0.;
00498     int count            = 0;
00499     void* data;
00500     int j = 0;  // index in iniVals
00501     while( iter != redEls.end() )
00502     {
00503         rval = mb.tag_iterate( tagElem, iter, redEls.end(), count, data );MB_CHK_ERR( rval );
00504         double* ptrTracer = (double*)data;
00505 
00506         rval = mb.tag_iterate( tagArea, iter, redEls.end(), count, data );MB_CHK_ERR( rval );
00507         double* ptrArea = (double*)data;
00508         for( int i = 0; i < count; i++, ++iter, ptrTracer++, ptrArea++, j++ )
00509         {
00510             // double area = *ptrArea;
00511             norm1 += fabs( *ptrTracer - iniVals[j] ) * ( *ptrArea );
00512         }
00513     }
00514 
00515     double total_norm1 = 0;
00516     int mpi_err        = MPI_Reduce( &norm1, &total_norm1, 1, MPI_DOUBLE, MPI_SUM, 0, MPI_COMM_WORLD );
00517     if( MPI_SUCCESS != mpi_err ) return 1;
00518     if( 0 == rank ) std::cout << " numSteps:" << numSteps << " 1-norm:" << total_norm1 << "\n";
00519     MPI_Finalize();
00520     return 0;
00521 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines