1 files changed, 155 insertions, 0 deletions
diff --git a/src/WaveToy.cc b/src/WaveToy.cc
new file mode 100644
index 0000000..b9621e7
--- /dev/null
+++ b/src/WaveToy.cc
@@ -0,0 +1,155 @@
+ /*@@
+   @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"
+
+
+ /*@@
+   @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; 
+
+	ApplySymmetry(cctkGH,"wavetoycxx::scalartmps"); 
+
+	if (CCTK_EQUALS(bound,"flat")) 
+	{
+		ierr = ApplyFlatBC(cctkGH,sw,"wavetoycxx::phi_new");
+	}
+	else if (CCTK_Equals(bound,"radiation")) 
+	{
+		ierr = ApplyRadiativeBC(cctkGH,0,1,sw,"wavetoycxx::phi_new","wavetoy::phi");
+	} 
+
+	if (ierr < 0) 
+	{
+		CCTK_WARN(0,"Boundary conditions not applied - giving up!");
+	}
+}
+
+#define val(gridfunc,i,j,k)  gridfunc[CCTK_GFINDEX3D(cctkGH,i,j,k)]
+
+
+ /*@@
+   @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];
+	}
+   }	
+
+}
+
+
+