/*@@
   @file      bnamr3dEvolve.C
   @date      Thu Apr  1 11:42:03 1999
   @author    Tom Goodale
   @desc 
   The Berger-Oliger Evolution stuff for GrACE
   @enddesc 
 @@*/

#include <stdio.h>

#include "cctk.h"
#include "cctk_Flesh.h"
#include "cctk_Groups.h"
#include "cctk_IO.h"
#include "rfrConstants.h"

#include "GrACE_extension.h"
#include "GrACEDefs.h"

//#include "bnoamr3d/bnoamr3d.h"
#include "cctk_parameters.h"

static char *rcsid = "$Header: /mnt/data2/cvs2svn/cvs-repositories/projects/CactusGrACE/GrACE/src/bnoamr3dEvolve.C,v 1.10 1999-11-28 23:11:10 goodale Exp $";

void grace_FillcGH(cGH* cgh,
		   const int iteration,
		   const int l,
		   const int c,
		   const int ident);

void W3dAMR_bnoRecursiveIntegrate(cGH *cgh,
				  GridHierarchy &GH,
				  INTEGER const Level,
				  GridFunction(3)<DOUBLE> &temp,
				  DOUBLE &dt, DOUBLE &h,
				  INTEGER &currentiter);

extern "C" int CCTK_StepGH(cGH *);

#define cfl 0.9

int BnoAMR_Evolve(tFleshConfig* Configdata)
{
  DECLARE_CCTK_PARAMETERS;

  cGH* cgh = Configdata->GH[0];

  GrACEGH* ggh = (GrACEGH *)cgh->extensions[GrACE_GHExtension];
  GridHierarchy& GH = *ggh->theGH;

  GridFunction(3)<DOUBLE> &tempGF = *((GridFunction(3)<DOUBLE> *)ggh->tempGF);

  DOUBLE h  = (xmax - xmin) / (DOUBLE)(global_nx-1);
  DOUBLE dt = h * cfl;

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

  /* GAB: Removed itfirst (zero in Cactus), and use cctk_itlast)*/
  INTEGER NumIters = cctk_itlast + 1;
  INTEGER currentiter = 0;

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

    W3dAMR_bnoRecursiveIntegrate(cgh,
				 GH, l, 
				 tempGF, 
				 dt, h,
				 currentiter);

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

    cout << "Iteration " << currentiter << endl << flush;
   foreachGF(GH,GF,3,DOUBLE)
      cout << "\t" <<  GF.GF_Name() << " - " 
           << "[MaxVal: " << MaxVal(GF,t,l,DAGH_Main) << "] "
           << "[MinVal: " << MinVal(GF,t,l,DAGH_Main) << "] "
           << "[Norm2: " << Norm2(GF,t,l,DAGH_Main) << "] "
           << endl << flush;
   end_foreachGF
  } // End loop over number of iterations

  return 0;

} // 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 W3dAMR_bnoRecursiveIntegrate(cGH *cgh,
				  GridHierarchy &GH,
				  INTEGER const Level,
				  GridFunction(3)<DOUBLE> &temp,
				  DOUBLE &dt, DOUBLE &h,
				  INTEGER &currentiter)
{
  DECLARE_CCTK_PARAMETERS;

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

  DOUBLE foo;

  //printf("Level %d, dt %f, h %f\n", Level, dt, h);


//---------------------------------------------------------------  
// 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 l = Level;
    INTEGER t = CurrentTime(GH,Level,DAGH_Main);
    INTEGER idt = 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
      //      W3dAMR_Lte(GH,t,l,A, B, C, D, temp);

      // Regrid using truncation error estimate
      //      W3dAMR_bnoRegridSystemAbove(GH,temp,Level,Threshold);

    } // End if statement Level...

    /*********** Take a step on the DAGH_Main hierarchy ***********/   
    DOUBLE dagh_h = DeltaX(GH,DAGH_X,l,DAGH_Main);
    DOUBLE dagh_dt = dagh_h * cfl;

    // Step on DAGH_Main
    forall(temp, t, l, c) {
      grace_FillcGH(cgh, currentiter, l, c, DAGH_Main);
      CCTK_StepGH(cgh);
      // FIXME - put this here temporarily to get some output.
      CCTK_OutputGH(cgh);
    } end_forall;

    /************* DAGH_Shadow hierarchy update *************/
    if (StepsTaken(GH,Level,DAGH_Main) % RefineFactor(GH) == 0) {

       DOUBLE dagh_h = DeltaX(GH,DAGH_X,l,DAGH_Shadow);
       DOUBLE dagh_dt = dagh_h * cfl;

       // First update DAGH_Shadow from DAGH_Main
       foreachGF(GH,GF,DIM,GFTYPE) {
	 forall(GF, t, l, c) {
	   GF(t,l,c,DAGH_Shadow) = GF(t,l,c,DAGH_Main);
	 } end_forall;
       } end_foreachGF;
       Sync(GH, t, l, DAGH_Shadow);       

       // Step on DAGH_Shadow
       forall(temp, t, l, c) {
	 grace_FillcGH(cgh, currentiter, l, c, DAGH_Shadow);
	 CCTK_StepGH(cgh);
	 // FIXME - put this here temporarily to get some output.
	 CCTK_OutputGH(cgh);
       } end_forall;
    } // Ends if statement for DAGH_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)) {
      // Recursive Integration of next level in hierarchy
      // With current level Level n+1, n, n-1 defined go and update children
      W3dAMR_bnoRecursiveIntegrate(cgh, GH, Level+1, temp, 
				   dt, h, currentiter);

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

    // Once children have been integrated move n to n-1 and n+1 to n
    foreachGF(GH,GF,DIM,GFTYPE) {
      CycleTimeLevels(GF,Level);
    } end_foreachGF;
  
    // Increment time step on current level
    IncrCurrentTime(GH,Level,DAGH_Main);
    
    if (Level < FineLevel(GH)) {  
      // Restriction on DAGH_Main
      foreachGF(GH,GF,DIM,GFTYPE) {
	W3dAMR_bnoUpdateLevelFromChildren(GF,DAGH_Main,GH,Level);
      } end_foreachGF;
    } // end if statement for level
  } // End For loop
} // End bnoRecursiveIntegrate

void W3dAMR_bnoUpdateLevelFromChildren(GridFunction(DIM)<GFTYPE>& u, 
				       INTEGER const IDENT, 
				       GridHierarchy &GH, 
				       INTEGER const Level)
{
  INTEGER t = CurrentTime(GH,Level);
  
  INTEGER myargc = 0; GFTYPE myargs[1];
  Restrict(u, t, Level+1, t, Level, myargs, myargc, IDENT);
  Sync(u, t, Level, IDENT);
}