Make the switchoff time configurable.

2 files changed, 21 insertions, 10 deletions

diff --git a/param.ccl b/param.ccl
index 4e2053a..63dc31c 100644
--- a/param.ccl
+++ b/param.ccl
@@ -26,6 +26,11 @@ CCTK_REAL scale_power "" STEERABLE=recover
     0: :: ""
 } 64.0
 
+CCTK_REAL switchoff_time "" STEERABLE=recover
+{
+   0: :: ""
+} 9999.0
+
 BOOLEAN export_coeffs "Export the coefficients of the spectral expansion in alpha_coeffs" STEERABLE=recover
 {
 } "no"
diff --git a/src/qms.c b/src/qms.c
index fd40e67..a40df08 100644
--- a/src/qms.c
+++ b/src/qms.c
@@ -264,6 +264,9 @@ void quasimaximal_slicing_axi_solve(CCTK_ARGUMENTS)
 
     time = cctkGH->cctk_time / ms->gh->cctk_delta_time;
 
+    if (cctkGH->cctk_time >= switchoff_time + 1.0)
+        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))
         return;
@@ -298,6 +301,8 @@ void quasimaximal_slicing_axi_eval(CCTK_ARGUMENTS)
     DECLARE_CCTK_ARGUMENTS;
     DECLARE_CCTK_PARAMETERS;
 
+    double time;
+
     int64_t expand_start;
 
     double *coeffs = NULL;
@@ -306,6 +311,11 @@ void quasimaximal_slicing_axi_eval(CCTK_ARGUMENTS)
     if (!qms_context)
         context_init(cctkGH);
 
+    time  = cctkGH->cctk_time;
+
+    if (time >= switchoff_time + 1.0)
+        return;
+
     ms = qms_context;
 
     cp = get_coord_patch(ms, x, y, z, scale_factor, scale_power);
@@ -324,22 +334,18 @@ void quasimaximal_slicing_axi_eval(CCTK_ARGUMENTS)
         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 time  = ms->gh->cctk_time;
 
-        //double fact;
+        double fact;
 
-        //if (time > 2.0)
-        //    fact = 1.0;
-        //else if (time < 0.1)
-        //    fact = 0.0;
-        //else
-        //    fact = (1.0 - exp(-pow((time - 0.0) / 0.25, 4.0)));
-        //fact = 1.0;
+        if (time < switchoff_time)
+            fact = 1.0;
+        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);
+            coeffs[i] = (coeffs1[i] * (time - time0) / (time1 - time0) + coeffs0[i] * (time - time1) / (time0 - time1)) * fact;
 
     }
 #endif