From 1424606c79dd3fa71125a7fda37e34372104ac32 Mon Sep 17 00:00:00 2001 From: Anton Khirnov Date: Fri, 13 Sep 2019 20:41:28 +0200 Subject: Correct past solutions to prevent wild oscillations. Same correction as in MaximalSlicingAxiMG. --- src/qms.c | 29 +++++++++++++++++++++++++---- 1 file 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 */ -- cgit v1.2.3