Drop the use of the pseudospectral solver.

It does not produce good enough results.

1 files 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);