Implement wrapper for radiation boundaries which allows all available options.

Options are now specified via:

key			value
---			-----
PREVIOUS TIME LEVEL	name or index of gv holding previous time level
			(CCTK_STRING or CCTK_INT respectively)
			(to be deprecated, please use proper time levels)
LIMIT			value of field at infinity (CCTK_REAL)
SPEED			wave speed (CCTK_REAL)
STENCIL WIDTH		width of stencil in each dimension (CCTK_ARRAY of d
			CCTK_INTs, for a d-dimensional gv)


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

1 files changed, 125 insertions, 50 deletions

diff --git a/src/RadiationBoundary.c b/src/RadiationBoundary.c
index 4eb2e5f..7f83635 100644
--- a/src/RadiationBoundary.c
+++ b/src/RadiationBoundary.c
@@ -71,6 +71,9 @@
 #include <string.h>
 
 #include "cctk.h"
+#include "cctk_Faces.h"
+#include "util_Table.h"
+#include "util_ErrorCodes.h"
 #include "cctk_FortranString.h"
 #include "cctk_Parameters.h"
 
@@ -81,6 +84,16 @@
 static const char *rcsid = "$Header$";
 CCTK_FILEVERSION(CactusBase_Boundary_RadiationBoundary_c);
 
+static int ApplyBndRadiative (const cGH *GH,
+                              int stencil_dir,
+                              const int *stencil_alldirs,
+                              int dir,
+                              CCTK_REAL var0,
+                              CCTK_REAL speed,
+                              int first_var_to,
+                              int first_var_from,
+                              int num_vars);
+
 /********************************************************************
  ********************    External Routines   ************************
  ********************************************************************/
@@ -93,11 +106,7 @@ CCTK_FILEVERSION(CactusBase_Boundary_RadiationBoundary_c);
                Top level function which is registered as handling 
 	       this boundary condition
    @enddesc
-   @calls      BndRadiativeVI, for the moment
-   @history
-               Perliminary version for testing, simply calls BndRadiativeVI 
-	       with hardwired arguments
-   @endhistory
+   @calls      ApplyBndRadiative
    @var        GH
    @vdesc      Pointer to CCTK grid hierarchy
    @vtype      const cGH *
@@ -114,6 +123,11 @@ CCTK_FILEVERSION(CactusBase_Boundary_RadiationBoundary_c);
    @vtype      int *
    @vio        in
    @endvar
+   @var        faces
+   @vdesc      array of set of faces to which to apply the bc
+   @vtype      int
+   @vio        in
+   @endvar
    @var        table_handles
    @vdesc      array of table handles which hold extra arguments
    @vtype      int
@@ -121,51 +135,133 @@ CCTK_FILEVERSION(CactusBase_Boundary_RadiationBoundary_c);
    @endvar
    @returntype int
    @returndesc
-               return code of @seeroutine BndRadiativeVI <BR>
-	       |= -16 if valid table handle passed
+               return code of @seeroutine ApplyBndRadiative <BR>
    @endreturndesc
 @@*/
 
-int BndRadiative(const cGH *GH, int num_vars, int *var_indices, 
-		 int *table_handles) 
+int BndRadiative(const cGH *GH, int num_vars, int *vars, int *faces,
+		 int *tables) 
 {
-  int i, retval=0;
-  int sw[3] = {1,1,1};
+  int i, j, k, value_type, value_size, err, prev_time_level,
+    gdim, max_gdim, retval;
+  char *prev_time_level_name;
+
+  /* variables to pass to ApplyBndRadiative */
+  /*int stencil_dir;*/ /* width of stencil in direction dir, not needed */
+  int *stencil_alldirs; /* width of stencil in all directions */
+  int dir; /* direction in which to apply bc */
+  CCTK_REAL limit, speed;
 
 #ifdef DEBUG
   printf("BndRadiative(): got passed GH=%p, num_vars=%d, var_indices[0]=%d, table_handles[0]=%d\n", (const void *) GH, num_vars, var_indices[0], table_handles[0]);
 #endif
 
-  /* loop through each variable (Is there some way to take advantage
-     of knowing the entire list of variables ahead of time?) */
-  for (i=0; i<num_vars; ++i) {
+  retval = 0; stencil_alldirs = NULL;
+
+  /* loop through variables, j at a time */
+  for (i=0; i<num_vars; i+=j) {
 
-    /* Check the input arguments */
-    if (table_handles[i]>=0) {
+    /* find other adjacent vars which are selected for identical bcs */
+    j=1;
+    while (i+j<num_vars && vars[i+j]==vars[i]+j && tables[i+j]==tables[i] && 
+	   faces[i+j]==faces[i])
+    {
+      ++j;
+    }
+
+    /* Check to see if faces specification is valid */
+    if (faces[i] != CCTK_ALL_FACES)
+    {
       CCTK_VWarn(1, __LINE__, __FILE__, CCTK_THORNSTRING,
-		 "Ignoring table handle %d, because parameters are simply "
-		 "hardwired at the moment.  Please complain.", i);
-      retval |= -16;
+		 "Faces specification %d for radiative boundary conditions on "
+		 "%s is not implemented yet.  "
+		 "Applying radiative bcs to all (external) faces.", faces[i],
+		 CCTK_VarName(vars[i]));
     }
+    dir = 0; /* apply bc to all faces */
 
-    /* hardwiring: stencil width
-                   constant var0 = 0.0
-	       	   wave speed = 1.0
-		   same var for from and to */
-    if ((retval = BndRadiativeVI(GH, sw, 0., 1., var_indices[i], 
-			      var_indices[i])) < 0) {
+    /* Set up default arguments for ApplyBndRadiative */
+    /* Allocate memory and populate stencil width array */
+    gdim = CCTK_GroupDimFromVarI(vars[i]);
+    if (!stencil_alldirs)
+    {
+      stencil_alldirs = (int *) malloc(gdim*sizeof(int));
+      max_gdim = gdim;
+    } else if (gdim > max_gdim)
+    {
+      realloc(stencil_alldirs, gdim*sizeof(int));
+      max_gdim = gdim;
+    }
+    for (k=0; k<gdim; ++k)
+    {
+      stencil_alldirs[k] = 1;
+    }
+    /* Defaults for remainder of arguments */
+    limit = 0.;
+    prev_time_level = vars[i];
+    speed = 1.;
+
+    /* Look on table for possible non-default arguments
+     * (If any of these table look-ups fail, the value will be unchanged 
+     * from its default value)
+     */
+    /* Asymptotic value of function at infinity */
+    err = Util_TableGetReal(tables[i], &limit, "LIMIT");
+    if (err == UTIL_ERROR_BAD_HANDLE)
+    {
+      CCTK_VWarn(3, __LINE__, __FILE__, CCTK_THORNSTRING,
+		 "Invalid table handle passed for Radiative boundary "
+		 "conditions for %s.  Using all default values.", 
+		 CCTK_VarName(vars[i]));
+    } else
+    { 
+      /* Previous time level (for GVs which don't use Cactus time levels) 
+	 (to be deprecated) */
+      if (Util_TableQueryValueInfo(tables[i], &value_type, &value_size, 
+				   "PREVIOUS TIME LEVEL"))
+      {
+	/* Key for previous time level is set */
+	if (value_type==CCTK_VARIABLE_STRING) 
+	{
+	  prev_time_level_name = (char *) malloc(value_size*sizeof(char));
+	  Util_TableGetString(tables[i], value_size, prev_time_level_name,
+			      "PREVIOUS TIME LEVEL");
+	  prev_time_level = CCTK_VarIndex(prev_time_level_name);
+	} else if (value_type==CCTK_VARIABLE_INT)
+	{
+	  Util_TableGetInt(tables[i], &prev_time_level,
+			   "PREVIOUS TIME LEVEL");
+	} else
+	{
+	  CCTK_VWarn(1, __LINE__, __FILE__, CCTK_THORNSTRING,
+		     "Invalid data type for key \"PREVIOUS TIME LEVEL\" "
+		     "Please use CCTK_STRING for the variable name, "
+		     "or CCTK_INT for the variable index.");
+	}
+      }
+
+      /* Stencil width array */
+      Util_TableGetIntArray(tables[i], gdim, stencil_alldirs, "STENCIL WIDTH");
+    
+      /* Wave speed */
+      Util_TableGetReal(tables[i], &speed, "SPEED");
+    }
+
+    if ((retval = ApplyBndRadiative(GH, 0, stencil_alldirs, dir, limit, speed,
+				    vars[i], prev_time_level, j)) < 0)
+    {
       CCTK_VWarn(1, __LINE__, __FILE__, CCTK_THORNSTRING,
-		 "BndRadiativeVI() returned %d", retval);
+		 "ApplyBndRadiative() returned %d", retval);
     }
   }
 #ifdef DEBUG
   printf("BndRadiative(): returning %d\n",retval);
 #endif
+  free(stencil_alldirs);
 
   return retval;
 }
 
-
 /* prototypes for external C routines are declared in header Boundary.h
    here only follow the fortran wrapper prototypes */
 void CCTK_FCALL CCTK_FNAME (BndRadiativeDirGI)
@@ -234,21 +330,9 @@ void CCTK_FCALL CCTK_FNAME (BndRadiativeVN)
                             TWO_FORTSTRING_ARG);
 
 
-
 /********************************************************************
  ********************    Internal Routines   ************************
  ********************************************************************/
-static int ApplyBndRadiative (const cGH *GH,
-                              int stencil_dir,
-                              const int *stencil_alldirs,
-                              int dir,
-                              CCTK_REAL var0,
-                              CCTK_REAL speed,
-                              int first_var_to,
-                              int first_var_from,
-                              int num_vars);
-
-
 /*@@
    @routine    BndRadiativeDirGI
    @date       Sat Jan 20 2001
@@ -1440,7 +1524,7 @@ static int ApplyBndRadiative (const cGH *GH,
                               int first_var_from,
                               int num_vars)
 {
-  DECLARE_CCTK_PARAMETERS
+  DECLARE_CCTK_PARAMETERS;
   int i, gdim, indx;
   int var_to, var_from;
   int timelvl_to, timelvl_from;
@@ -1490,11 +1574,6 @@ static int ApplyBndRadiative (const cGH *GH,
 
   /* Use next time level, if available */
   timelvl_to = 0;
-  /*  if (timelvl_to < 0)
-  {
-    timelvl_to = 0;
-    }*/
-  /* Use current time level, if available */
   if (CCTK_NumTimeLevelsFromVarI(first_var_from) > 1)
   {
     timelvl_from = 1;
@@ -1504,11 +1583,6 @@ static int ApplyBndRadiative (const cGH *GH,
     timelvl_from = 0;
   }
 
-  /*if (timelvl_from < 0)
-  {
-    timelvl_from = 0;
-    }*/
-
   /* Find Courant parameters. */
   dtv      = speed * GH->cctk_delta_time;
   dtvh     = 0.5 * dtv;
@@ -1538,6 +1612,7 @@ static int ApplyBndRadiative (const cGH *GH,
     offset[i] = i == 0 ? 1 : offset[i-1] * GH->cctk_lsh[i-1];
   }
 
+  /* Append r grid variable to end of coords[] array */
   sprintf (coord_system_name, "spher%dd", gdim);
   indx = CCTK_CoordIndex (-1, "r", coord_system_name);
   if (indx < 0)