mg2d: failer earlier when relaxation does not converge

3 files changed, 18 insertions, 1 deletions

diff --git a/libmg2d.v b/libmg2d.v
index 96151a5..f30680f 100644
--- a/libmg2d.v
+++ b/libmg2d.v
@@ -1,4 +1,4 @@
-LIBMG2D_5 {
+LIBMG2D_6 {
     global: mg2d_*;
     local: *;
 };
diff --git a/mg2d.c b/mg2d.c
index 061ae91..1f5d86e 100644
--- a/mg2d.c
+++ b/mg2d.c
@@ -17,6 +17,7 @@
  */
 
 #include <errno.h>
+#include <float.h>
 #include <math.h>
 #include <stdlib.h>
 #include <stdint.h>
@@ -256,6 +257,7 @@ static int mg_relax_step(MG2DLevel *level)
 
 static int mg_solve_subgrid(MG2DContext *ctx, MG2DLevel *level)
 {
+    double res_old, res_new;
     int ret;
 
     /* on the refined levels, use zero as the initial guess for the
@@ -270,6 +272,8 @@ static int mg_solve_subgrid(MG2DContext *ctx, MG2DLevel *level)
             return ret;
     }
 
+    res_old = level->solver->residual_max;
+
     /* handle coarsest grid */
     if (!level->child) {
 
@@ -333,6 +337,15 @@ static int mg_solve_subgrid(MG2DContext *ctx, MG2DLevel *level)
     }
 
 finish:
+    res_new = level->solver->residual_max;
+    if (!isfinite(res_new) ||
+        (res_new > DBL_EPSILON && res_old / res_new <= 1e-1)) {
+        mg2di_log(&ctx->priv->logger, MG2D_LOG_ERROR,
+                  "The relaxation step at level %d has diverged: %g -> %g\n",
+                  level->depth, res_old, res_new);
+        return MG2D_ERR_DIVERGE;
+    }
+
     return 0;
 }
 
diff --git a/mg2d.h b/mg2d.h
index a587a18..58308fd 100644
--- a/mg2d.h
+++ b/mg2d.h
@@ -22,6 +22,10 @@
 #include <stddef.h>
 #include <stdint.h>
 
+enum MG2DError {
+    MG2D_ERR_DIVERGE = -0xff00,
+};
+
 /**
  * Type of the boundary condition on a given boundary.
  */