Fix and enable 4th order phi derivatives for 4th order solver.

1 files changed, 13 insertions, 11 deletions

diff --git a/src/maximal_slicing_axi_mg.c b/src/maximal_slicing_axi_mg.c
index 14bd2c8..aa7b2e1 100644
--- a/src/maximal_slicing_axi_mg.c
+++ b/src/maximal_slicing_axi_mg.c
@@ -361,8 +361,9 @@ static void context_free(MSMGContext **pms)
     *pms = NULL;
 }
 
-static void fill_eq_coeffs(cGH *gh, CoordPatch *cp)
+static void fill_eq_coeffs(MSMGContext *ctx, CoordPatch *cp)
 {
+    cGH *gh = ctx->gh;
     int ret;
 
     const ptrdiff_t stride_z = CCTK_GFINDEX3D(gh, 0, 0, 1);
@@ -413,14 +414,6 @@ static void fill_eq_coeffs(cGH *gh, CoordPatch *cp)
             const double  Xtx =  a_Xtx[idx_src];
             const double  Xtz =  a_Xtz[idx_src];
 
-#if 1
-            const double phi_dx = (a_phi[idx_src + 1] - a_phi[idx_src - 1]) / (2.0 * dx);
-            const double phi_dz = (a_phi[idx_src + stride_z] - a_phi[idx_src - stride_z]) / (2.0 * dz);
-#else
-            const double phi_dx = (-1.0 * a_phi[idx_src + 2] + 8.0 * a_phi[idx_src + 1] - 8.0 * a_phi[idx_src - 1] + a_phi[idx_src - 1]) / (12.0 * dx);
-            const double phi_dz = (-1.0 * a_phi[idx_src + 2 * stride_z] + 8.0 * a_phi[idx_src + 1 * stride_z] - 8.0 * a_phi[idx_src - 1 * stride_z] + a_phi[idx_src - 1 * stride_z]) / (12.0 * dz);
-#endif
-
             const double det = gtxx * gtyy * gtzz + 2 * gtxy * gtyz * gtxz - gtzz * SQR(gtxy) - SQR(gtxz) * gtyy - gtxx * SQR(gtyz);
 
             const double At[3][3] = {
@@ -430,9 +423,18 @@ static void fill_eq_coeffs(cGH *gh, CoordPatch *cp)
 
 
             double Xx, Xz, k2;
+            double phi_dx, phi_dz;
 
             double Am[3][3], gtu[3][3];
 
+            if (ctx->fd_stencil == 1) {
+                phi_dx = (a_phi[idx_src + 1] - a_phi[idx_src - 1]) / (2.0 * dx);
+                phi_dz = (a_phi[idx_src + stride_z] - a_phi[idx_src - stride_z]) / (2.0 * dz);
+            } else {
+                phi_dx = (-1.0 * a_phi[idx_src + 2] + 8.0 * a_phi[idx_src + 1] - 8.0 * a_phi[idx_src - 1] + a_phi[idx_src - 2]) / (12.0 * dx);
+                phi_dz = (-1.0 * a_phi[idx_src + 2 * stride_z] + 8.0 * a_phi[idx_src + 1 * stride_z] - 8.0 * a_phi[idx_src - 1 * stride_z] + a_phi[idx_src - 2 * stride_z]) / (12.0 * dz);
+            }
+
 			gtu[0][0] =  (gtyy * gtzz - SQR(gtyz)) / det;
 			gtu[1][1] =  (gtxx * gtzz - SQR(gtxz)) / det;
 			gtu[2][2] =  (gtxx * gtyy - SQR(gtxy)) / det;
@@ -565,7 +567,7 @@ void msa_mg_eval(CCTK_ARGUMENTS)
         CoordPatch *cp = get_coord_patch(ms, reflevel);
 
         start = gettime();
-        fill_eq_coeffs(ms->gh, cp);
+        fill_eq_coeffs(ms, cp);
         ms->time_fine_fill_coeffs += gettime() - start;
 
         LOGDEBUG( "extrapolating BCs from t1=%g (%g) and t0=%g (%g)\n", lapse_prev1_time[reflevel], fact1, lapse_prev0_time[reflevel], fact0);
@@ -661,7 +663,7 @@ void msa_mg_solve(CCTK_ARGUMENTS)
     /* fill in the equation coefficients */
     cp = get_coord_patch(ms, reflevel);
     start = gettime();
-    fill_eq_coeffs(cctkGH, cp);
+    fill_eq_coeffs(ms, cp);
     ms->time_solve_fill_coeffs += gettime() - start;
 
     start = gettime();