Merge branch 'nextRelease' of ssh://git.code.sf.net/p/foam-extend/foam-extend-3.2 into nextRelease

ListOfContributors

@@ -36,14 +36,17 @@ Contents:
     Pierre-Olivier Dallaire
     Carsten Dehning
     Norman Del Puppo
+    Ilaria De Dominicis
     Eugene De Villiers
     Jovani Favero
     Andreas Feymark
     Flavio Galeazzo
     Maria Garcia Camprubi
+    Inno Gatin
     Christoph Goniva
     Chris Greenshields
     Bernhard Gschaider
+    Cesare Guardino
     Andy Heather
     David Hill
     Peter Janas
@@ -51,14 +54,17 @@ Contents:
     Klas Jareteg
     Franjo Juretic
     Aleksandar Karac
+    Robert Keser
     Dennis Kingsley
     Nikola Kornev
     Hannes Kroger
     Michael Leonard
     Niels Linnemann
+    Christian Lucas
     Tommaso Lucchini
     Luca Mangani
     Dubravko Matijasevic
+    Alexey Matveichev
     David McAuliffe
     Sandeep Menon
     Håkan Nilsson
@@ -74,19 +80,15 @@ Contents:
     David Schmidt
     Daniel Schmode
     Hua Shan
+    Vanja Skuric
     Richard Smith (Symscape)
     Hilary Spencer
     Darrin Stephens
     Gavin Tabor
     Tian Tang
     Zeljko Tukovic
-    Niklas Wikstrom
-    Vanja Skuric
-    Alexander Vakhrushev
-    Inno Gatin
-    Alexey Matveichev
-    Vuko Vukcevic
-    Robert Keser
-    Cesare Guardino
-    Ilaria De Dominicis
     Tessa Uroic
+    Alexander Vakhrushev
+    Vuko Vukcevic
+    Niklas Wikstrom
+

README_RealGasBranch

@@ -0,0 +1,63 @@
+Features
+
+
+/*************************************Libraries**************************************************/
+
+1. Included cubic equations of state for pure fluid.
+   --> Redlich Kwong --> class = redlichKwong
+   --> Peng Robinson --> class = pengRobinson
+   --> Soave Redlich Kwong --> class = soaveRedlichKwong
+   --> Aungier Redlich Kwong --> class = aungierRedlichKwong
+
+2. Included cubic equations of state for mixtures.
+   TODO: connect the classes to the reaction models
+   -->class = mixtureRedlichKwong
+   -->class = mixturePengRobinson
+   -->class = mixtureSoaveRedlichKwong
+   -->class = aungierSoaveRedlichKwong
+
+3. Included a new Thermo model
+   --> enthalpy calculated using NASA heat capacity polynomical
+       Added real gas correction in class
+       class = nasaHeatCapacityPolynomial
+
+4. Minor changed in thermo related models so that they can be used with real gas classes
+   --> class = realGasSpecieThermo
+   --> class = sutherland
+   --> class = const
+   --> class = basicMixture
+   --> class = basicPsiThermo
+
+5. Added enthalpy and internal energy thermodynamic class based on basicPsiThermo for real gases
+   --> enthalpy based thermo model, class = realGasHThermo
+   --> internal energy based thermo model, class = realGasEThermo
+
+6. Added high precision water properties (IAPWS97) based on freeSteam (external program, not included)
+   Classes will not be compiled by Allwmake file
+   Classes will compile without freeSteam. Freesteam is loaded at runTime (see Tutorial)
+   --> classes = externalMedia/IAPWS_Waterproperties
+
+/*************************************Solver**************************************************/
+ 
+1. Changed pressure equation of rhoPisoFoam. Orginal pressure equation assumes perfect gas 
+   (linear relationship between pressure and density)
+   --> new solver = realFluidPisoFoam
+
+/*************************************Tutorials**************************************************/
+
+1. pure fluid real gas mixture tutorial
+   --> realFluidPisoFoam/ras/backStep
+
+2. IAPWS97 water properties tutorial
+   freeSteam must be installed for this tutorial
+   additional classes loaded in controlDict
+   --> realFluidPisoFoam/ras/cavity_IAPWS97
+
+
+/*************************************Change log**************************************************/
+
+git commit: "add branch ReadMe file"
+   --> added this file
+
+
+   

applications/solvers/compressible/realFluidPisoFoam/Make/files

@@ -0,0 +1,3 @@
+realFluidPisoFoam.C
+
+EXE = $(FOAM_APPBIN)/realFluidPisoFoam

applications/solvers/compressible/realFluidPisoFoam/Make/options

@@ -0,0 +1,14 @@
+EXE_INC = \
+    -I$(LIB_SRC)/thermophysicalModels/specie/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
+    -I$(LIB_SRC)/turbulenceModels/compressible/turbulenceModel \
+    -I$(LIB_SRC)/finiteVolume/lnInclude
+
+EXE_LIBS = \
+    -lfiniteVolume \
+    -lmeshTools \
+    -lbasicThermophysicalModels \
+    -lspecie \
+    -lcompressibleTurbulenceModel \
+    -lcompressibleRASModels \
+    -lcompressibleLESModels 

applications/solvers/compressible/realFluidPisoFoam/UEqn.H

@@ -0,0 +1,11 @@
+    fvVectorMatrix UEqn
+    (
+        fvm::ddt(rho, U)
+      + fvm::div(phi, U)
+      + turbulence->divDevRhoReff()
+    );
+
+    if (piso.momentumPredictor())
+    {
+        solve(UEqn == -fvc::grad(p));
+    }

applications/solvers/compressible/realFluidPisoFoam/createFields.H

@@ -0,0 +1,59 @@
+    Info<< "Reading thermophysical properties\n" << endl;
+
+    autoPtr<basicPsiThermo> pThermo
+    (
+        basicPsiThermo::New(mesh)
+    );
+    basicPsiThermo& thermo = pThermo();
+
+    volScalarField& p = thermo.p();
+    volScalarField& h = thermo.h();
+    const volScalarField& psi = thermo.psi();
+    const volScalarField& drhodh = thermo.drhodh();
+
+    volScalarField rho
+    (
+        IOobject
+        (
+            "rho",
+            runTime.timeName(),
+            mesh,
+            IOobject::NO_READ,
+            IOobject::AUTO_WRITE
+        ),
+        thermo.rho()
+    );
+
+    Info<< "\nReading field U\n" << endl;
+    volVectorField U
+    (
+        IOobject
+        (
+            "U",
+            runTime.timeName(),
+            mesh,
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh
+    );
+
+#   include "compressibleCreatePhi.H"
+
+    Info<< "Creating turbulence model\n" << endl;
+    autoPtr<compressible::turbulenceModel> turbulence
+    (
+        compressible::turbulenceModel::New
+        (
+            rho,
+            U,
+            phi,
+            thermo
+        )
+    );
+
+    Info<< "Creating field DpDt\n" << endl;
+    volScalarField DpDt =
+        fvc::DDt(surfaceScalarField("phiU", phi/fvc::interpolate(rho)), p);
+
+    mesh.schemesDict().setFluxRequired(p.name());

applications/solvers/compressible/realFluidPisoFoam/hEqn.H

@@ -0,0 +1,12 @@
+{
+    solve
+    (
+        fvm::ddt(rho, h)
+      + fvm::div(phi, h)
+      - fvm::laplacian(turbulence->alphaEff(), h)
+     ==
+        DpDt
+    );
+
+    thermo.correct();
+}

applications/solvers/compressible/realFluidPisoFoam/pEqn.H

@@ -0,0 +1,37 @@
+rho = thermo.rho();
+
+volScalarField rUA = 1.0/UEqn.A();
+U = rUA*UEqn.H();
+
+phi =
+    fvc::interpolate(rho)*
+    (
+        (fvc::interpolate(U) & mesh.Sf())
+      + fvc::ddtPhiCorr(rUA, rho, U, phi)
+    );
+
+while (piso.correctNonOrthogonal())
+{
+    fvScalarMatrix pEqn
+    (
+        psi*fvm::ddt(p)
+      + drhodh*fvc::ddt(h)
+      + fvc::div(phi)
+      - fvm::laplacian(rho*rUA, p)
+    );
+
+    pEqn.solve();
+
+    if (piso.finalNonOrthogonalIter())
+    {
+        phi += pEqn.flux();
+    }
+}
+
+#include "rhoEqn.H"
+#include "compressibleContinuityErrs.H"
+
+U -= rUA*fvc::grad(p);
+U.correctBoundaryConditions();
+
+DpDt = fvc::DDt(surfaceScalarField("phiU", phi/fvc::interpolate(rho)), p);

applications/solvers/compressible/realFluidPisoFoam/realFluidPisoFoam.C

@@ -0,0 +1,100 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright held by original author
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+Application
+    realFluidPisoFoam
+
+Description
+    Transient PISO solver for compressible, laminar or turbulent flow
+    of real fluids e.g. real gases (cubic equations of state)
+
+    Solver cannot be used with perfect gas library
+
+\*---------------------------------------------------------------------------*/
+
+#include "fvCFD.H"
+#include "basicPsiThermo.H"
+#include "turbulenceModel.H"
+#include "pisoControl.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+int main(int argc, char *argv[])
+{
+    #include "setRootCase.H"
+
+    #include "createTime.H"
+    #include "createMesh.H"
+
+    pisoControl piso(mesh);
+
+    #include "createFields.H"
+    #include "initContinuityErrs.H"
+    #include "createTimeControls.H"
+    #include "compressibleCourantNo.H"
+    #include "setInitialDeltaT.H"
+
+
+    // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+    Info<< "\nStarting time loop\n" << endl;
+
+    while (runTime.run())
+    {
+        #include "readTimeControls.H"
+        #include "compressibleCourantNo.H"
+        #include "setDeltaT.H"
+
+        runTime++;
+
+        Info<< "Time = " << runTime.timeName() << nl << endl;
+
+        #include "rhoEqn.H"
+        #include "UEqn.H"
+
+        // --- PISO loop
+        while (piso.correct())
+        {
+            #include "hEqn.H"
+            #include "pEqn.H"
+        }
+
+        turbulence->correct();
+
+        rho = thermo.rho();
+
+        runTime.write();
+
+        Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
+            << "  ClockTime = " << runTime.elapsedClockTime() << " s"
+            << nl << endl;
+    }
+
+    Info<< "End\n" << endl;
+
+    return 0;
+}
+
+
+// ************************************************************************* //

applications/solvers/compressible/steadyUniversalFoam/UEqn.H

@@ -10,9 +10,4 @@
 
     UEqn.relax();
 
-    eqnResidual = solve
-    (
-        UEqn == -fvc::grad(p)
-    ).initialResidual();
-
-    maxResidual = max(eqnResidual, maxResidual);
+    solve(UEqn == -fvc::grad(p));

applications/solvers/compressible/steadyUniversalFoam/convergenceCheck.H

@@ -1,9 +0,0 @@
-// check convergence
-
-if (maxResidual < convergenceCriterion)
-{
-    Info<< "reached convergence criterion: " << convergenceCriterion << endl;
-    runTime.writeAndEnd();
-    Info<< "latestTime = " << runTime.timeName() << endl;
-}
-

applications/solvers/compressible/steadyUniversalFoam/createFields.H

@@ -39,3 +39,5 @@
             thermo
         )
     );
+
+    mesh.schemesDict().setFluxRequired(p.name());

applications/solvers/compressible/steadyUniversalFoam/hEqn.H

@@ -23,8 +23,7 @@
 
     hEqn.relax();
 
-    eqnResidual = hEqn.solve().initialResidual();
-    maxResidual = max(eqnResidual, maxResidual);
+    hEqn.solve();
 
     // Bounding of enthalpy taken out
     thermo.correct();

applications/solvers/compressible/steadyUniversalFoam/initConvergenceCheck.H

@@ -1,7 +0,0 @@
-// initialize values for convergence checks
-
-    scalar eqnResidual = 1, maxResidual = 0;
-    scalar convergenceCriterion = 0;
-
-    pimple.readIfPresent("convergence", convergenceCriterion);
-

applications/solvers/compressible/steadyUniversalFoam/pEqn.H

@@ -7,7 +7,7 @@
     // Needs to be outside of loop since p is changing, but psi and rho are not
     surfaceScalarField rhoReff = rhof - psisf*fvc::interpolate(p);
 
-    for (int corr = 0; corr < nCorr; corr++)
+    while (pimple.correct())
     {
         U = rUA*UEqn.H();
 
@@ -20,7 +20,7 @@
 
         p.storePrevIter();
 
-        for (int nonOrth = 0; nonOrth <= nNonOrthCorr; nonOrth++)
+        while (pimple.correctNonOrthogonal())
         {
             fvScalarMatrix pEqn
             (
@@ -30,16 +30,10 @@
               - fvm::laplacian(rho*rUA, p)
             );
 
-            // Retain the residual from the first pressure solution
-            eqnResidual = pEqn.solve().initialResidual();
-
-            if (corr == 0 && nonOrth == 0)
-            {
-                maxResidual = max(eqnResidual, maxResidual);
-            }
+            pEqn.solve();
 
             // Calculate the flux
-            if (nonOrth == nNonOrthCorr)
+            if (pimple.finalNonOrthogonalIter())
             {
                 phi = phid2 + pEqn.flux();
             }

applications/solvers/compressible/steadyUniversalFoam/steadyUniversalFoam.C

@@ -36,6 +36,7 @@ Author
 #include "basicPsiThermo.H"
 #include "basicRhoThermo.H"
 #include "RASModel.H"
+#include "pimpleControl.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -46,9 +47,11 @@ int main(int argc, char *argv[])
 
 #   include "createTime.H"
 #   include "createMesh.H"
+
+    pimpleControl pimple(mesh);
+
 #   include "createThermo.H"
 #   include "createFields.H"
-#   include "readPIMPLEControls.H"
 #   include "initContinuityErrs.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -59,11 +62,8 @@ int main(int argc, char *argv[])
     {
         Info<< "Time = " << runTime.timeName() << nl << endl;
 
-#       include "readPIMPLEControls.H"
 #       include "readFieldBounds.H"
 
-#       include "initConvergenceCheck.H"
-
 #       include "UEqn.H"
 #       include "pEqn.H"
 
@@ -79,8 +79,6 @@ int main(int argc, char *argv[])
         Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
             << "  ClockTime = " << runTime.elapsedClockTime() << " s"
             << nl << endl;
-
-#       include "convergenceCheck.H"
     }
 
 #   include "clearThermo.H"

applications/solvers/compressible/steadyUniversalMRFFoam/UEqn.H

@@ -13,9 +13,4 @@
 
     UEqn.relax();
 
-    eqnResidual = solve
-    (
-        UEqn == -fvc::grad(p)
-    ).initialResidual();
-
-    maxResidual = max(eqnResidual, maxResidual);
+    solve(UEqn == -fvc::grad(p));

applications/solvers/compressible/steadyUniversalMRFFoam/convergenceCheck.H

@@ -1,9 +0,0 @@
-// check convergence
-
-if (maxResidual < convergenceCriterion)
-{
-    Info<< "reached convergence criterion: " << convergenceCriterion << endl;
-    runTime.writeAndEnd();
-    Info<< "latestTime = " << runTime.timeName() << endl;
-}
-

applications/solvers/compressible/steadyUniversalMRFFoam/createFields.H

@@ -80,3 +80,5 @@
         mesh 
     );
     i == h - 0.5*(magSqr(Urot) - magSqr(Urel));
+
+    mesh.schemesDict().setFluxRequired(p.name());

applications/solvers/compressible/steadyUniversalMRFFoam/iEqn.H

@@ -7,7 +7,7 @@
     Urot == U - Urel;
 
     fvScalarMatrix iEqn
-    (      
+    (
 	      fvm::ddt(rho, i)
         + fvm::div(phi, i)
         - fvm::laplacian(turbulence->alphaEff(), i)
@@ -18,8 +18,7 @@
 
     iEqn.relax();
 
-    eqnResidual = iEqn.solve().initialResidual();
-    maxResidual = max(eqnResidual, maxResidual);
+    iEqn.solve();
 
     // From rothalpy, calculate enthalpy after solution of rothalpy equation
     h = i + 0.5*(magSqr(Urot) - magSqr(Urel));

applications/solvers/compressible/steadyUniversalMRFFoam/initConvergenceCheck.H

@@ -1,7 +0,0 @@
-// initialize values for convergence checks
-
-    scalar eqnResidual = 1, maxResidual = 0;
-    scalar convergenceCriterion = 0;
-
-    pimple.readIfPresent("convergence", convergenceCriterion);
-

applications/solvers/compressible/steadyUniversalMRFFoam/pEqn.H

@@ -7,7 +7,7 @@
     // Needs to be outside of loop since p is changing, but psi and rho are not
     surfaceScalarField rhoReff = rhof - psisf*fvc::interpolate(p);
 
-    for (int corr = 0; corr < nCorr; corr++)
+    while (pimple.correct())
     {
         U = rUA*UEqn.H();
 
@@ -20,7 +20,7 @@
 
         p.storePrevIter();
 
-        for (int nonOrth = 0; nonOrth <= nNonOrthCorr; nonOrth++)
+        while (pimple.correctNonOrthogonal())
         {
             fvScalarMatrix pEqn
             (
@@ -30,16 +30,10 @@
               - fvm::laplacian(rho*rUA, p)
             );
 
-            // Retain the residual from the first pressure solution
-            eqnResidual = pEqn.solve().initialResidual();
-
-            if (corr == 0 && nonOrth == 0)
-            {
-                maxResidual = max(eqnResidual, maxResidual);
-            }
+            pEqn.solve();
 
             // Calculate the flux
-            if (nonOrth == nNonOrthCorr)
+            if (pimple.finalNonOrthogonalIter())
             {
                 phi = phid2 + pEqn.flux();
             }

applications/solvers/compressible/steadyUniversalMRFFoam/steadyUniversalMRFFoam.C

@@ -43,6 +43,7 @@ GE CONFIDENTIAL INFORMATION 2016 General Electric Company. All Rights Reserved
 #include "basicRhoThermo.H"
 #include "RASModel.H"
 #include "MRFZones.H"
+#include "pimpleControl.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -53,9 +54,11 @@ int main(int argc, char *argv[])
 
 #   include "createTime.H"
 #   include "createMesh.H"
+
+    pimpleControl pimple(mesh);
+
 #   include "createThermo.H"
 #   include "createFields.H"
-#   include "readPIMPLEControls.H"
 #   include "initContinuityErrs.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -66,11 +69,8 @@ int main(int argc, char *argv[])
     {
         Info<< "Time = " << runTime.timeName() << nl << endl;
 
-#       include "readPIMPLEControls.H"
 #       include "readFieldBounds.H"
 
-#       include "initConvergenceCheck.H"
-
 #       include "UEqn.H"
 #       include "pEqn.H"
 
@@ -86,8 +86,6 @@ int main(int argc, char *argv[])
         Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
             << "  ClockTime = " << runTime.elapsedClockTime() << " s"
             << nl << endl;
-
-#       include "convergenceCheck.H"
     }
 
 #   include "clearThermo.H"

applications/solvers/finiteArea/liquidFilmFoam/createFaFields.H

@@ -154,3 +154,5 @@
         aMesh,
         dimensionedScalar("one", dimless, 0.01)
     );
+
+    aMesh.schemesDict().setFluxRequired("h");

applications/solvers/multiphase/multiphaseInterFoam/multiphaseMixture/multiphaseMixture.C

@@ -125,6 +125,7 @@ Foam::multiphaseMixture::multiphaseMixture
     forAllIter(PtrDictionary<phase>, phases_, iter)
     {
         alphaTable_.add(iter());
+        mesh_.schemesDict().setFluxRequired(iter().volScalarField::name());
     }
 }
 

applications/solvers/solidMechanics/elasticAcpSolidFoam/elasticAcpSolidFoam.C

@@ -102,7 +102,7 @@ int main(int argc, char *argv[])
         //#     include "waveCourantNo.H"
 
         int iCorr = 0;
-        lduMatrix::solverPerformance solverPerf;
+        lduSolverPerformance solverPerf;
         scalar initialResidual = 0;
         scalar relativeResidual = 1;
         //scalar forceResidual = 1;

applications/solvers/solidMechanics/elasticIncrAcpSolidFoam/elasticIncrAcpSolidFoam.C

@@ -86,7 +86,7 @@ int main(int argc, char *argv[])
         //#     include "waveCourantNo.H"
 
         int iCorr = 0;
-        lduMatrix::solverPerformance solverPerf;
+        lduSolverPerformance solverPerf;
         scalar initialResidual = 0;
         scalar relativeResidual = 1;
         //scalar forceResidual = 1;

applications/solvers/solidMechanics/elasticIncrSolidFoam/elasticIncrSolidFoam.C

@@ -70,7 +70,7 @@ int main(int argc, char *argv[])
         int iCorr = 0;
         scalar initialResidual = 0;
         scalar relativeResidual = 1.0;
-        lduMatrix::solverPerformance solverPerf;
+        lduSolverPerformance solverPerf;
         lduMatrix::debug = 0;
 
         do

applications/solvers/solidMechanics/elasticNonLinIncrTLSolidFoam/elasticNonLinIncrTLSolidFoam.C

@@ -60,7 +60,7 @@ int main(int argc, char *argv[])
 
         int iCorr = 0;
         scalar initialResidual = 0;
-        lduMatrix::solverPerformance solverPerf;
+        lduSolverPerformance solverPerf;
         scalar relativeResidual = 1.0;
         lduMatrix::debug=0;
 

applications/solvers/solidMechanics/elasticNonLinTLSolidFoam/elasticNonLinTLSolidFoam.C

@@ -62,7 +62,7 @@ int main(int argc, char *argv[])
 
         int iCorr = 0;
         scalar initialResidual = 0;
-        lduMatrix::solverPerformance solverPerf;
+        lduSolverPerformance solverPerf;
         scalar relativeResidual = 1.0;
         lduMatrix::debug = 0;
 

applications/solvers/solidMechanics/elasticNonLinULSolidFoam/elasticNonLinULSolidFoam.C

@@ -73,7 +73,7 @@ int main(int argc, char *argv[])
 #       include "readSolidMechanicsControls.H"
 
         int iCorr = 0;
-        lduMatrix::solverPerformance solverPerf;
+        lduSolverPerformance solverPerf;
         scalar initialResidual = 1.0;
         scalar relativeResidual = 1.0;
         lduMatrix::debug = 0;

applications/solvers/solidMechanics/elasticOrthoAcpSolidFoam/elasticOrthoAcpSolidFoam.C

@@ -103,7 +103,7 @@ int main(int argc, char *argv[])
         //#     include "waveCourantNo.H"
 
         int iCorr = 0;
-        lduMatrix::solverPerformance solverPerf;
+        lduSolverPerformance solverPerf;
         scalar initialResidual = 0;
         scalar relativeResidual = 1;
         //scalar forceResidual = 1;

applications/solvers/solidMechanics/elasticOrthoNonLinULSolidFoam/elasticOrthoNonLinULSolidFoam.C

@@ -78,7 +78,7 @@ int main(int argc, char *argv[])
 #       include "readSolidMechanicsControls.H"
 
         int iCorr = 0;
-        lduMatrix::solverPerformance solverPerf;
+        lduSolverPerformance solverPerf;
         scalar initialResidual = 1.0;
         lduMatrix::debug = 0;
 

applications/solvers/solidMechanics/elasticOrthoSolidFoam/elasticOrthoSolidFoam.C

@@ -67,7 +67,7 @@ int main(int argc, char *argv[])
 #       include "readSolidMechanicsControls.H"
 
         int iCorr = 0;
-        lduMatrix::solverPerformance solverPerf;
+        lduSolverPerformance solverPerf;
         scalar initialResidual = 1.0;
         scalar relativeResidual = 1.0;
         lduMatrix::debug = 0;

applications/solvers/solidMechanics/elasticPlasticNonLinTLSolidFoam/elasticPlasticNonLinTLSolidFoam.C

@@ -64,7 +64,7 @@ int main(int argc, char *argv[])
 
         int iCorr = 0;
         scalar initialResidual = 0;
-        lduMatrix::solverPerformance solverPerf;
+        lduSolverPerformance solverPerf;
         scalar relativeResidual = GREAT;
         lduMatrix::debug = 0;
 

applications/solvers/solidMechanics/elasticPlasticNonLinULSolidFoam/elasticPlasticNonLinULSolidFoam.C

@@ -80,7 +80,7 @@ int main(int argc, char *argv[])
 #       include "readSolidMechanicsControls.H"
 
         int iCorr = 0;
-        lduMatrix::solverPerformance solverPerf;
+        lduSolverPerformance solverPerf;
         scalar initialResidual = 0;
         scalar relativeResidual = 1.0;
         lduMatrix::debug = 0;

applications/solvers/solidMechanics/elasticPlasticSolidFoam/elasticPlasticSolidFoam.C

@@ -65,7 +65,7 @@ int main(int argc, char *argv[])
 #       include "readSolidMechanicsControls.H"
 
         int iCorr = 0;
-        lduMatrix::solverPerformance solverPerf;
+        lduSolverPerformance solverPerf;
         scalar initialResidual = 1.0;
         scalar relativeResidual = 1.0;
         scalar plasticResidual = 1.0;

applications/solvers/solidMechanics/elasticSolidFoam/elasticSolidFoam.C

@@ -69,7 +69,7 @@ int main(int argc, char *argv[])
 #       include "readSolidMechanicsControls.H"
 
         int iCorr = 0;
-        lduMatrix::solverPerformance solverPerf;
+        lduSolverPerformance solverPerf;
         scalar initialResidual = 1.0;
         scalar relativeResidual = 1.0;
 //         lduMatrix::debug = 0;

applications/solvers/solidMechanics/elasticThermalSolidFoam/elasticThermalSolidFoam.C

@@ -66,8 +66,8 @@ int main(int argc, char *argv[])
         scalar initialResidual = 1.0;
         scalar relResT = 1.0;
         scalar relResU = 1.0;
-        lduMatrix::solverPerformance solverPerfU;
-        lduMatrix::solverPerformance solverPerfT;
+        lduSolverPerformance solverPerfU;
+        lduSolverPerformance solverPerfT;
         lduMatrix::debug = 0;
 
         // solve energy equation for temperature

applications/solvers/solidMechanics/icoFsiElasticNonLinULSolidFoam/solveSolidBackward.H

@@ -2,7 +2,7 @@
 #       include "readSolidMechanicsControls.H"
 
         int iCorr = 0;
-        lduMatrix::solverPerformance solverPerf;
+        lduSolverPerformance solverPerf;
         scalar initialResidual = 0;
 
         lduMatrix::debug = 0;

applications/solvers/solidMechanics/icoFsiElasticNonLinULSolidFoam/solveSolidEuler.H

@@ -2,7 +2,7 @@
 #       include "readSolidMechanicsControls.H"
 
         int iCorr = 0;
-        lduMatrix::solverPerformance solverPerf;
+        lduSolverPerformance solverPerf;
         scalar initialResidual = 0;
 
         lduMatrix::debug = 0;

applications/solvers/solidMechanics/viscoElasticSolidFoam/calcInitialState.H

@@ -11,7 +11,7 @@
         Info << "lambda = " << average(lambda.internalField()) << endl;
 
         int iCorr = 0;
-        lduMatrix::solverPerformance solverPerf;
+        lduSolverPerformance solverPerf;
         scalar initialResidual = 0;
         scalar err = GREAT;
 

applications/solvers/solidMechanics/viscoElasticSolidFoam/viscoElasticSolidFoam.C

@@ -73,7 +73,7 @@ int main(int argc, char *argv[])
         Info << "average lambda = " << average(lambdaf.internalField()) << endl;
 
         int iCorr = 0;
-        lduMatrix::solverPerformance solverPerf;
+        lduSolverPerformance solverPerf;
         scalar initialResidual = 1.0;
         scalar residual = 1.0;
         surfaceSymmTensorField DSigmaCorrf = fvc::interpolate(DSigmaCorr);

applications/utilities/parallelProcessing/decomposePar/faMeshDecomposition.C

@@ -61,27 +61,87 @@ void faMeshDecomposition::distributeFaces()
             )
         );
 
-        labelHashSet faceProcAddressingHash
-        (
-            labelIOList
-            (
-                IOobject
-                (
-                    "faceProcAddressing",
-                    "constant",
-                    procMesh.meshSubDir,
-                    procMesh,
-                    IOobject::MUST_READ,
-                    IOobject::NO_WRITE
-                )
-            )
-        );
-
-        forAll (faceLabels(), faceI)
+        // If faMesh's fvPatch is a part of the global face zones, faces of that
+        // patch will be present on all processors. Because of that, looping
+        // through faceProcAddressing will decompose global faMesh faces to the
+        // very last processor regardless of where fvPatch is really decomposed.
+        // Since global faces which do not belong to specific processor are
+        // located at the end of the faceProcAddressing, cutting it at
+        // i = owner.size() will correctly decompose faMesh faces.
+        // Vanja Skuric, 2016-04-21
+        if (decompositionDict_.found("globalFaceZones"))
         {
-            if (faceProcAddressingHash.found(faceLabels()[faceI] + 1))
+            labelList faceProcAddressing
+            (
+                labelIOList
+                (
+                    IOobject
+                    (
+                        "faceProcAddressing",
+                        "constant",
+                        procMesh.meshSubDir,
+                        procMesh,
+                        IOobject::MUST_READ,
+                        IOobject::NO_WRITE
+                    )
+                )
+            );
+
+            const label ownerSize =
+            (
+                labelIOList
+                (
+                    IOobject
+                    (
+                        "owner",
+                        "constant",
+                        procMesh.meshSubDir,
+                        procMesh,
+                        IOobject::MUST_READ,
+                        IOobject::NO_WRITE
+                    )
+                )
+            ).size();
+
+            labelHashSet faceProcAddressingHash(ownerSize);
+
+            for (int i = 0; i < ownerSize; ++i)
             {
-                faceToProc_[faceI] = procI;
+                faceProcAddressingHash.insert(faceProcAddressing[i]);
+            }
+
+            forAll (faceLabels(), faceI)
+            {
+                if (faceProcAddressingHash.found(faceLabels()[faceI] + 1))
+                {
+                    faceToProc_[faceI] = procI;
+                }
+            }
+        }
+        else
+        {
+            labelHashSet faceProcAddressingHash
+            (
+                labelIOList
+                (
+                    IOobject
+                    (
+                        "faceProcAddressing",
+                        "constant",
+                        procMesh.meshSubDir,
+                        procMesh,
+                        IOobject::MUST_READ,
+                        IOobject::NO_WRITE
+                    )
+                )
+            );
+
+            forAll (faceLabels(), faceI)
+            {
+                if (faceProcAddressingHash.found(faceLabels()[faceI] + 1))
+                {
+                    faceToProc_[faceI] = procI;
+                }
             }
         }
     }

src/dynamicMesh/meshMotion/tetMotionSolver/tetMotionSolver/laplace/laplaceTetMotionSolver.H

@@ -93,7 +93,7 @@ protected:
         }
 
         // Holds recent solver performance
-        lduMatrix::solverPerformance solverPerf_;
+        lduSolverPerformance solverPerf_;
 
 
 public:
@@ -125,7 +125,7 @@ public:
         virtual void updateMesh(const mapPolyMesh&);
 
         //- Return recent solver performance
-        const lduMatrix::solverPerformance& solverPerformance() const
+        const lduSolverPerformance& solverPerformance() const
         {
             return solverPerf_;
         }

src/finiteArea/finiteArea/faSchemes/faSchemes.C

@@ -21,31 +21,45 @@ License
     You should have received a copy of the GNU General Public License
     along with foam-extend.  If not, see <http://www.gnu.org/licenses/>.
 
-Description
-
 \*---------------------------------------------------------------------------*/
 
-#include "error.H"
-#include "objectRegistry.H"
-#include "foamTime.H"
 #include "faSchemes.H"
+#include "objectRegistry.H"
 
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-namespace Foam
-{
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
 Foam::debug::debugSwitch
-faSchemes::debug
+Foam::faSchemes::debug
 (
     "faSchemes",
     false
 );
 
+// * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * * //
+
+void Foam::faSchemes::clear()
+{
+    ddtSchemes_.clear();
+    defaultDdtScheme_.clear();
+    d2dt2Schemes_.clear();
+    defaultD2dt2Scheme_.clear();
+    interpolationSchemes_.clear();
+    defaultInterpolationScheme_.clear();
+    divSchemes_.clear(); // optional
+    defaultDivScheme_.clear();
+    gradSchemes_.clear(); // optional
+    defaultGradScheme_.clear();
+    lnGradSchemes_.clear();
+    defaultLnGradScheme_.clear();
+    laplacianSchemes_.clear(); // optional
+    defaultLaplacianScheme_.clear();
+    fluxRequired_.clear();
+    defaultFluxRequired_ = false;
+}
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
-faSchemes::faSchemes(const objectRegistry& obr)
+Foam::faSchemes::faSchemes(const objectRegistry& obr)
 :
     IOdictionary
     (
@@ -59,21 +73,106 @@ faSchemes::faSchemes(const objectRegistry& obr)
             IOobject::NO_WRITE
         )
     ),
-    ddtSchemes_(ITstream("ddtSchemes", tokenList())()),
-    defaultDdtScheme_("default", tokenList()),
-    d2dt2Schemes_(ITstream("d2dt2Schemes", tokenList())()),
-    defaultD2dt2Scheme_("default", tokenList()),
-    interpolationSchemes_(ITstream("interpolationSchemes", tokenList())()),
-    defaultInterpolationScheme_("default", tokenList()),
-    divSchemes_(ITstream("divSchemes", tokenList())()),
-    defaultDivScheme_("default", tokenList()),
-    gradSchemes_(ITstream("gradSchemes", tokenList())()),
-    defaultGradScheme_("default", tokenList()),
-    lnGradSchemes_(ITstream("lnGradSchemes", tokenList())()),
-    defaultLnGradScheme_("default", tokenList()),
-    laplacianSchemes_(ITstream("laplacianSchemes", tokenList())()),
-    defaultLaplacianScheme_("default", tokenList()),
-    fluxRequired_(ITstream("fluxRequired", tokenList())())
+    ddtSchemes_
+    (
+        ITstream
+        (
+            objectPath() + "::ddtSchemes",
+            tokenList()
+        )()
+    ),
+    defaultDdtScheme_
+    (
+        ddtSchemes_.name() + "::default",
+        tokenList()
+    ),
+    d2dt2Schemes_
+    (
+        ITstream
+        (
+            objectPath() + "::d2dt2Schemes",
+            tokenList()
+        )()
+    ),
+    defaultD2dt2Scheme_
+    (
+        d2dt2Schemes_.name() + "::default",
+        tokenList()
+    ),
+    interpolationSchemes_
+    (
+        ITstream
+        (
+            objectPath() + "::interpolationSchemes",
+            tokenList()
+        )()
+    ),
+    defaultInterpolationScheme_
+    (
+        interpolationSchemes_.name() + "::default",
+        tokenList()
+    ),
+    divSchemes_
+    (
+        ITstream
+        (
+            objectPath() + "::divSchemes",
+            tokenList()
+        )()
+    ),
+    defaultDivScheme_
+    (
+        divSchemes_.name() + "::default",
+        tokenList()
+    ),
+    gradSchemes_
+    (
+        ITstream
+        (
+            objectPath() + "::gradSchemes",
+            tokenList()
+        )()
+    ),
+    defaultGradScheme_
+    (
+        gradSchemes_.name() + "::default",
+        tokenList()
+    ),
+    lnGradSchemes_
+    (
+        ITstream
+        (
+            objectPath() + "::snGradSchemes",
+            tokenList()
+        )()
+    ),
+    defaultLnGradScheme_
+    (
+        lnGradSchemes_.name() + "::default",
+        tokenList()
+    ),
+    laplacianSchemes_
+    (
+        ITstream
+        (
+            objectPath() + "::laplacianSchemes",
+            tokenList()
+        )()
+    ),
+    defaultLaplacianScheme_
+    (
+        laplacianSchemes_.name() + "::default",
+        tokenList()
+    ),
+    fluxRequired_
+    (
+        ITstream
+        (
+            objectPath() + "::fluxRequired",
+            tokenList()
+        )()
+    ),
+    defaultFluxRequired_(false)
 {
     if (!headerOk())
     {
@@ -85,6 +184,8 @@ faSchemes::faSchemes(const objectRegistry& obr)
             )   << "faSchemes dictionary not found.  Creating default."
                 << endl;
         }
+
+        regIOobject::write();
     }
 
     read();
@@ -93,27 +194,58 @@ faSchemes::faSchemes(const objectRegistry& obr)
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-bool faSchemes::read()
+bool Foam::faSchemes::read()
 {
-    bool readOk = false;
-
-    if (headerOk())
-    {
-        readOk = regIOobject::read();
-    }
-
-    if (readOk)
+    if (regIOobject::read())
     {
         const dictionary& dict = schemesDict();
 
-        // ddt Schemes
+        // persistent settings across reads is incorrect
+        clear();
+
         if (dict.found("ddtSchemes"))
         {
             ddtSchemes_ = dict.subDict("ddtSchemes");
         }
+        else if (dict.found("timeScheme"))
+        {
+            // For backward compatibility.
+            // The timeScheme will be deprecated with warning or removed
+            WarningIn("fvSchemes::read()")
+                << "Using deprecated 'timeScheme' instead of 'ddtSchemes'"
+                << nl << endl;
+
+            word schemeName(dict.lookup("timeScheme"));
+
+            if (schemeName == "EulerImplicit")
+            {
+                schemeName = "Euler";
+            }
+            else if (schemeName == "BackwardDifferencing")
+            {
+                schemeName = "backward";
+            }
+            else if (schemeName == "SteadyState")
+            {
+                schemeName = "steadyState";
+            }
+            else
+            {
+                FatalIOErrorIn("fvSchemes::read()", dict.lookup("timeScheme"))
+                    << "\n    Only EulerImplicit, BackwardDifferencing and "
+                       "SteadyState\n    are supported by the old timeScheme "
+                       "specification.\n    Please use ddtSchemes instead."
+                    << exit(FatalIOError);
+            }
+
+            ddtSchemes_.set("default", schemeName);
+
+            ddtSchemes_.lookup("default")[0].lineNumber() =
+                dict.lookup("timeScheme").lineNumber();
+        }
         else
         {
-            ddtSchemes_.add("default", "none");
+            ddtSchemes_.set("default", "none");
         }
 
         if
@@ -126,7 +258,6 @@ bool faSchemes::read()
         }
 
 
-        // d2dt2 schemes
         if (dict.found("d2dt2Schemes"))
         {
             d2dt2Schemes_ = dict.subDict("d2dt2Schemes");
@@ -134,29 +265,30 @@ bool faSchemes::read()
         else if (dict.found("timeScheme"))
         {
             // For backward compatibility.
-            // The timeScheme will be deprecated with warning or removed in 2.4.
+            // The timeScheme will be deprecated with warning or removed
+            WarningIn("fvSchemes::read()")
+                << "Using deprecated 'timeScheme' instead of 'd2dt2Schemes'"
+                << nl << endl;
 
-            word timeSchemeName(dict.lookup("timeScheme"));
+            word schemeName(dict.lookup("timeScheme"));
 
-            if (timeSchemeName == "EulerImplicit")
+            if (schemeName == "EulerImplicit")
             {
-                timeSchemeName = "Euler";
+                schemeName = "Euler";
             }
-            else if (timeSchemeName == "SteadyState")
+            else if (schemeName == "SteadyState")
             {
-                timeSchemeName = "steadyState";
+                schemeName = "steadyState";
             }
 
-            if (d2dt2Schemes_.found("default"))
-            {
-                d2dt2Schemes_.remove("default");
-            }
+            d2dt2Schemes_.set("default", schemeName);
 
-            d2dt2Schemes_.add("default", timeSchemeName);
+            d2dt2Schemes_.lookup("default")[0].lineNumber() =
+                dict.lookup("timeScheme").lineNumber();
         }
         else
         {
-            d2dt2Schemes_.add("default", "none");
+            d2dt2Schemes_.set("default", "none");
         }
 
         if
@@ -169,7 +301,6 @@ bool faSchemes::read()
         }
 
 
-        // Interpolation schemes
         if (dict.found("interpolationSchemes"))
         {
             interpolationSchemes_ = dict.subDict("interpolationSchemes");
@@ -190,15 +321,9 @@ bool faSchemes::read()
         }
 
 
-        // Div schemes
         if (dict.found("divSchemes"))
         {
             divSchemes_ = dict.subDict("divSchemes");
-
-        }
-        else
-        {
-            divSchemes_.add("default", "none");
         }
 
         if
@@ -210,16 +335,9 @@ bool faSchemes::read()
             defaultDivScheme_ = divSchemes_.lookup("default");
         }
 
-
-        // Grad schemes
         if (dict.found("gradSchemes"))
         {
             gradSchemes_ = dict.subDict("gradSchemes");
-
-        }
-        else
-        {
-            gradSchemes_.add("default", "none");
         }
 
         if
@@ -232,12 +350,11 @@ bool faSchemes::read()
         }
 
 
-        // lnGrad schemes
         if (dict.found("lnGradSchemes"))
         {
             lnGradSchemes_ = dict.subDict("lnGradSchemes");
         }
-        else
+        else if (!lnGradSchemes_.found("default"))
         {
             lnGradSchemes_.add("default", "corrected");
         }
@@ -252,15 +369,9 @@ bool faSchemes::read()
         }
 
 
-        // laplacian schemes
         if (dict.found("laplacianSchemes"))
         {
             laplacianSchemes_ = dict.subDict("laplacianSchemes");
-
-        }
-        else
-        {
-            laplacianSchemes_.add("default", "none");
         }
 
         if
@@ -273,18 +384,30 @@ bool faSchemes::read()
         }
 
 
-        // Flux required
         if (dict.found("fluxRequired"))
         {
             fluxRequired_ = dict.subDict("fluxRequired");
-        }
-    }
 
-    return readOk;
+            if
+            (
+                fluxRequired_.found("default")
+             && word(fluxRequired_.lookup("default")) != "none"
+            )
+            {
+                defaultFluxRequired_ = Switch(fluxRequired_.lookup("default"));
+            }
+        }
+
+        return true;
+    }
+    else
+    {
+        return false;
+    }
 }
 
 
-const dictionary& faSchemes::schemesDict() const
+const Foam::dictionary& Foam::faSchemes::schemesDict() const
 {
     if (found("select"))
     {
@@ -297,18 +420,14 @@ const dictionary& faSchemes::schemesDict() const
 }
 
 
-ITstream& faSchemes::ddtScheme(const word& name) const
+Foam::ITstream& Foam::faSchemes::ddtScheme(const word& name) const
 {
     if (debug)
     {
         Info<< "Lookup ddtScheme for " << name << endl;
     }
 
-    if
-    (
-        ddtSchemes_.found(name)
-     || !defaultDdtScheme_.size()
-    )
+    if (ddtSchemes_.found(name) || defaultDdtScheme_.empty())
     {
         return ddtSchemes_.lookup(name);
     }
@@ -320,18 +439,14 @@ ITstream& faSchemes::ddtScheme(const word& name) const
 }
 
 
-ITstream& faSchemes::d2dt2Scheme(const word& name) const
+Foam::ITstream& Foam::faSchemes::d2dt2Scheme(const word& name) const
 {
     if (debug)
     {
         Info<< "Lookup d2dt2Scheme for " << name << endl;
     }
 
-    if
-    (
-        d2dt2Schemes_.found(name)
-     || !defaultD2dt2Scheme_.size()
-    )
+    if (d2dt2Schemes_.found(name) || defaultD2dt2Scheme_.empty())
     {
         return d2dt2Schemes_.lookup(name);
     }
@@ -343,12 +458,17 @@ ITstream& faSchemes::d2dt2Scheme(const word& name) const
 }
 
 
-ITstream& faSchemes::interpolationScheme(const word& name) const
+Foam::ITstream& Foam::faSchemes::interpolationScheme(const word& name) const
 {
+    if (debug)
+    {
+        Info<< "Lookup interpolationScheme for " << name << endl;
+    }
+
     if
     (
         interpolationSchemes_.found(name)
-     || !defaultInterpolationScheme_.size()
+     || defaultInterpolationScheme_.empty()
     )
     {
         return interpolationSchemes_.lookup(name);
@@ -361,9 +481,14 @@ ITstream& faSchemes::interpolationScheme(const word& name) const
 }
 
 
-ITstream& faSchemes::divScheme(const word& name) const
+Foam::ITstream& Foam::faSchemes::divScheme(const word& name) const
 {
-    if (divSchemes_.found(name) || !defaultDivScheme_.size())
+    if (debug)
+    {
+        Info<< "Lookup divScheme for " << name << endl;
+    }
+
+    if (divSchemes_.found(name) || defaultDivScheme_.empty())
     {
         return divSchemes_.lookup(name);
     }
@@ -375,9 +500,14 @@ ITstream& faSchemes::divScheme(const word& name) const
 }
 
 
-ITstream& faSchemes::gradScheme(const word& name) const
+Foam::ITstream& Foam::faSchemes::gradScheme(const word& name) const
 {
-    if (gradSchemes_.found(name) || !defaultGradScheme_.size())
+    if (debug)
+    {
+        Info<< "Lookup gradScheme for " << name << endl;
+    }
+
+    if (gradSchemes_.found(name) || defaultGradScheme_.empty())
     {
         return gradSchemes_.lookup(name);
     }
@@ -389,9 +519,14 @@ ITstream& faSchemes::gradScheme(const word& name) const
 }
 
 
-ITstream& faSchemes::lnGradScheme(const word& name) const
+Foam::ITstream& Foam::faSchemes::lnGradScheme(const word& name) const
 {
-    if (lnGradSchemes_.found(name) || !defaultLnGradScheme_.size())
+    if (debug)
+    {
+        Info<< "Lookup snGradScheme for " << name << endl;
+    }
+
+    if (lnGradSchemes_.found(name) || defaultLnGradScheme_.empty())
     {
         return lnGradSchemes_.lookup(name);
     }
@@ -403,9 +538,14 @@ ITstream& faSchemes::lnGradScheme(const word& name) const
 }
 
 
-ITstream& faSchemes::laplacianScheme(const word& name) const
+Foam::ITstream& Foam::faSchemes::laplacianScheme(const word& name) const
 {
-    if (laplacianSchemes_.found(name) || !defaultLaplacianScheme_.size())
+    if (debug)
+    {
+        Info<< "Lookup laplacianScheme for " << name << endl;
+    }
+
+    if (laplacianSchemes_.found(name) || defaultLaplacianScheme_.empty())
     {
         return laplacianSchemes_.lookup(name);
     }
@@ -417,13 +557,24 @@ ITstream& faSchemes::laplacianScheme(const word& name) const
 }
 
 
-bool faSchemes::fluxRequired(const word& name) const
+void Foam::faSchemes::setFluxRequired(const word& name) const
+{
+    if (debug)
+    {
+        Info<< "Setting fluxRequired for " << name << endl;
+    }
+
+    fluxRequired_.add(name, true, true);
+}
+
+
+bool Foam::faSchemes::fluxRequired(const word& name) const
 {
     return fluxRequired_.found(name);
 }
 
 
-bool faSchemes::writeData(Ostream& os) const
+bool Foam::faSchemes::writeData(Ostream& os) const
 {
     // Write dictionaries
     os << nl << "ddtSchemes";
@@ -454,8 +605,4 @@ bool faSchemes::writeData(Ostream& os) const
 }
 
 
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-} // End namespace Foam
-
 // ************************************************************************* //

src/finiteArea/finiteArea/faSchemes/faSchemes.H

@@ -53,8 +53,6 @@ class faSchemes
 :
     public IOdictionary
 {
-private:
-
     // Private data
 
         dictionary ddtSchemes_;
@@ -78,13 +76,19 @@ private:
         dictionary laplacianSchemes_;
         ITstream defaultLaplacianScheme_;
 
-        dictionary fluxRequired_;
+        mutable dictionary fluxRequired_;
+        bool defaultFluxRequired_;
 
 
     // Private Member Functions
 
-        //- Disallow default bitwise copy construct and assignment
+        //- Clear the dictionaries and streams before reading
+        void clear();
+
+        //- Disallow default bitwise copy construct
         faSchemes(const faSchemes&);
+
+        //- Disallow default bitwise assignment
         void operator=(const faSchemes&);
 
 
@@ -120,6 +124,8 @@ public:
 
             ITstream& laplacianScheme(const word& name) const;
 
+            void setFluxRequired(const word& name) const;
+
             bool fluxRequired(const word& name) const;
 
 

src/finiteArea/finiteArea/fam/famLaplacian.C

@@ -329,6 +329,40 @@ laplacian
 }
 
 
+template<class Type>
+tmp<faMatrix<Type> >
+laplacian
+(
+    const edgeTensorField& gamma,
+    const GeometricField<Type, faPatchField, areaMesh>& vf,
+    const word& name
+)
+{
+    const faMesh& mesh = vf.mesh();
+
+    return fam::laplacian
+    (
+        (mesh.Le() & gamma & mesh.Le())/sqr(mesh.magLe()),
+        vf,
+        name
+    );
+}
+
+template<class Type>
+tmp<faMatrix<Type> >
+laplacian
+(
+    const tmp<edgeTensorField>& tgamma,
+    const GeometricField<Type, faPatchField, areaMesh>& vf,
+    const word& name
+)
+{
+    tmp<faMatrix<Type> > Laplacian = fam::laplacian(tgamma(), vf, name);
+    tgamma.clear();
+    return Laplacian;
+}
+
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 } // End namespace fam

src/finiteArea/finiteArea/fam/famLaplacian.H

@@ -180,6 +180,23 @@ namespace fam
         const tmp<edgeTensorField>&,
         const GeometricField<Type, faPatchField, areaMesh>&
     );
+
+
+    template<class Type>
+    tmp<faMatrix<Type> > laplacian
+    (
+        const edgeTensorField&,
+        const GeometricField<Type, faPatchField, areaMesh>&,
+        const word& name
+    );
+
+    template<class Type>
+    tmp<faMatrix<Type> > laplacian
+    (
+        const tmp<edgeTensorField>&,
+        const GeometricField<Type, faPatchField, areaMesh>&,
+        const word& name
+    );
 }
 
 

src/finiteVolume/cfdTools/general/solutionControl/pimpleControl/pimpleControl.C

@@ -63,7 +63,7 @@ bool Foam::pimpleControl::criteriaSatisfied()
     bool achieved = true;
     bool checked = false;    // safety that some checks were indeed performed
 
-    const dictionary& solverDict = mesh_.solverPerformanceDict();
+    const dictionary& solverDict = mesh_.solutionDict().solverPerformanceDict();
     forAllConstIter(dictionary, solverDict, iter)
     {
         const word& variableName = iter().keyword();
@@ -190,7 +190,7 @@ bool Foam::pimpleControl::loop()
         }
 
         corr_ = 0;
-        mesh_.data::remove("finalIteration");
+        mesh_.solutionDict().remove("finalIteration");
         return false;
     }
 
@@ -202,7 +202,7 @@ bool Foam::pimpleControl::loop()
             Info<< algorithmName_ << ": converged in " << corr_ - 1
                 << " iterations" << endl;
 
-            mesh_.data::remove("finalIteration");
+            mesh_.solutionDict().remove("finalIteration");
             corr_ = 0;
             converged_ = false;
 
@@ -213,7 +213,7 @@ bool Foam::pimpleControl::loop()
             Info<< algorithmName_ << ": iteration " << corr_ << endl;
             storePrevIterFields();
 
-            mesh_.data::add("finalIteration", true);
+            mesh_.solutionDict().add("finalIteration", true);
             converged_ = true;
         }
     }
@@ -221,7 +221,7 @@ bool Foam::pimpleControl::loop()
     {
         if (finalIter())
         {
-            mesh_.data::add("finalIteration", true);
+            mesh_.solutionDict().add("finalIteration", true);
         }
 
         if (corr_ <= nCorrPIMPLE_)

src/finiteVolume/cfdTools/general/solutionControl/simpleControl/simpleControl.C

@@ -52,7 +52,7 @@ bool Foam::simpleControl::criteriaSatisfied()
     bool achieved = true;
     bool checked = false;    // safety that some checks were indeed performed
 
-    const dictionary& solverDict = mesh_.solverPerformanceDict();
+    const dictionary& solverDict = mesh_.solutionDict().solverPerformanceDict();
     forAllConstIter(dictionary, solverDict, iter)
     {
         const word& variableName = iter().keyword();

src/finiteVolume/finiteVolume/fvSchemes/fvSchemes.C

@@ -28,7 +28,6 @@ License
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
-
 Foam::debug::debugSwitch
 Foam::fvSchemes::debug
 (
@@ -586,5 +585,35 @@ bool Foam::fvSchemes::fluxRequired(const word& name) const
     }
 }
 
+bool Foam::fvSchemes::writeData(Ostream& os) const
+{
+    // Write dictionaries
+    os << nl << "ddtSchemes";
+    ddtSchemes_.write(os, true);
+
+    os << nl << "d2dt2Schemes";
+    d2dt2Schemes_.write(os, true);
+
+    os << nl << "interpolationSchemes";
+    interpolationSchemes_.write(os, true);
+
+    os << nl << "divSchemes";
+    divSchemes_.write(os, true);
+
+    os << nl << "gradSchemes";
+    gradSchemes_.write(os, true);
+
+    os << nl << "snGradSchemes";
+    snGradSchemes_.write(os, true);
+
+    os << nl << "laplacianSchemes";
+    laplacianSchemes_.write(os, true);
+
+    os << nl << "fluxRequired";
+    fluxRequired_.write(os, true);
+
+    return true;
+}
+
 
 // ************************************************************************* //

src/finiteVolume/finiteVolume/fvSchemes/fvSchemes.H

@@ -95,7 +95,7 @@ class fvSchemes
 public:
 
     //- Debug switch
-        static debug::debugSwitch debug;
+    static debug::debugSwitch debug;
 
 
     // Constructors
@@ -184,6 +184,12 @@ public:
 
             //- Read the fvSchemes
             bool read();
+
+
+        // Write
+
+            //- WriteData function required for regIOobject write operation
+            virtual bool writeData(Ostream&) const;
 };
 
 

src/finiteVolume/fvMatrices/fvMatrix/fvMatrix.C

@@ -1390,7 +1390,7 @@ Foam::tmp<Foam::fvMatrix<Type> > Foam::relax(const fvMatrix<Type>& m)
 
 
 template<class Type>
-Foam::lduMatrix::solverPerformance Foam::solve
+Foam::lduSolverPerformance Foam::solve
 (
     fvMatrix<Type>& fvm,
     const dictionary& solverControls
@@ -1400,13 +1400,13 @@ Foam::lduMatrix::solverPerformance Foam::solve
 }
 
 template<class Type>
-Foam::lduMatrix::solverPerformance Foam::solve
+Foam::lduSolverPerformance Foam::solve
 (
     const tmp<fvMatrix<Type> >& tfvm,
     const dictionary& solverControls
 )
 {
-    lduMatrix::solverPerformance solverPerf =
+    lduSolverPerformance solverPerf =
         const_cast<fvMatrix<Type>&>(tfvm()).solve(solverControls);
 
     tfvm.clear();
@@ -1416,18 +1416,18 @@ Foam::lduMatrix::solverPerformance Foam::solve
 
 
 template<class Type>
-Foam::lduMatrix::solverPerformance Foam::solve(fvMatrix<Type>& fvm)
+Foam::lduSolverPerformance Foam::solve(fvMatrix<Type>& fvm)
 {
     return fvm.solve();
 }
 
 template<class Type>
-Foam::lduMatrix::solverPerformance Foam::solve
+Foam::lduSolverPerformance Foam::solve
 (
     const tmp<fvMatrix<Type> >& tfvm
 )
 {
-    lduMatrix::solverPerformance solverPerf =
+    lduSolverPerformance solverPerf =
         const_cast<fvMatrix<Type>&>(tfvm()).solve();
 
     tfvm.clear();

src/finiteVolume/fvMatrices/fvMatrix/fvMatrix.H

@@ -363,7 +363,7 @@ public:
 
             //- Solve returning the solution statistics.
             //  Use the given solver controls
-            lduMatrix::solverPerformance solve(const dictionary&);
+            lduSolverPerformance solve(const dictionary&);
 
             //- Solve returning the solution statistics.
             //  Solver controls read from fvSolution

src/finiteVolume/fvMatrices/fvMatrix/fvMatrixSolve.C

@@ -52,7 +52,7 @@ void Foam::fvMatrix<Type>::setComponentReference
 
 
 template<class Type>
-Foam::lduMatrix::solverPerformance Foam::fvMatrix<Type>::solve
+Foam::lduSolverPerformance Foam::fvMatrix<Type>::solve
 (
     const dictionary& solverControls
 )
@@ -140,7 +140,7 @@ Foam::lduMatrix::solverPerformance Foam::fvMatrix<Type>::solve
             cmpt
         );
 
-        lduMatrix::solverPerformance solverPerf;
+        lduSolverPerformance solverPerf;
 
         // Solver call
         solverPerf = lduMatrix::solver::New
@@ -170,14 +170,14 @@ Foam::lduMatrix::solverPerformance Foam::fvMatrix<Type>::solve
 
     psi.correctBoundaryConditions();
 
-    psi_.mesh().setSolverPerformance(psi_.name(), solverPerfVec);
+    psi_.mesh().solutionDict().setSolverPerformance(psi_.name(), solverPerfVec);
 
     return solverPerfVec;
 }
 
 
 template<class Type>
-Foam::lduMatrix::solverPerformance Foam::fvMatrix<Type>::solve()
+Foam::lduSolverPerformance Foam::fvMatrix<Type>::solve()
 {
     return solve(psi_.mesh().solutionDict().solverDict(psi_.name()));
 }

src/finiteVolume/fvMatrices/fvScalarMatrix/fvScalarMatrix.C

@@ -55,7 +55,7 @@ void Foam::fvMatrix<Foam::scalar>::setComponentReference
 
 
 template<>
-Foam::lduMatrix::solverPerformance Foam::fvMatrix<Foam::scalar>::solve
+Foam::lduSolverPerformance Foam::fvMatrix<Foam::scalar>::solve
 (
     const dictionary& solverControls
 )
@@ -104,7 +104,7 @@ Foam::lduMatrix::solverPerformance Foam::fvMatrix<Foam::scalar>::solve
 
     psi.correctBoundaryConditions();
 
-    psi_.mesh().setSolverPerformance(psi_.name(), solverPerf);
+    psi_.mesh().solutionDict().setSolverPerformance(psi_.name(), solverPerf);
 
     return solverPerf;
 }

src/finiteVolume/fvMatrices/fvScalarMatrix/fvScalarMatrix.H

@@ -56,7 +56,7 @@ void fvMatrix<scalar>::setComponentReference
 );
 
 template<>
-lduMatrix::solverPerformance fvMatrix<scalar>::solve
+lduSolverPerformance fvMatrix<scalar>::solve
 (
     const dictionary&
 );

src/finiteVolume/fvMesh/fvMesh.C

@@ -103,7 +103,6 @@ Foam::fvMesh::fvMesh(const IOobject& io)
 :
     polyMesh(io),
     surfaceInterpolation(*this),
-    data(static_cast<const objectRegistry&>(*this)),
     boundary_(*this, boundaryMesh()),
     lduPtr_(NULL),
     curTimeIndex_(time().timeIndex()),
@@ -206,7 +205,6 @@ Foam::fvMesh::fvMesh
 :
     polyMesh(io, points, faces, allOwner, allNeighbour, syncPar),
     surfaceInterpolation(*this),
-    data(static_cast<const objectRegistry&>(*this)),
     boundary_(*this),
     lduPtr_(NULL),
     curTimeIndex_(time().timeIndex()),
@@ -237,7 +235,6 @@ Foam::fvMesh::fvMesh
 :
     polyMesh(io, points, faces, cells, syncPar),
     surfaceInterpolation(*this),
-    data(static_cast<const objectRegistry&>(*this)),
     boundary_(*this),
     lduPtr_(NULL),
     curTimeIndex_(time().timeIndex()),
@@ -285,7 +282,6 @@ Foam::fvMesh::fvMesh
         syncPar
     ),
     surfaceInterpolation(*this),
-    data(static_cast<const objectRegistry&>(*this)),
     boundary_(*this),
     lduPtr_(NULL),
     curTimeIndex_(time().timeIndex()),

src/finiteVolume/fvMesh/fvMesh.H

@@ -52,7 +52,6 @@ SourceFiles
 #include "primitiveMesh.H"
 #include "fvBoundaryMesh.H"
 #include "surfaceInterpolation.H"
-#include "data.H"
 #include "DimensionedField.H"
 #include "volFieldsFwd.H"
 #include "surfaceFieldsFwd.H"
@@ -77,8 +76,7 @@ class fvMesh
 :
     public polyMesh,
     public lduMesh,
-    public surfaceInterpolation,
-    public data
+    public surfaceInterpolation
 {
     // Private data
 

src/foam/Make/files

@@ -242,7 +242,6 @@ lduMatrix = matrices/lduMatrix
 $(lduMatrix)/lduMatrix/lduMatrix.C
 $(lduMatrix)/lduMatrix/lduMatrixOperations.C
 $(lduMatrix)/lduMatrix/lduMatrixATmul.C
-$(lduMatrix)/lduMatrix/lduMatrixTests.C
 $(lduMatrix)/lduMatrix/lduMatrixUpdateMatrixInterfaces.C
 $(lduMatrix)/lduMatrix/lduMatrixSolver.C
 $(lduMatrix)/lduMatrix/lduMatrixSmoother.C
@@ -641,8 +640,6 @@ $(writers)/gnuplotGraph/gnuplotGraph.C
 $(writers)/xmgrGraph/xmgrGraph.C
 $(writers)/jplotGraph/jplotGraph.C
 
-meshes/data/data.C
-
 primitives/BlockCoeff/blockCoeffBase.C
 primitives/BlockCoeff/scalarBlockCoeff.C
 primitives/BlockCoeff/sphericalTensorBlockCoeff.C
@@ -733,6 +730,8 @@ $(BlockLduSmoothers)/BlockILUCpSmoother/blockILUCpSmoothers.C
 BlockLduSolvers = matrices/blockLduMatrix/BlockLduSolvers
 $(BlockLduSolvers)/blockVectorNSolvers.C
 $(BlockLduSolvers)/BlockLduSolver/blockLduSolvers.C
+$(BlockLduSolvers)/BlockLduSolver/BlockSolverPerformances.C
+$(BlockLduSolvers)/BlockLduSolver/BlockSolverPerformanceVectorN.C
 $(BlockLduSolvers)/BlockDiagonal/blockDiagonalSolvers.C
 $(BlockLduSolvers)/BlockGaussSeidel/blockGaussSeidelSolvers.C
 $(BlockLduSolvers)/BlockCG/blockCGSolvers.C

src/foam/fields/GeometricFields/GeometricField/GeometricField.C

@@ -27,7 +27,6 @@ License
 #include "foamTime.H"
 #include "demandDrivenData.H"
 #include "dictionary.H"
-#include "data.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 

src/foam/matrices/blockLduMatrix/BlockLduSolvers/BlockLduSolver/BlockSolverPerformance.C

@@ -34,7 +34,7 @@ bool Foam::BlockSolverPerformance<Type>::checkConvergence
     const scalar RelTolerance
 )
 {
-    if (BlockLduMatrix<Type>::debug >= 2)
+    if (debug >= 2)
     {
         Info<< solverName_
             << ":  Iteration " << nIterations_
@@ -92,7 +92,7 @@ bool Foam::BlockSolverPerformance<Type>::checkSingularity
 template<class Type>
 void Foam::BlockSolverPerformance<Type>::print() const
 {
-    if (BlockLduMatrix<Type>::debug)
+    if (debug)
     {
         Info<< solverName_ << ":  Solving for " << fieldName_;
 
@@ -111,4 +111,65 @@ void Foam::BlockSolverPerformance<Type>::print() const
 }
 
 
+template<class Type>
+bool Foam::BlockSolverPerformance<Type>::operator!=
+(
+    const BlockSolverPerformance<Type>& bsp
+) const
+{
+    return
+    (
+        solverName()      != bsp.solverName()
+     || fieldName()       != bsp.fieldName()
+     || initialResidual() != bsp.initialResidual()
+     || finalResidual()   != bsp.finalResidual()
+     || nIterations()     != bsp.nIterations()
+     || converged()       != bsp.converged()
+     || singular()        != bsp.singular()
+    );
+}
+
+
+template<class Type>
+Foam::Istream& Foam::operator>>
+(
+    Istream& is,
+    typename Foam::BlockSolverPerformance<Type>& bsp
+)
+{
+    is.readBeginList("BlockSolverPerformance<Type>");
+    is  >> bsp.solverName_
+        >> bsp.fieldName_
+        >> bsp.initialResidual_
+        >> bsp.finalResidual_
+        >> bsp.nIterations_
+        >> bsp.converged_
+        >> bsp.singular_;
+    is.readEndList("BlockSolverPerformance<Type>");
+
+    return is;
+}
+
+
+template<class Type>
+Foam::Ostream& Foam::operator<<
+(
+    Ostream& os,
+    const typename Foam::BlockSolverPerformance<Type>& bsp
+)
+{
+    os  << token::BEGIN_LIST
+        << bsp.solverName_ << token::SPACE
+        << bsp.fieldName_ << token::SPACE
+        << bsp.initialResidual_ << token::SPACE
+        << bsp.finalResidual_ << token::SPACE
+        << bsp.nIterations_ << token::SPACE
+        << bsp.converged_ << token::SPACE
+        << bsp.singular_ << token::SPACE
+        << token::END_LIST;
+
+    return os;
+}
+
+
 // ************************************************************************* //

src/foam/matrices/blockLduMatrix/BlockLduSolvers/BlockLduSolver/BlockSolverPerformance.H

@@ -51,6 +51,9 @@ namespace Foam
 template<class Type>
 class BlockSolverPerformance;
 
+template<class Type>
+Istream& operator>>(Istream&, BlockSolverPerformance<Type>&);
+
 template<class Type>
 Ostream& operator<<(Ostream&, const BlockSolverPerformance<Type>&);
 
@@ -88,8 +91,24 @@ class BlockSolverPerformance
 
 public:
 
+    // Static data
+
+        // Declare name of the class and its debug switch
+        ClassName("BlockSolverPerformance");
+
+
     // Constructors
 
+        //- Construct without solver and field name
+        BlockSolverPerformance()
+        :
+            initialResidual_(pTraits<Type>::zero),
+            finalResidual_(pTraits<Type>::zero),
+            nIterations_(0),
+            converged_(false),
+            singular_(false)
+        {}
+
         //- Construct with solver and field name
         BlockSolverPerformance
         (
@@ -120,6 +139,18 @@ public:
             return solverName_;
         }
 
+        //- Return solver name
+        word& solverName()
+        {
+            return solverName_;
+        }
+
+        //- Return field name
+        const word& fieldName() const
+        {
+            return fieldName_;
+        }
+
         //- Return initial residual
         const Type& initialResidual() const
         {
@@ -192,6 +223,26 @@ public:
 
         //- Print summary of solver performance
         void print() const;
+
+
+    // Member Operators
+
+        bool operator!=(const BlockSolverPerformance<Type>&) const;
+
+
+    // Friend functions
+
+        friend Istream& operator>> <Type>
+        (
+            Istream&,
+            BlockSolverPerformance<Type>&
+        );
+
+        friend Ostream& operator<< <Type>
+        (
+            Ostream&,
+            const BlockSolverPerformance<Type>&
+        );
 };
 
 

src/foam/matrices/blockLduMatrix/BlockLduSolvers/BlockLduSolver/BlockSolverPerformanceVectorN.C

@@ -0,0 +1,53 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | foam-extend: Open Source CFD
+   \\    /   O peration     |
+    \\  /    A nd           | For copyright notice see file Copyright
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of foam-extend.
+
+    foam-extend is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation, either version 3 of the License, or (at your
+    option) any later version.
+
+    foam-extend is distributed in the hope that it will be useful, but
+    WITHOUT ANY WARRANTY; without even the implied warranty of
+    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+    General Public License for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with foam-extend.  If not, see <http://www.gnu.org/licenses/>.
+
+Description
+    BlockSolverPerformance static data members
+
+Author
+    Vanja Skuric, FMENA Zagreb. All rights reserved.
+
+\*---------------------------------------------------------------------------*/
+
+#include "coeffFields.H"
+#include "BlockSolverPerformanceVectorN.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+#define makeNamedTemplateTypeNameAndDebug(type, Type, args...)                \
+    defineNamedTemplateTypeNameAndDebug(BlockSolverPerformance##Type, 1);
+
+forAllVectorNTypes(makeNamedTemplateTypeNameAndDebug);
+
+#undef makeNamedTemplateTypeNameAndDebug
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// ************************************************************************* //

src/foam/matrices/blockLduMatrix/BlockLduSolvers/BlockLduSolver/BlockSolverPerformanceVectorN.H

@@ -0,0 +1,66 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | foam-extend: Open Source CFD
+   \\    /   O peration     |
+    \\  /    A nd           | For copyright notice see file Copyright
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of foam-extend.
+
+    foam-extend is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation, either version 3 of the License, or (at your
+    option) any later version.
+
+    foam-extend is distributed in the hope that it will be useful, but
+    WITHOUT ANY WARRANTY; without even the implied warranty of
+    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+    General Public License for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with foam-extend.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    BlockSolverPerformance
+
+Description
+    Typedefs for BlockSolverPerformance
+
+Author
+    Vanja Skuric, FMENA Zagreb. All rights reserved.
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef BlockSolverPerformanceVectorN_H
+#define BlockSolverPerformanceVectorN_H
+
+#include "BlockSolverPerformances.H"
+#include "VectorTensorNFieldsFwd.H"
+#include "ExpandTensorNField.H"
+#include "VectorNFieldTypes.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#define makeTypedef(type, Type, args...)                                      \
+    typedef BlockSolverPerformance<type> BlockSolverPerformance##Type;
+
+forAllVectorNTypes(makeTypedef)
+
+#undef makeTypedef
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/foam/matrices/blockLduMatrix/BlockLduSolvers/BlockLduSolver/BlockSolverPerformances.C

@@ -0,0 +1,52 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | foam-extend: Open Source CFD
+   \\    /   O peration     |
+    \\  /    A nd           | For copyright notice see file Copyright
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of foam-extend.
+
+    foam-extend is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation, either version 3 of the License, or (at your
+    option) any later version.
+
+    foam-extend is distributed in the hope that it will be useful, but
+    WITHOUT ANY WARRANTY; without even the implied warranty of
+    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+    General Public License for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with foam-extend.  If not, see <http://www.gnu.org/licenses/>.
+
+Description
+    BlockSolverPerformance static data members
+
+Author
+    Vanja Skuric, FMENA Zagreb. All rights reserved.
+
+\*---------------------------------------------------------------------------*/
+
+#include "BlockSolverPerformances.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+defineNamedTemplateTypeNameAndDebug(BlockSolverPerformanceScalar, 1);
+defineNamedTemplateTypeNameAndDebug(BlockSolverPerformanceVector, 1);
+defineNamedTemplateTypeNameAndDebug(BlockSolverPerformanceSphericalTensor, 1);
+defineNamedTemplateTypeNameAndDebug(BlockSolverPerformanceSymmTensor, 1);
+defineNamedTemplateTypeNameAndDebug(BlockSolverPerformanceTensor, 1);
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// ************************************************************************* //

src/foam/matrices/blockLduMatrix/BlockLduSolvers/BlockLduSolver/BlockSolverPerformances.H

@@ -0,0 +1,64 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | foam-extend: Open Source CFD
+   \\    /   O peration     |
+    \\  /    A nd           | For copyright notice see file Copyright
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of foam-extend.
+
+    foam-extend is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation, either version 3 of the License, or (at your
+    option) any later version.
+
+    foam-extend is distributed in the hope that it will be useful, but
+    WITHOUT ANY WARRANTY; without even the implied warranty of
+    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+    General Public License for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with foam-extend.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    BlockSolverPerformance
+
+Description
+    Typedefs for BlockSolverPerformance
+
+Author
+    Vanja Skuric, FMENA Zagreb. All rights reserved.
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef BlockSolverPerformances_H
+#define BlockSolverPerformances_H
+
+#include "BlockSolverPerformance.H"
+#include "fieldTypes.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+typedef BlockSolverPerformance<scalar> BlockSolverPerformanceScalar;
+typedef BlockSolverPerformance<vector> BlockSolverPerformanceVector;
+typedef BlockSolverPerformance<sphericalTensor>
+    BlockSolverPerformanceSphericalTensor;
+typedef BlockSolverPerformance<symmTensor> BlockSolverPerformanceSymmTensor;
+typedef BlockSolverPerformance<tensor> BlockSolverPerformanceTensor;
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/foam/matrices/lduMatrix/lduMatrix/lduMatrix.H

@@ -42,7 +42,6 @@ SourceFiles
     lduMatrixOperations.C
     lduMatrixSolver.C
     lduMatrixPreconditioner.C
-    lduMatrixTests.C
     lduMatrixUpdateMatrixInterfaces.C
     lduMatrixBufferedUpdateMatrixInterfaces.C
 
@@ -60,6 +59,7 @@ SourceFiles
 #include "autoPtr.H"
 #include "runTimeSelectionTables.H"
 #include "profilingTrigger.H"
+#include "lduSolverPerformance.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -93,181 +93,6 @@ class lduMatrix
 
 public:
 
-    //- Class returned by the solver, containing performance statistics
-    class solverPerformance
-    {
-        // Private Data
-
-            //- Name of solver
-            word solverName_;
-
-            //- Name of field
-            word fieldName_;
-
-            //- Initial residual
-            scalar initialResidual_;
-
-            //- Final residual
-            scalar finalResidual_;
-
-            //- Number of iterations
-            label noIterations_;
-
-            //- Has the solver converged?
-            bool converged_;
-
-            //- Has the solver indicated a singular matrix?
-            bool singular_;
-
-
-    public:
-
-        // Constructors
-
-            //- Construct null
-            solverPerformance()
-            :
-                initialResidual_(0),
-                finalResidual_(0),
-                noIterations_(0),
-                converged_(false),
-                singular_(false)
-            {}
-
-
-            //- Construct from components
-            solverPerformance
-            (
-                const word& solverName,
-                const word& fieldName,
-                const scalar iRes = 0,
-                const scalar fRes = 0,
-                const label nIter = 0,
-                const bool converged = false,
-                const bool singular = false
-            )
-            :
-                solverName_(solverName),
-                fieldName_(fieldName),
-                initialResidual_(iRes),
-                finalResidual_(fRes),
-                noIterations_(nIter),
-                converged_(converged),
-                singular_(singular)
-            {}
-
-
-        // Member functions
-
-            //- Return solver name
-            const word& solverName() const
-            {
-                return solverName_;
-            }
-
-            //- Return access to solver name
-            word& solverName()
-            {
-                return solverName_;
-            }
-
-            //- Return field name
-            const word& fieldName() const
-            {
-                return fieldName_;
-            }
-
-            //- Return access to field name
-            word& fieldName()
-            {
-                return fieldName_;
-            }
-
-            //- Return initial residual
-            scalar initialResidual() const
-            {
-                return initialResidual_;
-            }
-
-            //- Return initial residual
-            scalar& initialResidual()
-            {
-                return initialResidual_;
-            }
-
-
-            //- Return final residual
-            scalar finalResidual() const
-            {
-                return finalResidual_;
-            }
-
-            //- Return final residual
-            scalar& finalResidual()
-            {
-                return finalResidual_;
-            }
-
-
-            //- Return number of iterations
-            label nIterations() const
-            {
-                return noIterations_;
-            }
-
-            //- Return number of iterations
-            label& nIterations()
-            {
-                return noIterations_;
-            }
-
-
-            //- Has the solver converged?
-            bool converged() const
-            {
-                return converged_;
-            }
-
-            //- Is the matrix singular?
-            bool singular() const
-            {
-                return singular_;
-            }
-
-            //- Convergence test
-            bool checkConvergence
-            (
-                const scalar tolerance,
-                const scalar relTolerance
-            );
-
-            //- Singularity test
-            bool checkSingularity(const scalar residual);
-
-            //- Print summary of solver performance
-            void print() const;
-
-        // Member Operators
-
-            bool operator!=(const solverPerformance&) const;
-
-        // Ostream Operator
-
-            friend Istream& operator>>
-            (
-                Istream&,
-                solverPerformance&
-            );
-
-            friend Ostream& operator<<
-            (
-                Ostream&,
-                const solverPerformance&
-            );
-
-    };
-
-
     //- Abstract base-class for lduMatrix solvers
     class solver
     {
@@ -322,10 +147,10 @@ public:
             }
 
             //- Is the stop criterion reached
-            bool stop(lduMatrix::solverPerformance& solverPerf) const;
+            bool stop(lduSolverPerformance& solverPerf) const;
 
             //- Is the converrgence criterion reached
-            bool converged(lduMatrix::solverPerformance& solverPerf) const;
+            bool converged(lduSolverPerformance& solverPerf) const;
 
             //- Read the control parameters from the dictionary
             virtual void readControls();
@@ -491,7 +316,7 @@ public:
             //- Read and reset the solver parameters from the given stream
             virtual void read(const dictionary&);
 
-            virtual solverPerformance solve
+            virtual lduSolverPerformance solve
             (
                 scalarField& x,
                 const scalarField& b,
@@ -1048,7 +873,6 @@ public:
 
 // Helper typedefs
 typedef lduMatrix::solver lduSolver;
-typedef lduMatrix::solverPerformance lduSolverPerformance;
 typedef lduMatrix::preconditioner lduPreconditioner;
 typedef lduMatrix::smoother lduSmoother;
 

src/foam/matrices/lduMatrix/lduMatrix/lduMatrixSolver.C

@@ -207,7 +207,7 @@ Foam::scalar Foam::lduMatrix::solver::normFactor
 
 bool Foam::lduMatrix::solver::stop
 (
-    lduMatrix::solverPerformance& solverPerf
+    lduSolverPerformance& solverPerf
 ) const
 {
     if (solverPerf.nIterations() < minIter_)
@@ -221,7 +221,7 @@ bool Foam::lduMatrix::solver::stop
 
 bool Foam::lduMatrix::solver::converged
 (
-    lduMatrix::solverPerformance& solverPerf
+    lduSolverPerformance& solverPerf
 ) const
 {
     if

src/foam/matrices/lduMatrix/lduMatrix/lduMatrixTests.C

@@ -1,165 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | foam-extend: Open Source CFD
-   \\    /   O peration     | Version:     3.2
-    \\  /    A nd           | Web:         http://www.foam-extend.org
-     \\/     M anipulation  | For copyright notice see file Copyright
--------------------------------------------------------------------------------
-License
-    This file is part of foam-extend.
-
-    foam-extend is free software: you can redistribute it and/or modify it
-    under the terms of the GNU General Public License as published by the
-    Free Software Foundation, either version 3 of the License, or (at your
-    option) any later version.
-
-    foam-extend is distributed in the hope that it will be useful, but
-    WITHOUT ANY WARRANTY; without even the implied warranty of
-    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
-    General Public License for more details.
-
-    You should have received a copy of the GNU General Public License
-    along with foam-extend.  If not, see <http://www.gnu.org/licenses/>.
-
-Description
-    Convergence and singularity tests for solvers.
-
-\*---------------------------------------------------------------------------*/
-
-#include "lduMatrix.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-bool Foam::lduMatrix::solverPerformance::checkConvergence
-(
-    const scalar Tolerance,
-    const scalar RelTolerance
-)
-{
-    if (debug >= 2)
-    {
-        Info<< solverName_
-            << ":  Iteration " << noIterations_
-            << " residual = " << finalResidual_
-            << endl;
-    }
-
-    if
-    (
-        finalResidual_ < Tolerance   // Abs. tolerance
-     || (
-            RelTolerance > SMALL
-         && finalResidual_ <= RelTolerance*initialResidual_ // Rel. tolerance
-        )
-  // || (solverName == "symSolve" && iter == 0)
-    )
-    {
-        converged_ = true;
-    }
-    else
-    {
-        converged_ = false;
-    }
-
-    return converged_;
-}
-
-
-bool Foam::lduMatrix::solverPerformance::checkSingularity
-(
-    const scalar residual
-)
-{
-    if (residual > VSMALL)
-    {
-        singular_ = false;
-    }
-    else
-    {
-        singular_ = true;
-    }
-
-    return singular_;
-}
-
-
-void Foam::lduMatrix::solverPerformance::print() const
-{
-    if (debug)
-    {
-        Info<< solverName_ << ":  Solving for " << fieldName_;
-
-        if (singular())
-        {
-            Info<< ":  solution singularity" << endl;
-        }
-        else
-        {
-            Info<< ", Initial residual = " << initialResidual_
-                << ", Final residual = " << finalResidual_
-                << ", No Iterations " << noIterations_
-                << endl;
-        }
-    }
-}
-
-
-bool Foam::lduSolverPerformance::operator!=
-(
-    const lduSolverPerformance& sp
-) const
-{
-    return
-    (
-        solverName()      != sp.solverName()
-     || fieldName()       != sp.fieldName()
-     || initialResidual() != sp.initialResidual()
-     || finalResidual()   != sp.finalResidual()
-     || nIterations()     != sp.nIterations()
-     || converged()       != sp.converged()
-     || singular()        != sp.singular()
-    );
-}
-
-
-Foam::Istream& Foam::operator>>
-(
-    Istream& is,
-    lduSolverPerformance& sp
-)
-{
-    is.readBeginList("SolverPerformance<Type>");
-    is  >> sp.solverName_
-        >> sp.fieldName_
-        >> sp.initialResidual_
-        >> sp.finalResidual_
-        >> sp.noIterations_
-        >> sp.converged_
-        >> sp.singular_;
-    is.readEndList("SolverPerformance<Type>");
-
-    return is;
-}
-
-
-Foam::Ostream& Foam::operator<<
-(
-    Ostream& os,
-    const lduSolverPerformance& sp
-)
-{
-    os  << token::BEGIN_LIST
-        << sp.solverName_ << token::SPACE
-        << sp.fieldName_ << token::SPACE
-        << sp.initialResidual_ << token::SPACE
-        << sp.finalResidual_ << token::SPACE
-        << sp.noIterations_ << token::SPACE
-        << sp.converged_ << token::SPACE
-        << sp.singular_ << token::SPACE
-        << token::END_LIST;
-
-    return os;
-}
-
-
-// ************************************************************************* //

src/foam/matrices/lduMatrix/lduMatrix/lduSolverPerformance.H

@@ -0,0 +1,49 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | foam-extend: Open Source CFD
+   \\    /   O peration     | Version:     3.2
+    \\  /    A nd           | Web:         http://www.foam-extend.org
+     \\/     M anipulation  | For copyright notice see file Copyright
+-------------------------------------------------------------------------------
+License
+    This file is part of foam-extend.
+
+    foam-extend is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation, either version 3 of the License, or (at your
+    option) any later version.
+
+    foam-extend is distributed in the hope that it will be useful, but
+    WITHOUT ANY WARRANTY; without even the implied warranty of
+    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+    General Public License for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with foam-extend.  If not, see <http://www.gnu.org/licenses/>.
+
+Typedef
+    Foam::lduSolverPerformance
+
+Description
+    BlockSolverPerformance instantiated for a scalar and lduMatrix
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef lduSolverPerformance_H
+#define lduSolverPerformance_H
+
+#include "BlockSolverPerformance.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    typedef BlockSolverPerformance<scalar> lduSolverPerformance;
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/foam/matrices/lduMatrix/solvers/AMG/GAMGSolverSolve.C

@@ -30,7 +30,7 @@ License
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-Foam::lduMatrix::solverPerformance Foam::GAMGSolver::solve
+Foam::lduSolverPerformance Foam::GAMGSolver::solve
 (
     scalarField& x,
     const scalarField& b,

src/foam/matrices/lduMatrix/solvers/diagonalSolver/diagonalSolver.C

@@ -59,7 +59,7 @@ Foam::diagonalSolver::diagonalSolver
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-Foam::lduMatrix::solverPerformance Foam::diagonalSolver::solve
+Foam::lduSolverPerformance Foam::diagonalSolver::solve
 (
     scalarField& x,
     const scalarField& b,

src/foam/matrices/solution/solution.C

@@ -146,7 +146,8 @@ Foam::solution::solution(const objectRegistry& obr, const fileName& dictName)
     eqnRelaxDict_(dictionary::null),
     fieldRelaxDefault_(0),
     eqnRelaxDefault_(0),
-    solvers_(dictionary::null)
+    solvers_(dictionary::null),
+    prevTimeIndex_(0)
 {
     if (!headerOk())
     {
@@ -162,6 +163,8 @@ Foam::solution::solution(const objectRegistry& obr, const fileName& dictName)
     }
 
     read(solutionDict());
+
+    set("solverPerformance", dictionary());
 }
 
 
@@ -408,5 +411,46 @@ bool Foam::solution::writeData(Ostream& os) const
     return true;
 }
 
+const Foam::dictionary& Foam::solution::solverPerformanceDict() const
+{
+    return subDict("solverPerformance");
+}
+
+
+void Foam::solution::setSolverPerformance
+(
+    const word& name,
+    const lduSolverPerformance& sp
+) const
+{
+    dictionary& dict = const_cast<dictionary&>(solverPerformanceDict());
+
+    List<lduSolverPerformance> perfs;
+
+    if (prevTimeIndex_ != this->time().timeIndex())
+    {
+        // Reset solver performance between iterations
+        prevTimeIndex_ = this->time().timeIndex();
+        dict.clear();
+    }
+    else
+    {
+        dict.readIfPresent(name, perfs);
+    }
+
+    // Append to list
+    perfs.setSize(perfs.size()+1, sp);
+
+    dict.set(name, perfs);
+}
+
+
+void Foam::solution::setSolverPerformance
+(
+    const lduSolverPerformance& sp
+) const
+{
+    setSolverPerformance(sp.fieldName(), sp);
+}
 
 // ************************************************************************* //

src/foam/matrices/solution/solution.H

@@ -37,6 +37,7 @@ SourceFiles
 
 #include "IOdictionary.H"
 #include "debugSwitch.H"
+#include "lduSolverPerformance.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -74,6 +75,9 @@ class solution
         //- Dictionary of solver parameters for all the fields
         dictionary solvers_;
 
+        //- Previously used time-index, used for reset between iterations
+        mutable label prevTimeIndex_;
+
 
     // Private Member Functions
 
@@ -144,6 +148,24 @@ public:
             //- Return the solver controls dictionary for the given field
             const dictionary& solver(const word& name) const;
 
+            //- Return the dictionary of solver performance data
+            //  which includes initial and final residuals for convergence
+            //  checking
+            const dictionary& solverPerformanceDict() const;
+
+            //- Add/set the solverPerformance entry for the named field
+            void setSolverPerformance
+            (
+                const word& name,
+                const lduSolverPerformance&
+            ) const;
+
+            //- Add/set the solverPerformance entry, using its fieldName
+            void setSolverPerformance
+            (
+                const lduSolverPerformance&
+            ) const;
+
 
         // Edit
 

src/foam/meshes/data/data.C

@@ -1,99 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | foam-extend: Open Source CFD
-   \\    /   O peration     | Version:     3.2
-    \\  /    A nd           | Web:         http://www.foam-extend.org
-     \\/     M anipulation  | For copyright notice see file Copyright
--------------------------------------------------------------------------------
-License
-    This file is part of foam-extend.
-
-    foam-extend is free software: you can redistribute it and/or modify it
-    under the terms of the GNU General Public License as published by the
-    Free Software Foundation, either version 3 of the License, or (at your
-    option) any later version.
-
-    foam-extend is distributed in the hope that it will be useful, but
-    WITHOUT ANY WARRANTY; without even the implied warranty of
-    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
-    General Public License for more details.
-
-    You should have received a copy of the GNU General Public License
-    along with foam-extend.  If not, see <http://www.gnu.org/licenses/>.
-
-\*---------------------------------------------------------------------------*/
-
-#include "data.H"
-#include "foamTime.H"
-
-// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
-
-int Foam::data::debug(Foam::debug::debugSwitchFromDict("data", false));
-
-// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
-
-Foam::data::data(const objectRegistry& obr)
-:
-    IOdictionary
-    (
-        IOobject
-        (
-            "data",
-            obr.time().system(),
-            obr,
-            IOobject::NO_READ,
-            IOobject::NO_WRITE
-        )
-    ),
-    prevTimeIndex_(0)
-{
-    set("solverPerformance", dictionary());
-}
-
-
-// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
-
-const Foam::dictionary& Foam::data::solverPerformanceDict() const
-{
-    return subDict("solverPerformance");
-}
-
-
-void Foam::data::setSolverPerformance
-(
-    const word& name,
-    const lduSolverPerformance& sp
-) const
-{
-    dictionary& dict = const_cast<dictionary&>(solverPerformanceDict());
-
-    List<lduSolverPerformance> perfs;
-
-    if (prevTimeIndex_ != this->time().timeIndex())
-    {
-        // Reset solver performance between iterations
-        prevTimeIndex_ = this->time().timeIndex();
-        dict.clear();
-    }
-    else
-    {
-        dict.readIfPresent(name, perfs);
-    }
-
-    // Append to list
-    perfs.setSize(perfs.size()+1, sp);
-
-    dict.set(name, perfs);
-}
-
-
-void Foam::data::setSolverPerformance
-(
-    const lduSolverPerformance& sp
-) const
-{
-    setSolverPerformance(sp.fieldName(), sp);
-}
-
-
-// ************************************************************************* //

src/foam/meshes/data/data.H

@@ -1,116 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | foam-extend: Open Source CFD
-   \\    /   O peration     | Version:     3.2
-    \\  /    A nd           | Web:         http://www.foam-extend.org
-     \\/     M anipulation  | For copyright notice see file Copyright
--------------------------------------------------------------------------------
-License
-    This file is part of foam-extend.
-
-    foam-extend is free software: you can redistribute it and/or modify it
-    under the terms of the GNU General Public License as published by the
-    Free Software Foundation, either version 3 of the License, or (at your
-    option) any later version.
-
-    foam-extend is distributed in the hope that it will be useful, but
-    WITHOUT ANY WARRANTY; without even the implied warranty of
-    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
-    General Public License for more details.
-
-    You should have received a copy of the GNU General Public License
-    along with foam-extend.  If not, see <http://www.gnu.org/licenses/>.
-
-Class
-    Foam::data
-
-Description
-    Database for solution data, solver performance and other reduced data.
-
-    fvMesh is derived from data so that all fields have access to the data from
-    the mesh reference they hold.
-
-SourceFiles
-    data.C
-
-\*---------------------------------------------------------------------------*/
-
-#ifndef data_H
-#define data_H
-
-#include "IOdictionary.H"
-#include "lduMatrix.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-namespace Foam
-{
-
-/*---------------------------------------------------------------------------*\
-                            Class data Declaration
-\*---------------------------------------------------------------------------*/
-
-class data
-:
-    public IOdictionary
-{
-    // Private data
-
-        //- Previously used time-index, used for reset between iterations
-        mutable label prevTimeIndex_;
-
-
-    // Private Member Functions
-
-        //- Disallow default bitwise copy construct
-        data(const data&);
-
-        //- Disallow default bitwise assignment
-        void operator=(const data&);
-
-
-public:
-
-    //- Debug switch
-    static int debug;
-
-
-    // Constructors
-
-        //- Construct for objectRegistry
-        data(const objectRegistry& obr);
-
-
-    // Member Functions
-
-        // Access
-
-            //- Return the dictionary of solver performance data
-            //  which includes initial and final residuals for convergence
-            //  checking
-            const dictionary& solverPerformanceDict() const;
-
-            //- Add/set the solverPerformance entry for the named field
-            void setSolverPerformance
-            (
-                const word& name,
-                const lduSolverPerformance&
-            ) const;
-
-            //- Add/set the solverPerformance entry, using its fieldName
-            void setSolverPerformance
-            (
-                const lduSolverPerformance&
-            ) const;
-};
-
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-} // End namespace Foam
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-#endif
-
-// ************************************************************************* //

src/thermophysicalModels/basic/IAPWS_Waterproperties/IAPWSThermo/IAPWS-IF97.C

@@ -0,0 +1,599 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright held by original author
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+Author
+Christian Lucas
+Institut für Thermodynamik
+Technische Universität Braunschweig
+Germany
+
+\*---------------------------------------------------------------------------*/
+
+#include "IAPWS-IF97.H"
+
+//CL: calculated all (minimal) needed properties for a given pressure and enthalpy
+void Foam::calculateProperties_ph
+(
+    scalar &p,
+    scalar &h,
+    scalar &T,
+    scalar &rho,
+    scalar &psi,
+    scalar &drhodh,
+    scalar &mu,
+    scalar &alpha
+)
+{
+    SteamState S;
+
+    // CL: vapor mass fraction is also calculated in calculateProperties_h
+    // CL: in this fuction, x is a dummy variable and x is not return to IAPWSThermo.C
+    scalar x;
+
+    S=freesteam_set_ph(p,h);
+    calculateProperties_h(S,p,h,T,rho,psi,drhodh,mu,alpha,x);
+}
+
+
+//CL: calculated all (minimal) needed properties + the vapor mass fraction for a given pressure and enthalpy
+void Foam::calculateProperties_ph
+(
+    scalar &p,
+    scalar &h,
+    scalar &T,
+    scalar &rho,
+    scalar &psi,
+    scalar &drhodh,
+    scalar &mu,
+    scalar &alpha,
+    scalar &x
+)
+{
+    SteamState S;
+
+    S=freesteam_set_ph(p,h);
+    calculateProperties_h(S,p,h,T,rho,psi,drhodh,mu,alpha,x);
+}
+
+
+//CL: calculated all (minimal) needed properties for a given pressure and temperature
+void Foam::calculateProperties_pT
+(
+    scalar &p,
+    scalar &T,
+    scalar &h,
+    scalar &rho,
+    scalar &psi,
+    scalar &drhodh,
+    scalar &mu,
+    scalar &alpha
+)
+{
+    SteamState S;
+
+    // CL: vapor mass fraction is also calculated in calculateProperties_h
+    // CL: in this fuction, x is a dummy variable and x is not return to IAPWSThermo.
+    scalar x;
+
+    S=freesteam_set_pT(p,T);
+    calculateProperties_h(S,p,h,T,rho,psi,drhodh,mu,alpha,x);
+}
+
+
+//CL: calculated all (minimal) needed properties + the vapor mass fraction for a given pressure and temperature
+void Foam::calculateProperties_pT
+(
+    scalar &p,
+    scalar &T,
+    scalar &h,
+    scalar &rho,
+    scalar &psi,
+    scalar &drhodh,
+    scalar &mu,
+    scalar &alpha,
+    scalar &x
+)
+{
+    SteamState S;
+
+    S=freesteam_set_pT(p,T);
+    calculateProperties_h(S,p,h,T,rho,psi,drhodh,mu,alpha,x);
+}
+
+
+//CL: calculated the properties --> this function is called by the functions above
+//CL: does not calulated the internal energy, if this is needed e.g. for sonicFoam
+//CL: the function has to be changed a little bit
+void Foam::calculateProperties_h
+(
+    SteamState S,
+    scalar &p,
+    scalar &h,
+    scalar &T,
+    scalar &rho,
+    scalar &psi,
+    scalar &drhodh,
+    scalar &mu,
+    scalar &alpha,
+    scalar &x
+)
+{
+    label region;
+    scalar kappa,lambda,cp,beta;
+
+    region=freesteam_region(S);
+
+    if (region==1)
+    {
+        p=S.R1.p;
+        T=S.R1.T;
+        rho=1/freesteam_region1_v_pT(S.R1.p,S.R1.T);
+        h=freesteam_region1_h_pT(S.R1.p,S.R1.T);
+        x=0;
+
+        //Cl: note: in FreeStream, beta=1/V*(dV/dP)_P=const is called alphaV (in this region)
+        //Cl: note: in FreeStream, kappa=1/V*(dV/dP)_T=const is called kappaT (in this region)
+        kappa=freesteam_region1_kappaT_pT(S.R1.p,S.R1.T);
+        beta=freesteam_region1_alphav_pT(S.R1.p,S.R1.T);
+        cp=freesteam_region1_cp_pT(S.R1.p,S.R1.T);
+
+        //CL: getting derivatives using Bridgmans table
+        //CL: psi=(drho/dp)_h=const
+        //CL: drhodh=(drho/dh)_p=const
+        psi=-((T*beta*beta-beta)/cp-kappa*rho);
+        drhodh=-rho*beta/cp;
+
+        //CL: getting transport properties
+        mu=freesteam_mu_rhoT(rho, T);
+        lambda=freesteam_k_rhoT(rho,T);
+        alpha=lambda/cp; //Cl: Important info -->alpha= thermal diffusivity time density
+    }
+    else if (region==2)
+    {
+        p=S.R2.p;
+        T=S.R2.T;
+        rho=1/freesteam_region2_v_pT(S.R2.p,S.R2.T);
+        h=freesteam_region2_h_pT(S.R2.p,S.R2.T);
+        x=1;
+
+        //Cl: note: in FreeStream, beta=1/V*(dV/dP)_P=const is called alphaV (in this region)
+        //Cl: note: in FreeStream, kappa=1/V*(dV/dP)_T=const is called kappaT (in this region)
+        kappa=freesteam_region2_kappaT_pT(S.R2.p,S.R2.T);
+        beta=freesteam_region2_alphav_pT(S.R2.p,S.R2.T);
+        cp=freesteam_region2_cp_pT(S.R2.p,S.R2.T);
+
+        //CL: getting derivatives using Bridgmans table
+        //CL: psi=(drho/dp)_h=const
+        //CL: drhodh=(drho/dh)_p=const
+        psi=-((T*beta*beta-beta)/cp-kappa*rho);
+        drhodh=-rho*beta/cp;
+
+        //CL: getting transport properties
+        mu=freesteam_mu_rhoT(rho, T);
+        lambda=freesteam_k_rhoT(rho,T);
+        alpha=lambda/cp; //Cl: Important info -->alpha= thermal diffusivity time density
+    }
+    else if (region==3)
+    {
+        scalar gamma,cv;
+
+        rho=S.R3.rho;
+        T=S.R3.T;
+        p=freesteam_region3_p_rhoT(S.R3.rho,S.R3.T);
+        h=freesteam_region3_h_rhoT(S.R3.rho,S.R3.T);
+
+        //CL= when h<h @ critical point -->x=0 else x=1
+        if (h<2084256.263)
+        {
+            x=0;
+        }
+        else
+        {
+            x=1;
+        }
+
+        //Cl: note: beta=1/V*(dV/dP)_P=const
+        //Cl: note: kappa=1/V*(dV/dP)_T=const
+        //Cl: note: in FreeStream, gamma=1/p*(dp/dT)_v=const is called alphap (in this region)
+        gamma=freesteam_region3_alphap_rhoT(S.R3.rho,S.R3.T);
+        cp=freesteam_region3_cp_rhoT(S.R3.rho,S.R3.T);
+        cv=freesteam_region3_cv_rhoT(S.R3.rho,S.R3.T);
+        beta=(cp-cv)/(S.R3.T/S.R3.rho*p*gamma);
+        kappa=(cp-cv)/(S.R3.T/S.R3.rho*p*p*gamma*gamma);
+
+        //CL: getting derivatives using Bridgmans table
+        //CL: psi=(drho/dp)_h=const
+        //CL: drhodh=(drho/dh)_p=const
+        psi=-((T*beta*beta-beta)/cp-kappa*rho);
+        drhodh=-rho*beta/cp;
+
+
+        //CL: getting transport properties
+        mu=freesteam_mu_rhoT(rho, T);
+        lambda=freesteam_k_rhoT(rho,T);
+        alpha=lambda/cp; //Cl: Important info -->alpha= thermal diffusivity time density
+
+    }
+    else if (region==4)
+    {
+        scalar rhov,rhol,betav,betal,kappav,kappal,vv,vl,cpl,hl,hv,cp;
+        scalar dvldp,dvvdp,dhldp,dhvdp;
+        scalar dpdT,dvdh,dvdp,dxdp;
+
+        SteamState Sl,Sv;
+
+        x=S.R4.x;
+        T=S.R4.T;
+        rho=1/freesteam_region4_v_Tx(S.R4.T,S.R4.x);
+        h=freesteam_region4_h_Tx(S.R4.T,S.R4.x);
+        p=freesteam_region4_psat_T(S.R4.T);
+        cp=freesteam_region4_cp_Tx(S.R4.T,S.R4.x);
+
+
+        //CL: Getting density on the vapour and liquid lines
+        rhov=freesteam_region4_rhog_T(S.R4.T);
+        rhol=freesteam_region4_rhof_T(S.R4.T);
+        vv=1/rhov;
+        vl=1/rhol;
+
+        //CL: getting derivatives --> this is a bit tricky inside the vapor dome
+
+        dpdT=freesteam_region4_dpsatdT_T(S.R4.T);
+
+        // getting the states outside the vapour dome
+        Sl=freesteam_set_pv(p,vl-0.0000001);  //inside region 1
+        Sv=freesteam_set_pv(p,vv+0.0000001);  //inside region 2
+
+        kappal=freesteam_region1_kappaT_pT(Sl.R1.p,Sl.R1.T);
+        kappav=freesteam_region2_kappaT_pT(Sv.R2.p,Sv.R2.T);
+
+        betal=freesteam_region1_alphav_pT(Sl.R1.p,Sl.R1.T);
+        betav=freesteam_region2_alphav_pT(Sv.R2.p,Sv.R2.T);
+
+        cpl=freesteam_region1_cp_pT(Sl.R1.p,Sl.R1.T);
+        //cpv=freesteam_region2_cp_pT(Sv.R2.p,Sv.R2.T);
+
+        hl=freesteam_region1_h_pT(Sl.R1.p,Sl.R1.T);
+        hv=freesteam_region2_h_pT(Sv.R2.p,Sv.R2.T);
+
+
+        //calculation derviatives on liquid and vapour line
+        dvldp=betal*vl/dpdT-kappal*vl;
+        dvvdp=betav*vv/dpdT-kappav*vv;
+
+        dhldp=vl*(1-betal*Sl.R1.T)+cpl/dpdT;
+        dhvdp=vv*(1-betav*Sv.R2.T)+cpl/dpdT;
+
+        dxdp=-dhldp/(hv-hl)
+                 +(h-hl)/((hv-hl)*(hv-hl))
+                     *(dhvdp-dhldp);
+
+        //CL: psi=(drho/dp)_h=const
+        dvdp=dvldp+(dvvdp-dvldp)*x+(vv-vl)*dxdp;
+        psi=-rho*rho*dvdp;
+
+        //CL: drhodh=(drho/dh)_p=const
+        dvdh=(vv-vl)/(hv-hl);
+        drhodh=-rho*rho*dvdh;
+
+        //CL: getting transport properties
+        mu=freesteam_mu_rhoT(rho, T);
+        lambda=freesteam_k_rhoT(rho,T);
+        alpha=lambda/cp; //Cl: Important info -->alpha= thermal diffusivity time density
+    }
+    else
+    {
+        std::cout<<"IAPWS-IF97 error, outside the regions 1-4"<<std::endl;
+    }
+}
+
+
+//CL: returns density for given pressure and temperature
+Foam::scalar Foam::rho_pT(scalar p,scalar T)
+{
+    return 1/freesteam_v(freesteam_set_pT(p,T));
+}
+
+
+//CL: returns density for given pressure and enthalpy
+Foam::scalar Foam::rho_ph(scalar p,scalar h)
+{
+    return 1/freesteam_v(freesteam_set_ph(p,h));
+}
+
+
+//CL: returns Cp(heat capacity @ contant pressure) for given pressure and temperature
+Foam::scalar Foam::cp_pT(scalar p,scalar T)
+{
+    return freesteam_cp(freesteam_set_pT(p,T));
+}
+
+
+//CL: returns Cp(heat capacity @ contant pressure) for given pressure and enthalpy
+Foam::scalar Foam::cp_ph(scalar p,scalar h)
+{
+    return freesteam_cp(freesteam_set_ph(p,h));
+}
+
+//CL: returns Cv (heat capacity @ contant volume) for given pressure and temperature
+Foam::scalar Foam::cv_pT(scalar p,scalar T)
+{
+    return freesteam_cv(freesteam_set_pT(p,T));
+}
+
+
+//CL: returns Cv (heat capacity @ contant volume) for given pressure and enthalpy
+Foam::scalar Foam::cv_ph(scalar p,scalar h)
+{
+    return freesteam_cv(freesteam_set_ph(p,h));
+}
+
+
+//CL: returns enthalpy for given pressure and temperature
+Foam::scalar Foam::h_pT(scalar p,scalar T)
+{
+    return freesteam_h(freesteam_set_pT(p,T));
+}
+
+
+//CL: returns temperature for given pressure and enthalpy
+Foam::scalar Foam::T_ph(scalar p,scalar h)
+{
+    return freesteam_T(freesteam_set_ph(p,h));
+}
+
+
+//CL: psiH=(drho/dp)_h=const
+Foam::scalar Foam::psiH_pT(scalar p,scalar T)
+{
+    return psiH(freesteam_set_pT(p,T));
+}
+
+
+//CL: psiH=(drho/dp)_h=const
+Foam::scalar Foam::psiH_ph(scalar p,scalar h)
+{
+    return psiH(freesteam_set_ph(p,h));
+}
+
+
+//CL: psiH=(drho/dp)_h=const
+Foam::scalar Foam::psiH(SteamState S)
+{
+    label region;
+    scalar kappa,cp,beta,psiH,rho;
+
+    region=freesteam_region(S);
+
+    if (region==1)
+    {
+        //Cl: note: in FreeStream, beta=1/V*(dV/dP)_P=const is called alphaV (in this region)
+        //Cl: note: in FreeStream, kappa=1/V*(dV/dP)_T=const is called kappaT (in this region)
+        kappa=freesteam_region1_kappaT_pT(S.R1.p,S.R1.T);
+        beta=freesteam_region1_alphav_pT(S.R1.p,S.R1.T);
+        cp=freesteam_region1_cp_pT(S.R1.p,S.R1.T);
+        rho=1/freesteam_region1_v_pT(S.R1.p,S.R1.T);
+
+        //CL: getting derivatives using Bridgmans table
+        //CL: psiH=(drho/dp)_h=const
+        psiH=-((S.R1.T*beta*beta-beta)/cp-kappa*rho);
+    }
+    else if (region==2)
+    {
+        //Cl: note: in FreeStream, beta=1/V*(dV/dP)_P=const is called alphaV (in this region)
+        //Cl: note: in FreeStream, kappa=1/V*(dV/dP)_T=const is called kappaT (in this region)
+        kappa=freesteam_region2_kappaT_pT(S.R2.p,S.R2.T);
+        beta=freesteam_region2_alphav_pT(S.R2.p,S.R2.T);
+        cp=freesteam_region2_cp_pT(S.R2.p,S.R2.T);
+        rho=1/freesteam_region2_v_pT(S.R2.p,S.R2.T);
+
+        //CL: getting derivatives using Bridgmans table
+        //CL: psiH=(drho/dp)_h=const
+        psiH=-((S.R2.T*beta*beta-beta)/cp-kappa*rho);
+    }
+    else if (region==3)
+    {
+
+        scalar gamma,cv,p;
+
+        rho=S.R3.rho;
+        p=freesteam_region3_p_rhoT(S.R3.rho,S.R3.T);
+
+        //Cl: note: beta=1/V*(dV/dP)_P=const
+        //Cl: note: kappa=1/V*(dV/dP)_T=const
+        //Cl: note: in FreeStream, gamma=1/p*(dp/dT)_v=const is called alphap (in this region
+        gamma=freesteam_region3_alphap_rhoT(S.R3.rho,S.R3.T);
+        cp=freesteam_region3_cp_rhoT(S.R3.rho,S.R3.T);
+        cv=freesteam_region3_cv_rhoT(S.R3.rho,S.R3.T);
+        beta=(cp-cv)/(S.R3.T/S.R3.rho*p*gamma);
+        kappa=(cp-cv)/(S.R3.T/S.R3.rho*p*p*gamma*gamma);
+
+        //CL: getting derivatives using Bridgmans table
+        //CL: psiH=(drho/dp)_h=const
+        psiH=-((S.R3.T*beta*beta-beta)/cp-kappa*rho);
+
+    }
+    else if (region==4)
+    {
+        scalar rhov,rhol,betav,betal,kappav,kappal,vv,vl,cpl,hl,hv,h,p;
+        scalar dvldp,dvvdp,dhldp,dhvdp;
+        scalar dpdT,dvdp,dxdp;
+
+        SteamState Sl,Sv;
+
+        rho=1/freesteam_region4_v_Tx(S.R4.T,S.R4.x);
+        h=freesteam_region4_h_Tx(S.R4.T,S.R4.x);
+        p=freesteam_region4_psat_T(S.R4.T);
+
+        //CL: Getting density on the vapour and liquid lines
+        rhov=freesteam_region4_rhog_T(S.R4.T);
+        rhol=freesteam_region4_rhof_T(S.R4.T);
+        vv=1/rhov;
+        vl=1/rhol;
+
+        //CL: getting derivatives --> this is a bit tricky in the vapor dome
+
+        dpdT=freesteam_region4_dpsatdT_T(S.R4.T);
+
+        // getting the states outside the vapour dome
+        Sl=freesteam_set_pv(p,vl-0.0000001);  //inside region 1
+        Sv=freesteam_set_pv(p,vv+0.0000001);  //inside region 2
+
+        kappal=freesteam_region1_kappaT_pT(Sl.R1.p,Sl.R1.T);
+        kappav=freesteam_region2_kappaT_pT(Sv.R2.p,Sv.R2.T);
+
+        betal=freesteam_region1_alphav_pT(Sl.R1.p,Sl.R1.T);
+        betav=freesteam_region2_alphav_pT(Sv.R2.p,Sv.R2.T);
+
+        cpl=freesteam_region1_cp_pT(Sl.R1.p,Sl.R1.T);
+        //cpv=freesteam_region2_cp_pT(Sv.R2.p,Sv.R2.T);
+
+        hl=freesteam_region1_h_pT(Sl.R1.p,Sl.R1.T);
+        hv=freesteam_region2_h_pT(Sv.R2.p,Sv.R2.T);
+
+        //calculation derviatives on liquid and vapour line
+        dvldp=betal*vl/dpdT-kappal*vl;
+        dvvdp=betav*vv/dpdT-kappav*vv;
+
+        dhldp=vl*(1-betal*Sl.R1.T)+cpl/dpdT;
+        dhvdp=vv*(1-betav*Sv.R2.T)+cpl/dpdT;
+
+        dxdp=-dhldp/(hv-hl)
+                 +(h-hl)/((hv-hl)*(hv-hl))
+                     *(dhvdp-dhldp);
+
+        //CL: psiH=(drho/dp)_h=const
+        dvdp=dvldp+(dvvdp-dvldp)*S.R4.x+(vv-vl)*dxdp;
+        psiH=-rho*rho*dvdp;
+    }
+    else
+    {
+        Info<<"IAPWS-IF97.C error, outside the regions 1-4"<<endl;
+
+        //Keep compiler happy
+        psiH = 0;
+    }
+
+    return psiH;
+}
+
+
+//CL: drhodh=(drho/dh)_p=const
+Foam::scalar Foam::drhodh_pT(scalar p,scalar T)
+{
+    return drhodh(freesteam_set_pT(p,T));
+}
+
+
+//CL: drhodh=(drho/dh)_p=const
+Foam::scalar Foam::drhodh_ph(scalar p,scalar h)
+{
+    return drhodh(freesteam_set_ph(p,h));
+}
+
+
+//CL: drhodh=(drho/dh)_p=const
+Foam::scalar Foam::drhodh(SteamState S)
+{
+    label region;
+    scalar cp,beta,drhodh,rho;
+
+    region=freesteam_region(S);
+
+    if (region==1)
+    {
+        rho=1/freesteam_region1_v_pT(S.R1.p,S.R1.T);
+
+        //Cl: note: in FreeStream, beta=1/V*(dV/dP)_P=const is called alphaV (in this region)
+        beta=freesteam_region1_alphav_pT(S.R1.p,S.R1.T);
+        cp=freesteam_region1_cp_pT(S.R1.p,S.R1.T);
+
+        //CL: getting derivatives using Bridgmans table
+        //CL: drhodh=(drho/dh)_p=const
+        drhodh=-rho*beta/cp;
+    }
+    else if (region==2)
+    {
+        rho=1/freesteam_region2_v_pT(S.R2.p,S.R2.T);
+
+        //Cl: note: in FreeStream, beta=1/V*(dV/dP)_P=const is called alphaV (in this region)
+        //Cl: note: in FreeStream, kappa=1/V*(dV/dP)_T=const is called kappaT (in this region)
+        beta=freesteam_region2_alphav_pT(S.R2.p,S.R2.T);
+        cp=freesteam_region2_cp_pT(S.R2.p,S.R2.T);
+
+        //CL: getting derivatives using Bridgmans table
+        //CL: drhodh=(drho/dh)_p=const
+        drhodh=-rho*beta/cp;
+    }
+    else if (region==3)
+    {
+        scalar gamma,cv,p;
+
+        p=freesteam_region3_p_rhoT(S.R3.rho,S.R3.T);
+
+        //Cl: note: beta=1/V*(dV/dP)_P=const
+        //Cl: note: in FreeStream, gamma=1/p*(dp/dT)_v=const is called alphap (in this region
+        gamma=freesteam_region3_alphap_rhoT(S.R3.rho,S.R3.T);
+        cp=freesteam_region3_cp_rhoT(S.R3.rho,S.R3.T);
+        cv=freesteam_region3_cv_rhoT(S.R3.rho,S.R3.T);
+        beta=(cp-cv)/(S.R3.T/S.R3.rho*p*gamma);
+
+        //CL: getting derivatives using Bridgmans table
+        //CL: drhodh=(drho/dh)_p=const
+        drhodh=-S.R3.rho*beta/cp;
+    }
+    else if (region==4)
+    {
+        scalar vv,vl,hl,hv,p;
+        SteamState Sl,Sv;
+
+        rho=1/freesteam_region4_v_Tx(S.R4.T,S.R4.x);
+        p=freesteam_region4_psat_T(S.R4.T);
+
+        //CL: Getting density on the vapour and liquid lines
+        vv=1/freesteam_region4_rhog_T(S.R4.T);
+        vl=1/freesteam_region4_rhof_T(S.R4.T);
+
+        // getting the states outside the vapour dome
+        Sl=freesteam_set_pv(p,vl-0.0000001);  //inside region 1
+        Sv=freesteam_set_pv(p,vv+0.0000001);  //inside region 2
+
+        hl=freesteam_region1_h_pT(Sl.R1.p,Sl.R1.T);
+        hv=freesteam_region2_h_pT(Sv.R2.p,Sv.R2.T);
+
+        //CL: drhodh=(drho/dh)_p=const
+        drhodh=-rho*rho*(vv-vl)/(hv-hl);
+    }
+    else
+    {
+         Info<<"IAPWS-IF97.C error, outside the regions 1-4"<<endl;
+         
+         // Keep compiler happy
+         drhodh = 0;
+    }
+
+    return drhodh;
+}
+
+

src/thermophysicalModels/basic/IAPWS_Waterproperties/IAPWSThermo/IAPWS-IF97.H

@@ -0,0 +1,154 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright held by original author
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+Description
+    IAPWS-IF97 (water) based thermodynamic class. Water properties calculated by freeSteam.
+
+    This code connects OpenFoam with freeSteam and provides the basic functions needed in OpenFOAM
+
+    For more information about freeSteam and its authors have a look @ http://freesteam.sourceforge.net/example.php
+
+SourceFiles
+  IAPWS-IF97.C
+
+Author
+Christian Lucas
+Institut für Thermodynamik
+Technische Universität Braunschweig
+Germany
+
+\*---------------------------------------------------------------------------*/
+#ifndef IAPWSIF97_H
+#define IAPWSIF97_H
+
+#include "basicPsiThermo.H"
+#include "steam.H"
+
+namespace Foam
+{
+    //CL: Functions to caluculate all fluid properties
+    void calculateProperties_h
+    (
+        SteamState S,
+        scalar &rho,
+        scalar &h,
+        scalar &T,
+        scalar &p,
+        scalar &psi,
+        scalar &drhodh,
+        scalar &mu,
+        scalar &alpha,
+        scalar &x
+    );
+
+    //CL: This functions returns all (minimal) needed propeties (p,T,h,rho,psi,drhodh,mu and alpha) for given p and T
+    void calculateProperties_pT
+    (
+        scalar &p,
+        scalar &T,
+        scalar &h,
+        scalar &rho,
+        scalar &psi,
+        scalar &drhodh,
+        scalar &mu,
+        scalar &alpha
+    );
+
+    //CL: This function returns the same values as the function above for given p and T
+    //CL: Additionally, the vapor mass fraction x is return
+    //CL: NOTE: This function is only included to have the possibility to update x at the fixedValue (Temperature) BC
+    //CL: can only return x=0 and x=1 because it is not possible to describe the vapour dome with p and T
+    void calculateProperties_pT
+    (
+        scalar &p,
+        scalar &T,
+        scalar &h,
+        scalar &rho,
+        scalar &psi,
+        scalar &drhodh,
+        scalar &mu,
+        scalar &alpha,
+        scalar &x
+    );
+
+    //CL: This functions returns all (minimal) needed properties (p,T,h,rho,psi,drhodh,mu and alpha) for given p and h
+    void calculateProperties_ph
+    (
+        scalar &p,
+        scalar &h,
+        scalar &T,
+        scalar &rho,
+        scalar &psi,
+        scalar &drhodh,
+        scalar &mu,
+        scalar &alpha
+    );
+
+    //CL: This function returns the same values as the function above for given p and h
+    //CL: Additionally, the vapor mass fraction x is return
+    void calculateProperties_ph
+    (
+        scalar &p,
+        scalar &h,
+        scalar &T,
+        scalar &rho,
+        scalar &psi,
+        scalar &drhodh,
+        scalar &mu,
+        scalar &alpha,
+        scalar &x
+    );
+
+
+    //CL: Return density for given pT or ph;
+    scalar rho_pT(scalar p,scalar T);
+    scalar rho_ph(scalar p,scalar h);
+
+    //CL: Return cp for given pT or ph;
+    scalar cp_pT(scalar p,scalar T);
+    scalar cp_ph(scalar p,scalar h);
+
+    //CL: Return cv for given pT or ph;
+    scalar cv_pT(scalar p,scalar T);
+    scalar cv_ph(scalar p,scalar h);
+
+    //CL: Return enthalpy for given pT;
+    scalar h_pT(scalar p,scalar T);
+
+    //CL: Return temperature for given ph;
+    scalar T_ph(scalar p,scalar T);
+
+    //CL: Return psiH=(drho/dp)_h=constant for given pT or ph;
+    scalar psiH_pT(scalar p,scalar T);
+    scalar psiH_ph(scalar p,scalar h);
+    scalar psiH(SteamState S);
+
+    //CL: Return drhodh=(drho/dh)_p=constant for given pT or ph;
+    scalar drhodh_pT(scalar p,scalar T);
+    scalar drhodh_ph(scalar p,scalar h);
+    scalar drhodh(SteamState S);
+
+}
+
+#endif //IAPWSIF97_C_

src/thermophysicalModels/basic/IAPWS_Waterproperties/IAPWSThermo/IAPWSThermo.C

@@ -0,0 +1,477 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright held by original author
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+Author
+Christian Lucas
+Institut für Thermodynamik
+Technische Universität Braunschweig
+Germany
+
+\*---------------------------------------------------------------------------*/
+
+#include "IAPWSThermo.H"
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+
+void Foam::IAPWSThermo::calculate()
+{
+    scalarField& hCells = h_.internalField();
+    scalarField& pCells = this->p_.internalField();
+    scalarField& TCells = this->T_.internalField();
+    scalarField& rhoCells= this->rho_.internalField();
+    scalarField& psiCells = this->psi_.internalField();
+    scalarField& drhodhCells = this->drhodh_.internalField();
+    scalarField& muCells = this->mu_.internalField();
+    scalarField& alphaCells = this->alpha_.internalField();
+
+    //CL: Updating all cell properties
+    //CL: loop through all cells
+    forAll(TCells, celli)
+    {
+        //CL: see IAPWAS-IF97.H
+        calculateProperties_ph
+        (
+            pCells[celli],
+            hCells[celli],
+            TCells[celli],
+            rhoCells[celli],
+            psiCells[celli],
+            drhodhCells[celli],
+            muCells[celli],
+            alphaCells[celli]
+        );
+    }
+
+    //CL: loop through all patches
+    forAll(T_.boundaryField(), patchi)
+    {
+        fvPatchScalarField& pp = this->p_.boundaryField()[patchi];
+        fvPatchScalarField& pT = this->T_.boundaryField()[patchi];
+        fvPatchScalarField& ppsi = this->psi_.boundaryField()[patchi];
+        fvPatchScalarField& pdrhodh = this->drhodh_.boundaryField()[patchi];
+        fvPatchScalarField& prho = this->rho_.boundaryField()[patchi];
+        fvPatchScalarField& ph = this->h_.boundaryField()[patchi];
+        fvPatchScalarField& pmu = this->mu_.boundaryField()[patchi];
+        fvPatchScalarField& palpha = this->alpha_.boundaryField()[patchi];
+
+        //CL: Updating the patch properties for patches with fixed temperature BC's
+        if (pT.fixesValue())
+        {
+            forAll(pT, facei)
+            {
+                //CL: see IAPWAS-IF97.H
+                calculateProperties_pT
+                (
+                    pp[facei],
+                    pT[facei],
+                    ph[facei],
+                    prho[facei],
+                    ppsi[facei],
+                    pdrhodh[facei],
+                    pmu[facei],
+                    palpha[facei]
+                );
+            }
+        }
+        //CL: Updating the patch properties for patches without fixed temperature BC's
+        else
+        {
+            forAll(pT, facei)
+            {
+                //CL: see IAPWAS-IF97.H
+                calculateProperties_ph
+                (
+                    pp[facei],
+                    ph[facei],
+                    pT[facei],
+                    prho[facei],
+                    ppsi[facei],
+                    pdrhodh[facei],
+                    pmu[facei],
+                    palpha[facei]
+                );
+            }
+        }
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+
+Foam::IAPWSThermo::IAPWSThermo
+(
+    const fvMesh& mesh,
+    const objectRegistry& obj
+)
+:
+    basicPsiThermo(mesh, obj),
+
+    h_
+    (
+        IOobject
+        (
+            "h",
+            mesh.time().timeName(),
+            obj,
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        ),
+        mesh,
+        dimensionSet(0, 2, -2, 0, 0),
+        this->hBoundaryTypes()
+    ),
+
+    rho_
+    (
+        IOobject
+        (
+            "rhoThermo",
+            mesh.time().timeName(),
+            obj,
+            IOobject::READ_IF_PRESENT,
+            IOobject::NO_WRITE
+        ),
+        mesh,
+        dimDensity
+    ),
+
+    drhodh_
+    (
+        IOobject
+        (
+            "drhodh",
+            mesh.time().timeName(),
+            obj,
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        ),
+        mesh,
+        dimensionSet(1, -5, 2, 0, 0)
+    )
+{
+
+    scalarField& hCells = h_.internalField();
+    scalarField& TCells = this->T_.internalField();
+    scalarField& pCells =this->p_.internalField();
+    scalarField& rhoCells =this->rho_.internalField();
+
+    forAll(hCells, celli)
+    {
+        hCells[celli] = h_pT(pCells[celli],TCells[celli]);
+    }
+
+    forAll(h_.boundaryField(), patchi)
+    {
+        h_.boundaryField()[patchi] ==
+            h(this->T_.boundaryField()[patchi], patchi);
+    }
+
+    forAll(rhoCells, celli)
+    {
+        rhoCells[celli] = rho_pT(pCells[celli],TCells[celli]);
+    }
+
+    forAll(rho_.boundaryField(), patchi)
+    {
+        rho_.boundaryField()[patchi] ==
+            rho(this->p_.boundaryField()[patchi] ,this->h_.boundaryField()[patchi], patchi);
+    }
+
+    hBoundaryCorrection(h_);
+
+    calculate();
+
+    // Switch on saving old time
+    this->psi_.oldTime();
+}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::IAPWSThermo::~IAPWSThermo()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+void Foam::IAPWSThermo::correct()
+{
+    if (debug)
+    {
+        Info<< "entering IAPWSThermo::correct()" << endl;
+    }
+
+    // force the saving of the old-time values
+    this->psi_.oldTime();
+
+    calculate();
+
+    if (debug)
+    {
+        Info<< "exiting IAPWSThermo::correct()" << endl;
+    }
+}
+
+
+Foam::tmp<Foam::scalarField> Foam::IAPWSThermo::h
+(
+    const scalarField& T,
+    const labelList& cells
+) const
+{
+    //getting pressure field
+    const scalarField& pCells = this->p_.internalField();
+
+    tmp<scalarField> th(new scalarField(T.size()));
+    scalarField& h = th();
+
+    forAll(T, celli)
+    {
+        h[celli] =  h_pT(pCells[cells[celli]],T[celli]);
+    }
+
+    return th;
+}
+
+
+Foam::tmp<Foam::scalarField> Foam::IAPWSThermo::h
+(
+    const scalarField& T,
+    const label patchi
+) const
+{
+    // getting pressure at the patch
+    const fvPatchScalarField& pp = this->p_.boundaryField()[patchi];
+
+    tmp<scalarField> th(new scalarField(T.size()));
+    scalarField& h = th();
+
+    forAll(T, facei)
+    {
+        h[facei] =  h_pT(pp[facei], T[facei]);
+    }
+
+    return th;
+}
+
+
+//CL: Calculates rho at patch
+Foam::tmp<Foam::scalarField> Foam::IAPWSThermo::rho
+(
+    const scalarField& p,
+    const scalarField& h,
+    const label patchi
+) const
+{
+    tmp<scalarField> trho(new scalarField(h.size()));
+    scalarField& rho = trho();
+
+    forAll(h, facei)
+    {
+        rho[facei] = rho_ph(p[facei], h[facei]);
+    }
+
+    return trho;
+}
+
+
+Foam::tmp<Foam::scalarField> Foam::IAPWSThermo::Cp
+(
+    const scalarField& T,
+    const label patchi
+) const
+{
+    // getting pressure at the patch
+    const fvPatchScalarField& pp = this->p_.boundaryField()[patchi];
+
+    tmp<scalarField> tCp(new scalarField(T.size()));
+    scalarField& cp = tCp();
+
+    forAll(T, facei)
+    {
+        cp[facei] = cp_ph(pp[facei],h_pT(pp[facei], T[facei]));
+    }
+
+    return tCp;
+}
+
+
+Foam::tmp<Foam::volScalarField> Foam::IAPWSThermo::Cp() const
+{
+    const fvMesh& mesh = this->T_.mesh();
+
+    tmp<volScalarField> tCp
+    (
+        new volScalarField
+        (
+            IOobject
+            (
+                "Cp",
+                mesh.time().timeName(),
+                this->T_.db(),
+                IOobject::NO_READ,
+                IOobject::NO_WRITE
+            ),
+            mesh,
+            dimensionSet(0, 2, -2, -1, 0)
+        )
+    );
+
+    volScalarField& cp = tCp();
+
+    forAll(this->T_, celli)
+    {
+        cp[celli] = cp_ph(this->p_[celli], this->h_[celli]);
+    }
+
+    forAll(this->T_.boundaryField(), patchi)
+    {
+        const fvPatchScalarField& pp = this->p_.boundaryField()[patchi];
+        const fvPatchScalarField& ph = this->h_.boundaryField()[patchi];
+        fvPatchScalarField& pCp = cp.boundaryField()[patchi];
+
+        forAll(ph, facei)
+        {
+            pCp[facei] = cp_ph(pp[facei], ph[facei]);
+        }
+    }
+
+    return tCp;
+}
+
+
+//CL: Returns an updated field for rho
+Foam::tmp<Foam::volScalarField> Foam::IAPWSThermo::rho() const
+{
+    const fvMesh& mesh = this->p_.mesh();
+
+    tmp<volScalarField> prho
+    (
+        new volScalarField
+        (
+            IOobject
+            (
+                "rhoThermo2",
+                mesh.time().timeName(),
+                this->T_.db(),
+                IOobject::NO_READ,
+                IOobject::NO_WRITE
+            ),
+            mesh,
+            dimDensity
+        )
+    );
+
+    volScalarField& rho = prho();
+
+    forAll(this->p_, celli)
+    {
+        rho[celli] = rho_ph(this->p_[celli], this->h_[celli]);
+    }
+
+    forAll(this->T_.boundaryField(), patchi)
+    {
+        const fvPatchScalarField& pp = this->p_.boundaryField()[patchi];
+        const fvPatchScalarField& ph = this->h_.boundaryField()[patchi];
+        fvPatchScalarField& prho = rho.boundaryField()[patchi];
+
+        forAll(ph, facei)
+        {
+            prho[facei] = rho_ph(pp[facei], ph[facei]);
+        }
+    }
+
+    return prho;
+}
+
+
+Foam::tmp<Foam::scalarField> Foam::IAPWSThermo::Cv
+(
+    const scalarField& T,
+    const label patchi
+) const
+{
+    // getting pressure at the patch
+    const fvPatchScalarField& pp = this->p_.boundaryField()[patchi];
+
+    tmp<scalarField> tCv(new scalarField(T.size()));
+    scalarField& cv = tCv();
+
+    forAll(T, facei)
+    {
+        cv[facei] = cv_ph(pp[facei], h_pT(pp[facei], T[facei]));
+    }
+
+    return tCv;
+}
+
+
+Foam::tmp<Foam::volScalarField> Foam::IAPWSThermo::Cv() const
+{
+    const fvMesh& mesh = this->T_.mesh();
+
+    tmp<volScalarField> tCv
+    (
+        new volScalarField
+        (
+            IOobject
+            (
+                "Cv",
+                mesh.time().timeName(),
+                this->T_.db(),
+                IOobject::NO_READ,
+                IOobject::NO_WRITE
+            ),
+            mesh,
+            dimensionSet(0, 2, -2, -1, 0)
+        )
+    );
+
+    volScalarField& cv = tCv();
+
+    forAll(this->h_, celli)
+    {
+        cv[celli] = cv_ph(this->p_[celli], this->h_[celli]);
+    }
+
+    forAll(this->T_.boundaryField(), patchi)
+    {
+        cv.boundaryField()[patchi] =
+        Cv(this->T_.boundaryField()[patchi], patchi);
+    }
+
+    return tCv;
+}
+
+
+bool Foam::IAPWSThermo::read()
+{
+    basicPsiThermo::read();
+    return true;
+}
+
+
+// ************************************************************************* //

src/thermophysicalModels/basic/IAPWS_Waterproperties/IAPWSThermo/IAPWSThermo.H

@@ -0,0 +1,201 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright held by original author
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+Class
+    Foam::IAPWSThermo
+
+Description:
+
+    Waterproperties based on the IAPWS 97 tables
+    The water properties are caluclated using freeSteam (http://freesteam.sourceforge.net/example.php)
+
+    General paper decribing the water tables:
+
+    "Revised Release on the IAPWS Industrial Formulation 1997 for the Thermodynamic Properties of Water and Steam"
+
+SourceFiles
+    IAPWSThermo.C
+
+Author
+Christian Lucas
+Institut für Thermodynamik
+Technische Universität Braunschweig
+Germany
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef IAPWSThermo_H
+#define IAPWSThermo_H
+
+#include "basicPsiThermo.H"
+#include "IAPWS-IF97.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+/*---------------------------------------------------------------------------*\
+                         Class IAPWSTHERMO Declaration
+\*---------------------------------------------------------------------------*/
+
+class IAPWSThermo
+:
+    public basicPsiThermo
+{
+    // Private data
+
+        //- Enthalpy field
+        volScalarField h_;
+
+        //- DensityField
+        volScalarField rho_;
+
+        // CL:drhodh Field
+        // CL:needed for pressure equation
+        volScalarField drhodh_;
+
+
+    // Private member functions
+
+        //- Calculate the thermo variables
+        void calculate();
+
+        //- Construct as copy (not implemented)
+        IAPWSThermo(const IAPWSThermo&);
+
+
+public:
+
+    //- Runtime type information
+    TypeName("IAPWSThermo");
+
+    // Constructors
+
+        //- Construct from mesh
+        IAPWSThermo(const fvMesh&, const objectRegistry& obj);
+
+
+    //- Destructor
+    virtual ~IAPWSThermo();
+
+    // Member functions
+
+
+
+        //- Update properties
+        virtual void correct();
+
+
+        // Access to thermodynamic state variables
+
+            //- Enthalpy [J/kg]
+            //  Non-const access allowed for transport equations
+            virtual volScalarField& h()
+            {
+                return h_;
+            }
+
+            //- Enthalpy [J/kg]
+            virtual const volScalarField& h() const
+            {
+                return h_;
+            }
+
+        // Fields derived from thermodynamic state variables
+
+            //- Enthalpy for cell-set [J/kg]
+            virtual tmp<scalarField> h
+            (
+                const scalarField& T,
+                const labelList& cells
+            ) const;
+
+
+            //- Enthalpy for patch [J/kg]
+            virtual tmp<scalarField> h
+            (
+                const scalarField& T,
+                const label patchi
+            ) const;
+
+
+            //- Density for patch [J/kg]
+            virtual tmp<scalarField> rho
+            (
+                const scalarField& p,
+                const scalarField& h,
+                const label patchi
+            ) const;
+
+            //- Heat capacity at constant pressure for patch [J/kg/K]
+            // dummy function needed for BC
+            virtual tmp<scalarField> Cp
+            (
+                const scalarField& T,
+                const label patchi
+            ) const;
+
+            //- Heat capacity at constant pressure [J/kg/K]
+            virtual tmp<volScalarField> Cp() const;
+
+            //- Heat capacity at constant volume for patch [J/kg/K]
+            virtual tmp<scalarField> Cv
+            (
+                const scalarField& T,
+                const label patchi
+            ) const;
+
+            //- Heat capacity at constant volume [J/kg/K]
+            virtual tmp<volScalarField> Cv() const;
+
+            //- Gradient drhodh @ constant pressure
+            virtual const volScalarField& drhodh() const
+            {
+                return drhodh_;
+            }
+
+             //- Density [kg/m^3] - uses current value of pressure
+             virtual tmp<volScalarField> rho() const;
+
+
+        //- Read thermophysicalProperties dictionary
+        virtual bool read();
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+#   include "IAPWSThermo.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/thermophysicalModels/basic/IAPWS_Waterproperties/IAPWSThermo/IAPWSThermos.C

@@ -0,0 +1,53 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           |  Copyright held by original author
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+Author
+Christian Lucas
+Institut für Thermodynamik
+Technische Universität Braunschweig
+Germany
+
+\*---------------------------------------------------------------------------*/
+
+#include "makeBasicPsiThermo.H"
+#include "IAPWSThermo.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/* * * * * * * * * * * * * * * private static data * * * * * * * * * * * * * */
+
+makeBasicExternalLibraryBasedThermo
+(
+    IAPWSThermo
+);
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// ************************************************************************* //

src/thermophysicalModels/basic/IAPWS_Waterproperties/IAPWSThermo/steam.H

@@ -0,0 +1,135 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright held by original author
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+
+Description:
+
+Includes the definition of a struct from freestream
+Code copied from freestream --> steam.h
+
+\*---------------------------------------------------------------------------*/
+
+#ifdef __cplusplus
+    #define EXTERN extern "C"
+#else
+    #define EXTERN extern
+#endif
+
+typedef struct SteamState_R1_struct{
+    double p, T;
+} SteamState_R1;
+
+typedef struct SteamState_R2_struct{
+    double p, T;
+} SteamState_R2;
+
+typedef struct SteamState_R3_struct{
+    double rho, T;
+} SteamState_R3;
+
+typedef struct SteamState_R4_struct{
+    double T, x;
+} SteamState_R4;
+
+typedef struct SteamState_struct{
+    char region;
+    union{
+        SteamState_R1 R1;
+        SteamState_R2 R2;
+        SteamState_R3 R3;
+        SteamState_R4 R4;
+    };
+} SteamState;
+
+//CL:
+EXTERN double freesteam_p(SteamState S);
+EXTERN double freesteam_T(SteamState S);
+EXTERN double freesteam_rho(SteamState S);
+EXTERN double freesteam_v(SteamState S);
+EXTERN double freesteam_u(SteamState S);
+EXTERN double freesteam_h(SteamState S);
+EXTERN double freesteam_s(SteamState S);
+EXTERN double freesteam_cp(SteamState S);
+EXTERN double freesteam_cv(SteamState S);
+EXTERN double freesteam_w(SteamState S);
+EXTERN double freesteam_x(SteamState S);
+EXTERN double freesteam_mu(SteamState S);
+EXTERN double freesteam_k(SteamState S);
+
+//CL: getting SteamState for two given properties e.g. pressure and temperatur
+EXTERN SteamState freesteam_set_pv(double,double);
+EXTERN SteamState freesteam_set_pu(double,double);
+EXTERN SteamState freesteam_set_pT(double,double);
+EXTERN SteamState freesteam_set_ph(double,double);
+
+//CL: getting region of the SteamState
+EXTERN int freesteam_region(SteamState);
+
+//CL: transport properties
+EXTERN double freesteam_mu_rhoT(double,double);
+EXTERN double freesteam_k_rhoT(double,double);
+
+//CL: Region 1 --> see region1.h (freesteam)
+EXTERN double freesteam_region1_v_pT(double,double);
+EXTERN double freesteam_region1_h_pT(double,double);
+EXTERN double freesteam_region1_kappaT_pT(double,double);
+EXTERN double freesteam_region1_alphav_pT(double,double);
+EXTERN double freesteam_region1_cp_pT(double,double);
+EXTERN double freesteam_region1_u_pT(double,double);
+EXTERN double freesteam_region1_s_pT(double,double);
+EXTERN double freesteam_region1_cv_pT(double,double);
+
+//CL: Region 2 --> see region2.h (freesteam)
+EXTERN double freesteam_region2_v_pT(double,double);
+EXTERN double freesteam_region2_u_pT(double,double);
+EXTERN double freesteam_region2_s_pT(double,double);
+EXTERN double freesteam_region2_h_pT(double,double);
+EXTERN double freesteam_region2_cp_pT(double,double);
+EXTERN double freesteam_region2_cv_pT(double,double);
+EXTERN double freesteam_region2_alphav_pT(double,double);
+EXTERN double freesteam_region2_kappaT_pT(double,double);
+
+//CL: Region 3 --> see region3.h (freesteam)
+EXTERN double freesteam_region3_p_rhoT(double,double);
+EXTERN double freesteam_region3_u_rhoT(double,double);
+EXTERN double freesteam_region3_s_rhoT(double,double);
+EXTERN double freesteam_region3_h_rhoT(double,double);
+EXTERN double freesteam_region3_cp_rhoT(double,double);
+EXTERN double freesteam_region3_cv_rhoT(double,double);
+EXTERN double freesteam_region3_alphap_rhoT(double,double);
+EXTERN double freesteam_region3_betap_rhoT(double,double);
+
+//CL: Region 4 --> see region4.h (freesteam)
+EXTERN double freesteam_region4_psat_T(double);
+EXTERN double freesteam_region4_Tsat_p(double);
+EXTERN double freesteam_region4_rhof_T(double);
+EXTERN double freesteam_region4_rhog_T(double);
+EXTERN double freesteam_region4_v_Tx(double,double);
+EXTERN double freesteam_region4_u_Tx(double,double);
+EXTERN double freesteam_region4_h_Tx(double,double);
+EXTERN double freesteam_region4_s_Tx(double,double);
+EXTERN double freesteam_region4_cp_Tx(double,double);
+EXTERN double freesteam_region4_cv_Tx(double,double);
+EXTERN double freesteam_region4_dpsatdT_T(double);
+

src/thermophysicalModels/basic/Make/files

@@ -8,12 +8,17 @@ psiThermo/basicPsiThermo/newBasicPsiThermo.C
 psiThermo/hPsiThermo/hPsiThermos.C
 psiThermo/hsPsiThermo/hsPsiThermos.C
 psiThermo/ePsiThermo/ePsiThermos.C
+psiThermo/realGasHThermo/realGasHThermos.C
+psiThermo/realGasEThermo/realGasEThermos.C
 
 rhoThermo/basicRhoThermo/basicRhoThermo.C
 rhoThermo/basicRhoThermo/newBasicRhoThermo.C
 rhoThermo/hRhoThermo/hRhoThermos.C
 rhoThermo/hsRhoThermo/hsRhoThermos.C
 
+IAPWS_Waterproperties/IAPWSThermo/IAPWS-IF97.C
+IAPWS_Waterproperties/IAPWSThermo/IAPWSThermos.C
+
 derivedFvPatchFields/fixedEnthalpy/fixedEnthalpyFvPatchScalarField.C
 derivedFvPatchFields/gradientEnthalpy/gradientEnthalpyFvPatchScalarField.C
 derivedFvPatchFields/mixedEnthalpy/mixedEnthalpyFvPatchScalarField.C

src/thermophysicalModels/basic/Make/options

@@ -1,7 +1,9 @@
 EXE_INC = \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
-    -I$(LIB_SRC)/thermophysicalModels/specie/lnInclude
+    -I$(LIB_SRC)/thermophysicalModels/specie/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/thermophysicalFunctions/lnInclude
 
 LIB_LIBS = \
     -lspecie \
+    -lthermophysicalFunctions \
     -lfiniteVolume

src/thermophysicalModels/basic/mixtures/basicMixture/basicMixtures.C

@@ -1,25 +1,26 @@
 /*---------------------------------------------------------------------------*\
   =========                 |
-  \\      /  F ield         | foam-extend: Open Source CFD
-   \\    /   O peration     | Version:     3.2
-    \\  /    A nd           | Web:         http://www.foam-extend.org
-     \\/     M anipulation  | For copyright notice see file Copyright
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright held by original author
+     \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
-    This file is part of foam-extend.
+    This file is part of OpenFOAM.
 
-    foam-extend is free software: you can redistribute it and/or modify it
+    OpenFOAM is free software; you can redistribute it and/or modify it
     under the terms of the GNU General Public License as published by the
-    Free Software Foundation, either version 3 of the License, or (at your
+    Free Software Foundation; either version 2 of the License, or (at your
     option) any later version.
 
-    foam-extend is distributed in the hope that it will be useful, but
-    WITHOUT ANY WARRANTY; without even the implied warranty of
-    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
-    General Public License for more details.
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
 
     You should have received a copy of the GNU General Public License
-    along with foam-extend.  If not, see <http://www.gnu.org/licenses/>.
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
 
 Description
     Mixture instantiation
@@ -32,12 +33,19 @@ Description
 #include "makeBasicMixture.H"
 
 #include "perfectGas.H"
+#include "redlichKwong.H"
+#include "pengRobinson.H"
+#include "aungierRedlichKwong.H"
+#include "soaveRedlichKwong.H"
 
 #include "eConstThermo.H"
 
 #include "hConstThermo.H"
 #include "janafThermo.H"
+#include "nasaHeatCapacityPolynomial.H"
+#include "constantHeatCapacity.H"
 #include "specieThermo.H"
+#include "realGasSpecieThermo.H"
 
 #include "constTransport.H"
 #include "sutherlandTransport.H"
@@ -101,6 +109,150 @@ makeBasicMixturePhys
     icoPoly8ThermoPhysics
 );
 
+makeBasicRealFluidMixture
+(
+    pureMixture,
+    sutherlandTransport,
+    realGasSpecieThermo,
+    nasaHeatCapacityPolynomial,
+    redlichKwong
+);
+
+makeBasicRealFluidMixture
+(
+    pureMixture,
+    sutherlandTransport,
+    realGasSpecieThermo,
+    nasaHeatCapacityPolynomial,
+    pengRobinson
+);
+
+makeBasicRealFluidMixture
+(
+    pureMixture,
+    sutherlandTransport,
+    realGasSpecieThermo,
+    nasaHeatCapacityPolynomial,
+    aungierRedlichKwong
+);
+
+makeBasicRealFluidMixture
+(
+    pureMixture,
+    sutherlandTransport,
+    realGasSpecieThermo,
+    nasaHeatCapacityPolynomial,
+    soaveRedlichKwong
+);
+
+makeBasicRealFluidMixture
+(
+    pureMixture,
+    constTransport,
+    realGasSpecieThermo,
+    nasaHeatCapacityPolynomial,
+    redlichKwong
+);
+
+makeBasicRealFluidMixture
+(
+    pureMixture,
+    constTransport,
+    realGasSpecieThermo,
+    nasaHeatCapacityPolynomial,
+    pengRobinson
+);
+
+makeBasicRealFluidMixture
+(
+    pureMixture,
+    constTransport,
+    realGasSpecieThermo,
+    nasaHeatCapacityPolynomial,
+    aungierRedlichKwong
+);
+
+makeBasicRealFluidMixture
+(
+    pureMixture,
+    constTransport,
+    realGasSpecieThermo,
+    nasaHeatCapacityPolynomial,
+    soaveRedlichKwong
+);
+
+makeBasicRealFluidMixture
+(
+    pureMixture,
+    sutherlandTransport,
+    realGasSpecieThermo,
+    constantHeatCapacity,
+    redlichKwong
+);
+
+makeBasicRealFluidMixture
+(
+    pureMixture,
+    sutherlandTransport,
+    realGasSpecieThermo,
+    constantHeatCapacity,
+    pengRobinson
+);
+
+makeBasicRealFluidMixture
+(
+    pureMixture,
+    sutherlandTransport,
+    realGasSpecieThermo,
+    constantHeatCapacity,
+    aungierRedlichKwong
+);
+
+makeBasicRealFluidMixture
+(
+    pureMixture,
+    sutherlandTransport,
+    realGasSpecieThermo,
+    constantHeatCapacity,
+    soaveRedlichKwong
+);
+
+makeBasicRealFluidMixture
+(
+    pureMixture,
+    constTransport,
+    realGasSpecieThermo,
+    constantHeatCapacity,
+    redlichKwong
+);
+
+makeBasicRealFluidMixture
+(
+    pureMixture,
+    constTransport,
+    realGasSpecieThermo,
+    constantHeatCapacity,
+    pengRobinson
+);
+
+makeBasicRealFluidMixture
+(
+    pureMixture,
+    constTransport,
+    realGasSpecieThermo,
+    constantHeatCapacity,
+    aungierRedlichKwong
+);
+
+makeBasicRealFluidMixture
+(
+    pureMixture,
+    constTransport,
+    realGasSpecieThermo,
+    constantHeatCapacity,
+    soaveRedlichKwong
+);
+
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 

src/thermophysicalModels/basic/mixtures/basicMixture/makeBasicMixture.H

@@ -56,6 +56,26 @@ defineTemplateTypeNameAndDebugWithName                                        \
     0                                                                         \
 );
 
+#define makeBasicRealFluidMixture(Mixture,Transport,SpecieThermo,Thermo,EqnOfState) \
+                                                                              \
+typedef Mixture<Transport<SpecieThermo<Thermo<EqnOfState> > > >               \
+    Mixture##Transport##SpecieThermo##Thermo##EqnOfState;                     \
+                                                                              \
+defineTemplateTypeNameAndDebugWithName                                        \
+    (Mixture##Transport##SpecieThermo##Thermo##EqnOfState,                    \
+    #Mixture"<"#Transport"<"#SpecieThermo"<"#Thermo"<"#EqnOfState">>>>", 0)
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#define makeBasicRealFluidMixture(Mixture,Transport,SpecieThermo,Thermo,EqnOfState)                 \
+                                                                              \
+typedef Mixture<Transport<SpecieThermo<Thermo<EqnOfState> > > >               \
+    Mixture##Transport##SpecieThermo##Thermo##EqnOfState;                     \
+                                                                              \
+defineTemplateTypeNameAndDebugWithName                                        \
+    (Mixture##Transport##SpecieThermo##Thermo##EqnOfState,                    \
+    #Mixture"<"#Transport"<"#SpecieThermo"<"#Thermo"<"#EqnOfState">>>>", 0)
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 

src/thermophysicalModels/basic/psiThermo/basicPsiThermo/basicPsiThermo.C

@@ -1,9 +1,9 @@
 /*---------------------------------------------------------------------------*\
   =========                 |
   \\      /  F ield         | foam-extend: Open Source CFD
-   \\    /   O peration     | Version:     3.2
-    \\  /    A nd           | Web:         http://www.foam-extend.org
-     \\/     M anipulation  | For copyright notice see file Copyright
+   \\    /   O peration     |
+    \\  /    A nd           | For copyright notice see file Copyright
+     \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
     This file is part of foam-extend.
@@ -48,4 +48,19 @@ Foam::basicPsiThermo::~basicPsiThermo()
 {}
 
 
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+const Foam::volScalarField& Foam::basicPsiThermo::drhodh() const
+{
+    notImplemented("basicPsiThermo::drhodh()");
+    return const_cast<volScalarField&>(volScalarField::null());
+}
+
+const Foam::volScalarField& Foam::basicPsiThermo::drhode() const
+{
+    notImplemented("basicPsiThermo::drhode()");
+    return const_cast<volScalarField&>(volScalarField::null());
+}
+
+
 // ************************************************************************* //

src/thermophysicalModels/basic/psiThermo/basicPsiThermo/basicPsiThermo.H

@@ -108,6 +108,12 @@ public:
             {
                 return p_*psi();
             }
+
+            //CL: drhodh needed for pressure equation of the real gas solver
+            virtual const volScalarField& drhodh() const;
+
+            //CL: drhode needed for pressure equation of the real gas solver
+            virtual const volScalarField& drhode() const;
 };
 
 

src/thermophysicalModels/basic/psiThermo/basicPsiThermo/makeBasicPsiThermo.H

@@ -58,6 +58,46 @@ addToRunTimeSelectionTable                                                    \
 )
 
 
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#define makeBasicRealGasThermo(Cthermo,Mixture,Transport,SpecieThermo,Thermo,EqnOfState) \
+                                                                              \
+typedef Cthermo<Mixture<Transport<SpecieThermo<Thermo<EqnOfState> > > > >     \
+    Cthermo##Mixture##Transport##SpecieThermo##Thermo##EqnOfState;            \
+                                                                              \
+defineTemplateTypeNameAndDebugWithName                                        \
+(                                                                             \
+    Cthermo##Mixture##Transport##SpecieThermo##Thermo##EqnOfState,            \
+    #Cthermo                                                                  \
+        "<"#Mixture"<"#Transport"<"#SpecieThermo"<"#Thermo"<"#EqnOfState">>>>>", \
+    0                                                                         \
+);                                                                            \
+                                                                              \
+addToRunTimeSelectionTable                                                    \
+(                                                                             \
+    basicPsiThermo,                                                           \
+    Cthermo##Mixture##Transport##SpecieThermo##Thermo##EqnOfState,            \
+    fvMesh                                                                    \
+)
+
+#define makeBasicExternalLibraryBasedThermo(Cthermo)                          \
+                                                                              \
+typedef Cthermo                                                               \
+    Cthermo;                                                                  \
+                                                                              \
+defineTypeNameAndDebug                                                        \
+(                                                                             \
+    Cthermo,                                                                  \
+    0                                                                         \
+);                                                                            \
+                                                                              \
+addToRunTimeSelectionTable                                                    \
+(                                                                             \
+    basicPsiThermo,                                                           \
+    Cthermo,                                                                  \
+    fvMesh                                                                    \
+)
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 #endif

src/thermophysicalModels/basic/psiThermo/hPsiThermo/hPsiThermo.C

@@ -111,7 +111,7 @@ Foam::hPsiThermo<MixtureType>::hPsiThermo
         (
             "h",
             mesh.time().timeName(),
-            mesh,
+            obj,
             IOobject::NO_READ,
             IOobject::NO_WRITE
         ),

src/thermophysicalModels/basic/psiThermo/realGasEThermo/realGasEThermo.C

@@ -0,0 +1,505 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           |
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+Author
+Christian Lucas
+Institut für Thermodynamik
+Technische Universität Braunschweig
+Germany
+
+\*---------------------------------------------------------------------------*/
+
+#include "realGasEThermo.H"
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+template<class MixtureType>
+void Foam::realGasEThermo<MixtureType>::calculate()
+{
+    const scalarField& eCells = e_.internalField();
+    const scalarField& pCells = this->p_.internalField();
+
+    scalarField& TCells = this->T_.internalField();
+    scalarField& rhoCells= this->rho_.internalField();
+    scalarField& psiCells = this->psi_.internalField();
+    scalarField& drhodeCells = this->drhode_.internalField();
+    scalarField& muCells = this->mu_.internalField();
+    scalarField& alphaCells = this->alpha_.internalField();
+
+    forAll(TCells, celli)
+    {
+        const typename MixtureType::thermoType& mixture_ =
+            this->cellMixture(celli);
+
+        mixture_.TE(eCells[celli], TCells[celli], pCells[celli], rhoCells[celli]);
+        psiCells[celli]=mixture_.psiE(rhoCells[celli], TCells[celli]);
+        drhodeCells[celli]=mixture_.drhodE(rhoCells[celli], TCells[celli]);
+        muCells[celli] = mixture_.mu(TCells[celli]);
+        alphaCells[celli] = mixture_.alpha(rhoCells[celli], TCells[celli]);
+    }
+
+
+    forAll(T_.boundaryField(), patchi)
+    {
+        fvPatchScalarField& pp = this->p_.boundaryField()[patchi];
+        fvPatchScalarField& pT = this->T_.boundaryField()[patchi];
+        fvPatchScalarField& ppsi = this->psi_.boundaryField()[patchi];
+        fvPatchScalarField& pdrhode = this->drhode_.boundaryField()[patchi];
+        fvPatchScalarField& prho = this->rho_.boundaryField()[patchi];
+        fvPatchScalarField& pe = e_.boundaryField()[patchi];
+        fvPatchScalarField& pmu = this->mu_.boundaryField()[patchi];
+        fvPatchScalarField& palpha = this->alpha_.boundaryField()[patchi];
+
+        if (pT.fixesValue())
+        {
+            forAll(pT, facei)
+            {
+                const typename MixtureType::thermoType& mixture_ =
+                    this->patchFaceMixture(patchi, facei);
+
+                prho[facei] = mixture_.rho(pp[facei], pT[facei],prho[facei]);
+                ppsi[facei]=mixture_.psiE(prho[facei],pT[facei]);
+                pdrhode[facei]=mixture_.drhodE(prho[facei],pT[facei]);
+                pe[facei] = mixture_.E(prho[facei], pT[facei]);
+                pmu[facei] = mixture_.mu(pT[facei]);
+                palpha[facei] = mixture_.alpha(prho[facei],pT[facei]);
+            }
+        }
+        else
+        {
+            forAll(pT, facei)
+            {
+                const typename MixtureType::thermoType& mixture_ =
+                    this->patchFaceMixture(patchi, facei);
+                mixture_.TE(pe[facei], pT[facei],pp[facei],prho[facei]);
+                pmu[facei] = mixture_.mu(pT[facei]);
+                ppsi[facei]=mixture_.psiE(prho[facei],pT[facei]);
+                pdrhode[facei]=mixture_.drhodE(prho[facei],pT[facei]);
+                palpha[facei] = mixture_.alpha(prho[facei],pT[facei]);
+            }
+        }
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+template<class MixtureType>
+Foam::realGasEThermo<MixtureType>::realGasEThermo
+(
+    const fvMesh& mesh,
+    const objectRegistry& obj
+)
+:
+    basicPsiThermo(mesh, obj),
+    MixtureType(*this, mesh, obj),
+
+    e_
+    (
+        IOobject
+        (
+            "e",
+            mesh.time().timeName(),
+            obj,
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        ),
+        mesh,
+        dimensionSet(0, 2, -2, 0, 0),
+        this->eBoundaryTypes()
+    ),
+    rho_
+    (
+        IOobject
+        (
+            "rhoThermo",
+            mesh.time().timeName(),
+            obj,
+            IOobject::READ_IF_PRESENT,
+            IOobject::NO_WRITE
+        ),
+        mesh,
+        dimDensity
+    ),
+    drhode_
+    (
+        IOobject
+        (
+            "drhode",
+            mesh.time().timeName(),
+            obj,
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        ),
+        mesh,
+        dimensionSet(1, -5, 2, 0, 0)
+    )
+{
+    scalarField& eCells = e_.internalField();
+    const scalarField& TCells = this->T_.internalField();
+    const scalarField& pCells =this->p_.internalField();
+    scalarField& rhoCells =this->rho_.internalField();
+
+
+    forAll(rhoCells, celli)
+    {
+     rhoCells[celli]=this->cellMixture(celli).rho(pCells[celli],TCells[celli]);
+    }
+
+    forAll(rho_.boundaryField(), patchi)
+    {
+        rho_.boundaryField()[patchi] ==
+            rho(this->T_.boundaryField()[patchi], patchi);
+    }
+
+    forAll(eCells, celli)
+    {
+        eCells[celli] = this->cellMixture(celli).E(rhoCells[celli],TCells[celli]);
+    }
+
+    forAll(e_.boundaryField(), patchi)
+    {
+        e_.boundaryField()[patchi] ==
+            e(this->T_.boundaryField()[patchi], patchi);
+    }
+
+    this->eBoundaryCorrection(e_);
+
+    calculate();
+
+    // Switch on saving old time
+    this->psi_.oldTime();
+}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+template<class MixtureType>
+Foam::realGasEThermo<MixtureType>::~realGasEThermo()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+template<class MixtureType>
+void Foam::realGasEThermo<MixtureType>::correct()
+{
+    if (debug)
+    {
+        Info<< "entering realGasEThermo<MixtureType>::correct()" << endl;
+    }
+
+    // force the saving of the old-time values
+    this->psi_.oldTime();
+
+    calculate();
+
+    if (debug)
+    {
+        Info<< "exiting realGasEThermo<MixtureType>::correct()" << endl;
+    }
+}
+
+
+template<class MixtureType>
+Foam::tmp<Foam::scalarField> Foam::realGasEThermo<MixtureType>::e
+(
+    const scalarField& T,
+    const labelList& cells
+) const
+{
+    //CL: need the pressure of the internal field to calculate the realGas internal energy
+    //CL: this is done this way to assure compatibility to old OF Thermo-Versions
+    const scalarField& pCells = this->p_.internalField();
+
+    tmp<scalarField> te(new scalarField(T.size()));
+    scalarField& e = te();
+
+    forAll(T, celli)
+    {
+        e[celli] = this->cellMixture(cells[celli]).E(this->cellMixture(cells[celli]).rho(pCells[cells[celli]],T[celli]),T[celli]);
+    }
+
+    return te;
+}
+
+template<class MixtureType>
+Foam::tmp<Foam::scalarField> Foam::realGasEThermo<MixtureType>::e
+(
+    const scalarField& T,
+    const label patchi
+) const
+{
+    //CL: need the pressure at the patch to calculate the realGas internal energy
+    //CL: this is done this way to assure compatibility to old OF Thermo-Versions
+    const fvPatchScalarField& pp = this->p_.boundaryField()[patchi];
+
+    tmp<scalarField> te(new scalarField(T.size()));
+    scalarField& e = te();
+
+    forAll(T, facei)
+    {
+        e[facei] = this->patchFaceMixture(patchi, facei).E(this->patchFaceMixture(patchi, facei).rho(pp[facei], T[facei]),T[facei]);
+    }
+
+    return te;
+}
+
+
+template<class MixtureType>
+Foam::tmp<Foam::scalarField> Foam::realGasEThermo<MixtureType>::rho
+(
+    const scalarField& T,
+    const label patchi
+) const
+{
+    //CL: need the pressure at the patch to calculate the realGas density
+    //CL: this is done this way to assure compatibility to old OF Thermo-Versions
+    const fvPatchScalarField& pp = this->p_.boundaryField()[patchi];
+
+    tmp<scalarField> trho(new scalarField(T.size()));
+    scalarField& rho = trho();
+
+    forAll(T, facei)
+    {
+        rho[facei] = this->patchFaceMixture(patchi, facei).rho(pp[facei], T[facei]);
+    }
+
+    return trho;
+}
+
+
+template<class MixtureType>
+Foam::tmp<Foam::scalarField> Foam::realGasEThermo<MixtureType>::Cp
+(
+    const scalarField& T,
+    const label patchi
+) const
+{
+    //CL: need the pressure at the patch to calculate the realGas cp
+    //CL: this is done this way to assure compatibility to old OF Thermo-Versions
+    const fvPatchScalarField& pp = this->p_.boundaryField()[patchi];
+
+    tmp<scalarField> tCp(new scalarField(T.size()));
+    scalarField& cp = tCp();
+
+    forAll(T, facei)
+    {
+        cp[facei] = this->patchFaceMixture(patchi, facei).Cp(this->patchFaceMixture(patchi, facei).rho(pp[facei], T[facei]),T[facei]);
+    }
+
+    return tCp;
+}
+
+
+template<class MixtureType>
+Foam::tmp<Foam::volScalarField> Foam::realGasEThermo<MixtureType>::Cp() const
+{
+    const fvMesh& mesh = this->T_.mesh();
+
+    tmp<volScalarField> tCp
+    (
+        new volScalarField
+        (
+            IOobject
+            (
+                "Cp",
+                mesh.time().timeName(),
+                T_.db(),
+                IOobject::NO_READ,
+                IOobject::NO_WRITE
+            ),
+            mesh,
+            dimensionSet(0, 2, -2, -1, 0)
+        )
+    );
+
+    volScalarField& cp = tCp();
+
+    forAll(this->T_, celli)
+    {
+        cp[celli] = this->cellMixture(celli).Cp(this->rho_[celli], this->T_[celli]);
+    }
+
+    forAll(this->T_.boundaryField(), patchi)
+    {
+        const fvPatchScalarField& pT = this->T_.boundaryField()[patchi];
+        const fvPatchScalarField& prho = this->rho_.boundaryField()[patchi];
+        fvPatchScalarField& pCp = cp.boundaryField()[patchi];
+
+        forAll(pT, facei)
+        {
+            pCp[facei] = this->patchFaceMixture(patchi, facei).Cp(prho[facei], pT[facei]);
+        }
+    }
+
+    return tCp;
+}
+
+
+template<class MixtureType>
+Foam::tmp<Foam::scalarField> Foam::realGasEThermo<MixtureType>::Cv
+(
+    const scalarField& T,
+    const label patchi
+) const
+{
+    //CL: need the pressure at the patch to calculate the realGas internal energy
+    //CL: this is done this way to assure compatibility to old OF Thermo-Versions
+    const fvPatchScalarField& pp = this->p_.boundaryField()[patchi];
+
+    tmp<scalarField> tCv(new scalarField(T.size()));
+    scalarField& cv = tCv();
+
+    forAll(T, facei)
+    {
+        cv[facei] = this->patchFaceMixture(patchi, facei).Cv(this->patchFaceMixture(patchi, facei).rho(pp[facei], T[facei]), T[facei]);
+    }
+
+    return tCv;
+}
+
+
+// CL: Maybe this function should be changed so that it is not "const" function anymore
+template<class MixtureType>
+Foam::tmp<Foam::volScalarField> Foam::realGasEThermo<MixtureType>::rho()  const
+{
+
+    const fvMesh& mesh = this->T_.mesh();
+
+    //CL: create an rho Field, which will be return
+    //CL: the problem is that this function is "const",
+    //CL: so a new variabel is needed
+      tmp<volScalarField> trho
+      (
+        new volScalarField
+        (
+            IOobject
+            (
+                "rhoFunctionThermo",
+                mesh.time().timeName(),
+                T_.db(),
+                IOobject::NO_READ,
+                IOobject::NO_WRITE
+            ),
+            mesh,
+            dimDensity
+        )
+    );
+
+   //CL: copy "old" rho value onto the new rho field as start point
+   //CL: for the newton solver used in this->TE( ... )
+   trho()=rho_;
+
+   volScalarField& rho = trho();
+
+    const scalarField& eCells = e_.internalField();
+    const scalarField& pCells = this->p_.internalField();
+    scalarField TCells = this->T_.internalField();
+
+    forAll(pCells, celli)
+    {
+        const typename MixtureType::thermoType& mixture_ =
+            this->cellMixture(celli);
+
+        // getting the new rho Field
+        mixture_.TE(eCells[celli], TCells[celli], pCells[celli], rho[celli]);
+    }
+
+    forAll(p_.boundaryField(), patchi)
+    {
+        fvPatchScalarField pp = this->p_.boundaryField()[patchi];
+        fvPatchScalarField pe = e_.boundaryField()[patchi];
+        fvPatchScalarField pT = this->T_.boundaryField()[patchi];
+
+        fvPatchScalarField& prho_ = rho.boundaryField()[patchi];
+
+        forAll(pp, facei)
+        {
+            const typename MixtureType::thermoType& mixture_ =
+                this->patchFaceMixture(patchi, facei);
+
+            // getting the new rho patch Field
+            mixture_.TE(pe[facei], pT[facei],pp[facei],prho_[facei]);
+        }
+    }
+    return trho;
+}
+
+
+template<class MixtureType>
+Foam::tmp<Foam::volScalarField> Foam::realGasEThermo<MixtureType>::Cv() const
+{
+    const fvMesh& mesh = this->T_.mesh();
+
+    tmp<volScalarField> tCv
+    (
+        new volScalarField
+        (
+            IOobject
+            (
+                "Cv",
+                mesh.time().timeName(),
+                T_.db(),
+                IOobject::NO_READ,
+                IOobject::NO_WRITE
+            ),
+            mesh,
+            dimensionSet(0, 2, -2, -1, 0)
+        )
+    );
+
+    volScalarField& cv = tCv();
+
+    forAll(this->T_, celli)
+    {
+        cv[celli] = this->cellMixture(celli).Cv(this->rho_[celli], this->T_[celli]);
+    }
+
+    forAll(this->T_.boundaryField(), patchi)
+    {
+        cv.boundaryField()[patchi] =
+            Cv(this->T_.boundaryField()[patchi], patchi);
+    }
+
+    return tCv;
+}
+
+
+template<class MixtureType>
+bool Foam::realGasEThermo<MixtureType>::read()
+{
+    if (basicPsiThermo::read())
+    {
+        MixtureType::read(*this);
+        return true;
+    }
+    else
+    {
+        return false;
+    }
+}
+
+
+// ************************************************************************* //

src/thermophysicalModels/basic/psiThermo/realGasEThermo/realGasEThermo.H

@@ -0,0 +1,205 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | .
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+Class
+    Foam::realGasEThermo
+
+Description
+    Internal energy for a real gas fluid libary
+
+
+SourceFiles
+    realGasEThermo.C
+
+Author
+Christian Lucas
+Institut für Thermodynamik
+Technische Universität Braunschweig
+Germany
+
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef realGasEThermo_H
+#define realGasEThermo_H
+
+#include "basicPsiThermo.H"
+#include "basicMixture.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                         Class realGasEThermo Declaration
+\*---------------------------------------------------------------------------*/
+
+template<class MixtureType>
+class realGasEThermo
+:
+    public basicPsiThermo,
+    public MixtureType
+{
+    // Private data
+
+        //- Enthalpy field
+        volScalarField e_;
+
+        //- Density field
+        volScalarField rho_;
+
+        //- drhode_Field
+        volScalarField drhode_;
+
+    // Private member functions
+
+        //- Calculate the thermo variables
+        void calculate();
+
+        //- Construct as copy (not implemented)
+        realGasEThermo(const realGasEThermo<MixtureType>&);
+
+
+public:
+
+    //- Runtime type information
+    TypeName("realGasEThermo");
+
+
+    // Constructors
+
+        //- Construct from mesh and object registry
+        realGasEThermo(const fvMesh&, const objectRegistry&);
+
+    //- Destructor
+    virtual ~realGasEThermo();
+
+
+    // Member functions
+
+        //- Return the compostion of the mixture
+        virtual basicMixture& composition()
+        {
+            return *this;
+        }
+
+        //- Return the compostion of the mixture
+        virtual const basicMixture& composition() const
+        {
+            return *this;
+        }
+
+        //- Update properties
+        virtual void correct();
+
+
+        // Access to thermodynamic state variables
+
+            //- Internal energy [J/kg]
+            //  Non-const access allowed for transport equations
+            virtual volScalarField& e()
+            {
+                return e_;
+            }
+
+            //- Internal energy [J/kg]
+            virtual const volScalarField& e() const
+            {
+                return e_;
+            }
+
+            //CL: drhode needed for pressure equation of the real gas solver
+            virtual const volScalarField& drhode() const
+            {
+                return drhode_;
+            }
+
+        // Fields derived from thermodynamic state variables
+
+            //- Internal energy for cell-set [J/kg]
+            virtual tmp<scalarField> e
+            (
+                const scalarField& T,
+                const labelList& cells
+            ) const;
+
+            //- Internal energy for patch [J/kg]
+            virtual tmp<scalarField> e
+            (
+                const scalarField& T,
+                const label patchi
+            ) const;
+
+            //- Density for patch [J/kg]
+            virtual tmp<scalarField> rho
+            (
+                const scalarField& T,
+                const label patchi
+            ) const;
+
+            //- Heat capacity at constant pressure for patch [J/kg/K]
+            virtual tmp<scalarField> Cp
+            (
+                const scalarField& T,
+                const label patchi
+            ) const;
+
+            //- Heat capacity at constant pressure [J/kg/K]
+            virtual tmp<volScalarField> Cp() const;
+
+            //- Heat capacity at constant volume for patch [J/kg/K]
+            virtual tmp<scalarField> Cv
+            (
+                const scalarField& T,
+                const label patchi
+            ) const;
+
+            //- Heat capacity at constant volume [J/kg/K]
+            virtual tmp<volScalarField> Cv() const;
+
+            //- Density [kg/m^3] - uses current value of pressure
+            virtual tmp<volScalarField> rho() const;
+
+
+        //- Read thermophysicalProperties dictionary
+        virtual bool read();
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
+
+#ifdef NoRepository
+#   include "realGasEThermo.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/thermophysicalModels/basic/psiThermo/realGasEThermo/realGasEThermos.C

@@ -0,0 +1,221 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           |
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+Author
+Christian Lucas
+Institut für Thermodynamik
+Technische Universität Braunschweig
+Germany
+
+\*---------------------------------------------------------------------------*/
+
+#include "makeBasicPsiThermo.H"
+
+#include "redlichKwong.H"
+#include "pengRobinson.H"
+#include "aungierRedlichKwong.H"
+#include "soaveRedlichKwong.H"
+#include "nasaHeatCapacityPolynomial.H"
+#include "realGasSpecieThermo.H"
+#include "constTransport.H"
+#include "sutherlandTransport.H"
+#include "constantHeatCapacity.H"
+#include "pureMixture.H"
+#include "realGasEThermo.H"
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/* * * * * * * * * * * * * * * private static data * * * * * * * * * * * * * */
+
+makeBasicRealGasThermo
+(
+    realGasEThermo,
+    pureMixture,
+    sutherlandTransport,
+    realGasSpecieThermo,
+    nasaHeatCapacityPolynomial,
+    pengRobinson
+);
+
+makeBasicRealGasThermo
+(
+    realGasEThermo,
+    pureMixture,
+    sutherlandTransport,
+    realGasSpecieThermo,
+    nasaHeatCapacityPolynomial,
+    aungierRedlichKwong
+);
+
+makeBasicRealGasThermo
+(
+    realGasEThermo,
+    pureMixture,
+    sutherlandTransport,
+    realGasSpecieThermo,
+    nasaHeatCapacityPolynomial,
+    redlichKwong
+);
+
+makeBasicRealGasThermo
+(
+    realGasEThermo,
+    pureMixture,
+    sutherlandTransport,
+    realGasSpecieThermo,
+    nasaHeatCapacityPolynomial,
+    soaveRedlichKwong
+);
+
+makeBasicRealGasThermo
+(
+    realGasEThermo,
+    pureMixture,
+    constTransport,
+    realGasSpecieThermo,
+    nasaHeatCapacityPolynomial,
+    pengRobinson
+);
+
+makeBasicRealGasThermo
+(
+    realGasEThermo,
+    pureMixture,
+    constTransport,
+    realGasSpecieThermo,
+    nasaHeatCapacityPolynomial,
+    aungierRedlichKwong
+);
+
+makeBasicRealGasThermo
+(
+    realGasEThermo,
+    pureMixture,
+    constTransport,
+    realGasSpecieThermo,
+    nasaHeatCapacityPolynomial,
+    redlichKwong
+);
+
+makeBasicRealGasThermo
+(
+    realGasEThermo,
+    pureMixture,
+    constTransport,
+    realGasSpecieThermo,
+    nasaHeatCapacityPolynomial,
+    soaveRedlichKwong
+);
+
+makeBasicRealGasThermo
+(
+    realGasEThermo,
+    pureMixture,
+    sutherlandTransport,
+    realGasSpecieThermo,
+    constantHeatCapacity,
+    pengRobinson
+);
+
+
+makeBasicRealGasThermo
+(
+    realGasEThermo,
+    pureMixture,
+    sutherlandTransport,
+    realGasSpecieThermo,
+    constantHeatCapacity,
+    aungierRedlichKwong
+);
+
+makeBasicRealGasThermo
+(
+    realGasEThermo,
+    pureMixture,
+    sutherlandTransport,
+    realGasSpecieThermo,
+    constantHeatCapacity,
+    redlichKwong
+);
+
+makeBasicRealGasThermo
+(
+    realGasEThermo,
+    pureMixture,
+    sutherlandTransport,
+    realGasSpecieThermo,
+    constantHeatCapacity,
+    soaveRedlichKwong
+);
+
+makeBasicRealGasThermo
+(
+    realGasEThermo,
+    pureMixture,
+    constTransport,
+    realGasSpecieThermo,
+    constantHeatCapacity,
+    pengRobinson
+);
+
+makeBasicRealGasThermo
+(
+    realGasEThermo,
+    pureMixture,
+    constTransport,
+    realGasSpecieThermo,
+    constantHeatCapacity,
+    aungierRedlichKwong
+);
+
+makeBasicRealGasThermo
+(
+    realGasEThermo,
+    pureMixture,
+    constTransport,
+    realGasSpecieThermo,
+    constantHeatCapacity,
+    redlichKwong
+);
+
+makeBasicRealGasThermo
+(
+    realGasEThermo,
+    pureMixture,
+    constTransport,
+    realGasSpecieThermo,
+    constantHeatCapacity,
+    soaveRedlichKwong
+);
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// ************************************************************************* //

src/thermophysicalModels/basic/psiThermo/realGasHThermo/realGasHThermo.C

@@ -0,0 +1,505 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           |
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+Author
+Christian Lucas
+Institut für Thermodynamik
+Technische Universität Braunschweig
+Germany
+
+\*---------------------------------------------------------------------------*/
+
+#include "realGasHThermo.H"
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+template<class MixtureType>
+void Foam::realGasHThermo<MixtureType>::calculate()
+{
+    const scalarField& hCells = h_.internalField();
+    const scalarField& pCells = this->p_.internalField();
+
+    scalarField& TCells = this->T_.internalField();
+    scalarField& rhoCells= this->rho_.internalField();
+    scalarField& psiCells = this->psi_.internalField();
+    scalarField& drhodhCells = this->drhodh_.internalField();
+    scalarField& muCells = this->mu_.internalField();
+    scalarField& alphaCells = this->alpha_.internalField();
+
+    forAll(TCells, celli)
+    {
+        const typename MixtureType::thermoType& mixture_ =
+            this->cellMixture(celli);
+
+        mixture_.TH(hCells[celli], TCells[celli], pCells[celli], rhoCells[celli]);
+        psiCells[celli] = mixture_.psi(rhoCells[celli], TCells[celli]);
+        drhodhCells[celli] = mixture_.drhodH(rhoCells[celli], TCells[celli]);
+        muCells[celli] = mixture_.mu(TCells[celli]);
+        alphaCells[celli] = mixture_.alpha(rhoCells[celli], TCells[celli]);
+    }
+
+    forAll(T_.boundaryField(), patchi)
+    {
+        fvPatchScalarField& pp = this->p_.boundaryField()[patchi];
+        fvPatchScalarField& pT = this->T_.boundaryField()[patchi];
+        fvPatchScalarField& ppsi = this->psi_.boundaryField()[patchi];
+        fvPatchScalarField& pdrhodh = this->drhodh_.boundaryField()[patchi];
+        fvPatchScalarField& prho = this->rho_.boundaryField()[patchi];
+        fvPatchScalarField& ph = h_.boundaryField()[patchi];
+        fvPatchScalarField& pmu = this->mu_.boundaryField()[patchi];
+        fvPatchScalarField& palpha = this->alpha_.boundaryField()[patchi];
+
+        if (pT.fixesValue())
+        {
+            forAll(pT, facei)
+            {
+                const typename MixtureType::thermoType& mixture_ =
+                    this->patchFaceMixture(patchi, facei);
+
+                prho[facei] = mixture_.rho(pp[facei], pT[facei],prho[facei]);
+                ppsi[facei]=mixture_.psi(prho[facei],pT[facei]);
+                pdrhodh[facei]=mixture_.drhodH(prho[facei],pT[facei]);
+                ph[facei] = mixture_.H(prho[facei], pT[facei]);
+                pmu[facei] = mixture_.mu(pT[facei]);
+                palpha[facei] = mixture_.alpha(prho[facei],pT[facei]);
+            }
+        }
+        else
+        {
+            forAll(pT, facei)
+            {
+                const typename MixtureType::thermoType& mixture_ =
+                    this->patchFaceMixture(patchi, facei);
+
+                mixture_.TH(ph[facei], pT[facei],pp[facei],prho[facei]);
+                pmu[facei] = mixture_.mu(pT[facei]);
+                ppsi[facei]=mixture_.psi(prho[facei],pT[facei]);
+                pdrhodh[facei]=mixture_.drhodH(prho[facei],pT[facei]);
+                palpha[facei] = mixture_.alpha(prho[facei],pT[facei]);
+            }
+        }
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+template<class MixtureType>
+Foam::realGasHThermo<MixtureType>::realGasHThermo
+(
+    const fvMesh& mesh,
+    const objectRegistry& obj
+)
+:
+    basicPsiThermo(mesh, obj),
+    MixtureType(*this, mesh, obj),
+
+    h_
+    (
+        IOobject
+        (
+            "h",
+            mesh.time().timeName(),
+            obj,
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        ),
+        mesh,
+        dimensionSet(0, 2, -2, 0, 0),
+        this->hBoundaryTypes()
+    ),
+    rho_
+    (
+        IOobject
+        (
+            "rhoThermo",
+            mesh.time().timeName(),
+            obj,
+            IOobject::READ_IF_PRESENT,
+            IOobject::NO_WRITE
+        ),
+        mesh,
+        dimDensity
+    ),
+    drhodh_
+    (
+        IOobject
+        (
+            "drhodh",
+            mesh.time().timeName(),
+            obj,
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        ),
+        mesh,
+        dimensionSet(1, -5, 2, 0, 0)
+    )
+{
+    scalarField& hCells = h_.internalField();
+    const scalarField& TCells = this->T_.internalField();
+    const scalarField& pCells =this->p_.internalField();
+    scalarField& rhoCells =this->rho_.internalField();
+
+    forAll(rhoCells, celli)
+    {
+        rhoCells[celli]=this->cellMixture(celli).rho(pCells[celli],TCells[celli]);
+    }
+
+    forAll(rho_.boundaryField(), patchi)
+    {
+        rho_.boundaryField()[patchi] ==
+            rho(this->T_.boundaryField()[patchi], patchi);
+    }
+
+    forAll(hCells, celli)
+    {
+        hCells[celli] = this->cellMixture(celli).H(rhoCells[celli],TCells[celli]);
+    }
+
+    forAll(h_.boundaryField(), patchi)
+    {
+        h_.boundaryField()[patchi] ==
+            h(this->T_.boundaryField()[patchi], patchi);
+    }
+
+    hBoundaryCorrection(h_);
+
+    calculate();
+
+    // Switch on saving old time
+    this->psi_.oldTime();
+}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+template<class MixtureType>
+Foam::realGasHThermo<MixtureType>::~realGasHThermo()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+template<class MixtureType>
+void Foam::realGasHThermo<MixtureType>::correct()
+{
+    if (debug)
+    {
+        Info<< "entering realGasHThermo<MixtureType>::correct()" << endl;
+    }
+
+    // force the saving of the old-time values
+    this->psi_.oldTime();
+
+    calculate();
+
+    if (debug)
+    {
+        Info<< "exiting realGasHThermo<MixtureType>::correct()" << endl;
+    }
+}
+
+
+template<class MixtureType>
+Foam::tmp<Foam::scalarField> Foam::realGasHThermo<MixtureType>::h
+(
+    const scalarField& T,
+    const labelList& cells
+) const
+{
+    //CL: need the pressure of the internal field to calculate the realGas enthalpy
+    //CL: this is done this way to assure compatibility to old OF Thermo-Versions
+    const scalarField& pCells = this->p_.internalField();
+
+    tmp<scalarField> th(new scalarField(T.size()));
+    scalarField& h = th();
+
+    forAll(T, celli)
+    {
+        h[celli] = this->cellMixture(cells[celli]).H(this->cellMixture(cells[celli]).rho(pCells[cells[celli]],T[celli]),T[celli]);
+    }
+
+    return th;
+}
+
+
+template<class MixtureType>
+Foam::tmp<Foam::scalarField> Foam::realGasHThermo<MixtureType>::h
+(
+    const scalarField& T,
+    const label patchi
+) const
+{
+    //CL: need the pressure at the patch to calculate the realGas enthalpy
+    //CL: this is done this way to assure compatibility to old OF Thermo-Versions
+    const fvPatchScalarField& pp = this->p_.boundaryField()[patchi];
+
+    tmp<scalarField> th(new scalarField(T.size()));
+    scalarField& h = th();
+
+    forAll(T, facei)
+    {
+        h[facei] = this->patchFaceMixture(patchi, facei).H(this->patchFaceMixture(patchi, facei).rho(pp[facei], T[facei]),T[facei]);
+    }
+
+    return th;
+}
+
+
+template<class MixtureType>
+Foam::tmp<Foam::scalarField> Foam::realGasHThermo<MixtureType>::rho
+(
+    const scalarField& T,
+    const label patchi
+) const
+{
+    //CL: need the pressure at the patch to calculate the realGas enthalpy
+    //CL: this is done this way to assure compatibility to old OF Thermo-Versions
+    const fvPatchScalarField& pp = this->p_.boundaryField()[patchi];
+
+    tmp<scalarField> trho(new scalarField(T.size()));
+    scalarField& rho = trho();
+
+    forAll(T, facei)
+    {
+        rho[facei] = this->patchFaceMixture(patchi, facei).rho(pp[facei], T[facei]);
+    }
+
+    return trho;
+}
+
+
+template<class MixtureType>
+Foam::tmp<Foam::scalarField> Foam::realGasHThermo<MixtureType>::Cp
+(
+    const scalarField& T,
+    const label patchi
+) const
+{
+    //CL: need the pressure at the patch to calculate the realGas enthalpy
+    //CL: this is done this way to assure compatibility to old OF Thermo-Versions
+    const fvPatchScalarField& pp = this->p_.boundaryField()[patchi];
+
+    tmp<scalarField> tCp(new scalarField(T.size()));
+    scalarField& cp = tCp();
+
+    forAll(T, facei)
+    {
+        cp[facei] = this->patchFaceMixture(patchi, facei).Cp(this->patchFaceMixture(patchi, facei).rho(pp[facei], T[facei]),T[facei]);
+    }
+
+    return tCp;
+}
+
+
+template<class MixtureType>
+Foam::tmp<Foam::volScalarField> Foam::realGasHThermo<MixtureType>::Cp() const
+{
+    const fvMesh& mesh = this->T_.mesh();
+
+    tmp<volScalarField> tCp
+    (
+        new volScalarField
+        (
+            IOobject
+            (
+                "Cp",
+                mesh.time().timeName(),
+                this->T_.db(),
+                IOobject::NO_READ,
+                IOobject::NO_WRITE
+            ),
+            mesh,
+            dimensionSet(0, 2, -2, -1, 0)
+        )
+    );
+
+    volScalarField& cp = tCp();
+
+    forAll(this->T_, celli)
+    {
+        cp[celli] = this->cellMixture(celli).Cp(this->rho_[celli], this->T_[celli]);
+    }
+
+    forAll(this->T_.boundaryField(), patchi)
+    {
+        const fvPatchScalarField& pT = this->T_.boundaryField()[patchi];
+        const fvPatchScalarField& prho = this->rho_.boundaryField()[patchi];
+        fvPatchScalarField& pCp = cp.boundaryField()[patchi];
+
+        forAll(pT, facei)
+        {
+            pCp[facei] = this->patchFaceMixture(patchi, facei).Cp(prho[facei], pT[facei]);
+        }
+    }
+
+    return tCp;
+}
+
+
+template<class MixtureType>
+Foam::tmp<Foam::scalarField> Foam::realGasHThermo<MixtureType>::Cv
+(
+    const scalarField& T,
+    const label patchi
+) const
+{
+    //CL: need the pressure at the patch to calculate the realGas enthalpy
+    //CL: this is done this way to assure compatibility to old OF Thermo-Versions
+    const fvPatchScalarField& pp = this->p_.boundaryField()[patchi];
+
+    tmp<scalarField> tCv(new scalarField(T.size()));
+    scalarField& cv = tCv();
+
+    forAll(T, facei)
+    {
+        cv[facei] = this->patchFaceMixture(patchi, facei).Cv(this->patchFaceMixture(patchi, facei).rho(pp[facei], T[facei]), T[facei]);
+    }
+
+    return tCv;
+}
+
+
+// CL: Maybe this function should be changed so that it is not "const" fucntion anymore
+template<class MixtureType>
+Foam::tmp<Foam::volScalarField> Foam::realGasHThermo<MixtureType>::rho()  const
+{
+
+    const fvMesh& mesh = this->T_.mesh();
+
+    //CL: create a rho Field, which will be return
+    //CL: the problem is that this function is "const",
+    //CL: so a new variabel is needed
+      tmp<volScalarField> trho
+      (
+        new volScalarField
+        (
+            IOobject
+            (
+                "rhoFunctionThermo",
+                mesh.time().timeName(),
+                this->T_.db(),
+                IOobject::NO_READ,
+                IOobject::NO_WRITE
+            ),
+            mesh,
+            dimDensity
+        )
+    );
+
+   //CL: copy "old" rho value onto the new rho field as start point
+   //CL: for the newton solver used in this->TH( ... )
+   trho()=rho_;
+
+   volScalarField& rho = trho();
+
+    const scalarField& hCells = h_.internalField();
+    const scalarField& pCells = this->p_.internalField();
+    scalarField TCells = this->T_.internalField();
+
+    forAll(pCells, celli)
+    {
+        const typename MixtureType::thermoType& mixture_ =
+            this->cellMixture(celli);
+
+        // getting the new rho Field
+        mixture_.TH(hCells[celli], TCells[celli], pCells[celli], rho[celli]);
+    }
+
+    forAll(p_.boundaryField(), patchi)
+    {
+        fvPatchScalarField pp = this->p_.boundaryField()[patchi];
+        fvPatchScalarField ph = h_.boundaryField()[patchi];
+        fvPatchScalarField pT = this->T_.boundaryField()[patchi];
+
+        fvPatchScalarField& prho_ = rho.boundaryField()[patchi];
+
+        forAll(pp, facei)
+        {
+            const typename MixtureType::thermoType& mixture_ =
+                this->patchFaceMixture(patchi, facei);
+
+            // getting the new rho patch Field
+            mixture_.TH(ph[facei], pT[facei],pp[facei],prho_[facei]);
+        }
+    }
+    return trho;
+}
+
+
+template<class MixtureType>
+Foam::tmp<Foam::volScalarField> Foam::realGasHThermo<MixtureType>::Cv() const
+{
+    const fvMesh& mesh = this->T_.mesh();
+
+    tmp<volScalarField> tCv
+    (
+        new volScalarField
+        (
+            IOobject
+            (
+                "Cv",
+                mesh.time().timeName(),
+                this->T_.db(),
+                IOobject::NO_READ,
+                IOobject::NO_WRITE
+            ),
+            mesh,
+            dimensionSet(0, 2, -2, -1, 0)
+        )
+    );
+
+    volScalarField& cv = tCv();
+
+    forAll(this->T_, celli)
+    {
+        cv[celli] = this->cellMixture(celli).Cv(this->rho_[celli], this->T_[celli]);
+    }
+
+    forAll(this->T_.boundaryField(), patchi)
+    {
+        cv.boundaryField()[patchi] =
+            Cv(this->T_.boundaryField()[patchi], patchi);
+    }
+
+    return tCv;
+}
+
+
+template<class MixtureType>
+bool Foam::realGasHThermo<MixtureType>::read()
+{
+    if (basicPsiThermo::read())
+    {
+        MixtureType::read(*this);
+        return true;
+    }
+    else
+    {
+        return false;
+    }
+}
+
+
+// ************************************************************************* //

src/thermophysicalModels/basic/psiThermo/realGasHThermo/realGasHThermo.H

@@ -0,0 +1,205 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | .
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+Class
+    Foam::realGasHThermo
+
+Description
+    Enthalpy for a real gas fluid libary
+
+
+SourceFiles
+    realGasHThermo.C
+
+Author
+Christian Lucas
+Institut für Thermodynamik
+Technische Universität Braunschweig
+Germany
+
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef realGasHThermo_H
+#define realGasHThermo_H
+
+#include "basicPsiThermo.H"
+#include "basicMixture.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                         Class realGasHThermo Declaration
+\*---------------------------------------------------------------------------*/
+
+template<class MixtureType>
+class realGasHThermo
+:
+    public basicPsiThermo,
+    public MixtureType
+{
+    // Private data
+
+        //- Enthalpy field
+        volScalarField h_;
+
+        //- DensityField
+        volScalarField rho_;
+
+        //- drhodh_Field
+        volScalarField drhodh_;
+
+    // Private member functions
+
+        //- Calculate the thermo variables
+        void calculate();
+
+        //- Construct as copy (not implemented)
+        realGasHThermo(const realGasHThermo<MixtureType>&);
+
+
+public:
+
+    //- Runtime type information
+    TypeName("realGasHThermo");
+
+
+    // Constructors
+
+        //- Construct from mesh
+        realGasHThermo(const fvMesh&, const objectRegistry& obj);
+
+    //- Destructor
+    virtual ~realGasHThermo();
+
+
+    // Member functions
+
+        //- Return the compostion of the mixture
+        virtual basicMixture& composition()
+        {
+            return *this;
+        }
+
+        //- Return the compostion of the mixture
+        virtual const basicMixture& composition() const
+        {
+            return *this;
+        }
+
+        //- Update properties
+        virtual void correct();
+
+
+        // Access to thermodynamic state variables
+
+            //- Enthalpy [J/kg]
+            //  Non-const access allowed for transport equations
+            virtual volScalarField& h()
+            {
+                return h_;
+            }
+
+            //- Enthalpy [J/kg]
+            virtual const volScalarField& h() const
+            {
+                return h_;
+            }
+
+            //CL: drhodh needed for pressure equation of the real gas solver
+            virtual const volScalarField& drhodh() const
+            {
+                return drhodh_;
+            }
+
+        // Fields derived from thermodynamic state variables
+
+            //- Enthalpy for cell-set [J/kg]
+            virtual tmp<scalarField> h
+            (
+                const scalarField& T,
+                const labelList& cells
+            ) const;
+
+            //- Enthalpy for patch [J/kg]
+            virtual tmp<scalarField> h
+            (
+                const scalarField& T,
+                const label patchi
+            ) const;
+
+            //- Density for patch [J/kg]
+            virtual tmp<scalarField> rho
+            (
+                const scalarField& T,
+                const label patchi
+            ) const;
+
+            //- Heat capacity at constant pressure for patch [J/kg/K]
+            virtual tmp<scalarField> Cp
+            (
+                const scalarField& T,
+                const label patchi
+            ) const;
+
+            //- Heat capacity at constant pressure [J/kg/K]
+            virtual tmp<volScalarField> Cp() const;
+
+            //- Heat capacity at constant volume for patch [J/kg/K]
+            virtual tmp<scalarField> Cv
+            (
+                const scalarField& T,
+                const label patchi
+            ) const;
+
+            //- Heat capacity at constant volume [J/kg/K]
+            virtual tmp<volScalarField> Cv() const;
+
+            //- Density [kg/m^3] - uses current value of pressure
+            virtual tmp<volScalarField> rho() const;
+
+
+        //- Read thermophysicalProperties dictionary
+        virtual bool read();
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
+
+#ifdef NoRepository
+#   include "realGasHThermo.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/thermophysicalModels/basic/psiThermo/realGasHThermo/realGasHThermos.C

@@ -0,0 +1,221 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           |
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+Author
+Christian Lucas
+Institut für Thermodynamik
+Technische Universität Braunschweig
+Germany
+
+\*---------------------------------------------------------------------------*/
+
+#include "makeBasicPsiThermo.H"
+
+#include "redlichKwong.H"
+#include "pengRobinson.H"
+#include "aungierRedlichKwong.H"
+#include "soaveRedlichKwong.H"
+#include "nasaHeatCapacityPolynomial.H"
+#include "realGasSpecieThermo.H"
+#include "constTransport.H"
+#include "sutherlandTransport.H"
+#include "constantHeatCapacity.H"
+#include "pureMixture.H"
+#include "realGasHThermo.H"
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/* * * * * * * * * * * * * * * private static data * * * * * * * * * * * * * */
+
+makeBasicRealGasThermo
+(
+    realGasHThermo,
+    pureMixture,
+    sutherlandTransport,
+    realGasSpecieThermo,
+    nasaHeatCapacityPolynomial,
+    pengRobinson
+);
+
+makeBasicRealGasThermo
+(
+    realGasHThermo,
+    pureMixture,
+    sutherlandTransport,
+    realGasSpecieThermo,
+    nasaHeatCapacityPolynomial,
+    aungierRedlichKwong
+);
+
+makeBasicRealGasThermo
+(
+    realGasHThermo,
+    pureMixture,
+    sutherlandTransport,
+    realGasSpecieThermo,
+    nasaHeatCapacityPolynomial,
+    redlichKwong
+);
+
+makeBasicRealGasThermo
+(
+    realGasHThermo,
+    pureMixture,
+    sutherlandTransport,
+    realGasSpecieThermo,
+    nasaHeatCapacityPolynomial,
+    soaveRedlichKwong
+);
+
+makeBasicRealGasThermo
+(
+    realGasHThermo,
+    pureMixture,
+    constTransport,
+    realGasSpecieThermo,
+    nasaHeatCapacityPolynomial,
+    pengRobinson
+);
+
+makeBasicRealGasThermo
+(
+    realGasHThermo,
+    pureMixture,
+    constTransport,
+    realGasSpecieThermo,
+    nasaHeatCapacityPolynomial,
+    aungierRedlichKwong
+);
+
+makeBasicRealGasThermo
+(
+    realGasHThermo,
+    pureMixture,
+    constTransport,
+    realGasSpecieThermo,
+    nasaHeatCapacityPolynomial,
+    redlichKwong
+);
+
+makeBasicRealGasThermo
+(
+    realGasHThermo,
+    pureMixture,
+    constTransport,
+    realGasSpecieThermo,
+    nasaHeatCapacityPolynomial,
+    soaveRedlichKwong
+);
+
+makeBasicRealGasThermo
+(
+    realGasHThermo,
+    pureMixture,
+    sutherlandTransport,
+    realGasSpecieThermo,
+    constantHeatCapacity,
+    pengRobinson
+);
+
+
+makeBasicRealGasThermo
+(
+    realGasHThermo,
+    pureMixture,
+    sutherlandTransport,
+    realGasSpecieThermo,
+    constantHeatCapacity,
+    aungierRedlichKwong
+);
+
+makeBasicRealGasThermo
+(
+    realGasHThermo,
+    pureMixture,
+    sutherlandTransport,
+    realGasSpecieThermo,
+    constantHeatCapacity,
+    redlichKwong
+);
+
+makeBasicRealGasThermo
+(
+    realGasHThermo,
+    pureMixture,
+    sutherlandTransport,
+    realGasSpecieThermo,
+    constantHeatCapacity,
+    soaveRedlichKwong
+);
+
+makeBasicRealGasThermo
+(
+    realGasHThermo,
+    pureMixture,
+    constTransport,
+    realGasSpecieThermo,
+    constantHeatCapacity,
+    pengRobinson
+);
+
+makeBasicRealGasThermo
+(
+    realGasHThermo,
+    pureMixture,
+    constTransport,
+    realGasSpecieThermo,
+    constantHeatCapacity,
+    aungierRedlichKwong
+);
+
+makeBasicRealGasThermo
+(
+    realGasHThermo,
+    pureMixture,
+    constTransport,
+    realGasSpecieThermo,
+    constantHeatCapacity,
+    redlichKwong
+);
+
+makeBasicRealGasThermo
+(
+    realGasHThermo,
+    pureMixture,
+    constTransport,
+    realGasSpecieThermo,
+    constantHeatCapacity,
+    soaveRedlichKwong
+);
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// ************************************************************************* //

src/thermophysicalModels/specie/Make/files

@@ -8,6 +8,10 @@ $(atomicWeights)/atomicWeights.C
 $(specie)/specie.C
 $(speciesTable)/speciesTable.C
 $(equationOfState)/perfectGas/perfectGas.C
+$(equationOfState)/cubicEquationOfState/redlichKwong/redlichKwong.C
+$(equationOfState)/cubicEquationOfState/aungierRedlichKwong/aungierRedlichKwong.C
+$(equationOfState)/cubicEquationOfState/pengRobinson/pengRobinson.C
+$(equationOfState)/cubicEquationOfState/soaveRedlichKwong/soaveRedlichKwong.C
 $(reactions)/makeChemkinReactions.C
 $(reactions)/makeReactionThermoReactions.C
 $(reactions)/makeLangmuirHinshelwoodReactions.C

src/thermophysicalModels/specie/equationOfState/cubicEquationOfState/aungierRedlichKwong/aungierRedlichKwong.C

@@ -0,0 +1,137 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright held by original author
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+Description
+    Aungier Redlich Kwong equation of state.
+
+Author
+Christian Lucas
+Institut für Thermodynamik
+Technische Universität Braunschweig
+Germany
+
+\*---------------------------------------------------------------------------*/
+
+#include "aungierRedlichKwong.H"
+#include "IOstreams.H"
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::aungierRedlichKwong::aungierRedlichKwong(Istream& is)
+:
+    specie(is),
+    pcrit_(readScalar(is)),
+    Tcrit_(readScalar(is)),
+    rhocrit_(readScalar(is)),
+    azentricFactor_(readScalar(is)),
+    a0_(0.42747*pow(this->RR(), 2)*pow(Tcrit_, 2)/pcrit_),
+    b_(0.08664*this->RR()*Tcrit_/pcrit_),
+    c_(this->RR()*Tcrit_/(pcrit_+(a0_/(this->W()/rhocrit_*(this->W()/rhocrit_ + b_)))) + b_ - this->W()/rhocrit_),
+    n_(0.4986 + 1.2735*azentricFactor_ + 0.4754*pow(azentricFactor_, 2)),
+    b2_(b_*b_),
+    b3_(b2_*b_),
+    b4_(b3_*b_),
+    b5_(b4_*b_),
+    c2_(c_*c_),
+    //CL: Only uses the default values
+    rhoMin_(1e-3),
+    rhoMax_(1500),
+    // Starting GUESS for the density by ideal gas law
+    rhostd_(this->rho(this->Pstd(), this->Tstd(), this->Pstd()/(this->Tstd()*this->R()))),
+    aSave(0.0),
+    daSave(0.0),
+    d2aSave(0.0),
+    TSave(0.0)
+{
+    is.check("aungierRedlichKwong::aungierRedlichKwong(Istream& is)");
+}
+
+//CL: Constructed needed in OpenFOAM 2.x.x
+//CL: Code works fine, but compiling problem in OpenFOAM 1.6.ext
+//CL:  because specie has no constructor using dict
+/*
+Foam::aungierRedlichKwong::aungierRedlichKwong(const dictionary& dict)
+:
+    specie(dict),
+    pcrit_(readScalar(dict.subDict("equationOfState").lookup("pCritical"))),
+    Tcrit_(readScalar(dict.subDict("equationOfState").lookup("TCritical"))),
+    rhocrit_(readScalar(dict.subDict("equationOfState").lookup("rhoCritical"))),
+    azentricFactor_(readScalar(dict.subDict("equationOfState").lookup("azentricFactor"))),
+    a0_(0.42747*pow(this->RR(), 2)*pow(Tcrit_, 2)/pcrit_),
+    b_(0.08664*this->RR()*Tcrit_/pcrit_),
+    c_(this->RR()*Tcrit_/(pcrit_+(a0_/(this->W()/rhocrit_*(this->W()/rhocrit_+b_))))+b_-this->W()/rhocrit_),
+    n_(0.4986+1.2735*azentricFactor_+0.4754*pow(azentricFactor_,2)),
+    b2_(b_*b_),
+    b3_(b2_*b_),
+    b4_(b3_*b_),
+    b5_(b4_*b_),
+    c2_(pow(c_,2)),
+    //CL: rhoMin and rhoMax are only used as boundaries for the bisection method (see rho function)
+    //CL: important: rhoMin and rhoMax are not used as boundary for the newton solver
+    //CL: therefore, rho can be larger than rhoMax and smaller than rhoMin
+    rhoMin_(dict.subDict("equationOfState").lookupOrDefault("rhoMin",1e-3)),
+    rhoMax_(dict.subDict("equationOfState").lookupOrDefault("rhoMax",1500)),
+    // Starting GUESS for the density by ideal gas law
+    rhostd_(this->rho(this->Pstd(), this->Tstd(), this->Pstd()/(this->Tstd()*this->R()))),
+    aSave(0.0),
+    daSave(0.0),
+    d2aSave(0.0),
+    TSave(0.0)
+{ 
+    is.check("aungierRedlichKwong::aungierRedlichKwong(Istream& is)");
+}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+void Foam::aungierRedlichKwong::write(Ostream& os) const
+{
+    specie::write(os);
+
+    dictionary dict("equationOfState");
+    dict.add("pCritical", pcrit_);
+    dict.add("TCritical", Tcrit_);
+    dict.add("azentricFactor", azentricFactor_);
+    dict.add("rhoCritical", rhocrit_);
+    dict.add("rhoMin", rhoMin_);
+    dict.add("rhoMax", rhoMax_);
+
+    os  << indent << dict.dictName() << dict;
+}
+*/
+
+
+// * * * * * * * * * * * * * * * Ostream Operator  * * * * * * * * * * * * * //
+
+Foam::Ostream& Foam::operator<<(Ostream& os, const aungierRedlichKwong& ark)
+{
+    os  << static_cast<const specie&>(ark)<< token::SPACE
+        << ark.pcrit_ << tab<< ark.Tcrit_<< tab<<ark.azentricFactor_<< tab<<ark.rhocrit_;
+
+    os.check("Ostream& operator<<(Ostream& os, const aungierRedlichKwong& st)");
+    return os;
+}
+
+
+// ************************************************************************* //