ENH: turbulence: better restart behaviour. Fixes #1172.

src/TurbulenceModels/incompressible/turbulentTransportModels/RAS/LamBremhorstKE/LamBremhorstKE.C

@@ -3,7 +3,7 @@
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
     \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
-     \\/     M anipulation  |
+     \\/     M anipulation  | Copyright (C) 2019 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -254,7 +254,8 @@ void LamBremhorstKE::correct()
     solve(kEqn);
     bound(k_, kMin_);
 
-    correctNut(fMu);
+    // Update nut with latest available k,epsilon
+    correctNut();
 }
 
 

src/TurbulenceModels/turbulenceModels/RAS/SpalartAllmaras/SpalartAllmaras.C

@@ -3,7 +3,7 @@
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
     \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
-     \\/     M anipulation  |
+     \\/     M anipulation  | Copyright (C) 2019 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -353,39 +353,42 @@ void SpalartAllmaras<BasicTurbulenceModel>::correct()
         return;
     }
 
-    // Local references
-    const alphaField& alpha = this->alpha_;
-    const rhoField& rho = this->rho_;
-    const surfaceScalarField& alphaRhoPhi = this->alphaRhoPhi_;
-    fv::options& fvOptions(fv::options::New(this->mesh_));
-
-    eddyViscosity<RASModel<BasicTurbulenceModel>>::correct();
+    {
+        // Local references
+        const alphaField& alpha = this->alpha_;
+        const rhoField& rho = this->rho_;
+        const surfaceScalarField& alphaRhoPhi = this->alphaRhoPhi_;
+        fv::options& fvOptions(fv::options::New(this->mesh_));
 
-    const volScalarField chi(this->chi());
-    const volScalarField fv1(this->fv1(chi));
+        eddyViscosity<RASModel<BasicTurbulenceModel>>::correct();
 
-    const volScalarField Stilda(this->Stilda(chi, fv1));
+        const volScalarField chi(this->chi());
+        const volScalarField fv1(this->fv1(chi));
 
-    tmp<fvScalarMatrix> nuTildaEqn
-    (
-        fvm::ddt(alpha, rho, nuTilda_)
-      + fvm::div(alphaRhoPhi, nuTilda_)
-      - fvm::laplacian(alpha*rho*DnuTildaEff(), nuTilda_)
-      - Cb2_/sigmaNut_*alpha*rho*magSqr(fvc::grad(nuTilda_))
-     ==
-        Cb1_*alpha*rho*Stilda*nuTilda_
-      - fvm::Sp(Cw1_*alpha*rho*fw(Stilda)*nuTilda_/sqr(y_), nuTilda_)
-      + fvOptions(alpha, rho, nuTilda_)
-    );
+        const volScalarField Stilda(this->Stilda(chi, fv1));
 
-    nuTildaEqn.ref().relax();
-    fvOptions.constrain(nuTildaEqn.ref());
-    solve(nuTildaEqn);
-    fvOptions.correct(nuTilda_);
-    bound(nuTilda_, dimensionedScalar(nuTilda_.dimensions(), Zero));
-    nuTilda_.correctBoundaryConditions();
+        tmp<fvScalarMatrix> nuTildaEqn
+        (
+            fvm::ddt(alpha, rho, nuTilda_)
+          + fvm::div(alphaRhoPhi, nuTilda_)
+          - fvm::laplacian(alpha*rho*DnuTildaEff(), nuTilda_)
+          - Cb2_/sigmaNut_*alpha*rho*magSqr(fvc::grad(nuTilda_))
+         ==
+            Cb1_*alpha*rho*Stilda*nuTilda_
+          - fvm::Sp(Cw1_*alpha*rho*fw(Stilda)*nuTilda_/sqr(y_), nuTilda_)
+          + fvOptions(alpha, rho, nuTilda_)
+        );
+
+        nuTildaEqn.ref().relax();
+        fvOptions.constrain(nuTildaEqn.ref());
+        solve(nuTildaEqn);
+        fvOptions.correct(nuTilda_);
+        bound(nuTilda_, dimensionedScalar(nuTilda_.dimensions(), Zero));
+        nuTilda_.correctBoundaryConditions();
+    }
 
-    correctNut(fv1);
+    // Update nut with latest available k,epsilon
+    correctNut();
 }
 
 

src/TurbulenceModels/turbulenceModels/derivedFvPatchFields/wallFunctions/nutWallFunctions/nutUBlendedWallFunction/nutUBlendedWallFunctionFvPatchScalarField.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2016-2017 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2016-2019 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -67,6 +67,10 @@ Note
     The full 'automatic wall treatment' description also requires use of the
     Foam::omegaWallFunction with the \c blended flag set to 'on'
 
+    Suffers from non-exact restart since correctNut() (called through
+    turbulence->validate) returns a slightly different value every time
+    it is called. See nutUSpaldingWallFunctionFvPatchScalarField.C
+
 SeeAlso
     Foam::nutWallFunctionFvPatchScalarField
     Foam::omegaWallFunctionFvPatchScalarField

src/TurbulenceModels/turbulenceModels/derivedFvPatchFields/wallFunctions/nutWallFunctions/nutURoughWallFunction/nutURoughWallFunctionFvPatchScalarField.H

@@ -3,7 +3,7 @@
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
     \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
-     \\/     M anipulation  |
+     \\/     M anipulation  | Copyright (C) 2019 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -50,6 +50,11 @@ Usage
     }
     \endverbatim
 
+Note
+    Suffers from non-exact restart since correctNut() (called through
+    turbulence->validate) returns a slightly different value every time
+    it is called. See nutUSpaldingWallFunctionFvPatchScalarField.C
+
 See also
     Foam::nutWallFunctionFvPatchScalarField
 

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

@@ -3,7 +3,7 @@
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
     \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
-     \\/     M anipulation  |
+     \\/     M anipulation  | Copyright (C) 2019 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -92,6 +92,8 @@ tmp<scalarField> nutUSpaldingWallFunctionFvPatchScalarField::calcUTau
     forAll(uTau, facei)
     {
         scalar ut = sqrt((nutw[facei] + nuw[facei])*magGradU[facei]);
+        // Note: for exact restart seed with laminar viscosity only:
+        //scalar ut = sqrt(nuw[facei]*magGradU[facei]);
 
         if (ut > ROOTVSMALL)
         {
@@ -117,9 +119,19 @@ tmp<scalarField> nutUSpaldingWallFunctionFvPatchScalarField::calcUTau
                 err = mag((ut - uTauNew)/ut);
                 ut = uTauNew;
 
+                if (debug)
+                {
+                    iterations_++;
+                }
+
             } while (ut > ROOTVSMALL && err > 0.01 && ++iter < 10);
 
             uTau[facei] = max(0.0, ut);
+
+            if (debug)
+            {
+                invocations_++;
+            }
         }
     }
 
@@ -136,7 +148,9 @@ nutUSpaldingWallFunctionFvPatchScalarField
     const DimensionedField<scalar, volMesh>& iF
 )
 :
-    nutWallFunctionFvPatchScalarField(p, iF)
+    nutWallFunctionFvPatchScalarField(p, iF),
+    invocations_(0),
+    iterations_(0)
 {}
 
 
@@ -149,7 +163,9 @@ nutUSpaldingWallFunctionFvPatchScalarField
     const fvPatchFieldMapper& mapper
 )
 :
-    nutWallFunctionFvPatchScalarField(ptf, p, iF, mapper)
+    nutWallFunctionFvPatchScalarField(ptf, p, iF, mapper),
+    invocations_(0),
+    iterations_(0)
 {}
 
 
@@ -161,7 +177,9 @@ nutUSpaldingWallFunctionFvPatchScalarField
     const dictionary& dict
 )
 :
-    nutWallFunctionFvPatchScalarField(p, iF, dict)
+    nutWallFunctionFvPatchScalarField(p, iF, dict),
+    invocations_(0),
+    iterations_(0)
 {}
 
 
@@ -171,7 +189,9 @@ nutUSpaldingWallFunctionFvPatchScalarField
     const nutUSpaldingWallFunctionFvPatchScalarField& wfpsf
 )
 :
-    nutWallFunctionFvPatchScalarField(wfpsf)
+    nutWallFunctionFvPatchScalarField(wfpsf),
+    invocations_(wfpsf.invocations_),
+    iterations_(wfpsf.iterations_)
 {}
 
 
@@ -182,10 +202,29 @@ nutUSpaldingWallFunctionFvPatchScalarField
     const DimensionedField<scalar, volMesh>& iF
 )
 :
-    nutWallFunctionFvPatchScalarField(wfpsf, iF)
+    nutWallFunctionFvPatchScalarField(wfpsf, iF),
+    invocations_(0),
+    iterations_(0)
 {}
 
 
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+nutUSpaldingWallFunctionFvPatchScalarField::
+~nutUSpaldingWallFunctionFvPatchScalarField()
+{
+    if (debug)
+    {
+        Info<< "nutUSpaldingWallFunctionFvPatchScalarField :"
+            << " total invocations:"
+            << returnReduce(invocations_, sumOp<label>())
+            << " total iterations:"
+            << returnReduce(iterations_, sumOp<label>())
+            << endl;
+    }
+}
+
+
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
 tmp<scalarField> nutUSpaldingWallFunctionFvPatchScalarField::yPlus() const

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

@@ -3,7 +3,7 @@
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
     \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
-     \\/     M anipulation  |
+     \\/     M anipulation  | Copyright (C) 2019 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -57,6 +57,13 @@ Usage
 See also
     Foam::nutWallFunctionFvPatchScalarField
 
+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.
+
 SourceFiles
     nutUSpaldingWallFunctionFvPatchScalarField.C
 
@@ -66,6 +73,7 @@ SourceFiles
 #define nutUSpaldingWallFunctionFvPatchScalarField_H
 
 #include "nutWallFunctionFvPatchScalarField.H"
+#include "uint64.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -82,6 +90,12 @@ class nutUSpaldingWallFunctionFvPatchScalarField
 {
 protected:
 
+    // Protected data
+
+        mutable uint64_t invocations_;
+        mutable uint64_t iterations_;
+
+
     // Protected Member Functions
 
         //- Calculate the turbulence viscosity
@@ -160,6 +174,10 @@ public:
         }
 
 
+    //- Destructor
+    virtual ~nutUSpaldingWallFunctionFvPatchScalarField();
+
+
     // Member functions
 
         // Evaluation functions

src/TurbulenceModels/turbulenceModels/derivedFvPatchFields/wallFunctions/nutWallFunctions/nutUWallFunction/nutUWallFunctionFvPatchScalarField.H

@@ -3,7 +3,7 @@
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
     \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
-     \\/     M anipulation  |
+     \\/     M anipulation  | Copyright (C) 2019 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -43,6 +43,11 @@ Usage
 See also
     Foam::nutWallFunctionFvPatchScalarField
 
+Note
+    Suffers from non-exact restart since correctNut() (called through
+    turbulence->validate) returns a slightly different value every time
+    it is called. See nutUSpaldingWallFunctionFvPatchScalarField.C
+
 SourceFiles
     nutUWallFunctionFvPatchScalarField.C