Merge branch 'develop' of develop.openfoam.com:Development/OpenFOAM-plus into develop

applications/solvers/combustion/PDRFoam/createFields.H

@@ -58,23 +58,6 @@ 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
 (
@@ -227,4 +210,8 @@ fields.add(thermo.he());
 fields.add(thermo.heu());
 flameWrinkling->addXi(fields);
 
+#include "createDpdt.H"
+
+#include "createK.H"
+
 #include "createMRF.H"

applications/solvers/combustion/XiFoam/createFields.H

@@ -58,21 +58,6 @@ 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
@@ -139,4 +124,8 @@ fields.add(b);
 fields.add(thermo.he());
 fields.add(thermo.heu());
 
+#include "createDpdt.H"
+
+#include "createK.H"
+
 #include "createMRF.H"

applications/solvers/combustion/coldEngineFoam/createFields.H

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

applications/solvers/combustion/fireFoam/createFields.H

@@ -145,22 +145,9 @@ volScalarField Qdot
     dimensionedScalar("Qdot", dimEnergy/dimVolume/dimTime, 0.0)
 );
 
+#include "createDpdt.H"
 
-Info<< "Creating field dpdt\n" << endl;
+#include "createK.H"
-volScalarField dpdt
-(
-    IOobject
-    (
-        "dpdt",
-        runTime.timeName(),
-        mesh
-    ),
-    mesh,
-    dimensionedScalar("dpdt", p.dimensions()/dimTime, 0)
-);
-
-Info<< "Creating field kinetic energy K\n" << endl;
-volScalarField K("K", 0.5*magSqr(U));
 
 #include "createClouds.H"
 #include "createSurfaceFilmModel.H"

applications/solvers/combustion/reactingFoam/createFields.H

@@ -71,24 +71,6 @@ 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)
@@ -111,4 +93,8 @@ volScalarField Qdot
     dimensionedScalar("Qdot", dimEnergy/dimVolume/dimTime, 0.0)
 );
 
+#include "createDpdt.H"
+
+#include "createK.H"
+
 #include "createMRF.H"

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

@@ -92,23 +92,6 @@ 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)
@@ -131,4 +114,8 @@ volScalarField Qdot
     dimensionedScalar("Qdot", dimEnergy/dimVolume/dimTime, 0.0)
 );
 
+#include "createDpdt.H"
+
+#include "createK.H"
+
 #include "createMRF.H"

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

@@ -72,24 +72,6 @@ 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)
@@ -112,6 +94,10 @@ volScalarField Qdot
     dimensionedScalar("Qdot", dimEnergy/dimVolume/dimTime, 0.0)
 );
 
+#include "createDpdt.H"
+
+#include "createK.H"
+
 #include "createMRF.H"
 
 dimensionedScalar rhoMax

applications/solvers/compressible/rhoPimpleFoam/createFields.H

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

applications/solvers/compressible/sonicFoam/createFields.H

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

applications/solvers/heatTransfer/buoyantPimpleFoam/createFields.H

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

applications/solvers/heatTransfer/thermoFoam/createFields.H

@@ -37,20 +37,8 @@ volVectorField U
 
 #include "setAlphaEff.H"
 
-Info<< "Creating field dpdt\n" << endl;
+#include "createDpdt.H"
-volScalarField dpdt
-(
-    IOobject
-    (
-        "dpdt",
-        runTime.timeName(),
-        mesh
-    ),
-    mesh,
-    dimensionedScalar("dpdt", p.dimensions()/dimTime, 0)
-);
 
-Info<< "Creating field kinetic energy K\n" << endl;
+#include "createK.H"
-volScalarField K("K", 0.5*magSqr(U));
 
 #include "createRadiationModel.H"

applications/solvers/lagrangian/coalChemistryFoam/createFields.H

@@ -115,22 +115,6 @@ autoPtr<compressible::turbulenceModel> turbulence
 // Set the turbulence into the combustion model
 combustion->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));
-
 volScalarField Qdot
 (
     IOobject
@@ -145,6 +129,10 @@ volScalarField Qdot
     dimensionedScalar("Qdot", dimEnergy/dimVolume/dimTime, 0.0)
 );
 
+#include "createDpdt.H"
+
+#include "createK.H"
+
 #include "createMRF.H"
 #include "createClouds.H"
 #include "createRadiationModel.H"

applications/solvers/lagrangian/reactingParcelFilmFoam/createFields.H

@@ -69,23 +69,6 @@ autoPtr<compressible::turbulenceModel> turbulence
 // Set the turbulence into the combustion model
 combustion->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));
-
-
 #include "readGravitationalAcceleration.H"
 #include "readhRef.H"
 #include "gh.H"
@@ -148,6 +131,10 @@ volScalarField Qdot
     dimensionedScalar("Qdot", dimEnergy/dimVolume/dimTime, 0.0)
 );
 
+#include "createDpdt.H"
+
+#include "createK.H"
+
 #include "createMRF.H"
 #include "createClouds.H"
 #include "createRadiationModel.H"

applications/solvers/lagrangian/reactingParcelFoam/createFields.H

@@ -96,22 +96,6 @@ autoPtr<compressible::turbulenceModel> turbulence
 // Set the turbulence into the combustion model
 combustion->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));
-
 Info<< "Creating multi-variate interpolation scheme\n" << endl;
 multivariateSurfaceInterpolationScheme<scalar>::fieldTable fields;
 
@@ -135,6 +119,10 @@ volScalarField Qdot
     dimensionedScalar("Qdot", dimEnergy/dimVolume/dimTime, 0.0)
 );
 
+#include "createDpdt.H"
+
+#include "createK.H"
+
 #include "createMRF.H"
 #include "createRadiationModel.H"
 #include "createClouds.H"

applications/solvers/lagrangian/sprayFoam/createFields.H

@@ -94,21 +94,6 @@ autoPtr<compressible::turbulenceModel> turbulence
 // Set the turbulence into the combustion model
 combustion->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)
@@ -131,6 +116,10 @@ volScalarField Qdot
     dimensionedScalar("Qdot", dimEnergy/dimVolume/dimTime, 0.0)
 );
 
+#include "createDpdt.H"
+
+#include "createK.H"
+
 #include "createMRF.H"
 #include "createClouds.H"
 #include "createRadiationModel.H"

applications/solvers/multiphase/compressibleInterFoam/createFields.H

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

applications/solvers/multiphase/compressibleMultiphaseInterFoam/createFields.H

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

applications/solvers/multiphase/interCondensingEvaporatingFoam/createFields.H

@@ -124,9 +124,6 @@ volScalarField kappaEff
     thermo->kappa()
 );
 
-Info<< "Creating field kinetic energy K\n" << endl;
-volScalarField K("K", 0.5*magSqr(U));
-
 Info<< "Creating field pDivU\n" << endl;
 volScalarField pDivU
 (

applications/solvers/multiphase/interCondensingEvaporatingFoam/pEqn.H

@@ -47,7 +47,6 @@
             U = HbyA + rAU*fvc::reconstruct((phig + p_rghEqn.flux())/rAUf);
             U.correctBoundaryConditions();
             fvOptions.correct(U);
-            K = 0.5*magSqr(U);
         }
     }
 

applications/utilities/postProcessing/dataConversion/foamToVTK/foamVtkLagrangianWriter.C

@@ -100,7 +100,7 @@ void Foam::vtk::lagrangianWriter::writeVertsLegacy()
 
     for (label i=0; i < nParcels_; ++i)
     {
-        format().write(1);  // Number of vertices for this cell (==1)
+        format().write(label(1)); // Number of vertices for this cell (==1)
         format().write(i);
     }
     format().flush();

bin/tools/foamCreateBashCompletions

@@ -29,40 +29,56 @@
 #     Create bash completions for OpenFOAM applications
 #
 #------------------------------------------------------------------------------
-#set -x
-
 usage() {
     exec 1>&2
     while [ "$#" -ge 1 ]; do echo "$1"; shift; done
     cat<<USAGE
 
-Usage: $Script [OPTION] <file>
+Usage: ${0##*/} [OPTION] <outputFile>
+options:
+  -d dir | -dir dir       Directory to process
+  -u | -user              Add \$FOAM_USER_APPBIN to the search directories
+  -h | -help              Print the usage
 
-* Create bash completions for OpenFOAM applications and write to <file>.
+Create bash completions for OpenFOAM applications and write to <outputFile>.
-  By default searches directories \$FOAM_APPBIN and \$FOAM_USER_APPBIN
+By default searches \$FOAM_APPBIN only.
-
-
-Options:
-  -d | -directory   Directory to process
-  -h | -help        Print the usage
 
 USAGE
     exit 1
 }
 
+# Report error and exit
+die()
+{
+    exec 1>&2
+    echo
+    echo "Error encountered:"
+    while [ "$#" -ge 1 ]; do echo "    $1"; shift; done
+    echo
+    echo "See '${0##*/} -help' for usage"
+    echo
+    exit 1
+}
+
+#-------------------------------------------------------------------------------
+#set -x
+
 unset outFile
-searchDirs="$FOAM_APPBIN $FOAM_USER_APPBIN"
+searchDirs="$FOAM_APPBIN"
 while [ "$#" -gt 0 ]
 do
     case "$1" in
         -h | -help)
             usage
             ;;
-        -d | -directory)
+        -d | -dir)
             searchDirs="$2"
             [ -d $searchDirs ] || usage "directory not found '$searchDirs'"
             shift
             ;;
+        -u | -user)
+            searchDirs="$searchDirs $FOAM_USER_APPBIN"
+            ;;
         -*)
             usage "unknown option: '$1'"
             ;;
@@ -74,14 +90,14 @@ do
     shift
 done
 
-[ -n "$outFile" ] || usage
+[ -n "$outFile" ] || usage "No output file specified"
-\rm -f $outFile
-touch $outFile
 
 
-writeFilterFunction()
+# Generate header
-{
+cat << HEADER > $outFile
-    cat<< WRITEFILTER >> $1
+#----------------------------------*-sh-*--------------------------------------
+# Bash completions for OpenFOAM applications
+
 unset -f _of_filter_opts 2>/dev/null
 _of_filter_opts()
 {
@@ -92,72 +108,74 @@ _of_filter_opts()
     done
 }
 
-WRITEFILTER
+HEADER
-}
 
+#------------------------------------------------------------------------------
 
+#
+# Produce contents for switch for common options
+#
 commonOptions()
 {
-    local options=$@
+    local indent1="    "
-    local indent1="        "
+    local indent2="        "
-    local indent2="            "
+    for opt
-    for o in ${options[@]}; do
+    do
-        case $o in
+        case $opt in
-            -case)
+        -case)
-                echo "${indent1}-case)"
+            echo "${indent1}-case)"
-                echo "${indent2}COMPREPLY=(\$(compgen -d -- \${cur}))"
+            echo "${indent2}COMPREPLY=(\$(compgen -d -- \${cur}))"
-                echo "${indent2};;"
+            echo "${indent2};;"
-                ;;
+            ;;
-            -srcDoc|-help)
+        -srcDoc|-help)
-                echo "${indent1}-srcDoc|-help)"
+            echo "${indent1}-srcDoc|-help)"
-                echo "${indent2}COMPREPLY=()"
+            echo "${indent2}COMPREPLY=()"
-                echo "${indent2};;"
+            echo "${indent2};;"
-                ;;
+            ;;
-            -time)
+        -time)
-                echo "${indent1}-time)"
+            echo "${indent1}-time)"
-                echo "${indent2}COMPREPLY=(\$(compgen -d -X '![-0-9]*' -- \${cur}))"
+            echo "${indent2}COMPREPLY=(\$(compgen -d -X '![-0-9]*' -- \${cur}))"
-                echo "${indent2};;"
+            echo "${indent2};;"
-                ;;
+            ;;
-            -region)
+        -region)
-                echo "${indent1}-region)"
+            echo "${indent1}-region)"
-                echo "${indent2}local regions=\$(sed 's#/##g' <<< \$([ -d system ] && (\cd system && (\ls -d */ 2>/dev/null))))"
+            echo "${indent2}local regions=\$(sed 's#/##g' <<< \$([ -d system ] && (\cd system && (\ls -d */ 2>/dev/null))))"
-                echo "${indent2}COMPREPLY=(\$(compgen -W \"\$regions\" -- \${cur}))"
+            echo "${indent2}COMPREPLY=(\$(compgen -W \"\$regions\" -- \${cur}))"
-                echo "${indent2};;"
+            echo "${indent2};;"
-                ;;
+            ;;
-            *Dict)
+        *Dict)
-                echo "${indent1}*Dict)"
+            echo "${indent1}*Dict)"
-#                echo "${indent2}local dirs=\$(\ls -d s*/)"
+#           echo "${indent2}local dirs=\$(\ls -d s*/)"
-#                echo "${indent2}local files=\$(\ls -f | grep Dict)"
+#           echo "${indent2}local files=\$(\ls -f | grep Dict)"
-#                echo "${indent2}COMPREPLY=(\$(compgen -W \"\$dirs \$files\" -- \${cur}))"
+#           echo "${indent2}COMPREPLY=(\$(compgen -W \"\$dirs \$files\" -- \${cur}))"
-                echo "${indent2}COMPREPLY=(\$(compgen -f -- \${cur}))"
+            echo "${indent2}COMPREPLY=(\$(compgen -f -- \${cur}))"
-                echo "${indent2};;"
+            echo "${indent2};;"
-                ;;
+            ;;
         esac
     done
 }
 
 
-# Generate header
+#------------------------------------------------------------------------------
-cat << HEADER > $outFile
-#----------------------------------*-sh-*--------------------------------------
-# Bash completions for OpenFOAM applications
-
-HEADER
-
-writeFilterFunction $outFile
 
 for dir in ${searchDirs}
 do
-    echo "Processing directory $dir" 1>&2
+    if [ -d "$dir" ]
+    then
+        echo "Processing directory $dir" 1>&2
+    else
+        echo "No such directory: $dir" 1>&2
+        continue
+    fi
 
     # Sort with ignore-case
-    apps=($(\ls $dir | sort -f))
+    set -- $(\ls $dir | sort -f)
-    for appName in "${apps[@]}"
+    for appName
     do
         [ -f "$dir/$appName" -a -x "$dir/$appName" ] || continue
         appHelp=$($appName -help)
 
-        echo "Processing $appName" 1>&2
+        echo "    $appName" 1>&2
 
         # Options with args
         optsWithArgs=($(awk '/^ {0,4}-[a-z]/ && /</ {print $1}' <<< "$appHelp"))
@@ -165,7 +183,7 @@ do
         # Options without args
         opts=($(awk '/^ {0,4}-[a-z]/ && !/</ {print $1}' <<< "$appHelp"))
 
-        cat<<WRITECOMPLETION >> $outFile
+        cat << WRITECOMPLETION >> $outFile
 unset -f _of_${appName} 2>/dev/null
 _of_${appName}()
 {
@@ -177,19 +195,19 @@ _of_${appName}()
 
     case \${prev} in
 $(commonOptions ${optsWithArgs[@]})
-        *)
+    *)
-            if [ "\${optsWithArgs/\${prev} /}" != "\${optsWithArgs}" ]; then
+        if [ "\${optsWithArgs/\${prev} /}" != "\${optsWithArgs}" ]
-                # Unknown what type of arg follows - set to files
+        then
-                # not always correct but at least can still navigate path if
+            # Unknown type of arg follows - set to files.
-                # needed...
+            # Not always correct but can still navigate path if needed...
-                COMPREPLY=(\$(compgen -f -- \${cur}))
+            COMPREPLY=(\$(compgen -f -- \${cur}))
-            else
+        else
-                # Catch-all - present all remaining options
+            # Catch-all - present all remaining options
-                opts=\$(_of_filter_opts "\${opts}" "\${COMP_LINE}")
+            opts=\$(_of_filter_opts "\${opts}" "\${COMP_LINE}")
-                optsWithArgs=\$(_of_filter_opts "\${optsWithArgs}" "\${COMP_LINE}")
+            optsWithArgs=\$(_of_filter_opts "\${optsWithArgs}" "\${COMP_LINE}")
-                COMPREPLY=(\$(compgen -W "\${opts} \${optsWithArgs}" -- \${cur}))
+            COMPREPLY=(\$(compgen -W "\${opts} \${optsWithArgs}" -- \${cur}))
-            fi
+        fi
-            ;;
+        ;;
     esac
 
     return 0

etc/bashrc

@@ -169,12 +169,16 @@ export PATH LD_LIBRARY_PATH MANPATH
 # Source project setup files
 # ~~~~~~~~~~~~~~~~~~~~~~~~~~
 _foamEtc config.sh/settings
-_foamEtc config.sh/aliases
 
-# Bash completions
+if /usr/bin/tty -s 2>/dev/null  # Interactive shell
-if command -v complete > /dev/null 2>&1
 then
-    _foamEtc config.sh/bashcompletion
+    _foamEtc config.sh/aliases
+
+    # Bash completions
+    if command -v complete > /dev/null 2>&1
+    then
+        _foamEtc config.sh/bashcompletion
+    fi
 fi
 
 

etc/config.csh/mpi

@@ -49,7 +49,7 @@ case SYSTEMOPENMPI:
     breaksw
 
 case OPENMPI:
-    setenv FOAM_MPI openmpi-2.0.2
+    setenv FOAM_MPI openmpi-2.1.1
     _foamEtc config.csh/openmpi                 # <- Adjustments (optional)
 
     setenv MPI_ARCH_PATH $WM_THIRD_PARTY_DIR/platforms/$WM_ARCH$WM_COMPILER/$FOAM_MPI

etc/config.csh/vtk

@@ -36,15 +36,15 @@
 #     equivalent config.sh version too, since that is the one which will
 #     be used during the build process.
 #
-# It is recommended to use the VTK sources from ParaView 5.0.1.
+# It is recommended to use VTK sources from ParaView (5.0.1 or later)
 #
-# See the BuildIssues.txt about problems that can be encountered when using
+# See BuildIssues.txt about problems that can be encountered when using
 # the 'plain' VTK sources.
 #------------------------------------------------------------------------------
 # USER EDITABLE PART: Changes made here may be lost with the next upgrade
 
-set vtk_version=VTK-7.1.0
+set vtk_version=VTK-8.1.0
-set mesa_version=mesa-13.0.1
+set mesa_version=mesa-17.1.1
 
 setenv VTK_DIR $WM_THIRD_PARTY_DIR/platforms/$WM_ARCH$WM_COMPILER/$vtk_version
 setenv MESA_ARCH_PATH $WM_THIRD_PARTY_DIR/platforms/$WM_ARCH$WM_COMPILER/$mesa_version

etc/config.sh/mpi

@@ -55,7 +55,7 @@ SYSTEMOPENMPI)
     ;;
 
 OPENMPI)
-    export FOAM_MPI=openmpi-2.0.2
+    export FOAM_MPI=openmpi-2.1.1
     _foamEtc config.sh/openmpi                  # <- Adjustments (optional)
 
     export MPI_ARCH_PATH=$WM_THIRD_PARTY_DIR/platforms/$WM_ARCH$WM_COMPILER/$FOAM_MPI

etc/config.sh/vtk

@@ -36,15 +36,15 @@
 #         - the vtk_version, mesa_version variables are retained.
 #           (for future integration into ThirdParty build)
 #
-# It is recommended to use the VTK sources from ParaView 5.0.1 or later
+# It is recommended to use VTK sources from ParaView (5.0.1 or later)
 #
-# See the BuildIssues.txt about problems that can be encountered when using
+# See BuildIssues.txt about problems that can be encountered when using
 # the 'plain' VTK sources.
 #------------------------------------------------------------------------------
 # USER EDITABLE PART: Changes made here may be lost with the next upgrade
 
-vtk_version=VTK-7.1.0
+vtk_version=VTK-8.1.0
-mesa_version=mesa-13.0.1
+mesa_version=mesa-17.1.1
 
 export VTK_DIR=$WM_THIRD_PARTY_DIR/platforms/$WM_ARCH$WM_COMPILER/$vtk_version
 export MESA_ARCH_PATH=$WM_THIRD_PARTY_DIR/platforms/$WM_ARCH$WM_COMPILER/$mesa_version

etc/cshrc

@@ -216,7 +216,11 @@ if ( $status == 0 ) setenv MANPATH $cleaned
 # Source project setup files
 # ~~~~~~~~~~~~~~~~~~~~~~~~~~
 _foamEtc config.csh/settings
-_foamEtc config.csh/aliases
+
+if ($?prompt) then  # Interactive shell
+    _foamEtc config.csh/aliases
+endif
+
 
 # Source user setup files for optional packages
 # ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

src/OpenFOAM/fields/GeometricFields/GeometricField/GeometricField.C

@@ -770,10 +770,11 @@ void Foam::GeometricField<Type, PatchField, GeoMesh>::storeOldTimes() const
     )
     {
         storeOldTime();
+        timeIndex_ = this->time().timeIndex();
     }
 
     // Correct time index
-    timeIndex_ = this->time().timeIndex();
+    //timeIndex_ = this->time().timeIndex();
 }
 
 
@@ -835,6 +836,17 @@ Foam::GeometricField<Type, PatchField, GeoMesh>::oldTime() const
             ),
             *this
         );
+
+        if (debug)
+        {
+            InfoInFunction
+                << "created old time field " << field0Ptr_->info() << endl;
+
+            if (debug&2)
+            {
+                error::printStack(Info);
+            }
+        }
     }
     else
     {

src/finiteVolume/cfdTools/compressible/createDpdt.H

@@ -0,0 +1,18 @@
+Info<< "Creating field dpdt\n" << endl;
+volScalarField dpdt
+(
+    IOobject
+    (
+        "dpdt",
+        runTime.timeName(),
+        mesh,
+        IOobject::NO_READ,
+        IOobject::NO_WRITE
+    ),
+    fvc::ddt(p)
+);
+
+if (!thermo.dpdt())
+{
+    dpdt == dimensionedScalar("0", dpdt.dimensions(), 0);
+}

src/finiteVolume/cfdTools/general/include/createK.H

@@ -0,0 +1,16 @@
+Info<< "Creating field kinetic energy K\n" << endl;
+volScalarField K("K", 0.5*magSqr(U));
+
+if (U.nOldTimes())
+{
+    volVectorField* Uold = &U.oldTime();
+    volScalarField* Kold = &K.oldTime();
+    *Kold == 0.5*magSqr(*Uold);
+
+    while (Uold->nOldTimes())
+    {
+        Uold = &Uold->oldTime();
+        Kold = &Kold->oldTime();
+        *Kold == 0.5*magSqr(*Uold);
+    }
+}

src/finiteVolume/finiteVolume/ddtSchemes/backwardDdtScheme/backwardDdtScheme.C

@@ -3,7 +3,7 @@
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
     \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
-     \\/     M anipulation  |
+     \\/     M anipulation  | Copyright (C) 2017 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -58,7 +58,7 @@ template<class Type>
 template<class GeoField>
 scalar backwardDdtScheme<Type>::deltaT0_(const GeoField& vf) const
 {
-    if (vf.nOldTimes() < 2)
+    if (mesh().time().timeIndex() < 2)
     {
         return GREAT;
     }

src/finiteVolume/finiteVolume/ddtSchemes/backwardDdtScheme/backwardDdtScheme.H

@@ -92,27 +92,13 @@ public:
         backwardDdtScheme(const fvMesh& mesh)
         :
             ddtScheme<Type>(mesh)
-        {
+        {}
-            // Ensure the old-old-time cell volumes are available
-            // for moving meshes
-            if (mesh.moving())
-            {
-                mesh.V00();
-            }
-        }
 
         //- Construct from mesh and Istream
         backwardDdtScheme(const fvMesh& mesh, Istream& is)
         :
             ddtScheme<Type>(mesh, is)
-        {
+        {}
-            // Ensure the old-old-time cell volumes are available
-            // for moving meshes
-            if (mesh.moving())
-            {
-                mesh.V00();
-            }
-        }
 
 
     // Member Functions

src/finiteVolume/fvMesh/fvMesh.C

@@ -198,8 +198,8 @@ void Foam::fvMesh::storeOldVol(const scalarField& V)
                     time().timeName(),
                     *this,
                     IOobject::NO_READ,
-                    IOobject::NO_WRITE,
+                    IOobject::AUTO_WRITE,
-                    false
+                    true
                 ),
                 *this,
                 dimVolume
@@ -281,8 +281,8 @@ Foam::fvMesh::fvMesh(const IOobject& io)
                 time().timeName(),
                 *this,
                 IOobject::MUST_READ,
-                IOobject::NO_WRITE,
+                IOobject::AUTO_WRITE,
-                false
+                true
             ),
             *this
         );

src/finiteVolume/fvMesh/fvMeshGeometry.C

@@ -260,11 +260,12 @@ const Foam::volScalarField::Internal& Foam::fvMesh::V00() const
                 *this,
                 IOobject::NO_READ,
                 IOobject::NO_WRITE,
-                false
+                true
             ),
             V0()
         );
 
+
         // If V00 is used then V0 should be stored for restart
         V0Ptr_->writeOpt() = IOobject::AUTO_WRITE;
     }

src/meshTools/coordinateSystems/coordinateRotation/EulerCoordinateRotation.C

@@ -3,7 +3,7 @@
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
     \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
-     \\/     M anipulation  |
+     \\/     M anipulation  | Copyright (C) 2017 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -195,7 +195,17 @@ void Foam::EulerCoordinateRotation::calcTransform
 Foam::EulerCoordinateRotation::EulerCoordinateRotation()
 :
     R_(sphericalTensor::I),
-    Rtr_(R_)
+    Rtr_(sphericalTensor::I)
+{}
+
+
+Foam::EulerCoordinateRotation::EulerCoordinateRotation
+(
+    const EulerCoordinateRotation& r
+)
+:
+    R_(r.R_),
+    Rtr_(r.Rtr_)
 {}
 
 
@@ -206,7 +216,7 @@ Foam::EulerCoordinateRotation::EulerCoordinateRotation
 )
 :
     R_(sphericalTensor::I),
-    Rtr_(R_)
+    Rtr_(sphericalTensor::I)
 {
     calcTransform
     (
@@ -227,7 +237,7 @@ Foam::EulerCoordinateRotation::EulerCoordinateRotation
 )
 :
     R_(sphericalTensor::I),
-    Rtr_(R_)
+    Rtr_(sphericalTensor::I)
 {
     calcTransform(phiAngle, thetaAngle, psiAngle, inDegrees);
 }
@@ -239,9 +249,9 @@ Foam::EulerCoordinateRotation::EulerCoordinateRotation
 )
 :
     R_(sphericalTensor::I),
-    Rtr_(R_)
+    Rtr_(sphericalTensor::I)
 {
-    vector rotation(dict.lookup("rotation"));
+    const vector rotation(dict.lookup("rotation"));
 
     calcTransform
     (
@@ -259,36 +269,8 @@ Foam::EulerCoordinateRotation::EulerCoordinateRotation
     const objectRegistry&
 )
 :
-    R_(sphericalTensor::I),
+    EulerCoordinateRotation(dict)
-    Rtr_(R_)
-{
-    vector rotation(dict.lookup("rotation"));
-
-    calcTransform
-    (
-        rotation.component(vector::X),
-        rotation.component(vector::Y),
-        rotation.component(vector::Z),
-        dict.lookupOrDefault("degrees", true)
-    );
-}
-
-
-Foam::EulerCoordinateRotation::EulerCoordinateRotation
-(
-    const EulerCoordinateRotation& r
-)
-:
-    R_(r.R_),
-    Rtr_(r.Rtr_)
 {}
 
 
-void Foam::EulerCoordinateRotation::write(Ostream& os) const
-{
-     os.writeKeyword("e1") << e1() << token::END_STATEMENT << nl;
-     os.writeKeyword("e2") << e2() << token::END_STATEMENT << nl;
-     os.writeKeyword("e3") << e3() << token::END_STATEMENT << nl;
-}
-
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

src/meshTools/coordinateSystems/coordinateRotation/EulerCoordinateRotation.H

@@ -3,7 +3,7 @@
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
     \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
-     \\/     M anipulation  |
+     \\/     M anipulation  | Copyright (C) 2017 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -91,7 +91,7 @@ class EulerCoordinateRotation
             const scalar phiAngle,
             const scalar thetaAngle,
             const scalar psiAngle,
-            const bool inDegrees=true
+            const bool inDegrees
         );
 
 
@@ -106,11 +106,14 @@ public:
         //- Construct null
         EulerCoordinateRotation();
 
+        //- Construct as copy
+        EulerCoordinateRotation(const EulerCoordinateRotation& r);
+
         //- Construct from rotation vector
         EulerCoordinateRotation
         (
             const vector& phiThetaPsi,
-            const bool inDegrees=true
+            const bool inDegrees
         );
 
         //- Construct from components of rotation vector
@@ -119,17 +122,18 @@ public:
             const scalar phiAngle,
             const scalar thetaAngle,
             const scalar psiAngle,
-            const bool inDegrees=true
+            const bool inDegrees
         );
 
         //- Construct from dictionary
-        EulerCoordinateRotation(const dictionary&);
+        explicit EulerCoordinateRotation(const dictionary& dict);
 
-        //- Construct from dictionary and mesh
+        //- Construct from dictionary and a registry (typically a mesh)
-        EulerCoordinateRotation(const dictionary&, const objectRegistry&);
+        EulerCoordinateRotation
-
+        (
-        //- Construct as copy
+            const dictionary& dict,
-        EulerCoordinateRotation(const EulerCoordinateRotation&);
+            const objectRegistry& unused
+        );
 
         //- Return clone
         autoPtr<coordinateRotation> clone() const
@@ -226,11 +230,6 @@ public:
         // symmetrical tensor
         virtual symmTensor transformVector(const vector& st) const;
 
-
-    // Write
-
-        //- Write
-        virtual void write(Ostream&) const;
 };
 
 

src/meshTools/coordinateSystems/coordinateRotation/STARCDCoordinateRotation.C

@@ -3,7 +3,7 @@
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
     \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
-     \\/     M anipulation  |
+     \\/     M anipulation  | Copyright (C) 2017 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -147,6 +147,7 @@ Foam::symmTensor Foam::STARCDCoordinateRotation::transformVector
     return transformPrincipal(R_, st);
 }
 
+
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
 void Foam::STARCDCoordinateRotation::calcTransform
@@ -195,7 +196,17 @@ void Foam::STARCDCoordinateRotation::calcTransform
 Foam::STARCDCoordinateRotation::STARCDCoordinateRotation()
 :
     R_(sphericalTensor::I),
-    Rtr_(R_)
+    Rtr_(sphericalTensor::I)
+{}
+
+
+Foam::STARCDCoordinateRotation::STARCDCoordinateRotation
+(
+    const STARCDCoordinateRotation& r
+)
+:
+    R_(r.R_),
+    Rtr_(r.Rtr_)
 {}
 
 
@@ -206,7 +217,7 @@ Foam::STARCDCoordinateRotation::STARCDCoordinateRotation
 )
 :
     R_(sphericalTensor::I),
-    Rtr_(R_)
+    Rtr_(sphericalTensor::I)
 {
     calcTransform
     (
@@ -227,7 +238,7 @@ Foam::STARCDCoordinateRotation::STARCDCoordinateRotation
 )
 :
     R_(sphericalTensor::I),
-    Rtr_(R_)
+    Rtr_(sphericalTensor::I)
 {
     calcTransform(rotZ, rotX, rotY, inDegrees);
 }
@@ -239,9 +250,9 @@ Foam::STARCDCoordinateRotation::STARCDCoordinateRotation
 )
 :
     R_(sphericalTensor::I),
-    Rtr_(R_)
+    Rtr_(sphericalTensor::I)
 {
-    vector rotation(dict.lookup("rotation"));
+    const vector rotation(dict.lookup("rotation"));
 
     calcTransform
     (
@@ -258,34 +269,9 @@ Foam::STARCDCoordinateRotation::STARCDCoordinateRotation
     const dictionary& dict,
     const objectRegistry&
 )
-{
-    vector rotation(dict.lookup("rotation"));
-
-    calcTransform
-    (
-        rotation.component(vector::X),
-        rotation.component(vector::Y),
-        rotation.component(vector::Z),
-        dict.lookupOrDefault("degrees", true)
-    );
-}
-
-
-Foam::STARCDCoordinateRotation::STARCDCoordinateRotation
-(
-    const STARCDCoordinateRotation& r
-)
 :
-    R_(r.R_),
+    STARCDCoordinateRotation(dict)
-    Rtr_(r.Rtr_)
 {}
 
 
-void Foam::STARCDCoordinateRotation::write(Ostream& os) const
-{
-     os.writeKeyword("e1") << e1() << token::END_STATEMENT << nl;
-     os.writeKeyword("e2") << e2() << token::END_STATEMENT << nl;
-     os.writeKeyword("e3") << e3() << token::END_STATEMENT << nl;
-}
-
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

src/meshTools/coordinateSystems/coordinateRotation/STARCDCoordinateRotation.H

@@ -3,7 +3,7 @@
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
     \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
-     \\/     M anipulation  |
+     \\/     M anipulation  | Copyright (C) 2017 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -83,7 +83,7 @@ class STARCDCoordinateRotation
             const scalar rotZ,
             const scalar rotX,
             const scalar rotY,
-            const bool inDegrees=true
+            const bool inDegrees
         );
 
 
@@ -98,11 +98,14 @@ public:
         //- Construct null
         STARCDCoordinateRotation();
 
+        //- Construct as copy
+        STARCDCoordinateRotation(const STARCDCoordinateRotation& r);
+
         //- Construct from rotation vector
         STARCDCoordinateRotation
         (
             const vector& rotZrotXrotY,
-            const bool inDegrees=true
+            const bool inDegrees
         );
 
         //- Construct from components of rotation vector
@@ -111,17 +114,19 @@ public:
             const scalar rotZ,
             const scalar rotX,
             const scalar rotY,
-            const bool inDegrees=true
+            const bool inDegrees
         );
 
         //- Construct from dictionary
-        STARCDCoordinateRotation(const dictionary&);
+        explicit STARCDCoordinateRotation(const dictionary& dict);
 
-        //- Construct from dictionary and mesh
+        //- Construct from dictionary and a registry (typically a mesh)
-        STARCDCoordinateRotation(const dictionary&, const objectRegistry&);
+        STARCDCoordinateRotation
+        (
+            const dictionary& dict,
+            const objectRegistry& unused
+        );
 
-         //- Construct as copy
-        STARCDCoordinateRotation(const STARCDCoordinateRotation&);
 
         //- Return clone
         autoPtr<coordinateRotation> clone() const
@@ -218,11 +223,6 @@ public:
         // symmetrical tensor
         virtual symmTensor transformVector(const vector& st) const;
 
-
-    // Write
-
-        //- Write
-        virtual void write(Ostream&) const;
 };
 
 

src/meshTools/coordinateSystems/coordinateRotation/axesRotation.C

@@ -3,7 +3,7 @@
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
     \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
-     \\/     M anipulation  |
+     \\/     M anipulation  | Copyright (C) 2017 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -32,7 +32,12 @@ License
 namespace Foam
 {
     defineTypeNameAndDebug(axesRotation, 0);
-    addToRunTimeSelectionTable(coordinateRotation, axesRotation, dictionary);
+    addToRunTimeSelectionTable
+    (
+        coordinateRotation,
+        axesRotation,
+        dictionary
+    );
     addToRunTimeSelectionTable
     (
         coordinateRotation,
@@ -41,119 +46,13 @@ namespace Foam
     );
 }
 
-// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
-
-void Foam::axesRotation::calcTransform
-(
-    const vector& axis1,
-    const vector& axis2,
-    const axisOrder& order
-)
-{
-    vector a = axis1/mag(axis1);
-    vector b = axis2;
-
-    b = b - (b & a)*a;
-
-    if (mag(b) < SMALL)
-    {
-        FatalErrorInFunction
-            << "axis1, axis2 appear co-linear: "
-            << axis1 << ", " << axis2 << endl
-            << abort(FatalError);
-    }
-
-    b = b/mag(b);
-    vector c = a^b;
-
-    tensor Rtr;
-    switch (order)
-    {
-        case e1e2:
-        {
-            Rtr = tensor(a, b, c);
-            break;
-        }
-        case e2e3:
-        {
-            Rtr = tensor(c, a, b);
-            break;
-        }
-        case e3e1:
-        {
-            Rtr = tensor(b, c, a);
-            break;
-        }
-        default:
-        {
-            FatalErrorInFunction
-                << "Unhandled axes specifictation" << endl
-                << abort(FatalError);
-
-            Rtr = Zero;
-            break;
-        }
-    }
-
-    // Global->local transformation
-    Rtr_ = Rtr;
-
-    // Local->global transformation
-    R_ = Rtr.T();
-}
-
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 Foam::axesRotation::axesRotation()
 :
     R_(sphericalTensor::I),
-    Rtr_(R_)
+    Rtr_(sphericalTensor::I)
-{}
-
-
-Foam::axesRotation::axesRotation
-(
-    const vector& axis,
-    const vector& dir
-)
-:
-    R_(sphericalTensor::I),
-    Rtr_(R_)
-{
-    calcTransform(axis, dir, e3e1);
-}
-
-
-Foam::axesRotation::axesRotation
-(
-    const dictionary& dict
-)
-:
-    R_(sphericalTensor::I),
-    Rtr_(R_)
-{
-    operator=(dict);
-}
-
-
-Foam::axesRotation::axesRotation
-(
-    const dictionary& dict,
-    const objectRegistry& obr
-)
-:
-    R_(sphericalTensor::I),
-    Rtr_(R_)
-{
-    operator=(dict);
-}
-
-
-Foam::axesRotation::axesRotation(const tensor& R)
-:
-    R_(R),
-    Rtr_(R_.T())
 {}
 
 
@@ -164,9 +63,145 @@ Foam::axesRotation::axesRotation(const axesRotation& r)
 {}
 
 
+Foam::axesRotation::axesRotation(const tensor& R)
+:
+    R_(R),
+    Rtr_(R_.T())
+{}
+
+
+Foam::axesRotation::axesRotation
+(
+    const vector& axis,
+    const vector& dir,
+    const axisOrder& order
+)
+:
+    R_(sphericalTensor::I),
+    Rtr_(sphericalTensor::I)
+{
+    setTransform(axis, dir, order);
+}
+
+
+Foam::axesRotation::axesRotation
+(
+    const vector& axis
+)
+:
+    R_(sphericalTensor::I),
+    Rtr_(sphericalTensor::I)
+{
+    direction maxCmpt = 0, dirCmpt = 1;
+
+    scalar maxVal = mag(axis[maxCmpt]);
+    bool negative = (axis[maxCmpt] < 0);
+
+    for (direction cmpt = 1; cmpt < vector::nComponents; ++cmpt)
+    {
+        const scalar val = mag(axis[cmpt]);
+
+        if (maxVal < val)
+        {
+            maxVal  = val;
+            maxCmpt = cmpt;
+            dirCmpt = maxCmpt+1;
+            negative = (axis[cmpt] < 0);
+
+            if (dirCmpt >= vector::nComponents)
+            {
+                dirCmpt = 0;
+            }
+        }
+    }
+
+    vector dir = Zero;
+    dir.component(dirCmpt) = (negative ? -1 : 1);
+
+    setTransform(axis, dir, E3_E1);
+}
+
+
+Foam::axesRotation::axesRotation
+(
+    const dictionary& dict
+)
+:
+    R_(sphericalTensor::I),
+    Rtr_(sphericalTensor::I)
+{
+    operator=(dict);
+}
+
+
+Foam::axesRotation::axesRotation
+(
+    const dictionary& dict,
+    const objectRegistry&
+)
+:
+    axesRotation(dict)
+{}
+
 
 // * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
 
+void Foam::axesRotation::setTransform
+(
+    const vector& axis1,
+    const vector& axis2,
+    const axisOrder& order
+)
+{
+    const vector a = axis1/mag(axis1);
+    vector b = axis2;
+
+    b = b - (b & a)*a;
+
+    if (mag(b) < SMALL)
+    {
+        FatalErrorInFunction
+            << "axis1, axis2 appear to be co-linear: "
+            << axis1 << ", " << axis2 << endl
+            << abort(FatalError);
+    }
+
+    b = b/mag(b);
+    const vector c = a^b;
+
+    // Global->local transformation
+    switch (order)
+    {
+        case E1_E2:
+        {
+            Rtr_ = tensor(a, b, c);
+            break;
+        }
+        case E2_E3:
+        {
+            Rtr_ = tensor(c, a, b);
+            break;
+        }
+        case E3_E1:
+        {
+            Rtr_ = tensor(b, c, a);
+            break;
+        }
+        default:
+        {
+            FatalErrorInFunction
+                << "Unhandled axes specification" << endl
+                << abort(FatalError);
+
+            break;
+        }
+    }
+
+    // Local->global transformation
+    R_ = Rtr_.T();
+}
+
+
 const Foam::tensorField& Foam::axesRotation::Tr() const
 {
     NotImplemented;
@@ -263,34 +298,28 @@ Foam::symmTensor Foam::axesRotation::transformVector
 
 void Foam::axesRotation::operator=(const dictionary& dict)
 {
-    if (debug)
-    {
-        Pout<< "axesRotation::operator=(const dictionary&) : "
-            << "assign from " << dict << endl;
-    }
-
     vector axis1, axis2;
-    axisOrder order(e3e1);
 
     if (dict.readIfPresent("e1", axis1) && dict.readIfPresent("e2", axis2))
     {
-        order = e1e2;
+        setTransform(axis1, axis2, E1_E2);
     }
-    else if (dict.readIfPresent("e2", axis1)&& dict.readIfPresent("e3", axis2))
+    else if (dict.readIfPresent("e2", axis1) && dict.readIfPresent("e3", axis2))
     {
-        order = e2e3;
+        setTransform(axis1, axis2, E2_E3);
     }
-    else if (dict.readIfPresent("e3", axis1)&& dict.readIfPresent("e1", axis2))
+    else if (dict.readIfPresent("e3", axis1) && dict.readIfPresent("e1", axis2))
     {
-        order = e3e1;
+        setTransform(axis1, axis2, E3_E1);
     }
     else if (dict.found("axis") || dict.found("direction"))
     {
         // Both "axis" and "direction" are required
         // If one is missing the appropriate error message will be generated
-        order = e3e1;
         dict.lookup("axis") >> axis1;
         dict.lookup("direction") >> axis2;
+
+        setTransform(axis1, axis2, E3_E1);
     }
     else
     {
@@ -299,16 +328,6 @@ void Foam::axesRotation::operator=(const dictionary& dict)
             << "found "
             << exit(FatalError);
     }
-
-    calcTransform(axis1, axis2, order);
-}
-
-
-void Foam::axesRotation::write(Ostream& os) const
-{
-     os.writeKeyword("e1") << e1() << token::END_STATEMENT << nl;
-     os.writeKeyword("e2") << e2() << token::END_STATEMENT << nl;
-     os.writeKeyword("e3") << e3() << token::END_STATEMENT << nl;
 }
 
 

src/meshTools/coordinateSystems/coordinateRotation/axesRotation.H

@@ -3,7 +3,7 @@
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
     \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
-     \\/     M anipulation  |
+     \\/     M anipulation  | Copyright (C) 2017 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -25,11 +25,12 @@ Class
     Foam::axesRotation
 
 Description
-    A coordinate rotation specified using global axis
+    A coordinate rotation specified using global axes
 
     The rotation is defined by a combination of vectors (e1/e2), (e2/e3)
     or (e3/e1). Any nonorthogonality will be absorbed into the second
-    vector.
+    vector. In terms of cylindrical coordinates, the 'axis' would
+    correspond to the \a z-axis and the 'direction' to the \a r-axis.
 
     \verbatim
     axesRotation
@@ -66,6 +67,20 @@ class axesRotation
 :
     public coordinateRotation
 {
+public:
+
+    //- The order/combination of local axes for the axes-rotation definition
+    //  Note that these follow the right-hand rule.
+    enum axisOrder
+    {
+        E1_E2,  //!< The axis is X-dominant, the direction is Y-dominant
+        E2_E3,  //!< The axis is Y-dominant, the direction is Z-dominant
+        E3_E1   //!< The axis is Z-dominant, the direction is X-dominant
+    };
+
+
+private:
+
     // Private data
 
         //- Local-to-Global transformation tensor
@@ -74,24 +89,6 @@ class axesRotation
         //- Global-to-Local transformation tensor
         tensor Rtr_;
 
-        //- The combination of local axes to be used
-        enum axisOrder
-        {
-            e1e2,
-            e2e3,
-            e3e1
-        };
-
-    // Private Member Functions
-
-        //- Calculate transformation tensor
-        void calcTransform
-        (
-            const vector& axis1,
-            const vector& axis2,
-            const axisOrder& order = e3e1
-        );
-
 
 public:
 
@@ -103,20 +100,34 @@ public:
         //- Construct null
         axesRotation();
 
-        //- Construct from 2 axes
+        //- Construct as copy
-        axesRotation(const vector& axis, const vector& dir);
+        axesRotation(const axesRotation& r);
+
+        //- Construct from local to global rotation matrix
+        explicit axesRotation(const tensor& R);
+
+        //- Construct from two axes (axis and direction)
+        axesRotation
+        (
+            const vector& axis,
+            const vector& dir,
+            const axisOrder& order = E3_E1
+        );
+
+        //- Construct from a single axis using a best-guess for the second axis
+        //  For the best-guess, the largest component value and sign of the
+        //  axis determines the direction orientation.
+        explicit axesRotation(const vector& axis);
 
         //- Construct from dictionary
-        axesRotation(const dictionary&);
+        explicit axesRotation(const dictionary& dict);
 
-        //- Construct from components
+        //- Construct from dictionary and a registry (typically a mesh)
-        axesRotation(const tensor& R);
+        axesRotation
-
+        (
-        //- Construct from dictionary and mesh
+            const dictionary& dict,
-        axesRotation(const dictionary&, const objectRegistry&);
+            const objectRegistry& unused
-
+        );
-        //- Construct as copy
-        axesRotation(const axesRotation&);
 
         //- Return clone
         autoPtr<coordinateRotation> clone() const
@@ -139,6 +150,14 @@ public:
             Rtr_ = sphericalTensor::I;
         }
 
+        //- Set the transformation tensors from two axes (axis and direction)
+        void setTransform
+        (
+            const vector& axis1,
+            const vector& axis2,
+            const axisOrder& order = E3_E1
+        );
+
         //- Update the rotation for a list of cells
         virtual void updateCells(const polyMesh&, const labelList&)
         {}
@@ -216,13 +235,8 @@ public:
     // Member Operators
 
         //- Assign from dictionary
-        void operator=(const dictionary&);
+        void operator=(const dictionary& dict);
 
-
-    // Write
-
-        //- Write
-        virtual void write(Ostream&) const;
 };
 
 

src/meshTools/coordinateSystems/coordinateRotation/coordinateRotation.C

@@ -3,7 +3,7 @@
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
     \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
-     \\/     M anipulation  |
+     \\/     M anipulation  | Copyright (C) 2017 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -74,4 +74,13 @@ Foam::symmTensor Foam::coordinateRotation::transformPrincipal
 
 }
 
+
+void Foam::coordinateRotation::write(Ostream& os) const
+{
+     os.writeEntry("e1", e1());
+     os.writeEntry("e2", e2());
+     os.writeEntry("e3", e3());
+}
+
+
 // ************************************************************************* //

src/meshTools/coordinateSystems/coordinateRotation/coordinateRotation.H

@@ -123,7 +123,8 @@ public:
         //- Select constructed from dictionary and objectRegistry
         static autoPtr<coordinateRotation> New
         (
-            const dictionary& dict, const objectRegistry& obr
+            const dictionary& dict,
+            const objectRegistry& obr
         );
 
         //- Select constructed from dictionary
@@ -216,8 +217,8 @@ public:
 
     // Write
 
-            //- Write
+        //- Write coordinateRotation as e1,e2,e3 vectors
-            virtual void write(Ostream&) const = 0;
+        virtual void write(Ostream& os) const;
 
 };
 

src/meshTools/coordinateSystems/coordinateRotation/coordinateRotationNew.C

@@ -3,7 +3,7 @@
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
     \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
-     \\/     M anipulation  |
+     \\/     M anipulation  | Copyright (C) 2017 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -34,26 +34,15 @@ Foam::autoPtr<Foam::coordinateRotation> Foam::coordinateRotation::New
     const objectRegistry& obr
 )
 {
-    if (debug)
-    {
-        Pout<< "coordinateRotation::New"
-            "(const dictionary&, const objectRegistry&) : "
-            << "constructing coordinateRotation"
-            << endl;
-    }
-
     const word rotType = dict.lookup("type");
 
-    objectRegistryConstructorTable::iterator cstrIter =
+    auto cstrIter = objectRegistryConstructorTablePtr_->cfind(rotType);
-        objectRegistryConstructorTablePtr_->find(rotType);
 
     if (!cstrIter.found())
     {
         FatalIOErrorInFunction(dict)
-            << "Unknown coordinateRotation type "
+            << "Unknown coordinateRotation type " << rotType << nl << nl
-            << rotType << nl << nl
             << "Valid coordinateRotation types are :" <<  nl
-            << "[default: axes ]"
             << objectRegistryConstructorTablePtr_->sortedToc()
             << exit(FatalIOError);
     }
@@ -67,23 +56,14 @@ Foam::autoPtr<Foam::coordinateRotation> Foam::coordinateRotation::New
     const dictionary& dict
 )
 {
-    if (debug)
-    {
-        Pout<< "coordinateRotation::New(const dictionary&) : "
-            << "constructing coordinateRotation"
-            << endl;
-    }
-
     const word rotType = dict.lookup("type");
 
-    dictionaryConstructorTable::iterator cstrIter =
+    auto cstrIter = dictionaryConstructorTablePtr_->cfind(rotType);
-        dictionaryConstructorTablePtr_->find(rotType);
 
     if (!cstrIter.found())
     {
         FatalIOErrorInFunction(dict)
-            << "Unknown coordinateRotation type "
+            << "Unknown coordinateRotation type " << rotType << nl << nl
-            << rotType << nl << nl
             << "Valid coordinateRotation types are :" <<  nl
             << dictionaryConstructorTablePtr_->sortedToc()
             << exit(FatalIOError);

src/meshTools/coordinateSystems/coordinateRotation/cylindrical.C

@@ -67,7 +67,7 @@ void Foam::cylindrical::init
         tensorField& R = Rptr_();
         forAll(cells, i)
         {
-            label celli = cells[i];
+            const label celli = cells[i];
             vector dir = cc[celli] - origin_;
             dir /= mag(dir) + VSMALL;
 
@@ -92,6 +92,40 @@ void Foam::cylindrical::init
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
+Foam::cylindrical::cylindrical(const cylindrical& r)
+:
+    Rptr_(r.Rptr_, false),  // clone
+    origin_(r.origin_),
+    e3_(r.e3_)
+{}
+
+
+Foam::cylindrical::cylindrical(const tensorField& R)
+:
+    Rptr_(),
+    origin_(Zero),
+    e3_(Zero)
+{
+    Rptr_() = R;
+}
+
+
+Foam::cylindrical::cylindrical(const dictionary& dict)
+:
+    Rptr_(),
+    origin_(),
+    e3_()
+{
+    FatalErrorInFunction
+        << " cylindrical can not be constructed from dictionary "
+        << " use the construtctor : "
+           "("
+           "    const dictionary&, const objectRegistry&"
+           ")"
+        << exit(FatalIOError);
+}
+
+
 Foam::cylindrical::cylindrical
 (
     const dictionary& dict,
@@ -103,12 +137,9 @@ Foam::cylindrical::cylindrical
     e3_(Zero)
 {
     // If origin is specified in the coordinateSystem
-    if (dict.parent().found("origin"))
+    dict.parent().readIfPresent("origin", origin_);
-    {
-        dict.parent().lookup("origin") >> origin_;
-    }
 
-    // rotation axis
+    // Rotation axis
     dict.lookup("e3") >> e3_;
 
     init(obr);
@@ -146,40 +177,6 @@ Foam::cylindrical::cylindrical
 }
 
 
-Foam::cylindrical::cylindrical(const dictionary& dict)
-:
-    Rptr_(),
-    origin_(),
-    e3_()
-{
-    FatalErrorInFunction
-        << " cylindrical can not be constructed from dictionary "
-        << " use the construtctor : "
-           "("
-           "    const dictionary&, const objectRegistry&"
-           ")"
-        << exit(FatalIOError);
-}
-
-
-Foam::cylindrical::cylindrical(const tensorField& R)
-:
-    Rptr_(),
-    origin_(Zero),
-    e3_(Zero)
-{
-    Rptr_() = R;
-}
-
-
-Foam::cylindrical::cylindrical(const cylindrical& r)
-:
-    Rptr_(r.Rptr_, false),  // clone
-    origin_(r.origin_),
-    e3_(r.e3_)
-{}
-
-
 
 // * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
 
@@ -360,7 +357,7 @@ Foam::symmTensor Foam::cylindrical::transformVector
 
 void Foam::cylindrical::write(Ostream& os) const
 {
-     os.writeKeyword("e3") << e3() << token::END_STATEMENT << nl;
+     os.writeEntry("e3", e3());
 }
 
 

src/meshTools/coordinateSystems/coordinateRotation/cylindrical.H

@@ -60,7 +60,7 @@ namespace Foam
 {
 
 /*---------------------------------------------------------------------------*\
-                    Class cylindrical Declaration
+                         Class cylindrical Declaration
 \*---------------------------------------------------------------------------*/
 
 class cylindrical
@@ -96,13 +96,22 @@ public:
 
     // Constructors
 
+        //- Construct as copy
+        cylindrical(const cylindrical& r);
+
+        //- Construct from tensor Field
+        explicit cylindrical(const tensorField& R);
+
+        //- Construct from dictionary - for API compatibility only
+        explicit cylindrical(const dictionary& dict);
+
         //- Construct from dictionary and objectRegistry
-        cylindrical(const dictionary&, const objectRegistry&);
+        cylindrical(const dictionary& dict, const objectRegistry& obr);
 
         //- Construct from components for all cells
         cylindrical
         (
-            const objectRegistry&,
+            const objectRegistry& obr,
             const vector& axis,
             const point& origin
         );
@@ -110,21 +119,12 @@ public:
         //- Construct from components for list of cells
         cylindrical
         (
-            const objectRegistry&,
+            const objectRegistry& obr,
             const vector& axis,
             const point& origin,
             const List<label>& cells
         );
 
-        //- Construct from dictionary
-        cylindrical(const dictionary&);
-
-        //- Construct from tensor Field
-        cylindrical(const tensorField&);
-
-        //- Construct as copy
-        cylindrical(const cylindrical&);
-
         //- Return clone
         autoPtr<coordinateRotation> clone() const
         {

src/meshTools/coordinateSystems/coordinateSystemNew.C

@@ -35,10 +35,9 @@ Foam::autoPtr<Foam::coordinateSystem> Foam::coordinateSystem::New
 )
 {
     const dictionary& coordDict = dict.subDict(typeName_());
-    word coordType = coordDict.lookup("type");
+    const word coordType = coordDict.lookup("type");
 
-    dictionaryConstructorTable::iterator cstrIter =
+    auto cstrIter = dictionaryConstructorTablePtr_->cfind(coordType);
-        dictionaryConstructorTablePtr_->find(coordType);
 
     if (!cstrIter.found())
     {

src/thermophysicalModels/basic/heThermo/heThermo.C

@@ -3,7 +3,7 @@
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
     \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
-     \\/     M anipulation  |
+     \\/     M anipulation  | Copyright (C) 2015-2017 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -52,11 +52,16 @@ heBoundaryCorrection(volScalarField& h)
 
 
 template<class BasicThermo, class MixtureType>
-void Foam::heThermo<BasicThermo, MixtureType>::init()
+void Foam::heThermo<BasicThermo, MixtureType>::init
+(
+    const volScalarField& p,
+    const volScalarField& T,
+    volScalarField& he
+)
 {
-    scalarField& heCells = he_.primitiveFieldRef();
+    scalarField& heCells = he.primitiveFieldRef();
-    const scalarField& pCells = this->p_;
+    const scalarField& pCells = p.primitiveField();
-    const scalarField& TCells = this->T_;
+    const scalarField& TCells = T.primitiveField();
 
     forAll(heCells, celli)
     {
@@ -64,19 +69,25 @@ void Foam::heThermo<BasicThermo, MixtureType>::init()
             this->cellMixture(celli).HE(pCells[celli], TCells[celli]);
     }
 
-    volScalarField::Boundary& heBf = he_.boundaryFieldRef();
+    volScalarField::Boundary& heBf = he.boundaryFieldRef();
 
     forAll(heBf, patchi)
     {
-        heBf[patchi] == he
+        heBf[patchi] == this->he
         (
-            this->p_.boundaryField()[patchi],
+            p.boundaryField()[patchi],
-            this->T_.boundaryField()[patchi],
+            T.boundaryField()[patchi],
             patchi
         );
     }
 
-    this->heBoundaryCorrection(he_);
+    this->heBoundaryCorrection(he);
+
+    // Note: T does not have oldTime
+    if (p.nOldTimes() > 0)
+    {
+        init(p.oldTime(), T.oldTime(), he.oldTime());
+    }
 }
 
 
@@ -112,7 +123,7 @@ Foam::heThermo<BasicThermo, MixtureType>::heThermo
         this->heBoundaryBaseTypes()
     )
 {
-    init();
+    init(this->p_, this->T_, he_);
 }
 
 
@@ -146,7 +157,7 @@ Foam::heThermo<BasicThermo, MixtureType>::heThermo
         this->heBoundaryBaseTypes()
     )
 {
-    init();
+    init(this->p_, this->T_, he_);
 }
 
 

src/thermophysicalModels/basic/heThermo/heThermo.H

@@ -3,7 +3,7 @@
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
     \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
-     \\/     M anipulation  |
+     \\/     M anipulation  | Copyright (C) 2015-2017 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -77,7 +77,12 @@ private:
 
 
         //- Initialize heThermo
-        void init();
+        void init
+        (
+            const volScalarField& p,
+            const volScalarField& T,
+            volScalarField& he
+        );
 
 
 public:

src/thermophysicalModels/basic/psiThermo/hePsiThermo.C

@@ -3,7 +3,7 @@
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
     \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
-     \\/     M anipulation  |
+     \\/     M anipulation  | Copyright (C) 2015-2017 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -28,15 +28,40 @@ License
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
 template<class BasicPsiThermo, class MixtureType>
-void Foam::hePsiThermo<BasicPsiThermo, MixtureType>::calculate()
+void Foam::hePsiThermo<BasicPsiThermo, MixtureType>::calculate
+(
+    const volScalarField& p,
+    volScalarField& T,
+    volScalarField& he,
+    volScalarField& psi,
+    volScalarField& mu,
+    volScalarField& alpha,
+    const bool doOldTimes
+)
 {
-    const scalarField& hCells = this->he_;
+    // Note: update oldTimes before current time so that if T.oldTime() is
-    const scalarField& pCells = this->p_;
+    // created from T, it starts from the unconverted T
+    if (doOldTimes && (p.nOldTimes() || T.nOldTimes()))
+    {
+        calculate
+        (
+            p.oldTime(),
+            T.oldTime(),
+            he.oldTime(),
+            psi.oldTime(),
+            mu.oldTime(),
+            alpha.oldTime(),
+            true
+        );
+    }
 
-    scalarField& TCells = this->T_.primitiveFieldRef();
+    const scalarField& hCells = he.primitiveField();
-    scalarField& psiCells = this->psi_.primitiveFieldRef();
+    const scalarField& pCells = p.primitiveField();
-    scalarField& muCells = this->mu_.primitiveFieldRef();
+
-    scalarField& alphaCells = this->alpha_.primitiveFieldRef();
+    scalarField& TCells = T.primitiveFieldRef();
+    scalarField& psiCells = psi.primitiveFieldRef();
+    scalarField& muCells = mu.primitiveFieldRef();
+    scalarField& alphaCells = alpha.primitiveFieldRef();
 
     forAll(TCells, celli)
     {
@@ -56,27 +81,16 @@ void Foam::hePsiThermo<BasicPsiThermo, MixtureType>::calculate()
         alphaCells[celli] = mixture_.alphah(pCells[celli], TCells[celli]);
     }
 
-    volScalarField::Boundary& pBf =
+    const volScalarField::Boundary& pBf = p.boundaryField();
-        this->p_.boundaryFieldRef();
+    volScalarField::Boundary& TBf = T.boundaryFieldRef();
+    volScalarField::Boundary& psiBf = psi.boundaryFieldRef();
+    volScalarField::Boundary& heBf = he.boundaryFieldRef();
+    volScalarField::Boundary& muBf = mu.boundaryFieldRef();
+    volScalarField::Boundary& alphaBf = alpha.boundaryFieldRef();
 
-    volScalarField::Boundary& TBf =
+    forAll(pBf, patchi)
-        this->T_.boundaryFieldRef();
-
-    volScalarField::Boundary& psiBf =
-        this->psi_.boundaryFieldRef();
-
-    volScalarField::Boundary& heBf =
-        this->he().boundaryFieldRef();
-
-    volScalarField::Boundary& muBf =
-        this->mu_.boundaryFieldRef();
-
-    volScalarField::Boundary& alphaBf =
-        this->alpha_.boundaryFieldRef();
-
-    forAll(this->T_.boundaryField(), patchi)
     {
-        fvPatchScalarField& pp = pBf[patchi];
+        const fvPatchScalarField& pp = pBf[patchi];
         fvPatchScalarField& pT = TBf[patchi];
         fvPatchScalarField& ppsi = psiBf[patchi];
         fvPatchScalarField& phe = heBf[patchi];
@@ -126,10 +140,16 @@ Foam::hePsiThermo<BasicPsiThermo, MixtureType>::hePsiThermo
 :
     heThermo<BasicPsiThermo, MixtureType>(mesh, phaseName)
 {
-    calculate();
+    calculate
-
+    (
-    // Switch on saving old time
+        this->p_,
-    this->psi_.oldTime();
+        this->T_,
+        this->he_,
+        this->psi_,
+        this->mu_,
+        this->alpha_,
+        true                    // Create old time fields
+    );
 }
 
 
@@ -145,20 +165,20 @@ Foam::hePsiThermo<BasicPsiThermo, MixtureType>::~hePsiThermo()
 template<class BasicPsiThermo, class MixtureType>
 void Foam::hePsiThermo<BasicPsiThermo, MixtureType>::correct()
 {
-    if (debug)
+    DebugInFunction << endl;
-    {
-        InfoInFunction << endl;
-    }
 
-    // force the saving of the old-time values
+    calculate
-    this->psi_.oldTime();
+    (
+        this->p_,
+        this->T_,
+        this->he_,
+        this->psi_,
+        this->mu_,
+        this->alpha_,
+        false           // No need to update old times
+    );
 
-    calculate();
+    DebugInFunction << "Finished" << endl;
-
-    if (debug)
-    {
-        Info<< "    Finished" << endl;
-    }
 }
 
 

src/thermophysicalModels/basic/psiThermo/hePsiThermo.H

@@ -3,7 +3,7 @@
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
     \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
-     \\/     M anipulation  |
+     \\/     M anipulation  | Copyright (C) 2015-2017 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -55,11 +55,21 @@ class hePsiThermo
     // Private Member Functions
 
         //- Calculate the thermo variables
-        void calculate();
+        void calculate
+        (
+            const volScalarField& p,
+            volScalarField& T,
+            volScalarField& he,
+            volScalarField& psi,
+            volScalarField& mu,
+            volScalarField& alpha,
+            const bool doOldTimes
+        );
 
         //- Construct as copy (not implemented)
         hePsiThermo(const hePsiThermo<BasicPsiThermo, MixtureType>&);
 
+
 public:
 
     //- Runtime type information

src/thermophysicalModels/basic/rhoThermo/heRhoThermo.C

@@ -3,7 +3,7 @@
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
     \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
-     \\/     M anipulation  |
+     \\/     M anipulation  | Copyright (C) 2015-2017 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -28,16 +28,43 @@ License
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
 template<class BasicPsiThermo, class MixtureType>
-void Foam::heRhoThermo<BasicPsiThermo, MixtureType>::calculate()
+void Foam::heRhoThermo<BasicPsiThermo, MixtureType>::calculate
+(
+    const volScalarField& p,
+    volScalarField& T,
+    volScalarField& he,
+    volScalarField& psi,
+    volScalarField& rho,
+    volScalarField& mu,
+    volScalarField& alpha,
+    const bool doOldTimes
+)
 {
-    const scalarField& hCells = this->he();
+    // Note: update oldTimes before current time so that if T.oldTime() is
-    const scalarField& pCells = this->p_;
+    // created from T, it starts from the unconverted T
+    if (doOldTimes && (p.nOldTimes() || T.nOldTimes()))
+    {
+        calculate
+        (
+            p.oldTime(),
+            T.oldTime(),
+            he.oldTime(),
+            psi.oldTime(),
+            rho.oldTime(),
+            mu.oldTime(),
+            alpha.oldTime(),
+            true
+        );
+    }
 
-    scalarField& TCells = this->T_.primitiveFieldRef();
+    const scalarField& hCells = he.primitiveField();
-    scalarField& psiCells = this->psi_.primitiveFieldRef();
+    const scalarField& pCells = p.primitiveField();
-    scalarField& rhoCells = this->rho_.primitiveFieldRef();
+
-    scalarField& muCells = this->mu_.primitiveFieldRef();
+    scalarField& TCells = T.primitiveFieldRef();
-    scalarField& alphaCells = this->alpha_.primitiveFieldRef();
+    scalarField& psiCells = psi.primitiveFieldRef();
+    scalarField& rhoCells = rho.primitiveFieldRef();
+    scalarField& muCells = mu.primitiveFieldRef();
+    scalarField& alphaCells = alpha.primitiveFieldRef();
 
     forAll(TCells, celli)
     {
@@ -58,30 +85,17 @@ void Foam::heRhoThermo<BasicPsiThermo, MixtureType>::calculate()
         alphaCells[celli] = mixture_.alphah(pCells[celli], TCells[celli]);
     }
 
-    volScalarField::Boundary& pBf =
+    const volScalarField::Boundary& pBf = p.boundaryField();
-        this->p_.boundaryFieldRef();
+    volScalarField::Boundary& TBf = T.boundaryFieldRef();
+    volScalarField::Boundary& psiBf = psi.boundaryFieldRef();
+    volScalarField::Boundary& rhoBf = rho.boundaryFieldRef();
+    volScalarField::Boundary& heBf = he.boundaryFieldRef();
+    volScalarField::Boundary& muBf = mu.boundaryFieldRef();
+    volScalarField::Boundary& alphaBf = alpha.boundaryFieldRef();
 
-    volScalarField::Boundary& TBf =
+    forAll(pBf, patchi)
-        this->T_.boundaryFieldRef();
-
-    volScalarField::Boundary& psiBf =
-        this->psi_.boundaryFieldRef();
-
-    volScalarField::Boundary& rhoBf =
-        this->rho_.boundaryFieldRef();
-
-    volScalarField::Boundary& heBf =
-        this->he().boundaryFieldRef();
-
-    volScalarField::Boundary& muBf =
-        this->mu_.boundaryFieldRef();
-
-    volScalarField::Boundary& alphaBf =
-        this->alpha_.boundaryFieldRef();
-
-    forAll(this->T_.boundaryField(), patchi)
     {
-        fvPatchScalarField& pp = pBf[patchi];
+        const fvPatchScalarField& pp = pBf[patchi];
         fvPatchScalarField& pT = TBf[patchi];
         fvPatchScalarField& ppsi = psiBf[patchi];
         fvPatchScalarField& prho = rhoBf[patchi];
@@ -134,7 +148,17 @@ Foam::heRhoThermo<BasicPsiThermo, MixtureType>::heRhoThermo
 :
     heThermo<BasicPsiThermo, MixtureType>(mesh, phaseName)
 {
-    calculate();
+    calculate
+    (
+        this->p_,
+        this->T_,
+        this->he_,
+        this->psi_,
+        this->rho_,
+        this->mu_,
+        this->alpha_,
+        true                    // Create old time fields
+    );
 }
 
 
@@ -150,17 +174,21 @@ Foam::heRhoThermo<BasicPsiThermo, MixtureType>::~heRhoThermo()
 template<class BasicPsiThermo, class MixtureType>
 void Foam::heRhoThermo<BasicPsiThermo, MixtureType>::correct()
 {
-    if (debug)
+    DebugInFunction << endl;
-    {
-        InfoInFunction << endl;
-    }
 
-    calculate();
+    calculate
+    (
+        this->p_,
+        this->T_,
+        this->he_,
+        this->psi_,
+        this->rho_,
+        this->mu_,
+        this->alpha_,
+        false           // No need to update old times
+    );
 
-    if (debug)
+    DebugInFunction << "Finished" << endl;
-    {
-        Info<< "    Finished" << endl;
-    }
 }
 
 

src/thermophysicalModels/basic/rhoThermo/heRhoThermo.H

@@ -3,7 +3,7 @@
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
     \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
-     \\/     M anipulation  |
+     \\/     M anipulation  | Copyright (C) 2015-2017 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -55,7 +55,17 @@ class heRhoThermo
     // Private Member Functions
 
         //- Calculate the thermo variables
-        void calculate();
+        void calculate
+        (
+            const volScalarField& p,
+            volScalarField& T,
+            volScalarField& he,
+            volScalarField& psi,
+            volScalarField& rho,
+            volScalarField& mu,
+            volScalarField& alpha,
+            const bool doOldTimes
+        );
 
         //- Construct as copy (not implemented)
         heRhoThermo(const heRhoThermo<BasicPsiThermo, MixtureType>&);