42 #include "TestUtil.hpp"
44 const char BRIEF_DESC[] =
"Simulate a transport problem in a semi-Lagrangian formulation\n";
81 assert( count == (
int)connecVerts.
size() );
82 double* ptr_DP = (
double*)data;
90 double params[] = { M_PI, M_PI / 3, M_PI, -M_PI / 3, 0.1, 0.9, 1., 0.5 };
116 double params[] = { p1[0], p1[1], p1[2], p2[0], p2[1], p2[2], 1, 5. };
134 double params[] = { M_PI, M_PI / 3, M_PI, -M_PI / 3, 0.1, 0.9, 0.5 };
154 double local_mass = 0.;
155 while( iter != polygons.
end() )
158 double* ptr = (
double*)data;
161 double* ptrArea = (
double*)data;
162 for(
int i = 0; i < count; i++, ++iter, ptr++, ptrArea++ )
167 if( num_nodes == 0 )
return MB_FAILURE;
168 std::vector< double > nodeVals( num_nodes );
171 for(
int j = 0; j < num_nodes; j++ )
172 average += nodeVals[j];
173 average /= num_nodes;
177 std::vector< double > coords;
178 coords.resize( 3 * num_nodes );
186 local_mass += *ptrArea * average;
189 double total_mass = 0.;
190 int mpi_err = MPI_Reduce( &local_mass, &total_mass, 1, MPI_DOUBLE, MPI_SUM, 0, MPI_COMM_WORLD );
191 if( MPI_SUCCESS != mpi_err )
return MB_FAILURE;
193 if( rank == 0 ) std::cout <<
"initial total mass:" << total_mass <<
"\n";
214 assert( count == (
int)connecVerts.
size() );
215 double* ptr_velo = (
double*)data;
227 ptr_velo[0] = velo[0];
228 ptr_velo[1] = velo[1];
229 ptr_velo[2] = velo[2];
235 std::stringstream velos;
236 velos <<
"Tracer" << rank <<
"_" << tStep <<
".vtk";
286 std::stringstream departureMesh;
287 departureMesh <<
"Departure" << rank <<
"_" << tStep <<
".vtk";
290 std::stringstream newTracer;
291 newTracer <<
"Tracer" << rank <<
"_" << tStep <<
".vtk";
294 std::stringstream lagr_cover;
295 lagr_cover <<
"Cover" << rank <<
"_" << tStep <<
".vtk";
305 std::stringstream intx_mesh;
306 intx_mesh <<
"Intx" << rank <<
"_" << tStep <<
".vtk";
315 std::stringstream resTrace;
317 <<
"_" << tStep - 1 <<
".h5m";
318 rval =
mb->
write_file( resTrace.str().c_str(), 0,
"PARALLEL=WRITE_PART", &euler_set, 1, &tagElem, 1 );
MB_CHK_ERR( rval );
324 std::stringstream resTrace;
326 <<
"_" << tStep <<
".h5m";
327 rval =
mb->
write_file( resTrace.str().c_str(), 0,
"PARALLEL=WRITE_PART", &euler_set, 1, &tagElem, 1 );
332 std::stringstream newIntx;
333 newIntx <<
"newIntx" << rank <<
"_" << tStep <<
".vtk";
360 if( rank == 0 ) std::cout <<
" step: " << tStep <<
"\n";
364 int main(
int argc,
char** argv )
366 MPI_Init( &argc, &argv );
367 LONG_DESC <<
"This program simulates a transport problem on a sphere"
368 " according to a benchmark from a Nair & Lauritzen paper.\n"
369 <<
"It starts with a partitioned mesh on a sphere, add a tracer, and steps through.\n"
370 <<
"The flow reverses after half time, and it should return to original "
371 "configuration, if the integration was exact. ";
375 std::string fileN = TestDir +
"unittest/mbcslam/fine4.h5m";
376 const char* filename_mesh1 = fileN.c_str();
378 opts.
addOpt<
double >(
"gtolerance,g",
"geometric absolute tolerance (used for point concidence on the sphere)",
382 opts.
addOpt< std::string >(
"input_file,i",
"input mesh file, partitioned", &
input_file );
383 std::string extra_read_opts;
384 opts.
addOpt< std::string >(
"extra_read_options,O",
"extra read options ", &extra_read_opts );
386 opts.
addOpt<
int >(
"field_type,f",
"field type-- 1: quasi-smooth; 2: smooth; 3: slotted cylinders (non-smooth)",
389 opts.
addOpt<
int >(
"num_steps,n",
"number of steps ", &
numSteps );
392 opts.
addOpt<
void >(
"write_debug_files,w",
"write debugging files during simulation ", &
writeFiles );
394 opts.
addOpt<
void >(
"write_velocity_files,v",
"Reorder mesh to group entities by partition", &
velocity );
396 opts.
addOpt<
void >(
"write_result_in_parallel,p",
"write tracer result files", &
parallelWrite );
403 std::string optsRead = std::string(
"PARALLEL=READ_PART;PARTITION=PARALLEL_PARTITION" ) +
404 std::string(
";PARALLEL_RESOLVE_SHARED_ENTS" ) + extra_read_opts;
421 std::cout <<
" case 1: use -gtol " <<
gtol <<
" -R " <<
radius <<
" -input " << filename_mesh1 <<
" -f "
423 std::cout <<
" write debug results: " << (
writeFiles ?
"yes" :
"no" ) <<
"\n";
424 std::cout <<
" write tracer in parallel: " << (
parallelWrite ?
"yes" :
"no" ) <<
"\n";
425 std::cout <<
" output velocity: " << (
velocity ?
"yes" :
"no" ) <<
"\n";
430 std::string tag_name(
"Tracer" );
435 std::string tag_name2(
"TracerAverage" );
441 std::string tag_name4(
"Area" );
450 std::vector< double > iniVals( redEls.
size() );
454 std::string tag_name3(
"Case1" );
469 worker.set_parallel_comm( pcomm );
474 rval = worker.build_processor_euler_boxes( euler_set, local_verts );
MB_CHK_ERR( rval );
479 for(
int i = 1; i <
numSteps + 1; i++ )
501 while( iter != redEls.
end() )
504 double* ptrTracer = (
double*)data;
507 double* ptrArea = (
double*)data;
508 for(
int i = 0; i < count; i++, ++iter, ptrTracer++, ptrArea++, j++ )
511 norm1 += fabs( *ptrTracer - iniVals[j] ) * ( *ptrArea );
515 double total_norm1 = 0;
516 int mpi_err = MPI_Reduce( &norm1, &total_norm1, 1, MPI_DOUBLE, MPI_SUM, 0, MPI_COMM_WORLD );
517 if( MPI_SUCCESS != mpi_err )
return 1;
518 if( 0 == rank ) std::cout <<
" numSteps:" <<
numSteps <<
" 1-norm:" << total_norm1 <<
"\n";