Rewrite using a new method.

Based on 'Adaptive mesh refinement for coupled elliptic-hyperbolic
problems' by Pretorius and Choptuik, 2006.

3 files changed, 326 insertions, 143 deletions

diff --git a/interface.ccl b/interface.ccl
index 43e512c..124823e 100644
--- a/interface.ccl
+++ b/interface.ccl
@@ -13,3 +13,11 @@ USES FUNCTION MoLNumIntegratorSubsteps
 REQUIRES FUNCTION MoLRegisterConstrained
 REQUIRES FUNCTION MoLRegisterSaveAndRestore
 REQUIRES FUNCTION MoLRegisterSaveAndRestoreGroup
+
+public:
+REAL lapse_mg TYPE=GF TIMELEVELS=2
+REAL lapse_mg_eval TYPE=GF tags='Prolongation="None"'
+REAL lapse_prev0 TYPE=GF tags='Prolongation="None"'
+REAL lapse_prev1 TYPE=GF tags='Prolongation="None"'
+REAL lapse_prev0_time TYPE=ARRAY DIM=1 SIZE=32
+REAL lapse_prev1_time TYPE=ARRAY DIM=1 SIZE=32
diff --git a/schedule.ccl b/schedule.ccl
index 246aa63..0b992b9 100644
--- a/schedule.ccl
+++ b/schedule.ccl
@@ -1,10 +1,26 @@
 # Schedule definitions for thorn MaximalSlicingAxiMG
 #
 if (CCTK_Equals(lapse_evolution_method, "maximal_axi_mg")) {
+    SCHEDULE msa_mg_inithist in CCTK_INITIAL {
+        LANG: C
+    } ""
 
-    SCHEDULE maximal_slicing_axi_mg IN ML_BSSN_evolCalcGroup BEFORE ML_BSSN_RHS {
+    SCHEDULE msa_mg_init IN CCTK_BASEGRID {
         LANG: C
-    } "Maximal slicing in axisymmetry"
+    } ""
+
+    SCHEDULE msa_mg_eval IN MoL_CalcRHS BEFORE ML_BSSN_evolCalcGroup {
+        LANG: C
+    } ""
+
+    SCHEDULE msa_mg_sync IN CCTK_POSTSTEP BEFORE msa_mg_solve {
+        SYNC: lapse_mg
+        LANG: C
+    } ""
+
+    SCHEDULE msa_mg_solve IN CCTK_POSTSTEP {
+        LANG: C
+    } ""
 
     #SCHEDULE maximal_slicing_axi IN MoL_PostStepModify {
     #    LANG: C
@@ -21,4 +37,11 @@ if (CCTK_Equals(lapse_evolution_method, "maximal_axi_mg")) {
     SCHEDULE maximal_slicing_axi_mg_register_mol IN MoL_Register {
         LANG: C
     } ""
+
+    STORAGE: lapse_mg[2]
+    STORAGE: lapse_mg_eval
+    STORAGE: lapse_prev0
+    STORAGE: lapse_prev1
+    STORAGE: lapse_prev0_time
+    STORAGE: lapse_prev1_time
 }
diff --git a/src/maximal_slicing_axi_mg.c b/src/maximal_slicing_axi_mg.c
index 974ef86..0838d6a 100644
--- a/src/maximal_slicing_axi_mg.c
+++ b/src/maximal_slicing_axi_mg.c
@@ -21,12 +21,14 @@
 #define SQR(x) ((x) * (x))
 #define ARRAY_ELEMS(x) (sizeof(x) / sizeof(*x))
 
+#define CPINDEX(cp, i, j, k) ((k * cp->grid_size[1] + j) * cp->grid_size[0] + i)
+
+typedef struct MaximalSlicingContext MaximalSlicingContext;
+
 /* precomputed values for a given refined grid */
 typedef struct CoordPatch {
-    double origin[3];
-    int    delta[3];
-    int    size[3];
     int    level;
+    ptrdiff_t grid_size[3];
 
     MG2DContext *solver;
 
@@ -37,21 +39,6 @@ typedef struct CoordPatch {
     ptrdiff_t offset_left[2];
     ptrdiff_t offset_right[2];
 
-    /* boundary interpolator parameters */
-    double *interp_coords[3];
-
-    double *boundary_interp_vals;
-    size_t nb_boundary_interp_vals;
-
-    int interp_values_codes[1];
-    int interp_var_indices[1];
-    int interp_operation_codes[1];
-    int interp_operation_indices[1];
-
-    int coord_system;
-    int interp_operator;
-    int interp_params;
-
     /**
      * pointers into boundary_interp_vals
      * the first row are the boundary points for the solver
@@ -60,9 +47,7 @@ typedef struct CoordPatch {
      * [1] is the upper-x patch
      *     rows run from z=0 to z=<upper physical boundary>
      */
-    double   *boundary_val[2];
     ptrdiff_t boundary_val_stride[2];
-
 } CoordPatch;
 
 typedef struct MSMGContext {
@@ -70,8 +55,17 @@ typedef struct MSMGContext {
 
     CoordPatch *patches;
     int nb_patches;
+
+    MaximalSlicingContext *ps_solver;
 } MSMGContext;
 
+int msa_context_init(cGH *gh, MaximalSlicingContext **ctx, int basis_order0, int basis_order1, double scale_factor, double filter_power);
+void msa_context_free(MaximalSlicingContext **ctx);
+void msa_lapse_solve(MaximalSlicingContext *ctx);
+void msa_lapse_eval(MaximalSlicingContext *ctx, double *dst);
+
+static MSMGContext *ms;
+
 static int ctz(int a)
 {
     int ret = 0;
@@ -92,10 +86,11 @@ static void coord_patch_free(CoordPatch *cp)
     mg2d_solver_free(&cp->solver);
 }
 
-static CoordPatch *get_coord_patch(MSMGContext *ms)
+static CoordPatch *get_coord_patch(MSMGContext *ms, int level)
 {
     cGH *gh = ms->gh;
 
+    const size_t grid_size = gh->cctk_lsh[2] * gh->cctk_lsh[1] * gh->cctk_lsh[0];
     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");
@@ -108,15 +103,7 @@ static CoordPatch *get_coord_patch(MSMGContext *ms)
     for (int i = 0; i < ms->nb_patches; i++) {
         cp = &ms->patches[i];
 
-        if (cp->origin[0] == ms->gh->cctk_origin_space[0] &&
-            cp->origin[1] == ms->gh->cctk_origin_space[1] &&
-            cp->origin[2] == ms->gh->cctk_origin_space[2] &&
-            cp->size[0]   == ms->gh->cctk_lsh[0]          &&
-            cp->size[1]   == ms->gh->cctk_lsh[1]          &&
-            cp->size[2]   == ms->gh->cctk_lsh[2]          &&
-            cp->delta[0]  == ms->gh->cctk_levfac[0]  &&
-            cp->delta[1]  == ms->gh->cctk_levfac[1]  &&
-            cp->delta[2]  == ms->gh->cctk_levfac[2])
+        if (cp->level == level)
             return cp;
     }
 
@@ -126,17 +113,17 @@ static CoordPatch *get_coord_patch(MSMGContext *ms)
 
     memset(cp, 0, sizeof(*cp));
 
-    memcpy(cp->origin, ms->gh->cctk_origin_space, sizeof(cp->origin));
-    memcpy(cp->size,   ms->gh->cctk_lsh,          sizeof(cp->size));
-    memcpy(cp->delta,  ms->gh->cctk_levfac,  sizeof(cp->delta));
     cp->level = ctz(ms->gh->cctk_levfac[0]);
+    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 < cp->size[1]; i++)
+    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 == cp->size[1])
+    if (i == gh->cctk_lsh[1])
         CCTK_WARN(0, "The grid does not include y==0");
 
     for (int i = 0; i < 2; i++) {
@@ -146,6 +133,11 @@ static CoordPatch *get_coord_patch(MSMGContext *ms)
         /* account for grid tapering on refined levels */
         if (cp->level > 0)
             offset_right *= integrator_substeps * 2;
+        ///* the outer boundary layer overlaps with the first ghost zone */
+        //if (cp->level > 0)
+        //    offset_right++;
+        if (cp->level > 0)
+            offset_right--;
 
         domain_size[i] = ms->gh->cctk_lsh[2 * i] - offset_left - offset_right;
         /* the outer boundary layer overlaps with the first ghost zone */
@@ -182,98 +174,9 @@ static CoordPatch *get_coord_patch(MSMGContext *ms)
     for (int i = 0; i < 4; i++)
         memset(cp->solver->boundaries[i]->val, 0, sizeof(double) * cp->solver->boundaries[i]->val_len);
 
-    /* initialize the interpolation state */
     if (cp->level > 0) {
         cp->boundary_val_stride[0] = gh->cctk_lsh[0] - cp->offset_left[0];
         cp->boundary_val_stride[1] = cp->solver->domain_size;
-
-        cp->nb_boundary_interp_vals = cp->boundary_val_stride[0] * cp->offset_right[1] +
-                                      cp->boundary_val_stride[1] * cp->offset_right[0];
-
-        ret = posix_memalign((void**)&cp->boundary_interp_vals, 32,
-                             cp->nb_boundary_interp_vals * sizeof(*cp->boundary_interp_vals));
-        if (ret != 0)
-            CCTK_WARN(0, "Error allocating arrays for boundary values");
-        cp->boundary_val[0] = cp->boundary_interp_vals;
-        cp->boundary_val[1] = cp->boundary_interp_vals + cp->boundary_val_stride[0] * cp->offset_right[1];
-
-        for (int i = 0; i < 3; i++) {
-            ret = posix_memalign((void**)&cp->interp_coords[i], 32,
-                                 cp->nb_boundary_interp_vals * sizeof(*cp->interp_coords[i]));
-            if (ret != 0)
-                CCTK_WARN(0, "Error allocating interp coords");
-        }
-
-        for (ptrdiff_t idx_row = 0; idx_row < cp->offset_right[1]; idx_row++)
-            for (ptrdiff_t i = 0; i < cp->boundary_val_stride[0]; i++) {
-                ptrdiff_t idx_dst = idx_row * cp->boundary_val_stride[0] + i;
-                ptrdiff_t idx_src_x = CCTK_GFINDEX3D(gh, cp->offset_left[0] + i, 0, 0);
-                ptrdiff_t idx_src_z = CCTK_GFINDEX3D(gh, 0, 0, gh->cctk_lsh[2] - cp->offset_right[1] + idx_row);
-
-                cp->interp_coords[0][idx_dst] = a_x[idx_src_x];
-                cp->interp_coords[1][idx_dst] = 0.0;
-                cp->interp_coords[2][idx_dst] = a_z[idx_src_z];
-            }
-
-        for (ptrdiff_t idx_row = 0; idx_row < cp->offset_right[0]; idx_row++)
-            for (ptrdiff_t i = 0; i < cp->boundary_val_stride[1]; i++) {
-                ptrdiff_t idx_dst = (cp->boundary_val[1] - cp->boundary_val[0]) +
-                                    idx_row * cp->boundary_val_stride[1] + i;
-                ptrdiff_t idx_src_x = CCTK_GFINDEX3D(gh, gh->cctk_lsh[0] - cp->offset_right[0] + idx_row,
-                                                     0, 0);
-                ptrdiff_t idx_src_z = CCTK_GFINDEX3D(gh, 0, 0, cp->offset_left[1] + i);
-
-                cp->interp_coords[0][idx_dst] = a_x[idx_src_x];
-                cp->interp_coords[1][idx_dst] = 0.0;
-                cp->interp_coords[2][idx_dst] = a_z[idx_src_z];
-            }
-
-#if 0
-        for (int i = 0; i < gh->cctk_lsh[0]; i++) {
-            cp->interp_coords[0][i] = a_x[CCTK_GFINDEX3D(gh, i, 0, 0)];
-            cp->interp_coords[1][i] = 0.0;
-            cp->interp_coords[2][i] = a_z[CCTK_GFINDEX3D(gh, 0, 0, gh->cctk_lsh[2] - gh->cctk_nghostzones[2] - 1)];
-        }
-        for (int i = 0; i < gh->cctk_lsh[2]; i++) {
-            int idx = gh->cctk_lsh[0] + i;
-            cp->interp_coords[0][idx] = a_x[CCTK_GFINDEX3D(gh, gh->cctk_lsh[0] - gh->cctk_nghostzones[0] - 1, 0, 0)];
-            cp->interp_coords[1][idx] = 0.0;
-            cp->interp_coords[2][idx] = a_z[CCTK_GFINDEX3D(gh, 0, 0, i)];
-        }
-#endif
-
-        cp->interp_values_codes[0] = CCTK_VARIABLE_REAL;
-        cp->interp_var_indices[0]  = CCTK_VarIndex("ML_BSSN::alpha");
-        if (cp->interp_var_indices[0] < 0)
-            CCTK_WARN(0, "Error getting the index of lapse variable");
-        cp->interp_operation_codes[0] = 0;
-        cp->interp_operation_indices[0] = 0;
-
-        cp->coord_system = CCTK_CoordSystemHandle("cart3d");
-        if (cp->coord_system < 0)
-            CCTK_WARN(0, "Error getting the coordinate system");
-
-        cp->interp_operator = CCTK_InterpHandle("Lagrange polynomial interpolation (tensor product)");
-        if (cp->interp_operator < 0)
-            CCTK_WARN(0, "Error getting the interpolation operator");
-
-        cp->interp_params = Util_TableCreateFromString("order=4 want_global_mode=1");
-        if (cp->interp_params < 0)
-            CCTK_WARN(0, "Error creating interpolation parameters table");
-
-        ret = Util_TableSetInt(cp->interp_params, cp->level, "max_reflevel");
-        if (ret < 0)
-            CCTK_WARN(0, "Error setting maximum reflevel for interpolation");
-
-        ret = Util_TableSetIntArray(cp->interp_params, ARRAY_ELEMS(cp->interp_operation_codes),
-                                    cp->interp_operation_codes, "operation_codes");
-        if (ret < 0)
-            CCTK_WARN(0, "Error setting operation codes");
-
-        ret = Util_TableSetIntArray(cp->interp_params, ARRAY_ELEMS(cp->interp_operation_indices),
-                                    cp->interp_operation_indices, "operand_indices");
-        if (ret < 0)
-            CCTK_WARN(0, "Error setting operand indices");
     }
 
     ms->nb_patches++;
@@ -363,8 +266,13 @@ 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);
 
@@ -413,39 +321,259 @@ static void fill_eq_coeffs(cGH *gh, CoordPatch *cp)
             cp->solver->diff_coeffs[MG2D_DIFF_COEFF_00][idx_dc] = -k2;
             cp->solver->rhs[idx_rhs]                            = k2;
         }
+}
 
-    if (cp->level > 0) {
-        ret = CCTK_InterpGridArrays(gh, 3, cp->interp_operator, cp->interp_params,
-                                    cp->coord_system, cp->nb_boundary_interp_vals, CCTK_VARIABLE_REAL,
-                                    (const void * const *)cp->interp_coords, ARRAY_ELEMS(cp->interp_var_indices),
-                                    cp->interp_var_indices, ARRAY_ELEMS(cp->interp_values_codes), cp->interp_values_codes,
-                                    (void * const *)&cp->boundary_interp_vals);
+static void solution_to_grid(CoordPatch *cp, double *dst)
+{
+    for (int j = 0; j < cp->solver->domain_size; j++)
+        for (int i = 0; i < cp->solver->domain_size; 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 * cp->solver->u_stride + i;
+            dst[idx_dst] = 1.0 + cp->solver->u[idx_src];
+        }
+
+    /* fill in the axial symmetry ghostpoints by mirroring */
+    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];
+        }
+    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];
+        }
+}
+
+static void solve_ps(MSMGContext *ctx, double *dst)
+{
+    if (!ctx->ps_solver)
+        msa_context_init(ctx->gh, &ctx->ps_solver, 60, 60, 8.0, 64.0);
+
+    msa_lapse_solve(ctx->ps_solver);
+    msa_lapse_eval(ctx->ps_solver, dst);
+}
+
+static int history_initialized(double *lapse_prev0_time, double *lapse_prev1_time, int nb_levels)
+{
+    for (int i = 0; i < nb_levels; i++) {
+        if (lapse_prev0_time[i] == DBL_MAX || lapse_prev1_time[i] == DBL_MAX)
+            return 0;
+    }
+    return 1;
+}
+
+void msa_mg_eval(CCTK_ARGUMENTS)
+{
+    DECLARE_CCTK_ARGUMENTS;
+    DECLARE_CCTK_PARAMETERS;
+
+    int nb_levels;
+    int nb_levels_type, ret;
+
+    const size_t grid_size = cctkGH->cctk_lsh[2] * cctkGH->cctk_lsh[1] * cctkGH->cctk_lsh[0];
+    const double t = cctkGH->cctk_time;
+    const int reflevel = ctz(cctkGH->cctk_levfac[0]);
+    double time_interp_step, fact0, fact1;
+
+    fprintf(stderr, "evaluating lapse at rl=%d, t=%g\n", reflevel, t);
+
+    nb_levels = *(int*)CCTK_ParameterGet("num_levels_1", "CarpetRegrid2", &nb_levels_type);
+
+    if (!history_initialized(lapse_prev0_time, lapse_prev1_time, nb_levels)) {
+        solve_ps(ms, lapse_mg_eval);
+        memcpy(alpha, lapse_mg_eval, sizeof(*alpha) * grid_size);
+        fprintf(stderr, "pseudospectral solve\n");
+        return;
+    }
+
+    if (lapse_prev1_time[reflevel] <= lapse_prev0_time[reflevel] ||
+        fabs(lapse_prev1_time[reflevel] - lapse_prev0_time[reflevel]) < 1e-14)
+        CCTK_WARN(0, "Invalid past value times");
+
+    time_interp_step = lapse_prev1_time[reflevel] - lapse_prev0_time[reflevel];
+    fact0 = (lapse_prev1_time[reflevel] - t) / time_interp_step;
+    fact1 = (t - lapse_prev0_time[reflevel]) / time_interp_step;
+
+    if (reflevel < nb_levels - 1) {
+        /* on coarse levels use extrapolated past solutions */
+        fprintf(stderr, "extrapolating from t1=%g (%g) and t0=%g (%g)\n", lapse_prev1_time[reflevel], fact1, lapse_prev0_time[reflevel], fact0);
+
+        for (size_t i = 0; i < grid_size; i++)
+            lapse_mg_eval[i] = fact0 * lapse_prev0[i] + fact1 * lapse_prev1[i];
+    } else {
+        /* on the finest level, use the extrapolation only for the boundary
+         * values and solve for the interior */
+        CoordPatch *cp = get_coord_patch(ms, reflevel);
+
+        fill_eq_coeffs(ms->gh, cp);
+
+        fprintf(stderr, "extrapolating BCs from t1=%g (%g) and t0=%g (%g)\n", lapse_prev1_time[reflevel], fact1, lapse_prev0_time[reflevel], fact0);
+
+        for (ptrdiff_t idx_z = 0; idx_z < cp->offset_right[1]; idx_z++)
+            for (ptrdiff_t idx_x = 0; idx_x < cp->boundary_val_stride[0]; idx_x++) {
+                const ptrdiff_t idx = CCTK_GFINDEX3D(cctkGH, idx_x + cp->offset_left[0], cp->y_idx, cctkGH->cctk_lsh[2] - cp->offset_right[1] + idx_z);
+                lapse_mg_eval[idx] = fact0 * lapse_prev0[idx] + fact1 * lapse_prev1[idx];
+            }
+        for (int idx_z = 0; idx_z < cp->boundary_val_stride[1]; idx_z++)
+            for (int idx_x = 0; idx_x < cp->offset_right[0]; idx_x++) {
+                const ptrdiff_t idx = CCTK_GFINDEX3D(cctkGH, cctkGH->cctk_lsh[0] - cp->offset_right[0] + idx_x, cp->y_idx,
+                                                     cp->offset_left[1] + idx_z);
+                lapse_mg_eval[idx] = fact0 * lapse_prev0[idx] + fact1 * lapse_prev1[idx];
+            }
+
+        for (size_t i = 0; i < cp->solver->boundaries[MG2D_BOUNDARY_0U]->val_len; i++) {
+            const ptrdiff_t idx = CCTK_GFINDEX3D(cctkGH, i + cp->offset_left[0], cp->y_idx, cctkGH->cctk_lsh[2] - cp->offset_right[1]);
+            cp->solver->boundaries[MG2D_BOUNDARY_0U]->val[i] = lapse_mg_eval[idx] - 1.0;
+        }
+        for (size_t i = 0; i < cp->solver->boundaries[MG2D_BOUNDARY_1U]->val_len; i++) {
+            const ptrdiff_t idx = CCTK_GFINDEX3D(cctkGH, cctkGH->cctk_lsh[0] - cp->offset_right[0], cp->y_idx, cp->offset_left[1] + i);
+            cp->solver->boundaries[MG2D_BOUNDARY_1U]->val[i] = lapse_mg_eval[idx] - 1.0;
+        }
+
+        fprintf(stderr, "mg solve\n");
+
+        ret = mg2d_solve(cp->solver);
         if (ret < 0)
-            CCTK_WARN(0, "Error interpolating");
+            CCTK_WARN(0, "Error solving the maximal slicing equation");
 
-        for (size_t i = 0; i < cp->solver->boundaries[MG2D_BOUNDARY_0U]->val_len; i++)
-            cp->solver->boundaries[MG2D_BOUNDARY_0U]->val[i] = cp->boundary_val[0][i] - 1.0;
-        for (size_t i = 0; i < cp->solver->boundaries[MG2D_BOUNDARY_1U]->val_len; i++)
-            cp->solver->boundaries[MG2D_BOUNDARY_1U]->val[i] = cp->boundary_val[1][i] - 1.0;
+        solution_to_grid(cp, lapse_mg_eval);
     }
+    memcpy(alpha, lapse_mg_eval, sizeof(*alpha) * grid_size);
 }
 
-void maximal_slicing_axi_mg(CCTK_ARGUMENTS)
+void msa_mg_solve(CCTK_ARGUMENTS)
 {
-    static MSMGContext *ms;
+    DECLARE_CCTK_ARGUMENTS;
+    DECLARE_CCTK_PARAMETERS;
 
     CoordPatch *cp;
 
+    int nb_levels, reflevel_top, ts_tmp;
+    int nb_levels_type, ret;
+
+    const size_t grid_size = cctkGH->cctk_lsh[2] * cctkGH->cctk_lsh[1] * cctkGH->cctk_lsh[0];
+    const double t = cctkGH->cctk_time;
+    const int timestep = lrint(t * cctkGH->cctk_levfac[0] / cctkGH->cctk_delta_time);
+    const int reflevel = ctz(cctkGH->cctk_levfac[0]);
+
+    double *lapse_mg1;
+
+    fprintf(stderr, "solve lapse at rl=%d, t=%g; step %d\n", reflevel, t, timestep);
+
+    nb_levels = *(int*)CCTK_ParameterGet("num_levels_1", "CarpetRegrid2", &nb_levels_type);
+
+    reflevel_top = reflevel;
+    ts_tmp = timestep;
+    while (reflevel_top > 0 && !(ts_tmp % 2)) {
+        ts_tmp /= 2;
+        reflevel_top--;
+    }
+
+    if (!history_initialized(lapse_prev0_time, lapse_prev1_time, nb_levels)) {
+        if (reflevel == 0 || !(timestep % 2)) {
+            if (lapse_prev1_time[reflevel] != DBL_MAX &&
+                reflevel_top == 0) {
+                fprintf(stderr, "copy %g\n", lapse_prev1_time[reflevel]);
+                memcpy(lapse_prev0, lapse_prev1, grid_size * sizeof(*lapse_prev1));
+                lapse_prev0_time[reflevel] = lapse_prev1_time[reflevel];
+            }
+
+            fprintf(stderr, "solve ps\n");
+            solve_ps(ms, alpha);
+            memcpy(lapse_prev1, alpha, grid_size * sizeof(*lapse_prev0));
+            lapse_prev1_time[reflevel] = t;
+        }
+        return;
+    }
+
+    fprintf(stderr, "mg solve cur %d top %d\n", reflevel, reflevel_top);
+
+    /* fill in the equation coefficients */
+    cp = get_coord_patch(ms, reflevel);
+    fill_eq_coeffs(cctkGH, cp);
+
+    if (reflevel == 0 || timestep % 2) {
+        /* outer-most level for this solve, use extrapolated lapse as the
+         * boundary condition */
+        double time_interp_step = lapse_prev1_time[reflevel] - lapse_prev0_time[reflevel];
+        double fact0 = (lapse_prev1_time[reflevel] - t) / time_interp_step;
+        double fact1 = (t - lapse_prev0_time[reflevel]) / time_interp_step;
+
+        fprintf(stderr, "extrapolating BCs from t1=%g (%g) and t0=%g (%g)\n", lapse_prev1_time[reflevel], fact1, lapse_prev0_time[reflevel], fact0);
+
+        for (size_t i = 0; i < cp->solver->boundaries[MG2D_BOUNDARY_0U]->val_len; i++) {
+            const ptrdiff_t idx = CCTK_GFINDEX3D(cctkGH, i + cp->offset_left[0], cp->y_idx, cctkGH->cctk_lsh[2] - cp->offset_right[1]);
+            lapse_mg[idx] = fact0 * lapse_prev0[idx] + fact1 * lapse_prev1[idx];
+        }
+        for (size_t i = 0; i < cp->solver->boundaries[MG2D_BOUNDARY_1U]->val_len; i++) {
+            const ptrdiff_t idx = CCTK_GFINDEX3D(cctkGH, cctkGH->cctk_lsh[0] - cp->offset_right[0], cp->y_idx, cp->offset_left[1] + i);
+            lapse_mg[idx] = fact0 * lapse_prev0[idx] + fact1 * lapse_prev1[idx];
+        }
+    }
+    /* if the condition above was false, then lapse_mg should be filled by
+     * prolongation from the coarser level */
+
+    /* fill the solver boundary conditions */
+    for (size_t i = 0; i < cp->solver->boundaries[MG2D_BOUNDARY_0U]->val_len; i++) {
+        const ptrdiff_t idx = CCTK_GFINDEX3D(cctkGH, i + cp->offset_left[0], cp->y_idx, cctkGH->cctk_lsh[2] - cp->offset_right[1]);
+        cp->solver->boundaries[MG2D_BOUNDARY_0U]->val[i] = lapse_mg[idx] - 1.0;
+    }
+    for (size_t i = 0; i < cp->solver->boundaries[MG2D_BOUNDARY_1U]->val_len; i++) {
+        const ptrdiff_t idx = CCTK_GFINDEX3D(cctkGH, cctkGH->cctk_lsh[0] - cp->offset_right[0], cp->y_idx, cp->offset_left[1] + i);
+        cp->solver->boundaries[MG2D_BOUNDARY_1U]->val[i] = lapse_mg[idx] - 1.0;
+    }
+
+    /* do the elliptic solve */
+    ret = mg2d_solve(cp->solver);
+    if (ret < 0)
+        CCTK_WARN(0, "Error solving the maximal slicing equation");
+
+    /* export the result */
+    solution_to_grid(cp, lapse_mg);
+
+    lapse_mg1 = CCTK_VarDataPtr(cctkGH, 1, "MaximalSlicingAxiMG::lapse_mg");
+    memcpy(lapse_mg1, lapse_mg, grid_size * sizeof(*lapse_mg1));
+
+    /* update lapse history for extrapolation */
+    if (reflevel == 0 || !(timestep % 2)) {
+        const double vel_fact = 1.0 / (1 << (reflevel - reflevel_top));
+
+        for (size_t j = 0; j < grid_size; j++) {
+            const double sol_new = lapse_mg[j];
+            const double delta = sol_new - lapse_mg_eval[j];
+            lapse_prev0[j] = lapse_prev1[j] + delta - delta * vel_fact;
+            lapse_prev1[j] = sol_new;
+            alpha[j]       = sol_new;
+        }
+        lapse_prev0_time[reflevel] = lapse_prev1_time[reflevel];
+        lapse_prev1_time[reflevel] = cctkGH->cctk_time;
+    }
+}
+
+void msa_mg_sync(CCTK_ARGUMENTS)
+{
+    DECLARE_CCTK_ARGUMENTS;
+    DECLARE_CCTK_PARAMETERS;
+    const int reflevel = ctz(cctkGH->cctk_levfac[0]);
+    fprintf(stderr, "\nsync %g %d\n\n", cctkGH->cctk_time, reflevel);
+
+    return;
+}
+
+#if 0
+static void maximal_slicing_axi_mg_old(CCTK_ARGUMENTS)
+{
+    CoordPatch *cp;
+
     DECLARE_CCTK_ARGUMENTS;
     DECLARE_CCTK_PARAMETERS;
 
     double *coeffs = NULL;
     int i, ret;
 
-    /* on the first run, init the solver */
-    if (!ms)
-        context_init(cctkGH, &ms);
-
     fprintf(stderr, "ms mg solve: level %d time %g\n", ctz(ms->gh->cctk_levfac[0]), ms->gh->cctk_time);
 
     cp = get_coord_patch(ms);
@@ -499,3 +627,27 @@ void maximal_slicing_axi_mg(CCTK_ARGUMENTS)
         }
     CCTK_TimerStop("MaximalSlicingAxiMG_Copy");
 }
+#endif
+
+void msa_mg_init(CCTK_ARGUMENTS)
+{
+    DECLARE_CCTK_ARGUMENTS;
+    DECLARE_CCTK_PARAMETERS;
+
+    context_init(cctkGH, &ms);
+}
+
+void msa_mg_inithist(CCTK_ARGUMENTS)
+{
+    DECLARE_CCTK_ARGUMENTS;
+    DECLARE_CCTK_PARAMETERS;
+    const size_t grid_size = cctkGH->cctk_lsh[2] * cctkGH->cctk_lsh[1] * cctkGH->cctk_lsh[0];
+
+    for (size_t i = 0; i < grid_size; i++)
+        lapse_mg[i] = 1.0;
+
+    for (int i = 0; i < 32; i++) {
+        lapse_prev0_time[i] = DBL_MAX;
+        lapse_prev1_time[i] = DBL_MAX;
+    }
+}