Add 8th order FD implementation.

Disabled for now.

1 files changed, 103 insertions, 50 deletions

diff --git a/src/qms.c b/src/qms.c
index d66a614..0634929 100644
--- a/src/qms.c
+++ b/src/qms.c
@@ -497,6 +497,59 @@ static double diff4_dxz(const double *arr, ptrdiff_t stride, double dx_inv, doub
             +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;
 }
 
+#define FD8(arr, stride, h_inv) \
+    ((-3.0 * arr[4 * stride] + 32.0 * arr[3 * stride] - 168.0 * arr[2 * stride] + 672.0 * arr[1 * stride] \
+      - 672.0 * arr[-1 * stride] + 168.0 * arr[-2 * stride] - 32.0 * arr[-3 * stride] + 3.0 * arr[-4 * stride]) * h_inv / 840.0)
+#define FD28(arr, stride, h_inv) \
+    ((-9.0 * arr[ 4 * stride] + 128.0 * arr[ 3 * stride] - 1008.0 * arr[ 2 * stride] + 8064.0 * arr[ 1 * stride]    \
+      -9.0 * arr[-4 * stride] + 128.0 * arr[-3 * stride] - 1008.0 * arr[-2 * stride] + 8064.0 * arr[-1 * stride]    \
+      - 14350.0 * arr[0]) * SQR(h_inv) / 5040.0)
+
+static double diff8_dx(const double *arr, ptrdiff_t stride, double h_inv)
+{
+    return FD8(arr, 1, h_inv);
+}
+static double diff8_dz(const double *arr, ptrdiff_t stride, double h_inv)
+{
+    return FD8(arr, stride, h_inv);
+}
+
+static double diff8_dxx(const double *arr, ptrdiff_t stride, double h_inv)
+{
+    return FD28(arr, 1, h_inv);
+}
+static double diff8_dzz(const double *arr, ptrdiff_t stride, double h_inv)
+{
+    return FD28(arr, stride, h_inv);
+}
+static double diff8_dxz(const double *arr, ptrdiff_t stride, double dx_inv, double dz_inv)
+{
+    static const double fd_coeffs[] = {  3.0, -32.0,   168.0, -672.0,      0.0,  672.0,  -168.0,  32.0, -3.0 };
+    static const int stencil = ARRAY_ELEMS(fd_coeffs) / 2;
+    double val = 0.0;
+    for (int j = -stencil; j <= stencil; j++) {
+        double tmp = 0.0;
+        for (int i = -stencil; i <= stencil; i++)
+            tmp += fd_coeffs[i + stencil] * arr[j * stride + i];
+        val += fd_coeffs[j + stencil] * tmp;
+    }
+    return val * dx_inv * dz_inv / SQR(840.0);
+}
+
+#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
+#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
+#endif
+
 static void fill_eq_coeffs(QMSMGContext *ctx, CoordPatch *cp, MG2DContext *solver)
 {
     cGH *gh = ctx->gh;
@@ -571,15 +624,15 @@ static void fill_eq_coeffs(QMSMGContext *ctx, CoordPatch *cp, MG2DContext *solve
                                      { gtxy, gtyy, gtyz },
                                      { gtxz, gtyz, gtzz }};
 
-            const double dx_gt11 = diff4_dx(&a_gtxx[idx_src], 1, dx_inv);
-            const double dx_gt22 = diff4_dx(&a_gtyy[idx_src], 1, dx_inv);
-            const double dx_gt33 = diff4_dx(&a_gtzz[idx_src], 1, dx_inv);
-            const double dx_gt13 = diff4_dx(&a_gtxz[idx_src], 1, dx_inv);
+            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 dz_gt11 = diff4_dz(&a_gtxx[idx_src], gh->cctk_lsh[0], dz_inv);
-            const double dz_gt22 = diff4_dz(&a_gtyy[idx_src], gh->cctk_lsh[0], dz_inv);
-            const double dz_gt33 = diff4_dz(&a_gtzz[idx_src], gh->cctk_lsh[0], dz_inv);
-            const double dz_gt13 = diff4_dz(&a_gtxz[idx_src], gh->cctk_lsh[0], dz_inv);
+            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 dgt[3][3][3] = {
                 {
@@ -599,20 +652,20 @@ static void fill_eq_coeffs(QMSMGContext *ctx, CoordPatch *cp, MG2DContext *solve
                 },
             };
 
-            const double dxx_gt11 = diff4_dxx(&a_gtxx[idx_src], 1, dx_inv);
-            const double dxx_gt22 = diff4_dxx(&a_gtyy[idx_src], 1, dx_inv);
-            const double dxx_gt33 = diff4_dxx(&a_gtzz[idx_src], 1, dx_inv);
-            const double dxx_gt13 = diff4_dxx(&a_gtxz[idx_src], 1, dx_inv);
+            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 dxz_gt11 = diff4_dxz(&a_gtxx[idx_src], gh->cctk_lsh[0], dx_inv, dz_inv);
-            const double dxz_gt22 = diff4_dxz(&a_gtyy[idx_src], gh->cctk_lsh[0], dx_inv, dz_inv);
-            const double dxz_gt33 = diff4_dxz(&a_gtzz[idx_src], gh->cctk_lsh[0], dx_inv, dz_inv);
-            const double dxz_gt13 = diff4_dxz(&a_gtxz[idx_src], gh->cctk_lsh[0], dx_inv, dz_inv);
+            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 dzz_gt11 = diff4_dzz(&a_gtxx[idx_src], gh->cctk_lsh[0], dz_inv);
-            const double dzz_gt22 = diff4_dzz(&a_gtyy[idx_src], gh->cctk_lsh[0], dz_inv);
-            const double dzz_gt33 = diff4_dzz(&a_gtzz[idx_src], gh->cctk_lsh[0], dz_inv);
-            const double dzz_gt13 = diff4_dzz(&a_gtxz[idx_src], gh->cctk_lsh[0], dz_inv);
+            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 d2gt[3][3][3][3] = {
                 {
@@ -688,15 +741,15 @@ static void fill_eq_coeffs(QMSMGContext *ctx, CoordPatch *cp, MG2DContext *solve
                 { Atxy, Atyy, Atyz },
                 { Atxz, Atyz, Atzz }};
 
-            const double dx_At11 = diff4_dx(&a_Atxx[idx_src], 1, dx_inv);
-            const double dx_At22 = diff4_dx(&a_Atyy[idx_src], 1, dx_inv);
-            const double dx_At33 = diff4_dx(&a_Atzz[idx_src], 1, dx_inv);
-            const double dx_At13 = diff4_dx(&a_Atxz[idx_src], 1, dx_inv);
+            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 dz_At11 = diff4_dz(&a_Atxx[idx_src], gh->cctk_lsh[0], dz_inv);
-            const double dz_At22 = diff4_dz(&a_Atyy[idx_src], gh->cctk_lsh[0], dz_inv);
-            const double dz_At33 = diff4_dz(&a_Atzz[idx_src], gh->cctk_lsh[0], dz_inv);
-            const double dz_At13 = diff4_dz(&a_Atxz[idx_src], gh->cctk_lsh[0], dz_inv);
+            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 dAt[3][3][3] = {
                 {
@@ -718,14 +771,14 @@ static void fill_eq_coeffs(QMSMGContext *ctx, CoordPatch *cp, MG2DContext *solve
 
             const double phi    = a_phi[idx_src];
 
-            const double dx_phi = diff4_dx(&a_phi[idx_src], 1, dx_inv);
-            const double dz_phi = diff4_dz(&a_phi[idx_src], gh->cctk_lsh[0], dz_inv);
+            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 dxx_phi = diff4_dxx(&a_phi[idx_src], gh->cctk_lsh[0], dx_inv);
-            const double dzz_phi = diff4_dzz(&a_phi[idx_src], gh->cctk_lsh[0], dz_inv);
-            const double dxz_phi = diff4_dxz(&a_phi[idx_src], gh->cctk_lsh[0], dx_inv, dz_inv);
+            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 d2phi[3][3] = {
                 { dxx_phi, 0.0, dxz_phi },
@@ -735,8 +788,8 @@ static void fill_eq_coeffs(QMSMGContext *ctx, CoordPatch *cp, MG2DContext *solve
 
             const double trK     = a_trK[idx_src];
 
-            const double dx_trk = diff4_dx(&a_trK[idx_src], 1, dx_inv);
-            const double dz_trk = diff4_dz(&a_trK[idx_src], gh->cctk_lsh[0], dz_inv);
+            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 dtrK[3] = { dx_trk, 0.0, dz_trk };
 
@@ -744,14 +797,14 @@ static void fill_eq_coeffs(QMSMGContext *ctx, CoordPatch *cp, MG2DContext *solve
             const double  Xtz  = a_Xtz[idx_src];
 
             const double alpha     = a_alpha[idx_src];
-            const double dx_alpha  = diff4_dx(&a_alpha[idx_src], 1, dx_inv);
-            const double dz_alpha  = diff4_dz(&a_alpha[idx_src], gh->cctk_lsh[0], dz_inv);
+            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 dalpha[3] = { dx_alpha, 0.0, dz_alpha };
 
-            const double dxx_alpha  = diff4_dxx(&a_alpha[idx_src], 1, dx_inv);
-            const double dzz_alpha  = diff4_dzz(&a_alpha[idx_src], gh->cctk_lsh[0], dz_inv);
-            const double dxz_alpha  = diff4_dxz(&a_alpha[idx_src], gh->cctk_lsh[0], dx_inv, dz_inv);
+            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 d2alpha[3][3] = {{ dxx_alpha,                                0, dxz_alpha },
                                           {         0, on_axis ? dxx_alpha : dx_alpha / x,         0 },
@@ -762,23 +815,23 @@ static void fill_eq_coeffs(QMSMGContext *ctx, CoordPatch *cp, MG2DContext *solve
 
             const double beta[3] = { betax, 0.0, betaz };
 
-            const double dx_betax  = diff4_dx(&a_betax[idx_src], 1, dx_inv);
-            const double dz_betax  = diff4_dz(&a_betax[idx_src], gh->cctk_lsh[0], dz_inv);
+            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_betaz  = diff4_dx(&a_betaz[idx_src], 1, dx_inv);
-            const double dz_betaz  = diff4_dz(&a_betaz[idx_src], gh->cctk_lsh[0], dz_inv);
+            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 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 = diff4_dxx(&a_betax[idx_src], gh->cctk_lsh[0], dx_inv);
-            const double dxz_betax = diff4_dxz(&a_betax[idx_src], gh->cctk_lsh[0], dx_inv, dz_inv);
-            const double dzz_betax = diff4_dzz(&a_betax[idx_src], gh->cctk_lsh[0], dz_inv);
+            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_betaz = diff4_dxx(&a_betaz[idx_src], gh->cctk_lsh[0], dx_inv);
-            const double dxz_betaz = diff4_dxz(&a_betaz[idx_src], gh->cctk_lsh[0], dx_inv, dz_inv);
-            const double dzz_betaz = diff4_dzz(&a_betaz[idx_src], gh->cctk_lsh[0], dz_inv);
+            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);
 
             const double d2beta[3][3][3] = {
                 {