Clean up the code:

Replace several magic numbers by named constants DIM, MAX_RANK,
MAX_TIME_LEVELS, and MAX_SPATIAL_DERIV_ORDER

Handle the case where time derivatives are calculated while
interpolating: Introduce a new constant MAX_TIME_DERIV_ORDER.  Add a
new rank to the variables thetensor and thevars.  Take the option
table key time_deriv_order into account.


git-svn-id: http://svn.cactuscode.org/arrangements/CactusNumerical/RotatingSymmetry90/trunk@34 c3c03602-0f4f-0410-b3fa-d2c81c8a7dc5

1 files changed, 149 insertions, 78 deletions

diff --git a/src/interpolate.c b/src/interpolate.c
index be75f3e..2bdb183 100644
--- a/src/interpolate.c
+++ b/src/interpolate.c
@@ -1,7 +1,6 @@
 /* $Header$ */
 
 #include <assert.h>
-#include <math.h>
 #include <stdio.h>
 #include <stdlib.h>
 #include <string.h>
@@ -14,6 +13,21 @@
 
 
 
+/* This is pretty hard coded into all the tensor types and cannot be
+   changed easily.  */
+#define DIM 3 /* spatial dimension */
+
+
+
+/* These can be increased if necessary, but not all necessary tensor
+   types may be defined .  */
+#define MAX_RANK                9 /* maximum tensor rank */
+#define MAX_TIME_LEVELS         3 /* maximum number of time levels */
+#define MAX_SPATIAL_DERIV_ORDER 2 /* maximum spatial derivative order */
+#define MAX_TIME_DERIV_ORDER    2 /* maximum time derivative order */
+
+
+
 /* A tensor description */
 struct tensor
 {
@@ -209,6 +223,36 @@ static struct tensor symmtensor4 = {
 
 
 
+static int
+ipow (int base, int expo)
+{
+  int res;
+  assert (expo >= 0);
+  res = expo & 1 ? base : 1;
+  while (expo >>= 1)
+  {
+    base *= base;
+    if (expo & 1) res *= base;
+  }
+  return res;
+}
+
+
+
+static int
+ilog (int res, int const base)
+{
+  int expo;
+  assert (base > 0);
+  assert (res >= 0);
+  for (expo = 0; res > 0; ++ expo) {
+    res /= base;
+  }
+  return expo;
+}
+
+
+
 /* Ensure that all tensor declarations are internally consistent */
 static void
 CheckTensorType (struct tensor const * restrict const atensor)
@@ -218,7 +262,7 @@ CheckTensorType (struct tensor const * restrict const atensor)
   assert (atensor->dim>=0);
   assert (atensor->rank>=0);
   assert (atensor->ncomps>=0);
-  assert (atensor->ncomps == floor(pow(atensor->dim, atensor->rank) + 0.5));
+  assert (atensor->ncomps == ipow (atensor->dim, atensor->rank));
   assert (atensor->nvars>=0 && atensor->nvars<=atensor->ncomps);
   assert (atensor->vars);
   for (i=0; i<atensor->ncomps; ++i) {
@@ -307,20 +351,20 @@ Rot90_CheckTensorTypes (CCTK_ARGUMENTS)
                || CCTK_EQUALS (tensortypealias, "d"))
     {
       /* vector */
-      assert (numvars == 3);
+      assert (numvars == DIM);
     } else if (CCTK_EQUALS (tensortypealias, "4u")
                || CCTK_EQUALS (tensortypealias, "4d"))
     {
       /* 4-vector */
-      assert (numvars == 4);
+      assert (numvars == DIM+1);
     } else if (CCTK_EQUALS (tensortypealias, "uu_sym")
                || CCTK_EQUALS (tensortypealias, "dd_sym")) {
       /* symmetric tensor */
-      assert (numvars == 6);
+      assert (numvars == DIM*(DIM+1)/2);
     } else if (CCTK_EQUALS (tensortypealias, "4uu_sym")
                || CCTK_EQUALS (tensortypealias, "4dd_sym")) {
       /* symmetric 4-tensor */
-      assert (numvars == 10);
+      assert (numvars == (DIM+1)*(DIM+2)/2);
     } else {
       char * groupname = CCTK_GroupName(gi);
       assert (groupname);
@@ -354,25 +398,28 @@ Rot90_SymmetryInterpolate (CCTK_POINTER_TO_CONST const cctkGH_,
 {
   cGH const * restrict const cctkGH = cctkGH_;
   
-  CCTK_POINTER new_interp_coords[3];
+  CCTK_POINTER new_interp_coords[DIM];
   CCTK_POINTER * restrict new_output_arrays;
   CCTK_INT new_faces;
   CCTK_INT * restrict input_array_time_levels;
   CCTK_INT * restrict operand_indices;
   CCTK_INT * restrict operation_codes;
+  CCTK_INT * restrict time_deriv_order;
   CCTK_INT * restrict output_array_indices;
   
-  struct tensor const * restrict * restrict thetensor[3][4];
-  int * restrict * restrict thevars[3][4];
+  struct tensor const * restrict * restrict thetensor
+    [MAX_TIME_DERIV_ORDER+1][MAX_SPATIAL_DERIV_ORDER+1][MAX_TIME_LEVELS+1];
+  int * restrict * restrict thevars
+    [MAX_TIME_DERIV_ORDER+1][MAX_SPATIAL_DERIV_ORDER+1][MAX_TIME_LEVELS+1];
   int * restrict thebase;
   int * restrict thevar;
   
   int m;                        /* output array */
   int n;                        /* point */
   int i;                        /* var */
-  int f;                        /* face */
   int d;                        /* dimension */
-  int q;                        /* number of derivatives */
+  int p;                        /* number of time derivatives */
+  int q;                        /* number of spatial derivatives */
   int tl;                       /* time level */
   int r;                        /* rank */
   
@@ -380,10 +427,10 @@ Rot90_SymmetryInterpolate (CCTK_POINTER_TO_CONST const cctkGH_,
   int ierr;
   
   /* Check arguments */
-  assert (N_dims == 3);
+  assert (N_dims == DIM);
   assert (N_interp_points >= 0);
   assert (interp_coords_type >= 0);
-  for (d=0; d<3; ++d) {
+  for (d=0; d<DIM; ++d) {
     assert (N_interp_points == 0 || interp_coords[d]);
   }
   assert (N_output_arrays >= 0);
@@ -406,8 +453,8 @@ Rot90_SymmetryInterpolate (CCTK_POINTER_TO_CONST const cctkGH_,
   
   
   /* Copy coordinates */
-  for (d=0; d<3; ++d) {
-    new_interp_coords[d] = (CCTK_REAL*)malloc (N_interp_points * sizeof(CCTK_REAL));
+  for (d=0; d<DIM; ++d) {
+    new_interp_coords[d] = malloc (N_interp_points * sizeof(CCTK_REAL));
     assert (new_interp_coords[d]);
   }
   
@@ -457,7 +504,7 @@ Rot90_SymmetryInterpolate (CCTK_POINTER_TO_CONST const cctkGH_,
   
   
   /* Free coordinates */
-  for (d=0; d<3; ++d) {
+  for (d=0; d<DIM; ++d) {
     free (new_interp_coords[d]);
   }
   
@@ -496,13 +543,33 @@ Rot90_SymmetryInterpolate (CCTK_POINTER_TO_CONST const cctkGH_,
   ierr = Util_TableGetIntArray
     (param_table_handle, N_output_arrays, operation_codes, "operation_codes");
   if (ierr == UTIL_ERROR_TABLE_NO_SUCH_KEY) {
-    assert (N_output_arrays == N_input_arrays);
+    assert (N_output_arrays == N_input_arrays); /* why? */
     for (m=0; m<N_output_arrays; ++m) {
-      operation_codes[m] = 0;     /* do not take derivatives */
+      operation_codes[m] = 0;     /* do not take spatial derivatives */
     }
   } else {
     assert (ierr == N_output_arrays);
   }
+  for (m=0; m<N_output_arrays; ++m) {
+    assert (operation_codes[m] >= 0
+            && operation_codes[m] < ipow (10, MAX_SPATIAL_DERIV_ORDER));
+  }
+  
+  time_deriv_order = malloc (N_output_arrays * sizeof *time_deriv_order);
+  assert (time_deriv_order);
+  ierr = Util_TableGetIntArray
+    (param_table_handle, N_output_arrays, time_deriv_order, "time_deriv_order");
+  if (ierr == UTIL_ERROR_TABLE_NO_SUCH_KEY) {
+    for (m=0; m<N_output_arrays; ++m) {
+      time_deriv_order[m] = 0;  /* do not take time derivatives */
+    }
+  } else {
+    assert (ierr == N_output_arrays);
+  }
+  for (m=0; m<N_output_arrays; ++m) {
+    assert (time_deriv_order[m] >= 0
+            && time_deriv_order[m] <= MAX_TIME_DERIV_ORDER);
+  }
   
   output_array_indices
     = malloc (N_output_arrays * sizeof *output_array_indices);
@@ -517,15 +584,18 @@ Rot90_SymmetryInterpolate (CCTK_POINTER_TO_CONST const cctkGH_,
   
   
   /* Map Cactus variables to tensor objects */
-  for (q=0; q<3; ++q) {
-    for (tl=0; tl<4; ++tl) {
-      thetensor[q][tl] = malloc (CCTK_NumVars() * sizeof *thetensor[q][tl]);
-      assert (thetensor[q][tl]);
-      thevars[q][tl] = malloc (CCTK_NumVars() * sizeof *thevars[q][tl]);
-      assert (thevars[q][tl]);
-      for (n=0; n<CCTK_NumVars(); ++n) {
-        thetensor[q][tl][n] = NULL;
-        thevars[q][tl][n] = NULL;
+  for (p=0; p<=MAX_TIME_DERIV_ORDER; ++p) {
+    for (q=0; q<=MAX_SPATIAL_DERIV_ORDER; ++q) {
+      for (tl=0; tl<=MAX_TIME_LEVELS; ++tl) {
+        thetensor[p][q][tl]
+          = malloc (CCTK_NumVars() * sizeof *thetensor[p][q][tl]);
+        assert (thetensor[p][q][tl]);
+        thevars[p][q][tl] = malloc (CCTK_NumVars() * sizeof *thevars[p][q][tl]);
+        assert (thevars[p][q][tl]);
+        for (n=0; n<CCTK_NumVars(); ++n) {
+          thetensor[p][q][tl][n] = NULL;
+          thevars[p][q][tl][n] = NULL;
+        }
       }
     }
   }
@@ -548,8 +618,9 @@ Rot90_SymmetryInterpolate (CCTK_POINTER_TO_CONST const cctkGH_,
     struct tensor const * tensortype;
     int basevar;
     int var;
-    int indices[10];
+    int indices[MAX_RANK+1];
     int oldrank;
+    int num_time_derivs;
     int num_derivs;
     int time_level;
     
@@ -619,7 +690,7 @@ Rot90_SymmetryInterpolate (CCTK_POINTER_TO_CONST const cctkGH_,
                || CCTK_EQUALS (tensortypealias, "d"))
     {
       /* vector */
-      assert (numvars == 3);
+      assert (numvars == DIM);
       tensortype = &vector;
       basevar = firstvar;
       var = vi - basevar;
@@ -627,7 +698,7 @@ Rot90_SymmetryInterpolate (CCTK_POINTER_TO_CONST const cctkGH_,
                || CCTK_EQUALS (tensortypealias, "4d"))
     {
       /* 4-vector */
-      assert (numvars == 4);
+      assert (numvars == DIM+1);
       if (vi == firstvar) {
         /* temporal component */
         int const off = 0;
@@ -644,21 +715,21 @@ Rot90_SymmetryInterpolate (CCTK_POINTER_TO_CONST const cctkGH_,
     } else if (CCTK_EQUALS (tensortypealias, "uu_sym")
                || CCTK_EQUALS (tensortypealias, "dd_sym")) {
       /* symmetric tensor */
-      assert (numvars == 6);
+      assert (numvars == DIM*(DIM+1)/2);
       tensortype = &symmtensor;
       basevar = firstvar;
       var = vi - basevar;
     } else if (CCTK_EQUALS (tensortypealias, "4uu_sym")
                || CCTK_EQUALS (tensortypealias, "4dd_sym")) {
       /* symmetric 4-tensor */
-      assert (numvars == 10);
+      assert (numvars == (DIM+1)*(DIM+2)/2);
       if (vi == firstvar) {
         /* temporal-temporal component */
         int const off = 0;
         tensortype = &scalar;
         basevar = firstvar + off;
         var = vi - basevar;
-      } else if (vi < firstvar+4) {
+      } else if (vi < firstvar+DIM+1) {
         /* temporal-spatial components */
         int const off = 1;
         tensortype = &vector;
@@ -666,7 +737,7 @@ Rot90_SymmetryInterpolate (CCTK_POINTER_TO_CONST const cctkGH_,
         var = vi - basevar;
       } else {
         /* spatial-spatial components */
-        int const off = 4;
+        int const off = DIM+1;
         tensortype = &symmtensor;
         basevar = firstvar + off;
         var = vi - basevar;
@@ -686,7 +757,7 @@ Rot90_SymmetryInterpolate (CCTK_POINTER_TO_CONST const cctkGH_,
       int component;
       assert (var>=0 && var<tensortype->nvars);
       component = tensortype->comps[var];
-      assert (tensortype->rank<10);
+      assert (tensortype->rank>=0 && tensortype->rank<=MAX_RANK);
       for (r=tensortype->rank-1; r>=0; --r) {
         indices[r] = component % tensortype->dim;
         assert (indices[r]>=0 && indices[r]<tensortype->dim);
@@ -696,20 +767,17 @@ Rot90_SymmetryInterpolate (CCTK_POINTER_TO_CONST const cctkGH_,
       oldrank = tensortype->rank;
     }
     
-    /* Take operation code (i.e., derivatives) into account */
-    {
-      int code = operation_codes[m];
-      num_derivs
-        = code == 0 ? 0
-        : code < 10 ? 1
-        : code < 100 ? 2
-        : -1;
-      assert (num_derivs>=0);
-    }
+    /* Take time derivative order (i.e., time derivatives) into account */
+    num_time_derivs = time_deriv_order[m];
+    assert (num_time_derivs>=0 && num_time_derivs<=MAX_TIME_DERIV_ORDER);
+    
+    /* Take operation code (i.e., spatial derivatives) into account */
+    num_derivs = ilog (operation_codes[m], 10);
+    assert (num_derivs>=0 && num_derivs<=MAX_SPATIAL_DERIV_ORDER);
     
     /* Get the time level */
     time_level = input_array_time_levels[operand_indices[m]];
-    assert (time_level>=0 && time_level<4);
+    assert (time_level>=0 && time_level<=MAX_TIME_LEVELS);
     
     if (tensortype == &scalar) {
       switch (num_derivs) {
@@ -739,7 +807,7 @@ Rot90_SymmetryInterpolate (CCTK_POINTER_TO_CONST const cctkGH_,
     /* Find the additional indices */
     {
       int code;
-      assert (tensortype->rank<10);
+      assert (tensortype->rank>=0 && tensortype->rank<=MAX_RANK);
       assert (tensortype->rank >= oldrank);
       code = operation_codes[m];
       for (r=oldrank; r<tensortype->rank; ++r) {
@@ -763,25 +831,28 @@ Rot90_SymmetryInterpolate (CCTK_POINTER_TO_CONST const cctkGH_,
     }
     
     /* Create or cross-check the tensor object */
-    if (! thetensor[num_derivs][time_level][basevar]) {
-      thetensor[num_derivs][time_level][basevar] = tensortype;
-      assert (! thevars[num_derivs][time_level][basevar]);
-      thevars[num_derivs][time_level][basevar]
+    if (! thetensor[num_time_derivs][num_derivs][time_level][basevar]) {
+      thetensor[num_time_derivs][num_derivs][time_level][basevar] = tensortype;
+      assert (! thevars[num_time_derivs][num_derivs][time_level][basevar]);
+      thevars[num_time_derivs][num_derivs][time_level][basevar]
         = malloc (tensortype->nvars
-                  * sizeof *thevars[num_derivs][time_level][basevar]);
-      assert (thevars[num_derivs][time_level][basevar]);
+                  * sizeof *thevars[num_time_derivs][num_derivs][time_level][basevar]);
+      assert (thevars[num_time_derivs][num_derivs][time_level][basevar]);
       for (i=0; i<tensortype->nvars; ++i) {
-        thevars[num_derivs][time_level][basevar][i] = -1;
+        thevars[num_time_derivs][num_derivs][time_level][basevar][i] = -1;
       }
     }
-    assert (thetensor[num_derivs][time_level][basevar] == tensortype);
-    assert (thevars[num_derivs][time_level][basevar]);
+    assert (thetensor[num_time_derivs][num_derivs][time_level][basevar]
+            == tensortype);
+    assert (thevars[num_time_derivs][num_derivs][time_level][basevar]);
     /* This does not hold if the caller requests the same
        interpolation to be done into different output arrays.  This
        may happen e.g. when CarpetInterp needs to differentiate in
        time.  This is arguably a performance bug in CarpetInterp.  */
-    /* assert (thevars[num_derivs][time_level][basevar][var] == -1); */
-    thevars[num_derivs][time_level][basevar][var] = m;
+    /* See whether this goes away now.  */
+    assert (thevars[num_time_derivs][num_derivs][time_level][basevar][var]
+            == -1);
+    thevars[num_time_derivs][num_derivs][time_level][basevar][var] = m;
     
   } /* for m */
   
@@ -790,35 +861,33 @@ Rot90_SymmetryInterpolate (CCTK_POINTER_TO_CONST const cctkGH_,
   /* Loop over all output arrays */
   for (m=0; m<N_output_arrays; ++m) {
     
+    int num_time_derivs;
     int num_derivs;
     int time_level;
     int basevar;
     struct tensor const * restrict tensortype;
     int const * restrict vars;
     int var;
-    int indices[10];
+    int indices[MAX_RANK+1];
     
-    /* Take operation code (i.e. derivatives) into account */
-    {
-      int code = operation_codes[m];
-      num_derivs
-        = code == 0 ? 0
-        : code < 10 ? 1
-        : code < 100 ? 2
-        : -1;
-      assert (num_derivs>=0);
-    }
+    /* Take time derivative order (i.e., time derivatives) into account */
+    num_time_derivs = time_deriv_order[m];
+    assert (num_time_derivs>=0 && num_time_derivs<=MAX_TIME_DERIV_ORDER);
+    
+    /* Take operation code (i.e., spatial derivatives) into account */
+    num_derivs = ilog (operation_codes[m], 10);
+    assert (num_derivs>=0 && num_derivs<=MAX_SPATIAL_DERIV_ORDER);
     
     /* Get the time level */
     time_level = input_array_time_levels[operand_indices[m]];
-    assert (time_level>=0 && time_level<4);
+    assert (time_level>=0 && time_level<=MAX_TIME_LEVELS);
     
     /* Get the tensor type */
     basevar = thebase[m];
     assert (basevar>=0 && basevar<=CCTK_NumVars());
-    tensortype = thetensor[num_derivs][time_level][basevar];
+    tensortype = thetensor[num_time_derivs][num_derivs][time_level][basevar];
     assert (tensortype);
-    vars = thevars[num_derivs][time_level][basevar];
+    vars = thevars[num_time_derivs][num_derivs][time_level][basevar];
     assert (vars);
     var = thevar[m];
     assert (var>=0 && var<tensortype->nvars);
@@ -826,7 +895,7 @@ Rot90_SymmetryInterpolate (CCTK_POINTER_TO_CONST const cctkGH_,
     /* Transform into indices */
     {
       int component = tensortype->comps[var];
-      assert (tensortype->rank<10);
+      assert (tensortype->rank>=0 && tensortype->rank<=MAX_RANK);
       for (r=tensortype->rank-1; r>=0; --r) {
         indices[r] = component % tensortype->dim;
         assert (indices[r]>=0 && indices[r]<tensortype->dim);
@@ -840,7 +909,7 @@ Rot90_SymmetryInterpolate (CCTK_POINTER_TO_CONST const cctkGH_,
       CCTK_REAL x = ((CCTK_REAL const *)interp_coords[0])[n];
       CCTK_REAL y = ((CCTK_REAL const *)interp_coords[1])[n];
       
-      int myindices[10];
+      int myindices[MAX_RANK+1];
       int parity;
       
       for (r=0; r<tensortype->rank; ++r) {
@@ -916,10 +985,12 @@ Rot90_SymmetryInterpolate (CCTK_POINTER_TO_CONST const cctkGH_,
   free (operation_codes);
   free (output_array_indices);
   
-  for (q=0; q<3; ++q) {
-    for (tl=0; tl<4; ++tl) {
-      free ((void *) (thetensor[q][tl]));
-      free ((void *) (thevars[q][tl]));
+  for (p=0; p<=MAX_TIME_DERIV_ORDER; ++p) {
+    for (q=0; q<=MAX_SPATIAL_DERIV_ORDER; ++q) {
+      for (tl=0; tl<=MAX_TIME_LEVELS; ++tl) {
+        free ((void *) (thetensor[p][q][tl]));
+        free ((void *) (thevars[p][q][tl]));
+      }
     }
   }
   free (thebase);