mg2d: account for effect of the ~u term on the relaxation factor

The maximum allowed time step in the presence of a -K * u term goes like
2 (dx ** 2) / (4 + K (dx ** 2))

3 files changed, 37 insertions, 0 deletions

diff --git a/ell_grid_solve.c b/ell_grid_solve.c
index 370deab..20de6bc 100644
--- a/ell_grid_solve.c
+++ b/ell_grid_solve.c
@@ -751,6 +751,9 @@ int mg2di_egs_init(EGSContext *ctx)
         else
             priv->r.relax_factor = r->relax_factor;
         priv->r.relax_factor *= ctx->step[0] * ctx->step[0];
+
+        if (r->relax_multiplier > 0.0)
+            priv->r.relax_factor *= r->relax_multiplier;
     }
 
     priv->fd_factors[MG2D_DIFF_COEFF_00] = 1.0 /  fd_denoms[ctx->fd_stencil - 1][MG2D_DIFF_COEFF_00];
diff --git a/ell_grid_solve.h b/ell_grid_solve.h
index a12e9c3..4ee6cb4 100644
--- a/ell_grid_solve.h
+++ b/ell_grid_solve.h
@@ -65,10 +65,24 @@ typedef struct EGSInternal EGSInternal;
 
 typedef struct EGSRelaxContext {
     /**
+     * The "time step" for relaxation is calculated as
+     * Δt = r_m * r_f * step[0] * step[1]
+     *
+     * Where
+     *  - r_m is relax_multiplier if specified, 1.0 otherwise
+     *  - r_f is relax_factor if specified, a default fd_stencil-dependent value
+     *    otherwise
+     */
+    /**
      * The time stepping factor in relaxation.
      */
     double relax_factor;
 
+    /**
+     * Multiplier for the time stepping factor.
+     */
+    double relax_multiplier;
+
     int64_t count_correct;
     int64_t time_correct;
 } EGSRelaxContext;
diff --git a/mg2d.c b/mg2d.c
index 003b88d..cbb2797 100644
--- a/mg2d.c
+++ b/mg2d.c
@@ -64,6 +64,8 @@ struct MG2DInternal {
 
     int cpuflags;
 
+    double fac00_max;
+
     int64_t time_solve;
     int64_t count_solve;
 };
@@ -407,11 +409,26 @@ static void restrict_diff_coeffs(MG2DContext *ctx, MG2DLevel *dst, MG2DLevel *sr
     }
 }
 
+static double findmax(double *arr, size_t size[2], ptrdiff_t stride)
+{
+    double ret = 0.0;
+    for (size_t y = 0; y < size[1]; y++)
+        for (size_t x = 0; x < size[0]; x++) {
+            double val = fabs(arr[y * stride + x]);
+            if (val > ret)
+                ret = val;
+        }
+    return ret;
+}
+
 static int mg_levels_init(MG2DContext *ctx)
 {
     MG2DInternal *priv = ctx->priv;
     MG2DLevel *cur, *prev;
 
+    priv->fac00_max = findmax(ctx->diff_coeffs[MG2D_DIFF_COEFF_00], priv->root->solver->domain_size,
+                              ctx->diff_coeffs_stride);
+
     cur  = priv->root;
     prev = NULL;
     while (cur) {
@@ -454,6 +471,9 @@ static int mg_levels_init(MG2DContext *ctx)
         if (cur->solver->solver_type == EGS_SOLVER_RELAXATION) {
             EGSRelaxContext *r = cur->solver->solver_data;
             r->relax_factor = ctx->cfl_factor;
+
+            r->relax_multiplier = 1.0 / (1.0 + cur->solver->step[0] * cur->solver->step[1] *
+                                               priv->fac00_max / 8.0);
         }
 
         prev = cur;