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))

1 files changed, 20 insertions, 0 deletions

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;