tportseq.C

//----------------------------------------------------------------------
//
// DAGH Driver for the Transport Application
// Unigrid Sequential Operation
//
//----------------------------------------------------------------------
#include <DAGH.h>
#include "tportfortran.h"

// If you have xgraph running on your system insert the path name
// to its directory. Redefine H3eDISPLAY for the machine you are
// running your program. Otherwise leave the comments in place.
/*
#define H3eXGRAPH "/u5/parashar/bin/xgraph"
#define H3eDISPLAY "indra.ticam.utexas.edu:0.0"
*/

void main(INTEGER argc, CHARACTER *argv[])  
{

  const INTEGER nx = 32;
  const INTEGER ny = 32;

  const INTEGER iterations = 2*nx;

  // Set up the Grid Structure
  INTEGER shape[2];
  DOUBLE bb[2*2];

  shape[0] = nx;
  bb[0] = 0.0;
  bb[1] = 1.0;

  shape[1] = ny;
  bb[2] = 0.0;
  bb[3] = 1.0;

  GridHierarchy GH(2,DAGHVertexCentered,1);
  SetBaseGrid (GH,bb,shape);
  SetBoundaryWidth(GH,1);
  ComposeHierarchy(GH);

  // Setup the GridFunctions....
  INTEGER t_sten = 1, s_sten = 1;
  GridFunction(2)<DOUBLE> u("u",t_sten,s_sten,GH);
  SetTimeAlias(u,1,-1);
  SetBoundaryType(u,DAGHBoundaryShift);

  INTEGER me = MY_PROC(GH);
  INTEGER num = NUM_PROC(GH);

  INTEGER Level = 0;
  INTEGER Time = CurrentTime(GH,Level);
  INTEGER TStep = TimeStep(GH,Level);

  DOUBLE dt, dx, dy, tzero = 0.0;
  const DOUBLE dtfac = 0.2;
  const DOUBLE sqrt3 = 1.7321;

  dx = 1.0/(nx - 1.0);
  dy = 1.0/(ny - 1.0);
  dt = dtfac *dx / sqrt3;

  forall(u,Time,Level,c)
    f_initial(FA(u(Time,Level,c,DAGH_Main)),&tzero,&dx,&dy);
  end_forall

  INTEGER currentiter = 0;

  cout << me << ": Iteration: " << currentiter
             << " Max: " << MaxVal (u,Time,Level,DAGH_Main)
             << " Min: " << MinVal (u,Time,Level,DAGH_Main)
             << " Norm: " << Norm2 (u,Time,Level,DAGH_Main) 
	     << endl;

// Remove the comments, if you have xgraph running.
/*
  BBox viewbbX (2, 0, ny/2, nx-1, ny/2, 1);
  BBox viewbbY (2, nx/2, 0, nx/2, ny-1, 1);

  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);
*/

  for (currentiter=1; currentiter < iterations; currentiter++) {

      cout << me << ": Iteration " << currentiter 
      << endl;

      // The evolved level will be at the new time
      forall(u,Time,Level,c)
        u(Time-TStep,Level,c,DAGH_Main) = u(Time,Level,c,DAGH_Main);
      end_forall

      forall(u,Time,Level,c)
        f_evolveP(FA(u(Time,Level,c,DAGH_Main)),
                  FA(u(Time-TStep,Level,c,DAGH_Main)),
                  &dt,&dx,&dy);
      end_forall

      Sync(u,Time,Level,DAGH_Main);

      forall(u,Time,Level,c)
        f_evolveC(FA(u(Time,Level,c,DAGH_Main)),
                  FA(u(Time-TStep,Level,c,DAGH_Main)),
                  &dt,&dx,&dy);
      end_forall

      Sync(u,Time,Level,DAGH_Main);

      BoundaryUpdate(u, Time, Level, DAGH_Main);
      cout << me << ": Iteration: " << currentiter
                 << " Max: " << MaxVal (u,Time,Level,DAGH_Main)
                 << " Min: " << MinVal (u,Time,Level,DAGH_Main)
                 << " Norm: " << Norm2 (u,Time,Level,DAGH_Main) 
		 << endl; 

// Remove the comments if you have xgraph running.
/*
      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);
        }
*/

  }


  cout << "Bye !" << endl;
}