From 6c01d9a6a8d80bebb442f5cf77eb12d5b07039ac Mon Sep 17 00:00:00 2001
From: Hrvoje Jasak <h.jasak@wikki.co.uk>
Date: Mon, 6 Feb 2017 11:47:00 +0000
Subject: [PATCH] Added rescue procedure on block solver divergence. 
 Experimental

---
 .../coupled/MRFPorousFoam/MRFPorousFoam.C     | 37 +++++++++++++++----
 .../MRFPorousFoam/initConvergenceCheck.H      |  3 +-
 2 files changed, 32 insertions(+), 8 deletions(-)

diff --git a/applications/solvers/coupled/MRFPorousFoam/MRFPorousFoam.C b/applications/solvers/coupled/MRFPorousFoam/MRFPorousFoam.C
index 4545f7093..8e96be5a4 100644
--- a/applications/solvers/coupled/MRFPorousFoam/MRFPorousFoam.C
+++ b/applications/solvers/coupled/MRFPorousFoam/MRFPorousFoam.C
@@ -81,17 +81,40 @@ int main(int argc, char *argv[])
 #       include "couplingTerms.H"
 
         // Solve the block matrix
-        maxResidual = cmptMax(UpEqn.solve().initialResidual());
+        residual = UpEqn.solve();
+        maxResidual = cmptMax(residual.initialResidual());
 
-        // Retrieve solution
-        UpEqn.retrieveSolution(0, U.internalField());
-        UpEqn.retrieveSolution(3, p.internalField());
+        // Check for divergence
+        if (mag(residual.finalResidual() - residual.initialResidual()) > 1)
+        {
+            Info<< "DIVERGENCE: rescue step" << endl;
 
-        U.correctBoundaryConditions();
-        p.correctBoundaryConditions();
+            // Solving potential flow equation and correcting velocities
+            phi = (fvc::interpolate(U) & mesh.Sf());
+            
+            solve
+            (
+                tpEqn
+             ==
+              - fvc::div(U)
+            );
 
-        phi = (fvc::interpolate(U) & mesh.Sf()) + tpEqn().flux() + tpresSource;
+            phi += tpEqn().flux() + tpresSource;
 
+            U = fvc::reconstruct(phi);
+        }
+        else
+        {
+            // Retrieve solution
+            UpEqn.retrieveSolution(0, U.internalField());
+            UpEqn.retrieveSolution(3, p.internalField());
+
+            U.correctBoundaryConditions();
+            p.correctBoundaryConditions();
+
+            phi = (fvc::interpolate(U) & mesh.Sf()) + tpEqn().flux() + tpresSource;
+        }
+        
         // Make flux relative in rotating zones
         mrfZones.relativeFlux(phi);
 
diff --git a/applications/solvers/coupled/MRFPorousFoam/initConvergenceCheck.H b/applications/solvers/coupled/MRFPorousFoam/initConvergenceCheck.H
index 482bc421d..6b4622db6 100644
--- a/applications/solvers/coupled/MRFPorousFoam/initConvergenceCheck.H
+++ b/applications/solvers/coupled/MRFPorousFoam/initConvergenceCheck.H
@@ -1,4 +1,5 @@
-// initialize values for convergence checks
+    // initialize values for convergence checks
+    BlockSolverPerformance<vector4> residual;
 
     scalar maxResidual = 0;
     scalar convergenceCriterion = 0;