MOAB: Mesh Oriented datABase
(version 5.2.1)
|
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 <math.h> // 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 >ol ); 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", ¶llelWrite ); 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 }