chtMultiRegionFoam: Added support for residual convergence controls

Patch contributed by Tobias Holzmann
Resolves patch request https://bugs.openfoam.org/view.php?id=2524

applications/solvers/heatTransfer/chtMultiRegionFoam/checkResidualControls.H

@@ -0,0 +1,66 @@
+// Residual control used?
+bool resControlUsed = false;
+int nFluidControlled = fluidRegions.size();
+int nSolidControlled = solidRegions.size();
+
+// Check wheater there is a single regions that uses residual control
+forAll(fluidRegions, i)
+{
+    if (residualControlUsedFluid[i])
+    {
+        resControlUsed = true;
+    }
+}
+
+forAll(solidRegions, i)
+{
+    if(residualControlUsedSolid[i])
+    {
+        resControlUsed = true;
+    }
+}
+
+if (resControlUsed)
+{
+    int nFluidConv = 0;
+    int nSolidConv = 0;
+
+    // Sum of all converged regions (Note: if no residual control is used
+    // the residualReached* flag is already set to true)
+    forAll(fluidRegions, i)
+    {
+        if (residualReachedFluid[i])
+        {
+            nFluidConv++;
+        }
+    }
+
+    forAll(solidRegions, i)
+    {
+        if (residualReachedSolid[i])
+        {
+            nSolidConv++;
+        }
+    }
+
+    if (nFluidConv == nFluidControlled && nSolidConv == nSolidControlled)
+    {
+        // Activate flag to go to the 'Final' loop using the 'Final'
+        // relaxation factors
+        allRegionsConverged = true;
+    }
+}
+
+if (finalIter && resControlUsed && !allRegionsConverged)
+{
+    Info<< "\nRegions not converged after " << nOuterCorr
+        << " outer correctors" << endl;
+}
+else if (finalIter && resControlUsed && allRegionsConverged)
+{
+    Info<< "\nRegions converged after " << oCorr
+        << " outer correctors" << endl;
+
+    // Leave PIMPLE loop
+    break;
+}

applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionFoam.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -86,11 +86,18 @@ int main(int argc, char *argv[])
             }
         }
 
+        bool allRegionsConverged = false;
+        bool finalIter = false;
 
         // --- PIMPLE loop
         for (int oCorr=0; oCorr<nOuterCorr; oCorr++)
         {
-            bool finalIter = oCorr == nOuterCorr-1;
+            Info<< "Pimple iteration " << oCorr << "\n";
+
+            if (oCorr == nOuterCorr-1 || allRegionsConverged)
+            {
+                finalIter = true;
+            }
 
             forAll(fluidRegions, i)
             {
@@ -98,7 +105,9 @@ int main(int argc, char *argv[])
                     << fluidRegions[i].name() << endl;
                 #include "setRegionFluidFields.H"
                 #include "readFluidMultiRegionPIMPLEControls.H"
+                #include "readFluidMultiRegionResidualControls.H"
                 #include "solveFluid.H"
+                #include "residualControlsFluid.H"
             }
 
             forAll(solidRegions, i)
@@ -107,9 +116,12 @@ int main(int argc, char *argv[])
                     << solidRegions[i].name() << endl;
                 #include "setRegionSolidFields.H"
                 #include "readSolidMultiRegionPIMPLEControls.H"
+                #include "readSolidMultiRegionResidualControls.H"
                 #include "solveSolid.H"
+                #include "residualControlsSolid.H"
             }
 
+            #include "checkResidualControls.H"
         }
 
         runTime.write();

applications/solvers/heatTransfer/chtMultiRegionFoam/createFields.H

@@ -1,2 +1,7 @@
 #include "createFluidFields.H"
 #include "createSolidFields.H"
+
+SolverPerformance<vector> solvPerfU;
+SolverPerformance<scalar> solvPerfE;
+SolverPerformance<scalar> solvPerfh;
+SolverPerformance<scalar> solvPerfp_rgh;

applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/EEqn.H

@@ -26,7 +26,7 @@
 
     fvOptions.constrain(EEqn);
 
-    EEqn.solve(mesh.solver(he.select(finalIter)));
+    solvPerfE = EEqn.solve(mesh.solver(he.select(finalIter)));
 
     fvOptions.correct(he);
 

applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/UEqn.H

@@ -18,19 +18,20 @@
 
     if (momentumPredictor)
     {
-        solve
+        solvPerfU =
-        (
+            solve
-            UEqn
-          ==
-            fvc::reconstruct
             (
+                UEqn
+              ==
+                fvc::reconstruct
                 (
-                  - ghf*fvc::snGrad(rho)
+                    (
-                  - fvc::snGrad(p_rgh)
+                      - ghf*fvc::snGrad(rho)
-                )*mesh.magSf()
+                      - fvc::snGrad(p_rgh)
-            ),
+                    )*mesh.magSf()
-            mesh.solver(U.select(finalIter))
+                ),
-        );
+                mesh.solver(U.select(finalIter))
+            );
 
         fvOptions.correct(U);
         K = 0.5*magSqr(U);

applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/createFluidFields.H

@@ -15,6 +15,8 @@ PtrList<volScalarField> dpdtFluid(fluidRegions.size());
 
 List<scalar> initialMassFluid(fluidRegions.size());
 List<bool> frozenFlowFluid(fluidRegions.size(), false);
+List<bool> residualReachedFluid(fluidRegions.size(), true);
+List<bool> residualControlUsedFluid(fluidRegions.size(), false);
 
 PtrList<IOMRFZoneList> MRFfluid(fluidRegions.size());
 PtrList<fv::options> fluidFvOptions(fluidRegions.size());

applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/pEqn.H

@@ -45,16 +45,16 @@
               - fvm::laplacian(rhorAUf, p_rgh)
             );
 
-            p_rghEqn.solve
+            solvPerfp_rgh = p_rghEqn.solve
             (
                 mesh.solver
                 (
                     p_rgh.select
                     (
                         (
-                           oCorr == nOuterCorr-1
+                            oCorr == nOuterCorr-1
-                        && corr == nCorr-1
+                         && corr == nCorr-1
-                        && nonOrth == nNonOrthCorr
+                         && nonOrth == nNonOrthCorr
                         )
                     )
                 )

applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/readFluidMultiRegionResidualControls.H

@@ -0,0 +1,35 @@
+    const dictionary& residualControl =
+        mesh.solutionDict().subDict("PIMPLE").subOrEmptyDict("residualControl");
+
+    scalar UTol = -1.;
+    scalar ETol = -1.;
+    scalar p_rghTol = -1.;
+
+    if (!residualControl.empty())
+    {
+        if (!residualControl.subOrEmptyDict("U").empty())
+        {
+            UTol = readScalar(residualControl.subDict("U").lookup("tolerance"));
+        }
+
+        if (!residualControl.subOrEmptyDict("p_rgh").empty())
+        {
+            p_rghTol =
+                readScalar
+                (
+                    residualControl.subDict("p_rgh").lookup("tolerance")
+                );
+        }
+
+        if (!residualControl.subOrEmptyDict("h").empty())
+        {
+            ETol = readScalar(residualControl.subDict("h").lookup("tolerance"));
+        }
+
+        // Residual control used?
+        if (UTol != -1 || ETol != -1 || p_rghTol != -1)
+        {
+            residualControlUsed = true;
+            resReachedFluid = false;
+        }
+    }

applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/residualControlsFluid.H

@@ -0,0 +1,56 @@
+// Residual control used
+if (residualControlUsed)
+{
+    bool UConv = false;
+    bool p_rghConv = false;
+    bool EConv = false;
+
+    // Check which field is not used for control
+    {
+        if (UTol == -1)
+        {
+            UConv =  true;
+        }
+
+        if (p_rghTol == -1)
+        {
+            p_rghConv = true;
+        }
+
+        if (ETol == -1)
+        {
+            EConv = true;
+        }
+    }
+
+    // Get the last initial residual of the solvers
+    if (momentumPredictor && !UConv)
+    {
+        if (UTol > cmptMax(solvPerfU.initialResidual()))
+        {
+            UConv = true;
+        }
+    }
+
+    if (!p_rghConv)
+    {
+        if (p_rghTol > solvPerfp_rgh.initialResidual())
+        {
+            p_rghConv = true;
+        }
+    }
+
+    if (!EConv)
+    {
+        if (ETol > solvPerfE.initialResidual())
+        {
+            EConv = true;
+        }
+    }
+
+    // Check if each field is converged
+    if (UConv && p_rghConv && EConv)
+    {
+        resReachedFluid = true;
+    }
+}

applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/setRegionFluidFields.H

@@ -33,3 +33,6 @@
     );
 
     const bool frozenFlow = frozenFlowFluid[i];
+
+    bool& resReachedFluid = residualReachedFluid[i];
+    bool& residualControlUsed = residualControlUsedFluid[i];

applications/solvers/heatTransfer/chtMultiRegionFoam/solid/createSolidFields.H

@@ -6,6 +6,9 @@
     PtrList<volScalarField> betavSolid(solidRegions.size());
     PtrList<volSymmTensorField> aniAlphas(solidRegions.size());
 
+    List<bool> residualReachedSolid(solidRegions.size(), false);
+    List<bool> residualControlUsedSolid(solidRegions.size(), false);
+
     // Populate solid field pointer lists
     forAll(solidRegions, i)
     {

applications/solvers/heatTransfer/chtMultiRegionFoam/solid/readSolidMultiRegionResidualControls.H

@@ -0,0 +1,25 @@
+    const dictionary& residualControl =
+        mesh.solutionDict().subDict("PIMPLE").subOrEmptyDict("residualControl");
+
+    scalar hTol = -1.;
+
+    if (!residualControl.empty())
+    {
+        if (!residualControl.subOrEmptyDict("h").empty())
+        {
+            hTol = readScalar(residualControl.subDict("h").lookup("tolerance"));
+
+            // Used residual control for actual solid region
+            if (hTol != -1)
+            {
+                residualControlUsed = true;
+            }
+        }
+
+        // Residual control used?
+        if (hTol != -1)
+        {
+            residualControlUsed = true;
+            resReachedSolid = false;
+        }
+    }

applications/solvers/heatTransfer/chtMultiRegionFoam/solid/residualControlsSolid.H

@@ -0,0 +1,8 @@
+// Residual control used
+if (residualControlUsed)
+{
+    if (hTol > solvPerfh.initialResidual())
+    {
+        resReachedSolid = true;
+    }
+}

applications/solvers/heatTransfer/chtMultiRegionFoam/solid/setRegionSolidFields.H

@@ -30,3 +30,6 @@ volScalarField& h = thermo.he();
 const volScalarField& betav = betavSolid[i];
 
 fv::options& fvOptions = solidHeatSources[i];
+
+bool& resReachedSolid = residualReachedSolid[i];
+bool& residualControlUsed = residualControlUsedSolid[i];

applications/solvers/heatTransfer/chtMultiRegionFoam/solid/solidRegionDiffusionNo.H

@@ -2,7 +2,7 @@ scalar DiNum = -GREAT;
 
 forAll(solidRegions, i)
 {
-    //- Note: do not use setRegionSolidFields.H to avoid double registering Cp
+    // Note: do not use setRegionSolidFields.H to avoid double registering Cp
     //#include "setRegionSolidFields.H"
     const solidThermo& thermo = thermos[i];
 

applications/solvers/heatTransfer/chtMultiRegionFoam/solid/solveSolid.H

@@ -22,7 +22,7 @@ if (finalIter)
 
         fvOptions.constrain(hEqn);
 
-        hEqn.solve(mesh.solver(h.select(finalIter)));
+        solvPerfh = hEqn.solve(mesh.solver(h.select(finalIter)));
 
         fvOptions.correct(h);
     }