(Changes from Ian Hawke.)

Possibly use ENO interpolation. Useful e.g. for hydro where shocks may appear.

To make Cartoon use the ENO interpolator you need an additional tags table entry. The same tags are used as for Carpet. Either

tags='Prolongation="TVD"'

or

tags='Prolongation="ENO"'

will work.

The default (i.e., with no tag) is to use the standard Lagrange polynomials. The tags "Lagrange" and "None" will also do this. Also, any direct call to the Cartoon functions (BndCartoon2DVI etc.) will use Lagrange interpolation.

Code from Burkhard Zink for interpolation; tags table stuff from me.


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

5 files changed, 240 insertions, 44 deletions

diff --git a/src/ApplyCartoon.c b/src/ApplyCartoon.c
index a203105..61424e8 100644
--- a/src/ApplyCartoon.c
+++ b/src/ApplyCartoon.c
@@ -25,6 +25,7 @@ CCTK_FILEVERSION(BetaThorns_Cartoon2D_ApplyCartoon_c);
  ********************************************************************/
 /* #define DEBUG 1 */
 #define TENSORTYPE_BUFF_SIZE 7
+#define PROLONG_BUFF_SIZE 1000
 
 /********************************************************************
  ********************* Local Routine Prototypes *********************
@@ -77,11 +78,15 @@ void Cartoon_ApplyBoundaries(CCTK_ARGUMENTS);
    @endreturndesc
 @@*/
 
-void Cartoon_ApplyBoundaries(CCTK_ARGUMENTS) {
+void Cartoon_ApplyBoundaries(CCTK_ARGUMENTS) 
+{
+
   DECLARE_CCTK_ARGUMENTS;
   int num_vars, err, i, gi, group_tags_table;
   CCTK_INT * vars;
   char tensortype[TENSORTYPE_BUFF_SIZE];
+  char prolongtype[PROLONG_BUFF_SIZE];
+  int prolongmethod;
 
   /* Allocate memory to hold selected bcs */
   num_vars = Boundary_SelectedGVs(cctkGH, 0, NULL, NULL, NULL, NULL, NULL);
@@ -145,6 +150,49 @@ void Cartoon_ApplyBoundaries(CCTK_ARGUMENTS) {
 #ifdef DEBUG
     printf("Cartoon_ApplyBoundaries: tensor type is %s\n",tensortype);
 #endif
+
+    /* Get prolongation type from group tags table */
+
+    err = Util_TableGetString(group_tags_table, PROLONG_BUFF_SIZE,
+                              prolongtype, "Prolongation");
+
+    if (err == UTIL_ERROR_TABLE_NO_SUCH_KEY) 
+    {
+      /* Use the default */
+
+      prolongmethod = PROLONG_LAGRANGE;
+    }
+    else if (err < 0)
+    {
+      CCTK_VWarn(0, __LINE__, __FILE__, CCTK_THORNSTRING,
+		 "Error (%d) in TAGS table for variable group %s", err,
+		 CCTK_GroupName(gi));
+    }
+    else
+    {
+      if (CCTK_Equals(prolongtype, "None")) 
+      {
+        prolongmethod = PROLONG_LAGRANGE; /* But why? */
+      } 
+      else if (CCTK_Equals(prolongtype, "Lagrange")) 
+      {
+        prolongmethod = PROLONG_LAGRANGE;
+      } 
+      else if (CCTK_Equals(prolongtype, "TVD")) 
+      {
+        prolongmethod = PROLONG_ENO;
+      } 
+      else if (CCTK_Equals(prolongtype, "ENO")) 
+      {
+        prolongmethod = PROLONG_ENO;
+      } 
+      else 
+      {
+        CCTK_VWarn(0, __LINE__, __FILE__, CCTK_THORNSTRING,
+                   "Error in TAGS table for variable group %s",
+                   CCTK_GroupName(gi));
+      }
+    }    
     
     /* Here one should check that the group's size is correct for the
      * specified tensortype.
@@ -154,13 +202,13 @@ void Cartoon_ApplyBoundaries(CCTK_ARGUMENTS) {
        macro */
     if (CCTK_Equals(tensortype, "scalar"))
     {
-      BndCartoon2DVI(cctkGH, TENSORTYPE_SCALAR, vars[i]);
+      BndCartoon2DVI(cctkGH, TENSORTYPE_SCALAR, prolongmethod, vars[i]);
     } else if (CCTK_Equals(tensortype, "u"))
     {
-      BndCartoon2DVI(cctkGH, TENSORTYPE_U, vars[i]);
+      BndCartoon2DVI(cctkGH, TENSORTYPE_U, prolongmethod, vars[i]);
     } else if (CCTK_Equals(tensortype, "dd_sym"))
     {
-      BndCartoon2DVI(cctkGH, TENSORTYPE_DDSYM, vars[i]);
+      BndCartoon2DVI(cctkGH, TENSORTYPE_DDSYM, prolongmethod, vars[i]);
     } else 
     {
       CCTK_VWarn(1, __LINE__, __FILE__, CCTK_THORNSTRING,
diff --git a/src/Cartoon2D.h b/src/Cartoon2D.h
index 930c647..476e985 100644
--- a/src/Cartoon2D.h
+++ b/src/Cartoon2D.h
@@ -3,7 +3,7 @@
 extern "C" {
 #endif
 
-int BndCartoon2DVI(const cGH *GH, int tensortype, int vi);
+int BndCartoon2DVI(const cGH *GH, int tensortype, int prolongtype, int vi);
 
 #ifdef __cplusplus 
 }
diff --git a/src/Cartoon2DBC.c b/src/Cartoon2DBC.c
index 840c042..2added6 100644
--- a/src/Cartoon2DBC.c
+++ b/src/Cartoon2DBC.c
@@ -19,6 +19,7 @@
 #include <math.h>
 
 #include "cctk.h"
+#include "cctk_Arguments.h"
 #include "cctk_Parameters.h"
 #include "cctk_FortranString.h"
 
@@ -30,11 +31,16 @@ CCTK_FILEVERSION(BetaThorns_Cartoon2D_Cartoon2DBC_c);
 
 int BndCartoon2DGN(const cGH *GH, int tensortype, const char *group);
 int BndCartoon2DVN(const cGH *GH, int tensortype, const char *var);
-int BndCartoon2DVI(const cGH *GH, int tensortype, int vi);
+int BndCartoon2DVI(const cGH *GH, int tensortype, int prolongtype, int vi);
 static CCTK_REAL Cartoon2DInterp(const cGH *GH, CCTK_REAL *f, int i, int di,
                              int ijk, CCTK_REAL x);
+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 interpolate_eno(CCTK_REAL x0, CCTK_REAL dx,
+                          CCTK_REAL y[], CCTK_REAL x);
 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,
@@ -63,10 +69,11 @@ void CCTK_FCALL CCTK_FNAME(BndCartoon2DGN)
    group, using the tensortype argument to determine how many there
    are and how to treat them. */
 
-int BndCartoon2DVI(const cGH *GH, int tensortype, int vi)
+int BndCartoon2DVI(const cGH *cctkGH, int tensortype, int prolongtype, int vi)
 
 {
-  DECLARE_CCTK_PARAMETERS
+  DECLARE_CCTK_ARGUMENTS;
+  DECLARE_CCTK_PARAMETERS;
   /*** FIXME: can CCTK_SyncGroup() have a 'const cGH *' parameter ?? ***/
   union
   {
@@ -74,14 +81,14 @@ int BndCartoon2DVI(const cGH *GH, int tensortype, int vi)
     cGH *non_const_ptr;
   } GH_fake_const;
 
-  static int xindx=-1;
+/*   static int xindx=-1; */
 
   int i, j, k, di, dj, ijk;
   int jj, n, s;
   int lnx, lnz, ny;
   int n_tensortype;
   int timelevel;
-  CCTK_REAL *xp;
+/*   CCTK_REAL *xp; */
   CCTK_REAL x, y, rho;
   CCTK_REAL lx0, dx0, dy0;
   CCTK_REAL rxx, rxy, ryx, ryy;
@@ -105,18 +112,26 @@ int BndCartoon2DVI(const cGH *GH, int tensortype, int vi)
     n_tensortype = 0;
   }
     
-  s = GH->cctk_nghostzones[0];
-  lnx = GH->cctk_lsh[0];
-  lnz = GH->cctk_lsh[2];
-  ny  = GH->cctk_gsh[1];
+  s = cctkGH->cctk_nghostzones[0];
+  lnx = cctkGH->cctk_lsh[0];
+  lnz = cctkGH->cctk_lsh[2];
+  ny  = cctkGH->cctk_gsh[1];
 
+  /*
   xindx = (xindx<0)?CCTK_CoordIndex(-1,"x","cart3d"):xindx;
-  xp    = GH->data[xindx][0];
+  xp    = cctkGH->data[xindx][0];
+  */
 
   /*  xp  = GH->data[CCTK_CoordIndex(-1,"x","cart3d")][0];*/
+  /*
   lx0 = xp[0];
-  dx0 = GH->cctk_delta_space[0]/GH->cctk_levfac[0];
-  dy0 = GH->cctk_delta_space[1]/GH->cctk_levfac[1];
+  dx0 = cctkGH->cctk_delta_space[0]/cctkGH->cctk_levfac[0];
+  dy0 = cctkGH->cctk_delta_space[1]/cctkGH->cctk_levfac[1];
+  */
+  
+  lx0 = CCTK_ORIGIN_SPACE(0);
+  dx0 = CCTK_DELTA_SPACE(0);
+  dy0 = CCTK_DELTA_SPACE(1);
 
   timelevel = 0; 
 
@@ -131,7 +146,7 @@ int BndCartoon2DVI(const cGH *GH, int tensortype, int vi)
 
   /* this union helps us to avoid compiler warning about discarding
      the const qualifier from a pointer target type */
-  GH_fake_const.const_ptr = GH;
+  GH_fake_const.const_ptr = cctkGH;
   /* make sure that the input data is synced */
   CCTK_SyncGroup(GH_fake_const.non_const_ptr, CCTK_GroupNameFromVarI(vi)); 
 
@@ -145,7 +160,7 @@ int BndCartoon2DVI(const cGH *GH, int tensortype, int vi)
   for (i = -s; i < lnx; i++) {
 
     /* index into linear memory */
-    ijk = CCTK_GFINDEX3D(GH,i,j,k);
+    ijk = CCTK_GFINDEX3D(cctkGH,i,j,k);
    
     /* what a neat way to hack in the zombie for x < 0, y = 0 */
     if (i < 0) {i += s; ijk += s; dj = 0;}
@@ -171,19 +186,32 @@ int BndCartoon2DVI(const cGH *GH, int tensortype, int vi)
     syy = rxx;
 
     /* interpolate grid functions */
-    for (n = 0; n < n_tensortype; n++)
-      f[n] = Cartoon2DInterp(GH, GH->data[vi+n][timelevel], i, di, ijk, rho); 
+    switch (prolongtype)
+    {
+      case PROLONG_LAGRANGE:
+        for (n = 0; n < n_tensortype; n++)
+          f[n] = Cartoon2DInterp(cctkGH, cctkGH->data[vi+n][timelevel], 
+                                 i, di, ijk, rho); 
+        break;
+      case PROLONG_ENO:
+        for (n = 0; n < n_tensortype; n++)
+          f[n] = Cartoon2DInterp_ENO(cctkGH, cctkGH->data[vi+n][timelevel], i, 
+                                     di, ijk, rho); 
+        break;
+      default:
+        CCTK_WARN(0, "Prolongation type not recognized by Cartoon2D.");
+    }
 
     /* rotate grid tensor by matrix multiplication */
     if (tensortype == TENSORTYPE_SCALAR) {
-      t = GH->data[vi][timelevel];
+      t = cctkGH->data[vi][timelevel];
       t[ijk+dj] = f[0];
       t[ijk-dj] = f[0];
     }
     else if (tensortype == TENSORTYPE_U) {
-      tx = GH->data[vi][timelevel];
-      ty = GH->data[vi+1][timelevel];
-      tz = GH->data[vi+2][timelevel];
+      tx = cctkGH->data[vi][timelevel];
+      ty = cctkGH->data[vi+1][timelevel];
+      tz = cctkGH->data[vi+2][timelevel];
       fx = f[0];
       fy = f[1];
       fz = f[2];
@@ -197,12 +225,12 @@ int BndCartoon2DVI(const cGH *GH, int tensortype, int vi)
       tz[ijk-dj] = fz;
     }
     else if (tensortype == TENSORTYPE_DDSYM) {
-      txx = GH->data[vi][timelevel];
-      txy = GH->data[vi+1][timelevel];
-      txz = GH->data[vi+2][timelevel];
-      tyy = GH->data[vi+3][timelevel];
-      tyz = GH->data[vi+4][timelevel];
-      tzz = GH->data[vi+5][timelevel];
+      txx = cctkGH->data[vi][timelevel];
+      txy = cctkGH->data[vi+1][timelevel];
+      txz = cctkGH->data[vi+2][timelevel];
+      tyy = cctkGH->data[vi+3][timelevel];
+      tyz = cctkGH->data[vi+4][timelevel];
+      tzz = cctkGH->data[vi+5][timelevel];
       fxx = f[0];
       fxy = f[1];
       fxz = f[2];
@@ -234,7 +262,7 @@ int BndCartoon2DVI(const cGH *GH, int tensortype, int vi)
 
   /* this union helps us to avoid compiler warning about discarding
      the const qualifier from a pointer target type */
-  GH_fake_const.const_ptr = GH;
+  GH_fake_const.const_ptr = cctkGH;
   /* syncs needed after interpolation (actually only for x direction) */
   CCTK_SyncGroup(GH_fake_const.non_const_ptr, CCTK_GroupNameFromVarI(vi)); 
 
@@ -244,7 +272,7 @@ int BndCartoon2DVI(const cGH *GH, int tensortype, int vi)
 void CCTK_FCALL CCTK_FNAME(BndCartoon2DVI)
                           (int *retval, const cGH **GH, int *tensortype, int *vi)
 {
-  *retval = BndCartoon2DVI(*GH, *tensortype, *vi);
+  *retval = BndCartoon2DVI(*GH, *tensortype, PROLONG_LAGRANGE, *vi);
 }
 
 int BndCartoon2DVN(const cGH *GH, int tensortype, const char *inpvarname)
@@ -256,7 +284,7 @@ int BndCartoon2DVN(const cGH *GH, int tensortype, const char *inpvarname)
   if (verbose) printf("Cartoon2D called for %s\n", inpvarname);
   vi = CCTK_VarIndex(inpvarname);
 
-  return(BndCartoon2DVI(GH, tensortype, vi));
+  return(BndCartoon2DVI(GH, tensortype, PROLONG_LAGRANGE, vi));
 }
 
 void CCTK_FCALL CCTK_FNAME(BndCartoon2DVN)
@@ -299,7 +327,8 @@ int BndCartoon2DGN(const cGH *GH, int tensortype, const char *group)
                 "%s should have only %d component(s).", group, n_tensortype);
   }
   
-  return(BndCartoon2DVI(GH, tensortype, CCTK_FirstVarIndex(group)));
+  return(BndCartoon2DVI(GH, tensortype, PROLONG_LAGRANGE, 
+                        CCTK_FirstVarIndex(group)));
 }
 
 void CCTK_FCALL CCTK_FNAME(BndCartoon2DGN)
@@ -311,31 +340,38 @@ void CCTK_FCALL CCTK_FNAME(BndCartoon2DGN)
 }
 
 /* interpolation on x-axis */
-static CCTK_REAL Cartoon2DInterp(const cGH *GH, CCTK_REAL *f, int i, int di, 
-              int ijk, CCTK_REAL x)
+static CCTK_REAL Cartoon2DInterp(const cGH *cctkGH, CCTK_REAL *f, 
+                                 int i, int di, int ijk, CCTK_REAL x)
 {
-  DECLARE_CCTK_PARAMETERS
+  DECLARE_CCTK_ARGUMENTS;
+  DECLARE_CCTK_PARAMETERS;
 
-  static int xindx=-1;
+/*   static int xindx=-1; */
 
   CCTK_REAL x0;
   CCTK_REAL lx0, dx0;
   CCTK_REAL y[6], r;
-  CCTK_REAL *xp;
+/*   CCTK_REAL *xp; */
   int lnx;
   int n, offset;
 
   (void) (di + 0);
 
-  lnx = GH->cctk_lsh[0];
+  lnx = cctkGH->cctk_lsh[0];
 
+  /*
   xindx = (xindx<0)?CCTK_CoordIndex(-1,"x","cart3d"):xindx;
-  xp    = GH->data[xindx][0];
-
+  xp    = cctkGH->data[xindx][0];
+  */
 
 /*  xp = GH->data[CCTK_CoordIndex(-1,"x","cart3d")][0];*/
+/*
   lx0 = xp[0];
-  dx0 = GH->cctk_delta_space[0]/GH->cctk_levfac[0];
+  dx0 = cctkGH->cctk_delta_space[0]/cctkGH->cctk_levfac[0];
+*/
+    
+  lx0 = CCTK_ORIGIN_SPACE(0);
+  dx0 = CCTK_DELTA_SPACE(0);
 
   /* find i such that x(i) < x <= x(i+1)
      for rotation on entry always x > x(i), but sometimes also x > x(i+1) */
@@ -366,3 +402,56 @@ static CCTK_REAL Cartoon2DInterp(const cGH *GH, CCTK_REAL *f, int i, int di,
   
   return r;
 }
+
+static CCTK_REAL Cartoon2DInterp_ENO(const cGH *cctkGH, 
+                                    CCTK_REAL *f, int i, int di, 
+                                    int ijk, CCTK_REAL x)
+{
+  
+  DECLARE_CCTK_PARAMETERS;
+  DECLARE_CCTK_ARGUMENTS;
+  
+/*   static int xindx=-1; */
+  
+  CCTK_REAL x0;
+  CCTK_REAL lx0, dx0;
+  CCTK_REAL y[5], r;
+/*   CCTK_REAL *xp; */
+  int lnx;
+  int n, offset;
+  
+  lnx = cctkGH->cctk_lsh[0];
+  
+  lx0 = CCTK_ORIGIN_SPACE(0);
+  dx0 = CCTK_DELTA_SPACE(0);
+  
+  /* find i such that x(i) < x <= x(i+1)
+     for rotation on entry always x > x(i), but sometimes also x > x(i+1) */
+  while (x > lx0 + dx0 * (i+1)) {i++; ijk++;}
+  
+  /* first attempt to interface to JT's interpolator */
+  
+  /* offset of stencil, note that rotation leads to x close to x0 
+     for large x */
+  offset = 2;
+  
+  /* shift stencil at boundaries */
+  /* note that for simplicity we don't distinguish between true boundaries
+     and decomposition boundaries: the sync fixes that! */
+  if (i-offset < 0) offset = i;
+  if (i-offset+4 >= lnx) offset = i+5-lnx;
+  
+  /* fill in info */
+  /* fills in old data near axis, but as long as stencil at axis is 
+     centered, no interpolation happens anyway */  
+  x0 = lx0 + dx0 * (i-offset);
+  for (n = 0; n <= 4; n++) 
+  {
+    y[n] = f[ijk-offset+n];
+  }
+  
+  /* call interpolator */
+  r = interpolate_eno(x0, dx0, y, x);
+  
+  return r;
+}
diff --git a/src/Cartoon2D_tensors.h b/src/Cartoon2D_tensors.h
index 4d0092b..29c72ec 100644
--- a/src/Cartoon2D_tensors.h
+++ b/src/Cartoon2D_tensors.h
@@ -10,3 +10,7 @@
 #define N_TENSORTYPE_SCALAR 1
 #define N_TENSORTYPE_U 3
 #define N_TENSORTYPE_DDSYM 6
+
+/* Prolongation methods */
+#define PROLONG_LAGRANGE 1
+#define PROLONG_ENO 2
diff --git a/src/interpolate.c b/src/interpolate.c
index 8a39c44..d48da6e 100644
--- a/src/interpolate.c
+++ b/src/interpolate.c
@@ -16,6 +16,7 @@
 */
 
 #include "cctk.h"
+#include "math.h"
 
 #define PRIVATE 
 #define PUBLIC 
@@ -40,6 +41,9 @@ PRIVATE CCTK_REAL interpolate_local_order5(CCTK_REAL x0, CCTK_REAL dx,
 					CCTK_REAL y[],
 					CCTK_REAL x);
 
+PRIVATE CCTK_REAL interpolate_eno(CCTK_REAL x0, CCTK_REAL dx,
+                                  CCTK_REAL y[],
+                                  CCTK_REAL x);
 
 /*****************************************************************************/
 
@@ -225,3 +229,54 @@ CCTK_REAL xr = (x - xc) / dx;
 return(  c0 + xr*(c1 + xr*(c2 + xr*(c3 + xr*(c4 + xr*c5))))  );
 }
 
+
+/*****************************************************************************/
+
+CCTK_REAL interpolate_eno(CCTK_REAL x0, CCTK_REAL dx,
+					CCTK_REAL y[],
+					CCTK_REAL x)
+{
+
+  CCTK_REAL diff[4];
+  CCTK_INT seed = 0;
+  CCTK_INT j,r;  
+
+  /* Find seed index */
+  while (x > x0 + dx * ((float)seed+0.5)) seed++;
+
+  if (seed!=2) {
+    /* Not enough stencil, only perform linear interpolation */   
+    seed = 0;
+    while (x > x0 + dx * ((float)(seed+1))) seed++;
+    if (seed==4) seed=3;
+    return y[seed] + (y[seed+1]-y[seed])/dx * (x-x0-(float)seed*dx);
+  }
+
+  /*  Calculate the undivided differences */
+  for (j=0; j<=3; j++)
+    diff[j] = y[j+1] - y[j];
+
+  /*  Find the stencil */
+  if ( fabs(diff[1]) < fabs(diff[2]) ) {
+    if ( fabs(diff[1]-diff[0]) < fabs(diff[2]-diff[1]) )
+      r = 0;
+    else
+      r = 1;
+  } else {
+    if ( fabs(diff[2]-diff[1]) < fabs(diff[3]-diff[2]) )
+      r = 1;
+    else
+      r = 2;
+  }
+
+  /* Interpolate second order */
+  double c0 = y[r+1];
+  double c1 = (-1.0/2.0)*y[r] + (+1.0/2.0)*y[r+2];
+  double c2 = (+1.0/2.0)*y[r] - y[r+1] + (+1.0/2.0)*y[r+2];
+
+  double xc = x0 + dx * (double)(r+1);
+  double xr = (x - xc) / dx;
+
+  return(  c0 + xr*(c1 + xr*c2)  );  
+
+}