Merge branch 'master' of /home/noisy3/OpenFOAM/OpenFOAM-dev

applications/utilities/parallelProcessing/reconstructPar/reconstructPar.C

@@ -179,23 +179,8 @@ int main(int argc, char *argv[])
         // Get list of objects from processor0 database
         IOobjectList objects(procMeshes.meshes()[0], databases[0].timeName());
 
-
-        // If there are any FV fields, reconstruct them
-
-        if
-        (
-            objects.lookupClass(volScalarField::typeName).size()
-         || objects.lookupClass(volVectorField::typeName).size()
-         || objects.lookupClass(volSphericalTensorField::typeName).size()
-         || objects.lookupClass(volSymmTensorField::typeName).size()
-         || objects.lookupClass(volTensorField::typeName).size()
-         || objects.lookupClass(surfaceScalarField::typeName).size()
-         || objects.lookupClass(surfaceVectorField::typeName).size()
-         || objects.lookupClass(surfaceSphericalTensorField::typeName).size()
-         || objects.lookupClass(surfaceSymmTensorField::typeName).size()
-         || objects.lookupClass(surfaceTensorField::typeName).size()
-        )
         {
+            // If there are any FV fields, reconstruct them
             Info<< "Reconstructing FV fields" << nl << endl;
 
             fvFieldReconstructor fvReconstructor
@@ -207,6 +192,32 @@ int main(int argc, char *argv[])
                 procMeshes.boundaryProcAddressing()
             );
 
+            fvReconstructor.reconstructFvVolumeInternalFields<scalar>
+            (
+                objects,
+                selectedFields
+            );
+            fvReconstructor.reconstructFvVolumeInternalFields<vector>
+            (
+                objects,
+                selectedFields
+            );
+            fvReconstructor.reconstructFvVolumeInternalFields<sphericalTensor>
+            (
+                objects,
+                selectedFields
+            );
+            fvReconstructor.reconstructFvVolumeInternalFields<symmTensor>
+            (
+                objects,
+                selectedFields
+            );
+            fvReconstructor.reconstructFvVolumeInternalFields<tensor>
+            (
+                objects,
+                selectedFields
+            );
+
             fvReconstructor.reconstructFvVolumeFields<scalar>
             (
                 objects,
@@ -258,22 +269,13 @@ int main(int argc, char *argv[])
                 objects,
                 selectedFields
             );
-        }
+
-        else
+            if (fvReconstructor.nReconstructed() == 0)
-        {
+            {
-            Info<< "No FV fields" << nl << endl;
+                Info<< "No FV fields" << nl << endl;
+            }
         }
 
-
-        // If there are any point fields, reconstruct them
-        if
-        (
-            objects.lookupClass(pointScalarField::typeName).size()
-         || objects.lookupClass(pointVectorField::typeName).size()
-         || objects.lookupClass(pointSphericalTensorField::typeName).size()
-         || objects.lookupClass(pointSymmTensorField::typeName).size()
-         || objects.lookupClass(pointTensorField::typeName).size()
-        )
         {
             Info<< "Reconstructing point fields" << nl << endl;
 
@@ -318,10 +320,11 @@ int main(int argc, char *argv[])
                 objects,
                 selectedFields
             );
-        }
+
-        else
+            if (pointReconstructor.nReconstructed() == 0)
-        {
+            {
-            Info<< "No point fields" << nl << endl;
+                Info<< "No point fields" << nl << endl;
+            }
         }
 
 

src/parallel/reconstruct/reconstruct/fvFieldReconstructor.C

@@ -40,7 +40,8 @@ Foam::fvFieldReconstructor::fvFieldReconstructor
     procMeshes_(procMeshes),
     faceProcAddressing_(faceProcAddressing),
     cellProcAddressing_(cellProcAddressing),
-    boundaryProcAddressing_(boundaryProcAddressing)
+    boundaryProcAddressing_(boundaryProcAddressing),
+    nReconstructed_(0)
 {}
 
 

src/parallel/reconstruct/reconstruct/fvFieldReconstructor.H

@@ -71,6 +71,9 @@ class fvFieldReconstructor
         //- List of processor boundary addressing lists
         const PtrList<labelIOList>& boundaryProcAddressing_;
 
+        //- Number of fields reconstructed
+        label nReconstructed_;
+
 
     // Private Member Functions
 
@@ -134,20 +137,33 @@ public:
 
     // Member Functions
 
+        //- Return number of fields reconstructed
+        label nReconstructed() const
+        {
+            return nReconstructed_;
+        }
+
+        //- Reconstruct volume internal field
+        template<class Type>
+        tmp<DimensionedField<Type, volMesh> >
+        reconstructFvVolumeInternalField(const IOobject& fieldIoObject);
+
         //- Reconstruct volume field
         template<class Type>
         tmp<GeometricField<Type, fvPatchField, volMesh> >
-        reconstructFvVolumeField
+        reconstructFvVolumeField(const IOobject& fieldIoObject);
-        (
-            const IOobject& fieldIoObject
-        );
 
         //- Reconstruct surface field
         template<class Type>
         tmp<GeometricField<Type, fvsPatchField, surfaceMesh> >
-        reconstructFvSurfaceField
+        reconstructFvSurfaceField(const IOobject& fieldIoObject);
+
+        //- Reconstruct and write all/selected volume internal fields
+        template<class Type>
+        void reconstructFvVolumeInternalFields
         (
-            const IOobject& fieldIoObject
+            const IOobjectList& objects,
+            const HashSet<word>& selectedFields
         );
 
         //- Reconstruct and write all/selected volume fields
@@ -158,7 +174,7 @@ public:
             const HashSet<word>& selectedFields
         );
 
-        //- Reconstruct and write all/selected volume fields
+        //- Reconstruct and write all/selected surface fields
         template<class Type>
         void reconstructFvSurfaceFields
         (

src/parallel/reconstruct/reconstruct/fvFieldReconstructorReconstructFields.C

@@ -33,6 +33,75 @@ License
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
+template<class Type>
+Foam::tmp<Foam::DimensionedField<Type, Foam::volMesh> >
+Foam::fvFieldReconstructor::reconstructFvVolumeInternalField
+(
+    const IOobject& fieldIoObject
+)
+{
+    // Read the field for all the processors
+    PtrList<DimensionedField<Type, volMesh> > procFields
+    (
+        procMeshes_.size()
+    );
+
+    forAll(procMeshes_, procI)
+    {
+        procFields.set
+        (
+            procI,
+            new DimensionedField<Type, volMesh>
+            (
+                IOobject
+                (
+                    fieldIoObject.name(),
+                    procMeshes_[procI].time().timeName(),
+                    procMeshes_[procI],
+                    IOobject::MUST_READ,
+                    IOobject::NO_WRITE
+                ),
+                procMeshes_[procI]
+            )
+        );
+    }
+
+
+    // Create the internalField
+    Field<Type> internalField(mesh_.nCells());
+
+    forAll(procMeshes_, procI)
+    {
+        const DimensionedField<Type, volMesh>& procField = procFields[procI];
+
+        // Set the cell values in the reconstructed field
+        internalField.rmap
+        (
+            procField.field(),
+            cellProcAddressing_[procI]
+        );
+    }
+
+    return tmp<DimensionedField<Type, volMesh> >
+    (
+        new DimensionedField<Type, volMesh>
+        (
+            IOobject
+            (
+                fieldIoObject.name(),
+                mesh_.time().timeName(),
+                mesh_,
+                IOobject::NO_READ,
+                IOobject::NO_WRITE
+            ),
+            mesh_,
+            procFields[0].dimensions(),
+            internalField
+        )
+    );
+}
+
+
 template<class Type>
 Foam::tmp<Foam::GeometricField<Type, Foam::fvPatchField, Foam::volMesh> >
 Foam::fvFieldReconstructor::reconstructFvVolumeField
@@ -451,7 +520,41 @@ Foam::fvFieldReconstructor::reconstructFvSurfaceField
 }
 
 
-// Reconstruct and write all/selected volume fields
+template<class Type>
+void Foam::fvFieldReconstructor::reconstructFvVolumeInternalFields
+(
+    const IOobjectList& objects,
+    const HashSet<word>& selectedFields
+)
+{
+    const word& fieldClassName = DimensionedField<Type, volMesh>::typeName;
+
+    IOobjectList fields = objects.lookupClass(fieldClassName);
+
+    if (fields.size())
+    {
+        Info<< "    Reconstructing " << fieldClassName << "s\n" << endl;
+
+        forAllConstIter(IOobjectList, fields, fieldIter)
+        {
+            if
+            (
+                !selectedFields.size()
+             || selectedFields.found(fieldIter()->name())
+            )
+            {
+                Info<< "        " << fieldIter()->name() << endl;
+
+                reconstructFvVolumeInternalField<Type>(*fieldIter())().write();
+
+                nReconstructed_++;
+            }
+        }
+        Info<< endl;
+    }
+}
+
+
 template<class Type>
 void Foam::fvFieldReconstructor::reconstructFvVolumeFields
 (
@@ -479,13 +582,15 @@ void Foam::fvFieldReconstructor::reconstructFvVolumeFields
                 Info<< "        " << fieldIter()->name() << endl;
 
                 reconstructFvVolumeField<Type>(*fieldIter())().write();
+
+                nReconstructed_++;
             }
         }
         Info<< endl;
     }
 }
 
-// Reconstruct and write all/selected surface fields
+
 template<class Type>
 void Foam::fvFieldReconstructor::reconstructFvSurfaceFields
 (
@@ -513,6 +618,8 @@ void Foam::fvFieldReconstructor::reconstructFvSurfaceFields
                 Info<< "        " << fieldIter()->name() << endl;
 
                 reconstructFvSurfaceField<Type>(*fieldIter())().write();
+
+                nReconstructed_++;
             }
         }
         Info<< endl;

src/parallel/reconstruct/reconstruct/pointFieldReconstructor.C

@@ -39,7 +39,8 @@ Foam::pointFieldReconstructor::pointFieldReconstructor
     procMeshes_(procMeshes),
     pointProcAddressing_(pointProcAddressing),
     boundaryProcAddressing_(boundaryProcAddressing),
-    patchPointAddressing_(procMeshes.size())
+    patchPointAddressing_(procMeshes.size()),
+    nReconstructed_(0)
 {
     // Inverse-addressing of the patch point labels.
     labelList pointMap(mesh_.size(), -1);

src/parallel/reconstruct/reconstruct/pointFieldReconstructor.H

@@ -68,6 +68,9 @@ class pointFieldReconstructor
         //- Point patch addressing
         labelListListList patchPointAddressing_;
 
+        //- Number of fields reconstructed
+        label nReconstructed_;
+
 
     // Private Member Functions
 
@@ -130,6 +133,12 @@ public:
 
     // Member Functions
 
+        //- Return number of fields reconstructed
+        label nReconstructed() const
+        {
+            return nReconstructed_;
+        }
+
         //- Reconstruct field
         template<class Type>
         tmp<GeometricField<Type, pointPatchField, pointMesh> >

src/thermophysicalModels/chemistryModel/chemistryModel/ODEChemistryModel/ODEChemistryModel.C

@@ -93,7 +93,8 @@ Foam::ODEChemistryModel<CompType, ThermoType>::~ODEChemistryModel()
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
 template<class CompType, class ThermoType>
-Foam::scalarField Foam::ODEChemistryModel<CompType, ThermoType>::omega
+Foam::tmp<Foam::scalarField>
+Foam::ODEChemistryModel<CompType, ThermoType>::omega
 (
     const scalarField& c,
     const scalar T,
@@ -103,7 +104,8 @@ Foam::scalarField Foam::ODEChemistryModel<CompType, ThermoType>::omega
     scalar pf, cf, pr, cr;
     label lRef, rRef;
 
-    scalarField om(nEqns(), 0.0);
+    tmp<scalarField> tom(new scalarField(nEqns(), 0.0));
+    scalarField& om = tom();
 
     forAll(reactions_, i)
     {
@@ -116,20 +118,20 @@ Foam::scalarField Foam::ODEChemistryModel<CompType, ThermoType>::omega
 
         forAll(R.lhs(), s)
         {
-            label si = R.lhs()[s].index;
+            const label si = R.lhs()[s].index;
-            scalar sl = R.lhs()[s].stoichCoeff;
+            const scalar sl = R.lhs()[s].stoichCoeff;
             om[si] -= sl*omegai;
         }
 
         forAll(R.rhs(), s)
         {
-            label si = R.rhs()[s].index;
+            const label si = R.rhs()[s].index;
-            scalar sr = R.rhs()[s].stoichCoeff;
+            const scalar sr = R.rhs()[s].stoichCoeff;
             om[si] += sr*omegai;
         }
     }
 
-    return om;
+    return tom;
 }
 
 
@@ -149,13 +151,13 @@ Foam::scalar Foam::ODEChemistryModel<CompType, ThermoType>::omega
 ) const
 {
     scalarField c2(nSpecie_, 0.0);
-    for (label i=0; i<nSpecie_; i++)
+    for (label i = 0; i < nSpecie_; i++)
     {
         c2[i] = max(0.0, c[i]);
     }
 
-    scalar kf = R.kf(T, p, c2);
+    const scalar kf = R.kf(T, p, c2);
-    scalar kr = R.kr(kf, T, p, c2);
+    const scalar kr = R.kr(kf, T, p, c2);
 
     pf = 1.0;
     pr = 1.0;
@@ -167,28 +169,28 @@ Foam::scalar Foam::ODEChemistryModel<CompType, ThermoType>::omega
     lRef = R.lhs()[slRef].index;
 
     pf = kf;
-    for (label s=1; s<Nl; s++)
+    for (label s = 1; s < Nl; s++)
     {
-        label si = R.lhs()[s].index;
+        const label si = R.lhs()[s].index;
 
         if (c[si] < c[lRef])
         {
-            scalar exp = R.lhs()[slRef].exponent;
+            const scalar exp = R.lhs()[slRef].exponent;
             pf *= pow(max(0.0, c[lRef]), exp);
             lRef = si;
             slRef = s;
         }
         else
         {
-            scalar exp = R.lhs()[s].exponent;
+            const scalar exp = R.lhs()[s].exponent;
             pf *= pow(max(0.0, c[si]), exp);
         }
     }
     cf = max(0.0, c[lRef]);
 
     {
-        scalar exp = R.lhs()[slRef].exponent;
+        const scalar exp = R.lhs()[slRef].exponent;
-        if (exp<1.0)
+        if (exp < 1.0)
         {
             if (cf > SMALL)
             {
@@ -210,27 +212,27 @@ Foam::scalar Foam::ODEChemistryModel<CompType, ThermoType>::omega
 
     // find the matrix element and element position for the rhs
     pr = kr;
-    for (label s=1; s<Nr; s++)
+    for (label s = 1; s < Nr; s++)
     {
-        label si = R.rhs()[s].index;
+        const label si = R.rhs()[s].index;
         if (c[si] < c[rRef])
         {
-            scalar exp = R.rhs()[srRef].exponent;
+            const scalar exp = R.rhs()[srRef].exponent;
             pr *= pow(max(0.0, c[rRef]), exp);
             rRef = si;
             srRef = s;
         }
         else
         {
-            scalar exp = R.rhs()[s].exponent;
+            const scalar exp = R.rhs()[s].exponent;
             pr *= pow(max(0.0, c[si]), exp);
         }
     }
     cr = max(0.0, c[rRef]);
 
     {
-        scalar exp = R.rhs()[srRef].exponent;
+        const scalar exp = R.rhs()[srRef].exponent;
-        if (exp<1.0)
+        if (exp < 1.0)
         {
             if (cr>SMALL)
             {
@@ -259,8 +261,8 @@ void Foam::ODEChemistryModel<CompType, ThermoType>::derivatives
     scalarField& dcdt
 ) const
 {
-    scalar T = c[nSpecie_];
+    const scalar T = c[nSpecie_];
-    scalar p = c[nSpecie_ + 1];
+    const scalar p = c[nSpecie_ + 1];
 
     dcdt = omega(c, T, p);
 
@@ -268,37 +270,37 @@ void Foam::ODEChemistryModel<CompType, ThermoType>::derivatives
     // dT/dt = ...
     scalar rho = 0.0;
     scalar cSum = 0.0;
-    for (label i=0; i<nSpecie_; i++)
+    for (label i = 0; i < nSpecie_; i++)
     {
-        scalar W = specieThermo_[i].W();
+        const scalar W = specieThermo_[i].W();
         cSum += c[i];
         rho += W*c[i];
     }
-    scalar mw = rho/cSum;
+    const scalar mw = rho/cSum;
     scalar cp = 0.0;
     for (label i=0; i<nSpecie_; i++)
     {
-        scalar cpi = specieThermo_[i].cp(T);
+        const scalar cpi = specieThermo_[i].cp(T);
-        scalar Xi = c[i]/rho;
+        const scalar Xi = c[i]/rho;
         cp += Xi*cpi;
     }
     cp /= mw;
 
     scalar dT = 0.0;
-    for (label i=0; i<nSpecie_; i++)
+    for (label i = 0; i < nSpecie_; i++)
     {
-        scalar hi = specieThermo_[i].h(T);
+        const scalar hi = specieThermo_[i].h(T);
         dT += hi*dcdt[i];
     }
     dT /= rho*cp;
 
     // limit the time-derivative, this is more stable for the ODE
     // solver when calculating the allowed time step
-    scalar dtMag = min(500.0, mag(dT));
+    const scalar dTLimited = min(500.0, mag(dT));
-    dcdt[nSpecie_] = -dT*dtMag/(mag(dT) + 1.0e-10);
+    dcdt[nSpecie_] = -dT*dTLimited/(mag(dT) + 1.0e-10);
 
     // dp/dt = ...
-    dcdt[nSpecie_+1] = 0.0;
+    dcdt[nSpecie_ + 1] = 0.0;
 }
 
 
@@ -311,11 +313,11 @@ void Foam::ODEChemistryModel<CompType, ThermoType>::jacobian
     scalarSquareMatrix& dfdc
 ) const
 {
-    scalar T = c[nSpecie_];
+    const scalar T = c[nSpecie_];
-    scalar p = c[nSpecie_ + 1];
+    const scalar p = c[nSpecie_ + 1];
 
     scalarField c2(nSpecie_, 0.0);
-    for (label i=0; i<nSpecie_; i++)
+    forAll(c2, i)
     {
         c2[i] = max(c[i], 0.0);
     }
@@ -335,22 +337,22 @@ void Foam::ODEChemistryModel<CompType, ThermoType>::jacobian
     {
         const Reaction<ThermoType>& R = reactions_[ri];
 
-        scalar kf0 = R.kf(T, p, c2);
+        const scalar kf0 = R.kf(T, p, c2);
-        scalar kr0 = R.kr(T, p, c2);
+        const scalar kr0 = R.kr(T, p, c2);
 
         forAll(R.lhs(), j)
         {
-            label sj = R.lhs()[j].index;
+            const label sj = R.lhs()[j].index;
             scalar kf = kf0;
             forAll(R.lhs(), i)
             {
-                label si = R.lhs()[i].index;
+                const label si = R.lhs()[i].index;
-                scalar el = R.lhs()[i].exponent;
+                const scalar el = R.lhs()[i].exponent;
                 if (i == j)
                 {
                     if (el < 1.0)
                     {
-                        if (c2[si]>SMALL)
+                        if (c2[si] > SMALL)
                         {
                             kf *= el*pow(c2[si] + VSMALL, el - 1.0);
                         }
@@ -372,31 +374,31 @@ void Foam::ODEChemistryModel<CompType, ThermoType>::jacobian
 
             forAll(R.lhs(), i)
             {
-                label si = R.lhs()[i].index;
+                const label si = R.lhs()[i].index;
-                scalar sl = R.lhs()[i].stoichCoeff;
+                const scalar sl = R.lhs()[i].stoichCoeff;
                 dfdc[si][sj] -= sl*kf;
             }
             forAll(R.rhs(), i)
             {
-                label si = R.rhs()[i].index;
+                const label si = R.rhs()[i].index;
-                scalar sr = R.rhs()[i].stoichCoeff;
+                const scalar sr = R.rhs()[i].stoichCoeff;
                 dfdc[si][sj] += sr*kf;
             }
         }
 
         forAll(R.rhs(), j)
         {
-            label sj = R.rhs()[j].index;
+            const label sj = R.rhs()[j].index;
             scalar kr = kr0;
             forAll(R.rhs(), i)
             {
-                label si = R.rhs()[i].index;
+                const label si = R.rhs()[i].index;
-                scalar er = R.rhs()[i].exponent;
+                const scalar er = R.rhs()[i].exponent;
-                if (i==j)
+                if (i == j)
                 {
-                    if (er<1.0)
+                    if (er < 1.0)
                     {
-                        if (c2[si]>SMALL)
+                        if (c2[si] > SMALL)
                         {
                             kr *= er*pow(c2[si] + VSMALL, er - 1.0);
                         }
@@ -418,25 +420,25 @@ void Foam::ODEChemistryModel<CompType, ThermoType>::jacobian
 
             forAll(R.lhs(), i)
             {
-                label si = R.lhs()[i].index;
+                const label si = R.lhs()[i].index;
-                scalar sl = R.lhs()[i].stoichCoeff;
+                const scalar sl = R.lhs()[i].stoichCoeff;
                 dfdc[si][sj] += sl*kr;
             }
             forAll(R.rhs(), i)
             {
-                label si = R.rhs()[i].index;
+                const label si = R.rhs()[i].index;
-                scalar sr = R.rhs()[i].stoichCoeff;
+                const scalar sr = R.rhs()[i].stoichCoeff;
                 dfdc[si][sj] -= sr*kr;
             }
         }
     }
 
     // calculate the dcdT elements numerically
-    scalar delta = 1.0e-8;
+    const scalar delta = 1.0e-8;
-    scalarField dcdT0 = omega(c2, T - delta, p);
+    const scalarField dcdT0 = omega(c2, T - delta, p);
-    scalarField dcdT1 = omega(c2, T + delta, p);
+    const scalarField dcdT1 = omega(c2, T + delta, p);
 
-    for (label i=0; i<nEqns(); i++)
+    for (label i = 0; i < nEqns(); i++)
     {
         dfdc[i][nSpecie()] = 0.5*(dcdT1[i] - dcdT0[i])/delta;
     }
@@ -485,7 +487,7 @@ Foam::ODEChemistryModel<CompType, ThermoType>::tc() const
 
     scalarField& tc = ttc();
 
-    label nReaction = reactions_.size();
+    const label nReaction = reactions_.size();
 
     if (this->chemistry_)
     {
@@ -550,6 +552,7 @@ Foam::ODEChemistryModel<CompType, ThermoType>::Sh() const
         )
     );
 
+
     if (this->chemistry_)
     {
         scalarField& Sh = tSh();
@@ -558,7 +561,7 @@ Foam::ODEChemistryModel<CompType, ThermoType>::Sh() const
         {
             forAll(Sh, cellI)
             {
-                scalar hi = specieThermo_[i].Hc();
+                const scalar hi = specieThermo_[i].Hc();
                 Sh[cellI] -= hi*RR_[i][cellI];
             }
         }
@@ -626,34 +629,31 @@ void Foam::ODEChemistryModel<CompType, ThermoType>::calculate()
         this->thermo().rho()
     );
 
-    for (label i=0; i<nSpecie_; i++)
+    if (this->mesh().changing())
     {
-        RR_[i].setSize(rho.size());
+        for (label i=0; i<nSpecie_; i++)
+        {
+            RR_[i].setSize(rho.size());
+            RR_[i] = 0.0;
+        }
     }
 
     if (this->chemistry_)
     {
         forAll(rho, celli)
         {
+            const scalar rhoi = rho[celli];
+            const scalar Ti = this->thermo().T()[celli];
+            const scalar pi = this->thermo().p()[celli];
+
+            scalarField c(nSpecie_, 0.0);
             for (label i=0; i<nSpecie_; i++)
             {
-                RR_[i][celli] = 0.0;
+                const scalar Yi = Y_[i][celli];
-            }
-
-            scalar rhoi = rho[celli];
-            scalar Ti = this->thermo().T()[celli];
-            scalar pi = this->thermo().p()[celli];
-
-            scalarField c(nSpecie_);
-            scalarField dcdt(nEqns(), 0.0);
-
-            for (label i=0; i<nSpecie_; i++)
-            {
-                scalar Yi = Y_[i][celli];
                 c[i] = rhoi*Yi/specieThermo_[i].W();
             }
 
-            dcdt = omega(c, Ti, pi);
+            const scalarField dcdt = omega(c, Ti, pi);
 
             for (label i=0; i<nSpecie_; i++)
             {
@@ -671,6 +671,8 @@ Foam::scalar Foam::ODEChemistryModel<CompType, ThermoType>::solve
     const scalar deltaT
 )
 {
+    scalar deltaTMin = GREAT;
+
     const volScalarField rho
     (
         IOobject
@@ -685,35 +687,33 @@ Foam::scalar Foam::ODEChemistryModel<CompType, ThermoType>::solve
         this->thermo().rho()
     );
 
-    for (label i=0; i<nSpecie_; i++)
+    if (this->mesh().changing())
     {
-        RR_[i].setSize(rho.size());
+        for (label i = 0; i < nSpecie_; i++)
+        {
+            RR_[i].setSize(this->mesh().nCells());
+            RR_[i] = 0.0;
+        }
     }
 
     if (!this->chemistry_)
     {
-        return GREAT;
+        return deltaTMin;
     }
 
-    scalar deltaTMin = GREAT;
 
     tmp<volScalarField> thc = this->thermo().hc();
     const scalarField& hc = thc();
 
     forAll(rho, celli)
     {
-        for (label i=0; i<nSpecie_; i++)
+        const scalar rhoi = rho[celli];
-        {
+        const scalar hi = this->thermo().hs()[celli] + hc[celli];
-            RR_[i][celli] = 0.0;
+        const scalar pi = this->thermo().p()[celli];
-        }
-
-        scalar rhoi = rho[celli];
         scalar Ti = this->thermo().T()[celli];
-        scalar hi = this->thermo().hs()[celli] + hc[celli];
-        scalar pi = this->thermo().p()[celli];
 
-        scalarField c(nSpecie_);
+        scalarField c(nSpecie_, 0.0);
-        scalarField c0(nSpecie_);
+        scalarField c0(nSpecie_, 0.0);
         scalarField dc(nSpecie_, 0.0);
 
         for (label i=0; i<nSpecie_; i++)
@@ -722,21 +722,20 @@ Foam::scalar Foam::ODEChemistryModel<CompType, ThermoType>::solve
         }
         c0 = c;
 
+        // initialise timing parameters
         scalar t = t0;
         scalar tauC = this->deltaTChem_[celli];
         scalar dt = min(deltaT, tauC);
         scalar timeLeft = deltaT;
 
         // calculate the chemical source terms
-        scalar cTot = 0.0;
-
         while (timeLeft > SMALL)
         {
             tauC = solver().solve(c, Ti, pi, t, dt);
             t += dt;
 
             // update the temperature
-            cTot = sum(c);
+            const scalar cTot = sum(c);
             ThermoType mixture(0.0*specieThermo_[0]);
             for (label i=0; i<nSpecie_; i++)
             {
@@ -746,19 +745,11 @@ Foam::scalar Foam::ODEChemistryModel<CompType, ThermoType>::solve
 
             timeLeft -= dt;
             this->deltaTChem_[celli] = tauC;
-            dt = min(timeLeft, tauC);
+            dt = max(SMALL, min(timeLeft, tauC));
-            dt = max(dt, SMALL);
         }
         deltaTMin = min(tauC, deltaTMin);
 
         dc = c - c0;
-        scalar WTot = 0.0;
-        for (label i=0; i<nSpecie_; i++)
-        {
-            WTot += c[i]*specieThermo_[i].W();
-        }
-        WTot /= cTot;
-
         for (label i=0; i<nSpecie_; i++)
         {
             RR_[i][celli] = dc[i]*specieThermo_[i].W()/deltaT;

src/thermophysicalModels/chemistryModel/chemistryModel/ODEChemistryModel/ODEChemistryModel.H

@@ -141,7 +141,7 @@ public:
         inline const chemistrySolver<CompType, ThermoType>& solver() const;
 
         //- dc/dt = omega, rate of change in concentration, for each species
-        virtual scalarField omega
+        virtual tmp<scalarField> omega
         (
             const scalarField& c,
             const scalar T,

src/thermophysicalModels/chemistryModel/chemistrySolver/EulerImplicit/EulerImplicit.C

@@ -39,7 +39,7 @@ Foam::EulerImplicit<CompType, ThermoType>::EulerImplicit
     chemistrySolver<CompType, ThermoType>(model, modelName),
     coeffsDict_(model.subDict(modelName + "Coeffs")),
     cTauChem_(readScalar(coeffsDict_.lookup("cTauChem"))),
-    equil_(coeffsDict_.lookup("equilibriumRateLimiter"))
+    eqRateLimiter_(coeffsDict_.lookup("equilibriumRateLimiter"))
 {}
 
 
@@ -65,12 +65,12 @@ Foam::scalar Foam::EulerImplicit<CompType, ThermoType>::solve
     const label nSpecie = this->model_.nSpecie();
     simpleMatrix<scalar> RR(nSpecie, 0, 0);
 
-    for (label i=0; i<nSpecie; i++)
+    for (label i = 0; i < nSpecie; i++)
     {
         c[i] = max(0.0, c[i]);
     }
 
-    for (label i=0; i<nSpecie; i++)
+    for (label i = 0; i < nSpecie; i++)
     {
         RR.source()[i] = c[i]/dt;
     }
@@ -88,9 +88,9 @@ Foam::scalar Foam::EulerImplicit<CompType, ThermoType>::solve
         );
 
         scalar corr = 1.0;
-        if (equil_)
+        if (eqRateLimiter_)
         {
-            if (omegai<0.0)
+            if (omegai < 0.0)
             {
                 corr = 1.0/(1.0 + pr*dt);
             }
@@ -100,31 +100,31 @@ Foam::scalar Foam::EulerImplicit<CompType, ThermoType>::solve
             }
         }
 
-        forAll(R.lhs(), s)
+        forAll(R.lhs(), specieI)
         {
-            label si = R.lhs()[s].index;
+            const label id = R.lhs()[specieI].index;
-            scalar sl = R.lhs()[s].stoichCoeff;
+            const scalar sc = R.lhs()[specieI].stoichCoeff;
-            RR[si][rRef] -= sl*pr*corr;
+            RR[id][rRef] -= sc*pr*corr;
-            RR[si][lRef] += sl*pf*corr;
+            RR[id][lRef] += sc*pf*corr;
         }
 
-        forAll(R.rhs(), s)
+        forAll(R.rhs(), specieI)
         {
-            label si = R.rhs()[s].index;
+            const label id = R.rhs()[specieI].index;
-            scalar sr = R.rhs()[s].stoichCoeff;
+            const scalar sc = R.rhs()[specieI].stoichCoeff;
-            RR[si][lRef] -= sr*pf*corr;
+            RR[id][lRef] -= sc*pf*corr;
-            RR[si][rRef] += sr*pr*corr;
+            RR[id][rRef] += sc*pr*corr;
         }
     }
 
 
-    for (label i=0; i<nSpecie; i++)
+    for (label i = 0; i < nSpecie; i++)
     {
         RR[i][i] += 1.0/dt;
     }
 
     c = RR.LUsolve();
-    for (label i=0; i<nSpecie; i++)
+    for (label i = 0; i < nSpecie; i++)
     {
         c[i] = max(0.0, c[i]);
     }
@@ -134,7 +134,7 @@ Foam::scalar Foam::EulerImplicit<CompType, ThermoType>::solve
     const label nEqns = this->model_.nEqns();
     scalarField c1(nEqns, 0.0);
 
-    for (label i=0; i<nSpecie; i++)
+    for (label i = 0; i < nSpecie; i++)
     {
         c1[i] = c[i];
     }
@@ -144,9 +144,9 @@ Foam::scalar Foam::EulerImplicit<CompType, ThermoType>::solve
     scalarField dcdt(nEqns, 0.0);
     this->model_.derivatives(0.0, c1, dcdt);
 
-    scalar sumC = sum(c);
+    const scalar sumC = sum(c);
 
-    for (label i=0; i<nSpecie; i++)
+    for (label i = 0; i < nSpecie; i++)
     {
         scalar d = dcdt[i];
         if (d < -SMALL)

src/thermophysicalModels/chemistryModel/chemistrySolver/EulerImplicit/EulerImplicit.H

@@ -57,12 +57,17 @@ class EulerImplicit
 {
     // Private data
 
+        //- Coefficients dictionary
         dictionary coeffsDict_;
 
+
         // Model constants
 
+            //- Chemistry timescale
             scalar cTauChem_;
-            Switch equil_;
+
+            //- Equilibrium rate limiter flag (on/off)
+            Switch eqRateLimiter_;
 
 
 public:

src/thermophysicalModels/chemistryModel/chemistrySolver/chemistrySolver/makeChemistrySolvers.C

@@ -29,6 +29,8 @@ License
 #include "psiChemistryModel.H"
 #include "rhoChemistryModel.H"
 
+#include "noChemistrySolver.H"
+
 #include "EulerImplicit.H"
 #include "ode.H"
 #include "sequential.H"
@@ -38,12 +40,24 @@ License
 namespace Foam
 {
     makeChemistrySolver(psiChemistryModel, gasThermoPhysics)
+    makeChemistrySolverType
+    (
+        noChemistrySolver,
+        psiChemistryModel,
+        gasThermoPhysics
+    )
     makeChemistrySolverType(EulerImplicit, psiChemistryModel, gasThermoPhysics)
     makeChemistrySolverType(ode, psiChemistryModel, gasThermoPhysics)
     makeChemistrySolverType(sequential, psiChemistryModel, gasThermoPhysics)
 
     makeChemistrySolver(psiChemistryModel, icoPoly8ThermoPhysics)
     makeChemistrySolverType
+    (
+        noChemistrySolver,
+        psiChemistryModel,
+        icoPoly8ThermoPhysics
+    )
+    makeChemistrySolverType
     (
         EulerImplicit,
         psiChemistryModel,
@@ -58,12 +72,24 @@ namespace Foam
     )
 
     makeChemistrySolver(rhoChemistryModel, gasThermoPhysics)
+    makeChemistrySolverType
+    (
+        noChemistrySolver,
+        rhoChemistryModel,
+        gasThermoPhysics
+    )
     makeChemistrySolverType(EulerImplicit, rhoChemistryModel, gasThermoPhysics)
     makeChemistrySolverType(ode, rhoChemistryModel, gasThermoPhysics)
     makeChemistrySolverType(sequential, rhoChemistryModel, gasThermoPhysics)
 
     makeChemistrySolver(rhoChemistryModel, icoPoly8ThermoPhysics)
     makeChemistrySolverType
+    (
+        noChemistrySolver,
+        rhoChemistryModel,
+        icoPoly8ThermoPhysics
+    )
+    makeChemistrySolverType
     (
         EulerImplicit,
         rhoChemistryModel,

src/thermophysicalModels/chemistryModel/chemistrySolver/noChemistrySolver/noChemistrySolver.C

@@ -0,0 +1,66 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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 "noChemistrySolver.H"
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+template<class CompType, class ThermoType>
+Foam::noChemistrySolver<CompType, ThermoType>::noChemistrySolver
+(
+    ODEChemistryModel<CompType, ThermoType>& model,
+    const word& modelName
+)
+:
+    chemistrySolver<CompType, ThermoType>(model, modelName)
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+template<class CompType, class ThermoType>
+Foam::noChemistrySolver<CompType, ThermoType>::~noChemistrySolver()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+template<class CompType, class ThermoType>
+Foam::scalar Foam::noChemistrySolver<CompType, ThermoType>::solve
+(
+    scalarField&,
+    const scalar,
+    const scalar,
+    const scalar,
+    const scalar
+) const
+{
+    return GREAT;
+}
+
+
+// ************************************************************************* //

src/thermophysicalModels/chemistryModel/chemistrySolver/noChemistrySolver/noChemistrySolver.H

@@ -0,0 +1,110 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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
+
+Class
+    Foam::noChemistry
+
+Description
+    Dummy chemistry solver for 'none' option
+
+SourceFiles
+    noChemistrySolver.H
+    noChemistrySolver.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef noChemistySolver_H
+#define noChemistySolver_H
+
+#include "chemistrySolver.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// Forward declaration of classes
+template<class CompType, class ThermoType>
+class noChemistrySolver;
+
+/*---------------------------------------------------------------------------*\
+                      Class noChemistrySolver Declaration
+\*---------------------------------------------------------------------------*/
+
+template<class CompType, class ThermoType>
+class noChemistrySolver
+:
+    public chemistrySolver<CompType, ThermoType>
+{
+
+public:
+
+    //- Runtime type information
+    TypeName("none");
+
+
+    // Constructors
+
+
+        //- Construct from components
+        noChemistrySolver
+        (
+            ODEChemistryModel<CompType, ThermoType>& model,
+            const word& modelName
+        );
+
+
+    //- Destructor
+    virtual ~noChemistrySolver();
+
+
+    // Member Functions
+
+        //- Update the concentrations and return the chemical time
+        scalar solve
+        (
+            scalarField &c,
+            const scalar T,
+            const scalar p,
+            const scalar t0,
+            const scalar dt
+        ) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+#   include "noChemistrySolver.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/thermophysicalModels/chemistryModel/chemistrySolver/ode/ode.C

@@ -39,8 +39,7 @@ Foam::ode<CompType, ThermoType>::ode
     coeffsDict_(model.subDict(modelName + "Coeffs")),
     solverName_(coeffsDict_.lookup("ODESolver")),
     odeSolver_(ODESolver::New(solverName_, model)),
-    eps_(readScalar(coeffsDict_.lookup("eps"))),
+    eps_(readScalar(coeffsDict_.lookup("eps")))
-    scale_(readScalar(coeffsDict_.lookup("scale")))
 {}
 
 
@@ -67,7 +66,7 @@ Foam::scalar Foam::ode<CompType, ThermoType>::solve
     scalarField c1(this->model_.nEqns(), 0.0);
 
     // copy the concentration, T and P to the total solve-vector
-    for (label i=0; i<nSpecie; i++)
+    for (label i = 0; i < nSpecie; i++)
     {
         c1[i] = c[i];
     }

src/thermophysicalModels/chemistryModel/chemistrySolver/ode/ode.H

@@ -65,7 +65,6 @@ class ode
         // Model constants
 
             scalar eps_;
-            scalar scale_;
 
 
 public:

src/thermophysicalModels/chemistryModel/chemistrySolver/sequential/sequential.C

@@ -38,7 +38,7 @@ Foam::sequential<CompType, ThermoType>::sequential
     chemistrySolver<CompType, ThermoType>(model, modelName),
     coeffsDict_(model.subDict(modelName + "Coeffs")),
     cTauChem_(readScalar(coeffsDict_.lookup("cTauChem"))),
-    equil_(coeffsDict_.lookup("equilibriumRateLimiter"))
+    eqRateLimiter_(coeffsDict_.lookup("equilibriumRateLimiter"))
 {}
 
 
@@ -70,45 +70,41 @@ Foam::scalar Foam::sequential<CompType, ThermoType>::solve
     {
         const Reaction<ThermoType>& R = this->model_.reactions()[i];
 
-        scalar om0 = this->model_.omega
+        scalar omega = this->model_.omega
         (
             R, c, T, p, pf, cf, lRef, pb, cb, rRef
         );
 
-        scalar omeg = 0.0;
+        if (eqRateLimiter_)
-        if (!equil_)
         {
-            omeg = om0;
+            if (omega < 0.0)
-        }
-        else
-        {
-            if (om0<0.0)
             {
-                omeg = om0/(1.0 + pb*dt);
+                omega /= 1.0 + pb*dt;
             }
             else
             {
-                omeg = om0/(1.0 + pf*dt);
+                omega /= 1.0 + pf*dt;
             }
         }
-        tChemInv = max(tChemInv, mag(omeg));
+
+        tChemInv = max(tChemInv, mag(omega));
 
 
         // update species
-        forAll(R.lhs(), s)
+        forAll(R.lhs(), specieI)
         {
-            label si = R.lhs()[s].index;
+            const label id = R.lhs()[specieI].index;
-            scalar sl = R.lhs()[s].stoichCoeff;
+            const scalar sc = R.lhs()[specieI].stoichCoeff;
-            c[si] -= dt*sl*omeg;
+            c[id] -= dt*sc*omega;
-            c[si] = max(0.0, c[si]);
+            c[id] = max(0.0, c[id]);
         }
 
-        forAll(R.rhs(), s)
+        forAll(R.rhs(), specieI)
         {
-            label si = R.rhs()[s].index;
+            const label id = R.rhs()[specieI].index;
-            scalar sr = R.rhs()[s].stoichCoeff;
+            const scalar sc = R.rhs()[specieI].stoichCoeff;
-            c[si] += dt*sr*omeg;
+            c[id] += dt*sc*omega;
-            c[si] = max(0.0, c[si]);
+            c[id] = max(0.0, c[id]);
         }
     }
 

src/thermophysicalModels/chemistryModel/chemistrySolver/sequential/sequential.H

@@ -59,12 +59,17 @@ class sequential
 {
     // Private data
 
+        //- Coefficients dictionary
         dictionary coeffsDict_;
 
+
         // Model constants
 
+            //- Chemistry time scale
             scalar cTauChem_;
-            Switch equil_;
+
+            //- Equilibrium rate limiter flag (on/off)
+            Switch eqRateLimiter_;
 
 
 public:
@@ -75,7 +80,6 @@ public:
 
     // Constructors
 
-
         //- Construct from components
         sequential
         (