//---------------------------------------------------------------------- 
//----------------------------------------------------------------------
//
// AMR Driver for the 2-D Wave Application
//
//----------------------------------------------------------------------
// ----------------------------------------------------------------------
// waveamr2d.C
// Program to evolve 2D wave equation using DAGH structures
// Mijan Huq and Manish Parashar
// University of Texas at Austin
//----------------------------------------------------------------------
// Program uses wave equation in first order form. That is, there are a
// set of three equations and three grid functions.
// A = D_t phi ,  B = D_x phi , C = D_y phi
// Wave equation D_t^2 phi = D_x^2 phi + D_y^2 phi is rewritten as
// D_t A = D_x B + D_y C
// D_t B = D_x A
// D_t C = D_x A
// and solved using a Leapfrog scheme.
//----------------------------------------------------------------------


#include "DAGH.h"
#include "DAGHCluster.h"
#include "DAGHIO.h"
#include "waveamr2d.h"
#include "wavefortran2d.h"
#include "bbhutil.h"

#include "math.h"

#define VizServer 0
#define WAVEXGRAPH  "/u5/parashar/xgraph/xgraph"
#define WAVEDISPLAY "zeus.ticam.utexas.edu:0.0"


char space[41];
INTEGER MaxLev;

DOUBLE smoothval = 0.5;

// AMR parameters
DOUBLE thresh;    // Threshold for regridding

// Clustering routine prototypes
//BBoxList Cluster2(GridData(2)<short> &flag,
//                  double Efficiency,
//                  int MinWidth,
//                  int BufferWidth);

static char *indent(int const l)
{
  assert(l>=0);

  for (int i=0;i<l;i++) {
   space[2*i] = ' ';
   space[2*i+1] = ' ';
  }

  space[2*i] = '\0';

  return space;
}

void main(INTEGER argc, CHARACTER *argv[])  
{
   MPI_Init( &argc, &argv ); // Initialize MPI

// Declarations

   INTEGER  nx, ny;

   DOUBLE  h, dt;         // Grid spacings
   DOUBLE const stp = 1.0e-8;
   INTEGER const maxsteps = 25;
   DOUBLE xmin, xmax, ymin, ymax; // Grid extent
   DOUBLE cfl, amplitude, width, vx, vy;  // Initial data parameters

   
   INTEGER niters = 0;
   INTEGER ml ;

// Read grid properties and initial data parameters from file
   cout << "****Reading in parameters****" << endl << flush;
   f_readinput(&nx, &ny, &xmin, &xmax, &ymin, &ymax, &cfl,
     &amplitude, &width, &vx, &vy, &niters, &ml,&thresh);

   MaxLev = ml;
   cout << "Using MaxLev = " << MaxLev << endl;
// Calculate h and dt from inputs

   h  = (xmax - xmin) / (DOUBLE)(nx-1);

   dt = h * cfl;

// Write out properties

   cout << "****Grid Properties****" << endl << flush;
   cout << "nx = " << nx << " ny =" << ny << endl << flush;
   cout << "xmin = " << xmin << " xmax = " << xmax << endl << flush;
   cout << "ymin = " << ymin << " ymax = " << ymax << endl << flush;
   cout << "h = " << h << " dt = " << dt << endl << flush;
   cout << "MaxLev = " << ml <<  flush;

   cout << "****Initial Data: Gaussian wave packet with****" << endl << flush;
   cout << "amplitude = " << amplitude << endl << flush;
   cout << "width = " << width << endl << flush;
   cout << "Velocity (vx,vy) = (" << vx << "," << vy << ") " << endl << flush;
   cout << "threshold = " << thresh << endl << flush;

// Create grid functions for use
  
   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 ;

// Setup grid hierarchy
   GridHierarchy GH(DIM,NON_CELL_CENTERED,MaxLev);
   SetBaseGrid (GH,bb,shape);
   SetRefineFactor(GH,2);
   SetRegExtBndry(GH,1);
   ComposeHierarchy(GH);
   cout << "Maximum level = " <<MaxLevel(GH) << endl << flush;

   INTEGER  t_sten0 = 0, t_sten = 1;
   INTEGER  s_sten = 1;

//  Declare grid functions
   GridFunction(DIM)<GFTYPE> A("a", t_sten, s_sten, GH, COMM, HAS_SHADOW);
   GridFunction(DIM)<GFTYPE> B("b", t_sten, s_sten, GH, COMM, HAS_SHADOW);
   GridFunction(DIM)<GFTYPE> C("c", t_sten, s_sten, GH, COMM, HAS_SHADOW);
   GridFunction(DIM)<GFTYPE> phi("phi", t_sten, s_sten, GH, COMM, HAS_SHADOW);
   GridFunction(DIM)<GFTYPE> temp("temp", t_sten, s_sten, GH, COMM, HAS_SHADOW);

// Additional variables needed for evolution
   DOUBLE anorm, bnorm, cnorm, phinorm;

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

  // Initialize Level = 0 (Base Grid)
  INTEGER Level = 0;

  // Tell DAGH which prolongation and restriction to use
  forallGF(GH,GF,DIM,GFTYPE)
    SetProlongFunction(GF, (void *) &f_prolong_amr);
    SetRestrictFunction(GF, (void *) &f_restrict_amr);
  end_forallGF

  INTEGER t = CurrentTime(GH,Level);
  INTEGER idt = TimeStep(GH,Level);
  INTEGER l = Level;

// Set up initial data

   if(me == 0) cout << "****Setting up Initial data on MAIN ****" 
                    << endl << flush;

   // grid spacings at current level on MAIN
   DOUBLE dagh_h = DeltaX(GH,X_AXIS,l,MAIN);
   DOUBLE dagh_dt = dagh_h * cfl;

   // Load in initial data 
   forall(phi, t, l, c)
      f_initial_data(FBA(phi(t,l,c,MAIN)), FDA(phi(t,l,c,MAIN)),
         FDA(A(t,l,c,MAIN)), FDA(B(t,l,c,MAIN)), FDA(C(t,l,c,MAIN)),
         &dagh_h, &dagh_dt, &amplitude, &width, &vx, &vy, &xmin, &ymin);
   end_forall

   // Initialize the first level on the MAIN hierarchy
   if(me == 0) cout << "****Initializing first step on MAIN****" 
                    << endl << flush;

   dagh_h = DeltaX(GH,X_AXIS,l,MAIN);
   dagh_dt = dagh_h * cfl;

   // Set up the n-1 level by iteration
   initial_step (MAIN, GH, A, B, C, phi, dagh_h, dagh_dt, stp, maxsteps,l); 

   if(me == 0) cout << "****Setting up Initial data on SHADOW ****" 
                    << endl << flush;

   // grid spacings at current level on SHADOW
   dagh_h = DeltaX(GH,X_AXIS,l,SHADOW);
   dagh_dt = dagh_h * cfl;

   // Load in initial data on SHADOW
   forall(phi, t, l, c)
      f_initial_data(FBA(phi(t,l,c,SHADOW)), FDA(phi(t,l,c,SHADOW)),
         FDA(A(t,l,c,SHADOW)), FDA(B(t,l,c,SHADOW)), FDA(C(t,l,c,SHADOW)),
         &dagh_h, &dagh_dt, &amplitude, &width, &vx, &vy, &xmin, &ymin);
   end_forall

   // Initialize the first level on the SHADOW hierarchy
   if(me == 0) cout << "****Initializing first step on SHADOW****" 
                    << endl << flush;

   dagh_h = DeltaX(GH,X_AXIS,l,SHADOW);
   dagh_dt = dagh_h * cfl;

   // Set up the n-1 level by iteration on SHADOW
   initial_step (SHADOW, GH, A, B, C, phi, dagh_h, dagh_dt, stp, maxsteps,l); 


// *************************Begin Evolution**************************

  INTEGER currentiter = 0;

       // Dump out the initial data
       char  fname[40];
       idt = TimeStep(GH,l);
       sprintf(fname,"SlicePhi%d.dat",currentiter);
       writeslice(MAIN, fname, GH, phi, t, l, VizServer);

       idt = TimeStep(GH,l);
       sprintf(fname,"SliceShadowPhi%d.dat",currentiter);
       writeslice(SHADOW, fname, GH, phi, t, l, VizServer);


  if(me == 0)cout << "Starting main loop" << endl << flush;
  

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

      // Carry out evolution steps on all grids on all levels

      // At this stage only timelevels n and n-1 are defined.

      t = CurrentTime(GH,Level);
      idt = TimeStep(GH,Level);

      bno_RecursiveIntegrate(GH, Level, A, B, C, phi, temp, dt, h, cfl);


      /******* MAIN *********/

      // Calculate norms of A, B, C, phi on MAIN hierarchy
      anorm = Norm2(A, t, l, MAIN);
      bnorm = Norm2(B, t, l, MAIN);
      cnorm = Norm2(C, t, l, MAIN);

      phinorm = Norm2(phi,t, l, MAIN);


      // Loop through existing levels  and dump out data (1 processor only)
      for(INTEGER curlev = 0; curlev <= FineLevel(GH); curlev++){
         if(me == 0)cout << "Dumping Level" << curlev << endl << flush;
         INTEGER myidt = TimeStep(GH,curlev);
         INTEGER myt = CurrentTime(GH,curlev);
         sprintf(fname,"SliceLevel%dPhi%d.dat",curlev,currentiter);
         writeslice(MAIN, fname, GH, phi, myt,curlev, VizServer);
      }
      t = CurrentTime(GH,Level);
      idt = TimeStep(GH,Level);

      if(me == 0){
         cout << "MAIN norms " << anorm << " " << bnorm << " " << cnorm << " " 
              << phinorm << endl << flush;
      }

      /******* SHADOW *********/
      if(StepsTaken(GH,Level) % RefineFactor(GH) == 0){

        // Calculate norms of A, B, C, phi on SHADOW hierarchy
        anorm = Norm2(A, t, l, SHADOW);
        bnorm = Norm2(B, t, l, SHADOW);
        cnorm = Norm2(C, t, l, SHADOW);

        phinorm = Norm2(phi,t, l, SHADOW);

        // Loop through existing levels  and dump out data (1 processor only)
        for(INTEGER curlev = 0; curlev <= FineLevel(GH); curlev++){
         if(me == 0)cout << "Dumping Level" << curlev << endl << flush;
         INTEGER myidt = TimeStep(GH,curlev);
         INTEGER myt = CurrentTime(GH,curlev);
         sprintf(fname,"SliceSHLevel%dPhi%d.dat",curlev,currentiter);
         writeslice(SHADOW, fname, GH, phi, myt,curlev, VizServer);
        }

        if(me == 0){
           cout << "SHADOW norms " << anorm << " " << bnorm << " " << cnorm 
                << " " << phinorm << endl << flush;
        }
      } // End if statement for SHADOW

  } // End loop over number of iterations

  DAGHMPI_Finalize(); 
} // End main

//----------------------------------------------------------------------
// 
// Recursive Integrate assumes:
//   if Level == 0 then timelevels n and n-1 are defined and n+1 is not.
//   if Level > 0 then timelevels n+1, n and n-1 are defined on level Level-1. 
//                     This allows for time interpolations for boundary 
//                     updates in the UpdateExternalBoundary call.
//                        


void bno_RecursiveIntegrate(
                GridHierarchy &GH,
                INTEGER const Level,
                GridFunction(DIM)<GFTYPE> &A,
                GridFunction(DIM)<GFTYPE> &B,
                GridFunction(DIM)<GFTYPE> &C,
                GridFunction(DIM)<GFTYPE> &phi,
                GridFunction(DIM)<GFTYPE> &temp,
                DOUBLE &dt, DOUBLE &h, DOUBLE const &cfl)
{

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

  INTEGER l = Level;
  INTEGER t = CurrentTime(GH,Level);
  INTEGER idt = TimeStep(GH,Level);

  INTEGER NoIterations = 0;

  // Select number of iterations to run based on current grid level
  if (Level==0)
      NoIterations = 1;
  else
     NoIterations = RefineFactor(GH);

  for (INTEGER i=0;i<NoIterations;i++) {
     
    /* Is it regridding time? */


    if(MY_PROC(GH) == 0)cout << "Steps taken " <<  StepsTaken(GH,Level) <<
       " at level " <<  Level << endl << flush;
    if (Level < MaxLevel(GH) && StepsTaken(GH,Level)%REGRID_EVERY == 0 &&
        StepsTaken(GH,Level) != 0) {

       comm_service::log() << indent(Level) 
         << "********* Regriding at Level " << Level << "*********" << endl ;

       cout << "********* Regriding at Level" << Level <<"*********" <<
                      endl << flush ;

       // Truncation error estimation using MAIN
       wave_trunc_est(GH,t,l,A, B, C, temp);

       // Regrid using truncation error estimate
       bno_RegridSystemAbove(GH,temp,Level,thresh);

       cout << "*** Components after regridding: " 
            << phi.len() << endl << flush;

//     At this point make sure that MAIN and SHADOW are synchronized !!!!

  	 idt = TimeStep(GH,Level);

    } // End if statement Level...

      if(me == 0){
         INTEGER fl = FineLevel(GH); 
         INTEGER ft = CurrentTime(GH,fl);
         INTEGER fidt = TimeStep(GH,fl);
         cout << "MaxVal on MAIN(t) " << MaxVal(A, ft, fl, MAIN) << " "
                                   << MaxVal(B, ft, fl, MAIN) << " "
                                   << MaxVal(C, ft, fl, MAIN) << " "
                                   << MaxVal(phi, ft, fl, MAIN) << " "
                                   << endl << flush;
         cout << "MaxVal on SHADOW(t) " << MaxVal(A, ft, fl, SHADOW) << " "
                                   << MaxVal(B, ft, fl, SHADOW) << " "
                                   << MaxVal(C, ft, fl, SHADOW) << " "
                                   << MaxVal(phi, ft, fl, SHADOW) << " "
                                   << endl << flush;
         cout << "MaxVal on MAIN(t-1) " << MaxVal(A, ft-fidt, fl, MAIN) << " "
                                   << MaxVal(B, ft-fidt, fl, MAIN) << " "
                                   << MaxVal(C, ft-fidt, fl, MAIN) << " "
                                   << MaxVal(phi, ft-fidt, fl, MAIN) << " "
                                   << endl << flush;
         cout << "MaxVal on SHADOW(t-1) " << MaxVal(A, ft-fidt, fl, SHADOW) << " "
                                   << MaxVal(B, ft-fidt, fl, SHADOW) << " "
                                   << MaxVal(C, ft-fidt, fl, SHADOW) << " "
                                   << MaxVal(phi, ft-fidt, fl, SHADOW) << " "
                                   << endl << flush;
      }
    /*********** Take a step on the MAIN hierarchy ***********/   

      comm_service::log() << indent(Level) 
                          << "********* Stepping Main *********" << endl ;



       l = Level;
       t = CurrentTime(GH,Level);
       idt = TimeStep(GH,Level);

      DOUBLE dagh_h = DeltaX(GH,X_AXIS,l,MAIN);
      DOUBLE dagh_dt = dagh_h * cfl;

      if(MY_PROC(GH) == 0 )cout << "Leapfrog on MAIN " << endl << flush;

      // Leapfrog on MAIN

      {
      // Calculate norms of A, B, C, phi on MAIN hierarchy
      DOUBLE anorm = Norm2(A, t, l, MAIN);
      DOUBLE bnorm = Norm2(B, t, l, MAIN);
      DOUBLE cnorm = Norm2(C, t, l, MAIN);

      DOUBLE phinorm = Norm2(phi,t, l, MAIN);

      if(me == 0){
         cout << "n: At level " << Level << " BEFORE LEAPFROG " << endl << flush;
         cout << "MAIN norms " << anorm << " " << bnorm << " " << cnorm << " " 
              << phinorm << endl << flush;
      }
      // Calculate norms of A, B, C, phi on MAIN hierarchy
       anorm = Norm2(A, t-idt, l, MAIN);
       bnorm = Norm2(B, t-idt, l, MAIN);
       cnorm = Norm2(C, t-idt, l, MAIN);

       phinorm = Norm2(phi,t-idt, l, MAIN);

      if(me == 0){
         cout << "n-1: At level " << Level <<" BEFORE LEAPFROG " <<  endl << flush;
         cout << "MAIN norms " << anorm << " " << bnorm << " " << cnorm << " " 
              << phinorm << endl << flush;
      }
      }
      forall(phi, t, l, c)
         f_leapfrog(FBA(A(t,l,c,MAIN)),
           FDA(A(t+idt,l,c,MAIN)), FDA(A(t-idt,l,c,MAIN)), FDA(A(t,l,c,MAIN)),
           FDA(B(t+idt,l,c,MAIN)), FDA(B(t-idt,l,c,MAIN)), FDA(B(t,l,c,MAIN)),
           FDA(C(t+idt,l,c,MAIN)), FDA(C(t-idt,l,c,MAIN)), FDA(C(t,l,c,MAIN)),
           FDA(phi(t+idt,l,c,MAIN)), FDA(phi(t-idt,l,c,MAIN)),
           FDA(phi(t,l,c,MAIN)), &dagh_dt, &dagh_h);
      end_forall

      Sync(A, t+idt, l, MAIN);
      Sync(B, t+idt, l, MAIN);
      Sync(C, t+idt, l, MAIN);

      Sync(phi, t+idt, l, MAIN);

       l = Level;
       t = CurrentTime(GH,Level);
       idt = TimeStep(GH,Level);

      // Smooth to clean up boundary effects due to large gradients
//      if(StepsTaken(GH,Level) % REGRID_EVERY == 0 ){
      if(MY_PROC(GH) == 0)cout << "Smoothing on MAIN" << endl << flush ; 
      DOUBLE wght = 0.5;
      if(MY_PROC(GH) == 0)cout << "Smoothing phi" << endl << flush ; 
      forall(phi,t+idt,l,c)
         f_smooth2(FBA(phi(t+idt,l,c)), FDA(phi(t+idt,l,c,MAIN)), &wght);
      end_forall
      if(MY_PROC(GH) == 0)cout << "Smoothing A" << endl << flush ; 
      forall(A,t+idt,l,c)
         f_smooth2(FBA(A(t+idt,l,c)), FDA(A(t+idt,l,c,MAIN)), &wght);
      end_forall
      if(MY_PROC(GH) == 0)cout << "Smoothing B" << endl << flush ; 
      forall(B,t+idt,l,c)
         f_smooth2(FBA(B(t+idt,l,c)), FDA(B(t+idt,l,c,MAIN)), &wght);
      end_forall
      if(MY_PROC(GH) == 0)cout << "Smoothing C" << endl << flush ; 
      forall(C,t+idt,l,c)
         f_smooth2(FBA(C(t+idt,l,c)), FDA(C(t+idt,l,c,MAIN)), &wght);
      end_forall

      Sync(A, t+idt, l, MAIN);
      Sync(B, t+idt, l, MAIN);
      Sync(C, t+idt, l, MAIN);

      Sync(phi, t+idt, l, MAIN);
      if(MY_PROC(GH) == 0)cout << "Done smoothing on MAIN" << endl << flush ; 
//      }
      

      // At this point t+idt, t, t-idt are defined. We do not move the
      // stencil up until we are down with children.

      // ******* Update Boundaries on MAIN ******

      if(MY_PROC(GH) == 0)cout << "Updating boundaries on MAIN" << endl << flush ; 
      UpdateExternalBoundary(phi,t+idt,l,MAIN);
      UpdateExternalBoundary(A,t+idt,l,MAIN);
      UpdateExternalBoundary(B,t+idt,l,MAIN);
      UpdateExternalBoundary(C,t+idt,l,MAIN);

      Sync(A, t+idt, l, MAIN);
      Sync(B, t+idt, l, MAIN);
      Sync(C, t+idt, l, MAIN);

      Sync(phi, t+idt, l, MAIN);

      // Calculate norms of A, B, C, phi on MAIN hierarchy
      DOUBLE anorm = Norm2(A, t+idt, l, MAIN);
      DOUBLE bnorm = Norm2(B, t+idt, l, MAIN);
      DOUBLE cnorm = Norm2(C, t+idt, l, MAIN);

      DOUBLE phinorm = Norm2(phi,t+idt, l, MAIN);

      if(me == 0){
         cout << "n+1: At level " << Level << endl << flush;
         cout << "MAIN norms " << anorm << " " << bnorm << " " << cnorm << " " 
              << phinorm << endl << flush;
      }
      // Calculate norms of A, B, C, phi on MAIN hierarchy
       anorm = Norm2(A, t-idt, l, MAIN);
       bnorm = Norm2(B, t-idt, l, MAIN);
       cnorm = Norm2(C, t-idt, l, MAIN);

       phinorm = Norm2(phi,t-idt, l, MAIN);

      if(me == 0){
         cout << "n-1: At level " << Level << endl << flush;
         cout << "MAIN norms " << anorm << " " << bnorm << " " << cnorm << " " 
              << phinorm << endl << flush;
      }


/************* SHADOW hierarchy update *************/

      if(StepsTaken(GH,Level) % RefineFactor(GH) == 0){

         DOUBLE dagh_h = DeltaX(GH,X_AXIS,l,SHADOW);
         DOUBLE dagh_dt = dagh_h * cfl;

         if(MY_PROC(GH) == 0 )cout << "Leapfrog on SHADOW " << endl << flush;

        // First update SHADOW from MAIN
         forall(A, t, l, c)
           A(t,l,c,SHADOW) = A(t,l,c,MAIN);
         end_forall
         forall(B, t, l, c)
           B(t,l,c,SHADOW) = B(t,l,c,MAIN);
         end_forall
         forall(C, t, l, c)
           C(t,l,c,SHADOW) = C(t,l,c,MAIN);
         end_forall

         forall(phi, t, l, c)
           phi(t,l,c,SHADOW) = phi(t,l,c,MAIN);
         end_forall
 
         Sync(A, t, l, SHADOW);
         Sync(B, t, l, SHADOW);
         Sync(C, t, l, SHADOW);

         Sync(phi, t, l, SHADOW);

         // Leapfrog on SHADOW

         forall(phi, t, l, c)
            f_leapfrog(FBA(A(t,l,c,SHADOW)),
              FDA(A(t+idt,l,c,SHADOW)), FDA(A(t-idt,l,c,SHADOW)), 
                FDA(A(t,l,c,SHADOW)),
              FDA(B(t+idt,l,c,SHADOW)), FDA(B(t-idt,l,c,SHADOW)), 
                FDA(B(t,l,c,SHADOW)),
              FDA(C(t+idt,l,c,SHADOW)), FDA(C(t-idt,l,c,SHADOW)), 
                FDA(C(t,l,c,SHADOW)),
              FDA(phi(t+idt,l,c,SHADOW)), FDA(phi(t-idt,l,c,SHADOW)),
              FDA(phi(t,l,c,SHADOW)), &dagh_dt, &dagh_h);
         end_forall

         Sync(A, t+idt, l, SHADOW);
         Sync(B, t+idt, l, SHADOW);
         Sync(C, t+idt, l, SHADOW);

         Sync(phi, t+idt, l, SHADOW);

         // Calculate norms of A, B, C, phi on SHADOW hierarchy
         DOUBLE anorm = Norm2(A, t+idt, l, SHADOW);
         DOUBLE bnorm = Norm2(B, t+idt, l, SHADOW);
         DOUBLE cnorm = Norm2(C, t+idt, l, SHADOW);
   
         DOUBLE phinorm = Norm2(phi,t+idt, l, SHADOW);

         if(me == 0){
            cout << "n+1: At level " << Level << endl << flush;
            cout << "SHADOW norms " << anorm << " " << bnorm << " " << cnorm << " " 
                 << phinorm << endl << flush;
         }
         // Calculate norms of A, B, C, phi on MAIN hierarchy
          anorm = Norm2(A, t-idt, l, SHADOW);
          bnorm = Norm2(B, t-idt, l, SHADOW);
          cnorm = Norm2(C, t-idt, l, SHADOW);
   
          phinorm = Norm2(phi,t-idt, l, SHADOW);
   
         if(me == 0){
            cout << "n-1: At level " << Level << endl << flush;
            cout << "SHADOW norms " << anorm << " " << bnorm << " " << cnorm << " " 
                 << phinorm << endl << flush;
         }

         // At this point t+idt, t, t-idt are defined. We do not move the
         // stencil up until we are down with children.

/**************We update SHADOW n,n-1 from MAIN at every SHADOW timestep***/

         /******* Update Boundaries on SHADOW ******/
//        if(MY_PROC(GH) == 0)cout << "Updating boundaries on SHADOW" << endl << flush ; 

//           UpdateExternalBoundary(phi,t+idt,l,SHADOW);
//           UpdateExternalBoundary(A,t+idt,l,SHADOW);
//           UpdateExternalBoundary(B,t+idt,l,SHADOW);
//           UpdateExternalBoundary(C,t+idt,l,SHADOW);

//         Sync(A, t+idt, l, SHADOW);
//         Sync(B, t+idt, l, SHADOW);
//         Sync(C, t+idt, l, SHADOW);

//         Sync(phi, t+idt, l, SHADOW);


      } // Ends if statement for SHADOW
    

    // If we are not at the finest level then go to next level
    // Go in with current level n+1 and n defined

    if (Level < FineLevel(GH)) {

      comm_service::log() << indent(Level) 
             << "********* Recursive Integrate Level+1 *********" << endl ;

       if(MY_PROC(GH) == 0 ){
            cout << "######################Integration of Level " 
                 << Level+1 << endl << flush ;
       }

       // Recursive Integration of next level in hierarchy

       // With current level Level n+1, n, n-1 defined go and update children

       bno_RecursiveIntegrate( GH, Level+1, A, B, C, phi, temp, dt, h, cfl);

       if(MY_PROC(GH) == 0 ){
            cout << "######################Done with Integration of Level " 
                 << Level+1 << endl << flush ;
       }

    } // End if (Level < FineLevel(GH))

    // Once children have been integrated move n to n-1 and n+1 to n

       l = Level;
       t = CurrentTime(GH,Level);
       idt = TimeStep(GH,Level);

    MoveStencilUp(MAIN, A, t, l,idt);
    MoveStencilUp(MAIN, B, t, l,idt);
    MoveStencilUp(MAIN, C, t, l,idt);
    MoveStencilUp(MAIN, phi, t, l,idt);

    if(StepsTaken(GH,Level) % RefineFactor(GH) == 0){
      MoveStencilUp(SHADOW, A, t, l,idt);
      MoveStencilUp(SHADOW, B, t, l,idt);
      MoveStencilUp(SHADOW, C, t, l,idt);
      MoveStencilUp(SHADOW, phi, t, l,idt);
    }

    // Increment time step on current level
    IncrCurrentTime(GH,Level);

    if (Level < FineLevel(GH)) {
      /* a sanity check */
      assert (CurrentTime(GH,Level) == CurrentTime(GH,Level+1)); 

      comm_service::log() << indent(Level) 
          << "********* Update Main From Children *********" << endl ;

      // Restriction on MAIN
      bno_UpdateLevelFromChildren(A,MAIN,GH,Level);
      bno_UpdateLevelFromChildren(B,MAIN,GH,Level);
      bno_UpdateLevelFromChildren(C,MAIN,GH,Level);
      bno_UpdateLevelFromChildren(phi,MAIN,GH,Level);

      // Restriction on SHADOW
//      if(StepsTaken(GH,Level) % RefineFactor(GH) == 0){
//        bno_UpdateLevelFromChildren(A,SHADOW,GH,Level);
//       bno_UpdateLevelFromChildren(B,SHADOW,GH,Level);
//      bno_UpdateLevelFromChildren(C,SHADOW,GH,Level);
//     bno_UpdateLevelFromChildren(phi,SHADOW,GH,Level);
//      }


    } // end if statement for level
    cout << "Steps taken = " << StepsTaken(GH,Level) << endl << flush;
  } // End For loop

  comm_service::log() << indent(Level) 
                << "********* ******************* *********" << endl << endl;



} // End bnoRecursiveIntegrate



void bno_UpdateLevelFromChildren(GridFunction(DIM)<GFTYPE>& u, 
        INTEGER const IDENT, GridHierarchy &GH, INTEGER const Level)
{
  INTEGER t = CurrentTime(GH,Level);

  Restrict(u, t, Level+1, Level, (GFTYPE *)0,0,IDENT);
  Sync(u, t, Level, IDENT);
}

void bno_RegridSystemAbove(GridHierarchy &GH, 
       GridFunction(DIM)<GFTYPE>& u, INTEGER const Level, GFTYPE const thresh)
{
    BBoxList bblist, tmpbblist;
    INTEGER flev = FineLevel(GH);

    INTEGER t = CurrentTime(GH,Level);
    INTEGER idt = TimeStep(GH,Level);
    INTEGER l = 0;

    DOUBLE newthresh;

    if (flev == 0) l = 0;
    else if (flev == MaxLev-1) l = flev-1;
    else l = flev;

    for (;l>=Level;l--) {

      t = CurrentTime(GH,l);
      idt = TimeStep(GH,l);

      int factor = 1 << (l+1);
//      newthresh = thresh / (1.0*factor) ;
//      newthresh = thresh / pow(2.0,(double)(factor));
//      if(MY_PROC(GH) == 0)cout << "Level " << l << "Threshold " << newthresh << endl << flush;
      Cluster2d(u, t, l,
              BLOCK_WIDTH, BUFFER_WIDTH,
              MIN_EFFICIENCY, thresh,
              bblist, bblist, MAIN);

      comm_service::log() << indent(Level) 
              << "********* From Clustering Routine *********" << endl << flush;
      comm_service::log() << indent(Level) << bblist << flush;

      if (1) {
         static int cnt = 0;
         char  fname[40];
         sprintf(fname,"GridSlice%d.dat",cnt);
         PlotBBox(MAIN, u, t, l, bblist, fname, VizServer);
         cnt++;
      }

      Refine(GH,bblist,l);

    } // Ends for loop over levels
    /*
        static mycnt2 = 0;
        if (1 && FineLevel(GH) != CoarseLevel(GH)) {
           char  fname[40];
           INTEGER idt = TimeStep(GH,FineLevel(GH));
           sprintf(fname,"SliceShadowFinePhiBefore%d.dat",mycnt2++);
           writeslice(SHADOW, fname, GH, u, t, FineLevel(GH), VizServer);
           sprintf(fname,"SliceFinePhiBefore%d.dat",mycnt2++);
           writeslice(MAIN, fname, GH, u, t, FineLevel(GH), VizServer);
        }
    */

    GH.DAGH_RecomposeHierarchy();

    /*
        static mycnt = 0;
        if (1 && FineLevel(GH) != CoarseLevel(GH)) {
           char  fname[40];
           INTEGER idt = TimeStep(GH,FineLevel(GH));
           sprintf(fname,"SliceShadowFinePhiAfter%d.dat",mycnt++);
           writeslice(SHADOW, fname, GH, u, t, FineLevel(GH), VizServer);
           sprintf(fname,"SliceFinePhiAfter%d.dat",mycnt++);
           writeslice(MAIN, fname, GH, u, t, FineLevel(GH), VizServer);
        }
    */

} // bno_RegridSystemAbove




void writeslice(int const IDENT,
                char const *fname,
                GridHierarchy &GH,
                GridFunction(DIM)<GFTYPE> &gf,
                INTEGER const t, INTEGER const l,
                INTEGER const viz_server)
  {
    FILE *dfile;
    INTEGER me = MY_PROC(GH);

    BBox bb;
    Coords s;
    forall(gf,t,l,c)
      bb += bbox(gf(t,l,c,IDENT));
      s = stepsize(gf(t,l,c,IDENT));
    end_forall
 
    bb.stepsize() = s;

    GridData(DIM)<GFTYPE> gd;
    gf.GF_gather_one(t,l,bb,bb,gd,viz_server,IDENT);
    if (me == viz_server) {
      if (!(dfile = fopen(fname, "w"))) exit(-1);

      for(INTEGER i=bb.lower(0);i<=bb.upper(0);i+=bb.stepsize(0)) {

        for(INTEGER j=bb.lower(1);j<=bb.upper(1);j+=bb.stepsize(1)) {
          fprintf(dfile,"%d\t%d\t%le\n",i,j,gd(i,j));
        } //end for over j

        fprintf(dfile,"\n");

      } //end for over i

      fclose(dfile);
     }
   }

void writeslicehdf(int const IDENT,
                char const *fname,
                GridHierarchy &GH,
                GridFunction(DIM)<GFTYPE> &gf,
                INTEGER const t, INTEGER const l,
                INTEGER const viz_server)
  {
    FILE *dfile;
    INTEGER me = MY_PROC(GH);
    INTEGER rank = DIM;
    INTEGER gft_rc;
 

    char *cnames[3];


    // Label coordinates with numbers 1 through DIM
    cnames[0] = (char *)malloc(sizeof(char)*3);
    cnames[1] = (char *)malloc(sizeof(char)*3);
    cnames[2] = (char *)malloc(sizeof(char)*3);

    sprintf(cnames[0],"1");
    sprintf(cnames[1],"2");
    sprintf(cnames[2],"3");

    

    BBox bb;
    Coords s;
    forall(gf,t,l,c)
      bb += bbox(gf(t,l,c,IDENT));
      s = stepsize(gf(t,l,c,IDENT));
    end_forall

    INTEGER *shape = bb.extents() ;

    DOUBLE *xy;

    xy = (DOUBLE *)malloc(sizeof(DOUBLE)*(shape[0]+shape[1]));
 
    bb.stepsize() = s;

    GridData(DIM)<GFTYPE> gd;
    gf.GF_gather_one(t,l,bb,bb,gd,viz_server,IDENT);
    if (me == viz_server) {
      if (!(dfile = fopen(fname, "w"))) exit(-1);

      // Set up coordinates
      for(INTEGER i=bb.lower(0);i<=bb.upper(0);i+=bb.stepsize(0)) {
          xy[i] = (DOUBLE)i;
      } // end for over i

      for(INTEGER j=bb.lower(1);j<=bb.upper(1);j+=bb.stepsize(1)) {
          xy[i + bb.upper(0)] = (DOUBLE)j;
      } //end for over j
      

    /****** DUMP HDF Data ******/
        gft_rc = gft_write_full(fname, t, bb.extents(),cnames,DIM, xy, gd.data());
        if(gft_rc == 0){
          cout << "Error in gft_write_full: Did not write to file " << fname 
               << endl << flush;
        }

     }
     free(cnames);
     free(xy);
   }


void SyncNext(INTEGER const IDENT, 
              GridFunction(DIM)<GFTYPE> &gf, 
              GridFunction(DIM)<GFTYPE> &temp, 
              INTEGER const t, INTEGER const l, INTEGER const idt)
  {
   Sync(gf,t+idt,l,IDENT);
  }

void SyncPrev(INTEGER const IDENT, 
              GridFunction(DIM)<GFTYPE> &gf, 
              GridFunction(DIM)<GFTYPE> &temp, 
              INTEGER const t, INTEGER const l, INTEGER const idt)
  {
   Sync(gf,t-idt,l,IDENT);
  }

void MoveStencilUp(INTEGER const IDENT,
                   GridFunction(DIM)<GFTYPE> &gf,
                   INTEGER const t, INTEGER const l, INTEGER const idt)
  {
   forall(gf, t, l, c)
      gf(t-idt,l,c,IDENT) = gf(t,l,c,IDENT);   // copy n to n-1
      gf(t,l,c,IDENT)  = gf(t+idt,l,c,IDENT);   // copy n+1 to n
   end_forall
  }

void PlotBBox(INTEGER const IDENT,
              GridFunction(DIM)<GFTYPE> &gf,
              INTEGER const t, INTEGER const l,
              BBoxList &bbl,
              char const *fname, INTEGER const viz_server)
  {
    BBox bb;
    Coords s;
    forall(gf,t,l,c)
      bb += bbox(gf(t,l,c,IDENT));
      s = stepsize(gf(t,l,c,IDENT));
    end_forall
 
    bb.stepsize() = s;
    Array(2)<INTEGER> flag_array(bb);
    flag_array = 0;

    BBox *_b = bbl.first();
    for (;_b;_b=bbl.next()) {
      flag_array.equals(1,*_b);
    }
             
    {
      FILE *dfile;
      INTEGER me = MY_PROC(GH);
      BBox viewbb = flag_array.bbox();
      if (me == viz_server) {
        if (!(dfile = fopen(fname, "w"))) exit(-1);

        for(INTEGER i=viewbb.lower(0);i<=viewbb.upper(0);i+=viewbb.stepsize(0))
        {
         for(INTEGER j=viewbb.lower(1);j<=viewbb.upper(1);j+=viewbb.stepsize(1))
          {
            fprintf(dfile,"%d\t%d\t%d\n",i,j,flag_array(i,j));
          } // End for over j
          fprintf(dfile,"\n");
        } // End for over i
        fclose(dfile);
      } 
    }
} //End PlotBBox