Schedule boundary conditions at postrestrict.

git-svn-id: http://svn.cactuscode.org/arrangements/CactusWave/WaveToyC/trunk@93 90eea020-d82d-4da5-bf6e-4ee79ff7632f

4 files changed, 213 insertions, 29 deletions

diff --git a/interface.ccl b/interface.ccl
index d963b48..1463592 100644
--- a/interface.ccl
+++ b/interface.ccl
@@ -14,7 +14,7 @@ cctk_real scalarevolve type = GF Timelevels=3 tags='tensortypealias="Scalar"'
 } "The evolved scalar field"
 
 
-CCTK_INT FUNCTION Boundary_SelectVarForBC(CCTK_POINTER IN GH, \
+CCTK_INT FUNCTION Boundary_SelectVarForBC(CCTK_POINTER_TO_CONST IN GH, \
   CCTK_INT IN faces, CCTK_INT IN boundary_width, CCTK_INT IN table_handle, \
   CCTK_STRING IN var_name, CCTK_STRING IN bc_name)
 USES FUNCTION Boundary_SelectVarForBC
diff --git a/param.ccl b/param.ccl
index a382ecb..87ea382 100644
--- a/param.ccl
+++ b/param.ccl
@@ -14,4 +14,22 @@ KEYWORD bound "Type of boundary condition to use"
 } "none"
 
 
- 
\ No newline at end of file
+
+private:
+
+INT repetitions "Number of timing repetions"
+{
+  0 : * :: "1 is the normal case"
+} 1
+
+INT level1_stride [3] "Stride for level 1 loops"
+{
+  -1    :: "off"
+  1 : * :: "on"
+} -1
+
+INT level2_stride [3] "Stride for level 2 loops"
+{
+  -1    :: "off"
+  1 : * :: "on: should be an integer multiple of level1_stride"
+} -1
diff --git a/schedule.ccl b/schedule.ccl
index bb18e8f..ae86396 100644
--- a/schedule.ccl
+++ b/schedule.ccl
@@ -11,6 +11,7 @@ schedule WaveToyC_Startup at STARTUP
 schedule WaveToyC_InitSymBound at BASEGRID
 {
   LANG: C
+  OPTIONS: global
 } "Schedule symmetries"
 
 schedule WaveToyC_Evolution as WaveToy_Evolution at EVOL
@@ -22,8 +23,19 @@ schedule WaveToyC_Evolution as WaveToy_Evolution at EVOL
 schedule WaveToyC_Boundaries as WaveToy_Boundaries at EVOL AFTER WaveToy_Evolution
 {
   LANG: C
+  # TODO: split this into local and level mode routines
 } "Boundaries of 3D wave equation"
   
 schedule GROUP ApplyBCs as WaveToyC_ApplyBCs at EVOL after WaveToy_Boundaries
 {
 } "Apply boundary conditions"
+
+schedule WaveToyC_Boundaries as WaveToy_Boundaries at POSTRESTRICT
+{
+  LANG: C
+  # TODO: split this into local and level mode routines
+} "Boundaries of 3D wave equation"
+
+schedule GROUP ApplyBCs as WaveToyC_ApplyBCs at POSTRESTRICT after WaveToy_Boundaries
+{
+} "Apply boundary conditions"
diff --git a/src/WaveToy.c b/src/WaveToy.c
index 7c98b1e..5ecbdd8 100644
--- a/src/WaveToy.c
+++ b/src/WaveToy.c
@@ -8,8 +8,15 @@
    @version $Header$
  @@*/
 
-  
-#include "cctk.h"
+#include <assert.h>
+
+#if 0
+#include <sys/time.h>
+#include <sys/resource.h>
+#include <unistd.h>
+#endif
+
+#include "cctk.h" 
 #include "cctk_Arguments.h"
 #include "cctk_Parameters.h"
 
@@ -20,7 +27,44 @@ static const char *rcsid = "$Header$";
 CCTK_FILEVERSION(CactusWave_WaveToyC_WaveToy_c);
 
 void WaveToyC_Boundaries(CCTK_ARGUMENTS);
-void  WaveToyC_Evolution(CCTK_ARGUMENTS);
+void WaveToyC_Evolution(CCTK_ARGUMENTS);
+
+#define DEBUG 0
+
+
+
+static inline int min (int const a, int const b)
+{
+  return a<b ? a : b;
+}
+
+static inline int max (int const a, int const b)
+{
+  return a>b ? a : b;
+}
+    
+#if 0
+static inline long long timeval2longlong (const struct timeval * const tv)
+{
+  return tv->tv_sec * 1000000LL + tv->tv_usec;
+}
+
+static long long currenttime (void)
+{
+  struct rusage ru;
+  getrusage (RUSAGE_SELF, &ru);
+  return timeval2longlong (ru.ru_utime) + timeval2longlong (ru.ru_stime);
+}
+#endif
+
+static long long currenttime (void)
+{
+  unsigned long eax, edx;
+  asm volatile ("rdtsc" : "=a" (eax), "=d" (edx));
+  return (unsigned long long)edx << 32 | (unsigned long long)eax;
+}
+
+
 
  /*@@
    @routine    WaveToyC_Evolution
@@ -37,14 +81,14 @@ void  WaveToyC_Evolution(CCTK_ARGUMENTS);
 
 void  WaveToyC_Evolution(CCTK_ARGUMENTS)
 {
-
   DECLARE_CCTK_ARGUMENTS;
+  DECLARE_CCTK_PARAMETERS;
       
   int i,j,k;
   int index;
   int istart, jstart, kstart, iend, jend, kend;
   CCTK_REAL dx,dy,dz,dt,dx2,dy2,dz2,dt2;
-  CCTK_REAL dx2i,dy2i,dz2i;
+  CCTK_REAL dt2dx2i,dt2dy2i,dt2dz2i;
 
   CCTK_REAL factor;
 
@@ -59,9 +103,9 @@ void  WaveToyC_Evolution(CCTK_ARGUMENTS)
   dz2 = dz*dz;
   dt2 = dt*dt;
 
-  dx2i = 1.0/dx2;
-  dy2i = 1.0/dy2;
-  dz2i = 1.0/dz2;
+  dt2dx2i = dt2/dx2;
+  dt2dy2i = dt2/dy2;
+  dt2dz2i = dt2/dz2;
 
   istart = 1;
   jstart = 1;
@@ -72,29 +116,139 @@ void  WaveToyC_Evolution(CCTK_ARGUMENTS)
   kend = cctk_lsh[2]-1;
 
   /* Do the evolution */
-  factor = 2*(1 - (dt2)*(dx2i + dy2i + dz2i));
+  factor = 2 * (1 - (dt2dx2i + dt2dy2i + dt2dz2i));
 
-  for (k=kstart; k<kend; k++)
   {
-    for (j=jstart; j<jend; j++)
+    /* user-defined values */
+    const char * const loopthorn = CCTK_THORNSTRING;
+    const char * const loopfile = __FILE__;
+    int          const loopline = __LINE__;
+    const char * const loopname = "Evolution";
+    
+    int const dim = 3;
+    
+    int const imin = istart;
+    int const jmin = jstart;
+    int const kmin = kstart;
+    int const imax = iend;
+    int const jmax = jend;
+    int const kmax = kend;
+    int const ishape = cctk_lsh[0];
+    int const jshape = cctk_lsh[1];
+    int const kshape = cctk_lsh[2];
+    
+    /* tuning parameters */
+    if (DEBUG) printf ("level1_stride = [%d,%d,%d]\n", level1_stride[0], level1_stride[1], level1_stride[2]);
+    if (DEBUG) printf ("level2_stride = [%d,%d,%d]\n", level2_stride[0], level2_stride[1], level2_stride[2]);
+    
+    int const i1step = level1_stride[0]<0 ? imax-imin : level1_stride[0];
+    int const j1step = level1_stride[1]<0 ? jmax-jmin : level1_stride[1];
+    int const k1step = level1_stride[2]<0 ? kmax-kmin : level1_stride[2];
+    
+    int const i2step = level2_stride[0]<0 ? imax-imin : level2_stride[0];
+    int const j2step = level2_stride[1]<0 ? jmax-jmin : level2_stride[1];
+    int const k2step = level2_stride[2]<0 ? kmax-kmin : level2_stride[2];
+    
+    /* internal stuff */
+    int rep;
+    for (rep = 0; rep < repetitions; ++rep) {
+    
+    long long l0time  = 0, l1time  = 0, l2time  = 0;
+    long long l0count = 0, l1count = 0, l2count = 0;
+    
+    /* level 2 loop */
+    l2time = currenttime() - l2time;
+    ++ l2count;
+    
+    int const i2min = imin;
+    int const j2min = jmin;
+    int const k2min = kmin;
+    int const i2max = imax;
+    int const j2max = jmax;
+    int const k2max = kmax;
+    
+    int i2, j2, k2;
+    
+    for (k2 = k2min; k2 < k2max; k2 += k2step) {
+    for (j2 = j2min; j2 < j2max; j2 += j2step) {
+    for (i2 = i2min; i2 < i2max; i2 += i2step) {
+    
+    /* level 1 loop */
+    l1time = currenttime() - l1time;
+    ++ l1count;
+    
+    int const i1min = i2;
+    int const j1min = j2;
+    int const k1min = k2;
+    int const i1max = min (imax, i2 + i2step);
+    int const j1max = min (jmax, j2 + j2step);
+    int const k1max = min (kmax, k2 + k2step);
+    
+    int i1, j1, k1;
+    
+    for (k1 = k1min; k1 < k1max; k1 += k1step) {
+    for (j1 = j1min; j1 < j1max; j1 += j1step) {
+    for (i1 = i1min; i1 < i1max; i1 += i1step) {
+    
+    /* level 0 loop */
+    l0time = currenttime() - l0time;
+    ++ l0count;
+    
+    int const i0min = i1;
+    int const j0min = j1;
+    int const k0min = k1;
+    int const i0max = min (imax, i1 + i1step);
+    int const j0max = min (jmax, j1 + j1step);
+    int const k0max = min (kmax, k1 + k1step);
+    
+    int i0, j0, k0;
+    
+    for (k0 = k0min; k0 < k0max; ++k0) {
+    for (j0 = j0min; j0 < j0max; ++j0) {
+    for (i0 = i0min; i0 < i0max; ++i0) {
+    
+    int const i = i0;
+    int const j = j0;
+    int const k = k0;
+    
+    int const di = 1;
+    int const dj = ishape;
+    int const dk = ishape * jshape;
+    
+    int const index = i *di + j * dj + k * dk;
+    
     {
-      for (i=istart; i<iend; i++)
-      {
-        index = CCTK_GFINDEX3D(cctkGH,i,j,k);
-        
-        phi[index] = factor*
-	             phi_p[index] - phi_p_p[index] 
-                     + (dt2) * 
-                   ( ( phi_p[CCTK_GFINDEX3D(cctkGH,i+1,j  ,k  )]
-                      +phi_p[CCTK_GFINDEX3D(cctkGH,i-1,j  ,k  )] )*dx2i
-                    +( phi_p[CCTK_GFINDEX3D(cctkGH,i  ,j+1,k  )]
-                      +phi_p[CCTK_GFINDEX3D(cctkGH,i  ,j-1,k  )] )*dy2i
-                    +( phi_p[CCTK_GFINDEX3D(cctkGH,i  ,j  ,k+1)]
-                      +phi_p[CCTK_GFINDEX3D(cctkGH,i  ,j,  k-1)] )*dz2i);
-      } 
-    }
+        phi[index] =   factor * phi_p[index] - phi_p_p[index] 
+                     + (phi_p[index+di] + phi_p[index-di]) * dt2dx2i
+                     + (phi_p[index+dj] + phi_p[index-dj]) * dt2dy2i
+                     + (phi_p[index+dk] + phi_p[index-dk]) * dt2dz2i;
+    } 
+    
+    } } } /* end level 0 loop */
+    l0time = currenttime() - l0time;
+    
+    } } } /* end level 1 loop */
+    l1time = currenttime() - l1time;
+    
+    } } } /* end level 2 loop */
+    l2time = currenttime() - l2time;
+    
+    if (DEBUG)
+      CCTK_VInfo (CCTK_THORNSTRING,
+                  "Timing information for loop %s::%s\n"
+                  "   (file \"%s\", line %d):\n"
+                  "   Level 0: shape %d %d %d, count %lld, time %lld\n"
+                  "   Level 1: shape %d %d %d, count %lld, time %lld\n"
+                  "   Level 2: shape %d %d %d, count %lld, time %lld",
+                  loopthorn, loopname, loopfile, loopline,
+                  i1step   , j1step   , k1step   , l0count, l0time,
+                  i2step   , j2step   , k2step   , l1count, l1time,
+                  imax-imin, jmax-jmin, kmax-kmin, l2count, l2time);
+    
+    } /* end repetition loop */
+    
   }
-
+    
   return;
 }