Merge commit 'OpenCFD/master' into olesenm

.gitignore

@@ -41,6 +41,10 @@ SunOS*Gcc*/
 # reinstate wmake/rules that might look like build folders
 !wmake/rules/*/
 
+# but do continue to ignore the derived wmake files
+wmake/rules/*/dirToString
+wmake/rules/*/wmkdep
+
 # doxygen generated documentation
 doc/[Dd]oxygen/html
 doc/[Dd]oxygen/latex

applications/solvers/combustion/PDRFoam/PDRFoam.C

@@ -75,7 +75,7 @@ int main(int argc, char *argv[])
     #include "createTime.H"
     #include "createMesh.H"
     #include "readCombustionProperties.H"
-    #include "readEnvironmentalProperties.H"
+    #include "readGravitationalAcceleration.H"
     #include "createFields.H"
     #include "initContinuityErrs.H"
     #include "readTimeControls.H"

applications/solvers/combustion/PDRFoam/PDRFoamAutoRefine.C

@@ -77,7 +77,7 @@ int main(int argc, char *argv[])
     #include "createTime.H"
     #include "createDynamicFvMesh.H"
     #include "readCombustionProperties.H"
-#   include "readEnvironmentalProperties.H"
+    #include "readGravitationalAcceleration.H"
     #include "createFields.H"
     #include "readPISOControls.H"
     #include "initContinuityErrs.H"
@@ -86,7 +86,7 @@ int main(int argc, char *argv[])
 
     scalar StCoNum = 0.0;
 
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+    // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
     Info<< "\nStarting time loop\n" << endl;
 

applications/solvers/combustion/XiFoam/XiFoam.C

@@ -66,7 +66,7 @@ int main(int argc, char *argv[])
     #include "createTime.H"
     #include "createMesh.H"
     #include "readCombustionProperties.H"
-    #include "readEnvironmentalProperties.H"
+    #include "readGravitationalAcceleration.H"
     #include "createFields.H"
     #include "initContinuityErrs.H"
     #include "readTimeControls.H"

applications/solvers/combustion/dieselEngineFoam/createSpray.H

@@ -23,7 +23,7 @@ spray dieselSpray
     composition,
     gasProperties,
     thermo,
-    environmentalProperties
+    g
 );
 
 scalar gasMass0 = fvc::domainIntegrate(rho).value();

applications/solvers/combustion/dieselEngineFoam/dieselEngineFoam.C

@@ -52,7 +52,7 @@ int main(int argc, char *argv[])
     #include "createEngineTime.H"
     #include "createEngineMesh.H"
     #include "createFields.H"
-    #include "readEnvironmentalProperties.H"
+    #include "readGravitationalAcceleration.H"
     #include "readCombustionProperties.H"
     #include "createSpray.H"
     #include "initContinuityErrs.H"
@@ -124,7 +124,10 @@ int main(int argc, char *argv[])
 
         rho = thermo.rho();
 
-        runTime.write();
+        if (runTime.write())
+        {
+            chemistry.dQ()().write();
+        }
 
         Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
             << "  ClockTime = " << runTime.elapsedClockTime() << " s"

applications/solvers/combustion/dieselFoam/dieselFoam.C

@@ -50,7 +50,7 @@ int main(int argc, char *argv[])
     #include "createTime.H"
     #include "createMesh.H"
     #include "createFields.H"
-    #include "readEnvironmentalProperties.H"
+    #include "readGravitationalAcceleration.H"
     #include "readCombustionProperties.H"
     #include "createSpray.H"
     #include "initContinuityErrs.H"
@@ -114,7 +114,10 @@ int main(int argc, char *argv[])
 
         rho = thermo.rho();
 
-        runTime.write();
+        if (runTime.write())
+        {
+            chemistry.dQ()().write();
+        }
 
         Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
             << "  ClockTime = " << runTime.elapsedClockTime() << " s"

applications/solvers/combustion/reactingFoam/reactingFoam.C

@@ -45,7 +45,7 @@ int main(int argc, char *argv[])
     #include "createTime.H"
     #include "createMesh.H"
     #include "readChemistryProperties.H"
-    #include "readEnvironmentalProperties.H"
+    #include "readGravitationalAcceleration.H"
     #include "createFields.H"
     #include "initContinuityErrs.H"
     #include "readTimeControls.H"
@@ -86,7 +86,10 @@ int main(int argc, char *argv[])
 
         turbulence->correct();
 
-        rho = thermo.rho();
+        if (runTime.write())
+        {
+            chemistry.dQ()().write();
+        }
 
         runTime.write();
 

applications/solvers/combustion/rhoReactingFoam/Make/files

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

applications/solvers/combustion/rhoReactingFoam/Make/options

@@ -0,0 +1,19 @@
+EXE_INC = \
+    -I../XiFoam \
+    -I$(LIB_SRC)/turbulenceModels/compressible/turbulenceModel \
+    -I$(LIB_SRC)/thermophysicalModels/specie/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/reactionThermo/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/chemistryModel/lnInclude \
+    -I$(LIB_SRC)/ODE/lnInclude \
+    -I$(LIB_SRC)/finiteVolume/lnInclude
+
+EXE_LIBS = \
+    -lcompressibleRASModels \
+    -lcompressibleLESModels \
+    -lreactionThermophysicalModels \
+    -lspecie \
+    -lbasicThermophysicalModels \
+    -lchemistryModel \
+    -lODE \
+    -lfiniteVolume

applications/solvers/combustion/rhoReactingFoam/YEqn.H

@@ -0,0 +1,43 @@
+tmp<fv::convectionScheme<scalar> > mvConvection
+(
+    fv::convectionScheme<scalar>::New
+    (
+        mesh,
+        fields,
+        phi,
+        mesh.divScheme("div(phi,Yi_h)")
+    )
+);
+
+{
+    label inertIndex = -1;
+    volScalarField Yt = 0.0*Y[0];
+
+    for (label i=0; i<Y.size(); i++)
+    {
+        if (Y[i].name() != inertSpecie)
+        {
+            volScalarField& Yi = Y[i];
+
+            solve
+            (
+                fvm::ddt(rho, Yi)
+              + mvConvection->fvmDiv(phi, Yi)
+              - fvm::laplacian(turbulence->muEff(), Yi)
+             ==
+                kappa*chemistry.RR(i),
+                mesh.solver("Yi")
+            );
+
+            Yi.max(0.0);
+            Yt += Yi;
+        }
+        else
+        {
+            inertIndex = i;
+        }
+    }
+
+    Y[inertIndex] = scalar(1) - Yt;
+    Y[inertIndex].max(0.0);
+}

applications/solvers/combustion/rhoReactingFoam/chemistry.H

@@ -0,0 +1,24 @@
+{
+    Info << "Solving chemistry" << endl;
+
+    chemistry.solve
+    (
+        runTime.value() - runTime.deltaT().value(),
+        runTime.deltaT().value()
+    );
+
+    // turbulent time scale
+    if (turbulentReaction)
+    {
+        volScalarField tk =
+                Cmix*sqrt(turbulence->muEff()/rho/turbulence->epsilon());
+        volScalarField tc = chemistry.tc();
+
+        // Chalmers PaSR model
+        kappa = (runTime.deltaT() + tc)/(runTime.deltaT() + tc + tk);
+    }
+    else
+    {
+        kappa = 1.0;
+    }
+}

applications/solvers/combustion/rhoReactingFoam/createFields.H

@@ -0,0 +1,85 @@
+Info<< nl << "Reading thermophysicalProperties" << endl;
+autoPtr<rhoChemistryModel> pChemistry
+(
+    rhoChemistryModel::New(mesh)
+);
+rhoChemistryModel& chemistry = pChemistry();
+
+hReactionThermo& thermo = chemistry.thermo();
+
+basicMultiComponentMixture& composition = thermo.composition();
+PtrList<volScalarField>& Y = composition.Y();
+
+word inertSpecie(thermo.lookup("inertSpecie"));
+
+volScalarField rho
+(
+    IOobject
+    (
+        "rho",
+        runTime.timeName(),
+        mesh
+    ),
+    thermo.rho()
+);
+
+Info<< "Reading field U\n" << endl;
+volVectorField U
+(
+    IOobject
+    (
+        "U",
+        runTime.timeName(),
+        mesh,
+        IOobject::MUST_READ,
+        IOobject::AUTO_WRITE
+    ),
+    mesh
+);
+
+
+volScalarField& p = thermo.p();
+const volScalarField& psi = thermo.psi();
+volScalarField& h = thermo.h();
+
+
+#include "compressibleCreatePhi.H"
+
+volScalarField kappa
+(
+    IOobject
+    (
+        "kappa",
+        runTime.timeName(),
+        mesh,
+        IOobject::NO_READ,
+        IOobject::AUTO_WRITE
+    ),
+    mesh,
+    dimensionedScalar("zero", dimless, 0.0)
+);
+
+Info << "Creating turbulence model.\n" << nl;
+autoPtr<compressible::turbulenceModel> turbulence
+(
+    compressible::turbulenceModel::New
+    (
+        rho,
+        U,
+        phi,
+        thermo
+    )
+);
+
+Info<< "Creating field DpDt\n" << endl;
+volScalarField DpDt =
+    fvc::DDt(surfaceScalarField("phiU", phi/fvc::interpolate(rho)), p);
+
+multivariateSurfaceInterpolationScheme<scalar>::fieldTable fields;
+
+forAll (Y, i)
+{
+    fields.add(Y[i]);
+}
+fields.add(h);
+

applications/solvers/combustion/rhoReactingFoam/pEqn.H

@@ -0,0 +1,93 @@
+{
+    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 rUA = 1.0/UEqn.A();
+    U = rUA*UEqn.H();
+
+    if (transonic)
+    {
+        surfaceScalarField phiv =
+            (fvc::interpolate(U) & mesh.Sf())
+          + fvc::ddtPhiCorr(rUA, rho, U, phi);
+
+        phi = fvc::interpolate(rho)*phiv;
+
+        surfaceScalarField phid
+        (
+            "phid",
+            fvc::interpolate(thermo.psi())*phiv
+        );
+
+        for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
+        {
+            fvScalarMatrix pEqn
+            (
+                fvc::ddt(rho) + fvc::div(phi)
+              + correction(fvm::ddt(psi, p) + fvm::div(phid, p))
+              - fvm::laplacian(rho*rUA, p)
+            );
+
+            if (ocorr == nOuterCorr && corr == nCorr && nonOrth == nNonOrthCorr)
+            {
+                pEqn.solve(mesh.solver(p.name() + "Final"));
+            }
+            else
+            {
+                pEqn.solve();
+            }
+
+            if (nonOrth == nNonOrthCorr)
+            {
+                phi += pEqn.flux();
+            }
+        }
+    }
+    else
+    {
+        phi =
+            fvc::interpolate(rho)
+           *(
+                (fvc::interpolate(U) & mesh.Sf())
+              + fvc::ddtPhiCorr(rUA, rho, U, phi)
+            );
+
+        for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
+        {
+            fvScalarMatrix pEqn
+            (
+                fvc::ddt(rho) + psi*correction(fvm::ddt(p))
+              + fvc::div(phi)
+              - fvm::laplacian(rho*rUA, p)
+            );
+
+            if (ocorr == nOuterCorr && corr == nCorr && nonOrth == nNonOrthCorr)
+            {
+                pEqn.solve(mesh.solver(p.name() + "Final"));
+            }
+            else
+            {
+                pEqn.solve();
+            }
+
+            if (nonOrth == nNonOrthCorr)
+            {
+                phi += pEqn.flux();
+            }
+        }
+    }
+
+    // Second part of thermodynamic density update
+    thermo.rho() += psi*p;
+
+    #include "rhoEqn.H"
+    #include "compressibleContinuityErrs.H"
+
+    U -= rUA*fvc::grad(p);
+    U.correctBoundaryConditions();
+
+    DpDt = fvc::DDt(surfaceScalarField("phiU", phi/fvc::interpolate(rho)), p);
+}

applications/solvers/combustion/rhoReactingFoam/readChemistryProperties.H

@@ -0,0 +1,23 @@
+Info<< "Reading chemistry properties\n" << endl;
+
+IOdictionary chemistryProperties
+(
+    IOobject
+    (
+        "chemistryProperties",
+        runTime.constant(),
+        mesh,
+        IOobject::MUST_READ,
+        IOobject::NO_WRITE,
+        false
+    )
+);
+
+Switch turbulentReaction(chemistryProperties.lookup("turbulentReaction"));
+
+dimensionedScalar Cmix("Cmix", dimless, 1.0);
+
+if (turbulentReaction)
+{
+    chemistryProperties.lookup("Cmix") >> Cmix;
+}

applications/solvers/combustion/rhoReactingFoam/rhoReactingFoam.C

@@ -0,0 +1,105 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2009-2009 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 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, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+Application
+    rhoReactingFoam
+
+Description
+    Chemical reaction code using density based thermodynamics package.
+
+\*---------------------------------------------------------------------------*/
+
+#include "fvCFD.H"
+#include "hReactionThermo.H"
+#include "turbulenceModel.H"
+#include "rhoChemistryModel.H"
+#include "chemistrySolver.H"
+#include "multivariateScheme.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+int main(int argc, char *argv[])
+{
+    #include "setRootCase.H"
+    #include "createTime.H"
+    #include "createMesh.H"
+    #include "readChemistryProperties.H"
+    #include "readGravitationalAcceleration.H"
+    #include "createFields.H"
+    #include "initContinuityErrs.H"
+    #include "readTimeControls.H"
+    #include "compressibleCourantNo.H"
+    #include "setInitialDeltaT.H"
+
+    // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+    Info << "\nStarting time loop\n" << endl;
+
+    while (runTime.run())
+    {
+        #include "readTimeControls.H"
+        #include "readPISOControls.H"
+        #include "compressibleCourantNo.H"
+        #include "setDeltaT.H"
+
+        runTime++;
+        Info<< "Time = " << runTime.timeName() << nl << endl;
+
+        #include "chemistry.H"
+        #include "rhoEqn.H"
+        #include "UEqn.H"
+
+        for (label ocorr=1; ocorr <= nOuterCorr; ocorr++)
+        {
+            #include "YEqn.H"
+            #include "hEqn.H"
+
+            // --- PISO loop
+            for (int corr=1; corr<=nCorr; corr++)
+            {
+                #include "pEqn.H"
+            }
+        }
+
+        turbulence->correct();
+
+        rho = thermo.rho();
+
+        if (runTime.write())
+        {
+            chemistry.dQ()().write();
+        }
+
+        Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
+            << "  ClockTime = " << runTime.elapsedClockTime() << " s"
+            << nl << endl;
+    }
+
+    Info<< "End\n" << endl;
+
+    return 0;
+}
+
+
+// ************************************************************************* //

applications/solvers/compressible/rhoCentralFoam/Make/options

@@ -8,5 +8,4 @@ EXE_LIBS = \
     -lfiniteVolume \
     -lbasicThermophysicalModels \
     -lspecie \
-    -L$(FOAM_USER_LIBBIN) \
     -lrhoCentralFoam

applications/solvers/compressible/rhoPimpleFoam/pEqn.H

@@ -58,8 +58,6 @@ else
         //+ fvc::ddtPhiCorr(rUA, rho, U, phi)
         );
 
-    //bool closedVolume = adjustPhi(phi, U, p);
-
     for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
     {
         // Pressure corrector

applications/solvers/heatTransfer/buoyantBoussinesqPisoFoam/buoyantBoussinesqPisoFoam.C

@@ -57,7 +57,7 @@ int main(int argc, char *argv[])
     #include "setRootCase.H"
     #include "createTime.H"
     #include "createMesh.H"
-    #include "readEnvironmentalProperties.H"
+    #include "readGravitationalAcceleration.H"
     #include "createFields.H"
     #include "initContinuityErrs.H"
     #include "readTimeControls.H"

applications/solvers/heatTransfer/buoyantBoussinesqSimpleFoam/buoyantBoussinesqSimpleFoam.C

@@ -57,7 +57,7 @@ int main(int argc, char *argv[])
     #include "setRootCase.H"
     #include "createTime.H"
     #include "createMesh.H"
-    #include "readEnvironmentalProperties.H"
+    #include "readGravitationalAcceleration.H"
     #include "createFields.H"
     #include "initContinuityErrs.H"
 

applications/solvers/heatTransfer/buoyantBoussinesqSimpleFoam/pEqn.H

@@ -7,6 +7,7 @@
 
     phi = fvc::interpolate(U) & mesh.Sf();
     adjustPhi(phi, U, p);
+
     surfaceScalarField buoyancyPhi =
         rUAf*fvc::interpolate(rhok)*(g & mesh.Sf());
     phi += buoyancyPhi;

applications/solvers/heatTransfer/buoyantPisoFoam/buoyantPisoFoam.C

@@ -44,7 +44,7 @@ int main(int argc, char *argv[])
     #include "setRootCase.H"
     #include "createTime.H"
     #include "createMesh.H"
-    #include "readEnvironmentalProperties.H"
+    #include "readGravitationalAcceleration.H"
     #include "createFields.H"
     #include "initContinuityErrs.H"
     #include "readTimeControls.H"

applications/solvers/heatTransfer/buoyantSimpleFoam/buoyantSimpleFoam.C

@@ -42,7 +42,7 @@ int main(int argc, char *argv[])
     #include "setRootCase.H"
     #include "createTime.H"
     #include "createMesh.H"
-    #include "readEnvironmentalProperties.H"
+    #include "readGravitationalAcceleration.H"
     #include "createFields.H"
     #include "initContinuityErrs.H"
 

applications/solvers/heatTransfer/buoyantSimpleFoam/pEqn.H

@@ -9,6 +9,7 @@
 
     phi = fvc::interpolate(rho)*(fvc::interpolate(U) & mesh.Sf());
     bool closedVolume = adjustPhi(phi, U, p);
+
     surfaceScalarField buoyancyPhi =
         rhorUAf*fvc::interpolate(rho)*(g & mesh.Sf());
     phi += buoyancyPhi;

applications/solvers/heatTransfer/buoyantSimpleRadiationFoam/buoyantSimpleRadiationFoam.C

@@ -44,7 +44,7 @@ int main(int argc, char *argv[])
     #include "setRootCase.H"
     #include "createTime.H"
     #include "createMesh.H"
-    #include "readEnvironmentalProperties.H"
+    #include "readGravitationalAcceleration.H"
     #include "createFields.H"
     #include "createRadiationModel.H"
     #include "initContinuityErrs.H"

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

@@ -4,6 +4,7 @@
     PtrList<volScalarField> KFluid(fluidRegions.size());
     PtrList<volVectorField> UFluid(fluidRegions.size());
     PtrList<surfaceScalarField> phiFluid(fluidRegions.size());
+    PtrList<uniformDimensionedVectorField> gFluid(fluidRegions.size());
     PtrList<compressible::turbulenceModel> turbulence(fluidRegions.size());
     PtrList<volScalarField> DpDtFluid(fluidRegions.size());
 
@@ -95,6 +96,23 @@
             )
         );
 
+        Info<< "    Adding to gFluid\n" << endl;
+        gFluid.set
+        (
+            i,
+            new uniformDimensionedVectorField
+            (
+                IOobject
+                (
+                    "g",
+                    runTime.constant(),
+                    fluidRegions[i],
+                    IOobject::MUST_READ,
+                    IOobject::NO_WRITE
+                )
+            )
+        );
+
         Info<< "    Adding to turbulence\n" << endl;
         turbulence.set
         (

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

@@ -19,27 +19,4 @@
                 )
             )
         );
-
-        // Force calculation of geometric properties to prevent it being done
-        // later in e.g. some boundary evaluation
-        //(void)fluidRegions[i].weights();
-        //(void)fluidRegions[i].deltaCoeffs();
-
-        // Attach environmental properties to each region
-        autoPtr<IOdictionary> environmentalProperties
-        (
-            new IOdictionary
-            (
-                IOobject
-                (
-                    "environmentalProperties",
-                    runTime.constant(),
-                    fluidRegions[i],
-                    IOobject::MUST_READ,
-                    IOobject::NO_WRITE
-                )
-            )
-        );
-
-        environmentalProperties.ptr()->store();
     }

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

@@ -5,6 +5,8 @@
     volScalarField& K = KFluid[i];
     volVectorField& U = UFluid[i];
     surfaceScalarField& phi = phiFluid[i];
+    const dimensionedVector& g = gFluid[i];
+
     compressible::turbulenceModel& turb = turbulence[i];
     volScalarField& DpDt = DpDtFluid[i];
 
@@ -13,9 +15,3 @@
     volScalarField& h = thermo.h();
 
     const dimensionedScalar massIni("massIni", dimMass, initialMassFluid[i]);
-
-    const dictionary& environmentalProperties =
-        fluidRegions[i].lookupObject<IOdictionary>
-        ("environmentalProperties");
-
-    const dimensionedVector g(environmentalProperties.lookup("g"));

applications/solvers/incompressible/pimpleDyMFoam/correctPhi.H

@@ -1,5 +1,9 @@
 {
-    wordList pcorrTypes(p.boundaryField().types());
+    wordList pcorrTypes
+    (
+        p.boundaryField().size(),
+        zeroGradientFvPatchScalarField::typeName
+    );
 
     for (label i=0; i<p.boundaryField().size(); i++)
     {

applications/solvers/incompressible/shallowWaterFoam/readGravitationalAcceleration.H

@@ -1,10 +1,10 @@
-Info<< "\nReading environmentalProperties" << endl;
+Info<< "\nReading gravitationalProperties" << endl;
 
-IOdictionary environmentalProperties
+IOdictionary gravitationalProperties
 (
     IOobject
     (
-        "environmentalProperties",
+        "gravitationalProperties",
         runTime.constant(),
         mesh,
         IOobject::MUST_READ,
@@ -12,10 +12,10 @@ IOdictionary environmentalProperties
     )
 );
 
-const dimensionedVector g(environmentalProperties.lookup("g"));
+const dimensionedVector g(gravitationalProperties.lookup("g"));
-const Switch rotating(environmentalProperties.lookup("rotating"));
+const Switch rotating(gravitationalProperties.lookup("rotating"));
 const dimensionedVector Omega =
-    rotating ? environmentalProperties.lookup("Omega")
+    rotating ? gravitationalProperties.lookup("Omega")
              : dimensionedVector("Omega", -dimTime, vector(0,0,0));
 const dimensionedScalar magg = mag(g);
 const dimensionedVector gHat = g/magg;

applications/solvers/incompressible/shallowWaterFoam/shallowWaterFoam.C

@@ -41,7 +41,7 @@ int main(int argc, char *argv[])
     #include "setRootCase.H"
     #include "createTime.H"
     #include "createMesh.H"
-    #include "readEnvironmentalProperties.H"
+    #include "readGravitationalAcceleration.H"
     #include "createFields.H"
 
     // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

applications/solvers/lagrangian/coalChemistryFoam/coalChemistryFoam.C

@@ -51,7 +51,7 @@ int main(int argc, char *argv[])
     #include "createTime.H"
     #include "createMesh.H"
     #include "readChemistryProperties.H"
-    #include "readEnvironmentalProperties.H"
+    #include "readGravitationalAcceleration.H"
     #include "createFields.H"
     #include "createClouds.H"
     #include "createRadiationModel.H"
@@ -114,7 +114,10 @@ int main(int argc, char *argv[])
 
         rho = thermo.rho();
 
-        runTime.write();
+        if (runTime.write())
+        {
+            chemistry.dQ()().write();
+        }
 
         Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
             << "  ClockTime = " << runTime.elapsedClockTime() << " s"

applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/UEqn.H

@@ -9,35 +9,10 @@
       + parcels.SU()
     );
 
-    UEqn.relax();
-
-    tmp<volScalarField> trAU;
-    tmp<volTensorField> trTU;
-
-    if (pressureImplicitPorosity)
-    {
-        tmp<volTensorField> tTU = tensor(I)*UEqn.A();
-        pZones.addResistance(UEqn, tTU());
-        trTU = inv(tTU());
-        trTU().rename("rAU");
-
-        volVectorField gradp = fvc::grad(p);
-
-        for (int UCorr=0; UCorr<nUCorr; UCorr++)
-        {
-            U = trTU() & (UEqn.H() - gradp);
-        }
-        U.correctBoundaryConditions();
-    }
-    else
-    {
     pZones.addResistance(UEqn);
 
-        trAU = 1.0/UEqn.A();
+    if (momentumPredictor)
-        trAU().rename("rAU");
+    {
-
+        solve(UEqn == -fvc::grad(p));
-        solve
-        (
-            UEqn == -fvc::grad(p)
-        );
     }
+

applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/YEqn.H

@@ -43,5 +43,4 @@ tmp<fv::convectionScheme<scalar> > mvConvection
 
     Y[inertIndex] = scalar(1) - Yt;
     Y[inertIndex].max(0.0);
-
 }

applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/createFields.H

@@ -73,13 +73,7 @@
         )
     );
 
-    Info<< "Creating field DpDt\n" << endl;
+    Info<< "Creating multi-variate interpolation scheme\n" << endl;
-    volScalarField DpDt
-    (
-        "DpDt",
-        fvc::DDt(surfaceScalarField("phiU", phi/fvc::interpolate(rho)), p)
-    );
-
     multivariateSurfaceInterpolationScheme<scalar>::fieldTable fields;
 
     forAll (Y, i)

applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/createPorousZones.H

@@ -1,21 +1,3 @@
     Info<< "Creating porous zones" << nl << endl;
 
     porousZones pZones(mesh);
-    Switch pressureImplicitPorosity(false);
-
-    label nUCorr = 0;
-    if (pZones.size())
-    {
-        // nUCorrectors for pressureImplicitPorosity
-        if (mesh.solutionDict().subDict("PISO").found("nUCorrectors"))
-        {
-            mesh.solutionDict().subDict("PISO").lookup("nUCorrectors")
-                >> nUCorr;
-        }
-
-        if (nUCorr > 0)
-        {
-            pressureImplicitPorosity = true;
-        }
-    }
-

applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/hEqn.H

@@ -1,20 +1,51 @@
 {
-    fvScalarMatrix hEqn
+    tmp<volScalarField> pWork
+    (
+        new volScalarField
+        (
+            IOobject
+            (
+                "pWork",
+                runTime.timeName(),
+                mesh,
+                IOobject::NO_READ,
+                IOobject::NO_WRITE
+            ),
+            mesh,
+            dimensionedScalar("zero", dimensionSet(1, -1, -3, 0, 0), 0.0)
+        )
+    );
+
+    if (dpdt)
+    {
+        pWork() += fvc::ddt(p);
+    }
+    if (eWork)
+    {
+        pWork() = -p*fvc::div(phi/fvc::interpolate(rho));
+    }
+    if (hWork)
+    {
+        pWork() += fvc::div(phi/fvc::interpolate(rho)*fvc::interpolate(p));
+    }
+
+    {
+        solve
         (
             fvm::ddt(rho, h)
           + mvConvection->fvmDiv(phi, h)
           - fvm::laplacian(turbulence->alphaEff(), h)
          ==
-        DpDt
+            pWork()
           + parcels.Sh()
           + radiation->Sh(thermo)
         );
 
-    hEqn.relax();
-
-    hEqn.solve();
-
         thermo.correct();
 
         radiation->correct();
+
+        Info<< "T gas min/max   = " << min(T).value() << ", "
+            << max(T).value() << endl;
+    }
 }

applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/pEqn.H

@@ -5,65 +5,13 @@
     // pressure solution - done in 2 parts. Part 1:
     thermo.rho() -= psi*p;
 
-    if (pressureImplicitPorosity)
+    volScalarField rAU = 1.0/UEqn.A();
-    {
+    U = rAU*UEqn.H();
-        U = trTU()&UEqn.H();
-    }
-    else
-    {
-        U = trAU()*UEqn.H();
-    }
 
-    if (transonic)
+    if (pZones.size() > 0)
     {
-        surfaceScalarField phiv = fvc::interpolate(U) & mesh.Sf();
+        // ddtPhiCorr not well defined for cases with porosity
-
+        phi = fvc::interpolate(rho)*(fvc::interpolate(U) & mesh.Sf());
-        phi = fvc::interpolate(rho)*phiv;
-
-        surfaceScalarField phid("phid", fvc::interpolate(psi)*phiv);
-
-        for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
-        {
-            tmp<fvScalarMatrix> lapTerm;
-
-            if (pressureImplicitPorosity)
-            {
-                lapTerm = fvm::laplacian(rho*trTU(), p);
-            }
-            else
-            {
-                lapTerm = fvm::laplacian(rho*trAU(), p);
-            }
-
-            fvScalarMatrix pEqn
-            (
-                fvc::ddt(rho) + fvc::div(phi)
-              + correction(psi*fvm::ddt(p) + fvm::div(phid, p))
-              - lapTerm()
-            ==
-                parcels.Srho()
-              + pointMassSources.Su()
-            );
-
-            if
-            (
-                oCorr == nOuterCorr-1
-             && corr == nCorr-1
-             && nonOrth == nNonOrthCorr
-            )
-            {
-                pEqn.solve(mesh.solver("pFinal"));
-            }
-            else
-            {
-                pEqn.solve();
-            }
-
-            if (nonOrth == nNonOrthCorr)
-            {
-                phi == pEqn.flux();
-            }
-        }
     }
     else
     {
@@ -71,38 +19,23 @@
             fvc::interpolate(rho)
            *(
                 (fvc::interpolate(U) & mesh.Sf())
-//            + fvc::ddtPhiCorr(trAU(), rho, U, phi)
+              + fvc::ddtPhiCorr(rAU, rho, U, phi)
             );
+    }
 
     for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
     {
-            tmp<fvScalarMatrix> lapTerm;
-
-            if (pressureImplicitPorosity)
-            {
-                lapTerm = fvm::laplacian(rho*trTU(), p);
-            }
-            else
-            {
-                lapTerm = fvm::laplacian(rho*trAU(), p);
-            }
-
         fvScalarMatrix pEqn
         (
             fvc::ddt(rho) + psi*correction(fvm::ddt(p))
           + fvc::div(phi)
-              - lapTerm()
+          - fvm::laplacian(rho*rAU, p)
         ==
             parcels.Srho()
           + pointMassSources.Su()
         );
 
-            if
+        if (corr == nCorr-1 && nonOrth == nNonOrthCorr)
-            (
-                oCorr == nOuterCorr-1
-             && corr == nCorr-1
-             && nonOrth == nNonOrthCorr
-            )
         {
             pEqn.solve(mesh.solver("pFinal"));
         }
@@ -116,7 +49,6 @@
             phi += pEqn.flux();
         }
     }
-    }
 
     // Second part of thermodynamic density update
     thermo.rho() += psi*p;
@@ -124,15 +56,6 @@
     #include "rhoEqn.H"
     #include "compressibleContinuityErrs.H"
 
-    if (pressureImplicitPorosity)
+    U -= rAU*fvc::grad(p);
-    {
-        U -= trTU()&fvc::grad(p);
-    }
-    else
-    {
-        U -= trAU()*fvc::grad(p);
-    }
     U.correctBoundaryConditions();
-
-    DpDt = fvc::DDt(surfaceScalarField("phiU", phi/fvc::interpolate(rho)), p);
 }

applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/porousExplicitSourceReactingParcelFoam.C

@@ -23,15 +23,16 @@ License
     Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
 
 Application
-    trackedReactingParcelFoam
+    porousExplicitSourceReactingParcelFoam
 
 Description
-    - reacting parcel cloud tracking
+    - reacting parcel cloud
     - porous media
     - point mass sources
     - polynomial based, incompressible thermodynamics (f(T))
 
-    Note: ddtPhiCorr not used here - not well defined for porous calcs
+    Note: ddtPhiCorr not used here when porous zones are active
+    - not well defined for porous calcs
 
 \*---------------------------------------------------------------------------*/
 
@@ -55,7 +56,7 @@ int main(int argc, char *argv[])
     #include "createTime.H"
     #include "createMesh.H"
     #include "readChemistryProperties.H"
-    #include "readEnvironmentalProperties.H"
+    #include "readGravitationalAcceleration.H"
     #include "createFields.H"
     #include "createRadiationModel.H"
     #include "createClouds.H"
@@ -74,6 +75,7 @@ int main(int argc, char *argv[])
     {
         #include "readTimeControls.H"
         #include "readPISOControls.H"
+        #include "readAdditionalSolutionControls.H"
         #include "compressibleCourantNo.H"
         #include "setDeltaT.H"
 
@@ -88,28 +90,23 @@ int main(int argc, char *argv[])
         #include "chemistry.H"
         #include "rhoEqn.H"
         #include "UEqn.H"
-
-        // --- PIMPLE loop
-        for (int oCorr=1; oCorr<=nOuterCorr; oCorr++)
-        {
         #include "YEqn.H"
         #include "hEqn.H"
 
         // --- PISO loop
-            for (int corr=1; corr<=nCorr; corr++)
+        for (int corr=0; corr<nCorr; corr++)
         {
             #include "pEqn.H"
         }
 
-            Info<< "T gas min/max   = " << min(T).value() << ", "
-                << max(T).value() << endl;
-        }
-
         turbulence->correct();
 
         rho = thermo.rho();
 
-        runTime.write();
+        if (runTime.write())
+        {
+            chemistry.dQ()().write();
+        }
 
         Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
             << "  ClockTime = " << runTime.elapsedClockTime() << " s"

applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/readAdditionalSolutionControls.H

@@ -0,0 +1,20 @@
+dictionary additional = mesh.solutionDict().subDict("additional");
+
+bool dpdt = true;
+if (additional.found("dpdt"))
+{
+    additional.lookup("dpdt") >> dpdt;
+}
+
+bool eWork = true;
+if (additional.found("eWork"))
+{
+    additional.lookup("eWork") >> eWork;
+}
+
+bool hWork = true;
+if (additional.found("hWork"))
+{
+    additional.lookup("hWork") >> hWork;
+}
+

applications/solvers/lagrangian/reactingParcelFoam/reactingParcelFoam.C

@@ -49,7 +49,7 @@ int main(int argc, char *argv[])
     #include "createTime.H"
     #include "createMesh.H"
     #include "readChemistryProperties.H"
-    #include "readEnvironmentalProperties.H"
+    #include "readGravitationalAcceleration.H"
     #include "createFields.H"
     #include "createClouds.H"
     #include "createRadiationModel.H"
@@ -101,7 +101,10 @@ int main(int argc, char *argv[])
 
         rho = thermo.rho();
 
-        runTime.write();
+        if (runTime.write())
+        {
+            chemistry.dQ()().write();
+        }
 
         Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
             << "  ClockTime = " << runTime.elapsedClockTime() << " s"

applications/solvers/lagrangian/uncoupledKinematicParcelFoam/uncoupledKinematicParcelFoam.C

@@ -45,7 +45,7 @@ int main(int argc, char *argv[])
     #include "setRootCase.H"
     #include "createTime.H"
     #include "createMesh.H"
-    #include "readEnvironmentalProperties.H"
+    #include "readGravitationalAcceleration.H"
     #include "createFields.H"
     #include "compressibleCourantNo.H"
 

applications/solvers/multiphase/bubbleFoam/alphaEqn.H

@@ -46,3 +46,5 @@
         << "  Max(alpha) = " << max(alpha).value()
         << endl;
 }
+
+rho = alpha*rhoa + beta*rhob;

applications/solvers/multiphase/bubbleFoam/bubbleFoam.C

@@ -44,7 +44,7 @@ int main(int argc, char *argv[])
 
     #include "createTime.H"
     #include "createMesh.H"
-    #include "readEnvironmentalProperties.H"
+    #include "readGravitationalAcceleration.H"
     #include "createFields.H"
     #include "initContinuityErrs.H"
 

applications/solvers/multiphase/bubbleFoam/createFields.H

@@ -156,6 +156,16 @@
       + fvc::interpolate(beta)*phib
     );
 
+    volScalarField rho
+    (
+        IOobject
+        (
+            "rho",
+            runTime.timeName(),
+            mesh
+        ),
+        alpha*rhoa + beta*rhob
+    );
 
     #include "createRASTurbulence.H"
 

applications/solvers/multiphase/compressibleInterDyMFoam/compressibleInterDyMFoam.C

@@ -50,7 +50,7 @@ int main(int argc, char *argv[])
     #include "setRootCase.H"
     #include "createTime.H"
     #include "createDynamicFvMesh.H"
-    #include "readEnvironmentalProperties.H"
+    #include "readGravitationalAcceleration.H"
     #include "readControls.H"
     #include "initContinuityErrs.H"
     #include "createFields.H"
@@ -104,12 +104,6 @@ int main(int argc, char *argv[])
             #include "meshCourantNo.H"
         }
 
-        if (mesh.changing())
-        {
-            gh = g & mesh.C();
-            ghf = g & mesh.Cf();
-        }
-
         turbulence->correct();
 
         // --- Outer-corrector loop

applications/solvers/multiphase/compressibleInterDyMFoam/createFields.H

@@ -46,11 +46,6 @@
     #include "createPhi.H"
 
 
-    Info<< "Calculating field g.h\n" << endl;
-    volScalarField gh("gh", g & mesh.C());
-    surfaceScalarField ghf("ghf", g & mesh.Cf());
-
-
     Info<< "Reading transportProperties\n" << endl;
     twoPhaseMixture twoPhaseProperties(U, phi);
 
@@ -134,7 +129,11 @@
     );
 
 
-    wordList pcorrTypes(p.boundaryField().types());
+    wordList pcorrTypes
+    (
+        p.boundaryField().size(),
+        zeroGradientFvPatchScalarField::typeName
+    );
 
     for (label i=0; i<p.boundaryField().size(); i++)
     {

applications/solvers/multiphase/compressibleInterFoam/compressibleInterFoam.C

@@ -49,7 +49,7 @@ int main(int argc, char *argv[])
     #include "setRootCase.H"
     #include "createTime.H"
     #include "createMesh.H"
-    #include "readEnvironmentalProperties.H"
+    #include "readGravitationalAcceleration.H"
     #include "readControls.H"
     #include "initContinuityErrs.H"
     #include "createFields.H"

applications/solvers/multiphase/interDyMFoam/createFields.H

@@ -92,7 +92,11 @@
         incompressible::turbulenceModel::New(U, phi, twoPhaseProperties)
     );
 
-    wordList pcorrTypes(p.boundaryField().types());
+    wordList pcorrTypes
+    (
+        p.boundaryField().size(),
+        zeroGradientFvPatchScalarField::typeName
+    );
 
     for (label i=0; i<p.boundaryField().size(); i++)
     {

applications/solvers/multiphase/interDyMFoam/interDyMFoam.C

@@ -48,7 +48,7 @@ int main(int argc, char *argv[])
     #include "setRootCase.H"
     #include "createTime.H"
     #include "createDynamicFvMesh.H"
-    #include "readEnvironmentalProperties.H"
+    #include "readGravitationalAcceleration.H"
     #include "readPISOControls.H"
     #include "initContinuityErrs.H"
     #include "createFields.H"
@@ -86,9 +86,6 @@ int main(int argc, char *argv[])
                 << " s" << endl;
         }
 
-        volScalarField gh("gh", g & mesh.C());
-        surfaceScalarField ghf("ghf", g & mesh.Cf());
-
         if (mesh.changing() && correctPhi)
         {
             #include "correctPhi.H"

applications/solvers/multiphase/interDyMFoam/pEqn.H

@@ -5,19 +5,19 @@
     U = rAU*UEqn.H();
     surfaceScalarField phiU("phiU", (fvc::interpolate(U) & mesh.Sf()));
 
+    if (p.needReference())
+    {
+        fvc::makeRelative(phiU, U);
+        adjustPhi(phiU, U, p);
+        fvc::makeAbsolute(phiU, U);
+    }
+
     phi = phiU +
     (
         fvc::interpolate(interface.sigmaK())*fvc::snGrad(alpha1)*mesh.magSf()
       + fvc::interpolate(rho)*(g & mesh.Sf())
     )*rAUf;
 
-    if (p.needReference())
-    {
-        fvc::makeRelative(phi, U);
-        adjustPhi(phi, U, p);
-        fvc::makeAbsolute(phi, U);
-    }
-
     for(int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
     {
         fvScalarMatrix pEqn

applications/solvers/multiphase/interFoam/correctPhi.H

@@ -1,7 +1,11 @@
 {
     #include "continuityErrs.H"
 
-    wordList pcorrTypes(p.boundaryField().types());
+    wordList pcorrTypes
+    (
+        p.boundaryField().size(),
+        zeroGradientFvPatchScalarField::typeName
+    );
 
     for (label i=0; i<p.boundaryField().size(); i++)
     {

applications/solvers/multiphase/interFoam/interFoam.C

@@ -51,7 +51,7 @@ int main(int argc, char *argv[])
     #include "setRootCase.H"
     #include "createTime.H"
     #include "createMesh.H"
-    #include "readEnvironmentalProperties.H"
+    #include "readGravitationalAcceleration.H"
     #include "readPISOControls.H"
     #include "initContinuityErrs.H"
     #include "createFields.H"

applications/solvers/multiphase/interFoam/pEqn.H

@@ -7,17 +7,18 @@
     surfaceScalarField phiU
     (
         "phiU",
-        (fvc::interpolate(U) & mesh.Sf()) + fvc::ddtPhiCorr(rUA, rho, U, phi)
+        (fvc::interpolate(U) & mesh.Sf())
+      + fvc::ddtPhiCorr(rUA, rho, U, phi)
     );
 
+    adjustPhi(phiU, U, p);
+
     phi = phiU +
         (
             fvc::interpolate(interface.sigmaK())*fvc::snGrad(alpha1)*mesh.magSf()
           + fvc::interpolate(rho)*(g & mesh.Sf())
         )*rUAf;
 
-    adjustPhi(phi, U, p);
-
     for(int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
     {
         fvScalarMatrix pEqn

applications/solvers/multiphase/interPhaseChangeFoam/correctPhi.H

@@ -1,7 +1,11 @@
 {
     #include "continuityErrs.H"
 
-    wordList pcorrTypes(p.boundaryField().types());
+    wordList pcorrTypes
+    (
+        p.boundaryField().size(),
+        zeroGradientFvPatchScalarField::typeName
+    );
 
     for (label i=0; i<p.boundaryField().size(); i++)
     {

applications/solvers/multiphase/interPhaseChangeFoam/interPhaseChangeFoam.C

@@ -53,7 +53,7 @@ int main(int argc, char *argv[])
     #include "setRootCase.H"
     #include "createTime.H"
     #include "createMesh.H"
-    #include "readEnvironmentalProperties.H"
+    #include "readGravitationalAcceleration.H"
     #include "readPISOControls.H"
     #include "initContinuityErrs.H"
     #include "createFields.H"

applications/solvers/multiphase/interPhaseChangeFoam/pEqn.H

@@ -11,14 +11,14 @@
       + fvc::ddtPhiCorr(rUA, rho, U, phi)
     );
 
+    adjustPhi(phiU, U, p);
+
     phi = phiU +
         (
             fvc::interpolate(interface.sigmaK())*fvc::snGrad(alpha1)*mesh.magSf()
           + fvc::interpolate(rho)*(g & mesh.Sf())
         )*rUAf;
 
-    adjustPhi(phi, U, p);
-
     Pair<tmp<volScalarField> > vDotP = twoPhaseProperties->vDotP();
     const volScalarField& vDotcP = vDotP[0]();
     const volScalarField& vDotvP = vDotP[1]();

applications/solvers/multiphase/multiphaseInterFoam/multiphaseInterFoam.C

@@ -44,7 +44,7 @@ int main(int argc, char *argv[])
     #include "setRootCase.H"
     #include "createTime.H"
     #include "createMesh.H"
-    #include "readEnvironmentalProperties.H"
+    #include "readGravitationalAcceleration.H"
     #include "readPISOControls.H"
     #include "initContinuityErrs.H"
     #include "createFields.H"

applications/solvers/multiphase/multiphaseInterFoam/pEqn.H

@@ -10,14 +10,14 @@
         (fvc::interpolate(U) & mesh.Sf()) + fvc::ddtPhiCorr(rUA, rho, U, phi)
     );
 
+    adjustPhi(phiU, U, p);
+
     phi = phiU +
         (
             mixture.surfaceTensionForce()*mesh.magSf()
           + fvc::interpolate(rho)*(g & mesh.Sf())
         )*rUAf;
 
-    adjustPhi(phi, U, p);
-
     for(int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
     {
         fvScalarMatrix pEqn

applications/solvers/multiphase/settlingFoam/settlingFoam.C

@@ -47,7 +47,7 @@ int main(int argc, char *argv[])
 
     #include "createTime.H"
     #include "createMesh.H"
-    #include "readEnvironmentalProperties.H"
+    #include "readGravitationalAcceleration.H"
     #include "createFields.H"
     #include "initContinuityErrs.H"
 

applications/solvers/multiphase/twoLiquidMixingFoam/twoLiquidMixingFoam.C

@@ -43,7 +43,7 @@ int main(int argc, char *argv[])
     #include "setRootCase.H"
     #include "createTime.H"
     #include "createMesh.H"
-    #include "readEnvironmentalProperties.H"
+    #include "readGravitationalAcceleration.H"
     #include "readPIMPLEControls.H"
     #include "initContinuityErrs.H"
     #include "createFields.H"

applications/solvers/multiphase/twoPhaseEulerFoam/alphaEqn.H

@@ -52,3 +52,5 @@
             << endl;
     }
 }
+
+rho = alpha*rhoa + beta*rhob;

applications/solvers/multiphase/twoPhaseEulerFoam/createFields.H

@@ -145,6 +145,16 @@
         fvc::interpolate(alpha)*phia + fvc::interpolate(beta)*phib
     );
 
+    volScalarField rho
+    (
+        IOobject
+        (
+            "rho",
+            runTime.timeName(),
+            mesh
+        ),
+        alpha*rhoa + beta*rhob
+    );
 
     IOdictionary RASProperties
     (

applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseEulerFoam.C

@@ -51,7 +51,7 @@ int main(int argc, char *argv[])
 
     #include "createTime.H"
     #include "createMesh.H"
-    #include "readEnvironmentalProperties.H"
+    #include "readGravitationalAcceleration.H"
     #include "createFields.H"
     #include "readPPProperties.H"
     #include "initContinuityErrs.H"

applications/utilities/mesh/generation/blockMesh/blockMeshApp.C

@@ -348,7 +348,8 @@ int main(int argc, char *argv[])
             << exit(FatalError);
     }
 
-    Info<< nl << "end" << endl;
+    Info<< nl << "End" << endl;
+
     return 0;
 }
 

applications/utilities/mesh/generation/extrudeMesh/extrudeMesh.C

@@ -23,8 +23,9 @@ License
     Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
 
 Description
-    Extrude mesh from existing patch (flipped so has inwards pointing
+    Extrude mesh from existing patch (by default outwards facing normals;
-    normals) or from patch read from file.
+    optional flips faces)
+    or from patch read from file.
     Note: Merges close points so be careful.
 
     Type of extrusion prescribed by run-time selectable model.
@@ -78,6 +79,7 @@ int main(int argc, char *argv[])
     if (sourceType == "patch")
     {
         fileName sourceCasePath(dict.lookup("sourceCase"));
+        sourceCasePath.expand();
         fileName sourceRootDir = sourceCasePath.path();
         fileName sourceCaseDir = sourceCasePath.name();
         word patchName(dict.lookup("sourcePatch"));
@@ -107,14 +109,11 @@ int main(int argc, char *argv[])
 
         const polyPatch& pp = mesh.boundaryMesh()[patchID];
         fMesh.reset(new faceMesh(pp.localFaces(), pp.localPoints()));
-        fMesh().flip();
 
         {
-            fileName surfName(patchName + ".sMesh");
+            fileName surfName(runTime.path()/patchName + ".sMesh");
-
+            Info<< "Writing patch as surfaceMesh to "
-            Info<< "Writing (flipped) patch as surfaceMesh to "
                 << surfName << nl << endl;
-
             OFstream os(surfName);
             os << fMesh() << nl;
         }
@@ -141,6 +140,20 @@ int main(int argc, char *argv[])
             << "patch or surface." << exit(FatalError);
     }
 
+    Switch flipNormals(dict.lookup("flipNormals"));
+
+    if (flipNormals)
+    {
+        Info<< "Flipping faces." << nl << endl;
+
+        faceList faces(fMesh().size());
+        forAll(faces, i)
+        {
+            faces[i] = fMesh()[i].reverseFace();
+        }
+        fMesh.reset(new faceMesh(faces, fMesh().localPoints()));
+    }
+
 
     Info<< "Extruding patch with :" << nl
             << "    points     : " << fMesh().points().size() << nl

applications/utilities/mesh/generation/extrudeMesh/extrudeModel/linearNormal/linearNormal.C

@@ -47,7 +47,14 @@ linearNormal::linearNormal(const dictionary& dict)
 :
     extrudeModel(typeName, dict),
     thickness_(readScalar(coeffDict_.lookup("thickness")))
-{}
+{
+    if (thickness_ <= 0)
+    {
+        FatalErrorIn("linearNormal(const dictionary&)")
+            << "thickness should be positive : " << thickness_
+            << exit(FatalError);
+    }
+}
 
 
 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //

applications/utilities/mesh/generation/extrudeMesh/extrudeProperties

@@ -19,31 +19,40 @@ FoamFile
 constructFrom patch;    //surface;
 
 // If construct from (flipped) patch
-sourceCase "../cavity";
+sourceCase "$FOAM_RUN/icoFoam/cavity";
 sourcePatch movingWall;
 
+// Flip surface normals before usage.
+flipNormals false;
+
 // If construct from surface
 surface "movingWall.sMesh";
 
 
-// Do front and back need to be merged?
+// Do front and back need to be merged? Usually only makes sense for 360
-mergeFaces false;
+// degree wedges.
+mergeFaces true;
+
 
 //- Linear extrusion in point-normal direction
 //extrudeModel        linearNormal;
+
 //- Wedge extrusion. If nLayers is 1 assumes symmetry around plane.
 extrudeModel        wedge;
+
 //- Extrudes into sphere around (0 0 0)
 //extrudeModel        linearRadial;
+
+//- Extrudes into sphere with grading according to pressure (atmospherics)
 //extrudeModel        sigmaRadial;
 
-nLayers             6;
+nLayers             20;
 
 wedgeCoeffs
 {
     axisPt      (0 0.1 0);
-    axis        (1 0 0);
+    axis        (-1 0 0);
-    angle       90.0;   // For nLayers=1 assume symmetry so angle/2 on each side
+    angle       360;  // For nLayers=1 assume symmetry so angle/2 on each side
 }
 
 linearNormalCoeffs

applications/utilities/mesh/generation/extrudeMesh/extrudedMesh/extrudedMesh.C

@@ -109,15 +109,33 @@ Foam::Xfer<Foam::faceList> Foam::extrudedMesh::extrudedFaces
         label currentLayerOffset = layer*surfacePoints.size();
         label nextLayerOffset = currentLayerOffset + surfacePoints.size();
 
-        // Side faces from layer to layer+1
+        // Vertical faces from layer to layer+1
-        for (label i=0; i<nInternalEdges; i++)
+        for (label edgeI=0; edgeI<nInternalEdges; edgeI++)
         {
-            quad[0] = surfaceEdges[i][1] + currentLayerOffset;
+            const edge& e = surfaceEdges[edgeI];
-            quad[1] = surfaceEdges[i][0] + currentLayerOffset;
+            const labelList& edgeFaces = extrudePatch.edgeFaces()[edgeI];
-            quad[2] = surfaceEdges[i][0] + nextLayerOffset;
-            quad[3] = surfaceEdges[i][1] + nextLayerOffset;
 
-            eFaces[facei++] = face(quad).reverseFace();
+            face& f = eFaces[facei++];
+            f.setSize(4);
+
+            if
+            (
+                (edgeFaces[0] < edgeFaces[1])
+             == sameOrder(surfaceFaces[edgeFaces[0]], e)
+            )
+            {
+                f[0] = e[0] + currentLayerOffset;
+                f[1] = e[1] + currentLayerOffset;
+                f[2] = e[1] + nextLayerOffset;
+                f[3] = e[0] + nextLayerOffset;
+            }
+            else
+            {
+                f[0] = e[1] + currentLayerOffset;
+                f[1] = e[0] + currentLayerOffset;
+                f[2] = e[0] + nextLayerOffset;
+                f[3] = e[1] + nextLayerOffset;
+            }
         }
 
         // Faces between layer and layer+1
@@ -128,9 +146,9 @@ Foam::Xfer<Foam::faceList> Foam::extrudedMesh::extrudedFaces
                 eFaces[facei++] =
                     face
                     (
-                        surfaceFaces[i].reverseFace()
+                        surfaceFaces[i] //.reverseFace()
                       + nextLayerOffset
-                    ).reverseFace();
+                    );
             }
         }
     }
@@ -142,40 +160,46 @@ Foam::Xfer<Foam::faceList> Foam::extrudedMesh::extrudedFaces
         label nextLayerOffset = currentLayerOffset + surfacePoints.size();
 
         // Side faces across layer
-        for (label i=nInternalEdges; i<surfaceEdges.size(); i++)
+        for (label edgeI=nInternalEdges; edgeI<surfaceEdges.size(); edgeI++)
         {
-            const edge& e = surfaceEdges[i];
+            const edge& e = surfaceEdges[edgeI];
-            quad[0] = e[1] + currentLayerOffset;
+            const labelList& edgeFaces = extrudePatch.edgeFaces()[edgeI];
-            quad[1] = e[0] + currentLayerOffset;
-            quad[2] = e[0] + nextLayerOffset;
-            quad[3] = e[1] + nextLayerOffset;
 
-            label ownerFace = extrudePatch.edgeFaces()[i][0];
+            face& f = eFaces[facei++];
+            f.setSize(4);
 
-            if (sameOrder(surfaceFaces[ownerFace], e))
+            if (sameOrder(surfaceFaces[edgeFaces[0]], e))
             {
-                reverse(quad);
+                f[0] = e[0] + currentLayerOffset;
+                f[1] = e[1] + currentLayerOffset;
+                f[2] = e[1] + nextLayerOffset;
+                f[3] = e[0] + nextLayerOffset;
             }
-
+            else
-            eFaces[facei++] = face(quad);
+            {
+                f[0] = e[1] + currentLayerOffset;
+                f[1] = e[0] + currentLayerOffset;
+                f[2] = e[0] + nextLayerOffset;
+                f[3] = e[1] + nextLayerOffset;
             }
         }
-
-    // Top faces
-    forAll(surfaceFaces, i)
-    {
-        eFaces[facei++] = face(surfaceFaces[i]).reverseFace();
     }
 
     // Bottom faces
     forAll(surfaceFaces, i)
+    {
+        eFaces[facei++] = face(surfaceFaces[i]).reverseFace();
+    }
+
+    // Top faces
+    forAll(surfaceFaces, i)
     {
         eFaces[facei++] =
             face
             (
-                surfaceFaces[i].reverseFace()
+                surfaceFaces[i]
               + nLayers*surfacePoints.size()
-            ).reverseFace();
+            );
     }
 
     // return points for transferring

applications/utilities/mesh/generation/extrudeMesh/setRoots.H

@@ -1,16 +0,0 @@
-    argList::validArgs.clear();
-    argList::noParallel();
-
-    argList::validOptions.insert("sourceCase", "source case");
-    argList::validOptions.insert("sourcePatch", "source patch");
-
-    argList::validOptions.insert("surface", "surface file");
-
-    argList::validOptions.insert("mergeFaces", "");
-
-    argList args(argc, argv);
-
-    if (!args.check())
-    {
-         FatalError.exit();
-    }

applications/utilities/mesh/generation/snappyHexMesh/snappyHexMeshDict

@@ -209,6 +209,8 @@ snapControls
 // Settings for the layer addition.
 addLayersControls
 {
+    // Are the thickness parameters below relative to the undistorted
+    // size of the refined cell outside layer (true) or absolute sizes (false).
     relativeSizes true;
 
     // Per final patch (so not geometry!) the layer information
@@ -230,14 +232,13 @@ addLayersControls
 
 
     //- Wanted thickness of final added cell layer. If multiple layers
-    //  is the
+    //  is the thickness of the layer furthest away from the wall.
-    //  thickness of the layer furthest away from the wall.
+    //  See relativeSizes parameter.
-    //  Relative to undistorted size of cell outside layer.
     finalLayerThickness 0.3;
 
     //- Minimum thickness of cell layer. If for any reason layer
     //  cannot be above minThickness do not add layer.
-    //  Relative to undistorted size of cell outside layer.
+    //  See relativeSizes parameter.
     minThickness 0.25;
 
     //- If points get not extruded do nGrow layers of connected faces that are

applications/utilities/mesh/manipulation/checkMesh/checkTopology.C

@@ -67,8 +67,13 @@ Foam::label Foam::checkTopology
     {
         forAll(mesh.faceZones(), zoneI)
         {
-            if (mesh.faceZones()[zoneI].checkParallelSync(true))
+            if (mesh.faceZones()[zoneI].checkParallelSync(false))
             {
+                Info<< " ***FaceZone " << mesh.faceZones()[zoneI].name()
+                    << " is not correctly synchronised"
+                    << " acrosss coupled boundaries."
+                    << " (coupled faces both"
+                    << " present in set but with opposite flipmap)" << endl;
                 noFailedChecks++;
             }
         }

applications/utilities/postProcessing/dataConversion/foamToEnsight/foamToEnsight.C

@@ -492,6 +492,8 @@ int main(int argc, char *argv[])
         delete ensightCaseFilePtr;
     }
 
+    Info<< "End\n" << endl;
+
     return 0;
 }
 

applications/utilities/postProcessing/dataConversion/foamToGMV/foamToGMV.C

@@ -88,7 +88,7 @@ int main(int argc, char *argv[])
         }
     }
 
-    Info << nl << "Translation Complete." << nl;
+    Info<< "\nEnd\n" << endl;
 
     return 0;
 }

applications/utilities/postProcessing/foamCalc/foamCalc.C

@@ -84,6 +84,8 @@ int main(int argc, char *argv[])
 
     utility().tryPostCalc(args, runTime, mesh);
 
+    Info<< "End\n" << endl;
+
     return 0;
 }
 

applications/utilities/postProcessing/graphics/PV3FoamReader/PV3FoamReader/CMakeLists.txt

@@ -23,7 +23,7 @@ INCLUDE_DIRECTORIES(
 )
 
 ADD_DEFINITIONS(
-    -D$ENV{WM_PRECISION_OPTION}
+    -DWM_$ENV{WM_PRECISION_OPTION}
 )
 
 # Set output library destination to plugin folder

applications/utilities/postProcessing/graphics/PVFoamReader/PVFoamReader/CMakeLists.txt

@@ -22,7 +22,7 @@ INCLUDE_DIRECTORIES(
     $ENV{ParaView_INST_DIR}/include
     ../vtkFoam/lnInclude
 )
-ADD_DEFINITIONS(-D$ENV{WM_PRECISION_OPTION})
+ADD_DEFINITIONS(-DWM_$ENV{WM_PRECISION_OPTION})
 
 #-----------------------------------------------------------------------------
 # Most users should not need to change anything below this line.

applications/utilities/postProcessing/lagrangian/particleTracks/particleTracks.C

@@ -38,7 +38,7 @@ Description
 #include "Time.H"
 #include "timeSelector.H"
 #include "OFstream.H"
-#include "passiveParticle.H"
+#include "passiveParticleCloud.H"
 
 using namespace Foam;
 
@@ -61,7 +61,8 @@ int main(int argc, char *argv[])
     fileName vtkPath(runTime.path()/"VTK");
     mkDir(vtkPath);
 
-    Info<< "Scanning times to determine track data" << nl << endl;
+    Info<< "Scanning times to determine track data for cloud " << cloudName
+        << nl << endl;
 
     labelList maxIds(Pstream::nProcs(), -1);
     forAll(timeDirs, timeI)
@@ -69,23 +70,12 @@ int main(int argc, char *argv[])
         runTime.setTime(timeDirs[timeI], timeI);
         Info<< "Time = " << runTime.timeName() << endl;
 
-        IOobject positionsHeader
-        (
-            "positions",
-            runTime.timeName(),
-            cloud::prefix/cloudName,
-            mesh,
-            IOobject::MUST_READ,
-            IOobject::NO_WRITE,
-            false
-        );
-
-        if (positionsHeader.headerOk())
-        {
         Info<< "    Reading particle positions" << endl;
-            Cloud<passiveParticle> myCloud(mesh, cloudName, false);
+        passiveParticleCloud myCloud(mesh, cloudName);
+        Info<< "    Read " << returnReduce(myCloud.size(), sumOp<label>())
+            << " particles" << endl;
 
-            forAllConstIter(Cloud<passiveParticle>, myCloud, iter)
+        forAllConstIter(passiveParticleCloud, myCloud, iter)
         {
             label origId = iter().origId();
             label origProc = iter().origProc();
@@ -93,11 +83,20 @@ int main(int argc, char *argv[])
             maxIds[origProc] = max(maxIds[origProc], origId);
         }
     }
-    }
+    Pstream::listCombineGather(maxIds, maxEqOp<label>());
-    Pstream::listCombineGather(maxIds, maxOp<label>());
     Pstream::listCombineScatter(maxIds);
+
     labelList numIds = maxIds + 1;
 
+    Info<< nl << "Particle statistics:" << endl;
+    forAll(maxIds, procI)
+    {
+        Info<< "    Found " << numIds[procI] << " particles originating"
+            << " from processor " << procI << endl;
+    }
+    Info<< nl << endl;
+
+
     // calc starting ids for particles on each processor
     List<label> startIds(numIds.size(), 0);
     for (label i = 0; i < numIds.size()-1; i++)
@@ -106,56 +105,49 @@ int main(int argc, char *argv[])
     }
     label nParticle = startIds[startIds.size()-1] + numIds[startIds.size()-1];
 
+
+
     // number of tracks to generate
     label nTracks = nParticle/sampleFrequency;
 
     // storage for all particle tracks
     List<DynamicList<vector> > allTracks(nTracks);
 
-    Info<< "\nGenerating " << nTracks << " particle tracks" << nl << endl;
+    Info<< "\nGenerating " << nTracks << " particle tracks for cloud "
+        << cloudName << nl << endl;
+
     forAll(timeDirs, timeI)
     {
         runTime.setTime(timeDirs[timeI], timeI);
         Info<< "Time = " << runTime.timeName() << endl;
 
-        IOobject positionsHeader
-        (
-            "positions",
-            runTime.timeName(),
-            cloud::prefix/cloudName,
-            mesh,
-            IOobject::MUST_READ,
-            IOobject::NO_WRITE,
-            false
-        );
-
-        if (positionsHeader.headerOk())
-        {
-            Info<< "    Reading particle positions" << endl;
-            Cloud<passiveParticle> myCloud(mesh, cloudName, false);
-
-            // collect the track data on the master processor
         List<pointField> allPositions(Pstream::nProcs());
+        List<labelField> allOrigIds(Pstream::nProcs());
+        List<labelField> allOrigProcs(Pstream::nProcs());
+
+        // Read particles. Will be size 0 if no particles.
+        Info<< "    Reading particle positions" << endl;
+        passiveParticleCloud myCloud(mesh, cloudName);
+
+        // collect the track data on all processors that have positions
         allPositions[Pstream::myProcNo()].setSize
         (
             myCloud.size(),
             point::zero
         );
-
-            List<labelField> allOrigIds(Pstream::nProcs());
         allOrigIds[Pstream::myProcNo()].setSize(myCloud.size(), 0);
-
-            List<labelField> allOrigProcs(Pstream::nProcs());
         allOrigProcs[Pstream::myProcNo()].setSize(myCloud.size(), 0);
 
         label i = 0;
-            forAllConstIter(Cloud<passiveParticle>, myCloud, iter)
+        forAllConstIter(passiveParticleCloud, myCloud, iter)
         {
             allPositions[Pstream::myProcNo()][i] = iter().position();
             allOrigIds[Pstream::myProcNo()][i] = iter().origId();
             allOrigProcs[Pstream::myProcNo()][i] = iter().origProc();
             i++;
         }
+
+        // collect the track data on the master processor
         Pstream::gatherList(allPositions);
         Pstream::gatherList(allOrigIds);
         Pstream::gatherList(allOrigProcs);
@@ -186,18 +178,14 @@ int main(int argc, char *argv[])
             }
         }
     }
-        else
-        {
-            Info<< "    No particles read" << nl << endl;
-        }
-    }
 
     if (Pstream::master())
     {
-        Info<< "\nWriting particle tracks" << nl << endl;
-
         OFstream vtkTracks(vtkPath/"particleTracks.vtk");
 
+        Info<< "\nWriting particle tracks to " << vtkTracks.name()
+            << nl << endl;
+
         // Total number of points in tracks + 1 per track
         label nPoints = 0;
         forAll(allTracks, trackI)

applications/utilities/postProcessing/velocityField/Co/Co.C

@@ -179,6 +179,8 @@ void Foam::calc(const argList& args, const Time& runTime, const fvMesh& mesh)
     {
         Info<< "    No phi" << endl;
     }
+
+    Info<< "\nEnd\n" << endl;
 }
 
 // ************************************************************************* //

applications/utilities/postProcessing/velocityField/Lambda2/Lambda2.C

@@ -78,6 +78,8 @@ void Foam::calc(const argList& args, const Time& runTime, const fvMesh& mesh)
     {
         Info<< "    No U" << endl;
     }
+
+    Info<< "\nEnd\n" << endl;
 }
 
 

applications/utilities/postProcessing/velocityField/Mach/Mach.C

@@ -134,6 +134,8 @@ void Foam::calc(const argList& args, const Time& runTime, const fvMesh& mesh)
     {
         Info<< "    Missing U or T" << endl;
     }
+
+    Info<< "\nEnd\n" << endl;
 }
 
 

applications/utilities/postProcessing/velocityField/Pe/Pe.C

@@ -361,6 +361,9 @@ void Foam::calc(const argList& args, const Time& runTime, const fvMesh& mesh)
     {
         Info<< "    No phi" << endl;
     }
+
+    Info<< "\nEnd\n" << endl;
 }
 
+
 // ************************************************************************* //

applications/utilities/postProcessing/velocityField/Q/Q.C

@@ -105,6 +105,9 @@ void Foam::calc(const argList& args, const Time& runTime, const fvMesh& mesh)
     {
         Info<< "    No U" << endl;
     }
+
+    Info<< "\nEnd\n" << endl;
 }
 
+
 // ************************************************************************* //

applications/utilities/postProcessing/velocityField/enstrophy/enstrophy.C

@@ -80,6 +80,8 @@ void Foam::calc(const argList& args, const Time& runTime, const fvMesh& mesh)
     {
         Info<< "    No U" << endl;
     }
+
+    Info<< "\nEnd\n" << endl;
 }
 
 

applications/utilities/postProcessing/velocityField/flowType/flowType.C

@@ -87,6 +87,8 @@ void Foam::calc(const argList& args, const Time& runTime, const fvMesh& mesh)
     {
         Info<< "    No U" << endl;
     }
+
+    Info<< "\nEnd\n" << endl;
 }
 
 

applications/utilities/postProcessing/velocityField/streamFunction/streamFunction.C

@@ -460,7 +460,7 @@ int main(int argc, char *argv[])
         }
     }
 
-    Info<< "End\n" << endl;
+    Info<< "\nEnd\n" << endl;
 
     return 0;
 }

applications/utilities/postProcessing/velocityField/uprime/uprime.C

@@ -80,6 +80,9 @@ void Foam::calc(const argList& args, const Time& runTime, const fvMesh& mesh)
     {
         Info<< "    No k" << endl;
     }
+
+    Info<< "\nEnd\n" << endl;
 }
 
+
 // ************************************************************************* //

applications/utilities/postProcessing/velocityField/vorticity/vorticity.C

@@ -93,6 +93,9 @@ void Foam::calc(const argList& args, const Time& runTime, const fvMesh& mesh)
     {
         Info<< "    No U" << endl;
     }
+
+    Info<< "\nEnd\n" << endl;
 }
 
+
 // ************************************************************************* //