Merge branch 'master' into particleInteractions

Conflicts:
	src/lagrangian/basic/Cloud/Cloud.H
	src/lagrangian/intermediate/clouds/Templates/ThermoCloud/ThermoCloud.C
	src/parallel/decompositionMethods/Make/options

applications/solvers/basic/potentialFoam/potentialFoam.C

@@ -38,11 +38,9 @@ Description
 
 int main(int argc, char *argv[])
 {
-
-    argList::validOptions.insert("writep", "");
+    argList::addBoolOption("writep", "write the final pressure field");
 
 #   include "setRootCase.H"
-
 #   include "createTime.H"
 #   include "createMesh.H"
 #   include "createFields.H"

applications/solvers/combustion/PDRFoam/Make/options

@@ -15,8 +15,7 @@ EXE_INC = \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
     -I$(LIB_SRC)/dynamicFvMesh/lnInclude \
     -I$(LIB_SRC)/dynamicMesh/lnInclude \
-    -I$(LIB_SRC)/meshTools/lnInclude \
-    -I$(LIB_SRC)/triSurface/lnInclude
+    -I$(LIB_SRC)/meshTools/lnInclude
 
 EXE_LIBS = \
     -lengine \

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

@@ -124,11 +124,6 @@ public:
 
     // Member Functions
 
-        tmp<volScalarField> mut() const
-        {
-            return mut_;
-        }
-
         //- Return the effective diffusivity for k
         tmp<volScalarField> DkEff() const
         {
@@ -147,41 +142,44 @@ public:
             );
         }
 
-        //- Return the effective turbulent thermal diffusivity
-        tmp<volScalarField> alphaEff() const
+        //- Return the turbulence viscosity
+        virtual tmp<volScalarField> mut() const
         {
-            return tmp<volScalarField>
-            (
-                new volScalarField("alphaEff", alphat_ + alpha())
-            );
+            return mut_;
+        }
+
+        //- Return the turbulence thermal diffusivity
+        virtual tmp<volScalarField> alphat() const
+        {
+            return alphat_;
         }
 
         //- Return the turbulence kinetic energy
-        tmp<volScalarField> k() const
+        virtual tmp<volScalarField> k() const
         {
             return k_;
         }
 
         //- Return the turbulence kinetic energy dissipation rate
-        tmp<volScalarField> epsilon() const
+        virtual tmp<volScalarField> epsilon() const
         {
             return epsilon_;
         }
 
         //- Return the Reynolds stress tensor
-        tmp<volSymmTensorField> R() const;
+        virtual tmp<volSymmTensorField> R() const;
 
         //- Return the effective stress tensor including the laminar stress
-        tmp<volSymmTensorField> devRhoReff() const;
+        virtual tmp<volSymmTensorField> devRhoReff() const;
 
         //- Return the source term for the momentum equation
-        tmp<fvVectorMatrix> divDevRhoReff(volVectorField& U) const;
+        virtual tmp<fvVectorMatrix> divDevRhoReff(volVectorField& U) const;
 
         //- Solve the turbulence equations and correct the turbulence viscosity
-        void correct();
+        virtual void correct();
 
         //- Read turbulenceProperties dictionary
-        bool read();
+        virtual bool read();
 };
 
 

applications/solvers/combustion/PDRFoam/XiEqns

@@ -22,7 +22,7 @@
         volScalarField GEta = GEtaCoef/tauEta;
         volScalarField XiEqEta = 1.0 + XiCoef*sqrt(up/(Su + SuMin))*Reta;
 
-        volScalarField R = 
+        volScalarField R =
             GEta*XiEqEta/(XiEqEta - 0.999) + GIn*XiIn/(XiIn - 0.999);
 
         volScalarField XiEqStar = R/(R - GEta - GIn);
@@ -42,7 +42,7 @@
         volScalarField GEta = GEtaCoef/tauEta;
         volScalarField XiEqEta = 1.0 + XiCoef*sqrt(up/(Su + SuMin))*Reta;
 
-        volScalarField R = 
+        volScalarField R =
             GEta*XiEqEta/(XiEqEta - 0.999) + GIn*XiIn/(XiIn - 0.999);
 
         volScalarField XiEqStar = R/(R - GEta - GIn);

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

@@ -91,7 +91,7 @@ bool Foam::XiModels::algebraic::read(const dictionary& XiProperties)
     XiModel::read(XiProperties);
 
     XiModelCoeffs_.lookup("XiShapeCoef") >> XiShapeCoef;
-    
+
     return true;
 }
 

applications/solvers/combustion/PDRFoam/laminarFlameSpeed/SCOPE/SCOPELaminarFlameSpeed.C

@@ -117,7 +117,7 @@ inline Foam::scalar Foam::laminarFlameSpeedModels::SCOPE::polyPhi
 {
     scalar x = phi - 1.0;
 
-    return 
+    return
         a[0]
        *(
            scalar(1)

applications/solvers/combustion/PDRFoam/setDeltaT.H

@@ -45,7 +45,7 @@ if (adjustTimeStep)
             maxDeltaT
         )
     );
-    
+
     Info<< "deltaT = " <<  runTime.deltaTValue() << endl;
 }
 

applications/solvers/combustion/dieselEngineFoam/createFields.H

@@ -6,7 +6,7 @@ autoPtr<psiChemistryModel> pChemistry
 );
 psiChemistryModel& chemistry = pChemistry();
 
-hCombustionThermo& thermo = chemistry.thermo();
+hsCombustionThermo& thermo = chemistry.thermo();
 
 basicMultiComponentMixture& composition = thermo.composition();
 PtrList<volScalarField>& Y = composition.Y();
@@ -50,7 +50,7 @@ volVectorField U
 volScalarField& p = thermo.p();
 const volScalarField& psi = thermo.psi();
 const volScalarField& T = thermo.T();
-volScalarField& h = thermo.h();
+volScalarField& hs = thermo.hs();
 
 
 #include "compressibleCreatePhi.H"
@@ -88,8 +88,22 @@ volScalarField DpDt =
 
 multivariateSurfaceInterpolationScheme<scalar>::fieldTable fields;
 
-forAll (Y, i)
+forAll(Y, i)
 {
     fields.add(Y[i]);
 }
-fields.add(h);
+fields.add(hs);
+
+DimensionedField<scalar, volMesh> chemistrySh
+(
+    IOobject
+    (
+        "chemistry::Sh",
+        runTime.timeName(),
+        mesh,
+        IOobject::NO_READ,
+        IOobject::NO_WRITE
+    ),
+    mesh,
+    dimensionedScalar("chemistrySh", dimEnergy/dimTime/dimVolume, 0.0)
+);

applications/solvers/combustion/dieselEngineFoam/dieselEngineFoam.C

@@ -23,7 +23,7 @@ License
     Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
 
 Application
-    dieselFoam
+    dieselEngineFoam
 
 Description
     Solver for diesel engine spray and combustion.
@@ -103,13 +103,15 @@ int main(int argc, char *argv[])
             kappa = (runTime.deltaT() + tc)/(runTime.deltaT() + tc + tk);
         }
 
+        chemistrySh = kappa*chemistry.Sh()();
+
         #include "rhoEqn.H"
         #include "UEqn.H"
 
         for (label ocorr=1; ocorr <= nOuterCorr; ocorr++)
         {
             #include "YEqn.H"
-            #include "hEqn.H"
+            #include "hsEqn.H"
 
             // --- PISO loop
             for (int corr=1; corr<=nCorr; corr++)

applications/solvers/combustion/dieselEngineFoam/hEqn.H

@@ -1,13 +0,0 @@
-{
-    solve
-    (
-        fvm::ddt(rho, h)
-      + mvConvection->fvmDiv(phi, h)
-      - fvm::laplacian(turbulence->alphaEff(), h)
-     ==
-       DpDt
-     + dieselSpray.heatTransferSource()
-    );
-
-    thermo.correct();
-}

applications/solvers/combustion/dieselEngineFoam/hsEqn.H

@@ -0,0 +1,14 @@
+{
+    solve
+    (
+        fvm::ddt(rho, hs)
+      + mvConvection->fvmDiv(phi, hs)
+      - fvm::laplacian(turbulence->alphaEff(), hs)
+     ==
+       DpDt
+     + dieselSpray.heatTransferSource()().dimensionedInternalField()
+     + chemistrySh
+    );
+
+    thermo.correct();
+}

applications/solvers/combustion/dieselFoam/dieselFoam.C

@@ -100,7 +100,7 @@ int main(int argc, char *argv[])
         for (label ocorr=1; ocorr <= nOuterCorr; ocorr++)
         {
             #include "YEqn.H"
-            #include "hEqn.H"
+            #include "hsEqn.H"
 
             // --- PISO loop
             for (int corr=1; corr<=nCorr; corr++)

applications/solvers/combustion/fireFoam/Allwclean

@@ -0,0 +1,8 @@
+#!/bin/sh
+cd ${0%/*} || exit 1    # run from this directory
+set -x
+
+wclean libso combustionModels
+wclean
+
+# ----------------------------------------------------------------- end-of-file

applications/solvers/combustion/fireFoam/Allwmake

@@ -0,0 +1,8 @@
+#!/bin/sh
+cd ${0%/*} || exit 1    # run from this directory
+set -x
+
+wmake libso combustionModels
+wmake
+
+# ----------------------------------------------------------------- end-of-file

applications/solvers/combustion/fireFoam/Make/files

@@ -0,0 +1,4 @@
+fireFoam.C
+
+EXE = $(FOAM_APPBIN)/fireFoam
+

applications/solvers/combustion/fireFoam/Make/options

@@ -0,0 +1,19 @@
+EXE_INC = \
+    -I./combustionModels/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/specie/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/reactionThermo/lnInclude \
+    -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/turbulenceModels/compressible/turbulenceModel \
+    -I$(LIB_SRC)/thermophysicalModels/radiation/lnInclude
+
+
+EXE_LIBS = \
+    -lcombustionModels \
+    -lspecie \
+    -lreactionThermophysicalModels \
+    -lbasicThermophysicalModels \
+    -lfiniteVolume \
+    -lcompressibleLESModels \
+    -lcompressibleRASModels \
+    -lradiation

applications/solvers/combustion/fireFoam/UEqn.H

@@ -0,0 +1,36 @@
+fvVectorMatrix UEqn
+(
+    fvm::ddt(rho, U)
+  + fvm::div(phi, U)
+  + turbulence->divDevRhoReff(U)
+);
+
+UEqn.relax();
+
+if (oCorr == nOuterCorr - 1)
+{
+    solve
+    (
+        UEqn
+      ==
+        fvc::reconstruct
+        (
+            fvc::interpolate(rho)*(g & mesh.Sf())
+          - fvc::snGrad(p)*mesh.magSf()
+        ),
+        mesh.solver("UFinal")
+    );
+}
+else
+{
+    solve
+    (
+        UEqn
+      ==
+        fvc::reconstruct
+        (
+            fvc::interpolate(rho)*(g & mesh.Sf())
+          - fvc::snGrad(p)*mesh.magSf()
+        )
+    );
+}

applications/solvers/combustion/fireFoam/combustionModels/Make/files

@@ -0,0 +1,9 @@
+combustionModel/combustionModel.C
+combustionModel/newCombustionModel.C
+
+infinitelyFastChemistry/infinitelyFastChemistry.C
+
+noCombustion/noCombustion.C
+
+LIB = $(FOAM_LIBBIN)/libcombustionModels
+

applications/solvers/combustion/fireFoam/combustionModels/Make/options

@@ -0,0 +1,11 @@
+EXE_INC = \
+    -I../sensibleEnthalpyCombustionThermophysicalModels/basic/lnInclude \
+    -I../sensibleEnthalpyCombustionThermophysicalModels/reactionThermo/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/reactionThermo/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/specie/lnInclude \
+    -I$(FOAM_SRC)/turbulenceModels/compressible/turbulenceModel \
+    -I$(FOAM_SRC)/finiteVolume/lnInclude
+
+LIB_LIBS = \
+    -lfiniteVolume

applications/solvers/combustion/fireFoam/combustionModels/combustionModel/combustionModel.C

@@ -0,0 +1,109 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-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
+
+\*---------------------------------------------------------------------------*/
+
+#include "combustionModel.H"
+#include "fvm.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    defineTypeNameAndDebug(combustionModel, 0);
+    defineRunTimeSelectionTable(combustionModel, dictionary);
+};
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::combustionModel::combustionModel
+(
+    const dictionary& combustionProperties,
+    const hsCombustionThermo& thermo,
+    const compressible::turbulenceModel& turbulence,
+    const surfaceScalarField& phi,
+    const volScalarField& rho
+)
+:
+    combustionModelCoeffs_
+    (
+        combustionProperties.subDict
+        (
+            word(combustionProperties.lookup("combustionModel")) + "Coeffs"
+        )
+    ),
+    thermo_(thermo),
+    turbulence_(turbulence),
+    mesh_(phi.mesh()),
+    phi_(phi),
+    rho_(rho),
+    stoicRatio_(thermo.lookup("stoichiometricAirFuelMassRatio")),
+    s_(thermo.lookup("stoichiometricOxygenFuelMassRatio")),
+    qFuel_(thermo_.lookup("qFuel")),
+    composition_(thermo.composition())
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructors * * * * * * * * * * * * * * * //
+
+Foam::combustionModel::~combustionModel()
+{}
+
+
+// * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
+
+Foam::tmp<Foam::fvScalarMatrix>
+Foam::combustionModel::combustionModel::R(volScalarField& fu) const
+{
+    const basicMultiComponentMixture& composition = thermo_.composition();
+    const volScalarField& ft = composition.Y("ft");
+    volScalarField fres = composition.fres(ft, stoicRatio_.value());
+    volScalarField wFuelNorm = this->wFuelNorm()*pos(fu - fres);
+
+    return wFuelNorm*fres - fvm::Sp(wFuelNorm, fu);
+}
+
+
+Foam::tmp<Foam::volScalarField> Foam::combustionModel::combustionModel::dQ
+(
+    const fvScalarMatrix& Rfu
+) const
+{
+    const basicMultiComponentMixture& composition = thermo_.composition();
+    const volScalarField& fu = composition.Y("fu");
+
+    return (-qFuel_)*(Rfu & fu);
+}
+
+
+bool Foam::combustionModel::read(const dictionary& combustionProperties)
+{
+    combustionModelCoeffs_ = combustionProperties.subDict(type() + "Coeffs");
+
+    return true;
+}
+
+
+// ************************************************************************* //

applications/solvers/combustion/fireFoam/combustionModels/combustionModel/combustionModel.H

@@ -0,0 +1,210 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-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
+
+Class
+    Foam::combustionModel
+
+Description
+    Base class for all non-premixed combustion models.
+
+SourceFiles
+    combustionModel.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef combustionModel_H
+#define combustionModel_H
+
+#include "IOdictionary.H"
+#include "hsCombustionThermo.H"
+#include "turbulenceModel.H"
+#include "multivariateSurfaceInterpolationScheme.H"
+#include "runTimeSelectionTables.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                        Class combustionModel Declaration
+\*---------------------------------------------------------------------------*/
+
+class combustionModel
+{
+
+protected:
+
+    // Protected data
+
+        //- Dictionary of coefficients for the particular model
+        dictionary combustionModelCoeffs_;
+
+        //- Reference to the thermodynamic
+        const hsCombustionThermo& thermo_;
+
+        //- Reference to the turbulence model
+        const compressible::turbulenceModel& turbulence_;
+
+        //- Reference to the mesh database
+        const fvMesh& mesh_;
+
+        //- Reference to mass-flux field
+        const surfaceScalarField& phi_;
+
+        //- Reference to the density field
+        const volScalarField& rho_;
+
+        //- Stoichiometric air-fuel mass ratio
+        dimensionedScalar stoicRatio_;
+
+        //- Stoichiometric oxygen-fuel mass ratio
+        dimensionedScalar s_;
+
+        //- Heat of combustion (J/Kg)
+        dimensionedScalar qFuel_;
+
+
+private:
+
+    // Private Member Functions
+
+        //- Disallow copy construct
+        combustionModel(const combustionModel&);
+
+        //- Disallow default bitwise assignment
+        void operator=(const combustionModel&);
+
+        const basicMultiComponentMixture& composition_;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("combustionModel");
+
+
+    // Declare run-time constructor selection table
+
+        declareRunTimeSelectionTable
+        (
+            autoPtr,
+            combustionModel,
+            dictionary,
+            (
+                const dictionary& combustionProperties,
+                const hsCombustionThermo& thermo,
+                const compressible::turbulenceModel& turbulence,
+                const surfaceScalarField& phi,
+                const volScalarField& rho
+            ),
+            (
+                combustionProperties,
+                thermo,
+                turbulence,
+                phi,
+                rho
+            )
+        );
+
+
+    // Selectors
+
+        //- Return a reference to the selected combustion model
+        static autoPtr<combustionModel> New
+        (
+            const dictionary& combustionProperties,
+            const hsCombustionThermo& thermo,
+            const compressible::turbulenceModel& turbulence,
+            const surfaceScalarField& phi,
+            const volScalarField& rho
+        );
+
+
+    // Constructors
+
+        //- Construct from components
+        combustionModel
+        (
+            const dictionary& combustionProperties,
+            const hsCombustionThermo& thermo,
+            const compressible::turbulenceModel& turbulence,
+            const surfaceScalarField& phi,
+            const volScalarField& rho
+        );
+
+
+    //- Destructor
+    virtual ~combustionModel();
+
+
+    // Member Functions
+
+        // Access functions
+
+            //- Access composition
+            const basicMultiComponentMixture& composition() const
+            {
+                return composition_;
+            }
+
+            //- Access combustion dictionary
+            const dictionary combustionModelCoeffs() const
+            {
+                return combustionModelCoeffs_;
+            }
+
+            //- Access heat of combustion
+            const dimensionedScalar qFuel() const
+            {
+                return qFuel_;
+            }
+
+        //- Return normalised consumption rate of (fu - fres)
+        virtual tmp<volScalarField> wFuelNorm() const = 0;
+
+        //- Fuel consumption rate matrix i.e. source-term for the fuel equation
+        virtual tmp<fvScalarMatrix> R(volScalarField& fu) const;
+
+        //- Heat-release rate calculated from the given
+        //  fuel consumption rate matrix
+        virtual tmp<volScalarField> dQ(const fvScalarMatrix& Rfu) const;
+
+        //- Correct combustion rate
+        virtual void correct() = 0;
+
+        //- Update properties from given dictionary
+        virtual bool read(const dictionary& combustionProperties) = 0;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

applications/solvers/combustion/fireFoam/combustionModels/combustionModel/newCombustionModel.C

@@ -0,0 +1,67 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-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
+
+\*---------------------------------------------------------------------------*/
+
+#include "combustionModel.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+Foam::autoPtr<Foam::combustionModel> Foam::combustionModel::New
+(
+    const dictionary& combustionProperties,
+    const hsCombustionThermo& thermo,
+    const compressible::turbulenceModel& turbulence,
+    const surfaceScalarField& phi,
+    const volScalarField& rho
+)
+{
+    word combustionModelTypeName = combustionProperties.lookup
+    (
+        "combustionModel"
+    );
+
+    Info<< "Selecting combustion model " << combustionModelTypeName << endl;
+
+    dictionaryConstructorTable::iterator cstrIter =
+        dictionaryConstructorTablePtr_->find(combustionModelTypeName);
+
+    if (cstrIter == dictionaryConstructorTablePtr_->end())
+    {
+        FatalErrorIn
+        (
+            "combustionModel::New"
+        )   << "Unknown combustionModel type "
+            << combustionModelTypeName << endl << endl
+            << "Valid  combustionModels are : " << endl
+            << dictionaryConstructorTablePtr_->toc()
+            << exit(FatalError);
+    }
+
+    return autoPtr<combustionModel>
+        (cstrIter()(combustionProperties, thermo, turbulence, phi, rho));
+}
+
+
+// ************************************************************************* //

applications/solvers/combustion/fireFoam/combustionModels/infinitelyFastChemistry/infinitelyFastChemistry.C

@@ -24,68 +24,70 @@ License
 
 \*---------------------------------------------------------------------------*/
 
-#include "simpleSplineEdge.H"
+#include "infinitelyFastChemistry.H"
 #include "addToRunTimeSelectionTable.H"
 
-
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
 namespace Foam
 {
-    defineTypeNameAndDebug(simpleSplineEdge, 0);
-    addToRunTimeSelectionTable(curvedEdge, simpleSplineEdge, Istream);
-}
+namespace combustionModels
+{
+    defineTypeNameAndDebug(infinitelyFastChemistry, 0);
+    addToRunTimeSelectionTable
+    (
+        combustionModel,
+        infinitelyFastChemistry,
+        dictionary
+    );
+};
+};
 
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
-Foam::simpleSplineEdge::simpleSplineEdge
+Foam::combustionModels::infinitelyFastChemistry::infinitelyFastChemistry
 (
-    const pointField& points,
-    const label start,
-    const label end,
-    const pointField& otherknots
+    const dictionary& combustionProperties,
+    const hsCombustionThermo& thermo,
+    const compressible::turbulenceModel& turbulence,
+    const surfaceScalarField& phi,
+    const volScalarField& rho
 )
 :
-    curvedEdge(points, start, end),
-    BSpline(appendEndPoints(points, start, end, otherknots))
+    combustionModel(combustionProperties, thermo, turbulence, phi, rho),
+    C_(readScalar(combustionModelCoeffs_.lookup("C")))
 {}
 
 
-Foam::simpleSplineEdge::simpleSplineEdge
-(
-    const pointField& points,
-    const label start,
-    const label end,
-    const pointField& otherknots,
-    const vector& fstend,
-    const vector& sndend
-)
-:
-    curvedEdge(points, start, end),
-    BSpline(appendEndPoints(points, start, end, otherknots), fstend, sndend)
+// * * * * * * * * * * * * * * * * Destructors * * * * * * * * * * * * * * * //
+
+Foam::combustionModels::infinitelyFastChemistry::~infinitelyFastChemistry()
 {}
 
 
-Foam::simpleSplineEdge::simpleSplineEdge(const pointField& points, Istream& is)
-:
-    curvedEdge(points, is),
-    BSpline(appendEndPoints(points, start_, end_, pointField(is)))
+// * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
+
+void Foam::combustionModels::infinitelyFastChemistry::correct()
 {}
 
 
-// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
-
-Foam::point Foam::simpleSplineEdge::position(const scalar mu) const
+Foam::tmp<Foam::volScalarField>
+Foam::combustionModels::infinitelyFastChemistry::wFuelNorm() const
 {
-    return BSpline::position(mu);
+    return rho_/(mesh_.time().deltaT()*C_);
 }
 
 
-Foam::scalar Foam::simpleSplineEdge::length() const
+bool Foam::combustionModels::infinitelyFastChemistry::read
+(
+    const dictionary& combustionProperties
+)
 {
-    notImplemented("simpleSplineEdge::length() const");
-    return 1.0;
+    combustionModel::read(combustionProperties);
+    combustionModelCoeffs_.lookup("C") >> C_ ;
+
+    return true;
 }
 
 

applications/solvers/combustion/fireFoam/combustionModels/infinitelyFastChemistry/infinitelyFastChemistry.H

@@ -23,96 +23,93 @@ License
     Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
 
 Class
-    Foam::simpleSplineEdge
+    Foam::combustionModels::infinitelyFastChemistry
 
 Description
-    A curvedEdge interface for B-splines.
+    Simple infinitely fast chemistry combustion model based on the principle
+    mixed is burnt. Additional parameter C is used to distribute the heat
+    release rate.in time
 
 SourceFiles
-    simpleSplineEdge.C
+    infinitelyFastChemistry.C
 
 \*---------------------------------------------------------------------------*/
 
-#ifndef simpleSplineEdge_H
-#define simpleSplineEdge_H
+#ifndef infinitelyFastChemistry_H
+#define infinitelyFastChemistry_H
 
-#include "curvedEdge.H"
-#include "BSpline.H"
+#include "fvc.H"
+#include "combustionModel.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 namespace Foam
 {
+namespace combustionModels
+{
 
 /*---------------------------------------------------------------------------*\
-                      Class simpleSplineEdge Declaration
+                Class infinitelyFastChemistry Declaration
 \*---------------------------------------------------------------------------*/
 
-class simpleSplineEdge
+class infinitelyFastChemistry
 :
-    public curvedEdge,
-    public BSpline
+    public combustionModel
 {
+    // Private data
+
+        //- Model constant
+        scalar C_;
+
     // Private Member Functions
 
-        //- Disallow default bitwise copy construct
-        simpleSplineEdge(const simpleSplineEdge&);
+        //- Disallow copy construct
+        infinitelyFastChemistry(const infinitelyFastChemistry&);
 
         //- Disallow default bitwise assignment
-        void operator=(const simpleSplineEdge&);
+        void operator=(const infinitelyFastChemistry&);
 
 
 public:
 
     //- Runtime type information
-    TypeName("simpleSpline");
+    TypeName("infinitelyFastChemistry");
 
 
     // Constructors
 
         //- Construct from components
-        simpleSplineEdge
+        infinitelyFastChemistry
         (
-            const pointField&,
-            const label start,
-            const label end,
-            const pointField& otherKnots
+            const dictionary& combustionProperties,
+            const hsCombustionThermo& thermo,
+            const compressible::turbulenceModel& turbulence,
+            const surfaceScalarField& phi,
+            const volScalarField& rho
         );
 
-        //- Construct from components
-        simpleSplineEdge
-        (
-            const pointField&,
-            const label start,
-            const label end,
-            const pointField& otherKnots,
-            const vector& fstend,
-            const vector& sndend
-        );
-
-        //- Construct from Istream setting pointsList
-        simpleSplineEdge(const pointField&, Istream&);
-
 
     // Destructor
 
-        virtual ~simpleSplineEdge()
-        {}
+        virtual ~infinitelyFastChemistry();
 
 
     // Member Functions
 
-        //- Return the point position corresponding to the curve parameter
-        //  0 <= lambda <= 1
-        point position(const scalar) const;
+        //- Update properties from given dictionary
+        virtual bool read(const dictionary& combustionProperties);
 
-        //- Return the length of the simple spline curve
-        scalar length() const;
+        //- Correct combustion rate
+        virtual void correct();
+
+        //- Return normalised consumption rate of (fu - fres)
+        virtual tmp<volScalarField> wFuelNorm() const;
 };
 
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
+} // End namespace combustionModels
 } // End namespace Foam
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

applications/solvers/combustion/fireFoam/combustionModels/noCombustion/noCombustion.C

@@ -0,0 +1,103 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-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
+
+\*---------------------------------------------------------------------------*/
+
+#include "noCombustion.H"
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace combustionModels
+{
+    defineTypeNameAndDebug(noCombustion, 0);
+    addToRunTimeSelectionTable
+    (
+        combustionModel,
+        noCombustion,
+        dictionary
+    );
+};
+};
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::combustionModels::noCombustion::noCombustion
+(
+    const dictionary& combustionProperties,
+    const hsCombustionThermo& thermo,
+    const compressible::turbulenceModel& turbulence,
+    const surfaceScalarField& phi,
+    const volScalarField& rho
+)
+:
+    combustionModel(combustionProperties, thermo, turbulence, phi, rho)
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructors * * * * * * * * * * * * * * * //
+
+Foam::combustionModels::noCombustion::~noCombustion()
+{}
+
+
+void Foam::combustionModels::noCombustion::correct()
+{}
+
+
+Foam::tmp<Foam::volScalarField>
+Foam::combustionModels::noCombustion::wFuelNorm() const
+{
+    return tmp<Foam::volScalarField>
+    (
+        new volScalarField
+        (
+            IOobject
+            (
+                "wFuelNorm",
+                mesh_.time().timeName(),
+                mesh_,
+                IOobject::NO_READ,
+                IOobject::NO_WRITE
+            ),
+            mesh_,
+            dimensionedScalar("wFuelNorm", dimMass/dimTime/pow3(dimLength), 0.0)
+        )
+    );
+}
+
+
+bool Foam::combustionModels::noCombustion::read
+(
+    const dictionary& combustionProperties
+)
+{
+    return combustionModel::read(combustionProperties);
+}
+
+
+// ************************************************************************* //

applications/solvers/combustion/fireFoam/combustionModels/noCombustion/noCombustion.H

@@ -23,91 +23,87 @@ License
     Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
 
 Class
-    Foam::polySplineEdge
+    Foam::combustionModel::noCombustion
 
 Description
-    A curvedEdge interface for B-splines.
+    No combustion
 
 SourceFiles
-    polySplineEdge.C
+    noCombustion.C
 
 \*---------------------------------------------------------------------------*/
 
-#ifndef polySplineEdge_H
-#define polySplineEdge_H
+#ifndef noCombustion_H
+#define noCombustion_H
 
-#include "curvedEdge.H"
-#include "polyLine.H"
+#include "combustionModel.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 namespace Foam
 {
+namespace combustionModels
+{
 
 /*---------------------------------------------------------------------------*\
-                       Class polySplineEdge Declaration
+                          Class noCombustion Declaration
 \*---------------------------------------------------------------------------*/
 
-class polySplineEdge
+class noCombustion
 :
-    public curvedEdge,
-    public polyLine
+    public combustionModel
 {
     // Private data
 
-        pointField otherKnots_;
 
+    // Private Member Functions
 
-    // Private member functions
+        //- Disallow copy construct
+        noCombustion(const noCombustion&);
+
+        //- Disallow default bitwise assignment
+        void operator=(const noCombustion&);
 
-        pointField intervening
-        (
-            const pointField& otherKnots,
-            const label nBetweenKnots,
-            const vector&,
-            const vector&
-        );
 
 public:
 
     //- Runtime type information
-    TypeName("polySpline");
+    TypeName("noCombustion");
 
 
     // Constructors
 
         //- Construct from components
-        polySplineEdge
+        noCombustion
         (
-            const pointField&,
-            const label start,
-            const label end,
-            const pointField& otherKnots,
-            const label nInterKnots = 20
+            const dictionary& combustionProperties,
+            const hsCombustionThermo& thermo,
+            const compressible::turbulenceModel& turbulence,
+            const surfaceScalarField& phi,
+            const volScalarField& rho
         );
 
-        //- Construct from Istream setting pointsList
-        polySplineEdge(const pointField&, Istream&);
 
+    // Destructor
 
-    //- Destructor
-        virtual ~polySplineEdge()
-        {}
-
+        virtual ~noCombustion();
 
     // Member Functions
 
-        //- Return the point position corresponding to the curve parameter
-        //  0 <= lambda <= 1
-        point position(const scalar mu) const;
+        //- Update properties from given dictionary
+        virtual bool read(const dictionary& combustionProperties);
 
-        //- Return the length of the curve
-        scalar length() const;
+        //- Correct combustion rate
+        virtual void correct();
+
+        //- Return normalised consumption rate of (fu - fres)
+        virtual tmp<volScalarField> wFuelNorm() const;
 };
 
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
+} // End namespace combustionModels
 } // End namespace Foam
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

applications/solvers/combustion/fireFoam/createFields.H

@@ -0,0 +1,130 @@
+Info<< "Reading thermophysical properties\n" << endl;
+
+autoPtr<hsCombustionThermo> pThermo
+(
+    hsCombustionThermo::New(mesh)
+);
+
+hsCombustionThermo& thermo = pThermo();
+
+basicMultiComponentMixture& composition = thermo.composition();
+
+volScalarField rho
+(
+    IOobject
+    (
+        "rho",
+        runTime.timeName(),
+        mesh,
+        IOobject::NO_READ,
+        IOobject::AUTO_WRITE
+    ),
+    thermo.rho()
+);
+
+dimensionedScalar stoicRatio
+(
+    thermo.lookup("stoichiometricAirFuelMassRatio")
+);
+
+volScalarField& p = thermo.p();
+volScalarField& hs = thermo.hs();
+
+const volScalarField& psi = thermo.psi();
+
+volScalarField& ft = composition.Y("ft");
+volScalarField& fu = composition.Y("fu");
+
+Info<< "Reading field U\n" << endl;
+
+volVectorField U
+(
+    IOobject
+    (
+        "U",
+        runTime.timeName(),
+        mesh,
+        IOobject::MUST_READ,
+        IOobject::AUTO_WRITE
+    ),
+    mesh
+);
+
+#include "compressibleCreatePhi.H"
+
+Info<< "Creating turbulence model\n" << endl;
+autoPtr<compressible::turbulenceModel> turbulence
+(
+    compressible::turbulenceModel::New(rho, U, phi, thermo)
+);
+
+IOdictionary combustionProperties
+(
+    IOobject
+    (
+        "combustionProperties",
+        runTime.constant(),
+        mesh,
+        IOobject::MUST_READ,
+        IOobject::NO_WRITE
+    )
+);
+
+Info<< "Creating combustion model\n" << endl;
+autoPtr<combustionModel::combustionModel> combustion
+(
+    combustionModel::combustionModel::New
+    (
+        combustionProperties,
+        thermo,
+        turbulence(),
+        phi,
+        rho
+    )
+);
+
+volScalarField dQ
+(
+    IOobject
+    (
+        "dQ",
+        runTime.timeName(),
+        mesh,
+        IOobject::NO_READ,
+        IOobject::AUTO_WRITE
+    ),
+    mesh,
+    dimensionedScalar("dQ", dimMass/pow3(dimTime)/dimLength, 0.0)
+);
+
+
+Info<< "Creating field DpDt\n" << endl;
+volScalarField DpDt =
+    fvc::DDt(surfaceScalarField("phiU", phi/fvc::interpolate(rho)), p);
+
+
+Info<< "Calculating field g.h\n" << endl;
+volScalarField gh("gh", g & mesh.C());
+
+surfaceScalarField ghf("gh", g & mesh.Cf());
+
+p += rho*gh;
+
+thermo.correct();
+
+dimensionedScalar initialMass = fvc::domainIntegrate(rho);
+
+
+multivariateSurfaceInterpolationScheme<scalar>::fieldTable fields;
+
+if (composition.contains("ft"))
+{
+    fields.add(composition.Y("ft"));
+}
+
+if (composition.contains("fu"))
+{
+    fields.add(composition.Y("fu"));
+}
+
+fields.add(hs);

applications/solvers/combustion/fireFoam/fireFoam.C

@@ -0,0 +1,103 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-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
+    fireFoam
+
+Description
+    Transient Solver for Fires and turbulent diffusion flames
+
+\*---------------------------------------------------------------------------*/
+
+#include "fvCFD.H"
+#include "hsCombustionThermo.H"
+#include "turbulenceModel.H"
+#include "combustionModel.H"
+
+#include "radiationModel.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+int main(int argc, char *argv[])
+{
+    #include "setRootCase.H"
+    #include "createTime.H"
+    #include "createMesh.H"
+    #include "readGravitationalAcceleration.H"
+    #include "initContinuityErrs.H"
+    #include "createFields.H"
+    #include "createRadiationModel.H"
+    #include "readTimeControls.H"
+    #include "compressibleCourantNo.H"
+    #include "setInitialDeltaT.H"
+
+    // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+    Info<< "\nStarting time loop\n" << endl;
+
+    while (runTime.run())
+    {
+        #include "readPISOControls.H"
+        #include "readTimeControls.H"
+        #include "compressibleCourantNo.H"
+        #include "setDeltaT.H"
+
+        runTime++;
+        Info<< "Time = " << runTime.timeName() << nl << endl;
+
+        #include "rhoEqn.H"
+
+        // --- Pressure-velocity PIMPLE corrector loop
+        for (int oCorr=0; oCorr<nOuterCorr; oCorr++)
+        {
+            #include "UEqn.H"
+
+            #include "ftEqn.H"
+            #include "fuhsEqn.H"
+
+            // --- PISO loop
+            for (int corr=0; corr<nCorr; corr++)
+            {
+                #include "pEqn.H"
+            }
+        }
+
+        turbulence->correct();
+
+        rho = thermo.rho();
+
+        runTime.write();
+
+        Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
+            << "  ClockTime = " << runTime.elapsedClockTime() << " s"
+            << nl << endl;
+
+    }
+
+    Info<< "End\n" << endl;
+
+    return 0;
+}
+
+// ************************************************************************* //

applications/solvers/combustion/fireFoam/ftEqn.H

@@ -0,0 +1,25 @@
+tmp<fv::convectionScheme<scalar> > mvConvection
+(
+    fv::convectionScheme<scalar>::New
+    (
+        mesh,
+        fields,
+        phi,
+        mesh.divScheme("div(phi,ft_b_h)")
+    )
+);
+
+{
+    fvScalarMatrix ftEqn
+    (
+        fvm::ddt(rho, ft)
+      + mvConvection->fvmDiv(phi, ft)
+      - fvm::laplacian(turbulence->alphaEff(), ft)
+    );
+
+    ftEqn.relax();
+    ftEqn.solve();
+}
+
+Info<< "max(ft) = " << max(ft).value() << endl;
+Info<< "min(ft) = " << min(ft).value() << endl;

applications/solvers/combustion/fireFoam/fuhsEqn.H

@@ -0,0 +1,47 @@
+{
+    // Solve fuel equation
+    // ~~~~~~~~~~~~~~~~~~~
+    fvScalarMatrix R = combustion->R(fu);
+
+    {
+        fvScalarMatrix fuEqn
+        (
+            fvm::ddt(rho, fu)
+          + mvConvection->fvmDiv(phi, fu)
+          - fvm::laplacian(turbulence->alphaEff(), fu)
+         ==
+            R
+        );
+
+        fuEqn.relax();
+        fuEqn.solve();
+    }
+
+    Info<< "max(fu) = " << max(fu).value() << endl;
+    Info<< "min(fu) = " << min(fu).value() << endl;
+
+
+    // Solve sensible enthalpy equation
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+    radiation->correct();
+    dQ = combustion->dQ(R);
+
+    {
+        fvScalarMatrix hsEqn
+        (
+            fvm::ddt(rho, hs)
+          + mvConvection->fvmDiv(phi,hs)
+          - fvm::laplacian(turbulence->alphaEff(), hs)
+         ==
+            DpDt
+          + dQ
+          + radiation->Shs(thermo)
+        );
+
+        hsEqn.relax();
+        hsEqn.solve();
+    }
+
+    thermo.correct();
+    combustion->correct();
+}

applications/solvers/combustion/fireFoam/pEqn.H

@@ -0,0 +1,64 @@
+bool closedVolume = false;
+
+rho = thermo.rho();
+
+volScalarField rUA = 1.0/UEqn.A();
+surfaceScalarField rhorUAf("(rho*(1|A(U)))", fvc::interpolate(rho*rUA));
+U = rUA*UEqn.H();
+
+surfaceScalarField phiU
+(
+    fvc::interpolate(rho)
+   *(
+        (fvc::interpolate(U) & mesh.Sf())
+      + fvc::ddtPhiCorr(rUA, rho, U, phi)
+    )
+);
+
+phi = phiU + rhorUAf*fvc::interpolate(rho)*(g & mesh.Sf());
+
+for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
+{
+       surfaceScalarField rhorUAf = fvc::interpolate(rho*rUA);
+
+        fvScalarMatrix pEqn
+        (
+            fvm::ddt(psi,p)
+          + fvc::div(phi)
+          - fvm::laplacian(rhorUAf, p)
+        );
+
+        closedVolume = p.needReference();
+
+        if (corr == nCorr-1 && nonOrth == nNonOrthCorr)
+        {
+            pEqn.solve(mesh.solver(p.name() + "Final"));
+        }
+        else
+        {
+            pEqn.solve(mesh.solver(p.name()));
+        }
+
+        if (nonOrth == nNonOrthCorr)
+        {
+            phi += pEqn.flux();
+        }
+}
+
+DpDt = fvc::DDt(surfaceScalarField("phiU", phi/fvc::interpolate(rho)), p);
+
+#include "rhoEqn.H"
+#include "compressibleContinuityErrs.H"
+
+U += rUA*fvc::reconstruct((phi - phiU)/rhorUAf);
+U.correctBoundaryConditions();
+
+// For closed-volume cases adjust the pressure and density levels
+// to obey overall mass continuity
+if (closedVolume)
+{
+    p +=
+        (initialMass - fvc::domainIntegrate(thermo.psi()*p))
+       /fvc::domainIntegrate(thermo.psi());
+    rho = thermo.rho();
+}

applications/solvers/combustion/reactingFoam/Make/options

@@ -1,5 +1,4 @@
 EXE_INC = \
-    -I../XiFoam \
     -I$(LIB_SRC)/turbulenceModels/compressible/turbulenceModel \
     -I$(LIB_SRC)/thermophysicalModels/specie/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/reactionThermo/lnInclude \

applications/solvers/combustion/reactingFoam/UEqn.H

@@ -0,0 +1,15 @@
+    fvVectorMatrix UEqn
+    (
+        fvm::ddt(rho, U)
+      + fvm::div(phi, U)
+      + turbulence->divDevRhoReff(U)
+     ==
+        rho*g
+    );
+
+    UEqn.relax();
+
+    if (momentumPredictor)
+    {
+        solve(UEqn == -fvc::grad(p));
+    }

applications/solvers/combustion/reactingFoam/chemistry.H

@@ -21,4 +21,6 @@
     {
         kappa = 1.0;
     }
+
+    chemistrySh = kappa*chemistry.Sh()();
 }

applications/solvers/combustion/reactingFoam/createFields.H

@@ -5,7 +5,7 @@ autoPtr<psiChemistryModel> pChemistry
 );
 psiChemistryModel& chemistry = pChemistry();
 
-hCombustionThermo& thermo = chemistry.thermo();
+hsCombustionThermo& thermo = chemistry.thermo();
 
 basicMultiComponentMixture& composition = thermo.composition();
 PtrList<volScalarField>& Y = composition.Y();
@@ -40,8 +40,8 @@ volVectorField U
 
 volScalarField& p = thermo.p();
 const volScalarField& psi = thermo.psi();
-volScalarField& h = thermo.h();
-
+volScalarField& hs = thermo.hs();
+const volScalarField& T = thermo.T();
 
 #include "compressibleCreatePhi.H"
 
@@ -77,9 +77,22 @@ volScalarField DpDt =
 
 multivariateSurfaceInterpolationScheme<scalar>::fieldTable fields;
 
-forAll (Y, i)
+forAll(Y, i)
 {
     fields.add(Y[i]);
 }
-fields.add(h);
+fields.add(hs);
 
+DimensionedField<scalar, volMesh> chemistrySh
+(
+    IOobject
+    (
+        "chemistry::Sh",
+        runTime.timeName(),
+        mesh,
+        IOobject::NO_READ,
+        IOobject::NO_WRITE
+    ),
+    mesh,
+    dimensionedScalar("chemistrySh", dimEnergy/dimTime/dimVolume, 0.0)
+);

applications/solvers/combustion/reactingFoam/hsEqn.H

@@ -0,0 +1,20 @@
+{
+    fvScalarMatrix hsEqn
+    (
+        fvm::ddt(rho, hs)
+      + mvConvection->fvmDiv(phi, hs)
+      - fvm::laplacian(turbulence->alphaEff(), hs)
+//      - fvm::laplacian(turbulence->muEff(), hs)  // unit lewis no.
+     ==
+        DpDt
+      + chemistrySh
+    );
+
+    hsEqn.relax();
+    hsEqn.solve();
+
+    thermo.correct();
+
+    Info<< "T gas min/max   = " << min(T).value() << ", "
+        << max(T).value() << endl;
+}

applications/solvers/combustion/reactingFoam/pEqn.H

@@ -0,0 +1,68 @@
+rho = thermo.rho();
+
+volScalarField rUA = 1.0/UEqn.A();
+U = rUA*UEqn.H();
+
+if (transonic)
+{
+    surfaceScalarField phid
+    (
+        "phid",
+        fvc::interpolate(psi)
+       *(
+            (fvc::interpolate(U) & mesh.Sf())
+          + fvc::ddtPhiCorr(rUA, rho, U, phi)
+        )
+    );
+
+    for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
+    {
+        fvScalarMatrix pEqn
+        (
+            fvm::ddt(psi, p)
+          + fvm::div(phid, p)
+          - fvm::laplacian(rho*rUA, p)
+        );
+
+        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
+        (
+            fvm::ddt(psi, p)
+          + fvc::div(phi)
+          - fvm::laplacian(rho*rUA, p)
+        );
+
+        pEqn.solve();
+
+        if (nonOrth == nNonOrthCorr)
+        {
+            phi += pEqn.flux();
+        }
+    }
+}
+
+#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/reactingFoam/reactingFoam.C

@@ -73,9 +73,7 @@ int main(int argc, char *argv[])
         {
             #include "UEqn.H"
             #include "YEqn.H"
-
-            #define Db turbulence->alphaEff()
-            #include "hEqn.H"
+            #include "hsEqn.H"
 
             // --- PISO loop
             for (int corr=1; corr<=nCorr; corr++)

applications/solvers/combustion/rhoReactingFoam/Make/options

@@ -1,5 +1,4 @@
 EXE_INC = \
-    -I../XiFoam \
     -I$(LIB_SRC)/turbulenceModels/compressible/turbulenceModel \
     -I$(LIB_SRC)/thermophysicalModels/specie/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/reactionThermo/lnInclude \

applications/solvers/combustion/rhoReactingFoam/UEqn.H

@@ -0,0 +1,15 @@
+    fvVectorMatrix UEqn
+    (
+        fvm::ddt(rho, U)
+      + fvm::div(phi, U)
+      + turbulence->divDevRhoReff(U)
+     ==
+        rho*g
+    );
+
+    UEqn.relax();
+
+    if (momentumPredictor)
+    {
+        solve(UEqn == -fvc::grad(p));
+    }

applications/solvers/combustion/rhoReactingFoam/chemistry.H

@@ -21,4 +21,6 @@
     {
         kappa = 1.0;
     }
+
+    chemistrySh = kappa*chemistry.Sh()();
 }

applications/solvers/combustion/rhoReactingFoam/createFields.H

@@ -5,7 +5,7 @@ autoPtr<rhoChemistryModel> pChemistry
 );
 rhoChemistryModel& chemistry = pChemistry();
 
-hReactionThermo& thermo = chemistry.thermo();
+hsReactionThermo& thermo = chemistry.thermo();
 
 basicMultiComponentMixture& composition = thermo.composition();
 PtrList<volScalarField>& Y = composition.Y();
@@ -40,7 +40,8 @@ volVectorField U
 
 volScalarField& p = thermo.p();
 const volScalarField& psi = thermo.psi();
-volScalarField& h = thermo.h();
+volScalarField& hs = thermo.hs();
+const volScalarField& T = thermo.T();
 
 
 #include "compressibleCreatePhi.H"
@@ -77,9 +78,22 @@ volScalarField DpDt =
 
 multivariateSurfaceInterpolationScheme<scalar>::fieldTable fields;
 
-forAll (Y, i)
+forAll(Y, i)
 {
     fields.add(Y[i]);
 }
-fields.add(h);
+fields.add(hs);
 
+DimensionedField<scalar, volMesh> chemistrySh
+(
+    IOobject
+    (
+        "chemistry::Sh",
+        runTime.timeName(),
+        mesh,
+        IOobject::NO_READ,
+        IOobject::NO_WRITE
+    ),
+    mesh,
+    dimensionedScalar("chemistrySh", dimEnergy/dimTime/dimVolume, 0.0)
+);

applications/solvers/combustion/rhoReactingFoam/hsEqn.H

@@ -0,0 +1,19 @@
+{
+    fvScalarMatrix hsEqn
+    (
+        fvm::ddt(rho, hs)
+      + mvConvection->fvmDiv(phi, hs)
+      - fvm::laplacian(turbulence->alphaEff(), hs)
+     ==
+        DpDt
+      + chemistrySh
+    );
+
+    hsEqn.relax();
+    hsEqn.solve();
+
+    thermo.correct();
+
+    Info<< "T gas min/max   = " << min(T).value() << ", "
+        << max(T).value() << endl;
+}

applications/solvers/combustion/rhoReactingFoam/rhoReactingFoam.C

@@ -74,7 +74,7 @@ int main(int argc, char *argv[])
         {
             #include "UEqn.H"
             #include "YEqn.H"
-            #include "hEqn.H"
+            #include "hsEqn.H"
 
             // --- PISO loop
             for (int corr=1; corr<=nCorr; corr++)

applications/solvers/compressible/rhoCentralFoam/Allwmake

@@ -2,7 +2,7 @@
 cd ${0%/*} || exit 1    # run from this directory
 set -x
 
-wmake libso BCs 
+wmake libso BCs
 wmake
 
 # ----------------------------------------------------------------- end-of-file

applications/solvers/compressible/rhoCentralFoam/BCs/Make/options

@@ -6,5 +6,5 @@ EXE_INC = \
 EXE_LIBS = \
     -lfiniteVolume \
     -lbasicThermophysicalModels \
-	-lspecie
+    -lspecie
 

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

@@ -22,8 +22,6 @@ License
     along with OpenFOAM; if not, write to the Free Software Foundation,
     Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
 
-Description
-
 \*---------------------------------------------------------------------------*/
 
 #include "maxwellSlipUFvPatchVectorField.H"

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

@@ -22,8 +22,6 @@ License
     along with OpenFOAM; if not, write to the Free Software Foundation,
     Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
 
-Description
-
 \*---------------------------------------------------------------------------*/
 
 #include "mixedFixedValueSlipFvPatchFields.H"

applications/solvers/compressible/rhoSonicFoam/rhoSonicFoam.C

@@ -84,7 +84,7 @@ int main(int argc, char *argv[])
         (
             fvm::ddt(rhoU)
           + fvm::div(phiv, rhoU)
-         == 
+         ==
           - fvc::grad(p)
         );
 

applications/solvers/compressible/rhopSonicFoam/Allwmake

@@ -2,7 +2,7 @@
 cd ${0%/*} || exit 1    # run from this directory
 set -x
 
-wmake libso BCs 
+wmake libso BCs
 wmake
 
 # ----------------------------------------------------------------- end-of-file

applications/solvers/compressible/rhopSonicFoam/BCs/p/inviscidWallPFvPatchScalarField.C

@@ -131,7 +131,7 @@ void inviscidWallPFvPatchScalarField::updateCoeffs()
     const fvPatchField<vector>& rhoUp =
         patch().lookupPatchField<volVectorField, vector>("rhoU");
 
-    const fvsPatchField<scalar>& phip = 
+    const fvsPatchField<scalar>& phip =
         patch().lookupPatchField<surfaceScalarField, scalar>("phi");
 
     const fvsPatchField<scalar>& rAp =
@@ -147,7 +147,7 @@ void inviscidWallPFvPatchScalarField::updateCoeffs()
 void inviscidWallPFvPatchScalarField::write(Ostream& os) const
 {
     fixedGradientFvPatchScalarField::write(os);
-    os.writeKeyword("fluxFraction") 
+    os.writeKeyword("fluxFraction")
         << fluxFraction_ << token::END_STATEMENT << nl;
     writeEntry("value", os);
 }

applications/solvers/compressible/rhopSonicFoam/BCs/rhoE/mixedRhoEFvPatchScalarField.C

@@ -145,7 +145,7 @@ void mixedRhoEFvPatchScalarField::updateCoeffs()
 
     const volScalarField& T = db().lookupObject<volScalarField>("T");
     const label patchi = patch().index();
-    fvPatchScalarField& Tp = 
+    fvPatchScalarField& Tp =
         const_cast<fvPatchScalarField&>(T.boundaryField()[patchi]);
 
     Tp.evaluate();
@@ -164,7 +164,7 @@ void mixedRhoEFvPatchScalarField::updateCoeffs()
     refGrad() =
         rhop*Cv.value()*Tp.snGrad()
       + (
-            refValue() 
+            refValue()
           - (0.5*rhop.patchInternalField()*
                 magSqr(rhoUp.patchInternalField()/rhop.patchInternalField()))
         )*patch().deltaCoeffs();

applications/solvers/compressible/rhopSonicFoam/BCs/rhoU/fixedRhoUFvPatchVectorField.H

@@ -66,7 +66,7 @@ public:
             const dictionary&
         );
 
-        //- Construct by mapping given fixedRhoUFvPatchVectorField 
+        //- Construct by mapping given fixedRhoUFvPatchVectorField
         //  onto a new patch
         fixedRhoUFvPatchVectorField
         (

applications/solvers/electromagnetics/electrostaticFoam/Make/files

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

applications/solvers/heatTransfer/chtMultiRegionFoam/Make/files

@@ -1,6 +1,10 @@
+regionProperties/regionProperties.C
+
 derivedFvPatchFields/solidWallHeatFluxTemperature/solidWallHeatFluxTemperatureFvPatchScalarField.C
+derivedFvPatchFields/solidWallMixedTemperatureCoupled/solidWallMixedTemperatureCoupledFvPatchScalarField.C
 
 fluid/compressibleCourantNo.C
+solid/solidRegionDiffNo.C
 
 chtMultiRegionFoam.C
 

applications/solvers/heatTransfer/chtMultiRegionFoam/Make/options

@@ -2,11 +2,11 @@ EXE_INC = \
     -Ifluid \
     -Isolid \
     -IregionProperties \
+    -Iinclude \
     -I$(LIB_SRC)/meshTools/lnInclude \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
-    -I$(LIB_SRC)/turbulenceModels/compressible/turbulenceModel \
-    -I$(LIB_SRC)/turbulenceModels/compressible/RAS/lnInclude
+    -I$(LIB_SRC)/turbulenceModels/compressible/turbulenceModel
 
 EXE_LIBS = \
     -lbasicThermophysicalModels \

applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionFoam.C

@@ -37,6 +37,7 @@ Description
 #include "fixedGradientFvPatchFields.H"
 #include "regionProperties.H"
 #include "compressibleCourantNo.H"
+#include "solidRegionDiffNo.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -56,23 +57,24 @@ int main(int argc, char *argv[])
     #include "initContinuityErrs.H"
 
     #include "readTimeControls.H"
+    #include "readSolidTimeControls.H"
+
+
+    #include "compressibleMultiRegionCourantNo.H"
+    #include "solidRegionDiffusionNo.H"
+    #include "setInitialMultiRegionDeltaT.H"
 
-    if (fluidRegions.size())
-    {
-        #include "compressibleMultiRegionCourantNo.H"
-        #include "setInitialDeltaT.H"
-    }
 
     while (runTime.run())
     {
         #include "readTimeControls.H"
+        #include "readSolidTimeControls.H"
         #include "readPIMPLEControls.H"
 
-        if (fluidRegions.size())
-        {
-            #include "compressibleMultiRegionCourantNo.H"
-            #include "setDeltaT.H"
-        }
+
+        #include "compressibleMultiRegionCourantNo.H"
+        #include "solidRegionDiffusionNo.H"
+        #include "setMultiRegionDeltaT.H"
 
         runTime++;
 
@@ -117,7 +119,7 @@ int main(int argc, char *argv[])
             << nl << endl;
     }
 
-    Info<< "End\n" << endl;
+    Info << "End\n" << endl;
 
     return 0;
 }

applications/solvers/heatTransfer/chtMultiRegionFoam/derivedFvPatchFields/solidWallMixedTemperatureCoupled/solidWallMixedTemperatureCoupledFvPatchScalarField.C

@@ -0,0 +1,382 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-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
+
+\*---------------------------------------------------------------------------*/
+
+#include "solidWallMixedTemperatureCoupledFvPatchScalarField.H"
+#include "addToRunTimeSelectionTable.H"
+#include "fvPatchFieldMapper.H"
+#include "volFields.H"
+#include "directMappedPatchBase.H"
+#include "regionProperties.H"
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+bool Foam::solidWallMixedTemperatureCoupledFvPatchScalarField::interfaceOwner
+(
+    const polyMesh& nbrRegion
+) const
+{
+    const fvMesh& myRegion = patch().boundaryMesh().mesh();
+
+    const regionProperties& props =
+        myRegion.objectRegistry::parent().lookupObject<regionProperties>
+        (
+            "regionProperties"
+        );
+
+    label myIndex = findIndex(props.fluidRegionNames(), myRegion.name());
+    if (myIndex == -1)
+    {
+        label i = findIndex(props.solidRegionNames(), myRegion.name());
+
+        if (i == -1)
+        {
+            FatalErrorIn
+            (
+                "solidWallMixedTemperatureCoupledFvPatchScalarField"
+                "::interfaceOwner(const polyMesh&) const"
+            )   << "Cannot find region " << myRegion.name()
+                << " neither in fluids " << props.fluidRegionNames()
+                << " nor in solids " << props.solidRegionNames()
+                << exit(FatalError);
+        }
+        myIndex = props.fluidRegionNames().size() + i;
+    }
+    label nbrIndex = findIndex(props.fluidRegionNames(), nbrRegion.name());
+    if (nbrIndex == -1)
+    {
+        label i = findIndex(props.solidRegionNames(), nbrRegion.name());
+
+        if (i == -1)
+        {
+            FatalErrorIn("coupleManager::interfaceOwner(const polyMesh&) const")
+                << "Cannot find region " << nbrRegion.name()
+                << " neither in fluids " << props.fluidRegionNames()
+                << " nor in solids " << props.solidRegionNames()
+                << exit(FatalError);
+        }
+        nbrIndex = props.fluidRegionNames().size() + i;
+    }
+
+    return myIndex < nbrIndex;
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::solidWallMixedTemperatureCoupledFvPatchScalarField::
+solidWallMixedTemperatureCoupledFvPatchScalarField
+(
+    const fvPatch& p,
+    const DimensionedField<scalar, volMesh>& iF
+)
+:
+    mixedFvPatchScalarField(p, iF),
+    neighbourFieldName_("undefined-neighbourFieldName"),
+    KName_("undefined-K")
+{
+    this->refValue() = 0.0;
+    this->refGrad() = 0.0;
+    this->valueFraction() = 1.0;
+    this->fixesValue_ = true;
+}
+
+
+Foam::solidWallMixedTemperatureCoupledFvPatchScalarField::
+solidWallMixedTemperatureCoupledFvPatchScalarField
+(
+    const solidWallMixedTemperatureCoupledFvPatchScalarField& ptf,
+    const fvPatch& p,
+    const DimensionedField<scalar, volMesh>& iF,
+    const fvPatchFieldMapper& mapper
+)
+:
+    mixedFvPatchScalarField(ptf, p, iF, mapper),
+    neighbourFieldName_(ptf.neighbourFieldName_),
+    KName_(ptf.KName_),
+    fixesValue_(ptf.fixesValue_)
+{}
+
+
+Foam::solidWallMixedTemperatureCoupledFvPatchScalarField::
+solidWallMixedTemperatureCoupledFvPatchScalarField
+(
+    const fvPatch& p,
+    const DimensionedField<scalar, volMesh>& iF,
+    const dictionary& dict
+)
+:
+    mixedFvPatchScalarField(p, iF),
+    neighbourFieldName_(dict.lookup("neighbourFieldName")),
+    KName_(dict.lookup("K"))
+{
+    if (!isA<directMappedPatchBase>(this->patch().patch()))
+    {
+        FatalErrorIn
+        (
+            "solidWallMixedTemperatureCoupledFvPatchScalarField::"
+            "solidWallMixedTemperatureCoupledFvPatchScalarField\n"
+            "(\n"
+            "    const fvPatch& p,\n"
+            "    const DimensionedField<scalar, volMesh>& iF,\n"
+            "    const dictionary& dict\n"
+            ")\n"
+        )   << "\n    patch type '" << p.type()
+            << "' not type '" << directMappedPatchBase::typeName << "'"
+            << "\n    for patch " << p.name()
+            << " of field " << dimensionedInternalField().name()
+            << " in file " << dimensionedInternalField().objectPath()
+            << exit(FatalError);
+    }
+
+    fvPatchScalarField::operator=(scalarField("value", dict, p.size()));
+
+    if (dict.found("refValue"))
+    {
+        // Full restart
+        refValue() = scalarField("refValue", dict, p.size());
+        refGrad() = scalarField("refGradient", dict, p.size());
+        valueFraction() = scalarField("valueFraction", dict, p.size());
+        fixesValue_ = readBool(dict.lookup("fixesValue"));
+    }
+    else
+    {
+        // Start from user entered data. Assume fixedValue.
+        refValue() = *this;
+        refGrad() = 0.0;
+        valueFraction() = 1.0;
+        fixesValue_ = true;
+    }
+}
+
+
+Foam::solidWallMixedTemperatureCoupledFvPatchScalarField::
+solidWallMixedTemperatureCoupledFvPatchScalarField
+(
+    const solidWallMixedTemperatureCoupledFvPatchScalarField& wtcsf,
+    const DimensionedField<scalar, volMesh>& iF
+)
+:
+    mixedFvPatchScalarField(wtcsf, iF),
+    neighbourFieldName_(wtcsf.neighbourFieldName_),
+    KName_(wtcsf.KName_),
+    fixesValue_(wtcsf.fixesValue_)
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+const Foam::fvPatchScalarField&
+Foam::solidWallMixedTemperatureCoupledFvPatchScalarField::K() const
+{
+    return this->patch().lookupPatchField<volScalarField, scalar>(KName_);
+}
+
+
+void Foam::solidWallMixedTemperatureCoupledFvPatchScalarField::updateCoeffs()
+{
+    if (updated())
+    {
+        return;
+    }
+
+    // Get the coupling information from the directMappedPatchBase
+    const directMappedPatchBase& mpp = refCast<const directMappedPatchBase>
+    (
+        patch().patch()
+    );
+    const polyMesh& nbrMesh = mpp.sampleMesh();
+    // Force recalculation of mapping and schedule
+    const mapDistribute& distMap = mpp.map();
+    (void)distMap.schedule();
+
+    tmp<scalarField> intFld = patchInternalField();
+
+    if (interfaceOwner(nbrMesh))
+    {
+        // Note: other side information could be cached - it only needs
+        // to be updated the first time round the iteration (i.e. when
+        // switching regions) but unfortunately we don't have this information.
+
+        const fvPatch& nbrPatch = refCast<const fvMesh>
+        (
+            nbrMesh
+        ).boundary()[mpp.samplePolyPatch().index()];
+
+
+        // Calculate the temperature by harmonic averaging
+        // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+        const solidWallMixedTemperatureCoupledFvPatchScalarField& nbrField =
+        refCast<const solidWallMixedTemperatureCoupledFvPatchScalarField>
+        (
+            nbrPatch.lookupPatchField<volScalarField, scalar>
+            (
+                neighbourFieldName_
+            )
+        );
+
+        // Swap to obtain full local values of neighbour internal field
+        scalarField nbrIntFld = nbrField.patchInternalField();
+        mapDistribute::distribute
+        (
+            Pstream::defaultCommsType,
+            distMap.schedule(),
+            distMap.constructSize(),
+            distMap.subMap(),           // what to send
+            distMap.constructMap(),     // what to receive
+            nbrIntFld
+        );
+
+        // Swap to obtain full local values of neighbour K*delta
+        scalarField nbrKDelta = nbrField.K()*nbrPatch.deltaCoeffs();
+        mapDistribute::distribute
+        (
+            Pstream::defaultCommsType,
+            distMap.schedule(),
+            distMap.constructSize(),
+            distMap.subMap(),           // what to send
+            distMap.constructMap(),     // what to receive
+            nbrKDelta
+        );
+
+
+        tmp<scalarField> myKDelta = K()*patch().deltaCoeffs();
+
+        // Calculate common wall temperature. Reuse *this to store common value.
+        scalarField Twall
+        (
+            (myKDelta()*intFld() + nbrKDelta*nbrIntFld)
+          / (myKDelta() + nbrKDelta)
+        );
+        // Assign to me
+        fvPatchScalarField::operator=(Twall);
+        // Distribute back and assign to neighbour
+        mapDistribute::distribute
+        (
+            Pstream::defaultCommsType,
+            distMap.schedule(),
+            nbrField.size(),
+            distMap.constructMap(),     // reverse : what to send
+            distMap.subMap(),
+            Twall
+        );
+        const_cast<solidWallMixedTemperatureCoupledFvPatchScalarField&>
+        (
+            nbrField
+        ).fvPatchScalarField::operator=(Twall);
+    }
+
+
+    // Switch between fixed value (of harmonic avg) or gradient
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+    label nFixed = 0;
+
+    // Like snGrad but bypass switching on refValue/refGrad.
+    tmp<scalarField> normalGradient = (*this-intFld())*patch().deltaCoeffs();
+
+    if (debug)
+    {
+        scalar Q = gSum(K()*patch().magSf()*normalGradient());
+
+        Info<< "solidWallMixedTemperatureCoupledFvPatchScalarField::"
+            << "updateCoeffs() :"
+            << " patch:" << patch().name()
+            << " heatFlux:" << Q
+            << " walltemperature "
+            << " min:" << gMin(*this)
+            << " max:" << gMax(*this)
+            << " avg:" << gAverage(*this)
+            << endl;
+    }
+
+    forAll(*this, i)
+    {
+        // if outgoing flux use fixed value.
+        if (normalGradient()[i] < 0.0)
+        {
+            this->refValue()[i] = operator[](i);
+            this->refGrad()[i] = 0.0;   // not used by me
+            this->valueFraction()[i] = 1.0;
+            nFixed++;
+        }
+        else
+        {
+            // Fixed gradient. Make sure to have valid refValue (even though
+            // I am not using it - other boundary conditions might)
+            this->refValue()[i] = operator[](i);
+            this->refGrad()[i] = normalGradient()[i];
+            this->valueFraction()[i] = 0.0;
+        }
+    }
+
+    reduce(nFixed, sumOp<label>());
+
+    fixesValue_ = (nFixed > 0);
+
+    if (debug)
+    {
+        label nTotSize = returnReduce(this->size(), sumOp<label>());
+
+        Info<< "solidWallMixedTemperatureCoupledFvPatchScalarField::"
+            << "updateCoeffs() :"
+            << " patch:" << patch().name()
+            << " out of:" << nTotSize
+            << " fixedBC:" << nFixed
+            << " gradient:" << nTotSize-nFixed << endl;
+    }
+
+    mixedFvPatchScalarField::updateCoeffs();
+}
+
+
+void Foam::solidWallMixedTemperatureCoupledFvPatchScalarField::write
+(
+    Ostream& os
+) const
+{
+    mixedFvPatchScalarField::write(os);
+    os.writeKeyword("neighbourFieldName")<< neighbourFieldName_
+        << token::END_STATEMENT << nl;
+    os.writeKeyword("K") << KName_ << token::END_STATEMENT << nl;
+    os.writeKeyword("fixesValue") << fixesValue_ << token::END_STATEMENT << nl;
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+makePatchTypeField
+(
+    fvPatchScalarField,
+    solidWallMixedTemperatureCoupledFvPatchScalarField
+);
+
+} // End namespace Foam
+
+// ************************************************************************* //

applications/solvers/heatTransfer/chtMultiRegionFoam/derivedFvPatchFields/solidWallMixedTemperatureCoupled/solidWallMixedTemperatureCoupledFvPatchScalarField.H

@@ -0,0 +1,193 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-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
+
+Class
+    solidWallMixedTemperatureCoupledFvPatchScalarField
+
+Description
+    Mixed boundary condition for temperature, to be used by the
+    conjugate heat transfer solver.
+    If my temperature is T1, neighbour is T2:
+
+    T1 > T2: my side becomes fixedValue T2 bc, other side becomes fixedGradient.
+    
+
+    Example usage:
+        myInterfacePatchName
+        {
+            type                solidWallMixedTemperatureCoupled;
+            neighbourFieldName  T;
+            K                   K;
+            value               uniform 300;
+        }
+
+    Needs to be on underlying directMapped(Wall)FvPatch.
+
+    Note: runs in parallel with arbitrary decomposition. Uses directMapped
+    functionality to calculate exchange.
+
+    Note: lags interface data so both sides use same data.
+    - problem: schedule to calculate average would interfere
+    with standard processor swaps.
+    - so: updateCoeffs sets both to same Twall. Only need to do
+    this for last outer iteration but don't have access to this.
+
+SourceFiles
+    solidWallMixedTemperatureCoupledFvPatchScalarField.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef solidWallMixedTemperatureCoupledFvPatchScalarField_H
+#define solidWallMixedTemperatureCoupledFvPatchScalarField_H
+
+#include "fvPatchFields.H"
+#include "mixedFvPatchFields.H"
+#include "fvPatch.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+        Class solidWallMixedTemperatureCoupledFvPatchScalarField Declaration
+\*---------------------------------------------------------------------------*/
+
+class solidWallMixedTemperatureCoupledFvPatchScalarField
+:
+    public mixedFvPatchScalarField
+{
+    // Private data
+
+        //- Name of field on the neighbour region
+        const word neighbourFieldName_;
+        
+        //- Name of thermal conductivity field
+        const word KName_;
+
+        bool fixesValue_;
+
+
+    // Private Member Functions
+
+        //- Am I or neighbour owner of interface
+        bool interfaceOwner(const polyMesh& nbrRegion) const;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("solidWallMixedTemperatureCoupled");
+
+
+    // Constructors
+
+        //- Construct from patch and internal field
+        solidWallMixedTemperatureCoupledFvPatchScalarField
+        (
+            const fvPatch&,
+            const DimensionedField<scalar, volMesh>&
+        );
+
+        //- Construct from patch, internal field and dictionary
+        solidWallMixedTemperatureCoupledFvPatchScalarField
+        (
+            const fvPatch&,
+            const DimensionedField<scalar, volMesh>&,
+            const dictionary&
+        );
+
+        //- Construct by mapping given
+        //  solidWallMixedTemperatureCoupledFvPatchScalarField onto a new patch
+        solidWallMixedTemperatureCoupledFvPatchScalarField
+        (
+            const solidWallMixedTemperatureCoupledFvPatchScalarField&,
+            const fvPatch&,
+            const DimensionedField<scalar, volMesh>&,
+            const fvPatchFieldMapper&
+        );
+
+        //- Construct and return a clone
+        virtual tmp<fvPatchScalarField> clone() const
+        {
+            return tmp<fvPatchScalarField>
+            (
+                new solidWallMixedTemperatureCoupledFvPatchScalarField(*this)
+            );
+        }
+
+        //- Construct as copy setting internal field reference
+        solidWallMixedTemperatureCoupledFvPatchScalarField
+        (
+            const solidWallMixedTemperatureCoupledFvPatchScalarField&,
+            const DimensionedField<scalar, volMesh>&
+        );
+
+        //- Construct and return a clone setting internal field reference
+        virtual tmp<fvPatchScalarField> clone
+        (
+            const DimensionedField<scalar, volMesh>& iF
+        ) const
+        {
+            return tmp<fvPatchScalarField>
+            (
+                new solidWallMixedTemperatureCoupledFvPatchScalarField
+                (
+                    *this,
+                    iF
+                )
+            );
+        }
+
+
+    // Member functions
+
+        //- Get corresponding K field
+        const fvPatchScalarField& K() const;
+
+        //- Return true if this patch field fixes a value.
+        //  Needed to check if a level has to be specified while solving
+        //  Poissons equations.
+        virtual bool fixesValue() const
+        {
+            return fixesValue_;
+        }
+
+        //- Update the coefficients associated with the patch field
+        virtual void updateCoeffs();
+
+        //- Write
+        virtual void write(Ostream&) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/compressibleCourantNo.C

@@ -47,10 +47,10 @@ Foam::scalar Foam::compressibleCourantNo
           / fvc::interpolate(rho);
 
         CoNum = max(SfUfbyDelta/mesh.magSf())
-            .value()*runTime.deltaTValue();
+            .value()*runTime.deltaT().value();
 
         meanCoNum = (sum(SfUfbyDelta)/sum(mesh.magSf()))
-            .value()*runTime.deltaTValue();
+            .value()*runTime.deltaT().value();
     }
 
     Info<< "Region: " << mesh.name() << " Courant Number mean: " << meanCoNum

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

@@ -1,15 +1,18 @@
     scalar CoNum = -GREAT;
-    forAll(fluidRegions, regionI)
+    if (fluidRegions.size())
     {
-        CoNum = max
-        (
-            compressibleCourantNo
+        forAll(fluidRegions, regionI)
+        {
+            CoNum = max
             (
-                fluidRegions[regionI],
-                runTime,
-                rhoFluid[regionI],
-                phiFluid[regionI]
-            ),
-            CoNum
-        );
+                compressibleCourantNo
+                (
+                    fluidRegions[regionI],
+                    runTime,
+                    rhoFluid[regionI],
+                    phiFluid[regionI]
+                ),
+                CoNum
+            );
+        }
     }

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

@@ -41,6 +41,24 @@
             )
         );
 
+        Info<< "    Adding to KFluid\n" << endl;
+        KFluid.set
+        (
+            i,
+            new volScalarField
+            (
+                IOobject
+                (
+                    "K",
+                    runTime.timeName(),
+                    fluidRegions[i],
+                    IOobject::NO_READ,
+                    IOobject::NO_WRITE
+                ),
+                thermoFluid[i].Cp()*thermoFluid[i].alpha()
+            )
+        );
+
         Info<< "    Adding to UFluid\n" << endl;
         UFluid.set
         (
@@ -111,24 +129,6 @@
             ).ptr()
         );
 
-        Info<< "    Adding to KFluid\n" << endl;
-        KFluid.set
-        (
-            i,
-            new volScalarField
-            (
-                IOobject
-                (
-                    "K",
-                    runTime.timeName(),
-                    fluidRegions[i],
-                    IOobject::NO_READ,
-                    IOobject::NO_WRITE
-                ),
-                thermoFluid[i].Cp()*turbulence[i].alphaEff()
-            )
-        );
-
         Info<< "    Adding to DpDtFluid\n" << endl;
         DpDtFluid.set
         (

applications/solvers/heatTransfer/chtMultiRegionFoam/include/setInitialMultiRegionDeltaT.H

@@ -0,0 +1,60 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2010-2010 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
+
+Global
+    setInitialDeltaT
+
+Description
+    Set the initial timestep for the CHT MultiRegion solver.
+
+\*---------------------------------------------------------------------------*/
+
+if (adjustTimeStep)
+{
+    if ((runTime.timeIndex() == 0) && ((CoNum > SMALL) || (DiNum > SMALL)))
+    {
+
+        if (CoNum == -GREAT)
+        {
+            CoNum = SMALL;
+        }
+
+        if (DiNum == -GREAT)
+        {
+            DiNum = SMALL;
+        }
+
+        runTime.setDeltaT
+        (
+            min
+            (
+                min(maxCo/CoNum, maxDi/DiNum)*runTime.deltaT().value(),
+                maxDeltaT
+            )
+        );
+        Info<< "deltaT = " <<  runTime.deltaT().value() << endl;
+    }
+}
+
+// ************************************************************************* //

applications/solvers/heatTransfer/chtMultiRegionFoam/include/setMultiRegionDeltaT.H

@@ -0,0 +1,69 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2010-2010 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
+
+Global
+    setMultiRegionDeltaT
+
+Description
+    Reset the timestep to maintain a constant maximum courant and diffusion
+    Numbers. Reduction of time-step is immediate, but increase is damped to avoid
+    unstable oscillations.
+
+\*---------------------------------------------------------------------------*/
+
+if (adjustTimeStep)
+{
+    if (CoNum == -GREAT)
+    {
+        CoNum = SMALL;
+    }
+
+    if (DiNum == -GREAT)
+    {
+        DiNum = SMALL;
+    }
+
+    scalar maxDeltaTFluid = maxCo/(CoNum + SMALL);
+    scalar maxDeltaTSolid = maxDi/(DiNum + SMALL);
+
+    scalar deltaTFluid =
+        min
+        (
+            min(maxDeltaTFluid, 1.0 + 0.1*maxDeltaTFluid),
+            1.2
+        );
+
+    runTime.setDeltaT
+    (
+        min
+        (
+            min(deltaTFluid, maxDeltaTSolid)*runTime.deltaT().value(),
+            maxDeltaT
+        )
+    );
+
+    Info<< "deltaT = " <<  runTime.deltaT().value() << endl;
+}
+
+// ************************************************************************* //

applications/solvers/heatTransfer/chtMultiRegionFoam/regionProperties/regionProperties.C

@@ -0,0 +1,69 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-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
+
+\*---------------------------------------------------------------------------*/
+
+#include "regionProperties.H"
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::regionProperties::regionProperties(const Time& runTime)
+:
+    IOdictionary
+    (
+        IOobject
+        (
+            "regionProperties",
+            runTime.time().constant(),
+            runTime.db(),
+            IOobject::MUST_READ,
+            IOobject::NO_WRITE
+        )
+    ),
+    fluidRegionNames_(lookup("fluidRegionNames")),
+    solidRegionNames_(lookup("solidRegionNames"))
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::regionProperties::~regionProperties()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+const Foam::List<Foam::word>& Foam::regionProperties::fluidRegionNames() const
+{
+    return fluidRegionNames_;
+}
+
+
+const Foam::List<Foam::word>& Foam::regionProperties::solidRegionNames() const
+{
+    return solidRegionNames_;
+}
+
+
+// ************************************************************************* //

applications/solvers/heatTransfer/chtMultiRegionFoam/regionProperties/regionProperties.H

@@ -23,20 +23,21 @@ License
     Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
 
 Class
-    Foam::BSpline
+    regionProperties
 
 Description
-    A cubic spline going through all the knots
+    Simple class to hold region information for coupled region simulations
 
 SourceFiles
-    BSpline.C
+    regionProperties.C
 
 \*---------------------------------------------------------------------------*/
 
-#ifndef BSpline_H
-#define BSpline_H
+#ifndef regionProperties_H
+#define regionProperties_H
 
-#include "spline.H"
+#include "IOdictionary.H"
+#include "Time.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -44,27 +45,29 @@ namespace Foam
 {
 
 /*---------------------------------------------------------------------------*\
-                           Class BSpline Declaration
+                       Class regionProperties Declaration
 \*---------------------------------------------------------------------------*/
 
-class BSpline
+class regionProperties
 :
-    public spline
+    public IOdictionary
 {
+    // Private data
+
+        //- List of the fluid region names
+        List<word> fluidRegionNames_;
+
+        //- List of the solid region names
+        List<word> solidRegionNames_;
+
+
     // Private Member Functions
 
-        pointField findKnots
-        (
-            const pointField&,
-            const vector& fstend,
-            const vector& sndend
-        );
-
         //- Disallow default bitwise copy construct
-        BSpline(const BSpline&);
+        regionProperties(const regionProperties&);
 
         //- Disallow default bitwise assignment
-        void operator=(const BSpline&);
+        void operator=(const regionProperties&);
 
 
 public:
@@ -72,26 +75,23 @@ public:
     // Constructors
 
         //- Construct from components
-        BSpline
-        (
-            const pointField& knots,
-            const vector& fstend = vector::zero,
-            const vector& sndend = vector::zero
-        );
+        regionProperties(const Time& runTime);
+
+
+    // Destructor
+
+        ~regionProperties();
 
 
     // Member Functions
 
-        //- Return the real point position corresponding to the curve parameter
-        //  0 <= lambda <= 1
-        point realPosition(const scalar lambda) const;
+        // Access
 
-        //- Return the point position corresponding to the curve parameter
-        //  0 <= lambda <= 1
-        point position(const scalar lambda) const;
+            //- Return const reference to the list of fluid region names
+            const List<word>& fluidRegionNames() const;
 
-        //- Return the length of the curve
-        scalar length() const;
+            //- Return const reference to the list of solid region names
+            const List<word>& solidRegionNames() const;
 };
 
 

applications/solvers/heatTransfer/chtMultiRegionFoam/solid/readSolidTimeControls.H

@@ -0,0 +1,35 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2010-2010 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
+
+Global
+    readSolidTimeControls
+
+Description
+    Read the control parameters used in the solid
+
+\*---------------------------------------------------------------------------*/
+
+scalar maxDi = runTime.controlDict().lookupOrDefault<scalar>("maxDi", 10.0);
+
+// ************************************************************************* //

applications/solvers/heatTransfer/chtMultiRegionFoam/solid/solidRegionDiffNo.C

@@ -0,0 +1,61 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2010-2010 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
+
+\*---------------------------------------------------------------------------*/
+
+#include "solidRegionDiffNo.H"
+#include "fvc.H"
+
+Foam::scalar Foam::solidRegionDiffNo
+(
+    const fvMesh& mesh,
+    const Time& runTime,
+    const volScalarField& Cprho,
+    const volScalarField& K
+)
+{
+    scalar DiNum = 0.0;
+    scalar meanDiNum = 0.0;
+
+    //- Can have fluid domains with 0 cells so do not test.
+    if (mesh.nInternalFaces())
+    {
+        surfaceScalarField KrhoCpbyDelta =
+            mesh.surfaceInterpolation::deltaCoeffs()
+          * fvc::interpolate(K)
+          / fvc::interpolate(Cprho);
+
+        DiNum = max(KrhoCpbyDelta.internalField())*runTime.deltaT().value();
+
+        meanDiNum = (average(KrhoCpbyDelta)).value()*runTime.deltaT().value();
+    }
+
+    Info<< "Region: " << mesh.name() << " Diffusion Number mean: " << meanDiNum
+        << " max: " << DiNum << endl;
+
+    return DiNum;
+}
+
+
+// ************************************************************************* //

applications/solvers/heatTransfer/chtMultiRegionFoam/solid/solidRegionDiffNo.H

@@ -0,0 +1,49 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2010-2010 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
+
+Description
+    Calculates and outputs the mean and maximum Diffusion Numbers for the solid
+    regions
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef solidRegionDiff_H
+#define solidRegionDiff_H
+
+#include "fvMesh.H"
+
+namespace Foam
+{
+    scalar solidRegionDiffNo
+    (
+        const fvMesh& mesh,
+        const Time& runTime,
+        const volScalarField& Cprho,
+        const volScalarField& K
+    );
+}
+
+#endif
+
+// ************************************************************************* //

applications/solvers/heatTransfer/chtMultiRegionFoam/solid/solidRegionDiffusionNo.H

@@ -0,0 +1,18 @@
+    scalar DiNum = -GREAT;
+    if (solidRegions.size())
+    {
+        forAll(solidRegions, regionI)
+        {
+            DiNum = max
+            (
+                solidRegionDiffNo
+                (
+                    solidRegions[regionI],
+                    runTime,
+                    rhosCps[regionI],
+                    Ks[regionI]
+                ),
+                DiNum
+            );
+        }
+    }

applications/solvers/heatTransfer/chtMultiRegionSimpleFoam/chtMultiRegionSimpleFoam.C

@@ -55,7 +55,7 @@ int main(int argc, char *argv[])
     #include "initContinuityErrs.H"
 
 
-    while (runTime.run())
+    while (runTime.loop())
     {
         Info<< "Time = " << runTime.timeName() << nl << endl;
 
@@ -81,8 +81,6 @@ int main(int argc, char *argv[])
             #include "convergenceCheck.H"
         }
 
-        runTime++;
-
         runTime.write();
 
         Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"

applications/solvers/incompressible/icoFoam/Make/files

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

applications/solvers/incompressible/nonNewtonianIcoFoam/Make/files

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

applications/solvers/incompressible/simpleFoam/Make/files

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

applications/solvers/incompressible/simpleFoam/pEqn.H

@@ -40,4 +40,3 @@
     // Momentum corrector
     U -= fvc::grad(p)/AU;
     U.correctBoundaryConditions();
-

applications/solvers/lagrangian/coalChemistryFoam/chemistry.H

@@ -22,4 +22,6 @@
     {
         kappa = 1.0;
     }
+
+    chemistrySh = kappa*chemistry.Sh()();
 }

applications/solvers/lagrangian/coalChemistryFoam/coalChemistryFoam.C

@@ -26,8 +26,13 @@ Application
     coalChemistryFoam
 
 Description
-    Transient solver for compressible, turbulent flow with coal and
-    limestone parcel injections, and combustion.
+    Transient solver for:
+    - compressible,
+    - turbulent flow,
+    with
+    - coal and limestone parcel injections,
+    - energy source, and
+    - combustion.
 
 \*---------------------------------------------------------------------------*/
 
@@ -38,7 +43,7 @@ Description
 #include "CoalCloud.H"
 #include "psiChemistryModel.H"
 #include "chemistrySolver.H"
-#include "timeActivatedExplicitCellSource.H"
+#include "timeActivatedExplicitSource.H"
 #include "radiationModel.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -79,16 +84,8 @@ int main(int argc, char *argv[])
 
         coalParcels.evolve();
 
-        coalParcels.info();
-
-        Info<< endl;
-
         limestoneParcels.evolve();
 
-        limestoneParcels.info();
-
-        Info<< endl;
-
         #include "chemistry.H"
         #include "rhoEqn.H"
 
@@ -97,16 +94,13 @@ int main(int argc, char *argv[])
         {
             #include "UEqn.H"
             #include "YEqn.H"
-            #include "hEqn.H"
+            #include "hsEqn.H"
 
             // --- PISO loop
             for (int corr=1; corr<=nCorr; corr++)
             {
                 #include "pEqn.H"
             }
-
-            Info<< "T gas min/max   = " << min(T).value() << ", "
-                << max(T).value() << endl;
         }
 
         turbulence->correct();

applications/solvers/lagrangian/coalChemistryFoam/createFields.H

@@ -6,7 +6,7 @@
     );
     psiChemistryModel& chemistry = pChemistry();
 
-    hCombustionThermo& thermo = chemistry.thermo();
+    hsCombustionThermo& thermo = chemistry.thermo();
 
     basicMultiComponentMixture& composition = thermo.composition();
     PtrList<volScalarField>& Y = composition.Y();
@@ -22,17 +22,17 @@
     }
 
     volScalarField& p = thermo.p();
-    volScalarField& h = thermo.h();
+    volScalarField& hs = thermo.hs();
     const volScalarField& T = thermo.T();
     const volScalarField& psi = thermo.psi();
 
     multivariateSurfaceInterpolationScheme<scalar>::fieldTable fields;
 
-    forAll (Y, i)
+    forAll(Y, i)
     {
         fields.add(Y[i]);
     }
-    fields.add(h);
+    fields.add(hs);
 
     volScalarField rho
     (
@@ -127,10 +127,25 @@
         fvc::DDt(surfaceScalarField("phiU", phi/fvc::interpolate(rho)), p)
     );
 
-    Info<< "\nConstructing explicit enthalpy cell source" << endl;
-    timeActivatedExplicitCellSource enthalpySource
+    Info<< "\nConstructing explicit enthalpy source" << endl;
+    scalarTimeActivatedExplicitSourceList enthalpySource
     (
-        "enthalpySource",
+        "energy",
         mesh,
-        h.dimensions()*phi.dimensions()/mesh.V().dimensions()
+        dimEnergy/dimTime/dimVolume,
+        "hs"
+    );
+
+    DimensionedField<scalar, volMesh> chemistrySh
+    (
+        IOobject
+        (
+            "chemistry::Sh",
+            runTime.timeName(),
+            mesh,
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        ),
+        mesh,
+        dimensionedScalar("chemistrySh", dimEnergy/dimTime/dimVolume, 0.0)
     );

applications/solvers/lagrangian/coalChemistryFoam/hEqn.H

@@ -1,22 +0,0 @@
-{
-    fvScalarMatrix hEqn
-    (
-        fvm::ddt(rho, h)
-      + mvConvection->fvmDiv(phi, h)
-      - fvm::laplacian(turbulence->alphaEff(), h)
-     ==
-        DpDt
-     +  coalParcels.Sh()
-     +  limestoneParcels.Sh()
-     +  enthalpySource.Su()
-     +  radiation->Sh(thermo)
-    );
-
-    hEqn.relax();
-
-    hEqn.solve();
-
-    thermo.correct();
-
-    radiation->correct();
-}

applications/solvers/lagrangian/coalChemistryFoam/hsEqn.H

@@ -0,0 +1,26 @@
+{
+    fvScalarMatrix hsEqn
+    (
+        fvm::ddt(rho, hs)
+      + mvConvection->fvmDiv(phi, hs)
+      - fvm::laplacian(turbulence->alphaEff(), hs)
+     ==
+        DpDt
+     +  coalParcels.Sh()
+     +  limestoneParcels.Sh()
+     +  enthalpySource.Su()
+     +  radiation->Shs(thermo)
+     +  chemistrySh
+    );
+
+    hsEqn.relax();
+
+    hsEqn.solve();
+
+    thermo.correct();
+
+    radiation->correct();
+
+    Info<< "T gas min/max   = " << min(T).value() << ", "
+        << max(T).value() << endl;
+}

applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/UEqn.H

@@ -7,6 +7,7 @@
      ==
         rho.dimensionedInternalField()*g
       + parcels.SU()
+      + momentumSource.Su()
     );
 
     pZones.addResistance(UEqn);

applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/YEqn.H

@@ -28,7 +28,7 @@ tmp<fv::convectionScheme<scalar> > mvConvection
               ==
                 parcels.Srho(i)
               + kappa*chemistry.RR(i)().dimensionedInternalField()
-              + pointMassSources.Su(i),
+              + massSource.Su(i),
                 mesh.solver("Yi")
             );
 

applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/chemistry.H

@@ -22,4 +22,6 @@
     {
         kappa = 1.0;
     }
+
+    chemistrySh = kappa*chemistry.Sh()();
 }

applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/createClouds.H

@@ -1,5 +1,5 @@
 Info<< "\nConstructing reacting cloud" << endl;
-icoPoly8ThermoReactingCloud parcels
+icoPoly8ThermoReactingMultiphaseCloud parcels
 (
     "reactingCloud1",
     rho,

applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/createExplicitSources.H

@@ -0,0 +1,27 @@
+Info<< "Creating mass source\n" << endl;
+scalarTimeActivatedExplicitSourceList massSource
+(
+    "mass",
+    mesh,
+    dimMass/dimTime/dimVolume,
+    composition.species()
+);
+
+
+Info<< "Creating momentum source\n" << endl;
+vectorTimeActivatedExplicitSourceList momentumSource
+(
+    "momentum",
+    mesh,
+    dimMass*dimVelocity/dimTime/dimVolume,
+    "U"
+);
+
+Info<< "Creating energy source\n" << endl;
+scalarTimeActivatedExplicitSourceList energySource
+(
+    "energy",
+    mesh,
+    dimEnergy/dimTime/dimVolume,
+    "h"
+);

applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/createFields.H

@@ -6,7 +6,7 @@
     );
     rhoChemistryModel& chemistry = pChemistry();
 
-    hReactionThermo& thermo = chemistry.thermo();
+    hsReactionThermo& thermo = chemistry.thermo();
 
     basicMultiComponentMixture& composition = thermo.composition();
     PtrList<volScalarField>& Y = composition.Y();
@@ -22,7 +22,7 @@
     }
 
     volScalarField& p = thermo.p();
-    volScalarField& h = thermo.h();
+    volScalarField& hs = thermo.hs();
     const volScalarField& T = thermo.T();
     const volScalarField& psi = thermo.psi();
 
@@ -84,8 +84,22 @@
     Info<< "Creating multi-variate interpolation scheme\n" << endl;
     multivariateSurfaceInterpolationScheme<scalar>::fieldTable fields;
 
-    forAll (Y, i)
+    forAll(Y, i)
     {
         fields.add(Y[i]);
     }
-    fields.add(h);
+    fields.add(hs);
+
+    DimensionedField<scalar, volMesh> chemistrySh
+    (
+        IOobject
+        (
+            "chemistry::Sh",
+            runTime.timeName(),
+            mesh,
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        ),
+        mesh,
+        dimensionedScalar("chemistrySh", dimEnergy/dimTime/dimVolume, 0.0)
+    );

applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/createMulticomponentPointSources.H

@@ -1,8 +0,0 @@
-Info<< "\nConstructing multi-component mass flow rate point sources" << endl;
-timeActivatedExplicitMulticomponentPointSource pointMassSources
-(
-    "pointMassSources",
-    mesh,
-    Y,
-    dimMass/dimVolume/dimTime
-);

applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/hsEqn.H

@@ -32,13 +32,15 @@
     {
         solve
         (
-            fvm::ddt(rho, h)
-          + mvConvection->fvmDiv(phi, h)
-          - fvm::laplacian(turbulence->alphaEff(), h)
+            fvm::ddt(rho, hs)
+          + mvConvection->fvmDiv(phi, hs)
+          - fvm::laplacian(turbulence->alphaEff(), hs)
          ==
             pWork()
           + parcels.Sh()
-          + radiation->Sh(thermo)
+          + radiation->Shs(thermo)
+          + energySource.Su()
+          + chemistrySh
         );
 
         thermo.correct();

applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/pEqn.H

@@ -32,7 +32,7 @@
           - fvm::laplacian(rho*rAU, p)
         ==
             parcels.Srho()
-          + pointMassSources.Su()
+          + massSource.SuTot()
         );
 
         if (corr == nCorr-1 && nonOrth == nNonOrthCorr)

applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/porousExplicitSourceReactingParcelFoam.C

@@ -27,12 +27,13 @@ Application
 
 Description
     Transient PISO solver for compressible, laminar or turbulent flow with
-    reacting Lagrangian parcels for porous media, including explicit sources
+    reacting multiphase Lagrangian parcels for porous media, including explicit
+    sources for mass, momentum and energy
 
     The solver includes:
-    - reacting parcel cloud
+    - reacting multiphase parcel cloud
     - porous media
-    - point mass sources
+    - mass, momentum and energy sources
     - polynomial based, incompressible thermodynamics (f(T))
 
     Note: ddtPhiCorr not used here when porous zones are active
@@ -43,12 +44,12 @@ Description
 #include "fvCFD.H"
 #include "hReactionThermo.H"
 #include "turbulenceModel.H"
-#include "BasicReactingCloud.H"
+#include "BasicReactingMultiphaseCloud.H"
 #include "rhoChemistryModel.H"
 #include "chemistrySolver.H"
 #include "radiationModel.H"
 #include "porousZones.H"
-#include "timeActivatedExplicitMulticomponentPointSource.H"
+#include "timeActivatedExplicitSource.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -63,7 +64,7 @@ int main(int argc, char *argv[])
     #include "createFields.H"
     #include "createRadiationModel.H"
     #include "createClouds.H"
-    #include "createMulticomponentPointSources.H"
+    #include "createExplicitSources.H"
     #include "createPorousZones.H"
     #include "initContinuityErrs.H"
     #include "readTimeControls.H"
@@ -88,13 +89,11 @@ int main(int argc, char *argv[])
 
         parcels.evolve();
 
-        parcels.info();
-
         #include "chemistry.H"
         #include "rhoEqn.H"
         #include "UEqn.H"
         #include "YEqn.H"
-        #include "hEqn.H"
+        #include "hsEqn.H"
 
         // --- PISO loop
         for (int corr=0; corr<nCorr; corr++)

applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/rhoEqn.H

@@ -37,7 +37,7 @@ Description
       + fvc::div(phi)
       ==
         parcels.Srho()
-      + pointMassSources.Su()
+      + massSource.SuTot()
     );
 
     rhoEqn.relax();

applications/solvers/lagrangian/reactingParcelFoam/chemistry.H

@@ -22,4 +22,6 @@
     {
         kappa = 1.0;
     }
+
+    chemistrySh = kappa*chemistry.Sh()();
 }

applications/solvers/lagrangian/reactingParcelFoam/createFields.H

@@ -6,7 +6,7 @@
     );
     psiChemistryModel& chemistry = pChemistry();
 
-    hCombustionThermo& thermo = chemistry.thermo();
+    hsCombustionThermo& thermo = chemistry.thermo();
 
     basicMultiComponentMixture& composition = thermo.composition();
     PtrList<volScalarField>& Y = composition.Y();
@@ -22,7 +22,7 @@
     }
 
     volScalarField& p = thermo.p();
-    volScalarField& h = thermo.h();
+    volScalarField& hs = thermo.hs();
     const volScalarField& T = thermo.T();
     const volScalarField& psi = thermo.psi();
 
@@ -90,8 +90,22 @@
 
     multivariateSurfaceInterpolationScheme<scalar>::fieldTable fields;
 
-    forAll (Y, i)
+    forAll(Y, i)
     {
         fields.add(Y[i]);
     }
-    fields.add(h);
+    fields.add(hs);
+
+    DimensionedField<scalar, volMesh> chemistrySh
+    (
+        IOobject
+        (
+            "chemistry::Sh",
+            runTime.timeName(),
+            mesh,
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        ),
+        mesh,
+        dimensionedScalar("chemistry::Sh", dimEnergy/dimTime/dimVolume, 0.0)
+    );

applications/solvers/lagrangian/reactingParcelFoam/hEqn.H

@@ -1,20 +0,0 @@
-{
-    fvScalarMatrix hEqn
-    (
-        fvm::ddt(rho, h)
-      + mvConvection->fvmDiv(phi, h)
-      - fvm::laplacian(turbulence->alphaEff(), h)
-     ==
-        DpDt
-     +  parcels.Sh()
-     +  radiation->Sh(thermo)
-    );
-
-    hEqn.relax();
-
-    hEqn.solve();
-
-    thermo.correct();
-
-    radiation->correct();
-}

applications/solvers/lagrangian/reactingParcelFoam/hsEqn.H

@@ -0,0 +1,24 @@
+{
+    fvScalarMatrix hEqn
+    (
+        fvm::ddt(rho, hs)
+      + mvConvection->fvmDiv(phi, hs)
+      - fvm::laplacian(turbulence->alphaEff(), hs)
+     ==
+        DpDt
+     +  parcels.Sh()
+     +  radiation->Shs(thermo)
+     +  chemistrySh
+    );
+
+    hEqn.relax();
+
+    hEqn.solve();
+
+    thermo.correct();
+
+    radiation->correct();
+
+    Info<< "T gas min/max   = " << min(T).value() << ", "
+        << max(T).value() << endl;
+}

applications/solvers/lagrangian/reactingParcelFoam/reactingParcelFoam.C

@@ -74,8 +74,6 @@ int main(int argc, char *argv[])
 
         parcels.evolve();
 
-        parcels.info();
-
         #include "chemistry.H"
         #include "rhoEqn.H"
 
@@ -88,12 +86,9 @@ int main(int argc, char *argv[])
             // --- PISO loop
             for (int corr=1; corr<=nCorr; corr++)
             {
-                #include "hEqn.H"
+                #include "hsEqn.H"
                 #include "pEqn.H"
             }
-
-            Info<< "T gas min/max   = " << min(T).value() << ", "
-                << max(T).value() << endl;
         }
 
         turbulence->correct();

applications/solvers/lagrangian/uncoupledKinematicParcelFoam/createFields.H

@@ -47,8 +47,10 @@
         )
     );
 
-    word kinematicCloudName("kinematicCloud");
-    args.optionReadIfPresent("cloudName", kinematicCloudName);
+    const word kinematicCloudName
+    (
+        args.optionLookupOrDefault<word>("cloudName", "kinematicCloud")
+    );
 
     Info<< "Constructing kinematicCloud " << kinematicCloudName << endl;
     basicKinematicCloud kinematicCloud

applications/solvers/lagrangian/uncoupledKinematicParcelFoam/uncoupledKinematicParcelFoam.C

@@ -42,7 +42,7 @@ Description
 
 int main(int argc, char *argv[])
 {
-    argList::validOptions.insert("cloudName", "cloud name");
+    argList::addOption("cloudName", "cloud name");
 
     #include "setRootCase.H"
     #include "createTime.H"
@@ -61,7 +61,6 @@ int main(int argc, char *argv[])
 
         Info<< "Evolving " << kinematicCloud.name() << endl;
         kinematicCloud.evolve();
-        kinematicCloud.info();
 
         runTime.write();
 

applications/solvers/multiphase/bubbleFoam/UEqns.H

@@ -41,7 +41,7 @@ fvVectorMatrix UbEqn(Ub, Ub.dimensions()*dimVol/dimTime);
        - fvc::interpolate(nuEffb)
         *mesh.magSf()*fvc::snGrad(beta)/fvc::interpolate(beta + scalar(0.001));
 
-    UbEqn = 
+    UbEqn =
     (
         (scalar(1) + Cvm*rhob*alpha/rhob)*
         (

applications/solvers/multiphase/cavitatingFoam/CourantNo.H

@@ -36,7 +36,7 @@ scalar acousticCoNum = 0.0;
 
 if (mesh.nInternalFaces())
 {
-    surfaceScalarField SfUfbyDelta = 
+    surfaceScalarField SfUfbyDelta =
         mesh.surfaceInterpolation::deltaCoeffs()*mag(phiv);
 
     CoNum = max(SfUfbyDelta/mesh.magSf())