MOAB: Mesh Oriented datABase  (version 5.2.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, Intx2MeshOnSphere& worker, EntityHandle euler_set, EntityHandle lagr_set,
00243                                 EntityHandle out_set, Tag& tagElem, Tag& tagArea, int rank, int tStep,
00244                                 Range& connecVerts )
00245 {
00246     ErrorCode rval;
00247     EntityHandle dum = 0;
00248     Tag corrTag;
00249     rval = mb->tag_get_handle( CORRTAGNAME, 1, MB_TYPE_HANDLE, corrTag, MB_TAG_DENSE, &dum );MB_CHK_ERR( rval );
00250 
00251     double t       = tStep * T / numSteps;  // numSteps is global; so is T
00252     double delta_t = T / numSteps;          // this is global too, actually
00253     Range polys;
00254     rval = mb->get_entities_by_dimension( euler_set, 2, polys );MB_CHK_ERR( rval );
00255 
00256     // change coordinates of lagr mesh vertices
00257     for( Range::iterator vit = connecVerts.begin(); vit != connecVerts.end(); ++vit )
00258     {
00259         EntityHandle oldV = *vit;
00260         CartVect posi;
00261         rval = mb->get_coords( &oldV, 1, &( posi[0] ) );MB_CHK_ERR( rval );
00262         // Intx utils, case 1
00263         CartVect newPos;
00264         IntxUtilsCSLAM::departure_point_case1( posi, t, delta_t, newPos );
00265         newPos = radius * newPos;  // do we need this? the radius should be 1
00266         EntityHandle new_vert;
00267         rval = mb->tag_get_data( corrTag, &oldV, 1, &new_vert );MB_CHK_ERR( rval );
00268         // set the new position for the new vertex
00269         rval = mb->set_coords( &new_vert, 1, &( newPos[0] ) );MB_CHK_ERR( rval );
00270     }
00271 
00272     // if in parallel, we have to move some elements to another proc, and receive other cells
00273     // from other procs
00274     // lagr and euler are preserved
00275     EntityHandle covering_set;
00276     rval = worker.create_departure_mesh_3rd_alg( lagr_set, covering_set );MB_CHK_ERR( rval );
00277     if( writeFiles )  // so if write
00278     {
00279         std::stringstream departureMesh;
00280         departureMesh << "Departure" << rank << "_" << tStep << ".vtk";
00281         rval = mb->write_file( departureMesh.str().c_str(), 0, 0, &lagr_set, 1 );MB_CHK_ERR( rval );
00282 
00283         std::stringstream newTracer;
00284         newTracer << "Tracer" << rank << "_" << tStep << ".vtk";
00285         rval = mb->write_file( newTracer.str().c_str(), 0, 0, &euler_set, 1 );MB_CHK_ERR( rval );
00286 
00287         std::stringstream lagr_cover;
00288         lagr_cover << "Cover" << rank << "_" << tStep << ".vtk";
00289         rval = mb->write_file( lagr_cover.str().c_str(), 0, 0, &covering_set, 1 );MB_CHK_ERR( rval );
00290     }
00291     // so we have now the departure at the previous time
00292     // intersect the 2 meshes (what about some checking of convexity?) for sufficient
00293     // small dt, it is not an issue;
00294 
00295     rval = worker.intersect_meshes( covering_set, euler_set, out_set );MB_CHK_ERR( rval );
00296     if( writeFiles )  // so if write
00297     {
00298         std::stringstream intx_mesh;
00299         intx_mesh << "Intx" << rank << "_" << tStep << ".vtk";
00300         rval = mb->write_file( intx_mesh.str().c_str(), 0, 0, &out_set, 1 );MB_CHK_ERR( rval );
00301     }
00302 
00303     // serially: lagr is the same order as euler;
00304     // we need to update now the tracer information on each element, based on
00305     // initial value and areas of each resulting polygons
00306     if( parallelWrite && tStep == 1 )
00307     {
00308         std::stringstream resTrace;
00309         resTrace << "Tracer"
00310                  << "_" << tStep - 1 << ".h5m";
00311         rval = mb->write_file( resTrace.str().c_str(), 0, "PARALLEL=WRITE_PART", &euler_set, 1, &tagElem, 1 );MB_CHK_ERR( rval );
00312     }
00313     rval = worker.update_tracer_data( out_set, tagElem, tagArea );MB_CHK_ERR( rval );
00314 
00315     if( parallelWrite )
00316     {
00317         std::stringstream resTrace;
00318         resTrace << "Tracer"
00319                  << "_" << tStep << ".h5m";
00320         rval = mb->write_file( resTrace.str().c_str(), 0, "PARALLEL=WRITE_PART", &euler_set, 1, &tagElem, 1 );
00321     }
00322 
00323     if( writeFiles )  // so if write
00324     {
00325         std::stringstream newIntx;
00326         newIntx << "newIntx" << rank << "_" << tStep << ".vtk";
00327         rval = mb->write_file( newIntx.str().c_str(), 0, 0, &out_set, 1 );MB_CHK_ERR( rval );
00328     }
00329     // delete now the polygons and the elements of out_set
00330     // also, all verts that are not in euler set or lagr_set
00331     Range allVerts;
00332     rval = mb->get_entities_by_dimension( 0, 0, allVerts );MB_CHK_ERR( rval );
00333 
00334     Range allElems;
00335     rval = mb->get_entities_by_dimension( 0, 2, allElems );MB_CHK_ERR( rval );
00336 
00337     // add to polys range the lagr polys
00338     // do not delete lagr set either, with its vertices
00339     rval = mb->get_entities_by_dimension( lagr_set, 2, polys );MB_CHK_ERR( rval );
00340     // add to the connecVerts range all verts, from all initial polys
00341     Range vertsToStay;
00342     rval = mb->get_connectivity( polys, vertsToStay );MB_CHK_ERR( rval );
00343 
00344     Range todeleteVerts = subtract( allVerts, vertsToStay );
00345 
00346     Range todeleteElem = subtract( allElems, polys );
00347 
00348     // empty the out mesh set
00349     rval = mb->clear_meshset( &out_set, 1 );MB_CHK_ERR( rval );
00350 
00351     rval = mb->delete_entities( todeleteElem );MB_CHK_ERR( rval );
00352     rval = mb->delete_entities( todeleteVerts );MB_CHK_ERR( rval );
00353     if( rank == 0 ) std::cout << " step: " << tStep << "\n";
00354     return rval;
00355 }
00356 
00357 int main( int argc, char** argv )
00358 {
00359     MPI_Init( &argc, &argv );
00360     LONG_DESC << "This program simulates a transport problem on a sphere"
00361                  " according to a benchmark from a Nair & Lauritzen paper.\n"
00362               << "It starts with a partitioned mesh on a sphere, add a tracer, and steps through.\n"
00363               << "The flow reverses after half time, and it should return to original "
00364                  "configuration, if the integration was exact. ";
00365     ProgOptions opts( LONG_DESC.str(), BRIEF_DESC );
00366 
00367     // read a homme file, partitioned in 16 so far
00368     std::string fileN          = TestDir + "/mbcslam/fine4.h5m";
00369     const char* filename_mesh1 = fileN.c_str();
00370 
00371     opts.addOpt< double >( "gtolerance,g", "geometric absolute tolerance (used for point concidence on the sphere)",
00372                            &gtol );
00373 
00374     std::string input_file;
00375     opts.addOpt< std::string >( "input_file,i", "input mesh file, partitioned", &input_file );
00376     std::string extra_read_opts;
00377     opts.addOpt< std::string >( "extra_read_options,O", "extra read options ", &extra_read_opts );
00378     // int field_type;
00379     opts.addOpt< int >( "field_type,f", "field type--  1: quasi-smooth; 2: smooth; 3: slotted cylinders (non-smooth)",
00380                         &field_type );
00381 
00382     opts.addOpt< int >( "num_steps,n", "number of  steps ", &numSteps );
00383 
00384     // bool reorder = false;
00385     opts.addOpt< void >( "write_debug_files,w", "write debugging files during simulation ", &writeFiles );
00386 
00387     opts.addOpt< void >( "write_velocity_files,v", "Reorder mesh to group entities by partition", &velocity );
00388 
00389     opts.addOpt< void >( "write_result_in_parallel,p", "write tracer result files", &parallelWrite );
00390 
00391     opts.parseCommandLine( argc, argv );
00392 
00393     if( !input_file.empty() ) filename_mesh1 = input_file.c_str();
00394 
00395     // read in parallel, in the "euler_set", the initial mesh
00396     std::string optsRead = std::string( "PARALLEL=READ_PART;PARTITION=PARALLEL_PARTITION" ) +
00397                            std::string( ";PARALLEL_RESOLVE_SHARED_ENTS" ) + extra_read_opts;
00398     Core moab;
00399     Interface& mb = moab;
00400     EntityHandle euler_set;
00401     ErrorCode rval;
00402     rval = mb.create_meshset( MESHSET_SET, euler_set );MB_CHK_ERR( rval );
00403 
00404     rval = mb.load_file( filename_mesh1, &euler_set, optsRead.c_str() );MB_CHK_ERR( rval );
00405 
00406     ParallelComm* pcomm = ParallelComm::get_pcomm( &mb, 0 );
00407 
00408     rval = pcomm->check_all_shared_handles();MB_CHK_ERR( rval );
00409 
00410     int rank = pcomm->proc_config().proc_rank();
00411 
00412     if( 0 == rank )
00413     {
00414         std::cout << " case 1: use -gtol " << gtol << " -R " << radius << " -input " << filename_mesh1 << " -f "
00415                   << field_type << " numSteps: " << numSteps << "\n";
00416         std::cout << " write debug results: " << ( writeFiles ? "yes" : "no" ) << "\n";
00417         std::cout << " write tracer in parallel: " << ( parallelWrite ? "yes" : "no" ) << "\n";
00418         std::cout << " output velocity: " << ( velocity ? "yes" : "no" ) << "\n";
00419     }
00420 
00421     // tagTracer is the value at nodes
00422     Tag tagTracer = 0;
00423     std::string tag_name( "Tracer" );
00424     rval = mb.tag_get_handle( tag_name.c_str(), 1, MB_TYPE_DOUBLE, tagTracer, MB_TAG_DENSE | MB_TAG_CREAT );MB_CHK_ERR( rval );
00425 
00426     // tagElem is the average computed at each element, from nodal values
00427     Tag tagElem = 0;
00428     std::string tag_name2( "TracerAverage" );
00429     rval = mb.tag_get_handle( tag_name2.c_str(), 1, MB_TYPE_DOUBLE, tagElem, MB_TAG_DENSE | MB_TAG_CREAT );MB_CHK_ERR( rval );
00430 
00431     // area of the euler element is fixed, store it; it is used to recompute the averages at each
00432     // time step
00433     Tag tagArea = 0;
00434     std::string tag_name4( "Area" );
00435     rval = mb.tag_get_handle( tag_name4.c_str(), 1, MB_TYPE_DOUBLE, tagArea, MB_TAG_DENSE | MB_TAG_CREAT );MB_CHK_ERR( rval );
00436 
00437     // add a field value, quasi smooth first
00438     rval = add_field_value( &mb, euler_set, rank, tagTracer, tagElem, tagArea );MB_CHK_ERR( rval );
00439 
00440     // iniVals are used for 1-norm error computation
00441     Range redEls;
00442     rval = mb.get_entities_by_dimension( euler_set, 2, redEls );MB_CHK_ERR( rval );
00443     std::vector< double > iniVals( redEls.size() );
00444     rval = mb.tag_get_data( tagElem, redEls, &iniVals[0] );MB_CHK_ERR( rval );
00445 
00446     Tag tagh = 0;
00447     std::string tag_name3( "Case1" );
00448     rval = mb.tag_get_handle( tag_name3.c_str(), 3, MB_TYPE_DOUBLE, tagh, MB_TAG_DENSE | MB_TAG_CREAT );MB_CHK_ERR( rval );
00449 
00450     EntityHandle out_set, lagr_set;
00451     rval = mb.create_meshset( MESHSET_SET, out_set );MB_CHK_ERR( rval );
00452     rval = mb.create_meshset( MESHSET_SET, lagr_set );MB_CHK_ERR( rval );
00453     // copy the initial mesh in the lagrangian set
00454     // initial vertices will be at the same position as euler;
00455 
00456     rval = IntxUtilsCSLAM::deep_copy_set( &mb, euler_set, lagr_set );MB_CHK_ERR( rval );
00457 
00458     Intx2MeshOnSphere worker( &mb );
00459     worker.set_radius_source_mesh( radius );
00460     worker.set_radius_destination_mesh( radius );
00461     worker.set_error_tolerance( gtol );
00462     worker.set_parallel_comm( pcomm );
00463 
00464     rval = worker.FindMaxEdges( lagr_set, euler_set );MB_CHK_ERR( rval );
00465     Range local_verts;
00466     // output also the local_verts
00467     rval = worker.build_processor_euler_boxes( euler_set, local_verts );MB_CHK_ERR( rval );
00468     // these stay fixed for one run
00469     // other things from intersection might need to change, like input blue set (departure set)
00470     // so we need also a method to clean memory
00471 
00472     for( int i = 1; i < numSteps + 1; i++ )
00473     {
00474         // time depends on i; t = i*T/numSteps: ( 0, T/numSteps, 2*T/numSteps, ..., T )
00475         // this is really just to create some plots; it is not really needed to proceed
00476         // the compute_tracer_case1 method actually computes the departure point position
00477         if( velocity )
00478         {
00479             rval = compute_velocity_case1( &mb, euler_set, tagh, rank, i );MB_CHK_ERR( rval );
00480         }
00481 
00482         // this is to actually compute concentrations at time step i, using the
00483         //  current concentrations
00484         rval =
00485             compute_tracer_case1( &mb, worker, euler_set, lagr_set, out_set, tagElem, tagArea, rank, i, local_verts );MB_CHK_ERR( rval );
00486     }
00487 
00488     // final vals and 1-norm
00489     Range::iterator iter = redEls.begin();
00490     double norm1         = 0.;
00491     int count            = 0;
00492     void* data;
00493     int j = 0;  // index in iniVals
00494     while( iter != redEls.end() )
00495     {
00496         rval = mb.tag_iterate( tagElem, iter, redEls.end(), count, data );MB_CHK_ERR( rval );
00497         double* ptrTracer = (double*)data;
00498 
00499         rval = mb.tag_iterate( tagArea, iter, redEls.end(), count, data );MB_CHK_ERR( rval );
00500         double* ptrArea = (double*)data;
00501         for( int i = 0; i < count; i++, ++iter, ptrTracer++, ptrArea++, j++ )
00502         {
00503             // double area = *ptrArea;
00504             norm1 += fabs( *ptrTracer - iniVals[j] ) * ( *ptrArea );
00505         }
00506     }
00507 
00508     double total_norm1 = 0;
00509     int mpi_err        = MPI_Reduce( &norm1, &total_norm1, 1, MPI_DOUBLE, MPI_SUM, 0, MPI_COMM_WORLD );
00510     if( MPI_SUCCESS != mpi_err ) return 1;
00511     if( 0 == rank ) std::cout << " numSteps:" << numSteps << " 1-norm:" << total_norm1 << "\n";
00512     MPI_Finalize();
00513     return 0;
00514 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines