Add tensor types to handle the Weyl scalars.

Check that the boundary specification is consistent between this thorn
and CoordBase, when CoordBase is used.


git-svn-id: http://svn.cactuscode.org/arrangements/CactusNumerical/ReflectionSymmetry/trunk@17 082bdb00-0f4f-0410-b49e-b1835e5f2039

3 files changed, 190 insertions, 9 deletions

diff --git a/interface.ccl b/interface.ccl
index ae0af7d..d1c90c4 100644
--- a/interface.ccl
+++ b/interface.ccl
@@ -68,3 +68,13 @@ CCTK_INT FUNCTION Boundary_SelectedGVs \
    CCTK_INT ARRAY OUT table_handles, \
    CCTK_STRING IN  bc_name)
 REQUIRES FUNCTION Boundary_SelectedGVs
+
+
+
+CCTK_INT FUNCTION GetBoundarySpecification \
+  (CCTK_INT IN size, \
+   CCTK_INT OUT ARRAY nboundaryzones, \
+   CCTK_INT OUT ARRAY is_internal, \
+   CCTK_INT OUT ARRAY is_staggered, \
+   CCTK_INT OUT ARRAY shiftout)
+USES FUNCTION GetBoundarySpecification
diff --git a/src/apply.c b/src/apply.c
index 6170234..50854d0 100644
--- a/src/apply.c
+++ b/src/apply.c
@@ -146,7 +146,7 @@ BndReflectVI (cGH const * restrict const cctkGH,
   
   int table;
   char tensortypealias[1000];
-  enum tensortype { UNKNOWN, SCALAR, VECTOR, TENSOR };
+  enum tensortype { UNKNOWN, SCALAR, VECTOR, TENSOR, WEYLSCALARS_REAL };
   enum tensortype ttype;
   CCTK_INT tensorparity;
   int tcomponent;
@@ -263,9 +263,8 @@ BndReflectVI (cGH const * restrict const cctkGH,
     /* scalar */
     ttype = SCALAR;
     tcomponent = 0;
-  }
-  else if (CCTK_EQUALS (tensortypealias, "u")
-           || CCTK_EQUALS (tensortypealias, "d"))
+  } else if (CCTK_EQUALS (tensortypealias, "u")
+             || CCTK_EQUALS (tensortypealias, "d"))
   {
     /* vector */
     assert (numvars == 3);
@@ -277,6 +276,13 @@ BndReflectVI (cGH const * restrict const cctkGH,
     assert (numvars == 6);
     ttype = TENSOR;
     tcomponent = vi - firstvar;
+  } else if (CCTK_EQUALS (tensortypealias, "weylscalars_real")) {
+    /* Weyl scalars, stored as 10 real values.  NOTE: This assumes
+       that Psi_0 comes first, which is NOT the default with
+       PsiKadelia.  */
+    assert (numvars == 10);
+    ttype = WEYLSCALARS_REAL;
+    tcomponent = vi - firstvar;
   } else {
     char * groupname = CCTK_GroupName(gi);
     assert (groupname);
@@ -297,6 +303,9 @@ BndReflectVI (cGH const * restrict const cctkGH,
   case TENSOR:
     assert (tcomponent>=0 && tcomponent<6);
     break;
+  case WEYLSCALARS_REAL:
+    assert (tcomponent>=0 && tcomponent<10);
+    break;
   default:
     assert (0);
   }
@@ -366,6 +375,22 @@ BndReflectVI (cGH const * restrict const cctkGH,
             default: assert (0);
             }
             break;
+          case WEYLSCALARS_REAL:
+            {
+              static int const weylparities[10][3] =
+                {{+1,+1,+1},
+                 {-1,-1,-1},
+                 {+1,+1,-1},
+                 {-1,-1,+1},
+                 {+1,+1,+1},
+                 {-1,-1,-1},
+                 {+1,+1,-1},
+                 {-1,-1,+1},
+                 {+1,+1,+1},
+                 {-1,-1,-1}};
+              parity *= weylparities[tcomponent][dir];
+            }
+            break;
           default:
             assert (0);
           }
@@ -563,6 +588,126 @@ BndReflectVI (cGH const * restrict const cctkGH,
 
 
 
+/* When CoordBase is used to specify the location of the boundary
+   points, then ensure that the CoordBase parameters and this thorn's
+   parameters are consistent.  */
+static
+void
+CheckBoundaryParameters (cGH const * restrict const cctkGH,
+                         int const vi,
+                         int const * restrict const stencil)
+{
+  DECLARE_CCTK_PARAMETERS;
+  
+  static int did_check = 0;
+  
+  int type;                     /* Parameter type */
+  void const * ptr;             /* Pointer to parameter value */
+  char const * coordtype;       /* CartGrid3D::type */
+  
+  int dim;                      /* Number of dimensions of vi */
+  
+  CCTK_INT * restrict nboundaryzones; /* CoordBase boundary location */
+  CCTK_INT * restrict is_internal;
+  CCTK_INT * restrict is_staggered;
+  CCTK_INT * restrict shiftout;
+  
+  int do_reflection[6];         /* This thorn's parameters */
+  int do_stagger[6];
+  
+  int d;
+  int ierr;
+  
+  
+  
+  /* Check only once to save time */
+  if (did_check) return;
+  
+  /* Check only for grid functions */
+  if (CCTK_GroupTypeFromVarI (vi) != CCTK_GF) return;
+  
+  did_check = 1;
+  
+  /* Check whether CartGrid3D is active */
+  if (! CCTK_IsThornActive ("CartGrid3D")) return;
+  
+  /* Check whether CoordBase is used */
+  ptr = CCTK_ParameterGet ("type", "CartGrid3D", & type);
+  assert (ptr != 0);
+  assert (type == PARAMETER_KEYWORD);
+  coordtype = * (char const * *) ptr;
+  if (! CCTK_EQUALS (coordtype, "coordbase")) return;
+  
+  /* Get the boundary specification */
+  dim = CCTK_GroupDimFromVarI (vi);
+  assert (dim >= 0);
+  nboundaryzones = malloc (2*dim * sizeof *nboundaryzones);
+  is_internal = malloc (2*dim * sizeof *is_internal);
+  is_staggered = malloc (2*dim * sizeof *is_staggered);
+  shiftout = malloc (2*dim * sizeof *shiftout);
+  ierr = GetBoundarySpecification
+    (2*dim, nboundaryzones, is_internal, is_staggered, shiftout);
+  assert (! ierr);
+  
+  /* Reflection symmetry information */
+  assert (dim == 3);
+  do_reflection[0] = reflection_x;
+  do_reflection[1] = reflection_upper_x;
+  do_reflection[2] = reflection_y;
+  do_reflection[3] = reflection_upper_y;
+  do_reflection[4] = reflection_z;
+  do_reflection[5] = reflection_upper_z;
+  
+  do_stagger[0] = avoid_origin_x;
+  do_stagger[1] = avoid_origin_upper_x;
+  do_stagger[2] = avoid_origin_y;
+  do_stagger[3] = avoid_origin_upper_y;
+  do_stagger[4] = avoid_origin_z;
+  do_stagger[5] = avoid_origin_upper_z;
+  
+  /* Check the boundary sizes */
+  for (d=0; d<6; ++d) {
+    if (do_reflection[d]) {
+      char const dir = "xyz"[d/2];
+      char const * const face = (d%2==0) ? "lower" : "upper";
+      if (stencil[d/2] != nboundaryzones[d]) {
+        CCTK_VWarn (CCTK_WARN_ABORT,
+                    __LINE__, __FILE__, CCTK_THORNSTRING,
+                    "The %s %c face is a symmetry boundary.  Since there are %d ghost zones in the %c direction, the corresponding CoordBase boundary width must also be %d.  The boundary width is currently %d.",
+                    face, dir,
+                    stencil[d/2], dir, stencil[d/2],
+                    (int) nboundaryzones[d]);
+      }
+      if (is_internal[d]) {
+        CCTK_VWarn (CCTK_WARN_ABORT,
+                    __LINE__, __FILE__, CCTK_THORNSTRING,
+                    "The %s %c face is a symmetry boundary.  The corresponding CoordBase boundary must not be internal.",
+                    face, dir);
+      }
+      if (do_stagger[d] != is_staggered[d]) {
+        CCTK_VWarn (CCTK_WARN_ABORT,
+                    __LINE__, __FILE__, CCTK_THORNSTRING,
+                    "The %s %c face is a symmetry boundary.  The symmetry condition and the corresponding CoordBase boundary must either be both staggered or both not staggered.",
+                    face, dir);
+      }
+      if ((do_stagger[d] ? 0 : 1) != shiftout[d]) {
+        CCTK_VWarn (CCTK_WARN_ABORT,
+                    __LINE__, __FILE__, CCTK_THORNSTRING,
+                    "The %s %c face is a symmetry boundary.  If the symmetry condition is staggered, then the corresponding CoordBase shiftout must be 0; otherwise it must be 1.",
+                    face, dir);
+      }
+    }
+  }
+  
+  /* Free memory */
+  free (nboundaryzones);
+  free (is_internal);
+  free (is_staggered);
+  free (shiftout);
+}
+
+
+
 void
 ReflectionSymmetry_Apply (CCTK_ARGUMENTS)
 {
@@ -655,6 +800,8 @@ ReflectionSymmetry_Apply (CCTK_ARGUMENTS)
       CCTK_WARN (0, "Internal error in CCTK_GroupnghostzonesVI");
     }
     
+    CheckBoundaryParameters (cctkGH, vi, stencil);
+    
     ierr = BndReflectVI (cctkGH, stencil, vi);
     if (ierr)
     {
diff --git a/src/interpolate.c b/src/interpolate.c
index 884abe6..7c8831e 100644
--- a/src/interpolate.c
+++ b/src/interpolate.c
@@ -178,7 +178,7 @@ ReflectionSymmetry_Interpolate (CCTK_POINTER_TO_CONST restrict const cctkGH_,
     
     int table;
     char tensortypealias[1000];
-    enum tensortype { UNKNOWN, SCALAR, VECTOR, TENSOR };
+    enum tensortype { UNKNOWN, SCALAR, VECTOR, TENSOR, WEYLSCALARS_REAL };
     enum tensortype ttype;
     CCTK_INT tensorparity;
     int tcomponent;
@@ -224,8 +224,7 @@ ReflectionSymmetry_Interpolate (CCTK_POINTER_TO_CONST restrict const cctkGH_,
     
     ttype = UNKNOWN;
     tcomponent = 0;
-    if (CCTK_EQUALS (tensortypealias, "scalar"))
-    {
+    if (CCTK_EQUALS (tensortypealias, "scalar")) {
       /* scalar */
       if (numvars != 1) {
         groupname = CCTK_GroupName(gi);
@@ -239,8 +238,7 @@ ReflectionSymmetry_Interpolate (CCTK_POINTER_TO_CONST restrict const cctkGH_,
       tcomponent = 0;
     }
     else if (CCTK_EQUALS (tensortypealias, "u")
-             || CCTK_EQUALS (tensortypealias, "d"))
-    {
+             || CCTK_EQUALS (tensortypealias, "d")) {
       /* vector */
       assert (numvars == 3);
       ttype = VECTOR;
@@ -251,6 +249,11 @@ ReflectionSymmetry_Interpolate (CCTK_POINTER_TO_CONST restrict const cctkGH_,
       assert (numvars == 6);
       ttype = TENSOR;
       tcomponent = vi - firstvar;
+    } else if (CCTK_EQUALS (tensortypealias, "weylscalars_real")) {
+      /* Weyl scalars, stored as 10 real numbers */
+      assert (numvars == 10);
+      ttype = WEYLSCALARS_REAL;
+      tcomponent = vi - firstvar;
     } else {
       char * groupname = CCTK_GroupName(gi);
       assert (groupname);
@@ -271,6 +274,9 @@ ReflectionSymmetry_Interpolate (CCTK_POINTER_TO_CONST restrict const cctkGH_,
     case TENSOR:
       assert (tcomponent>=0 && tcomponent<6);
       break;
+    case WEYLSCALARS_REAL:
+      assert (tcomponent>=0 && tcomponent<10);
+      break;
     default:
       assert (0);
     }
@@ -311,6 +317,24 @@ ReflectionSymmetry_Interpolate (CCTK_POINTER_TO_CONST restrict const cctkGH_,
       default: assert (0);
       }
       break;
+    case WEYLSCALARS_REAL:
+      {
+        static int const weylparities[10][3] =
+          {{+1,+1,+1},
+           {-1,-1,-1},
+           {+1,+1,-1},
+           {-1,-1,+1},
+           {+1,+1,+1},
+           {-1,-1,-1},
+           {+1,+1,-1},
+           {-1,-1,+1},
+           {+1,+1,+1},
+           {-1,-1,-1}};
+        for (dir=0; dir<3; ++dir) {
+          parities[dir] = weylparities[tcomponent][dir];
+        }
+      }
+      break;
     default:
       assert (0);
     }