Improve some comments.

Add more self-tests.

Apply extrapolation and radiative boundaries even in the ghost zones
of the outer boundary.  This means that it is not necessary any more
to synchronise afterwards.


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

2 files changed, 87 insertions, 9 deletions

diff --git a/src/extrap.cc b/src/extrap.cc
index b752aa0..b44d0ae 100644
--- a/src/extrap.cc
+++ b/src/extrap.cc
@@ -18,6 +18,8 @@ using namespace std;
 
 
 // Adapted from BSSN_MoL's files Init.F
+// Erik Schnetter: I assume this code was probably originally written
+// by Miguel Alcubierre.
 static
 void extrap_kernel (cGH const* restrict const cctkGH,
                     int const* restrict const bmin,
@@ -55,6 +57,8 @@ void extrap_kernel (cGH const* restrict const cctkGH,
     }
   }
   
+  // Note: These loops are not parallel, since each iteration may
+  // access grid points set by previous iterations.
   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]) {
@@ -102,9 +106,9 @@ void extrap_loop (cGH const* restrict const cctkGH,
      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
+  // Loop over faces first, then corners, and then 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) {
@@ -117,10 +121,15 @@ void extrap_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;
+            // at least one boundary face is a physical boundary
+            bool any_physbnd = false;
+            // all boundary faces are either physical boundaries or
+            // ghost zones
+            bool all_physbnd_or_ghostbnd = true;
             
             int bmin[3], bmax[3];
             for (int d=0; d<3; ++d) {
@@ -130,6 +139,9 @@ void extrap_loop (cGH const* restrict const cctkGH,
                 bmax[d] = imin[d];
                 have_bnd = true;
                 all_physbnd = all_physbnd and is_physbnd[2*d+0];
+                any_physbnd = any_physbnd or is_physbnd[2*d+0];
+                all_physbnd_or_ghostbnd = all_physbnd_or_ghostbnd and
+                  (is_physbnd[2*d+0] or is_ipbnd[2*d+0]);
                 break;
               case 0:
                 bmin[d] = imin[d];
@@ -140,11 +152,15 @@ void extrap_loop (cGH const* restrict const cctkGH,
                 bmax[d] = cctkGH->cctk_lssh[CCTK_LSSH_IDX(0,d)];
                 have_bnd = true;
                 all_physbnd = all_physbnd and is_physbnd[2*d+1];
+                any_physbnd = any_physbnd or is_physbnd[2*d+1];
+                all_physbnd_or_ghostbnd = all_physbnd_or_ghostbnd and
+                  (is_physbnd[2*d+1] or not is_ipbnd[2*d+1]);
                 break;
               }
             }
+            assert (have_bnd);  // must be true since nnz>0
             
-            if (have_bnd and all_physbnd) {
+            if (have_bnd and any_physbnd and all_physbnd_or_ghostbnd) {
               extrap_kernel (cctkGH, bmin, bmax, dir, var);
             }
             
diff --git a/src/newrad.cc b/src/newrad.cc
index f87e78b..c39c360 100644
--- a/src/newrad.cc
+++ b/src/newrad.cc
@@ -18,6 +18,8 @@ using namespace std;
 
 
 // Adapted from BSSN_MoL's files NewRad.F and newrad.h
+// Erik Schnetter: I assume this code was probably originally written
+// by Miguel Alcubierre.
 static
 void newrad_kernel (cGH const* restrict const cctkGH,
                     int const* restrict const bmin,
@@ -93,6 +95,22 @@ void newrad_kernel (cGH const* restrict const cctkGH,
         if (j==nj-1) assert (idir[1]>0);
         if (k==nk-1) assert (idir[2]>0);
         
+        if (si==0) {
+          assert (i-1>=0 and i+1<ni);
+        } else {
+          assert (i-2*si>=0 and i-2*si<ni);
+        }
+        if (sj==0) {
+          assert (j-1>=0 and j+1<nj);
+        } else {
+          assert (j-2*sj>=0 and j-2*sj<nj);
+        }
+        if (sk==0) {
+          assert (k-1>=0 and k+1<nk);
+        } else {
+          assert (k-2*sk>=0 and k-2*sk<nk);
+        }
+        
         {
           // The main part of the boundary condition assumes that we
           // have an outgoing radial wave with some speed v0:
@@ -105,6 +123,22 @@ void newrad_kernel (cGH const* restrict const cctkGH,
           //
           // where  vi = v0 xi/r
           
+          if (si==0) {
+            assert (i-1>=0 and i+1<ni);
+          } else {
+            assert (i-2*si>=0 and i-2*si<ni);
+          }
+          if (sj==0) {
+            assert (j-1>=0 and j+1<nj);
+          } else {
+            assert (j-2*sj>=0 and j-2*sj<nj);
+          }
+          if (sk==0) {
+            assert (k-1>=0 and k+1<nk);
+          } else {
+            assert (k-2*sk>=0 and k-2*sk<nk);
+          }
+          
           // Find local wave speeds
           CCTK_REAL const rp = r[ind];
           CCTK_REAL const rpi = 1.0/rp;
@@ -169,6 +203,22 @@ void newrad_kernel (cGH const* restrict const cctkGH,
           assert (jp>=0 and jp<nj);
           assert (kp>=0 and kp<nk);
           
+          if (si==0) {
+            assert (ip-1>=0 and ip+1<ni);
+          } else {
+            assert (ip-2*si>=0 and ip-2*si<ni);
+          }
+          if (sj==0) {
+            assert (jp-1>=0 and jp+1<nj);
+          } else {
+            assert (jp-2*sj>=0 and jp-2*sj<nj);
+          }
+          if (sk==0) {
+            assert (kp-1>=0 and kp+1<nk);
+          } else {
+            assert (kp-2*sk>=0 and kp-2*sk<nk);
+          }
+          
           // Find local wave speeds
           int const indp = CCTK_GFINDEX3D(cctkGH, ip,jp,kp);
           
@@ -244,9 +294,9 @@ void newrad_loop (cGH const* restrict const cctkGH,
      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
+  // Loop over faces first, then corners, and then 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) {
@@ -263,6 +313,11 @@ void newrad_loop (cGH const* restrict const cctkGH,
             bool have_bnd = false;
             // all boundary faces are physical boundaries
             bool all_physbnd = true;
+            // at least one boundary face is a physical boundary
+            bool any_physbnd = false;
+            // all boundary faces are either physical boundaries or
+            // ghost zones
+            bool all_physbnd_or_ghostbnd = true;
             
             int bmin[3], bmax[3];
             for (int d=0; d<3; ++d) {
@@ -272,6 +327,9 @@ void newrad_loop (cGH const* restrict const cctkGH,
                 bmax[d] = imin[d];
                 have_bnd = true;
                 all_physbnd = all_physbnd and is_physbnd[2*d+0];
+                any_physbnd = any_physbnd or is_physbnd[2*d+0];
+                all_physbnd_or_ghostbnd = all_physbnd_or_ghostbnd and
+                  (is_physbnd[2*d+0] or not is_ipbnd[2*d+0]);
                 break;
               case 0:
                 bmin[d] = imin[d];
@@ -282,11 +340,15 @@ void newrad_loop (cGH const* restrict const cctkGH,
                 bmax[d] = cctkGH->cctk_lssh[CCTK_LSSH_IDX(0,d)];
                 have_bnd = true;
                 all_physbnd = all_physbnd and is_physbnd[2*d+1];
+                any_physbnd = any_physbnd or is_physbnd[2*d+1];
+                all_physbnd_or_ghostbnd = all_physbnd_or_ghostbnd and
+                  (is_physbnd[2*d+1] or not is_ipbnd[2*d+1]);
                 break;
               }
             }
+            assert (have_bnd);  // must be true since nnz>0
             
-            if (have_bnd and all_physbnd) {
+            if (have_bnd and any_physbnd and all_physbnd_or_ghostbnd) {
               newrad_kernel (cctkGH, bmin, bmax, dir,
                              var, rhs, x,y,z,r, var0, v0, radpower);
             }