ENH: exact restart. Fixes #1172.

src/TurbulenceModels/turbulenceModels/derivedFvPatchFields/wallFunctions/nutWallFunctions/nutUSpaldingWallFunction/nutUSpaldingWallFunctionFvPatchScalarField.C

@@ -53,11 +53,39 @@ tmp<scalarField> nutUSpaldingWallFunctionFvPatchScalarField::calcNut() const
     const tmp<scalarField> tnuw = turbModel.nu(patchi);
     const scalarField& nuw = tnuw();
 
-    return max
+
+    // Calculate new viscosity
+    tmp<scalarField> tnutw
     (
-        scalar(0),
-        sqr(calcUTau(magGradU))/(magGradU + ROOTVSMALL) - nuw
+        max
+        (
+            scalar(0),
+            sqr(calcUTau(magGradU))/(magGradU + ROOTVSMALL) - nuw
+        )
     );
+
+    if (tolerance_ != 0.01)
+    {
+        // User-specified tolerance. Find out if current nut already satisfies
+        // eqns.
+
+        // Run ut for one iteration
+        scalarField err;
+        tmp<scalarField> UTau(calcUTau(magGradU, 1, err));
+
+        // Preserve nutw if the initial error (err) already lower than the
+        // tolerance.
+
+        scalarField& nutw = tnutw.ref();
+        forAll(err, facei)
+        {
+            if (err[facei] < tolerance_)
+            {
+                nutw[facei] = this->operator[](facei);
+            }
+        }
+    }
+    return tnutw;
 }
 
 
@@ -65,6 +93,18 @@ tmp<scalarField> nutUSpaldingWallFunctionFvPatchScalarField::calcUTau
 (
     const scalarField& magGradU
 ) const
+{
+    scalarField err;
+    return calcUTau(magGradU, maxIter_, err);
+}
+
+
+tmp<scalarField> nutUSpaldingWallFunctionFvPatchScalarField::calcUTau
+(
+    const scalarField& magGradU,
+    const label maxIter,
+    scalarField& err
+) const
 {
     const label patchi = patch().index();
 
@@ -89,6 +129,9 @@ tmp<scalarField> nutUSpaldingWallFunctionFvPatchScalarField::calcUTau
     tmp<scalarField> tuTau(new scalarField(patch().size(), Zero));
     scalarField& uTau = tuTau.ref();
 
+    err.setSize(uTau.size());
+    err = 0.0;
+
     forAll(uTau, facei)
     {
         scalar ut = sqrt((nutw[facei] + nuw[facei])*magGradU[facei]);
@@ -98,7 +141,6 @@ tmp<scalarField> nutUSpaldingWallFunctionFvPatchScalarField::calcUTau
         if (ut > ROOTVSMALL)
         {
             int iter = 0;
-            scalar err = GREAT;
 
             do
             {
@@ -116,12 +158,17 @@ tmp<scalarField> nutUSpaldingWallFunctionFvPatchScalarField::calcUTau
                   + 1/E_*kUu*fkUu/ut;
 
                 scalar uTauNew = ut + f/df;
-                err = mag((ut - uTauNew)/ut);
+                err[facei] = mag((ut - uTauNew)/ut);
                 ut = uTauNew;
 
                 //iterations_++;
 
-            } while (ut > ROOTVSMALL && err > tolerance_ && ++iter < maxIter_);
+            } while
+            (
+                ut > ROOTVSMALL
+             && err[facei] > tolerance_
+             && ++iter < maxIter
+            );
 
             uTau[facei] = max(0.0, ut);
 

src/TurbulenceModels/turbulenceModels/derivedFvPatchFields/wallFunctions/nutWallFunctions/nutUSpaldingWallFunction/nutUSpaldingWallFunctionFvPatchScalarField.H

@@ -71,9 +71,18 @@ Note
     Suffers from non-exact restart since correctNut() (called through
     turbulence->validate) returns a slightly different value every time
     it is called. This is since the seed for the Newton-Raphson iteration
-    uses the current value of *this (= nut). Can be avoided by seeding the
-    NR with e.g. the laminar viscosity or tightening the convergence tolerance
-    to e.g. 1e-7 and the max number of iterations to 100.
+    uses the current value of *this (= nut).
+
+    This can be avoided by overriding the tolerance. This also switches on
+    a pre-detection whether the current nut already satisfies the turbulence
+    conditions and if so does not change it at all. This means that the nut
+    only changes if it 'has really changed'. This probably should be used with
+    a tight tolerance, e.g.
+
+        maxIter     100;
+        tolerance   1e-7;
+
+    to make sure to kick every iteration.
 
 SourceFiles
     nutUSpaldingWallFunctionFvPatchScalarField.C
@@ -123,6 +132,15 @@ protected:
         //- Calculate the friction velocity
         virtual tmp<scalarField> calcUTau(const scalarField& magGradU) const;
 
+        //- Calculate the friction velocity and number of iterations for
+        //  convergence
+        virtual tmp<scalarField> calcUTau
+        (
+            const scalarField& magGradU,
+            const label maxIter,
+            scalarField& err
+        ) const;
+
         //- Write local wall function variables
         virtual void writeLocalEntries(Ostream&) const;