GIT: Resolved merge conflict when merging develop branch

.gitignore

@@ -62,6 +62,7 @@ doc/Doxygen/DTAGS
 /etc/prefs.sh
 /etc/config.csh/prefs.csh
 /etc/config.sh/prefs.sh
+/wmake/rules/General/mplibUSER*
 
 # Source packages - anywhere
 *.tar.bz2

Allwmake

@@ -5,8 +5,8 @@ cd ${0%/*} && wmakeCheckPwd "$WM_PROJECT_DIR" 2>/dev/null || {
     echo "    Check your OpenFOAM environment and installation"
     exit 1
 }
-[ -n "$FOAM_EXT_LIBBIN" ] || {
+[ -d "$WM_PROJECT_DIR" -a -f "$WM_PROJECT_DIR/etc/bashrc" ] || {
-    echo "Error (${0##*/}) : FOAM_EXT_LIBBIN not set"
+    echo "Error (${0##*/}) : WM_PROJECT_DIR appears to be incorrect"
     echo "    Check your OpenFOAM environment and installation"
     exit 1
 }
@@ -22,8 +22,8 @@ echo "  $WM_COMPILER $WM_COMPILER_TYPE compiler"
 echo "  ${WM_OPTIONS}, with ${WM_MPLIB} ${FOAM_MPI}"
 echo
 
-# Compile wmake support applications
+# Compile wmake tools
-(cd wmake/src && make)
+(cd "${WM_DIR:-wmake}/src" && make)
 
 # Compile ThirdParty libraries and applications
 if [ -d "$WM_THIRD_PARTY_DIR" ]

README.md

@@ -1,5 +1,5 @@
 # About OpenFOAM
-OpenFOAM is a free, open source CFD software [released and developed primarily by OpenCFD Ltd](http://www.openfoam.com) since 2004released and developed primarily by. It has a large user base across most areas of engineering and science, from both commercial and academic organisations. OpenFOAM has an extensive range of features to solve anything from complex fluid flows involving chemical reactions, turbulence and heat transfer, to acoustics, solid mechanics and electromagnetics.  [More...](http://www.openfoam.com/documentation)
+OpenFOAM is a free, open source CFD software [released and developed primarily by OpenCFD Ltd](http://www.openfoam.com) since 2004. It has a large user base across most areas of engineering and science, from both commercial and academic organisations. OpenFOAM has an extensive range of features to solve anything from complex fluid flows involving chemical reactions, turbulence and heat transfer, to acoustics, solid mechanics and electromagnetics.  [More...](http://www.openfoam.com/documentation)
 
 OpenFOAM+ is professionally released every six months to include customer sponsored developments and contributions from the community, including the OpenFOAM Foundation. Releases designated OpenFOAM+ contain several man years of client-sponsored developments of which much has been transferred to, but not released in the OpenFOAM Foundation branch.
 
@@ -23,4 +23,4 @@ Violations of the Trademark are continuously monitored, and will be duly prosecu
 - [OpenFOAM Community](http://www.openfoam.com/community/)
 - [Contacting OpenCFD](http://www.openfoam.com/contact/)
 
-Copyright 2016 OpenCFD Ltd
+Copyright 2016-2017 OpenCFD Ltd

applications/Allwmake

@@ -5,8 +5,8 @@ cd ${0%/*} && wmakeCheckPwd "$WM_PROJECT_DIR/applications" 2>/dev/null || {
     echo "    Check your OpenFOAM environment and installation"
     exit 1
 }
-[ -n "$FOAM_EXT_LIBBIN" ] || {
+[ -d "$WM_PROJECT_DIR" -a -f "$WM_PROJECT_DIR/etc/bashrc" ] || {
-    echo "Error (${0##*/}) : FOAM_EXT_LIBBIN not set"
+    echo "Error (${0##*/}) : WM_PROJECT_DIR appears to be incorrect"
     echo "    Check your OpenFOAM environment and installation"
     exit 1
 }

applications/solvers/basic/laplacianFoam/Make/options

@@ -4,4 +4,5 @@ EXE_INC = \
 
 EXE_LIBS = \
     -lfiniteVolume \
+    -lfvOptions \
     -lmeshTools

applications/solvers/basic/laplacianFoam/createFields.H

@@ -1,37 +1,39 @@
-    Info<< "Reading field T\n" << endl;
+Info<< "Reading field T\n" << endl;
 
-    volScalarField T
+volScalarField T
+(
+    IOobject
     (
-        IOobject
+        "T",
-        (
+        runTime.timeName(),
-            "T",
+        mesh,
-            runTime.timeName(),
+        IOobject::MUST_READ,
-            mesh,
+        IOobject::AUTO_WRITE
-            IOobject::MUST_READ,
+    ),
-            IOobject::AUTO_WRITE
+    mesh
-        ),
+);
-        mesh
-    );
 
 
-    Info<< "Reading transportProperties\n" << endl;
+Info<< "Reading transportProperties\n" << endl;
 
-    IOdictionary transportProperties
+IOdictionary transportProperties
+(
+    IOobject
     (
-        IOobject
+        "transportProperties",
-        (
+        runTime.constant(),
-            "transportProperties",
+        mesh,
-            runTime.constant(),
+        IOobject::MUST_READ_IF_MODIFIED,
-            mesh,
+        IOobject::NO_WRITE
-            IOobject::MUST_READ_IF_MODIFIED,
+    )
-            IOobject::NO_WRITE
+);
-        )
-    );
 
 
-    Info<< "Reading diffusivity DT\n" << endl;
+Info<< "Reading diffusivity DT\n" << endl;
 
-    dimensionedScalar DT
+dimensionedScalar DT
-    (
+(
-        transportProperties.lookup("DT")
+    "DT",
-    );
+    dimArea/dimTime,
+    transportProperties
+);

applications/solvers/basic/scalarTransportFoam/createFields.H

@@ -1,55 +1,57 @@
-    Info<< "Reading field T\n" << endl;
+Info<< "Reading field T\n" << endl;
 
-    volScalarField T
+volScalarField T
+(
+    IOobject
     (
-        IOobject
+        "T",
-        (
+        runTime.timeName(),
-            "T",
+        mesh,
-            runTime.timeName(),
+        IOobject::MUST_READ,
-            mesh,
+        IOobject::AUTO_WRITE
-            IOobject::MUST_READ,
+    ),
-            IOobject::AUTO_WRITE
+    mesh
-        ),
+);
-        mesh
-    );
 
 
-    Info<< "Reading field U\n" << endl;
+Info<< "Reading field U\n" << endl;
 
-    volVectorField U
+volVectorField U
+(
+    IOobject
     (
-        IOobject
+        "U",
-        (
+        runTime.timeName(),
-            "U",
+        mesh,
-            runTime.timeName(),
+        IOobject::MUST_READ,
-            mesh,
+        IOobject::AUTO_WRITE
-            IOobject::MUST_READ,
+    ),
-            IOobject::AUTO_WRITE
+    mesh
-        ),
+);
-        mesh
-    );
 
 
-    Info<< "Reading transportProperties\n" << endl;
+Info<< "Reading transportProperties\n" << endl;
 
-    IOdictionary transportProperties
+IOdictionary transportProperties
+(
+    IOobject
     (
-        IOobject
+        "transportProperties",
-        (
+        runTime.constant(),
-            "transportProperties",
+        mesh,
-            runTime.constant(),
+        IOobject::MUST_READ_IF_MODIFIED,
-            mesh,
+        IOobject::NO_WRITE
-            IOobject::MUST_READ_IF_MODIFIED,
+    )
-            IOobject::NO_WRITE
+);
-        )
-    );
 
 
-    Info<< "Reading diffusivity DT\n" << endl;
+Info<< "Reading diffusivity DT\n" << endl;
 
-    dimensionedScalar DT
+dimensionedScalar DT
-    (
+(
-        transportProperties.lookup("DT")
+    "DT",
-    );
+    dimArea/dimTime,
+    transportProperties
+);
 
-    #include "createPhi.H"
+#include "createPhi.H"

applications/solvers/combustion/PDRFoam/Make/options

@@ -15,8 +15,8 @@ EXE_INC = \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
     -I$(LIB_SRC)/dynamicFvMesh/lnInclude \
     -I$(LIB_SRC)/dynamicMesh/lnInclude \
-    -I$(LIB_SRC)/meshTools/lnInclude \
+    -I$(LIB_SRC)/surfMesh/lnInclude \
-    -I$(LIB_SRC)/triSurface/lnInclude \
+    -I$(LIB_SRC)/meshTools/lnInclude
 
 EXE_LIBS = \
     -lengine \

applications/solvers/combustion/PDRFoam/PDRModels/dragModels/PDRDragModel/PDRDragModel.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2015 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -79,7 +79,7 @@ Foam::PDRDragModel::~PDRDragModel()
 
 bool Foam::PDRDragModel::read(const dictionary& PDRProperties)
 {
-    PDRDragModelCoeffs_ = PDRProperties.subDict(type() + "Coeffs");
+    PDRDragModelCoeffs_ = PDRProperties.optionalSubDict(type() + "Coeffs");
 
     PDRDragModelCoeffs_.lookup("drag") >> on_;
 

applications/solvers/combustion/PDRFoam/PDRModels/dragModels/PDRDragModel/PDRDragModelNew.C

@@ -43,7 +43,7 @@ Foam::autoPtr<Foam::PDRDragModel> Foam::PDRDragModel::New
     dictionaryConstructorTable::iterator cstrIter =
         dictionaryConstructorTablePtr_->find(modelType);
 
-    if (cstrIter == dictionaryConstructorTablePtr_->end())
+    if (!cstrIter.found())
     {
         FatalErrorInFunction
             << "Unknown PDRDragModel type "

applications/solvers/combustion/PDRFoam/XiModels/XiEqModels/XiEqModel/XiEqModel.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -67,7 +67,7 @@ Foam::XiEqModel::~XiEqModel()
 
 bool Foam::XiEqModel::read(const dictionary& XiEqProperties)
 {
-    XiEqModelCoeffs_ = XiEqProperties.subDict(type() + "Coeffs");
+    XiEqModelCoeffs_ = XiEqProperties.optionalSubDict(type() + "Coeffs");
 
     return true;
 }

applications/solvers/combustion/PDRFoam/XiModels/XiEqModels/XiEqModel/XiEqModelNew.C

@@ -42,7 +42,7 @@ Foam::autoPtr<Foam::XiEqModel> Foam::XiEqModel::New
     dictionaryConstructorTable::iterator cstrIter =
         dictionaryConstructorTablePtr_->find(modelType);
 
-    if (cstrIter == dictionaryConstructorTablePtr_->end())
+    if (!cstrIter.found())
     {
         FatalErrorInFunction
             << "Unknown XiEqModel type "

applications/solvers/combustion/PDRFoam/XiModels/XiGModels/XiGModel/XiGModel.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2015 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -67,7 +67,7 @@ Foam::XiGModel::~XiGModel()
 
 bool Foam::XiGModel::read(const dictionary& XiGProperties)
 {
-    XiGModelCoeffs_ = XiGProperties.subDict(type() + "Coeffs");
+    XiGModelCoeffs_ = XiGProperties.optionalSubDict(type() + "Coeffs");
 
     return true;
 }

applications/solvers/combustion/PDRFoam/XiModels/XiGModels/XiGModel/XiGModelNew.C

@@ -42,7 +42,7 @@ Foam::autoPtr<Foam::XiGModel> Foam::XiGModel::New
     dictionaryConstructorTable::iterator cstrIter =
         dictionaryConstructorTablePtr_->find(modelType);
 
-    if (cstrIter == dictionaryConstructorTablePtr_->end())
+    if (!cstrIter.found())
     {
         FatalErrorInFunction
             << "Unknown XiGModel type "

applications/solvers/combustion/PDRFoam/XiModels/XiModel/XiModel.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2015 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -85,7 +85,7 @@ Foam::XiModel::~XiModel()
 
 bool Foam::XiModel::read(const dictionary& XiProperties)
 {
-    XiModelCoeffs_ = XiProperties.subDict(type() + "Coeffs");
+    XiModelCoeffs_ = XiProperties.optionalSubDict(type() + "Coeffs");
 
     return true;
 }

applications/solvers/combustion/PDRFoam/XiModels/XiModel/XiModelNew.C

@@ -45,7 +45,7 @@ Foam::autoPtr<Foam::XiModel> Foam::XiModel::New
     dictionaryConstructorTable::iterator cstrIter =
         dictionaryConstructorTablePtr_->find(modelType);
 
-    if (cstrIter == dictionaryConstructorTablePtr_->end())
+    if (!cstrIter.found())
     {
         FatalErrorInFunction
             << "Unknown XiModel type "

applications/solvers/combustion/PDRFoam/laminarFlameSpeed/SCOPE/SCOPELaminarFlameSpeed.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -80,7 +80,7 @@ Foam::laminarFlameSpeedModels::SCOPE::SCOPE
                   dict.lookup("fuelFile")
               )
           )()
-        ).subDict(typeName + "Coeffs")
+        ).optionalSubDict(typeName + "Coeffs")
     ),
     LFL_(readScalar(coeffsDict_.lookup("lowerFlamabilityLimit"))),
     UFL_(readScalar(coeffsDict_.lookup("upperFlamabilityLimit"))),

applications/solvers/combustion/chemFoam/createBaseFields.H

@@ -54,4 +54,3 @@ Info<< "Creating base fields for time " << runTime.timeName() << endl;
 
     T.write();
 }
-

applications/solvers/combustion/chemFoam/output.H

@@ -1,6 +1,6 @@
     runTime.write();
 
-    Info<< "Sh = " << Sh
+    Info<< "Qdot = " << Qdot
         << ", T = " << thermo.T()[0]
         << ", p = " << thermo.p()[0]
         << ", " << Y[0].name() << " = " << Y[0][0]
@@ -8,4 +8,3 @@
 
     post<< runTime.value() << token::TAB << thermo.T()[0] << token::TAB
         << thermo.p()[0] << endl;
-

applications/solvers/combustion/chemFoam/setDeltaT.H

@@ -3,4 +3,3 @@ if (adjustTimeStep)
     runTime.setDeltaT(min(dtChem, maxDeltaT));
     Info<< "deltaT = " <<  runTime.deltaT().value() << endl;
 }
-

applications/solvers/combustion/chemFoam/solveChemistry.H

@@ -1,3 +1,3 @@
-    dtChem = chemistry.solve(runTime.deltaT().value());
+dtChem = chemistry.solve(runTime.deltaT().value());
-    scalar Sh = chemistry.Sh()()[0]/rho[0];
+scalar Qdot = chemistry.Qdot()()[0]/rho[0];
-    integratedHeat += Sh*runTime.deltaT().value();
+integratedHeat += Qdot*runTime.deltaT().value();

applications/solvers/combustion/fireFoam/Make/options

@@ -13,10 +13,7 @@ EXE_INC = \
     -I$(LIB_SRC)/thermophysicalModels/chemistryModel/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/solidChemistryModel/lnInclude \
     -I$(LIB_SRC)/combustionModels/lnInclude \
-    -I$(LIB_SRC)/thermophysicalModels/properties/solidProperties/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/thermophysicalProperties/lnInclude \
-    -I$(LIB_SRC)/thermophysicalModels/properties/solidMixtureProperties/lnInclude \
-    -I$(LIB_SRC)/thermophysicalModels/properties/liquidProperties/lnInclude \
-    -I$(LIB_SRC)/thermophysicalModels/properties/liquidMixtureProperties/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/thermophysicalFunctions/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/reactionThermo/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/SLGThermo/lnInclude \
@@ -38,9 +35,6 @@ EXE_LIBS = \
     -lspecie \
     -lcompressibleTransportModels \
     -lfluidThermophysicalModels \
-    -lsolidProperties \
-    -lsolidMixtureProperties \
-    -lthermophysicalFunctions \
     -lreactionThermophysicalModels \
     -lSLGThermo \
     -lchemistryModel \

applications/solvers/combustion/fireFoam/YEEqn.H

@@ -10,7 +10,7 @@ tmp<fv::convectionScheme<scalar>> mvConvection
 );
 {
     combustion->correct();
-    dQ = combustion->dQ();
+    Qdot = combustion->Qdot();
     volScalarField Yt(0.0*Y[0]);
 
     forAll(Y, i)
@@ -67,8 +67,8 @@ tmp<fv::convectionScheme<scalar>> mvConvection
         )
       - fvm::laplacian(turbulence->alphaEff(), he)
      ==
-        combustion->Sh()
+        Qdot
-      + radiation->Sh(thermo)
+      + radiation->Sh(thermo, he)
       + parcels.Sh(he)
       + surfaceFilm.Sh()
       + fvOptions(rho, he)

applications/solvers/combustion/fireFoam/createFields.H

@@ -131,18 +131,18 @@ Switch solvePyrolysisRegion
     additionalControlsDict.lookupOrDefault<bool>("solvePyrolysisRegion", true)
 );
 
-volScalarField dQ
+volScalarField Qdot
 (
     IOobject
     (
-        "dQ",
+        "Qdot",
         runTime.timeName(),
         mesh,
-        IOobject::NO_READ,
+        IOobject::READ_IF_PRESENT,
         IOobject::AUTO_WRITE
     ),
     mesh,
-    dimensionedScalar("dQ", dimEnergy/dimTime, 0.0)
+    dimensionedScalar("Qdot", dimEnergy/dimVolume/dimTime, 0.0)
 );
 
 

applications/solvers/combustion/reactingFoam/EEqn.H

@@ -17,7 +17,7 @@
         )
       - fvm::laplacian(turbulence->alphaEff(), he)
      ==
-        reaction->Sh()
+        Qdot
       + fvOptions(rho, he)
     );
 

applications/solvers/combustion/reactingFoam/YEqn.H

@@ -11,7 +11,7 @@ tmp<fv::convectionScheme<scalar>> mvConvection
 
 {
     reaction->correct();
-    dQ = reaction->dQ();
+    Qdot = reaction->Qdot();
     volScalarField Yt(0.0*Y[0]);
 
     forAll(Y, i)

applications/solvers/combustion/reactingFoam/createFields.H

@@ -117,18 +117,18 @@ forAll(Y, i)
 }
 fields.add(thermo.he());
 
-volScalarField dQ
+volScalarField Qdot
 (
     IOobject
     (
-        "dQ",
+        "Qdot",
         runTime.timeName(),
         mesh,
-        IOobject::NO_READ,
+        IOobject::READ_IF_PRESENT,
         IOobject::AUTO_WRITE
     ),
     mesh,
-    dimensionedScalar("dQ", dimEnergy/dimTime, 0.0)
+    dimensionedScalar("Qdot", dimEnergy/dimVolume/dimTime, 0.0)
 );
 
 #include "createMRF.H"

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

@@ -117,18 +117,18 @@ forAll(Y, i)
 }
 fields.add(thermo.he());
 
-volScalarField dQ
+volScalarField Qdot
 (
     IOobject
     (
-        "dQ",
+        "Qdot",
         runTime.timeName(),
         mesh,
-        IOobject::NO_READ,
+        IOobject::READ_IF_PRESENT,
         IOobject::AUTO_WRITE
     ),
     mesh,
-    dimensionedScalar("dQ", dimEnergy/dimTime, 0.0)
+    dimensionedScalar("Qdot", dimEnergy/dimVolume/dimTime, 0.0)
 );
 
 #include "createMRF.H"

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

@@ -96,18 +96,18 @@ forAll(Y, i)
 }
 fields.add(thermo.he());
 
-volScalarField dQ
+volScalarField Qdot
 (
     IOobject
     (
-        "dQ",
+        "Qdot",
         runTime.timeName(),
         mesh,
-        IOobject::NO_READ,
+        IOobject::READ_IF_PRESENT,
         IOobject::AUTO_WRITE
     ),
     mesh,
-    dimensionedScalar("dQ", dimEnergy/dimTime, 0.0)
+    dimensionedScalar("Qdot", dimEnergy/dimVolume/dimTime, 0.0)
 );
 
 #include "createMRF.H"

applications/solvers/combustion/reactingFoam/setRDeltaT.H

@@ -43,13 +43,16 @@ License
     // Damping coefficient (1-0)
     scalar rDeltaTDampingCoeff
     (
-        pimpleDict.lookupOrDefault<scalar>("rDeltaTDampingCoeff", 1)
+        pimpleDict.lookupOrDefault<scalar>("rDeltaTDampingCoeff", 1.0)
     );
 
     // Maximum change in cell temperature per iteration
     // (relative to previous value)
     scalar alphaTemp(pimpleDict.lookupOrDefault("alphaTemp", 0.05));
 
+    // Maximum change in cell concentration per iteration
+    // (relative to reference value)
+    scalar alphaY(pimpleDict.lookupOrDefault("alphaY", 1.0));
 
     Info<< "Time scales min/max:" << endl;
 
@@ -68,34 +71,89 @@ License
         rDeltaT.max(1/maxDeltaT);
 
         Info<< "    Flow        = "
-            << gMin(1/rDeltaT.primitiveField()) << ", "
+            << 1/gMax(rDeltaT.primitiveField()) << ", "
-            << gMax(1/rDeltaT.primitiveField()) << endl;
+            << 1/gMin(rDeltaT.primitiveField()) << endl;
     }
 
-    // Reaction source time scale
+    // Heat release rate time scale
-    if (alphaTemp < 1.0)
+    if (alphaTemp < 1)
     {
         volScalarField::Internal rDeltaTT
         (
-            mag(reaction->Sh())/(alphaTemp*rho*thermo.Cp()*T)
+            mag(Qdot)/(alphaTemp*rho*thermo.Cp()*T)
         );
 
         Info<< "    Temperature = "
-            << gMin(1/(rDeltaTT.field() + VSMALL)) << ", "
+            << 1/(gMax(rDeltaTT.field()) + VSMALL) << ", "
-            << gMax(1/(rDeltaTT.field() + VSMALL)) << endl;
+            << 1/(gMin(rDeltaTT.field()) + VSMALL) << endl;
 
-        rDeltaT.ref() = max
+        rDeltaT.ref() = max(rDeltaT(), rDeltaTT);
+    }
+
+    // Reaction rate time scale
+    if (alphaY < 1)
+    {
+        dictionary Yref(pimpleDict.subDict("Yref"));
+
+        volScalarField::Internal rDeltaTY
         (
-            rDeltaT(),
+            IOobject
-            rDeltaTT
+            (
+                "rDeltaTY",
+                runTime.timeName(),
+                mesh
+            ),
+            mesh,
+            dimensionedScalar("rDeltaTY", rDeltaT.dimensions(), 0)
         );
+
+        bool foundY = false;
+
+        forAll(Y, i)
+        {
+            if (i != inertIndex && composition.active(i))
+            {
+                volScalarField& Yi = Y[i];
+
+                if (Yref.found(Yi.name()))
+                {
+                    foundY = true;
+                    scalar Yrefi = readScalar(Yref.lookup(Yi.name()));
+
+                    rDeltaTY.field() = max
+                    (
+                        mag
+                        (
+                            reaction->R(Yi)().source()
+                           /((Yrefi*alphaY)*(rho*mesh.V()))
+                        ),
+                        rDeltaTY
+                    );
+                }
+            }
+        }
+
+        if (foundY)
+        {
+            Info<< "    Composition = "
+                << 1/(gMax(rDeltaTY.field()) + VSMALL) << ", "
+                << 1/(gMin(rDeltaTY.field()) + VSMALL) << endl;
+
+            rDeltaT.ref() = max(rDeltaT(), rDeltaTY);
+        }
+        else
+        {
+            IOWarningIn(args.executable().c_str(), Yref)
+                << "Cannot find any active species in Yref " << Yref
+                << endl;
+        }
     }
 
     // Update tho boundary values of the reciprocal time-step
     rDeltaT.correctBoundaryConditions();
 
     // Spatially smooth the time scale field
-    if (rDeltaTSmoothingCoeff < 1.0)
+    if (rDeltaTSmoothingCoeff < 1)
     {
         fvc::smooth(rDeltaT, rDeltaTSmoothingCoeff);
     }
@@ -105,7 +163,7 @@ License
     // - only increase at a fraction of old time scale
     if
     (
-        rDeltaTDampingCoeff < 1.0
+        rDeltaTDampingCoeff < 1
      && runTime.timeIndex() > runTime.startTimeIndex() + 1
     )
     {
@@ -120,8 +178,8 @@ License
     rDeltaT.correctBoundaryConditions();
 
     Info<< "    Overall     = "
-        << gMin(1/rDeltaT.primitiveField())
+        << 1/gMax(rDeltaT.primitiveField())
-        << ", " << gMax(1/rDeltaT.primitiveField()) << endl;
+        << ", " << 1/gMin(rDeltaT.primitiveField()) << endl;
 }
 
 

applications/solvers/compressible/rhoCentralFoam/BCs/mixedFixedValueSlip/mixedFixedValueSlipFvPatchField.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -184,7 +184,8 @@ public:
             //- Evaluate the patch field
             virtual void evaluate
             (
-                const Pstream::commsTypes commsType=Pstream::blocking
+                const Pstream::commsTypes commsType=
+                    Pstream::commsTypes::blocking
             );
 
             //- Return face-gradient transform diagonal

applications/solvers/compressible/rhoPimpleFoam/createFields.H

@@ -2,11 +2,11 @@
 
 Info<< "Reading thermophysical properties\n" << endl;
 
-autoPtr<psiThermo> pThermo
+autoPtr<fluidThermo> pThermo
 (
-    psiThermo::New(mesh)
+    fluidThermo::New(mesh)
 );
-psiThermo& thermo = pThermo();
+fluidThermo& thermo = pThermo();
 thermo.validate(args.executable(), "h", "e");
 
 volScalarField& p = thermo.p();
@@ -40,27 +40,7 @@ volVectorField U
 
 #include "compressibleCreatePhi.H"
 
-dimensionedScalar rhoMax
+pressureControl pressureControl(p, rho, pimple.dict(), false);
-(
-    dimensionedScalar::lookupOrDefault
-    (
-        "rhoMax",
-        pimple.dict(),
-        dimDensity,
-        GREAT
-    )
-);
-
-dimensionedScalar rhoMin
-(
-    dimensionedScalar::lookupOrDefault
-    (
-        "rhoMin",
-        pimple.dict(),
-        dimDensity,
-        0
-    )
-);
 
 Info<< "Creating turbulence model\n" << endl;
 autoPtr<compressible::turbulenceModel> turbulence

applications/solvers/compressible/rhoPimpleFoam/pEqn.H

@@ -1,8 +1,3 @@
-rho = thermo.rho();
-rho = max(rho, rhoMin);
-rho = min(rho, rhoMax);
-rho.relax();
-
 volScalarField rAU(1.0/UEqn.A());
 surfaceScalarField rhorAUf("rhorAUf", fvc::interpolate(rho*rAU));
 volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p));
@@ -12,55 +7,52 @@ if (pimple.nCorrPISO() <= 1)
     tUEqn.clear();
 }
 
+surfaceScalarField phiHbyA
+(
+    "phiHbyA",
+    fvc::flux(rho*HbyA)
+  + rhorAUf*fvc::ddtCorr(rho, U, phi)
+);
+
+MRF.makeRelative(fvc::interpolate(rho), phiHbyA);
+
+// Update the pressure BCs to ensure flux consistency
+constrainPressure(p, rho, U, phiHbyA, rhorAUf, MRF);
+
 if (pimple.transonic())
 {
     surfaceScalarField phid
     (
         "phid",
-        fvc::interpolate(psi)
+        (fvc::interpolate(psi)/fvc::interpolate(rho))*phiHbyA
-       *(
-            fvc::flux(HbyA)
-          + rhorAUf*fvc::ddtCorr(rho, U, phi)/fvc::interpolate(rho)
-        )
     );
-
+    phiHbyA -= fvc::interpolate(p)*phid;
-    MRF.makeRelative(fvc::interpolate(psi), phid);
 
     while (pimple.correctNonOrthogonal())
     {
         fvScalarMatrix pEqn
         (
             fvm::ddt(psi, p)
+          + fvc::div(phiHbyA)
           + fvm::div(phid, p)
           - fvm::laplacian(rhorAUf, p)
          ==
             fvOptions(psi, p, rho.name())
         );
 
+        // Relax the pressure equation to ensure diagonal-dominance
+        pEqn.relax();
+
         pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter())));
 
         if (pimple.finalNonOrthogonalIter())
         {
-            phi == pEqn.flux();
+            phi = phiHbyA + pEqn.flux();
         }
     }
 }
 else
 {
-    surfaceScalarField phiHbyA
-    (
-        "phiHbyA",
-        (
-            fvc::flux(rho*HbyA)
-          + rhorAUf*fvc::ddtCorr(rho, U, phi)
-        )
-    );
-
-    MRF.makeRelative(fvc::interpolate(rho), phiHbyA);
-
-    // Update the pressure BCs to ensure flux consistency
-    constrainPressure(p, rho, U, phiHbyA, rhorAUf, MRF);
-
     while (pimple.correctNonOrthogonal())
     {
         fvScalarMatrix pEqn
@@ -87,19 +79,20 @@ else
 // Explicitly relax pressure for momentum corrector
 p.relax();
 
-// Recalculate density from the relaxed pressure
-rho = thermo.rho();
-rho = max(rho, rhoMin);
-rho = min(rho, rhoMax);
-rho.relax();
-Info<< "rho max/min : " << max(rho).value()
-    << " " << min(rho).value() << endl;
-
 U = HbyA - rAU*fvc::grad(p);
 U.correctBoundaryConditions();
 fvOptions.correct(U);
 K = 0.5*magSqr(U);
 
+pressureControl.limit(p);
+p.correctBoundaryConditions();
+rho = thermo.rho();
+
+if (!pimple.transonic())
+{
+    rho.relax();
+}
+
 if (thermo.dpdt())
 {
     dpdt = fvc::ddt(p);

applications/solvers/compressible/rhoPimpleFoam/pcEqn.H

@@ -1,8 +1,3 @@
-rho = thermo.rho();
-rho = max(rho, rhoMin);
-rho = min(rho, rhoMax);
-rho.relax();
-
 volScalarField rAU(1.0/UEqn.A());
 volScalarField rAtU(1.0/(1.0/rAU - UEqn.H1()));
 volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p));
@@ -12,72 +7,64 @@ if (pimple.nCorrPISO() <= 1)
     tUEqn.clear();
 }
 
+surfaceScalarField phiHbyA
+(
+    "phiHbyA",
+    (
+        fvc::flux(rho*HbyA)
+      + fvc::interpolate(rho*rAU)*fvc::ddtCorr(rho, U, phi)
+    )
+);
+
+MRF.makeRelative(fvc::interpolate(rho), phiHbyA);
+
+volScalarField rhorAtU("rhorAtU", rho*rAtU);
+
+// Update the pressure BCs to ensure flux consistency
+constrainPressure(p, rho, U, phiHbyA, rhorAtU, MRF);
+
 if (pimple.transonic())
 {
     surfaceScalarField phid
     (
         "phid",
-        fvc::interpolate(psi)
+        (fvc::interpolate(psi)/fvc::interpolate(rho))*phiHbyA
-       *(
-            fvc::flux(HbyA)
-          + fvc::interpolate(rho*rAU)*fvc::ddtCorr(rho, U, phi)
-           /fvc::interpolate(rho)
-        )
     );
 
-    MRF.makeRelative(fvc::interpolate(psi), phid);
+    phiHbyA +=
-
-    surfaceScalarField phic
-    (
-        "phic",
         fvc::interpolate(rho*(rAtU - rAU))*fvc::snGrad(p)*mesh.magSf()
-    );
+      - fvc::interpolate(p)*phid;
 
     HbyA -= (rAU - rAtU)*fvc::grad(p);
 
-    volScalarField rhorAtU("rhorAtU", rho*rAtU);
-
     while (pimple.correctNonOrthogonal())
     {
         fvScalarMatrix pEqn
         (
             fvm::ddt(psi, p)
+          + fvc::div(phiHbyA)
           + fvm::div(phid, p)
-          + fvc::div(phic)
           - fvm::laplacian(rhorAtU, p)
          ==
             fvOptions(psi, p, rho.name())
         );
 
+        // Relax the pressure equation to ensure diagonal-dominance
+        pEqn.relax();
+
         pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter())));
 
         if (pimple.finalNonOrthogonalIter())
         {
-            phi == phic + pEqn.flux();
+            phi = phiHbyA + pEqn.flux();
         }
     }
 }
 else
 {
-    surfaceScalarField phiHbyA
-    (
-        "phiHbyA",
-        (
-            fvc::flux(rho*HbyA)
-          + fvc::interpolate(rho*rAU)*fvc::ddtCorr(rho, U, phi)
-        )
-    );
-
-    MRF.makeRelative(fvc::interpolate(rho), phiHbyA);
-
     phiHbyA += fvc::interpolate(rho*(rAtU - rAU))*fvc::snGrad(p)*mesh.magSf();
     HbyA -= (rAU - rAtU)*fvc::grad(p);
 
-    volScalarField rhorAtU("rhorAtU", rho*rAtU);
-
-    // Update the pressure BCs to ensure flux consistency
-    constrainPressure(p, rho, U, phiHbyA, rhorAtU, MRF);
-
     while (pimple.correctNonOrthogonal())
     {
         fvScalarMatrix pEqn
@@ -109,19 +96,16 @@ U.correctBoundaryConditions();
 fvOptions.correct(U);
 K = 0.5*magSqr(U);
 
-if (thermo.dpdt())
+pressureControl.limit(p);
-{
+p.correctBoundaryConditions();
-    dpdt = fvc::ddt(p);
-}
-
-// Recalculate density from the relaxed pressure
 rho = thermo.rho();
-rho = max(rho, rhoMin);
-rho = min(rho, rhoMax);
 
 if (!pimple.transonic())
 {
     rho.relax();
 }
 
-Info<< "rho max/min : " << max(rho).value() << " " << min(rho).value() << endl;
+if (thermo.dpdt())
+{
+    dpdt = fvc::ddt(p);
+}

applications/solvers/compressible/rhoPimpleFoam/rhoPimpleDyMFoam/pEqn.H

@@ -1,8 +1,3 @@
-rho = thermo.rho();
-rho = max(rho, rhoMin);
-rho = min(rho, rhoMax);
-rho.relax();
-
 volScalarField rAU(1.0/UEqn.A());
 surfaceScalarField rhorAUf("rhorAUf", fvc::interpolate(rho*rAU));
 volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p));
@@ -12,55 +7,53 @@ if (pimple.nCorrPISO() <= 1)
     tUEqn.clear();
 }
 
+surfaceScalarField phiHbyA
+(
+    "phiHbyA",
+    fvc::flux(rho*HbyA)
+  + rhorAUf*fvc::ddtCorr(rho, U, rhoUf)
+);
+
+fvc::makeRelative(phiHbyA, rho, U);
+MRF.makeRelative(fvc::interpolate(rho), phiHbyA);
+
+// Update the pressure BCs to ensure flux consistency
+constrainPressure(p, rho, U, phiHbyA, rhorAUf, MRF);
+
 if (pimple.transonic())
 {
     surfaceScalarField phid
     (
         "phid",
-        fvc::interpolate(psi)
+        (fvc::interpolate(psi)/fvc::interpolate(rho))*phiHbyA
-       *(
-            fvc::flux(HbyA)
-          + rhorAUf*fvc::ddtCorr(rho, U, rhoUf)/fvc::interpolate(rho)
-        )
     );
-
+    phiHbyA -= fvc::interpolate(p)*phid;
-    fvc::makeRelative(phid, psi, U);
-    MRF.makeRelative(fvc::interpolate(psi), phid);
 
     while (pimple.correctNonOrthogonal())
     {
         fvScalarMatrix pEqn
         (
             fvm::ddt(psi, p)
+          + fvc::div(phiHbyA)
           + fvm::div(phid, p)
           - fvm::laplacian(rhorAUf, p)
          ==
             fvOptions(psi, p, rho.name())
         );
 
+        // Relax the pressure equation to ensure diagonal-dominance
+        pEqn.relax();
+
         pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter())));
 
         if (pimple.finalNonOrthogonalIter())
         {
-            phi == pEqn.flux();
+            phi = phiHbyA + pEqn.flux();
         }
     }
 }
 else
 {
-    surfaceScalarField phiHbyA
-    (
-        "phiHbyA",
-        fvc::flux(rho*HbyA)
-      + rhorAUf*fvc::ddtCorr(rho, U, rhoUf)
-    );
-
-    fvc::makeRelative(phiHbyA, rho, U);
-    MRF.makeRelative(fvc::interpolate(rho), phiHbyA);
-
-    // Update the pressure BCs to ensure flux consistency
-    constrainPressure(p, rho, U, phiHbyA, rhorAUf, MRF);
-
     while (pimple.correctNonOrthogonal())
     {
         // Pressure corrector
@@ -88,19 +81,20 @@ else
 // Explicitly relax pressure for momentum corrector
 p.relax();
 
-// Recalculate density from the relaxed pressure
-rho = thermo.rho();
-rho = max(rho, rhoMin);
-rho = min(rho, rhoMax);
-rho.relax();
-Info<< "rho max/min : " << max(rho).value()
-    << " " << min(rho).value() << endl;
-
 U = HbyA - rAU*fvc::grad(p);
 U.correctBoundaryConditions();
 fvOptions.correct(U);
 K = 0.5*magSqr(U);
 
+pressureControl.limit(p);
+p.correctBoundaryConditions();
+rho = thermo.rho();
+
+if (!pimple.transonic())
+{
+    rho.relax();
+}
+
 {
     rhoUf = fvc::interpolate(rho*U);
     surfaceVectorField n(mesh.Sf()/mesh.magSf());

applications/solvers/compressible/rhoPimpleFoam/rhoPimpleDyMFoam/rhoPimpleDyMFoam.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -22,7 +22,7 @@ License
     along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
 
 Application
-    rhoPimpleFoam
+    rhoPimpleDyMFoam
 
 Group
     grpCompressibleSolvers grpMovingMeshSolvers
@@ -38,10 +38,11 @@ Description
 
 #include "fvCFD.H"
 #include "dynamicFvMesh.H"
-#include "psiThermo.H"
+#include "fluidThermo.H"
 #include "turbulentFluidThermoModel.H"
 #include "bound.H"
 #include "pimpleControl.H"
+#include "pressureControl.H"
 #include "CorrectPhi.H"
 #include "fvOptions.H"
 #include "localEulerDdtScheme.H"

applications/solvers/compressible/rhoPimpleFoam/rhoPimpleFoam.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -37,10 +37,11 @@ Description
 \*---------------------------------------------------------------------------*/
 
 #include "fvCFD.H"
-#include "psiThermo.H"
+#include "fluidThermo.H"
 #include "turbulentFluidThermoModel.H"
 #include "bound.H"
 #include "pimpleControl.H"
+#include "pressureControl.H"
 #include "fvOptions.H"
 #include "localEulerDdtScheme.H"
 #include "fvcSmooth.H"

applications/solvers/compressible/rhoSimpleFoam/createFields.H

@@ -1,11 +1,14 @@
 Info<< "Reading thermophysical properties\n" << endl;
-autoPtr<rhoThermo> pThermo
+
+autoPtr<fluidThermo> pThermo
 (
-    rhoThermo::New(mesh)
+    fluidThermo::New(mesh)
 );
-rhoThermo& thermo = pThermo();
+fluidThermo& thermo = pThermo();
 thermo.validate(args.executable(), "h", "e");
 
+volScalarField& p = thermo.p();
+
 volScalarField rho
 (
     IOobject
@@ -19,8 +22,6 @@ volScalarField rho
     thermo.rho()
 );
 
-volScalarField& p = thermo.p();
-
 Info<< "Reading field U\n" << endl;
 volVectorField U
 (
@@ -37,35 +38,10 @@ volVectorField U
 
 #include "compressibleCreatePhi.H"
 
-
+pressureControl pressureControl(p, rho, simple.dict());
-label pRefCell = 0;
-scalar pRefValue = 0.0;
-setRefCell(p, simple.dict(), pRefCell, pRefValue);
 
 mesh.setFluxRequired(p.name());
 
-dimensionedScalar rhoMax
-(
-    dimensionedScalar::lookupOrDefault
-    (
-        "rhoMax",
-        simple.dict(),
-        dimDensity,
-        GREAT
-    )
-);
-
-dimensionedScalar rhoMin
-(
-    dimensionedScalar::lookupOrDefault
-    (
-        "rhoMin",
-        simple.dict(),
-        dimDensity,
-        0
-    )
-);
-
 Info<< "Creating turbulence model\n" << endl;
 autoPtr<compressible::turbulenceModel> turbulence
 (

applications/solvers/compressible/rhoSimpleFoam/pEqn.H

@@ -1,6 +1,4 @@
 {
-    //const volScalarField& psi = thermo.psi();
-
     volScalarField rAU(1.0/UEqn.A());
     surfaceScalarField rhorAUf("rhorAUf", fvc::interpolate(rho*rAU));
     volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p));
@@ -8,22 +6,27 @@
 
     bool closedVolume = false;
 
+    surfaceScalarField phiHbyA("phiHbyA", fvc::flux(rho*HbyA));
+    MRF.makeRelative(fvc::interpolate(rho), phiHbyA);
+
+    // Update the pressure BCs to ensure flux consistency
+    constrainPressure(p, rho, U, phiHbyA, rhorAUf, MRF);
+
     if (simple.transonic())
     {
         surfaceScalarField phid
         (
             "phid",
-            fvc::interpolate(psi)
+            (fvc::interpolate(psi)/fvc::interpolate(rho))*phiHbyA
-           *fvc::flux(HbyA)
         );
-
+        phiHbyA -= fvc::interpolate(p)*phid;
-        MRF.makeRelative(fvc::interpolate(psi), phid);
 
         while (simple.correctNonOrthogonal())
         {
             fvScalarMatrix pEqn
             (
-                fvm::div(phid, p)
+                fvc::div(phiHbyA)
+              + fvm::div(phid, p)
               - fvm::laplacian(rhorAUf, p)
               ==
                 fvOptions(psi, p, rho.name())
@@ -32,26 +35,24 @@
             // Relax the pressure equation to ensure diagonal-dominance
             pEqn.relax();
 
-            pEqn.setReference(pRefCell, pRefValue);
+            pEqn.setReference
+            (
+                pressureControl.refCell(),
+                pressureControl.refValue()
+            );
 
             pEqn.solve();
 
             if (simple.finalNonOrthogonalIter())
             {
-                phi == pEqn.flux();
+                phi = phiHbyA + pEqn.flux();
             }
         }
     }
     else
     {
-        surfaceScalarField phiHbyA("phiHbyA", fvc::flux(rho*HbyA));
-        MRF.makeRelative(fvc::interpolate(rho), phiHbyA);
-
         closedVolume = adjustPhi(phiHbyA, U, p);
 
-        // Update the pressure BCs to ensure flux consistency
-        constrainPressure(p, rho, U, phiHbyA, rhorAUf, MRF);
-
         while (simple.correctNonOrthogonal())
         {
             fvScalarMatrix pEqn
@@ -62,7 +63,11 @@
                 fvOptions(psi, p, rho.name())
             );
 
-            pEqn.setReference(pRefCell, pRefValue);
+            pEqn.setReference
+            (
+                pressureControl.refCell(),
+                pressureControl.refValue()
+            );
 
             pEqn.solve();
 
@@ -83,6 +88,8 @@
     U.correctBoundaryConditions();
     fvOptions.correct(U);
 
+    pressureControl.limit(p);
+
     // For closed-volume cases adjust the pressure and density levels
     // to obey overall mass continuity
     if (closedVolume)
@@ -91,19 +98,12 @@
             /fvc::domainIntegrate(psi);
     }
 
-    rho = thermo.rho();
+    p.correctBoundaryConditions();
-    rho = max(rho, rhoMin);
-    rho = min(rho, rhoMax);
 
-    thermo.rho() = max(thermo.rho(), rhoMin);
+    rho = thermo.rho();
-    thermo.rho() = min(thermo.rho(), rhoMax);
 
     if (!simple.transonic())
     {
         rho.relax();
     }
-
-    Info<< "rho max/min : "
-        << max(rho).value() << " "
-        << min(rho).value() << endl;
 }

applications/solvers/compressible/rhoSimpleFoam/pcEqn.H

@@ -5,33 +5,34 @@ tUEqn.clear();
 
 bool closedVolume = false;
 
+surfaceScalarField phiHbyA("phiHbyA", fvc::flux(rho*HbyA));
+MRF.makeRelative(fvc::interpolate(rho), phiHbyA);
+
+volScalarField rhorAtU("rhorAtU", rho*rAtU);
+
+// Update the pressure BCs to ensure flux consistency
+constrainPressure(p, rho, U, phiHbyA, rhorAtU, MRF);
+
 if (simple.transonic())
 {
     surfaceScalarField phid
     (
         "phid",
-        fvc::interpolate(psi)
+        (fvc::interpolate(psi)/fvc::interpolate(rho))*phiHbyA
-       *fvc::flux(HbyA)
     );
 
-    MRF.makeRelative(fvc::interpolate(psi), phid);
+    phiHbyA +=
-
-    surfaceScalarField phic
-    (
-        "phic",
         fvc::interpolate(rho*(rAtU - rAU))*fvc::snGrad(p)*mesh.magSf()
-    );
+      - fvc::interpolate(p)*phid;
 
     HbyA -= (rAU - rAtU)*fvc::grad(p);
 
-    volScalarField rhorAtU("rhorAtU", rho*rAtU);
-
     while (simple.correctNonOrthogonal())
     {
         fvScalarMatrix pEqn
         (
-            fvm::div(phid, p)
+            fvc::div(phiHbyA)
-          + fvc::div(phic)
+          + fvm::div(phid, p)
           - fvm::laplacian(rhorAtU, p)
          ==
             fvOptions(psi, p, rho.name())
@@ -40,31 +41,27 @@ if (simple.transonic())
         // Relax the pressure equation to maintain diagonal dominance
         pEqn.relax();
 
-        pEqn.setReference(pRefCell, pRefValue);
+        pEqn.setReference
+        (
+            pressureControl.refCell(),
+            pressureControl.refValue()
+        );
 
         pEqn.solve();
 
         if (simple.finalNonOrthogonalIter())
         {
-            phi == phic + pEqn.flux();
+            phi = phiHbyA + pEqn.flux();
         }
     }
 }
 else
 {
-    surfaceScalarField phiHbyA("phiHbyA", fvc::flux(rho*HbyA));
-    MRF.makeRelative(fvc::interpolate(rho), phiHbyA);
-
     closedVolume = adjustPhi(phiHbyA, U, p);
 
     phiHbyA += fvc::interpolate(rho*(rAtU - rAU))*fvc::snGrad(p)*mesh.magSf();
     HbyA -= (rAU - rAtU)*fvc::grad(p);
 
-    volScalarField rhorAtU("rhorAtU", rho*rAtU);
-
-    // Update the pressure BCs to ensure flux consistency
-    constrainPressure(p, rho, U, phiHbyA, rhorAtU, MRF);
-
     while (simple.correctNonOrthogonal())
     {
         fvScalarMatrix pEqn
@@ -75,7 +72,11 @@ else
             fvOptions(psi, p, rho.name())
         );
 
-        pEqn.setReference(pRefCell, pRefValue);
+        pEqn.setReference
+        (
+            pressureControl.refCell(),
+            pressureControl.refValue()
+        );
 
         pEqn.solve();
 
@@ -97,6 +98,8 @@ U = HbyA - rAtU*fvc::grad(p);
 U.correctBoundaryConditions();
 fvOptions.correct(U);
 
+pressureControl.limit(p);
+
 // For closed-volume cases adjust the pressure and density levels
 // to obey overall mass continuity
 if (closedVolume)
@@ -105,16 +108,12 @@ if (closedVolume)
         /fvc::domainIntegrate(psi);
 }
 
+p.correctBoundaryConditions();
+
 // Recalculate density from the relaxed pressure
 rho = thermo.rho();
-rho = max(rho, rhoMin);
-rho = min(rho, rhoMax);
-thermo.rho() = max(thermo.rho(), rhoMin);
-thermo.rho() = min(thermo.rho(), rhoMax);
 
 if (!simple.transonic())
 {
     rho.relax();
 }
-
-Info<< "rho max/min : " << max(rho).value() << " " << min(rho).value() << endl;

applications/solvers/compressible/rhoSimpleFoam/rhoPorousSimpleFoam/EEqn.H

@@ -1,26 +0,0 @@
-{
-    volScalarField& he = thermo.he();
-
-    fvScalarMatrix EEqn
-    (
-        fvm::div(phi, he)
-      + (
-            he.name() == "e"
-          ? fvc::div(phi, volScalarField("Ekp", 0.5*magSqr(U) + p/rho))
-          : fvc::div(phi, volScalarField("K", 0.5*magSqr(U)))
-        )
-      - fvm::laplacian(turbulence->alphaEff(), he)
-     ==
-        fvOptions(rho, he)
-    );
-
-    EEqn.relax();
-
-    fvOptions.constrain(EEqn);
-
-    EEqn.solve();
-
-    fvOptions.correct(he);
-
-    thermo.correct();
-}

applications/solvers/compressible/rhoSimpleFoam/rhoPorousSimpleFoam/createFields.H

@@ -1,84 +0,0 @@
-Info<< "Reading thermophysical properties\n" << endl;
-
-autoPtr<rhoThermo> pThermo
-(
-    rhoThermo::New(mesh)
-);
-rhoThermo& thermo = pThermo();
-thermo.validate(args.executable(), "h", "e");
-
-volScalarField rho
-(
-    IOobject
-    (
-        "rho",
-        runTime.timeName(),
-        mesh,
-        IOobject::READ_IF_PRESENT,
-        IOobject::AUTO_WRITE
-    ),
-    thermo.rho()
-);
-
-volScalarField& p = thermo.p();
-
-Info<< "Reading field U\n" << endl;
-volVectorField U
-(
-    IOobject
-    (
-        "U",
-        runTime.timeName(),
-        mesh,
-        IOobject::MUST_READ,
-        IOobject::AUTO_WRITE
-    ),
-    mesh
-);
-
-#include "compressibleCreatePhi.H"
-
-
-label pRefCell = 0;
-scalar pRefValue = 0.0;
-setRefCell(p, simple.dict(), pRefCell, pRefValue);
-
-mesh.setFluxRequired(p.name());
-
-dimensionedScalar rhoMax
-(
-    dimensionedScalar::lookupOrDefault
-    (
-        "rhoMax",
-        simple.dict(),
-        dimDensity,
-        GREAT
-    )
-);
-
-dimensionedScalar rhoMin
-(
-    dimensionedScalar::lookupOrDefault
-    (
-        "rhoMin",
-        simple.dict(),
-        dimDensity,
-        0
-    )
-);
-
-Info<< "Creating turbulence model\n" << endl;
-autoPtr<compressible::turbulenceModel> turbulence
-(
-    compressible::turbulenceModel::New
-    (
-        rho,
-        U,
-        phi,
-        thermo
-    )
-);
-
-dimensionedScalar initialMass = fvc::domainIntegrate(rho);
-
-#include "createMRF.H"

applications/solvers/compressible/rhoSimpleFoam/rhoPorousSimpleFoam/pEqn.H

@@ -48,7 +48,11 @@
 
         fvScalarMatrix& pEqn = tpEqn.ref();
 
-        pEqn.setReference(pRefCell, pRefValue);
+        pEqn.setReference
+        (
+            pressureControl.refCell(),
+            pressureControl.refValue()
+        );
 
         pEqn.solve();
 
@@ -75,6 +79,8 @@
     U.correctBoundaryConditions();
     fvOptions.correct(U);
 
+    pressureControl.limit(p);
+
     // For closed-volume cases adjust the pressure and density levels
     // to obey overall mass continuity
     if (closedVolume)
@@ -84,14 +90,5 @@
     }
 
     rho = thermo.rho();
-    rho = max(rho, rhoMin);
-    rho = min(rho, rhoMax);
-
-    thermo.rho() = max(thermo.rho(), rhoMin);
-    thermo.rho() = min(thermo.rho(), rhoMax);
-
     rho.relax();
-    Info<< "rho max/min : "
-        << max(rho).value() << " "
-        << min(rho).value() << endl;
 }

applications/solvers/compressible/rhoSimpleFoam/rhoPorousSimpleFoam/rhoPorousSimpleFoam.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -34,11 +34,12 @@ Description
 \*---------------------------------------------------------------------------*/
 
 #include "fvCFD.H"
-#include "rhoThermo.H"
+#include "fluidThermo.H"
 #include "turbulentFluidThermoModel.H"
+#include "simpleControl.H"
+#include "pressureControl.H"
 #include "fvOptions.H"
 #include "IOporosityModelList.H"
-#include "simpleControl.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 

applications/solvers/compressible/rhoSimpleFoam/rhoSimpleFoam.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -33,9 +33,10 @@ Description
 \*---------------------------------------------------------------------------*/
 
 #include "fvCFD.H"
-#include "rhoThermo.H"
+#include "fluidThermo.H"
 #include "turbulentFluidThermoModel.H"
 #include "simpleControl.H"
+#include "pressureControl.H"
 #include "fvOptions.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

applications/solvers/compressible/sonicFoam/sonicLiquidFoam/readThermodynamicProperties.H

@@ -1,37 +1,37 @@
-    Info<< "Reading thermodynamicProperties\n" << endl;
+Info<< "Reading thermodynamicProperties\n" << endl;
 
-    IOdictionary thermodynamicProperties
+IOdictionary thermodynamicProperties
+(
+    IOobject
     (
-        IOobject
+        "thermodynamicProperties",
-        (
+        runTime.constant(),
-            "thermodynamicProperties",
+        mesh,
-            runTime.constant(),
+        IOobject::MUST_READ_IF_MODIFIED,
-            mesh,
+        IOobject::NO_WRITE
-            IOobject::MUST_READ_IF_MODIFIED,
+    )
-            IOobject::NO_WRITE
+);
-        )
-    );
 
-    dimensionedScalar rho0
+dimensionedScalar rho0
-    (
+(
-        "rho0",
+    "rho0",
-        dimDensity,
+    dimDensity,
-        thermodynamicProperties
+    thermodynamicProperties
-    );
+);
 
-    dimensionedScalar p0
+dimensionedScalar p0
-    (
+(
-        "p0",
+    "p0",
-        dimPressure,
+    dimPressure,
-        thermodynamicProperties
+    thermodynamicProperties
-    );
+);
 
-    dimensionedScalar psi
+dimensionedScalar psi
-    (
+(
-        "psi",
+    "psi",
-        dimCompressibility,
+    dimCompressibility,
-        thermodynamicProperties
+    thermodynamicProperties
-    );
+);
 
-    // Density offset, i.e. the constant part of the density
+// Density offset, i.e. the constant part of the density
-    dimensionedScalar rhoO("rhoO", rho0 - psi*p0);
+dimensionedScalar rhoO("rhoO", rho0 - psi*p0);

applications/solvers/electromagnetics/mhdFoam/readBPISOControls.H

@@ -1,4 +1,3 @@
     const dictionary& Bpiso = mesh.solutionDict().subDict("BPISO");
 
     const int nBcorr = Bpiso.lookupOrDefault<int>("nCorrectors", 1);
-

applications/solvers/heatTransfer/buoyantPimpleFoam/EEqn.H

@@ -18,7 +18,7 @@
       - fvm::laplacian(turbulence->alphaEff(), he)
      ==
         rho*(U&g)
-      + radiation->Sh(thermo)
+      + radiation->Sh(thermo, he)
       + fvOptions(rho, he)
     );
 

applications/solvers/heatTransfer/buoyantSimpleFoam/EEqn.H

@@ -12,7 +12,7 @@
       - fvm::laplacian(turbulence->alphaEff(), he)
      ==
         rho*(U&g)
-      + radiation->Sh(thermo)
+      + radiation->Sh(thermo, he)
       + fvOptions(rho, he)
     );
 

applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/fluid/EEqn.H

@@ -12,7 +12,7 @@
       - fvm::laplacian(turb.alphaEff(), he)
      ==
         rho*(U&g)
-      + rad.Sh(thermo)
+      + rad.Sh(thermo, he)
       + fvOptions(rho, he)
     );
 

applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/fluid/readFluidMultiRegionSIMPLEControls.H

@@ -2,4 +2,3 @@
 
     const int nNonOrthCorr =
         simple.lookupOrDefault<int>("nNonOrthogonalCorrectors", 0);
-

applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/EEqn.H

@@ -18,7 +18,7 @@
       - fvm::laplacian(turb.alphaEff(), he)
      ==
         rho*(U&g)
-      + rad.Sh(thermo)
+      + rad.Sh(thermo, he)
       + fvOptions(rho, he)
     );
 

applications/solvers/heatTransfer/chtMultiRegionFoam/solid/solidRegionDiffNo.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -24,7 +24,9 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "solidRegionDiffNo.H"
-#include "fvc.H"
+#include "surfaceInterpolate.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 Foam::scalar Foam::solidRegionDiffNo
 (
@@ -34,21 +36,16 @@ Foam::scalar Foam::solidRegionDiffNo
     const volScalarField& kappa
 )
 {
-    scalar DiNum = 0.0;
-    scalar meanDiNum = 0.0;
-
-    //- Take care: can have fluid domains with 0 cells so do not test for
-    //  zero internal faces.
     surfaceScalarField kapparhoCpbyDelta
     (
-        mesh.surfaceInterpolation::deltaCoeffs()
+        sqr(mesh.surfaceInterpolation::deltaCoeffs())
-      * fvc::interpolate(kappa)
+       *fvc::interpolate(kappa)
-      / fvc::interpolate(Cprho)
+       /fvc::interpolate(Cprho)
     );
 
-    DiNum = max(kapparhoCpbyDelta).value()*runTime.deltaT().value();
+    const scalar DiNum = max(kapparhoCpbyDelta).value()*runTime.deltaTValue();
-
+    const scalar meanDiNum =
-    meanDiNum = (average(kapparhoCpbyDelta)).value()*runTime.deltaT().value();
+        average(kapparhoCpbyDelta).value()*runTime.deltaTValue();
 
     Info<< "Region: " << mesh.name() << " Diffusion Number mean: " << meanDiNum
         << " max: " << DiNum << endl;

applications/solvers/heatTransfer/thermoFoam/EEqn.H

@@ -17,7 +17,7 @@
         )
       - fvm::laplacian(alphaEff, he)
      ==
-        radiation->Sh(thermo)
+        radiation->Sh(thermo, he)
       + fvOptions(rho, he)
     );
 

applications/solvers/incompressible/adjointShapeOptimizationFoam/adjointShapeOptimizationFoam.C

@@ -93,8 +93,6 @@ int main(int argc, char *argv[])
     {
         Info<< "Time = " << runTime.timeName() << nl << endl;
 
-        laminarTransport.lookup("lambda") >> lambda;
-
         //alpha +=
         //    mesh.relaxationFactor("alpha")
         //   *(lambda*max(Ua & U, zeroSensitivity) - alpha);

applications/solvers/incompressible/adjointShapeOptimizationFoam/createFields.H

@@ -89,8 +89,19 @@ autoPtr<incompressible::turbulenceModel> turbulence
 dimensionedScalar zeroSensitivity("0", dimVelocity*dimVelocity, 0.0);
 dimensionedScalar zeroAlpha("0", dimless/dimTime, 0.0);
 
-dimensionedScalar lambda(laminarTransport.lookup("lambda"));
+dimensionedScalar lambda
-dimensionedScalar alphaMax(laminarTransport.lookup("alphaMax"));
+(
+    "lambda",
+    dimTime/sqr(dimLength),
+    laminarTransport
+);
+
+dimensionedScalar alphaMax
+(
+    "alphaMax",
+    dimless/dimTime,
+    laminarTransport
+);
 
 const labelList& inletCells = mesh.boundary()["inlet"].faceCells();
 //const labelList& outletCells = mesh.boundary()["outlet"].faceCells();

applications/solvers/incompressible/pimpleFoam/pimpleDyMFoam/pimpleDyMFoam.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -70,7 +70,6 @@ int main(int argc, char *argv[])
     {
         #include "readControls.H"
         #include "CourantNo.H"
-
         #include "setDeltaT.H"
 
         runTime++;

applications/solvers/lagrangian/DPMFoam/Allwclean

@@ -2,7 +2,10 @@
 cd ${0%/*} || exit 1    # Run from this directory
 
 wclean libso DPMTurbulenceModels
+
 wclean
 wclean MPPICFoam
+wclean DPMDyMFoam
+wclean DPMDyMFoam/MPPICDyMFoam
 
 #------------------------------------------------------------------------------

applications/solvers/lagrangian/DPMFoam/Allwmake

@@ -1,5 +1,5 @@
 #!/bin/sh
-cd ${0%/*} || exit 1
+cd ${0%/*} || exit 1    # Run from this directory
 
 # Parse arguments for library compilation
 . $WM_PROJECT_DIR/wmake/scripts/AllwmakeParseArguments
@@ -8,5 +8,7 @@ wmake $targetType DPMTurbulenceModels
 
 wmake $targetType
 wmake $targetType MPPICFoam
+wmake $targetType DPMDyMFoam
+wmake $targetType DPMDyMFoam/MPPICDyMFoam
 
 #------------------------------------------------------------------------------

applications/solvers/lagrangian/DPMFoam/DPMDyMFoam/DPMDyMFoam.C

@@ -0,0 +1,163 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2017 OpenFOAM Foundation
+     \\/     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 3 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, see <http://www.gnu.org/licenses/>.
+
+Application
+    DPMDyMFoam
+
+Description
+    Transient solver for the coupled transport of a single kinematic particle
+    cloud including the effect of the volume fraction of particles on the
+    continuous phase, with optional mesh motion and mesh topology changes.
+
+\*---------------------------------------------------------------------------*/
+
+#include "fvCFD.H"
+#include "dynamicFvMesh.H"
+#include "singlePhaseTransportModel.H"
+#include "PhaseIncompressibleTurbulenceModel.H"
+#include "pimpleControl.H"
+#include "CorrectPhi.H"
+
+#ifdef MPPIC
+    #include "basicKinematicMPPICCloud.H"
+    #define basicKinematicTypeCloud basicKinematicMPPICCloud
+#else
+    #include "basicKinematicCollidingCloud.H"
+    #define basicKinematicTypeCloud basicKinematicCollidingCloud
+#endif
+
+int main(int argc, char *argv[])
+{
+    argList::addOption
+    (
+        "cloudName",
+        "name",
+        "specify alternative cloud name. default is 'kinematicCloud'"
+    );
+
+    #include "postProcess.H"
+
+    #include "setRootCase.H"
+    #include "createTime.H"
+    #include "createDynamicFvMesh.H"
+    #include "createControls.H"
+    #include "createFields.H"
+    #include "createUcf.H"
+    #include "initContinuityErrs.H"
+
+    Info<< "\nStarting time loop\n" << endl;
+
+    while (runTime.run())
+    {
+        #include "readControls.H"
+        #include "CourantNo.H"
+        #include "setDeltaT.H"
+
+        runTime++;
+
+        Info<< "Time = " << runTime.timeName() << nl << endl;
+
+        mesh.update();
+
+        // Calculate absolute flux from the mapped surface velocity
+        phic = mesh.Sf() & Ucf;
+
+        if (mesh.changing() && correctPhi)
+        {
+            #include "correctPhic.H"
+        }
+
+        // Make the flux relative to the mesh motion
+        fvc::makeRelative(phic, Uc);
+
+        if (mesh.changing() && checkMeshCourantNo)
+        {
+            #include "meshCourantNo.H"
+        }
+
+        continuousPhaseTransport.correct();
+        muc = rhoc*continuousPhaseTransport.nu();
+
+        Info<< "Evolving " << kinematicCloud.name() << endl;
+        kinematicCloud.evolve();
+
+        // Update continuous phase volume fraction field
+        alphac = max(1.0 - kinematicCloud.theta(), alphacMin);
+        alphac.correctBoundaryConditions();
+        alphacf = fvc::interpolate(alphac);
+        alphaPhic = alphacf*phic;
+
+        fvVectorMatrix cloudSU(kinematicCloud.SU(Uc));
+        volVectorField cloudVolSUSu
+        (
+            IOobject
+            (
+                "cloudVolSUSu",
+                runTime.timeName(),
+                mesh
+            ),
+            mesh,
+            dimensionedVector
+            (
+                "0",
+                cloudSU.dimensions()/dimVolume,
+                Zero
+            ),
+            zeroGradientFvPatchVectorField::typeName
+        );
+
+        cloudVolSUSu.primitiveFieldRef() = -cloudSU.source()/mesh.V();
+        cloudVolSUSu.correctBoundaryConditions();
+        cloudSU.source() = Zero;
+
+        // --- Pressure-velocity PIMPLE corrector loop
+        while (pimple.loop())
+        {
+            #include "UcEqn.H"
+
+            // --- PISO loop
+            while (pimple.correct())
+            {
+                #include "pEqn.H"
+            }
+
+            if (pimple.turbCorr())
+            {
+                continuousPhaseTurbulence->correct();
+            }
+        }
+
+        runTime.write();
+
+        Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
+            << "  ClockTime = " << runTime.elapsedClockTime() << " s"
+            << nl << endl;
+    }
+
+    Info<< "End\n" << endl;
+
+    return 0;
+}
+
+
+// ************************************************************************* //

applications/solvers/lagrangian/DPMFoam/DPMDyMFoam/MPPICDyMFoam/MPPICDyMFoam.C

@@ -0,0 +1,41 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2017 OpenFOAM Foundation
+     \\/     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 3 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, see <http://www.gnu.org/licenses/>.
+
+Application
+    MPPICDyMFoam
+
+Description
+    Transient solver for the coupled transport of a single kinematic particle
+    cloud including the effect of the volume fraction of particles on the
+    continuous phase. Multi-Phase Particle In Cell (MPPIC) modeling is used to
+    represent collisions without resolving particle-particle interactions,
+    with optional mesh motion and mesh topology changes.
+
+\*---------------------------------------------------------------------------*/
+
+#define MPPIC
+
+#include "DPMDyMFoam.C"
+
+
+// ************************************************************************* //

applications/solvers/lagrangian/DPMFoam/DPMDyMFoam/MPPICDyMFoam/Make/files

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

applications/solvers/lagrangian/DPMFoam/DPMDyMFoam/MPPICDyMFoam/Make/options

@@ -0,0 +1,42 @@
+EXE_INC = \
+    -I.. \
+    -I../.. \
+    -I../DPMTurbulenceModels/lnInclude \
+    -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/meshTools/lnInclude \
+    -I$(LIB_SRC)/lagrangian/basic/lnInclude \
+    -I$(LIB_SRC)/lagrangian/intermediate/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/specie/lnInclude \
+    -I$(LIB_SRC)/transportModels/compressible/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/reactionThermo/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/radiation/lnInclude \
+    -I$(LIB_SRC)/transportModels \
+    -I$(LIB_SRC)/transportModels/incompressible/singlePhaseTransportModel \
+    -I$(LIB_SRC)/TurbulenceModels/turbulenceModels/lnInclude \
+    -I$(LIB_SRC)/TurbulenceModels/incompressible/lnInclude \
+    -I$(LIB_SRC)/TurbulenceModels/phaseIncompressible/lnInclude \
+    -I$(LIB_SRC)/regionModels/regionModel/lnInclude \
+    -I$(LIB_SRC)/regionModels/surfaceFilmModels/lnInclude \
+    -I$(LIB_SRC)/dynamicFvMesh/lnInclude \
+    -I$(LIB_SRC)/dynamicMesh/lnInclude
+
+EXE_LIBS = \
+    -lfiniteVolume \
+    -lfvOptions \
+    -lmeshTools \
+    -llagrangian \
+    -llagrangianIntermediate \
+    -llagrangianTurbulence \
+    -lspecie \
+    -lradiationModels \
+    -lincompressibleTransportModels \
+    -lturbulenceModels \
+    -lincompressibleTurbulenceModels \
+    -lDPMTurbulenceModels \
+    -lregionModels \
+    -lsurfaceFilmModels \
+    -lsampling \
+    -ldynamicFvMesh \
+    -ltopoChangerFvMesh \
+    -ldynamicMesh

applications/solvers/lagrangian/DPMFoam/DPMDyMFoam/Make/files

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

applications/solvers/lagrangian/DPMFoam/DPMDyMFoam/Make/options

@@ -0,0 +1,41 @@
+EXE_INC = \
+    -I.. \
+    -I../DPMTurbulenceModels/lnInclude \
+    -I$(LIB_SRC)/lagrangian/basic/lnInclude \
+    -I$(LIB_SRC)/lagrangian/intermediate/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/specie/lnInclude \
+    -I$(LIB_SRC)/transportModels/compressible/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/reactionThermo/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/radiation/lnInclude \
+    -I$(LIB_SRC)/transportModels \
+    -I$(LIB_SRC)/transportModels/incompressible/singlePhaseTransportModel \
+    -I$(LIB_SRC)/TurbulenceModels/turbulenceModels/lnInclude \
+    -I$(LIB_SRC)/TurbulenceModels/incompressible/lnInclude \
+    -I$(LIB_SRC)/TurbulenceModels/phaseIncompressible/lnInclude \
+    -I$(LIB_SRC)/regionModels/regionModel/lnInclude \
+    -I$(LIB_SRC)/regionModels/surfaceFilmModels/lnInclude \
+    -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/meshTools/lnInclude \
+    -I$(LIB_SRC)/dynamicFvMesh/lnInclude \
+    -I$(LIB_SRC)/dynamicMesh/lnInclude
+
+EXE_LIBS = \
+    -llagrangian \
+    -llagrangianIntermediate \
+    -llagrangianTurbulence \
+    -lspecie \
+    -lradiationModels \
+    -lincompressibleTransportModels \
+    -lturbulenceModels \
+    -lincompressibleTurbulenceModels \
+    -lDPMTurbulenceModels \
+    -lregionModels \
+    -lsurfaceFilmModels \
+    -lsampling \
+    -lfiniteVolume \
+    -lfvOptions \
+    -lmeshTools \
+    -ldynamicFvMesh \
+    -ltopoChangerFvMesh \
+    -ldynamicMesh

applications/solvers/lagrangian/DPMFoam/DPMDyMFoam/correctPhic.H

@@ -0,0 +1,11 @@
+CorrectPhi
+(
+    Uc,
+    phic,
+    p,
+    dimensionedScalar("rAUf", dimTime, 1),
+    geometricZeroField(),
+    pimple
+);
+
+#include "continuityErrs.H"

applications/solvers/lagrangian/DPMFoam/DPMDyMFoam/createControls.H

@@ -0,0 +1,12 @@
+#include "createControl.H"
+#include "createTimeControls.H"
+
+bool correctPhi
+(
+    pimple.dict().lookupOrDefault("correctPhi", false)
+);
+
+bool checkMeshCourantNo
+(
+    pimple.dict().lookupOrDefault("checkMeshCourantNo", false)
+);

applications/solvers/lagrangian/DPMFoam/DPMDyMFoam/createUcf.H

@@ -2,8 +2,8 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2013 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2017 OpenFOAM Foundation
-     \\/     M anipulation  | Copyright (C) 2016 OpenCFD Ltd.
+     \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -21,34 +21,29 @@ License
     You should have received a copy of the GNU General Public License
     along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
 
-\*---------------------------------------------------------------------------*/
+Global
+    createUcf
 
-#include "IMULES.H"
+Description
-#include "profiling.H"
+    Creates and initialises the velocity velocity field Ucf.
+
+\*---------------------------------------------------------------------------*/
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-void Foam::MULES::implicitSolve
+Info<< "Reading/calculating continuous phase face velocity Ucf\n" << endl;
+
+surfaceVectorField Ucf
 (
-    volScalarField& psi,
+    IOobject
-    const surfaceScalarField& phi,
-    surfaceScalarField& phiPsi,
-    const scalar psiMax,
-    const scalar psiMin
-)
-{
-    addProfiling(solve, "MULES::implicitSolve");
-
-    implicitSolve
     (
-        geometricOneField(),
+        "Ucf",
-        psi,
+        runTime.timeName(),
-        phi,
+        mesh,
-        phiPsi,
+        IOobject::READ_IF_PRESENT,
-        zeroField(), zeroField(),
+        IOobject::AUTO_WRITE
-        psiMax, psiMin
+    ),
-    );
+    fvc::interpolate(Uc)
-}
+);
-
 
 // ************************************************************************* //

applications/solvers/lagrangian/DPMFoam/DPMDyMFoam/pEqn.H

@@ -0,0 +1,62 @@
+{
+    volVectorField HbyA(constrainHbyA(rAUc*UcEqn.H(), Uc, p));
+
+    surfaceScalarField phiHbyA
+    (
+        "phiHbyA",
+        (
+           fvc::flux(HbyA)
+         + alphacf*rAUcf*fvc::ddtCorr(Uc, Ucf)
+        )
+    );
+
+    if (p.needReference())
+    {
+        fvc::makeRelative(phiHbyA, Uc);
+        adjustPhi(phiHbyA, Uc, p);
+        fvc::makeAbsolute(phiHbyA, Uc);
+    }
+
+    phiHbyA += phicForces;
+
+    // Update the pressure BCs to ensure flux consistency
+    constrainPressure(p, Uc, phiHbyA, rAUcf);
+
+    // Non-orthogonal pressure corrector loop
+    while (pimple.correctNonOrthogonal())
+    {
+        fvScalarMatrix pEqn
+        (
+            fvm::laplacian(alphacf*rAUcf, p)
+         ==
+            fvc::ddt(alphac) + fvc::div(alphacf*phiHbyA)
+        );
+
+        pEqn.setReference(pRefCell, pRefValue);
+
+        pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter())));
+
+        if (pimple.finalNonOrthogonalIter())
+        {
+            phic = phiHbyA - pEqn.flux()/alphacf;
+
+            p.relax();
+
+            Uc = HbyA
+               + rAUc
+                *fvc::reconstruct((phicForces - pEqn.flux()/alphacf)/rAUcf);
+            Uc.correctBoundaryConditions();
+
+            {
+                Ucf = fvc::interpolate(Uc);
+                surfaceVectorField n(mesh.Sf()/mesh.magSf());
+                Ucf += n*(phic/mesh.magSf() - (n & Ucf));
+            }
+
+            // Make the fluxes relative to the mesh motion
+            fvc::makeRelative(phic, Uc);
+        }
+    }
+}
+
+#include "continuityErrs.H"

applications/solvers/lagrangian/DPMFoam/DPMDyMFoam/readControls.H

@@ -0,0 +1,5 @@
+#include "readTimeControls.H"
+
+correctPhi = pimple.dict().lookupOrDefault("correctPhi", false);
+
+checkMeshCourantNo = pimple.dict().lookupOrDefault("checkMeshCourantNo", false);

applications/solvers/lagrangian/DPMFoam/MPPICFoam/Make/options

@@ -25,7 +25,6 @@ EXE_LIBS = \
     -llagrangian \
     -llagrangianIntermediate \
     -llagrangianTurbulence \
-    -lthermophysicalFunctions \
     -lspecie \
     -lradiationModels \
     -lincompressibleTransportModels \

applications/solvers/lagrangian/DPMFoam/Make/options

@@ -21,7 +21,6 @@ EXE_LIBS = \
     -llagrangian \
     -llagrangianIntermediate \
     -llagrangianTurbulence \
-    -lthermophysicalFunctions \
     -lspecie \
     -lradiationModels \
     -lincompressibleTransportModels \

applications/solvers/lagrangian/DPMFoam/pEqn.H

@@ -1,6 +1,5 @@
 {
-    volVectorField HbyA("HbyA", Uc);
+    volVectorField HbyA(constrainHbyA(rAUc*UcEqn.H(), Uc, p));
-    HbyA = rAUc*UcEqn.H();
 
     surfaceScalarField phiHbyA
     (
@@ -8,10 +7,16 @@
         (
            fvc::flux(HbyA)
          + alphacf*rAUcf*fvc::ddtCorr(Uc, phic)
-         + phicForces
         )
     );
 
+    if (p.needReference())
+    {
+        adjustPhi(phiHbyA, Uc, p);
+    }
+
+    phiHbyA += phicForces;
+
     // Update the pressure BCs to ensure flux consistency
     constrainPressure(p, Uc, phiHbyA, rAUcf);
 

applications/solvers/lagrangian/coalChemistryFoam/EEqn.H

@@ -18,10 +18,10 @@
       - fvm::laplacian(turbulence->alphaEff(), he)
      ==
         rho*(U&g)
-      + combustion->Sh()
+      + Qdot
       + coalParcels.Sh(he)
       + limestoneParcels.Sh(he)
-      + radiation->Sh(thermo)
+      + radiation->Sh(thermo, he)
       + fvOptions(rho, he)
     );
 

applications/solvers/lagrangian/coalChemistryFoam/Make/options

@@ -10,10 +10,7 @@ EXE_INC = \
     -I$(LIB_SRC)/thermophysicalModels/specie/lnInclude \
     -I$(LIB_SRC)/transportModels/compressible/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
-    -I$(LIB_SRC)/thermophysicalModels/properties/liquidProperties/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/thermophysicalProperties/lnInclude \
-    -I$(LIB_SRC)/thermophysicalModels/properties/liquidMixtureProperties/lnInclude \
-    -I$(LIB_SRC)/thermophysicalModels/properties/solidProperties/lnInclude \
-    -I$(LIB_SRC)/thermophysicalModels/properties/solidMixtureProperties/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/thermophysicalFunctions/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/reactionThermo/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/SLGThermo/lnInclude \
@@ -38,11 +35,7 @@ EXE_LIBS = \
     -lspecie \
     -lcompressibleTransportModels \
     -lfluidThermophysicalModels \
-    -lliquidProperties \
+    -lthermophysicalProperties \
-    -lliquidMixtureProperties \
-    -lsolidProperties \
-    -lsolidMixtureProperties \
-    -lthermophysicalFunctions \
     -lreactionThermophysicalModels \
     -lSLGThermo \
     -lchemistryModel \

applications/solvers/lagrangian/coalChemistryFoam/YEqn.H

@@ -12,7 +12,7 @@ tmp<fv::convectionScheme<scalar>> mvConvection
 
 {
     combustion->correct();
-    dQ = combustion->dQ();
+    Qdot = combustion->Qdot();
     volScalarField Yt(0.0*Y[0]);
 
     forAll(Y, i)

applications/solvers/lagrangian/coalChemistryFoam/createFields.H

@@ -131,18 +131,18 @@ volScalarField dpdt
 Info<< "Creating field kinetic energy K\n" << endl;
 volScalarField K("K", 0.5*magSqr(U));
 
-volScalarField dQ
+volScalarField Qdot
 (
     IOobject
     (
-        "dQ",
+        "Qdot",
         runTime.timeName(),
         mesh,
-        IOobject::NO_READ,
+        IOobject::READ_IF_PRESENT,
         IOobject::AUTO_WRITE
     ),
     mesh,
-    dimensionedScalar("dQ", dimEnergy/dimTime, 0.0)
+    dimensionedScalar("Qdot", dimEnergy/dimVolume/dimTime, 0.0)
 );
 
 #include "createMRF.H"

applications/solvers/lagrangian/coalChemistryFoam/setRDeltaT.H

@@ -80,7 +80,7 @@ License
             (
                (coalParcels.hsTrans() + limestoneParcels.hsTrans())
               /(mesh.V()*runTime.deltaT())
-             + combustion->Sh()()
+             + Qdot
             )
            /(
                alphaTemp

applications/solvers/lagrangian/icoUncoupledKinematicParcelFoam/Make/options

@@ -19,7 +19,6 @@ EXE_LIBS = \
     -llagrangian \
     -llagrangianIntermediate \
     -llagrangianTurbulence \
-    -lthermophysicalFunctions \
     -lcompressibleTransportModels \
     -lfluidThermophysicalModels \
     -lspecie \

applications/solvers/lagrangian/icoUncoupledKinematicParcelFoam/icoUncoupledKinematicParcelDyMFoam/Make/options

@@ -23,7 +23,6 @@ EXE_LIBS = \
     -llagrangian \
     -llagrangianIntermediate \
     -llagrangianTurbulence \
-    -lthermophysicalFunctions \
     -lcompressibleTransportModels \
     -lfluidThermophysicalModels \
     -lspecie \

applications/solvers/lagrangian/reactingParcelFilmFoam/EEqn.H

@@ -20,8 +20,8 @@
         rho*(U&g)
       + parcels.Sh(he)
       + surfaceFilm.Sh()
-      + radiation->Sh(thermo)
+      + radiation->Sh(thermo, he)
-      + combustion->Sh()
+      + Qdot
       + fvOptions(rho, he)
     );
 

applications/solvers/lagrangian/reactingParcelFilmFoam/Make/options

@@ -9,10 +9,7 @@ EXE_INC = \
     -I$(LIB_SRC)/thermophysicalModels/specie/lnInclude \
     -I$(LIB_SRC)/transportModels/compressible/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
-    -I$(LIB_SRC)/thermophysicalModels/properties/solidProperties/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/thermophysicalProperties/lnInclude \
-    -I$(LIB_SRC)/thermophysicalModels/properties/solidMixtureProperties/lnInclude \
-    -I$(LIB_SRC)/thermophysicalModels/properties/liquidProperties/lnInclude \
-    -I$(LIB_SRC)/thermophysicalModels/properties/liquidMixtureProperties/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/thermophysicalFunctions/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/reactionThermo/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/SLGThermo/lnInclude \
@@ -36,9 +33,6 @@ EXE_LIBS = \
     -lspecie \
     -lcompressibleTransportModels \
     -lfluidThermophysicalModels \
-    -lsolidProperties \
-    -lsolidMixtureProperties \
-    -lthermophysicalFunctions \
     -lreactionThermophysicalModels \
     -lSLGThermo \
     -lchemistryModel \

applications/solvers/lagrangian/reactingParcelFilmFoam/YEqn.H

@@ -12,7 +12,7 @@ tmp<fv::convectionScheme<scalar>> mvConvection
 
 {
     combustion->correct();
-    dQ = combustion->dQ();
+    Qdot = combustion->Qdot();
     volScalarField Yt(0.0*Y[0]);
 
     forAll(Y, i)

applications/solvers/lagrangian/reactingParcelFilmFoam/createFields.H

@@ -134,18 +134,18 @@ Switch solvePrimaryRegion
     additionalControlsDict.lookup("solvePrimaryRegion")
 );
 
-volScalarField dQ
+volScalarField Qdot
 (
     IOobject
     (
-        "dQ",
+        "Qdot",
         runTime.timeName(),
         mesh,
-        IOobject::NO_READ,
+        IOobject::READ_IF_PRESENT,
         IOobject::AUTO_WRITE
     ),
     mesh,
-    dimensionedScalar("dQ", dimEnergy/dimTime, 0.0)
+    dimensionedScalar("Qdot", dimEnergy/dimVolume/dimTime, 0.0)
 );
 
 #include "createMRF.H"

applications/solvers/lagrangian/reactingParcelFoam/EEqn.H

@@ -19,8 +19,8 @@
      ==
         rho*(U&g)
       + parcels.Sh(he)
-      + radiation->Sh(thermo)
+      + radiation->Sh(thermo, he)
-      + combustion->Sh()
+      + Qdot
       + fvOptions(rho, he)
     );
 
@@ -35,6 +35,6 @@
     thermo.correct();
     radiation->correct();
 
-    Info<< "T gas min/max   = " << min(T).value() << ", "
+    Info<< "T gas min/max   " << min(T).value() << ", "
         << max(T).value() << endl;
 }

applications/solvers/lagrangian/reactingParcelFoam/Make/options

@@ -11,10 +11,7 @@ EXE_INC = \
     -I$(LIB_SRC)/thermophysicalModels/specie/lnInclude \
     -I$(LIB_SRC)/transportModels/compressible/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
-    -I$(LIB_SRC)/thermophysicalModels/properties/liquidProperties/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/thermophysicalProperties/lnInclude \
-    -I$(LIB_SRC)/thermophysicalModels/properties/liquidMixtureProperties/lnInclude \
-    -I$(LIB_SRC)/thermophysicalModels/properties/solidProperties/lnInclude \
-    -I$(LIB_SRC)/thermophysicalModels/properties/solidMixtureProperties/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/thermophysicalFunctions/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/reactionThermo/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/SLGThermo/lnInclude \
@@ -39,11 +36,7 @@ EXE_LIBS = \
     -lspecie \
     -lcompressibleTransportModels \
     -lfluidThermophysicalModels \
-    -lliquidProperties \
+    -lthermophysicalProperties \
-    -lliquidMixtureProperties \
-    -lsolidProperties \
-    -lsolidMixtureProperties \
-    -lthermophysicalFunctions \
     -lreactionThermophysicalModels \
     -lSLGThermo \
     -lchemistryModel \

applications/solvers/lagrangian/reactingParcelFoam/YEqn.H

@@ -11,7 +11,7 @@ tmp<fv::convectionScheme<scalar>> mvConvection
 
 {
     combustion->correct();
-    dQ = combustion->dQ();
+    Qdot = combustion->Qdot();
     volScalarField Yt(0.0*Y[0]);
 
     forAll(Y, i)

applications/solvers/lagrangian/reactingParcelFoam/createFields.H

@@ -121,18 +121,18 @@ forAll(Y, i)
 }
 fields.add(thermo.he());
 
-volScalarField dQ
+volScalarField Qdot
 (
     IOobject
     (
-        "dQ",
+        "Qdot",
         runTime.timeName(),
         mesh,
-        IOobject::NO_READ,
+        IOobject::READ_IF_PRESENT,
         IOobject::AUTO_WRITE
     ),
     mesh,
-    dimensionedScalar("dQ", dimEnergy/dimTime, 0.0)
+    dimensionedScalar("Qdot", dimEnergy/dimVolume/dimTime, 0.0)
 );
 
 #include "createMRF.H"

applications/solvers/lagrangian/reactingParcelFoam/setRDeltaT.H

@@ -79,7 +79,7 @@ License
             mag
             (
                parcels.hsTrans()/(mesh.V()*runTime.deltaT())
-             + combustion->Sh()()
+             + Qdot
             )
            /(
                alphaTemp

applications/solvers/lagrangian/reactingParcelFoam/simpleReactingParcelFoam/EEqn.H

@@ -13,8 +13,8 @@
      ==
         rho*(U&g)
       + parcels.Sh(he)
-      + radiation->Sh(thermo)
+      + radiation->Sh(thermo, he)
-      + combustion->Sh()
+      + Qdot
       + fvOptions(rho, he)
     );
 

applications/solvers/lagrangian/reactingParcelFoam/simpleReactingParcelFoam/Make/options

@@ -11,10 +11,7 @@ EXE_INC = \
     -I$(LIB_SRC)/thermophysicalModels/specie/lnInclude \
     -I$(LIB_SRC)/transportModels/compressible/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
-    -I$(LIB_SRC)/thermophysicalModels/properties/liquidProperties/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/thermophysicalProperties/lnInclude \
-    -I$(LIB_SRC)/thermophysicalModels/properties/liquidMixtureProperties/lnInclude \
-    -I$(LIB_SRC)/thermophysicalModels/properties/solidProperties/lnInclude \
-    -I$(LIB_SRC)/thermophysicalModels/properties/solidMixtureProperties/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/thermophysicalFunctions/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/reactionThermo/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/SLGThermo/lnInclude \
@@ -39,11 +36,7 @@ EXE_LIBS = \
     -lspecie \
     -lcompressibleTransportModels \
     -lfluidThermophysicalModels \
-    -lliquidProperties \
+    -lthermophysicalProperties \
-    -lliquidMixtureProperties \
-    -lsolidProperties \
-    -lsolidMixtureProperties \
-    -lthermophysicalFunctions \
     -lreactionThermophysicalModels \
     -lSLGThermo \
     -lchemistryModel \

applications/solvers/lagrangian/reactingParcelFoam/simpleReactingParcelFoam/YEqn.H

@@ -11,7 +11,7 @@ tmp<fv::convectionScheme<scalar>> mvConvection
 
 {
     combustion->correct();
-    dQ = combustion->dQ();
+    Qdot = combustion->Qdot();
     volScalarField Yt(0.0*Y[0]);
 
     forAll(Y, i)

applications/solvers/lagrangian/reactingParcelFoam/simpleReactingParcelFoam/createFields.H

@@ -103,18 +103,18 @@ forAll(Y, i)
 }
 fields.add(thermo.he());
 
-volScalarField dQ
+volScalarField Qdot
 (
     IOobject
     (
-        "dQ",
+        "Qdot",
         runTime.timeName(),
         mesh,
-        IOobject::NO_READ,
+        IOobject::READ_IF_PRESENT,
         IOobject::AUTO_WRITE
     ),
     mesh,
-    dimensionedScalar("dQ", dimEnergy/dimTime, 0.0)
+    dimensionedScalar("Qdot", dimEnergy/dimVolume/dimTime, 0.0)
 );
 
 #include "createMRF.H"

applications/solvers/lagrangian/simpleCoalParcelFoam/EEqn.H

@@ -12,9 +12,9 @@
       - fvm::laplacian(turbulence->alphaEff(), he)
      ==
         rho*(U&g)
+      + Qdot
       + parcels.Sh(he)
-      + radiation->Sh(thermo)
+      + radiation->Sh(thermo, he)
-      + combustion->Sh()
       + fvOptions(rho, he)
     );
 
@@ -25,6 +25,7 @@
     EEqn.solve();
 
     fvOptions.correct(he);
+
     thermo.correct();
     radiation->correct();
 

applications/solvers/lagrangian/simpleCoalParcelFoam/Make/options

@@ -10,10 +10,7 @@ EXE_INC = \
     -I$(LIB_SRC)/thermophysicalModels/specie/lnInclude \
     -I$(LIB_SRC)/transportModels/compressible/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
-    -I$(LIB_SRC)/thermophysicalModels/properties/liquidProperties/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/thermophysicalProperties/lnInclude \
-    -I$(LIB_SRC)/thermophysicalModels/properties/liquidMixtureProperties/lnInclude \
-    -I$(LIB_SRC)/thermophysicalModels/properties/solidProperties/lnInclude \
-    -I$(LIB_SRC)/thermophysicalModels/properties/solidMixtureProperties/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/thermophysicalFunctions/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/reactionThermo/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/SLGThermo/lnInclude \
@@ -37,11 +34,6 @@ EXE_LIBS = \
     -llagrangianTurbulence \
     -lspecie \
     -lfluidThermophysicalModels \
-    -lliquidProperties \
-    -lliquidMixtureProperties \
-    -lsolidProperties \
-    -lsolidMixtureProperties \
-    -lthermophysicalFunctions \
     -lreactionThermophysicalModels \
     -lSLGThermo \
     -lchemistryModel \

applications/solvers/lagrangian/simpleCoalParcelFoam/YEqn.H

@@ -11,7 +11,7 @@ tmp<fv::convectionScheme<scalar>> mvConvection
 
 {
     combustion->correct();
-    dQ = combustion->dQ();
+    Qdot = combustion->Qdot();
     volScalarField Yt(0.0*Y[0]);
 
     forAll(Y, i)

applications/solvers/lagrangian/simpleCoalParcelFoam/createFields.H

@@ -103,18 +103,18 @@ forAll(Y, i)
 }
 fields.add(thermo.he());
 
-volScalarField dQ
+volScalarField Qdot
 (
     IOobject
     (
-        "dQ",
+        "Qdot",
         runTime.timeName(),
         mesh,
         IOobject::NO_READ,
         IOobject::AUTO_WRITE
     ),
     mesh,
-    dimensionedScalar("dQ", dimEnergy/dimTime, 0.0)
+    dimensionedScalar("Qdot", dimEnergy/dimTime, 0.0)
 );
 
 #include "createMRF.H"

applications/solvers/lagrangian/sprayFoam/Make/options

@@ -13,10 +13,7 @@ EXE_INC = \
     -I$(LIB_SRC)/thermophysicalModels/specie/lnInclude \
     -I$(LIB_SRC)/transportModels/compressible/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
-    -I$(LIB_SRC)/thermophysicalModels/properties/liquidProperties/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/thermophysicalProperties/lnInclude \
-    -I$(LIB_SRC)/thermophysicalModels/properties/liquidMixtureProperties/lnInclude \
-    -I$(LIB_SRC)/thermophysicalModels/properties/solidProperties/lnInclude \
-    -I$(LIB_SRC)/thermophysicalModels/properties/solidMixtureProperties/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/thermophysicalFunctions/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/reactionThermo/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/SLGThermo/lnInclude \
@@ -37,11 +34,7 @@ EXE_LIBS = \
     -lspecie \
     -lcompressibleTransportModels \
     -lfluidThermophysicalModels \
-    -lliquidProperties \
+    -lthermophysicalProperties \
-    -lliquidMixtureProperties \
-    -lsolidProperties \
-    -lsolidMixtureProperties \
-    -lthermophysicalFunctions \
     -lreactionThermophysicalModels \
     -lSLGThermo \
     -lchemistryModel \

applications/solvers/lagrangian/sprayFoam/createFields.H

@@ -117,18 +117,18 @@ forAll(Y, i)
 }
 fields.add(thermo.he());
 
-volScalarField dQ
+volScalarField Qdot
 (
     IOobject
     (
-        "dQ",
+        "Qdot",
         runTime.timeName(),
         mesh,
         IOobject::NO_READ,
         IOobject::AUTO_WRITE
     ),
     mesh,
-    dimensionedScalar("dQ", dimEnergy/dimTime, 0.0)
+    dimensionedScalar("Qdot", dimEnergy/dimVolume/dimTime, 0.0)
 );
 
 #include "createMRF.H"

applications/solvers/lagrangian/sprayFoam/sprayDyMFoam/Make/options

@@ -14,10 +14,7 @@ EXE_INC = \
     -I$(LIB_SRC)/thermophysicalModels/specie/lnInclude \
     -I$(LIB_SRC)/transportModels/compressible/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
-    -I$(LIB_SRC)/thermophysicalModels/properties/liquidProperties/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/thermophysicalProperties/lnInclude \
-    -I$(LIB_SRC)/thermophysicalModels/properties/liquidMixtureProperties/lnInclude \
-    -I$(LIB_SRC)/thermophysicalModels/properties/solidProperties/lnInclude \
-    -I$(LIB_SRC)/thermophysicalModels/properties/solidMixtureProperties/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/thermophysicalFunctions/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/reactionThermo/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/SLGThermo/lnInclude \
@@ -41,11 +38,7 @@ EXE_LIBS = \
     -lspecie \
     -lcompressibleTransportModels \
     -lfluidThermophysicalModels \
-    -lliquidProperties \
+    -lthermophysicalProperties \
-    -lliquidMixtureProperties \
-    -lsolidProperties \
-    -lsolidMixtureProperties \
-    -lthermophysicalFunctions \
     -lreactionThermophysicalModels \
     -lSLGThermo \
     -lchemistryModel \