Replace Cactus complex number type with C/C++ complex numbers

Map CCTK_COMPLEX to "double complex" in C, and "complex<double>" in
C++. (It is already mapped to "double complex" in Fortran.)

Update type definitions.

Re-implement Cactus complex number math functions by calling the
respective C functions.

Update thorn that access real and imaginary parts of complex numbers
to use standard-conforming methods instead.


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

1 files changed, 106 insertions, 47 deletions

diff --git a/src/apply.c b/src/apply.c
index f179c0e..0ef2674 100644
--- a/src/apply.c
+++ b/src/apply.c
@@ -56,20 +56,14 @@
           int const jj = jjoff + jdir * j;                              \
           int const kk = kkoff + kdir * k;                              \
           int const srcind = ii + ni * (jj + nj * kk);                  \
-          REAL(dstvar[dstind] RE = parity * srcvar[srcind] RE;)         \
-          IMAG(dstvar[dstind] IM = parity * srcvar[srcind] IM;)         \
+          dstvar[dstind] = parity * srcvar[srcind];                     \
         }                                                               \
       }                                                                 \
     }                                                                   \
 
-#define COPY_POST(VARTYPE)                      \
+#define COPY_POST(VARTYPE)                                              \
   }
 
-#define REAL(x) x
-#define IMAG(x) /* nothing */
-#define RE /* nothing */
-#define IM ERROR ERROR ERROR
-
 #ifdef HAVE_CCTK_INT1
 COPY_PRE(CCTK_INT1)
 #pragma omp parallel for
@@ -119,16 +113,6 @@ COPY_LOOP(CCTK_REAL16)
 COPY_POST(CCTK_REAL16)
 #endif
 
-#undef REAL
-#undef IMAG
-#undef RE
-#undef IM
-
-#define REAL(x) x
-#define IMAG(x) x
-#define RE .Re
-#define IM .Im
-
 #ifdef HAVE_CCTK_COMPLEX8
 COPY_PRE(CCTK_COMPLEX8)
 #pragma omp parallel for
@@ -150,11 +134,6 @@ COPY_LOOP(CCTK_COMPLEX32)
 COPY_POST(CCTK_COMPLEX32)
 #endif
 
-#undef REAL
-#undef IMAG
-#undef RE
-#undef IM
-
 #undef COPY_PRE
 #undef COPY_LOOP
 #undef COPY_POST
@@ -177,21 +156,35 @@ BndReflectVI (cGH const * restrict const cctkGH,
   
   int table;
   char tensortypealias[1000];
-  enum tensortype { UNKNOWN,
+  enum tensortype { TT_UNKNOWN,
                     SCALAR, VECTOR, SYMTENSOR, SYMTENSOR3, TENSOR,
                     WEYLSCALARS_REAL, MANUALCARTESIAN };
   enum tensortype ttype;
   CCTK_INT tensorparity;
   int tcomponent;
   
+  /* The stagger type is defined for grid variable groups, and
+     depending on this type, each variable in the group may be
+     staggered differently. This implies that groups with FACE or EDGE
+     staggering have well-defined tensor types. */
+  /* The stagger type is defined assuming a cell-centered grid, since
+     this seems to be the common case where different stagger types
+     are used. If the grid is vertex centered, then we currently abort
+     when staggered grid variables are encountered. (We determine the
+     centering of the grid via the avoid_origin_* parameters.) */
+  char staggertype[1000];
+  enum staggertype { ST_UNKNOWN, CELL, FACE, EDGE, VERTEX };
+  enum staggertype stype;
+  
   int do_reflection[6];
-  int do_stagger[6];
+  int do_stagger_grid[6];
   
   int dir, face;
   
-  int lsh[3], imin[3], imax[3], ioff[3], idir[3];
+  int ash[3], imin[3], imax[3], ioff[3], idir[3];
   
   int parity;
+  int do_stagger;
   int manual_parities[3];
 
   int d;
@@ -287,7 +280,7 @@ BndReflectVI (cGH const * restrict const cctkGH,
     free (groupname);
   }
   
-  ttype = UNKNOWN;
+  ttype = TT_UNKNOWN;
   tcomponent = 0;
   if (CCTK_EQUALS (tensortypealias, "scalar")) {
     /* scalar */
@@ -414,6 +407,45 @@ BndReflectVI (cGH const * restrict const cctkGH,
   
   
   
+  /* Get and check stagger type information */
+  ierr = Util_TableGetString
+    (table, sizeof staggertype, staggertype, "staggertype");
+  if (ierr == UTIL_ERROR_TABLE_NO_SUCH_KEY) {
+    /* assume cell centering */
+    strcpy (staggertype, "cell");
+  } else if (ierr<0) {
+    char * groupname = CCTK_GroupName(gi);
+    assert (groupname);
+    CCTK_VWarn (0, __LINE__, __FILE__, CCTK_THORNSTRING,
+                "Error in stagger type declaration for group \"%s\": %d",
+                groupname, ierr);
+    free (groupname);
+  }
+  
+  stype = ST_UNKNOWN;
+  if (CCTK_EQUALS (staggertype, "cell")) {
+    /* cell */
+    stype = CELL;
+  } else if (CCTK_EQUALS (staggertype, "face")) {
+    /* face */
+    stype = FACE;
+  } else if (CCTK_EQUALS (staggertype, "edge")) {
+    /* edge */
+    stype = EDGE;
+  } else if (CCTK_EQUALS (staggertype, "vertex")) {
+    /* vertex */
+    stype = VERTEX;
+  } else {
+    char * groupname = CCTK_GroupName(gi);
+    assert (groupname);
+    CCTK_VWarn (0, __LINE__, __FILE__, CCTK_THORNSTRING,
+                "Illegal stagger type \"%s\" for group \"%s\"",
+                staggertype, groupname);
+    free (groupname);
+  }
+  
+  
+  
   /* Reflection symmetry information */
   do_reflection[0] = reflection_x;
   do_reflection[1] = reflection_upper_x;
@@ -422,12 +454,12 @@ BndReflectVI (cGH const * restrict const cctkGH,
   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;
+  do_stagger_grid[0] = avoid_origin_x;
+  do_stagger_grid[1] = avoid_origin_upper_x;
+  do_stagger_grid[2] = avoid_origin_y;
+  do_stagger_grid[3] = avoid_origin_upper_y;
+  do_stagger_grid[4] = avoid_origin_z;
+  do_stagger_grid[5] = avoid_origin_upper_z;
   
   
   
@@ -512,9 +544,36 @@ BndReflectVI (cGH const * restrict const cctkGH,
             assert (0);
           }
           
+          /* Find staggering */
+          do_stagger = do_stagger_grid[2*dir+face];
+          switch (stype) {
+          case CELL:
+            /* do nothing */
+            break;
+          case FACE:
+            assert (face == 0);  /* assume lower face */
+            assert (do_stagger); /* assume cell-centered grid */
+            assert (ttype == VECTOR); /* TODO: support other tensor types */
+            if (dir == tcomponent) do_stagger = !do_stagger;
+            break;
+          case EDGE:
+            assert (face == 0);  /* assume lower face */
+            assert (do_stagger); /* assume cell-centered grid */
+            assert (ttype == VECTOR); /* TODO: support other tensor types */
+            if (dir != tcomponent) do_stagger = !do_stagger;
+            break;
+          case VERTEX:
+            assert (face == 0);  /* assume lower face */
+            assert (do_stagger); /* assume cell-centered grid */
+            do_stagger = !do_stagger;
+            break;
+          default:
+            assert (0);
+          }
+          
           /* Find region extent */
           for (d=0; d<3; ++d) {
-            lsh[d] = cctkGH->cctk_lsh[d];
+            ash[d] = cctkGH->cctk_ash[d];
             imin[d] = 0;
             imax[d] = cctkGH->cctk_lsh[d];
             ioff[d] = 0;
@@ -539,11 +598,11 @@ BndReflectVI (cGH const * restrict const cctkGH,
           if (face == 0) {
             imax[dir] = cctkGH->cctk_nghostzones[dir];
             ioff[dir] = (+ 2*cctkGH->cctk_nghostzones[dir] - 1 
-                         + (do_stagger[2*dir+face] ? 0 : 1));
+                         + (do_stagger ? 0 : 1));
             idir[dir] = -1;
           } else {
             imin[dir] = cctkGH->cctk_lsh[dir] - cctkGH->cctk_nghostzones[dir];
-            ioff[dir] = - 1 - (do_stagger[2*dir+face] ? 0 : 1); 
+            ioff[dir] = - 1 - (do_stagger ? 0 : 1); 
             idir[dir] = -1;
           }
           
@@ -554,7 +613,7 @@ BndReflectVI (cGH const * restrict const cctkGH,
             int const have_points = cctkGH->cctk_gsh[dir];
             int const need_points =
               3 * cctkGH->cctk_nghostzones[dir]
-              + !do_stagger[2*dir] + !do_stagger[2*dir+1];
+              + !do_stagger_grid[2*dir] + !do_stagger_grid[2*dir+1];
             if (need_points > have_points) {
               CCTK_VWarn (CCTK_WARN_ABORT, __LINE__, __FILE__, CCTK_THORNSTRING,
                           "Cannot apply symmetry boundary zones in the %s %c direction, since there seem to be more symmetry zones than interior zones",
@@ -567,7 +626,7 @@ BndReflectVI (cGH const * restrict const cctkGH,
           switch (group.vartype) {
             
 #define ARGS  varptr, varptr,                   \
-              lsh[0], lsh[1], lsh[2],           \
+              ash[0], ash[1], ash[2],           \
               imin[0], imin[1], imin[2],        \
               imax[0], imax[1], imax[2],        \
               ioff[0], ioff[1], ioff[2],        \
@@ -706,7 +765,7 @@ CheckBoundaryParameters (cGH const * restrict const cctkGH,
   CCTK_INT * restrict shiftout;
   
   int do_reflection[6];         /* This thorn's parameters */
-  int do_stagger[6];
+  int do_stagger_grid[6];
   
   int d;
   int ierr;
@@ -751,12 +810,12 @@ CheckBoundaryParameters (cGH const * restrict const cctkGH,
   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;
+  do_stagger_grid[0] = avoid_origin_x;
+  do_stagger_grid[1] = avoid_origin_upper_x;
+  do_stagger_grid[2] = avoid_origin_y;
+  do_stagger_grid[3] = avoid_origin_upper_y;
+  do_stagger_grid[4] = avoid_origin_z;
+  do_stagger_grid[5] = avoid_origin_upper_z;
   
   /* Check the boundary sizes */
   for (d=0; d<6; ++d) {
@@ -777,13 +836,13 @@ CheckBoundaryParameters (cGH const * restrict const cctkGH,
                     "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]) {
+      if (do_stagger_grid[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]) {
+      if ((do_stagger_grid[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.",