From fe99659bbdc43965546e53cd7c01adc9183d43aa Mon Sep 17 00:00:00 2001
From: schnetter <schnetter@16f5cafb-89ad-4b9f-9d0b-9d7a0a37d03b>
Date: Mon, 25 Jan 2010 16:26:09 +0000
Subject: Add assert statements with consistency checks.

Rename some variables for clarity.


git-svn-id: http://svn.einsteintoolkit.org/cactus/EinsteinEvolve/NewRad/trunk@5 16f5cafb-89ad-4b9f-9d0b-9d7a0a37d03b
---
 src/newrad.cc | 61 ++++++++++++++++++++++++++++++++---------------------------
 1 file changed, 33 insertions(+), 28 deletions(-)

diff --git a/src/newrad.cc b/src/newrad.cc
index 3845075..f87e78b 100644
--- a/src/newrad.cc
+++ b/src/newrad.cc
@@ -56,17 +56,30 @@ void newrad_kernel (cGH const* restrict const cctkGH,
   for (int d=0; d<3; ++d) {
     if (dir[d]<0) {
       // lower boundary
+      assert (bmin[d] >= 0);
+      assert (bmax[d] + 2 <= cctkGH->cctk_lsh[d]);
       imin[d] = bmax[d]-1;
       imax[d] = bmin[d]-1;
       idir[d] = -1;
+    } else if (dir[d]>0) {
+      // upper boundary
+      assert (bmin[d] - 2 >= 0);
+      assert (bmax[d] <= cctkGH->cctk_lsh[d]);
+      imin[d] = bmin[d];
+      imax[d] = bmax[d];
+      idir[d] = +1;
     } else {
-      // interior and upper boundary
+      // interior
+      assert (bmin[d] - 1 >= 0);
+      assert (bmax[d] + 1 <= cctkGH->cctk_lsh[d]);
       imin[d] = bmin[d];
       imax[d] = bmax[d];
       idir[d] = +1;
     }
   }
   
+  // Warning: these loops are not parallel, since previously
+  // calculated RHS are accessed for radpower>=0
   for (int k=imin[2]; k!=imax[2]; k+=idir[2]) {
     for (int j=imin[1]; j!=imax[1]; j+=idir[1]) {
       for (int i=imin[0]; i!=imax[0]; i+=idir[0]) {
@@ -100,35 +113,31 @@ void newrad_kernel (cGH const* restrict const cctkGH,
           CCTK_REAL const vy = v0*y[ind]*rpi;
           CCTK_REAL const vz = v0*z[ind]*rpi;
           
-          int const svi = dir[0];
-          int const svj = dir[1];
-          int const svk = dir[2];
-          
           // Find x derivative
           CCTK_REAL derivx;
-          if (svi==0) {
+          if (si==0) {
             derivx = 0.5*(var[ind+di]-var[ind-di])*idx;
           } else {
-            derivx = svi*0.5*(3*var[ind] - 4*var[ind-svi*di]
-                              + var[ind-2*svi*di])*idx;
+            derivx = si*0.5*(3*var[ind] - 4*var[ind-si*di]
+                             + var[ind-2*si*di])*idx;
           }
           
           // Find y derivative
           CCTK_REAL derivy;
-          if (svj==0) {
+          if (sj==0) {
             derivy = 0.5*(var[ind+dj]-var[ind-dj])*idy;
           } else {
-            derivy = svj*0.5*(3*var[ind] - 4*var[ind-svj*dj]
-                              + var[ind-2*svj*dj])*idy;
+            derivy = sj*0.5*(3*var[ind] - 4*var[ind-sj*dj]
+                             + var[ind-2*sj*dj])*idy;
           }
           
           // Find z derivative
           CCTK_REAL derivz;
-          if (svk==0) {
+          if (sk==0) {
             derivz = 0.5*(var[ind+dk]-var[ind-dk])*idz;
           } else {
-            derivz = svk*0.5*(3*var[ind] - 4*var[ind-svk*dk]
-                              + var[ind-2*svk*dk])*idz;
+            derivz = sk*0.5*(3*var[ind] - 4*var[ind-sk*dk]
+                             + var[ind-2*sk*dk])*idz;
           }
           
           // Calculate source term
@@ -170,35 +179,31 @@ void newrad_kernel (cGH const* restrict const cctkGH,
           CCTK_REAL const vy = v0*y[indp]*rpi;
           CCTK_REAL const vz = v0*z[indp]*rpi;
           
-          int const svi = dir[0];
-          int const svj = dir[1];
-          int const svk = dir[2];
-          
           // Find x derivative
           CCTK_REAL derivx;
-          if (svi==0) {
+          if (si==0) {
             derivx = 0.5*(var[indp+di]-var[indp-di])*idx;
           } else {
-            derivx = svi*0.5*(3*var[indp] - 4*var[indp-svi*di]
-                              + var[indp-2*svi*di])*idx;
+            derivx = si*0.5*(3*var[indp] - 4*var[indp-si*di]
+                             + var[indp-2*si*di])*idx;
           }
           
           // Find y derivative
           CCTK_REAL derivy;
-          if (svj==0) {
+          if (sj==0) {
             derivy = 0.5*(var[indp+dj]-var[indp-dj])*idy;
           } else {
-            derivy = svj*0.5*(3*var[indp] - 4*var[indp-svj*dj]
-                              + var[indp-2*svj*dj])*idy;
+            derivy = sj*0.5*(3*var[indp] - 4*var[indp-sj*dj]
+                             + var[indp-2*sj*dj])*idy;
           }
           
           // Find z derivative
           CCTK_REAL derivz;
-          if (svk==0) {
+          if (sk==0) {
             derivz = 0.5*(var[indp+dk]-var[indp-dk])*idz;
           } else {
-            derivz = svk*0.5*(3*var[indp] - 4*var[indp-svk*dk]
-                              + var[indp-2*svk*dk])*idz;
+            derivz = sk*0.5*(3*var[indp] - 4*var[indp-sk*dk]
+                             + var[indp-2*sk*dk])*idz;
           }
           
           // Find difference in sources
@@ -254,7 +259,7 @@ void newrad_loop (cGH const* restrict const cctkGH,
           }
           if (nnz == ifec) {
             
-            // one of tahe faces is a boundary
+            // one of the faces is a boundary
             bool have_bnd = false;
             // all boundary faces are physical boundaries
             bool all_physbnd = true;
-- 
cgit v1.2.3