Initial commit.

7 files changed, 573 insertions, 0 deletions

diff --git a/configuration.ccl b/configuration.ccl
new file mode 100644
index 0000000..f124283
--- /dev/null
+++ b/configuration.ccl
@@ -0,0 +1,2 @@
+# Configuration definition for thorn MaximalSlicingAxiMG
+
diff --git a/interface.ccl b/interface.ccl
new file mode 100644
index 0000000..43e512c
--- /dev/null
+++ b/interface.ccl
@@ -0,0 +1,15 @@
+# Interface definition for thorn MaximalSlicingAxiMG
+implements: MaximalSlicingAxiMG
+
+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)
+
+CCTK_INT FUNCTION MoLNumIntegratorSubsteps ()
+USES FUNCTION MoLNumIntegratorSubsteps
+
+REQUIRES FUNCTION MoLRegisterConstrained
+REQUIRES FUNCTION MoLRegisterSaveAndRestore
+REQUIRES FUNCTION MoLRegisterSaveAndRestoreGroup
diff --git a/param.ccl b/param.ccl
new file mode 100644
index 0000000..e6b39e2
--- /dev/null
+++ b/param.ccl
@@ -0,0 +1,10 @@
+# Parameter definitions for thorn MaximalSlicingAxiMG
+
+SHARES: ADMBase
+
+EXTENDS KEYWORD lapse_evolution_method
+{
+  "maximal_axi_mg" :: "Maximal slicing for an axisymmetric spacetime"
+}
+
+RESTRICTED:
diff --git a/schedule.ccl b/schedule.ccl
new file mode 100644
index 0000000..246aa63
--- /dev/null
+++ b/schedule.ccl
@@ -0,0 +1,24 @@
+# Schedule definitions for thorn MaximalSlicingAxiMG
+#
+if (CCTK_Equals(lapse_evolution_method, "maximal_axi_mg")) {
+
+    SCHEDULE maximal_slicing_axi_mg IN ML_BSSN_evolCalcGroup BEFORE ML_BSSN_RHS {
+        LANG: C
+    } "Maximal slicing in axisymmetry"
+
+    #SCHEDULE maximal_slicing_axi IN MoL_PostStepModify {
+    #    LANG: C
+    #} "Maximal slicing in axisymmetry"
+
+    #SCHEDULE maximal_slicing_axi IN MoL_PseudoEvolution {
+    #    LANG: C
+    #} "Maximal slicing in axisymmetry"
+
+    SCHEDULE maximal_slicing_axi_mg_register AT Startup {
+        LANG: C
+    } ""
+
+    SCHEDULE maximal_slicing_axi_mg_register_mol IN MoL_Register {
+        LANG: C
+    } ""
+}
diff --git a/src/make.code.defn b/src/make.code.defn
new file mode 100644
index 0000000..4bc0d9c
--- /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 = maximal_slicing_axi_mg.c register.c
+
+# Subdirectories containing source files
+SUBDIRS = 
diff --git a/src/maximal_slicing_axi_mg.c b/src/maximal_slicing_axi_mg.c
new file mode 100644
index 0000000..668a702
--- /dev/null
+++ b/src/maximal_slicing_axi_mg.c
@@ -0,0 +1,497 @@
+#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 <cblas.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 ARRAY_ELEMS(x) (sizeof(x) / sizeof(*x))
+
+/* precomputed values for a given refined grid */
+typedef struct CoordPatch {
+    double origin[3];
+    int    delta[3];
+    int    size[3];
+    int    level;
+
+    MG2DContext *solver;
+
+    ptrdiff_t y_idx;
+
+    /* number of x/z grid points by which the elliptic solver domain is offset
+     * from the cactus grid */
+    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
+     * [0] is the upper-z patch
+     *     rows run from x=0 to x[lsh[0] - 1]
+     * [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 {
+    cGH *gh;
+
+    CoordPatch *patches;
+    int nb_patches;
+} MSMGContext;
+
+static int ctz(int a)
+{
+    int ret = 0;
+
+    if (!a)
+        return INT_MAX;
+
+    while (!(a & 1)) {
+        a >>= 1;
+        ret++;
+    }
+
+    return ret;
+}
+
+static void coord_patch_free(CoordPatch *cp)
+{
+    mg2d_solver_free(&cp->solver);
+}
+
+static CoordPatch *get_coord_patch(MSMGContext *ms)
+{
+    cGH *gh = ms->gh;
+
+    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");
+
+    CoordPatch *cp;
+    size_t domain_size[2];
+    int integrator_substeps = MoLNumIntegratorSubsteps();
+    int i, ret;
+
+    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])
+            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));
+
+    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]);
+
+    for (i = 0; i < cp->size[1]; i++)
+        if (fabs(a_y[CCTK_GFINDEX3D(gh, 0, i, 0)]) < 1e-8) {
+            cp->y_idx = i;
+            break;
+        }
+    if (i == cp->size[1])
+        CCTK_WARN(0, "The grid does not include y==0");
+
+    for (int i = 0; i < 2; i++) {
+        ptrdiff_t offset_left = ms->gh->cctk_nghostzones[2 * i];
+        ptrdiff_t offset_right = offset_left;
+
+        /* account for grid tapering on refined levels */
+        if (cp->level > 0)
+            offset_right *= integrator_substeps * 2;
+
+        domain_size[i] = ms->gh->cctk_lsh[2 * i] - offset_left - offset_right;
+        /* the outer boundary layer overlaps with the first ghost zone */
+        if (cp->level > 0)
+            domain_size[i]++;
+
+        cp->offset_left[i]  = offset_left;
+        cp->offset_right[i] = offset_right;
+    }
+
+    if (domain_size[0] != domain_size[1])
+        CCTK_WARN(0, "The grid is non-square, only square grids are supported");
+
+    cp->solver = mg2d_solver_alloc(domain_size[0]);
+    if (!cp->solver)
+        CCTK_WARN(0, "Error allocating the solver");
+
+    cp->solver->step[0] = a_x[1] - a_x[0];
+    cp->solver->step[1] = cp->solver->step[0];
+
+    cp->solver->fd_stencil = 1;
+
+    cp->solver->boundaries[MG2D_BOUNDARY_0L]->type = MG2D_BC_TYPE_FIXDIFF;
+    cp->solver->boundaries[MG2D_BOUNDARY_1L]->type = MG2D_BC_TYPE_FIXDIFF;
+    cp->solver->boundaries[MG2D_BOUNDARY_0U]->type = MG2D_BC_TYPE_FIXVAL;
+    cp->solver->boundaries[MG2D_BOUNDARY_1U]->type = MG2D_BC_TYPE_FIXVAL;
+
+    cp->solver->maxiter   = 32;
+    cp->solver->tol       = 1e-10;
+    cp->solver->nb_cycles = 2;
+
+    /* 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++)
+        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_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++;
+    return cp;
+}
+
+static int context_init(cGH *gh, MSMGContext **ctx)
+{
+    MSMGContext *ms;
+    int ret;
+
+    ms = calloc(1, sizeof(*ms));
+    if (!ms)
+        return -ENOMEM;
+
+    ms->gh = gh;
+
+    *ctx = ms;
+
+    return 0;
+}
+
+static void context_free(MSMGContext **pms)
+{
+    MSMGContext *ms = *pms;
+
+    if (!ms)
+        return;
+
+    for (int i = 0; i < ms->nb_patches; i++)
+        coord_patch_free(&ms->patches[i]);
+    free(ms->patches);
+
+    free(ms);
+    *pms = NULL;
+}
+
+static void fill_eq_coeffs(cGH *gh, CoordPatch *cp)
+{
+    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];
+
+    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_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");
+
+    for (int idx_z = 0; idx_z < cp->solver->domain_size; idx_z++)
+        for (int idx_x = 0; idx_x < cp->solver->domain_size; idx_x++) {
+            const int idx_src = CCTK_GFINDEX3D(gh, idx_x + gh->cctk_nghostzones[0], cp->y_idx, idx_z + gh->cctk_nghostzones[2]);
+            const int idx_dc  = idx_z * cp->solver->diff_coeffs_stride + idx_x;
+            const int idx_rhs = idx_z * cp->solver->rhs_stride + idx_x;
+
+            const double x = a_x[idx_src];
+
+            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 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  phi =  a_phi[idx_src];
+            const double  Xtx =  a_Xtx[idx_src];
+            const double  Xtz =  a_Xtz[idx_src];
+
+            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);
+
+            const double det = gtxx * gtyy * gtzz + 2 * gtxy * gtyz * gtxz - gtzz * SQR(gtxy) - SQR(gtxz) * gtyy - gtxx * SQR(gtyz);
+
+            const double At[3][3] = {
+                { Atxx, Atxy, Atxz },
+                { Atxy, Atyy, Atyz },
+                { Atxz, Atyz, Atzz }};
+
+
+            double Xx, Xz, k2;
+
+            double Am[3][3], gtu[3][3];
+
+			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];
+
+			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 += gtu[j][l] * At[l][k];
+					Am[j][k] = val;
+				}
+
+			// K_{ij} K^{ij}
+			k2 = 0.0;
+			for (int j = 0; j < 3; j++)
+				for (int k = 0; k < 3; k++)
+					k2 += Am[j][k] * Am[k][j];
+
+            Xx = SQR(phi) * (Xtx + (phi_dx * gtu[0][0] + phi_dz * gtu[0][2]) / phi);
+            Xz = SQR(phi) * (Xtz + (phi_dx * gtu[0][2] + phi_dz * gtu[2][2]) / phi);
+
+            cp->solver->diff_coeffs[MG2D_DIFF_COEFF_20][idx_dc] = SQR(phi) * (gtu[0][0] + ((idx_x == 0) ? gtu[1][1] : 0.0));
+            cp->solver->diff_coeffs[MG2D_DIFF_COEFF_02][idx_dc] = SQR(phi) * gtu[2][2];
+            cp->solver->diff_coeffs[MG2D_DIFF_COEFF_11][idx_dc] = SQR(phi) * gtu[0][2] * 2;
+            cp->solver->diff_coeffs[MG2D_DIFF_COEFF_10][idx_dc] = -Xx + (idx_x ? SQR(phi) * gtu[1][1] / x : 0.0);
+            cp->solver->diff_coeffs[MG2D_DIFF_COEFF_01][idx_dc] = -Xz;
+            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);
+        if (ret < 0)
+            CCTK_WARN(0, "Error interpolating");
+
+        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;
+    }
+}
+
+void maximal_slicing_axi_mg(CCTK_ARGUMENTS)
+{
+    static MSMGContext *ms;
+
+    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);
+
+    fill_eq_coeffs(ms->gh, cp);
+
+    CCTK_TimerStart("MaximalSlicingAxiMG_Solve");
+    ret = mg2d_solve(cp->solver);
+    if (ret < 0)
+        CCTK_WARN(0, "Error solving the maximal slicing equation");
+    CCTK_TimerStop("MaximalSlicingAxiMG_Solve");
+
+
+    CCTK_TimerStart("MaximalSlicingAxiMG_Copy");
+//#pragma omp parallel for
+    for (int j = 0; j < cp->solver->domain_size; j++)
+        for (int i = 0; i < cp->solver->domain_size; i++) {
+            const int idx_dst = CCTK_GFINDEX3D(ms->gh, i + cp->offset_left[0], cp->y_idx, j + cp->offset_left[1]);
+            const int idx_src = j * cp->solver->u_stride + i;
+            alpha[idx_dst] = 1.0 + cp->solver->u[idx_src];
+        }
+    if (cp->level > 0) {
+        /* fill in the interpolated outer buffer points */
+        for (int idx_z = 0; idx_z < cp->offset_right[1]; idx_z++)
+            for (int idx_x = 0; idx_x < cp->boundary_val_stride[0]; idx_x++) {
+                const ptrdiff_t idx_src = idx_z * cp->boundary_val_stride[0] + idx_x;
+                const ptrdiff_t idx_dst = CCTK_GFINDEX3D(ms->gh, cp->offset_left[0] + idx_x, cp->y_idx, ms->gh->cctk_lsh[2] - cp->offset_right[1] + idx_z);
+                alpha[idx_dst] = cp->boundary_val[0][idx_src];
+            }
+        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_src = idx_x * cp->boundary_val_stride[1] + idx_z;
+                const ptrdiff_t idx_dst = CCTK_GFINDEX3D(ms->gh, ms->gh->cctk_lsh[0] - cp->offset_right[0] + idx_x, cp->y_idx,
+                                                         cp->offset_left[1] + idx_z);
+                alpha[idx_dst] = cp->boundary_val[1][idx_src];
+            }
+    }
+
+    /* fill in the axial symmetry ghostpoints by mirroring */
+    for (int idx_z = cp->offset_left[1]; idx_z < ms->gh->cctk_lsh[2]; idx_z++)
+        for (int idx_x = 0; idx_x < cp->offset_left[0]; idx_x++) {
+            const ptrdiff_t idx_dst = CCTK_GFINDEX3D(ms->gh, idx_x, cp->y_idx, idx_z);
+            const ptrdiff_t idx_src = CCTK_GFINDEX3D(ms->gh, 2 * cp->offset_left[0] - idx_x, cp->y_idx, idx_z);
+            alpha[idx_dst] = alpha[idx_src];
+        }
+    for (int idx_z = 0; idx_z < cp->offset_left[1]; idx_z++)
+        for (int idx_x = 0; idx_x < ms->gh->cctk_lsh[0]; idx_x++) {
+            const ptrdiff_t idx_dst = CCTK_GFINDEX3D(ms->gh, idx_x, cp->y_idx, idx_z);
+            const ptrdiff_t idx_src = CCTK_GFINDEX3D(ms->gh, idx_x, cp->y_idx, 2 * cp->offset_left[1] - idx_z);
+            alpha[idx_dst] = alpha[idx_src];
+        }
+    CCTK_TimerStop("MaximalSlicingAxiMG_Copy");
+}
diff --git a/src/register.c b/src/register.c
new file mode 100644
index 0000000..9bc64cd
--- /dev/null
+++ b/src/register.c
@@ -0,0 +1,18 @@
+
+#include "cctk.h"
+#include "cctk_Arguments.h"
+#include "Slicing.h"
+
+int maximal_slicing_axi_mg_register(void)
+{
+    Einstein_RegisterSlicing("maximal_axi_mg");
+    return 0;
+}
+
+void maximal_slicing_axi_mg_register_mol(CCTK_ARGUMENTS)
+{
+    MoLRegisterConstrained(CCTK_VarIndex("ADMBase::alp"));
+
+    MoLRegisterSaveAndRestoreGroup(CCTK_GroupIndex("ADMBase::metric"));
+    MoLRegisterSaveAndRestoreGroup(CCTK_GroupIndex("ADMBase::curv"));
+}