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

1 files changed, 212 insertions, 5 deletions

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;
+    }
+  }
+  
+}