Add an option to scale residual tolerance with dx**-2

2 files changed, 14 insertions, 3 deletions

diff --git a/param.ccl b/param.ccl
index 80f2c55..55b52a6 100644
--- a/param.ccl
+++ b/param.ccl
@@ -50,6 +50,12 @@ CCTK_REAL tol_residual "maximum absolute value of the residual" STEERABLE=always
 } 1e-12
 
 RESTRICTED:
+CCTK_REAL tol_residual_base "maximum absolute value of the residual for dx=1.0, scaled by 1/(dx ** 2)" STEERABLE=always
+{
+    0: :: ""
+} 0.0
+
+RESTRICTED:
 CCTK_REAL cfl_factor "" STEERABLE=always
 {
     0: :: ""
diff --git a/src/maximal_slicing_axi_mg.c b/src/maximal_slicing_axi_mg.c
index e7cbd52..bef2139 100644
--- a/src/maximal_slicing_axi_mg.c
+++ b/src/maximal_slicing_axi_mg.c
@@ -72,6 +72,7 @@ typedef struct MSMGContext {
     int nb_relax_post;
     int nb_relax_pre;
     double tol_residual;
+    double tol_residual_base;
     double cfl_factor;
 
     CoordPatch *patches;
@@ -219,7 +220,10 @@ static CoordPatch *get_coord_patch(MSMGContext *ms, int level)
     cp->solver->boundaries[MG2D_BOUNDARY_1U]->type = level ? MG2D_BC_TYPE_FIXVAL : MG2D_BC_TYPE_FALLOFF;
 
     cp->solver->maxiter   = ms->maxiter;
-    cp->solver->tol       = ms->tol_residual;
+    if (ms->tol_residual_base > 0.0)
+        cp->solver->tol = ms->tol_residual_base / SQR(cp->solver->step[0]);
+    else
+        cp->solver->tol       = ms->tol_residual;
     cp->solver->nb_cycles = ms->nb_cycles;
     cp->solver->nb_relax_post = ms->nb_relax_post;
     cp->solver->nb_relax_pre  = ms->nb_relax_pre;
@@ -294,7 +298,7 @@ static void print_stats(MSMGContext *ms)
 }
 
 static int context_init(cGH *gh, int fd_stencil, int maxiter, int exact_size, int nb_cycles,
-                        int nb_relax_pre, int nb_relax_post, double tol_residual,
+                        int nb_relax_pre, int nb_relax_post, double tol_residual, double tol_residual_base,
                         double cfl_factor,
                         const char *loglevel_str,
                         MSMGContext **ctx)
@@ -314,6 +318,7 @@ static int context_init(cGH *gh, int fd_stencil, int maxiter, int exact_size, in
     ms->nb_relax_pre  = nb_relax_pre;
     ms->nb_relax_post = nb_relax_post;
     ms->tol_residual  = tol_residual;
+    ms->tol_residual_base = tol_residual_base;
     ms->cfl_factor    = cfl_factor;
 
     for (int i = 0; i < ARRAY_ELEMS(log_levels); i++) {
@@ -839,7 +844,7 @@ void msa_mg_init(CCTK_ARGUMENTS)
     int ret;
 
     ret = context_init(cctkGH, fd_stencil, maxiter, exact_size, nb_cycles,
-                       nb_relax_pre, nb_relax_post, tol_residual,
+                       nb_relax_pre, nb_relax_post, tol_residual, tol_residual_base,
                        cfl_factor,
                        loglevel, &ms);
     if (ret < 0)