Zeroth attempt at doing excision with Cartoon. Now requires the

Spacemask thorn from CactusEinstein.

Initial tests suggest perfect second order convergence but not quite
perfect fourth order convergence, at least at reasonable
resolutions. Convergence is still better than 3rd order, though.

Excision with the new spacemask has been added but not tested.


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

8 files changed, 342 insertions, 16 deletions

diff --git a/configuration.ccl b/configuration.ccl
new file mode 100644
index 0000000..6582098
--- /dev/null
+++ b/configuration.ccl
@@ -0,0 +1,4 @@
+# Configuration definition for thorn Cartoon2D
+# $Header$
+
+REQUIRES SpaceMask
diff --git a/interface.ccl b/interface.ccl
index 6d97591..b432be3 100644
--- a/interface.ccl
+++ b/interface.ccl
@@ -6,6 +6,8 @@ implements: cartoon2d
 INCLUDES HEADER: Cartoon2D_tensors.h in Cartoon2D_tensors.h
 INCLUDES HEADER: Cartoon2D.h in Cartoon2D.h
 
+USES INCLUDE: SpaceMask.h
+
 # Function aliases:
 
 CCTK_INT FUNCTION SymmetryRegister (CCTK_STRING IN sym_name)
@@ -67,3 +69,14 @@ CCTK_INT FUNCTION Boundary_SelectedGVs(CCTK_POINTER_TO_CONST IN GH, \
   CCTK_INT ARRAY OUT faces, CCTK_INT ARRAY OUT boundary_widths, \
   CCTK_INT ARRAY OUT table_handles, CCTK_STRING IN bc_name)
 REQUIRES FUNCTION Boundary_SelectedGVs
+
+private:
+
+int excision_variables type=SCALAR
+{
+  excision_active
+  old_mask_vi
+  new_mask_vi
+  new_excision_field
+  new_excision_descriptor
+} "Internal variables to store information about excision"
diff --git a/param.ccl b/param.ccl
index 3f95f13..76c71a8 100644
--- a/param.ccl
+++ b/param.ccl
@@ -28,3 +28,31 @@ INT eno_order "The interpolation order applied to the ENO interpolator"
 BOOLEAN allow_grid_resize "Allow grid to be resized in a cartoon-compatible way"
 {
 } "no"
+
+BOOLEAN old_excision "Are we doing excision based on the old style mask?"
+{
+} "no"
+
+BOOLEAN new_excision "Are we doing excision based on the new style mask?"
+{
+} "no"
+
+STRING old_style_excision_var "The variable to be checked for old style excision"
+{
+  ".*" :: "Expected to be \'Spacemask::emask\'"
+} ""
+
+STRING new_style_excision_var "The variable to be checked for new style excision"
+{
+  ".*" :: "Expected to be \'Spacemask::space_mask\'"
+} ""
+
+STRING new_mask_field_name "The name of the field that describes excision for the new mask"
+{
+  ".*" :: "Could be anything"
+} ""
+
+STRING new_mask_excised_name "The name of the descriptor that says the point is excised for the new mask"
+{
+  ".*" :: "Could be anything"
+} ""
diff --git a/schedule.ccl b/schedule.ccl
index f0f19a4..55719b8 100644
--- a/schedule.ccl
+++ b/schedule.ccl
@@ -3,12 +3,20 @@
 
 if (cartoon_active)
 {
+  STORAGE: excision_variables
+
   schedule Cartoon2D_CheckTensorTypes at PARAMCHECK
   {
     LANG: C
     OPTIONS: meta
   } "Check tensor type definitions for consistency"
 
+  schedule Cartoon2D_InitExcisionVars at BASEGRID
+  {
+    LANG: C
+    OPTIONS: global
+  } "Initialize the excision variables"
+
   if (allow_grid_resize)
   {
     schedule Cartoon2D_SetGrid at CCTK_RECOVER_PARAMETERS
diff --git a/src/ApplyCartoon.c b/src/ApplyCartoon.c
index 87680cd..b1417d6 100644
--- a/src/ApplyCartoon.c
+++ b/src/ApplyCartoon.c
@@ -178,7 +178,7 @@ void Cartoon_ApplyBoundaries(CCTK_ARGUMENTS)
     {
       if (CCTK_Equals(prolongtype, "None")) 
       {
-        prolongmethod = PROLONG_LAGRANGE; /* But why? */
+        prolongmethod = PROLONG_NONE; /* But why? */
       } 
       else if (CCTK_Equals(prolongtype, "Lagrange")) 
       {
diff --git a/src/Cartoon2DBC.c b/src/Cartoon2DBC.c
index f74c30f..7653424 100644
--- a/src/Cartoon2DBC.c
+++ b/src/Cartoon2DBC.c
@@ -24,6 +24,7 @@
 #include "cctk_FortranString.h"
 
 #include "Cartoon2D_tensors.h"
+#include "SpaceMask.h"
 
 static const char *rcsid = "$Id$";
 
@@ -38,9 +39,11 @@ static CCTK_REAL Cartoon2DInterp_ENO(const cGH *cctkGH,
                                     CCTK_REAL *f, int i, int di, 
                                     int ijk, CCTK_REAL x);
 CCTK_REAL interpolate_local(int order, CCTK_REAL x0, CCTK_REAL dx, 
-			 CCTK_REAL y[], CCTK_REAL x);
+                            CCTK_REAL y[], CCTK_REAL x, 
+                            CCTK_INT excised_points[]);
 CCTK_REAL interpolate_eno(CCTK_INT order, CCTK_REAL x0, CCTK_REAL dx,
-                          CCTK_REAL y[], CCTK_REAL x);
+                          CCTK_REAL y[], CCTK_REAL x,
+                          CCTK_INT excised_points[]);
 void CCTK_FCALL CCTK_FNAME(BndCartoon2DVI)(int *retval, const cGH **GH, 
                 int *tensortype, int *vi);
 void CCTK_FCALL CCTK_FNAME(BndCartoon2DVN)(int *retval, const cGH **GH,
@@ -49,6 +52,9 @@ void CCTK_FCALL CCTK_FNAME(BndCartoon2DGN)
                           (int *retval, const cGH **GH, int *tensortype, ONE_FORTSTRING_ARG);
 
 
+void Cartoon2D_InitExcisionVars(CCTK_ARGUMENTS);
+void SetupExcisionVars(CCTK_ARGUMENTS);
+
 /* set boundaries of a grid tensor assuming axisymmetry 
    - handles lower boundary in x
    - does not touch other boundaries
@@ -90,7 +96,6 @@ int BndCartoon2DVI(const cGH *cctkGH, int tensortype, int prolongtype, int vi)
   CCTK_REAL f[100], fx, fy, fz, fxx, fxy, fxz, fyy, fyz, fzz;
   CCTK_REAL *t, *tx, *ty, *tz, *txx, *txy, *txz, *tyy, *tyz, *tzz;
 
-
   gi = CCTK_GroupIndexFromVarI (vi);
   assert (gi >= 0);
   n_vars = CCTK_NumVarsInGroupI (gi);
@@ -342,6 +347,59 @@ static CCTK_REAL Cartoon2DInterp(const cGH *cctkGH, CCTK_REAL *f,
   int lnx;
   int n, offset;
 
+  CCTK_REAL *old_mask;
+  CCTK_INT *new_mask;
+
+  CCTK_INT excised_points[6];
+
+  if (*excision_active < 0)
+  {
+    SetupExcisionVars((cGH*)cctkGH);
+  }
+
+  if (old_excision)
+  {
+    old_mask = (CCTK_REAL*)CCTK_VarDataPtrI(cctkGH, 0, *old_mask_vi);
+  }
+  else
+  {
+    old_mask = NULL;
+  }
+  
+  if (new_excision)
+  {
+    new_mask = (CCTK_INT*)CCTK_VarDataPtrI(cctkGH, 0, *new_mask_vi);
+  }
+  else
+  {
+    new_mask = NULL;
+  }
+
+  /* 
+     If this boundary point is excised then we return 
+     without doing anything 
+  */
+
+  if (*excision_active)
+  {
+    if ((*excision_active == 1)||(*excision_active == 3))
+    {
+      if (old_mask[ijk] < 0.75)
+      {
+        return f[ijk];
+      }
+    }
+    if ((*excision_active == 2)||(*excision_active == 3))
+    {
+      if (SpaceMask_CheckStateBits(new_mask, (ijk),
+                                   *new_excision_field,
+                                   *new_excision_descriptor))
+      {
+        return f[ijk];
+      }
+    }
+  }  
+
   lnx = cctkGH->cctk_lsh[0];
     
   dx0 = CCTK_DELTA_SPACE(0);
@@ -369,10 +427,25 @@ static CCTK_REAL Cartoon2DInterp(const cGH *cctkGH, CCTK_REAL *f,
   x0 = lx0 + dx0 * (i-offset);
   for (n = 0; n <= order; n++) {
     y[n] = f[ijk-offset+n];
+    excised_points[n] = 0;
+    if (*excision_active)
+    {
+      if ((*excision_active == 1)||(*excision_active == 3))
+      {
+        excised_points[n] = (old_mask[ijk-offset+n] > 0.75) ? 0 : 1;
+      }
+      if ((*excision_active == 2)||(*excision_active == 3))
+      {
+        excised_points[n] |= 
+          (SpaceMask_CheckStateBits(new_mask, (ijk-offset+n),
+                                    *new_excision_field,
+                                    *new_excision_descriptor));
+      }
+    }
   }
   
   /* call interpolator */
-  r = interpolate_local(order, x0, dx0, y, x);
+  r = interpolate_local(order, x0, dx0, y, x, excised_points);
   
   return r;
 }
@@ -391,6 +464,59 @@ static CCTK_REAL Cartoon2DInterp_ENO(const cGH *cctkGH,
   int lnx;
   int n, offset;
   int tmp_order = 2;
+
+  CCTK_REAL *old_mask;
+  CCTK_INT *new_mask;
+
+  CCTK_INT excised_points[11];
+
+  if (*excision_active < 0)
+  {
+    SetupExcisionVars((cGH*)cctkGH);
+  }
+
+  if (old_excision)
+  {
+    old_mask = (CCTK_REAL*)CCTK_VarDataPtrI(cctkGH, 0, *old_mask_vi);
+  }
+  else
+  {
+    old_mask = NULL;
+  }
+  
+  if (new_excision)
+  {
+    new_mask = (CCTK_INT*)CCTK_VarDataPtrI(cctkGH, 0, *new_mask_vi);
+  }
+  else
+  {
+    new_mask = NULL;
+  }
+
+  /* 
+     If this boundary point is excised then we return 
+     without doing anything 
+  */
+
+  if (*excision_active)
+  {
+    if ((*excision_active == 1)||(*excision_active == 3))
+    {
+      if (old_mask[ijk] < 0.75)
+      {
+        return f[ijk];
+      }
+    }
+    if ((*excision_active == 2)||(*excision_active == 3))
+    {
+      if (SpaceMask_CheckStateBits(new_mask, (ijk),
+                                   *new_excision_field,
+                                   *new_excision_descriptor))
+      {
+        return f[ijk];
+      }
+    }
+  }
   
   lnx = cctkGH->cctk_lsh[0];
   
@@ -424,11 +550,92 @@ static CCTK_REAL Cartoon2DInterp_ENO(const cGH *cctkGH,
   for (n = 0; n < 2 * eno_order + 1; ++n) 
   {
     y[n] = f[ijk-offset+n];
+    excised_points[n] = 0;
+    if (*excision_active)
+    {
+      if ((*excision_active == 1)||(*excision_active == 3))
+      {
+        excised_points[n] = (old_mask[ijk-offset+n] > 0.75) ? 1 : 0;
+      }
+      if ((*excision_active == 2)||(*excision_active == 3))
+      {
+        excised_points[n] |= 
+          (SpaceMask_CheckStateBits(new_mask, (ijk-offset+n),
+                                    *new_excision_field,
+                                    *new_excision_descriptor));
+      }
+    }
   }
   
   /* call interpolator */
 /*   r = interpolate_eno(tmp_order, x0, dx0, y, x); */
-  r = interpolate_eno(eno_order, x0, dx0, y, x);
+  r = interpolate_eno(eno_order, x0, dx0, y, x, excised_points);
   
   return r;
 }
+
+void Cartoon2D_InitExcisionVars(CCTK_ARGUMENTS)
+{
+  DECLARE_CCTK_ARGUMENTS;
+  DECLARE_CCTK_PARAMETERS;
+
+  *excision_active = -1;
+  *old_mask_vi = -1;
+  *new_mask_vi = -1;
+  *new_excision_field = -1;
+  *new_excision_descriptor = -1;
+
+}
+
+void SetupExcisionVars(CCTK_ARGUMENTS)
+{
+  DECLARE_CCTK_ARGUMENTS;
+  DECLARE_CCTK_PARAMETERS;
+
+  if (old_excision)
+  {
+    *old_mask_vi = CCTK_VarIndex(old_style_excision_var);
+  }
+  else
+  {
+    *old_mask_vi = -1;
+  }
+  if (new_excision)
+  {
+    *new_mask_vi = CCTK_VarIndex(new_style_excision_var);
+  }
+  else
+  {
+    *new_mask_vi = -1;
+  }
+
+  if (*old_mask_vi > -1)
+  {
+    if (*new_mask_vi > -1)
+    {
+      *excision_active = 3;
+      *new_excision_field = SpaceMask_GetTypeBits(new_mask_field_name);
+      *new_excision_descriptor = SpaceMask_GetStateBits(new_mask_field_name,
+                                                        new_mask_excised_name);
+    }
+    else
+    {
+      *excision_active = 1;
+    }
+  }
+  else
+  {
+    if (*new_mask_vi > -1)
+    {
+      *excision_active = 2;
+      *new_excision_field = SpaceMask_GetTypeBits(new_mask_field_name);
+      *new_excision_descriptor = SpaceMask_GetStateBits(new_mask_field_name,
+                                                        new_mask_excised_name);
+    }
+    else
+    {
+      *excision_active = 0;
+    }
+  }
+  
+}
diff --git a/src/SymInterp.c b/src/SymInterp.c
index 0efc291..24668ab 100644
--- a/src/SymInterp.c
+++ b/src/SymInterp.c
@@ -360,6 +360,13 @@ 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;
+      }
+    }
   }
   
   
diff --git a/src/interpolate.c b/src/interpolate.c
index 1d774a2..e207cb2 100644
--- a/src/interpolate.c
+++ b/src/interpolate.c
@@ -23,11 +23,13 @@
 
 CCTK_REAL interpolate_local(int order, CCTK_REAL x0, CCTK_REAL dx, 
                             CCTK_REAL y[],
-                            CCTK_REAL x);
+                            CCTK_REAL x,
+                            CCTK_INT excised_points[]);
 
 CCTK_REAL interpolate_eno(CCTK_INT order, CCTK_REAL x0, CCTK_REAL dx,
                           CCTK_REAL y[],
-                          CCTK_REAL x);
+                          CCTK_REAL x,
+                            CCTK_INT excised_points[]);
 
 /* prototypes for private functions defined in this file */
 static CCTK_REAL interpolate_local_order1(CCTK_REAL x0, CCTK_REAL dx,
@@ -75,14 +77,70 @@ static CCTK_REAL interpolate_local_order5(CCTK_REAL x0, CCTK_REAL dx,
 CCTK_REAL interpolate_local(int order,
                             CCTK_REAL x0, CCTK_REAL dx,
                             CCTK_REAL y[],
-                            CCTK_REAL x)
+                            CCTK_REAL x,
+                            CCTK_INT excised_points[])
 {
-  switch (order) {
-  case 1: return interpolate_local_order1(x0, dx, y, x);
-  case 2: return interpolate_local_order2(x0, dx, y, x);
-  case 3: return interpolate_local_order3(x0, dx, y, x);
-  case 4: return interpolate_local_order4(x0, dx, y, x);
-  case 5: return interpolate_local_order5(x0, dx, y, x);
+
+  /* 
+     We assume that the excised region has one of the following 4 forms:
+     
+     o o o o o         (No excision)
+     x o o o o         (left excision)
+     o o o o x         (right excision)
+     x o o o x         (left and right excision)
+
+     Anything else (e.g., x o x o x) will go horribly wrong.
+  */
+
+  int start, end, available_order;
+  
+  start = 0;
+  end   = order;
+  
+  while ((start < order + 1)&&(excised_points[start]))
+  {
+    ++start;
+  }
+  while ((end > -1)&&(excised_points[end]))
+  {
+    --end;
+  }
+  
+  if (end < start)
+  {
+    /* The whole block is excised */
+    return y[0];
+  }
+
+  /* The possible interpolation order that can be used */
+
+  available_order = end - start;
+
+  if (available_order == 0)
+  {
+    /* There is only one non-excised point */
+    return y[start];
+  }
+
+  if (dx * start + x0 > x)
+  {
+    /* We would have to extrapolate to the left */
+    return y[start];
+  }
+  
+  if (dx * end + x0 < x)
+  {
+    /* We would have to extrapolate to the right */
+    return y[end];
+  }
+
+  switch (available_order) 
+  {
+  case 1: return interpolate_local_order1(x0 + start * dx, dx, &y[start], x);
+  case 2: return interpolate_local_order2(x0 + start * dx, dx, &y[start], x);
+  case 3: return interpolate_local_order3(x0 + start * dx, dx, &y[start], x);
+  case 4: return interpolate_local_order4(x0 + start * dx, dx, &y[start], x);
+  case 5: return interpolate_local_order5(x0 + start * dx, dx, &y[start], x);
   default: CCTK_WARN(0, "Interpolation order not supported");
   }
 }
@@ -218,7 +276,8 @@ return(  c0 + xr*(c1 + xr*(c2 + xr*(c3 + xr*(c4 + xr*c5))))  );
 CCTK_REAL interpolate_eno(CCTK_INT order, 
                           CCTK_REAL x0, CCTK_REAL dx,
                           CCTK_REAL y[],
-                          CCTK_REAL x)
+                          CCTK_REAL x,
+                          CCTK_INT excised_points[])
 {
 
 /*   CCTK_REAL diff[4]; */