Added rescue procedure on block solver divergence. Experimental

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
+        // Check for divergence
-        UpEqn.retrieveSolution(0, U.internalField());
+        if (mag(residual.finalResidual() - residual.initialResidual()) > 1)
-        UpEqn.retrieveSolution(3, p.internalField());
+        {
+            Info<< "DIVERGENCE: rescue step" << endl;
 
-        U.correctBoundaryConditions();
+            // Solving potential flow equation and correcting velocities
-        p.correctBoundaryConditions();
+            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);
 

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;