Adding directional boundary conditions for everything apart from Robin.

See thorn documentation for details

Next to come, boundary conditions for grid functions which aren't CCTK_REALs.


git-svn-id: http://svn.cactuscode.org/arrangements/CactusBase/Boundary/trunk@131 6a38eb6e-646e-4a02-a296-d141613ad6c4

1 files changed, 517 insertions, 185 deletions

diff --git a/src/RadiationBoundary.c b/src/RadiationBoundary.c
index dcc853f..68ab210 100644
--- a/src/RadiationBoundary.c
+++ b/src/RadiationBoundary.c
@@ -8,7 +8,7 @@
    @version $Header$
  @@*/
 
-/*#define DEBUG_BOUND*/
+/*#define DEBUG_BOUNDARY*/
  
 #include <stdio.h>
 #include <assert.h>
@@ -34,10 +34,9 @@ CCTK_FILEVERSION(CactusBase_Boundary_RadiationBoundaryWrappers_c)
 /* Internal routine prototypes */
 
 static int BndApplyRadiative3Di(cGH *GH, 
-				int  gdim,
-				int *sw,  
-				int *doBC, 
+				int doBC, 
 				int *lssh,
+				int *stencil_size,  
 				CCTK_REAL *dxyz, 
 				CCTK_REAL dt, 
 				CCTK_REAL *var_n, 
@@ -50,12 +49,14 @@ static int BndApplyRadiative3Di(cGH *GH,
 				CCTK_REAL v0) ;
 
 static int ApplyBndRadiative(cGH *GH, 
-			     int *sw,
+			     int stencil,
+			     int *stencil_array,
 			     CCTK_REAL var0, 
 			     CCTK_REAL v0, 
 			     int var_current, 
 			     int var_previous,
-			     int num_vars);
+			     int num_vars,
+			     int dir);
 
 /* Local variables */
 
@@ -65,6 +66,46 @@ static int ApplyBndRadiative(cGH *GH,
  ********************************************************************/
 
 /*@@
+   @routine    BndRadiativeDirGI
+   @date       Sat Jan 20 2001
+   @author     Gabrielle Allen
+   @desc 
+      Aply radiative BC's by group index in given direction
+   @enddesc 
+   @calls     
+   @calledby   
+   @history 
+   @endhistory 
+
+@@*/
+
+int BndRadiativeDirGI(cGH *GH, 
+		   int stencil_size,
+		   int dir,
+		   CCTK_REAL var0, 
+		   CCTK_REAL speed, 
+		   int gi, 
+		   int gi_p) 
+{
+  int numvars, first_vi, first_vi_p; 
+
+  first_vi   = CCTK_FirstVarIndexI(gi);
+  first_vi_p = CCTK_FirstVarIndexI(gi_p);
+  numvars    = CCTK_NumVarsInGroupI(gi);
+
+  return ApplyBndRadiative(GH,stencil_size,NULL,var0,
+			   speed,first_vi,first_vi_p,numvars,dir);
+
+}
+
+void CCTK_FCALL CCTK_FNAME(BndRadiativeDirGI)
+     (int *ierr, cGH *GH, int *stencil_size, int *dir, CCTK_REAL *var0, 
+      CCTK_REAL *speed, int *gi, int *gi_p)
+{
+  *ierr = BndRadiativeDirGI(GH, *stencil_size, *dir, *var0, *speed, *gi, *gi_p);
+}
+
+/*@@
    @routine    BndRadiativeGI
    @date       Tue Jul 18 18:04:07 2000
    @author     Gerd Lanfermann
@@ -79,9 +120,9 @@ static int ApplyBndRadiative(cGH *GH,
 @@*/
 
 int BndRadiativeGI(cGH *GH, 
-		   int *sw,
+		   int *stencil_array,
 		   CCTK_REAL var0, 
-		   CCTK_REAL v0, 
+		   CCTK_REAL speed, 
 		   int gi, 
 		   int gi_p) 
 {
@@ -91,17 +132,72 @@ int BndRadiativeGI(cGH *GH,
   first_vi_p = CCTK_FirstVarIndexI(gi_p);
   numvars    = CCTK_NumVarsInGroupI(gi);
 
-  return ApplyBndRadiative(GH,sw,var0,v0,first_vi,first_vi_p,numvars);
+  return ApplyBndRadiative(GH,-1,stencil_array,var0,speed,first_vi,first_vi_p,numvars,0);
 
 }
 
 void CCTK_FCALL CCTK_FNAME(BndRadiativeGI)
-     (int *ierr, cGH *GH, int *sw, CCTK_REAL *var0, 
-      CCTK_REAL *v0, int *gi, int *gi_p)
+     (int *ierr, cGH *GH, int *stencil_array, CCTK_REAL *var0, 
+      CCTK_REAL *speed, int *gi, int *gi_p)
 {
-  *ierr = BndRadiativeGI(GH, sw, *var0, *v0, *gi, *gi_p);
+  *ierr = BndRadiativeGI(GH, stencil_array, *var0, *speed, *gi, *gi_p);
 }
 
+
+/* ===================================================================== */
+
+
+/*@@
+   @routine    BndRadiativeDirGN
+   @date       Sat Jan 20 2001
+   @author     Gabrielle Allen
+   @desc 
+      Apply radiative BC's by group name in given direction
+   @enddesc 
+   @calls     
+   @calledby   
+   @history 
+   @endhistory 
+
+@@*/
+
+int BndRadiativeDirGN(cGH *GH, 
+		      int stencil_size,
+		      int dir,
+		      CCTK_REAL var0, 
+		      CCTK_REAL speed, 
+		      const char *gn, 
+		      const char *gn_p)
+{
+  int retval;
+  int gi   = CCTK_GroupIndex(gn);
+  int gi_p = CCTK_GroupIndex(gn_p);
+  
+  if ((gi>-1)&&(gi_p>-1))
+    retval = BndRadiativeDirGI(GH, stencil_size, dir, var0, speed, gi, gi_p);
+  else 
+  {
+    CCTK_VWarn(1,__LINE__,__FILE__,"Boundary",
+	      "BndRadiativeDirGN: Grid variable groups %s or %s not found",
+	      gn, gn_p);
+    retval = -1;
+  }
+
+  return retval;
+}
+ 
+void CCTK_FCALL CCTK_FNAME(BndRadiativeDirGN)
+     (int *ierr, cGH *GH, int *stencil_size, int *dir, CCTK_REAL *var0, 
+      CCTK_REAL *speed, TWO_FORTSTRINGS_ARGS)
+{
+  TWO_FORTSTRINGS_CREATE(gn, gn_p)
+  *ierr = BndRadiativeDirGN(GH, *stencil_size, *dir, *var0, *speed, gn, gn_p);
+  free(gn);
+  free(gn_p);
+  return;
+}
+
+
 /*@@
    @routine    BndRadiativeGI
    @date       Tue Jul 18 18:04:07 2000
@@ -116,9 +212,9 @@ void CCTK_FCALL CCTK_FNAME(BndRadiativeGI)
 @@*/
 
 int BndRadiativeGN(cGH *GH, 
-		   int *sw, 
+		   int *stencil_array, 
 		   CCTK_REAL var0, 
-		   CCTK_REAL v0, 
+		   CCTK_REAL speed, 
 		   const char *gn, 
 		   const char *gn_p)
 {
@@ -127,7 +223,7 @@ int BndRadiativeGN(cGH *GH,
   int gi_p = CCTK_GroupIndex(gn_p);
   
   if ((gi>-1)&&(gi_p>-1))
-    retval = BndRadiativeGI(GH, sw, var0, v0, gi, gi_p);
+    retval = BndRadiativeGI(GH, stencil_array, var0, speed, gi, gi_p);
   else 
   {
     CCTK_VWarn(1,__LINE__,__FILE__,"Boundary",
@@ -140,18 +236,58 @@ int BndRadiativeGN(cGH *GH,
 }
  
 void CCTK_FCALL CCTK_FNAME(BndRadiativeGN)
-     (int *ierr, cGH *GH, int *sw, CCTK_REAL *var0, 
-      CCTK_REAL *v0, TWO_FORTSTRINGS_ARGS)
+     (int *ierr, cGH *GH, int *stencil_array, CCTK_REAL *var0, 
+      CCTK_REAL *speed, TWO_FORTSTRINGS_ARGS)
 {
   TWO_FORTSTRINGS_CREATE(gn, gn_p)
-  *ierr = BndRadiativeGN(GH, sw, *var0, *v0, gn, gn_p);
+  *ierr = BndRadiativeGN(GH, stencil_array, *var0, *speed, gn, gn_p);
   free(gn);
   free(gn_p);
   return;
 }
 
+
+/* ===================================================================== */
+
+
 /*@@
-   @routine    BndRadiativeGI
+   @routine    BndRadiativeDirVI
+   @date       Sat Jan 20 2001
+   @author     Gabrielle Allen
+   @desc 
+      Apply radiative BC's by variable index in given direction
+   @enddesc 
+   @calls     
+   @calledby   
+   @history 
+   @endhistory 
+
+@@*/
+
+int BndRadiativeDirVI(cGH *GH, 
+		      int stencil_size, 
+		      int dir,
+		      CCTK_REAL var0, 
+		      CCTK_REAL speed, 
+		      int vi, 
+		      int vi_p) 
+{
+  return ApplyBndRadiative(GH,stencil_size,NULL,var0,speed,vi,vi_p,1,dir);
+}
+
+void CCTK_FCALL CCTK_FNAME(BndRadiativeDirVI)
+     (int *ierr, cGH *GH, int *stencil_size, int *dir, CCTK_REAL *var0, 
+      CCTK_REAL *speed, int *vi, int *vi_p)
+{
+  *ierr = BndRadiativeDirVI(GH, *stencil_size, *dir, *var0, 
+			    *speed, *vi, *vi_p);
+  return;
+}
+
+
+
+/*@@
+   @routine    BndRadiativeVI
    @date       Tue Jul 18 18:04:07 2000
    @author     Gerd Lanfermann
    @desc 
@@ -165,23 +301,81 @@ void CCTK_FCALL CCTK_FNAME(BndRadiativeGN)
 @@*/
 
 int BndRadiativeVI(cGH *GH, 
-                    int *sw, 
+                    int *stencil_array, 
                     CCTK_REAL var0, 
-                    CCTK_REAL v0, 
+                    CCTK_REAL speed, 
                     int vi, 
                     int vi_p) 
 {
-  return ApplyBndRadiative(GH,sw,var0,v0,vi,vi_p,1);
+  return ApplyBndRadiative(GH,-1,stencil_array,var0,speed,vi,vi_p,1,0);
 }
 
 void CCTK_FCALL CCTK_FNAME(BndRadiativeVI)
-     (int *ierr, cGH *GH, int *sw, CCTK_REAL *var0, 
-      CCTK_REAL *v0, int *vi, int *vi_p)
+     (int *ierr, cGH *GH, int *stencil_array, CCTK_REAL *var0, 
+      CCTK_REAL *speed, int *vi, int *vi_p)
 {
-  *ierr = BndRadiativeVI(GH, sw, *var0, *v0, *vi, *vi_p);
+  *ierr = BndRadiativeVI(GH, stencil_array, *var0, *speed, *vi, *vi_p);
   return;
 }
 
+
+/* ======================================================================= */
+
+
+/*@@
+   @routine    BndRadiativeDirGI
+   @date       Sat Jan 20 2001
+   @author     Gabrielle Allen
+   @desc 
+      Apply radiative BC's by variable name in given direction
+   @enddesc 
+   @calls     
+   @calledby   
+   @history 
+   @endhistory 
+
+@@*/
+
+int BndRadiativeDirVN(cGH *GH, 
+                   int stencil_size,
+		   int dir,
+                   CCTK_REAL var0, 
+                   CCTK_REAL speed, 
+                   const char *vn, 
+                   const char *vn_p) 
+{
+  int retval;
+  int vi   = CCTK_VarIndex(vn);
+  int vi_p = CCTK_VarIndex(vn_p);
+  
+  if (vi>=0 && vi_p>=0)
+  {
+    retval = BndRadiativeDirVI(GH, stencil_size, dir, var0, speed, vi, vi_p);
+  }
+  else
+  {
+    CCTK_VWarn(1,__LINE__,__FILE__,"Boundary",
+	      "BndRadiativeDirVN: Grid variable %s or %s not found",
+	      vn, vn_p);
+    retval = -1;
+  }
+
+  return retval;
+
+}
+
+void CCTK_FCALL CCTK_FNAME(BndRadiativeDirVN)
+     (int *ierr, cGH *GH, int *stencil_size, int *dir, CCTK_REAL *var0, 
+      CCTK_REAL *speed, TWO_FORTSTRINGS_ARGS)
+{
+  TWO_FORTSTRINGS_CREATE(vn,vn_p)
+  *ierr = BndRadiativeDirVN(GH, *stencil_size, *dir, *var0, *speed, vn, vn_p);
+  free(vn);
+  free(vn_p);
+  return;
+}  
+
+
 /*@@
    @routine    BndRadiativeGI
    @date       Tue Jul 18 18:04:07 2000
@@ -197,9 +391,9 @@ void CCTK_FCALL CCTK_FNAME(BndRadiativeVI)
 @@*/
 
 int BndRadiativeVN(cGH *GH, 
-                   int *sw,
+                   int *stencil_array,
                    CCTK_REAL var0, 
-                   CCTK_REAL v0, 
+                   CCTK_REAL speed, 
                    const char *vn, 
                    const char *vn_p) 
 {
@@ -209,7 +403,7 @@ int BndRadiativeVN(cGH *GH,
   
   if (vi>=0 && vi_p>=0)
   {
-    retval = BndRadiativeVI(GH, sw, var0, v0, vi, vi_p);
+    retval = BndRadiativeVI(GH, stencil_array, var0, speed, vi, vi_p);
   }
   else
   {
@@ -224,11 +418,11 @@ int BndRadiativeVN(cGH *GH,
 }
 
 void CCTK_FCALL CCTK_FNAME(BndRadiativeVN)
-     (int *ierr, cGH *GH, int *sw, CCTK_REAL *var0, 
-      CCTK_REAL *v0, TWO_FORTSTRINGS_ARGS)
+     (int *ierr, cGH *GH, int *stencil_array, CCTK_REAL *var0, 
+      CCTK_REAL *speed, TWO_FORTSTRINGS_ARGS)
 {
   TWO_FORTSTRINGS_CREATE(vn,vn_p)
-  *ierr = BndRadiativeVN(GH, sw, *var0, *v0, vn, vn_p);
+  *ierr = BndRadiativeVN(GH, stencil_array, *var0, *speed, vn, vn_p);
   free(vn);
   free(vn_p);
   return;
@@ -240,7 +434,7 @@ void CCTK_FCALL CCTK_FNAME(BndRadiativeVN)
  *********************     Local Routines   *************************
  ********************************************************************/
 /*@@
-   @routine    BndRadiativeGI
+   @routine    ApplyBndRadiative
    @date       Tue Jul 18 18:04:07 2000
    @author     Gerd Lanfermann
    @desc 
@@ -254,44 +448,117 @@ void CCTK_FCALL CCTK_FNAME(BndRadiativeVN)
 @@*/
 
 static int ApplyBndRadiative(cGH *GH, 
-			     int *sw,
+			     int stencil,
+			     int *stencil_array,
 			     CCTK_REAL var0, 
-			     CCTK_REAL v0, 
+			     CCTK_REAL speed, 
 			     int var_current, 
 			     int var_previous,
-			     int num_vars) 
+			     int num_vars,
+			     int dir) 
 {
-  SymmetryGHex *sGHex;
   int xi=-1;
   int yi=-1;
   int zi=-1;
   int ri=-1;    
+  int count;
   int gi;
   int vi1, vi2;
   int symmetry_handle;
   int idim, dim;
   int time, time_p;
   int berr,ierr;
-  int *doBC, *dstag, *lssh;
+  int dirloop1,dirloop2;
+  int doBC;
+  int retval;
+  int *sw;
+  int *dstag;
+  int *lssh;
   CCTK_REAL dx[3], dt;
+  SymmetryGHex *sGHex;
+
+#ifdef DEBUG_BOUNDARY
+  printf("Input arguments for ApplyBndRadiative:\n");
+  printf("GH = %x\n",GH);
+  printf("stencil = %d\n",stencil);
+  printf("stencil_array = %x\n",stencil_array);
+  printf("var0 = %f\n",var0);
+  printf("speed = %f\n",speed);
+  printf("var_current = %d\n",var_current);
+  printf("var_previous = %d\n",var_previous);
+  printf("num_vars = %d\n",num_vars); 
+  printf("dir = %d\n",dir);
+  printf("-----------------------------------\n");
+#endif
+  
+  /* Get group dimension */
+  dim = CCTK_GroupDimFromVarI(var_current);
 
-  /* See if we have a symmetry array */
-  symmetry_handle = CCTK_GHExtensionHandle("Symmetry");
-  if (symmetry_handle < 0)
+  if (dim<3 || dim>3)
   {
-    sGHex = NULL;
+    CCTK_VWarn(1,__LINE__,__FILE__,CCTK_THORNSTRING,
+	       "ApplyBndRadiative: Invalid grid function dimension %d",
+	       dim);
+    retval=-1;
   }
   else
   {
-    sGHex  = (SymmetryGHex*)GH->extensions[symmetry_handle];
-  }
+    /* See if we have a symmetry array */
+    symmetry_handle = CCTK_GHExtensionHandle("Symmetry");
+    if (symmetry_handle < 0)
+    {
+      sGHex = NULL;
+    }
+    else
+    {
+      sGHex  = (SymmetryGHex*)GH->extensions[symmetry_handle];
+    }
 
-  /* Get group dimension */
-  dim = CCTK_GroupDimFromVarI(var_current);
+    /* Set up stencil width array if needed */
+    if (stencil_array && dir == 0)
+    {
+      sw = stencil_array;
+    }
+    else if (dir==0)
+    {
+      CCTK_WARN(0,"ApplyBndRadiative: Direction 0 invalid for directional BC");
+    }
+    else if (dir>dim||-dir>dim)
+    {
+      CCTK_VWarn(0,__LINE__,__FILE__,CCTK_THORNSTRING,
+		 "ApplyBndRadiative: Direction %d greater than dimension %d",
+		 dir,dim);
+    }
+    else
+    {
 
-  /* Radiative boundaries need the underlying Cartesian coordinates */
-  switch (dim) 
-  {
+      int i;
+      /* Should be single direction */
+      sw  = (int *)malloc(dim*sizeof(int));
+      if (sw)
+      {
+	for (i=0;i<dim;i++)
+	{
+	  sw[i] = 0;
+	}
+	if (dir > 0)
+	{
+	  sw[dir-1] = stencil;
+	}
+	else
+	{
+	  sw[-dir-1] = stencil;
+	}
+      }
+      else
+      {
+	CCTK_WARN(0,"ApplyBndRadiative: Malloc sw failed");
+      }
+    }
+
+    /* Radiative boundaries need the underlying Cartesian coordinates */
+    switch (dim) 
+    {
     case 1: 
       xi = CCTK_CoordIndex(-1,"x","cart1d");
       break;
@@ -306,130 +573,148 @@ static int ApplyBndRadiative(cGH *GH,
       ri = CCTK_CoordIndex(-1,"r","spher3d");
       break;
     default:
-      CCTK_VWarn(1,__LINE__,__FILE__,"Boundary",
-		 "Radiative boundary only for dim=3: grid varible '%s'",
+      CCTK_VWarn(1,__LINE__,__FILE__,CCTK_THORNSTRING,
+		 "ApplyBndRadiative: Radiative boundary only for dim=3: grid variable '%s'",
 		 CCTK_VarName(var_current));
       break;
-  }
+    }
 
+    /* Allocate temporary arrays */
+    dstag = (int *)malloc(dim*sizeof(int));
+    lssh  = (int *)malloc(dim*sizeof(int));
 
-  /* Allocate temporary arrays */
-  doBC  = (int *)malloc((2*dim)*sizeof(int)); 
-  dstag = (int *)malloc(dim*sizeof(int));
-  lssh  = (int *)malloc(dim*sizeof(int));
+    /* get the directional staggering of the group */
+    gi   = CCTK_GroupIndexFromVarI(var_current);
+    ierr = CCTK_GroupStaggerDirArrayGI(dstag, dim, gi);
 
-  /* get the directional staggering of the group */
-  gi   = CCTK_GroupIndexFromVarI(var_current);
-  ierr = CCTK_GroupStaggerDirArrayGI(dstag, dim, gi);
+    
+    /* Use next time level, if available */
+    time = CCTK_NumTimeLevelsFromVarI(var_current) - 1;
+    if (time < 0)
+    {
+      time = 0;
+    }
+    /* Use current time level, if available */
+    time_p = CCTK_NumTimeLevelsFromVarI(var_previous) - 2;
+    if (time_p < 0)
+    {
+      time_p = 0;
+    }
 
+    /* Treat all boundaries or just one? */
+    if (dir>0)
+    {
+      dirloop1 = 2*(dir-1)+1;
+      dirloop2 = dirloop1+1;
+    }
+    else if (dir<0)
+    {
+      dirloop1 = 2*(-dir-1);
+      dirloop2 = dirloop1+1;
+    }
+    else
+    {
+      dirloop1 = 0;
+      dirloop2 = 2*dim;
+    }
+  
+    for (vi1=var_current; vi1<var_current+num_vars; vi1++) 
+    {
+      
+      vi2 = (vi1-var_current+var_previous);
 
-  /* Use next time level, if available */
-  time = CCTK_NumTimeLevelsFromVarI(var_current) - 1;
-  if (time < 0)
-  {
-    time = 0;
-  }
-  /* Use current time level, if available */
-  time_p = CCTK_NumTimeLevelsFromVarI(var_previous) - 2;
-  if (time_p < 0)
-  {
-    time_p = 0;
-  }
+      /* Apply condition if:
+	 + boundary is not a symmetry boundary 
+	 (no symmetry or unset(=unsed))
+	 + boundary is a physical boundary
+	 + have enough grid points            */
 
-  for (vi1=var_current; vi1<var_current+num_vars; vi1++) 
-  {
-     
-    vi2 = (vi1-var_current+var_previous);
-
-    /* Apply condition if:
-       + boundary is not a symmetry boundary 
-       (no symmetry or unset(=unsed))
-       + boundary is a physical boundary
-       + have enough grid points
-    */
-    if (sGHex)
-    {
-      for (idim=0;idim<dim;idim++) 
+      for (idim=0;idim<dim;idim++)
       {
-        doBC[idim*2]   = (((sGHex->GFSym[vi1][idim*2]==GFSYM_NOSYM)||
-                           (sGHex->GFSym[vi1][idim*2]==GFSYM_UNSET)) &&
-                          GH->cctk_lsh[idim]>1 && GH->cctk_bbox[idim*2]);
-        doBC[idim*2+1] = (((sGHex->GFSym[vi1][idim*2+1]==GFSYM_NOSYM)||
-                           (sGHex->GFSym[vi1][idim*2+1]==GFSYM_UNSET)) &&
-                          GH->cctk_lsh[idim]>1 && GH->cctk_bbox[idim*2+1]);
 	/* FIXME: THIS LINE FAILS WAVE TESTSUITES */
         /*lssh[idim] = GH->cctk_lssh[CCTK_LSSH_IDX(dstag[idim],idim)];*/
 	lssh[idim] = GH->cctk_lsh[idim];
       }
-    }
-    else
-    {
-      for (idim=0;idim<dim;idim++) 
+
+      for (count=dirloop1;count<dirloop2;count++)
       {
-        doBC[idim*2]   = (GH->cctk_lsh[idim]>1 && GH->cctk_bbox[idim*2]);
-        doBC[idim*2+1] = (GH->cctk_lsh[idim]>1 && GH->cctk_bbox[idim*2+1]);
-        lssh[idim]     = GH->cctk_lssh[CCTK_LSSH_IDX(dstag[idim],idim)];
+	if (sGHex)
+	{
+	  doBC = (
+		  ( (sGHex->GFSym[vi1][count]==GFSYM_NOSYM) ||
+		    (sGHex->GFSym[vi1][count]==GFSYM_UNSET) ) &&
+		  GH->cctk_lsh[count/2]>1 && GH->cctk_bbox[count]) ?
+	    count : -1;
+	}
+	else
+	{
+	  doBC  = (GH->cctk_lsh[count/2]>1 && GH->cctk_bbox[count]) ?
+	    count : -1 ;
+	}
+	
+	if (doBC >=0)
+	{
+	  switch (dim)
+	  {
+	  case 1: 
+	    berr = -1; 
+	    CCTK_VWarn(1,__LINE__,__FILE__,CCTK_THORNSTRING,
+		       "ApplyBndRadiative: Radiative boundary conditions" 
+		       "not implemented for dimension = 1, grid variable %s",
+		       CCTK_VarName(vi1));
+	    break;
+	  case 2: 
+	    berr = -1; 
+	    CCTK_VWarn(1,__LINE__,__FILE__,CCTK_THORNSTRING,
+		       "ApplyBndRadiative: Radiative boundary conditions"
+		       "not implemented for dimension = 2, grid variable %s",
+		       CCTK_VarName(vi1));
+	    break;	  
+	  case 3:   
+	    /* According to the Cactus spec, the true dx and dt values for a
+	       grid are calculated as follows: */
+	    dx[0] = GH->cctk_delta_space[0] / GH->cctk_levfac[0];
+	    dx[1] = GH->cctk_delta_space[1] / GH->cctk_levfac[1];
+	    dx[2] = GH->cctk_delta_space[2] / GH->cctk_levfac[2];
+	    dt    = GH->cctk_delta_time;
+    
+	    berr = BndApplyRadiative3Di
+	      (GH,
+	       doBC,		      /* apply BC on which face */
+	       lssh, 	              /* local shape, respects staggering */
+	       sw,	              /* stencil width array */
+	       dx,	              /* dx,dy,dz */
+	       dt,		      /* dt */
+	       GH->data[vi1][time  ], /* start of data array GF[] */
+	       GH->data[vi2][time_p], /* start of prev data array */
+	       GH->data[xi][0],       /* pointer to the X field */
+	       GH->data[yi][0],       /* pointer to the Y field */
+	       GH->data[zi][0],       /* pointer to the Z field */
+	       GH->data[ri][0],       /* pointer to the R field */
+	       var0,		      /* asymptotic limit */
+	       speed);		      /* wave speed */
+	    break;
+	  default : 
+	    berr=-1; 
+	    CCTK_VWarn(1,__LINE__,__FILE__,CCTK_THORNSTRING,
+		       "ApplyBndRadiative: No Radiative BC for dim>3, grid variable %s",
+		       CCTK_VarName(vi1));
+	  }
+	}
       }
+      retval=(berr<0)?-1:0;
     }
+    free(dstag);
+    free(lssh);
 
-    switch (dim) 
+    if (!stencil_array)
     {
-      case 1: 
-	berr = -1; 
-	CCTK_VWarn(1,__LINE__,__FILE__,"Boundary",
-		   "ApplyBndRadiative: Radiative boundary conditions" 
-		   "not implemented for dimension = 1, grid variable '%s'",
-		   CCTK_VarName(vi1));
-	break;
-      case 2: 
-	berr = -1; 
-	CCTK_VWarn(1,__LINE__,__FILE__,"Boundary",
-		   "ApplyBndRadiative: Radiative boundary conditions"
-		   "not implemented for dimension = 2, grid variable '%s'",
-		   CCTK_VarName(vi1));
-	break;
-      case 3:   
-	/* According to the Cactus spec, the true dx and dt values for a
-		   grid are calculated as follows: */
-	dx[0] = GH->cctk_delta_space[0] / GH->cctk_levfac[0];
-	dx[1] = GH->cctk_delta_space[1] / GH->cctk_levfac[1];
-	dx[2] = GH->cctk_delta_space[2] / GH->cctk_levfac[2];
-	dt    = GH->cctk_delta_time;
-    
-	berr = BndApplyRadiative3Di
-	  (GH,
-	   dim,
-	   sw,		          /* stencil_width */
-	   doBC,		  /* flag applying low/up BCs */
-	   lssh, 	          /* local shape, respects staggering */
-	   dx,		          /* dx,dy,dz */
-	   dt,		          /* dt */
-	   GH->data[vi1][time  ], /* start of data array GF[] */
-	   GH->data[vi2][time_p], /* start of prev data array */
-	   GH->data[xi][0],       /* pointer to the X field */
-	   GH->data[yi][0],       /* pointer to the Y field */
-	   GH->data[zi][0],       /* pointer to the Z field */
-	   GH->data[ri][0],       /* pointer to the R field */
-	   var0,		  /* asymptotic limit */
-	   v0);		          /* wave speed */
-	break;
-      default : 
-	berr=-1; 
-	CCTK_VWarn(1,__LINE__,__FILE__,"Boundary",
-		   "No Radiative BC for dim>3, grid variable '%s'",
-		   CCTK_VarName(vi1));
+      free(sw);
     }
-    berr=(berr<0)?-1:0;
   }
-  
-  free(lssh);
-  free(dstag);
-  free(doBC);
-  
-  return ierr;
-}
-
 
+  return retval;
+}
 
 /*@@
    @routine BndApplyRadiative3Di
@@ -515,10 +800,9 @@ static int ApplyBndRadiative(cGH *GH,
 @@*/
 
 static int BndApplyRadiative3Di(cGH *GH,
-				int  gdim,
-				int *sw,  
-				int *doBC, 
+				int doBC, 
 				int *lssh,
+				int *sw,  
 				CCTK_REAL *dxyz, 
 				CCTK_REAL dt, 
 				CCTK_REAL *var_n, 
@@ -528,7 +812,7 @@ static int BndApplyRadiative3Di(cGH *GH,
 				CCTK_REAL *z, 
 				CCTK_REAL *r,
 				CCTK_REAL var0, 
-				CCTK_REAL v0) 
+				CCTK_REAL speed) 
 {
 
   DECLARE_CCTK_PARAMETERS
@@ -588,7 +872,7 @@ static int BndApplyRadiative3Di(cGH *GH,
 
   /* Find Courant parameters. */
 
-  dtv  = v0*dt;
+  dtv  = speed*dt;
   dtvh = half*dtv;
 
   dtvvar0 = dtv*var0;
@@ -610,11 +894,17 @@ static int BndApplyRadiative3Di(cGH *GH,
 
   /* Lower x-bound: x(2,:,:) --> x[xgp2] */
 
-  if (doBC[0]==1) {
-    for (k=0;k<lssh2;k++) {
-      for (j=0;j<lssh1;j++) {
-	for (i=sw0-1;i>=0;i--) {
-
+  if (doBC == 0) 
+  {
+#ifdef DEBUG_BOUNDARY
+    printf("Applying radiative boundary (lower x)\n");
+#endif
+    for (k=0;k<lssh2;k++) 
+    {
+      for (j=0;j<lssh1;j++) 
+      {
+	for (i=sw0-1;i>=0;i--) 
+	{
 	  xgp0 = CCTK_GFINDEX3D(GH,i  ,j,k);
 	  xgp1 = CCTK_GFINDEX3D(GH,i+1,j,k);
 	  xgp2 = CCTK_GFINDEX3D(GH,i+2,j,k);
@@ -633,11 +923,14 @@ static int BndApplyRadiative3Di(cGH *GH,
           ri0 = 1.0/r0;
           ri1 = 1.0/r1;
 
-          if (fac==0) {
+          if (fac==0) 
+	  {
 
              H = 0;
 
-          } else {
+          } 
+	  else 
+	  {
 
              H = 0;
 
@@ -670,10 +963,17 @@ static int BndApplyRadiative3Di(cGH *GH,
 
   /* Upper x-bound: x(nx,:,:) --> xgp[xgp0] */
 
-  if (doBC[1]==1){
-    for (k=0;k<lssh2;k++) {
-      for (j=0;j<lssh1;j++) {
-	for (i=lssh0-sw0;i<lssh0;i++) {
+  if (doBC == 1)
+  {
+#ifdef DEBUG_BOUNDARY
+    printf("Applying radiative boundary (upper x)\n");
+#endif
+    for (k=0;k<lssh2;k++) 
+    {
+      for (j=0;j<lssh1;j++) 
+      {
+	for (i=lssh0-sw0;i<lssh0;i++) 
+	{
 
 	  xgp0 = CCTK_GFINDEX3D(GH,i  ,j,k);
 	  xgp1 = CCTK_GFINDEX3D(GH,i-1,j,k);
@@ -693,11 +993,14 @@ static int BndApplyRadiative3Di(cGH *GH,
           ri0 = 1.0/r0;
           ri1 = 1.0/r1;
 
-          if (fac==0) {
+          if (fac==0) 
+	  {
 
              H = 0;
 
-          } else {
+          } 
+	  else 
+	  {
 
              H = 0;
 
@@ -730,8 +1033,11 @@ static int BndApplyRadiative3Di(cGH *GH,
 
   /* Lower y-bound */
 
-  if (doBC[2] == 1) 
+  if (doBC == 2) 
   {
+#ifdef DEBUG_BOUNDARY
+    printf("Applying radiative boundary (lower y)\n");
+#endif
     for (k=0;k<lssh2;k++) 
     {
       for (i=0;i<lssh0;i++) 
@@ -756,11 +1062,14 @@ static int BndApplyRadiative3Di(cGH *GH,
           ri0 = 1.0/r0;
           ri1 = 1.0/r1;
 
-          if (fac==0) {
+          if (fac==0) 
+	  {
 
              H = 0;
 
-          } else {
+          } 
+	  else 
+	  {
 
              H = 0;
 
@@ -793,8 +1102,11 @@ static int BndApplyRadiative3Di(cGH *GH,
 
   /* Upper y bound */
 
-  if (doBC[3] == 1) 
+  if (doBC == 3) 
   {
+#ifdef DEBUG_BOUNDARY
+    printf("Applying radiative boundary (upper y)\n");
+#endif
     for (k=0;k<lssh2;k++) 
     {
       for (i=0;i<lssh0;i++) 
@@ -860,10 +1172,17 @@ static int BndApplyRadiative3Di(cGH *GH,
 
   /* Lower z-bound */
 
-  if (doBC[4]==1) {
-    for (j=0;j<lssh1;j++) {
-      for (i=0;i<lssh0;i++) {
-	for (k=sw2-1;k>=0;k--) {
+  if (doBC==4) 
+  {
+#ifdef DEBUG_BOUNDARY
+    printf("Applying radiative boundary (lower z)\n");
+#endif
+    for (j=0;j<lssh1;j++) 
+    {
+      for (i=0;i<lssh0;i++) 
+      {
+	for (k=sw2-1;k>=0;k--) 
+	{
 
 	  zgp0 = CCTK_GFINDEX3D(GH,i,j,k  );
 	  zgp1 = CCTK_GFINDEX3D(GH,i,j,k+1);
@@ -883,11 +1202,14 @@ static int BndApplyRadiative3Di(cGH *GH,
           ri0 = 1.0/r0;
           ri1 = 1.0/r1;
 
-          if (fac==0) {
+          if (fac==0) 
+	  {
 
              H = 0;
 
-          } else {
+          } 
+	  else 
+	  {
 
              H = 0;
 
@@ -920,10 +1242,17 @@ static int BndApplyRadiative3Di(cGH *GH,
 
   /* Upper z-bound */
 
-  if (doBC[5] == 1) {
-    for (j=0;j<lssh1;j++) {
-      for (i=0;i<lssh0;i++) {
-	for (k=lssh2-sw2;k<lssh2;k++) {
+  if (doBC ==  5) 
+  {
+#ifdef DEBUG_BOUNDARY
+    printf("Applying radiative boundary (upper z)\n");
+#endif
+    for (j=0;j<lssh1;j++) 
+    {
+      for (i=0;i<lssh0;i++) 
+      {
+	for (k=lssh2-sw2;k<lssh2;k++) 
+	{
 
 	  zgp0 = CCTK_GFINDEX3D(GH,i,j,k  );
 	  zgp1 = CCTK_GFINDEX3D(GH,i,j,k-1);
@@ -943,11 +1272,14 @@ static int BndApplyRadiative3Di(cGH *GH,
           ri0 = 1.0/r0;
           ri1 = 1.0/r1;
 
-          if (fac==0) {
+          if (fac==0) 
+	  {
 
              H = 0;
 
-          } else {
+          } 
+	  else 
+	  {
 
              H = 0;