From dade53184986b1993b560d54ecbea14a13787f0b Mon Sep 17 00:00:00 2001 From: Anton Khirnov Date: Fri, 30 Jun 2017 10:49:39 +0200 Subject: Ensure the solution is continuous in time. --- src/qms.c | 88 ++++++++++++++++++++++++++++++++++++++++++--------------------- 1 file 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; -- cgit v1.2.3