//=============================================================================
// application interface functions for wave example
//=============================================================================
// PAMR example to solve 2D wave equation
// See http://laplace.physics.ubc.ca/People/arman/file/wave2d_pamr.tar.gz for source

#include <stdlib.h>
#include <stdio.h>
#include <pamr.h>
#include <amrd.h>
#include <math.h>
#include <mpi.h>
#include "wave.h"

#include <sys/timeb.h>
void elapsed_time(void);

//=============================================================================
// id parameters (time symmetric gaussian)
//=============================================================================

//Initialization parameters:
real amp, tdam, xc, xwid, yc, ywid;
//=====================================

//=============================================================================
// some convenient, "local" global variables
//=============================================================================

//The pointers to actual grid function data:
real *f_t_n, *f_t_np1, *f_n, *f_np1;
//=======================================

//The pointers to grid:
real *x, *y;
//===========================

//Parameters that describe the grid:
int shape[2]; // number of grid points in each dim, basically shape[0]=Nx, shape[1]=Ny
int ghost_width[4]; //side of ghost cells, 2*dim array, at x_min, x_max, y_min and so on
int Nx,Ny;
int phys_bdy[4]; // 1 or 0, wether or not boundary is physical, has size = dim*2
int lb_flag,rb_flag,tb_flag,bb_flag; //convenient variable for phys boundary indication
int size; //size=Nx*Ny*...
int dim; //the spatial dimension of the grid (only 1,2 or 3)
int g_rank; // To store the MPI rank of the execution

real base_bbox[4], bbox[4], dx, dy, dt;
// bbox stores the coordinate bounding box [x1,x2,y1...]
// base_bbox is the coordinate of the entire domain
int g_L; //to store the grid level (in AMR heirarchy)


//To store the grid function numbers, see set_gfns subroutine
int f_t_n_gfn, f_t_np1_gfn, f_n_gfn, f_np1_gfn;
//===========


//=============================================================================
// call after variables have been defined
//=============================================================================
void set_gfns(void)
{
      if ((f_t_np1_gfn = PAMR_get_gfn("f_t",PAMR_AMRH,1))<0) AMRD_stop("set_gnfs error",0);
      if ((f_t_n_gfn   = PAMR_get_gfn("f_t",PAMR_AMRH,2))<0) AMRD_stop("set_gnfs error",0);

      if ((f_np1_gfn   = PAMR_get_gfn("f",PAMR_AMRH,1))<0) AMRD_stop("set_gnfs error",0);
      if ((f_n_gfn     = PAMR_get_gfn("f",PAMR_AMRH,2))<0) AMRD_stop("set_gnfs error",0);

// get_gfn returns the grid function number:
// PAMR_get_gfn(A,B,C)
// A >> grid function name
// B >> in which heirarchy, two options: PAMR_AMRH: AMR heirarchy, PAMR_MGH: multigrid h.
// C >> time level if it is in AMRH, else 0, latest time ordered first.
}

//=============================================================================
// call with valid iter to set up globals:
//=============================================================================
void ldptr(void)
{
     real dx0[2];
     real *x0[2]; //array of size dim, pointing to local perimeter coordinate array
     real *gfs[PAMR_MAX_GFNS]; //array of size ngfs, pointers to local data array
   static int first=1;

   if (first)
   {
      first=0;
      set_gfns();
      PAMR_get_global_bbox(base_bbox); //AAK??
   }
//These are to get the information of the grid
   PAMR_get_g_dim(&dim); //saves the dimension of the grid in dim
   PAMR_get_g_rank(&g_rank); // gets the MPI rank
   PAMR_get_g_shape(shape); //gets the shape (Nx,Ny...)
   PAMR_get_g_bbox(bbox); //gets the boundary (x1,x2,y1,y2...)
   PAMR_get_g_ghost_width(ghost_width); //gets the ghots width at boundaries
   PAMR_get_g_level(&g_L); // gets the level number of the current grid
   PAMR_get_dxdt(g_L,dx0,&dt); //clear
   dx=dx0[0]; //convenient variables??
   dy=dx0[1];

//AAK: This is for finding the real boundary
   if ((bbox[0]-base_bbox[0])<dx/2) phys_bdy[0]=1; else phys_bdy[0]=0;
   if ((base_bbox[1]-bbox[1])<dx/2) phys_bdy[1]=1; else phys_bdy[1]=0;
   Nx=shape[0];

   if ((bbox[2]-base_bbox[2])<dy/2) phys_bdy[2]=1; else phys_bdy[2]=0;
      if ((base_bbox[3]-bbox[3])<dy/2) phys_bdy[3]=1; else phys_bdy[3]=0;
      Ny=shape[1];
      size=Nx*Ny;

   //setting the convenient variables
   lb_flag = phys_bdy[0];
   rb_flag = phys_bdy[1];
   bb_flag = phys_bdy[2];
   tb_flag = phys_bdy[3];
//=========================
   PAMR_get_g_x(x0); //gets the pointer to the coordinate array

   x=x0[0]; //the pointers to the array of actual grid data (size Nx*Ny..)
   y=x0[1];

   PAMR_get_g_gfs(gfs); //gets the pointer to local data of grid function

     //AAK:store the pointer in convenient variables
     f_n     = gfs[f_n_gfn-1]; //AAK?? what is -1, numbering starts from  1?
     f_np1   = gfs[f_np1_gfn-1];
     f_t_n   = gfs[f_t_n_gfn-1];
     f_t_np1 = gfs[f_t_np1_gfn-1];

}

//=============================================================================
// utility routines
//=============================================================================
void const_f(real *f, real c)
{
   int i;

     for (i=0; i<Nx*Ny; i++) f[i]=c;
}

void zero(real *f)
{
   const_f(f,0);
}

//=============================================================================
// Routines required by amrd:
//=============================================================================

//=============================================================================
// Returns 0 to use default mechanism, or is expected to calculate
// the correct initial hierarchy and return 1:
//=============================================================================
int wave_id(void)
{
   if( my_rank == 0 ) elapsed_time();
   return 0;
}

//=============================================================================
// Sets custom parameters, variables, etc. Split up into two segments,
// one called before the pamr context is initialized and standard
// parameters are read, and the other afterwards
//=============================================================================
void wave_var_pre_init(char *pfile)
{
   return;
}


//AAK: one of the hook functions passed to amrd, to read the parameters from file
void wave_var_post_init(char *pfile)
{
   int i,j;
   char buf[64]; //AAK??

   if (my_rank==0)
   {
      system("date > date.dat");
      printf("===================================================================\n");
      printf("Reading wave parameters:\n\n");
   }

  
     amp=tdam=xc=yc=0;
     xwid=ywid=1;
     //AAK: the file is in RNPL format
     // last argument is the size of the array, 1 for just a real number
   AMRD_real_param(pfile,"amp",&amp,1);
   AMRD_real_param(pfile,"tdam",&tdam,1);
   AMRD_real_param(pfile,"xc",&xc,1);
   AMRD_real_param(pfile,"yc",&yc,1);
   AMRD_real_param(pfile,"xwid",&xwid,1);
   AMRD_real_param(pfile,"ywid",&ywid,1);
  

   if (my_rank==0) printf("===================================================================\n");
   return;
}

//=============================================================================
// Sets all t=n variables to their 'zero' values:
//=============================================================================
void wave_AMRH_var_clear(void)
{
    ldptr();

//   zero(phi_n);

     zero(f_n);
     zero(f_t_n);
   return;
}

//=============================================================================
// Initial data for free fields: (at tn=2)
//
// currently, we only allow for time-symmetric initial data
//=============================================================================
void wave_free_data(void)
{
   ldptr();
//This initializes the grid, all the calculations are encapsulated in the routine.
     initializer0_(f_t_n,f_n,&Nx,&Ny,x,y,&amp,&tdam,&xc,&xwid,&yc,&ywid);


   return;
}  

//=============================================================================
// Initial constraint data --- called after each MG iteration.
//=============================================================================
void wave_t0_cnst_data(void)
{
   return;
}

//=============================================================================
// Calculations prior to saving info to disk.
//
// NOTE: at this point, the time sequence is: n,nm1,np1
//=============================================================================
void wave_pre_io_calc(void)
{
   return;
}

//=============================================================================
// Returns some norm of the residual for the evolution variables,
// called after an evolution iteration.
// We're using an explicit scheme to solve for phi, hence return 0
//=============================================================================
real wave_evo_residual(void)
{
   return 0;
}

//=============================================================================
// Returns some norm of the residual for the MG variables, *AND*
// stores the point-wise residual in "f_res" for each MG variable "f" (for
// computing new RHS's)
//=============================================================================
real wave_MG_residual(void)
{
   return 0;
}

//=============================================================================
// Performs 1 iteration of the evolution equations
//=============================================================================
void wave_evolve(int iter)
{
   ldptr();
//This updates the grid, all the calculations are encapsulated in the routine
#ifdef __L
    update0__(f_t_np1,f_t_n,f_np1,f_n,&Nx,&Ny,&dt,&dx,&dy,&bb_flag,&lb_flag,&rb_flag,&tb_flag);
#else
   update0_(f_t_np1,f_t_n,f_np1,f_n,&Nx,&Ny,&dt,&dx,&dy,&bb_flag,&lb_flag,&rb_flag,&tb_flag);
#endif

   return;
}

//=============================================================================
// sets excision mask (NO ITERATOR, SO DON'T LOAD POINTERS!!!)
//=============================================================================
void wave_fill_ex_mask(real *mask, int dim, int *shape, real *bbox, real excised)
{
}

//=============================================================================
void wave_fill_bh_bboxes(real *bbox, int *num, int max_num)
{
}

//=============================================================================
void wave_post_tstep(int L)
{
   return;
}

//=============================================================================
// Performs 1 relaxation sweep of the MG equations, and returns an estimate
// of the norm of the residual.
//=============================================================================
real wave_MG_relax(void)
{
   return 0;
}

//=============================================================================
// Computes the differential operator L acting on the current grid,
// in a region specified by the grid function "mask". Stores the result
// in "f_lop" for each MG variable "f"
//=============================================================================
void wave_L_op(void)
{
   return;
}

//=============================================================================
// Called after calculating the TRE for all variables
//=============================================================================
void wave_scale_tre(void)
{
   return;
}

//=============================================================================
// post-regrid initialization of constant functions
//=============================================================================
void wave_post_regrid(void)
{
}

//=============================================================================
int main(int argc, char **argv)
{
   amrd(argc,argv,&wave_id,&wave_var_pre_init,
        &wave_var_post_init, &wave_AMRH_var_clear,
        &wave_free_data, &wave_t0_cnst_data,
        &wave_evo_residual, &wave_MG_residual,
        &wave_evolve, &wave_MG_relax, &wave_L_op,
        &wave_pre_io_calc, &wave_scale_tre,
        &wave_post_regrid, &wave_post_tstep,
        &wave_fill_ex_mask, &wave_fill_bh_bboxes);
   if (my_rank==0) elapsed_time();
}

//=============================================================================
// Maintains/reports elapsed wall-clock time.
//=============================================================================
void elapsed_time(void) {
   static int    first = 1;
   struct        timeb t;
   static double msinit;
   double        mscurr, mselapsed;

   ftime(&t);
   mscurr = 1000.0 * t.time + t.millitm;
   if( first ) {
      msinit = mscurr;
      first = 0;
   }
   mselapsed = mscurr - msinit;
   printf("elapsed_time: Seconds since initial call: %12.3f\n",
         mselapsed / 1000.0);
}