Port Erik's changes from RotatingSymmetry90 over to Cartoon. This

ensures that if the interpolator is taking time derivatives that the
symmetry boundary conditions are suitably aware of what's going on.


git-svn-id: http://svn.cactuscode.org/arrangements/CactusNumerical/Cartoon2D/trunk@98 eec4d7dc-71c2-46d6-addf-10296150bf52

1 files changed, 302 insertions, 69 deletions

diff --git a/src/SymInterp.c b/src/SymInterp.c
index 24668ab..44a3ddd 100644
--- a/src/SymInterp.c
+++ b/src/SymInterp.c
@@ -31,6 +31,22 @@
 
 
 
+
+/* 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
 {
@@ -225,6 +241,35 @@ 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
@@ -252,6 +297,9 @@ CheckTensorType (struct tensor const * restrict const atensor)
 void
 Cartoon2D_CheckTensorTypes (CCTK_ARGUMENTS)
 {
+  int gi;
+  
+  /* Check all internal tensor type definitions */
   CheckTensorType (&scalar);
   CheckTensorType (&vector);
   CheckTensorType (&tensor);
@@ -259,6 +307,92 @@ Cartoon2D_CheckTensorTypes (CCTK_ARGUMENTS)
   CheckTensorType (&symmtensor3a);
   CheckTensorType (&symmtensor3b);
   CheckTensorType (&symmtensor4);
+  
+  /* Check tensor types of all groups */
+  for (gi=0; gi<CCTK_NumGroups(); ++gi) {
+    
+    char tensortypealias[100];
+    int numvars, firstvar;
+    int table;
+    int ierr;
+    
+    numvars = CCTK_NumVarsInGroupI(gi);
+    if (numvars == 0) continue;
+    assert (numvars>0);
+    firstvar = CCTK_FirstVarIndexI(gi);
+    assert (firstvar>=0);
+    table = CCTK_GroupTagsTableI(gi);
+    assert (table>=0);
+    
+    ierr = Util_TableGetString
+      (table, sizeof tensortypealias, tensortypealias, "tensortypealias");
+    if (ierr == UTIL_ERROR_TABLE_NO_SUCH_KEY) {
+      char * groupname = CCTK_GroupName(gi);
+      assert (groupname);
+      CCTK_VWarn (2, __LINE__, __FILE__, CCTK_THORNSTRING,
+                  "Tensor type alias not declared for group \"%s\" -- assuming a scalar",
+                  groupname);
+      free (groupname);
+      strcpy (tensortypealias, "");
+    } else if (ierr<0) {
+      char * groupname = CCTK_GroupName(gi);
+      assert (groupname);
+      CCTK_VWarn (0, __LINE__, __FILE__, CCTK_THORNSTRING,
+                  "Error in tensor type alias declaration for group \"%s\"",
+                  groupname);
+      free (groupname);
+    }
+    
+    if (CCTK_EQUALS (tensortypealias, "")) {
+      /* do nothing */
+    } else if (CCTK_EQUALS (tensortypealias, "scalar")) {
+      /* scalar */
+      if (numvars != 1) {
+        char * groupname = CCTK_GroupName(gi);
+        assert (groupname);
+        CCTK_VWarn (3, __LINE__, __FILE__, CCTK_THORNSTRING,
+                    "Group \"%s\" has the tensor type alias \"scalar\", but contains more than 1 element",
+                    groupname);
+        free (groupname);
+      }
+    } else if (CCTK_EQUALS (tensortypealias, "4scalar")) {
+      /* 4-scalar */
+      if (numvars != 1) {
+        char * groupname = CCTK_GroupName(gi);
+        assert (groupname);
+        CCTK_VWarn (3, __LINE__, __FILE__, CCTK_THORNSTRING,
+                    "Group \"%s\" has the tensor type alias \"4scalar\", but contains more than 1 element",
+                    groupname);
+        free (groupname);
+      }
+    } else if (CCTK_EQUALS (tensortypealias, "u")
+               || CCTK_EQUALS (tensortypealias, "d"))
+    {
+      /* vector */
+      assert (numvars == DIM);
+    } else if (CCTK_EQUALS (tensortypealias, "4u")
+               || CCTK_EQUALS (tensortypealias, "4d"))
+    {
+      /* 4-vector */
+      assert (numvars == DIM+1);
+    } else if (CCTK_EQUALS (tensortypealias, "uu_sym")
+               || CCTK_EQUALS (tensortypealias, "dd_sym")) {
+      /* symmetric tensor */
+      assert (numvars == DIM*(DIM+1)/2);
+    } else if (CCTK_EQUALS (tensortypealias, "4uu_sym")
+               || CCTK_EQUALS (tensortypealias, "4dd_sym")) {
+      /* symmetric 4-tensor */
+      assert (numvars == (DIM+1)*(DIM+2)/2);
+    } else {
+      char * groupname = CCTK_GroupName(gi);
+      assert (groupname);
+      CCTK_VWarn (0, __LINE__, __FILE__, CCTK_THORNSTRING,
+                  "Illegal tensor type alias for group \"%s\"",
+                  groupname);
+      free (groupname);
+    }
+    
+  }
 }
   
 
@@ -282,15 +416,19 @@ Cartoon2D_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];
-  int * restrict * restrict thevars[3];
+  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;
   
@@ -299,17 +437,19 @@ Cartoon2D_SymmetryInterpolate (CCTK_POINTER_TO_CONST const cctkGH_,
   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 */
   
   int iret;                     /* interpolator return value */
   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);
@@ -334,7 +474,7 @@ Cartoon2D_SymmetryInterpolate (CCTK_POINTER_TO_CONST const cctkGH_,
   
   
   /* Copy coordinates */
-  for (d=0; d<3; ++d) {
+  for (d=0; d<DIM; ++d) {
     new_interp_coords[d] = malloc (N_interp_points * sizeof(CCTK_REAL));
     assert (new_interp_coords[d]);
   }
@@ -360,13 +500,6 @@ Cartoon2D_SymmetryInterpolate (CCTK_POINTER_TO_CONST const cctkGH_,
   for (m=0; m<N_output_arrays; ++m) {
     new_output_arrays[m] = malloc (N_interp_points * sizeof(CCTK_REAL));
     assert (new_output_arrays[m]);
-    {
-      CCTK_INT mm;
-      for (mm = 0; mm < N_interp_points; ++mm)
-      {
-        ((CCTK_REAL*)new_output_arrays[m])[mm] = 0;
-      }
-    }
   }
   
   
@@ -383,13 +516,27 @@ Cartoon2D_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]);
   }
   
   
   
   /* Find output variable indices */
+  input_array_time_levels
+    = malloc (N_input_arrays * sizeof *input_array_time_levels);
+  assert (input_array_time_levels);
+  ierr = Util_TableGetIntArray
+    (param_table_handle, N_input_arrays,
+     input_array_time_levels, "input_array_time_levels");
+  if (ierr == UTIL_ERROR_TABLE_NO_SUCH_KEY) {
+    for (m=0; m<N_input_arrays; ++m) {
+      input_array_time_levels[m] = 0; /* time level is 0 */
+    }
+  } else {
+    assert (ierr == N_input_arrays);
+  }
+  
   operand_indices = malloc (N_output_arrays * sizeof *operand_indices);
   assert (operand_indices);
   ierr = Util_TableGetIntArray
@@ -408,14 +555,34 @@ Cartoon2D_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 */
     }
   } 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);
   assert (output_array_indices);
@@ -429,14 +596,19 @@ Cartoon2D_SymmetryInterpolate (CCTK_POINTER_TO_CONST const cctkGH_,
   
   
   /* Map Cactus variables to tensor objects */
-  for (q=0; q<3; ++q) {
-    thetensor[q] = malloc (CCTK_NumVars() * sizeof *thetensor[q]);
-    assert (thetensor[q]);
-    thevars[q] = malloc (CCTK_NumVars() * sizeof *thevars[q]);
-    assert (thevars[q]);
-    for (n=0; n<CCTK_NumVars(); ++n) {
-      thetensor[q][n] = NULL;
-      thevars[q][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;
+        }
+      }
     }
   }
   
@@ -458,9 +630,11 @@ Cartoon2D_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;
     
     /* Get some variable information */
     vi = output_array_indices[m];
@@ -513,16 +687,57 @@ Cartoon2D_SymmetryInterpolate (CCTK_POINTER_TO_CONST const cctkGH_,
       var = 0;
     } else if ((CCTK_EQUALS (tensortypealias, "u")) ||(CCTK_EQUALS (tensortypealias, "d"))) {
       /* vector */
-      assert (numvars == 3);
+      assert (numvars == DIM);
       tensortype = &vector;
       basevar = firstvar;
       var = vi - firstvar;
+    } else if (CCTK_EQUALS (tensortypealias, "4u")
+               || CCTK_EQUALS (tensortypealias, "4d"))
+    {
+      /* 4-vector */
+      assert (numvars == DIM+1);
+      if (vi == firstvar) {
+        /* temporal component */
+        int const off = 0;
+        tensortype = &scalar;
+        basevar = firstvar + off;
+        var = vi - basevar;
+      } else {
+        /* spatial components */
+        int const off = 1;
+        tensortype = &vector;
+        basevar = firstvar + off;
+        var = vi - basevar;
+      }
     } else if (CCTK_EQUALS (tensortypealias, "dd_sym")) {
       /* symmetric tensor */
-      assert (numvars == 6);
+      assert (numvars == DIM*(DIM+1)/2);
       tensortype = &symmtensor;
       basevar = firstvar;
       var = vi - firstvar;
+    } else if (CCTK_EQUALS (tensortypealias, "4uu_sym")
+               || CCTK_EQUALS (tensortypealias, "4dd_sym")) {
+      /* symmetric 4-tensor */
+      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+DIM+1) {
+        /* temporal-spatial components */
+        int const off = 1;
+        tensortype = &vector;
+        basevar = firstvar + off;
+        var = vi - basevar;
+      } else {
+        /* spatial-spatial components */
+        int const off = DIM+1;
+        tensortype = &symmtensor;
+        basevar = firstvar + off;
+        var = vi - basevar;
+      }
     } else {
       char * groupname = CCTK_GroupName(gi);
       assert (groupname);
@@ -538,7 +753,7 @@ Cartoon2D_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);
@@ -547,17 +762,18 @@ Cartoon2D_SymmetryInterpolate (CCTK_POINTER_TO_CONST const cctkGH_,
       assert (component == 0);
       oldrank = tensortype->rank;
     }
+        
+    /* 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. 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 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<=MAX_TIME_LEVELS);
     
     if (tensortype == &scalar) {
       switch (num_derivs) {
@@ -587,7 +803,7 @@ Cartoon2D_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) {
@@ -611,20 +827,28 @@ Cartoon2D_SymmetryInterpolate (CCTK_POINTER_TO_CONST const cctkGH_,
     }
     
     /* Create or cross-check the tensor object */
-    if (! thetensor[num_derivs][basevar]) {
-      thetensor[num_derivs][basevar] = tensortype;
-      assert (! thevars[num_derivs][basevar]);
-      thevars[num_derivs][basevar]
-        = malloc (tensortype->nvars * sizeof *thevars[num_derivs][basevar]);
-      assert (thevars[num_derivs][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_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][basevar][i] = -1;
+        thevars[num_time_derivs][num_derivs][time_level][basevar][i] = -1;
       }
     }
-    assert (thetensor[num_derivs][basevar] == tensortype);
-    assert (thevars[num_derivs][basevar]);
-    assert (thevars[num_derivs][basevar][var] == -1);
-    thevars[num_derivs][basevar][var] = m;
+    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.  */
+    /* 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 */
   
@@ -633,30 +857,33 @@ Cartoon2D_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<=MAX_TIME_LEVELS);
+        
     /* Get the tensor type */
     basevar = thebase[m];
     assert (basevar>=0 && basevar<=CCTK_NumVars());
-    tensortype = thetensor[num_derivs][basevar];
+    tensortype = thetensor[num_time_derivs][num_derivs][time_level][basevar];
     assert (tensortype);
-    vars = thevars[num_derivs][basevar];
+    vars = thevars[num_time_derivs][num_derivs][time_level][basevar];
     assert (vars);
     var = thevar[m];
     assert (var>=0 && var<tensortype->nvars);
@@ -664,7 +891,7 @@ Cartoon2D_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);
@@ -702,11 +929,11 @@ Cartoon2D_SymmetryInterpolate (CCTK_POINTER_TO_CONST const cctkGH_,
       
       int thecomponent;
       for (thecomponent=0; thecomponent<tensortype->ncomps; ++thecomponent) {
-        int theindices[10];
+        int theindices[MAX_RANK+1];
         int thevar;
         int theparity;
         int mm;
-        assert (tensortype->rank<10);
+        assert (tensortype->rank<MAX_RANK+1);
         {
           int r;
           int component = thecomponent;
@@ -748,13 +975,19 @@ Cartoon2D_SymmetryInterpolate (CCTK_POINTER_TO_CONST const cctkGH_,
   }
   free (new_output_arrays);
   
+  free (input_array_time_levels);
   free (operand_indices);
   free (operation_codes);
   free (output_array_indices);
   
-  for (q=0; q<3; ++q) {
-    free ((void *) (thetensor[q]));
-    free ((void *) (thevars[q]));
+  
+  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);
   free (thevar);