Correct errors.

Add extrapolator to set up the boundary conditions Gamma for initial
data.


git-svn-id: http://svn.einsteintoolkit.org/cactus/EinsteinEvolve/NewRad/trunk@3 16f5cafb-89ad-4b9f-9d0b-9d7a0a37d03b

4 files changed, 296 insertions, 66 deletions

diff --git a/interface.ccl b/interface.ccl
index 57f6416..234262d 100644
--- a/interface.ccl
+++ b/interface.ccl
@@ -5,12 +5,17 @@ IMPLEMENTS: NewRad
 USES INCLUDE HEADER: GenericFD.h
 
 CCTK_INT FUNCTION                         \
+    ExtrapolateGammas                     \
+        (CCTK_POINTER_TO_CONST IN cctkGH, \
+         CCTK_REAL ARRAY INOUT var)
+PROVIDES FUNCTION ExtrapolateGammas WITH ExtrapolateGammas1 LANGUAGE C
+
+CCTK_INT FUNCTION                         \
     NewRad_Apply                          \
         (CCTK_POINTER_TO_CONST IN cctkGH, \
          CCTK_REAL ARRAY IN var,          \
          CCTK_REAL ARRAY INOUT rhs,       \
          CCTK_REAL IN var0,               \
          CCTK_REAL IN v0,                 \
-         CCTK_REAL IN radpower,           \
-         CCTK_INT IN width)
+         CCTK_INT IN radpower)
 PROVIDES FUNCTION NewRad_Apply WITH NewRad_Apply1 LANGUAGE C
diff --git a/src/extrap.cc b/src/extrap.cc
new file mode 100644
index 0000000..b752aa0
--- /dev/null
+++ b/src/extrap.cc
@@ -0,0 +1,186 @@
+#include <cassert>
+#include <cmath>
+
+#include <cctk.h>
+
+#ifdef CCTK_CXX_RESTRICT
+#  undef restrict
+#  define restrict CCTK_CXX_RESTRICT
+#endif
+
+#define KRANC_C
+extern "C" {
+#include <GenericFD.h>
+}
+
+using namespace std;
+
+
+
+// Adapted from BSSN_MoL's files Init.F
+static
+void extrap_kernel (cGH const* restrict const cctkGH,
+                    int const* restrict const bmin,
+                    int const* restrict const bmax,
+                    int const* restrict const dir,
+                    CCTK_REAL* restrict const var)
+{
+  int const ni = cctkGH->cctk_lsh[0];
+  int const nj = cctkGH->cctk_lsh[1];
+  int const nk = cctkGH->cctk_lsh[2];
+  
+  int const si = -dir[0];
+  int const sj = -dir[1];
+  int const sk = -dir[2];
+  
+  int const di = 1;
+  int const dj = ni;
+  int const dk = ni*nj;
+  int const np = ni*nj*nk;
+  
+  int const dind = si*di + sj*dj + sk*dk;
+  
+  int imin[3], imax[3], idir[3];
+  for (int d=0; d<3; ++d) {
+    if (dir[d]<0) {
+      // lower boundary
+      imin[d] = bmax[d]-1;
+      imax[d] = bmin[d]-1;
+      idir[d] = -1;
+    } else {
+      // interior and upper boundary
+      imin[d] = bmin[d];
+      imax[d] = bmax[d];
+      idir[d] = +1;
+    }
+  }
+  
+  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]) {
+        int const ind = CCTK_GFINDEX3D(cctkGH, i,j,k);
+        
+        // Test looping directions
+        if (i==0) assert (idir[0]<0);
+        if (j==0) assert (idir[1]<0);
+        if (k==0) assert (idir[2]<0);
+        if (i==ni-1) assert (idir[0]>0);
+        if (j==nj-1) assert (idir[1]>0);
+        if (k==nk-1) assert (idir[2]>0);
+        
+        // Apply boundary conditions to get e.g. Gammas on physical
+        // boundaries.  Notice that we only want to apply boundary
+        // conditions to the Gammas, all other variables have been
+        // defined point-wise all the way to the boundaries from the
+        // original ADM quantities (which SHOULD be correctly defined
+        // all the way to the boundaries by the initial data
+        // routines).  Here I use cubic extrapolation (though rational
+        // extrapolation might be better).
+        
+        assert (ind       >=0 and ind       <np);
+        assert (ind+4*dind>=0 and ind+4*dind<np);
+        var[ind] = (4*var[ind+dind] - 6*var[ind+2*dind] + 4*var[ind+3*dind]
+                    - var[ind+4*dind]);
+        
+      } // for i j k
+    }
+  }
+}
+
+
+
+// Adapted from Kranc's KrancNumericalTools/GenericFD's file
+// GenericFD.c
+static
+void extrap_loop (cGH const* restrict const cctkGH,
+                  CCTK_REAL* restrict const var)
+{
+  int imin[3], imax[3], is_symbnd[6], is_physbnd[6], is_ipbnd[6];
+  GenericFD_GetBoundaryInfo
+    (cctkGH, cctkGH->cctk_lsh, cctkGH->cctk_lssh, cctkGH->cctk_bbox,
+     cctkGH->cctk_nghostzones, 
+     imin, imax, is_symbnd, is_physbnd, is_ipbnd);
+  
+  // Loop over all faces:
+  // Loop over faces first, then corners, athen edges, so that the
+  // stencil only sees points that have already been treated
+  // ifec means: interior-face-edge-corner
+  for (int ifec=1; ifec<=3; ++ifec) {
+    for (int dir2=-1; dir2<=+1; ++dir2) {
+      for (int dir1=-1; dir1<=+1; ++dir1) {
+        for (int dir0=-1; dir0<=+1; ++dir0) {
+          int const dir[3] = { dir0, dir1, dir2 };
+          
+          int nnz = 0;
+          for (int d=0; d<3; ++d) {
+            if (dir[d]) ++nnz;
+          }
+          if (nnz == ifec) {
+            
+            // one of tahe faces is a boundary
+            bool have_bnd = false;
+            // all boundary faces are physical boundaries
+            bool all_physbnd = true;
+            
+            int bmin[3], bmax[3];
+            for (int d=0; d<3; ++d) {
+              switch (dir[d]) {
+              case -1:
+                bmin[d] = 0;
+                bmax[d] = imin[d];
+                have_bnd = true;
+                all_physbnd = all_physbnd and is_physbnd[2*d+0];
+                break;
+              case 0:
+                bmin[d] = imin[d];
+                bmax[d] = imax[d];
+                break;
+              case +1:
+                bmin[d] = imax[d];
+                bmax[d] = cctkGH->cctk_lssh[CCTK_LSSH_IDX(0,d)];
+                have_bnd = true;
+                all_physbnd = all_physbnd and is_physbnd[2*d+1];
+                break;
+              }
+            }
+            
+            if (have_bnd and all_physbnd) {
+              extrap_kernel (cctkGH, bmin, bmax, dir, var);
+            }
+            
+          }
+        } // for dir0 dir1 dir2
+      }
+    }
+  }
+}
+
+
+
+extern "C"
+CCTK_INT ExtrapolateGammas1 (CCTK_POINTER_TO_CONST const cctkGH_,
+                             CCTK_REAL* restrict const var)
+{
+  cGH const* restrict const cctkGH = static_cast<cGH const*> (cctkGH_);
+  if (not cctkGH) {
+    CCTK_WARN (CCTK_WARN_ABORT,
+               "cctkGH is NULL");
+  }
+  
+#if 0
+  CCTK_REAL* restrict const var = CCTK_VarDataPtr (cctkGH, 0, varname);
+  if (not var) {
+    CCTK_VWarn (CCTK_WARN_ABORT, __LINE__, __FILE__, CCTK_THORNSTRING,
+                "Cannot access variable \"%s\"", varname);
+  }
+#endif
+  
+  if (not var) {
+    CCTK_WARN (CCTK_WARN_ABORT,
+               "Pointer to variable is NULL");
+  }
+  
+  extrap_loop (cctkGH, var);
+  
+  return 0;
+}
diff --git a/src/make.code.defn b/src/make.code.defn
index 244b721..40d1780 100644
--- a/src/make.code.defn
+++ b/src/make.code.defn
@@ -1,7 +1,7 @@
 # Main make.code.defn file for thorn NewRad
 
 # Source files in this directory
-SRCS = newrad.cc
+SRCS = extrap.cc newrad.cc
 
 # Subdirectories containing source files
 SUBDIRS = 
diff --git a/src/newrad.cc b/src/newrad.cc
index a1edcf0..3845075 100644
--- a/src/newrad.cc
+++ b/src/newrad.cc
@@ -1,4 +1,5 @@
-#include <math.h>
+#include <cassert>
+#include <cmath>
 
 #include <cctk.h>
 
@@ -12,6 +13,8 @@ extern "C" {
 #include <GenericFD.h>
 }
 
+using namespace std;
+
 
 
 // Adapted from BSSN_MoL's files NewRad.F and newrad.h
@@ -28,15 +31,15 @@ void newrad_kernel (cGH const* restrict const cctkGH,
                     CCTK_REAL const* restrict const r,
                     CCTK_REAL const& var0,
                     CCTK_REAL const& v0,
-                    CCTK_REAL const& radpower)
+                    int const radpower)
 {
   int const ni = cctkGH->cctk_lsh[0];
   int const nj = cctkGH->cctk_lsh[1];
   int const nk = cctkGH->cctk_lsh[2];
   
-  int const si = -dir[0];
-  int const sj = -dir[1];
-  int const sk = -dir[2];
+  int const si = dir[0];
+  int const sj = dir[1];
+  int const sk = dir[2];
   
   int const di = 1;
   int const dj = ni;
@@ -49,11 +52,34 @@ void newrad_kernel (cGH const* restrict const cctkGH,
   CCTK_REAL const idy = 1.0/dy;
   CCTK_REAL const idz = 1.0/dz;
   
-  for (int k=bmin[2]; k<bmax[2]; ++k) {
-    for (int j=bmin[1]; j<bmax[1]; ++j) {
-      for (int i=bmin[0]; i<bmax[0]; ++i) {
+  int imin[3], imax[3], idir[3];
+  for (int d=0; d<3; ++d) {
+    if (dir[d]<0) {
+      // lower boundary
+      imin[d] = bmax[d]-1;
+      imax[d] = bmin[d]-1;
+      idir[d] = -1;
+    } else {
+      // interior and upper boundary
+      imin[d] = bmin[d];
+      imax[d] = bmax[d];
+      idir[d] = +1;
+    }
+  }
+  
+  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]) {
         int const ind = CCTK_GFINDEX3D(cctkGH, i,j,k);
         
+        // Test looping directions
+        if (i==0) assert (idir[0]<0);
+        if (j==0) assert (idir[1]<0);
+        if (k==0) assert (idir[2]<0);
+        if (i==ni-1) assert (idir[0]>0);
+        if (j==nj-1) assert (idir[1]>0);
+        if (k==nk-1) assert (idir[2]>0);
+        
         {
           // The main part of the boundary condition assumes that we
           // have an outgoing radial wave with some speed v0:
@@ -74,13 +100,13 @@ 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 = i==0 ? +1 : i==ni-1 ? -1 : vx>0 ? +1 : -1;
-          int const svj = j==0 ? +1 : j==nj-1 ? -1 : vy>0 ? +1 : -1;
-          int const svk = k==0 ? +1 : k==nk-1 ? -1 : vz>0 ? +1 : -1;
+          int const svi = dir[0];
+          int const svj = dir[1];
+          int const svk = dir[2];
           
           // Find x derivative
           CCTK_REAL derivx;
-          if (i>0 and i<ni-1) {
+          if (svi==0) {
             derivx = 0.5*(var[ind+di]-var[ind-di])*idx;
           } else {
             derivx = svi*0.5*(3*var[ind] - 4*var[ind-svi*di]
@@ -89,7 +115,7 @@ void newrad_kernel (cGH const* restrict const cctkGH,
           
           // Find y derivative
           CCTK_REAL derivy;
-          if (j>0 and j<nj-1) {
+          if (svj==0) {
             derivy = 0.5*(var[ind+dj]-var[ind-dj])*idy;
           } else {
             derivy = svj*0.5*(3*var[ind] - 4*var[ind-svj*dj]
@@ -98,7 +124,7 @@ void newrad_kernel (cGH const* restrict const cctkGH,
           
           // Find z derivative
           CCTK_REAL derivz;
-          if (k>0 and k<nk-1) {
+          if (svk==0) {
             derivz = 0.5*(var[ind+dk]-var[ind-dk])*idz;
           } else {
             derivz = svk*0.5*(3*var[ind] - 4*var[ind-svk*dk]
@@ -130,6 +156,9 @@ void newrad_kernel (cGH const* restrict const cctkGH,
           int const ip = i-si;
           int const jp = j-sj;
           int const kp = k-sk;
+          assert (ip>=0 and ip<ni);
+          assert (jp>=0 and jp<nj);
+          assert (kp>=0 and kp<nk);
           
           // Find local wave speeds
           int const indp = CCTK_GFINDEX3D(cctkGH, ip,jp,kp);
@@ -141,13 +170,13 @@ 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 = i==0 ? +1 : i==ni-1 ? -1 : vx>0 ? +1 : -1;
-          int const svj = j==0 ? +1 : j==nj-1 ? -1 : vy>0 ? +1 : -1;
-          int const svk = k==0 ? +1 : k==nk-1 ? -1 : vz>0 ? +1 : -1;
+          int const svi = dir[0];
+          int const svj = dir[1];
+          int const svk = dir[2];
           
           // Find x derivative
           CCTK_REAL derivx;
-          if (ip>0 and ip<ni-1) {
+          if (svi==0) {
             derivx = 0.5*(var[indp+di]-var[indp-di])*idx;
           } else {
             derivx = svi*0.5*(3*var[indp] - 4*var[indp-svi*di]
@@ -156,7 +185,7 @@ void newrad_kernel (cGH const* restrict const cctkGH,
           
           // Find y derivative
           CCTK_REAL derivy;
-          if (jp>0 and jp<nj-1) {
+          if (svj==0) {
             derivy = 0.5*(var[indp+dj]-var[indp-dj])*idy;
           } else {
             derivy = svj*0.5*(3*var[indp] - 4*var[indp-svj*dj]
@@ -165,7 +194,7 @@ void newrad_kernel (cGH const* restrict const cctkGH,
           
           // Find z derivative
           CCTK_REAL derivz;
-          if (kp>0 and kp<nk-1) {
+          if (svk==0) {
             derivz = 0.5*(var[indp+dk]-var[indp-dk])*idz;
           } else {
             derivz = svk*0.5*(3*var[indp] - 4*var[indp-svk*dk]
@@ -201,8 +230,7 @@ void newrad_loop (cGH const* restrict const cctkGH,
                   CCTK_REAL const* restrict const r,
                   CCTK_REAL const& var0,
                   CCTK_REAL const& v0,
-                  CCTK_REAL const& radpower,
-                  int const width)
+                  int const radpower)
 {
   int imin[3], imax[3], is_symbnd[6], is_physbnd[6], is_ipbnd[6];
   GenericFD_GetBoundaryInfo
@@ -210,45 +238,57 @@ void newrad_loop (cGH const* restrict const cctkGH,
      cctkGH->cctk_nghostzones, 
      imin, imax, is_symbnd, is_physbnd, is_ipbnd);
   
-  // Loop over all faces
-  for (int dir2=-1; dir2<=+1; ++dir2) {
-    for (int dir1=-1; dir1<=+1; ++dir1) {
-      for (int dir0=-1; dir0<=+1; ++dir0) {
-        int const dir[3] = { dir0, dir1, dir2 };
-        
-        // one of tahe faces is a boundary
-        bool have_bnd = false;
-        // all boundary faces are physical boundaries
-        bool all_physbnd = true;
-        
-        int bmin[3], bmax[3];
-        for (int d=0; d<3; ++d) {
-          switch (dir[d]) {
-          case -1:
-            bmin[d] = 0;
-            bmax[d] = imin[d];
-            have_bnd = true;
-            all_physbnd = all_physbnd and is_physbnd[2*d+0];
-            break;
-          case 0:
-            bmin[d] = imin[d];
-            bmax[d] = imax[d];
-            break;
-          case +1:
-            bmin[d] = imax[d];
-            bmax[d] = cctkGH->cctk_lssh[CCTK_LSSH_IDX(0,d)];
-            have_bnd = true;
-            all_physbnd = all_physbnd and is_physbnd[2*d+1];
-            break;
+  // Loop over all faces:
+  // Loop over faces first, then corners, athen edges, so that the
+  // stencil only sees points that have already been treated
+  // ifec means: interior-face-edge-corner
+  for (int ifec=1; ifec<=3; ++ifec) {
+    for (int dir2=-1; dir2<=+1; ++dir2) {
+      for (int dir1=-1; dir1<=+1; ++dir1) {
+        for (int dir0=-1; dir0<=+1; ++dir0) {
+          int const dir[3] = { dir0, dir1, dir2 };
+          
+          int nnz = 0;
+          for (int d=0; d<3; ++d) {
+            if (dir[d]) ++nnz;
           }
-        }
-        
-        if (have_bnd and all_physbnd) {
-          newrad_kernel (cctkGH, bmin, bmax, dir,
-                         var, rhs, x,y,z,r, var0, v0, radpower);
-        }
-        
-      } // for dir0 dir1 dir2
+          if (nnz == ifec) {
+            
+            // one of tahe faces is a boundary
+            bool have_bnd = false;
+            // all boundary faces are physical boundaries
+            bool all_physbnd = true;
+            
+            int bmin[3], bmax[3];
+            for (int d=0; d<3; ++d) {
+              switch (dir[d]) {
+              case -1:
+                bmin[d] = 0;
+                bmax[d] = imin[d];
+                have_bnd = true;
+                all_physbnd = all_physbnd and is_physbnd[2*d+0];
+                break;
+              case 0:
+                bmin[d] = imin[d];
+                bmax[d] = imax[d];
+                break;
+              case +1:
+                bmin[d] = imax[d];
+                bmax[d] = cctkGH->cctk_lssh[CCTK_LSSH_IDX(0,d)];
+                have_bnd = true;
+                all_physbnd = all_physbnd and is_physbnd[2*d+1];
+                break;
+              }
+            }
+            
+            if (have_bnd and all_physbnd) {
+              newrad_kernel (cctkGH, bmin, bmax, dir,
+                             var, rhs, x,y,z,r, var0, v0, radpower);
+            }
+            
+          }
+        } // for dir0 dir1 dir2
+      }
     }
   }
 }
@@ -261,8 +301,7 @@ CCTK_INT NewRad_Apply1 (CCTK_POINTER_TO_CONST const cctkGH_,
                         CCTK_REAL      * restrict const rhs,
                         CCTK_REAL const var0,
                         CCTK_REAL const v0,
-                        CCTK_REAL const radpower,
-                        CCTK_INT const width)
+                        CCTK_INT const radpower)
 {
   cGH const* restrict const cctkGH = static_cast<cGH const*> (cctkGH_);
   if (not cctkGH) {
@@ -305,7 +344,7 @@ CCTK_INT NewRad_Apply1 (CCTK_POINTER_TO_CONST const cctkGH_,
                "Cannot access coordinate variables x, y, z, and r");
   }
   
-  newrad_loop (cctkGH, var, rhs, x,y,z,r, var0, v0, radpower, width);
+  newrad_loop (cctkGH, var, rhs, x,y,z,r, var0, v0, radpower);
   
   return 0;
 }