From 80dff3e38de86f886952ebf133e03551cdbef541 Mon Sep 17 00:00:00 2001
From: Anton Khirnov <anton@khirnov.net>
Date: Thu, 12 Sep 2019 12:23:18 +0200
Subject: Allow specifying maximum level on which to solve for W.

---
 param.ccl |  5 +++++
 src/qms.c | 19 ++++++++++++++-----
 2 files changed, 19 insertions(+), 5 deletions(-)

diff --git a/param.ccl b/param.ccl
index b7d9bcf..22c4d21 100644
--- a/param.ccl
+++ b/param.ccl
@@ -12,6 +12,11 @@ CCTK_INT solve_level "" STEERABLE=recover
     0: :: ""
 } 0
 
+CCTK_INT solve_level_max "" STEERABLE=recover
+{
+    0: :: ""
+} 0
+
 CCTK_REAL switchoff_time "" STEERABLE=recover
 {
    : :: ""
diff --git a/src/qms.c b/src/qms.c
index c22cc97..7f0f3e1 100644
--- a/src/qms.c
+++ b/src/qms.c
@@ -82,6 +82,7 @@ typedef struct QMSMGContext {
     double tol_residual_base;
     double cfl_factor;
     int solve_level;
+    int solve_level_max;
     int boundary_offset;
 
     CoordPatch **patches;
@@ -1292,6 +1293,9 @@ void qms_mg_solve(CCTK_ARGUMENTS)
 
     grid_size = cctkGH->cctk_lsh[0] * cctkGH->cctk_lsh[1] * cctkGH->cctk_lsh[2];
 
+    if (reflevel < ms->solve_level_max)
+        goto skip_solve;
+
     start = gettime();
 
     fill_eq_coeffs(ms, cp, cp->solver);
@@ -1300,7 +1304,7 @@ void qms_mg_solve(CCTK_ARGUMENTS)
     ms->count_fill++;
 
     start = gettime();
-    if (reflevel > 0) {
+    if (reflevel > ms->solve_level_max) {
         /* outer-most level for this solve, use extrapolated values as boundary condition */
         if (fabs(time - W_val1_time[reflevel - 1]) > 1e-13) {
             const double fact0 = (W_val1_time[reflevel] - time) / dt;
@@ -1353,8 +1357,10 @@ void qms_mg_solve(CCTK_ARGUMENTS)
     if (ret < 0)
         CCTK_WARN(0, "Error solving the quasi-maximal slicing equation");
 
+skip_solve:
+
     start = gettime();
-    {
+    if (reflevel >= ms->solve_level_max) {
         double *W_val_tl1 = CCTK_VarDataPtr(cctkGH, 1, "QuasiMaximalSlicingMG::W_val");
         solution_to_grid(cp, cp->solver, W_val);
         memcpy(W_val_tl1, W_val, grid_size * sizeof(*W_val_tl1));
@@ -1528,7 +1534,7 @@ void qms_mg_eval(CCTK_ARGUMENTS)
     ms->count_eval++;
 }
 
-static int context_init(cGH *gh, int solve_level, int fd_stencil, int maxiter, int exact_size, int nb_cycles,
+static int context_init(cGH *gh, int solve_level, int solve_level_max, int fd_stencil, int maxiter, int exact_size, int nb_cycles,
                         int nb_relax_pre, int nb_relax_post, double tol_residual_base,
                         double cfl_factor, const char *loglevel_str, int boundary_offset,
                         QMSMGContext **ctx)
@@ -1544,6 +1550,7 @@ static int context_init(cGH *gh, int solve_level, int fd_stencil, int maxiter, i
     qms->gh = gh;
 
     qms->solve_level = solve_level;
+    qms->solve_level_max = solve_level_max;
     qms->fd_stencil    = fd_stencil;
     qms->maxiter       = maxiter;
     qms->max_exact_size = exact_size;
@@ -1591,16 +1598,18 @@ void qms_mg_init(CCTK_ARGUMENTS)
     DECLARE_CCTK_ARGUMENTS;
     DECLARE_CCTK_PARAMETERS;
 
+    const int reflevel = ctz(cctk_levfac[0]);
+
     int use_tapered_grids = *(int*)CCTK_ParameterGet("use_tapered_grids", "Carpet", &use_tapered_grids);
     int ret;
 
-    fprintf(stderr, "%d init\n", ctz(cctk_levfac[0]));
+    fprintf(stderr, "%d init\n", reflevel);
 
     if (!use_tapered_grids)
         CCTK_WARN(0, "MaximalSlicingAxiMG only works with use_tapered_grids=1");
 
     if (!qms_context) {
-        ret = context_init(cctkGH, solve_level, fd_stencil, maxiter, exact_size, nb_cycles,
+        ret = context_init(cctkGH, solve_level, solve_level_max, fd_stencil, maxiter, exact_size, nb_cycles,
                            nb_relax_pre, nb_relax_post, tol_residual_base,
                            cfl_factor, loglevel, boundary_offset, &qms_context);
         if (ret < 0)
-- 
cgit v1.2.3