Added ability to use the boundary normal vector in a calculation

2 files changed, 205 insertions, 69 deletions

diff --git a/Auxiliary/Cactus/KrancNumericalTools/GenericFD/src/GenericFD.c b/Auxiliary/Cactus/KrancNumericalTools/GenericFD/src/GenericFD.c
index e4630e8..c34c1d7 100644
--- a/Auxiliary/Cactus/KrancNumericalTools/GenericFD/src/GenericFD.c
+++ b/Auxiliary/Cactus/KrancNumericalTools/GenericFD/src/GenericFD.c
@@ -33,6 +33,8 @@
 #include "cctk_Parameters.h"                  
 #include "util_Table.h"
 #include <assert.h>
+#include <stdlib.h>
+#include <math.h>
             
 #include "Symmetry.h"                         
 
@@ -154,84 +156,103 @@ void GenericFD_LoopOverEverything(cGH *cctkGH, Kranc_Calculation calc)
   CCTK_INT   dir = 0;
   CCTK_INT   face = 0;
   CCTK_REAL  normal[] = {0,0,0};
-  CCTK_REAL  tangent1[] = {0,0,0};
-  CCTK_REAL  tangent2[] = {0,0,0};
+  CCTK_REAL  tangentA[] = {0,0,0};
+  CCTK_REAL  tangentB[] = {0,0,0};
   CCTK_INT   bmin[] = {0,0,0};
   CCTK_INT   bmax[] = {cctk_lsh[0], cctk_lsh[1], cctk_lsh[2]};
 
-  calc(cctkGH, dir, face, normal, tangent1, tangent2, bmin, bmax, 0, NULL);
+  calc(cctkGH, dir, face, normal, tangentA, tangentB, bmin, bmax, 0, NULL);
   return;
 }
 
-static void basis_on_boundary(int dir, int face, CCTK_REAL normal[3], CCTK_REAL tangent1[3], 
-                              CCTK_REAL tangent2[3])
-{
-  for (int i = 0; i < 3; i++)
-  {
-    normal[i] = 0.0;
-    tangent1[i] = 0.0;
-    tangent2[i] = 0.0;
-  }
-
-  // Do we need these to be chosen in a more intelligent way?  For
-  // example, do they need to form a right handed set, or a set of
-  // constant handedness?
-  normal[dir] = (face == 0) ? -1 : 1;
-  tangent1[(dir + 1) % 3] = 1;
-  tangent2[(dir + 2) % 3] = 1;
-}
-
 
 void GenericFD_LoopOverBoundary(cGH *cctkGH, Kranc_Calculation calc)
 {
   DECLARE_CCTK_ARGUMENTS
 
-  CCTK_INT   dir = 0;
-  CCTK_INT   face = 0;
-  CCTK_REAL  normal[] = {0,0,0};
-  CCTK_REAL  tangent1[] = {0,0,0};
-  CCTK_REAL  tangent2[] = {0,0,0};
-  CCTK_INT   bmin[] = {0,0,0};
-  CCTK_INT   bmax[] = {cctk_lsh[0], cctk_lsh[1], cctk_lsh[2]};
+  CCTK_INT   dir1, dir2, dir3;
+  CCTK_INT   dir[3];
+  CCTK_REAL  normal[3];
+  CCTK_REAL  tangentA[3];
+  CCTK_REAL  tangentB[3];
+  CCTK_INT   bmin[3];
+  CCTK_INT   bmax[3];
+  CCTK_INT   here_is_physbnd;
+  CCTK_INT   d;
 
   CCTK_INT   is_symbnd[6], is_physbnd[6], is_ipbnd[6];
   CCTK_INT   imin[3], imax[3];
+  int        old_dir = 0;
+  int        old_face = 0;
 
   GenericFD_GetBoundaryInfo(cctkGH, cctk_lsh, cctk_bbox, cctk_nghostzones, 
                             imin, imax, is_symbnd, is_physbnd, is_ipbnd);
 
  /* Loop over all faces */
-  for (dir = 0; dir < 3; dir++)
+  for (dir3 = -1; dir3 <= +1; dir3++)
   {
-    for (face = 0; face < 2; face++)
+    for (dir2 = -1; dir2 <= +1; dir2++)
     {
-      /* Start by looping over the whole grid, minus the NON-PHYSICAL
-         boundary points, which are set by synchronization.  */
-      bmin[0] = is_physbnd[0*2+0] ? 0 : imin[0];
-      bmin[1] = is_physbnd[1*2+0] ? 0 : imin[1];
-      bmin[2] = is_physbnd[2*2+0] ? 0 : imin[2];
-      bmax[0] = is_physbnd[0*2+1] ? cctk_lsh[0] : imax[0];
-      bmax[1] = is_physbnd[1*2+1] ? cctk_lsh[1] : imax[1];
-      bmax[2] = is_physbnd[2*2+1] ? cctk_lsh[2] : imax[2];
-
-      /* Now restrict to only the boundary points on the current face */
-      switch(face)
+      for (dir1 = -1; dir1 <= +1; dir1++)
       {
-      case 0:
-        bmax[dir] = imin[dir];
-        bmin[dir] = 0;
-        break;
-      case 1:
-        bmin[dir] = imax[dir];
-        bmax[dir] = cctk_lsh[dir];
-        break;
-      }
+        dir[0] = dir1;
+        dir[1] = dir2;
+        dir[2] = dir3;
+
+        here_is_physbnd = 0;
+
+        /* Start by looping over the whole grid, minus the NON-PHYSICAL
+           boundary points, which are set by synchronization.  */
+        for (d = 0; d < 3; d++)
+        {
+          bmin[d] = is_physbnd[d*2+0] ? 0 : imin[d];
+          bmax[d] = is_physbnd[d*2+1] ? cctk_lsh[d] : imax[d];
+
+          /* Now restrict to only the boundary points on the current face */
+          switch(dir[d])
+          {
+          case -1:
+            bmin[d] = 0;
+            bmax[d] = imin[d];
+            here_is_physbnd = here_is_physbnd || is_physbnd[2*d+0];
+            break;
+          case 0:
+            /* do nothing */
+            break;
+          case +1:
+            bmin[d] = imax[d];
+            bmax[d] = cctk_lsh[d];
+            here_is_physbnd = here_is_physbnd || is_physbnd[2*d+1];
+            break;
+          }
+
+          /* Choose a basis */
+          normal[d] = dir[d];
+          tangentA[d] = dir[d]; // FIXME
+          tangentB[d] = dir[d]; // FIXME
+        }
+
+        /* Ensure the normal vector is normalized */
+        CCTK_REAL normal_norm = 0;
+        for (int i = 0; i < 3; i++)
+        {
+          normal_norm += pow(normal[i],2);
+        }
+        normal_norm = sqrt(normal_norm);
+
+        if (fabs(normal_norm) > 1e-10)
+        {
+          for (int i = 0; i < 3; i++)
+          {
+            normal[i] /= normal_norm;
+          }
+        }
+
+        if (here_is_physbnd)
+        {
+          calc(cctkGH, old_dir, old_face, normal, tangentA, tangentB, bmin, bmax, 0, NULL);
+        }
 
-      basis_on_boundary(dir, face, normal, tangent1, tangent2);
-
-      if (is_physbnd[dir * 2 + face])
-      {
-        calc(cctkGH, dir, face, normal, tangent1, tangent2, bmin, bmax, 0, NULL);
       }
     }
   }
@@ -243,19 +264,19 @@ void GenericFD_LoopOverInterior(cGH *cctkGH, Kranc_Calculation calc)
 {
   DECLARE_CCTK_ARGUMENTS
 
-  CCTK_INT   dir = 0;
-  CCTK_INT   face = 0;
   CCTK_REAL  normal[] = {0,0,0};
-  CCTK_REAL  tangent1[] = {0,0,0};
-  CCTK_REAL  tangent2[] = {0,0,0};
+  CCTK_REAL  tangentA[] = {0,0,0};
+  CCTK_REAL  tangentB[] = {0,0,0};
 
   CCTK_INT   is_symbnd[6], is_physbnd[6], is_ipbnd[6];
   CCTK_INT   imin[3], imax[3];
+  int        dir = 0;
+  int        face = 0;
 
   GenericFD_GetBoundaryInfo(cctkGH, cctk_lsh, cctk_bbox, cctk_nghostzones, 
                             imin, imax, is_symbnd, is_physbnd, is_ipbnd);
 
-  calc(cctkGH, dir, face, normal, tangent1, tangent2, imin, imax, 0, NULL);
+  calc(cctkGH, dir, face, normal, tangentA, tangentB, imin, imax, 0, NULL);
   
   return;
 }
@@ -274,12 +295,14 @@ void GenericFD_PenaltyPrim2Char(cGH *cctkGH, CCTK_INT const dir,
   DECLARE_CCTK_ARGUMENTS
 
   CCTK_REAL  normal[] = {0,0,0};
-  CCTK_REAL  tangent1[] = {0,0,0};
-  CCTK_REAL  tangent2[] = {0,0,0};
+  CCTK_REAL  tangentA[] = {0,0,0};
+  CCTK_REAL  tangentB[] = {0,0,0};
   CCTK_INT   bmin[] = {0,0,0};
   CCTK_INT   bmax[] = {cctk_lsh[0], cctk_lsh[1], cctk_lsh[2]};
   CCTK_REAL  **all_vars;
   int        i = 0;
+  int        dir = 0;
+  int        face = 0;
 
   all_vars = malloc(num_modes*2*sizeof(CCTK_REAL *));
   assert(all_vars != NULL);
@@ -292,11 +315,11 @@ void GenericFD_PenaltyPrim2Char(cGH *cctkGH, CCTK_INT const dir,
 
   for (int d=0; d<3; ++d) {
     normal[d] = base[d];        /* A covector, index down */
-    tangent1[d] = base[d+3];    /* A vector, index up */
-    tangent2[d] = base[d+6];    /* A vector, index up */
+    tangentA[d] = base[d+3];    /* A vector, index up */
+    tangentB[d] = base[d+6];    /* A vector, index up */
   }
 
-  calc(cctkGH, dir, face, normal, tangent1, tangent2, bmin, bmax, num_modes * 2, all_vars);
+  calc(cctkGH, dir, face, normal, tangentA, tangentB, bmin, bmax, num_modes * 2, all_vars);
 
   free(all_vars);
   
diff --git a/Auxiliary/Cactus/KrancNumericalTools/GenericFD/src/GenericFD.h b/Auxiliary/Cactus/KrancNumericalTools/GenericFD/src/GenericFD.h
index cb88b46..2bd3aa9 100644
--- a/Auxiliary/Cactus/KrancNumericalTools/GenericFD/src/GenericFD.h
+++ b/Auxiliary/Cactus/KrancNumericalTools/GenericFD/src/GenericFD.h
@@ -66,7 +66,7 @@
 #endif
 
 #if defined(KRANC_C)
-#define IMAX(int1, int2) (int1 > int2) ? int1 : int2
+#define IMAX(int1, int2) ((int1) > (int2) ? (int1) : (int2))
 #endif
 
 
@@ -78,6 +78,18 @@
 /* add method argument to shorthands */
 /*                                   */
 
+/* third derivatives */
+#define D111x(gf)  string(D111,gf)
+#define D211x(gf)  string(D211,gf)
+#define D311x(gf)  string(D311,gf)
+#define D221x(gf)  string(D221,gf)
+#define D321x(gf)  string(D321,gf)
+#define D331x(gf)  string(D331,gf)
+#define D222x(gf)  string(D222,gf)
+#define D322x(gf)  string(D322,gf)
+#define D332x(gf)  string(D332,gf)
+#define D333x(gf)  string(D333,gf)
+
 /* second derivatives */
 #define D11x(gf)  string(D11,gf)
 #define D22x(gf)  string(D22,gf)
@@ -93,6 +105,18 @@
 
 #ifdef PRECOMPUTE
 
+/* third derivatives */
+#define D111(gf)  string(D111,gf)
+#define D211(gf)  string(D211,gf)
+#define D311(gf)  string(D311,gf)
+#define D221(gf)  string(D221,gf)
+#define D321(gf)  string(D321,gf)
+#define D331(gf)  string(D331,gf)
+#define D222(gf)  string(D222,gf)
+#define D322(gf)  string(D322,gf)
+#define D332(gf)  string(D332,gf)
+#define D333(gf)  string(D333,gf)
+
 /* second derivatives */
 #define D11(gf,i,j,k)  string(D11,gf)
 #define D22(gf,i,j,k)  string(D22,gf)
@@ -108,6 +132,18 @@
 
 #else
 
+/* third derivatives */
+#define D111(gf)  D111gf(gf,i,j,k)
+#define D211(gf)  D211gf(gf,i,j,k)
+#define D311(gf)  D311gf(gf,i,j,k)
+#define D221(gf)  D221gf(gf,i,j,k)
+#define D321(gf)  D321gf(gf,i,j,k)
+#define D331(gf)  D331gf(gf,i,j,k)
+#define D222(gf)  D222gf(gf,i,j,k)
+#define D322(gf)  D322gf(gf,i,j,k)
+#define D332(gf)  D332gf(gf,i,j,k)
+#define D333(gf)  D333gf(gf,i,j,k)
+
 /* second derivatives */
 #define D11(gf,i,j,k)  D11gf(gf,i,j,k)
 #define D22(gf,i,j,k)  D22gf(gf,i,j,k)
@@ -125,6 +161,18 @@
 
 
 #ifdef FD_C0
+/* third derivatives */
+#define D111gf(gf,i,j,k)  D111_c0(gf,i,j,k)
+#define D211gf(gf,i,j,k)  D211_c0(gf,i,j,k)
+#define D311gf(gf,i,j,k)  D311_c0(gf,i,j,k)
+#define D221gf(gf,i,j,k)  D221_c0(gf,i,j,k)
+#define D321gf(gf,i,j,k)  D321_c0(gf,i,j,k)
+#define D331gf(gf,i,j,k)  D331_c0(gf,i,j,k)
+#define D222gf(gf,i,j,k)  D222_c0(gf,i,j,k)
+#define D322gf(gf,i,j,k)  D322_c0(gf,i,j,k)
+#define D332gf(gf,i,j,k)  D332_c0(gf,i,j,k)
+#define D333gf(gf,i,j,k)  D333_c0(gf,i,j,k)
+
 /* second derivatives */
 #define D11gf(gf,i,j,k)  D11_c0(gf,i,j,k)
 #define D22gf(gf,i,j,k)  D22_c0(gf,i,j,k)
@@ -142,6 +190,18 @@
 
 
 #ifdef FD_C2
+/* third derivatives */
+#define D111gf(gf,i,j,k)  D111_c2(gf,i,j,k)
+#define D211gf(gf,i,j,k)  D211_c2(gf,i,j,k)
+#define D311gf(gf,i,j,k)  D311_c2(gf,i,j,k)
+#define D221gf(gf,i,j,k)  D221_c2(gf,i,j,k)
+#define D321gf(gf,i,j,k)  D321_c2(gf,i,j,k)
+#define D331gf(gf,i,j,k)  D331_c2(gf,i,j,k)
+#define D222gf(gf,i,j,k)  D222_c2(gf,i,j,k)
+#define D322gf(gf,i,j,k)  D322_c2(gf,i,j,k)
+#define D332gf(gf,i,j,k)  D332_c2(gf,i,j,k)
+#define D333gf(gf,i,j,k)  D333_c2(gf,i,j,k)
+
 /* second derivatives */
 #define D11gf(gf,i,j,k)  D11_c2(gf,i,j,k)
 #define D22gf(gf,i,j,k)  D22_c2(gf,i,j,k)
@@ -159,6 +219,18 @@
 
 
 #ifdef FD_C4
+/* third derivatives */
+#define D111gf(gf,i,j,k)  D111_c4(gf,i,j,k)
+#define D211gf(gf,i,j,k)  D211_c4(gf,i,j,k)
+#define D311gf(gf,i,j,k)  D311_c4(gf,i,j,k)
+#define D221gf(gf,i,j,k)  D221_c4(gf,i,j,k)
+#define D321gf(gf,i,j,k)  D321_c4(gf,i,j,k)
+#define D331gf(gf,i,j,k)  D331_c4(gf,i,j,k)
+#define D222gf(gf,i,j,k)  D222_c4(gf,i,j,k)
+#define D322gf(gf,i,j,k)  D322_c4(gf,i,j,k)
+#define D332gf(gf,i,j,k)  D332_c4(gf,i,j,k)
+#define D333gf(gf,i,j,k)  D333_c4(gf,i,j,k)
+
 /* second derivatives */
 #define D11gf(gf,i,j,k)  D11_c4(gf,i,j,k) 
 #define D22gf(gf,i,j,k)  D22_c4(gf,i,j,k)
@@ -175,6 +247,18 @@
 
 
 #ifdef FD_C2C4
+/* third derivatives */
+#define D111gf(gf,i,j,k)  D111_c2c4(gf,i,j,k)
+#define D211gf(gf,i,j,k)  D211_c2c4(gf,i,j,k)
+#define D311gf(gf,i,j,k)  D311_c2c4(gf,i,j,k)
+#define D221gf(gf,i,j,k)  D221_c2c4(gf,i,j,k)
+#define D321gf(gf,i,j,k)  D321_c2c4(gf,i,j,k)
+#define D331gf(gf,i,j,k)  D331_c2c4(gf,i,j,k)
+#define D222gf(gf,i,j,k)  D222_c2c4(gf,i,j,k)
+#define D322gf(gf,i,j,k)  D322_c2c4(gf,i,j,k)
+#define D332gf(gf,i,j,k)  D332_c2c4(gf,i,j,k)
+#define D333gf(gf,i,j,k)  D333_c2c4(gf,i,j,k)
+
 /* second derivatives */
 #define D11gf(gf,i,j,k)  D11_c2c4(gf,i,j,k)
 #define D22gf(gf,i,j,k)  D22_c2c4(gf,i,j,k)
@@ -203,6 +287,19 @@
 /*  set all derivatives = 0         */
 /*  for debugging and benchmarking  */
 
+/* third derivatives */
+
+#define D111_c0(gf,i,j,k)   0.
+#define D211_c0(gf,i,j,k)   0.
+#define D311_c0(gf,i,j,k)   0.
+#define D221_c0(gf,i,j,k)   0.
+#define D321_c0(gf,i,j,k)   0.
+#define D331_c0(gf,i,j,k)   0.
+#define D222_c0(gf,i,j,k)   0.
+#define D322_c0(gf,i,j,k)   0.
+#define D332_c0(gf,i,j,k)   0.
+#define D333_c0(gf,i,j,k)   0.
+
 /*  second derivatives              */
 
 #define D11_c0(gf,i,j,k)   0.
@@ -219,13 +316,29 @@
 #define D3_c0(gf,i,j,k)    0.
 
 
+
+#ifndef KRANC_C
+
 /*  c2                  */
 /*                      */
 /*  2nd order centered  */
 /*                      */
+
+/* third derivatives, centered, 2nd order */
+
+#define D111_c2(gf,i,j,k)   ((- gf(i+2,j,k) + 2*gf(i+1,j,k) - 2*gf(i-1,j,k) + gf(i-2,j,k)) * dxi*dxi*dxi * (1.0/2.0))
+#define D211_c2(gf,i,j,k)   ((gf(i+1,j+1,k) - 2*gf(i,j+1,k) + gf(i-1,j+1,k) - gf(i+1,j-1,k) + 2*gf(i,j-1,k) - gf(i-1,j-1,k)) * dxi*dxi*dyi * (1.0/2.0))
+#define D311_c2(gf,i,j,k)   ((gf(i+1,j,k+1) - 2*gf(i,j,k+1) + gf(i-1,j,k+1) - gf(i+1,j,k-1) + 2*gf(i,j,k-1) - gf(i-1,j,k-1)) * dxi*dxi*dzi * (1.0/2.0))
+#define D221_c2(gf,i,j,k)   ((gf(i+1,j+1,k) - 2*gf(i+1,j,k) + gf(i+1,j-1,k) - gf(i-1,j+1,k) + 2*gf(i-1,j,k) - gf(i-1,j-1,k)) * dxi*dyi*dyi * (1.0/2.0))
+#define D321_c2(gf,i,j,k)   ((gf(i+1,j+1,k+1) - gf(i-1,j+1,k+1) - gf(i+1,j-1,k+1) + gf(i-1,j-1,k+1) - gf(i+1,j+1,k-1) + gf(i-1,j+1,k-1) + gf(i+1,j-1,k-1) - gf(i-1,j-1,k-1)) * dxi*dyi*dzi * (1.0/8.0))
+#define D331_c2(gf,i,j,k)   ((gf(i+1,j,k+1) - 2*gf(i+1,j,k) + gf(i+1,j,k-1) - gf(i-1,j,k+1) + 2*gf(i-1,j,k) - gf(i-1,j,k-1)) * dxi*dzi*dzi * (1.0/2.0))
+#define D222_c2(gf,i,j,k)   ((- gf(i,j+2,k) + 2*gf(i,j+1,k) - 2*gf(i,j-1,k) + gf(i,j-2,k)) * dyi*dyi*dyi * (1.0/2.0))
+#define D322_c2(gf,i,j,k)   ((gf(i,j+1,k+1) - 2*gf(i,j,k+1) + gf(i,j-1,k+1) - gf(i,j+1,k-1) + 2*gf(i,j,k-1) - gf(i,j-1,k-1)) * dyi*dyi*dzi * (1.0/2.0))
+#define D332_c2(gf,i,j,k)   ((gf(i,j+1,k+1) - 2*gf(i,j+1,k) + gf(i,j+1,k-1) - gf(i,j-1,k+1) + 2*gf(i,j-1,k) - gf(i,j-1,k-1)) * dyi*dzi*dzi * (1.0/2.0))
+#define D333_c2(gf,i,j,k)   ((- gf(i,j,k+2) + 2*gf(i,j,k+1) - 2*gf(i,j,k-1) + gf(i,j,k-2)) * dzi*dzi*dzi * (1.0/2.0))
+
 /* second derivatives, centered, 2nd order */
 
-#ifndef KRANC_C
 #define D11_c2(gf,i,j,k)                         \
 	 ((   gf(i+1,j,k) \
 	 - 2.*gf(i,  j,k) \
@@ -613,8 +726,8 @@ static inline CCTK_REAL sbp_deriv_z(int i, int j, int k,
 
 typedef void(*Kranc_Calculation)(cGH *cctkGH, CCTK_INT dir, CCTK_INT   face,
                                  CCTK_REAL  normal[3],
-                                 CCTK_REAL  tangent1[3],
-                                 CCTK_REAL  tangent2[3],
+                                 CCTK_REAL  tangentA[3],
+                                 CCTK_REAL  tangentB[3],
                                  CCTK_INT   min[3],
                                  CCTK_INT   max[3], 
                                  CCTK_INT   n_subblock_gfs,