We add the adaptive mesh refinement to the program in this chapter. The driver (tportamr.C) is more complex than before. The user must define clustering, mesh refinement, error estimation, etc. The prolongation and restriction routines in this example are found in grid.f. To cluster the flagged points we need the routines Cluster1.C and Cluster3.C. The file tportutil.f has subroutine readinput that reads the input parameters from a file input.par. The truncation error estimate code is in the file truncation.C. The same makefile with a few additions can be used to obtain the executable code and the same command to run on the program on several processors using MPI.
We will define a flag to indicate that processor number 0 will write to stdout.
#define VizServer 0
#define H3eXGRAPH "/u5/parashar/xgraph/xgraph"
#define H3eDISPLAY "indra.ticam.utexas.edu:0.0"
We will give the pre-processor directive to include the DAGH header files.
#include "DAGH.h"
#include "DAGHCluster.h"
The application specific header files are also included.
#include "tportamr.h"
#include "tportfortran.h"
We will set the default AMR parameters. These parameters will be read in from the input file later on. The maximum level of refinement is set at 1, indicating a single grid. A buffer zone is added around the fine grid. The larger the buffer zone the longer the time interval over which the grids are stable and the time between regridding is lengthened. But greater the buffer width, the more work is needed to integrate the extra points in the buffer zone. The block width defines the minimum size of the cluster. The minimum efficiency of clustering is given by the ratio of the number of flagged points to the number of grid points.
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; // Minimum efficiency of clustering
GFTYPE Thresh = 0.0; // Truncation error threshold for regriding
INTEGER RegridEvery = 4; // How often do I regrid?
INTEGER NumIters = 0; // Number of timesteps on the base grid
We start the main function and initialize MPI.
void main (INTEGER argc, CHARACTER *argv[])
{
cout << "Initializing MPI \n";
MPI_Init (&argc, &argv); // Initialize MPI
We will declare or define the local variables used.
INTEGER nx, ny; // Number of grid points
DOUBLE dx, dy, dt; // Grid spacings
DOUBLE xmin, xmax, ymin, ymax; // Grid extent
DOUBLE cfl; // Courant factor
INTEGER niters = 0; // Number of iterations
INTEGER ml; // Maximum number of levels
INTEGER regride; // Regridding time interval
INTEGER buffw; // Buffer width
INTEGER blkw; // Minimum block width
GFTYPE thresh = 0; // Error threshold
DOUBLE mineff = 0; // Minimum efficiency
We will call the fortran subroutine readinput to read the input parameters from the file input.par. After reading the input file the AMR parameters get assigned.
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;
We will calculate the coarse grid spacing and the time step.
dx = (xmax - xmin) / (DOUBLE) (nx - 1);
dy = (ymax - ymin) / (DOUBLE) (ny - 1);
dt = dx * cfl;
We will set up a 2-dimensional grid hierarchy GH. The constant DIM is set to 2 in the header files. The arrays giving the shape of the grid and the bounding boxes are declared and initialized. The refinement factor is set at 2, which implies that the mesh spacings are halved at each finer level. The DAGHNoBoundary is a directive not to do anything on the external boundary.
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);
The data structures need to be distributed across processors.
ComposeHierarchy (GH);
MPI sets up a virtual machine having a network of processors. NUM_PROC returns the number of processors in this network. MY_PROC returns the local processor number.
INTEGER me = MY_PROC (GH);
INTEGER num = NUM_PROC (GH);
A ghost region is maintained around the distributed grid blocks for inter-grid communication and updating. The space stencil defines the width of the ghost region. The space stencil is set to 1. The time stencil is also set to 1. Storage is maintained for 3 (= 2 * t_sten + 1) time levels numbered from -1 to +1.
INTEGER t_sten = 1;
INTEGER s_sten = 1;
The GridFunction associates the values from the function u with grid points in the GridHierarchy. It defines storage for these values and establishes communication. In our example we want communication on the face of the cells only. We also want a shadow hierarchy to be used for error estimation.
GridFunction (DIM)<GFTYPE> u("u", t_sten, s_sten, GH, DAGHCommFaceOnly, DAGHHasShadow);
To reduce the amount of storage we will share storage of the values of the function u. Storage is shared between time step +1 from the current time and time step -1 from current time.
SetTimeAlias (u, 1, -1);
The boundary values are treated differently. We want to shift the values on the inner cells of the boundary to the boundary.
SetBoundaryType (u, DAGHBoundaryShift);
A similar GridFunction is defined for the function temp .
GridFunction (DIM)<GFTYPE> temp("temp", t_sten, s_sten, GH, DAGHCommFaceOnly, DAGHHasShadow);
SetTimeAlias (temp, 1, -1);
The prolongation function maps the values at the coarse grid points to the next finer level of grid points. The restriction function maps the values at the fine grid points to the next level of coarser grid points. We do this for all the distributed grid structures in the hierarchy.
foreachGF (GH, GF, DIM, GFTYPE)
SetProlongFunction (GF, (void *) &f_prolong_amr);
SetRestrictFunction (GF, (void *) &f_restrict_amr);
end_foreachGF
We will declare and initialize the local variables for the base or coarse grid.
INTEGER Level = 0;
INTEGER Time = CurrentTime (GH, Level, DAGH_Main);
INTEGER TStep = TimeStep (GH, Level, DAGH_Main);
We will initialize the main grid structure and the shadow using the Fortran initialization routine.
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
If you have xgraph you can see the result of the initialization. If you do not have xgraph leave this section commented out.
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);
We will start the main iterative loop. This loop carries out the evolution of the transport equation over all the grids at all levels. A recursive function bno_RecursiveIntegrate is used to integrate at time levels n and n-1.
INTEGER currentiter = 0;
for (currentiter=1; currentiter<=NumIters; currentiter++) {
bno_RecursiveIntegrate (GH, Level, u, temp, dt, h, cfl);
The built-in functions MaxVal and MinVal return the maximum and minimum values. To obtain the result we need to define two local variables for the current time and the time step.
Time = CurrentTime (GH, Level, DAGH_Main);
TStep = TimeStep (GH, Level, DAGH_Main);
cout << me << " -> Iteration: " << currentiter
<< " MaxVal: " << MaxVal (u, Time, Level, DAGH_Main)
<< " MinVal: " << MinVal (u, Time, Level, DAGH_Main) << endl;
If you have xgraph running you can display the results of the evolution, otherwise leave this section commented out.
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);
This is the end of the loop. We will close the MPI environment, and that will be the end of the main function.
DAGHMPI_Finalize();
void bno_RecursiveIntegrate (
GridHierarchy &GH,
INTEGER const Level,
GridFunction (DIM) <GFTYPE> &u,
GridFunction (DIM) <GFTYPE> &temp,
DOUBLE &dt, DOUBLE const &cfl)
We need to get the processor information. me is the id of the current processor and num is the number of processors that the data structures are distributed over.
INTEGER me = MY_PROC (GH);
INTEGER num = NUM_PROC (GH);
We select the number of iterations to run based on the current grid level.
INTEGER NoIterations = 0;
if (Level == 0)
NoIterations = 1;
else
NoIterations = RefineFactor (GH);
We take the recursive steps. First we declare and initialize the local variables.
for (INTEGER i = 0; i < NoIterations; i++) {
INTEGER Time = CurrentTime (GH, Level, DAGH_Main);
INTEGER TStep = TimeStep (GH, Level, DAGH_Main);
We will check if it is time for re-gridding. We get the truncation error from the user supplied error estimation routine wave_trunc_est and use the error threshold read in to determine if it is re-gridding time. This check and regridding is done in the function bno_RegridSystemAbove.
if (Level < MaxLevel (GH) &&
StepsTaken (GH, Level, DAGH_Main) % RegridEvery == 0 &&
StepsTaken (GH, Level, DAGH_Main) != 0){
wave_trunc_est (GH, t, l, u, temp);
bno_RegridSystemAbove (GH, temp, Level, Thresh);
}
We will take one step on DAGH_Main hierarchy. The sequence of steps are similar to the ones we saw for the uni-grid sequential code. We begin by defining the grid spacing and the time step. We call the predictor subroutine. We update the boundary values of the various grids using the values in the ghost regions. We call the corrector subroutine. We do another update of the boundary using the values in the ghost region. Finally we shift the values just interior to the boundary to the boundary.
DOUBLE dagh_dx = DeltaX (GH, DAGH_X, Level, DAGH_Main);
DOUBLE dagh_dy = DeltaY (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);
We take a step on the shadow hierarchy. The shadow hierarchy is used to estimate the local truncation error. We define the grid spacings and time step. Then we update the shadow hierarchy from the main. The next sequence of steps are similar to the steps we took on the main hierarchy as discussed above.
if (StepsTaken(GH,Level,DAGH_Main) % RefineFactor(GH) == 0) {
DOUBLE dagh_dx = DeltaX (GH, DAGH_X, Level, DAGH_Shadow);
DOUBLE dagh_dy = DeltaY (GH, DAGH_Y, Level, DAGH_Shadow);
DOUBLE dagh_dt = dagh_dx * cfl;
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);
If we are not at the finest level go to the next higher level and do a recursive integration at that level.
if (Level < FineLevel (GH)) {
bno_RecursiveIntegrate (GH, Level+1, u, temp, dt, cfl);
}
We will move from the current time to the previous time.
CycleTimeLevels (u, Level);
The time step at the current level needs to be incremented.
IncrCurrentTime (GH, Level, DAGH_Main)
Finally we will map the function values at the fine grid points to the next level of coarser grid points.
if (Level < FineLevel (GH)) {
Time = CurrentTime (GH, Level, DAGH_Main);
INTEGER myargc = 0; GFTYPE myargs[1];
Restrict (u, Time, Level+1, t, Level, myargs, myargc, DAGH_Main);
Sync (u, Time, Level, DAGH_Main);
}
void bno_RegridSystemAbove (
GridHierarchy &GH,
GridFunction (DIM) <GFTYPE> &u,
INTEGER const Level, GFTYPE const thresh) {
INTEGER me = MY_PROC (GH);
INTEGER num = NUM_PROC (GH);
We will compute the levels to regrid from the finest level down to the current level.
INTEGER flev = FineLevel (GH);
INTEGER lev = 0;
if (flev == 0) lev = 0;
else if (flev == MaxLev - 1) lev = flev - 1;
else lev = flev;
We declare a local bounding box list. Then we loop from the the finest level to the current level. The flagged points are clustered using a 3-dimensional clustering algorithm. A refined grid is overlaid over the flagged points that have been clustered.
BBoxList bblist;
for (; lev >= Level; lev--) {
INTEGER Time = CurrentTime (GH, lev, DAGH_Main);
INTEGER TStep = TimeStep (GH, lev, DAGH_Main);
DAGHCluster3d (u, Time, lev,
BlockWidth, BufferWidth,
MinEfficiency, thresh,
bblist, bblist, DAGH_Main);
Refine (GH, bblist, lev);
}
Finally the coarse and fine grid blocks are redistributed over the various processors.
RecomposeHierarchy (GH);