mg2d: take into account the factor in front of ∂_xx when computing the CFL factor

1 files changed, 9 insertions, 5 deletions

diff --git a/mg2d.c b/mg2d.c
index 0a17793..5b94022 100644
--- a/mg2d.c
+++ b/mg2d.c
@@ -64,8 +64,6 @@ struct MG2DInternal {
 
     int cpuflags;
 
-    double fac00_max;
-
     int64_t time_solve;
     int64_t count_solve;
 };
@@ -425,9 +423,15 @@ static int mg_levels_init(MG2DContext *ctx)
 {
     MG2DInternal *priv = ctx->priv;
     MG2DLevel *cur, *prev;
+    double diff0_max, diff2_max, tmp;
 
-    priv->fac00_max = findmax(ctx->diff_coeffs[MG2D_DIFF_COEFF_00], priv->root->solver->domain_size,
+    diff0_max = findmax(ctx->diff_coeffs[MG2D_DIFF_COEFF_00], priv->root->solver->domain_size,
+                        ctx->diff_coeffs_stride);
+    diff2_max = findmax(ctx->diff_coeffs[MG2D_DIFF_COEFF_20], priv->root->solver->domain_size,
+                        ctx->diff_coeffs_stride);
+    tmp = findmax(ctx->diff_coeffs[MG2D_DIFF_COEFF_02], priv->root->solver->domain_size,
                               ctx->diff_coeffs_stride);
+    diff2_max = MAX(diff2_max, tmp);
 
     cur  = priv->root;
     prev = NULL;
@@ -473,8 +477,8 @@ static int mg_levels_init(MG2DContext *ctx)
             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);
+            r->relax_multiplier = 1.0 / (diff2_max + cur->solver->step[0] * cur->solver->step[1] *
+                                                     diff0_max / 8.0);
         }
 
         prev = cur;