/*@@
   @file      WaveToy.c
   @date      
   @author    Tom Goodale
   @desc 
              Evolution routines for the wave equation solver
   @enddesc 
 @@*/
 
#include "cctk_WarnLevel.h" 
#include "cctk.h" 
#include "cctk_parameters.h"
#include "cctk_arguments.h"

#include "cctk_Flesh.h"
#include "cctk_Misc.h"
#include "cctk_Comm.h"
#include "cctk_Groups.h"

#include "CactusBase/Boundary/src/Boundary.h"
#include "CactusBase/CartGrid3D/src/Symmetry.h"

static char *rcsid = "$Header$";

#define val(gridfunc,i,j,k)  gridfunc[CCTK_GFINDEX3D(cctkGH,i,j,k)]

static void WaveToyC_Boundaries(CCTK_CARGUMENTS);

 /*@@
   @routine    WaveToyC_Evolution
   @date       
   @author     Tom Goodale
   @desc 
               Evolution for the wave equation
   @enddesc 
   @calls      CCTK_SyncGroup, WaveToyC_Boundaries
   @calledby   
   @history 
 
   @endhistory 

@@*/

extern "C" void  WaveToyCXX_Evolution(CCTK_CARGUMENTS)
{

  DECLARE_CCTK_CARGUMENTS

  // Set up shorthands 

  CCTK_REAL dx = CCTK_DELTA_SPACE(0),
            dy = CCTK_DELTA_SPACE(1),
            dz = CCTK_DELTA_SPACE(2),
            dt = CCTK_DELTA_TIME,
      
            dx2=dx*dx,
	    dy2=dy*dy,
	    dz2=dz*dz,
	    dt2=dt*dt,

            factor = 2*(1 - (dt2)*(1/dx2 + 1/dy2 + 1/dz2));

  int   iend = cctk_lsh[0]-1,
        jend = cctk_lsh[1]-1,
        kend = cctk_lsh[2]-1; 

  //
  // Do the evolution 
  // 

  for (int k=1; k<kend; k++)
    for (int j=1; j<jend; j++)
      for (int i=1; i<iend; i++)
      {
        val(phi_new, i,j,k) =
		     factor* val( phi,i,j,k ) - val( phi_old,i,j,k) 
              + dt2 * 
                   ( ( val( phi, i+1,j  ,k  ) + val( phi, i-1,j  ,k)   )/dx2
                    +( val( phi, i  ,j+1,k  ) + val( phi, i  ,j-1,k)   )/dy2
                    +( val( phi, i  ,j  ,k+1) + val( phi, i  ,j,  k-1) )/dz2
		   );
      } 


  // 
  // Synchronize before applying boundary conditions 
  // 
  CCTK_SyncGroup(cctkGH,"wavetoycxx::scalartmps"); 

  //
  // Apply boundary conditions 
  //
  WaveToyC_Boundaries(CCTK_PASS_CTOC); 

  //
  // Update timeslices 
  //
  {	
    for (int k=0; k<cctk_lsh[2]; k++)
      for (int j=0; j<cctk_lsh[1]; j++)
	for (int i=0; i<cctk_lsh[0]; i++)
	{
          int     index = CCTK_GFINDEX3D(cctkGH,i,j,k); 
                  phi_old[index]  = phi[index]; 
		  phi[index]      = phi_new[index];
	}
   }	

  USE_CCTK_CARGUMENTS

}

 /*@@
   @routine    WaveToyC_Boundaries
   @date       
   @author     Tom Goodale
   @desc 
               Boundary conditions for the wave equation
   @enddesc 
   @calls      ApplySymmetry,ApplyFlatBC,ApplyRadiativeBC 
   @calledby   
   @history 
 
   @endhistory 

@@*/

static void WaveToyC_Boundaries(CCTK_CARGUMENTS)
{
  DECLARE_CCTK_CARGUMENTS
  DECLARE_CCTK_PARAMETERS

  int ierr;
  int sw[3]; 

  /* Set the stencil width */ 
  sw[0]=1; 
  sw[1]=1; 
  sw[2]=1; 

  ierr = CartSymBCGroup(cctkGH,"wavetoycxx::scalartmps");

  if (CCTK_Equals(bound,"flat")) 
  {
    ierr = FlatBCVar(cctkGH,sw,"wavetoycxx::phi_new");
  }
  else if (CCTK_Equals(bound,"radiation")) 
  {
    ierr = RadiativeBCVar(cctkGH,0,1,sw,"wavetoycxx::phi_new","wavetoy::phi");
  } 

  if (ierr < 0) 
  {
    CCTK_WARN(0,"Boundary conditions not applied - giving up!");
  }

  USE_CCTK_PARAMETERS
  USE_CCTK_CARGUMENTS

}