INT: solver function objects: add outer-loop convergence checks

Co-authored-by: Kutalmis Bercin <kutalmis.bercin@esi-group.com>

src/functionObjects/solvers/electricPotential/electricPotential.C

@@ -208,6 +208,7 @@ Foam::functionObjects::electricPotential::electricPotential
         )
     ),
     fvOptions_(mesh_),
+    tol_(1),
     nCorr_(1),
     writeDerivedFields_(false),
     electricField_(false)
@@ -258,6 +259,7 @@ bool Foam::functionObjects::electricPotential::read(const dictionary& dict)
     dict.readIfPresent("sigma", sigma_);
     dict.readIfPresent("epsilonr", epsilonr_);
     dict.readIfPresent("nCorr", nCorr_);
+    dict.readIfPresent("tolerance", tol_);
     dict.readIfPresent("writeDerivedFields", writeDerivedFields_);
     dict.readIfPresent("electricField", electricField_);
 
@@ -346,6 +348,10 @@ bool Foam::functionObjects::electricPotential::execute()
 {
     Log << type() << " execute: " << name() << endl;
 
+    // Convergence monitor parameters
+    bool converged = false;
+    label iter = 0;
+
     tmp<volScalarField> tsigma = this->sigma();
     const auto& sigma = tsigma();
 
@@ -362,7 +368,9 @@ bool Foam::functionObjects::electricPotential::execute()
 
         fvOptions_.constrain(eVEqn);
 
-        eVEqn.solve();
+        ++iter;
+        converged = (eVEqn.solve().initialResidual() < tol_);
+        if (converged) break;
     }
 
     if (electricField_)
@@ -371,6 +379,14 @@ bool Foam::functionObjects::electricPotential::execute()
         E == -fvc::grad(eV);
     }
 
+    if (converged)
+    {
+        Log << type() << ": " << name() << ": "
+            << eV.name() << " is converged." << nl
+            << tab << "initial-residual tolerance: " << tol_ << nl
+            << tab << "outer iteration: " << iter << nl;
+    }
+
     Log << endl;
 
     return true;

src/functionObjects/solvers/electricPotential/electricPotential.H

@@ -124,6 +124,7 @@ Usage
         electricField         <bool>;
         E                     <word>;
         fvOptions             <dict>;
+        tolerance             <scalar>;
 
         // Inherited entries
         ...
@@ -143,6 +144,7 @@ Usage
       electricField | Flag to calculate electric field  | bool | no | false
       E    | Name of electric field           | word | no | electricPotential:E
       fvOptions | List of finite-volume options       | dict | no   | -
+      tolerance | Outer-loop initial-residual tolerance | scalar | no | 1
     \endtable
 
     The inherited entries are elaborated in:
@@ -218,6 +220,9 @@ class electricPotential
         //- Run-time selectable finite volume options
         fv::optionList fvOptions_;
 
+        //- Outer-loop initial-residual tolerance
+        scalar tol_;
+
         //- Number of corrector iterations
         int nCorr_;
 

src/functionObjects/solvers/energyTransport/energyTransport.C

@@ -192,6 +192,7 @@ Foam::functionObjects::energyTransport::energyTransport
     phiName_(dict.getOrDefault<word>("phi", "phi")),
     rhoName_(dict.getOrDefault<word>("rho", "rho")),
     nCorr_(0),
+    tol_(1),
     schemesField_("unknown-schemesField"),
     fvOptions_(mesh_),
     multiphaseThermo_(dict.subOrEmptyDict("phaseThermos")),
@@ -323,6 +324,7 @@ bool Foam::functionObjects::energyTransport::read(const dictionary& dict)
     schemesField_ = dict.getOrDefault("schemesField", fieldName_);
 
     dict.readIfPresent("nCorr", nCorr_);
+    dict.readIfPresent("tolerance", tol_);
 
     if (dict.found("fvOptions"))
     {
@@ -355,6 +357,10 @@ bool Foam::functionObjects::energyTransport::execute()
     scalar relaxCoeff = 0;
     mesh_.relaxEquation(schemesField_, relaxCoeff);
 
+    // Convergence monitor parameters
+    bool converged = false;
+    label iter = 0;
+
     if (phi.dimensions() == dimMass/dimTime)
     {
         rhoCp_ = rho()*Cp();
@@ -376,7 +382,9 @@ bool Foam::functionObjects::energyTransport::execute()
 
             fvOptions_.constrain(sEqn);
 
-            sEqn.solve(schemesField_);
+            ++iter;
+            converged = (sEqn.solve(schemesField_).initialResidual() < tol_);
+            if (converged) break;
         }
     }
     else if (phi.dimensions() == dimVolume/dimTime)
@@ -408,7 +416,9 @@ bool Foam::functionObjects::energyTransport::execute()
 
             fvOptions_.constrain(sEqn);
 
-            sEqn.solve(schemesField_);
+            ++iter;
+            converged = (sEqn.solve(schemesField_).initialResidual() < tol_);
+            if (converged) break;
         }
     }
     else
@@ -419,6 +429,14 @@ bool Foam::functionObjects::energyTransport::execute()
             << dimVolume/dimTime << exit(FatalError);
     }
 
+    if (converged)
+    {
+        Log << type() << ": " << name() << ": "
+            << s.name() << " is converged." << nl
+            << tab << "initial-residual tolerance: " << tol_ << nl
+            << tab << "outer iteration: " << iter << nl;
+    }
+
     Log << endl;
 
     return true;

src/functionObjects/solvers/energyTransport/energyTransport.H

@@ -163,6 +163,7 @@ Usage
         phi          | Name of flux field      | no          | phi
         rho          | Name of density field   | no          | rho
         nCorr        | Number of correctors    | no          | 0
+        tolerance    | Outer-loop initial-residual tolerance | no | 1
         schemesField | Name of field to specify schemes | no | field name
         fvOptions    | List of scalar sources  | no          |
         Cp           | Heat capacity for single phase | no   | 0
@@ -217,6 +218,9 @@ class energyTransport
         //- Number of corrector iterations (optional)
         label nCorr_;
 
+        //- Outer-loop initial-residual tolerance
+        scalar tol_;
+
         //- Name of field whose schemes are used (optional)
         word schemesField_;
 

src/functionObjects/solvers/scalarTransport/scalarTransport.C

@@ -183,6 +183,7 @@ Foam::functionObjects::scalarTransport::scalarTransport
     ),
     D_(0),
     constantD_(false),
+    tol_(1),
     nCorr_(0),
     resetOnStartUp_(false),
     schemesField_("unknown-schemesField"),
@@ -225,6 +226,7 @@ bool Foam::functionObjects::scalarTransport::read(const dictionary& dict)
     alphaD_ = dict.getOrDefault<scalar>("alphaD", 1);
     alphaDt_ = dict.getOrDefault<scalar>("alphaDt", 1);
 
+    dict.readIfPresent("tolerance", tol_);
     dict.readIfPresent("nCorr", nCorr_);
     dict.readIfPresent("resetOnStartUp", resetOnStartUp_);
 
@@ -256,6 +258,10 @@ bool Foam::functionObjects::scalarTransport::execute()
     scalar relaxCoeff = 0;
     mesh_.relaxEquation(schemesField_, relaxCoeff);
 
+    // Convergence monitor parameters
+    bool converged = false;
+    label iter = 0;
+
     // Two phase scalar transport
     if (phaseName_ != "none")
     {
@@ -285,9 +291,13 @@ bool Foam::functionObjects::scalarTransport::execute()
 
             sEqn.relax(relaxCoeff);
             fvOptions_.constrain(sEqn);
-            sEqn.solve(schemesField_);
+
+            ++iter;
+            converged = (sEqn.solve(schemesField_).initialResidual() < tol_);
 
             tTPhiUD = sEqn.flux();
+
+            if (converged) break;
         }
 
         if (bounded01_)
@@ -323,7 +333,9 @@ bool Foam::functionObjects::scalarTransport::execute()
 
             fvOptions_.constrain(sEqn);
 
-            sEqn.solve(schemesField_);
+            ++iter;
+            converged = (sEqn.solve(schemesField_).initialResidual() < tol_);
+            if (converged) break;
         }
     }
     else if (phi.dimensions() == dimVolume/dimTime)
@@ -343,7 +355,9 @@ bool Foam::functionObjects::scalarTransport::execute()
 
             fvOptions_.constrain(sEqn);
 
-            sEqn.solve(schemesField_);
+            ++iter;
+            converged = (sEqn.solve(schemesField_).initialResidual() < tol_);
+            if (converged) break;
         }
     }
     else
@@ -354,6 +368,14 @@ bool Foam::functionObjects::scalarTransport::execute()
             << dimVolume/dimTime << exit(FatalError);
     }
 
+    if (converged)
+    {
+        Log << type() << ": " << name() << ": "
+            << s.name() << " is converged." << nl
+            << tab << "initial-residual tolerance: " << tol_ << nl
+            << tab << "outer iteration: " << iter << nl;
+    }
+
     Log << endl;
 
     return true;

src/functionObjects/solvers/scalarTransport/scalarTransport.H

@@ -90,6 +90,7 @@ Usage
             field           s;
             bounded01       false;
             phase           alpha.water;
+            tolerance       1e-5;
 
             write           true;
 
@@ -123,6 +124,7 @@ Usage
         phi          | Name of flux field      | no          | phi
         rho          | Name of density field   | no          | rho
         phase        | Name of the phase       | no          | none
+        tolerance    | Outer-loop initial-residual tolerance | no | 1
         nut          | Name of the turbulence viscosity | no | none
         D            | Diffusion coefficient   | no          | auto generated
         nCorr        | Number of correctors    | no          | 0
@@ -195,6 +197,9 @@ class scalarTransport
         //- Turbulent diffusion coefficient (optional)
         scalar alphaDt_;
 
+        //- Outer-loop initial-residual tolerance
+        scalar tol_;
+
         //- Number of corrector iterations (optional)
         label nCorr_;
 
@@ -223,6 +228,7 @@ class scalarTransport
             const surfaceScalarField& phi
         ) const;
 
+
         //- No copy construct
         scalarTransport(const scalarTransport&) = delete;