Initial commit.

6 files changed, 1283 insertions, 0 deletions

diff --git a/configuration.ccl b/configuration.ccl
new file mode 100644
index 0000000..5fad107
--- /dev/null
+++ b/configuration.ccl
@@ -0,0 +1,3 @@
+# Configuration definition for thorn QuasiMaximalSlicingMG
+REQUIRES MPI
+
diff --git a/interface.ccl b/interface.ccl
new file mode 100644
index 0000000..37e039a
--- /dev/null
+++ b/interface.ccl
@@ -0,0 +1,26 @@
+# Interface definition for thorn QuasiMaximalSlicingMG
+implements: QuasiMaximalSlicingMG
+
+INHERITS: ADMBase grid CoordBase MethodOfLines ML_BSSN
+
+CCTK_INT FUNCTION MoLRegisterConstrained(CCTK_INT IN idx)
+CCTK_INT FUNCTION MoLRegisterSaveAndRestore(CCTK_INT IN idx)
+CCTK_INT FUNCTION MoLRegisterSaveAndRestoreGroup(CCTK_INT IN idx)
+
+REQUIRES FUNCTION MoLRegisterConstrained
+REQUIRES FUNCTION MoLRegisterSaveAndRestore
+REQUIRES FUNCTION MoLRegisterSaveAndRestoreGroup
+
+CCTK_INT FUNCTION MoLNumIntegratorSubsteps ()
+USES FUNCTION MoLNumIntegratorSubsteps
+
+public:
+REAL W_val0    TYPE=GF TIMELEVELS=1
+REAL W_val1    TYPE=GF TIMELEVELS=1
+REAL W_val0_time TYPE=ARRAY DIM=1 SIZE=32 DISTRIB=constant
+REAL W_val1_time TYPE=ARRAY DIM=1 SIZE=32 DISTRIB=constant
+
+REAL W_pred0 TYPE=GF tags='Prolongation="None"' TIMELEVELS=1
+REAL W_pred1 TYPE=GF tags='Prolongation="None"' TIMELEVELS=1
+REAL W_pred0_time TYPE=ARRAY DIM=1 SIZE=32 DISTRIB=constant
+REAL W_pred1_time TYPE=ARRAY DIM=1 SIZE=32 DISTRIB=constant
diff --git a/param.ccl b/param.ccl
new file mode 100644
index 0000000..a738832
--- /dev/null
+++ b/param.ccl
@@ -0,0 +1,83 @@
+# Parameter definitions for thorn QuasiMaximalSlicing
+
+shares: ML_BSSN
+EXTENDS CCTK_KEYWORD lapse_source "lapse_source"
+{
+  QMS_MG :: ""
+}
+
+PRIVATE:
+CCTK_INT solve_level "" STEERABLE=recover
+{
+    0: :: ""
+} 0
+
+CCTK_REAL switchoff_time "" STEERABLE=recover
+{
+   : :: ""
+} 9999.0
+
+CCTK_INT ccz4 "Use CCZ4 or BSSN" STEERABLE=recover
+{
+    0:1 :: ""
+} 0
+
+RESTRICTED:
+CCTK_INT fd_stencil "finite differencing stencil"
+{
+    1: :: ""
+} 1
+
+RESTRICTED:
+CCTK_INT maxiter "maximum number of solver iterations" STEERABLE=always
+{
+    1: :: ""
+} 64
+
+RESTRICTED:
+CCTK_INT nb_cycles "number of coarse-correction cycles per level" STEERABLE=always
+{
+    1: :: ""
+} 1
+
+RESTRICTED:
+CCTK_INT nb_relax_pre "number of relaxation steps before coarse-grid correction" STEERABLE=always
+{
+    0: :: ""
+} 2
+
+RESTRICTED:
+CCTK_INT nb_relax_post "number of relaxation steps after coarse-grid correction" STEERABLE=always
+{
+    0: :: ""
+} 2
+
+RESTRICTED:
+CCTK_INT exact_size "" STEERABLE=always
+{
+    0: :: ""
+} 5
+
+RESTRICTED:
+CCTK_REAL tol_residual_base "maximum absolute value of the residual for dx=1.0, scaled by 1/(dx ** 2)" STEERABLE=always
+{
+    0: :: ""
+} 0.0
+
+RESTRICTED:
+CCTK_REAL cfl_factor "" STEERABLE=always
+{
+    0: :: ""
+} 0.0
+
+RESTRICTED:
+CCTK_INT stats_every "print elliptic solver stats every <count> coarsest-level steps" STEERABLE=always
+{
+    0: :: ""
+} 0
+
+RESTRICTED:
+STRING loglevel "logging level" STEERABLE=always
+{
+    :: ""
+} "info"
diff --git a/schedule.ccl b/schedule.ccl
new file mode 100644
index 0000000..70e90c0
--- /dev/null
+++ b/schedule.ccl
@@ -0,0 +1,52 @@
+# Schedule definitions for thorn MaximalSlicingAxi
+#
+if (CCTK_EQUALS(lapse_source, "QMS_MG")) {
+    SCHEDULE msa_mg_init IN CCTK_BASEGRID AFTER TemporalSpacings {
+        LANG: C
+    } "QMS MG init"
+
+    SCHEDULE qms_mg_eval IN ML_BSSN_evolCalcGroup BEFORE ML_BSSN_lapse_evol {
+        LANG: C
+        SYNC: ML_lapse
+    } "Quasimaximal slicing eval W"
+    SCHEDULE qms_mg_eval IN ML_CCZ4_evolCalcGroup BEFORE ML_CCZ4_lapse_evol {
+        LANG: C
+        SYNC: ML_lapse
+    } "Quasimaximal slicing eval W"
+
+    #SCHEDULE quasimaximal_slicing_axi_solve IN ML_BSSN_evolCalcGroup BEFORE quasimaximal_slicing_axi_eval {
+    #SCHEDULE quasimaximal_slicing_axi_solve IN MoL_PostStep AFTER ML_BSSN_ApplyBCs {
+    SCHEDULE qms_mg_solve IN CCTK_POSTSTEP {
+        LANG: C
+        SYNC: W_pred0
+        SYNC: W_pred1
+        SYNC: W_val0
+        SYNC: W_val1
+    } "Quasimaximal slicing solve W"
+
+    SCHEDULE qms_mg_sync IN CCTK_POSTSTEP BEFORE qms_mg_solve {
+        SYNC: W_val1
+        LANG: C
+    } ""
+
+    SCHEDULE qms_mg_inithist IN CCTK_INITIAL {
+        LANG: C
+    } ""
+
+    SCHEDULE qms_mg_init IN CCTK_BASEGRID AFTER TemporalSpacings {
+        LANG: C
+    } ""
+
+    #SCHEDULE quasimaximal_slicing_axi IN MoL_PseudoEvolution {
+    #    LANG: C
+    #} "Quasimaximal slicing"
+
+    STORAGE: W_pred0
+    STORAGE: W_pred1
+    STORAGE: W_pred0_time
+    STORAGE: W_pred1_time
+    STORAGE: W_val0
+    STORAGE: W_val1
+    STORAGE: W_val0_time
+    STORAGE: W_val1_time
+}
diff --git a/src/make.code.defn b/src/make.code.defn
new file mode 100644
index 0000000..7f29e1a
--- /dev/null
+++ b/src/make.code.defn
@@ -0,0 +1,7 @@
+# Main make.code.defn file for thorn MaximalSlicingAxi
+
+# Source files in this directory
+SRCS = qms.c
+
+# Subdirectories containing source files
+SUBDIRS = 
diff --git a/src/qms.c b/src/qms.c
new file mode 100644
index 0000000..de2bb12
--- /dev/null
+++ b/src/qms.c
@@ -0,0 +1,1112 @@
+#include <ctype.h>
+#include <errno.h>
+#include <float.h>
+#include <inttypes.h>
+#include <limits.h>
+#include <math.h>
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+
+#include <mg2d.h>
+#include <mg2d_boundary.h>
+
+#include "cctk.h"
+#include "cctk_Arguments.h"
+#include "cctk_Parameters.h"
+#include "cctk_Timers.h"
+#include "util_Table.h"
+
+#define SQR(x) ((x) * (x))
+#define SGN(x) ((x) >= 0.0 ? 1.0 : -1.0)
+#define MAX(x, y) ((x) > (y) ? (x) : (y))
+#define MIN(x, y) ((x) > (y) ? (y) : (x))
+#define ABS(x) ((x >= 0) ? (x) : -(x))
+#define ARRAY_ELEMS(arr) (sizeof(arr) / sizeof(*arr))
+
+#define CPINDEX(cp, i, j, k) ((k * cp->grid_size[1] + j) * cp->grid_size[0] + i)
+
+/* precomputed values for a given refined grid */
+typedef struct CoordPatch {
+    int level;
+    ptrdiff_t grid_size[3];
+
+    MG2DContext *solver;
+
+    int y_idx;
+
+    ptrdiff_t offset_left[2];
+    int bnd_intercomp[2][2];
+} CoordPatch;
+
+typedef struct QMSMGContext {
+    cGH *gh;
+
+    int log_level;
+    int fd_stencil;
+    int maxiter;
+    int max_exact_size;
+    int nb_cycles;
+    int nb_relax_post;
+    int nb_relax_pre;
+    double tol_residual_base;
+    double cfl_factor;
+    int solve_level;
+
+    CoordPatch *patches;
+    int nb_patches;
+} QMSMGContext;
+
+static const struct {
+    const char *str;
+    enum MG2DLogLevel level;
+} log_levels[] = {
+    { "fatal",   MG2D_LOG_FATAL   },
+    { "error",   MG2D_LOG_ERROR   },
+    { "warning", MG2D_LOG_WARNING },
+    { "info",    MG2D_LOG_INFO    },
+    { "verbose", MG2D_LOG_VERBOSE },
+    { "debug",   MG2D_LOG_DEBUG   },
+};
+
+static QMSMGContext *qms_context;
+
+static int ctz(int a)
+{
+    int ret = 0;
+
+    if (!a)
+        return INT_MAX;
+
+    while (!(a & 1)) {
+        a >>= 1;
+        ret++;
+    }
+
+    return ret;
+}
+
+static void log_callback(const MG2DContext *ctx, int level,
+                         const char *fmt, va_list vl)
+{
+    QMSMGContext *ms = ctx->opaque;
+    int target_level = ms->log_level;
+    uint8_t buf[1024];
+    int ret;
+
+    if (level > target_level)
+        return;
+
+    ret = snprintf(buf, sizeof(buf), "[%d] t=%g ", ctz(ms->gh->cctk_levfac[0]), ms->gh->cctk_time);
+    if (ret >= sizeof(buf))
+        return;
+    vsnprintf(buf + ret, sizeof(buf) - ret, fmt, vl);
+    fputs(buf, stderr);
+}
+
+static CoordPatch *get_coord_patch(QMSMGContext *ms, int level)
+{
+    cGH *gh = ms->gh;
+
+    const char *omp_threads = getenv("OMP_NUM_THREADS");
+
+    const int integrator_substeps = MoLNumIntegratorSubsteps();
+
+    const double *a_x = CCTK_VarDataPtr(gh, 0, "grid::x");
+    const double *a_y = CCTK_VarDataPtr(gh, 0, "grid::y");
+    const double *a_z = CCTK_VarDataPtr(gh, 0, "grid::z");
+
+    size_t solver_size[2];
+
+    CoordPatch *cp;
+    MG2DContext *solver;
+    int i, ret;
+
+    for (int i = 0; i < ms->nb_patches; i++) {
+        cp = &ms->patches[i];
+
+        if (cp->level == level)
+            return cp;
+    }
+
+    /* create a new patch */
+    ms->patches = realloc(ms->patches, sizeof(*ms->patches) * (ms->nb_patches + 1));
+    cp = &ms->patches[ms->nb_patches];
+
+    memset(cp, 0, sizeof(*cp));
+
+    cp->level = level;
+    cp->grid_size[0] = gh->cctk_lsh[0];
+    cp->grid_size[1] = gh->cctk_lsh[1];
+    cp->grid_size[2] = gh->cctk_lsh[2];
+
+    for (i = 0; i < gh->cctk_lsh[1]; i++)
+        if (fabs(a_y[CCTK_GFINDEX3D(gh, 0, i, 0)]) < 1e-8) {
+            cp->y_idx = i;
+            break;
+        }
+    if (i == gh->cctk_lsh[1])
+        CCTK_WARN(0, "The grid does not include y==0");
+
+    for (int dir = 0; dir < 2; dir++) {
+        ptrdiff_t offset_right;
+
+        cp->offset_left[dir] = gh->cctk_nghostzones[2 * dir];
+
+        offset_right = gh->cctk_nghostzones[2 * dir];
+        if (cp->level > 0) {
+            offset_right *= integrator_substeps * 2;
+            offset_right -= 2;
+        }
+
+        solver_size[dir] = cp->grid_size[2 * dir] - cp->offset_left[dir] - offset_right;
+    }
+
+    if (solver_size[0] != solver_size[1])
+        CCTK_WARN(0, "Non-square domain");
+
+    solver = mg2d_solver_alloc(solver_size[0]);
+    if (!solver)
+        CCTK_WARN(0, "Error allocating the solver");
+
+    cp->solver = solver;
+
+    solver->step[0] = a_x[1] - a_x[0];
+    solver->step[1] = a_z[CCTK_GFINDEX3D(gh, 0, 0, 1)] - a_z[0];
+
+    solver->fd_stencil = ms->fd_stencil;
+
+    solver->boundaries[MG2D_BOUNDARY_0L]->type = MG2D_BC_TYPE_REFLECT;
+    solver->boundaries[MG2D_BOUNDARY_1L]->type = MG2D_BC_TYPE_REFLECT;
+    solver->boundaries[MG2D_BOUNDARY_0U]->type = MG2D_BC_TYPE_FIXVAL;
+    solver->boundaries[MG2D_BOUNDARY_1U]->type = MG2D_BC_TYPE_FIXVAL;
+
+    solver->maxiter        = ms->maxiter;
+    solver->tol            = ms->tol_residual_base / (solver->step[0] * solver->step[1]);
+    solver->nb_cycles      = ms->nb_cycles;
+    solver->nb_relax_pre   = ms->nb_relax_pre;
+    solver->nb_relax_post  = ms->nb_relax_post;
+    solver->cfl_factor     = ms->cfl_factor;
+    solver->max_exact_size = ms->max_exact_size;
+
+    solver->opaque       = ms;
+    solver->log_callback = log_callback;
+
+    if (omp_threads)
+        solver->nb_threads = strtol(omp_threads, NULL, 0);
+    if (solver->nb_threads <= 0)
+        solver->nb_threads = 1;
+
+    /* initialize boundary values to zero,
+     *      * non-zero values on the outer boundaries of refined levels are filled in elsewhere */
+    for (int i = 0; i < 4; i++) {
+        int ci = mg2d_bnd_coord_idx(i);
+        for (size_t j = 0; j < solver->fd_stencil; j++) {
+            memset(solver->boundaries[i]->val + j * solver->boundaries[i]->val_stride - j,
+                   0, sizeof(*solver->boundaries[i]->val) * (solver_size[!ci] + 2 * j));
+        }
+    }
+
+
+    ms->nb_patches++;
+    return cp;
+}
+
+#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)
+
+static double diff4_dx(const double *arr, ptrdiff_t stride, 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);
+}
+
+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)
+{
+    return FD24(arr, stride, h_inv);
+}
+static double diff4_dxz(const double *arr, ptrdiff_t stride, 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;
+}
+
+static void fill_eq_coeffs(QMSMGContext *ctx, CoordPatch *cp, MG2DContext *solver)
+{
+    cGH *gh = ctx->gh;
+    int ret;
+
+    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_kdotxx = CCTK_VarDataPtr(gh, 0, "ML_BSSN::Kdot11");
+    double *a_kdotyy = CCTK_VarDataPtr(gh, 0, "ML_BSSN::Kdot22");
+    double *a_kdotzz = CCTK_VarDataPtr(gh, 0, "ML_BSSN::Kdot33");
+    double *a_kdotxy = CCTK_VarDataPtr(gh, 0, "ML_BSSN::Kdot12");
+    double *a_kdotxz = CCTK_VarDataPtr(gh, 0, "ML_BSSN::Kdot13");
+    double *a_kdotyz = CCTK_VarDataPtr(gh, 0, "ML_BSSN::Kdot23");
+
+    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");
+
+    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;
+
+#pragma omp parallel for
+    for (int idx_z = 0; idx_z < solver->local_size[1]; idx_z++)
+        for (int idx_x = 0; idx_x < solver->local_size[0]; idx_x++) {
+            const int idx_src = CCTK_GFINDEX3D(gh, idx_x + cp->offset_left[0], cp->y_idx, idx_z + cp->offset_left[1]);
+            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 int on_axis = fabs(x) < 1e-13;
+
+            double K[3][3], Km[3][3], Ku[3][3], dK[3][3][3];
+            double gtu[3][3], g[3][3], gu[3][3];
+            double dg[3][3][3], dgu[3][3][3], d2g[3][3][3][3], gudot[3][3];
+            double G[3][3][3], dG[3][3][3][3], Gdot[3][3][3], X[3];
+            double Ric[3][3], Ricm[3][3];
+            double dgdot[3][3][3];
+            double k2, Kij_Dij_alpha, k_kdot, k3;
+            double Gammadot_term, beta_term;
+
+            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];
+
+            const double gt[3][3] = {{ gtxx, gtxy, gtxz },
+                                     { 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 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 dgt[3][3][3] = {
+                {
+                    { dx_gt11,     0.0, dx_gt13 },
+                    {     0.0, dx_gt22,     0.0 },
+                    { dx_gt13,     0.0, dx_gt33 },
+                },
+                {
+                    {     0.0, on_axis ? dx_gt11 - dx_gt22 : (gtxx - gtyy) / x, 0.0 },
+                    { on_axis ? dx_gt11 - dx_gt22 : (gtxx - gtyy) / x, 0.0, on_axis ? dx_gt13 : gtxz / x },
+                    { 0.0, on_axis ? dx_gt13 : gtxz / x, 0.0 },
+                },
+                {
+                    { dz_gt11,     0.0, dz_gt13 },
+                    {     0.0, dz_gt22,     0.0 },
+                    { dz_gt13,     0.0, dz_gt33 },
+                },
+            };
+
+            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 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 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 d2gt[3][3][3][3] = {
+                {
+                    {
+                        { dxx_gt11,      0.0, dxx_gt13 },
+                        { 0.0,      dxx_gt22,      0.0 },
+                        { dxx_gt13,      0.0, dxx_gt33 },
+                    },
+                    {
+                        {      0.0, on_axis ? 0.5 * (dxx_gt11 - dxx_gt22) : (dx_gt11 - dx_gt22) / x - (gtxx - gtyy) / SQR(x), 0.0 },
+                        { on_axis ? 0.5 * (dxx_gt11 - dxx_gt22) : (dx_gt11 - dx_gt22) / x - (gtxx - gtyy) / SQR(x), 0.0,
+                            on_axis ? 0.5 * dxx_gt13 : dx_gt13 / x - gtxz / SQR(x) },
+                        { 0.0, on_axis ? 0.5 * dxx_gt13 : dx_gt13 / x - gtxz / SQR(x), 0.0 },
+                    },
+                    {
+                        { dxz_gt11,     0.0, dxz_gt13 },
+                        { 0.0,     dxz_gt22,      0.0 },
+                        { dxz_gt13,     0.0, dxz_gt33 },
+                    },
+
+                },
+                {
+                    {
+                        { 0.0, on_axis ? 0.5 * (dxx_gt11 - dxx_gt22) : (dx_gt11 - dx_gt22) / x - (gtxx - gtyy) / SQR(x), 0.0 },
+                        { on_axis ? 0.5 * (dxx_gt11 - dxx_gt22) : (dx_gt11 - dx_gt22) / x - (gtxx - gtyy) / SQR(x), 0.0,
+                            on_axis ? 0.5 * dxx_gt13 : dx_gt13 / x - gtxz / SQR(x) },
+                        { 0.0, on_axis ? 0.5 * dxx_gt13 : dx_gt13 / x - gtxz / SQR(x), 0.0 },
+                    },
+                    {
+                        { on_axis ? dxx_gt22 : dx_gt11 / x - 2 * (gtxx - gtyy) / SQR(x), 0.0,
+                           on_axis ? 0.5 * dxx_gt13 : dx_gt13 / x - gtxz / SQR(x) },
+                        { 0.0, on_axis ? dxx_gt11 : dx_gt22 / x + 2.0 * (gtxx - gtyy) / SQR(x), 0.0 },
+                        { on_axis ? 0.5 * dxx_gt13 : dx_gt13 / x - gtxz / SQR(x), 0.0, on_axis ? dxx_gt33 : dx_gt33 / x },
+                    },
+                    {
+                        { 0.0, on_axis ? dxz_gt11 - dxz_gt22 : (dz_gt11 - dz_gt22) / x, 0.0 },
+                        { on_axis ? dxz_gt11 - dxz_gt22 : (dz_gt11 - dz_gt22) / x, 0.0,
+                           on_axis ? dxz_gt13 : dz_gt13 / x },
+                        { 0.0, on_axis ? dxz_gt13 : dz_gt13 / x, 0.0 },
+                    },
+
+                },
+                {
+                    {
+                        { dxz_gt11,      0.0, dxz_gt13 },
+                        {      0.0, dxz_gt22,      0.0 },
+                        { dxz_gt13,      0.0, dxz_gt33 },
+                    },
+                    {
+                        { 0.0, on_axis ? dxz_gt11 - dxz_gt22 : (dz_gt11 - dz_gt22) / x, 0.0 },
+                        { on_axis ? dxz_gt11 - dxz_gt22 : (dz_gt11 - dz_gt22) / x, 0.0,
+                           on_axis ? dxz_gt13 : dz_gt13 / x },
+                        { 0.0, on_axis ? dxz_gt13 : dz_gt13 / x, 0.0 },
+                    },
+                    {
+                        { dzz_gt11,      0.0, dzz_gt13 },
+                        {      0.0, dzz_gt22,      0.0 },
+                        { dzz_gt13,      0.0, dzz_gt33 },
+                    },
+
+                },
+            };
+
+            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 At[3][3] = {
+                { Atxx, Atxy, Atxz },
+                { 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 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 dAt[3][3][3] = {
+                {
+                    { dx_At11,     0.0, dx_At13 },
+                    {     0.0, dx_At22,     0.0 },
+                    { dx_At13,     0.0, dx_At33 },
+                },
+                {
+                    {     0.0, on_axis ? dx_At11 - dx_At22 : (Atxx - Atyy) / x, 0.0 },
+                    { on_axis ? dx_At11 - dx_At22 : (Atxx - Atyy) / x, 0.0, on_axis ? dx_At13 : Atxz / x },
+                    { 0.0, on_axis ? dx_At13 : Atxz / x, 0.0 },
+                },
+                {
+                    { dz_At11,     0.0, dz_At13 },
+                    {     0.0, dz_At22,     0.0 },
+                    { dz_At13,     0.0, dz_At33 },
+                },
+            };
+
+            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 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 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 },
+            };
+
+            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 dtrK[3] = { dx_trk, 0.0, dz_trk };
+
+            const double kdot_xx = a_kdotxx[idx_src];
+            const double kdot_yy = a_kdotyy[idx_src];
+            const double kdot_zz = a_kdotzz[idx_src];
+            const double kdot_xy = a_kdotxy[idx_src];
+            const double kdot_xz = a_kdotxz[idx_src];
+            const double kdot_yz = a_kdotyz[idx_src];
+
+            const double kdot[3][3] = {{ kdot_xx, kdot_xy, kdot_xz },
+                                       { kdot_xy, kdot_yy, kdot_yz },
+                                       { kdot_xz, kdot_yz, kdot_zz }};
+
+            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  = 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 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 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 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_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 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_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 d2beta[3][3][3] = {
+                {
+                    { dxx_betax, 0.0, dxx_betaz },
+                    { 0.0, on_axis ? 0.5 * dxx_betax : dx_betax / x - betax / SQR(x), 0.0 },
+                    { dxz_betax, 0.0, dxz_betaz },
+                },
+                {
+                    { 0.0, on_axis ? 0.5 * dxx_betax : dx_betax / x - betax / SQR(x), 0.0 },
+                    { on_axis ? 0.5 * dxx_betax : dx_betax / x - betax / SQR(x), 0.0, on_axis ? dxx_betaz : dx_betaz / x },
+                    { 0.0, on_axis ? dxz_betax : dz_betax / x, 0.0 },
+                },
+                {
+                    { dxz_betax, 0.0, dxz_betaz },
+                    { 0.0, on_axis ? dxz_betax : dz_betax / x, 0.0 },
+                    { dzz_betax, 0.0, dzz_betaz },
+                },
+            };
+
+            const double det = gtxx * gtyy * gtzz + 2 * gtxy * gtyz * gtxz - gtzz * SQR(gtxy) - SQR(gtxz) * gtyy - gtxx * SQR(gtyz);
+
+            // \tilde{γ}^{ij}
+            gtu[0][0] =  (gtyy * gtzz - SQR(gtyz)) / det;
+            gtu[1][1] =  (gtxx * gtzz - SQR(gtxz)) / det;
+            gtu[2][2] =  (gtxx * gtyy - SQR(gtxy)) / det;
+            gtu[0][1] = -(gtxy * gtzz - gtyz * gtxz) / det;
+            gtu[0][2] =  (gtxy * gtyz - gtyy * gtxz) / det;
+            gtu[1][2] = -(gtxx * gtyz - gtxy * gtxz) / det;
+            gtu[1][0] = gtu[0][1];
+            gtu[2][0] = gtu[0][2];
+            gtu[2][1] = gtu[1][2];
+
+            //　γ_{jk}/^{jk}
+            for (int j = 0; j < 3; j++)
+                for (int k = 0; k < 3; k++) {
+                    gu[j][k] = SQR(phi) * gtu[j][k];
+                    g[j][k]  = gt[j][k] / SQR(phi);
+                }
+
+            // β_j
+            for (int j = 0; j < 3; j++) {
+                double val = 0.0;
+                for (int k = 0; k < 3; k++)
+                    val += g[j][k] * beta[k];
+                betal[j] = val;
+            }
+
+            // ∂_j γ_{kl}
+            // ∂_j K_{kl}
+            for (int j = 0; j < 3; j++)
+                for (int k = 0; k < 3; k++)
+                    for (int l = 0; l < 3; l++) {
+                        dg[j][k][l] = -2.0 * dphi[j] * gt[k][l] / (phi * SQR(phi)) + dgt[j][k][l] / SQR(phi);
+                        dK[j][k][l] = -2.0 * dphi[j] * At[k][l] / (phi * SQR(phi)) + dAt[j][k][l] / SQR(phi) +
+                                      (1.0 / 3.0) * (dtrK[j] * g[k][l] + trK * dg[j][k][l]);
+                    }
+
+            // ∂_{jk} g_{lm}
+            for (int j = 0; j < 3; j++)
+                for (int k = 0; k < 3; k++)
+                    for (int l = 0; l < 3; l++)
+                        for (int m = 0; m < 3; m++) {
+                            d2g[j][k][l][m] = 6.0 *  gt      [l][m] * dphi[j] * dphi[k] / SQR(SQR(phi))    -
+                                              2.0 *  gt      [l][m] * d2phi[j][k]       / (phi * SQR(phi)) -
+                                              2.0 * dgt   [j][l][m] * dphi[k]           / (phi * SQR(phi)) -
+                                              2.0 * dgt   [k][l][m] * dphi[j]           / (phi * SQR(phi)) +
+                                                   d2gt[j][k][l][m]                     / SQR(phi);
+                        }
+
+            // ∂_j γ^{kl}
+            for (int j = 0; j < 3; j++)
+                for (int k = 0; k < 3; k++)
+                    for (int l = 0; l < 3; l++) {
+                        double val = 0.0;
+                        for (int m = 0; m < 3; m++)
+                            for (int n = 0; n < 3; n++)
+                                val += -gu[k][m] * gu[l][n] * dg[j][m][n];
+                        dgu[j][k][l] = val;
+                    }
+
+            // Γ^j_{kl}
+            for (int j = 0; j < 3; j++)
+                for (int k = 0; k < 3; k++)
+                    for (int l = 0; l < 3; l++) {
+                        double val = 0.0;
+                        for (int m = 0; m < 3; m++)
+                            val += 0.5 * gu[j][m] * (dg[k][l][m] + dg[l][k][m] - dg[m][k][l]);
+                        G[j][k][l] = val;
+                    }
+
+            // ∂_j Γ^k_{lm}
+            for (int j = 0; j < 3; j++)
+                for (int k = 0; k < 3; k++)
+                    for (int l = 0; l < 3; l++)
+                        for (int m = 0; m < 3; m++) {
+                            double val = 0.0;
+                            for (int n = 0; n < 3; n++) {
+                                val += dgu[j][k][n] * (dg    [l][m][n] +  dg   [m][l][n] -  dg   [n][l][m]) +
+                                        gu   [k][n] * (d2g[j][l][m][n] + d2g[j][m][l][n] - d2g[j][n][l][m]);
+                            }
+                            dG[j][k][l][m] = 0.5 * val;
+                        }
+
+            // Γ^j = γ^{kl} Γ_{kl}^j
+            for (int j = 0; j < 3; j++) {
+                double val = 0.0;
+                for (int k = 0; k < 3; k++)
+                    for (int l = 0; l < 3; l++)
+                        val += gu[k][l] * G[j][k][l];
+                X[j] = val;
+            }
+
+            // ∂_j β_k
+            for (int j = 0; j < 3; j++)
+                for (int k = 0; k < 3; k++) {
+                    double val = 0.0;
+                    for (int l = 0; l < 3; l++)
+                        val += g[k][l] * dbeta[j][l] + dg[j][k][l] * beta[l];
+                    dbetal[j][k] = val;
+                }
+
+            // ∂_{jk} β_l
+            for (int j = 0; j < 3; j++)
+                for (int k = 0; k < 3; k++)
+                    for (int l = 0; l < 3; l++) {
+                        double val = 0.0;
+                        for (int m = 0; m < 3; m++)
+                            val += d2g[j][k][l][m] * beta[m] + dg[k][l][m] * dbeta[j][m] + dg[j][l][m] * dbeta[k][m] +
+                                   g[l][m] * d2beta[j][k][m];
+                        d2betal[j][k][l] = val;
+                    }
+
+            // K_{ij}
+            for (int j = 0; j < 3; j++)
+                for (int k = 0; k < 3; k++)
+                    K[j][k] = At[j][k] / SQR(phi) + (1.0 / 3.0) * g[j][k] * trK;
+
+            for (int j = 0; j < 3; j++)
+                for (int k = 0; k < 3; k++) {
+                    double val = 0.0;
+                    for (int l = 0; l < 3; l++)
+                        val += gu[j][l] * K[l][k];
+                    Km[j][k] = val;
+                }
+
+            // K^{ij}
+            for (int j = 0; j < 3; j++)
+                for (int k = 0; k < 3; k++) {
+                    double val = 0.0;
+                    for (int l = 0; l < 3; l++)
+                        val += gu[j][l] * Km[k][l];
+                    Ku[j][k] = val;
+                }
+
+            // ∂_j \dot{γ_kl}
+            for (int j = 0; j < 3; j++)
+                for (int k = 0; k < 3; k++)
+                    for (int l = 0; l < 3; l++) {
+                        dgdot[j][k][l] = -2.0 * (dalpha[j] * K[k][l] + alpha * dK[j][k][l]) +
+                                         d2betal[j][k][l] + d2betal[j][l][k];
+                        for (int m = 0; m < 3; m++)
+                            dgdot[j][k][l] += -2.0 * (dG[j][m][k][l] * betal[m] + G[m][k][l] * dbetal[j][m]);
+                    }
+
+            // \dot{γ^{jk}}
+            for (int j = 0; j < 3; j++)
+                for (int k = 0; k < 3; k++) {
+                    gudot[j][k] = 2.0 * alpha * Ku[j][k];
+
+                    for (int l = 0; l < 3; l++) {
+                        gudot[j][k] += -gu[j][l] * dbeta[l][k] - gu[k][l] * dbeta[l][j];
+
+                        for (int m = 0; m < 3; m++)
+                            gudot[j][k] += -gu[j][l] * G[k][l][m] * beta[m] - gu[k][l] * G[j][l][m] * beta[m];
+                    }
+                }
+
+
+            // \dot{Γ}^j_{kl}
+            for (int j = 0; j < 3; j++)
+                for (int k = 0; k < 3; k++)
+                    for (int l = 0; l < 3; l++) {
+                        double val = 0.0;
+                        for (int m = 0; m < 3; m++)
+                            val += gudot[j][m] * (dg   [k][l][m] + dg   [l][k][m] - dg   [m][k][l]) +
+                                   gu   [j][m] * (dgdot[k][l][m] + dgdot[l][k][m] - dgdot[m][k][l]);
+
+                        Gdot[j][k][l] = 0.5 * val;
+                    }
+
+            Gammadot_term = 0.0;
+            for (int j = 0; j < 3; j++) {
+                double val = 0.0;
+                for (int k = 0; k < 3; k++)
+                    for (int l = 0; l < 3; l++) {
+                        val += gu[k][l] * Gdot[j][k][l];
+                    }
+
+                Gammadot_term += val * dalpha[j];
+            }
+
+            Kij_Dij_alpha = 0.0;
+            for (int j = 0; j < 3; j++)
+                for (int k = 0; k < 3; k++) {
+                    double val = 0.0;
+                    for (int l = 0; l < 3; l++)
+                        val += G[l][j][k] * dalpha[l];
+
+                    Kij_Dij_alpha += Ku[j][k] * (d2alpha[j][k] - val);
+                }
+
+            k_kdot = 0.0;
+            for (int j = 0; j < 3; j++)
+                for (int k = 0; k < 3; k++)
+                    k_kdot += kdot[j][k] * Ku[j][k];
+
+            k3 = 0.0;
+            for (int j = 0; j < 3; j++)
+                for (int k = 0; k < 3; k++) {
+                    double val = 0.0;
+                    for (int l = 0; l < 3; l++)
+                        val += Km[k][l] * Km[l][j];
+                    k3 += val * Km[j][k];
+                }
+
+            // K_{ij} K^{ij}
+            k2 = 0.0;
+            for (int j = 0; j < 3; j++)
+                for (int k = 0; k < 3; k++)
+                    k2 += Km[j][k] * Km[k][j];
+
+            // Ric_{jk}
+            for (int j = 0; j < 3; j++)
+                for (int k = 0; k < 3; k++) {
+                    double val = 0.0;
+                    for (int m = 0; m < 3; m++)
+                        val += dG[m][m][j][k] - dG[k][m][j][m];
+                    for (int m = 0; m < 3; m++)
+                        for (int l = 0; l < 3; l++)
+                            val += G[l][l][m] * G[m][j][k] - G[l][k][m] * G[m][j][l];
+                    Ric[j][k] = val;
+                }
+
+            // Ric^j_k
+            for (int j = 0; j < 3; j++)
+                for (int k = 0; k < 3; k++) {
+                    double val = 0.0;
+                    for (int l = 0; l < 3; l++)
+                        val += gu[j][l] * Ric[l][k];
+                    Ricm[j][k] = val;
+                }
+
+            beta_term = 0.0;
+            for (int j = 0; j < 3; j++)
+                for (int k = 0; k < 3; k++)
+                    for (int l = 0; l < 3; l++) {
+                        double D2beta = 0.0;
+
+                        D2beta += d2beta[j][k][l];
+
+                        for (int m = 0; m < 3; m++) {
+                            D2beta += dG[j][l][k][m] * beta[m] + G[l][k][m] * dbeta[j][m] + G[l][j][m] * dbeta[k][m] - G[m][j][k] * dbeta[m][l];
+
+                            for (int n = 0; n < 3; n++)
+                                D2beta += G[l][j][m] * G[m][k][n] * beta[n] - G[m][j][k] * G[l][m][n] * beta[n];
+                        }
+
+                        beta_term += -gu[j][k] * D2beta * dalpha[l];
+                    }
+            for (int j = 0; j < 3; j++)
+                for (int k = 0; k < 3; k++)
+                    beta_term += -beta[k] * Ricm[j][k] * dalpha[j];
+
+            solver->diff_coeffs[MG2D_DIFF_COEFF_20]->data[idx_dc] = gu[0][0] + (on_axis ? gu[1][1] : 0.0);
+            solver->diff_coeffs[MG2D_DIFF_COEFF_02]->data[idx_dc] = gu[2][2];
+            solver->diff_coeffs[MG2D_DIFF_COEFF_11]->data[idx_dc] = 2.0 * gu[0][2];
+            solver->diff_coeffs[MG2D_DIFF_COEFF_10]->data[idx_dc] = -X[0] + (on_axis ? 0.0 :  gu[1][1] / x);
+            solver->diff_coeffs[MG2D_DIFF_COEFF_01]->data[idx_dc] = -X[2];
+            solver->diff_coeffs[MG2D_DIFF_COEFF_00]->data[idx_dc] = -k2;
+
+            solver->rhs[idx_rhs] = -2.0 * alpha * Kij_Dij_alpha + Gammadot_term +
+                                    2.0 * alpha * (k_kdot + 2.0 * alpha * k3) + beta_term;
+
+            if (on_axis) {
+                solver->diff_coeffs[MG2D_DIFF_COEFF_20]->boundaries[MG2D_BOUNDARY_0L].val[idx_z] = gu[1][1];
+                solver->diff_coeffs[MG2D_DIFF_COEFF_10]->boundaries[MG2D_BOUNDARY_0L].val[idx_z] = gu[1][1];
+            }
+        }
+}
+
+static void solution_to_grid(CoordPatch *cp, MG2DContext *solver, double *dst)
+{
+#pragma omp parallel for
+    for (int j = 0; j < solver->local_size[1]; j++)
+        for (int i = 0; i < solver->local_size[0]; i++) {
+            const ptrdiff_t idx_dst = CPINDEX(cp, i + cp->offset_left[0], cp->y_idx, j + cp->offset_left[1]);
+            const int idx_src = j * solver->u_stride + i;
+            dst[idx_dst] = solver->u[idx_src];
+        }
+
+    /* fill in the axial symmetry ghostpoints by mirroring */
+    if (!cp->bnd_intercomp[0][0]) {
+#pragma omp parallel for
+        for (int idx_z = cp->offset_left[1]; idx_z < cp->grid_size[2]; idx_z++)
+            for (int idx_x = 0; idx_x < cp->offset_left[0]; idx_x++) {
+                const ptrdiff_t idx_dst = CPINDEX(cp, idx_x, cp->y_idx, idx_z);
+                const ptrdiff_t idx_src = CPINDEX(cp, 2 * cp->offset_left[0] - idx_x, cp->y_idx, idx_z);
+                dst[idx_dst] = dst[idx_src];
+            }
+    }
+    if (!cp->bnd_intercomp[1][0]) {
+#pragma omp parallel for
+        for (int idx_z = 0; idx_z < cp->offset_left[1]; idx_z++)
+            for (int idx_x = 0; idx_x < cp->grid_size[0]; idx_x++) {
+                const ptrdiff_t idx_dst = CPINDEX(cp, idx_x, cp->y_idx, idx_z);
+                const ptrdiff_t idx_src = CPINDEX(cp, idx_x, cp->y_idx, 2 * cp->offset_left[1] - idx_z);
+                dst[idx_dst] = dst[idx_src];
+            }
+    }
+}
+
+void qms_mg_solve(CCTK_ARGUMENTS)
+{
+    QMSMGContext *ms = qms_context;
+    CoordPatch *cp;
+
+    DECLARE_CCTK_ARGUMENTS;
+    DECLARE_CCTK_PARAMETERS;
+
+    const int reflevel = ctz(cctkGH->cctk_levfac[0]);
+    const double time  = ms->gh->cctk_time;
+
+    int64_t grid_size;
+    int ret;
+
+
+    if (cctkGH->cctk_time >= switchoff_time + 1.0)
+        return;
+
+    if (reflevel < ms->solve_level)
+        return;
+    if (reflevel > ms->solve_level &&
+        fabs(time - W_val1_time[ms->solve_level]) > 1e-13)
+        return;
+
+    cp = get_coord_patch(ms, reflevel);
+
+    grid_size = cctkGH->cctk_lsh[0] * cctkGH->cctk_lsh[1] * cctkGH->cctk_lsh[2];
+
+    fill_eq_coeffs(ms, cp, cp->solver);
+
+    if (reflevel > ms->solve_level) {
+        /* fill the solver boundary conditions */
+        if (cp->solver->local_start[1] + cp->solver->local_size[1] == cp->solver->domain_size) {
+#pragma omp parallel for
+            for (int j = 0; j < cp->solver->fd_stencil; j++) {
+                MG2DBoundary *bnd = cp->solver->boundaries[MG2D_BOUNDARY_1U];
+                double *dst = bnd->val + j * bnd->val_stride;
+                for (ptrdiff_t i = -j; i < (ptrdiff_t)cp->solver->local_size[0] + j; i++) {
+                    const ptrdiff_t idx = CCTK_GFINDEX3D(cctkGH, ABS(i) + cp->offset_left[0], cp->y_idx, cp->offset_left[1] + cp->solver->local_size[1] - 1 + j);
+                    dst[i] = W_val1[idx];
+                }
+            }
+        }
+        if (cp->solver->local_start[0] + cp->solver->local_size[0] == cp->solver->domain_size) {
+#pragma omp parallel for
+            for (int j = 0; j < cp->solver->fd_stencil; j++) {
+                MG2DBoundary *bnd = cp->solver->boundaries[MG2D_BOUNDARY_0U];
+                double *dst = bnd->val + j * bnd->val_stride;
+                for (ptrdiff_t i = -j; i < (ptrdiff_t)cp->solver->local_size[1] + j; i++) {
+                    const ptrdiff_t idx = CCTK_GFINDEX3D(cctkGH, cp->offset_left[1] + cp->solver->local_size[0] - 1 + j, cp->y_idx, cp->offset_left[1] + ABS(i));
+                    dst[i] = W_val1[idx];
+                }
+            }
+        }
+    }
+
+    fprintf(stderr, "%d qms solve: time %g\n", reflevel, time);
+    ret = mg2d_solve(cp->solver);
+    if (ret < 0)
+        CCTK_WARN(0, "Error solving the quasi-maximal slicing equation");
+
+    /* add the solution to the list of past solutions */
+    {
+        memcpy(W_val0, W_val1, sizeof(*W_val0) * grid_size);
+        W_val0_time[reflevel] = W_val1_time[reflevel];
+
+        solution_to_grid(cp, cp->solver, W_val1);
+        W_val1_time[reflevel] = ms->gh->cctk_time;
+    }
+
+    /* linearly extrapolate the past two solution to predict the next one */
+    {
+        double *u0 = W_val0;
+        double *u1 = W_val1;
+        double time0 = W_val0_time[reflevel];
+        double time1 = W_val1_time[reflevel];
+        double time  = time1 + ms->gh->cctk_delta_time / (1.0 * (1 << ms->solve_level));
+
+        double fact0, fact1;
+
+        if (time0 == DBL_MAX && time1 == DBL_MAX) {
+            fact0 = 0.0;
+            fact1 = 0.0;
+        } else if (time0 == DBL_MAX) {
+            fact0 = 0.0;
+            fact1 = 0.0;
+        } else {
+            fact0 = (time - time1) / (time0 - time1);
+            fact1 = (time - time0) / (time1 - time0);
+        }
+
+        memcpy(W_pred0, W_pred1, sizeof(*W_pred0) * grid_size);
+        W_pred0_time[reflevel] = W_pred1_time[reflevel];
+
+        for (int i = 0; i < grid_size; i++)
+            W_pred1[i] = (u1[i] * fact1 + u0[i] * fact0);
+
+        W_pred1_time[reflevel] = time;
+    }
+}
+
+void qms_mg_sync(CCTK_ARGUMENTS)
+{
+}
+
+void qms_mg_eval(CCTK_ARGUMENTS)
+{
+    QMSMGContext *ms;
+
+    DECLARE_CCTK_ARGUMENTS;
+    DECLARE_CCTK_PARAMETERS;
+
+    const int      reflevel = ctz(cctkGH->cctk_levfac[0]);
+    const double       time = cctkGH->cctk_time;
+    const int64_t grid_size = cctkGH->cctk_lsh[0] * cctkGH->cctk_lsh[1] * cctkGH->cctk_lsh[2];
+
+    int i, ret;
+
+    ms = qms_context;
+
+    if (reflevel < ms->solve_level || time >= switchoff_time + 1.0)
+        return;
+
+    {
+        double *u0 = W_pred0;
+        double *u1 = W_pred1;
+        double time0 = W_pred0_time[reflevel];
+        double time1 = W_pred1_time[reflevel];
+
+        double fact, fact0, fact1;
+
+        if (time0 == DBL_MAX && time1 == DBL_MAX) {
+            fact0 = 0.0;
+            fact1 = 0.0;
+        } else if (time0 == DBL_MAX) {
+            fact0 = 0.0;
+            fact1 = 1.0;
+        } else {
+            fact0 = (time - time1) / (time0 - time1);
+            fact1 = (time - time0) / (time1 - time0);
+        }
+
+        if (time < switchoff_time)
+            fact = 1.0;
+        else
+            fact = exp(-36.0 * pow((time - switchoff_time), 4.0));
+
+#pragma omp parallel for
+        for (int i = 0; i < grid_size; i++)
+            W[i] = (u1[i] * fact1 + u0[i] * fact0) * fact;
+
+    }
+}
+
+static int context_init(cGH *gh, int solve_level, int fd_stencil, int maxiter, int exact_size, int nb_cycles,
+                        int nb_relax_pre, int nb_relax_post, double tol_residual_base,
+                        double cfl_factor, const char *loglevel_str,
+                        QMSMGContext **ctx)
+{
+    QMSMGContext *qms;
+    int loglevel = -1;
+    int ret;
+
+    qms = calloc(1, sizeof(*qms));
+    if (!qms)
+        return -ENOMEM;
+
+    qms->gh = gh;
+
+    qms->solve_level = solve_level;
+    qms->fd_stencil    = fd_stencil;
+    qms->maxiter       = maxiter;
+    qms->max_exact_size = exact_size;
+    qms->nb_cycles     = nb_cycles;
+    qms->nb_relax_pre  = nb_relax_pre;
+    qms->nb_relax_post = nb_relax_post;
+    qms->tol_residual_base = tol_residual_base;
+    qms->cfl_factor    = cfl_factor;
+
+   for (int i = 0; i < ARRAY_ELEMS(log_levels); i++) {
+        if (!strcmp(loglevel_str, log_levels[i].str)) {
+            loglevel = log_levels[i].level;
+            break;
+        }
+    }
+    if (loglevel < 0) {
+        fprintf(stderr, "Invalid loglevel: %s\n", loglevel_str);
+        return -EINVAL;
+    }
+    qms->log_level = loglevel;
+
+    *ctx = qms;
+
+    return 0;
+}
+
+void qms_mg_init(CCTK_ARGUMENTS)
+{
+    DECLARE_CCTK_ARGUMENTS;
+    DECLARE_CCTK_PARAMETERS;
+
+    int use_tapered_grids = *(int*)CCTK_ParameterGet("use_tapered_grids", "Carpet", &use_tapered_grids);
+    int ret;
+
+    fprintf(stderr, "%d init\n", ctz(cctk_levfac[0]));
+
+    if (!use_tapered_grids)
+        CCTK_WARN(0, "MaximalSlicingAxiMG only works with use_tapered_grids=1");
+
+    if (!qms_context) {
+        ret = context_init(cctkGH, solve_level, fd_stencil, maxiter, exact_size, nb_cycles,
+                           nb_relax_pre, nb_relax_post, tol_residual_base,
+                           cfl_factor, loglevel, &qms_context);
+        if (ret < 0)
+            CCTK_WARN(0, "Error initializing the solver context");
+    }
+}
+
+void qms_mg_inithist(CCTK_ARGUMENTS)
+{
+    DECLARE_CCTK_ARGUMENTS;
+    DECLARE_CCTK_PARAMETERS;
+
+    size_t grid_size = cctk_lsh[0] * cctk_lsh[1] * cctk_lsh[2];
+
+    fprintf(stderr, "%d inithist\n", ctz(cctk_levfac[0]));
+
+    for (int i = 0; i < 32; i++) {
+        W_pred0_time[i] = DBL_MAX;
+        W_pred1_time[i] = DBL_MAX;
+        W_val0_time[i]  = DBL_MAX;
+        W_val1_time[i]  = DBL_MAX;
+    }
+
+    memset(Kdot11, 0, sizeof(double) * grid_size);
+    memset(Kdot22, 0, sizeof(double) * grid_size);
+    memset(Kdot33, 0, sizeof(double) * grid_size);
+    memset(Kdot12, 0, sizeof(double) * grid_size);
+    memset(Kdot13, 0, sizeof(double) * grid_size);
+    memset(Kdot23, 0, sizeof(double) * grid_size);
+
+    memset(Xtdot1, 0, sizeof(double) * grid_size);
+    memset(Xtdot2, 0, sizeof(double) * grid_size);
+    memset(Xtdot3, 0, sizeof(double) * grid_size);
+
+    memset(W_pred0, 0, sizeof(double) * grid_size);
+    memset(W_pred1, 0, sizeof(double) * grid_size);
+    memset(W_val0, 0, sizeof(double) * grid_size);
+    memset(W_val1, 0, sizeof(double) * grid_size);
+}