Modify the ENO interpolators so that it takes an "order" parameter.

Add knowledge of the PROLONG_NONE type to lower level bits of the
code, in case that case propagates.


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

4 files changed, 173 insertions, 43 deletions

diff --git a/param.ccl b/param.ccl
index f334526..3f95f13 100644
--- a/param.ccl
+++ b/param.ccl
@@ -20,6 +20,11 @@ INT order "Cartoon's interpolation order"
   1:5 :: "From linear to fifth order."
 } 4
 
+INT eno_order "The interpolation order applied to the ENO interpolator"
+{
+  1:5 :: "From linear to fifth order."
+} 4
+
 BOOLEAN allow_grid_resize "Allow grid to be resized in a cartoon-compatible way"
 {
-} "no"
\ No newline at end of file
+} "no"
diff --git a/src/Cartoon2DBC.c b/src/Cartoon2DBC.c
index ba995d0..f74c30f 100644
--- a/src/Cartoon2DBC.c
+++ b/src/Cartoon2DBC.c
@@ -39,7 +39,7 @@ static CCTK_REAL Cartoon2DInterp_ENO(const cGH *cctkGH,
                                     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 interpolate_eno(CCTK_INT order, 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);
@@ -198,6 +198,10 @@ int BndCartoon2DVI(const cGH *cctkGH, int tensortype, int prolongtype, int vi)
     /* interpolate grid functions */
     switch (prolongtype)
     {
+      case PROLONG_NONE:
+        /* no-op */
+        return (0);
+        break;
       case PROLONG_LAGRANGE:
         for (n = 0; n < n_vars; n++)
           f[n] = Cartoon2DInterp(cctkGH, cctkGH->data[var0+n][timelevel], 
@@ -383,9 +387,10 @@ static CCTK_REAL Cartoon2DInterp_ENO(const cGH *cctkGH,
   
   CCTK_REAL x0;
   CCTK_REAL lx0, dx0;
-  CCTK_REAL y[5], r;
+  CCTK_REAL y[11], r;
   int lnx;
   int n, offset;
+  int tmp_order = 2;
   
   lnx = cctkGH->cctk_lsh[0];
   
@@ -400,25 +405,30 @@ static CCTK_REAL Cartoon2DInterp_ENO(const cGH *cctkGH,
   
   /* offset of stencil, note that rotation leads to x close to x0 
      for large x */
-  offset = 2;
+  offset = eno_order;
+/*   offset = tmp_order; */
   
   /* 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;
+/*   if (i-offset+4 >= lnx) offset = i+5-lnx; */
+/*   if (i-offset+2*tmp_order >= lnx) offset = i+2*tmp_order+1-lnx; */
+  if (i-offset+2*eno_order >= lnx) offset = i+2*eno_order+1-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++) 
+/*   for (n = 0; n < 2 * tmp_order + 1; ++n)  */
+  for (n = 0; n < 2 * eno_order + 1; ++n) 
   {
     y[n] = f[ijk-offset+n];
   }
   
   /* call interpolator */
-  r = interpolate_eno(x0, dx0, y, x);
+/*   r = interpolate_eno(tmp_order, x0, dx0, y, x); */
+  r = interpolate_eno(eno_order, x0, dx0, y, x);
   
   return r;
 }
diff --git a/src/Cartoon2D_tensors.h b/src/Cartoon2D_tensors.h
index bbb4fa8..7670eab 100644
--- a/src/Cartoon2D_tensors.h
+++ b/src/Cartoon2D_tensors.h
@@ -9,3 +9,4 @@
 /* Prolongation methods */
 #define PROLONG_LAGRANGE 1
 #define PROLONG_ENO 2
+#define PROLONG_NONE 3
diff --git a/src/interpolate.c b/src/interpolate.c
index d990cf5..1d774a2 100644
--- a/src/interpolate.c
+++ b/src/interpolate.c
@@ -17,6 +17,7 @@
  */
 
 #include <math.h>
+#include <stdio.h>
 
 #include "cctk.h"
 
@@ -24,7 +25,7 @@ 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 interpolate_eno(CCTK_INT order, CCTK_REAL x0, CCTK_REAL dx,
                           CCTK_REAL y[],
                           CCTK_REAL x);
 
@@ -214,58 +215,171 @@ return(  c0 + xr*(c1 + xr*(c2 + xr*(c3 + xr*(c4 + xr*c5))))  );
 
 /*****************************************************************************/
 
-CCTK_REAL interpolate_eno(CCTK_REAL x0, CCTK_REAL dx,
+CCTK_REAL interpolate_eno(CCTK_INT order, 
+                          CCTK_REAL x0, CCTK_REAL dx,
                           CCTK_REAL y[],
                           CCTK_REAL x)
 {
 
-  CCTK_REAL diff[4];
-  CCTK_INT seed = 0;
-  CCTK_INT j,r;  
+/*   CCTK_REAL diff[4]; */
+/*   CCTK_INT seed = 0; */
+/*   CCTK_INT j,r;   */
 
-  CCTK_REAL c0;
-  CCTK_REAL c1;
-  CCTK_REAL c2;
+/*   CCTK_REAL c0; */
+/*   CCTK_REAL c1; */
+/*   CCTK_REAL c2; */
 
-  CCTK_REAL xc;
-  CCTK_REAL xr;
+/*   CCTK_REAL xc; */
+/*   CCTK_REAL xr; */
+
+  CCTK_REAL result;
+  CCTK_REAL undiv_diff [13][6];
+  CCTK_INT diff, i, ii, r;
+  CCTK_REAL fp_i, fp_ii, x_rel;
 
   /* Find seed index */
-  while (x > x0 + dx * ((CCTK_REAL)seed+0.5)) seed++;
+/*   while (x > x0 + dx * ((CCTK_REAL)seed+0.5)) seed++; */
 
-  if (seed!=2) {
+/*   if (seed!=2) { */
     /* Not enough stencil, only perform linear interpolation */   
-    seed = 0;
-    while (x > x0 + dx * ((CCTK_REAL)(seed+1))) seed++;
-    if (seed==4) seed=3;
-    return y[seed] + (y[seed+1]-y[seed])/dx * (x-x0-(CCTK_REAL)seed*dx);
-  }
+/*     seed = 0; */
+/*     while (x > x0 + dx * ((CCTK_REAL)(seed+1))) seed++; */
+/*     if (seed==4) seed=3; */
+/*     return y[seed] + (y[seed+1]-y[seed])/dx * (x-x0-(CCTK_REAL)seed*dx); */
+/*   } */
 
   /*  Calculate the undivided differences */
-  for (j=0; j<=3; j++)
-    diff[j] = y[j+1] - y[j];
+/*   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;
-  }
+/*   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 */
-  c0 = y[r+1];
-  c1 = (-1.0/2.0)*y[r] + (+1.0/2.0)*y[r+2];
-  c2 = (+1.0/2.0)*y[r] - y[r+1] + (+1.0/2.0)*y[r+2];
+/*   c0 = y[r+1]; */
+/*   c1 = (-1.0/2.0)*y[r] + (+1.0/2.0)*y[r+2]; */
+/*   c2 = (+1.0/2.0)*y[r] - y[r+1] + (+1.0/2.0)*y[r+2]; */
 
-  xc = x0 + dx * (CCTK_REAL)(r+1);
-  xr = (x - xc) / dx;
+/*   xc = x0 + dx * (CCTK_REAL)(r+1); */
+/*   xr = (x - xc) / dx; */
+/*   result = (  c0 + xr*(c1 + xr*c2)  ); */
 
-  return(  c0 + xr*(c1 + xr*c2)  );  
+  for (i = 1; i < 2 * order + 2; ++i)
+  {
+    undiv_diff[i][0] = y[i-1];
+  }
+  for (i = 0; i < 6; ++i)
+  {
+    undiv_diff[0][i] = 1e10;
+    undiv_diff[2 * order + 2][i] = 1e10;
+  }
+
+  for (diff = 1; diff < order + 1; ++diff)
+  {
+    for (i = 1; i < 2 * order + 2; ++i)
+    {
+      undiv_diff[i][diff] = undiv_diff[i][diff-1] - undiv_diff[i-1][diff-1];
+    }
+    undiv_diff[0][diff] = -10 * undiv_diff[1][diff];
+    undiv_diff[2 * order + 2][diff]= -10 * undiv_diff[2 * order + 1][diff];
+  }
+  
+  fp_i = (x - x0) / dx;
+  fp_ii = floor(fp_i);
+  x_rel = fp_i - fp_ii;
+  
+  ii = (CCTK_INT)(fp_ii) + 1;
+  
+  if (ii < 1)
+  {
+    ii = 1;
+    x_rel = fp_i - (CCTK_REAL)(ii) + 1.0;
+  }
+  else if (ii > 2 * order + 1)
+  {
+    ii = 2 * order + 1;
+    x_rel = fp_i - (CCTK_REAL)(ii) + 1.0;
+  }
+  
+  r = 0;
+  for (diff = 1; diff < order + 1; ++diff)
+  {
+    if (fabs(undiv_diff[ii-r][diff]) < fabs(undiv_diff[ii-r+1][diff]) )
+    {
+      ++r;
+    }
+    
+    if (ii - r < diff + 1)
+    {
+      --r;
+    }
+    else if (ii - r + diff > 2 * order + 1)
+    {
+      ++r;
+    }
+  }
+
+  result = undiv_diff[ii-r][0];
+
+  switch (order)
+  {
+    case 1:
+      result += (x_rel + r - 0.0) / 1.0 * undiv_diff[ii-r+1][1];
+      break;
+    case 2:
+      result += (x_rel + r - 0.0) / 1.0 * (undiv_diff[ii-r+1][1] +
+                (x_rel + r - 1.0) / 2.0 * (undiv_diff[ii-r+2][2]));
+      break;
+    case 3:
+      result += (x_rel + r - 0.0) / 1.0 * (undiv_diff[ii-r+1][1] +
+                (x_rel + r - 1.0) / 2.0 * (undiv_diff[ii-r+2][2] +
+                (x_rel + r - 2.0) / 3.0 * (undiv_diff[ii-r+3][3])));
+      break;
+    case 4:
+      result += (x_rel + r - 0.0) / 1.0 * (undiv_diff[ii-r+1][1] +
+                (x_rel + r - 1.0) / 2.0 * (undiv_diff[ii-r+2][2] +
+                (x_rel + r - 2.0) / 3.0 * (undiv_diff[ii-r+3][3] +
+                (x_rel + r - 3.0) / 4.0 * (undiv_diff[ii-r+4][4]))));
+      break;
+    case 5:
+      result += (x_rel + r - 0.0) / 1.0 * (undiv_diff[ii-r+1][1] +
+                (x_rel + r - 1.0) / 2.0 * (undiv_diff[ii-r+2][2] +
+                (x_rel + r - 2.0) / 3.0 * (undiv_diff[ii-r+3][3] +
+                (x_rel + r - 3.0) / 4.0 * (undiv_diff[ii-r+4][4] +
+                (x_rel + r - 4.0) / 5.0 * (undiv_diff[ii-r+5][5])))));
+      break;
+  }
+  
+/* #define CARTOON_ENO_DBG 1 */
+#ifdef CARTOON_ENO_DBG
+
+  printf("Cartoon ENO interpolation.\n"
+         "Input: order %d x0 %g dx %g x %g.\n",
+         order, x0, dx, x);
+  printf("Data is\n");
+  for (i = 0; i < 2 * order + 1; ++i)
+  {
+    printf(" %g", y[i]);
+  }
+  printf("\nEntries used were\n");
+  printf("0: %g   1: %g   2: %g\n",
+         undiv_diff[ii-r][0], undiv_diff[ii-r+1][1], undiv_diff[ii-r+2][2]);
+  printf("Parameters are fp_i %g fp_ii %g x_rel %g ii %d r %d\n",
+         fp_i, fp_ii, x_rel, ii, r);
+  printf("Result is %g\n", result);
+  
+#endif
+
+  return result;  
 
 }