From 5e053d3453bfc4d1c75de237bb830725e87d3130 Mon Sep 17 00:00:00 2001
From: Anton Khirnov <anton@khirnov.net>
Date: Tue, 4 Jun 2019 16:26:53 +0200
Subject: mg2d: use the all-comp global residual norm for divergence testing

Residual may grow dramatically in certain components while globally
decreasing.
---
 mg2d.c | 33 ++++++++++++++++-----------------
 1 file changed, 16 insertions(+), 17 deletions(-)

diff --git a/mg2d.c b/mg2d.c
index 578dfe0..9313ce3 100644
--- a/mg2d.c
+++ b/mg2d.c
@@ -220,6 +220,18 @@ static int mg_prolong_u(MG2DContext *ctx, MG2DLevel *l)
     return 0;
 }
 
+static double mg_resmax_get(MG2DLevel *l)
+{
+    double ret = l->solver->residual_max;
+
+    if (l->dg->nb_components > 1) {
+        MPI_Allreduce(MPI_IN_PLACE, &ret, 1,
+                      MPI_DOUBLE, MPI_MAX, l->mpi_comm);
+    }
+
+    return ret;
+}
+
 static int mg_solve_subgrid(MG2DContext *ctx, MG2DLevel *level, int allow_exact)
 {
     enum EGSType solve_type = (allow_exact && level->dg->nb_components == 1 &&
@@ -244,7 +256,7 @@ static int mg_solve_subgrid(MG2DContext *ctx, MG2DLevel *level, int allow_exact)
         level->egs_init_flags |= EGS_INIT_FLAG_SAME_DIFF_COEFFS;
     }
 
-    res_old = level->solver->residual_max;
+    res_old = mg_resmax_get(level);
 
     /* handle exact solve */
     if (solve_type == EGS_SOLVE_EXACT) {
@@ -315,7 +327,7 @@ static int mg_solve_subgrid(MG2DContext *ctx, MG2DLevel *level, int allow_exact)
     }
 
 finish:
-    res_new = level->solver->residual_max;
+    res_new = mg_resmax_get(level);
     if (!isfinite(res_new) ||
         (res_new > 1e2 * DBL_EPSILON && res_old / res_new <= 1e-1)) {
         mg2di_log(&ctx->priv->logger, MG2D_LOG_ERROR,
@@ -843,19 +855,6 @@ static int threadpool_init(MG2DContext *ctx)
 
 static int mg_levels_alloc(MG2DContext *ctx);
 
-static double mg_resmax_get(MG2DContext *ctx)
-{
-    MG2DInternal *priv = ctx->priv;
-    double ret = priv->root->solver->residual_max;
-
-    if (priv->root->dg->nb_components > 1) {
-        MPI_Allreduce(MPI_IN_PLACE, &ret, 1,
-                      MPI_DOUBLE, MPI_MAX, priv->mpi_comm);
-    }
-
-    return ret;
-}
-
 int mg2d_solve(MG2DContext *ctx)
 {
     MG2DInternal *priv = ctx->priv;
@@ -895,7 +894,7 @@ int mg2d_solve(MG2DContext *ctx)
 
     root->egs_init_flags |= EGS_INIT_FLAG_SAME_DIFF_COEFFS;
 
-    res_orig = mg_resmax_get(ctx);
+    res_orig = mg_resmax_get(root);
     res_prev = res_orig;
 
     for (int i = 0; i < ctx->maxiter; i++) {
@@ -903,7 +902,7 @@ int mg2d_solve(MG2DContext *ctx)
         if (ret < 0)
             goto fail;
 
-        res_cur = mg_resmax_get(ctx);
+        res_cur = mg_resmax_get(root);
 
         if (res_cur < ctx->tol) {
             mg2di_log(&priv->logger, MG2D_LOG_INFO, "converged on iteration %d, residual %g\n",
-- 
cgit v1.2.3