diffusion.cpp File Reference

#include <iostream>
#include <sstream>
#include <ctime>
#include <cstdlib>
#include <cstdio>
#include <cstring>
#include <cmath>
#include "moab/Core.hpp"
#include "moab/Interface.hpp"
#include "moab/IntxMesh/Intx2MeshOnSphere.hpp"
#include "moab/ProgOptions.hpp"
#include "MBTagConventions.hpp"
#include "moab/ParallelComm.hpp"
#include "moab/IntxMesh/IntxUtils.hpp"
#include "IntxUtilsCSLAM.hpp"
#include "TestUtil.hpp"

Include dependency graph for diffusion.cpp:

Go to the source code of this file.

Functions
ErrorCode	add_field_value (Interface *mb, EntityHandle euler_set, int rank, Tag &tagTracer, Tag &tagElem, Tag &tagArea)
ErrorCode	compute_velocity_case1 (Interface *mb, EntityHandle euler_set, Tag &tagh, int rank, int tStep)
ErrorCode	compute_tracer_case1 (Interface *mb, Intx2MeshOnSphere &worker, EntityHandle euler_set, EntityHandle lagr_set, EntityHandle out_set, Tag &tagElem, Tag &tagArea, int rank, int tStep, Range &connecVerts)
int	main (int argc, char **argv)

Variables
const char	BRIEF_DESC [] = "Simulate a transport problem in a semi-Lagrangian formulation\n"
std::ostringstream	LONG_DESC
double	gtol = 1.e-9
double	radius = 1.
int	numSteps = 3
double	T = 1
int	case_number = 1
bool	writeFiles = false
bool	parallelWrite = false
bool	velocity = false
int	field_type = 1

Function Documentation

add_field_value()

ErrorCode add_field_value	(	Interface *	mb,
		EntityHandle	euler_set,
		int	rank,
		Tag &	tagTracer,
		Tag &	tagElem,
		Tag &	tagArea
	)

Definition at line 61 of file diffusion.cpp.

{
 /*
 * get all plys first, then vertices, then move them on the surface of the sphere
 * radius is 1., most of the time
 *
 */
 Range polygons;
 ErrorCode rval = mb->get_entities_by_dimension( euler_set, 2, polygons );
 if( MB_SUCCESS != rval ) return rval;
 
 Range connecVerts;
 rval = mb->get_connectivity( polygons, connecVerts );
 if( MB_SUCCESS != rval ) return rval;
 
 void* data; // pointer to the LOC in memory, for each vertex
 int count;
 
 rval = mb->tag_iterate( tagTracer, connecVerts.begin(), connecVerts.end(), count, data );MB_CHK_ERR( rval );
 // here we are checking contiguity
 assert( count == (int)connecVerts.size() );
 double* ptr_DP = (double*)data;
 // lambda is for longitude, theta for latitude
 // param will be: (la1, te1), (la2, te2), b, c; hmax=1, r=1/2
 // nondivergent flow, page 5, case 1, (la1, te1) = (M_PI, M_PI/3)
 // (la2, te2) = (M_PI, -M_PI/3)
 // la1, te1 la2 te2 b c hmax r
 if( field_type == 1 ) // quasi smooth
 {
 double params[] = { M_PI, M_PI / 3, M_PI, -M_PI / 3, 0.1, 0.9, 1., 0.5 };
 for( Range::iterator vit = connecVerts.begin(); vit != connecVerts.end(); ++vit )
 {
 EntityHandle oldV = *vit;
 CartVect posi;
 rval = mb->get_coords( &oldV, 1, &( posi[0] ) );MB_CHK_ERR( rval );
 
 moab::IntxUtils::SphereCoords sphCoord = moab::IntxUtils::cart_to_spherical( posi );
 
 ptr_DP[0] = IntxUtilsCSLAM::quasi_smooth_field( sphCoord.lon, sphCoord.lat, params );
 ;
 
 ptr_DP++; // increment to the next node
 }
 }
 else if( 2 == field_type ) // smooth
 {
 CartVect p1, p2;
 moab::IntxUtils::SphereCoords spr;
 spr.R = 1;
 spr.lat = M_PI / 3;
 spr.lon = M_PI;
 p1 = moab::IntxUtils::spherical_to_cart( spr );
 spr.lat = -M_PI / 3;
 p2 = moab::IntxUtils::spherical_to_cart( spr );
 // x1, y1, z1, x2, y2, z2, h_max, b0
 double params[] = { p1[0], p1[1], p1[2], p2[0], p2[1], p2[2], 1, 5. };
 for( Range::iterator vit = connecVerts.begin(); vit != connecVerts.end(); ++vit )
 {
 EntityHandle oldV = *vit;
 CartVect posi;
 rval = mb->get_coords( &oldV, 1, &( posi[0] ) );MB_CHK_ERR( rval );
 
 moab::IntxUtils::SphereCoords sphCoord = moab::IntxUtils::cart_to_spherical( posi );
 
 ptr_DP[0] = IntxUtilsCSLAM::smooth_field( sphCoord.lon, sphCoord.lat, params );
 ;
 
 ptr_DP++; // increment to the next node
 }
 }
 else if( 3 == field_type ) // slotted
 {
 // la1, te1, la2, te2, b, c, r
 double params[] = { M_PI, M_PI / 3, M_PI, -M_PI / 3, 0.1, 0.9, 0.5 }; // no h_max
 for( Range::iterator vit = connecVerts.begin(); vit != connecVerts.end(); ++vit )
 {
 EntityHandle oldV = *vit;
 CartVect posi;
 rval = mb->get_coords( &oldV, 1, &( posi[0] ) );MB_CHK_ERR( rval );
 
 moab::IntxUtils::SphereCoords sphCoord = moab::IntxUtils::cart_to_spherical( posi );
 
 ptr_DP[0] = IntxUtilsCSLAM::slotted_cylinder_field( sphCoord.lon, sphCoord.lat, params );
 ;
 
 ptr_DP++; // increment to the next node
 }
 }
 
 // add average value for quad/polygon (average corners)
 // do some averages
 
 Range::iterator iter = polygons.begin();
 double local_mass = 0.; // this is total mass on one proc
 while( iter != polygons.end() )
 {
 rval = mb->tag_iterate( tagElem, iter, polygons.end(), count, data );MB_CHK_ERR( rval );
 double* ptr = (double*)data;
 
 rval = mb->tag_iterate( tagArea, iter, polygons.end(), count, data );MB_CHK_ERR( rval );
 double* ptrArea = (double*)data;
 for( int i = 0; i < count; i++, ++iter, ptr++, ptrArea++ )
 {
 const moab::EntityHandle* conn = NULL;
 int num_nodes = 0;
 rval = mb->get_connectivity( *iter, conn, num_nodes );MB_CHK_ERR( rval );
 if( num_nodes == 0 ) return MB_FAILURE;
 std::vector< double > nodeVals( num_nodes );
 double average = 0.;
 rval = mb->tag_get_data( tagTracer, conn, num_nodes, &nodeVals[0] );MB_CHK_ERR( rval );
 for( int j = 0; j < num_nodes; j++ )
 average += nodeVals[j];
 average /= num_nodes;
 *ptr = average;
 
 // now get area
 std::vector< double > coords;
 coords.resize( 3 * num_nodes );
 rval = mb->get_coords( conn, num_nodes, &coords[0] );MB_CHK_ERR( rval );
 
 moab::IntxAreaUtils sphAreaUtils;
 *ptrArea = sphAreaUtils.area_spherical_polygon( &coords[0], num_nodes, radius );
 
 // we should have used some
 // total mass:
 local_mass += *ptrArea * average;
 }
 }
 double total_mass = 0.;
 int mpi_err = MPI_Reduce( &local_mass, &total_mass, 1, MPI_DOUBLE, MPI_SUM, 0, MPI_COMM_WORLD );
 if( MPI_SUCCESS != mpi_err ) return MB_FAILURE;
 
 if( rank == 0 ) std::cout << "initial total mass:" << total_mass << "\n";
 
 // now we can delete the tags? not yet
 return MB_SUCCESS;
}

References moab::IntxAreaUtils::area_spherical_polygon(), moab::Range::begin(), moab::IntxUtils::cart_to_spherical(), moab::Range::end(), ErrorCode, field_type, moab::Core::get_connectivity(), moab::Core::get_coords(), moab::Core::get_entities_by_dimension(), moab::IntxUtils::SphereCoords::lat, moab::IntxUtils::SphereCoords::lon, mb, MB_CHK_ERR, MB_SUCCESS, IntxUtilsCSLAM::quasi_smooth_field(), moab::IntxUtils::SphereCoords::R, radius, moab::Range::size(), IntxUtilsCSLAM::slotted_cylinder_field(), IntxUtilsCSLAM::smooth_field(), moab::IntxUtils::spherical_to_cart(), moab::Core::tag_get_data(), and moab::Core::tag_iterate().

Referenced by main().

compute_tracer_case1()

ErrorCode compute_tracer_case1	(	Interface *	mb,
		Intx2MeshOnSphere &	worker,
		EntityHandle	euler_set,
		EntityHandle	lagr_set,
		EntityHandle	out_set,
		Tag &	tagElem,
		Tag &	tagArea,
		int	rank,
		int	tStep,
		Range &	connecVerts
	)

Definition at line 242 of file diffusion.cpp.

{
 ErrorCode rval;
 EntityHandle dum = 0;
 Tag corrTag;
 rval = mb->tag_get_handle( CORRTAGNAME, 1, MB_TYPE_HANDLE, corrTag, MB_TAG_DENSE, &dum );MB_CHK_ERR( rval );
 
 double t = tStep * T / numSteps; // numSteps is global; so is T
 double delta_t = T / numSteps; // this is global too, actually
 Range polys;
 rval = mb->get_entities_by_dimension( euler_set, 2, polys );MB_CHK_ERR( rval );
 
 // change coordinates of lagr mesh vertices
 for( Range::iterator vit = connecVerts.begin(); vit != connecVerts.end(); ++vit )
 {
 EntityHandle oldV = *vit;
 CartVect posi;
 rval = mb->get_coords( &oldV, 1, &( posi[0] ) );MB_CHK_ERR( rval );
 // Intx utils, case 1
 CartVect newPos;
 IntxUtilsCSLAM::departure_point_case1( posi, t, delta_t, newPos );
 newPos = radius * newPos; // do we need this? the radius should be 1
 EntityHandle new_vert;
 rval = mb->tag_get_data( corrTag, &oldV, 1, &new_vert );MB_CHK_ERR( rval );
 // set the new position for the new vertex
 rval = mb->set_coords( &new_vert, 1, &( newPos[0] ) );MB_CHK_ERR( rval );
 }
 
 // if in parallel, we have to move some elements to another proc, and receive other cells
 // from other procs
 // lagr and euler are preserved
 EntityHandle covering_set;
 rval = worker.create_departure_mesh_3rd_alg( lagr_set, covering_set );MB_CHK_ERR( rval );
 if( writeFiles ) // so if write
 {
 std::stringstream departureMesh;
 departureMesh << "Departure" << rank << "_" << tStep << ".vtk";
 rval = mb->write_file( departureMesh.str().c_str(), 0, 0, &lagr_set, 1 );MB_CHK_ERR( rval );
 
 std::stringstream newTracer;
 newTracer << "Tracer" << rank << "_" << tStep << ".vtk";
 rval = mb->write_file( newTracer.str().c_str(), 0, 0, &euler_set, 1 );MB_CHK_ERR( rval );
 
 std::stringstream lagr_cover;
 lagr_cover << "Cover" << rank << "_" << tStep << ".vtk";
 rval = mb->write_file( lagr_cover.str().c_str(), 0, 0, &covering_set, 1 );MB_CHK_ERR( rval );
 }
 // so we have now the departure at the previous time
 // intersect the 2 meshes (what about some checking of convexity?) for sufficient
 // small dt, it is not an issue;
 
 rval = worker.intersect_meshes( covering_set, euler_set, out_set );MB_CHK_ERR( rval );
 if( writeFiles ) // so if write
 {
 std::stringstream intx_mesh;
 intx_mesh << "Intx" << rank << "_" << tStep << ".vtk";
 rval = mb->write_file( intx_mesh.str().c_str(), 0, 0, &out_set, 1 );MB_CHK_ERR( rval );
 }
 
 // serially: lagr is the same order as euler;
 // we need to update now the tracer information on each element, based on
 // initial value and areas of each resulting polygons
 if( parallelWrite && tStep == 1 )
 {
 std::stringstream resTrace;
 resTrace << "Tracer"
 << "_" << tStep - 1 << ".h5m";
 rval = mb->write_file( resTrace.str().c_str(), 0, "PARALLEL=WRITE_PART", &euler_set, 1, &tagElem, 1 );MB_CHK_ERR( rval );
 }
 rval = worker.update_tracer_data( out_set, tagElem, tagArea );MB_CHK_ERR( rval );
 
 if( parallelWrite )
 {
 std::stringstream resTrace;
 resTrace << "Tracer"
 << "_" << tStep << ".h5m";
 rval = mb->write_file( resTrace.str().c_str(), 0, "PARALLEL=WRITE_PART", &euler_set, 1, &tagElem, 1 );
 }
 
 if( writeFiles ) // so if write
 {
 std::stringstream newIntx;
 newIntx << "newIntx" << rank << "_" << tStep << ".vtk";
 rval = mb->write_file( newIntx.str().c_str(), 0, 0, &out_set, 1 );MB_CHK_ERR( rval );
 }
 // delete now the polygons and the elements of out_set
 // also, all verts that are not in euler set or lagr_set
 Range allVerts;
 rval = mb->get_entities_by_dimension( 0, 0, allVerts );MB_CHK_ERR( rval );
 
 Range allElems;
 rval = mb->get_entities_by_dimension( 0, 2, allElems );MB_CHK_ERR( rval );
 
 // add to polys range the lagr polys
 // do not delete lagr set either, with its vertices
 rval = mb->get_entities_by_dimension( lagr_set, 2, polys );MB_CHK_ERR( rval );
 // add to the connecVerts range all verts, from all initial polys
 Range vertsToStay;
 rval = mb->get_connectivity( polys, vertsToStay );MB_CHK_ERR( rval );
 
 Range todeleteVerts = subtract( allVerts, vertsToStay );
 
 Range todeleteElem = subtract( allElems, polys );
 
 // empty the out mesh set
 rval = mb->clear_meshset( &out_set, 1 );MB_CHK_ERR( rval );
 
 rval = mb->delete_entities( todeleteElem );MB_CHK_ERR( rval );
 rval = mb->delete_entities( todeleteVerts );MB_CHK_ERR( rval );
 if( rank == 0 ) std::cout << " step: " << tStep << "\n";
 return rval;
}

References moab::Range::begin(), moab::Core::clear_meshset(), CORRTAGNAME, moab::Intx2Mesh::create_departure_mesh_3rd_alg(), moab::Core::delete_entities(), delta_t, IntxUtilsCSLAM::departure_point_case1(), moab::dum, moab::Range::end(), ErrorCode, moab::Core::get_connectivity(), moab::Core::get_coords(), moab::Core::get_entities_by_dimension(), moab::Intx2Mesh::intersect_meshes(), mb, MB_CHK_ERR, MB_TAG_DENSE, MB_TYPE_HANDLE, numSteps, parallelWrite, radius, moab::Core::set_coords(), moab::subtract(), t, T, moab::Core::tag_get_data(), moab::Core::tag_get_handle(), moab::Intx2MeshOnSphere::update_tracer_data(), moab::Core::write_file(), and writeFiles.

Referenced by main().

compute_velocity_case1()

ErrorCode compute_velocity_case1	(	Interface *	mb,
		EntityHandle	euler_set,
		Tag &	tagh,
		int	rank,
		int	tStep
	)

Definition at line 199 of file diffusion.cpp.

{
 Range polygons;
 ErrorCode rval = mb->get_entities_by_dimension( euler_set, 2, polygons );
 if( MB_SUCCESS != rval ) return rval;
 
 Range connecVerts;
 rval = mb->get_connectivity( polygons, connecVerts );
 if( MB_SUCCESS != rval ) return rval;
 
 void* data; // pointer to the velo in memory, for each vertex
 int count;
 
 rval = mb->tag_iterate( tagh, connecVerts.begin(), connecVerts.end(), count, data );MB_CHK_ERR( rval );
 // here we are checking contiguity
 assert( count == (int)connecVerts.size() );
 double* ptr_velo = (double*)data;
 // lambda is for longitude, theta for latitude
 
 for( Range::iterator vit = connecVerts.begin(); vit != connecVerts.end(); ++vit )
 {
 EntityHandle oldV = *vit;
 CartVect posi;
 rval = mb->get_coords( &oldV, 1, &( posi[0] ) );MB_CHK_ERR( rval );
 CartVect velo;
 double t = T * tStep / numSteps; //
 IntxUtilsCSLAM::velocity_case1( posi, t, velo );
 
 ptr_velo[0] = velo[0];
 ptr_velo[1] = velo[1];
 ptr_velo[2] = velo[2];
 
 // increment to the next node
 ptr_velo += 3; // to next velocity
 }
 
 std::stringstream velos;
 velos << "Tracer" << rank << "_" << tStep << ".vtk";
 rval = mb->write_file( velos.str().c_str(), 0, 0, &euler_set, 1, &tagh, 1 );MB_CHK_ERR( rval );
 
 return MB_SUCCESS;
}

References moab::Range::begin(), moab::Range::end(), ErrorCode, moab::Core::get_connectivity(), moab::Core::get_coords(), moab::Core::get_entities_by_dimension(), mb, MB_CHK_ERR, MB_SUCCESS, numSteps, moab::Range::size(), t, T, moab::Core::tag_iterate(), IntxUtilsCSLAM::velocity_case1(), and moab::Core::write_file().

Referenced by main().

main()

int main	(	int	argc,
		char **	argv
	)

Definition at line 364 of file diffusion.cpp.

{
 MPI_Init( &argc, &argv );
 LONG_DESC << "This program simulates a transport problem on a sphere"
 " according to a benchmark from a Nair & Lauritzen paper.\n"
 << "It starts with a partitioned mesh on a sphere, add a tracer, and steps through.\n"
 << "The flow reverses after half time, and it should return to original "
 "configuration, if the integration was exact. ";
 ProgOptions opts( LONG_DESC.str(), BRIEF_DESC );
 
 // read a homme file, partitioned in 16 so far
 std::string fileN = TestDir + "unittest/mbcslam/fine4.h5m";
 const char* filename_mesh1 = fileN.c_str();
 
 opts.addOpt< double >( "gtolerance,g", "geometric absolute tolerance (used for point concidence on the sphere)",
 &gtol );
 
 std::string input_file;
 opts.addOpt< std::string >( "input_file,i", "input mesh file, partitioned", &input_file );
 std::string extra_read_opts;
 opts.addOpt< std::string >( "extra_read_options,O", "extra read options ", &extra_read_opts );
 // int field_type;
 opts.addOpt< int >( "field_type,f", "field type-- 1: quasi-smooth; 2: smooth; 3: slotted cylinders (non-smooth)",
 &field_type );
 
 opts.addOpt< int >( "num_steps,n", "number of steps ", &numSteps );
 
 // bool reorder = false;
 opts.addOpt< void >( "write_debug_files,w", "write debugging files during simulation ", &writeFiles );
 
 opts.addOpt< void >( "write_velocity_files,v", "Reorder mesh to group entities by partition", &velocity );
 
 opts.addOpt< void >( "write_result_in_parallel,p", "write tracer result files", &parallelWrite );
 
 opts.parseCommandLine( argc, argv );
 
 if( !input_file.empty() ) filename_mesh1 = input_file.c_str();
 
 // read in parallel, in the "euler_set", the initial mesh
 std::string optsRead = std::string( "PARALLEL=READ_PART;PARTITION=PARALLEL_PARTITION" ) +
 std::string( ";PARALLEL_RESOLVE_SHARED_ENTS" ) + extra_read_opts;
 Core moab;
 Interface& mb = moab;
 EntityHandle euler_set;
 ErrorCode rval;
 rval = mb.create_meshset( MESHSET_SET, euler_set );MB_CHK_ERR( rval );
 
 rval = mb.load_file( filename_mesh1, &euler_set, optsRead.c_str() );MB_CHK_ERR( rval );
 
 ParallelComm* pcomm = ParallelComm::get_pcomm( &mb, 0 );
 
 rval = pcomm->check_all_shared_handles();MB_CHK_ERR( rval );
 
 int rank = pcomm->proc_config().proc_rank();
 
 if( 0 == rank )
 {
 std::cout << " case 1: use -gtol " << gtol << " -R " << radius << " -input " << filename_mesh1 << " -f "
 << field_type << " numSteps: " << numSteps << "\n";
 std::cout << " write debug results: " << ( writeFiles ? "yes" : "no" ) << "\n";
 std::cout << " write tracer in parallel: " << ( parallelWrite ? "yes" : "no" ) << "\n";
 std::cout << " output velocity: " << ( velocity ? "yes" : "no" ) << "\n";
 }
 
 // tagTracer is the value at nodes
 Tag tagTracer = 0;
 std::string tag_name( "Tracer" );
 rval = mb.tag_get_handle( tag_name.c_str(), 1, MB_TYPE_DOUBLE, tagTracer, MB_TAG_DENSE | MB_TAG_CREAT );MB_CHK_ERR( rval );
 
 // tagElem is the average computed at each element, from nodal values
 Tag tagElem = 0;
 std::string tag_name2( "TracerAverage" );
 rval = mb.tag_get_handle( tag_name2.c_str(), 1, MB_TYPE_DOUBLE, tagElem, MB_TAG_DENSE | MB_TAG_CREAT );MB_CHK_ERR( rval );
 
 // area of the euler element is fixed, store it; it is used to recompute the averages at each
 // time step
 Tag tagArea = 0;
 std::string tag_name4( "Area" );
 rval = mb.tag_get_handle( tag_name4.c_str(), 1, MB_TYPE_DOUBLE, tagArea, MB_TAG_DENSE | MB_TAG_CREAT );MB_CHK_ERR( rval );
 
 // add a field value, quasi smooth first
 rval = add_field_value( &mb, euler_set, rank, tagTracer, tagElem, tagArea );MB_CHK_ERR( rval );
 
 // iniVals are used for 1-norm error computation
 Range redEls;
 rval = mb.get_entities_by_dimension( euler_set, 2, redEls );MB_CHK_ERR( rval );
 std::vector< double > iniVals( redEls.size() );
 rval = mb.tag_get_data( tagElem, redEls, &iniVals[0] );MB_CHK_ERR( rval );
 
 Tag tagh = 0;
 std::string tag_name3( "Case1" );
 rval = mb.tag_get_handle( tag_name3.c_str(), 3, MB_TYPE_DOUBLE, tagh, MB_TAG_DENSE | MB_TAG_CREAT );MB_CHK_ERR( rval );
 
 EntityHandle out_set, lagr_set;
 rval = mb.create_meshset( MESHSET_SET, out_set );MB_CHK_ERR( rval );
 rval = mb.create_meshset( MESHSET_SET, lagr_set );MB_CHK_ERR( rval );
 // copy the initial mesh in the lagrangian set
 // initial vertices will be at the same position as euler;
 
 rval = IntxUtilsCSLAM::deep_copy_set( &mb, euler_set, lagr_set );MB_CHK_ERR( rval );
 
 Intx2MeshOnSphere worker( &mb );
 worker.set_radius_source_mesh( radius );
 worker.set_radius_destination_mesh( radius );
 worker.set_error_tolerance( gtol );
 worker.set_parallel_comm( pcomm );
 
 rval = worker.FindMaxEdges( lagr_set, euler_set );MB_CHK_ERR( rval );
 Range local_verts;
 // output also the local_verts
 rval = worker.build_processor_euler_boxes( euler_set, local_verts );MB_CHK_ERR( rval );
 // these stay fixed for one run
 // other things from intersection might need to change, like input blue set (departure set)
 // so we need also a method to clean memory
 
 for( int i = 1; i < numSteps + 1; i++ )
 {
 // time depends on i; t = i*T/numSteps: ( 0, T/numSteps, 2*T/numSteps, ..., T )
 // this is really just to create some plots; it is not really needed to proceed
 // the compute_tracer_case1 method actually computes the departure point position
 if( velocity )
 {
 rval = compute_velocity_case1( &mb, euler_set, tagh, rank, i );MB_CHK_ERR( rval );
 }
 
 // this is to actually compute concentrations at time step i, using the
 // current concentrations
 rval =
 compute_tracer_case1( &mb, worker, euler_set, lagr_set, out_set, tagElem, tagArea, rank, i, local_verts );MB_CHK_ERR( rval );
 }
 
 // final vals and 1-norm
 Range::iterator iter = redEls.begin();
 double norm1 = 0.;
 int count = 0;
 void* data;
 int j = 0; // index in iniVals
 while( iter != redEls.end() )
 {
 rval = mb.tag_iterate( tagElem, iter, redEls.end(), count, data );MB_CHK_ERR( rval );
 double* ptrTracer = (double*)data;
 
 rval = mb.tag_iterate( tagArea, iter, redEls.end(), count, data );MB_CHK_ERR( rval );
 double* ptrArea = (double*)data;
 for( int i = 0; i < count; i++, ++iter, ptrTracer++, ptrArea++, j++ )
 {
 // double area = *ptrArea;
 norm1 += fabs( *ptrTracer - iniVals[j] ) * ( *ptrArea );
 }
 }
 
 double total_norm1 = 0;
 int mpi_err = MPI_Reduce( &norm1, &total_norm1, 1, MPI_DOUBLE, MPI_SUM, 0, MPI_COMM_WORLD );
 if( MPI_SUCCESS != mpi_err ) return 1;
 if( 0 == rank ) std::cout << " numSteps:" << numSteps << " 1-norm:" << total_norm1 << "\n";
 MPI_Finalize();
 return 0;
}

References add_field_value(), ProgOptions::addOpt(), moab::Range::begin(), BRIEF_DESC, moab::ParallelComm::check_all_shared_handles(), compute_tracer_case1(), compute_velocity_case1(), moab::Core::create_meshset(), IntxUtilsCSLAM::deep_copy_set(), moab::Range::end(), ErrorCode, field_type, moab::Intx2Mesh::FindMaxEdges(), moab::Core::get_entities_by_dimension(), moab::ParallelComm::get_pcomm(), gtol, input_file, moab::Core::load_file(), LONG_DESC, mb, MB_CHK_ERR, MB_TAG_CREAT, MB_TAG_DENSE, MB_TYPE_DOUBLE, MESHSET_SET, numSteps, parallelWrite, ProgOptions::parseCommandLine(), moab::ParallelComm::proc_config(), moab::ProcConfig::proc_rank(), radius, moab::Intx2Mesh::set_error_tolerance(), moab::Intx2MeshOnSphere::set_radius_destination_mesh(), moab::Intx2MeshOnSphere::set_radius_source_mesh(), moab::Range::size(), moab::Core::tag_get_data(), moab::Core::tag_get_handle(), moab::Core::tag_iterate(), velocity, and writeFiles.

Variable Documentation

BRIEF_DESC

const char BRIEF_DESC[] = "Simulate a transport problem in a semi-Lagrangian formulation\n"

Definition at line 44 of file diffusion.cpp.

Referenced by main().

case_number

int case_number = 1

Definition at line 54 of file diffusion.cpp.

field_type

int field_type = 1

Definition at line 58 of file diffusion.cpp.

Referenced by add_field_value(), and main().

gtol

double gtol = 1.e-9

Definition at line 49 of file diffusion.cpp.

Referenced by main().

LONG_DESC

std::ostringstream LONG_DESC

Definition at line 45 of file diffusion.cpp.

Referenced by main().

numSteps

int numSteps = 3

Definition at line 52 of file diffusion.cpp.

Referenced by compute_tracer_case1(), compute_velocity_case1(), and main().

parallelWrite

bool parallelWrite = false

Definition at line 56 of file diffusion.cpp.

Referenced by compute_tracer_case1(), and main().

radius

double radius = 1.

Definition at line 50 of file diffusion.cpp.

Referenced by add_field_value(), compute_tracer_case1(), and main().

T

double T = 1

Definition at line 53 of file diffusion.cpp.

Referenced by moab::SysUtil::byteswap(), moab::SmoothCurve::compute_control_points_on_boundary_edges(), moab::SmoothCurve::compute_tangents_for_each_edge(), moab::SmoothFace::compute_tangents_for_each_edge(), compute_tracer_case1(), compute_velocity_case1(), moab::connectivity_match(), IntxUtilsCSLAM::departure_point_case1(), get_vartag_data(), ProgOptions::getReqArg(), moab::SmoothFace::init_bezier_edge(), moab::ScdNCHelper::kji_to_jik(), moab::UcdNCHelper::kji_to_jik_stride(), moab::TagTypeEqual< T >::operator()(), moab::TagTypeLess< T >::operator()(), moab::TagOneTypeEqual< T >::operator()(), moab::TagOneTypeLess< T >::operator()(), moab::TagVarTypeEqual< T >::operator()(), moab::TagVarTypeLess< T >::operator()(), moab::PACK(), parse_value(), parse_values_typed(), moab::permute_this(), moab::common_tree::print_vector(), moab::ParallelComm::reduce(), moab::WriteUtilIface::reorder(), moab::rev_permute_this(), tensor_c3(), tensor_r3(), moab::UNPACK(), IntxUtilsCSLAM::velocity_case1(), moab::WriteVtk::write_data(), and moab::WriteVtk::write_tag().

velocity

bool velocity = false

Definition at line 57 of file diffusion.cpp.

Referenced by main().

writeFiles

bool writeFiles = false

Definition at line 55 of file diffusion.cpp.

Referenced by compute_tracer_case1(), and main().