Correct past solutions to prevent wild oscillations.

Same correction as in MaximalSlicingAxiMG.

1 files changed, 25 insertions, 4 deletions

diff --git a/src/qms.c b/src/qms.c
index 7f0f3e1..e64c560 100644
--- a/src/qms.c
+++ b/src/qms.c
@@ -1272,7 +1272,9 @@ void qms_mg_solve(CCTK_ARGUMENTS)
 
     const int reflevel = ctz(cctkGH->cctk_levfac[0]);
     const double time  = ms->gh->cctk_time;
+    const int timestep = lrint(time * cctkGH->cctk_levfac[0] / cctkGH->cctk_delta_time);
 
+    int reflevel_top, ts_tmp;
     int64_t grid_size;
     int64_t total_start, start;
     int ret;
@@ -1293,6 +1295,13 @@ void qms_mg_solve(CCTK_ARGUMENTS)
 
     grid_size = cctkGH->cctk_lsh[0] * cctkGH->cctk_lsh[1] * cctkGH->cctk_lsh[2];
 
+    reflevel_top = reflevel;
+    ts_tmp = timestep;
+    while (reflevel_top > 0 && !(ts_tmp % 2)) {
+        ts_tmp /= 2;
+        reflevel_top--;
+    }
+
     if (reflevel < ms->solve_level_max)
         goto skip_solve;
 
@@ -1368,11 +1377,23 @@ skip_solve:
 
     /* 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];
+        //memcpy(W_val0, W_val1, sizeof(*W_val0) * grid_size);
+        //memcpy(W_val1, W_val, sizeof(*W_val1) * grid_size);
+
+        const double vel_fact = 1.0 / (1 << (reflevel - reflevel_top));
+
+        qms_assert(fabs(time - W_pred1_time[reflevel]) < 1e-13);
 
-        memcpy(W_val1, W_val, sizeof(*W_val1) * grid_size);
-        W_val1_time[reflevel] = ms->gh->cctk_time;
+#pragma omp parallel for
+        for (int i = 0; i < grid_size; i++) {
+            const double sol_new = W_val[i];
+            const double delta   = sol_new - W_pred1[i];
+
+            W_val0[i] = W_val1[i] + delta - delta * vel_fact;
+            W_val1[i] = sol_new;
+        }
+        W_val0_time[reflevel] = W_val1_time[reflevel];
+        W_val1_time[reflevel] = time;
     }
 
     /* linearly extrapolate the past two solution to predict the next one */