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

.gitmodules

@@ -1,6 +0,0 @@
-[submodule "cfmesh"]
-	path = modules/cfmesh
-	url = https://develop.openfoam.com/Community/integration-cfmesh.git
-[submodule "avalanche"]
-	path = modules/avalanche
-	url = https://develop.openfoam.com/Community/avalanche.git

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" ]
@@ -43,15 +38,6 @@ echo "Compile OpenFOAM applications"
 echo
 applications/Allwmake $targetType $*
 
-# Additional components/modules
-if [ -d "$WM_PROJECT_DIR/modules" ]
-then
-    echo "========================================"
-    echo "Compile OpenFOAM modules"
-    echo
-    (cd $WM_PROJECT_DIR/modules 2>/dev/null && wmake -all)
-fi
-
 # Some summary information
 echo
 date "+%Y-%m-%d %H:%M:%S %z" 2>/dev/null || echo "date is unknown"

BuildIssues.txt

@@ -1,23 +1,8 @@
-OpenFOAM-1712
+OpenFOAM-v1612+
 ==================
 Known Build Issues
 ==================
 
----------------------
-Intel MPI (Gcc/Clang)
----------------------
-
-  Either I_MPI_ROOT or MPI_ROOT can be used to specify the Intel-MPI
-  installation directory path.
-
-  The ThirdParty build of ptscotch uses `mpiicc` for Intel-MPI
-  instead of the usual `mpicc`.
-  When gcc or clang are used, it is highly likely that the
-  I_MPI_CC environment variable also needs to be set accordingly.
-
-  See `mpiicc -help` for more information about environment variables.
-
-
 --------------
 Intel Compiler
 --------------
@@ -37,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:
@@ -74,48 +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.
 
-
----------------------------------
-ThirdParty clang without gmp/mpfr
----------------------------------
-
-If using ThirdParty clang without gmp/mpfr, the ThirdParty makeCGAL
-script will need to be run manually and specify that there is no
-gmp/mpfr. Eg,
-
-    cd $WM_THIRD_PARTY_DIR
-    ./makeCGAL gmp-none mpfr-none
-
-Subequent compilation with Allwmake will now run largely without any
-problems, except that the components linking against CGAL
-(foamyMesh and surfaceBooleanFeatures) will also try to link against
-a nonexistent mpfr library. As a workaround, the link-dependency can
-be removed in wmake/rules/General/CGAL :
-
-    CGAL_LIBS = \
-        -L$(BOOST_ARCH_PATH)/lib \
-        -L$(BOOST_ARCH_PATH)/lib$(WM_COMPILER_LIB_ARCH) \
-        -L$(CGAL_ARCH_PATH)/lib \
-        -L$(CGAL_ARCH_PATH)/lib$(WM_COMPILER_LIB_ARCH) \
-        -lCGAL
-
-This is a temporary inconvenience until a more robust solution is found.
-
-
--------------------------
-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

@@ -1,8 +1,7 @@
 # About OpenFOAM
 OpenFOAM is a free, open source CFD software [released and developed primarily by OpenCFD Ltd](http://www.openfoam.com) since 2004. It has a large user base across most areas of engineering and science, from both commercial and academic organisations. OpenFOAM has an extensive range of features to solve anything from complex fluid flows involving chemical reactions, turbulence and heat transfer, to acoustics, solid mechanics and electromagnetics.  [More...](http://www.openfoam.com/documentation)
 
-
-OpenFOAM is professionally released every six months to include customer sponsored developments and contributions from the community - individual and group contributors, fork re-integrations including from FOAM-extend and OpenFOAM Foundation Ltd - in this Official Release sanctioned by the OpenFOAM Worldwide Trademark Owner aiming towards one OpenFOAM.
+OpenFOAM+ is professionally released every six months to include customer sponsored developments and contributions from the community, including the OpenFOAM Foundation. Releases designated OpenFOAM+ contain several man years of client-sponsored developments of which much has been transferred to, but not released in the OpenFOAM Foundation branch.
 
 
 # Copyright
@@ -10,11 +9,11 @@ OpenFOAM is free software: you can redistribute it and/or modify it under the te
 
 
 # OpenFOAM Trademark
-OpenCFD Ltd grants use of its OpenFOAM trademark by Third Parties on a licence basis. ESI Group and OpenFOAM Foundation Ltd are currently permitted to use the Name and agreed Domain Name. For information on trademark use, please refer to the [trademark policy guidelines](http://www.openfoam.com/legal/trademark-policy.php).
+OpenCFD Ltd grants use of the OpenFOAM trademark by Third Parties on a licence basis. ESI Group and the OpenFOAM Foundation Ltd are currently permitted to use the Name and agreed Domain Name. For information on trademark use, please refer to the [trademark policy guidelines](http://www.openfoam.com/legal/trademark-policy.php).
 
 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/potentialFoam/overPotentialFoam/Make/files

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

applications/solvers/basic/potentialFoam/overPotentialFoam/Make/options

@@ -1,12 +0,0 @@
-EXE_INC = \
-    -I$(LIB_SRC)/finiteVolume/lnInclude \
-    -I$(LIB_SRC)/dynamicFvMesh/lnInclude \
-    -I$(LIB_SRC)/overset/lnInclude \
-    -I$(LIB_SRC)/meshTools/lnInclude \
-    -I$(LIB_SRC)/sampling/lnInclude
-
-EXE_LIBS = \
-    -lfiniteVolume \
-    -lmeshTools \
-    -lsampling \
-    -loverset

applications/solvers/basic/potentialFoam/overPotentialFoam/createControls.H

@@ -1,9 +0,0 @@
-const dictionary& potentialFlow
-(
-    mesh.solutionDict().subDict("potentialFlow")
-);
-
-const int nNonOrthCorr
-(
-    potentialFlow.lookupOrDefault<int>("nNonOrthogonalCorrectors", 0)
-);

applications/solvers/basic/potentialFoam/overPotentialFoam/createFields.H

@@ -1,146 +0,0 @@
-Info<< "Reading velocity field U\n" << endl;
-volVectorField U
-(
-    IOobject
-    (
-        "U",
-        runTime.timeName(),
-        mesh,
-        IOobject::MUST_READ,
-        IOobject::AUTO_WRITE
-    ),
-    mesh
-);
-
-// Initialise the velocity internal field to zero
-U = dimensionedVector("0", U.dimensions(), Zero);
-
-surfaceScalarField phi
-(
-    IOobject
-    (
-        "phi",
-        runTime.timeName(),
-        mesh,
-        IOobject::NO_READ,
-        IOobject::AUTO_WRITE
-    ),
-    fvc::flux(U)
-);
-
-if (args.optionFound("initialiseUBCs"))
-{
-    U.correctBoundaryConditions();
-    phi = fvc::flux(U);
-}
-
-
-// Construct a pressure field
-// If it is available read it otherwise construct from the velocity BCs
-// converting fixed-value BCs to zero-gradient and vice versa.
-word pName("p");
-
-// Update name of the pressure field from the command-line option
-args.optionReadIfPresent("pName", pName);
-
-// Infer the pressure BCs from the velocity
-wordList pBCTypes
-(
-    U.boundaryField().size(),
-    fixedValueFvPatchScalarField::typeName
-);
-
-forAll(U.boundaryField(), patchi)
-{
-    if (U.boundaryField()[patchi].fixesValue())
-    {
-        pBCTypes[patchi] = zeroGradientFvPatchScalarField::typeName;
-    }
-}
-
-// Note that registerObject is false for the pressure field. The pressure
-// field in this solver doesn't have a physical value during the solution.
-// It shouldn't be looked up and used by sub models or boundary conditions.
-Info<< "Constructing pressure field " << pName << nl << endl;
-volScalarField p
-(
-    IOobject
-    (
-        pName,
-        runTime.timeName(),
-        mesh,
-        IOobject::READ_IF_PRESENT,
-        IOobject::NO_WRITE,
-        false
-    ),
-    mesh,
-    dimensionedScalar(pName, sqr(dimVelocity), 0),
-    pBCTypes
-);
-
-// Infer the velocity potential BCs from the pressure
-wordList PhiBCTypes
-(
-    p.boundaryField().size(),
-    zeroGradientFvPatchScalarField::typeName
-);
-
-forAll(p.boundaryField(), patchi)
-{
-    if (p.boundaryField()[patchi].fixesValue())
-    {
-        PhiBCTypes[patchi] = fixedValueFvPatchScalarField::typeName;
-    }
-}
-
-Info<< "Constructing velocity potential field Phi\n" << endl;
-volScalarField Phi
-(
-    IOobject
-    (
-        "Phi",
-        runTime.timeName(),
-        mesh,
-        IOobject::READ_IF_PRESENT,
-        IOobject::NO_WRITE
-    ),
-    mesh,
-    dimensionedScalar("Phi", dimLength*dimVelocity, 0),
-    PhiBCTypes
-);
-
-label PhiRefCell = 0;
-scalar PhiRefValue = 0;
-setRefCell
-(
-    Phi,
-    potentialFlow.dict(),
-    PhiRefCell,
-    PhiRefValue
-);
-mesh.setFluxRequired(Phi.name());
-
-#include "createMRF.H"
-
-// Add overset specific interpolations
-{
-    dictionary oversetDict;
-    oversetDict.add("Phi", true);
-    oversetDict.add("U", true);
-
-    const_cast<dictionary&>
-    (
-        mesh.schemesDict()
-    ).add
-    (
-        "oversetInterpolationRequired",
-        oversetDict,
-        true
-    );
-}
-
-// Mask field for zeroing out contributions on hole cells
-#include "createCellMask.H"
-
-// Create bool field with interpolated cells
-#include "createInterpolatedCells.H"

applications/solvers/basic/potentialFoam/overPotentialFoam/potentialFoam.C

@@ -1,260 +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
-    potentialFoam
-
-Group
-    grpBasicSolvers
-
-Description
-    Potential flow solver which solves for the velocity potential, to
-    calculate the flux-field, from which the velocity field is obtained by
-    reconstructing the flux.
-
-    \heading Solver details
-    The potential flow solution is typically employed to generate initial fields
-    for full Navier-Stokes codes.  The flow is evolved using the equation:
-
-    \f[
-        \laplacian \Phi = \div(\vec{U})
-    \f]
-
-    Where:
-    \vartable
-        \Phi      | Velocity potential [m2/s]
-        \vec{U}   | Velocity [m/s]
-    \endvartable
-
-    The corresponding pressure field could be calculated from the divergence
-    of the Euler equation:
-
-    \f[
-        \laplacian p + \div(\div(\vec{U}\otimes\vec{U})) = 0
-    \f]
-
-    but this generates excessive pressure variation in regions of large
-    velocity gradient normal to the flow direction.  A better option is to
-    calculate the pressure field corresponding to velocity variation along the
-    stream-lines:
-
-    \f[
-        \laplacian p + \div(\vec{F}\cdot\div(\vec{U}\otimes\vec{U})) = 0
-    \f]
-    where the flow direction tensor \f$\vec{F}\f$ is obtained from
-    \f[
-        \vec{F} = \hat{\vec{U}}\otimes\hat{\vec{U}}
-    \f]
-
-    \heading Required fields
-    \plaintable
-        U         | Velocity [m/s]
-    \endplaintable
-
-    \heading Optional fields
-    \plaintable
-        p         | Kinematic pressure [m2/s2]
-        Phi       | Velocity potential [m2/s]
-                  | Generated from p (if present) or U if not present
-    \endplaintable
-
-    \heading Options
-    \plaintable
-        -writep   | write the Euler pressure
-        -writePhi | Write the final velocity potential
-        -initialiseUBCs | Update the velocity boundaries before solving for Phi
-    \endplaintable
-
-\*---------------------------------------------------------------------------*/
-
-#include "fvCFD.H"
-#include "pisoControl.H"
-#include "dynamicFvMesh.H"
-#include "cellCellStencilObject.H"
-#include "localMin.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-int main(int argc, char *argv[])
-{
-    argList::addOption
-    (
-        "pName",
-        "pName",
-        "Name of the pressure field"
-    );
-
-    argList::addBoolOption
-    (
-        "initialiseUBCs",
-        "Initialise U boundary conditions"
-    );
-
-    argList::addBoolOption
-    (
-        "writePhi",
-        "Write the final velocity potential field"
-    );
-
-    argList::addBoolOption
-    (
-        "writep",
-        "Calculate and write the Euler pressure field"
-    );
-
-    argList::addBoolOption
-    (
-        "withFunctionObjects",
-        "execute functionObjects"
-    );
-
-    #include "setRootCase.H"
-    #include "createTime.H"
-    #include "createNamedDynamicFvMesh.H"
-
-    pisoControl potentialFlow(mesh, "potentialFlow");
-
-    #include "createFields.H"
-
-    // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-    Info<< nl << "Calculating potential flow" << endl;
-
-    mesh.update();
-
-    surfaceScalarField faceMask(localMin<scalar>(mesh).interpolate(cellMask));
-
-    // Since solver contains no time loop it would never execute
-    // function objects so do it ourselves
-    runTime.functionObjects().start();
-
-    MRF.makeRelative(phi);
-    adjustPhi(phi, U, p);
-
-    // Non-orthogonal velocity potential corrector loop
-    while (potentialFlow.correct())
-    {
-        phi = fvc::flux(U);
-
-        while (potentialFlow.correctNonOrthogonal())
-        {
-            fvScalarMatrix PhiEqn
-            (
-                fvm::laplacian(faceMask, Phi)
-            ==
-                fvc::div(phi)
-            );
-
-            PhiEqn.setReference(PhiRefCell, PhiRefValue);
-            PhiEqn.solve();
-
-            if (potentialFlow.finalNonOrthogonalIter())
-            {
-                phi -= PhiEqn.flux();
-            }
-        }
-
-        MRF.makeAbsolute(phi);
-
-        Info<< "Continuity error = "
-            << mag(fvc::div(phi))().weightedAverage(mesh.V()).value()
-            << endl;
-
-        U = fvc::reconstruct(phi);
-        U.correctBoundaryConditions();
-
-        Info<< "Interpolated velocity error = "
-            << (sqrt(sum(sqr(fvc::flux(U) - phi)))/sum(mesh.magSf())).value()
-            << endl;
-    }
-
-    // Write U and phi
-    U.write();
-    phi.write();
-
-    // Optionally write Phi
-    if (args.optionFound("writePhi"))
-    {
-        Phi.write();
-    }
-
-    // Calculate the pressure field from the Euler equation
-    if (args.optionFound("writep"))
-    {
-        Info<< nl << "Calculating approximate pressure field" << endl;
-
-        label pRefCell = 0;
-        scalar pRefValue = 0.0;
-        setRefCell
-        (
-            p,
-            potentialFlow.dict(),
-            pRefCell,
-            pRefValue
-        );
-
-        // Calculate the flow-direction filter tensor
-        volScalarField magSqrU(magSqr(U));
-        volSymmTensorField F(sqr(U)/(magSqrU + SMALL*average(magSqrU)));
-
-        // Calculate the divergence of the flow-direction filtered div(U*U)
-        // Filtering with the flow-direction generates a more reasonable
-        // pressure distribution in regions of high velocity gradient in the
-        // direction of the flow
-        volScalarField divDivUU
-        (
-            fvc::div
-            (
-                F & fvc::div(phi, U),
-                "div(div(phi,U))"
-            )
-        );
-
-        // Solve a Poisson equation for the approximate pressure
-        while (potentialFlow.correctNonOrthogonal())
-        {
-            fvScalarMatrix pEqn
-            (
-                fvm::laplacian(p) + divDivUU
-            );
-
-            pEqn.setReference(pRefCell, pRefValue);
-            pEqn.solve();
-        }
-
-        p.write();
-    }
-
-    runTime.functionObjects().end();
-
-    Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
-        << "  ClockTime = " << runTime.elapsedClockTime() << " s"
-        << nl << endl;
-
-    Info<< "End\n" << endl;
-
-    return 0;
-}
-
-
-// ************************************************************************* //

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

@@ -40,14 +40,15 @@ Foam::autoPtr<Foam::PDRDragModel> Foam::PDRDragModel::New
 
     Info<< "Selecting flame-wrinkling model " << modelType << endl;
 
-    auto cstrIter = dictionaryConstructorTablePtr_->cfind(modelType);
+    dictionaryConstructorTable::iterator cstrIter =
+        dictionaryConstructorTablePtr_->find(modelType);
 
-    if (!cstrIter.found())
+    if (cstrIter == dictionaryConstructorTablePtr_->end())
     {
         FatalErrorInFunction
             << "Unknown PDRDragModel type "
             << modelType << nl << nl
-            << "Valid PDRDragModel types :" << endl
+            << "Valid  PDRDragModels are : " << endl
             << dictionaryConstructorTablePtr_->sortedToc()
             << exit(FatalError);
     }

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

@@ -39,14 +39,15 @@ Foam::autoPtr<Foam::XiEqModel> Foam::XiEqModel::New
 
     Info<< "Selecting flame-wrinkling model " << modelType << endl;
 
-    auto cstrIter = dictionaryConstructorTablePtr_->cfind(modelType);
+    dictionaryConstructorTable::iterator cstrIter =
+        dictionaryConstructorTablePtr_->find(modelType);
 
-    if (!cstrIter.found())
+    if (cstrIter == dictionaryConstructorTablePtr_->end())
     {
         FatalErrorInFunction
             << "Unknown XiEqModel type "
             << modelType << nl << nl
-            << "Valid XiEqModel types :" << endl
+            << "Valid XiEqModels are : " << endl
             << dictionaryConstructorTablePtr_->sortedToc()
             << exit(FatalError);
     }

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

@@ -39,14 +39,15 @@ Foam::autoPtr<Foam::XiGModel> Foam::XiGModel::New
 
     Info<< "Selecting flame-wrinkling model " << modelType << endl;
 
-    auto cstrIter = dictionaryConstructorTablePtr_->cfind(modelType);
+    dictionaryConstructorTable::iterator cstrIter =
+        dictionaryConstructorTablePtr_->find(modelType);
 
-    if (!cstrIter.found())
+    if (cstrIter == dictionaryConstructorTablePtr_->end())
     {
         FatalErrorInFunction
             << "Unknown XiGModel type "
             << modelType << nl << nl
-            << "Valid XiGModel types :" << endl
+            << "Valid XiGModels are : " << endl
             << dictionaryConstructorTablePtr_->sortedToc()
             << exit(FatalError);
     }

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

@@ -42,14 +42,15 @@ Foam::autoPtr<Foam::XiModel> Foam::XiModel::New
 
     Info<< "Selecting flame-wrinkling model " << modelType << endl;
 
-    auto cstrIter = dictionaryConstructorTablePtr_->cfind(modelType);
+    dictionaryConstructorTable::iterator cstrIter =
+        dictionaryConstructorTablePtr_->find(modelType);
 
-    if (!cstrIter.found())
+    if (cstrIter == dictionaryConstructorTablePtr_->end())
     {
         FatalErrorInFunction
             << "Unknown XiModel type "
             << modelType << nl << nl
-            << "Valid XiModel types :" << endl
+            << "Valid XiModels are : " << endl
             << dictionaryConstructorTablePtr_->sortedToc()
             << exit(FatalError);
     }

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/chemFoam.C

@@ -42,7 +42,7 @@ Description
 #include "OFstream.H"
 #include "thermoPhysicsTypes.H"
 #include "basicMultiComponentMixture.H"
-#include "cellModel.H"
+#include "cellModeller.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/createSingleCellMesh.H

@@ -13,7 +13,8 @@ points[5] = vector(1, 0, 1);
 points[6] = vector(1, 1, 1);
 points[7] = vector(0, 1, 1);
 
-faceList faces = cellModel::ref(cellModel::HEX).modelFaces();
+const cellModel& hexa = *(cellModeller::lookup("hex"));
+faceList faces = hexa.modelFaces();
 
 fvMesh mesh
 (
@@ -24,7 +25,7 @@ fvMesh mesh
         runTime,
         IOobject::READ_IF_PRESENT
     ),
-    xferMove<pointField>(points),
+    xferMove<Field<vector>>(points),
     faces.xfer(),
     owner.xfer(),
     neighbour.xfer()

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/createFieldRefs.H

@@ -1,4 +1,4 @@
 const volScalarField& psi = thermo.psi();
 const volScalarField& T = thermo.T();
-regionModels::surfaceFilmModel& surfaceFilm = tsurfaceFilm();
+filmModelType& surfaceFilm = tsurfaceFilm();
 const label inertIndex(composition.species()[inertSpecie]);

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/fireFoam/createSurfaceFilmModel.H

@@ -1,6 +1,5 @@
 Info<< "\nConstructing surface film model" << endl;
 
-autoPtr<regionModels::surfaceFilmModel> tsurfaceFilm
-(
-    regionModels::surfaceFilmModel::New(mesh, g)
-);
+typedef regionModels::surfaceFilmModels::surfaceFilmModel filmModelType;
+
+autoPtr<filmModelType> tsurfaceFilm(filmModelType::New(mesh, g));

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/T/smoluchowskiJumpTFvPatchScalarField.C

@@ -209,14 +209,16 @@ void Foam::smoluchowskiJumpTFvPatchScalarField::write(Ostream& os) const
 {
     fvPatchScalarField::write(os);
 
-    os.writeEntryIfDifferent<word>("U", "U", UName_);
-    os.writeEntryIfDifferent<word>("rho", "rho", rhoName_);
-    os.writeEntryIfDifferent<word>("psi", "thermo:psi", psiName_);
-    os.writeEntryIfDifferent<word>("mu", "thermo:mu", muName_);
+    writeEntryIfDifferent<word>(os, "U", "U", UName_);
+    writeEntryIfDifferent<word>(os, "rho", "rho", rhoName_);
+    writeEntryIfDifferent<word>(os, "psi", "thermo:psi", psiName_);
+    writeEntryIfDifferent<word>(os, "mu", "thermo:mu", muName_);
 
-    os.writeEntry("accommodationCoeff", accommodationCoeff_);
+    os.writeKeyword("accommodationCoeff")
+        << accommodationCoeff_ << token::END_STATEMENT << nl;
     Twall_.writeEntry("Twall", os);
-    os.writeEntry("gamma", gamma_);
+    os.writeKeyword("gamma")
+        << gamma_ << token::END_STATEMENT << nl;
     writeEntry("value", os);
 }
 

applications/solvers/compressible/rhoCentralFoam/BCs/T/smoluchowskiJumpTFvPatchScalarField.H

@@ -32,8 +32,8 @@ SourceFiles
 
 \*---------------------------------------------------------------------------*/
 
-#ifndef smoluchowskiJumpTFvPatchScalarField_H
-#define smoluchowskiJumpTFvPatchScalarField_H
+#ifndef smoluchowskiJumpTFvPatchScalarFields_H
+#define smoluchowskiJumpTFvPatchScalarFields_H
 
 #include "mixedFvPatchFields.H"
 
@@ -43,7 +43,7 @@ namespace Foam
 {
 
 /*---------------------------------------------------------------------------*\
-            Class smoluchowskiJumpTFvPatchScalarField Declaration
+                 Class smoluchowskiJumpTFvPatch Declaration
 \*---------------------------------------------------------------------------*/
 
 class smoluchowskiJumpTFvPatchScalarField
@@ -74,7 +74,6 @@ class smoluchowskiJumpTFvPatchScalarField
         //- Heat capacity ratio (default 1.4)
         scalar gamma_;
 
-
 public:
 
     //- Runtime type information

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,23 +193,25 @@ void Foam::maxwellSlipUFvPatchVectorField::updateCoeffs()
         refValue() -= C1/prho*transform(I - n*n, (n & ptauMC));
     }
 
-    partialSlipFvPatchVectorField::updateCoeffs();
+    mixedFixedValueSlipFvPatchVectorField::updateCoeffs();
 }
 
 
 void Foam::maxwellSlipUFvPatchVectorField::write(Ostream& os) const
 {
     fvPatchVectorField::write(os);
-    os.writeEntryIfDifferent<word>("T", "T", TName_);
-    os.writeEntryIfDifferent<word>("rho", "rho", rhoName_);
-    os.writeEntryIfDifferent<word>("psi", "thermo:psi", psiName_);
-    os.writeEntryIfDifferent<word>("mu", "thermo:mu", muName_);
-    os.writeEntryIfDifferent<word>("tauMC", "tauMC", tauMCName_);
+    writeEntryIfDifferent<word>(os, "T", "T", TName_);
+    writeEntryIfDifferent<word>(os, "rho", "rho", rhoName_);
+    writeEntryIfDifferent<word>(os, "psi", "thermo:psi", psiName_);
+    writeEntryIfDifferent<word>(os, "mu", "thermo:mu", muName_);
+    writeEntryIfDifferent<word>(os, "tauMC", "tauMC", tauMCName_);
 
-    os.writeEntry("accommodationCoeff", accommodationCoeff_);
+    os.writeKeyword("accommodationCoeff")
+        << accommodationCoeff_ << token::END_STATEMENT << nl;
     Uwall_.writeEntry("Uwall", os);
-    os.writeEntry("thermalCreep", thermalCreep_);
-    os.writeEntry("curvature", curvature_);
+    os.writeKeyword("thermalCreep")
+        << thermalCreep_ << token::END_STATEMENT << nl;
+    os.writeKeyword("curvature") << curvature_ << token::END_STATEMENT << nl;
 
     refValue().writeEntry("refValue", os);
     valueFraction().writeEntry("valueFraction", os);

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

@@ -2,10 +2,8 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           |
+    \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
      \\/     M anipulation  |
--------------------------------------------------------------------------------
-                            | Copyright (C) 2016-2017 Wikki Ltd
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -24,23 +22,22 @@ License
     along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
 
 Class
-    Foam::mixedFaPatchField
+    Foam::mixedFixedValueSlipFvPatchField
 
 Description
-
-Author
-    Zeljko Tukovic, FMENA
-    Hrvoje Jasak, Wikki Ltd.
+    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
-    mixedFaPatchField.C
+    mixedFixedValueSlipFvPatchField.C
 
 \*---------------------------------------------------------------------------*/
 
-#ifndef mixedFaPatchField_H
-#define mixedFaPatchField_H
+#ifndef mixedFixedValueSlipFvPatchField_H
+#define mixedFixedValueSlipFvPatchField_H
 
-#include "faPatchField.H"
+#include "transformFvPatchField.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -48,106 +45,115 @@ namespace Foam
 {
 
 /*---------------------------------------------------------------------------*\
-                      Class mixedFaPatchField Declaration
+                      Class mixedFixedValueSlipFvPatch Declaration
 \*---------------------------------------------------------------------------*/
 
 template<class Type>
-class mixedFaPatchField
+class mixedFixedValueSlipFvPatchField
 :
-    public faPatchField<Type>
+    public transformFvPatchField<Type>
 {
     // Private data
 
-        //- Value field
+        //- Value field used for boundary condition
         Field<Type> refValue_;
 
-        //- Normal gradient field
-        Field<Type> refGrad_;
-
         //- Fraction (0-1) of value used for boundary condition
         scalarField valueFraction_;
 
-
 public:
 
     //- Runtime type information
-    TypeName("mixed");
+    TypeName("mixedFixedValueSlip");
 
 
     // Constructors
 
         //- Construct from patch and internal field
-        mixedFaPatchField
+        mixedFixedValueSlipFvPatchField
         (
-            const faPatch&,
-            const DimensionedField<Type, areaMesh>&
+            const fvPatch&,
+            const DimensionedField<Type, volMesh>&
         );
 
         //- Construct from patch, internal field and dictionary
-        mixedFaPatchField
+        mixedFixedValueSlipFvPatchField
         (
-            const faPatch&,
-            const DimensionedField<Type, areaMesh>&,
+            const fvPatch&,
+            const DimensionedField<Type, volMesh>&,
             const dictionary&
         );
 
-        //- Construct by mapping the given mixedFaPatchField onto a new patch
-        mixedFaPatchField
+        //- Construct by mapping given mixedFixedValueSlipFvPatchField
+        //- Onto a new patch
+        mixedFixedValueSlipFvPatchField
         (
-            const mixedFaPatchField<Type>&,
-            const faPatch&,
-            const DimensionedField<Type, areaMesh>&,
-            const faPatchFieldMapper&
+            const mixedFixedValueSlipFvPatchField<Type>&,
+            const fvPatch&,
+            const DimensionedField<Type, volMesh>&,
+            const fvPatchFieldMapper&
         );
 
         //- Construct as copy
-        mixedFaPatchField
+        mixedFixedValueSlipFvPatchField
         (
-            const mixedFaPatchField<Type>&
+            const mixedFixedValueSlipFvPatchField<Type>&
         );
 
         //- Construct and return a clone
-        virtual tmp<faPatchField<Type>> clone() const
+        virtual tmp<fvPatchField<Type>> clone() const
         {
-            return tmp<faPatchField<Type>>
+            return tmp<fvPatchField<Type>>
             (
-                new mixedFaPatchField<Type>(*this)
+                new mixedFixedValueSlipFvPatchField<Type>(*this)
             );
         }
 
         //- Construct as copy setting internal field reference
-        mixedFaPatchField
+        mixedFixedValueSlipFvPatchField
         (
-            const mixedFaPatchField<Type>&,
-            const DimensionedField<Type, areaMesh>&
+            const mixedFixedValueSlipFvPatchField<Type>&,
+            const DimensionedField<Type, volMesh>&
         );
 
         //- Construct and return a clone setting internal field reference
-        virtual tmp<faPatchField<Type>> clone
+        virtual tmp<fvPatchField<Type>> clone
         (
-            const DimensionedField<Type, areaMesh>& iF
+            const DimensionedField<Type, volMesh>& iF
         ) const
         {
-            return tmp<faPatchField<Type>>
+            return tmp<fvPatchField<Type>>
             (
-                new mixedFaPatchField<Type>(*this, iF)
+                new mixedFixedValueSlipFvPatchField<Type>(*this, iF)
             );
         }
 
-
     // Member functions
 
         // Access
 
-            //- Return true if this patch field fixes a value.
-            //  Needed to check if a level has to be specified while solving
-            //  Poissons equations.
-            virtual bool fixesValue() const
+            //- Return false: this patch field is not altered by assignment
+            virtual bool assignable() const
             {
-                return true;
+                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()
@@ -160,16 +166,6 @@ public:
                 return refValue_;
             }
 
-            virtual Field<Type>& refGrad()
-            {
-                return refGrad_;
-            }
-
-            virtual const Field<Type>& refGrad() const
-            {
-                return refGrad_;
-            }
-
             virtual scalarField& valueFraction()
             {
                 return valueFraction_;
@@ -180,23 +176,6 @@ public:
                 return valueFraction_;
             }
 
-
-        // Mapping functions
-
-            //- Map (and resize as needed) from self given a mapping object
-            virtual void autoMap
-            (
-                const faPatchFieldMapper&
-            );
-
-            //- Reverse map the given faPatchField onto this faPatchField
-            virtual void rmap
-            (
-                const faPatchField<Type>&,
-                const labelList&
-            );
-
-
         // Evaluation functions
 
             //- Return gradient at boundary
@@ -205,30 +184,11 @@ public:
             //- Evaluate the patch field
             virtual void evaluate
             (
-                const Pstream::commsTypes commsType = Pstream::commsTypes::blocking
+                const Pstream::commsTypes commsType=Pstream::blocking
             );
 
-            //- Return the matrix diagonal coefficients corresponding to the
-            //  evaluation of the value of this patchField with given weights
-            virtual tmp<Field<Type>> valueInternalCoeffs
-            (
-                const tmp<scalarField>&
-            ) const;
-
-            //- Return the matrix source coefficients corresponding to the
-            //  evaluation of the value of this patchField with given weights
-            virtual tmp<Field<Type>> valueBoundaryCoeffs
-            (
-                const tmp<scalarField>&
-            ) const;
-
-            //- Return the matrix diagonal coefficients corresponding to the
-            //  evaluation of the gradient of this patchField
-            virtual tmp<Field<Type>> gradientInternalCoeffs() const;
-
-            //- Return the matrix source coefficients corresponding to the
-            //  evaluation of the gradient of this patchField
-            virtual tmp<Field<Type>> gradientBoundaryCoeffs() const;
+            //- Return face-gradient transform diagonal
+            virtual tmp<Field<Type>> snGradTransformDiag() const;
 
 
         //- Write
@@ -239,11 +199,11 @@ public:
 
         virtual void operator=(const UList<Type>&) {}
 
-        virtual void operator=(const faPatchField<Type>&) {}
-        virtual void operator+=(const faPatchField<Type>&) {}
-        virtual void operator-=(const faPatchField<Type>&) {}
-        virtual void operator*=(const faPatchField<scalar>&) {}
-        virtual void operator/=(const faPatchField<scalar>&) {}
+        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>&) {}
@@ -266,7 +226,7 @@ public:
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 #ifdef NoRepository
-    #include "mixedFaPatchField.C"
+    #include "mixedFixedValueSlipFvPatchField.C"
 #endif
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

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

@@ -23,9 +23,9 @@ License
 
 \*---------------------------------------------------------------------------*/
 
-#include "genericFaPatchFields.H"
+#include "mixedFixedValueSlipFvPatchFields.H"
 #include "addToRunTimeSelectionTable.H"
-#include "areaFields.H"
+#include "volFields.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -34,7 +34,7 @@ namespace Foam
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
-makeFaPatchFields(generic);
+makePatchFields(mixedFixedValueSlip);
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 

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

@@ -23,10 +23,10 @@ License
 
 \*---------------------------------------------------------------------------*/
 
-#ifndef genericFaPatchFields_H
-#define genericFaPatchFields_H
+#ifndef mixedFixedValueSlipFvPatchFields_H
+#define mixedFixedValueSlipFvPatchFields_H
 
-#include "genericFaPatchField.H"
+#include "mixedFixedValueSlipFvPatchField.H"
 #include "fieldTypes.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -36,7 +36,7 @@ namespace Foam
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-makeFaPatchTypeFieldTypedefs(generic);
+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) 2016-2017 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -23,10 +23,9 @@ License
 
 \*---------------------------------------------------------------------------*/
 
-#ifndef oversetFvPatchFields_H
-#define oversetFvPatchFields_H
+#ifndef mixedFixedValueSlipFvPatchFieldsFwd_H
+#define mixedFixedValueSlipFvPatchFieldsFwd_H
 
-#include "oversetFvPatchField.H"
 #include "fieldTypes.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -36,7 +35,9 @@ namespace Foam
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-makePatchTypeFieldTypedefs(overset);
+template<class Type> class mixedFixedValueSlipFvPatchField;
+
+makePatchTypeFieldTypedefs(mixedFixedValueSlip);
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 

applications/solvers/compressible/rhoCentralFoam/BCs/rho/fixedRhoFvPatchScalarField.C

@@ -117,8 +117,8 @@ void Foam::fixedRhoFvPatchScalarField::write(Ostream& os) const
 {
     fvPatchScalarField::write(os);
 
-    os.writeEntryIfDifferent<word>("p", "p", pName_);
-    os.writeEntryIfDifferent<word>("psi", "thermo:psi", psiName_);
+    writeEntryIfDifferent<word>(os, "p", "p", this->pName_);
+    writeEntryIfDifferent<word>(os, "psi", "thermo:psi", psiName_);
     writeEntry("value", os);
 }
 

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

@@ -113,9 +113,6 @@ int main(int argc, char *argv[])
             phiv_pos -= mesh.phi();
             phiv_neg -= mesh.phi();
         }
-        // Note: extracted out the orientation so becomes unoriented
-        phiv_pos.setOriented(false);
-        phiv_neg.setOriented(false);
 
         volScalarField c("c", sqrt(thermo.Cp()/thermo.Cv()*rPsi));
         surfaceScalarField cSf_pos
@@ -168,12 +165,11 @@ int main(int argc, char *argv[])
 
         phi = aphiv_pos*rho_pos + aphiv_neg*rho_neg;
 
-        surfaceVectorField phiU(aphiv_pos*rhoU_pos + aphiv_neg*rhoU_neg);
-        // Note: reassembled orientation from the pos and neg parts so becomes
-        // oriented
-        phiU.setOriented(true);
-
-        surfaceVectorField phiUp(phiU + (a_pos*p_pos + a_neg*p_neg)*mesh.Sf());
+        surfaceVectorField phiUp
+        (
+            (aphiv_pos*rhoU_pos + aphiv_neg*rhoU_neg)
+          + (a_pos*p_pos + a_neg*p_neg)*mesh.Sf()
+        );
 
         surfaceScalarField phiEp
         (
@@ -186,10 +182,7 @@ int main(int argc, char *argv[])
         // Make flux for pressure-work absolute
         if (mesh.moving())
         {
-            surfaceScalarField phia(a_pos*p_pos + a_neg*p_neg);
-            phia.setOriented(true);
-
-            phiEp += mesh.phi()*phia;
+            phiEp += mesh.phi()*(a_pos*p_pos + a_neg*p_neg);
         }
 
         volScalarField muEff("muEff", turbulence->muEff());
@@ -226,7 +219,7 @@ int main(int argc, char *argv[])
                 fvc::interpolate(muEff)*mesh.magSf()*fvc::snGrad(U)
               + fvc::dotInterpolate(mesh.Sf(), tauMC)
             )
-          & (a_pos*U_pos + a_neg*U_neg)
+            & (a_pos*U_pos + a_neg*U_neg)
         );
 
         solve
@@ -276,5 +269,4 @@ int main(int argc, char *argv[])
     return 0;
 }
 
-
 // ************************************************************************* //

applications/solvers/compressible/rhoCentralFoam/rhoCentralFoam.C

@@ -93,10 +93,7 @@ int main(int argc, char *argv[])
         surfaceScalarField p_neg("p_neg", rho_neg*rPsi_neg);
 
         surfaceScalarField phiv_pos("phiv_pos", U_pos & mesh.Sf());
-        // Note: extracted out the orientation so becomes unoriented
-        phiv_pos.setOriented(false);
         surfaceScalarField phiv_neg("phiv_neg", U_neg & mesh.Sf());
-        phiv_neg.setOriented(false);
 
         volScalarField c("c", sqrt(thermo.Cp()/thermo.Cv()*rPsi));
         surfaceScalarField cSf_pos
@@ -104,7 +101,6 @@ int main(int argc, char *argv[])
             "cSf_pos",
             interpolate(c, pos, T.name())*mesh.magSf()
         );
-
         surfaceScalarField cSf_neg
         (
             "cSf_neg",
@@ -116,7 +112,6 @@ int main(int argc, char *argv[])
             "ap",
             max(max(phiv_pos + cSf_pos, phiv_neg + cSf_neg), v_zero)
         );
-
         surfaceScalarField am
         (
             "am",
@@ -165,12 +160,11 @@ int main(int argc, char *argv[])
 
         phi = aphiv_pos*rho_pos + aphiv_neg*rho_neg;
 
-        surfaceVectorField phiU(aphiv_pos*rhoU_pos + aphiv_neg*rhoU_neg);
-        // Note: reassembled orientation from the pos and neg parts so becomes
-        // oriented
-        phiU.setOriented(true);
-
-        surfaceVectorField phiUp(phiU + (a_pos*p_pos + a_neg*p_neg)*mesh.Sf());
+        surfaceVectorField phiUp
+        (
+            (aphiv_pos*rhoU_pos + aphiv_neg*rhoU_neg)
+          + (a_pos*p_pos + a_neg*p_neg)*mesh.Sf()
+        );
 
         surfaceScalarField phiEp
         (

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/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/createFields.H

@@ -1,117 +0,0 @@
-Info<< "Reading thermophysical properties\n" << endl;
-
-autoPtr<psiThermo> pThermo
-(
-    psiThermo::New(mesh)
-);
-psiThermo& thermo = pThermo();
-thermo.validate(args.executable(), "h", "e");
-
-volScalarField& p = thermo.p();
-const volScalarField& psi = thermo.psi();
-
-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
-);
-
-#include "compressibleCreatePhi.H"
-
-dimensionedScalar rhoMax
-(
-    dimensionedScalar::lookupOrDefault
-    (
-        "rhoMax",
-        pimple.dict(),
-        dimDensity,
-        GREAT
-    )
-);
-
-dimensionedScalar rhoMin
-(
-    dimensionedScalar::lookupOrDefault
-    (
-        "rhoMin",
-        pimple.dict(),
-        dimDensity,
-        0
-    )
-);
-
-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));
-
-
-//- 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"
-
-Info<< "Creating turbulence model\n" << endl;
-autoPtr<compressible::turbulenceModel> turbulence
-(
-    compressible::turbulenceModel::New
-    (
-        rho,
-        U,
-        phi,
-        thermo
-    )
-);

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", U);
-HbyA = constrainHbyA(cellMask*rAU*UEqn.H(), U, p);
-
-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/overRhoSimpleFoam/Make/files

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

applications/solvers/compressible/rhoSimpleFoam/overRhoSimpleFoam/Make/options

@@ -1,27 +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)/sampling/lnInclude \
-    -I$(LIB_SRC)/dynamicFvMesh/lnInclude \
-    -I$(LIB_SRC)/dynamicMesh/lnInclude \
-    -I$(LIB_SRC)/meshTools/lnInclude \
-    -I$(LIB_SRC)/overset/lnInclude
-
-EXE_LIBS = \
-    -lcompressibleTransportModels \
-    -lfluidThermophysicalModels \
-    -lspecie \
-    -lturbulenceModels \
-    -lcompressibleTurbulenceModels \
-    -lfiniteVolume \
-    -lsampling \
-    -lmeshTools \
-    -lfvOptions \
-    -loverset \
-    -ldynamicFvMesh \
-    -ltopoChangerFvMesh

applications/solvers/compressible/rhoSimpleFoam/overRhoSimpleFoam/UEqn.H

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

applications/solvers/compressible/rhoSimpleFoam/overRhoSimpleFoam/createFieldRefs.H

@@ -1 +0,0 @@
-const volScalarField& psi = thermo.psi();

applications/solvers/compressible/rhoSimpleFoam/overRhoSimpleFoam/createUpdatedDynamicFvMesh.H

@@ -1,22 +0,0 @@
-    Info<< "Create dynamic mesh for time = "
-        << runTime.timeName() << nl << endl;
-
-    autoPtr<dynamicFvMesh> meshPtr
-    (
-        dynamicFvMesh::New
-        (
-            IOobject
-            (
-                dynamicFvMesh::defaultRegion,
-                runTime.timeName(),
-                runTime,
-                IOobject::MUST_READ
-            )
-        )
-    );
-
-    dynamicFvMesh& mesh = meshPtr();
-
-    // Calculate initial mesh-to-mesh mapping. Note that this should be
-    // done under the hood, e.g. as a MeshObject
-    mesh.update();

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

@@ -1,123 +0,0 @@
-{
-    surfaceScalarField faceMask(localMin<scalar>(mesh).interpolate(cellMask));
-
-    volScalarField rAU(1.0/UEqn.A());
-    surfaceScalarField rhorAUf("rhorAUf", faceMask*fvc::interpolate(rho*rAU));
-    volVectorField HbyA("HbyA", U);
-    HbyA = constrainHbyA(cellMask*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
-        );
-
-        phiHbyA -= fvc::interpolate(psi*p)*phiHbyA/fvc::interpolate(rho);
-
-        while (simple.correctNonOrthogonal())
-        {
-            fvScalarMatrix pEqn
-            (
-                fvc::div(phiHbyA)
-              + fvm::div(phid, p)
-              - fvm::laplacian(rhorAUf, p)
-            ==
-                fvOptions(psi, p, rho.name())
-            );
-
-            // Relax the pressure equation to ensure diagonal-dominance
-            pEqn.relax();
-
-            pEqn.setReference
-            (
-                pressureControl.refCell(),
-                pressureControl.refValue()
-            );
-
-            pEqn.solve();
-
-            if (simple.finalNonOrthogonalIter())
-            {
-                phi = phiHbyA + pEqn.flux();
-            }
-        }
-    }
-    else
-    {
-        closedVolume = adjustPhi(phiHbyA, U, p);
-
-        if (adjustFringe)
-        {
-            oversetAdjustPhi(phiHbyA, U);
-        }
-
-        while (simple.correctNonOrthogonal())
-        {
-            fvScalarMatrix pEqn
-            (
-                fvc::div(phiHbyA)
-              - fvm::laplacian(rhorAUf, p)
-            ==
-                fvOptions(psi, p, rho.name())
-            );
-
-            pEqn.setReference
-            (
-                pressureControl.refCell(),
-                pressureControl.refValue()
-            );
-
-            pEqn.solve();
-
-            if (simple.finalNonOrthogonalIter())
-            {
-                phi = phiHbyA + pEqn.flux();
-            }
-        }
-    }
-
-    #include "incompressible/continuityErrs.H"
-
-    // Explicitly relax pressure for momentum corrector
-    p.relax();
-
-    volVectorField gradP(fvc::grad(p));
-
-    U = HbyA - rAU*cellMask*gradP;
-    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();
-    }
-}

applications/solvers/compressible/rhoSimpleFoam/overRhoSimpleFoam/rhoSimpleFoam.C

@@ -1,92 +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
-    overRhoSimpleFoam
-
-Group
-    grpCompressibleSolvers
-
-Description
-    Overset steady-state solver for turbulent flow of compressible fluids.
-
-\*---------------------------------------------------------------------------*/
-
-#include "fvCFD.H"
-#include "dynamicFvMesh.H"
-#include "fluidThermo.H"
-#include "turbulentFluidThermoModel.H"
-#include "simpleControl.H"
-#include "pressureControl.H"
-#include "fvOptions.H"
-#include "cellCellStencilObject.H"
-#include "localMin.H"
-#include "oversetAdjustPhi.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-int main(int argc, char *argv[])
-{
-    #define CREATE_MESH createUpdatedDynamicFvMesh.H
-    #include "postProcess.H"
-
-    #include "setRootCase.H"
-    #include "createTime.H"
-    #include "createUpdatedDynamicFvMesh.H"
-    #include "createControl.H"
-    #include "createFields.H"
-    #include "createFieldRefs.H"
-    #include "createFvOptions.H"
-    #include "initContinuityErrs.H"
-
-    turbulence->validate();
-
-    // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-    Info<< "\nStarting time loop\n" << endl;
-
-    while (simple.loop())
-    {
-        Info<< "Time = " << runTime.timeName() << nl << endl;
-
-        // Pressure-velocity SIMPLE corrector
-        #include "UEqn.H"
-        #include "EEqn.H"
-        #include "pEqn.H"
-
-        turbulence->correct();
-
-        runTime.write();
-
-        Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
-            << "  ClockTime = " << runTime.elapsedClockTime() << " s"
-            << nl << endl;
-    }
-
-    Info<< "End\n" << endl;
-
-    return 0;
-}
-
-
-// ************************************************************************* //