Make recovery from checkpoints work seamlessly.

3 files changed, 28 insertions, 37 deletions

diff --git a/interface.ccl b/interface.ccl
index 324fc5f..7216b74 100644
--- a/interface.ccl
+++ b/interface.ccl
@@ -13,4 +13,5 @@ REQUIRES FUNCTION MoLRegisterSaveAndRestoreGroup
 
 public:
 CCTK_REAL W TYPE=GF
-CCTK_REAL w_coeffs TYPE=array DIM=2 SIZE=basis_order_z,basis_order_r DISTRIB=constant
+CCTK_REAL W_coeffs TYPE=array DIM=3 SIZE=2,basis_order_z,basis_order_r DISTRIB=constant
+CCTK_REAL W_pred   TYPE=array DIM=3 SIZE=2,basis_order_z,basis_order_r DISTRIB=constant
diff --git a/schedule.ccl b/schedule.ccl
index b4fcda2..d8eb719 100644
--- a/schedule.ccl
+++ b/schedule.ccl
@@ -26,7 +26,5 @@ SCHEDULE qms_init IN ADMBase_InitialData {
 } ""
 
 STORAGE: W
-
-if (export_coeffs) {
-    STORAGE: w_coeffs
-}
+STORAGE: W_pred
+STORAGE: W_coeffs
diff --git a/src/qms.c b/src/qms.c
index 3d3fb64..4be70c6 100644
--- a/src/qms.c
+++ b/src/qms.c
@@ -256,24 +256,16 @@ static int context_init(cGH *cctkGH)
         return -ENOMEM;
 
     for (int i = 0; i < ARRAY_ELEMS(qms->solutions); i++) {
-        ret = posix_memalign((void**)&qms->solutions[i].coeffs, 32,
-                             basis_order_r * basis_order_z * sizeof(*qms->solutions[i].coeffs));
-        if (ret)
-            return -ENOMEM;
-
-        memset(qms->solutions[i].coeffs, 0, basis_order_r * basis_order_z * sizeof(*qms->solutions[i].coeffs));
-        qms->solutions[i].time = -i * cctkGH->cctk_delta_time;
+        qms->solutions[i].coeffs = W_coeffs + i * basis_order_r * basis_order_z;
+        //memset(qms->solutions[i].coeffs, 0, basis_order_r * basis_order_z * sizeof(*qms->solutions[i].coeffs));
+        qms->solutions[i].time = cctkGH->cctk_time - (2 - i) * cctkGH->cctk_delta_time / (1.0 * (1 << solve_level));
         fprintf(stderr, "init time %d %f\n", i, qms->solutions[i].time);
     }
 
     for (int i = 0; i < ARRAY_ELEMS(qms->pred); i++) {
-        ret = posix_memalign((void**)&qms->pred[i].coeffs, 32,
-                             basis_order_r * basis_order_z * sizeof(*qms->pred[i].coeffs));
-        if (ret)
-            return -ENOMEM;
-
-        memset(qms->pred[i].coeffs, 0, basis_order_r * basis_order_z * sizeof(*qms->pred[i].coeffs));
-        qms->pred[i].time = -i * cctkGH->cctk_delta_time;
+        qms->pred[i].coeffs = W_pred + i * basis_order_r * basis_order_z;
+        //memset(qms->pred[i].coeffs, 0, basis_order_r * basis_order_z * sizeof(*qms->pred[i].coeffs));
+        qms->pred[i].time = cctkGH->cctk_time - (1 - i) * cctkGH->cctk_delta_time / (1.0 * (1 << solve_level));
         fprintf(stderr, "init time %d %f\n", i, qms->pred[i].time);
     }
 
@@ -316,30 +308,33 @@ void quasimaximal_slicing_axi_solve(CCTK_ARGUMENTS)
 
     /* add the solution to the list of past solutions */
     {
-        double *tmp = ms->solutions[0].coeffs;
+        int N = ms->solver->nb_coeffs[0] * ms->solver->nb_coeffs[1];
 
-        memmove(ms->solutions, ms->solutions + 1, sizeof(ms->solutions[0]) * (ARRAY_ELEMS(ms->solutions) - 1));
-        ms->solutions[ARRAY_ELEMS(ms->solutions) - 1].coeffs = ms->solver->coeffs;
-        ms->solutions[ARRAY_ELEMS(ms->solutions) - 1].time   = ms->gh->cctk_time;
+        memcpy(ms->solutions[0].coeffs, ms->solutions[1].coeffs, sizeof(*ms->solutions[0].coeffs) * N);
+        ms->solutions[0].time = ms->solutions[1].time;
 
-        ms->solver->coeffs = tmp;
+        memcpy(ms->solutions[1].coeffs, ms->solver->coeffs, sizeof(*ms->solutions[1].coeffs) * N);
+        ms->solutions[1].time = ms->gh->cctk_time;
     }
 
     /* linearly extrapolate the past two solution to predict the next one */
     {
-        double time0 = ms->solutions[ARRAY_ELEMS(ms->solutions) - 2].time;
-        double time1 = ms->solutions[ARRAY_ELEMS(ms->solutions) - 1].time;
-        double time  = time1 + ms->gh->cctk_delta_time;
+        int N = ms->solver->nb_coeffs[0] * ms->solver->nb_coeffs[1];
 
-        double *coeffs0 = ms->solutions[ARRAY_ELEMS(ms->solutions) - 2].coeffs;;
-        double *coeffs1 = ms->solutions[ARRAY_ELEMS(ms->solutions) - 1].coeffs;;
-        double *pred = ms->pred[0].coeffs;
+        double time0 = ms->solutions[0].time;
+        double time1 = ms->solutions[1].time;
+        double time  = time1 + ms->gh->cctk_delta_time / (1.0 * (1 << solve_level));
 
-        for (int i = 0; i < ms->solver->nb_coeffs[0] * ms->solver->nb_coeffs[1]; i++)
-            pred[i] = (coeffs1[i] * (time - time0) / (time1 - time0) + coeffs0[i] * (time - time1) / (time0 - time1));
+        double *coeffs0 = ms->solutions[0].coeffs;
+        double *coeffs1 = ms->solutions[1].coeffs;
+        double *pred = ms->pred[1].coeffs;
+
+        memcpy(ms->pred[0].coeffs, ms->pred[1].coeffs, sizeof(*ms->pred[0].coeffs) * N);
+        ms->pred[0].time = ms->pred[1].time;
+
+        for (int i = 0; i < N; i++)
+            ms->pred[1].coeffs[i] = (coeffs1[i] * (time - time0) / (time1 - time0) + coeffs0[i] * (time - time1) / (time0 - time1));
 
-        memmove(ms->pred, ms->pred + 1, sizeof(ms->pred[0]) * (ARRAY_ELEMS(ms->pred) - 1));
-        ms->pred[1].coeffs = pred;
         ms->pred[1].time   = time;
     }
 }
@@ -397,9 +392,6 @@ void quasimaximal_slicing_axi_eval(CCTK_ARGUMENTS)
     }
 #endif
 
-    if (export_coeffs)
-        memcpy(w_coeffs, coeffs, sizeof(*w_coeffs) * ms->solver->nb_coeffs[0] * ms->solver->nb_coeffs[1]);
-
     CCTK_TimerStart("QuasiMaximalSlicing_Expand");
     expand_start = gettime();
 #if 0