Introduce cctk_ash, retire cctk_lssh

Introduce cctk_ash, describing the process-local array shape that has
been allocated. This may be larger than cctk_lsh, the process-local
shape that should be used.

Retire cctk_lssh and related infrastructure to handle staggered grid
functions.

6 files changed, 94 insertions, 70 deletions

diff --git a/ML_WaveToy/src/WT_Dirichlet.cc b/ML_WaveToy/src/WT_Dirichlet.cc
index daf0a6a..4629fda 100644
--- a/ML_WaveToy/src/WT_Dirichlet.cc
+++ b/ML_WaveToy/src/WT_Dirichlet.cc
@@ -18,10 +18,14 @@
 
 /* Define macros used in calculations */
 #define INITVALUE (42)
-#define QAD(x) (SQR(SQR(x)))
+#define ScalarINV(x) ((CCTK_REAL)1.0 / (x))
+#define ScalarSQR(x) ((x) * (x))
+#define ScalarCUB(x) ((x) * ScalarSQR(x))
+#define ScalarQAD(x) (ScalarSQR(ScalarSQR(x)))
 #define INV(x) (kdiv(ToReal(1.0),x))
 #define SQR(x) (kmul(x,x))
 #define CUB(x) (kmul(x,SQR(x)))
+#define QAD(x) (SQR(SQR(x)))
 
 extern "C" void WT_Dirichlet_SelectBCs(CCTK_ARGUMENTS)
 {
@@ -71,15 +75,15 @@ static void WT_Dirichlet_Body(cGH const * restrict const cctkGH, int const dir,
   CCTK_REAL_VEC const hdzi = kmul(ToReal(0.5), dzi);
   
   /* Initialize predefined quantities */
-  CCTK_REAL_VEC const p1o12dx = kmul(INV(dx),ToReal(0.0833333333333333333333333333333));
-  CCTK_REAL_VEC const p1o12dy = kmul(INV(dy),ToReal(0.0833333333333333333333333333333));
-  CCTK_REAL_VEC const p1o12dz = kmul(INV(dz),ToReal(0.0833333333333333333333333333333));
-  CCTK_REAL_VEC const p1o144dxdy = kmul(INV(kmul(dx,dy)),ToReal(0.00694444444444444444444444444444));
-  CCTK_REAL_VEC const p1o144dxdz = kmul(INV(kmul(dx,dz)),ToReal(0.00694444444444444444444444444444));
-  CCTK_REAL_VEC const p1o144dydz = kmul(INV(kmul(dy,dz)),ToReal(0.00694444444444444444444444444444));
-  CCTK_REAL_VEC const pm1o12dx2 = kmul(INV(SQR(dx)),ToReal(-0.0833333333333333333333333333333));
-  CCTK_REAL_VEC const pm1o12dy2 = kmul(INV(SQR(dy)),ToReal(-0.0833333333333333333333333333333));
-  CCTK_REAL_VEC const pm1o12dz2 = kmul(INV(SQR(dz)),ToReal(-0.0833333333333333333333333333333));
+  CCTK_REAL_VEC const p1o12dx = kdiv(ToReal(0.0833333333333333333333333333333),dx);
+  CCTK_REAL_VEC const p1o12dy = kdiv(ToReal(0.0833333333333333333333333333333),dy);
+  CCTK_REAL_VEC const p1o12dz = kdiv(ToReal(0.0833333333333333333333333333333),dz);
+  CCTK_REAL_VEC const p1o144dxdy = kdiv(ToReal(0.00694444444444444444444444444444),kmul(dy,dx));
+  CCTK_REAL_VEC const p1o144dxdz = kdiv(ToReal(0.00694444444444444444444444444444),kmul(dz,dx));
+  CCTK_REAL_VEC const p1o144dydz = kdiv(ToReal(0.00694444444444444444444444444444),kmul(dz,dy));
+  CCTK_REAL_VEC const pm1o12dx2 = kdiv(ToReal(-0.0833333333333333333333333333333),kmul(dx,dx));
+  CCTK_REAL_VEC const pm1o12dy2 = kdiv(ToReal(-0.0833333333333333333333333333333),kmul(dy,dy));
+  CCTK_REAL_VEC const pm1o12dz2 = kdiv(ToReal(-0.0833333333333333333333333333333),kmul(dz,dz));
   
   /* Assign local copies of arrays functions */
   
@@ -93,7 +97,7 @@ static void WT_Dirichlet_Body(cGH const * restrict const cctkGH, int const dir,
   #pragma omp parallel
   LC_LOOP3VEC(WT_Dirichlet,
     i,j,k, imin[0],imin[1],imin[2], imax[0],imax[1],imax[2],
-    cctk_lsh[0],cctk_lsh[1],cctk_lsh[2],
+    cctk_ash[0],cctk_ash[1],cctk_ash[2],
     CCTK_REAL_VEC_SIZE)
   {
     ptrdiff_t const index = di*i + dj*j + dk*k;
diff --git a/ML_WaveToy/src/WT_Energy.cc b/ML_WaveToy/src/WT_Energy.cc
index ba36b94..96355af 100644
--- a/ML_WaveToy/src/WT_Energy.cc
+++ b/ML_WaveToy/src/WT_Energy.cc
@@ -18,10 +18,14 @@
 
 /* Define macros used in calculations */
 #define INITVALUE (42)
-#define QAD(x) (SQR(SQR(x)))
+#define ScalarINV(x) ((CCTK_REAL)1.0 / (x))
+#define ScalarSQR(x) ((x) * (x))
+#define ScalarCUB(x) ((x) * ScalarSQR(x))
+#define ScalarQAD(x) (ScalarSQR(ScalarSQR(x)))
 #define INV(x) (kdiv(ToReal(1.0),x))
 #define SQR(x) (kmul(x,x))
 #define CUB(x) (kmul(x,SQR(x)))
+#define QAD(x) (SQR(SQR(x)))
 
 extern "C" void WT_Energy_SelectBCs(CCTK_ARGUMENTS)
 {
@@ -68,15 +72,15 @@ static void WT_Energy_Body(cGH const * restrict const cctkGH, int const dir, int
   CCTK_REAL_VEC const hdzi = kmul(ToReal(0.5), dzi);
   
   /* Initialize predefined quantities */
-  CCTK_REAL_VEC const p1o12dx = kmul(INV(dx),ToReal(0.0833333333333333333333333333333));
-  CCTK_REAL_VEC const p1o12dy = kmul(INV(dy),ToReal(0.0833333333333333333333333333333));
-  CCTK_REAL_VEC const p1o12dz = kmul(INV(dz),ToReal(0.0833333333333333333333333333333));
-  CCTK_REAL_VEC const p1o144dxdy = kmul(INV(kmul(dx,dy)),ToReal(0.00694444444444444444444444444444));
-  CCTK_REAL_VEC const p1o144dxdz = kmul(INV(kmul(dx,dz)),ToReal(0.00694444444444444444444444444444));
-  CCTK_REAL_VEC const p1o144dydz = kmul(INV(kmul(dy,dz)),ToReal(0.00694444444444444444444444444444));
-  CCTK_REAL_VEC const pm1o12dx2 = kmul(INV(SQR(dx)),ToReal(-0.0833333333333333333333333333333));
-  CCTK_REAL_VEC const pm1o12dy2 = kmul(INV(SQR(dy)),ToReal(-0.0833333333333333333333333333333));
-  CCTK_REAL_VEC const pm1o12dz2 = kmul(INV(SQR(dz)),ToReal(-0.0833333333333333333333333333333));
+  CCTK_REAL_VEC const p1o12dx = kdiv(ToReal(0.0833333333333333333333333333333),dx);
+  CCTK_REAL_VEC const p1o12dy = kdiv(ToReal(0.0833333333333333333333333333333),dy);
+  CCTK_REAL_VEC const p1o12dz = kdiv(ToReal(0.0833333333333333333333333333333),dz);
+  CCTK_REAL_VEC const p1o144dxdy = kdiv(ToReal(0.00694444444444444444444444444444),kmul(dy,dx));
+  CCTK_REAL_VEC const p1o144dxdz = kdiv(ToReal(0.00694444444444444444444444444444),kmul(dz,dx));
+  CCTK_REAL_VEC const p1o144dydz = kdiv(ToReal(0.00694444444444444444444444444444),kmul(dz,dy));
+  CCTK_REAL_VEC const pm1o12dx2 = kdiv(ToReal(-0.0833333333333333333333333333333),kmul(dx,dx));
+  CCTK_REAL_VEC const pm1o12dy2 = kdiv(ToReal(-0.0833333333333333333333333333333),kmul(dy,dy));
+  CCTK_REAL_VEC const pm1o12dz2 = kdiv(ToReal(-0.0833333333333333333333333333333),kmul(dz,dz));
   
   /* Assign local copies of arrays functions */
   
@@ -90,7 +94,7 @@ static void WT_Energy_Body(cGH const * restrict const cctkGH, int const dir, int
   #pragma omp parallel
   LC_LOOP3VEC(WT_Energy,
     i,j,k, imin[0],imin[1],imin[2], imax[0],imax[1],imax[2],
-    cctk_lsh[0],cctk_lsh[1],cctk_lsh[2],
+    cctk_ash[0],cctk_ash[1],cctk_ash[2],
     CCTK_REAL_VEC_SIZE)
   {
     ptrdiff_t const index = di*i + dj*j + dk*k;
@@ -110,7 +114,7 @@ static void WT_Energy_Body(cGH const * restrict const cctkGH, int const dir, int
     
     /* Calculate temporaries and grid functions */
     CCTK_REAL_VEC epsL = 
-      kmul(kadd(SQR(rhoL),kadd(SQR(PDstandardNth1u),kadd(SQR(PDstandardNth2u),SQR(PDstandardNth3u)))),ToReal(0.5));
+      kmul(kmadd(rhoL,rhoL,kmadd(PDstandardNth1u,PDstandardNth1u,kmadd(PDstandardNth2u,PDstandardNth2u,kmul(PDstandardNth3u,PDstandardNth3u)))),ToReal(0.5));
     
     /* Copy local copies back to grid functions */
     vec_store_partial_prepare(i,lc_imin,lc_imax);
diff --git a/ML_WaveToy/src/WT_EnergyBoundary.cc b/ML_WaveToy/src/WT_EnergyBoundary.cc
index be8978d..a1672f4 100644
--- a/ML_WaveToy/src/WT_EnergyBoundary.cc
+++ b/ML_WaveToy/src/WT_EnergyBoundary.cc
@@ -18,10 +18,14 @@
 
 /* Define macros used in calculations */
 #define INITVALUE (42)
-#define QAD(x) (SQR(SQR(x)))
+#define ScalarINV(x) ((CCTK_REAL)1.0 / (x))
+#define ScalarSQR(x) ((x) * (x))
+#define ScalarCUB(x) ((x) * ScalarSQR(x))
+#define ScalarQAD(x) (ScalarSQR(ScalarSQR(x)))
 #define INV(x) (kdiv(ToReal(1.0),x))
 #define SQR(x) (kmul(x,x))
 #define CUB(x) (kmul(x,SQR(x)))
+#define QAD(x) (SQR(SQR(x)))
 
 extern "C" void WT_EnergyBoundary_SelectBCs(CCTK_ARGUMENTS)
 {
@@ -68,15 +72,15 @@ static void WT_EnergyBoundary_Body(cGH const * restrict const cctkGH, int const
   CCTK_REAL_VEC const hdzi = kmul(ToReal(0.5), dzi);
   
   /* Initialize predefined quantities */
-  CCTK_REAL_VEC const p1o12dx = kmul(INV(dx),ToReal(0.0833333333333333333333333333333));
-  CCTK_REAL_VEC const p1o12dy = kmul(INV(dy),ToReal(0.0833333333333333333333333333333));
-  CCTK_REAL_VEC const p1o12dz = kmul(INV(dz),ToReal(0.0833333333333333333333333333333));
-  CCTK_REAL_VEC const p1o144dxdy = kmul(INV(kmul(dx,dy)),ToReal(0.00694444444444444444444444444444));
-  CCTK_REAL_VEC const p1o144dxdz = kmul(INV(kmul(dx,dz)),ToReal(0.00694444444444444444444444444444));
-  CCTK_REAL_VEC const p1o144dydz = kmul(INV(kmul(dy,dz)),ToReal(0.00694444444444444444444444444444));
-  CCTK_REAL_VEC const pm1o12dx2 = kmul(INV(SQR(dx)),ToReal(-0.0833333333333333333333333333333));
-  CCTK_REAL_VEC const pm1o12dy2 = kmul(INV(SQR(dy)),ToReal(-0.0833333333333333333333333333333));
-  CCTK_REAL_VEC const pm1o12dz2 = kmul(INV(SQR(dz)),ToReal(-0.0833333333333333333333333333333));
+  CCTK_REAL_VEC const p1o12dx = kdiv(ToReal(0.0833333333333333333333333333333),dx);
+  CCTK_REAL_VEC const p1o12dy = kdiv(ToReal(0.0833333333333333333333333333333),dy);
+  CCTK_REAL_VEC const p1o12dz = kdiv(ToReal(0.0833333333333333333333333333333),dz);
+  CCTK_REAL_VEC const p1o144dxdy = kdiv(ToReal(0.00694444444444444444444444444444),kmul(dy,dx));
+  CCTK_REAL_VEC const p1o144dxdz = kdiv(ToReal(0.00694444444444444444444444444444),kmul(dz,dx));
+  CCTK_REAL_VEC const p1o144dydz = kdiv(ToReal(0.00694444444444444444444444444444),kmul(dz,dy));
+  CCTK_REAL_VEC const pm1o12dx2 = kdiv(ToReal(-0.0833333333333333333333333333333),kmul(dx,dx));
+  CCTK_REAL_VEC const pm1o12dy2 = kdiv(ToReal(-0.0833333333333333333333333333333),kmul(dy,dy));
+  CCTK_REAL_VEC const pm1o12dz2 = kdiv(ToReal(-0.0833333333333333333333333333333),kmul(dz,dz));
   
   /* Assign local copies of arrays functions */
   
@@ -90,7 +94,7 @@ static void WT_EnergyBoundary_Body(cGH const * restrict const cctkGH, int const
   #pragma omp parallel
   LC_LOOP3VEC(WT_EnergyBoundary,
     i,j,k, imin[0],imin[1],imin[2], imax[0],imax[1],imax[2],
-    cctk_lsh[0],cctk_lsh[1],cctk_lsh[2],
+    cctk_ash[0],cctk_ash[1],cctk_ash[2],
     CCTK_REAL_VEC_SIZE)
   {
     ptrdiff_t const index = di*i + dj*j + dk*k;
diff --git a/ML_WaveToy/src/WT_Gaussian.cc b/ML_WaveToy/src/WT_Gaussian.cc
index cde55be..2c8f46a 100644
--- a/ML_WaveToy/src/WT_Gaussian.cc
+++ b/ML_WaveToy/src/WT_Gaussian.cc
@@ -18,10 +18,14 @@
 
 /* Define macros used in calculations */
 #define INITVALUE (42)
-#define QAD(x) (SQR(SQR(x)))
+#define ScalarINV(x) ((CCTK_REAL)1.0 / (x))
+#define ScalarSQR(x) ((x) * (x))
+#define ScalarCUB(x) ((x) * ScalarSQR(x))
+#define ScalarQAD(x) (ScalarSQR(ScalarSQR(x)))
 #define INV(x) (kdiv(ToReal(1.0),x))
 #define SQR(x) (kmul(x,x))
 #define CUB(x) (kmul(x,SQR(x)))
+#define QAD(x) (SQR(SQR(x)))
 
 static void WT_Gaussian_Body(cGH const * restrict const cctkGH, int const dir, int const face, CCTK_REAL const normal[3], CCTK_REAL const tangentA[3], CCTK_REAL const tangentB[3], int const imin[3], int const imax[3], int const n_subblock_gfs, CCTK_REAL * restrict const subblock_gfs[])
 {
@@ -56,15 +60,15 @@ static void WT_Gaussian_Body(cGH const * restrict const cctkGH, int const dir, i
   CCTK_REAL_VEC const hdzi = kmul(ToReal(0.5), dzi);
   
   /* Initialize predefined quantities */
-  CCTK_REAL_VEC const p1o12dx = kmul(INV(dx),ToReal(0.0833333333333333333333333333333));
-  CCTK_REAL_VEC const p1o12dy = kmul(INV(dy),ToReal(0.0833333333333333333333333333333));
-  CCTK_REAL_VEC const p1o12dz = kmul(INV(dz),ToReal(0.0833333333333333333333333333333));
-  CCTK_REAL_VEC const p1o144dxdy = kmul(INV(kmul(dx,dy)),ToReal(0.00694444444444444444444444444444));
-  CCTK_REAL_VEC const p1o144dxdz = kmul(INV(kmul(dx,dz)),ToReal(0.00694444444444444444444444444444));
-  CCTK_REAL_VEC const p1o144dydz = kmul(INV(kmul(dy,dz)),ToReal(0.00694444444444444444444444444444));
-  CCTK_REAL_VEC const pm1o12dx2 = kmul(INV(SQR(dx)),ToReal(-0.0833333333333333333333333333333));
-  CCTK_REAL_VEC const pm1o12dy2 = kmul(INV(SQR(dy)),ToReal(-0.0833333333333333333333333333333));
-  CCTK_REAL_VEC const pm1o12dz2 = kmul(INV(SQR(dz)),ToReal(-0.0833333333333333333333333333333));
+  CCTK_REAL_VEC const p1o12dx = kdiv(ToReal(0.0833333333333333333333333333333),dx);
+  CCTK_REAL_VEC const p1o12dy = kdiv(ToReal(0.0833333333333333333333333333333),dy);
+  CCTK_REAL_VEC const p1o12dz = kdiv(ToReal(0.0833333333333333333333333333333),dz);
+  CCTK_REAL_VEC const p1o144dxdy = kdiv(ToReal(0.00694444444444444444444444444444),kmul(dy,dx));
+  CCTK_REAL_VEC const p1o144dxdz = kdiv(ToReal(0.00694444444444444444444444444444),kmul(dz,dx));
+  CCTK_REAL_VEC const p1o144dydz = kdiv(ToReal(0.00694444444444444444444444444444),kmul(dz,dy));
+  CCTK_REAL_VEC const pm1o12dx2 = kdiv(ToReal(-0.0833333333333333333333333333333),kmul(dx,dx));
+  CCTK_REAL_VEC const pm1o12dy2 = kdiv(ToReal(-0.0833333333333333333333333333333),kmul(dy,dy));
+  CCTK_REAL_VEC const pm1o12dz2 = kdiv(ToReal(-0.0833333333333333333333333333333),kmul(dz,dz));
   
   /* Assign local copies of arrays functions */
   
@@ -78,7 +82,7 @@ static void WT_Gaussian_Body(cGH const * restrict const cctkGH, int const dir, i
   #pragma omp parallel
   LC_LOOP3VEC(WT_Gaussian,
     i,j,k, imin[0],imin[1],imin[2], imax[0],imax[1],imax[2],
-    cctk_lsh[0],cctk_lsh[1],cctk_lsh[2],
+    cctk_ash[0],cctk_ash[1],cctk_ash[2],
     CCTK_REAL_VEC_SIZE)
   {
     ptrdiff_t const index = di*i + dj*j + dk*k;
@@ -94,7 +98,7 @@ static void WT_Gaussian_Body(cGH const * restrict const cctkGH, int const dir, i
     
     /* Calculate temporaries and grid functions */
     CCTK_REAL_VEC uL = 
-      kmul(kexp(kmul(INV(SQR(ToReal(width))),kmul(SQR(rL),ToReal(-0.5)))),ToReal(amplitude));
+      kmul(kexp(kmul(kmul(kmul(rL,rL),ToReal(-0.5)),ToReal(ScalarINV(ScalarSQR(width))))),ToReal(amplitude));
     
     CCTK_REAL_VEC rhoL = ToReal(0);
     
diff --git a/ML_WaveToy/src/WT_RHS.cc b/ML_WaveToy/src/WT_RHS.cc
index 42c14a8..bddf28f 100644
--- a/ML_WaveToy/src/WT_RHS.cc
+++ b/ML_WaveToy/src/WT_RHS.cc
@@ -18,10 +18,14 @@
 
 /* Define macros used in calculations */
 #define INITVALUE (42)
-#define QAD(x) (SQR(SQR(x)))
+#define ScalarINV(x) ((CCTK_REAL)1.0 / (x))
+#define ScalarSQR(x) ((x) * (x))
+#define ScalarCUB(x) ((x) * ScalarSQR(x))
+#define ScalarQAD(x) (ScalarSQR(ScalarSQR(x)))
 #define INV(x) (kdiv(ToReal(1.0),x))
 #define SQR(x) (kmul(x,x))
 #define CUB(x) (kmul(x,SQR(x)))
+#define QAD(x) (SQR(SQR(x)))
 
 extern "C" void WT_RHS_SelectBCs(CCTK_ARGUMENTS)
 {
@@ -71,15 +75,15 @@ static void WT_RHS_Body(cGH const * restrict const cctkGH, int const dir, int co
   CCTK_REAL_VEC const hdzi = kmul(ToReal(0.5), dzi);
   
   /* Initialize predefined quantities */
-  CCTK_REAL_VEC const p1o12dx = kmul(INV(dx),ToReal(0.0833333333333333333333333333333));
-  CCTK_REAL_VEC const p1o12dy = kmul(INV(dy),ToReal(0.0833333333333333333333333333333));
-  CCTK_REAL_VEC const p1o12dz = kmul(INV(dz),ToReal(0.0833333333333333333333333333333));
-  CCTK_REAL_VEC const p1o144dxdy = kmul(INV(kmul(dx,dy)),ToReal(0.00694444444444444444444444444444));
-  CCTK_REAL_VEC const p1o144dxdz = kmul(INV(kmul(dx,dz)),ToReal(0.00694444444444444444444444444444));
-  CCTK_REAL_VEC const p1o144dydz = kmul(INV(kmul(dy,dz)),ToReal(0.00694444444444444444444444444444));
-  CCTK_REAL_VEC const pm1o12dx2 = kmul(INV(SQR(dx)),ToReal(-0.0833333333333333333333333333333));
-  CCTK_REAL_VEC const pm1o12dy2 = kmul(INV(SQR(dy)),ToReal(-0.0833333333333333333333333333333));
-  CCTK_REAL_VEC const pm1o12dz2 = kmul(INV(SQR(dz)),ToReal(-0.0833333333333333333333333333333));
+  CCTK_REAL_VEC const p1o12dx = kdiv(ToReal(0.0833333333333333333333333333333),dx);
+  CCTK_REAL_VEC const p1o12dy = kdiv(ToReal(0.0833333333333333333333333333333),dy);
+  CCTK_REAL_VEC const p1o12dz = kdiv(ToReal(0.0833333333333333333333333333333),dz);
+  CCTK_REAL_VEC const p1o144dxdy = kdiv(ToReal(0.00694444444444444444444444444444),kmul(dy,dx));
+  CCTK_REAL_VEC const p1o144dxdz = kdiv(ToReal(0.00694444444444444444444444444444),kmul(dz,dx));
+  CCTK_REAL_VEC const p1o144dydz = kdiv(ToReal(0.00694444444444444444444444444444),kmul(dz,dy));
+  CCTK_REAL_VEC const pm1o12dx2 = kdiv(ToReal(-0.0833333333333333333333333333333),kmul(dx,dx));
+  CCTK_REAL_VEC const pm1o12dy2 = kdiv(ToReal(-0.0833333333333333333333333333333),kmul(dy,dy));
+  CCTK_REAL_VEC const pm1o12dz2 = kdiv(ToReal(-0.0833333333333333333333333333333),kmul(dz,dz));
   
   /* Assign local copies of arrays functions */
   
@@ -93,7 +97,7 @@ static void WT_RHS_Body(cGH const * restrict const cctkGH, int const dir, int co
   #pragma omp parallel
   LC_LOOP3VEC(WT_RHS,
     i,j,k, imin[0],imin[1],imin[2], imax[0],imax[1],imax[2],
-    cctk_lsh[0],cctk_lsh[1],cctk_lsh[2],
+    cctk_ash[0],cctk_ash[1],cctk_ash[2],
     CCTK_REAL_VEC_SIZE)
   {
     ptrdiff_t const index = di*i + dj*j + dk*k;
diff --git a/ML_WaveToy/src/WT_Standing.cc b/ML_WaveToy/src/WT_Standing.cc
index 8a67e52..a3d0fa6 100644
--- a/ML_WaveToy/src/WT_Standing.cc
+++ b/ML_WaveToy/src/WT_Standing.cc
@@ -18,10 +18,14 @@
 
 /* Define macros used in calculations */
 #define INITVALUE (42)
-#define QAD(x) (SQR(SQR(x)))
+#define ScalarINV(x) ((CCTK_REAL)1.0 / (x))
+#define ScalarSQR(x) ((x) * (x))
+#define ScalarCUB(x) ((x) * ScalarSQR(x))
+#define ScalarQAD(x) (ScalarSQR(ScalarSQR(x)))
 #define INV(x) (kdiv(ToReal(1.0),x))
 #define SQR(x) (kmul(x,x))
 #define CUB(x) (kmul(x,SQR(x)))
+#define QAD(x) (SQR(SQR(x)))
 
 static void WT_Standing_Body(cGH const * restrict const cctkGH, int const dir, int const face, CCTK_REAL const normal[3], CCTK_REAL const tangentA[3], CCTK_REAL const tangentB[3], int const imin[3], int const imax[3], int const n_subblock_gfs, CCTK_REAL * restrict const subblock_gfs[])
 {
@@ -56,15 +60,15 @@ static void WT_Standing_Body(cGH const * restrict const cctkGH, int const dir, i
   CCTK_REAL_VEC const hdzi = kmul(ToReal(0.5), dzi);
   
   /* Initialize predefined quantities */
-  CCTK_REAL_VEC const p1o12dx = kmul(INV(dx),ToReal(0.0833333333333333333333333333333));
-  CCTK_REAL_VEC const p1o12dy = kmul(INV(dy),ToReal(0.0833333333333333333333333333333));
-  CCTK_REAL_VEC const p1o12dz = kmul(INV(dz),ToReal(0.0833333333333333333333333333333));
-  CCTK_REAL_VEC const p1o144dxdy = kmul(INV(kmul(dx,dy)),ToReal(0.00694444444444444444444444444444));
-  CCTK_REAL_VEC const p1o144dxdz = kmul(INV(kmul(dx,dz)),ToReal(0.00694444444444444444444444444444));
-  CCTK_REAL_VEC const p1o144dydz = kmul(INV(kmul(dy,dz)),ToReal(0.00694444444444444444444444444444));
-  CCTK_REAL_VEC const pm1o12dx2 = kmul(INV(SQR(dx)),ToReal(-0.0833333333333333333333333333333));
-  CCTK_REAL_VEC const pm1o12dy2 = kmul(INV(SQR(dy)),ToReal(-0.0833333333333333333333333333333));
-  CCTK_REAL_VEC const pm1o12dz2 = kmul(INV(SQR(dz)),ToReal(-0.0833333333333333333333333333333));
+  CCTK_REAL_VEC const p1o12dx = kdiv(ToReal(0.0833333333333333333333333333333),dx);
+  CCTK_REAL_VEC const p1o12dy = kdiv(ToReal(0.0833333333333333333333333333333),dy);
+  CCTK_REAL_VEC const p1o12dz = kdiv(ToReal(0.0833333333333333333333333333333),dz);
+  CCTK_REAL_VEC const p1o144dxdy = kdiv(ToReal(0.00694444444444444444444444444444),kmul(dy,dx));
+  CCTK_REAL_VEC const p1o144dxdz = kdiv(ToReal(0.00694444444444444444444444444444),kmul(dz,dx));
+  CCTK_REAL_VEC const p1o144dydz = kdiv(ToReal(0.00694444444444444444444444444444),kmul(dz,dy));
+  CCTK_REAL_VEC const pm1o12dx2 = kdiv(ToReal(-0.0833333333333333333333333333333),kmul(dx,dx));
+  CCTK_REAL_VEC const pm1o12dy2 = kdiv(ToReal(-0.0833333333333333333333333333333),kmul(dy,dy));
+  CCTK_REAL_VEC const pm1o12dz2 = kdiv(ToReal(-0.0833333333333333333333333333333),kmul(dz,dz));
   
   /* Assign local copies of arrays functions */
   
@@ -78,7 +82,7 @@ static void WT_Standing_Body(cGH const * restrict const cctkGH, int const dir, i
   #pragma omp parallel
   LC_LOOP3VEC(WT_Standing,
     i,j,k, imin[0],imin[1],imin[2], imax[0],imax[1],imax[2],
-    cctk_lsh[0],cctk_lsh[1],cctk_lsh[2],
+    cctk_ash[0],cctk_ash[1],cctk_ash[2],
     CCTK_REAL_VEC_SIZE)
   {
     ptrdiff_t const index = di*i + dj*j + dk*k;
@@ -95,9 +99,9 @@ static void WT_Standing_Body(cGH const * restrict const cctkGH, int const dir, i
     /* Precompute derivatives */
     
     /* Calculate temporaries and grid functions */
-    CCTK_REAL_VEC k = kmul(INV(ToReal(width)),ToReal(Pi));
+    CCTK_REAL_VEC k = ToReal(Pi*ScalarINV(width));
     
-    CCTK_REAL_VEC omega = ksqrt(kmul(SQR(k),ToReal(3)));
+    CCTK_REAL_VEC omega = ksqrt(kmul(kmul(k,k),ToReal(3)));
     
     CCTK_REAL_VEC uL = 
       kmul(kcos(kmul(xL,k)),kmul(kcos(kmul(yL,k)),kmul(kcos(kmul(zL,k)),kmul(kcos(kmul(omega,t)),ToReal(amplitude)))));