//----------------------------------------------------------------------
//----------------------------------------------------------------------
// AMR Driver for the transport equation
//----------------------------------------------------------------------
//----------------------------------------------------------------------
//----------------------------------------------------------------------
// Flags
//----------------------------------------------------------------------
#define VizServer 0 // Who will write to stdout
// Remove the comments if you have XGRAPH running.
/*
#define H3eXGRAPH "/u5/parashar/xgraph/xgraph"
#define H3eDISPLAY "telluride.ticam.utexas.edu:0.0"
*/
//----------------------------------------------------------------------
// DAGH includes
//----------------------------------------------------------------------
#include "DAGH.h"
#include "DAGHCluster.h"
//----------------------------------------------------------------------
// Application Includes
//----------------------------------------------------------------------
#include "tportamr.h"
#include "tportfortran.h"
//----------------------------------------------------------------------
// AMR parameters
//----------------------------------------------------------------------
INTEGER MaxLev = 1; // Maximum level of refinement
INTEGER BufferWidth = 1; // Buffer width used by clusterer
INTEGER BlockWidth = 1; // Minimum granularity used by clusterer
DOUBLE MinEfficiency = 0.7; // Minumum efficiency of clustering
GFTYPE Thresh = 0.0; // Truncation error threshold for regriding
INTEGER RegridEvery = 4; // How oftern do I regrid ?
INTEGER NumIters = 0; // Number of timesteps on the base grid
INTEGER OutEvery = 0; //
void main(INTEGER argc, CHARACTER *argv[])
{
//---------------------------------------------------------------
// If I using MPI initialize it...
//---------------------------------------------------------------
cout << "Initializing MPI\n";
MPI_Init( &argc, &argv ); // Initialize MPI
//---------------------------------------------------------------
// Local variable declarations
//---------------------------------------------------------------
INTEGER nx, ny;
DOUBLE dx, dy, dt; // Grid spacings
DOUBLE xmin, xmax, ymin, ymax; // Grid extent
DOUBLE cfl; // Initial data parameters
//---------------------------------------------------------------
// Read grid properties and initial data parameters from file
//---------------------------------------------------------------
INTEGER niters = 0;
INTEGER ml, regride, buffw, blkw, outevery;
GFTYPE thresh = 0;
DOUBLE mineff = 0;
INTEGER d = 0;
cout << "****Reading in parameters****" << endl << flush;
f_readinput(&nx, &ny, // Number of grid points
&xmin, &xmax, // X Coords range
&ymin, &ymax, // Y Coords range
&cfl, // Courant factor
&niters, // Number of iterations
&ml, // Max levels
&thresh, // Threshold
®ride, // Regrid every
&mineff, // Min efficiency
&buffw, // Buffer width
&blkw, // Min block width
&outevery); // Output every time step
MaxLev = ml;
RegridEvery = regride;
BufferWidth = buffw;
BlockWidth = blkw;
MinEfficiency = mineff;
Thresh = thresh;
NumIters = niters;
OutEvery = outevery;
//---------------------------------------------------------------
// Calculate dx, dy and dt from inputs
//---------------------------------------------------------------
dx = (xmax - xmin) / (DOUBLE)(nx-1);
dy = (ymax - ymin) / (DOUBLE)(ny-1);
dt = cfl*dx;
//---------------------------------------------------------------
// Setup grid hierarchy
//---------------------------------------------------------------
INTEGER shape[DIM];
DOUBLE bb[2*DIM];
shape[0] = nx-1;
shape[1] = ny-1;
bb[0] = xmin ;
bb[1] = xmax ;
bb[2] = ymin ;
bb[3] = ymax ;
GridHierarchy GH(DIM,DAGHVertexCentered,MaxLev);
SetBaseGrid (GH,bb,shape);
SetRefineFactor(GH,2);
SetBoundaryWidth(GH,1);
SetBoundaryType(GH,DAGHNoBoundary);
ComposeHierarchy(GH);
//---------------------------------------------------------------
// Processor information
//---------------------------------------------------------------
INTEGER me = MY_PROC(GH);
INTEGER num = NUM_PROC(GH);
//---------------------------------------------------------------
// Declare grid functions
//---------------------------------------------------------------
INTEGER t_sten = 1;
INTEGER s_sten = 1;
GridFunction(DIM)<GFTYPE> u("u", t_sten, s_sten, GH,
DAGHCommFaceOnly, DAGHHasShadow);
SetTimeAlias(u,1,-1);
SetBoundaryType(u,DAGHBoundaryShift);
GridFunction(DIM)<GFTYPE> temp("temp", t_sten, s_sten, GH,
DAGHCommFaceOnly, DAGHHasShadow);
SetTimeAlias(temp,1,-1);
foreachGF(GH,GF,DIM,GFTYPE)
SetProlongFunction(GF, (void *) &f_prolong_amr);
SetRestrictFunction(GF, (void *) &f_restrict_amr);
end_foreachGF
//---------------------------------------------------------------
// Initialize Level = 0 (Base Grid)
//---------------------------------------------------------------
INTEGER Level = 0;
INTEGER Time = CurrentTime(GH,Level,DAGH_Main); // Current Time at Level
INTEGER TStep = TimeStep(GH,Level,DAGH_Main); // Time Step at Level
//---------------------------------------------------------------
// Set up initial data on previous step (Main & Shadow)
//---------------------------------------------------------------
GFTYPE tzero = 0.0;
forall(u,Time,Level,c)
f_initial(FA(u(Time,Level,c,DAGH_Main)),&tzero,&dx,&dy);
end_forall
forall(u,Time,Level,c)
f_initial(FA(u(Time,Level,c,DAGH_Shadow)),&tzero,&dx,&dy);
end_forall
//---------------------------------------------------------------
// Dump out the initial data
//---------------------------------------------------------------
// Remove the comments if you are using XGRAPH
/*
const int s = StepSize(GH,Level,DAGH_Main);
const BBox gbb = GH.glbbbox();
const int snx = gbb.upper(0);
const int sny = gbb.upper(1);
BBox viewbbX(2,0,sny/2,snx,sny/2,s);
BBox viewbbY(2,snx/2,0,snx/2,sny,s);
View(u,Time,Level,DAGH_X,viewbbX,0,DAGHViz_XGRAPH,
H3eXGRAPH,H3eDISPLAY,DAGH_Main);
View(u,Time,Level,DAGH_Y,viewbbY,0,DAGHViz_XGRAPH,
H3eXGRAPH,H3eDISPLAY,DAGH_Main);
*/
//---------------------------------------------------------------
// Begin Evolution
//---------------------------------------------------------------
if (me == VizServer)
cout << "Starting main loop" << endl << flush;
INTEGER currentiter = 0;
for (currentiter=1; currentiter < NumIters; currentiter++) {
bno_RecursiveIntegrate(GH, Level, u, temp, dt, cfl);
Time = CurrentTime(GH,Level,DAGH_Main); // Current Time at Level
TStep = TimeStep(GH,Level,DAGH_Main); // Time Step at Level
cout << me << " -> Iteration: " << currentiter
<< " MaxVal: " << MaxVal(u,Time,Level,DAGH_Main)
<< " MinVal: " << MinVal(u,Time,Level,DAGH_Main)
<< endl;
// Remove the comments if you are running XGRAPH
/*
if (currentiter%10 == 0) {
View(u,Time,Level,DAGH_X,viewbbX,0,DAGHViz_XGRAPH,
H3eXGRAPH,H3eDISPLAY,DAGH_Main);
View(u,Time,Level,DAGH_Y,viewbbY,0,DAGHViz_XGRAPH,
H3eXGRAPH,H3eDISPLAY,DAGH_Main);
}
*/
} // End loop over number of iterations
cout << "Finalizing MPI\n";
DAGHMPI_Finalize();
} // End main
//----------------------------------------------------------------------
//
// Recursive Integrate
//
//----------------------------------------------------------------------
void bno_RecursiveIntegrate(
GridHierarchy &GH,
INTEGER const Level,
GridFunction(DIM)<GFTYPE> &u,
GridFunction(DIM)<GFTYPE> &temp,
DOUBLE &dt, DOUBLE const &cfl)
{
//---------------------------------------------------------------
// Processor Info
//---------------------------------------------------------------
INTEGER me = MY_PROC(GH);
INTEGER num = NUM_PROC(GH);
//---------------------------------------------------------------
// Select number of iterations to run based on current grid level
//---------------------------------------------------------------
INTEGER NoIterations = 0;
if (Level==0)
NoIterations = 1;
else
NoIterations = RefineFactor(GH);
//---------------------------------------------------------------
// Take recursive steps..
//---------------------------------------------------------------
for (INTEGER i = 0; i < NoIterations; i++) {
INTEGER Time = CurrentTime(GH,Level,DAGH_Main);
INTEGER TStep = TimeStep(GH,Level,DAGH_Main);
//---------------------------------------------------------------
// Is it regridding time?
//---------------------------------------------------------------
if (Level < MaxLevel(GH) &&
StepsTaken(GH,Level,DAGH_Main)%RegridEvery == 0 &&
StepsTaken(GH,Level,DAGH_Main) != 0) {
// Truncation error estimation using DAGH_Main
wave_trunc_est(GH,Time,Level,u,temp);
// Regrid using truncation error estimate
bno_RegridSystemAbove(GH,temp,Level,Thresh);
} // End if statement Level...
//---------------------------------------------------------------
// Take a step on the DAGH_Main hierarchy
//---------------------------------------------------------------
DOUBLE dagh_dx = DeltaX(GH,DAGH_X,Level,DAGH_Main);
DOUBLE dagh_dy = DeltaX(GH,DAGH_Y,Level,DAGH_Main);
DOUBLE dagh_dt = dagh_dx * cfl;
forall(u,Time,Level,c)
f_evolveP(FA(u(Time+TStep,Level,c,DAGH_Main)),
FA(u(Time,Level,c,DAGH_Main)),
&dagh_dt,&dagh_dx,&dagh_dy);
end_forall
Sync(u,Time,Level,DAGH_Main);
forall(u,Time,Level,c)
f_evolveC(FA(u(Time+TStep,Level,c,DAGH_Main)),
FA(u(Time,Level,c,DAGH_Main)),
&dagh_dt,&dagh_dx,&dagh_dy);
end_forall
Sync(u,Time,Level,DAGH_Main);
BoundaryUpdate(u,Time,Level,DAGH_Main);
Sync(u,Time,Level,DAGH_Main);
//---------------------------------------------------------------
// Take a step on the DAGH_Shadow hierarchy
//---------------------------------------------------------------
if (StepsTaken(GH,Level,DAGH_Main) % RefineFactor(GH) == 0){
dagh_dx = DeltaX(GH,DAGH_X,Level,DAGH_Shadow);
dagh_dy = DeltaX(GH,DAGH_Y,Level,DAGH_Shadow);
dagh_dt = dagh_dx * cfl;
// First update DAGH_Shadow from DAGH_Main
forall(u,Time,Level,c)
u(Time,Level,c,DAGH_Shadow) = u(Time,Level,c,DAGH_Main);
end_forall
Sync(u,Time,Level,DAGH_Shadow);
forall(u,Time,Level,c)
f_evolveP(FA(u(Time+TStep,Level,c,DAGH_Shadow)),
FA(u(Time,Level,c,DAGH_Shadow)),
&dagh_dt,&dagh_dx,&dagh_dy);
end_forall
Sync(u,Time,Level,DAGH_Shadow);
forall(u,Time,Level,c)
f_evolveC(FA(u(Time+TStep,Level,c,DAGH_Shadow)),
FA(u(Time,Level,c,DAGH_Shadow)),
&dagh_dt,&dagh_dx,&dagh_dy);
end_forall
Sync(u,Time,Level,DAGH_Shadow);
BoundaryUpdate(u,Time,Level,DAGH_Shadow);
Sync(u,Time,Level,DAGH_Shadow);
} // Ends if statement for DAGH_Shadow
//---------------------------------------------------------------
// If we are not at the finest level then go to next level
//---------------------------------------------------------------
if (Level < FineLevel(GH)) {
// Recursive Integration of next level in hierarchy
bno_RecursiveIntegrate( GH, Level+1, u, temp, dt, cfl );
} // End if (Level < FineLevel(GH))
//---------------------------------------------------------------
// Move from current time to previous.
//---------------------------------------------------------------
CycleTimeLevels(u,Level);
//---------------------------------------------------------------
// Increment time step on current level
//---------------------------------------------------------------
IncrCurrentTime(GH,Level,DAGH_Main);
//---------------------------------------------------------------
// Restriction on DAGH_Main
//---------------------------------------------------------------
if (Level < FineLevel(GH)) {
Time = CurrentTime(GH,Level,DAGH_Main);
INTEGER myargc = 0; GFTYPE myargs[1];
Restrict(u, Time, Level+1, Time, Level, myargs, myargc, DAGH_Main);
Sync(u, Time, Level, DAGH_Main);
} // end if statement for level
} // End For loop
} // End bnoRecursiveIntegrate
//---------------------------------------------------------------
// Regrid GridHierarchy at Level
//---------------------------------------------------------------
void bno_RegridSystemAbove(GridHierarchy &GH,
GridFunction(DIM)<GFTYPE>amp; u,
INTEGER const Level, GFTYPE const thresh)
{
INTEGER me = MY_PROC(GH);
INTEGER num = NUM_PROC(GH);
//---------------------------------------------------------------
// Compute levels to regrid at - from the finest down to Level
//---------------------------------------------------------------
INTEGER flev = FineLevel(GH);
INTEGER lev = 0;
if (flev == 0) lev = 0;
else if (flev == MaxLev-1) lev = flev-1;
else lev = flev;
//---------------------------------------------------------------
// Cluster at refine
//---------------------------------------------------------------
BBoxList bblist;
for ( ; lev > Level; lev--) {
INTEGER Time = CurrentTime(GH,lev,DAGH_Main);
INTEGER TStep = TimeStep(GH,lev,DAGH_Main);
DAGHCluster2d(u, Time, lev,
BlockWidth, BufferWidth,
MinEfficiency, thresh,
bblist, bblist, DAGH_Main);
Refine(GH,bblist,lev);
} // Ends for loop over levels
//---------------------------------------------------------------
// Recompose the GridHierarchy
//---------------------------------------------------------------
RecomposeHierarchy(GH);
} // bno_RegridSystemAbove