From e4fd9b1001d9597abc47150475e2ba63fbf0a8e2 Mon Sep 17 00:00:00 2001 From: Anton Khirnov Date: Sat, 16 Feb 2019 11:34:06 +0100 Subject: Drop the use of the pseudospectral solver. It does not produce good enough results. --- src/maximal_slicing_axi_mg.c | 108 ++++++++++--------------------------------- 1 file changed, 24 insertions(+), 84 deletions(-) diff --git a/src/maximal_slicing_axi_mg.c b/src/maximal_slicing_axi_mg.c index c83c2fe..db57979 100644 --- a/src/maximal_slicing_axi_mg.c +++ b/src/maximal_slicing_axi_mg.c @@ -47,8 +47,6 @@ static inline int64_t gettime(void) return (int64_t)tv.tv_sec * 1000000 + tv.tv_usec; } -typedef struct MaximalSlicingContext MaximalSlicingContext; - /* precomputed values for a given refined grid */ typedef struct CoordPatch { int level; @@ -88,17 +86,12 @@ typedef struct MSMGContext { CoordPatch *patches; int nb_patches; - MaximalSlicingContext *ps_solver; - int log_level; /* timings */ int64_t time_solve; int64_t count_solve; - int64_t time_solve_ps; - int64_t count_solve_ps; - int64_t time_solve_mg; int64_t count_solve_mg; int64_t time_solve_fill_coeffs; @@ -110,9 +103,6 @@ typedef struct MSMGContext { int64_t time_eval; int64_t count_eval; - int64_t time_eval_ps; - int64_t count_eval_ps; - int64_t time_extrapolate; int64_t count_extrapolate; @@ -124,11 +114,6 @@ typedef struct MSMGContext { int64_t time_fine_export; } MSMGContext; -int msa_context_init(cGH *gh, MaximalSlicingContext **ctx, int basis_order0, int basis_order1, double outer_bound, double filter_power); -void msa_context_free(MaximalSlicingContext **ctx); -void msa_lapse_solve(MaximalSlicingContext *ctx); -void msa_lapse_eval(MaximalSlicingContext *ctx, double *dst); - static MSMGContext *ms; static int ctz(int a) @@ -282,9 +267,6 @@ static void print_stats(MSMGContext *ms) fprintf(stderr, "%2.2f%% eval: %ld runs; %g s total; %g ms avg per run\n", (double)ms->time_eval * 100.0 / total, ms->count_eval, ms->time_eval / 1e6, ms->time_eval / 1e3 / ms->count_eval); - fprintf(stderr, " %2.2f%% eval ps: %ld runs; %g s total; %g ms avg per run\n", - (double)ms->time_eval_ps * 100.0 / total, ms->count_eval_ps, - ms->time_eval_ps / 1e6, ms->time_eval_ps / 1e3 / ms->count_eval_ps); fprintf(stderr, " %2.2f%% eval extrapolate: %ld runs; %g s total; %g ms avg per run\n", (double)ms->time_extrapolate * 100.0 / total, ms->count_extrapolate, ms->time_extrapolate / 1e6, ms->time_extrapolate / 1e3 / ms->count_extrapolate); @@ -300,9 +282,6 @@ static void print_stats(MSMGContext *ms) fprintf(stderr, "%2.2f%% solve: %ld runs; %g s total; %g ms avg per run\n", (double)ms->time_solve * 100.0 / total, ms->count_solve, ms->time_solve / 1e6, ms->time_solve / 1e3 / ms->count_solve); - fprintf(stderr, " %2.2f%% solve ps: %ld runs; %g s total; %g ms avg per run\n", - (double)ms->time_solve_ps * 100.0 / total, ms->count_solve_ps, - ms->time_solve_ps / 1e6, ms->time_solve_ps / 1e3 / ms->count_solve_ps); fprintf(stderr, " %2.2f%% solve mg2d: %ld runs; %g s total; %g ms avg per run || " "%2.2f%% fill coeffs %2.2f%% boundaries %2.2f%% mg2d %2.2f%% export %2.2f%% history\n", (double)ms->time_solve_mg * 100.0 / total, ms->count_solve_mg, @@ -532,27 +511,6 @@ static void solution_to_grid(CoordPatch *cp, double *dst) } } -static void solve_ps(MSMGContext *ctx, double *dst) -{ - if (!ctx->ps_solver) { - const double *a_x = CCTK_VarDataPtr(ctx->gh, 0, "grid::x"); - const double outer_bound = a_x[CCTK_GFINDEX3D(ctx->gh, ctx->gh->cctk_lsh[0] - 1, 0, 0)]; - msa_context_init(ctx->gh, &ctx->ps_solver, 80, 80, outer_bound, 64.0); - } - - msa_lapse_solve(ctx->ps_solver); - msa_lapse_eval(ctx->ps_solver, dst); -} - -static int history_initialized(double *lapse_prev0_time, double *lapse_prev1_time, int nb_levels) -{ - for (int i = 0; i < nb_levels; i++) { - if (lapse_prev0_time[i] == DBL_MAX || lapse_prev1_time[i] == DBL_MAX) - return 0; - } - return 1; -} - void msa_mg_eval(CCTK_ARGUMENTS) { DECLARE_CCTK_ARGUMENTS; @@ -574,27 +532,20 @@ void msa_mg_eval(CCTK_ARGUMENTS) nb_levels = *(int*)CCTK_ParameterGet("num_levels_1", "CarpetRegrid2", &nb_levels_type); - if (!history_initialized(lapse_prev0_time, lapse_prev1_time, nb_levels)) { - int64_t eval_ps_start = gettime(); - - solve_ps(ms, lapse_mg_eval); - memcpy(alpha, lapse_mg_eval, sizeof(*alpha) * grid_size); - - ms->time_eval_ps += gettime() - eval_ps_start; - ms->count_eval_ps++; + time_interp_step = lapse_prev1_time[reflevel] - lapse_prev0_time[reflevel]; - LOGDEBUG( "pseudospectral solve\n"); - goto finish; + if (lapse_prev0_time[reflevel] == DBL_MAX && + lapse_prev1_time[reflevel] == DBL_MAX) { + fact0 = 0.0; + fact1 = 0.0; + } else if (lapse_prev0_time[reflevel] == DBL_MAX) { + fact0 = 0.0; + fact1 = 1.0; + } else { + fact0 = (lapse_prev1_time[reflevel] - t) / time_interp_step; + fact1 = (t - lapse_prev0_time[reflevel]) / time_interp_step; } - if (lapse_prev1_time[reflevel] <= lapse_prev0_time[reflevel] || - fabs(lapse_prev1_time[reflevel] - lapse_prev0_time[reflevel]) < 1e-14) - CCTK_WARN(0, "Invalid past value times"); - - time_interp_step = lapse_prev1_time[reflevel] - lapse_prev0_time[reflevel]; - fact0 = (lapse_prev1_time[reflevel] - t) / time_interp_step; - fact1 = (t - lapse_prev0_time[reflevel]) / time_interp_step; - if (reflevel < nb_levels - 1) { /* on coarse levels use extrapolated past solutions */ int64_t extrap_start = gettime(); @@ -704,28 +655,6 @@ void msa_mg_solve(CCTK_ARGUMENTS) reflevel_top--; } - if (!history_initialized(lapse_prev0_time, lapse_prev1_time, nb_levels)) { - if (reflevel == 0 || !(timestep % 2)) { - int64_t solve_ps_start = gettime(); - - if (lapse_prev1_time[reflevel] != DBL_MAX && - reflevel_top == 0) { - LOGDEBUG( "copy %g\n", lapse_prev1_time[reflevel]); - memcpy(lapse_prev0, lapse_prev1, grid_size * sizeof(*lapse_prev1)); - lapse_prev0_time[reflevel] = lapse_prev1_time[reflevel]; - } - - LOGDEBUG( "solve ps\n"); - solve_ps(ms, alpha); - memcpy(lapse_prev1, alpha, grid_size * sizeof(*lapse_prev0)); - lapse_prev1_time[reflevel] = t; - - ms->time_solve_ps += gettime() - solve_ps_start; - ms->count_solve_ps++; - } - goto finish; - } - mg_start = gettime(); LOGDEBUG( "mg solve cur %d top %d\n", reflevel, reflevel_top); @@ -741,8 +670,19 @@ void msa_mg_solve(CCTK_ARGUMENTS) /* outer-most level for this solve, use extrapolated lapse as the * boundary condition */ double time_interp_step = lapse_prev1_time[reflevel] - lapse_prev0_time[reflevel]; - double fact0 = (lapse_prev1_time[reflevel] - t) / time_interp_step; - double fact1 = (t - lapse_prev0_time[reflevel]) / time_interp_step; + double fact0, fact1; + + if (lapse_prev0_time[reflevel] == DBL_MAX && + lapse_prev1_time[reflevel] == DBL_MAX) { + fact0 = 0.0; + fact1 = 0.0; + } else if (lapse_prev0_time[reflevel] == DBL_MAX) { + fact0 = 0.0; + fact1 = 1.0; + } else { + fact0 = (lapse_prev1_time[reflevel] - t) / time_interp_step; + fact1 = (t - lapse_prev0_time[reflevel]) / time_interp_step; + } LOGDEBUG( "extrapolating BCs from t1=%g (%g) and t0=%g (%g)\n", lapse_prev1_time[reflevel], fact1, lapse_prev0_time[reflevel], fact0); -- cgit v1.2.3