Merge branch 'feature-secondOrderRestart' into 'develop'

Feature second order restart

See merge request !114

applications/solvers/combustion/PDRFoam/createFields.H

@@ -58,23 +58,6 @@ autoPtr<compressible::RASModel> turbulence
 );
 
 
-Info<< "Creating field dpdt\n" << endl;
-volScalarField dpdt
-(
-    IOobject
-    (
-        "dpdt",
-        runTime.timeName(),
-        mesh
-    ),
-    mesh,
-    dimensionedScalar("dpdt", p.dimensions()/dimTime, 0)
-);
-
-Info<< "Creating field kinetic energy K\n" << endl;
-volScalarField K("K", 0.5*magSqr(U));
-
-
 Info<< "Creating the unstrained laminar flame speed\n" << endl;
 autoPtr<laminarFlameSpeed> unstrainedLaminarFlameSpeed
 (
@@ -227,4 +210,8 @@ fields.add(thermo.he());
 fields.add(thermo.heu());
 flameWrinkling->addXi(fields);
 
+#include "createDpdt.H"
+
+#include "createK.H"
+
 #include "createMRF.H"

applications/solvers/combustion/XiFoam/createFields.H

@@ -58,21 +58,6 @@ autoPtr<compressible::turbulenceModel> turbulence
     )
 );
 
-Info<< "Creating field dpdt\n" << endl;
-volScalarField dpdt
-(
-    IOobject
-    (
-        "dpdt",
-        runTime.timeName(),
-        mesh
-    ),
-    mesh,
-    dimensionedScalar("dpdt", p.dimensions()/dimTime, 0)
-);
-
-Info<< "Creating field kinetic energy K\n" << endl;
-volScalarField K("K", 0.5*magSqr(U));
 
 Info<< "Creating field Xi\n" << endl;
 volScalarField Xi
@@ -139,4 +124,8 @@ fields.add(b);
 fields.add(thermo.he());
 fields.add(thermo.heu());
 
+#include "createDpdt.H"
+
+#include "createK.H"
+
 #include "createMRF.H"

applications/solvers/combustion/coldEngineFoam/createFields.H

@@ -51,20 +51,8 @@ autoPtr<compressible::turbulenceModel> turbulence
     )
 );
 
-Info<< "Creating field dpdt\n" << endl;
+#include "createDpdt.H"
-volScalarField dpdt
-(
-    IOobject
-    (
-        "dpdt",
-        runTime.timeName(),
-        mesh
-    ),
-    mesh,
-    dimensionedScalar("dpdt", p.dimensions()/dimTime, 0)
-);
 
-Info<< "Creating field kinetic energy K\n" << endl;
+#include "createK.H"
-volScalarField K("K", 0.5*magSqr(U));
 
 #include "createMRF.H"

applications/solvers/combustion/fireFoam/createFields.H

@@ -145,22 +145,9 @@ volScalarField Qdot
     dimensionedScalar("Qdot", dimEnergy/dimVolume/dimTime, 0.0)
 );
 
+#include "createDpdt.H"
 
-Info<< "Creating field dpdt\n" << endl;
+#include "createK.H"
-volScalarField dpdt
-(
-    IOobject
-    (
-        "dpdt",
-        runTime.timeName(),
-        mesh
-    ),
-    mesh,
-    dimensionedScalar("dpdt", p.dimensions()/dimTime, 0)
-);
-
-Info<< "Creating field kinetic energy K\n" << endl;
-volScalarField K("K", 0.5*magSqr(U));
 
 #include "createClouds.H"
 #include "createSurfaceFilmModel.H"

applications/solvers/combustion/reactingFoam/createFields.H

@@ -71,24 +71,6 @@ autoPtr<compressible::turbulenceModel> turbulence
 // Set the turbulence into the reaction model
 reaction->setTurbulence(turbulence());
 
-
-Info<< "Creating field dpdt\n" << endl;
-volScalarField dpdt
-(
-    IOobject
-    (
-        "dpdt",
-        runTime.timeName(),
-        mesh
-    ),
-    mesh,
-    dimensionedScalar("dpdt", p.dimensions()/dimTime, 0)
-);
-
-Info<< "Creating field kinetic energy K\n" << endl;
-volScalarField K("K", 0.5*magSqr(U));
-
-
 multivariateSurfaceInterpolationScheme<scalar>::fieldTable fields;
 
 forAll(Y, i)
@@ -111,4 +93,8 @@ volScalarField Qdot
     dimensionedScalar("Qdot", dimEnergy/dimVolume/dimTime, 0.0)
 );
 
+#include "createDpdt.H"
+
+#include "createK.H"
+
 #include "createMRF.H"

applications/solvers/combustion/reactingFoam/rhoReactingBuoyantFoam/createFields.H

@@ -92,23 +92,6 @@ volScalarField p_rgh
 // Force p_rgh to be consistent with p
 p_rgh = p - rho*gh;
 
-Info<< "Creating field dpdt\n" << endl;
-volScalarField dpdt
-(
-    IOobject
-    (
-        "dpdt",
-        runTime.timeName(),
-        mesh
-    ),
-    mesh,
-    dimensionedScalar("dpdt", p.dimensions()/dimTime, 0)
-);
-
-Info<< "Creating field kinetic energy K\n" << endl;
-volScalarField K("K", 0.5*magSqr(U));
-
-
 multivariateSurfaceInterpolationScheme<scalar>::fieldTable fields;
 
 forAll(Y, i)
@@ -131,4 +114,8 @@ volScalarField Qdot
     dimensionedScalar("Qdot", dimEnergy/dimVolume/dimTime, 0.0)
 );
 
+#include "createDpdt.H"
+
+#include "createK.H"
+
 #include "createMRF.H"

applications/solvers/combustion/reactingFoam/rhoReactingFoam/createFields.H

@@ -72,24 +72,6 @@ autoPtr<compressible::turbulenceModel> turbulence
 // Set the turbulence into the reaction model
 reaction->setTurbulence(turbulence());
 
-
-Info<< "Creating field dpdt\n" << endl;
-volScalarField dpdt
-(
-    IOobject
-    (
-        "dpdt",
-        runTime.timeName(),
-        mesh
-    ),
-    mesh,
-    dimensionedScalar("dpdt", p.dimensions()/dimTime, 0)
-);
-
-Info<< "Creating field kinetic energy K\n" << endl;
-volScalarField K("K", 0.5*magSqr(U));
-
-
 multivariateSurfaceInterpolationScheme<scalar>::fieldTable fields;
 
 forAll(Y, i)
@@ -112,6 +94,10 @@ volScalarField Qdot
     dimensionedScalar("Qdot", dimEnergy/dimVolume/dimTime, 0.0)
 );
 
+#include "createDpdt.H"
+
+#include "createK.H"
+
 #include "createMRF.H"
 
 dimensionedScalar rhoMax

applications/solvers/compressible/rhoPimpleFoam/createFields.H

@@ -56,21 +56,9 @@ autoPtr<compressible::turbulenceModel> turbulence
     )
 );
 
-Info<< "Creating field dpdt\n" << endl;
+#include "createDpdt.H"
-volScalarField dpdt
-(
-    IOobject
-    (
-        "dpdt",
-        runTime.timeName(),
-        mesh
-    ),
-    mesh,
-    dimensionedScalar("dpdt", p.dimensions()/dimTime, 0)
-);
 
-Info<< "Creating field kinetic energy K\n" << endl;
+#include "createK.H"
-volScalarField K("K", 0.5*magSqr(U));
 
 #include "createMRF.H"
 

applications/solvers/compressible/sonicFoam/createFields.H

@@ -50,7 +50,6 @@ autoPtr<compressible::turbulenceModel> turbulence
     )
 );
 
-Info<< "Creating field kinetic energy K\n" << endl;
+#include "createK.H"
-volScalarField K("K", 0.5*magSqr(U));
 
 #include "createMRF.H"

applications/solvers/heatTransfer/buoyantPimpleFoam/createFields.H

@@ -73,22 +73,6 @@ p_rgh = p - rho*gh;
 
 mesh.setFluxRequired(p_rgh.name());
 
-Info<< "Creating field dpdt\n" << endl;
-volScalarField dpdt
-(
-    IOobject
-    (
-        "dpdt",
-        runTime.timeName(),
-        mesh
-    ),
-    mesh,
-    dimensionedScalar("dpdt", p.dimensions()/dimTime, 0)
-);
-
-Info<< "Creating field kinetic energy K\n" << endl;
-volScalarField K("K", 0.5*magSqr(U));
-
 label pRefCell = 0;
 scalar pRefValue = 0.0;
 
@@ -113,5 +97,9 @@ if (p_rgh.needReference())
 
 dimensionedScalar initialMass("initialMass", fvc::domainIntegrate(rho));
 
+#include "createDpdt.H"
+
+#include "createK.H"
+
 #include "createMRF.H"
 #include "createRadiationModel.H"

applications/solvers/heatTransfer/thermoFoam/createFields.H

@@ -37,20 +37,8 @@ volVectorField U
 
 #include "setAlphaEff.H"
 
-Info<< "Creating field dpdt\n" << endl;
+#include "createDpdt.H"
-volScalarField dpdt
-(
-    IOobject
-    (
-        "dpdt",
-        runTime.timeName(),
-        mesh
-    ),
-    mesh,
-    dimensionedScalar("dpdt", p.dimensions()/dimTime, 0)
-);
 
-Info<< "Creating field kinetic energy K\n" << endl;
+#include "createK.H"
-volScalarField K("K", 0.5*magSqr(U));
 
 #include "createRadiationModel.H"

applications/solvers/lagrangian/coalChemistryFoam/createFields.H

@@ -115,22 +115,6 @@ autoPtr<compressible::turbulenceModel> turbulence
 // Set the turbulence into the combustion model
 combustion->setTurbulence(turbulence());
 
-Info<< "Creating field dpdt\n" << endl;
-volScalarField dpdt
-(
-    IOobject
-    (
-        "dpdt",
-        runTime.timeName(),
-        mesh
-    ),
-    mesh,
-    dimensionedScalar("dpdt", p.dimensions()/dimTime, 0)
-);
-
-Info<< "Creating field kinetic energy K\n" << endl;
-volScalarField K("K", 0.5*magSqr(U));
-
 volScalarField Qdot
 (
     IOobject
@@ -145,6 +129,10 @@ volScalarField Qdot
     dimensionedScalar("Qdot", dimEnergy/dimVolume/dimTime, 0.0)
 );
 
+#include "createDpdt.H"
+
+#include "createK.H"
+
 #include "createMRF.H"
 #include "createClouds.H"
 #include "createRadiationModel.H"

applications/solvers/lagrangian/reactingParcelFilmFoam/createFields.H

@@ -69,23 +69,6 @@ autoPtr<compressible::turbulenceModel> turbulence
 // Set the turbulence into the combustion model
 combustion->setTurbulence(turbulence());
 
-Info<< "Creating field dpdt\n" << endl;
-volScalarField dpdt
-(
-    IOobject
-    (
-        "dpdt",
-        runTime.timeName(),
-        mesh
-    ),
-    mesh,
-    dimensionedScalar("dpdt", p.dimensions()/dimTime, 0)
-);
-
-Info<< "Creating field kinetic energy K\n" << endl;
-volScalarField K("K", 0.5*magSqr(U));
-
-
 #include "readGravitationalAcceleration.H"
 #include "readhRef.H"
 #include "gh.H"
@@ -148,6 +131,10 @@ volScalarField Qdot
     dimensionedScalar("Qdot", dimEnergy/dimVolume/dimTime, 0.0)
 );
 
+#include "createDpdt.H"
+
+#include "createK.H"
+
 #include "createMRF.H"
 #include "createClouds.H"
 #include "createRadiationModel.H"

applications/solvers/lagrangian/reactingParcelFoam/createFields.H

@@ -96,22 +96,6 @@ autoPtr<compressible::turbulenceModel> turbulence
 // Set the turbulence into the combustion model
 combustion->setTurbulence(turbulence());
 
-Info<< "Creating field dpdt\n" << endl;
-volScalarField dpdt
-(
-    IOobject
-    (
-        "dpdt",
-        runTime.timeName(),
-        mesh
-    ),
-    mesh,
-    dimensionedScalar("dpdt", p.dimensions()/dimTime, 0)
-);
-
-Info<< "Creating field kinetic energy K\n" << endl;
-volScalarField K("K", 0.5*magSqr(U));
-
 Info<< "Creating multi-variate interpolation scheme\n" << endl;
 multivariateSurfaceInterpolationScheme<scalar>::fieldTable fields;
 
@@ -135,6 +119,10 @@ volScalarField Qdot
     dimensionedScalar("Qdot", dimEnergy/dimVolume/dimTime, 0.0)
 );
 
+#include "createDpdt.H"
+
+#include "createK.H"
+
 #include "createMRF.H"
 #include "createRadiationModel.H"
 #include "createClouds.H"

applications/solvers/lagrangian/sprayFoam/createFields.H

@@ -94,21 +94,6 @@ autoPtr<compressible::turbulenceModel> turbulence
 // Set the turbulence into the combustion model
 combustion->setTurbulence(turbulence());
 
-Info<< "Creating field dpdt\n" << endl;
-volScalarField dpdt
-(
-    IOobject
-    (
-        "dpdt",
-        runTime.timeName(),
-        mesh
-    ),
-    mesh,
-    dimensionedScalar("dpdt", p.dimensions()/dimTime, 0)
-);
-
-Info<< "Creating field kinetic energy K\n" << endl;
-volScalarField K("K", 0.5*magSqr(U));
 multivariateSurfaceInterpolationScheme<scalar>::fieldTable fields;
 
 forAll(Y, i)
@@ -131,6 +116,10 @@ volScalarField Qdot
     dimensionedScalar("Qdot", dimEnergy/dimVolume/dimTime, 0.0)
 );
 
+#include "createDpdt.H"
+
+#include "createK.H"
+
 #include "createMRF.H"
 #include "createClouds.H"
 #include "createRadiationModel.H"

applications/solvers/multiphase/compressibleInterFoam/createFields.H

@@ -98,7 +98,6 @@ autoPtr<compressible::turbulenceModel> turbulence
     compressible::turbulenceModel::New(rho, U, rhoPhi, mixture)
 );
 
-Info<< "Creating field kinetic energy K\n" << endl;
+#include "createK.H"
-volScalarField K("K", 0.5*magSqr(U));
 
 #include "createMRF.H"

applications/solvers/multiphase/compressibleMultiphaseInterFoam/createFields.H

@@ -65,5 +65,4 @@ autoPtr<compressible::turbulenceModel> turbulence
     )
 );
 
-Info<< "Creating field kinetic energy K\n" << endl;
+#include "createK.H"
-volScalarField K("K", 0.5*magSqr(U));

applications/solvers/multiphase/interCondensingEvaporatingFoam/createFields.H

@@ -124,9 +124,6 @@ volScalarField kappaEff
     thermo->kappa()
 );
 
-Info<< "Creating field kinetic energy K\n" << endl;
-volScalarField K("K", 0.5*magSqr(U));
-
 Info<< "Creating field pDivU\n" << endl;
 volScalarField pDivU
 (

applications/solvers/multiphase/interCondensingEvaporatingFoam/pEqn.H

@@ -47,7 +47,6 @@
             U = HbyA + rAU*fvc::reconstruct((phig + p_rghEqn.flux())/rAUf);
             U.correctBoundaryConditions();
             fvOptions.correct(U);
-            K = 0.5*magSqr(U);
         }
     }
 

src/OpenFOAM/fields/GeometricFields/GeometricField/GeometricField.C

@@ -770,10 +770,11 @@ void Foam::GeometricField<Type, PatchField, GeoMesh>::storeOldTimes() const
     )
     {
         storeOldTime();
+        timeIndex_ = this->time().timeIndex();
     }
 
     // Correct time index
-    timeIndex_ = this->time().timeIndex();
+    //timeIndex_ = this->time().timeIndex();
 }
 
 
@@ -835,6 +836,17 @@ Foam::GeometricField<Type, PatchField, GeoMesh>::oldTime() const
             ),
             *this
         );
+
+        if (debug)
+        {
+            InfoInFunction
+                << "created old time field " << field0Ptr_->info() << endl;
+
+            if (debug&2)
+            {
+                error::printStack(Info);
+            }
+        }
     }
     else
     {

src/finiteVolume/cfdTools/compressible/createDpdt.H

@@ -0,0 +1,18 @@
+Info<< "Creating field dpdt\n" << endl;
+volScalarField dpdt
+(
+    IOobject
+    (
+        "dpdt",
+        runTime.timeName(),
+        mesh,
+        IOobject::NO_READ,
+        IOobject::NO_WRITE
+    ),
+    fvc::ddt(p)
+);
+
+if (!thermo.dpdt())
+{
+    dpdt == dimensionedScalar("0", dpdt.dimensions(), 0);
+}

src/finiteVolume/cfdTools/general/include/createK.H

@@ -0,0 +1,16 @@
+Info<< "Creating field kinetic energy K\n" << endl;
+volScalarField K("K", 0.5*magSqr(U));
+
+if (U.nOldTimes())
+{
+    volVectorField* Uold = &U.oldTime();
+    volScalarField* Kold = &K.oldTime();
+    *Kold == 0.5*magSqr(*Uold);
+
+    while (Uold->nOldTimes())
+    {
+        Uold = &Uold->oldTime();
+        Kold = &Kold->oldTime();
+        *Kold == 0.5*magSqr(*Uold);
+    }
+}

src/finiteVolume/finiteVolume/ddtSchemes/backwardDdtScheme/backwardDdtScheme.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) 2017 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -58,7 +58,7 @@ template<class Type>
 template<class GeoField>
 scalar backwardDdtScheme<Type>::deltaT0_(const GeoField& vf) const
 {
-    if (vf.nOldTimes() < 2)
+    if (mesh().time().timeIndex() < 2)
     {
         return GREAT;
     }

src/finiteVolume/finiteVolume/ddtSchemes/backwardDdtScheme/backwardDdtScheme.H

@@ -92,27 +92,13 @@ public:
         backwardDdtScheme(const fvMesh& mesh)
         :
             ddtScheme<Type>(mesh)
-        {
+        {}
-            // Ensure the old-old-time cell volumes are available
-            // for moving meshes
-            if (mesh.moving())
-            {
-                mesh.V00();
-            }
-        }
 
         //- Construct from mesh and Istream
         backwardDdtScheme(const fvMesh& mesh, Istream& is)
         :
             ddtScheme<Type>(mesh, is)
-        {
+        {}
-            // Ensure the old-old-time cell volumes are available
-            // for moving meshes
-            if (mesh.moving())
-            {
-                mesh.V00();
-            }
-        }
 
 
     // Member Functions

src/finiteVolume/fvMesh/fvMesh.C

@@ -198,8 +198,8 @@ void Foam::fvMesh::storeOldVol(const scalarField& V)
                     time().timeName(),
                     *this,
                     IOobject::NO_READ,
-                    IOobject::NO_WRITE,
+                    IOobject::AUTO_WRITE,
-                    false
+                    true
                 ),
                 *this,
                 dimVolume
@@ -281,8 +281,8 @@ Foam::fvMesh::fvMesh(const IOobject& io)
                 time().timeName(),
                 *this,
                 IOobject::MUST_READ,
-                IOobject::NO_WRITE,
+                IOobject::AUTO_WRITE,
-                false
+                true
             ),
             *this
         );

src/finiteVolume/fvMesh/fvMeshGeometry.C

@@ -260,11 +260,12 @@ const Foam::volScalarField::Internal& Foam::fvMesh::V00() const
                 *this,
                 IOobject::NO_READ,
                 IOobject::NO_WRITE,
-                false
+                true
             ),
             V0()
         );
 
+
         // If V00 is used then V0 should be stored for restart
         V0Ptr_->writeOpt() = IOobject::AUTO_WRITE;
     }

src/thermophysicalModels/basic/heThermo/heThermo.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) 2015-2017 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -52,11 +52,16 @@ heBoundaryCorrection(volScalarField& h)
 
 
 template<class BasicThermo, class MixtureType>
-void Foam::heThermo<BasicThermo, MixtureType>::init()
+void Foam::heThermo<BasicThermo, MixtureType>::init
+(
+    const volScalarField& p,
+    const volScalarField& T,
+    volScalarField& he
+)
 {
-    scalarField& heCells = he_.primitiveFieldRef();
+    scalarField& heCells = he.primitiveFieldRef();
-    const scalarField& pCells = this->p_;
+    const scalarField& pCells = p.primitiveField();
-    const scalarField& TCells = this->T_;
+    const scalarField& TCells = T.primitiveField();
 
     forAll(heCells, celli)
     {
@@ -64,19 +69,25 @@ void Foam::heThermo<BasicThermo, MixtureType>::init()
             this->cellMixture(celli).HE(pCells[celli], TCells[celli]);
     }
 
-    volScalarField::Boundary& heBf = he_.boundaryFieldRef();
+    volScalarField::Boundary& heBf = he.boundaryFieldRef();
 
     forAll(heBf, patchi)
     {
-        heBf[patchi] == he
+        heBf[patchi] == this->he
         (
-            this->p_.boundaryField()[patchi],
+            p.boundaryField()[patchi],
-            this->T_.boundaryField()[patchi],
+            T.boundaryField()[patchi],
             patchi
         );
     }
 
-    this->heBoundaryCorrection(he_);
+    this->heBoundaryCorrection(he);
+
+    // Note: T does not have oldTime
+    if (p.nOldTimes() > 0)
+    {
+        init(p.oldTime(), T.oldTime(), he.oldTime());
+    }
 }
 
 
@@ -112,7 +123,7 @@ Foam::heThermo<BasicThermo, MixtureType>::heThermo
         this->heBoundaryBaseTypes()
     )
 {
-    init();
+    init(this->p_, this->T_, he_);
 }
 
 
@@ -146,7 +157,7 @@ Foam::heThermo<BasicThermo, MixtureType>::heThermo
         this->heBoundaryBaseTypes()
     )
 {
-    init();
+    init(this->p_, this->T_, he_);
 }
 
 

src/thermophysicalModels/basic/heThermo/heThermo.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) 2015-2017 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -77,7 +77,12 @@ private:
 
 
         //- Initialize heThermo
-        void init();
+        void init
+        (
+            const volScalarField& p,
+            const volScalarField& T,
+            volScalarField& he
+        );
 
 
 public:

src/thermophysicalModels/basic/psiThermo/hePsiThermo.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) 2015-2017 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -28,15 +28,40 @@ License
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
 template<class BasicPsiThermo, class MixtureType>
-void Foam::hePsiThermo<BasicPsiThermo, MixtureType>::calculate()
+void Foam::hePsiThermo<BasicPsiThermo, MixtureType>::calculate
+(
+    const volScalarField& p,
+    volScalarField& T,
+    volScalarField& he,
+    volScalarField& psi,
+    volScalarField& mu,
+    volScalarField& alpha,
+    const bool doOldTimes
+)
 {
-    const scalarField& hCells = this->he_;
+    // Note: update oldTimes before current time so that if T.oldTime() is
-    const scalarField& pCells = this->p_;
+    // created from T, it starts from the unconverted T
+    if (doOldTimes && (p.nOldTimes() || T.nOldTimes()))
+    {
+        calculate
+        (
+            p.oldTime(),
+            T.oldTime(),
+            he.oldTime(),
+            psi.oldTime(),
+            mu.oldTime(),
+            alpha.oldTime(),
+            true
+        );
+    }
 
-    scalarField& TCells = this->T_.primitiveFieldRef();
+    const scalarField& hCells = he.primitiveField();
-    scalarField& psiCells = this->psi_.primitiveFieldRef();
+    const scalarField& pCells = p.primitiveField();
-    scalarField& muCells = this->mu_.primitiveFieldRef();
+
-    scalarField& alphaCells = this->alpha_.primitiveFieldRef();
+    scalarField& TCells = T.primitiveFieldRef();
+    scalarField& psiCells = psi.primitiveFieldRef();
+    scalarField& muCells = mu.primitiveFieldRef();
+    scalarField& alphaCells = alpha.primitiveFieldRef();
 
     forAll(TCells, celli)
     {
@@ -56,27 +81,16 @@ void Foam::hePsiThermo<BasicPsiThermo, MixtureType>::calculate()
         alphaCells[celli] = mixture_.alphah(pCells[celli], TCells[celli]);
     }
 
-    volScalarField::Boundary& pBf =
+    const volScalarField::Boundary& pBf = p.boundaryField();
-        this->p_.boundaryFieldRef();
+    volScalarField::Boundary& TBf = T.boundaryFieldRef();
+    volScalarField::Boundary& psiBf = psi.boundaryFieldRef();
+    volScalarField::Boundary& heBf = he.boundaryFieldRef();
+    volScalarField::Boundary& muBf = mu.boundaryFieldRef();
+    volScalarField::Boundary& alphaBf = alpha.boundaryFieldRef();
 
-    volScalarField::Boundary& TBf =
+    forAll(pBf, patchi)
-        this->T_.boundaryFieldRef();
-
-    volScalarField::Boundary& psiBf =
-        this->psi_.boundaryFieldRef();
-
-    volScalarField::Boundary& heBf =
-        this->he().boundaryFieldRef();
-
-    volScalarField::Boundary& muBf =
-        this->mu_.boundaryFieldRef();
-
-    volScalarField::Boundary& alphaBf =
-        this->alpha_.boundaryFieldRef();
-
-    forAll(this->T_.boundaryField(), patchi)
     {
-        fvPatchScalarField& pp = pBf[patchi];
+        const fvPatchScalarField& pp = pBf[patchi];
         fvPatchScalarField& pT = TBf[patchi];
         fvPatchScalarField& ppsi = psiBf[patchi];
         fvPatchScalarField& phe = heBf[patchi];
@@ -126,10 +140,16 @@ Foam::hePsiThermo<BasicPsiThermo, MixtureType>::hePsiThermo
 :
     heThermo<BasicPsiThermo, MixtureType>(mesh, phaseName)
 {
-    calculate();
+    calculate
-
+    (
-    // Switch on saving old time
+        this->p_,
-    this->psi_.oldTime();
+        this->T_,
+        this->he_,
+        this->psi_,
+        this->mu_,
+        this->alpha_,
+        true                    // Create old time fields
+    );
 }
 
 
@@ -145,20 +165,20 @@ Foam::hePsiThermo<BasicPsiThermo, MixtureType>::~hePsiThermo()
 template<class BasicPsiThermo, class MixtureType>
 void Foam::hePsiThermo<BasicPsiThermo, MixtureType>::correct()
 {
-    if (debug)
+    DebugInFunction << endl;
-    {
-        InfoInFunction << endl;
-    }
 
-    // force the saving of the old-time values
+    calculate
-    this->psi_.oldTime();
+    (
+        this->p_,
+        this->T_,
+        this->he_,
+        this->psi_,
+        this->mu_,
+        this->alpha_,
+        false           // No need to update old times
+    );
 
-    calculate();
+    DebugInFunction << "Finished" << endl;
-
-    if (debug)
-    {
-        Info<< "    Finished" << endl;
-    }
 }
 
 

src/thermophysicalModels/basic/psiThermo/hePsiThermo.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) 2015-2017 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -55,11 +55,21 @@ class hePsiThermo
     // Private Member Functions
 
         //- Calculate the thermo variables
-        void calculate();
+        void calculate
+        (
+            const volScalarField& p,
+            volScalarField& T,
+            volScalarField& he,
+            volScalarField& psi,
+            volScalarField& mu,
+            volScalarField& alpha,
+            const bool doOldTimes
+        );
 
         //- Construct as copy (not implemented)
         hePsiThermo(const hePsiThermo<BasicPsiThermo, MixtureType>&);
 
+
 public:
 
     //- Runtime type information

src/thermophysicalModels/basic/rhoThermo/heRhoThermo.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) 2015-2017 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -28,16 +28,43 @@ License
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
 template<class BasicPsiThermo, class MixtureType>
-void Foam::heRhoThermo<BasicPsiThermo, MixtureType>::calculate()
+void Foam::heRhoThermo<BasicPsiThermo, MixtureType>::calculate
+(
+    const volScalarField& p,
+    volScalarField& T,
+    volScalarField& he,
+    volScalarField& psi,
+    volScalarField& rho,
+    volScalarField& mu,
+    volScalarField& alpha,
+    const bool doOldTimes
+)
 {
-    const scalarField& hCells = this->he();
+    // Note: update oldTimes before current time so that if T.oldTime() is
-    const scalarField& pCells = this->p_;
+    // created from T, it starts from the unconverted T
+    if (doOldTimes && (p.nOldTimes() || T.nOldTimes()))
+    {
+        calculate
+        (
+            p.oldTime(),
+            T.oldTime(),
+            he.oldTime(),
+            psi.oldTime(),
+            rho.oldTime(),
+            mu.oldTime(),
+            alpha.oldTime(),
+            true
+        );
+    }
 
-    scalarField& TCells = this->T_.primitiveFieldRef();
+    const scalarField& hCells = he.primitiveField();
-    scalarField& psiCells = this->psi_.primitiveFieldRef();
+    const scalarField& pCells = p.primitiveField();
-    scalarField& rhoCells = this->rho_.primitiveFieldRef();
+
-    scalarField& muCells = this->mu_.primitiveFieldRef();
+    scalarField& TCells = T.primitiveFieldRef();
-    scalarField& alphaCells = this->alpha_.primitiveFieldRef();
+    scalarField& psiCells = psi.primitiveFieldRef();
+    scalarField& rhoCells = rho.primitiveFieldRef();
+    scalarField& muCells = mu.primitiveFieldRef();
+    scalarField& alphaCells = alpha.primitiveFieldRef();
 
     forAll(TCells, celli)
     {
@@ -58,30 +85,17 @@ void Foam::heRhoThermo<BasicPsiThermo, MixtureType>::calculate()
         alphaCells[celli] = mixture_.alphah(pCells[celli], TCells[celli]);
     }
 
-    volScalarField::Boundary& pBf =
+    const volScalarField::Boundary& pBf = p.boundaryField();
-        this->p_.boundaryFieldRef();
+    volScalarField::Boundary& TBf = T.boundaryFieldRef();
+    volScalarField::Boundary& psiBf = psi.boundaryFieldRef();
+    volScalarField::Boundary& rhoBf = rho.boundaryFieldRef();
+    volScalarField::Boundary& heBf = he.boundaryFieldRef();
+    volScalarField::Boundary& muBf = mu.boundaryFieldRef();
+    volScalarField::Boundary& alphaBf = alpha.boundaryFieldRef();
 
-    volScalarField::Boundary& TBf =
+    forAll(pBf, patchi)
-        this->T_.boundaryFieldRef();
-
-    volScalarField::Boundary& psiBf =
-        this->psi_.boundaryFieldRef();
-
-    volScalarField::Boundary& rhoBf =
-        this->rho_.boundaryFieldRef();
-
-    volScalarField::Boundary& heBf =
-        this->he().boundaryFieldRef();
-
-    volScalarField::Boundary& muBf =
-        this->mu_.boundaryFieldRef();
-
-    volScalarField::Boundary& alphaBf =
-        this->alpha_.boundaryFieldRef();
-
-    forAll(this->T_.boundaryField(), patchi)
     {
-        fvPatchScalarField& pp = pBf[patchi];
+        const fvPatchScalarField& pp = pBf[patchi];
         fvPatchScalarField& pT = TBf[patchi];
         fvPatchScalarField& ppsi = psiBf[patchi];
         fvPatchScalarField& prho = rhoBf[patchi];
@@ -134,7 +148,17 @@ Foam::heRhoThermo<BasicPsiThermo, MixtureType>::heRhoThermo
 :
     heThermo<BasicPsiThermo, MixtureType>(mesh, phaseName)
 {
-    calculate();
+    calculate
+    (
+        this->p_,
+        this->T_,
+        this->he_,
+        this->psi_,
+        this->rho_,
+        this->mu_,
+        this->alpha_,
+        true                    // Create old time fields
+    );
 }
 
 
@@ -150,17 +174,21 @@ Foam::heRhoThermo<BasicPsiThermo, MixtureType>::~heRhoThermo()
 template<class BasicPsiThermo, class MixtureType>
 void Foam::heRhoThermo<BasicPsiThermo, MixtureType>::correct()
 {
-    if (debug)
+    DebugInFunction << endl;
-    {
-        InfoInFunction << endl;
-    }
 
-    calculate();
+    calculate
+    (
+        this->p_,
+        this->T_,
+        this->he_,
+        this->psi_,
+        this->rho_,
+        this->mu_,
+        this->alpha_,
+        false           // No need to update old times
+    );
 
-    if (debug)
+    DebugInFunction << "Finished" << endl;
-    {
-        Info<< "    Finished" << endl;
-    }
 }
 
 

src/thermophysicalModels/basic/rhoThermo/heRhoThermo.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) 2015-2017 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -55,7 +55,17 @@ class heRhoThermo
     // Private Member Functions
 
         //- Calculate the thermo variables
-        void calculate();
+        void calculate
+        (
+            const volScalarField& p,
+            volScalarField& T,
+            volScalarField& he,
+            volScalarField& psi,
+            volScalarField& rho,
+            volScalarField& mu,
+            volScalarField& alpha,
+            const bool doOldTimes
+        );
 
         //- Construct as copy (not implemented)
         heRhoThermo(const heRhoThermo<BasicPsiThermo, MixtureType>&);