fill_eq_coeffs(): refactor loading variables/calculating derivatives

Prepare the code for templatization/vectorization.

1 files changed, 202 insertions, 139 deletions

diff --git a/src/qms.c b/src/qms.c
index ab803c2..30e1c73 100644
--- a/src/qms.c
+++ b/src/qms.c
@@ -494,34 +494,54 @@ static void print_stats(QMSMGContext *ms)
     ms->log_level = orig_log_level;
 }
 
-#define FD4(arr, stride, h_inv) \
-    ((-1.0 * arr[2 * stride] + 8.0 * arr[1 * stride] - 8.0 * arr[-1 * stride] + 1.0 * arr[-2 * stride]) * h_inv / 12.0)
-#define FD24(arr, stride, h_inv) \
-    ((-1.0 * arr[2 * stride] + 16.0 * arr[1 * stride] - 30.0 * arr[0] + 16.0 * arr[-1 * stride] - 1.0 * arr[-2 * stride]) * SQR(h_inv) / 12.0)
+#define LOAD(arr, idx) (arr[idx])
 
-static double diff4_dx(const double *arr, ptrdiff_t stride, double h_inv)
+static inline double fd1_4(const double *arr, const ptrdiff_t stride, const double h_inv)
 {
-    return FD4(arr, 1, h_inv);
-}
-static double diff4_dz(const double *arr, ptrdiff_t stride, double h_inv)
-{
-    return FD4(arr, stride, h_inv);
+    const double ap2 = LOAD(arr,  2 * stride);
+    const double ap1 = LOAD(arr,  1 * stride);
+    const double am1 = LOAD(arr, -1 * stride);
+    const double am2 = LOAD(arr, -2 * stride);
+    return (-ap2 + 8.0 * ap1 - 8.0 * am1 + am2) * h_inv / 12.0;
 }
 
-static double diff4_dxx(const double *arr, ptrdiff_t stride, double h_inv)
-{
-    return FD24(arr, 1, h_inv);
-}
-static double diff4_dzz(const double *arr, ptrdiff_t stride, double h_inv)
+static inline double fd2_4(const double *arr, const ptrdiff_t stride, const double h_inv)
 {
-    return FD24(arr, stride, h_inv);
+    const double ap2 = LOAD(arr,  2 * stride);
+    const double ap1 = LOAD(arr,  1 * stride);
+    const double a0  = LOAD(arr,  0 * stride);
+    const double am1 = LOAD(arr, -1 * stride);
+    const double am2 = LOAD(arr, -2 * stride);
+    return (-ap2 + 16.0 * ap1 - 30.0 * a0 + 16.0 * am1 - am2) * SQR(h_inv) / 12.0;
 }
-static double diff4_dxz(const double *arr, ptrdiff_t stride, double dx_inv, double dz_inv)
+
+static inline double fd11_4(const double *arr, ptrdiff_t stride_x, ptrdiff_t stride_z,
+                            double dx_inv, double dz_inv)
 {
-    return (-1.0 * (-1.0 * arr[ 2 * stride + 2] + 8.0 * arr[ 2 * stride + 1] - 8.0 * arr[ 2 * stride - 1] + 1.0 * arr[ 2 * stride - 2])
-            +8.0 * (-1.0 * arr[ 1 * stride + 2] + 8.0 * arr[ 1 * stride + 1] - 8.0 * arr[ 1 * stride - 1] + 1.0 * arr[ 1 * stride - 2])
-            -8.0 * (-1.0 * arr[-1 * stride + 2] + 8.0 * arr[-1 * stride + 1] - 8.0 * arr[-1 * stride - 1] + 1.0 * arr[-1 * stride - 2])
-            +1.0 * (-1.0 * arr[-2 * stride + 2] + 8.0 * arr[-2 * stride + 1] - 8.0 * arr[-2 * stride - 1] + 1.0 * arr[-2 * stride - 2])) * dx_inv * dz_inv / 144.0;
+    const double ap2p2 = LOAD(arr,  2 * stride_x + 2 * stride_z);
+    const double ap2p1 = LOAD(arr,  2 * stride_x + 1 * stride_z);
+    const double ap2m1 = LOAD(arr,  2 * stride_x - 1 * stride_z);
+    const double ap2m2 = LOAD(arr,  2 * stride_x - 2 * stride_z);
+
+    const double ap1p2 = LOAD(arr,  1 * stride_x + 2 * stride_z);
+    const double ap1p1 = LOAD(arr,  1 * stride_x + 1 * stride_z);
+    const double ap1m1 = LOAD(arr,  1 * stride_x - 1 * stride_z);
+    const double ap1m2 = LOAD(arr,  1 * stride_x - 2 * stride_z);
+
+    const double am1p2 = LOAD(arr, -1 * stride_x + 2 * stride_z);
+    const double am1p1 = LOAD(arr, -1 * stride_x + 1 * stride_z);
+    const double am1m1 = LOAD(arr, -1 * stride_x - 1 * stride_z);
+    const double am1m2 = LOAD(arr, -1 * stride_x - 2 * stride_z);
+
+    const double am2p2 = LOAD(arr, -2 * stride_x + 2 * stride_z);
+    const double am2p1 = LOAD(arr, -2 * stride_x + 1 * stride_z);
+    const double am2m1 = LOAD(arr, -2 * stride_x - 1 * stride_z);
+    const double am2m2 = LOAD(arr, -2 * stride_x - 2 * stride_z);
+
+    return (-1.0 * (-1.0 * ap2p2 + 8.0 * ap1p2 - 8.0 * am1p2 + 1.0 * am2p2)
+            +8.0 * (-1.0 * ap2p1 + 8.0 * ap1p1 - 8.0 * am1p1 + 1.0 * am2p1)
+            -8.0 * (-1.0 * ap2m1 + 8.0 * ap1m1 - 8.0 * am1m1 + 1.0 * am2m1)
+            +1.0 * (-1.0 * ap2m2 + 8.0 * ap1m2 - 8.0 * am1m2 + 1.0 * am2m2)) * dx_inv * dz_inv / 144.0;
 }
 
 #define FD8(arr, stride, h_inv) \
@@ -564,19 +584,55 @@ static double diff8_dxz(const double *arr, ptrdiff_t stride, double dx_inv, doub
 }
 
 #if 0
-#define diff_dx diff8_dx
-#define diff_dz diff8_dz
-#define diff_dxx diff8_dxx
-#define diff_dzz diff8_dzz
-#define diff_dxz diff8_dxz
+#define diff1   fd1_8
+#define diff2   fd2_8
+#define diff11 fd11_8
 #else
-#define diff_dx diff4_dx
-#define diff_dz diff4_dz
-#define diff_dxx diff4_dxx
-#define diff_dzz diff4_dzz
-#define diff_dxz diff4_dxz
+#define diff1   fd1_4
+#define diff2   fd2_4
+#define diff11 fd11_4
 #endif
 
+#define diff_dx(arr, stride, dx_inv, dz_inv)    (diff1 (arr,      1,         dx_inv))
+#define diff_dz(arr, stride, dx_inv, dz_inv)    (diff1 (arr, stride,         dz_inv))
+#define diff_dxx(arr, stride, dx_inv, dz_inv)   (diff2 (arr,      1,         dx_inv))
+#define diff_dzz(arr, stride, dx_inv, dz_inv)   (diff2 (arr, stride,         dz_inv))
+#define diff_dxz(arr, stride, dx_inv, dz_inv)   (diff11(arr,      1, stride, dx_inv, dz_inv))
+
+enum RHSVar {
+    RHS_VAR_GTXX,
+    RHS_VAR_GTYY,
+    RHS_VAR_GTZZ,
+    RHS_VAR_GTXY,
+    RHS_VAR_GTXZ,
+    RHS_VAR_GTYZ,
+
+    RHS_VAR_ATXX,
+    RHS_VAR_ATYY,
+    RHS_VAR_ATZZ,
+    RHS_VAR_ATXY,
+    RHS_VAR_ATXZ,
+    RHS_VAR_ATYZ,
+
+    RHS_VAR_BETAX,
+    RHS_VAR_BETAY,
+    RHS_VAR_BETAZ,
+
+    RHS_VAR_XTX,
+    RHS_VAR_XTY,
+    RHS_VAR_XTZ,
+
+    RHS_VAR_TRK,
+    RHS_VAR_PHI,
+    RHS_VAR_ALPHA,
+
+    RHS_VAR_X,
+    RHS_VAR_Y,
+    RHS_VAR_Z,
+
+    RHS_VAR_NB,
+};
+
 static void fill_eq_coeffs(QMSMGContext *ctx, CoordPatch *cp, MG2DContext *solver)
 {
     cGH *gh = ctx->gh;
@@ -584,35 +640,36 @@ static void fill_eq_coeffs(QMSMGContext *ctx, CoordPatch *cp, MG2DContext *solve
 
     const ptrdiff_t stride_z = CCTK_GFINDEX3D(gh, 0, 0, 1);
 
-    double *a_x = CCTK_VarDataPtr(gh, 0, "grid::x");
-    double *a_z = CCTK_VarDataPtr(gh, 0, "grid::z");
-
-    const double dx = a_x[1] - a_x[0];
-    const double dz = a_z[stride_z] - a_z[0];
-    const double dx_inv = 1.0 / dx;
-    const double dz_inv = 1.0 / dz;
-
-    double *a_gtxx   = CCTK_VarDataPtr(gh, 0, "ML_BSSN::gt11");
-    double *a_gtyy   = CCTK_VarDataPtr(gh, 0, "ML_BSSN::gt22");
-    double *a_gtzz   = CCTK_VarDataPtr(gh, 0, "ML_BSSN::gt33");
-    double *a_gtxy   = CCTK_VarDataPtr(gh, 0, "ML_BSSN::gt12");
-    double *a_gtxz   = CCTK_VarDataPtr(gh, 0, "ML_BSSN::gt13");
-    double *a_gtyz   = CCTK_VarDataPtr(gh, 0, "ML_BSSN::gt23");
-    double *a_Atxx   = CCTK_VarDataPtr(gh, 0, "ML_BSSN::At11");
-    double *a_Atyy   = CCTK_VarDataPtr(gh, 0, "ML_BSSN::At22");
-    double *a_Atzz   = CCTK_VarDataPtr(gh, 0, "ML_BSSN::At33");
-    double *a_Atxy   = CCTK_VarDataPtr(gh, 0, "ML_BSSN::At12");
-    double *a_Atxz   = CCTK_VarDataPtr(gh, 0, "ML_BSSN::At13");
-    double *a_Atyz   = CCTK_VarDataPtr(gh, 0, "ML_BSSN::At23");
-    double *a_trK    = CCTK_VarDataPtr(gh, 0, "ML_BSSN::trK");
-    double  *a_phi   = CCTK_VarDataPtr(gh, 0, "ML_BSSN::phi");
-    double  *a_Xtx   = CCTK_VarDataPtr(gh, 0, "ML_BSSN::Xt1");
-    double  *a_Xtz   = CCTK_VarDataPtr(gh, 0, "ML_BSSN::Xt3");
-
-    double  *a_alpha = CCTK_VarDataPtr(gh, 0, "ML_BSSN::alpha");
-
-    double  *a_betax = CCTK_VarDataPtr(gh, 0, "ML_BSSN::beta1");
-    double  *a_betaz = CCTK_VarDataPtr(gh, 0, "ML_BSSN::beta3");
+    const double * const vars[] = {
+        [RHS_VAR_GTXX]   = CCTK_VarDataPtr(gh, 0, "ML_BSSN::gt11"),
+        [RHS_VAR_GTYY]   = CCTK_VarDataPtr(gh, 0, "ML_BSSN::gt22"),
+        [RHS_VAR_GTZZ]   = CCTK_VarDataPtr(gh, 0, "ML_BSSN::gt33"),
+        [RHS_VAR_GTXY]   = CCTK_VarDataPtr(gh, 0, "ML_BSSN::gt12"),
+        [RHS_VAR_GTXZ]   = CCTK_VarDataPtr(gh, 0, "ML_BSSN::gt13"),
+        [RHS_VAR_GTYZ]   = CCTK_VarDataPtr(gh, 0, "ML_BSSN::gt23"),
+
+        [RHS_VAR_ATXX]   = CCTK_VarDataPtr(gh, 0, "ML_BSSN::At11"),
+        [RHS_VAR_ATYY]   = CCTK_VarDataPtr(gh, 0, "ML_BSSN::At22"),
+        [RHS_VAR_ATZZ]   = CCTK_VarDataPtr(gh, 0, "ML_BSSN::At33"),
+        [RHS_VAR_ATXY]   = CCTK_VarDataPtr(gh, 0, "ML_BSSN::At12"),
+        [RHS_VAR_ATXZ]   = CCTK_VarDataPtr(gh, 0, "ML_BSSN::At13"),
+        [RHS_VAR_ATYZ]   = CCTK_VarDataPtr(gh, 0, "ML_BSSN::At23"),
+
+        [RHS_VAR_TRK]    = CCTK_VarDataPtr(gh, 0, "ML_BSSN::trK"),
+        [RHS_VAR_PHI]    = CCTK_VarDataPtr(gh, 0, "ML_BSSN::phi"),
+        [RHS_VAR_ALPHA]  = CCTK_VarDataPtr(gh, 0, "ML_BSSN::alpha"),
+
+        [RHS_VAR_XTX]    = CCTK_VarDataPtr(gh, 0, "ML_BSSN::Xt1"),
+        [RHS_VAR_XTZ]    = CCTK_VarDataPtr(gh, 0, "ML_BSSN::Xt3"),
+
+        [RHS_VAR_BETAX]  = CCTK_VarDataPtr(gh, 0, "ML_BSSN::beta1"),
+        [RHS_VAR_BETAZ]  = CCTK_VarDataPtr(gh, 0, "ML_BSSN::beta3"),
+        [RHS_VAR_X]      = CCTK_VarDataPtr(gh, 0, "grid::x"),
+        [RHS_VAR_Z]      = CCTK_VarDataPtr(gh, 0, "grid::z"),
+    };
+
+    const double dx_inv = 1.0 / (vars[RHS_VAR_X][1]        - vars[RHS_VAR_X][0]);
+    const double dz_inv = 1.0 / (vars[RHS_VAR_Z][stride_z] - vars[RHS_VAR_Z][0]);
 
     solver->diff_coeffs[MG2D_DIFF_COEFF_20]->boundaries[MG2D_BOUNDARY_0L].flags = MG2D_DC_FLAG_DISCONT;
     solver->diff_coeffs[MG2D_DIFF_COEFF_10]->boundaries[MG2D_BOUNDARY_0L].flags = MG2D_DC_FLAG_POLE;
@@ -624,7 +681,7 @@ static void fill_eq_coeffs(QMSMGContext *ctx, CoordPatch *cp, MG2DContext *solve
             const int idx_dc  = idx_z * solver->diff_coeffs[0]->stride + idx_x;
             const int idx_rhs = idx_z * solver->rhs_stride + idx_x;
 
-            const double x = a_x[idx_src];
+            const double x = vars[RHS_VAR_X][idx_src];
             const int on_axis = fabs(x) < 1e-13;
 
             double K[3][3], Km[3][3], Ku[3][3], dK[3][3][3];
@@ -640,26 +697,33 @@ static void fill_eq_coeffs(QMSMGContext *ctx, CoordPatch *cp, MG2DContext *solve
 
             double betal[3], dbetal[3][3], d2betal[3][3][3];
 
-            const double gtxx  = a_gtxx[idx_src];
-            const double gtyy  = a_gtyy[idx_src];
-            const double gtzz  = a_gtzz[idx_src];
-            const double gtxy  = a_gtxy[idx_src];
-            const double gtxz  = a_gtxz[idx_src];
-            const double gtyz  = a_gtyz[idx_src];
+#define LOAD_VAR(var)  (        vars[RHS_VAR_ ## var][idx_src])
+#define DX(var)        (diff1 (&vars[RHS_VAR_ ## var][idx_src],        1,           dx_inv        ))
+#define DZ(var)        (diff1 (&vars[RHS_VAR_ ## var][idx_src],           stride_z,         dz_inv))
+#define DXX(var)       (diff2 (&vars[RHS_VAR_ ## var][idx_src],        1,           dx_inv))
+#define DZZ(var)       (diff2 (&vars[RHS_VAR_ ## var][idx_src],           stride_z,         dz_inv))
+#define DXZ(var)       (diff11(&vars[RHS_VAR_ ## var][idx_src],        1, stride_z, dx_inv, dz_inv))
+
+            const double gtxx  = LOAD_VAR(GTXX);
+            const double gtyy  = LOAD_VAR(GTYY);
+            const double gtzz  = LOAD_VAR(GTZZ);
+            const double gtxy  = LOAD_VAR(GTXY);
+            const double gtxz  = LOAD_VAR(GTXZ);
+            const double gtyz  = LOAD_VAR(GTYZ);
 
             const double gt[3][3] = {{ gtxx, gtxy, gtxz },
                                      { gtxy, gtyy, gtyz },
                                      { gtxz, gtyz, gtzz }};
 
-            const double dx_gt11 = diff_dx(&a_gtxx[idx_src], 1, dx_inv);
-            const double dx_gt22 = diff_dx(&a_gtyy[idx_src], 1, dx_inv);
-            const double dx_gt33 = diff_dx(&a_gtzz[idx_src], 1, dx_inv);
-            const double dx_gt13 = diff_dx(&a_gtxz[idx_src], 1, dx_inv);
+            const double dx_gt11 = DX(GTXX);
+            const double dx_gt22 = DX(GTYY);
+            const double dx_gt33 = DX(GTZZ);
+            const double dx_gt13 = DX(GTXZ);
 
-            const double dz_gt11 = diff_dz(&a_gtxx[idx_src], gh->cctk_lsh[0], dz_inv);
-            const double dz_gt22 = diff_dz(&a_gtyy[idx_src], gh->cctk_lsh[0], dz_inv);
-            const double dz_gt33 = diff_dz(&a_gtzz[idx_src], gh->cctk_lsh[0], dz_inv);
-            const double dz_gt13 = diff_dz(&a_gtxz[idx_src], gh->cctk_lsh[0], dz_inv);
+            const double dz_gt11 = DZ(GTXX);
+            const double dz_gt22 = DZ(GTYY);
+            const double dz_gt33 = DZ(GTZZ);
+            const double dz_gt13 = DZ(GTXZ);
 
             const double dgt[3][3][3] = {
                 {
@@ -679,20 +743,20 @@ static void fill_eq_coeffs(QMSMGContext *ctx, CoordPatch *cp, MG2DContext *solve
                 },
             };
 
-            const double dxx_gt11 = diff_dxx(&a_gtxx[idx_src], 1, dx_inv);
-            const double dxx_gt22 = diff_dxx(&a_gtyy[idx_src], 1, dx_inv);
-            const double dxx_gt33 = diff_dxx(&a_gtzz[idx_src], 1, dx_inv);
-            const double dxx_gt13 = diff_dxx(&a_gtxz[idx_src], 1, dx_inv);
+            const double dxx_gt11 = DXX(GTXX);
+            const double dxx_gt22 = DXX(GTYY);
+            const double dxx_gt33 = DXX(GTZZ);
+            const double dxx_gt13 = DXX(GTXZ);
 
-            const double dxz_gt11 = diff_dxz(&a_gtxx[idx_src], gh->cctk_lsh[0], dx_inv, dz_inv);
-            const double dxz_gt22 = diff_dxz(&a_gtyy[idx_src], gh->cctk_lsh[0], dx_inv, dz_inv);
-            const double dxz_gt33 = diff_dxz(&a_gtzz[idx_src], gh->cctk_lsh[0], dx_inv, dz_inv);
-            const double dxz_gt13 = diff_dxz(&a_gtxz[idx_src], gh->cctk_lsh[0], dx_inv, dz_inv);
+            const double dxz_gt11 = DXZ(GTXX);
+            const double dxz_gt22 = DXZ(GTYY);
+            const double dxz_gt33 = DXZ(GTZZ);
+            const double dxz_gt13 = DXZ(GTXZ);
 
-            const double dzz_gt11 = diff_dzz(&a_gtxx[idx_src], gh->cctk_lsh[0], dz_inv);
-            const double dzz_gt22 = diff_dzz(&a_gtyy[idx_src], gh->cctk_lsh[0], dz_inv);
-            const double dzz_gt33 = diff_dzz(&a_gtzz[idx_src], gh->cctk_lsh[0], dz_inv);
-            const double dzz_gt13 = diff_dzz(&a_gtxz[idx_src], gh->cctk_lsh[0], dz_inv);
+            const double dzz_gt11 = DZZ(GTXX);
+            const double dzz_gt22 = DZZ(GTYY);
+            const double dzz_gt33 = DZZ(GTZZ);
+            const double dzz_gt13 = DZZ(GTXZ);
 
             const double d2gt[3][3][3][3] = {
                 {
@@ -756,27 +820,27 @@ static void fill_eq_coeffs(QMSMGContext *ctx, CoordPatch *cp, MG2DContext *solve
                 },
             };
 
-            const double Atxx  = a_Atxx[idx_src];
-            const double Atyy  = a_Atyy[idx_src];
-            const double Atzz  = a_Atzz[idx_src];
-            const double Atxy  = a_Atxy[idx_src];
-            const double Atxz  = a_Atxz[idx_src];
-            const double Atyz  = a_Atyz[idx_src];
+            const double Atxx  = LOAD_VAR(ATXX);
+            const double Atyy  = LOAD_VAR(ATYY);
+            const double Atzz  = LOAD_VAR(ATZZ);
+            const double Atxy  = LOAD_VAR(ATXY);
+            const double Atxz  = LOAD_VAR(ATXZ);
+            const double Atyz  = LOAD_VAR(ATYZ);
 
             const double At[3][3] = {
                 { Atxx, Atxy, Atxz },
                 { Atxy, Atyy, Atyz },
                 { Atxz, Atyz, Atzz }};
 
-            const double dx_At11 = diff_dx(&a_Atxx[idx_src], 1, dx_inv);
-            const double dx_At22 = diff_dx(&a_Atyy[idx_src], 1, dx_inv);
-            const double dx_At33 = diff_dx(&a_Atzz[idx_src], 1, dx_inv);
-            const double dx_At13 = diff_dx(&a_Atxz[idx_src], 1, dx_inv);
+            const double dx_At11 = DX(ATXX);
+            const double dx_At22 = DX(ATYY);
+            const double dx_At33 = DX(ATZZ);
+            const double dx_At13 = DX(ATXZ);
 
-            const double dz_At11 = diff_dz(&a_Atxx[idx_src], gh->cctk_lsh[0], dz_inv);
-            const double dz_At22 = diff_dz(&a_Atyy[idx_src], gh->cctk_lsh[0], dz_inv);
-            const double dz_At33 = diff_dz(&a_Atzz[idx_src], gh->cctk_lsh[0], dz_inv);
-            const double dz_At13 = diff_dz(&a_Atxz[idx_src], gh->cctk_lsh[0], dz_inv);
+            const double dz_At11 = DZ(ATXX);
+            const double dz_At22 = DZ(ATYY);
+            const double dz_At33 = DZ(ATZZ);
+            const double dz_At13 = DZ(ATXZ);
 
             const double dAt[3][3][3] = {
                 {
@@ -796,69 +860,68 @@ static void fill_eq_coeffs(QMSMGContext *ctx, CoordPatch *cp, MG2DContext *solve
                 },
             };
 
-            const double phi    = a_phi[idx_src];
-
-            const double dx_phi = diff_dx(&a_phi[idx_src], 1, dx_inv);
-            const double dz_phi = diff_dz(&a_phi[idx_src], gh->cctk_lsh[0], dz_inv);
-
-            const double dphi[3] = { dx_phi, 0.0, dz_phi };
+            const double phi     = LOAD_VAR(PHI);
+            const double dphi[3] = { DX(PHI), 0.0, DZ(PHI) };
 
-            const double dxx_phi = diff_dxx(&a_phi[idx_src], gh->cctk_lsh[0], dx_inv);
-            const double dzz_phi = diff_dzz(&a_phi[idx_src], gh->cctk_lsh[0], dz_inv);
-            const double dxz_phi = diff_dxz(&a_phi[idx_src], gh->cctk_lsh[0], dx_inv, dz_inv);
+            const double dxx_phi = DXX(PHI);
+            const double dzz_phi = DZZ(PHI);
+            const double dxz_phi = DXZ(PHI);
 
             const double d2phi[3][3] = {
-                { dxx_phi, 0.0, dxz_phi },
-                {     0.0, on_axis ? dxx_phi : dx_phi / x, 0.0 },
-                { dxz_phi,    0.0,  dzz_phi },
+                { dxx_phi,                             0.0, dxz_phi },
+                {     0.0, on_axis ? dxx_phi : dphi[0] / x,     0.0 },
+                { dxz_phi,                             0.0, dzz_phi },
             };
 
-            const double trK     = a_trK[idx_src];
+            const double trK     = LOAD_VAR(TRK);
+            const double dtrK[3] = { DX(TRK), 0.0, DZ(TRK) };
 
-            const double dx_trk = diff_dx(&a_trK[idx_src], 1, dx_inv);
-            const double dz_trk = diff_dz(&a_trK[idx_src], gh->cctk_lsh[0], dz_inv);
+            const double  Xtx  = LOAD_VAR(XTX);
+            const double  Xtz  = LOAD_VAR(XTZ);
 
-            const double dtrK[3] = { dx_trk, 0.0, dz_trk };
-
-            const double  Xtx  = a_Xtx[idx_src];
-            const double  Xtz  = a_Xtz[idx_src];
-
-            const double alpha     = a_alpha[idx_src];
-            const double dx_alpha  = diff_dx(&a_alpha[idx_src], 1, dx_inv);
-            const double dz_alpha  = diff_dz(&a_alpha[idx_src], gh->cctk_lsh[0], dz_inv);
+            const double alpha     = LOAD_VAR(ALPHA);
+            const double dx_alpha  = DX(ALPHA);
+            const double dz_alpha  = DZ(ALPHA);
 
             const double dalpha[3] = { dx_alpha, 0.0, dz_alpha };
 
-            const double dxx_alpha  = diff_dxx(&a_alpha[idx_src], 1, dx_inv);
-            const double dzz_alpha  = diff_dzz(&a_alpha[idx_src], gh->cctk_lsh[0], dz_inv);
-            const double dxz_alpha  = diff_dxz(&a_alpha[idx_src], gh->cctk_lsh[0], dx_inv, dz_inv);
+            const double dxx_alpha  = DXX(ALPHA);
+            const double dzz_alpha  = DZZ(ALPHA);
+            const double dxz_alpha  = DXZ(ALPHA);
 
             const double d2alpha[3][3] = {{ dxx_alpha,                                0, dxz_alpha },
                                           {         0, on_axis ? dxx_alpha : dx_alpha / x,         0 },
                                           { dxz_alpha,                                0, dzz_alpha }};
 
-            const double betax  = a_betax[idx_src];
-            const double betaz  = a_betaz[idx_src];
+            const double betax  = LOAD_VAR(BETAX);
+            const double betaz  = LOAD_VAR(BETAX);
 
             const double beta[3] = { betax, 0.0, betaz };
 
-            const double dx_betax  = diff_dx(&a_betax[idx_src], 1, dx_inv);
-            const double dz_betax  = diff_dz(&a_betax[idx_src], gh->cctk_lsh[0], dz_inv);
+            const double dx_betax  = DX(BETAX);
+            const double dz_betax  = DZ(BETAX);
 
-            const double dx_betaz  = diff_dx(&a_betaz[idx_src], 1, dx_inv);
-            const double dz_betaz  = diff_dz(&a_betaz[idx_src], gh->cctk_lsh[0], dz_inv);
+            const double dx_betaz  = DX(BETAZ);
+            const double dz_betaz  = DZ(BETAZ);
 
             const double dbeta[3][3] = {{ dx_betax,      0.0, dx_betaz },
                                         {      0.0, on_axis ? dx_betax : betax / x,      0.0 },
                                         { dz_betax,      0.0, dz_betaz }};
 
-            const double dxx_betax = diff_dxx(&a_betax[idx_src], gh->cctk_lsh[0], dx_inv);
-            const double dxz_betax = diff_dxz(&a_betax[idx_src], gh->cctk_lsh[0], dx_inv, dz_inv);
-            const double dzz_betax = diff_dzz(&a_betax[idx_src], gh->cctk_lsh[0], dz_inv);
+            const double dxx_betax = DXX(BETAX);
+            const double dxz_betax = DXZ(BETAX);
+            const double dzz_betax = DZZ(BETAX);
+
+            const double dxx_betaz = DXX(BETAZ);
+            const double dxz_betaz = DXZ(BETAZ);
+            const double dzz_betaz = DZZ(BETAZ);
 
-            const double dxx_betaz = diff_dxx(&a_betaz[idx_src], gh->cctk_lsh[0], dx_inv);
-            const double dxz_betaz = diff_dxz(&a_betaz[idx_src], gh->cctk_lsh[0], dx_inv, dz_inv);
-            const double dzz_betaz = diff_dzz(&a_betaz[idx_src], gh->cctk_lsh[0], dz_inv);
+#undef LOAD_VAR
+#undef DX
+#undef DZ
+#undef DXX
+#undef DZZ
+#undef DXZ
 
             const double d2beta[3][3][3] = {
                 {