BACKPORT: spack configuration patches (as per OpenFOAM-v1706)

etc/bashrc
  - improved robustness when sourcing
  - source top-level prefs.sh first (for sysadmin changes)

etc/config.*/settings
  - write job control information to the user directory
  - site/ directory under the OpenFOAM project dir, not its parent dir

etc/config.*/mpi
  - added USERMPI as place for spack mpi information

mgridgen, zoltan:
  - make location configurable

metis, scotch:
  - also check lib path (not just lib64)

.gitignore

@@ -39,12 +39,8 @@ solaris*Gcc*/
 SunOS*Gcc*/
 platforms/
 
-# Top-level build directories
-/build/
-/platforms/
-
-# Reinstate wmake rules that might look like build directories
-!/wmake/rules/*/
+# Reinstate wmake/rules that might look like build directories
+!wmake/rules/*/
 
 # doxygen generated documentation
 doc/Doxygen/html
@@ -62,12 +58,10 @@ 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
 *.tar.gz
-*.tar.xz
 *.tar
 *.tgz
 *.gtgz

Allwmake

@@ -1,29 +1,24 @@
 #!/bin/sh
-# Run from OPENFOAM top-level directory only
-cd ${0%/*} && wmakeCheckPwd "$WM_PROJECT_DIR" 2>/dev/null || {
-    echo "Error (${0##*/}) : not located in \$WM_PROJECT_DIR"
-    echo "    Check your OpenFOAM environment and installation"
-    exit 1
-}
-[ -d "$WM_PROJECT_DIR" -a -f "$WM_PROJECT_DIR/etc/bashrc" ] || {
-    echo "Error (${0##*/}) : WM_PROJECT_DIR appears to be incorrect"
-    echo "    Check your OpenFOAM environment and installation"
-    exit 1
-}
+cd ${0%/*} || exit 1    # Run from this directory
 
 # Parse arguments for library compilation
 . $WM_PROJECT_DIR/wmake/scripts/AllwmakeParseArguments
 
-#------------------------------------------------------------------------------
-echo "========================================"
-date "+%Y-%m-%d %H:%M:%S %z" 2>/dev/null || echo "date is unknown"
-echo "Starting ${WM_PROJECT_DIR##*/} ${0##*}"
-echo "  $WM_COMPILER $WM_COMPILER_TYPE compiler"
-echo "  ${WM_OPTIONS}, with ${WM_MPLIB} ${FOAM_MPI}"
-echo
+wmakeCheckPwd "$WM_PROJECT_DIR" || {
+    echo "Allwmake error: Current directory is not \$WM_PROJECT_DIR"
+    echo "    The environment variables are inconsistent with the installation."
+    echo "    Check the OpenFOAM entries in your dot-files and source them."
+    exit 1
+}
 
-# Compile wmake tools
-(cd "${WM_DIR:-wmake}/src" && make)
+[ -n "$FOAM_EXT_LIBBIN" ] || {
+    echo "Allwmake error: FOAM_EXT_LIBBIN not set"
+    echo "    Check the OpenFOAM entries in your dot-files and source them."
+    exit 1
+}
+
+# Compile wmake support applications
+(cd wmake/src && make)
 
 # Compile ThirdParty libraries and applications
 if [ -d "$WM_THIRD_PARTY_DIR" ]

BuildIssues.txt

@@ -1,4 +1,4 @@
-OpenFOAM-1706
+OpenFOAM-v1612+
 ==================
 Known Build Issues
 ==================
@@ -22,15 +22,16 @@ VTK
 ---
 
 If using the runTimePostProcessing to create on-the-fly images, you
-can simply just compile ParaView and these libraries will be used.
+can either simply just compile ParaView-5.0.1 and these libraries will
+be used.
 
 If you elect to use a separate VTK compilation (for example for
 off-screen rendering), it is advisable to reuse the VTK libraries that
-are provided with ParaView by making an appropriate symlink
+are provided with ParaView-5.0.1, by making an appropriate symlink
 prior to using makeVTK. This doesn't just reduce disk-space, but works
-much better than using the VTK tar file.
+much better than using the VTK-7.1.0.tar file.
 
-Using runTimePostProcessing with the 'plain' VTK libraries does
+Using runTimePostProcessing with the 'plain' VTK-7.1.0 libraries does
 generally work, but does not exit cleanly:
 
     symbol lookup error: .../linux64Gcc/VTK-7.1.0/lib/libvtkCommonExecutionModel-7.1.so.1:
@@ -59,21 +60,4 @@ there are also minimum gcc/g++ requirements there:
 If your system compiler is too old to build the minimum required gcc or
 clang/llvm, it is just simply too old.
 
-
--------------------------
-Building with spack
--------------------------
-
-If you are building with spack, note that the depends_on for paraview
-resolves poorly. The +qt dependency (for building the reader module)
-may need to be specified as a preference by including the following in
-your `~/.spack/packages.yaml` file:
-
-    packages:
-        paraview:
-            variants: +qt
-
-It appears that spack will otherwise ignore any paraview+qt version
-and attempt to install a paraview~qt version instead.
-
 --

META-INFO/.gitignore

@@ -0,0 +1,9 @@
+# Do not track build information
+build-info
+
+# Do not track time-stamp
+time-stamp
+
+# Do not track any manifest files
+Manifest.txt
+manifest.txt

META-INFO/README.md

@@ -0,0 +1,5 @@
+# META-INFO
+
+Meta-information is for OpenFOAM internal use only.
+
+(Backport from 1812)

META-INFO/api-info

@@ -0,0 +1,2 @@
+api=1612
+patch=180618

README.md

@@ -13,7 +13,7 @@ OpenCFD Ltd grants use of the OpenFOAM trademark by Third Parties on a licence b
 
 Please [contact OpenCFD](http://www.openfoam.com/contact) if you have any questions on the use of the OpenFOAM trademark.
 
-Violations of the Trademark are continuously monitored, and will be duly prosecuted.
+Violations of the Trademark are continuously monitored, and will be duly prosecuted. 
 
 
 # Useful Links

applications/Allwmake

@@ -1,20 +1,21 @@
 #!/bin/sh
-# Run from OPENFOAM applications/ directory only
-cd ${0%/*} && wmakeCheckPwd "$WM_PROJECT_DIR/applications" 2>/dev/null || {
-    echo "Error (${0##*/}) : not located in \$WM_PROJECT_DIR/applications"
-    echo "    Check your OpenFOAM environment and installation"
-    exit 1
-}
-[ -d "$WM_PROJECT_DIR" -a -f "$WM_PROJECT_DIR/etc/bashrc" ] || {
-    echo "Error (${0##*/}) : WM_PROJECT_DIR appears to be incorrect"
-    echo "    Check your OpenFOAM environment and installation"
-    exit 1
-}
+cd ${0%/*} || exit 1    # Run from this directory
 
 # Parse arguments for library compilation
 . $WM_PROJECT_DIR/wmake/scripts/AllwmakeParseArguments
 
-#------------------------------------------------------------------------------
+wmakeCheckPwd "$WM_PROJECT_DIR/applications" || {
+    echo "Allwmake error: Current directory is not \$WM_PROJECT_DIR/applications"
+    echo "    The environment variables are inconsistent with the installation."
+    echo "    Check the OpenFOAM entries in your dot-files and source them."
+    exit 1
+}
+
+[ -n "$FOAM_EXT_LIBBIN" ] || {
+    echo "Allwmake error: FOAM_EXT_LIBBIN not set"
+    echo "    Check the OpenFOAM entries in your dot-files and source them."
+    exit 1
+}
 
 wmake -all $targetType solvers
 wmake -all $targetType utilities

applications/solvers/basic/laplacianFoam/Make/options

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

applications/solvers/basic/laplacianFoam/createFields.H

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

applications/solvers/basic/laplacianFoam/overLaplacianDyMFoam/Make/files

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

applications/solvers/basic/laplacianFoam/overLaplacianDyMFoam/Make/options

@@ -1,8 +0,0 @@
-EXE_INC = \
-    -I$(LIB_SRC)/finiteVolume/lnInclude \
-    -I$(LIB_SRC)/dynamicFvMesh/lnInclude \
-    -I$(LIB_SRC)/overset/lnInclude \
-    -I$(LIB_SRC)/meshTools/lnInclude
-
-EXE_LIBS = \
-    -loverset

applications/solvers/basic/laplacianFoam/overLaplacianDyMFoam/createFields.H

@@ -1,53 +0,0 @@
-    Info<< "Reading field T\n" << endl;
-
-    volScalarField T
-    (
-        IOobject
-        (
-            "T",
-            runTime.timeName(),
-            mesh,
-            IOobject::MUST_READ,
-            IOobject::AUTO_WRITE
-        ),
-        mesh
-    );
-
-    // Add overset specific interpolations
-    {
-        dictionary oversetDict;
-        oversetDict.add("T", true);
-
-        const_cast<dictionary&>
-        (
-            mesh.schemesDict()
-        ).add
-        (
-            "oversetInterpolationRequired",
-            oversetDict,
-            true
-        );
-    }
-
-
-    Info<< "Reading transportProperties\n" << endl;
-
-    IOdictionary transportProperties
-    (
-        IOobject
-        (
-            "transportProperties",
-            runTime.constant(),
-            mesh,
-            IOobject::MUST_READ_IF_MODIFIED,
-            IOobject::NO_WRITE
-        )
-    );
-
-
-    Info<< "Reading diffusivity DT\n" << endl;
-
-    dimensionedScalar DT
-    (
-        transportProperties.lookup("DT")
-    );

applications/solvers/basic/laplacianFoam/overLaplacianDyMFoam/laplacianDyMFoam.C

@@ -1,110 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2015 OpenFOAM Foundation
-     \\/     M anipulation  | Copyright (C) 2016-2017 OpenCFD Ltd.
--------------------------------------------------------------------------------
-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
-    laplacianFoam
-
-Group
-    grpBasicSolvers
-
-Description
-    Laplace equation solver for a scalar quantity.
-
-    \heading Solver details
-    The solver is applicable to, e.g. for thermal diffusion in a solid.  The
-    equation is given by:
-
-    \f[
-        \ddt{T}  = \div \left( D_T \grad T \right)
-    \f]
-
-    Where:
-    \vartable
-        T     | Scalar field which is solved for, e.g. temperature
-        D_T   | Diffusion coefficient
-    \endvartable
-
-    \heading Required fields
-    \plaintable
-        T     | Scalar field which is solved for, e.g. temperature
-    \endplaintable
-
-\*---------------------------------------------------------------------------*/
-
-#include "fvCFD.H"
-#include "fvOptions.H"
-#include "simpleControl.H"
-#include "dynamicFvMesh.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-int main(int argc, char *argv[])
-{
-    #include "setRootCase.H"
-
-    #include "createTime.H"
-    #include "createNamedDynamicFvMesh.H"
-
-    simpleControl simple(mesh);
-
-    #include "createFields.H"
-    #include "createFvOptions.H"
-
-    // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-    Info<< "\nCalculating temperature distribution\n" << endl;
-
-    while (simple.loop())
-    {
-        Info<< "Time = " << runTime.timeName() << nl << endl;
-
-        mesh.update();
-
-        while (simple.correctNonOrthogonal())
-        {
-            fvScalarMatrix TEqn
-            (
-                fvm::ddt(T) - fvm::laplacian(DT, T)
-             ==
-                fvOptions(T)
-            );
-
-            fvOptions.constrain(TEqn);
-            TEqn.solve();
-            fvOptions.correct(T);
-        }
-
-        #include "write.H"
-
-        Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
-            << "  ClockTime = " << runTime.elapsedClockTime() << " s"
-            << nl << endl;
-    }
-
-    Info<< "End\n" << endl;
-
-    return 0;
-}
-
-
-// ************************************************************************* //

applications/solvers/basic/laplacianFoam/overLaplacianDyMFoam/write.H

@@ -1,46 +0,0 @@
-    if (runTime.outputTime())
-    {
-        volVectorField gradT(fvc::grad(T));
-
-        volScalarField gradTx
-        (
-            IOobject
-            (
-                "gradTx",
-                runTime.timeName(),
-                mesh,
-                IOobject::NO_READ,
-                IOobject::AUTO_WRITE
-            ),
-            gradT.component(vector::X)
-        );
-
-        volScalarField gradTy
-        (
-            IOobject
-            (
-                "gradTy",
-                runTime.timeName(),
-                mesh,
-                IOobject::NO_READ,
-                IOobject::AUTO_WRITE
-            ),
-            gradT.component(vector::Y)
-        );
-
-        volScalarField gradTz
-        (
-            IOobject
-            (
-                "gradTz",
-                runTime.timeName(),
-                mesh,
-                IOobject::NO_READ,
-                IOobject::AUTO_WRITE
-            ),
-            gradT.component(vector::Z)
-        );
-
-
-        runTime.write();
-    }

applications/solvers/basic/scalarTransportFoam/createFields.H

@@ -1,57 +1,55 @@
-Info<< "Reading field T\n" << endl;
+    Info<< "Reading field T\n" << endl;
 
-volScalarField T
-(
-    IOobject
+    volScalarField T
     (
-        "T",
-        runTime.timeName(),
-        mesh,
-        IOobject::MUST_READ,
-        IOobject::AUTO_WRITE
-    ),
-    mesh
-);
+        IOobject
+        (
+            "T",
+            runTime.timeName(),
+            mesh,
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh
+    );
 
 
-Info<< "Reading field U\n" << endl;
+    Info<< "Reading field U\n" << endl;
 
-volVectorField U
-(
-    IOobject
+    volVectorField U
     (
-        "U",
-        runTime.timeName(),
-        mesh,
-        IOobject::MUST_READ,
-        IOobject::AUTO_WRITE
-    ),
-    mesh
-);
+        IOobject
+        (
+            "U",
+            runTime.timeName(),
+            mesh,
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh
+    );
 
 
-Info<< "Reading transportProperties\n" << endl;
+    Info<< "Reading transportProperties\n" << endl;
 
-IOdictionary transportProperties
-(
-    IOobject
+    IOdictionary transportProperties
     (
-        "transportProperties",
-        runTime.constant(),
-        mesh,
-        IOobject::MUST_READ_IF_MODIFIED,
-        IOobject::NO_WRITE
-    )
-);
+        IOobject
+        (
+            "transportProperties",
+            runTime.constant(),
+            mesh,
+            IOobject::MUST_READ_IF_MODIFIED,
+            IOobject::NO_WRITE
+        )
+    );
 
 
-Info<< "Reading diffusivity DT\n" << endl;
+    Info<< "Reading diffusivity DT\n" << endl;
 
-dimensionedScalar DT
-(
-    "DT",
-    dimArea/dimTime,
-    transportProperties
-);
+    dimensionedScalar DT
+    (
+        transportProperties.lookup("DT")
+    );
 
-#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)/surfMesh/lnInclude \
-    -I$(LIB_SRC)/meshTools/lnInclude
+    -I$(LIB_SRC)/meshTools/lnInclude \
+    -I$(LIB_SRC)/triSurface/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-2017 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2015 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -79,7 +79,7 @@ Foam::PDRDragModel::~PDRDragModel()
 
 bool Foam::PDRDragModel::read(const dictionary& PDRProperties)
 {
-    PDRDragModelCoeffs_ = PDRProperties.optionalSubDict(type() + "Coeffs");
+    PDRDragModelCoeffs_ = PDRProperties.subDict(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.found())
+    if (cstrIter == dictionaryConstructorTablePtr_->end())
     {
         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-2017 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -67,7 +67,7 @@ Foam::XiEqModel::~XiEqModel()
 
 bool Foam::XiEqModel::read(const dictionary& XiEqProperties)
 {
-    XiEqModelCoeffs_ = XiEqProperties.optionalSubDict(type() + "Coeffs");
+    XiEqModelCoeffs_ = XiEqProperties.subDict(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.found())
+    if (cstrIter == dictionaryConstructorTablePtr_->end())
     {
         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-2017 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2015 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -67,7 +67,7 @@ Foam::XiGModel::~XiGModel()
 
 bool Foam::XiGModel::read(const dictionary& XiGProperties)
 {
-    XiGModelCoeffs_ = XiGProperties.optionalSubDict(type() + "Coeffs");
+    XiGModelCoeffs_ = XiGProperties.subDict(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.found())
+    if (cstrIter == dictionaryConstructorTablePtr_->end())
     {
         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-2017 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2015 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -85,7 +85,7 @@ Foam::XiModel::~XiModel()
 
 bool Foam::XiModel::read(const dictionary& XiProperties)
 {
-    XiModelCoeffs_ = XiProperties.optionalSubDict(type() + "Coeffs");
+    XiModelCoeffs_ = XiProperties.subDict(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.found())
+    if (cstrIter == dictionaryConstructorTablePtr_->end())
     {
         FatalErrorInFunction
             << "Unknown XiModel type "

applications/solvers/combustion/PDRFoam/createFields.H

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

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-2017 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -80,7 +80,7 @@ Foam::laminarFlameSpeedModels::SCOPE::SCOPE
                   dict.lookup("fuelFile")
               )
           )()
-        ).optionalSubDict(typeName + "Coeffs")
+        ).subDict(typeName + "Coeffs")
     ),
     LFL_(readScalar(coeffsDict_.lookup("lowerFlamabilityLimit"))),
     UFL_(readScalar(coeffsDict_.lookup("upperFlamabilityLimit"))),

applications/solvers/combustion/XiFoam/createFields.H

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

applications/solvers/combustion/chemFoam/createBaseFields.H

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

applications/solvers/combustion/chemFoam/output.H

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

applications/solvers/combustion/chemFoam/setDeltaT.H

@@ -3,3 +3,4 @@ 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());
-scalar Qdot = chemistry.Qdot()()[0]/rho[0];
-integratedHeat += Qdot*runTime.deltaT().value();
+    dtChem = chemistry.solve(runTime.deltaT().value());
+    scalar Sh = chemistry.Sh()()[0]/rho[0];
+    integratedHeat += Sh*runTime.deltaT().value();

applications/solvers/combustion/coldEngineFoam/createFields.H

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

applications/solvers/combustion/fireFoam/Make/options

@@ -13,7 +13,10 @@ EXE_INC = \
     -I$(LIB_SRC)/thermophysicalModels/chemistryModel/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/solidChemistryModel/lnInclude \
     -I$(LIB_SRC)/combustionModels/lnInclude \
-    -I$(LIB_SRC)/thermophysicalModels/thermophysicalProperties/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/properties/solidProperties/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 \
@@ -35,6 +38,9 @@ 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();
-    Qdot = combustion->Qdot();
+    dQ = combustion->dQ();
     volScalarField Yt(0.0*Y[0]);
 
     forAll(Y, i)
@@ -67,8 +67,8 @@ tmp<fv::convectionScheme<scalar>> mvConvection
         )
       - fvm::laplacian(turbulence->alphaEff(), he)
      ==
-        Qdot
-      + radiation->Sh(thermo, he)
+        combustion->Sh()
+      + radiation->Sh(thermo)
       + parcels.Sh(he)
       + surfaceFilm.Sh()
       + fvOptions(rho, he)

applications/solvers/combustion/fireFoam/createFields.H

@@ -131,23 +131,36 @@ Switch solvePyrolysisRegion
     additionalControlsDict.lookupOrDefault<bool>("solvePyrolysisRegion", true)
 );
 
-volScalarField Qdot
+volScalarField dQ
 (
     IOobject
     (
-        "Qdot",
+        "dQ",
         runTime.timeName(),
         mesh,
-        IOobject::READ_IF_PRESENT,
+        IOobject::NO_READ,
         IOobject::AUTO_WRITE
     ),
     mesh,
-    dimensionedScalar("Qdot", dimEnergy/dimVolume/dimTime, 0.0)
+    dimensionedScalar("dQ", dimEnergy/dimTime, 0.0)
 );
 
-#include "createDpdt.H"
 
-#include "createK.H"
+Info<< "Creating field dpdt\n" << endl;
+volScalarField dpdt
+(
+    IOobject
+    (
+        "dpdt",
+        runTime.timeName(),
+        mesh
+    ),
+    mesh,
+    dimensionedScalar("dpdt", p.dimensions()/dimTime, 0)
+);
+
+Info<< "Creating field kinetic energy K\n" << endl;
+volScalarField K("K", 0.5*magSqr(U));
 
 #include "createClouds.H"
 #include "createSurfaceFilmModel.H"

applications/solvers/combustion/reactingFoam/EEqn.H

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

applications/solvers/combustion/reactingFoam/YEqn.H

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

applications/solvers/combustion/reactingFoam/createFields.H

@@ -52,7 +52,27 @@ volScalarField& p = thermo.p();
 
 #include "compressibleCreatePhi.H"
 
-pressureControl pressureControl(p, rho, pimple.dict(), false);
+dimensionedScalar rhoMax
+(
+    dimensionedScalar::lookupOrDefault
+    (
+        "rhoMax",
+        pimple.dict(),
+        dimDensity,
+        GREAT
+    )
+);
+
+dimensionedScalar rhoMin
+(
+    dimensionedScalar::lookupOrDefault
+    (
+        "rhoMin",
+        pimple.dict(),
+        dimDensity,
+        0
+    )
+);
 
 mesh.setFluxRequired(p.name());
 
@@ -71,6 +91,24 @@ autoPtr<compressible::turbulenceModel> turbulence
 // Set the turbulence into the reaction model
 reaction->setTurbulence(turbulence());
 
+
+Info<< "Creating field dpdt\n" << endl;
+volScalarField dpdt
+(
+    IOobject
+    (
+        "dpdt",
+        runTime.timeName(),
+        mesh
+    ),
+    mesh,
+    dimensionedScalar("dpdt", p.dimensions()/dimTime, 0)
+);
+
+Info<< "Creating field kinetic energy K\n" << endl;
+volScalarField K("K", 0.5*magSqr(U));
+
+
 multivariateSurfaceInterpolationScheme<scalar>::fieldTable fields;
 
 forAll(Y, i)
@@ -79,22 +117,18 @@ forAll(Y, i)
 }
 fields.add(thermo.he());
 
-volScalarField Qdot
+volScalarField dQ
 (
     IOobject
     (
-        "Qdot",
+        "dQ",
         runTime.timeName(),
         mesh,
-        IOobject::READ_IF_PRESENT,
+        IOobject::NO_READ,
         IOobject::AUTO_WRITE
     ),
     mesh,
-    dimensionedScalar("Qdot", dimEnergy/dimVolume/dimTime, 0.0)
+    dimensionedScalar("dQ", dimEnergy/dimTime, 0.0)
 );
 
-#include "createDpdt.H"
-
-#include "createK.H"
-
 #include "createMRF.H"

applications/solvers/combustion/reactingFoam/pEqn.H

@@ -1,4 +1,7 @@
 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));
@@ -84,18 +87,19 @@ 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);
 
-if (pressureControl.limit(p))
-{
-    p.correctBoundaryConditions();
-}
-
-rho = thermo.rho();
-
 if (thermo.dpdt())
 {
     dpdt = fvc::ddt(p);

applications/solvers/combustion/reactingFoam/pcEqn.H

@@ -1,4 +1,7 @@
 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()));
@@ -106,13 +109,19 @@ U.correctBoundaryConditions();
 fvOptions.correct(U);
 K = 0.5*magSqr(U);
 
-if (pressureControl.limit(p))
-{
-    p.correctBoundaryConditions();
-    rho = thermo.rho();
-}
-
 if (thermo.dpdt())
 {
     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;

applications/solvers/combustion/reactingFoam/reactingFoam.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -37,7 +37,6 @@ Description
 #include "psiCombustionModel.H"
 #include "multivariateScheme.H"
 #include "pimpleControl.H"
-#include "pressureControl.H"
 #include "fvOptions.H"
 #include "localEulerDdtScheme.H"
 #include "fvcSmooth.H"
@@ -115,8 +114,6 @@ int main(int argc, char *argv[])
             }
         }
 
-        rho = thermo.rho();
-
         runTime.write();
 
         Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"

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

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

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

@@ -1,74 +1,76 @@
-rho = thermo.rho();
-
-// Thermodynamic density needs to be updated by psi*d(p) after the
-// pressure solution
-const volScalarField psip0(psi*p);
-
-volScalarField rAU(1.0/UEqn.A());
-surfaceScalarField rhorAUf("rhorAUf", fvc::interpolate(rho*rAU));
-volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p));
-
-surfaceScalarField phig(-rhorAUf*ghf*fvc::snGrad(rho)*mesh.magSf());
-
-surfaceScalarField phiHbyA
-(
-    "phiHbyA",
-    (
-        fvc::flux(rho*HbyA)
-      + rhorAUf*fvc::ddtCorr(rho, U, phi)
-    )
-  + phig
-);
-
-MRF.makeRelative(fvc::interpolate(rho), phiHbyA);
-
-// Update the pressure BCs to ensure flux consistency
-constrainPressure(p_rgh, rho, U, phiHbyA, rhorAUf, MRF);
-
-fvScalarMatrix p_rghDDtEqn
-(
-    fvc::ddt(rho) + psi*correction(fvm::ddt(p_rgh))
-  + fvc::div(phiHbyA)
-  ==
-    fvOptions(psi, p_rgh, rho.name())
-);
-
-while (pimple.correctNonOrthogonal())
 {
-    fvScalarMatrix p_rghEqn
+    rho = thermo.rho();
+
+    // Thermodynamic density needs to be updated by psi*d(p) after the
+    // pressure solution - done in 2 parts. Part 1:
+    thermo.rho() -= psi*p;
+
+    volScalarField rAU(1.0/UEqn.A());
+    surfaceScalarField rhorAUf("rhorAUf", fvc::interpolate(rho*rAU));
+    volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p));
+
+    surfaceScalarField phig(-rhorAUf*ghf*fvc::snGrad(rho)*mesh.magSf());
+
+    surfaceScalarField phiHbyA
     (
-        p_rghDDtEqn
-      - fvm::laplacian(rhorAUf, p_rgh)
+        "phiHbyA",
+        (
+            fvc::flux(rho*HbyA)
+          + rhorAUf*fvc::ddtCorr(rho, U, phi)
+        )
+      + phig
     );
 
-    p_rghEqn.solve(mesh.solver(p_rgh.select(pimple.finalInnerIter())));
+    MRF.makeRelative(fvc::interpolate(rho), phiHbyA);
 
-    if (pimple.finalNonOrthogonalIter())
+    // Update the pressure BCs to ensure flux consistency
+    constrainPressure(p_rgh, rho, U, phiHbyA, rhorAUf, MRF);
+
+    fvScalarMatrix p_rghDDtEqn
+    (
+        fvc::ddt(rho) + psi*correction(fvm::ddt(p_rgh))
+      + fvc::div(phiHbyA)
+      ==
+        fvOptions(psi, p_rgh, rho.name())
+    );
+
+    while (pimple.correctNonOrthogonal())
     {
-        // Calculate the conservative fluxes
-        phi = phiHbyA + p_rghEqn.flux();
+        fvScalarMatrix p_rghEqn
+        (
+            p_rghDDtEqn
+          - fvm::laplacian(rhorAUf, p_rgh)
+        );
 
-        // Explicitly relax pressure for momentum corrector
-        p_rgh.relax();
+        p_rghEqn.solve(mesh.solver(p_rgh.select(pimple.finalInnerIter())));
 
-        // Correct the momentum source with the pressure gradient flux
-        // calculated from the relaxed pressure
-        U = HbyA + rAU*fvc::reconstruct((phig + p_rghEqn.flux())/rhorAUf);
-        U.correctBoundaryConditions();
-        fvOptions.correct(U);
-        K = 0.5*magSqr(U);
+        if (pimple.finalNonOrthogonalIter())
+        {
+            // Calculate the conservative fluxes
+            phi = phiHbyA + p_rghEqn.flux();
+
+            // Explicitly relax pressure for momentum corrector
+            p_rgh.relax();
+
+            // Correct the momentum source with the pressure gradient flux
+            // calculated from the relaxed pressure
+            U = HbyA + rAU*fvc::reconstruct((phig + p_rghEqn.flux())/rhorAUf);
+            U.correctBoundaryConditions();
+            fvOptions.correct(U);
+            K = 0.5*magSqr(U);
+        }
     }
+
+    p = p_rgh + rho*gh;
+
+    // Second part of thermodynamic density update
+    thermo.rho() += psi*p;
+
+    if (thermo.dpdt())
+    {
+        dpdt = fvc::ddt(p);
+    }
+
+    #include "rhoEqn.H"
+    #include "compressibleContinuityErrs.H"
 }
-
-p = p_rgh + rho*gh;
-
-// Thermodynamic density update
-thermo.correctRho(psi*p - psip0);
-
-if (thermo.dpdt())
-{
-    dpdt = fvc::ddt(p);
-}
-
-#include "rhoEqn.H"
-#include "compressibleContinuityErrs.H"

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

@@ -52,8 +52,6 @@ volScalarField& p = thermo.p();
 
 #include "compressibleCreatePhi.H"
 
-pressureControl pressureControl(p, rho, pimple.dict(), false);
-
 mesh.setFluxRequired(p.name());
 
 
@@ -72,6 +70,24 @@ autoPtr<compressible::turbulenceModel> turbulence
 // Set the turbulence into the reaction model
 reaction->setTurbulence(turbulence());
 
+
+Info<< "Creating field dpdt\n" << endl;
+volScalarField dpdt
+(
+    IOobject
+    (
+        "dpdt",
+        runTime.timeName(),
+        mesh
+    ),
+    mesh,
+    dimensionedScalar("dpdt", p.dimensions()/dimTime, 0)
+);
+
+Info<< "Creating field kinetic energy K\n" << endl;
+volScalarField K("K", 0.5*magSqr(U));
+
+
 multivariateSurfaceInterpolationScheme<scalar>::fieldTable fields;
 
 forAll(Y, i)
@@ -80,45 +96,18 @@ forAll(Y, i)
 }
 fields.add(thermo.he());
 
-volScalarField Qdot
+volScalarField dQ
 (
     IOobject
     (
-        "Qdot",
+        "dQ",
         runTime.timeName(),
         mesh,
-        IOobject::READ_IF_PRESENT,
+        IOobject::NO_READ,
         IOobject::AUTO_WRITE
     ),
     mesh,
-    dimensionedScalar("Qdot", dimEnergy/dimVolume/dimTime, 0.0)
+    dimensionedScalar("dQ", dimEnergy/dimTime, 0.0)
 );
 
-#include "createDpdt.H"
-
-#include "createK.H"
-
 #include "createMRF.H"
-
-dimensionedScalar rhoMax
-(
-    dimensionedScalar::lookupOrDefault
-    (
-        "rhoMax",
-        pimple.dict(),
-        dimDensity,
-        GREAT
-    )
-);
-
-dimensionedScalar rhoMin
-(
-    dimensionedScalar::lookupOrDefault
-    (
-        "rhoMin",
-        pimple.dict(),
-        dimDensity,
-        0
-    )
-);
-

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

@@ -0,0 +1,109 @@
+{
+    rho = thermo.rho();
+
+    // Thermodynamic density needs to be updated by psi*d(p) after the
+    // pressure solution - done in 2 parts. Part 1:
+    thermo.rho() -= psi*p;
+
+    volScalarField rAU(1.0/UEqn.A());
+    surfaceScalarField rhorAUf("rhorAUf", fvc::interpolate(rho*rAU));
+    volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p));
+
+    if (pimple.transonic())
+    {
+        surfaceScalarField phiHbyA
+        (
+            "phiHbyA",
+            (
+                fvc::flux(HbyA)
+              + rhorAUf*fvc::ddtCorr(rho, U, phi)/fvc::interpolate(rho)
+            )
+        );
+
+        MRF.makeRelative(phiHbyA);
+
+        surfaceScalarField phid("phid", fvc::interpolate(thermo.psi())*phiHbyA);
+
+        phiHbyA *= fvc::interpolate(rho);
+
+        fvScalarMatrix pDDtEqn
+        (
+            fvc::ddt(rho) + fvc::div(phiHbyA)
+          + correction(psi*fvm::ddt(p) + fvm::div(phid, p))
+        );
+
+        while (pimple.correctNonOrthogonal())
+        {
+            fvScalarMatrix pEqn
+            (
+                pDDtEqn
+              - fvm::laplacian(rhorAUf, p)
+             ==
+                fvOptions(psi, p, rho.name())
+            );
+
+            pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter())));
+
+            if (pimple.finalNonOrthogonalIter())
+            {
+                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);
+
+        fvScalarMatrix pDDtEqn
+        (
+            fvc::ddt(rho) + psi*correction(fvm::ddt(p))
+          + fvc::div(phiHbyA)
+         ==
+            fvOptions(psi, p, rho.name())
+        );
+
+        while (pimple.correctNonOrthogonal())
+        {
+            fvScalarMatrix pEqn
+            (
+                pDDtEqn
+              - fvm::laplacian(rhorAUf, p)
+            );
+
+            pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter())));
+
+            if (pimple.finalNonOrthogonalIter())
+            {
+                phi = phiHbyA + pEqn.flux();
+            }
+        }
+    }
+
+    // Second part of thermodynamic density update
+    thermo.rho() += psi*p;
+
+    #include "rhoEqn.H"
+    #include "compressibleContinuityErrs.H"
+
+    U = HbyA - rAU*fvc::grad(p);
+    U.correctBoundaryConditions();
+    fvOptions.correct(U);
+    K = 0.5*magSqr(U);
+
+    if (thermo.dpdt())
+    {
+        dpdt = fvc::ddt(p);
+    }
+}

applications/solvers/combustion/reactingFoam/rhoReactingFoam/rhoReactingFoam.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -38,7 +38,6 @@ Description
 #include "turbulentFluidThermoModel.H"
 #include "multivariateScheme.H"
 #include "pimpleControl.H"
-#include "pressureControl.H"
 #include "fvOptions.H"
 #include "localEulerDdtScheme.H"
 #include "fvcSmooth.H"
@@ -101,14 +100,7 @@ int main(int argc, char *argv[])
             // --- Pressure corrector loop
             while (pimple.correct())
             {
-                if (pimple.consistent())
-                {
-                    #include "../../../compressible/rhoPimpleFoam/pcEqn.H"
-                }
-                else
-                {
-                    #include "../../../compressible/rhoPimpleFoam/pEqn.H"
-                }
+                #include "pEqn.H"
             }
 
             if (pimple.turbCorr())

applications/solvers/combustion/reactingFoam/setRDeltaT.H

@@ -43,16 +43,13 @@ License
     // Damping coefficient (1-0)
     scalar rDeltaTDampingCoeff
     (
-        pimpleDict.lookupOrDefault<scalar>("rDeltaTDampingCoeff", 1.0)
+        pimpleDict.lookupOrDefault<scalar>("rDeltaTDampingCoeff", 1)
     );
 
     // 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;
 
@@ -71,89 +68,34 @@ License
         rDeltaT.max(1/maxDeltaT);
 
         Info<< "    Flow        = "
-            << 1/gMax(rDeltaT.primitiveField()) << ", "
-            << 1/gMin(rDeltaT.primitiveField()) << endl;
+            << gMin(1/rDeltaT.primitiveField()) << ", "
+            << gMax(1/rDeltaT.primitiveField()) << endl;
     }
 
-    // Heat release rate time scale
-    if (alphaTemp < 1)
+    // Reaction source time scale
+    if (alphaTemp < 1.0)
     {
         volScalarField::Internal rDeltaTT
         (
-            mag(Qdot)/(alphaTemp*rho*thermo.Cp()*T)
+            mag(reaction->Sh())/(alphaTemp*rho*thermo.Cp()*T)
         );
 
         Info<< "    Temperature = "
-            << 1/(gMax(rDeltaTT.field()) + VSMALL) << ", "
-            << 1/(gMin(rDeltaTT.field()) + VSMALL) << endl;
+            << gMin(1/(rDeltaTT.field() + VSMALL)) << ", "
+            << gMax(1/(rDeltaTT.field() + VSMALL)) << endl;
 
-        rDeltaT.ref() = max(rDeltaT(), rDeltaTT);
-    }
-
-    // Reaction rate time scale
-    if (alphaY < 1)
-    {
-        dictionary Yref(pimpleDict.subDict("Yref"));
-
-        volScalarField::Internal rDeltaTY
+        rDeltaT.ref() = max
         (
-            IOobject
-            (
-                "rDeltaTY",
-                runTime.timeName(),
-                mesh
-            ),
-            mesh,
-            dimensionedScalar("rDeltaTY", rDeltaT.dimensions(), 0)
+            rDeltaT(),
+            rDeltaTT
         );
-
-        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)
+    if (rDeltaTSmoothingCoeff < 1.0)
     {
         fvc::smooth(rDeltaT, rDeltaTSmoothingCoeff);
     }
@@ -163,7 +105,7 @@ License
     // - only increase at a fraction of old time scale
     if
     (
-        rDeltaTDampingCoeff < 1
+        rDeltaTDampingCoeff < 1.0
      && runTime.timeIndex() > runTime.startTimeIndex() + 1
     )
     {
@@ -178,8 +120,8 @@ License
     rDeltaT.correctBoundaryConditions();
 
     Info<< "    Overall     = "
-        << 1/gMax(rDeltaT.primitiveField())
-        << ", " << 1/gMin(rDeltaT.primitiveField()) << endl;
+        << gMin(1/rDeltaT.primitiveField())
+        << ", " << gMax(1/rDeltaT.primitiveField()) << endl;
 }
 
 

applications/solvers/compressible/rhoCentralFoam/BCs/Make/files

@@ -1,3 +1,4 @@
+mixedFixedValueSlip/mixedFixedValueSlipFvPatchFields.C
 U/maxwellSlipUFvPatchVectorField.C
 T/smoluchowskiJumpTFvPatchScalarField.C
 rho/fixedRhoFvPatchScalarField.C

applications/solvers/compressible/rhoCentralFoam/BCs/U/maxwellSlipUFvPatchVectorField.C

@@ -38,7 +38,7 @@ Foam::maxwellSlipUFvPatchVectorField::maxwellSlipUFvPatchVectorField
     const DimensionedField<vector, volMesh>& iF
 )
 :
-    partialSlipFvPatchVectorField(p, iF),
+    mixedFixedValueSlipFvPatchVectorField(p, iF),
     TName_("T"),
     rhoName_("rho"),
     psiName_("thermo:psi"),
@@ -59,7 +59,7 @@ Foam::maxwellSlipUFvPatchVectorField::maxwellSlipUFvPatchVectorField
     const fvPatchFieldMapper& mapper
 )
 :
-    partialSlipFvPatchVectorField(mspvf, p, iF, mapper),
+    mixedFixedValueSlipFvPatchVectorField(mspvf, p, iF, mapper),
     TName_(mspvf.TName_),
     rhoName_(mspvf.rhoName_),
     psiName_(mspvf.psiName_),
@@ -79,7 +79,7 @@ Foam::maxwellSlipUFvPatchVectorField::maxwellSlipUFvPatchVectorField
     const dictionary& dict
 )
 :
-    partialSlipFvPatchVectorField(p, iF),
+    mixedFixedValueSlipFvPatchVectorField(p, iF),
     TName_(dict.lookupOrDefault<word>("T", "T")),
     rhoName_(dict.lookupOrDefault<word>("rho", "rho")),
     psiName_(dict.lookupOrDefault<word>("psi", "thermo:psi")),
@@ -132,7 +132,7 @@ Foam::maxwellSlipUFvPatchVectorField::maxwellSlipUFvPatchVectorField
     const DimensionedField<vector, volMesh>& iF
 )
 :
-    partialSlipFvPatchVectorField(mspvf, iF),
+    mixedFixedValueSlipFvPatchVectorField(mspvf, iF),
     TName_(mspvf.TName_),
     rhoName_(mspvf.rhoName_),
     psiName_(mspvf.psiName_),
@@ -193,7 +193,7 @@ void Foam::maxwellSlipUFvPatchVectorField::updateCoeffs()
         refValue() -= C1/prho*transform(I - n*n, (n & ptauMC));
     }
 
-    partialSlipFvPatchVectorField::updateCoeffs();
+    mixedFixedValueSlipFvPatchVectorField::updateCoeffs();
 }
 
 

applications/solvers/compressible/rhoCentralFoam/BCs/U/maxwellSlipUFvPatchVectorField.H

@@ -3,7 +3,7 @@
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
     \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
-     \\/     M anipulation  | Copyright (C) 2017 OpenCFD Ltd.
+     \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -36,7 +36,7 @@ SourceFiles
 #ifndef maxwellSlipUFvPatchVectorField_H
 #define maxwellSlipUFvPatchVectorField_H
 
-#include "partialSlipFvPatchFields.H"
+#include "mixedFixedValueSlipFvPatchFields.H"
 #include "Switch.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -50,7 +50,7 @@ namespace Foam
 
 class maxwellSlipUFvPatchVectorField
 :
-    public partialSlipFvPatchVectorField
+    public mixedFixedValueSlipFvPatchVectorField
 {
     // Private data
 

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

@@ -0,0 +1,196 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011-2016 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/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "mixedFixedValueSlipFvPatchField.H"
+#include "symmTransformField.H"
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+template<class Type>
+Foam::mixedFixedValueSlipFvPatchField<Type>::mixedFixedValueSlipFvPatchField
+(
+    const fvPatch& p,
+    const DimensionedField<Type, volMesh>& iF
+)
+:
+    transformFvPatchField<Type>(p, iF),
+    refValue_(p.size()),
+    valueFraction_(p.size(), 1.0)
+{}
+
+
+template<class Type>
+Foam::mixedFixedValueSlipFvPatchField<Type>::mixedFixedValueSlipFvPatchField
+(
+    const fvPatch& p,
+    const DimensionedField<Type, volMesh>& iF,
+    const dictionary& dict
+)
+:
+    transformFvPatchField<Type>(p, iF),
+    refValue_("refValue", dict, p.size()),
+    valueFraction_("valueFraction", dict, p.size())
+{}
+
+
+template<class Type>
+Foam::mixedFixedValueSlipFvPatchField<Type>::mixedFixedValueSlipFvPatchField
+(
+    const mixedFixedValueSlipFvPatchField<Type>& ptf,
+    const fvPatch& p,
+    const DimensionedField<Type, volMesh>& iF,
+    const fvPatchFieldMapper& mapper
+)
+:
+    transformFvPatchField<Type>(ptf, p, iF, mapper),
+    refValue_(ptf.refValue_, mapper),
+    valueFraction_(ptf.valueFraction_, mapper)
+{}
+
+
+template<class Type>
+Foam::mixedFixedValueSlipFvPatchField<Type>::mixedFixedValueSlipFvPatchField
+(
+    const mixedFixedValueSlipFvPatchField<Type>& ptf
+)
+:
+    transformFvPatchField<Type>(ptf),
+    refValue_(ptf.refValue_),
+    valueFraction_(ptf.valueFraction_)
+{}
+
+
+template<class Type>
+Foam::mixedFixedValueSlipFvPatchField<Type>::mixedFixedValueSlipFvPatchField
+(
+    const mixedFixedValueSlipFvPatchField<Type>& ptf,
+    const DimensionedField<Type, volMesh>& iF
+)
+:
+    transformFvPatchField<Type>(ptf, iF),
+    refValue_(ptf.refValue_),
+    valueFraction_(ptf.valueFraction_)
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+template<class Type>
+void Foam::mixedFixedValueSlipFvPatchField<Type>::autoMap
+(
+    const fvPatchFieldMapper& m
+)
+{
+    Field<Type>::autoMap(m);
+    refValue_.autoMap(m);
+    valueFraction_.autoMap(m);
+}
+
+
+template<class Type>
+void Foam::mixedFixedValueSlipFvPatchField<Type>::rmap
+(
+    const fvPatchField<Type>& ptf,
+    const labelList& addr
+)
+{
+    transformFvPatchField<Type>::rmap(ptf, addr);
+
+    const mixedFixedValueSlipFvPatchField<Type>& dmptf =
+        refCast<const mixedFixedValueSlipFvPatchField<Type>>(ptf);
+
+    refValue_.rmap(dmptf.refValue_, addr);
+    valueFraction_.rmap(dmptf.valueFraction_, addr);
+}
+
+
+template<class Type>
+Foam::tmp<Foam::Field<Type>>
+Foam::mixedFixedValueSlipFvPatchField<Type>::snGrad() const
+{
+    tmp<vectorField> nHat = this->patch().nf();
+    Field<Type> pif(this->patchInternalField());
+
+    return
+    (
+        valueFraction_*refValue_
+      + (1.0 - valueFraction_)*transform(I - sqr(nHat), pif) - pif
+    )*this->patch().deltaCoeffs();
+}
+
+
+template<class Type>
+void Foam::mixedFixedValueSlipFvPatchField<Type>::evaluate
+(
+    const Pstream::commsTypes
+)
+{
+    if (!this->updated())
+    {
+        this->updateCoeffs();
+    }
+
+    vectorField nHat(this->patch().nf());
+
+    Field<Type>::operator=
+    (
+        valueFraction_*refValue_
+      +
+        (1.0 - valueFraction_)
+       *transform(I - nHat*nHat, this->patchInternalField())
+    );
+
+    transformFvPatchField<Type>::evaluate();
+}
+
+
+template<class Type>
+Foam::tmp<Foam::Field<Type>>
+Foam::mixedFixedValueSlipFvPatchField<Type>::snGradTransformDiag() const
+{
+    vectorField nHat(this->patch().nf());
+    vectorField diag(nHat.size());
+
+    diag.replace(vector::X, mag(nHat.component(vector::X)));
+    diag.replace(vector::Y, mag(nHat.component(vector::Y)));
+    diag.replace(vector::Z, mag(nHat.component(vector::Z)));
+
+    return
+        valueFraction_*Type(pTraits<Type>::one)
+      + (1.0 - valueFraction_)
+       *transformFieldMask<Type>(pow<vector, pTraits<Type>::rank>(diag));
+}
+
+
+template<class Type>
+void Foam::mixedFixedValueSlipFvPatchField<Type>::write(Ostream& os) const
+{
+    transformFvPatchField<Type>::write(os);
+    refValue_.writeEntry("refValue", os);
+    valueFraction_.writeEntry("valueFraction", os);
+}
+
+
+// ************************************************************************* //

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

@@ -0,0 +1,236 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011-2016 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/>.
+
+Class
+    Foam::mixedFixedValueSlipFvPatchField
+
+Description
+    A mixed boundary type that blends between fixedValue and slip, as opposed
+    to the standard mixed condition that blends between fixedValue and
+    fixedGradient; required to implement maxwellSlipU condition.
+
+SourceFiles
+    mixedFixedValueSlipFvPatchField.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef mixedFixedValueSlipFvPatchField_H
+#define mixedFixedValueSlipFvPatchField_H
+
+#include "transformFvPatchField.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                      Class mixedFixedValueSlipFvPatch Declaration
+\*---------------------------------------------------------------------------*/
+
+template<class Type>
+class mixedFixedValueSlipFvPatchField
+:
+    public transformFvPatchField<Type>
+{
+    // Private data
+
+        //- Value field used for boundary condition
+        Field<Type> refValue_;
+
+        //- Fraction (0-1) of value used for boundary condition
+        scalarField valueFraction_;
+
+public:
+
+    //- Runtime type information
+    TypeName("mixedFixedValueSlip");
+
+
+    // Constructors
+
+        //- Construct from patch and internal field
+        mixedFixedValueSlipFvPatchField
+        (
+            const fvPatch&,
+            const DimensionedField<Type, volMesh>&
+        );
+
+        //- Construct from patch, internal field and dictionary
+        mixedFixedValueSlipFvPatchField
+        (
+            const fvPatch&,
+            const DimensionedField<Type, volMesh>&,
+            const dictionary&
+        );
+
+        //- Construct by mapping given mixedFixedValueSlipFvPatchField
+        //- Onto a new patch
+        mixedFixedValueSlipFvPatchField
+        (
+            const mixedFixedValueSlipFvPatchField<Type>&,
+            const fvPatch&,
+            const DimensionedField<Type, volMesh>&,
+            const fvPatchFieldMapper&
+        );
+
+        //- Construct as copy
+        mixedFixedValueSlipFvPatchField
+        (
+            const mixedFixedValueSlipFvPatchField<Type>&
+        );
+
+        //- Construct and return a clone
+        virtual tmp<fvPatchField<Type>> clone() const
+        {
+            return tmp<fvPatchField<Type>>
+            (
+                new mixedFixedValueSlipFvPatchField<Type>(*this)
+            );
+        }
+
+        //- Construct as copy setting internal field reference
+        mixedFixedValueSlipFvPatchField
+        (
+            const mixedFixedValueSlipFvPatchField<Type>&,
+            const DimensionedField<Type, volMesh>&
+        );
+
+        //- Construct and return a clone setting internal field reference
+        virtual tmp<fvPatchField<Type>> clone
+        (
+            const DimensionedField<Type, volMesh>& iF
+        ) const
+        {
+            return tmp<fvPatchField<Type>>
+            (
+                new mixedFixedValueSlipFvPatchField<Type>(*this, iF)
+            );
+        }
+
+    // Member functions
+
+        // Access
+
+            //- Return false: this patch field is not altered by assignment
+            virtual bool assignable() const
+            {
+                return false;
+            }
+
+
+        // Mapping functions
+
+            //- Map (and resize as needed) from self given a mapping object
+            virtual void autoMap
+            (
+                const fvPatchFieldMapper&
+            );
+
+            //- Reverse map the given fvPatchField onto this fvPatchField
+            virtual void rmap
+            (
+                const fvPatchField<Type>&,
+                const labelList&
+            );
+
+        // Return defining fields
+
+            virtual Field<Type>& refValue()
+            {
+                return refValue_;
+            }
+
+            virtual const Field<Type>& refValue() const
+            {
+                return refValue_;
+            }
+
+            virtual scalarField& valueFraction()
+            {
+                return valueFraction_;
+            }
+
+            virtual const scalarField& valueFraction() const
+            {
+                return valueFraction_;
+            }
+
+        // Evaluation functions
+
+            //- Return gradient at boundary
+            virtual tmp<Field<Type>> snGrad() const;
+
+            //- Evaluate the patch field
+            virtual void evaluate
+            (
+                const Pstream::commsTypes commsType=Pstream::blocking
+            );
+
+            //- Return face-gradient transform diagonal
+            virtual tmp<Field<Type>> snGradTransformDiag() const;
+
+
+        //- Write
+        virtual void write(Ostream&) const;
+
+
+    // Member operators
+
+        virtual void operator=(const UList<Type>&) {}
+
+        virtual void operator=(const fvPatchField<Type>&) {}
+        virtual void operator+=(const fvPatchField<Type>&) {}
+        virtual void operator-=(const fvPatchField<Type>&) {}
+        virtual void operator*=(const fvPatchField<scalar>&) {}
+        virtual void operator/=(const fvPatchField<scalar>&) {}
+
+        virtual void operator+=(const Field<Type>&) {}
+        virtual void operator-=(const Field<Type>&) {}
+
+        virtual void operator*=(const Field<scalar>&) {}
+        virtual void operator/=(const Field<scalar>&) {}
+
+        virtual void operator=(const Type&) {}
+        virtual void operator+=(const Type&) {}
+        virtual void operator-=(const Type&) {}
+        virtual void operator*=(const scalar) {}
+        virtual void operator/=(const scalar) {}
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+    #include "mixedFixedValueSlipFvPatchField.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

applications/solvers/compressible/rhoCentralFoam/BCs/mixedFixedValueSlip/mixedFixedValueSlipFvPatchFields.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2016-2017 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -23,7 +23,7 @@ License
 
 \*---------------------------------------------------------------------------*/
 
-#include "oversetFvPatchFields.H"
+#include "mixedFixedValueSlipFvPatchFields.H"
 #include "addToRunTimeSelectionTable.H"
 #include "volFields.H"
 
@@ -34,7 +34,7 @@ namespace Foam
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
-makePatchFields(overset);
+makePatchFields(mixedFixedValueSlip);
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 

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

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2016-2017 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -23,10 +23,10 @@ License
 
 \*---------------------------------------------------------------------------*/
 
-#ifndef oversetFvPatchFields_H
-#define oversetFvPatchFields_H
+#ifndef mixedFixedValueSlipFvPatchFields_H
+#define mixedFixedValueSlipFvPatchFields_H
 
-#include "oversetFvPatchField.H"
+#include "mixedFixedValueSlipFvPatchField.H"
 #include "fieldTypes.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -36,7 +36,7 @@ namespace Foam
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-makePatchTypeFieldTypedefs(overset);
+makePatchTypeFieldTypedefs(mixedFixedValueSlip);
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 

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

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2017 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -21,21 +21,12 @@ License
     You should have received a copy of the GNU General Public License
     along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
 
-Typedef
-    Foam::triadField
-
-Description
-    Forward declarations of the specialisation of Field\<T\> for triad.
-
-SourceFiles
-    triadField.C
-
 \*---------------------------------------------------------------------------*/
 
-#ifndef triadFieldFwd_H
-#define triadFieldFwd_H
+#ifndef mixedFixedValueSlipFvPatchFieldsFwd_H
+#define mixedFixedValueSlipFvPatchFieldsFwd_H
 
-#include "triad.H"
+#include "fieldTypes.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -44,8 +35,9 @@ namespace Foam
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-template<class Type> class Field;
-typedef Field<triad> triadField;
+template<class Type> class mixedFixedValueSlipFvPatchField;
+
+makePatchTypeFieldTypedefs(mixedFixedValueSlip);
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 

applications/solvers/compressible/rhoCentralFoam/rhoCentralDyMFoam/rhoCentralDyMFoam.C

@@ -120,14 +120,11 @@ int main(int argc, char *argv[])
             "cSf_pos",
             interpolate(c, pos, T.name())*mesh.magSf()
         );
-        cSf_pos.setOriented();
-
         surfaceScalarField cSf_neg
         (
             "cSf_neg",
             interpolate(c, neg, T.name())*mesh.magSf()
         );
-        cSf_neg.setOriented();
 
         surfaceScalarField ap
         (
@@ -272,5 +269,4 @@ int main(int argc, char *argv[])
     return 0;
 }
 
-
 // ************************************************************************* //

applications/solvers/compressible/rhoCentralFoam/rhoCentralFoam.C

@@ -101,14 +101,11 @@ int main(int argc, char *argv[])
             "cSf_pos",
             interpolate(c, pos, T.name())*mesh.magSf()
         );
-        cSf_pos.setOriented();
-
         surfaceScalarField cSf_neg
         (
             "cSf_neg",
             interpolate(c, neg, T.name())*mesh.magSf()
         );
-        cSf_neg.setOriented();
 
         surfaceScalarField ap
         (

applications/solvers/compressible/rhoPimpleAdiabaticFoam/EEqn.H

@@ -1,50 +0,0 @@
-{
-    volScalarField& he = thermo.he();
-
-    const tmp<volScalarField>& tCp = thermo.Cp();
-    const tmp<volScalarField>& tCv = thermo.Cv();
-
-    const volScalarField& Cp = tCp();
-    const volScalarField& Cv = tCv();
-    const scalar gamma = max(Cp/Cv).value();
-
-    if (mag(gamma - min(Cp/Cv).value()) > VSMALL)
-    {
-        notImplemented("gamma not constant in space");
-    }
-
-    const dictionary& thermoDict = thermo.subDict("mixture");
-
-    const dictionary& eosDict = thermoDict.subDict("equationOfState");
-
-    bool local = eosDict.lookupOrDefault<bool>("local", false);
-
-    // Evolve T as:
-    //
-    // T_1 = T_0 \frac{p}{p_0}^{\frac{\gamma - 1}{\gamma}}
-
-    if (!local)
-    {
-        const scalar T0 = readScalar(eosDict.lookup("T0"));
-        const scalar p0 = readScalar(eosDict.lookup("p0"));
-
-        he = thermo.he(p, pow(p/p0, (gamma - scalar(1))/gamma)*T0);
-    }
-    else
-    {
-        const volScalarField& T0 = T.oldTime();
-        const volScalarField& p0 = p.oldTime();
-
-        he = thermo.he(p, pow(p/p0, (gamma - scalar(1))/gamma)*T0);
-    }
-
-    thermo.correct();
-
-    psi = 1.0/((Cp - Cv)*T);
-
-    rho = thermo.rho();
-    rho.relax();
-
-    rho.writeMinMax(Info);
-}
-

applications/solvers/compressible/rhoPimpleAdiabaticFoam/Make/files

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

applications/solvers/compressible/rhoPimpleAdiabaticFoam/Make/options

@@ -1,20 +0,0 @@
-EXE_INC = \
-    -I$(LIB_SRC)/transportModels/compressible/lnInclude \
-    -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
-    -I$(LIB_SRC)/TurbulenceModels/turbulenceModels/lnInclude \
-    -I$(LIB_SRC)/TurbulenceModels/compressible/lnInclude \
-    -I$(LIB_SRC)/finiteVolume/cfdTools \
-    -I$(LIB_SRC)/finiteVolume/lnInclude \
-    -I$(LIB_SRC)/meshTools/lnInclude \
-    -I$(LIB_SRC)/sampling/lnInclude \
-
-EXE_LIBS = \
-    -lcompressibleTransportModels \
-    -lfluidThermophysicalModels \
-    -lspecie \
-    -lturbulenceModels \
-    -lcompressibleTurbulenceModels \
-    -lfiniteVolume \
-    -lmeshTools \
-    -lsampling \
-    -lfvOptions

applications/solvers/compressible/rhoPimpleAdiabaticFoam/UEqn.H

@@ -1,24 +0,0 @@
-// Solve the Momentum equation
-
-MRF.correctBoundaryVelocity(U);
-
-tmp<fvVectorMatrix> tUEqn
-(
-    fvm::ddt(rho, U) + fvm::div(phi, U)
-  + MRF.DDt(rho, U)
-  + turbulence->divDevRhoReff(U)
- ==
-    fvOptions(rho, U)
-);
-fvVectorMatrix& UEqn = tUEqn.ref();
-
-UEqn.relax();
-
-fvOptions.constrain(UEqn);
-
-if (pimple.momentumPredictor())
-{
-    solve(UEqn == -fvc::grad(p));
-
-    fvOptions.correct(U);
-}

applications/solvers/compressible/rhoPimpleAdiabaticFoam/createFields.H

@@ -1,102 +0,0 @@
-Info<< "Reading thermophysical properties\n" << endl;
-
-autoPtr<fluidThermo> pThermo
-(
-    fluidThermo::New(mesh)
-);
-fluidThermo& thermo = pThermo();
-thermo.validate(args.executable(), "h", "e");
-
-volScalarField& p = thermo.p();
-volScalarField& T = thermo.T();
-
-volScalarField rho
-(
-    IOobject
-    (
-        "rho",
-        runTime.timeName(),
-        mesh,
-        IOobject::READ_IF_PRESENT,
-        IOobject::AUTO_WRITE
-    ),
-    thermo.rho()
-);
-
-Info<< "Reading field U\n" << endl;
-volVectorField U
-(
-    IOobject
-    (
-        "U",
-        runTime.timeName(),
-        mesh,
-        IOobject::MUST_READ,
-        IOobject::AUTO_WRITE
-    ),
-    mesh
-);
-
-Info<< "Calculating face flux field phi\n" << endl;
-
-surfaceScalarField phi
-(
-    IOobject
-    (
-        "phi",
-        runTime.timeName(),
-        mesh,
-        IOobject::READ_IF_PRESENT,
-        IOobject::AUTO_WRITE
-    ),
-    linearInterpolate(rho)*linearInterpolate(U) & mesh.Sf()
-);
-
-Info<< "Calculating face flux field phiByRho\n" << endl;
-
-surfaceScalarField phiByRho
-(
-    IOobject
-    (
-        "phiByRho",
-        runTime.timeName(),
-        mesh,
-        IOobject::READ_IF_PRESENT,
-        IOobject::AUTO_WRITE
-    ),
-    phi/linearInterpolate(rho)
-);
-
-Info<< "Creating turbulence model\n" << endl;
-autoPtr<compressible::turbulenceModel> turbulence
-(
-    compressible::turbulenceModel::New
-    (
-        rho,
-        U,
-        phi,
-        thermo
-    )
-);
-
-mesh.setFluxRequired(p.name());
-
-Info<< "Creating field dpdt\n" << endl;
-volScalarField dpdt
-(
-    IOobject
-    (
-        "dpdt",
-        runTime.timeName(),
-        mesh
-    ),
-    mesh,
-    dimensionedScalar("dpdt", p.dimensions()/dimTime, 0)
-);
-
-#include "createMRF.H"
-
-Info<< "Creating compressibility field psi\n" << endl;
-volScalarField psi("psi", 1.0/((thermo.Cp() - thermo.Cv())*T));
-psi.oldTime() = 1.0/((thermo.Cp() - thermo.Cv())*T.oldTime());
-psi.oldTime().oldTime() = 1.0/((thermo.Cp()-thermo.Cv())*T.oldTime().oldTime());

applications/solvers/compressible/rhoPimpleAdiabaticFoam/pEqn.H

@@ -1,95 +0,0 @@
-{
-    volScalarField rAU(1.0/UEqn.A());
-    volVectorField HbyA("HbyA", U);
-    HbyA = rAU*UEqn.H();
-
-
-    // Define coefficients and pseudo-velocities for RCM interpolation
-    // M[U] = AU - H = -grad(p)
-    // U = H/A - 1/A grad(p)
-    // H/A = U + 1/A grad(p)
-    surfaceScalarField rhorAUf
-    (
-        "rhorAUf",
-        fvc::interpolate(rho)/fvc::interpolate(UEqn.A())
-    );
-
-    surfaceVectorField rhoHbyAf
-    (
-        "rhoHbyAf",
-        fvc::interpolate(rho)*fvc::interpolate(U)
-      + rhorAUf*fvc::interpolate(fvc::grad(p))
-    );
-
-    #include "resetBoundaries.H"
-
-    if (pimple.nCorrPISO() <= 1)
-    {
-        tUEqn.clear();
-    }
-
-    if (pimple.transonic())
-    {
-         FatalError
-             << "\nTransonic option not available for " << args.executable()
-             << exit(FatalError);
-    }
-    else
-    {
-        // Rhie & Chow interpolation (part 1)
-        surfaceScalarField phiHbyA
-        (
-            "phiHbyA",
-            (
-                (rhoHbyAf & mesh.Sf())
-              + rhorAUf*fvc::interpolate(rho)*fvc::ddtCorr(U, phiByRho)
-              + fvc::interpolate(rho)
-              * fvc::alphaCorr(U, phiByRho, pimple.finalInnerIter())
-            )
-        );
-
-        MRF.makeRelative(fvc::interpolate(rho), phiHbyA);
-
-        // Non-orthogonal pressure corrector loop
-        while (pimple.correctNonOrthogonal())
-        {
-            // Pressure corrector
-            fvScalarMatrix pEqn
-            (
-                fvm::ddt(psi, p)
-              + fvc::div(phiHbyA)
-              - fvm::laplacian(rhorAUf, p)
-              ==
-                fvOptions(psi, p, rho.name())
-            );
-
-            pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter())));
-
-            // Rhie & Chow interpolation (part 2)
-            if (pimple.finalNonOrthogonalIter())
-            {
-                phi = phiHbyA + pEqn.flux();
-            }
-        }
-    }
-
-    phiByRho = phi/fvc::interpolate(rho);
-
-    #include "rhoEqn.H"
-    #include "compressibleContinuityErrs.H"
-
-    // Explicitly relax pressure for momentum corrector
-    p.relax();
-
-    U = HbyA - rAU*fvc::grad(p);
-    U.correctBoundaryConditions();
-    fvOptions.correct(U);
-}
-
-rho = thermo.rho();
-
-if (thermo.dpdt())
-{
-    dpdt = fvc::ddt(p);
-}
-

applications/solvers/compressible/rhoPimpleAdiabaticFoam/resetBoundaries.H

@@ -1,25 +0,0 @@
-{
-    // Keep standard formulation on domain boundaries to ensure compatibility
-    // with existing boundary conditions
-    const Foam::FieldField<Foam::fvsPatchField, scalar> rhorAUf_orig
-    (
-        fvc::interpolate(rho.boundaryField()*rAU.boundaryField())
-    );
-
-    const Foam::FieldField<Foam::fvsPatchField, vector> rhoHbyA_orig
-    (
-        fvc::interpolate(rho.boundaryField()*HbyA.boundaryField())
-    );
-
-    surfaceScalarField::Boundary& rhorAUfbf = rhorAUf.boundaryFieldRef();
-    surfaceVectorField::Boundary& rhoHbyAfbf = rhoHbyAf.boundaryFieldRef();
-
-    forAll(U.boundaryField(), patchi)
-    {
-        if (!U.boundaryField()[patchi].coupled())
-        {
-            rhorAUfbf[patchi] = rhorAUf_orig[patchi];
-            rhoHbyAfbf[patchi] = rhoHbyA_orig[patchi];
-        }
-    }
-}

applications/solvers/compressible/rhoPimpleAdiabaticFoam/rhoPimpleAdiabaticFoam.C

@@ -1,129 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2017 OpenCFD Ltd.
-     \\/     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
-    rhoPimpleAdiabaticFoam
-
-Description
-    Transient solver for laminar or turbulent flow of weakly compressible
-    fluids for low Mach number aeroacoustic applications.
-
-    Uses the flexible PIMPLE (PISO-SIMPLE) solution for time-resolved and
-    pseudo-transient simulations. The RCM interpolation is used as in
-
-    \verbatim
-        Knacke, T. (2013).
-        Potential effects of Rhie & Chow type interpolations in airframe
-        noise simulations. In: Schram, C., Dénos, R., Lecomte E. (ed):
-        Accurate and efficient aeroacoustic prediction approaches for
-        airframe noise, VKI LS 2013-03.
-    \endverbatim
-
-
-    Contact:      info@cfd-berlin.com
-
-\*---------------------------------------------------------------------------*/
-
-#include "fvCFD.H"
-#include "fluidThermo.H"
-#include "turbulentFluidThermoModel.H"
-#include "bound.H"
-#include "pimpleControl.H"
-#include "fvOptions.H"
-
-#include "ddtScheme.H"
-#include "fvcCorrectAlpha.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-int main(int argc, char *argv[])
-{
-    #include "postProcess.H"
-
-    #include "setRootCase.H"
-    #include "createTime.H"
-    #include "createMesh.H"
-    #include "createControl.H"
-    #include "createTimeControls.H"
-
-    #include "createFields.H"
-    #include "createFvOptions.H"
-    #include "initContinuityErrs.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;
-
-        if (pimple.nCorrPIMPLE() <= 1)
-        {
-            #include "rhoEqn.H"
-        }
-
-        // --- Pressure-velocity PIMPLE corrector loop
-        while (pimple.loop())
-        {
-            U.storePrevIter();
-            rho.storePrevIter();
-            phi.storePrevIter();
-            phiByRho.storePrevIter();
-
-            #include "UEqn.H"
-
-            // --- Pressure corrector loop
-            while (pimple.correct())
-            {
-                #include "pEqn.H"
-            }
-
-            #include "EEqn.H"
-
-            if (pimple.turbCorr())
-            {
-                turbulence->correct();
-            }
-        }
-
-        runTime.write();
-
-        Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
-            << "  ClockTime = " << runTime.elapsedClockTime() << " s"
-            << nl << endl;
-    }
-
-    Info<< "End\n" << endl;
-
-    return 0;
-}
-
-
-// ************************************************************************* //

applications/solvers/compressible/rhoPimpleFoam/createFields.H

@@ -2,11 +2,11 @@
 
 Info<< "Reading thermophysical properties\n" << endl;
 
-autoPtr<fluidThermo> pThermo
+autoPtr<psiThermo> pThermo
 (
-    fluidThermo::New(mesh)
+    psiThermo::New(mesh)
 );
-fluidThermo& thermo = pThermo();
+psiThermo& thermo = pThermo();
 thermo.validate(args.executable(), "h", "e");
 
 volScalarField& p = thermo.p();
@@ -40,29 +40,6 @@ volVectorField U
 
 #include "compressibleCreatePhi.H"
 
-pressureControl pressureControl(p, rho, pimple.dict(), false);
-
-mesh.setFluxRequired(p.name());
-
-Info<< "Creating turbulence model\n" << endl;
-autoPtr<compressible::turbulenceModel> turbulence
-(
-    compressible::turbulenceModel::New
-    (
-        rho,
-        U,
-        phi,
-        thermo
-    )
-);
-
-#include "createDpdt.H"
-
-#include "createK.H"
-
-#include "createMRF.H"
-
-
 dimensionedScalar rhoMax
 (
     dimensionedScalar::lookupOrDefault
@@ -84,3 +61,35 @@ dimensionedScalar rhoMin
         0
     )
 );
+
+Info<< "Creating turbulence model\n" << endl;
+autoPtr<compressible::turbulenceModel> turbulence
+(
+    compressible::turbulenceModel::New
+    (
+        rho,
+        U,
+        phi,
+        thermo
+    )
+);
+
+mesh.setFluxRequired(p.name());
+
+Info<< "Creating field dpdt\n" << endl;
+volScalarField dpdt
+(
+    IOobject
+    (
+        "dpdt",
+        runTime.timeName(),
+        mesh
+    ),
+    mesh,
+    dimensionedScalar("dpdt", p.dimensions()/dimTime, 0)
+);
+
+Info<< "Creating field kinetic energy K\n" << endl;
+volScalarField K("K", 0.5*magSqr(U));
+
+#include "createMRF.H"

applications/solvers/compressible/rhoPimpleFoam/overRhoPimpleDyMFoam/Make/files

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

applications/solvers/compressible/rhoPimpleFoam/overRhoPimpleDyMFoam/Make/options

@@ -1,25 +0,0 @@
-EXE_INC = \
-    -I.. \
-    -I$(LIB_SRC)/transportModels/compressible/lnInclude \
-    -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
-    -I$(LIB_SRC)/TurbulenceModels/turbulenceModels/lnInclude \
-    -I$(LIB_SRC)/TurbulenceModels/compressible/lnInclude \
-    -I$(LIB_SRC)/finiteVolume/cfdTools \
-    -I$(LIB_SRC)/finiteVolume/lnInclude \
-    -I$(LIB_SRC)/meshTools/lnInclude \
-    -I$(LIB_SRC)/sampling/lnInclude \
-    -I$(LIB_SRC)/dynamicFvMesh/lnInclude \
-    -I$(LIB_SRC)/dynamicMesh/lnInclude \
-    -I$(LIB_SRC)/overset/lnInclude \
-    -I$(LIB_SRC)/meshTools/lnInclude \
-
-EXE_LIBS = \
-    -lcompressibleTransportModels \
-    -lfluidThermophysicalModels \
-    -lspecie \
-    -lturbulenceModels \
-    -lcompressibleTurbulenceModels \
-    -loverset \
-    -lfvOptions \
-    -ldynamicFvMesh \
-    -ltopoChangerFvMesh

applications/solvers/compressible/rhoPimpleFoam/overRhoPimpleDyMFoam/correctPhi.H

@@ -1,11 +0,0 @@
-CorrectPhi
-(
-    U,
-    phi,
-    p,
-    rho,
-    psi,
-    dimensionedScalar("rAUf", dimTime, 1),
-    divrhoU,
-    pimple
-);

applications/solvers/compressible/rhoPimpleFoam/overRhoPimpleDyMFoam/createControls.H

@@ -1,11 +0,0 @@
-#include "createTimeControls.H"
-
-bool correctPhi
-(
-    pimple.dict().lookupOrDefault<Switch>("correctPhi", true)
-);
-
-bool checkMeshCourantNo
-(
-    pimple.dict().lookupOrDefault<Switch>("checkMeshCourantNo", false)
-);

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

@@ -1,123 +0,0 @@
-rho = thermo.rho();
-rho = max(rho, rhoMin);
-rho = min(rho, rhoMax);
-rho.relax();
-
-surfaceScalarField faceMask(localMin<scalar>(mesh).interpolate(cellMask));
-
-volScalarField rAU(1.0/UEqn.A());
-surfaceScalarField rhorAUf("rhorAUf", faceMask*fvc::interpolate(rho*rAU));
-volVectorField HbyA("HbyA", constrainHbyA(rAU*UEqn.H(), U, p));
-//mesh.interpolate(HbyA);
-
-if (pimple.nCorrPISO() <= 1)
-{
-    tUEqn.clear();
-}
-
-if (pimple.transonic())
-{
-    surfaceScalarField phid
-    (
-        "phid",
-        fvc::interpolate(psi)
-       *(
-            fvc::flux(HbyA)
-          + rhorAUf*fvc::ddtCorr(rho, U, rhoUf)/fvc::interpolate(rho)
-        )
-    );
-
-    fvc::makeRelative(phid, psi, U);
-    MRF.makeRelative(fvc::interpolate(psi), phid);
-
-    while (pimple.correctNonOrthogonal())
-    {
-        fvScalarMatrix pEqn
-        (
-            fvm::ddt(psi, p)
-          + fvm::div(phid, p)
-          - fvm::laplacian(rhorAUf, p)
-         ==
-            fvOptions(psi, p, rho.name())
-        );
-
-        pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter())));
-
-        if (pimple.finalNonOrthogonalIter())
-        {
-            phi == 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
-        fvScalarMatrix pEqn
-        (
-            fvm::ddt(psi, p)
-          + fvc::div(phiHbyA)
-          - fvm::laplacian(rhorAUf, p)
-         ==
-            fvOptions(psi, p, rho.name())
-        );
-
-        pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter())));
-
-        if (pimple.finalNonOrthogonalIter())
-        {
-            phi = phiHbyA + pEqn.flux();
-        }
-    }
-}
-
-#include "rhoEqn.H"
-#include "compressibleContinuityErrs.H"
-
-// 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;
-
-volVectorField gradP(fvc::grad(p));
-//mesh.interpolate(gradP);
-U = HbyA - rAU*cellMask*gradP;
-U.correctBoundaryConditions();
-fvOptions.correct(U);
-K = 0.5*magSqr(U);
-
-{
-    rhoUf = fvc::interpolate(rho*U);
-    surfaceVectorField n(mesh.Sf()/mesh.magSf());
-    rhoUf += n*(fvc::absolute(phi, rho, U)/mesh.magSf() - (n & rhoUf));
-}
-
-if (thermo.dpdt())
-{
-    dpdt = fvc::ddt(p);
-
-    if (mesh.moving())
-    {
-        dpdt -= fvc::div(fvc::meshPhi(rho, U), p);
-    }
-}

applications/solvers/compressible/rhoPimpleFoam/overRhoPimpleDyMFoam/readControls.H

@@ -1,6 +0,0 @@
-#include "readTimeControls.H"
-
-correctPhi = pimple.dict().lookupOrDefault<Switch>("correctPhi", true);
-
-checkMeshCourantNo =
-    pimple.dict().lookupOrDefault<Switch>("checkMeshCourantNo", false);

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

@@ -1,171 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2015 OpenFOAM Foundation
-     \\/     M anipulation  | Copyright (C) 2016-2017 OpenCFD Ltd.
--------------------------------------------------------------------------------
-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
-    rhoPimpleFoam
-
-Group
-    grpCompressibleSolvers grpMovingMeshSolvers
-
-Description
-    Transient solver for laminar or turbulent flow of compressible fluids
-    for HVAC and similar applications.
-
-    Uses the flexible PIMPLE (PISO-SIMPLE) solution for time-resolved and
-    pseudo-transient simulations.
-
-\*---------------------------------------------------------------------------*/
-
-#include "fvCFD.H"
-#include "dynamicFvMesh.H"
-#include "psiThermo.H"
-#include "turbulentFluidThermoModel.H"
-#include "bound.H"
-#include "pimpleControl.H"
-#include "CorrectPhi.H"
-#include "fvOptions.H"
-#include "localEulerDdtScheme.H"
-#include "fvcSmooth.H"
-#include "cellCellStencilObject.H"
-#include "localMin.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-int main(int argc, char *argv[])
-{
-    #include "setRootCase.H"
-    #include "createTime.H"
-    #include "createDynamicFvMesh.H"
-    #include "createControl.H"
-    #include "createRDeltaT.H"
-    #include "initContinuityErrs.H"
-    #include "createFields.H"
-    #include "createMRF.H"
-    #include "createFvOptions.H"
-    #include "createRhoUf.H"
-    #include "createControls.H"
-
-    turbulence->validate();
-
-    if (!LTS)
-    {
-        #include "compressibleCourantNo.H"
-        #include "setInitialDeltaT.H"
-    }
-
-    // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-    Info<< "\nStarting time loop\n" << endl;
-
-    while (runTime.run())
-    {
-        #include "readControls.H"
-
-        {
-            // Store divrhoU from the previous mesh so that it can be mapped
-            // and used in correctPhi to ensure the corrected phi has the
-            // same divergence
-            volScalarField divrhoU
-            (
-                "divrhoU",
-                fvc::div(fvc::absolute(phi, rho, U))
-            );
-
-            if (LTS)
-            {
-                #include "setRDeltaT.H"
-            }
-            else
-            {
-                #include "compressibleCourantNo.H"
-                #include "setDeltaT.H"
-            }
-
-            runTime++;
-
-            Info<< "Time = " << runTime.timeName() << nl << endl;
-
-            // Store momentum to set rhoUf for introduced faces.
-            volVectorField rhoU("rhoU", rho*U);
-
-            // Do any mesh changes
-            mesh.update();
-
-            if (mesh.changing())
-            {
-                #include "setCellMask.H"
-            }
-
-            if (mesh.changing() && correctPhi)
-            {
-                // Calculate absolute flux from the mapped surface velocity
-                phi = mesh.Sf() & rhoUf;
-
-                #include "correctPhi.H"
-
-                // Make the fluxes relative to the mesh-motion
-                fvc::makeRelative(phi, rho, U);
-            }
-        }
-
-        if (mesh.changing() && checkMeshCourantNo)
-        {
-            #include "meshCourantNo.H"
-        }
-
-        #include "rhoEqn.H"
-        Info<< "rhoEqn max/min : " << max(rho).value()
-            << " " << min(rho).value() << endl;
-
-        // --- Pressure-velocity PIMPLE corrector loop
-        while (pimple.loop())
-        {
-            #include "UEqn.H"
-            #include "EEqn.H"
-
-            // --- Pressure corrector loop
-            while (pimple.correct())
-            {
-                #include "pEqn.H"
-            }
-
-            if (pimple.turbCorr())
-            {
-                turbulence->correct();
-            }
-        }
-
-        runTime.write();
-
-        Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
-            << "  ClockTime = " << runTime.elapsedClockTime() << " s"
-            << nl << endl;
-    }
-
-    Info<< "End\n" << endl;
-
-    return 0;
-}
-
-
-// ************************************************************************* //

applications/solvers/compressible/rhoPimpleFoam/pEqn.H

@@ -1,11 +1,7 @@
-if (!pimple.SIMPLErho())
-{
-    rho = thermo.rho();
-}
-
-// Thermodynamic density needs to be updated by psi*d(p) after the
-// pressure solution
-const volScalarField psip0(psi*p);
+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));
@@ -16,64 +12,65 @@ if (pimple.nCorrPISO() <= 1)
     tUEqn.clear();
 }
 
-surfaceScalarField phiHbyA
-(
-    "phiHbyA",
-    fvc::interpolate(rho)*fvc::flux(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(rho))*phiHbyA
+        fvc::interpolate(psi)
+       *(
+            fvc::flux(HbyA)
+          + rhorAUf*fvc::ddtCorr(rho, U, phi)/fvc::interpolate(rho)
+        )
     );
 
-    phiHbyA -= fvc::interpolate(psi*p)*phiHbyA/fvc::interpolate(rho);
-
-    fvScalarMatrix pDDtEqn
-    (
-        fvc::ddt(rho) + psi*correction(fvm::ddt(p))
-      + fvc::div(phiHbyA) + fvm::div(phid, p)
-     ==
-        fvOptions(psi, p, rho.name())
-    );
+    MRF.makeRelative(fvc::interpolate(psi), phid);
 
     while (pimple.correctNonOrthogonal())
     {
-        fvScalarMatrix pEqn(pDDtEqn - fvm::laplacian(rhorAUf, p));
-
-        // Relax the pressure equation to ensure diagonal-dominance
-        pEqn.relax();
+        fvScalarMatrix pEqn
+        (
+            fvm::ddt(psi, p)
+          + fvm::div(phid, p)
+          - fvm::laplacian(rhorAUf, p)
+         ==
+            fvOptions(psi, p, rho.name())
+        );
 
         pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter())));
 
         if (pimple.finalNonOrthogonalIter())
         {
-            phi = phiHbyA + pEqn.flux();
+            phi == pEqn.flux();
         }
     }
 }
 else
 {
-    fvScalarMatrix pDDtEqn
+    surfaceScalarField phiHbyA
     (
-        fvc::ddt(rho) + psi*correction(fvm::ddt(p))
-      + fvc::div(phiHbyA)
-     ==
-        fvOptions(psi, p, rho.name())
+        "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(pDDtEqn - fvm::laplacian(rhorAUf, p));
+        fvScalarMatrix pEqn
+        (
+            fvm::ddt(psi, p)
+          + fvc::div(phiHbyA)
+          - fvm::laplacian(rhorAUf, p)
+         ==
+            fvOptions(psi, p, rho.name())
+        );
 
         pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter())));
 
@@ -90,19 +87,19 @@ 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);
 
-if (pressureControl.limit(p))
-{
-    p.correctBoundaryConditions();
-}
-
-thermo.correctRho(psi*p - psip0, rhoMin, rhoMax) ;
-rho = thermo.rho();
-
 if (thermo.dpdt())
 {
     dpdt = fvc::ddt(p);

applications/solvers/compressible/rhoPimpleFoam/pcEqn.H

@@ -1,11 +1,7 @@
-if (!pimple.SIMPLErho())
-{
-    rho = thermo.rho();
-}
-
-// Thermodynamic density needs to be updated by psi*d(p) after the
-// pressure solution
-const volScalarField psip0(psi*p);
+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()));
@@ -16,75 +12,82 @@ if (pimple.nCorrPISO() <= 1)
     tUEqn.clear();
 }
 
-surfaceScalarField phiHbyA
-(
-    "phiHbyA",
-    (
-        fvc::interpolate(rho)*fvc::flux(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(rho))*phiHbyA
+        fvc::interpolate(psi)
+       *(
+            fvc::flux(HbyA)
+          + fvc::interpolate(rho*rAU)*fvc::ddtCorr(rho, U, phi)
+           /fvc::interpolate(rho)
+        )
     );
 
-    phiHbyA +=
+    MRF.makeRelative(fvc::interpolate(psi), phid);
+
+    surfaceScalarField phic
+    (
+        "phic",
         fvc::interpolate(rho*(rAtU - rAU))*fvc::snGrad(p)*mesh.magSf()
-      - fvc::interpolate(psi*p)*phiHbyA/fvc::interpolate(rho);
+    );
 
     HbyA -= (rAU - rAtU)*fvc::grad(p);
 
-    fvScalarMatrix pDDtEqn
-    (
-        fvc::ddt(rho) + psi*correction(fvm::ddt(p))
-      + fvc::div(phiHbyA) + fvm::div(phid, p)
-     ==
-        fvOptions(psi, p, rho.name())
-    );
+    volScalarField rhorAtU("rhorAtU", rho*rAtU);
 
     while (pimple.correctNonOrthogonal())
     {
-        fvScalarMatrix pEqn(pDDtEqn - fvm::laplacian(rhorAtU, p));
-
-        // Relax the pressure equation to ensure diagonal-dominance
-        pEqn.relax();
+        fvScalarMatrix pEqn
+        (
+            fvm::ddt(psi, p)
+          + fvm::div(phid, p)
+          + fvc::div(phic)
+          - fvm::laplacian(rhorAtU, p)
+         ==
+            fvOptions(psi, p, rho.name())
+        );
 
         pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter())));
 
         if (pimple.finalNonOrthogonalIter())
         {
-            phi = phiHbyA + pEqn.flux();
+            phi == phic + 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);
 
-    fvScalarMatrix pDDtEqn
-    (
-        fvc::ddt(rho) + psi*correction(fvm::ddt(p))
-      + fvc::div(phiHbyA)
-     ==
-        fvOptions(psi, p, rho.name())
-    );
+    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(pDDtEqn - fvm::laplacian(rhorAtU, p));
+        fvScalarMatrix pEqn
+        (
+            fvm::ddt(psi, p)
+          + fvc::div(phiHbyA)
+          - fvm::laplacian(rhorAtU, p)
+         ==
+            fvOptions(psi, p, rho.name())
+        );
 
         pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter())));
 
@@ -106,15 +109,19 @@ U.correctBoundaryConditions();
 fvOptions.correct(U);
 K = 0.5*magSqr(U);
 
-if (pressureControl.limit(p))
-{
-    p.correctBoundaryConditions();
-}
-
-thermo.correctRho(psi*p - psip0, rhoMin, rhoMax) ;
-rho = thermo.rho();
-
 if (thermo.dpdt())
 {
     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;

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

@@ -1,11 +1,7 @@
-if (!pimple.SIMPLErho())
-{
-    rho = thermo.rho();
-}
-
-// Thermodynamic density needs to be updated by psi*d(p) after the
-// pressure solution
-const volScalarField psip0(psi*p);
+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));
@@ -16,60 +12,61 @@ 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(rho))*phiHbyA
+        fvc::interpolate(psi)
+       *(
+            fvc::flux(HbyA)
+          + rhorAUf*fvc::ddtCorr(rho, U, rhoUf)/fvc::interpolate(rho)
+        )
     );
 
-    phiHbyA -= fvc::interpolate(psi*p)*phiHbyA/fvc::interpolate(rho);
+    fvc::makeRelative(phid, psi, U);
+    MRF.makeRelative(fvc::interpolate(psi), phid);
 
     while (pimple.correctNonOrthogonal())
     {
         fvScalarMatrix pEqn
         (
-            fvc::ddt(rho) + psi*correction(fvm::ddt(p))
-          + fvc::div(phiHbyA)
+            fvm::ddt(psi, p)
           + 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 = phiHbyA + pEqn.flux();
+            phi == 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
         fvScalarMatrix pEqn
         (
-            fvc::ddt(rho) + psi*correction(fvm::ddt(p))
+            fvm::ddt(psi, p)
           + fvc::div(phiHbyA)
           - fvm::laplacian(rhorAUf, p)
          ==
@@ -91,19 +88,19 @@ 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);
 
-if (pressureControl.limit(p))
-{
-    p.correctBoundaryConditions();
-}
-
-thermo.correctRho(psi*p - psip0, rhoMin, rhoMax);
-rho = thermo.rho();
-
 {
     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-2017 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -22,7 +22,7 @@ License
     along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
 
 Application
-    rhoPimpleDyMFoam
+    rhoPimpleFoam
 
 Group
     grpCompressibleSolvers grpMovingMeshSolvers
@@ -38,11 +38,10 @@ Description
 
 #include "fvCFD.H"
 #include "dynamicFvMesh.H"
-#include "fluidThermo.H"
+#include "psiThermo.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-2017 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -37,11 +37,10 @@ Description
 \*---------------------------------------------------------------------------*/
 
 #include "fvCFD.H"
-#include "fluidThermo.H"
+#include "psiThermo.H"
 #include "turbulentFluidThermoModel.H"
 #include "bound.H"
 #include "pimpleControl.H"
-#include "pressureControl.H"
 #include "fvOptions.H"
 #include "localEulerDdtScheme.H"
 #include "fvcSmooth.H"
@@ -122,8 +121,6 @@ int main(int argc, char *argv[])
             }
         }
 
-        rho = thermo.rho();
-
         runTime.write();
 
         Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"

applications/solvers/compressible/rhoSimpleFoam/createFields.H

@@ -1,14 +1,11 @@
 Info<< "Reading thermophysical properties\n" << endl;
-
-autoPtr<fluidThermo> pThermo
+autoPtr<rhoThermo> pThermo
 (
-    fluidThermo::New(mesh)
+    rhoThermo::New(mesh)
 );
-fluidThermo& thermo = pThermo();
+rhoThermo& thermo = pThermo();
 thermo.validate(args.executable(), "h", "e");
 
-volScalarField& p = thermo.p();
-
 volScalarField rho
 (
     IOobject
@@ -22,6 +19,8 @@ volScalarField rho
     thermo.rho()
 );
 
+volScalarField& p = thermo.p();
+
 Info<< "Reading field U\n" << endl;
 volVectorField U
 (
@@ -38,10 +37,35 @@ 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,110 +1,109 @@
-volScalarField rAU(1.0/UEqn.A());
-surfaceScalarField rhorAUf("rhorAUf", fvc::interpolate(rho*rAU));
-volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p));
-tUEqn.clear();
-
-bool closedVolume = false;
-
-surfaceScalarField phiHbyA("phiHbyA", fvc::interpolate(rho)*fvc::flux(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(rho))*phiHbyA
-    );
+    //const volScalarField& psi = thermo.psi();
 
-    phiHbyA -= fvc::interpolate(psi*p)*phiHbyA/fvc::interpolate(rho);
+    volScalarField rAU(1.0/UEqn.A());
+    surfaceScalarField rhorAUf("rhorAUf", fvc::interpolate(rho*rAU));
+    volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p));
+    tUEqn.clear();
 
-    while (simple.correctNonOrthogonal())
+    bool closedVolume = false;
+
+    if (simple.transonic())
     {
-        fvScalarMatrix pEqn
+        surfaceScalarField phid
         (
-            fvc::div(phiHbyA)
-          + fvm::div(phid, p)
-          - fvm::laplacian(rhorAUf, p)
-          ==
-            fvOptions(psi, p, rho.name())
+            "phid",
+            fvc::interpolate(psi)
+           *fvc::flux(HbyA)
         );
 
-        // Relax the pressure equation to ensure diagonal-dominance
-        pEqn.relax();
+        MRF.makeRelative(fvc::interpolate(psi), phid);
 
-        pEqn.setReference
-        (
-            pressureControl.refCell(),
-            pressureControl.refValue()
-        );
-
-        pEqn.solve();
-
-        if (simple.finalNonOrthogonalIter())
+        while (simple.correctNonOrthogonal())
         {
-            phi = phiHbyA + pEqn.flux();
+            fvScalarMatrix pEqn
+            (
+                fvm::div(phid, p)
+              - fvm::laplacian(rhorAUf, p)
+              ==
+                fvOptions(psi, p, rho.name())
+            );
+
+            // Relax the pressure equation to ensure diagonal-dominance
+            pEqn.relax();
+
+            pEqn.setReference(pRefCell, pRefValue);
+
+            pEqn.solve();
+
+            if (simple.finalNonOrthogonalIter())
+            {
+                phi == pEqn.flux();
+            }
         }
     }
-}
-else
-{
-    closedVolume = adjustPhi(phiHbyA, U, p);
-
-    while (simple.correctNonOrthogonal())
+    else
     {
-        fvScalarMatrix pEqn
-        (
-            fvc::div(phiHbyA)
-          - fvm::laplacian(rhorAUf, p)
-          ==
-            fvOptions(psi, p, rho.name())
-        );
+        surfaceScalarField phiHbyA("phiHbyA", fvc::flux(rho*HbyA));
+        MRF.makeRelative(fvc::interpolate(rho), phiHbyA);
 
-        pEqn.setReference
-        (
-            pressureControl.refCell(),
-            pressureControl.refValue()
-        );
+        closedVolume = adjustPhi(phiHbyA, U, p);
 
-        pEqn.solve();
+        // Update the pressure BCs to ensure flux consistency
+        constrainPressure(p, rho, U, phiHbyA, rhorAUf, MRF);
 
-        if (simple.finalNonOrthogonalIter())
+        while (simple.correctNonOrthogonal())
         {
-            phi = phiHbyA + pEqn.flux();
+            fvScalarMatrix pEqn
+            (
+                fvc::div(phiHbyA)
+              - fvm::laplacian(rhorAUf, p)
+              ==
+                fvOptions(psi, p, rho.name())
+            );
+
+            pEqn.setReference(pRefCell, pRefValue);
+
+            pEqn.solve();
+
+            if (simple.finalNonOrthogonalIter())
+            {
+                phi = phiHbyA + pEqn.flux();
+            }
         }
     }
-}
-
-#include "incompressible/continuityErrs.H"
-
-// Explicitly relax pressure for momentum corrector
-p.relax();
-
-U = HbyA - rAU*fvc::grad(p);
-U.correctBoundaryConditions();
-fvOptions.correct(U);
-
-bool pLimited = pressureControl.limit(p);
-
-// For closed-volume cases adjust the pressure and density levels
-// to obey overall mass continuity
-if (closedVolume)
-{
-    p += (initialMass - fvc::domainIntegrate(psi*p))
-        /fvc::domainIntegrate(psi);
-}
-
-if (pLimited || closedVolume)
-{
-    p.correctBoundaryConditions();
-}
-
-rho = thermo.rho();
-
-if (!simple.transonic())
-{
-    rho.relax();
+
+
+    #include "incompressible/continuityErrs.H"
+
+    // Explicitly relax pressure for momentum corrector
+    p.relax();
+
+    U = HbyA - rAU*fvc::grad(p);
+    U.correctBoundaryConditions();
+    fvOptions.correct(U);
+
+    // For closed-volume cases adjust the pressure and density levels
+    // to obey overall mass continuity
+    if (closedVolume)
+    {
+        p += (initialMass - fvc::domainIntegrate(psi*p))
+            /fvc::domainIntegrate(psi);
+    }
+
+    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/pcEqn.H

@@ -1,5 +1,3 @@
-rho = thermo.rho();
-
 volScalarField rAU(1.0/UEqn.A());
 volScalarField rAtU(1.0/(1.0/rAU - UEqn.H1()));
 volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p));
@@ -7,34 +5,33 @@ tUEqn.clear();
 
 bool closedVolume = false;
 
-surfaceScalarField phiHbyA("phiHbyA", fvc::interpolate(rho)*fvc::flux(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(rho))*phiHbyA
+        fvc::interpolate(psi)
+       *fvc::flux(HbyA)
     );
 
-    phiHbyA +=
+    MRF.makeRelative(fvc::interpolate(psi), phid);
+
+    surfaceScalarField phic
+    (
+        "phic",
         fvc::interpolate(rho*(rAtU - rAU))*fvc::snGrad(p)*mesh.magSf()
-      - fvc::interpolate(psi*p)*phiHbyA/fvc::interpolate(rho);
+    );
 
     HbyA -= (rAU - rAtU)*fvc::grad(p);
 
+    volScalarField rhorAtU("rhorAtU", rho*rAtU);
+
     while (simple.correctNonOrthogonal())
     {
         fvScalarMatrix pEqn
         (
-            fvc::div(phiHbyA)
-          + fvm::div(phid, p)
+            fvm::div(phid, p)
+          + fvc::div(phic)
           - fvm::laplacian(rhorAtU, p)
          ==
             fvOptions(psi, p, rho.name())
@@ -43,27 +40,31 @@ if (simple.transonic())
         // Relax the pressure equation to maintain diagonal dominance
         pEqn.relax();
 
-        pEqn.setReference
-        (
-            pressureControl.refCell(),
-            pressureControl.refValue()
-        );
+        pEqn.setReference(pRefCell, pRefValue);
 
         pEqn.solve();
 
         if (simple.finalNonOrthogonalIter())
         {
-            phi = phiHbyA + pEqn.flux();
+            phi == phic + 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
@@ -74,11 +75,7 @@ else
             fvOptions(psi, p, rho.name())
         );
 
-        pEqn.setReference
-        (
-            pressureControl.refCell(),
-            pressureControl.refValue()
-        );
+        pEqn.setReference(pRefCell, pRefValue);
 
         pEqn.solve();
 
@@ -100,8 +97,6 @@ U = HbyA - rAtU*fvc::grad(p);
 U.correctBoundaryConditions();
 fvOptions.correct(U);
 
-bool pLimited = pressureControl.limit(p);
-
 // For closed-volume cases adjust the pressure and density levels
 // to obey overall mass continuity
 if (closedVolume)
@@ -110,14 +105,16 @@ if (closedVolume)
         /fvc::domainIntegrate(psi);
 }
 
-if (pLimited || closedVolume)
-{
-    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

@@ -0,0 +1,26 @@
+{
+    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,15 +1,12 @@
 Info<< "Reading thermophysical properties\n" << endl;
 
-autoPtr<psiThermo> pThermo
+autoPtr<rhoThermo> pThermo
 (
-    psiThermo::New(mesh)
+    rhoThermo::New(mesh)
 );
-psiThermo& thermo = pThermo();
+rhoThermo& thermo = pThermo();
 thermo.validate(args.executable(), "h", "e");
 
-volScalarField& p = thermo.p();
-const volScalarField& psi = thermo.psi();
-
 volScalarField rho
 (
     IOobject
@@ -23,6 +20,8 @@ volScalarField rho
     thermo.rho()
 );
 
+volScalarField& p = thermo.p();
+
 Info<< "Reading field U\n" << endl;
 volVectorField U
 (
@@ -39,12 +38,19 @@ volVectorField U
 
 #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",
-        pimple.dict(),
+        simple.dict(),
         dimDensity,
         GREAT
     )
@@ -55,7 +61,7 @@ dimensionedScalar rhoMin
     dimensionedScalar::lookupOrDefault
     (
         "rhoMin",
-        pimple.dict(),
+        simple.dict(),
         dimDensity,
         0
     )
@@ -73,45 +79,6 @@ autoPtr<compressible::turbulenceModel> turbulence
     )
 );
 
-mesh.setFluxRequired(p.name());
+dimensionedScalar initialMass = fvc::domainIntegrate(rho);
 
-Info<< "Creating field dpdt\n" << endl;
-volScalarField dpdt
-(
-    IOobject
-    (
-        "dpdt",
-        runTime.timeName(),
-        mesh
-    ),
-    mesh,
-    dimensionedScalar("dpdt", p.dimensions()/dimTime, 0)
-);
-
-Info<< "Creating field kinetic energy K\n" << endl;
-volScalarField K("K", 0.5*magSqr(U));
-
-
-//- Overset specific
-
-// Add solver-specific interpolations
-{
-    dictionary oversetDict;
-    oversetDict.add("U", true);
-    oversetDict.add("p", true);
-    oversetDict.add("HbyA", true);
-    oversetDict.add("grad(p)", true);
-
-    const_cast<dictionary&>
-    (
-        mesh.schemesDict()
-    ).add
-    (
-        "oversetInterpolationRequired",
-        oversetDict,
-        true
-    );
-}
-
-// Mask field for zeroing out contributions on hole cells
-#include "createCellMask.H"
+#include "createMRF.H"

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

@@ -48,11 +48,7 @@
 
         fvScalarMatrix& pEqn = tpEqn.ref();
 
-        pEqn.setReference
-        (
-            pressureControl.refCell(),
-            pressureControl.refValue()
-        );
+        pEqn.setReference(pRefCell, pRefValue);
 
         pEqn.solve();
 
@@ -79,8 +75,6 @@
     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)
@@ -90,5 +84,14 @@
     }
 
     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-2017 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -34,12 +34,11 @@ Description
 \*---------------------------------------------------------------------------*/
 
 #include "fvCFD.H"
-#include "fluidThermo.H"
+#include "rhoThermo.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-2017 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -33,10 +33,9 @@ Description
 \*---------------------------------------------------------------------------*/
 
 #include "fvCFD.H"
-#include "fluidThermo.H"
+#include "rhoThermo.H"
 #include "turbulentFluidThermoModel.H"
 #include "simpleControl.H"
-#include "pressureControl.H"
 #include "fvOptions.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

applications/solvers/compressible/sonicFoam/createFields.H

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

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

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

applications/solvers/electromagnetics/magneticFoam/magnet.H

@@ -138,7 +138,9 @@ public:
                 >> m.orientation_;
             is.readEnd("magnet");
 
-            is.check(FUNCTION_NAME);
+            // Check state of Istream
+            is.check("operator>>(Istream&, magnet&)");
+
             return is;
         }
 

applications/solvers/electromagnetics/mhdFoam/readBPISOControls.H

@@ -1,3 +1,4 @@
     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, he)
+      + radiation->Sh(thermo)
       + fvOptions(rho, he)
     );
 

applications/solvers/heatTransfer/buoyantPimpleFoam/createFields.H

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

applications/solvers/heatTransfer/buoyantPimpleFoam/pEqn.H

@@ -1,108 +1,115 @@
-
-dimensionedScalar compressibility = fvc::domainIntegrate(psi);
-bool compressible = (compressibility.value() > SMALL);
-
-rho = thermo.rho();
-
-// Thermodynamic density needs to be updated by psi*d(p) after the
-// pressure solution
-const volScalarField psip0(psi*p);
-
-volScalarField rAU(1.0/UEqn.A());
-surfaceScalarField rhorAUf("rhorAUf", fvc::interpolate(rho*rAU));
-volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p_rgh));
-
-surfaceScalarField phig(-rhorAUf*ghf*fvc::snGrad(rho)*mesh.magSf());
-
-surfaceScalarField phiHbyA
-(
-    "phiHbyA",
-    (
-        fvc::flux(rho*HbyA)
-      + rhorAUf*fvc::ddtCorr(rho, U, phi)
-    )
-  + phig
-);
-
-MRF.makeRelative(fvc::interpolate(rho), phiHbyA);
-
-// Update the pressure BCs to ensure flux consistency
-constrainPressure(p_rgh, rho, U, phiHbyA, rhorAUf, MRF);
-
-fvScalarMatrix p_rghDDtEqn
-(
-    fvc::ddt(rho) + psi*correction(fvm::ddt(p_rgh))
-  + fvc::div(phiHbyA)
-  ==
-    fvOptions(psi, p_rgh, rho.name())
-);
-
-while (pimple.correctNonOrthogonal())
 {
-    fvScalarMatrix p_rghEqn
+    bool closedVolume = p_rgh.needReference();
+
+    dimensionedScalar compressibility = fvc::domainIntegrate(psi);
+    bool compressible = (compressibility.value() > SMALL);
+
+    rho = thermo.rho();
+
+    // Thermodynamic density needs to be updated by psi*d(p) after the
+    // pressure solution - done in 2 parts. Part 1:
+    thermo.rho() -= psi*p_rgh;
+
+    volScalarField rAU(1.0/UEqn.A());
+    surfaceScalarField rhorAUf("rhorAUf", fvc::interpolate(rho*rAU));
+    volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p_rgh));
+
+    surfaceScalarField phig(-rhorAUf*ghf*fvc::snGrad(rho)*mesh.magSf());
+
+    surfaceScalarField phiHbyA
     (
-        p_rghDDtEqn
-      - fvm::laplacian(rhorAUf, p_rgh)
-    );
-
-    p_rghEqn.setReference
-    (
-        pRefCell,
-        compressible ? getRefCellValue(p_rgh, pRefCell) : pRefValue
-    );
-
-    p_rghEqn.solve(mesh.solver(p_rgh.select(pimple.finalInnerIter())));
-
-    if (pimple.finalNonOrthogonalIter())
-    {
-        // Calculate the conservative fluxes
-        phi = phiHbyA + p_rghEqn.flux();
-
-        // Explicitly relax pressure for momentum corrector
-        p_rgh.relax();
-
-        // Correct the momentum source with the pressure gradient flux
-        // calculated from the relaxed pressure
-        U = HbyA + rAU*fvc::reconstruct((phig + p_rghEqn.flux())/rhorAUf);
-        U.correctBoundaryConditions();
-        fvOptions.correct(U);
-        K = 0.5*magSqr(U);
-    }
-}
-
-p = p_rgh + rho*gh;
-
-#include "rhoEqn.H"
-#include "compressibleContinuityErrs.H"
-
-if (p_rgh.needReference())
-{
-    if (!compressible)
-    {
-        p += dimensionedScalar
+        "phiHbyA",
         (
-            "p",
-            p.dimensions(),
-            pRefValue - getRefCellValue(p, pRefCell)
+            fvc::flux(rho*HbyA)
+          + rhorAUf*fvc::ddtCorr(rho, U, phi)
+        )
+      + phig
+    );
+
+    MRF.makeRelative(fvc::interpolate(rho), phiHbyA);
+
+    // Update the pressure BCs to ensure flux consistency
+    constrainPressure(p_rgh, rho, U, phiHbyA, rhorAUf, MRF);
+
+    tmp<fvScalarMatrix> p_rghDDtEqn
+    (
+        new fvScalarMatrix(p_rgh, dimMass/dimTime)
+    );
+
+    if (compressible)
+    {
+        p_rghDDtEqn =
+        (
+            fvc::ddt(rho) + psi*correction(fvm::ddt(p_rgh))
         );
     }
-    else
+
+    while (pimple.correctNonOrthogonal())
     {
-        p += (initialMass - fvc::domainIntegrate(psi*p))
-            /compressibility;
-        thermo.correctRho(psi*p - psip0);
-        rho = thermo.rho();
+        fvScalarMatrix p_rghEqn
+        (
+            p_rghDDtEqn()
+          + fvc::div(phiHbyA)
+          - fvm::laplacian(rhorAUf, p_rgh)
+        ==
+            fvOptions(psi, p_rgh, rho.name())
+        );
+
+        p_rghEqn.solve(mesh.solver(p_rgh.select(pimple.finalInnerIter())));
+
+        if (pimple.finalNonOrthogonalIter())
+        {
+            // Calculate the conservative fluxes
+            phi = phiHbyA + p_rghEqn.flux();
+
+            // Explicitly relax pressure for momentum corrector
+            p_rgh.relax();
+
+            // Correct the momentum source with the pressure gradient flux
+            // calculated from the relaxed pressure
+            U = HbyA + rAU*fvc::reconstruct((phig + p_rghEqn.flux())/rhorAUf);
+            U.correctBoundaryConditions();
+            fvOptions.correct(U);
+            K = 0.5*magSqr(U);
+        }
+    }
+
+    p = p_rgh + rho*gh;
+
+    // Second part of thermodynamic density update
+    thermo.rho() += psi*p_rgh;
+
+    if (thermo.dpdt())
+    {
+        dpdt = fvc::ddt(p);
+    }
+
+    if (compressible)
+    {
+        #include "rhoEqn.H"
+    }
+    #include "compressibleContinuityErrs.H"
+
+    if (closedVolume)
+    {
+        if (!compressible)
+        {
+            p += dimensionedScalar
+            (
+                "p",
+                p.dimensions(),
+                pRefValue - getRefCellValue(p, pRefCell)
+            );
+        }
+        else
+        {
+            p += (initialMass - fvc::domainIntegrate(thermo.rho()))
+                /compressibility;
+            rho = thermo.rho();
+        }
         p_rgh = p - rho*gh;
     }
-}
-else
-{
-    thermo.correctRho(psi*p - psip0);
-}
 
-rho = thermo.rho();
-
-if (thermo.dpdt())
-{
-    dpdt = fvc::ddt(p);
+    Info<< "rho max/min : " << max(rho).value() << " " << min(rho).value()
+        << endl;
 }

applications/solvers/heatTransfer/buoyantSimpleFoam/EEqn.H

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

applications/solvers/heatTransfer/buoyantSimpleFoam/createFields.H

@@ -88,26 +88,3 @@ dimensionedScalar totalVolume = sum(mesh.V());
 
 #include "createMRF.H"
 #include "createRadiationModel.H"
-
-dimensionedScalar rhoMax
-(
-    dimensionedScalar::lookupOrDefault
-    (
-        "rhoMax",
-        simple.dict(),
-        dimDensity,
-        GREAT
-    )
-);
-
-dimensionedScalar rhoMin
-(
-    dimensionedScalar::lookupOrDefault
-    (
-        "rhoMin",
-        simple.dict(),
-        dimDensity,
-        0
-    )
-);
-

applications/solvers/heatTransfer/buoyantSimpleFoam/pEqn.H

@@ -1,4 +1,7 @@
 {
+    rho = thermo.rho();
+    rho.relax();
+
     volScalarField rAU("rAU", 1.0/UEqn.A());
     surfaceScalarField rhorAUf("rhorAUf", fvc::interpolate(rho*rAU));
     volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p_rgh));
@@ -21,8 +24,6 @@
     // Update the pressure BCs to ensure flux consistency
     constrainPressure(p_rgh, rho, U, phiHbyA, rhorAUf, MRF);
 
-    dimensionedScalar compressibility = fvc::domainIntegrate(psi);
-    bool compressible = (compressibility.value() > SMALL);
 
     while (simple.correctNonOrthogonal())
     {
@@ -32,7 +33,6 @@
         );
 
         p_rghEqn.setReference(pRefCell, getRefCellValue(p_rgh, pRefCell));
-
         p_rghEqn.solve();
 
         if (simple.finalNonOrthogonalIter())
@@ -51,9 +51,12 @@
         }
     }
 
+    #include "continuityErrs.H"
+
     p = p_rgh + rho*gh;
 
-    #include "continuityErrs.H"
+    dimensionedScalar compressibility = fvc::domainIntegrate(psi);
+    bool compressible = (compressibility.value() > SMALL);
 
     // For closed-volume cases adjust the pressure level
     // to obey overall mass continuity
@@ -72,7 +75,6 @@
         {
             p += (initialMass - fvc::domainIntegrate(psi*p))
                 /fvc::domainIntegrate(psi);
-
         }
         p_rgh = p - rho*gh;
     }

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

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

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

@@ -1,4 +1,9 @@
 {
+    rho = thermo.rho();
+    rho = max(rho, rhoMin[i]);
+    rho = min(rho, rhoMax[i]);
+    rho.relax();
+
     volScalarField rAU("rAU", 1.0/UEqn.A());
     surfaceScalarField rhorAUf("rhorAUf", fvc::interpolate(rho*rAU));
     volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p_rgh));
@@ -32,7 +37,11 @@
             fvm::laplacian(rhorAUf, p_rgh) == fvc::div(phiHbyA)
         );
 
-        p_rghEqn.setReference(pRefCell, getRefCellValue(p_rgh, pRefCell));
+        p_rghEqn.setReference
+        (
+            pRefCell,
+            compressible ? getRefCellValue(p_rgh, pRefCell) : pRefValue
+        );
 
         p_rghEqn.solve();
 
@@ -58,22 +67,10 @@
 
     // For closed-volume cases adjust the pressure level
     // to obey overall mass continuity
-    if (closedVolume)
+    if (closedVolume && compressible)
     {
-        if (!compressible)
-        {
-            p += dimensionedScalar
-            (
-                "p",
-                p.dimensions(),
-                pRefValue - getRefCellValue(p, pRefCell)
-            );
-        }
-        else
-        {
-            p += (initialMass - fvc::domainIntegrate(psi*p))
-                /compressibility;
-        }
+        p += (initialMass - fvc::domainIntegrate(thermo.rho()))
+            /compressibility;
         p_rgh = p - rho*gh;
     }
 

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

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