chtMultiRegionSimpleFoam: Added residual convergence control

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

applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/checkResidualControls.H

@@ -0,0 +1,68 @@
+// 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;
+        break;
+    }
+}
+
+forAll(solidRegions, i)
+{
+    if(residualControlUsedSolid[i])
+    {
+        resControlUsed = true;
+        break;
+    }
+}
+
+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 (!runTime.run())
+{
+    Info<< "\nRegions not converged after " << runTime.timeName()
+        << " iterations" << endl;
+}
+else if (runTime.run() && resControlUsed && allRegionsConverged)
+{
+    Info<< "\nRegions converged after " << runTime.timeName()
+        << " iterations" << endl;
+
+    // Leave SIMPLE loop
+    break;
+}

applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/chtMultiRegionSimpleFoam.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
@@ -54,9 +54,13 @@ int main(int argc, char *argv[])
     #include "createFields.H"
     #include "initContinuityErrs.H"
 
+    // Residual control
+    bool allRegionsConverged = false;
 
-    while (runTime.loop())
+    while (runTime.run())
     {
+        runTime++;
+
         Info<< "Time = " << runTime.timeName() << nl << endl;
 
         forAll(fluidRegions, i)
@@ -65,7 +69,9 @@ int main(int argc, char *argv[])
                 << fluidRegions[i].name() << endl;
             #include "setRegionFluidFields.H"
             #include "readFluidMultiRegionSIMPLEControls.H"
+            #include "readFluidMultiRegionResidualControls.H"
             #include "solveFluid.H"
+            #include "residualControlsFluid.H"
         }
 
         forAll(solidRegions, i)
@@ -74,9 +80,13 @@ int main(int argc, char *argv[])
                 << solidRegions[i].name() << endl;
             #include "setRegionSolidFields.H"
             #include "readSolidMultiRegionSIMPLEControls.H"
+            #include "readSolidMultiRegionResidualControls.H"
             #include "solveSolid.H"
+            #include "residualControlsSolid.H"
         }
 
+        #include "checkResidualControls.H"
+
         runTime.write();
 
         Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"

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

@@ -20,7 +20,7 @@
 
     fvOptions.constrain(EEqn);
 
-    EEqn.solve();
+    solvPerfE = EEqn.solve();
 
     fvOptions.correct(he);
 

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

@@ -16,17 +16,20 @@
 
     fvOptions.constrain(UEqn);
 
-    solve
+    if (momentumPredictor)
-    (
+    {
-        UEqn
+        solvPerfU = solve
-     ==
-        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()
-    );
+            )
+        );
+    }
 
     fvOptions.correct(U);

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

@@ -15,6 +15,9 @@ List<scalar> initialMassFluid(fluidRegions.size());
 List<label> pRefCellFluid(fluidRegions.size(), 0);
 List<scalar> pRefValueFluid(fluidRegions.size(), 0.0);
 List<bool> frozenFlowFluid(fluidRegions.size(), false);
+List<bool> residualReachedFluid(fluidRegions.size(), true);
+List<bool> residualControlUsedFluid(fluidRegions.size(), false);
+List<bool> momentumPredictorFluid(fluidRegions.size(), false);
 
 PtrList<dimensionedScalar> rhoMax(fluidRegions.size());
 PtrList<dimensionedScalar> rhoMin(fluidRegions.size());

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

@@ -43,7 +43,7 @@
             compressible ? getRefCellValue(p_rgh, pRefCell) : pRefValue
         );
 
-        p_rghEqn.solve();
+        solvPerfp_rgh = p_rghEqn.solve();
 
         if (nonOrth == nNonOrthCorr)
         {

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

@@ -0,0 +1,35 @@
+    const dictionary& residualControl =
+        mesh.solutionDict().subDict("SIMPLE").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/chtMultiRegionSimpleFoam/fluid/residualControlsFluid.H

@@ -0,0 +1 @@
+../../fluid/residualControlsFluid.H

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

@@ -33,3 +33,8 @@
     const dimensionedVector& g = gFluid[i];
     const volScalarField& gh = ghFluid[i];
     const surfaceScalarField& ghf = ghfFluid[i];
+
+    const bool momentumPredictor = momentumPredictorFluid[i];
+
+    bool& resReachedFluid = residualReachedFluid[i];
+    bool& residualControlUsed = residualControlUsedFluid[i];

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

@@ -0,0 +1,19 @@
+    const dictionary& residualControl =
+        mesh.solutionDict().subDict("SIMPLE").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;
+                resReachedSolid = false;
+            }
+        }
+    }

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

@@ -16,7 +16,7 @@
 
         fvOptions.constrain(hEqn);
 
-        hEqn.solve();
+        solvPerfh = hEqn.solve();
 
         fvOptions.correct(h);
     }

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

@@ -7,7 +7,7 @@ if (residualControlUsed)
 
     // Check which field is not used for control
     {
-        if (UTol == -1)
+        if (UTol == -1 || !momentumPredictor)
         {
             UConv =  true;
         }