Ensure the solution is continuous in time.

1 files changed, 59 insertions, 29 deletions

diff --git a/src/qms.c b/src/qms.c
index a40df08..761c3ae 100644
--- a/src/qms.c
+++ b/src/qms.c
@@ -46,11 +46,13 @@ struct QMSContext {
     struct {
         double time;
         double *coeffs;
-    } solution_cache[8];
-    int nb_solutions;
+    } solutions[2], pred[2];
 
     double *coeffs_eval;
 
+    double max_radius;
+    int solve_level;
+
     uint64_t grid_expand_count;
     uint64_t grid_expand_time;
 
@@ -74,6 +76,7 @@ static CoordPatch *get_coord_patch(QMSContext *ms,
 
     CoordPatch *cp;
     int64_t grid_size;
+    int i;
 
     for (int i = 0; i < ms->nb_patches; i++) {
         cp = &ms->patches[i];
@@ -138,14 +141,14 @@ static CoordPatch *get_coord_patch(QMSContext *ms,
             //    }
             for (int k = 0; k < ms->solver->nb_coeffs[0]; k++) {
                 double dx = calc_basis_freq(ms->solver->basis[0], k);
-                double r0 = dx * scale_factor;
+                double r0 = MIN(ms->max_radius, dx * scale_factor);
                 double fact =  exp(-36.0 * pow(rr / r0, scale_power));
 
                 cp->transform_matrix[idx_grid + cp->size[0] * cp->size[2] * k] = ms->solver->basis[0]->eval(coord0, k) * fact;
             }
             for (int k = 0; k < ms->solver->nb_coeffs[1]; k++) {
                 double dx = calc_basis_freq(ms->solver->basis[1], k);
-                double r0 = dx * scale_factor;
+                double r0 = MIN(ms->max_radius, dx * scale_factor);
                 double fact = exp(-36.0 * pow(rr / r0, scale_power));
 
                 cp->transform_matrix1[idx_grid * ms->solver->nb_coeffs[1] + k] = ms->solver->basis[1]->eval(coord1, k) * fact;
@@ -226,6 +229,8 @@ static int context_init(cGH *cctkGH)
 
     qms->gh = cctkGH;
 
+    qms->max_radius = 60.0;
+
     ret = qms_solver_init(&qms->solver, cctkGH, basis_order_r, basis_order_z,
                           scale_factor, filter_power, 0.0);
     if (ret < 0)
@@ -236,11 +241,26 @@ static int context_init(cGH *cctkGH)
     if (ret)
         return -ENOMEM;
 
-    for (int i = 0; i < ARRAY_ELEMS(qms->solution_cache); i++) {
-        ret = posix_memalign((void**)&qms->solution_cache[i].coeffs, 32,
-                             basis_order_r * basis_order_z * sizeof(*qms->solution_cache[i].coeffs));
+    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;
+        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;
+        fprintf(stderr, "init time %d %f\n", i, qms->pred[i].time);
     }
 
     qms_context = qms;
@@ -268,7 +288,7 @@ void quasimaximal_slicing_axi_solve(CCTK_ARGUMENTS)
         return;
 
     if (ms->gh->cctk_levfac[0] != 1 || fabs(time - ceilf(time)) > 1e-8 ||
-        (ms->nb_solutions && ms->solution_cache[ms->nb_solutions - 1].time == cctkGH->cctk_time))
+        (ms->solutions[ARRAY_ELEMS(ms->solutions) - 1].time >= cctkGH->cctk_time))
         return;
 
     CCTK_TimerStart("QuasiMaximalSlicing_Solve");
@@ -276,20 +296,35 @@ void quasimaximal_slicing_axi_solve(CCTK_ARGUMENTS)
     CCTK_TimerStop("QuasiMaximalSlicing_Solve");
 
     fprintf(stderr, "%d qms solve: time %g %g %g\n", ms->gh->cctk_levfac[0], ms->gh->cctk_time, time, ms->solver->coeffs[0]);
-    if (1) {
-        double *tmp;
-        if (ms->nb_solutions == ARRAY_ELEMS(ms->solution_cache)) {
-            tmp = ms->solution_cache[0].coeffs;
-            memmove(ms->solution_cache, ms->solution_cache + 1, sizeof(ms->solution_cache[0]) * (ARRAY_ELEMS(ms->solution_cache) - 1));
-        } else {
-            ms->nb_solutions++;
-            tmp = ms->solution_cache[ms->nb_solutions - 1].coeffs;
-        }
-        ms->solution_cache[ms->nb_solutions - 1].coeffs = ms->solver->coeffs;
-        ms->solution_cache[ms->nb_solutions - 1].time   = ms->gh->cctk_time;
+
+    /* add the solution to the list of past solutions */
+    {
+        double *tmp = ms->solutions[0].coeffs;
+
+        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;
 
         ms->solver->coeffs = tmp;
     }
+
+    /* 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;
+
+        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;
+
+        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));
+
+        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;
+    }
 }
 
 void quasimaximal_slicing_axi_eval(CCTK_ARGUMENTS)
@@ -326,14 +361,11 @@ void quasimaximal_slicing_axi_eval(CCTK_ARGUMENTS)
 #else
     coeffs = ms->coeffs_eval;
 
-    if (cctkGH->cctk_levfac[0] < 1 || ms->nb_solutions < 2) {
-        memset(coeffs, 0, sizeof(*coeffs) * ms->solver->nb_coeffs[0] * ms->solver->nb_coeffs[1]);
-        //fprintf(stderr, "qms eval: time %g zero\n", ms->gh->cctk_time);
-    } else {
-        double *coeffs0 = ms->solution_cache[ms->nb_solutions - 2].coeffs;
-        double *coeffs1 = ms->solution_cache[ms->nb_solutions - 1].coeffs;
-        double time0 = ms->solution_cache[ms->nb_solutions - 2].time;
-        double time1 = ms->solution_cache[ms->nb_solutions - 1].time;
+    {
+        double *coeffs0 = ms->pred[0].coeffs;
+        double *coeffs1 = ms->pred[1].coeffs;
+        double time0 = ms->pred[0].time;
+        double time1 = ms->pred[1].time;
 
         double fact;
 
@@ -342,8 +374,6 @@ void quasimaximal_slicing_axi_eval(CCTK_ARGUMENTS)
         else
             fact = exp(-36.0 * pow((time - switchoff_time), 4.0));
 
-        //fprintf(stderr, "qms eval: time %g interp from %g %g %g\n", ms->gh->cctk_time, time0, time1, fact);
-
         for (int i = 0; i < ms->solver->nb_coeffs[0] * ms->solver->nb_coeffs[1]; i++)
             coeffs[i] = (coeffs1[i] * (time - time0) / (time1 - time0) + coeffs0[i] * (time - time1) / (time0 - time1)) * fact;