BUG: parse issues for nastran long format (closes #589)

- the regex was not being unset on assignment.
- drop unused recompile() method as being dubious

.gitignore

@@ -6,6 +6,8 @@
 .*~
 *.bak
 *.bak[0-9][0-9]
+*.orig
+*.orig[0-9][0-9]
 \#*\#
 
 # File-browser settings - anywhere
@@ -39,12 +41,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 +60,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
@@ -81,10 +77,9 @@ doc/Doxygen/DTAGS
 # Ignore .tags in the main directory
 /.tags
 
-# Ignore project files in the main directory
+# Ignore eclipse project files in the main directory
 /.cproject
 /.project
-/.dir-locals.el
 
 # Ignore the test directory
 /tutorialsTest

Allwmake

@@ -1,54 +1,46 @@
 #!/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
+# Parse arguments for library compilation without documentation by default
+genDoc=0
+targetType=libso
 . $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" ]
 then
     $WM_THIRD_PARTY_DIR/Allwmake
 else
-    echo "No ThirdParty directory found - skipping"
+    echo "Allwmake: no ThirdParty directory found - skipping"
 fi
 
-echo "========================================"
-echo "Compile OpenFOAM libraries"
-echo
+# Compile OpenFOAM libraries and applications
 src/Allwmake $targetType $*
 
-echo "========================================"
-echo "Compile OpenFOAM applications"
-echo
+# Compile OpenFOAM libraries and applications
 applications/Allwmake $targetType $*
 
-# Some summary information
-echo
-date "+%Y-%m-%d %H:%M:%S %z" 2>/dev/null || echo "date is unknown"
-echo "========================================"
-echo "  ${WM_PROJECT_DIR##*/}"
-echo "  $WM_COMPILER $WM_COMPILER_TYPE compiler"
-echo "  ${WM_OPTIONS}, with ${WM_MPLIB} ${FOAM_MPI}"
-echo
+# Optionally build OpenFOAM Doxygen documentation
+if [ $genDoc -eq 1 ]
+then
+    doc/Allwmake
+fi
+
+
 #------------------------------------------------------------------------------

BuildIssues.txt

@@ -1,79 +0,0 @@
-OpenFOAM-1706
-==================
-Known Build Issues
-==================
-
---------------
-Intel Compiler
---------------
-
-  Since OpenFOAM uses C++11, a fairly recent version is required.
-  The Intel compiler - icc (ICC) 17.0.1 20161005 is ok, but the
-  initial release - icc (ICC) 17.0.0 20160721 - has a bug that
-  will result in these types of error messages.
-
-    MatrixSpaceI.H(492): error: no instance of overloaded function
-    "Foam::MatrixSpace<Form, Cmpt, Mrows, Ncols>::Block<SubTensor,
-    BRowStart, BColStart>::operator=" matches the specified type
-
-
----
-VTK
----
-
-If using the runTimePostProcessing to create on-the-fly images, you
-can simply just compile ParaView 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
-prior to using makeVTK. This doesn't just reduce disk-space, but works
-much better than using the VTK tar file.
-
-Using runTimePostProcessing with the 'plain' VTK libraries does
-generally work, but does not exit cleanly:
-
-    symbol lookup error: .../linux64Gcc/VTK-7.1.0/lib/libvtkCommonExecutionModel-7.1.so.1:
-    undefined symbol: _ZN33vtkFilteringInformationKeyManager13ClassFinalizeEv
-
-    symbol lookup error: .../linux64Gcc/VTK-7.1.0/lib/libvtkCommonDataModel-7.1.so.1:
-    undefined symbol: _ZN49vtkInformationQuadratureSchemeDefinitionVectorKeyD1Ev
-
-This error appears to be suppressed if VTK is compiled with a Debug build-type.
-
-
--------------------------
-Building on older systems
--------------------------
-
-If the system gcc is too old for building OpenFOAM, a third-party gcc or
-clang/llvm installation can be used. If building clang/llvm, note that
-there are also minimum gcc/g++ requirements there:
-
-   Min gcc/g++
-   ===========   ==========
-    4.4          llvm-3.4.2
-    4.7          llvm-3.5.2 - llvm-3.7.0
-
-
-If your system compiler is too old to build the minimum required gcc or
-clang/llvm, it is just simply too old.
-
-
--------------------------
-Building with spack
--------------------------
-
-If you are building with spack, note that the depends_on for paraview
-resolves poorly. The +qt dependency (for building the reader module)
-may need to be specified as a preference by including the following in
-your `~/.spack/packages.yaml` file:
-
-    packages:
-        paraview:
-            variants: +qt
-
-It appears that spack will otherwise ignore any paraview+qt version
-and attempt to install a paraview~qt version instead.
-
---

README.md

@@ -1,26 +1,22 @@
 # 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 a free, open source CFD software package developed primarily by [OpenCFD](http://www.openfoam.com) since 2004 and is currently distributed by [ESI-OpenCFD](http://www.openfoam.com) and the [OpenFOAM Foundation](http://openfoam.org). It has a large user base across most areas of engineering and science, from both commercial and academic organisations. OpenFOAM has an extensive range of features to solve anything from complex fluid flows involving chemical reactions, turbulence and heat transfer, to acoustics, solid mechanics and electromagnetics.  [More...](http://www.openfoam.com/documentation)
 
-OpenFOAM+ is professionally released every six months to include customer sponsored developments and contributions from the community, including the OpenFOAM Foundation. Releases designated OpenFOAM+ contain several man years of client-sponsored developments of which much has been transferred to, but not released in the OpenFOAM Foundation branch.
+OpenFOAM+ uses the OpenFOAM Foundation version as a common code base, with wider functionality and platform support. Its purpose is to accelerate the public availability of new features which are sponsored by OpenCFD's customers and contributed by the OpenFOAM community.
 
 
 # Copyright
-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.  See the file COPYING in this directory or [http://www.gnu.org/licenses/](http://www.gnu.org/licenses), for a description of the GNU General Public License terms under which you can copy the files.
+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.  See the file COPYING in this directory or [http://www.gnu.org/licenses/](http://www.gnu.org/licenses/), for a description of the GNU General Public License terms under which you can copy the files.
 
 
 # OpenFOAM Trademark
-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.
+OpenCFD Limited, owner of the OpenFOAM Trademark, has granted the use of the Trademark to the OpenFOAM Foundation.  Violations of the Trademark are continuously monitored and will be duly prosecuted. Please [contact OpenCFD](http://www.openfoam.com/contact) if you have any questions on the use of the OpenFOAM trademark.
 
 
 # Useful Links
-- [Download and installation instructions](http://www.openfoam.com/download/)
+- [Download and installation instructions](http://www.openfoam.com/releases)
 - [Documentation](http://www.openfoam.com/documentation)
-- [Reporting bugs/issues/feature requests](http://www.openfoam.com/code/bug-reporting.php)
-- [OpenFOAM Community](http://www.openfoam.com/community/)
-- [Contacting OpenCFD](http://www.openfoam.com/contact/)
+- [Reporting bugs/issues (including bugs/suggestions/feature requests) in OpenFOAM+](http://www.openfoam.com/code/bug-reporting.php)
+- [Collaborative and Community-based Developments](http://www.openfoam.com/services/community-projects.php)
+- [Contacting OpenCFD](http://www.openfoam.com/contact)
 
-Copyright 2016-2017 OpenCFD Ltd
+Copyright 2016 OpenCFD Ltd

applications/Allwmake

@@ -1,22 +1,27 @@
 #!/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
+targetType=libso
 . $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
+}
+
+set -x
+
+wmake -all solvers $*
+wmake -all utilities $*
 
-wmake -all $targetType solvers
-wmake -all $targetType utilities
 
 #------------------------------------------------------------------------------

applications/solvers/DNS/dnsFoam/Allwmake

@@ -6,7 +6,9 @@ if [ -f "$FFTW_ARCH_PATH/include/fftw3.h" ] || \
 then
     wmake
 else
-    echo "==> skip dnsFoam solver (no FFTW)"
+    echo
+    echo "Skipping dnsFoam solver (no FFTW)"
+    echo
 fi
 
 #------------------------------------------------------------------------------

applications/solvers/DNS/dnsFoam/createFields.H

@@ -1,5 +1,3 @@
-#include "readTransportProperties.H"
-
 Info<< "Reading field p\n" << endl;
 volScalarField p
 (
@@ -31,5 +29,3 @@ volVectorField U
 #include "createPhi.H"
 
 mesh.setFluxRequired(p.name());
-
-#include "readTurbulenceProperties.H"

applications/solvers/DNS/dnsFoam/dnsFoam.C

@@ -28,7 +28,7 @@ Group
     grpDNSSolvers
 
 Description
-    Direct numerical simulation solver for boxes of isotropic turbulence.
+    Direct numerical simulation solver for boxes of isotropic turbulence
 
 \*---------------------------------------------------------------------------*/
 
@@ -44,13 +44,16 @@ Description
 
 int main(int argc, char *argv[])
 {
-    #include "postProcess.H"
-
     #include "setRootCase.H"
+
     #include "createTime.H"
     #include "createMeshNoClear.H"
-    #include "createControl.H"
+
+    pisoControl piso(mesh);
+
+    #include "readTransportProperties.H"
     #include "createFields.H"
+    #include "readTurbulenceProperties.H"
     #include "initContinuityErrs.H"
 
     // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -61,7 +64,7 @@ int main(int argc, char *argv[])
     {
         Info<< "Time = " << runTime.timeName() << nl << endl;
 
-        force.primitiveFieldRef() = ReImSum
+        force.internalField() = ReImSum
         (
             fft::reverseTransform
             (
@@ -116,7 +119,7 @@ int main(int argc, char *argv[])
 
         runTime.write();
 
-        if (runTime.writeTime())
+        if (runTime.outputTime())
         {
             calcEk(U, K).write
             (

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

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2015 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -52,7 +52,6 @@ Description
 \*---------------------------------------------------------------------------*/
 
 #include "fvCFD.H"
-#include "fvOptions.H"
 #include "simpleControl.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -67,7 +66,6 @@ int main(int argc, char *argv[])
     simpleControl simple(mesh);
 
     #include "createFields.H"
-    #include "createFvOptions.H"
 
     // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -79,16 +77,10 @@ int main(int argc, char *argv[])
 
         while (simple.correctNonOrthogonal())
         {
-            fvScalarMatrix TEqn
+            solve
             (
                 fvm::ddt(T) - fvm::laplacian(DT, T)
-             ==
-                fvOptions(T)
             );
-
-            fvOptions.constrain(TEqn);
-            TEqn.solve();
-            fvOptions.correct(T);
         }
 
         #include "write.H"

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/laplacianFoam/write.H

@@ -1,4 +1,4 @@
-    if (runTime.writeTime())
+    if (runTime.outputTime())
     {
         volVectorField gradT(fvc::grad(T));
 

applications/solvers/basic/potentialFoam/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/createFields.H

@@ -43,7 +43,7 @@ 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
+// Infer the pressure BCs from the velocity BCs
 wordList pBCTypes
 (
     U.boundaryField().size(),
@@ -78,37 +78,63 @@ volScalarField p
     pBCTypes
 );
 
-// Infer the velocity potential BCs from the pressure
-wordList PhiBCTypes
-(
-    p.boundaryField().size(),
-    zeroGradientFvPatchScalarField::typeName
-);
-
-forAll(p.boundaryField(), patchi)
+label pRefCell = 0;
+scalar pRefValue = 0.0;
+if (args.optionFound("writep"))
 {
-    if (p.boundaryField()[patchi].fixesValue())
-    {
-        PhiBCTypes[patchi] = fixedValueFvPatchScalarField::typeName;
-    }
+    setRefCell
+    (
+        p,
+        potentialFlow.dict(),
+        pRefCell,
+        pRefValue
+    );
 }
 
+
 Info<< "Constructing velocity potential field Phi\n" << endl;
-volScalarField Phi
+autoPtr<volScalarField> PhiPtr;
+
+IOobject io
 (
-    IOobject
-    (
-        "Phi",
-        runTime.timeName(),
-        mesh,
-        IOobject::READ_IF_PRESENT,
-        IOobject::NO_WRITE
-    ),
+    "Phi",
+    runTime.timeName(),
     mesh,
-    dimensionedScalar("Phi", dimLength*dimVelocity, 0),
-    PhiBCTypes
+    IOobject::MUST_READ,
+    IOobject::NO_WRITE
 );
 
+if (io.typeHeaderOk<volScalarField>())
+{
+    PhiPtr.reset(new volScalarField(io, mesh));
+}
+else
+{
+    // Cannot just use p.boundaryField().types() since does not initialise
+    // complex boundary types. Instead re-clone them from p.
+
+    io.readOpt() = IOobject::NO_READ;
+    PhiPtr.reset
+    (
+        new volScalarField
+        (
+            io,
+            mesh,
+            dimensionedScalar("Phi", dimLength*dimVelocity, 0),
+            p.boundaryField().types()
+        )
+    );
+
+    const volScalarField::GeometricBoundaryField& bp = p.boundaryField();
+    volScalarField::GeometricBoundaryField& bPhi = PhiPtr().boundaryField();
+
+    forAll(bp, patchI)
+    {
+        bPhi.set(patchI, bp[patchI].clone(PhiPtr().dimensionedInternalField()));
+    }
+}
+volScalarField& Phi = PhiPtr();
+
 label PhiRefCell = 0;
 scalar PhiRefValue = 0;
 setRefCell
@@ -119,5 +145,3 @@ setRefCell
     PhiRefValue
 );
 mesh.setFluxRequired(Phi.name());
-
-#include "createMRF.H"

applications/solvers/basic/potentialFoam/potentialFoam.C

@@ -28,9 +28,7 @@ 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.
+    Potential flow solver.
 
     \heading Solver details
     The potential flow solution is typically employed to generate initial fields
@@ -132,6 +130,7 @@ int main(int argc, char *argv[])
     pisoControl potentialFlow(mesh, "potentialFlow");
 
     #include "createFields.H"
+    #include "createMRF.H"
 
     // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 

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

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2015 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -87,16 +87,17 @@ Description
 
 int main(int argc, char *argv[])
 {
-    #include "postProcess.H"
-
     #include "setRootCase.H"
+
     #include "createTime.H"
     #include "createMesh.H"
-    #include "createControl.H"
+
+    pimpleControl pimple(mesh);
+
     #include "readCombustionProperties.H"
     #include "readGravitationalAcceleration.H"
     #include "createFields.H"
-    #include "createFieldRefs.H"
+    #include "createMRF.H"
     #include "createFvOptions.H"
     #include "initContinuityErrs.H"
     #include "createTimeControls.H"
@@ -112,7 +113,7 @@ int main(int argc, char *argv[])
 
     while (runTime.run())
     {
-        #include "readTimeControls.H"
+        #include "createTimeControls.H"
         #include "compressibleCourantNo.H"
         #include "setDeltaT.H"
 

applications/solvers/combustion/PDRFoam/PDRFoamAutoRefine.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2015 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -98,7 +98,7 @@ int main(int argc, char *argv[])
 
     while (runTime.run())
     {
-        #include "readTimeControls.H"
+        #include "createTimeControls.H"
         #include "compressibleCourantNo.H"
         #include "setDeltaT.H"
 
@@ -131,11 +131,11 @@ int main(int argc, char *argv[])
                 protectedCell = 0;
             }
 
-            forAll(betav, celli)
+            forAll(betav, cellI)
             {
-                if (betav[celli] < 0.99)
+                if (betav[cellI] < 0.99)
                 {
-                    protectedCell[celli] = 1;
+                    protectedCell[cellI] = 1;
                 }
             }
 

applications/solvers/combustion/PDRFoam/PDRModels/XiEqModels/basicXiSubXiEq/basicXiSubXiEq.C

@@ -105,7 +105,7 @@ Foam::tmp<Foam::volScalarField> Foam::XiEqModels::basicSubGrid::XiEq() const
         mesh,
         dimensionedScalar("zero", Nv.dimensions(), 0.0)
     );
-    N.primitiveFieldRef() = Nv.primitiveField()*Cw;
+    N.internalField() = Nv.internalField()*Cw;
 
     volSymmTensorField ns
     (
@@ -125,7 +125,7 @@ Foam::tmp<Foam::volScalarField> Foam::XiEqModels::basicSubGrid::XiEq() const
             Zero
         )
     );
-    ns.primitiveFieldRef() = nsv.primitiveField()*Cw;
+    ns.internalField() = nsv.internalField()*Cw;
 
     volScalarField n(max(N - (Uhat & ns & Uhat), scalar(1e-4)));
     volScalarField b((Uhat & B_ & Uhat)/sqrt(n));

applications/solvers/combustion/PDRFoam/PDRModels/XiGModels/basicXiSubG/basicXiSubG.C

@@ -73,7 +73,7 @@ Foam::tmp<Foam::volScalarField> Foam::XiGModels::basicSubGrid::G() const
     volScalarField& Gtot = tGtot.ref();
 
     const scalarField Cw = pow(Su_.mesh().V(), 2.0/3.0);
-    scalarField N(Nv.primitiveField()*Cw);
+    scalarField N(Nv.internalField()*Cw);
 
     forAll(N, celli)
     {

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

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

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

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

applications/solvers/combustion/PDRFoam/PDRModels/turbulence/PDRkEpsilon/PDRkEpsilon.C

@@ -130,7 +130,7 @@ void PDRkEpsilon::correct()
     tgradU.clear();
 
     // Update espsilon and G at the wall
-    epsilon_.boundaryFieldRef().updateCoeffs();
+    epsilon_.boundaryField().updateCoeffs();
 
     // Add the blockage generation term so that it is included consistently
     // in both the k and epsilon equations
@@ -163,7 +163,7 @@ void PDRkEpsilon::correct()
 
     epsEqn.ref().relax();
 
-    epsEqn.ref().boundaryManipulate(epsilon_.boundaryFieldRef());
+    epsEqn.ref().boundaryManipulate(epsilon_.boundaryField());
 
     solve(epsEqn);
     bound(epsilon_, epsilonMin_);

applications/solvers/combustion/PDRFoam/StCourantNo.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -36,8 +36,8 @@ Description
     {
        scalarField sumPhi
         (
-            fvc::surfaceSum(mag(phiSt))().primitiveField()
-          / rho.primitiveField()
+            fvc::surfaceSum(mag(phiSt))().internalField()
+          / rho.internalField()
         );
 
         StCoNum = 0.5*gMax(sumPhi/mesh.V().field())*runTime.deltaTValue();

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

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2015 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -74,7 +74,7 @@ Foam::tmp<Foam::volScalarField> Foam::XiEqModels::Gulder::XiEq() const
 
     if (subGridSchelkin_)
     {
-        up.primitiveFieldRef() += calculateSchelkinEffect(uPrimeCoef_);
+        up.internalField() += calculateSchelkinEffect(uPrimeCoef_);
     }
 
     volScalarField tauEta(sqrt(mag(thermo_.muu()/(thermo_.rhou()*epsilon))));

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

@@ -82,7 +82,7 @@ Foam::tmp<Foam::volScalarField> Foam::XiEqModels::SCOPEXiEq::XiEq() const
     volScalarField up(sqrt((2.0/3.0)*k));
     if (subGridSchelkin_)
     {
-        up.primitiveFieldRef() += calculateSchelkinEffect(uPrimeCoef_);
+        up.internalField() += calculateSchelkinEffect(uPrimeCoef_);
     }
 
     volScalarField l(lCoef_*sqrt(3.0/2.0)*up*k/epsilon);
@@ -119,11 +119,9 @@ Foam::tmp<Foam::volScalarField> Foam::XiEqModels::SCOPEXiEq::XiEq() const
         }
     }
 
-    volScalarField::Boundary& xieqBf = xieq.boundaryFieldRef();
-
     forAll(xieq.boundaryField(), patchi)
     {
-        scalarField& xieqp = xieqBf[patchi];
+        scalarField& xieqp = xieq.boundaryField()[patchi];
         const scalarField& Kp = K.boundaryField()[patchi];
         const scalarField& Map = Ma.boundaryField()[patchi];
         const scalarField& upBySup = upBySu.boundaryField()[patchi];

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;
 }
@@ -114,7 +114,7 @@ Foam::XiEqModel::calculateSchelkinEffect(const scalar uPrimeCoef) const
         )
     );
     volScalarField& N = tN.ref();
-    N.primitiveFieldRef() = Nv.primitiveField()*pow(mesh.V(), 2.0/3.0);
+    N.internalField() = Nv.internalField()*pow(mesh.V(), 2.0/3.0);
 
     volSymmTensorField ns
     (
@@ -134,7 +134,7 @@ Foam::XiEqModel::calculateSchelkinEffect(const scalar uPrimeCoef) const
             Zero
         )
     );
-    ns.primitiveFieldRef() = nsv.primitiveField()*pow(mesh.V(), 2.0/3.0);
+    ns.internalField() = nsv.internalField()*pow(mesh.V(), 2.0/3.0);
 
     const volVectorField Uhat
     (
@@ -150,7 +150,7 @@ Foam::XiEqModel::calculateSchelkinEffect(const scalar uPrimeCoef) const
     const scalarField deltaUp(upLocal*(max(scalar(1.0), pow(nr, 0.5)) - 1.0));
 
     // Re use tN
-    N.primitiveFieldRef() = upLocal*(max(scalar(1.0), pow(nr, 0.5)) - 1.0);
+    N.internalField() = upLocal*(max(scalar(1.0), pow(nr, 0.5)) - 1.0);
 
     return tN;
 }

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

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

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

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

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

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

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

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

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

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

applications/solvers/combustion/PDRFoam/createFieldRefs.H

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

applications/solvers/combustion/PDRFoam/createFields.H

@@ -23,6 +23,7 @@ volScalarField rho
 );
 
 volScalarField& p = thermo.p();
+const volScalarField& psi = thermo.psi();
 
 volScalarField& b = composition.Y("b");
 Info<< "min(b) = " << min(b).value() << endl;
@@ -58,6 +59,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
 (
@@ -209,9 +227,3 @@ fields.add(b);
 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"))),
@@ -266,11 +266,9 @@ Foam::tmp<Foam::volScalarField> Foam::laminarFlameSpeedModels::SCOPE::Su0pTphi
         Su0[celli] = Su0pTphi(p[celli], Tu[celli], phi);
     }
 
-    volScalarField::Boundary& Su0Bf = Su0.boundaryFieldRef();
-
-    forAll(Su0Bf, patchi)
+    forAll(Su0.boundaryField(), patchi)
     {
-        scalarField& Su0p = Su0Bf[patchi];
+        scalarField& Su0p = Su0.boundaryField()[patchi];
         const scalarField& pp = p.boundaryField()[patchi];
         const scalarField& Tup = Tu.boundaryField()[patchi];
 
@@ -315,11 +313,9 @@ Foam::tmp<Foam::volScalarField> Foam::laminarFlameSpeedModels::SCOPE::Su0pTphi
         Su0[celli] = Su0pTphi(p[celli], Tu[celli], phi[celli]);
     }
 
-    volScalarField::Boundary& Su0Bf = Su0.boundaryFieldRef();
-
-    forAll(Su0Bf, patchi)
+    forAll(Su0.boundaryField(), patchi)
     {
-        scalarField& Su0p = Su0Bf[patchi];
+        scalarField& Su0p = Su0.boundaryField()[patchi];
         const scalarField& pp = p.boundaryField()[patchi];
         const scalarField& Tup = Tu.boundaryField()[patchi];
         const scalarField& phip = phi.boundaryField()[patchi];
@@ -369,11 +365,9 @@ Foam::tmp<Foam::volScalarField> Foam::laminarFlameSpeedModels::SCOPE::Ma
         ma[celli] = Ma(phi[celli]);
     }
 
-    volScalarField::Boundary& maBf = ma.boundaryFieldRef();
-
-    forAll(maBf, patchi)
+    forAll(ma.boundaryField(), patchi)
     {
-        scalarField& map = maBf[patchi];
+        scalarField& map = ma.boundaryField()[patchi];
         const scalarField& phip = phi.boundaryField()[patchi];
 
         forAll(map, facei)

applications/solvers/combustion/XiFoam/XiDyMFoam/XiDyMFoam.C

@@ -3,7 +3,7 @@
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
     \\  /    A nd           | Copyright (C) 2011-2015 OpenFOAM Foundation
-     \\/     M anipulation  | Copyright (C) 2015-2016 OpenCFD Ltd.
+     \\/     M anipulation  | Copyright (C) 2015 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -67,26 +67,25 @@ Description
 
 int main(int argc, char *argv[])
 {
-    #include "postProcess.H"
-
     #include "setRootCase.H"
+
     #include "createTime.H"
     #include "createDynamicFvMesh.H"
-    #include "createControl.H"
+    #include "initContinuityErrs.H"
+
+    pimpleControl pimple(mesh);
+
     #include "readCombustionProperties.H"
     #include "readGravitationalAcceleration.H"
     #include "createFields.H"
-    #include "createFieldRefs.H"
+    #include "createMRF.H"
     #include "createFvOptions.H"
-    #include "initContinuityErrs.H"
     #include "createRhoUf.H"
     #include "createControls.H"
     #include "initContinuityErrs.H"
     #include "compressibleCourantNo.H"
     #include "setInitialDeltaT.H"
 
-    turbulence->validate();
-
     // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
     Info<< "\nStarting time loop\n" << endl;

applications/solvers/combustion/XiFoam/XiDyMFoam/readControls.H

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

applications/solvers/combustion/XiFoam/XiFoam.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2015 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -65,15 +65,16 @@ Description
 
 int main(int argc, char *argv[])
 {
-    #include "postProcess.H"
-
     #include "setRootCase.H"
+
     #include "createTime.H"
     #include "createMesh.H"
-    #include "createControl.H"
+
+    pimpleControl pimple(mesh);
+
     #include "readCombustionProperties.H"
     #include "createFields.H"
-    #include "createFieldRefs.H"
+    #include "createMRF.H"
     #include "createFvOptions.H"
     #include "initContinuityErrs.H"
     #include "createTimeControls.H"
@@ -88,7 +89,7 @@ int main(int argc, char *argv[])
 
     while (runTime.run())
     {
-        #include "readTimeControls.H"
+        #include "createTimeControls.H"
         #include "compressibleCourantNo.H"
         #include "setDeltaT.H"
 

applications/solvers/combustion/XiFoam/createFieldRefs.H

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

applications/solvers/combustion/XiFoam/createFields.H

@@ -23,6 +23,7 @@ volScalarField rho
 );
 
 volScalarField& p = thermo.p();
+const volScalarField& psi = thermo.psi();
 
 volScalarField& b = composition.Y("b");
 Info<< "min(b) = " << min(b).value() << endl;
@@ -58,6 +59,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
@@ -123,9 +139,3 @@ if (composition.contains("ft"))
 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

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2015 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -28,10 +28,12 @@ Group
     grpCombustionSolvers
 
 Description
-    Solver for chemistry problems, designed for use on single cell cases to
-    provide comparison against other chemistry solvers, that uses a single cell
-    mesh, and fields created from the initial conditions.
+    Solver for chemistry problems designed for use on single cell cases to
+    provide comparison against other chemistry solvers
 
+    Note:
+    - single cell mesh created on-the-fly
+    - fields created on the fly from the initial conditions
 
 \*---------------------------------------------------------------------------*/
 
@@ -50,15 +52,10 @@ int main(int argc, char *argv[])
 {
     argList::noParallel();
 
-    #define CREATE_MESH createSingleCellMesh.H
-    #define NO_CONTROL
-    #include "postProcess.H"
-
     #include "setRootCase.H"
     #include "createTime.H"
     #include "createSingleCellMesh.H"
     #include "createFields.H"
-    #include "createFieldRefs.H"
     #include "readInitialConditions.H"
     #include "createControls.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/createFieldRefs.H

@@ -1,4 +0,0 @@
-scalar dtChem = refCast<const psiChemistryModel>(chemistry).deltaTChem()[0];
-basicMultiComponentMixture& composition = thermo.composition();
-PtrList<volScalarField>& Y = composition.Y();
-volScalarField& p = thermo.p();

applications/solvers/combustion/chemFoam/createFields.H

@@ -12,7 +12,7 @@
         (
             "initialConditions",
             runTime.constant(),
-            mesh,
+            runTime,
             IOobject::MUST_READ_IF_MODIFIED,
             IOobject::NO_WRITE
         )
@@ -27,30 +27,36 @@
     autoPtr<psiChemistryModel> pChemistry(psiChemistryModel::New(mesh));
 
     psiChemistryModel& chemistry = pChemistry();
+    scalar dtChem = refCast<const psiChemistryModel>(chemistry).deltaTChem()[0];
 
     psiReactionThermo& thermo = chemistry.thermo();
     thermo.validate(args.executable(), "h");
 
+    basicMultiComponentMixture& composition = thermo.composition();
+    PtrList<volScalarField>& Y = composition.Y();
+
     volScalarField rho
     (
         IOobject
         (
             "rho",
             runTime.timeName(),
-            mesh,
+            runTime,
             IOobject::NO_READ,
             IOobject::AUTO_WRITE
         ),
         thermo.rho()
     );
 
+    volScalarField& p = thermo.p();
+
     volScalarField Rspecific
     (
         IOobject
         (
             "Rspecific",
             runTime.timeName(),
-            mesh,
+            runTime,
             IOobject::NO_READ,
             IOobject::AUTO_WRITE
         ),
@@ -69,7 +75,7 @@
         (
             "U",
             runTime.timeName(),
-            mesh,
+            runTime,
             IOobject::NO_READ,
             IOobject::NO_WRITE
         ),

applications/solvers/combustion/chemFoam/output.H

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

applications/solvers/combustion/chemFoam/setDeltaT.H

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

applications/solvers/combustion/chemFoam/solveChemistry.H

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

applications/solvers/combustion/coldEngineFoam/Make/options

@@ -1,5 +1,4 @@
 EXE_INC = \
-    -I. \
     -I../engineFoam \
     -I../XiFoam \
     -I../../compressible/rhoPimpleFoam \

applications/solvers/combustion/coldEngineFoam/coldEngineFoam.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2015 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -45,16 +45,15 @@ Description
 
 int main(int argc, char *argv[])
 {
-    #define CREATE_TIME createEngineTime.H
-    #define CREATE_MESH createEngineMesh.H
-    #include "postProcess.H"
-
     #include "setRootCase.H"
+
     #include "createEngineTime.H"
     #include "createEngineMesh.H"
-    #include "createControl.H"
+
+    pimpleControl pimple(mesh);
+
     #include "createFields.H"
-    #include "createFieldRefs.H"
+    #include "createMRF.H"
     #include "createFvOptions.H"
     #include "createRhoUf.H"
     #include "initContinuityErrs.H"

applications/solvers/combustion/coldEngineFoam/createFieldRefs.H

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

applications/solvers/combustion/coldEngineFoam/createFields.H

@@ -1,58 +1,71 @@
-Info<< "Reading thermophysical properties\n" << endl;
+    Info<< "Reading thermophysical properties\n" << endl;
 
-autoPtr<psiThermo> pThermo
-(
-    psiThermo::New(mesh)
-);
-psiThermo& thermo = pThermo();
-thermo.validate(args.executable(), "h", "e");
-
-volScalarField rho
-(
-    IOobject
+    autoPtr<psiThermo> pThermo
     (
-        "rho",
-        runTime.timeName(),
+        psiThermo::New(mesh)
+    );
+    psiThermo& thermo = pThermo();
+    thermo.validate(args.executable(), "h", "e");
+
+    volScalarField rho
+    (
+        IOobject
+        (
+            "rho",
+            runTime.timeName(),
+            mesh,
+            IOobject::NO_READ,
+            IOobject::AUTO_WRITE
+        ),
+        thermo.rho()
+    );
+
+    volScalarField& p = thermo.p();
+    const volScalarField& psi = thermo.psi();
+    const volScalarField& T = thermo.T();
+
+
+    Info<< "\nReading field U\n" << endl;
+    volVectorField U
+    (
+        IOobject
+        (
+            "U",
+            runTime.timeName(),
+            mesh,
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh
+    );
+
+    #include "compressibleCreatePhi.H"
+
+
+    Info<< "Creating turbulence model\n" << endl;
+    autoPtr<compressible::turbulenceModel> turbulence
+    (
+        compressible::turbulenceModel::New
+        (
+            rho,
+            U,
+            phi,
+            thermo
+        )
+    );
+
+    Info<< "Creating field dpdt\n" << endl;
+    volScalarField dpdt
+    (
+        IOobject
+        (
+            "dpdt",
+            runTime.timeName(),
+            mesh
+        ),
         mesh,
-        IOobject::NO_READ,
-        IOobject::AUTO_WRITE
-    ),
-    thermo.rho()
-);
+        dimensionedScalar("dpdt", p.dimensions()/dimTime, 0)
+    );
 
-volScalarField& p = thermo.p();
-
-Info<< "\nReading field U\n" << endl;
-volVectorField U
-(
-    IOobject
-    (
-        "U",
-        runTime.timeName(),
-        mesh,
-        IOobject::MUST_READ,
-        IOobject::AUTO_WRITE
-    ),
-    mesh
-);
-
-#include "compressibleCreatePhi.H"
-
-
-Info<< "Creating turbulence model\n" << endl;
-autoPtr<compressible::turbulenceModel> turbulence
-(
-    compressible::turbulenceModel::New
-    (
-        rho,
-        U,
-        phi,
-        thermo
-    )
-);
-
-#include "createDpdt.H"
-
-#include "createK.H"
-
-#include "createMRF.H"
+    Info<< "Creating field kinetic energy K\n" << endl;
+    volScalarField K("K", 0.5*magSqr(U));

applications/solvers/combustion/engineFoam/engineFoam.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2015 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -68,17 +68,16 @@ Description
 
 int main(int argc, char *argv[])
 {
-    #define CREATE_TIME createEngineTime.H
-    #define CREATE_MESH createEngineMesh.H
-    #include "postProcess.H"
-
     #include "setRootCase.H"
+
     #include "createEngineTime.H"
     #include "createEngineMesh.H"
-    #include "createControl.H"
+
+    pimpleControl pimple(mesh);
+
     #include "readCombustionProperties.H"
     #include "createFields.H"
-    #include "createFieldRefs.H"
+    #include "createMRF.H"
     #include "createFvOptions.H"
     #include "createRhoUf.H"
     #include "initContinuityErrs.H"

applications/solvers/combustion/fireFoam/Make/options

@@ -13,11 +13,15 @@ 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 \
     -I$(LIB_SRC)/thermophysicalModels/radiation/lnInclude \
+    -I$(LIB_SRC)/combustionModels/lnInclude \
     -I$(LIB_SRC)/regionModels/regionModel/lnInclude \
     -I$(LIB_SRC)/regionModels/surfaceFilmModels/lnInclude \
     -I$(LIB_SRC)/regionModels/pyrolysisModels/lnInclude \
@@ -25,6 +29,7 @@ EXE_INC = \
     -I$(LIB_SRC)/lagrangian/intermediate/lnInclude \
     -I$(LIB_SRC)/ODE/lnInclude
 
+
 EXE_LIBS = \
     -lfiniteVolume \
     -lfvOptions \
@@ -35,6 +40,9 @@ EXE_LIBS = \
     -lspecie \
     -lcompressibleTransportModels \
     -lfluidThermophysicalModels \
+    -lsolidProperties \
+    -lsolidMixtureProperties \
+    -lthermophysicalFunctions \
     -lreactionThermophysicalModels \
     -lSLGThermo \
     -lchemistryModel \

applications/solvers/combustion/fireFoam/YEEqn.H

@@ -9,13 +9,15 @@ tmp<fv::convectionScheme<scalar>> mvConvection
     )
 );
 {
+    radiation->correct();
     combustion->correct();
-    Qdot = combustion->Qdot();
+    dQ = combustion->dQ();
+    label inertIndex = -1;
     volScalarField Yt(0.0*Y[0]);
 
     forAll(Y, i)
     {
-        if (i != inertIndex && composition.active(i))
+        if (Y[i].name() != inertSpecie)
         {
             volScalarField& Yi = Y[i];
 
@@ -42,6 +44,10 @@ tmp<fv::convectionScheme<scalar>> mvConvection
             Yi.max(0.0);
             Yt += Yi;
         }
+        else
+        {
+            inertIndex = i;
+        }
     }
 
     Y[inertIndex] = scalar(1) - Yt;
@@ -67,8 +73,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 +0,0 @@
-const volScalarField& psi = thermo.psi();
-const volScalarField& T = thermo.T();
-filmModelType& surfaceFilm = tsurfaceFilm();
-const label inertIndex(composition.species()[inertSpecie]);

applications/solvers/combustion/fireFoam/createFields.H

@@ -19,13 +19,6 @@ basicMultiComponentMixture& composition = thermo.composition();
 PtrList<volScalarField>& Y = composition.Y();
 
 const word inertSpecie(thermo.lookup("inertSpecie"));
-if (!composition.species().found(inertSpecie))
-{
-    FatalIOErrorIn(args.executable().c_str(), thermo)
-        << "Inert specie " << inertSpecie << " not found in available species "
-        << composition.species()
-        << exit(FatalIOError);
-}
 
 Info<< "Creating field rho\n" << endl;
 volScalarField rho
@@ -42,6 +35,8 @@ volScalarField rho
 );
 
 volScalarField& p = thermo.p();
+const volScalarField& T = thermo.T();
+const volScalarField& psi = thermo.psi();
 
 Info<< "\nReading field U\n" << endl;
 volVectorField U
@@ -59,9 +54,6 @@ volVectorField U
 
 #include "compressibleCreatePhi.H"
 
-#include "createMRF.H"
-
-
 Info<< "Creating turbulence model\n" << endl;
 autoPtr<compressible::turbulenceModel> turbulence
 (
@@ -77,11 +69,42 @@ autoPtr<compressible::turbulenceModel> turbulence
 // Set the turbulence into the combustion model
 combustion->setTurbulence(turbulence());
 
+volScalarField dQ
+(
+    IOobject
+    (
+        "dQ",
+        runTime.timeName(),
+        mesh,
+        IOobject::NO_READ,
+        IOobject::AUTO_WRITE
+    ),
+    mesh,
+    dimensionedScalar("dQ", dimEnergy/dimTime, 0.0)
+);
+
+
+Info<< "Creating field dpdt\n" << endl;
+volScalarField dpdt
+(
+    IOobject
+    (
+        "dpdt",
+        runTime.timeName(),
+        mesh
+    ),
+    mesh,
+    dimensionedScalar("dpdt", p.dimensions()/dimTime, 0)
+);
+
+Info<< "Creating field kinetic energy K\n" << endl;
+volScalarField K("K", 0.5*magSqr(U));
+
 
 #include "readGravitationalAcceleration.H"
 #include "readhRef.H"
 #include "gh.H"
-#include "readpRef.H"
+
 
 volScalarField p_rgh
 (
@@ -96,11 +119,11 @@ volScalarField p_rgh
     mesh
 );
 
+// Force p_rgh to be consistent with p
+p_rgh = p - rho*gh;
+
 mesh.setFluxRequired(p_rgh.name());
 
-#include "phrghEqn.H"
-
-
 multivariateSurfaceInterpolationScheme<scalar>::fieldTable fields;
 
 forAll(Y, i)
@@ -125,31 +148,3 @@ Switch solvePrimaryRegion
 (
     additionalControlsDict.lookup("solvePrimaryRegion")
 );
-
-Switch solvePyrolysisRegion
-(
-    additionalControlsDict.lookupOrDefault<bool>("solvePyrolysisRegion", true)
-);
-
-volScalarField Qdot
-(
-    IOobject
-    (
-        "Qdot",
-        runTime.timeName(),
-        mesh,
-        IOobject::READ_IF_PRESENT,
-        IOobject::AUTO_WRITE
-    ),
-    mesh,
-    dimensionedScalar("Qdot", dimEnergy/dimVolume/dimTime, 0.0)
-);
-
-#include "createDpdt.H"
-
-#include "createK.H"
-
-#include "createClouds.H"
-#include "createSurfaceFilmModel.H"
-#include "createPyrolysisModel.H"
-#include "createRadiationModel.H"

applications/solvers/combustion/fireFoam/createSurfaceFilmModel.H

@@ -3,3 +3,4 @@ Info<< "\nConstructing surface film model" << endl;
 typedef regionModels::surfaceFilmModels::surfaceFilmModel filmModelType;
 
 autoPtr<filmModelType> tsurfaceFilm(filmModelType::New(mesh, g));
+filmModelType& surfaceFilm = tsurfaceFilm();

applications/solvers/combustion/fireFoam/fireFoam.C

@@ -28,8 +28,8 @@ Group
     grpCombustionSolvers
 
 Description
-    Transient solver for fires and turbulent diffusion flames with reacting
-    particle clouds, surface film and pyrolysis modelling.
+    Transient PIMPLE solver for fires and turbulent diffusion flames with
+    reacting Lagrangian parcels, surface film and pyrolysis modelling.
 
 \*---------------------------------------------------------------------------*/
 
@@ -49,15 +49,20 @@ Description
 
 int main(int argc, char *argv[])
 {
-    #include "postProcess.H"
-
     #include "setRootCase.H"
+
     #include "createTime.H"
     #include "createMesh.H"
-    #include "createControl.H"
+
+    pimpleControl pimple(mesh);
+
     #include "createFields.H"
-    #include "createFieldRefs.H"
+    #include "createMRF.H"
     #include "createFvOptions.H"
+    #include "createClouds.H"
+    #include "createSurfaceFilmModel.H"
+    #include "createPyrolysisModel.H"
+    #include "createRadiationModel.H"
     #include "initContinuityErrs.H"
     #include "createTimeControls.H"
     #include "compressibleCourantNo.H"
@@ -72,7 +77,7 @@ int main(int argc, char *argv[])
 
     while (runTime.run())
     {
-        #include "readTimeControls.H"
+        #include "createTimeControls.H"
         #include "compressibleCourantNo.H"
         #include "solidRegionDiffusionNo.H"
         #include "setMultiRegionDeltaT.H"
@@ -86,10 +91,7 @@ int main(int argc, char *argv[])
 
         surfaceFilm.evolve();
 
-        if(solvePyrolysisRegion)
-        {
-            pyrolysis.evolve();
-        }
+        pyrolysis.evolve();
 
         if (solvePrimaryRegion)
         {

applications/solvers/combustion/fireFoam/pEqn.H

@@ -4,7 +4,7 @@ volScalarField rAU(1.0/UEqn.A());
 surfaceScalarField rhorAUf("rhorAUf", fvc::interpolate(rho*rAU));
 volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p));
 
-surfaceScalarField phig("phig", -rhorAUf*ghf*fvc::snGrad(rho)*mesh.magSf());
+surfaceScalarField phig(-rhorAUf*ghf*fvc::snGrad(rho)*mesh.magSf());
 
 surfaceScalarField phiHbyA
 (
@@ -25,10 +25,9 @@ while (pimple.correctNonOrthogonal())
 {
     fvScalarMatrix p_rghEqn
     (
-        fvm::ddt(psi, p_rgh)
-      + fvc::ddt(psi, rho)*gh
-      + fvc::ddt(psi)*pRef
+        fvc::ddt(psi, rho)*gh
       + fvc::div(phiHbyA)
+      + fvm::ddt(psi, p_rgh)
       - fvm::laplacian(rhorAUf, p_rgh)
      ==
         parcels.Srho()
@@ -47,7 +46,7 @@ while (pimple.correctNonOrthogonal())
     }
 }
 
-p = p_rgh + rho*gh + pRef;
+p = p_rgh + rho*gh;
 
 #include "rhoEqn.H"
 #include "compressibleContinuityErrs.H"

applications/solvers/combustion/fireFoam/phrghEqn.H

@@ -1,62 +0,0 @@
-if (pimple.dict().lookupOrDefault<bool>("hydrostaticInitialization", false))
-{
-    volScalarField& ph_rgh = regIOobject::store
-    (
-        new volScalarField
-        (
-            IOobject
-            (
-                "ph_rgh",
-                "0",
-                mesh,
-                IOobject::MUST_READ,
-                IOobject::NO_WRITE
-            ),
-            mesh
-        )
-    );
-
-    if (equal(runTime.value(), 0))
-    {
-        p = ph_rgh + rho*gh + pRef;
-        thermo.correct();
-        rho = thermo.rho();
-
-        label nCorr
-        (
-            pimple.dict().lookupOrDefault<label>("nHydrostaticCorrectors", 5)
-        );
-
-        for (label i=0; i<nCorr; i++)
-        {
-            surfaceScalarField rhof("rhof", fvc::interpolate(rho));
-
-            surfaceScalarField phig
-            (
-                "phig",
-               -rhof*ghf*fvc::snGrad(rho)*mesh.magSf()
-            );
-
-            // Update the pressure BCs to ensure flux consistency
-            constrainPressure(ph_rgh, rho, U, phig, rhof);
-
-            fvScalarMatrix ph_rghEqn
-            (
-                fvm::laplacian(rhof, ph_rgh) == fvc::div(phig)
-            );
-
-            ph_rghEqn.solve();
-
-            p = ph_rgh + rho*gh + pRef;
-            thermo.correct();
-            rho = thermo.rho();
-
-            Info<< "Hydrostatic pressure variation "
-                << (max(ph_rgh) - min(ph_rgh)).value() << endl;
-        }
-
-        ph_rgh.write();
-
-        p_rgh = ph_rgh;
-    }
-}

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,12 +11,13 @@ tmp<fv::convectionScheme<scalar>> mvConvection
 
 {
     reaction->correct();
-    Qdot = reaction->Qdot();
+    dQ = reaction->dQ();
+    label inertIndex = -1;
     volScalarField Yt(0.0*Y[0]);
 
     forAll(Y, i)
     {
-        if (i != inertIndex && composition.active(i))
+        if (Y[i].name() != inertSpecie)
         {
             volScalarField& Yi = Y[i];
 
@@ -41,6 +42,10 @@ tmp<fv::convectionScheme<scalar>> mvConvection
             Yi.max(0.0);
             Yt += Yi;
         }
+        else
+        {
+            inertIndex = i;
+        }
     }
 
     Y[inertIndex] = scalar(1) - Yt;

applications/solvers/combustion/reactingFoam/createFieldRefs.H

@@ -1,3 +0,0 @@
-const volScalarField& psi = thermo.psi();
-const volScalarField& T = thermo.T();
-const label inertIndex(composition.species()[inertSpecie]);

applications/solvers/combustion/reactingFoam/createFields.H

@@ -1,5 +1,3 @@
-#include "createRDeltaT.H"
-
 Info<< "Creating reaction model\n" << endl;
 
 autoPtr<combustionModels::psiCombustionModel> reaction
@@ -13,14 +11,7 @@ thermo.validate(args.executable(), "h", "e");
 basicMultiComponentMixture& composition = thermo.composition();
 PtrList<volScalarField>& Y = composition.Y();
 
-const word inertSpecie(thermo.lookup("inertSpecie"));
-if (!composition.species().found(inertSpecie))
-{
-    FatalIOErrorIn(args.executable().c_str(), thermo)
-        << "Inert specie " << inertSpecie << " not found in available species "
-        << composition.species()
-        << exit(FatalIOError);
-}
+word inertSpecie(thermo.lookup("inertSpecie"));
 
 volScalarField rho
 (
@@ -49,10 +40,32 @@ volVectorField U
 
 
 volScalarField& p = thermo.p();
+const volScalarField& psi = thermo.psi();
+const volScalarField& T = thermo.T();
 
 #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 +84,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 +110,16 @@ 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-2015 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"
@@ -46,16 +45,17 @@ Description
 
 int main(int argc, char *argv[])
 {
-    #include "postProcess.H"
-
     #include "setRootCase.H"
     #include "createTime.H"
     #include "createMesh.H"
-    #include "createControl.H"
+
+    pimpleControl pimple(mesh);
+
     #include "createTimeControls.H"
+    #include "createRDeltaT.H"
     #include "initContinuityErrs.H"
     #include "createFields.H"
-    #include "createFieldRefs.H"
+    #include "createMRF.H"
     #include "createFvOptions.H"
 
     turbulence->validate();
@@ -115,8 +115,6 @@ int main(int argc, char *argv[])
             }
         }
 
-        rho = thermo.rho();
-
         runTime.write();
 
         Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"

applications/solvers/combustion/reactingFoam/rhoReactingBuoyantFoam/Make/options

@@ -1,5 +1,4 @@
 EXE_INC = \
-    -I. \
     -I$(FOAM_SOLVERS)/combustion/reactingFoam \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
     -I$(LIB_SRC)/meshTools/lnInclude \

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

@@ -1,3 +0,0 @@
-const volScalarField& psi = thermo.psi();
-const volScalarField& T = thermo.T();
-const label inertIndex(composition.species()[inertSpecie]);

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

@@ -1,5 +1,3 @@
-#include "createRDeltaT.H"
-
 Info<< "Creating reaction model\n" << endl;
 
 autoPtr<combustionModels::rhoCombustionModel> reaction
@@ -13,14 +11,7 @@ thermo.validate(args.executable(), "h", "e");
 basicMultiComponentMixture& composition = thermo.composition();
 PtrList<volScalarField>& Y = composition.Y();
 
-const word inertSpecie(thermo.lookup("inertSpecie"));
-if (!composition.species().found(inertSpecie))
-{
-    FatalIOErrorIn(args.executable().c_str(), thermo)
-        << "Inert specie " << inertSpecie << " not found in available species "
-        << composition.species()
-        << exit(FatalIOError);
-}
+word inertSpecie(thermo.lookup("inertSpecie"));
 
 volScalarField rho
 (
@@ -49,6 +40,8 @@ volVectorField U
 
 
 volScalarField& p = thermo.p();
+const volScalarField& psi = thermo.psi();
+const volScalarField& T = thermo.T();
 
 
 #include "compressibleCreatePhi.H"
@@ -92,6 +85,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 +110,16 @@ 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/rhoReactingBuoyantFoam/rhoReactingBuoyantFoam.C

@@ -28,8 +28,8 @@ Group
     grpCombustionSolvers
 
 Description
-    Solver for combustion with chemical reactions using a density based
-    thermodynamics package with enhanced buoyancy treatment.
+    Solver for combustion with chemical reactions using density based
+    thermodynamics package, using enahanced buoyancy treatment.
 
 \*---------------------------------------------------------------------------*/
 
@@ -46,16 +46,17 @@ Description
 
 int main(int argc, char *argv[])
 {
-    #include "postProcess.H"
-
     #include "setRootCase.H"
     #include "createTime.H"
     #include "createMesh.H"
-    #include "createControl.H"
+
+    pimpleControl pimple(mesh);
+
     #include "createTimeControls.H"
+    #include "createRDeltaT.H"
     #include "initContinuityErrs.H"
     #include "createFields.H"
-    #include "createFieldRefs.H"
+    #include "createMRF.H"
     #include "createFvOptions.H"
 
     turbulence->validate();

applications/solvers/combustion/reactingFoam/rhoReactingFoam/Make/options

@@ -1,5 +1,4 @@
 EXE_INC = \
-    -I. \
     -I$(FOAM_SOLVERS)/combustion/reactingFoam \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
     -I$(LIB_SRC)/meshTools/lnInclude \

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

@@ -1,3 +0,0 @@
-const volScalarField& psi = thermo.psi();
-const volScalarField& T = thermo.T();
-const label inertIndex(composition.species()[inertSpecie]);

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

@@ -1,5 +1,3 @@
-#include "createRDeltaT.H"
-
 Info<< "Creating reaction model\n" << endl;
 
 autoPtr<combustionModels::rhoCombustionModel> reaction
@@ -13,14 +11,7 @@ thermo.validate(args.executable(), "h", "e");
 basicMultiComponentMixture& composition = thermo.composition();
 PtrList<volScalarField>& Y = composition.Y();
 
-const word inertSpecie(thermo.lookup("inertSpecie"));
-if (!composition.species().found(inertSpecie))
-{
-    FatalIOErrorIn(args.executable().c_str(), thermo)
-        << "Inert specie " << inertSpecie << " not found in available species "
-        << composition.species()
-        << exit(FatalIOError);
-}
+word inertSpecie(thermo.lookup("inertSpecie"));
 
 volScalarField rho
 (
@@ -49,11 +40,11 @@ volVectorField U
 
 
 volScalarField& p = thermo.p();
+const volScalarField& psi = thermo.psi();
+const volScalarField& T = thermo.T();
 
 #include "compressibleCreatePhi.H"
 
-pressureControl pressureControl(p, rho, pimple.dict(), false);
-
 mesh.setFluxRequired(p.name());
 
 
@@ -72,6 +63,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 +89,16 @@ 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-2015 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"
@@ -47,16 +46,17 @@ Description
 
 int main(int argc, char *argv[])
 {
-    #include "postProcess.H"
-
     #include "setRootCase.H"
     #include "createTime.H"
     #include "createMesh.H"
-    #include "createControl.H"
+
+    pimpleControl pimple(mesh);
+
     #include "createTimeControls.H"
+    #include "createRDeltaT.H"
     #include "initContinuityErrs.H"
     #include "createFields.H"
-    #include "createFieldRefs.H"
+    #include "createMRF.H"
     #include "createFvOptions.H"
 
     turbulence->validate();
@@ -101,14 +101,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;
 
@@ -61,99 +58,44 @@ License
 
     // Flow time scale
     {
-        rDeltaT.ref() =
+        rDeltaT.dimensionedInternalField() =
         (
-            fvc::surfaceSum(mag(phi))()()
-           /((2*maxCo)*mesh.V()*rho())
+            fvc::surfaceSum(mag(phi))().dimensionedInternalField()
+           /((2*maxCo)*mesh.V()*rho.dimensionedInternalField())
         );
 
         // Limit the largest time scale
         rDeltaT.max(1/maxDeltaT);
 
         Info<< "    Flow        = "
-            << 1/gMax(rDeltaT.primitiveField()) << ", "
-            << 1/gMin(rDeltaT.primitiveField()) << endl;
+            << gMin(1/rDeltaT.internalField()) << ", "
+            << gMax(1/rDeltaT.internalField()) << endl;
     }
 
-    // Heat release rate time scale
-    if (alphaTemp < 1)
+    // Reaction source time scale
+    if (alphaTemp < 1.0)
     {
-        volScalarField::Internal rDeltaTT
+        volScalarField::DimensionedInternalField 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.dimensionedInternalField() = max
         (
-            IOobject
-            (
-                "rDeltaTY",
-                runTime.timeName(),
-                mesh
-            ),
-            mesh,
-            dimensionedScalar("rDeltaTY", rDeltaT.dimensions(), 0)
+            rDeltaT.dimensionedInternalField(),
+            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.internalField())
+        << ", " << gMax(1/rDeltaT.internalField()) << endl;
 }
 
 

applications/solvers/compressible/rhoCentralFoam/Allwclean

@@ -1,5 +1,6 @@
 #!/bin/sh
 cd ${0%/*} || exit 1    # Run from this directory
+set -x
 
 wclean libso BCs
 wclean

applications/solvers/compressible/rhoCentralFoam/Allwmake

@@ -2,8 +2,10 @@
 cd ${0%/*} || exit 1    # Run from this directory
 
 # Parse arguments for library compilation
+targetType=libso
 . $WM_PROJECT_DIR/wmake/scripts/AllwmakeParseArguments
+set -x
 
-(wmake $targetType BCs && wmake $targetType && wmake $targetType rhoCentralDyMFoam)
+(wmake $targetType BCs && wmake && wmake rhoCentralDyMFoam)
 
 #------------------------------------------------------------------------------

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

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

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

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

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

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

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

@@ -0,0 +1,196 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "mixedFixedValueSlipFvPatchField.H"
+#include "symmTransformField.H"
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+template<class Type>
+Foam::mixedFixedValueSlipFvPatchField<Type>::mixedFixedValueSlipFvPatchField
+(
+    const fvPatch& p,
+    const DimensionedField<Type, volMesh>& iF
+)
+:
+    transformFvPatchField<Type>(p, iF),
+    refValue_(p.size()),
+    valueFraction_(p.size(), 1.0)
+{}
+
+
+template<class Type>
+Foam::mixedFixedValueSlipFvPatchField<Type>::mixedFixedValueSlipFvPatchField
+(
+    const fvPatch& p,
+    const DimensionedField<Type, volMesh>& iF,
+    const dictionary& dict
+)
+:
+    transformFvPatchField<Type>(p, iF),
+    refValue_("refValue", dict, p.size()),
+    valueFraction_("valueFraction", dict, p.size())
+{}
+
+
+template<class Type>
+Foam::mixedFixedValueSlipFvPatchField<Type>::mixedFixedValueSlipFvPatchField
+(
+    const mixedFixedValueSlipFvPatchField<Type>& ptf,
+    const fvPatch& p,
+    const DimensionedField<Type, volMesh>& iF,
+    const fvPatchFieldMapper& mapper
+)
+:
+    transformFvPatchField<Type>(ptf, p, iF, mapper),
+    refValue_(ptf.refValue_, mapper),
+    valueFraction_(ptf.valueFraction_, mapper)
+{}
+
+
+template<class Type>
+Foam::mixedFixedValueSlipFvPatchField<Type>::mixedFixedValueSlipFvPatchField
+(
+    const mixedFixedValueSlipFvPatchField<Type>& ptf
+)
+:
+    transformFvPatchField<Type>(ptf),
+    refValue_(ptf.refValue_),
+    valueFraction_(ptf.valueFraction_)
+{}
+
+
+template<class Type>
+Foam::mixedFixedValueSlipFvPatchField<Type>::mixedFixedValueSlipFvPatchField
+(
+    const mixedFixedValueSlipFvPatchField<Type>& ptf,
+    const DimensionedField<Type, volMesh>& iF
+)
+:
+    transformFvPatchField<Type>(ptf, iF),
+    refValue_(ptf.refValue_),
+    valueFraction_(ptf.valueFraction_)
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+template<class Type>
+void Foam::mixedFixedValueSlipFvPatchField<Type>::autoMap
+(
+    const fvPatchFieldMapper& m
+)
+{
+    Field<Type>::autoMap(m);
+    refValue_.autoMap(m);
+    valueFraction_.autoMap(m);
+}
+
+
+template<class Type>
+void Foam::mixedFixedValueSlipFvPatchField<Type>::rmap
+(
+    const fvPatchField<Type>& ptf,
+    const labelList& addr
+)
+{
+    transformFvPatchField<Type>::rmap(ptf, addr);
+
+    const mixedFixedValueSlipFvPatchField<Type>& dmptf =
+        refCast<const mixedFixedValueSlipFvPatchField<Type>>(ptf);
+
+    refValue_.rmap(dmptf.refValue_, addr);
+    valueFraction_.rmap(dmptf.valueFraction_, addr);
+}
+
+
+template<class Type>
+Foam::tmp<Foam::Field<Type>>
+Foam::mixedFixedValueSlipFvPatchField<Type>::snGrad() const
+{
+    tmp<vectorField> nHat = this->patch().nf();
+    Field<Type> pif(this->patchInternalField());
+
+    return
+    (
+        valueFraction_*refValue_
+      + (1.0 - valueFraction_)*transform(I - sqr(nHat), pif) - pif
+    )*this->patch().deltaCoeffs();
+}
+
+
+template<class Type>
+void Foam::mixedFixedValueSlipFvPatchField<Type>::evaluate
+(
+    const Pstream::commsTypes
+)
+{
+    if (!this->updated())
+    {
+        this->updateCoeffs();
+    }
+
+    vectorField nHat(this->patch().nf());
+
+    Field<Type>::operator=
+    (
+        valueFraction_*refValue_
+      +
+        (1.0 - valueFraction_)
+       *transform(I - nHat*nHat, this->patchInternalField())
+    );
+
+    transformFvPatchField<Type>::evaluate();
+}
+
+
+template<class Type>
+Foam::tmp<Foam::Field<Type>>
+Foam::mixedFixedValueSlipFvPatchField<Type>::snGradTransformDiag() const
+{
+    vectorField nHat(this->patch().nf());
+    vectorField diag(nHat.size());
+
+    diag.replace(vector::X, mag(nHat.component(vector::X)));
+    diag.replace(vector::Y, mag(nHat.component(vector::Y)));
+    diag.replace(vector::Z, mag(nHat.component(vector::Z)));
+
+    return
+        valueFraction_*Type(pTraits<Type>::one)
+      + (1.0 - valueFraction_)
+       *transformFieldMask<Type>(pow<vector, pTraits<Type>::rank>(diag));
+}
+
+
+template<class Type>
+void Foam::mixedFixedValueSlipFvPatchField<Type>::write(Ostream& os) const
+{
+    transformFvPatchField<Type>::write(os);
+    refValue_.writeEntry("refValue", os);
+    valueFraction_.writeEntry("valueFraction", os);
+}
+
+
+// ************************************************************************* //

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

@@ -0,0 +1,236 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::mixedFixedValueSlipFvPatchField
+
+Description
+    A mixed boundary type that blends between fixedValue and slip, as opposed
+    to the standard mixed condition that blends between fixedValue and
+    fixedGradient; required to implement maxwellSlipU condition.
+
+SourceFiles
+    mixedFixedValueSlipFvPatchField.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef mixedFixedValueSlipFvPatchField_H
+#define mixedFixedValueSlipFvPatchField_H
+
+#include "transformFvPatchField.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                      Class mixedFixedValueSlipFvPatch Declaration
+\*---------------------------------------------------------------------------*/
+
+template<class Type>
+class mixedFixedValueSlipFvPatchField
+:
+    public transformFvPatchField<Type>
+{
+    // Private data
+
+        //- Value field used for boundary condition
+        Field<Type> refValue_;
+
+        //- Fraction (0-1) of value used for boundary condition
+        scalarField valueFraction_;
+
+public:
+
+    //- Runtime type information
+    TypeName("mixedFixedValueSlip");
+
+
+    // Constructors
+
+        //- Construct from patch and internal field
+        mixedFixedValueSlipFvPatchField
+        (
+            const fvPatch&,
+            const DimensionedField<Type, volMesh>&
+        );
+
+        //- Construct from patch, internal field and dictionary
+        mixedFixedValueSlipFvPatchField
+        (
+            const fvPatch&,
+            const DimensionedField<Type, volMesh>&,
+            const dictionary&
+        );
+
+        //- Construct by mapping given mixedFixedValueSlipFvPatchField
+        //- Onto a new patch
+        mixedFixedValueSlipFvPatchField
+        (
+            const mixedFixedValueSlipFvPatchField<Type>&,
+            const fvPatch&,
+            const DimensionedField<Type, volMesh>&,
+            const fvPatchFieldMapper&
+        );
+
+        //- Construct as copy
+        mixedFixedValueSlipFvPatchField
+        (
+            const mixedFixedValueSlipFvPatchField<Type>&
+        );
+
+        //- Construct and return a clone
+        virtual tmp<fvPatchField<Type>> clone() const
+        {
+            return tmp<fvPatchField<Type>>
+            (
+                new mixedFixedValueSlipFvPatchField<Type>(*this)
+            );
+        }
+
+        //- Construct as copy setting internal field reference
+        mixedFixedValueSlipFvPatchField
+        (
+            const mixedFixedValueSlipFvPatchField<Type>&,
+            const DimensionedField<Type, volMesh>&
+        );
+
+        //- Construct and return a clone setting internal field reference
+        virtual tmp<fvPatchField<Type>> clone
+        (
+            const DimensionedField<Type, volMesh>& iF
+        ) const
+        {
+            return tmp<fvPatchField<Type>>
+            (
+                new mixedFixedValueSlipFvPatchField<Type>(*this, iF)
+            );
+        }
+
+    // Member functions
+
+        // Access
+
+            //- Return false: this patch field is not altered by assignment
+            virtual bool assignable() const
+            {
+                return false;
+            }
+
+
+        // Mapping functions
+
+            //- Map (and resize as needed) from self given a mapping object
+            virtual void autoMap
+            (
+                const fvPatchFieldMapper&
+            );
+
+            //- Reverse map the given fvPatchField onto this fvPatchField
+            virtual void rmap
+            (
+                const fvPatchField<Type>&,
+                const labelList&
+            );
+
+        // Return defining fields
+
+            virtual Field<Type>& refValue()
+            {
+                return refValue_;
+            }
+
+            virtual const Field<Type>& refValue() const
+            {
+                return refValue_;
+            }
+
+            virtual scalarField& valueFraction()
+            {
+                return valueFraction_;
+            }
+
+            virtual const scalarField& valueFraction() const
+            {
+                return valueFraction_;
+            }
+
+        // Evaluation functions
+
+            //- Return gradient at boundary
+            virtual tmp<Field<Type>> snGrad() const;
+
+            //- Evaluate the patch field
+            virtual void evaluate
+            (
+                const Pstream::commsTypes commsType=Pstream::blocking
+            );
+
+            //- Return face-gradient transform diagonal
+            virtual tmp<Field<Type>> snGradTransformDiag() const;
+
+
+        //- Write
+        virtual void write(Ostream&) const;
+
+
+    // Member operators
+
+        virtual void operator=(const UList<Type>&) {}
+
+        virtual void operator=(const fvPatchField<Type>&) {}
+        virtual void operator+=(const fvPatchField<Type>&) {}
+        virtual void operator-=(const fvPatchField<Type>&) {}
+        virtual void operator*=(const fvPatchField<scalar>&) {}
+        virtual void operator/=(const fvPatchField<scalar>&) {}
+
+        virtual void operator+=(const Field<Type>&) {}
+        virtual void operator-=(const Field<Type>&) {}
+
+        virtual void operator*=(const Field<scalar>&) {}
+        virtual void operator/=(const Field<scalar>&) {}
+
+        virtual void operator=(const Type&) {}
+        virtual void operator+=(const Type&) {}
+        virtual void operator-=(const Type&) {}
+        virtual void operator*=(const scalar) {}
+        virtual void operator/=(const scalar) {}
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+    #include "mixedFixedValueSlipFvPatchField.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

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

@@ -0,0 +1,43 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011 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 "mixedFixedValueSlipFvPatchFields.H"
+#include "addToRunTimeSelectionTable.H"
+#include "volFields.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+makePatchFields(mixedFixedValueSlip);
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// ************************************************************************* //

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

@@ -0,0 +1,49 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011 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/>.
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef mixedFixedValueSlipFvPatchFields_H
+#define mixedFixedValueSlipFvPatchFields_H
+
+#include "mixedFixedValueSlipFvPatchField.H"
+#include "fieldTypes.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+makePatchTypeFieldTypedefs(mixedFixedValueSlip);
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

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

@@ -0,0 +1,50 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011 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/>.
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef mixedFixedValueSlipFvPatchFieldsFwd_H
+#define mixedFixedValueSlipFvPatchFieldsFwd_H
+
+#include "fieldTypes.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+template<class Type> class mixedFixedValueSlipFvPatchField;
+
+makePatchTypeFieldTypedefs(mixedFixedValueSlip);
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

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

@@ -44,16 +44,17 @@ Description
     \endvartable
 
 
-Usage
+    \heading Patch usage
+
     \table
-        Property     | Description                | Required    | Default value
-        p            | Pressure field name        | no          | p
-        psi          | Compressibility field name | no          | thermo:psi
+        Property     | Description             | Required    | Default value
+        pName        | Pressure field name     | no          | p
+        psiName      | Compressibility field name | no       | thermo:psi
     \endtable
 
     Example of the boundary condition specification:
     \verbatim
-    <patchName>
+    myPatch
     {
         type            fixedRho;
     }

applications/solvers/compressible/rhoCentralFoam/centralCourantNo.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2015 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -31,7 +31,7 @@ Description
 
 if (mesh.nInternalFaces())
 {
-    scalarField sumAmaxSf(fvc::surfaceSum(amaxSf)().primitiveField());
+    scalarField sumAmaxSf(fvc::surfaceSum(amaxSf)().internalField());
 
     CoNum = 0.5*gMax(sumAmaxSf/mesh.V().field())*runTime.deltaTValue();
 

applications/solvers/compressible/rhoCentralFoam/createFieldRefs.H

@@ -1,10 +0,0 @@
-volScalarField& p = thermo.p();
-const volScalarField& T = thermo.T();
-const volScalarField& psi = thermo.psi();
-const volScalarField& mu = thermo.mu();
-
-bool inviscid(true);
-if (max(mu.primitiveField()) > 0.0)
-{
-    inviscid = false;
-}

applications/solvers/compressible/rhoCentralFoam/createFields.H

@@ -1,5 +1,3 @@
-#include "createRDeltaT.H"
-
 Info<< "Reading thermophysical properties\n" << endl;
 
 autoPtr<psiThermo> pThermo
@@ -8,7 +6,17 @@ autoPtr<psiThermo> pThermo
 );
 psiThermo& thermo = pThermo();
 
+volScalarField& p = thermo.p();
 volScalarField& e = thermo.he();
+const volScalarField& T = thermo.T();
+const volScalarField& psi = thermo.psi();
+const volScalarField& mu = thermo.mu();
+
+bool inviscid(true);
+if (max(mu.internalField()) > 0.0)
+{
+    inviscid = false;
+}
 
 Info<< "Reading field U\n" << endl;
 volVectorField U

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

@@ -29,7 +29,7 @@ Group
 
 Description
     Density-based compressible flow solver based on central-upwind schemes of
-    Kurganov and Tadmor with support for mesh-motion and topology changes.
+    Kurganov and Tadmor with support for mesh-motion and topology changes
 
 \*---------------------------------------------------------------------------*/
 
@@ -39,22 +39,16 @@ Description
 #include "turbulentFluidThermoModel.H"
 #include "fixedRhoFvPatchScalarField.H"
 #include "directionInterpolate.H"
-#include "localEulerDdtScheme.H"
-#include "fvcSmooth.H"
 #include "motionSolver.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 int main(int argc, char *argv[])
 {
-    #define NO_CONTROL
-    #include "postProcess.H"
-
     #include "setRootCase.H"
     #include "createTime.H"
     #include "createDynamicFvMesh.H"
     #include "createFields.H"
-    #include "createFieldRefs.H"
     #include "createTimeControls.H"
 
     turbulence->validate();
@@ -120,14 +114,11 @@ int main(int argc, char *argv[])
             "cSf_pos",
             interpolate(c, pos, T.name())*mesh.magSf()
         );
-        cSf_pos.setOriented();
-
         surfaceScalarField cSf_neg
         (
             "cSf_neg",
             interpolate(c, neg, T.name())*mesh.magSf()
         );
-        cSf_neg.setOriented();
 
         surfaceScalarField ap
         (
@@ -197,11 +188,11 @@ int main(int argc, char *argv[])
         // --- Solve momentum
         solve(fvm::ddt(rhoU) + fvc::div(phiUp));
 
-        U.ref() =
-            rhoU()
-           /rho();
+        U.dimensionedInternalField() =
+            rhoU.dimensionedInternalField()
+           /rho.dimensionedInternalField();
         U.correctBoundaryConditions();
-        rhoU.boundaryFieldRef() == rho.boundaryField()*U.boundaryField();
+        rhoU.boundaryField() == rho.boundaryField()*U.boundaryField();
 
         if (!inviscid)
         {
@@ -235,7 +226,7 @@ int main(int argc, char *argv[])
         e = rhoE/rho - 0.5*magSqr(U);
         e.correctBoundaryConditions();
         thermo.correct();
-        rhoE.boundaryFieldRef() ==
+        rhoE.boundaryField() ==
             rho.boundaryField()*
             (
                 e.boundaryField() + 0.5*magSqr(U.boundaryField())
@@ -252,11 +243,11 @@ int main(int argc, char *argv[])
             rhoE = rho*(e + 0.5*magSqr(U));
         }
 
-        p.ref() =
-            rho()
-           /psi();
+        p.dimensionedInternalField() =
+            rho.dimensionedInternalField()
+           /psi.dimensionedInternalField();
         p.correctBoundaryConditions();
-        rho.boundaryFieldRef() == psi.boundaryField()*p.boundaryField();
+        rho.boundaryField() == psi.boundaryField()*p.boundaryField();
 
         turbulence->correct();
 
@@ -272,5 +263,4 @@ int main(int argc, char *argv[])
     return 0;
 }
 
-
 // ************************************************************************* //

applications/solvers/compressible/rhoCentralFoam/rhoCentralFoam.C

@@ -29,7 +29,7 @@ Group
 
 Description
     Density-based compressible flow solver based on central-upwind schemes of
-    Kurganov and Tadmor.
+    Kurganov and Tadmor
 
 \*---------------------------------------------------------------------------*/
 
@@ -45,15 +45,13 @@ Description
 
 int main(int argc, char *argv[])
 {
-    #define NO_CONTROL
-    #include "postProcess.H"
-
     #include "setRootCase.H"
+
     #include "createTime.H"
     #include "createMesh.H"
     #include "createFields.H"
-    #include "createFieldRefs.H"
     #include "createTimeControls.H"
+    #include "createRDeltaT.H"
 
     turbulence->validate();
 
@@ -101,14 +99,11 @@ int main(int argc, char *argv[])
             "cSf_pos",
             interpolate(c, pos, T.name())*mesh.magSf()
         );
-        cSf_pos.setOriented();
-
         surfaceScalarField cSf_neg
         (
             "cSf_neg",
             interpolate(c, neg, T.name())*mesh.magSf()
         );
-        cSf_neg.setOriented();
 
         surfaceScalarField ap
         (
@@ -186,11 +181,11 @@ int main(int argc, char *argv[])
         // --- Solve momentum
         solve(fvm::ddt(rhoU) + fvc::div(phiUp));
 
-        U.ref() =
-            rhoU()
-           /rho();
+        U.dimensionedInternalField() =
+            rhoU.dimensionedInternalField()
+           /rho.dimensionedInternalField();
         U.correctBoundaryConditions();
-        rhoU.boundaryFieldRef() == rho.boundaryField()*U.boundaryField();
+        rhoU.boundaryField() == rho.boundaryField()*U.boundaryField();
 
         if (!inviscid)
         {
@@ -224,7 +219,7 @@ int main(int argc, char *argv[])
         e = rhoE/rho - 0.5*magSqr(U);
         e.correctBoundaryConditions();
         thermo.correct();
-        rhoE.boundaryFieldRef() ==
+        rhoE.boundaryField() ==
             rho.boundaryField()*
             (
                 e.boundaryField() + 0.5*magSqr(U.boundaryField())
@@ -241,11 +236,11 @@ int main(int argc, char *argv[])
             rhoE = rho*(e + 0.5*magSqr(U));
         }
 
-        p.ref() =
-            rho()
-           /psi();
+        p.dimensionedInternalField() =
+            rho.dimensionedInternalField()
+           /psi.dimensionedInternalField();
         p.correctBoundaryConditions();
-        rho.boundaryFieldRef() == psi.boundaryField()*p.boundaryField();
+        rho.boundaryField() == psi.boundaryField()*p.boundaryField();
 
         turbulence->correct();