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

applications/solvers/combustion/fireFoam/pEqn.H

@@ -4,6 +4,9 @@ volScalarField rAU(1.0/UEqn.A());
 surfaceScalarField rhorAUf("Dp", fvc::interpolate(rho*rAU));
 volVectorField HbyA("HbyA", U);
 HbyA = rAU*UEqn.H();
+phi.boundaryField() =
+    fvc::interpolate(rho.boundaryField()*U.boundaryField())
+  & mesh.Sf().boundaryField();
 
 surfaceScalarField phig(-rhorAUf*ghf*fvc::snGrad(rho)*mesh.magSf());
 

applications/solvers/multiphase/bubbleFoam/pEqn.H

@@ -38,16 +38,18 @@
     (
         (fvc::interpolate(HbyA1) & mesh.Sf())
       + fvc::ddtPhiCorr(rAU1, U1, phi1)
-      + phiDrag1
     );
 
     surfaceScalarField phiHbyA2
     (
         (fvc::interpolate(HbyA2) & mesh.Sf())
       + fvc::ddtPhiCorr(rAU2, U2, phi2)
-      + phiDrag2
     );
 
+    phi = alpha1f*phiHbyA1 + alpha2f*phiHbyA2;
+
+    phiHbyA1 += phiDrag1;
+    phiHbyA2 += phiDrag2;
     surfaceScalarField phiHbyA("phiHbyA", alpha1f*phiHbyA1 + alpha2f*phiHbyA2);
 
     surfaceScalarField Dp

applications/solvers/multiphase/cavitatingFoam/UEqn.H

@@ -1,17 +1,8 @@
-    surfaceScalarField muEff
-    (
-        "muEff",
-        twoPhaseProperties.muf()
-      + fvc::interpolate(rho*turbulence->nut())
-    );
-
     fvVectorMatrix UEqn
     (
         fvm::ddt(rho, U)
       + fvm::div(phi, U)
-      - fvm::laplacian(muEff, U)
+      + turbulence->divDevRhoReff(rho, U)
-    //- (fvc::grad(U) & fvc::grad(muf))
-      - fvc::div(muEff*(fvc::interpolate(dev(fvc::grad(U))) & mesh.Sf()))
     );
 
     UEqn.relax();

applications/solvers/multiphase/compressibleInterFoam/UEqn.H

@@ -1,17 +1,8 @@
-    surfaceScalarField muEff
-    (
-        "muEff",
-        twoPhaseProperties.muf()
-      + fvc::interpolate(rho*turbulence->nut())
-    );
-
     fvVectorMatrix UEqn
     (
         fvm::ddt(rho, U)
       + fvm::div(rhoPhi, U)
-      - fvm::laplacian(muEff, U)
+      + turbulence->divDevRhoReff(rho, U)
-      - (fvc::grad(U) & fvc::grad(muEff))
-    //- fvc::div(muEff*(fvc::interpolate(dev(fvc::grad(U))) & mesh.Sf()))
     );
 
     UEqn.relax();

applications/solvers/multiphase/compressibleInterFoam/pEqn.H

@@ -36,6 +36,7 @@
         (fvc::interpolate(HbyA) & mesh.Sf())
       + fvc::ddtPhiCorr(rAU, rho, U, phi)
     );
+    phi = phiHbyA;
 
     surfaceScalarField phig
     (

applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/pEqn.H

@@ -2,11 +2,6 @@
     rho1 = rho10 + psi1*p;
     rho2 = rho20 + psi2*p;
 
-    mrfZones.absoluteFlux(phi1.oldTime());
-    mrfZones.absoluteFlux(phi1);
-    mrfZones.absoluteFlux(phi2.oldTime());
-    mrfZones.absoluteFlux(phi2);
-
     surfaceScalarField alpha1f(fvc::interpolate(alpha1));
     surfaceScalarField alpha2f(scalar(1) - alpha1f);
 
@@ -22,6 +17,11 @@
     volVectorField HbyA2("HbyA2", U2);
     HbyA2 = rAU2*U2Eqn.H();
 
+    mrfZones.absoluteFlux(phi1.oldTime());
+    mrfZones.absoluteFlux(phi1);
+    mrfZones.absoluteFlux(phi2.oldTime());
+    mrfZones.absoluteFlux(phi2);
+
     surfaceScalarField ppDrag("ppDrag", 0.0*phi1);
 
     if (g0.value() > 0.0)
@@ -41,18 +41,29 @@
         "phiHbyA1",
         (fvc::interpolate(HbyA1) & mesh.Sf())
       + fvc::ddtPhiCorr(rAU1, alpha1, U1, phi1)
-      + fvc::interpolate((1.0/rho1)*rAU1*dragCoeff)*phi2
-      + ppDrag
-      + rAlphaAU1f*(g & mesh.Sf())
     );
-    mrfZones.relativeFlux(phiHbyA1);
 
     surfaceScalarField phiHbyA2
     (
         "phiHbyA2",
         (fvc::interpolate(HbyA2) & mesh.Sf())
       + fvc::ddtPhiCorr(rAU2, alpha2, U2, phi2)
-      + fvc::interpolate((1.0/rho2)*rAU2*dragCoeff)*phi1
+    );
+
+    phi = alpha1f*phiHbyA1 + alpha2f*phiHbyA2;
+    mrfZones.relativeFlux(phi);
+
+    phiHbyA1 +=
+    (
+        fvc::interpolate((1.0/rho1)*rAU1*dragCoeff)*phi2
+      + ppDrag
+      + rAlphaAU1f*(g & mesh.Sf())
+    );
+    mrfZones.relativeFlux(phiHbyA1);
+
+    phiHbyA2 +=
+    (
+        fvc::interpolate((1.0/rho2)*rAU2*dragCoeff)*phi1
       + rAlphaAU2f*(g & mesh.Sf())
     );
     mrfZones.relativeFlux(phiHbyA2);

applications/solvers/multiphase/interFoam/MRFInterFoam/pEqn.H

@@ -11,9 +11,9 @@
         (fvc::interpolate(HbyA) & mesh.Sf())
       + fvc::ddtPhiCorr(rAU, rho, U, phi)
     );
-    mrfZones.relativeFlux(phiHbyA);
-
     adjustPhi(phiHbyA, U, p_rgh);
+    mrfZones.relativeFlux(phiHbyA);
+    phi = phiHbyA;
 
     surfaceScalarField phig
     (

applications/solvers/multiphase/interFoam/UEqn.H

@@ -1,17 +1,8 @@
-    surfaceScalarField muEff
-    (
-        "muEff",
-        twoPhaseProperties.muf()
-      + fvc::interpolate(rho*turbulence->nut())
-    );
-
     fvVectorMatrix UEqn
     (
         fvm::ddt(rho, U)
       + fvm::div(rhoPhi, U)
-      - fvm::laplacian(muEff, U)
+      + turbulence->divDevRhoReff(rho, U)
-      - (fvc::grad(U) & fvc::grad(muEff))
-    //- fvc::div(muEff*(fvc::interpolate(dev(fvc::grad(U))) & mesh.Sf()))
     );
 
     UEqn.relax();

applications/solvers/multiphase/interFoam/interDyMFoam/pEqn.H

@@ -19,6 +19,8 @@
         fvc::makeAbsolute(phiHbyA, U);
     }
 
+    phiAbs = phiHbyA;
+
     surfaceScalarField phig
     (
         (

applications/solvers/multiphase/interFoam/pEqn.H

@@ -13,6 +13,7 @@
     );
 
     adjustPhi(phiHbyA, U, p_rgh);
+    phi = phiHbyA;
 
     surfaceScalarField phig
     (

applications/solvers/multiphase/interPhaseChangeFoam/UEqn.H

@@ -1,18 +1,9 @@
-    surfaceScalarField muEff
-    (
-        "muEff",
-        twoPhaseProperties->muf()
-      + fvc::interpolate(rho*turbulence->nut())
-    );
-
     fvVectorMatrix UEqn
     (
         fvm::ddt(rho, U)
       + fvm::div(rhoPhi, U)
       - fvm::Sp(fvc::ddt(rho) + fvc::div(rhoPhi), U)
-      - fvm::laplacian(muEff, U)
+      + turbulence->divDevRhoReff(rho, U)
-      - (fvc::grad(U) & fvc::grad(muEff))
-    //- fvc::div(muEff*(fvc::interpolate(dev2(fvc::grad(U))) & mesh.Sf()))
     );
 
     UEqn.relax();

applications/solvers/multiphase/interPhaseChangeFoam/pEqn.H

@@ -11,8 +11,8 @@
         (fvc::interpolate(HbyA) & mesh.Sf())
       + fvc::ddtPhiCorr(rAU, rho, U, phi)
     );
-
     adjustPhi(phiHbyA, U, p_rgh);
+    phi = phiHbyA;
 
     surfaceScalarField phig
     (

applications/solvers/multiphase/multiphaseEulerFoam/pEqn.H

@@ -58,6 +58,8 @@
         dimensionedScalar("phiHbyA", dimArea*dimVelocity, 0)
     );
 
+    phi = dimensionedScalar("phi", phi.dimensions(), 0);
+
     phasei = 0;
     forAllIter(PtrDictionary<phaseModel>, fluid.phases(), iter)
     {
@@ -74,15 +76,18 @@
         (
             (fvc::interpolate(HbyAs[phasei]) & mesh.Sf())
           + fvc::ddtPhiCorr(rAUs[phasei], alpha, phase.U(), phase.phi())
-          + rAlphaAUfs[phasei]
+        );
+        mrfZones.relativeFlux(phiHbyAs[phasei]);
+
+        phi += alphafs[phasei]*phiHbyAs[phasei];
+
+        phiHbyAs[phasei] +=
+            rAlphaAUfs[phasei]
            *(
                fluid.surfaceTension(phase)*mesh.magSf()/phase.rho()
              + (g & mesh.Sf())
-            )
             );
 
-        mrfZones.relativeFlux(phiHbyAs[phasei]);
-
         multiphaseSystem::dragModelTable::const_iterator dmIter =
             fluid.dragModels().begin();
         multiphaseSystem::dragCoeffFields::const_iterator dcIter =

applications/solvers/multiphase/multiphaseInterFoam/MRFMultiphaseInterFoam/pEqn.H

@@ -11,9 +11,9 @@
         (fvc::interpolate(HbyA) & mesh.Sf())
       + fvc::ddtPhiCorr(rAU, rho, U, phi)
     );
-    mrfZones.relativeFlux(phiHbyA);
-
     adjustPhi(phiHbyA, U, p_rgh);
+    mrfZones.relativeFlux(phiHbyA);
+    phi = phiHbyA;
 
     surfaceScalarField phig
     (

applications/solvers/multiphase/multiphaseInterFoam/UEqn.H

@@ -1,17 +1,8 @@
-    surfaceScalarField muEff
-    (
-        "muEff",
-        mixture.muf()
-      + fvc::interpolate(rho*turbulence->nut())
-    );
-
     fvVectorMatrix UEqn
     (
         fvm::ddt(rho, U)
       + fvm::div(mixture.rhoPhi(), U)
-      - fvm::laplacian(muEff, U)
+      + turbulence->divDevRhoReff(rho, U)
-      - (fvc::grad(U) & fvc::grad(muEff))
-    //- fvc::div(muEff*(fvc::interpolate(dev(fvc::grad(U))) & mesh.Sf()))
     );
 
     UEqn.relax();

applications/solvers/multiphase/multiphaseInterFoam/pEqn.H

@@ -11,8 +11,8 @@
         (fvc::interpolate(HbyA) & mesh.Sf())
       + fvc::ddtPhiCorr(rAU, rho, U, phi)
     );
-
     adjustPhi(phiHbyA, U, p_rgh);
+    phi = phiHbyA;
 
     surfaceScalarField phig
     (

applications/solvers/multiphase/settlingFoam/UEqn.H

@@ -11,7 +11,6 @@
         )
       - fvm::laplacian(muEff, U, "laplacian(muEff,U)")
       - (fvc::grad(U) & fvc::grad(muEff))
-    //- fvc::div(muEff*dev2(T(fvc::grad(U))))
     );
 
     UEqn.relax();

applications/solvers/multiphase/settlingFoam/pEqn.H

@@ -15,6 +15,7 @@
          + fvc::ddtPhiCorr(rAU, rho, U, phi)
         )
     );
+    phi = phiHbyA;
 
     surfaceScalarField phig
     (

applications/solvers/multiphase/twoLiquidMixingFoam/UEqn.H

@@ -1,17 +1,8 @@
-    surfaceScalarField muEff
-    (
-        "muEff",
-        twoPhaseProperties.muf()
-      + fvc::interpolate(rho*turbulence->nut())
-    );
-
     fvVectorMatrix UEqn
     (
         fvm::ddt(rho, U)
       + fvm::div(rhoPhi, U)
-      - fvm::laplacian(muEff, U)
+      + turbulence->divDevRhoReff(rho, U)
-      - (fvc::grad(U) & fvc::grad(muEff))
-    //- fvc::div(muEff*(fvc::interpolate(dev(fvc::grad(U))) & mesh.Sf()))
     );
 
     UEqn.relax();

applications/solvers/multiphase/twoLiquidMixingFoam/pEqn.H

@@ -11,8 +11,8 @@
         (fvc::interpolate(HbyA) & mesh.Sf())
       + fvc::ddtPhiCorr(rAU, rho, U, phi)
     );
-
     adjustPhi(phiHbyA, U, p_rgh);
+    phi = phiHbyA;
 
     surfaceScalarField phig
     (

applications/solvers/multiphase/twoPhaseEulerFoam/pEqn.H

@@ -36,18 +36,29 @@
         "phiHbyA1",
         (fvc::interpolate(HbyA1) & mesh.Sf())
       + fvc::ddtPhiCorr(rAU1, U1, phi1)
-      + fvc::interpolate(alpha2/rho1*K*rAU1)*phi2
-      + ppDrag
-      + rAU1f*(g & mesh.Sf())
     );
-    mrfZones.relativeFlux(phiHbyA1);
 
     surfaceScalarField phiHbyA2
     (
         "phiHbyA2",
         (fvc::interpolate(HbyA2) & mesh.Sf())
       + fvc::ddtPhiCorr(rAU2, U2, phi2)
-      + fvc::interpolate(alpha1/rho2*K*rAU2)*phi1
+    );
+
+    phi = alpha1f*phiHbyA1 + alpha2f*phiHbyA2;
+    mrfZones.relativeFlux(phi);
+
+    phiHbyA1 +=
+    (
+        fvc::interpolate(alpha2/rho1*K*rAU1)*phi2
+      + ppDrag
+      + rAU1f*(g & mesh.Sf())
+    );
+    mrfZones.relativeFlux(phiHbyA1);
+
+    phiHbyA2 +=
+    (
+        fvc::interpolate(alpha1/rho2*K*rAU2)*phi1
       + rAU2f*(g & mesh.Sf())
     );
     mrfZones.relativeFlux(phiHbyA2);

applications/solvers/multiphase/twoPhaseEulerFoam/pEqn.H.old

@@ -1,118 +0,0 @@
-{
-    surfaceScalarField alpha1f(fvc::interpolate(alpha1));
-    surfaceScalarField alpha2f(scalar(1) - alpha1f);
-
-    volScalarField rAU1(1.0/U1Eqn.A());
-    volScalarField rAU2(1.0/U2Eqn.A());
-
-    rAU1f = fvc::interpolate(rAU1);
-    surfaceScalarField rAU2f(fvc::interpolate(rAU2));
-
-    volVectorField HbyA1("HbyA1", U1);
-    HbyA1 = rAU1*U1Eqn.H();
-
-    volVectorField HbyA2("HbyA2", U2);
-    HbyA2 = rAU2*U2Eqn.H();
-
-    mrfZones.absoluteFlux(phi1.oldTime());
-    mrfZones.absoluteFlux(phi1);
-
-    mrfZones.absoluteFlux(phi2.oldTime());
-    mrfZones.absoluteFlux(phi2);
-
-    surfaceScalarField phiDrag1
-    (
-        fvc::interpolate(alpha2/rho1*K*rAU1)*phi2 + rAU1f*(g & mesh.Sf())
-    );
-
-    if (g0.value() > 0.0)
-    {
-        phiDrag1 -= ppMagf*fvc::snGrad(alpha1)*mesh.magSf();
-    }
-
-    if (kineticTheory.on())
-    {
-        phiDrag1 -= rAU1f*fvc::snGrad(kineticTheory.pa()/rho1)*mesh.magSf();
-    }
-
-
-    surfaceScalarField phiDrag2
-    (
-        fvc::interpolate(alpha1/rho2*K*rAU2)*phi1 + rAU2f*(g & mesh.Sf())
-    );
-
-    // Fix for gravity on outlet boundary.
-    forAll(p.boundaryField(), patchi)
-    {
-        if (isA<zeroGradientFvPatchScalarField>(p.boundaryField()[patchi]))
-        {
-            phiDrag1.boundaryField()[patchi] = 0.0;
-            phiDrag2.boundaryField()[patchi] = 0.0;
-        }
-    }
-
-    surfaceScalarField phiHbyA1
-    (
-        "phiHbyA1",
-        (fvc::interpolate(HbyA1) & mesh.Sf())
-      + fvc::ddtPhiCorr(rAU1, U1, phi1)
-      + phiDrag1
-    );
-    mrfZones.relativeFlux(phiHbyA1);
-
-    surfaceScalarField phiHbyA2
-    (
-        "phiHbyA2",
-        (fvc::interpolate(HbyA2) & mesh.Sf())
-      + fvc::ddtPhiCorr(rAU2, U2, phi2)
-      + phiDrag2
-    );
-    mrfZones.relativeFlux(phiHbyA2);
-
-    mrfZones.relativeFlux(phi1.oldTime());
-    mrfZones.relativeFlux(phi1);
-    mrfZones.relativeFlux(phi2.oldTime());
-    mrfZones.relativeFlux(phi2);
-
-    surfaceScalarField phiHbyA("phiHbyA", alpha1f*phiHbyA1 + alpha2f*phiHbyA2);
-
-    surfaceScalarField Dp
-    (
-        "Dp",
-        alpha1f*rAU1f/rho1 + alpha2f*rAU2f/rho2
-    );
-
-    while (pimple.correctNonOrthogonal())
-    {
-        fvScalarMatrix pEqn
-        (
-            fvm::laplacian(Dp, p) == fvc::div(phiHbyA)
-        );
-
-        pEqn.setReference(pRefCell, pRefValue);
-
-        pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter())));
-
-        if (pimple.finalNonOrthogonalIter())
-        {
-            surfaceScalarField SfGradp(pEqn.flux()/Dp);
-
-            phi1 = phiHbyA1 - rAU1f*SfGradp/rho1;
-            phi2 = phiHbyA2 - rAU2f*SfGradp/rho2;
-            phi = alpha1f*phi1 + alpha2f*phi2;
-
-            p.relax();
-            SfGradp = pEqn.flux()/Dp;
-
-            U1 = HbyA1 + fvc::reconstruct(phiDrag1 - rAU1f*SfGradp/rho1);
-            U1.correctBoundaryConditions();
-
-            U2 = HbyA2 + fvc::reconstruct(phiDrag2 - rAU2f*SfGradp/rho2);
-            U2.correctBoundaryConditions();
-
-            U = alpha1*U1 + alpha2*U2;
-        }
-    }
-}
-
-#include "continuityErrs.H"

applications/solvers/stressAnalysis/solidDisplacementFoam/readMechanicalProperties.H

@@ -161,7 +161,7 @@
             )
         );
     }
-    else if(nuType == "field")
+    else if (nuType == "field")
     {
         nuIO.readOpt() = IOobject::MUST_READ;
         nuPtr.reset

applications/solvers/stressAnalysis/solidDisplacementFoam/readThermalProperties.H

@@ -81,7 +81,7 @@ if (thermalStress)
         );
 
     }
-    else if(CType == "field")
+    else if (CType == "field")
     {
         CIO.readOpt() = IOobject::MUST_READ;
 
@@ -228,5 +228,4 @@ if (thermalStress)
 
     threeKalpha = threeK*alpha;
     DT = k/C;
-
 }

applications/utilities/mesh/manipulation/polyDualMesh/polyDualMeshApp.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -56,7 +56,6 @@ Usage
 
 #include "argList.H"
 #include "Time.H"
-#include "timeSelector.H"
 #include "fvMesh.H"
 #include "unitConversion.H"
 #include "polyTopoChange.H"
@@ -354,7 +353,6 @@ int main(int argc, char *argv[])
 {
 #   include "addOverwriteOption.H"
     argList::noParallel();
-    timeSelector::addOptions(true, false);
 
     argList::validArgs.append("featureAngle [0-180]");
     argList::addBoolOption
@@ -376,9 +374,6 @@ int main(int argc, char *argv[])
 
 #   include "setRootCase.H"
 #   include "createTime.H"
-
-    instantList timeDirs = timeSelector::select0(runTime, args);
-
 #   include "createMesh.H"
 
     const word oldInstance = mesh.pointsInstance();

applications/utilities/mesh/manipulation/topoSet/topoSetDict

@@ -171,13 +171,13 @@ FoamFile
 //        insidePoint (1 2 3);    // point inside region to select
 //    }
 //
-//    // Cells underneath plane such that volume is reached. Can be used
+//    // Cells underneath plane such that volume is reached. E.g. for use
-//    // in setFields.
+//    // in setFields to set the level given a wanted volume.
 //    source targetVolumeToCell;
 //    sourceInfo
 //    {
 //        volume  2e-05;
-//        normal  (1 1 1);
+//        normal  (0 1 0);          // usually in direction of gravity
 //    }
 //
 //

applications/utilities/postProcessing/dataConversion/foamToEnsight/checkMeshMoving.H

@@ -1,21 +1,27 @@
-// check for "points" in all of the result directories
+// check for "points" in any of the result directories
 
-bool meshMoving = true;
+bool meshMoving = false;
-if (Times.size() > 2)
+if (Times.size() > 1)
 {
-    for (label n1=2; n1<Times.size() && meshMoving; ++n1)
+    // We already loaded a mesh (usually from constant). See if any other
+    // points files
+    forAll(Times, timeI)
     {
-        meshMoving = IOobject
+        if (Times[timeI].name() != mesh.pointsInstance())
+        {
+            IOobject io
             (
                 "points",
-            Times[n1].name(),
+                Times[timeI].name(),
                 polyMesh::meshSubDir,
                 mesh,
                 IOobject::NO_READ
-        ).headerOk();
+            );
+            if (io.headerOk())
+            {
+                meshMoving = true;
+                break;
+            }
+        }
     }
 }
-else
-{
-    meshMoving = false;
-}

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

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -1014,6 +1014,7 @@ void Foam::ensightMesh::write
     const fileName& postProcPath,
     const word& prepend,
     const label timeIndex,
+    const bool meshMoving,
     Ostream& ensightCaseFile
 ) const
 {
@@ -1027,7 +1028,7 @@ void Foam::ensightMesh::write
     {
         timeFile += "000.";
     }
-    else if (mesh_.moving())
+    else if (meshMoving)
     {
         timeFile += itoa(timeIndex) + '.';
     }

applications/utilities/postProcessing/dataConversion/foamToEnsight/ensightMesh.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -380,6 +380,7 @@ public:
             const fileName& postProcPath,
             const word& prepend,
             const label timeIndex,
+            const bool meshMoving,
             Ostream& ensightCaseFile
         ) const;
 

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

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -228,6 +228,13 @@ int main(int argc, char *argv[])
 
 #   include "checkMeshMoving.H"
 
+    if (meshMoving)
+    {
+        Info<< "Detected a moving mesh (multiple polyMesh/points files)."
+            << " Writing meshes for every timestep." << endl;
+    }
+
+
     wordHashSet allCloudNames;
     if (Pstream::master())
     {
@@ -346,13 +353,14 @@ int main(int argc, char *argv[])
             eMesh.correct();
         }
 
-        if (timeIndex == 0 || (meshState != polyMesh::UNCHANGED))
+        if (timeIndex == 0 || meshMoving)
         {
             eMesh.write
             (
                 ensightDir,
                 prepend,
                 timeIndex,
+                meshMoving,
                 ensightCaseFile
             );
         }

applications/utilities/preProcessing/setFields/setFields.C

@@ -27,7 +27,6 @@ Description
 \*---------------------------------------------------------------------------*/
 
 #include "argList.H"
-#include "timeSelector.H"
 #include "Time.H"
 #include "fvMesh.H"
 #include "topoSetSource.H"
@@ -360,14 +359,8 @@ public:
 
 int main(int argc, char *argv[])
 {
-    timeSelector::addOptions();
-
 #   include "setRootCase.H"
 #   include "createTime.H"
-
-    // Get times list
-    instantList timeDirs = timeSelector::select0(runTime, args);
-
 #   include "createMesh.H"
 
     Info<< "Reading setFieldsDict\n" << endl;

applications/utilities/surface/surfaceMeshConvert/surfaceMeshConvert.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -58,7 +58,6 @@ Note
 \*---------------------------------------------------------------------------*/
 
 #include "argList.H"
-#include "timeSelector.H"
 #include "Time.H"
 
 #include "MeshedSurfaces.H"

applications/utilities/surface/surfaceMeshConvertTesting/surfaceMeshConvertTesting.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -52,7 +52,6 @@ Note
 \*---------------------------------------------------------------------------*/
 
 #include "argList.H"
-#include "timeSelector.H"
 #include "Time.H"
 #include "polyMesh.H"
 #include "triSurface.H"

applications/utilities/surface/surfaceMeshExport/surfaceMeshExport.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -59,7 +59,6 @@ Note
 \*---------------------------------------------------------------------------*/
 
 #include "argList.H"
-#include "timeSelector.H"
 #include "Time.H"
 
 #include "MeshedSurfaces.H"

applications/utilities/surface/surfaceMeshImport/surfaceMeshImport.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -59,7 +59,6 @@ Note
 \*---------------------------------------------------------------------------*/
 
 #include "argList.H"
-#include "timeSelector.H"
 #include "Time.H"
 
 #include "MeshedSurfaces.H"

applications/utilities/surface/surfaceMeshInfo/surfaceMeshInfo.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -56,7 +56,6 @@ Note
 \*---------------------------------------------------------------------------*/
 
 #include "argList.H"
-#include "timeSelector.H"
 #include "Time.H"
 
 #include "UnsortedMeshedSurfaces.H"

src/OpenFOAM/db/IOobjects/IOdictionary/IOdictionaryIO.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -95,7 +95,9 @@ void Foam::IOdictionary::readFile(const bool masterOnly)
             }
 
             // Note: use ASCII for now - binary IO of dictionaries is
-            // not currently supported
+            // not currently supported or rather the primitiveEntries of
+            // the dictionary think they are in binary form whereas they are
+            // not. Could reset all the ITstreams to ascii?
             IPstream fromAbove
             (
                 Pstream::scheduled,

src/OpenFOAM/db/dictionary/dictionaryIO.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -80,6 +80,12 @@ Foam::autoPtr<Foam::dictionary> Foam::dictionary::New(Istream& is)
 
 bool Foam::dictionary::read(Istream& is, const bool keepHeader)
 {
+    // Check for empty dictionary
+    if (is.eof())
+    {
+        return true;
+    }
+
     if (!is.good())
     {
         FatalIOErrorIn("dictionary::read(Istream&, bool)", is)

src/OpenFOAM/meshes/polyMesh/globalMeshData/globalMeshData.H

@@ -323,8 +323,6 @@ class globalMeshData
             // its own master. Maybe store as well?
 
             void calcGlobalCoPointSlaves() const;
-            const labelListList& globalCoPointSlaves() const;
-            const mapDistribute& globalCoPointSlavesMap() const;
 
 
         //- Disallow default bitwise copy construct
@@ -547,6 +545,11 @@ public:
                 //- Is my edge same orientation as master edge
                 const PackedBoolList& globalEdgeOrientation() const;
 
+            // Collocated point to collocated point
+
+                const labelListList& globalCoPointSlaves() const;
+                const mapDistribute& globalCoPointSlavesMap() const;
+
             // Coupled point to boundary faces. These are uncoupled boundary
             // faces only but include empty patches.
 

src/OpenFOAM/meshes/polyMesh/globalMeshData/globalPoints.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -339,17 +339,15 @@ bool Foam::globalPoints::storeInitialInfo
 }
 
 
-void Foam::globalPoints::printProcPoints
+void Foam::globalPoints::printProcPoint
 (
     const labelList& patchToMeshPoint,
-    const labelPairList& pointInfo
+    const labelPair& pointInfo
 ) const
 {
-    forAll(pointInfo, i)
+    label procI = globalIndexAndTransform::processor(pointInfo);
-    {
+    label index = globalIndexAndTransform::index(pointInfo);
-        label procI = globalIndexAndTransform::processor(pointInfo[i]);
+    label trafoI = globalIndexAndTransform::transformIndex(pointInfo);
-        label index = globalIndexAndTransform::index(pointInfo[i]);
-        label trafoI = globalIndexAndTransform::transformIndex(pointInfo[i]);
 
     Pout<< "    proc:" << procI;
     Pout<< " localpoint:";
@@ -363,7 +361,18 @@ void Foam::globalPoints::printProcPoints
         label meshPointI = localToMeshPoint(patchToMeshPoint, index);
         Pout<< " at:" <<  mesh_.points()[meshPointI];
     }
+}
 
+
+void Foam::globalPoints::printProcPoints
+(
+    const labelList& patchToMeshPoint,
+    const labelPairList& pointInfo
+) const
+{
+    forAll(pointInfo, i)
+    {
+        printProcPoint(patchToMeshPoint, pointInfo[i]);
         Pout<< endl;
     }
 }

src/OpenFOAM/meshes/polyMesh/globalMeshData/globalPoints.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -206,6 +206,12 @@ class globalPoints
         );
 
         //- Debug printing
+        void printProcPoint
+        (
+            const labelList& patchToMeshPoint,
+            const labelPair& pointInfo
+        ) const;
+
         void printProcPoints
         (
             const labelList& patchToMeshPoint,

src/OpenFOAM/meshes/polyMesh/polyMeshClear.C

@@ -97,6 +97,16 @@ void Foam::polyMesh::clearAddressing()
     geometricD_ = Vector<label>::zero;
     solutionD_ = Vector<label>::zero;
 
+    // Update zones
+    pointZones_.clearAddressing();
+    faceZones_.clearAddressing();
+    cellZones_.clearAddressing();
+
+    // Remove the stored tet base points
+    tetBasePtIsPtr_.clear();
+    // Remove the cell tree
+    cellTreePtr_.clear();
+
     pointMesh::Delete(*this);
 }
 

src/OpenFOAM/primitives/functions/DataEntry/DataEntry/DataEntryFwd.H

@@ -0,0 +1,52 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2012 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef DataEntryFws_H
+#define DataEntryFws_H
+
+#include "DataEntry.H"
+#include "vector.H"
+#include "symmTensor.H"
+#include "sphericalTensor.H"
+#include "tensor.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    typedef DataEntry<label> labelDataEntry;
+    typedef DataEntry<scalar> scalarDataEntry;
+    typedef DataEntry<vector> vectorDataEntry;
+    typedef DataEntry<symmTensor> symmTensorDataEntry;
+    typedef DataEntry<sphericalTensor> sphericalTensorDataEntry;
+    typedef DataEntry<tensor> tensorDataEntry;
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/OpenFOAM/primitives/functions/DataEntry/polynomial/polynomial.C

@@ -25,14 +25,14 @@ License
 
 #include "polynomial.H"
 #include "Time.H"
+#include "addToRunTimeSelectionTable.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
 namespace Foam
 {
     defineTypeNameAndDebug(polynomial, 0);
-    DataEntry<scalar>::adddictionaryConstructorToTable<polynomial>
+    addToRunTimeSelectionTable(scalarDataEntry, polynomial, dictionary);
-        addpolynomialConstructorToTable_;
 }
 
 
@@ -40,7 +40,7 @@ namespace Foam
 
 Foam::polynomial::polynomial(const word& entryName, const dictionary& dict)
 :
-    DataEntry<scalar>(entryName),
+    scalarDataEntry(entryName),
     coeffs_(),
     canIntegrate_(true),
     dimensions_(dimless)
@@ -59,8 +59,10 @@ Foam::polynomial::polynomial(const word& entryName, const dictionary& dict)
 
     if (!coeffs_.size())
     {
-        FatalErrorIn("Foam::polynomial::polynomial(const word&, dictionary&)")
+        FatalErrorIn
-            << "polynomial coefficients for entry " << this->name_
+        (
+            "Foam::polynomial::polynomial(const word&, const dictionary&)"
+        )   << "polynomial coefficients for entry " << this->name_
             << " are invalid (empty)" << nl << exit(FatalError);
     }
 
@@ -77,8 +79,10 @@ Foam::polynomial::polynomial(const word& entryName, const dictionary& dict)
     {
         if (!canIntegrate_)
         {
-            WarningIn("Foam::polynomial::polynomial(const word&, dictionary&)")
+            WarningIn
-                << "Polynomial " << this->name_ << " cannot be integrated"
+            (
+                "Foam::polynomial::polynomial(const word&, const dictionary&)"
+            )   << "Polynomial " << this->name_ << " cannot be integrated"
                 << endl;
         }
     }
@@ -91,7 +95,7 @@ Foam::polynomial::polynomial
     const List<Tuple2<scalar, scalar> >& coeffs
 )
 :
-    DataEntry<scalar>(entryName),
+    scalarDataEntry(entryName),
     coeffs_(coeffs),
     canIntegrate_(true),
     dimensions_(dimless)
@@ -101,7 +105,7 @@ Foam::polynomial::polynomial
         FatalErrorIn
         (
             "Foam::polynomial::polynomial"
-            "(const word&, const List<Tuple2<scalar, scalar> >&&)"
+            "(const word&, const List<Tuple2<scalar, scalar> >&)"
         )   << "polynomial coefficients for entry " << this->name_
             << " are invalid (empty)" << nl << exit(FatalError);
     }
@@ -122,7 +126,7 @@ Foam::polynomial::polynomial
             WarningIn
             (
                 "Foam::polynomial::polynomial"
-                "(const word&, const List<Tuple2<scalar, scalar> >&&)"
+                "(const word&, const List<Tuple2<scalar, scalar> >&)"
             )   << "Polynomial " << this->name_ << " cannot be integrated"
                 << endl;
         }
@@ -132,7 +136,7 @@ Foam::polynomial::polynomial
 
 Foam::polynomial::polynomial(const polynomial& poly)
 :
-    DataEntry<scalar>(poly),
+    scalarDataEntry(poly),
     coeffs_(poly.coeffs_),
     canIntegrate_(poly.canIntegrate_),
     dimensions_(poly.dimensions_)
@@ -201,7 +205,8 @@ Foam::dimensioned<Foam::scalar> Foam::polynomial::dimValue
 
 Foam::dimensioned<Foam::scalar> Foam::polynomial::dimIntegrate
 (
-    const scalar x1, const scalar x2
+    const scalar x1,
+    const scalar x2
 ) const
 {
     return dimensioned<scalar>
@@ -212,4 +217,5 @@ Foam::dimensioned<Foam::scalar> Foam::polynomial::dimIntegrate
     );
 }
 
+
 // ************************************************************************* //

src/OpenFOAM/primitives/functions/DataEntry/polynomial/polynomial.H

@@ -48,6 +48,7 @@ SourceFiles
 #include "DataEntry.H"
 #include "Tuple2.H"
 #include "dimensionSet.H"
+#include "DataEntryFwd.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -70,7 +71,7 @@ Ostream& operator<<
 
 class polynomial
 :
-    public DataEntry<scalar>
+    public scalarDataEntry
 {
     // Private data
 
@@ -107,9 +108,9 @@ public:
         polynomial(const polynomial& poly);
 
         //- Construct and return a clone
-        virtual tmp<DataEntry<scalar> > clone() const
+        virtual tmp<scalarDataEntry> clone() const
         {
-            return tmp<DataEntry<scalar> >(new polynomial(*this));
+            return tmp<scalarDataEntry>(new polynomial(*this));
         }
 
 

src/Pstream/mpi/UIPread.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -107,12 +107,7 @@ Foam::UIPstream::UIPstream
 
         if (!messageSize_)
         {
-            FatalErrorIn
+            setEof();
-            (
-                "UIPstream::UIPstream(const commsTypes, const int, "
-                "DynamicList<char>&, streamFormat, versionNumber)"
-            )   << "read failed"
-                << Foam::abort(FatalError);
         }
     }
 }
@@ -199,11 +194,7 @@ Foam::UIPstream::UIPstream(const int fromProcNo, PstreamBuffers& buffers)
 
         if (!messageSize_)
         {
-            FatalErrorIn
+            setEof();
-            (
-                "UIPstream::UIPstream(const int, PstreamBuffers&)"
-            )   << "read failed"
-                << Foam::abort(FatalError);
         }
     }
 }

src/combustionModels/FSD/FSD.C

@@ -272,7 +272,7 @@ void FSD<CombThermoType, ThermoType>::calculateSourceNorm()
 
     forAll(ft_, cellI)
     {
-        if(ft_[cellI] < ftStoich)
+        if (ft_[cellI] < ftStoich)
         {
             pc[cellI] = ft_[cellI]*(YprodTotal/ftStoich);
         }

src/combustionModels/Make/files

@@ -10,6 +10,7 @@ rhoCombustionModel/rhoCombustionModel/rhoCombustionModelNew.C
 rhoCombustionModel/rhoThermoCombustion/rhoThermoCombustion.C
 rhoCombustionModel/rhoChemistryCombustion/rhoChemistryCombustion.C
 
+diffusion/diffusions.C
 infinitelyFastChemistry/infinitelyFastChemistrys.C
 
 PaSR/PaSRs.C

src/combustionModels/PaSR/PaSR.C

@@ -138,7 +138,7 @@ void Foam::combustionModels::PaSR<Type>::correct()
 
 template<class Type>
 Foam::tmp<Foam::fvScalarMatrix>
-Foam::combustionModels::PaSR<Type>::R(const volScalarField& Y) const
+Foam::combustionModels::PaSR<Type>::R(volScalarField& Y) const
 {
     tmp<fvScalarMatrix> tSu(new fvScalarMatrix(Y, dimMass/dimTime));
 

src/combustionModels/PaSR/PaSR.H

@@ -103,7 +103,7 @@ public:
             virtual void correct();
 
             //- Fuel consumption rate matrix.
-            virtual tmp<fvScalarMatrix> R(const volScalarField& Y) const;
+            virtual tmp<fvScalarMatrix> R(volScalarField& Y) const;
 
             //- Heat release rate calculated from fuel consumption rate matrix
             virtual tmp<volScalarField> dQ() const;

src/combustionModels/combustionModel/combustionModel.H

@@ -135,7 +135,7 @@ public:
         virtual void correct() = 0;
 
         //- Fuel consumption rate matrix, i.e. source term for fuel equation
-        virtual tmp<fvScalarMatrix> R(const volScalarField& Y) const = 0;
+        virtual tmp<fvScalarMatrix> R(volScalarField& Y) const = 0;
 
         //- Heat release rate calculated from fuel consumption rate matrix
         virtual tmp<volScalarField> dQ() const = 0;

src/combustionModels/diffusion/diffusion.C

@@ -0,0 +1,106 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2012 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "diffusion.H"
+#include "fvcGrad.H"
+
+namespace Foam
+{
+namespace combustionModels
+{
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+template<class CombThermoType, class ThermoType>
+diffusion<CombThermoType, ThermoType>::diffusion
+(
+    const word& modelType, const fvMesh& mesh
+)
+:
+    singleStepCombustion<CombThermoType, ThermoType>(modelType, mesh),
+    C_(readScalar(this->coeffs().lookup("C"))),
+    oxidantName_(this->coeffs().template lookupOrDefault<word>("oxidant", "O2"))
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructors * * * * * * * * * * * * * * * //
+
+template<class CombThermoType, class ThermoType>
+diffusion<CombThermoType, ThermoType>::~diffusion()
+{}
+
+
+// * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
+
+template<class CombThermoType, class ThermoType>
+void diffusion<CombThermoType, ThermoType>::correct()
+{
+    this->wFuel_ ==
+        dimensionedScalar("zero", dimMass/pow3(dimLength)/dimTime, 0.0);
+
+    if (this->active())
+    {
+        this->singleMixturePtr_->fresCorrect();
+
+        const label fuelI = this->singleMixturePtr_->fuelIndex();
+
+        const volScalarField& YFuel =
+            this->thermoPtr_->composition().Y()[fuelI];
+
+        if (this->thermoPtr_->composition().contains(oxidantName_))
+        {
+            const volScalarField& YO2 =
+                this->thermoPtr_->composition().Y(oxidantName_);
+
+            this->wFuel_ ==
+                C_*this->turbulence().muEff()
+               *mag(fvc::grad(YFuel) & fvc::grad(YO2))
+               *pos(YFuel)*pos(YO2);
+        }
+    }
+}
+
+
+template<class CombThermoType, class ThermoType>
+bool diffusion<CombThermoType, ThermoType>::read()
+{
+    if (singleStepCombustion<CombThermoType, ThermoType>::read())
+    {
+        this->coeffs().lookup("C") >> C_ ;
+        this->coeffs().readIfPresent("oxidant", oxidantName_);
+        return true;
+    }
+    else
+    {
+        return false;
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace combustionModels
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

src/combustionModels/diffusion/diffusion.H

@@ -0,0 +1,122 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2012 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::combustionModels::diffusion
+
+Description
+    Simple diffusion-based combustion model based on the principle mixed is
+    burnt. Additional parameter C is used to distribute the heat release rate
+    in time.
+
+SourceFiles
+    diffusion.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef diffusion_H
+#define diffusion_H
+
+#include "singleStepCombustion.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace combustionModels
+{
+
+/*---------------------------------------------------------------------------*\
+                         Class diffusion Declaration
+\*---------------------------------------------------------------------------*/
+
+template<class CombThermoType, class ThermoType>
+class diffusion
+:
+    public singleStepCombustion<CombThermoType, ThermoType>
+{
+    // Private data
+
+        //- Model constant
+        scalar C_;
+
+        //- Name of oxidant - default is "O2"
+        word oxidantName_;
+
+
+    // Private Member Functions
+
+        //- Disallow copy construct
+        diffusion(const diffusion&);
+
+        //- Disallow default bitwise assignment
+        void operator=(const diffusion&);
+
+
+public:
+
+    //- Runtime type information
+    TypeName("diffusion");
+
+
+    // Constructors
+
+        //- Construct from components
+        diffusion(const word& modelType, const fvMesh& mesh);
+
+
+    //- Destructor
+    virtual ~diffusion();
+
+
+    // Member Functions
+
+        // Evolution
+
+            //- Correct combustion rate
+            virtual void correct();
+
+
+        // I-O
+
+            //- Update properties
+            virtual bool read();
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace combustionModels
+} // End namespace Foam
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+#   include "diffusion.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/combustionModels/diffusion/diffusions.C

@@ -0,0 +1,74 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2012 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "makeCombustionTypes.H"
+
+#include "thermoPhysicsTypes.H"
+#include "psiThermoCombustion.H"
+#include "rhoThermoCombustion.H"
+#include "diffusion.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace combustionModels
+{
+    makeCombustionTypesThermo
+    (
+        diffusion,
+        psiThermoCombustion,
+        gasThermoPhysics,
+        psiCombustionModel
+    );
+
+    makeCombustionTypesThermo
+    (
+        diffusion,
+        psiThermoCombustion,
+        constGasThermoPhysics,
+        psiCombustionModel
+    );
+
+    makeCombustionTypesThermo
+    (
+        diffusion,
+        rhoThermoCombustion,
+        gasThermoPhysics,
+        rhoCombustionModel
+    );
+
+    makeCombustionTypesThermo
+    (
+        diffusion,
+        rhoThermoCombustion,
+        constGasThermoPhysics,
+        rhoCombustionModel
+    );
+}
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

src/combustionModels/noCombustion/noCombustion.C

@@ -59,7 +59,7 @@ template<class CombThermoType>
 Foam::tmp<Foam::fvScalarMatrix>
 Foam::combustionModels::noCombustion<CombThermoType>::R
 (
-    const volScalarField& Y
+    volScalarField& Y
 ) const
 {
     tmp<fvScalarMatrix> tSu

src/combustionModels/noCombustion/noCombustion.H

@@ -84,8 +84,8 @@ public:
             //- Correct combustion rate
             virtual void correct();
 
-            //- Fuel consumption rate matrix.
+            //- Fuel consumption rate matrix
-            virtual tmp<fvScalarMatrix> R(const volScalarField& Y) const;
+            virtual tmp<fvScalarMatrix> R(volScalarField& Y) const;
 
             //- Heat release rate calculated from fuel consumption rate matrix
             virtual tmp<volScalarField> dQ() const;

src/combustionModels/singleStepCombustion/singleStepCombustion.C

@@ -53,7 +53,8 @@ singleStepCombustion<CombThermoType, ThermoType>::singleStepCombustion
         ),
         this->mesh(),
         dimensionedScalar("zero", dimMass/dimVolume/dimTime, 0.0)
-    )
+    ),
+    semiImplicit_(readBool(this->coeffs_.lookup("semiImplicit")))
 {
     if (isA<singleStepReactingMixture<ThermoType> >(this->thermo()))
     {
@@ -79,6 +80,15 @@ singleStepCombustion<CombThermoType, ThermoType>::singleStepCombustion
             << "Please select a thermo package based on "
             << "singleStepReactingMixture" << exit(FatalError);
     }
+
+    if (semiImplicit_)
+    {
+        Info<< "Combustion mode: semi-implicit" << endl;
+    }
+    else
+    {
+        Info<< "Combustion mode: explicit" << endl;
+    }
 }
 
 
@@ -95,17 +105,28 @@ template<class CombThermoType, class ThermoType>
 Foam::tmp<Foam::fvScalarMatrix>
 singleStepCombustion<CombThermoType, ThermoType>::R
 (
-    const volScalarField& Y
+    volScalarField& Y
 ) const
 {
     const label specieI = this->thermoPtr_->composition().species()[Y.name()];
 
-    const volScalarField wSpecie
+    volScalarField wSpecie
     (
         wFuel_*singleMixturePtr_->specieStoichCoeffs()[specieI]
     );
 
+    if (semiImplicit_)
+    {
+        const label fNorm = singleMixturePtr_->specieProd()[specieI];
+        const volScalarField fres(singleMixturePtr_->fres(specieI));
+        wSpecie /= max(fNorm*(Y - fres), 1e-2);
+
+        return -fNorm*wSpecie*fres + fNorm*fvm::Sp(wSpecie, Y);
+    }
+    else
+    {
         return wSpecie + fvm::Sp(0.0*wSpecie, Y);
+    }
 }
 
 
@@ -114,7 +135,8 @@ Foam::tmp<Foam::volScalarField>
 singleStepCombustion< CombThermoType, ThermoType>::Sh() const
 {
     const label fuelI = singleMixturePtr_->fuelIndex();
-    const volScalarField& YFuel = this->thermoPtr_->composition().Y(fuelI);
+    volScalarField& YFuel =
+        const_cast<volScalarField&>(this->thermoPtr_->composition().Y(fuelI));
 
     return -singleMixturePtr_->qFuel()*(R(YFuel) & YFuel);
 }

src/combustionModels/singleStepCombustion/singleStepCombustion.H

@@ -71,6 +71,9 @@ protected:
         //- Fuel consumption rate
         volScalarField wFuel_;
 
+        //- Semi-implicit (true) or explicit (false) treatment
+        bool semiImplicit_;
+
 
 public:
 
@@ -89,7 +92,7 @@ public:
         // Evolution
 
             //- Fuel consumption rate matrix
-            virtual tmp<fvScalarMatrix> R(const volScalarField& Y) const;
+            virtual tmp<fvScalarMatrix> R(volScalarField& Y) const;
 
             //- Heat release rate calculated from fuel consumption rate matrix
             virtual tmp<volScalarField> dQ() const;

src/fieldSources/basicSource/interRegionHeatTransferModel/interRegionHeatTransferModel/interRegionHeatTransferModel.C

@@ -68,17 +68,11 @@ void Foam::interRegionHeatTransferModel::check()
         }
     }
 
-    if(!secSourceFound)
+    if (!secSourceFound)
     {
         FatalErrorIn
         (
-            "constantHeatTransfer::interRegionHeatTransferModel"
+            "constantHeatTransfer::interRegionHeatTransferModel::check()"
-            "("
-            "   const word& name,"
-            "   const word& modelType,"
-            "   const dictionary& dict,"
-            "   const fvMesh& mesh"
-            ")"
         )   << "Secondary source name not found" << secondarySourceName_
             << " in region " << secondaryMesh.name()
             << nl
@@ -125,9 +119,12 @@ Foam::interRegionHeatTransferModel::interRegionHeatTransferModel
 }
 
 
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
 Foam::interRegionHeatTransferModel::~interRegionHeatTransferModel()
 {}
 
+
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
 
@@ -222,7 +219,7 @@ void Foam::interRegionHeatTransferModel::addSup
                     << endl;
             }
         }
-        else if(h.dimensions() == dimTemperature)
+        else if (h.dimensions() == dimTemperature)
         {
             eEqn += htc_*Tmapped - fvm::Sp(htc_, h);
 
@@ -269,4 +266,6 @@ bool Foam::interRegionHeatTransferModel::read(const dictionary& dict)
         return false;
     }
 }
+
+
 // ************************************************************************* //

src/finiteVolume/fields/fvPatchFields/basic/calculated/calculatedFvPatchField.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -98,20 +98,19 @@ calculatedFvPatchField<Type>::calculatedFvPatchField
 
 
 template<class Type>
-template<class Type2>
 tmp<fvPatchField<Type> > fvPatchField<Type>::NewCalculatedType
 (
-    const fvPatchField<Type2>& pf
+    const fvPatch& p
 )
 {
     typename patchConstructorTable::iterator patchTypeCstrIter =
-        patchConstructorTablePtr_->find(pf.patch().type());
+        patchConstructorTablePtr_->find(p.type());
 
     if (patchTypeCstrIter != patchConstructorTablePtr_->end())
     {
         return patchTypeCstrIter()
         (
-            pf.patch(),
+            p,
             DimensionedField<Type, volMesh>::null()
         );
     }
@@ -121,7 +120,7 @@ tmp<fvPatchField<Type> > fvPatchField<Type>::NewCalculatedType
         (
             new calculatedFvPatchField<Type>
             (
-                pf.patch(),
+                p,
                 DimensionedField<Type, volMesh>::null()
             )
         );
@@ -129,6 +128,17 @@ tmp<fvPatchField<Type> > fvPatchField<Type>::NewCalculatedType
 }
 
 
+template<class Type>
+template<class Type2>
+tmp<fvPatchField<Type> > fvPatchField<Type>::NewCalculatedType
+(
+    const fvPatchField<Type2>& pf
+)
+{
+    return NewCalculatedType(pf.patch());
+}
+
+
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
 template<class Type>

src/finiteVolume/fields/fvPatchFields/fvPatchField/fvPatchField.H

@@ -257,6 +257,13 @@ public:
             const dictionary&
         );
 
+        //- Return a pointer to a new calculatedFvPatchField created on
+        //  freestore without setting patchField values
+        static tmp<fvPatchField<Type> > NewCalculatedType
+        (
+            const fvPatch&
+        );
+
         //- Return a pointer to a new calculatedFvPatchField created on
         //  freestore without setting patchField values
         template<class Type2>

src/finiteVolume/fields/fvsPatchFields/basic/calculated/calculatedFvsPatchField.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -97,21 +97,20 @@ calculatedFvsPatchField<Type>::calculatedFvsPatchField
 
 
 template<class Type>
-template<class Type2>
 tmp<fvsPatchField<Type> > fvsPatchField<Type>::NewCalculatedType
 (
-    const fvsPatchField<Type2>& pf
+    const fvPatch& p
 )
 {
     typename patchConstructorTable::iterator patchTypeCstrIter =
-        patchConstructorTablePtr_->find(pf.patch().type());
+        patchConstructorTablePtr_->find(p.type());
 
     if (patchTypeCstrIter != patchConstructorTablePtr_->end())
     {
         return patchTypeCstrIter()
         (
-            pf.patch(),
+            p,
-            Field<Type>::null()
+            DimensionedField<Type, surfaceMesh>::null()
         );
     }
     else
@@ -120,14 +119,25 @@ tmp<fvsPatchField<Type> > fvsPatchField<Type>::NewCalculatedType
         (
             new calculatedFvsPatchField<Type>
             (
-                pf.patch(),
+                p,
-                Field<Type>::null()
+                DimensionedField<Type, surfaceMesh>::null()
             )
         );
     }
 }
 
 
+template<class Type>
+template<class Type2>
+tmp<fvsPatchField<Type> > fvsPatchField<Type>::NewCalculatedType
+(
+    const fvsPatchField<Type2>& pf
+)
+{
+    return NewCalculatedType(pf.patch());
+}
+
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 } // End namespace Foam

src/finiteVolume/fields/fvsPatchFields/fvsPatchField/fvsPatchField.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -247,6 +247,13 @@ public:
             const dictionary&
         );
 
+        //- Return a pointer to a new calculatedFvsPatchField created on
+        //  freestore without setting patchField values
+        static tmp<fvsPatchField<Type> > NewCalculatedType
+        (
+            const fvPatch&
+        );
+
         //- Return a pointer to a new calculatedFvsPatchField created on
         //  freestore without setting patchField values
         template<class Type2>

src/finiteVolume/interpolation/surfaceInterpolation/surfaceInterpolation/surfaceInterpolate.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -302,6 +302,43 @@ interpolate
 }
 
 
+template<class Type>
+tmp<FieldField<fvsPatchField, Type> > interpolate
+(
+    const FieldField<fvPatchField, Type>& fvpff
+)
+{
+    FieldField<fvsPatchField, Type>* fvspffPtr
+    (
+        new FieldField<fvsPatchField, Type>(fvpff.size())
+    );
+
+    forAll(*fvspffPtr, patchi)
+    {
+        fvspffPtr->set
+        (
+            patchi,
+            fvsPatchField<Type>::NewCalculatedType(fvpff[patchi].patch()).ptr()
+        );
+        (*fvspffPtr)[patchi] = fvpff[patchi];
+    }
+
+    return tmp<FieldField<fvsPatchField, Type> >(fvspffPtr);
+}
+
+
+template<class Type>
+tmp<FieldField<fvsPatchField, Type> > interpolate
+(
+    const tmp<FieldField<fvPatchField, Type> >& tfvpff
+)
+{
+    tmp<FieldField<fvsPatchField, Type> > tfvspff = interpolate(tfvpff());
+    tfvpff.clear();
+    return tfvspff;
+}
+
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 } // End namespace fvc

src/finiteVolume/interpolation/surfaceInterpolation/surfaceInterpolation/surfaceInterpolate.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -163,12 +163,26 @@ namespace fvc
         const tmp<GeometricField<Type, fvPatchField, volMesh> >& tvf
     );
 
-    //- Interpolate tmp field onto faces using 'interpolate(\<name\>)'
+    //- Interpolate field onto faces using 'interpolate(\<name\>)'
     template<class Type>
     static tmp<GeometricField<Type, fvsPatchField, surfaceMesh> > interpolate
     (
         const GeometricField<Type, fvPatchField, volMesh>& tvf
     );
+
+    //- Interpolate boundary field onto faces (simply a type conversion)
+    template<class Type>
+    static tmp<FieldField<fvsPatchField, Type> > interpolate
+    (
+        const FieldField<fvPatchField, Type>& fvpff
+    );
+
+    //- Interpolate boundary field onto faces (simply a type conversion)
+    template<class Type>
+    static tmp<FieldField<fvsPatchField, Type> > interpolate
+    (
+        const tmp<FieldField<fvPatchField, Type> >& tfvpff
+    );
 }
 
 

src/finiteVolume/interpolation/volPointInterpolation/volPointInterpolate.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -51,8 +51,8 @@ void volPointInterpolation::syncUntransformedData
     const indirectPrimitivePatch& cpp = gmd.coupledPatch();
     const labelList& meshPoints = cpp.meshPoints();
 
-    const mapDistribute& slavesMap = gmd.globalPointSlavesMap();
+    const mapDistribute& slavesMap = gmd.globalCoPointSlavesMap();
-    const labelListList& slaves = gmd.globalPointSlaves();
+    const labelListList& slaves = gmd.globalCoPointSlaves();
 
     List<Type> elems(slavesMap.constructSize());
     forAll(meshPoints, i)
@@ -105,8 +105,8 @@ void volPointInterpolation::pushUntransformedData
     const indirectPrimitivePatch& cpp = gmd.coupledPatch();
     const labelList& meshPoints = cpp.meshPoints();
 
-    const mapDistribute& slavesMap = gmd.globalPointSlavesMap();
+    const mapDistribute& slavesMap = gmd.globalCoPointSlavesMap();
-    const labelListList& slaves = gmd.globalPointSlaves();
+    const labelListList& slaves = gmd.globalCoPointSlaves();
 
     List<Type> elems(slavesMap.constructSize());
     forAll(meshPoints, i)
@@ -155,14 +155,14 @@ void volPointInterpolation::addSeparated
             refCast<coupledPointPatchField<Type> >
                 (pf.boundaryField()[patchI]).initSwapAddSeparated
                 (
-                    Pstream::blocking,  //Pstream::nonBlocking,
+                    Pstream::nonBlocking,
                     pf.internalField()
                 );
         }
     }
 
     // Block for any outstanding requests
-    //Pstream::waitRequests();
+    Pstream::waitRequests();
 
     forAll(pf.boundaryField(), patchI)
     {
@@ -171,7 +171,7 @@ void volPointInterpolation::addSeparated
             refCast<coupledPointPatchField<Type> >
                 (pf.boundaryField()[patchI]).swapAddSeparated
                 (
-                    Pstream::blocking,  //Pstream::nonBlocking,
+                    Pstream::nonBlocking,
                     pf.internalField()
                 );
         }
@@ -306,7 +306,6 @@ void volPointInterpolation::interpolateBoundaryField
     }
 
     // Sum collocated contributions
-    //mesh().globalData().syncPointData(pfi, plusEqOp<Type>());
     syncUntransformedData(pfi, plusEqOp<Type>());
 
     // And add separated contributions
@@ -314,9 +313,7 @@ void volPointInterpolation::interpolateBoundaryField
 
     // Push master data to slaves. It is possible (not sure how often) for
     // a coupled point to have its master on a different patch so
-    // to make sure just push master data to slaves. Reuse the syncPointData
+    // to make sure just push master data to slaves.
-    // structure.
-    //mesh().globalData().syncPointData(pfi, nopEqOp<Type>());
     pushUntransformedData(pfi);
 
 

src/finiteVolume/interpolation/volPointInterpolation/volPointInterpolation.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -113,10 +113,8 @@ void volPointInterpolation::calcBoundaryAddressing()
     {
         boolList oldData(isPatchPoint_);
 
-        //mesh().globalData().syncPointData(isPatchPoint_, orEqOp<bool>());
         syncUntransformedData(isPatchPoint_, orEqOp<bool>());
 
-
         forAll(isPatchPoint_, pointI)
         {
             if (isPatchPoint_[pointI] != oldData[pointI])
@@ -272,7 +270,7 @@ void volPointInterpolation::makeWeights()
     makeInternalWeights(sumWeights);
 
 
-    // Create boundary weights; add to sumWeights
+    // Create boundary weights; override sumWeights
     makeBoundaryWeights(sumWeights);
 
 
@@ -292,7 +290,6 @@ void volPointInterpolation::makeWeights()
 
 
     // Sum collocated contributions
-    //mesh().globalData().syncPointData(sumWeights, plusEqOp<scalar>());
     syncUntransformedData(sumWeights, plusEqOp<scalar>());
 
     // And add separated contributions
@@ -302,7 +299,6 @@ void volPointInterpolation::makeWeights()
     // a coupled point to have its master on a different patch so
     // to make sure just push master data to slaves. Reuse the syncPointData
     // structure.
-    //mesh().globalData().syncPointData(sumWeights, nopEqOp<scalar>());
     pushUntransformedData(sumWeights);
 
 

src/lagrangian/coalCombustion/Make/files

@@ -1,4 +1,5 @@
 /* Coal parcel and sub-models */
 coalParcel/makeCoalParcelSubmodels.C
+coalCloudList/coalCloudList.C
 
 LIB = $(FOAM_LIBBIN)/libcoalCombustion

src/lagrangian/coalCombustion/coalCloudList/coalCloudList.C

@@ -0,0 +1,92 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2012 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "coalCloudList.H"
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::coalCloudList::coalCloudList
+(
+    const volScalarField& rho,
+    const volVectorField& U,
+    const dimensionedVector& g,
+    const SLGThermo& slgThermo
+)
+:
+    PtrList<coalCloud>(),
+    mesh_(rho.mesh())
+{
+    IOdictionary props
+    (
+        IOobject
+        (
+            "coalCloudList",
+            mesh_.time().constant(),
+            mesh_,
+            IOobject::MUST_READ
+        )
+    );
+
+    const wordHashSet cloudNames(wordList(props.lookup("clouds")));
+
+    setSize(cloudNames.size());
+
+    label i = 0;
+    forAllConstIter(wordHashSet, cloudNames, iter)
+    {
+        const word& name = iter.key();
+
+        Info<< "creating cloud: " << name << endl;
+
+        set
+        (
+            i++,
+            new coalCloud
+            (
+                name,
+                rho,
+                U,
+                g,
+                slgThermo
+            )
+        );
+
+        Info<< endl;
+    }
+}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+void Foam::coalCloudList::evolve()
+{
+    forAll(*this, i)
+    {
+        operator[](i).evolve();
+    }
+}
+
+
+// ************************************************************************* //

src/lagrangian/coalCombustion/coalCloudList/coalCloudList.H

@@ -0,0 +1,126 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2012 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef coalCloudList_H
+#define coalCloudList_H
+
+#include "coalCloud.H"
+#include "volFieldsFwd.H"
+#include "fvMatricesFwd.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                       Class coalCloudList Declaration
+\*---------------------------------------------------------------------------*/
+
+class coalCloudList
+:
+    public PtrList<coalCloud>
+{
+private:
+
+    //- Reference to the mesh
+    const fvMesh& mesh_;
+
+
+public:
+
+    // Constructor
+
+        coalCloudList
+        (
+            const volScalarField& rho,
+            const volVectorField& U,
+            const dimensionedVector& g,
+            const SLGThermo& slgThermo
+        );
+
+
+    // Member Functions
+
+        // Evolution
+
+            //- Evolve the cloud collection
+            void evolve();
+
+
+        // Source terms
+
+            //- Return const reference to momentum source
+            inline tmp<DimensionedField<vector, volMesh> > UTrans() const;
+
+            //- Return tmp momentum source term
+            inline tmp<fvVectorMatrix> SU(volVectorField& U) const;
+
+            //- Sensible enthalpy transfer [J/kg]
+            inline tmp<DimensionedField<scalar, volMesh> > hsTrans() const;
+
+            //- Return sensible enthalpy source term [J/kg/m3/s]
+            inline tmp<fvScalarMatrix> Sh(volScalarField& hs) const;
+
+
+            //- Return mass source term for specie i - specie eqn
+            inline tmp<fvScalarMatrix> SYi
+            (
+                const label i,
+                volScalarField& Yi
+            ) const;
+
+            //- Return total mass transfer [kg/m3]
+            inline tmp<DimensionedField<scalar, volMesh> > rhoTrans() const;
+
+            //- Return tmp total mass source for carrier phase
+            //  - fully explicit
+            inline tmp<DimensionedField<scalar, volMesh> > Srho() const;
+
+            //- Return tmp total mass source for carrier phase specie i
+            //  - fully explicit
+            inline tmp<DimensionedField<scalar, volMesh> > Srho
+            (
+                const label i
+            ) const;
+
+            //- Return total mass source term [kg/m3/s]
+            inline tmp<fvScalarMatrix> Srho(volScalarField& rho) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#include "coalCloudListI.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/lagrangian/coalCombustion/coalCloudList/coalCloudListI.H

@@ -0,0 +1,262 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2012 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "fvMatrices.H"
+#include "volFields.H"
+#include "DimensionedField.H"
+
+Foam::tmp<Foam::DimensionedField<Foam::vector, Foam::volMesh> >
+Foam::coalCloudList::UTrans() const
+{
+    tmp<DimensionedField<vector, volMesh> > tfld
+    (
+        new DimensionedField<vector, volMesh>
+        (
+            IOobject
+            (
+                "UTransEff",
+                mesh_.time().timeName(),
+                mesh_,
+                IOobject::NO_READ,
+                IOobject::NO_WRITE
+            ),
+            mesh_,
+            dimensionedVector("zero", dimMass*dimVelocity, vector::zero)
+        )
+    );
+
+    DimensionedField<vector, volMesh>& fld = tfld();
+
+    forAll(*this, i)
+    {
+        fld += operator[](i).UTrans();
+    }
+
+    return tfld;
+}
+
+
+Foam::tmp<Foam::fvVectorMatrix> Foam::coalCloudList::SU
+(
+    volVectorField& U
+) const
+{
+    tmp<fvVectorMatrix> tfvm(new fvVectorMatrix(U, dimForce));
+    fvVectorMatrix& fvm = tfvm();
+
+    forAll(*this, i)
+    {
+        fvm += operator[](i).SU(U);
+    }
+
+    return tfvm;
+}
+
+
+Foam::tmp<Foam::DimensionedField<Foam::scalar, Foam::volMesh> >
+Foam::coalCloudList::hsTrans() const
+{
+    tmp<DimensionedField<scalar, volMesh> > tfld
+    (
+        new DimensionedField<scalar, volMesh>
+        (
+            IOobject
+            (
+                "hsTransEff",
+                mesh_.time().timeName(),
+                mesh_,
+                IOobject::NO_READ,
+                IOobject::NO_WRITE
+            ),
+            mesh_,
+            dimensionedScalar("zero", dimEnergy, 0.0)
+        )
+    );
+
+    DimensionedField<scalar, volMesh>& fld = tfld();
+
+    forAll(*this, i)
+    {
+        fld += operator[](i).hsTrans();
+    }
+
+    return tfld;
+}
+
+
+Foam::tmp<Foam::fvScalarMatrix> Foam::coalCloudList::Sh
+(
+    volScalarField& hs
+) const
+{
+    tmp<fvScalarMatrix> tfvm(new fvScalarMatrix(hs, dimEnergy/dimTime));
+    fvScalarMatrix& fvm = tfvm();
+
+    forAll(*this, i)
+    {
+        fvm += operator[](i).Sh(hs);
+    }
+
+    return tfvm;
+}
+
+
+Foam::tmp<Foam::fvScalarMatrix> Foam::coalCloudList::SYi
+(
+    const label ii,
+    volScalarField& Yi
+) const
+{
+    tmp<fvScalarMatrix> tfvm(new fvScalarMatrix(Yi, dimMass/dimTime));
+    fvScalarMatrix& fvm = tfvm();
+
+    forAll(*this, i)
+    {
+        fvm += operator[](i).SYi(ii, Yi);
+    }
+
+    return tfvm;
+}
+
+
+Foam::tmp<Foam::DimensionedField<Foam::scalar, Foam::volMesh> >
+Foam::coalCloudList::rhoTrans() const
+{
+    tmp<DimensionedField<scalar, volMesh> > tfld
+    (
+        new DimensionedField<scalar, volMesh>
+        (
+            IOobject
+            (
+                "rhoTransEff",
+                mesh_.time().timeName(),
+                mesh_,
+                IOobject::NO_READ,
+                IOobject::NO_WRITE
+            ),
+            mesh_,
+            dimensionedScalar("zero", dimMass, 0.0)
+        )
+    );
+
+    DimensionedField<scalar, volMesh>& fld = tfld();
+
+    forAll(*this, i)
+    {
+        forAll(operator[](i).rhoTrans(), j)
+        {
+            fld += operator[](i).rhoTrans()[j];
+        }
+    }
+
+    return tfld;
+}
+
+
+
+
+Foam::tmp<Foam::DimensionedField<Foam::scalar, Foam::volMesh> >
+Foam::coalCloudList::Srho() const
+{
+    tmp<DimensionedField<scalar, volMesh> > tfld
+    (
+        new DimensionedField<scalar, volMesh>
+        (
+            IOobject
+            (
+                "rhoTransEff",
+                mesh_.time().timeName(),
+                mesh_,
+                IOobject::NO_READ,
+                IOobject::NO_WRITE
+            ),
+            mesh_,
+            dimensionedScalar("zero", dimDensity/dimTime, 0.0)
+        )
+    );
+
+    DimensionedField<scalar, volMesh>& fld = tfld();
+
+    forAll(*this, i)
+    {
+        fld += operator[](i).Srho();
+    }
+
+    return tfld;
+}
+
+
+Foam::tmp<Foam::DimensionedField<Foam::scalar, Foam::volMesh> >
+Foam::coalCloudList::Srho
+(
+    const label i
+) const
+{
+    tmp<DimensionedField<scalar, volMesh> > tfld
+    (
+        new DimensionedField<scalar, volMesh>
+        (
+            IOobject
+            (
+                "rhoTransEff",
+                mesh_.time().timeName(),
+                mesh_,
+                IOobject::NO_READ,
+                IOobject::NO_WRITE
+            ),
+            mesh_,
+            dimensionedScalar("zero", dimDensity/dimTime, 0.0)
+        )
+    );
+
+    DimensionedField<scalar, volMesh>& fld = tfld();
+
+    forAll(*this, j)
+    {
+        fld += operator[](j).Srho(i);
+    }
+
+    return tfld;
+}
+
+
+Foam::tmp<Foam::fvScalarMatrix> Foam::coalCloudList::Srho
+(
+    volScalarField& rho
+) const
+{
+    tmp<fvScalarMatrix> tfvm(new fvScalarMatrix(rho, dimMass/dimTime));
+    fvScalarMatrix& fvm = tfvm();
+
+    forAll(*this, i)
+    {
+        fvm += operator[](i).Srho(rho);
+    }
+
+    return tfvm;
+}
+
+
+// ************************************************************************* //

src/lagrangian/coalCombustion/include/makeCoalParcelSurfaceReactionModels.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -31,6 +31,7 @@ License
 #include "NoSurfaceReaction.H"
 #include "COxidationDiffusionLimitedRate.H"
 #include "COxidationKineticDiffusionLimitedRate.H"
+#include "COxidationHurtMitchell.H"
 #include "COxidationMurphyShaddix.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -43,6 +44,7 @@ License
         COxidationKineticDiffusionLimitedRate,                                \
         CloudType                                                             \
     );                                                                        \
+    makeSurfaceReactionModelType(COxidationHurtMitchell, CloudType);          \
     makeSurfaceReactionModelType(COxidationMurphyShaddix, CloudType);
 
 

src/lagrangian/coalCombustion/submodels/surfaceReactionModel/COxidationHurtMitchell/COxidationHurtMitchell.C

@@ -0,0 +1,206 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2012 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "COxidationHurtMitchell.H"
+#include "mathematicalConstants.H"
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+template<class CloudType>
+Foam::COxidationHurtMitchell<CloudType>::COxidationHurtMitchell
+(
+    const dictionary& dict,
+    CloudType& owner
+)
+:
+    SurfaceReactionModel<CloudType>(dict, owner, typeName),
+    Sb_(readScalar(this->coeffDict().lookup("Sb"))),
+    CsLocalId_(-1),
+    ashLocalId_(-1),
+    O2GlobalId_(owner.composition().globalCarrierId("O2")),
+    CO2GlobalId_(owner.composition().globalCarrierId("CO2")),
+    WC_(0.0),
+    WO2_(0.0),
+    HcCO2_(0.0),
+    heatOfReaction_(-1.0)
+{
+    // Determine Cs and ash ids
+    label idSolid = owner.composition().idSolid();
+    CsLocalId_ = owner.composition().localId(idSolid, "C");
+    ashLocalId_ = owner.composition().localId(idSolid, "ash", true);
+
+    // Set local copies of thermo properties
+    WO2_ = owner.thermo().carrier().W(O2GlobalId_);
+    const scalar WCO2 = owner.thermo().carrier().W(CO2GlobalId_);
+    WC_ = WCO2 - WO2_;
+
+    HcCO2_ = owner.thermo().carrier().Hc(CO2GlobalId_);
+
+    const scalar YCloc = owner.composition().Y0(idSolid)[CsLocalId_];
+    const scalar YSolidTot = owner.composition().YMixture0()[idSolid];
+    Info<< "    C(s): particle mass fraction = " << YCloc*YSolidTot << endl;
+
+    if (this->coeffDict().readIfPresent("heatOfReaction", heatOfReaction_))
+    {
+        Info<< "    Using user specified heat of reaction: "
+            << heatOfReaction_ << " [J/kg]" << endl;
+    }
+}
+
+
+template<class CloudType>
+Foam::COxidationHurtMitchell<CloudType>::COxidationHurtMitchell
+(
+    const COxidationHurtMitchell<CloudType>& srm
+)
+:
+    SurfaceReactionModel<CloudType>(srm),
+    Sb_(srm.Sb_),
+    CsLocalId_(srm.CsLocalId_),
+    ashLocalId_(srm.ashLocalId_),
+    O2GlobalId_(srm.O2GlobalId_),
+    CO2GlobalId_(srm.CO2GlobalId_),
+    WC_(srm.WC_),
+    WO2_(srm.WO2_),
+    HcCO2_(srm.HcCO2_),
+    heatOfReaction_(srm.heatOfReaction_)
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+template<class CloudType>
+Foam::COxidationHurtMitchell<CloudType>::~COxidationHurtMitchell()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+template<class CloudType>
+Foam::scalar Foam::COxidationHurtMitchell<CloudType>::calculate
+(
+    const scalar dt,
+    const label cellI,
+    const scalar d,
+    const scalar T,
+    const scalar Tc,
+    const scalar pc,
+    const scalar rhoc,
+    const scalar mass,
+    const scalarField& YGas,
+    const scalarField& YLiquid,
+    const scalarField& YSolid,
+    const scalarField& YMixture,
+    const scalar N,
+    scalarField& dMassGas,
+    scalarField& dMassLiquid,
+    scalarField& dMassSolid,
+    scalarField& dMassSRCarrier
+) const
+{
+    const label idGas = CloudType::parcelType::GAS;
+    const label idSolid = CloudType::parcelType::SLD;
+    const scalar Ychar = YMixture[idSolid]*YSolid[CsLocalId_];
+
+    // Surface combustion until combustible fraction is consumed
+    if (Ychar < SMALL)
+    {
+        return 0.0;
+    }
+
+    const SLGThermo& thermo = this->owner().thermo();
+
+    // Local mass fraction of O2 in the carrier phase
+    const scalar YO2 = thermo.carrier().Y(O2GlobalId_)[cellI];
+
+    // No combustion if no oxygen present
+    if (YO2 < SMALL)
+    {
+        return 0.0;
+    }
+
+    // Conversion from [g/cm^2) to [kg/m^2]
+    const scalar convSI = 1000.0/10000.0;
+
+    // Universal gas constant in [kcal/mol/K]
+    const scalar RRcal = 1985.877534;
+
+    // Dry mass fraction
+    scalar Ydaf = YMixture[idGas] + YMixture[idSolid];
+    if (ashLocalId_ != -1)
+    {
+        Ydaf -= YMixture[idSolid]*YSolid[ashLocalId_];
+    }
+
+    // Char percentage
+    const scalar charPrc = Ychar/Ydaf*100.0;
+
+    // Particle surface area
+    const scalar Ap = constant::mathematical::pi*sqr(d);
+
+    // Far field partial pressure O2 [Pa]
+    // Note: Should really use the surface partial pressure
+    const scalar ppO2 = max(0.0, rhoc*YO2/WO2_*specie::RR*Tc);
+
+    // Activation energy [kcal/mol]
+    const scalar E = -5.94 + 0.355*charPrc;
+
+    // Pre-exponential factor [g/(cm^2.s.atm^0.5)]
+    const scalar lnK1750 = 2.8 - 0.0758*charPrc;
+    const scalar A = exp(lnK1750 + E/RRcal/1750.0);
+
+    // Kinetic rate of char oxidation [g/(cm^2.s.atm^0.5)]
+    const scalar Rk = A*exp(-E/(RRcal*T));
+
+    // Molar reaction rate per unit surface area [kmol/(m^2.s)]
+    const scalar qCsLim = mass*Ychar/(WC_*Ap*dt);
+    const scalar qCs = min(convSI*Rk*Foam::sqrt(ppO2/101325.0), qCsLim);
+
+    // Calculate the number of molar units reacted [kmol]
+    const scalar dOmega = qCs*Ap*dt;
+
+    // Add to carrier phase mass transfer
+    dMassSRCarrier[O2GlobalId_] += -dOmega*Sb_*WO2_;
+    dMassSRCarrier[CO2GlobalId_] += dOmega*(WC_ + Sb_*WO2_);
+
+    // Add to particle mass transfer
+    dMassSolid[CsLocalId_] += dOmega*WC_;
+
+
+    // Return the heat of reaction [J]
+    // note: carrier sensible enthalpy exchange handled via change in mass
+    if (heatOfReaction_ < 0)
+    {
+        const scalar HsC = thermo.solids().properties()[CsLocalId_].Hs(T);
+        return dOmega*(WC_*HsC - (WC_ + Sb_*WO2_)*HcCO2_);
+    }
+    else
+    {
+        return dOmega*WC_*heatOfReaction_;
+    }
+}
+
+
+// ************************************************************************* //

src/lagrangian/coalCombustion/submodels/surfaceReactionModel/COxidationHurtMitchell/COxidationHurtMitchell.H

@@ -0,0 +1,184 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2012 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    COxidationHurtMitchell
+
+Description
+    Char oxidation model given by Hurt and Mitchell:
+
+    Based on the reference:
+        Hurt R. and Mitchell R., "Unified high-temperature char combustion
+        kinetics for a suite of coals of various rank", 24th Symposium in
+        Combustion, The Combustion Institute, 1992, p 1243-1250
+
+    Model specifies the rate of char combustion.
+
+        C(s) + Sb*O2 -> CO2
+
+    where Sb is the stoichiometry of the reaction
+
+    Model validity:
+        Gas temperature: Tc > 1500 K
+        Particle sizes:  75 um -> 200 um
+        Pox > 0.3 atm        
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef COxidationHurtMitchell_H
+#define COxidationHurtMitchell_H
+
+#include "SurfaceReactionModel.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// Forward class declarations
+template<class CloudType>
+class COxidationHurtMitchell;
+
+/*---------------------------------------------------------------------------*\
+                  Class COxidationHurtMitchell Declaration
+\*---------------------------------------------------------------------------*/
+
+template<class CloudType>
+class COxidationHurtMitchell
+:
+    public SurfaceReactionModel<CloudType>
+{
+    // Private data
+
+        // Model constants
+
+            //- Stoichiometry of reaction
+            const scalar Sb_;
+
+
+        // Addressing
+
+            //- Cs position in global/local lists
+            label CsLocalId_;
+
+            //- Ash position in global/local lists
+            label ashLocalId_;
+
+            //- O2 position in global list
+            label O2GlobalId_;
+
+            //- CO2 positions in global list
+            label CO2GlobalId_;
+
+
+        // Local copies of thermo properties
+
+            //- Molecular weight of C [kg/kmol]
+            scalar WC_;
+
+            //- Molecular weight of O2 [kg/kmol]
+            scalar WO2_;
+
+            //- Formation enthalpy for CO2 [J/kg]
+            scalar HcCO2_;
+
+            //- Heat of reaction [J/kg] (optional)
+            scalar heatOfReaction_;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("COxidationHurtMitchell");
+
+
+    // Constructors
+
+        //- Construct from dictionary
+        COxidationHurtMitchell
+        (
+            const dictionary& dict,
+            CloudType& owner
+        );
+
+        //- Construct copy
+        COxidationHurtMitchell
+        (
+            const COxidationHurtMitchell<CloudType>& srm
+        );
+
+        //- Construct and return a clone
+        virtual autoPtr<SurfaceReactionModel<CloudType> > clone() const
+        {
+            return autoPtr<SurfaceReactionModel<CloudType> >
+            (
+                new COxidationHurtMitchell<CloudType>(*this)
+            );
+        }
+
+
+    //- Destructor
+    virtual ~COxidationHurtMitchell();
+
+
+    // Member Functions
+
+        //- Update surface reactions
+        virtual scalar calculate
+        (
+            const scalar dt,
+            const label cellI,
+            const scalar d,
+            const scalar T,
+            const scalar Tc,
+            const scalar pc,
+            const scalar rhoc,
+            const scalar mass,
+            const scalarField& YGas,
+            const scalarField& YLiquid,
+            const scalarField& YSolid,
+            const scalarField& YMixture,
+            const scalar N,
+            scalarField& dMassGas,
+            scalarField& dMassLiquid,
+            scalarField& dMassSolid,
+            scalarField& dMassSRCarrier
+        ) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+    #include "COxidationHurtMitchell.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/lagrangian/dsmc/submodels/BinaryCollisionModel/NoBinaryCollision/NoBinaryCollision.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -68,15 +68,13 @@ Foam::scalar Foam::NoBinaryCollision<CloudType>::sigmaTcR
     (
         "Foam::scalar Foam::NoBinaryCollision<CloudType>::sigmaTcR"
         "("
-            "label typeIdP,"
+            "const typename CloudType::parcelType&, "
-            "label typeIdQ,"
+            "const typename CloudType::parcelType"
-            "const vector& UP,"
-            "const vector& UQ"
         ") const"
     )
         << "sigmaTcR called on NoBinaryCollision model, this should "
         << "not happen, this is not an actual model." << nl
-        << "Enclose calls to sigmaTcR within a if(binaryCollision().active()) "
+        << "Enclose calls to sigmaTcR within a if (binaryCollision().active()) "
         << " check."
         << abort(FatalError);
 

src/lagrangian/intermediate/clouds/Templates/KinematicCloud/KinematicCloud.C

@@ -670,6 +670,159 @@ void  Foam::KinematicCloud<CloudType>::motion(TrackData& td)
 }
 
 
+template<class CloudType>
+void Foam::KinematicCloud<CloudType>::patchData
+(
+    const parcelType& p,
+    const polyPatch& pp,
+    const scalar trackFraction,
+    const tetIndices& tetIs,
+    vector& nw,
+    vector& Up
+) const
+{
+    label patchI = pp.index();
+    label patchFaceI = pp.whichFace(p.face());
+
+    vector n = tetIs.faceTri(mesh_).normal();
+    n /= mag(n);
+
+    vector U = U_.boundaryField()[patchI][patchFaceI];
+
+    // Unless the face is rotating, the required normal is n;
+    nw = n;
+
+    if (!mesh_.moving())
+    {
+        // Only wall patches may have a non-zero wall velocity from
+        // the velocity field when the mesh is not moving.
+
+        if (isA<wallPolyPatch>(pp))
+        {
+            Up = U;
+        }
+        else
+        {
+            Up = vector::zero;
+        }
+    }
+    else
+    {
+        vector U00 = U_.oldTime().boundaryField()[patchI][patchFaceI];
+
+        vector n00 = tetIs.oldFaceTri(mesh_).normal();
+
+        // Difference in normal over timestep
+        vector dn = vector::zero;
+
+        if (mag(n00) > SMALL)
+        {
+            // If the old normal is zero (for example in layer
+            // addition) then use the current normal, meaning that the
+            // motion can only be translational, and dn remains zero,
+            // otherwise, calculate dn:
+
+            n00 /= mag(n00);
+
+            dn = n - n00;
+        }
+
+        // Total fraction through the timestep of the motion,
+        // including stepFraction before the current tracking step
+        // and the current trackFraction
+        // i.e.
+        // let s = stepFraction, t = trackFraction
+        // Motion of x in time:
+        // |-----------------|---------|---------|
+        // x00               x0        xi        x
+        //
+        // where xi is the correct value of x at the required
+        // tracking instant.
+        //
+        // x0 = x00 + s*(x - x00) = s*x + (1 - s)*x00
+        //
+        // i.e. the motion covered by previous tracking portions
+        // within this timestep, and
+        //
+        // xi = x0 + t*(x - x0)
+        //    = t*x + (1 - t)*x0
+        //    = t*x + (1 - t)*(s*x + (1 - s)*x00)
+        //    = (s + t - s*t)*x + (1 - (s + t - s*t))*x00
+        //
+        // let m = (s + t - s*t)
+        //
+        // xi = m*x + (1 - m)*x00 = x00 + m*(x - x00);
+        //
+        // In the same form as before.
+
+        scalar m =
+            p.stepFraction()
+          + trackFraction
+          - (p.stepFraction()*trackFraction);
+
+        // When the mesh is moving, the velocity field on wall patches
+        // will contain the velocity associated with the motion of the
+        // mesh, in which case it is interpolated in time using m.
+        // For other patches the face velocity will need to be
+        // reconstructed from the face centre motion.
+
+        const vector& Cf = mesh_.faceCentres()[p.face()];
+
+        vector Cf00 = mesh_.faces()[p.face()].centre(mesh_.oldPoints());
+
+        if (isA<wallPolyPatch>(pp))
+        {
+            Up = U00 + m*(U - U00);
+        }
+        else
+        {
+            Up = (Cf - Cf00)/mesh_.time().deltaTValue();
+        }
+
+        if (mag(dn) > SMALL)
+        {
+            // Rotational motion, nw requires interpolation and a
+            // rotational velocity around face centre correction to Up
+            // is required.
+
+            nw = n00 + m*dn;
+
+            // Cf at tracking instant
+            vector Cfi = Cf00 + m*(Cf - Cf00);
+
+            // Normal vector cross product
+            vector omega = (n00 ^ n);
+
+            scalar magOmega = mag(omega);
+
+            // magOmega = sin(angle between unit normals)
+            // Normalise omega vector by magOmega, then multiply by
+            // angle/dt to give the correct angular velocity vector.
+            omega *= Foam::asin(magOmega)/(magOmega*mesh_.time().deltaTValue());
+
+            // Project position onto face and calculate this position
+            // relative to the face centre.
+            vector facePos =
+                p.position()
+              - ((p.position() - Cfi) & nw)*nw
+              - Cfi;
+
+            Up += (omega ^ facePos);
+        }
+
+        // No further action is required if the motion is
+        // translational only, nw and Up have already been set.
+    }
+}
+
+
+template<class CloudType>
+void Foam::KinematicCloud<CloudType>::updateMesh()
+{
+    injectors_.updateMesh();
+}
+
+
 template<class CloudType>
 void Foam::KinematicCloud<CloudType>::autoMap(const mapPolyMesh& mapper)
 {
@@ -678,6 +831,8 @@ void Foam::KinematicCloud<CloudType>::autoMap(const mapPolyMesh& mapper)
     tdType td(*this);
 
     Cloud<parcelType>::template autoMap<tdType>(td, mapper);
+
+    updateMesh();
 }
 
 

src/lagrangian/intermediate/clouds/Templates/KinematicCloud/KinematicCloud.H

@@ -554,9 +554,24 @@ public:
             template<class TrackData>
             void motion(TrackData& td);
 
+            //- Calculate the patch normal and velocity to interact with,
+            //  accounting for patch motion if required.
+            void patchData
+            (
+                const parcelType& p,
+                const polyPatch& pp,
+                const scalar trackFraction,
+                const tetIndices& tetIs,
+                vector& normal,
+                vector& Up
+            ) const;
+
 
         // Mapping
 
+            //- Update mesh
+            void updateMesh();
+
             //- Remap the cells of particles corresponding to the
             //  mesh topology change with a default tracking data object
             virtual void autoMap(const mapPolyMesh&);

src/lagrangian/intermediate/clouds/Templates/ReactingCloud/ReactingCloud.C

@@ -87,8 +87,6 @@ void Foam::ReactingCloud<CloudType>::cloudReset(ReactingCloud<CloudType>& c)
 
     compositionModel_.reset(c.compositionModel_.ptr());
     phaseChangeModel_.reset(c.phaseChangeModel_.ptr());
-
-    dMassPhaseChange_ = c.dMassPhaseChange_;
 }
 
 
@@ -111,8 +109,7 @@ Foam::ReactingCloud<CloudType>::ReactingCloud
     constProps_(this->particleProperties(), this->solution().active()),
     compositionModel_(NULL),
     phaseChangeModel_(NULL),
-    rhoTrans_(thermo.carrier().species().size()),
+    rhoTrans_(thermo.carrier().species().size())
-    dMassPhaseChange_(0.0)
 {
     if (this->solution().active())
     {
@@ -167,8 +164,7 @@ Foam::ReactingCloud<CloudType>::ReactingCloud
     constProps_(c.constProps_),
     compositionModel_(c.compositionModel_->clone()),
     phaseChangeModel_(c.phaseChangeModel_->clone()),
-    rhoTrans_(c.rhoTrans_.size()),
+    rhoTrans_(c.rhoTrans_.size())
-    dMassPhaseChange_(c.dMassPhaseChange_)
 {
     forAll(c.rhoTrans_, i)
     {
@@ -209,8 +205,7 @@ Foam::ReactingCloud<CloudType>::ReactingCloud
     compositionModel_(c.compositionModel_->clone()),
 //    compositionModel_(NULL),
     phaseChangeModel_(NULL),
-    rhoTrans_(0),
+    rhoTrans_(0)
-    dMassPhaseChange_(0.0)
 {}
 
 
@@ -350,6 +345,8 @@ void Foam::ReactingCloud<CloudType>::autoMap(const mapPolyMesh& mapper)
     tdType td(*this);
 
     Cloud<parcelType>::template autoMap<tdType>(td, mapper);
+
+    this->updateMesh();
 }
 
 

src/lagrangian/intermediate/clouds/Templates/ReactingCloud/ReactingCloud.H

@@ -121,12 +121,6 @@ protected:
             PtrList<DimensionedField<scalar, volMesh> > rhoTrans_;
 
 
-        // Check
-
-            //- Total mass transferred to continuous phase via phase change
-            scalar dMassPhaseChange_;
-
-
     // Protected Member Functions
 
         // New parcel helper functions

src/lagrangian/intermediate/clouds/Templates/ReactingMultiphaseCloud/ReactingMultiphaseCloud.C

@@ -264,6 +264,8 @@ void Foam::ReactingMultiphaseCloud<CloudType>::autoMap
     tdType td(*this);
 
     Cloud<parcelType>::template autoMap<tdType>(td, mapper);
+
+    this->updateMesh();
 }
 
 

src/lagrangian/intermediate/clouds/Templates/ThermoCloud/ThermoCloud.C

@@ -51,7 +51,69 @@ void Foam::ThermoCloud<CloudType>::setModels()
         ).ptr()
     );
 
+    if (this->solution().coupled())
+    {
         this->subModelProperties().lookup("radiation") >> radiation_;
+    }
+
+    if (radiation_)
+    {
+        radAreaP_.reset
+        (
+            new DimensionedField<scalar, volMesh>
+            (
+                IOobject
+                (
+                    this->name() + "::radAreaP",
+                    this->db().time().timeName(),
+                    this->db(),
+                    IOobject::READ_IF_PRESENT,
+                    IOobject::AUTO_WRITE
+                ),
+                this->mesh(),
+                dimensionedScalar("zero", dimArea, 0.0)
+            )
+        );
+
+        radT4_.reset
+        (
+            new DimensionedField<scalar, volMesh>
+            (
+                IOobject
+                (
+                    this->name() + "::radT4",
+                    this->db().time().timeName(),
+                    this->db(),
+                    IOobject::READ_IF_PRESENT,
+                    IOobject::AUTO_WRITE
+                ),
+                this->mesh(),
+                dimensionedScalar("zero", pow4(dimTemperature), 0.0)
+            )
+        );
+
+        radAreaPT4_.reset
+        (
+            new DimensionedField<scalar, volMesh>
+            (
+                IOobject
+                (
+                    this->name() + "::radAreaPT4",
+                    this->db().time().timeName(),
+                    this->db(),
+                    IOobject::READ_IF_PRESENT,
+                    IOobject::AUTO_WRITE
+                ),
+                this->mesh(),
+                dimensionedScalar
+                (
+                    "zero",
+                    sqr(dimLength)*pow4(dimTemperature),
+                    0.0
+                )
+            )
+        );
+    }
 }
 
 
@@ -98,13 +160,16 @@ Foam::ThermoCloud<CloudType>::ThermoCloud
     heatTransferModel_(NULL),
     TIntegrator_(NULL),
     radiation_(false),
+    radAreaP_(NULL),
+    radT4_(NULL),
+    radAreaPT4_(NULL),
     hsTrans_
     (
         new DimensionedField<scalar, volMesh>
         (
             IOobject
             (
-                this->name() + "hsTrans",
+                this->name() + "::hsTrans",
                 this->db().time().timeName(),
                 this->db(),
                 IOobject::READ_IF_PRESENT,
@@ -120,7 +185,7 @@ Foam::ThermoCloud<CloudType>::ThermoCloud
         (
             IOobject
             (
-                this->name() + "hsCoeff",
+                this->name() + "::hsCoeff",
                 this->db().time().timeName(),
                 this->db(),
                 IOobject::READ_IF_PRESENT,
@@ -165,13 +230,16 @@ Foam::ThermoCloud<CloudType>::ThermoCloud
     heatTransferModel_(c.heatTransferModel_->clone()),
     TIntegrator_(c.TIntegrator_->clone()),
     radiation_(c.radiation_),
+    radAreaP_(NULL),
+    radT4_(NULL),
+    radAreaPT4_(NULL),
     hsTrans_
     (
         new DimensionedField<scalar, volMesh>
         (
             IOobject
             (
-                this->name() + "hsTrans",
+                this->name() + "::hsTrans",
                 this->db().time().timeName(),
                 this->db(),
                 IOobject::NO_READ,
@@ -187,7 +255,7 @@ Foam::ThermoCloud<CloudType>::ThermoCloud
         (
             IOobject
             (
-                this->name() + "hsCoeff",
+                this->name() + "::hsCoeff",
                 this->db().time().timeName(),
                 this->db(),
                 IOobject::NO_READ,
@@ -197,7 +265,61 @@ Foam::ThermoCloud<CloudType>::ThermoCloud
             c.hsCoeff()
         )
     )
-{}
+{
+    if (radiation_)
+    {
+        radAreaP_.reset
+        (
+            new DimensionedField<scalar, volMesh>
+            (
+                IOobject
+                (
+                    this->name() + "::radAreaP",
+                    this->db().time().timeName(),
+                    this->db(),
+                    IOobject::NO_READ,
+                    IOobject::NO_WRITE,
+                    false
+                ),
+                c.radAreaP()
+            )
+        );
+
+        radT4_.reset
+        (
+            new DimensionedField<scalar, volMesh>
+            (
+                IOobject
+                (
+                    this->name() + "::radT4",
+                    this->db().time().timeName(),
+                    this->db(),
+                    IOobject::NO_READ,
+                    IOobject::NO_WRITE,
+                    false
+                ),
+                c.radT4()
+            )
+        );
+
+        radAreaPT4_.reset
+        (
+            new DimensionedField<scalar, volMesh>
+            (
+                IOobject
+                (
+                    this->name() + "::radAreaPT4",
+                    this->db().time().timeName(),
+                    this->db(),
+                    IOobject::NO_READ,
+                    IOobject::NO_WRITE,
+                    false
+                ),
+                c.radAreaPT4()
+            )
+        );
+    }
+}
 
 
 template<class CloudType>
@@ -218,6 +340,9 @@ Foam::ThermoCloud<CloudType>::ThermoCloud
     heatTransferModel_(NULL),
     TIntegrator_(NULL),
     radiation_(false),
+    radAreaP_(NULL),
+    radT4_(NULL),
+    radAreaPT4_(NULL),
     hsTrans_(NULL),
     hsCoeff_(NULL)
 {}
@@ -285,6 +410,13 @@ void Foam::ThermoCloud<CloudType>::resetSourceTerms()
     CloudType::resetSourceTerms();
     hsTrans_->field() = 0.0;
     hsCoeff_->field() = 0.0;
+
+    if (radiation_)
+    {
+        radAreaP_->field() = 0.0;
+        radT4_->field() = 0.0;
+        radAreaPT4_->field() = 0.0;
+    }
 }
 
 
@@ -298,6 +430,13 @@ void Foam::ThermoCloud<CloudType>::relaxSources
 
     this->relax(hsTrans_(), cloudOldTime.hsTrans(), "h");
     this->relax(hsCoeff_(), cloudOldTime.hsCoeff(), "h");
+
+    if (radiation_)
+    {
+        this->relax(radAreaP_(), cloudOldTime.radAreaP(), "radiation");
+        this->relax(radT4_(), cloudOldTime.radT4(), "radiation");
+        this->relax(radAreaPT4_(), cloudOldTime.radAreaPT4(), "radiation");
+    }
 }
 
 
@@ -308,6 +447,13 @@ void Foam::ThermoCloud<CloudType>::scaleSources()
 
     this->scale(hsTrans_(), "h");
     this->scale(hsCoeff_(), "h");
+
+    if (radiation_)
+    {
+        this->scale(radAreaP_(), "radiation");
+        this->scale(radT4_(), "radiation");
+        this->scale(radAreaPT4_(), "radiation");
+    }
 }
 
 
@@ -341,6 +487,8 @@ void Foam::ThermoCloud<CloudType>::autoMap(const mapPolyMesh& mapper)
     tdType td(*this);
 
     Cloud<parcelType>::template autoMap<tdType>(td, mapper);
+
+    this->updateMesh();
 }
 
 

src/lagrangian/intermediate/clouds/Templates/ThermoCloud/ThermoCloud.H

@@ -132,6 +132,15 @@ protected:
             //- Include radiation
             Switch radiation_;
 
+            //- Radiation sum of parcel projected areas
+            autoPtr<DimensionedField<scalar, volMesh> > radAreaP_;
+
+            //- Radiation sum of parcel temperature^4
+            autoPtr<DimensionedField<scalar, volMesh> > radT4_;
+
+            //- Radiation sum of parcel projected areas * temperature^4
+            autoPtr<DimensionedField<scalar, volMesh> > radAreaPT4_;
+
 
         // Sources
 
@@ -244,6 +253,26 @@ public:
                  //- Radiation flag
                 inline bool radiation() const;
 
+                //- Radiation sum of parcel projected areas [m2]
+                inline DimensionedField<scalar, volMesh>& radAreaP();
+
+                //- Radiation sum of parcel projected areas [m2]
+                inline const DimensionedField<scalar, volMesh>&
+                    radAreaP() const;
+
+                //- Radiation sum of parcel temperature^4 [K4]
+                inline DimensionedField<scalar, volMesh>& radT4();
+
+                //- Radiation sum of parcel temperature^4 [K4]
+                inline const DimensionedField<scalar, volMesh>& radT4() const;
+
+                //- Radiation sum of parcel projected area*temperature^4 [m2K4]
+                inline DimensionedField<scalar, volMesh>& radAreaPT4();
+
+                //- Radiation sum of parcel temperature^4 [m2K4]
+                inline const DimensionedField<scalar, volMesh>&
+                    radAreaPT4() const;
+
 
             // Sources
 

src/lagrangian/intermediate/clouds/Templates/ThermoCloud/ThermoCloudI.H

@@ -89,6 +89,114 @@ inline bool Foam::ThermoCloud<CloudType>::radiation() const
 }
 
 
+template<class CloudType>
+inline Foam::DimensionedField<Foam::scalar, Foam::volMesh>&
+Foam::ThermoCloud<CloudType>::radAreaP()
+{
+    if (!radiation_)
+    {
+        FatalErrorIn
+        (
+            "inline Foam::DimensionedField<Foam::scalar, Foam::volMesh> "
+            "Foam::ThermoCloud<CloudType>::radAreaP()"
+        )   << "Radiation field requested, but radiation model not active"
+            << abort(FatalError);
+    }
+
+    return radAreaP_();
+}
+
+
+template<class CloudType>
+inline const Foam::DimensionedField<Foam::scalar, Foam::volMesh>&
+Foam::ThermoCloud<CloudType>::radAreaP() const
+{
+    if (!radiation_)
+    {
+        FatalErrorIn
+        (
+            "inline Foam::DimensionedField<Foam::scalar, Foam::volMesh> "
+            "Foam::ThermoCloud<CloudType>::radAreaP()"
+        )   << "Radiation field requested, but radiation model not active"
+            << abort(FatalError);
+    }
+
+    return radAreaP_();
+}
+
+
+template<class CloudType>
+inline Foam::DimensionedField<Foam::scalar, Foam::volMesh>&
+Foam::ThermoCloud<CloudType>::radT4()
+{
+    if (!radiation_)
+    {
+        FatalErrorIn
+        (
+            "inline Foam::DimensionedField<Foam::scalar, Foam::volMesh> "
+            "Foam::ThermoCloud<CloudType>::radT4()"
+        )   << "Radiation field requested, but radiation model not active"
+            << abort(FatalError);
+    }
+
+    return radT4_();
+}
+
+
+template<class CloudType>
+inline const Foam::DimensionedField<Foam::scalar, Foam::volMesh>&
+Foam::ThermoCloud<CloudType>::radT4() const
+{
+    if (!radiation_)
+    {
+        FatalErrorIn
+        (
+            "inline Foam::DimensionedField<Foam::scalar, Foam::volMesh> "
+            "Foam::ThermoCloud<CloudType>::radT4()"
+        )   << "Radiation field requested, but radiation model not active"
+            << abort(FatalError);
+    }
+
+    return radT4_();
+}
+
+
+template<class CloudType>
+inline Foam::DimensionedField<Foam::scalar, Foam::volMesh>&
+Foam::ThermoCloud<CloudType>::radAreaPT4()
+{
+    if (!radiation_)
+    {
+        FatalErrorIn
+        (
+            "inline Foam::DimensionedField<Foam::scalar, Foam::volMesh> "
+            "Foam::ThermoCloud<CloudType>::radAreaPT4()"
+        )   << "Radiation field requested, but radiation model not active"
+            << abort(FatalError);
+    }
+
+    return radAreaPT4_();
+}
+
+
+template<class CloudType>
+inline const Foam::DimensionedField<Foam::scalar, Foam::volMesh>&
+Foam::ThermoCloud<CloudType>::radAreaPT4() const
+{
+    if (!radiation_)
+    {
+        FatalErrorIn
+        (
+            "inline Foam::DimensionedField<Foam::scalar, Foam::volMesh> "
+            "Foam::ThermoCloud<CloudType>::radAreaPT4()"
+        )   << "Radiation field requested, but radiation model not active"
+            << abort(FatalError);
+    }
+
+    return radAreaPT4_();
+}
+
+
 template<class CloudType>
 inline Foam::DimensionedField<Foam::scalar, Foam::volMesh>&
 Foam::ThermoCloud<CloudType>::hsTrans()
@@ -182,21 +290,15 @@ inline Foam::tmp<Foam::volScalarField> Foam::ThermoCloud<CloudType>::Ep() const
         )
     );
 
-    // Need to check if coupled as field is created on-the-fly
+    if (radiation_)
-    if (radiation_ && this->solution().coupled())
     {
         scalarField& Ep = tEp().internalField();
+        const scalar dt = this->db().time().deltaTValue();
         const scalarField& V = this->mesh().V();
         const scalar epsilon = constProps_.epsilon0();
+        const scalarField& sumAreaPT4 = radAreaPT4_->field();
 
-        forAllConstIter(typename ThermoCloud<CloudType>, *this, iter)
+        Ep = sumAreaPT4*epsilon*physicoChemical::sigma.value()/V/dt;
-        {
-            const parcelType& p = iter();
-            const label cellI = p.cell();
-            Ep[cellI] += p.nParticle()*p.areaP()*pow4(p.T());
-        }
-
-        Ep *= epsilon*physicoChemical::sigma.value()/V;
     }
 
     return tEp;
@@ -224,21 +326,15 @@ inline Foam::tmp<Foam::volScalarField> Foam::ThermoCloud<CloudType>::ap() const
         )
     );
 
-    // Need to check if coupled as field is created on-the-fly
+    if (radiation_)
-    if (radiation_ && this->solution().coupled())
     {
         scalarField& ap = tap().internalField();
+        const scalar dt = this->db().time().deltaTValue();
         const scalarField& V = this->mesh().V();
         const scalar epsilon = constProps_.epsilon0();
+        const scalarField& sumAreaP = radAreaP_->field();
 
-        forAllConstIter(typename ThermoCloud<CloudType>, *this, iter)
+        ap = sumAreaP*epsilon/V/dt;
-        {
-            const parcelType& p = iter();
-            const label cellI = p.cell();
-            ap[cellI] += p.nParticle()*p.areaP();
-        }
-
-        ap *= epsilon/V;
     }
 
     return tap;
@@ -267,23 +363,16 @@ Foam::ThermoCloud<CloudType>::sigmap() const
         )
     );
 
-    // Need to check if coupled as field is created on-the-fly
+    if (radiation_)
-    if (radiation_ && this->solution().coupled())
     {
         scalarField& sigmap = tsigmap().internalField();
-
+        const scalar dt = this->db().time().deltaTValue();
         const scalarField& V = this->mesh().V();
         const scalar epsilon = constProps_.epsilon0();
         const scalar f = constProps_.f0();
+        const scalarField& sumAreaP = radAreaP_->field();
 
-        forAllConstIter(typename ThermoCloud<CloudType>, *this, iter)
+        sigmap *= sumAreaP*(1.0 - f)*(1.0 - epsilon)/V/dt;
-        {
-            const parcelType& p = iter();
-            const label cellI = p.cell();
-            sigmap[cellI] += p.nParticle()*p.areaP();
-        }
-
-        sigmap *= (1.0 - f)*(1.0 - epsilon)/V;
     }
 
     return tsigmap;

src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcel.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -363,7 +363,7 @@ bool Foam::KinematicParcel<ParcelType>::hitPatch
         static_cast<typename TrackData::cloudType::parcelType&>(*this);
 
     // Invoke post-processing model
-    td.cloud().functions().postPatch(p, patchI, pp.whichFace(p.face()));
+    td.cloud().functions().postPatch(p, pp, trackFraction, tetIs);
 
     // Invoke surface film model
     if (td.cloud().surfaceFilm().transferParcel(p, pp, td.keepParticle))

src/lagrangian/intermediate/parcels/Templates/ReactingMultiphaseParcel/ReactingMultiphaseParcel.C

@@ -268,12 +268,15 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::calc
     (
         td,
         dt,
+        this->age_,
         Ts,
         d0,
         T0,
         mass0,
         this->mass0_,
         YMix[GAS]*YGas_,
+        YMix[LIQ]*YLiquid_,
+        YMix[SLD]*YSolid_,
         canCombust_,
         dMassDV,
         Sh,
@@ -466,6 +469,16 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::calc
         // Update sensible enthalpy transfer
         td.cloud().hsTrans()[cellI] += np0*dhsTrans;
         td.cloud().hsCoeff()[cellI] += np0*Sph;
+
+        // Update radiation fields
+        if (td.cloud().radiation())
+        {
+            const scalar ap = this->areaP();
+            const scalar T4 = pow4(this->T_);
+            td.cloud().radAreaP()[cellI] += dt*np0*ap;
+            td.cloud().radT4()[cellI] += dt*np0*T4;
+            td.cloud().radAreaP()[cellI] += dt*np0*ap*T4;
+        }
     }
 }
 
@@ -476,12 +489,15 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::calcDevolatilisation
 (
     TrackData& td,
     const scalar dt,
+    const scalar age,
     const scalar Ts,
     const scalar d,
     const scalar T,
     const scalar mass,
     const scalar mass0,
     const scalarField& YGasEff,
+    const scalarField& YLiquidEff,
+    const scalarField& YSolidEff,
     bool& canCombust,
     scalarField& dMassDV,
     scalar& Sh,
@@ -510,10 +526,13 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::calcDevolatilisation
     td.cloud().devolatilisation().calculate
     (
         dt,
+        age,
         mass0,
         mass,
         T,
         YGasEff,
+        YLiquidEff,
+        YSolidEff,
         canCombust,
         dMassDV
     );
@@ -638,7 +657,8 @@ Foam::ReactingMultiphaseParcel<ParcelType>::ReactingMultiphaseParcel
     ParcelType(p),
     YGas_(p.YGas_),
     YLiquid_(p.YLiquid_),
-    YSolid_(p.YSolid_)
+    YSolid_(p.YSolid_),
+    canCombust_(p.canCombust_)
 {}
 
 
@@ -652,7 +672,8 @@ Foam::ReactingMultiphaseParcel<ParcelType>::ReactingMultiphaseParcel
     ParcelType(p, mesh),
     YGas_(p.YGas_),
     YLiquid_(p.YLiquid_),
-    YSolid_(p.YSolid_)
+    YSolid_(p.YSolid_),
+    canCombust_(p.canCombust_)
 {}
 
 

src/lagrangian/intermediate/parcels/Templates/ReactingMultiphaseParcel/ReactingMultiphaseParcel.H

@@ -196,12 +196,15 @@ protected:
         (
             TrackData& td,
             const scalar dt,           // timestep
-            const scalar Ts,           // Surface temperature
+            const scalar age,          // age
+            const scalar Ts,           // surface temperature
             const scalar d,            // diameter
             const scalar T,            // temperature
             const scalar mass,         // mass
             const scalar mass0,        // mass (initial on injection)
-            const scalarField& YGasEff,// Gas component mass fractions
+            const scalarField& YGasEff,// gas component mass fractions
+            const scalarField& YLiquidEff,// liquid component mass fractions
+            const scalarField& YSolidEff,// solid component mass fractions
             bool& canCombust,          // 'can combust' flag
             scalarField& dMassDV,      // mass transfer - local to particle
             scalar& Sh,                // explicit particle enthalpy source

src/lagrangian/intermediate/parcels/Templates/ReactingParcel/ReactingParcel.C

@@ -458,6 +458,16 @@ void Foam::ReactingParcel<ParcelType>::calc
         // Update sensible enthalpy transfer
         td.cloud().hsTrans()[cellI] += np0*dhsTrans;
         td.cloud().hsCoeff()[cellI] += np0*Sph;
+
+        // Update radiation fields
+        if (td.cloud().radiation())
+        {
+            const scalar ap = this->areaP();
+            const scalar T4 = pow4(this->T_);
+            td.cloud().radAreaP()[cellI] += dt*np0*ap;
+            td.cloud().radT4()[cellI] += dt*np0*T4;
+            td.cloud().radAreaP()[cellI] += dt*np0*ap*T4;
+        }
     }
 }
 

src/lagrangian/intermediate/parcels/Templates/ThermoParcel/ThermoParcel.C

@@ -144,7 +144,7 @@ void Foam::ThermoParcel<ParcelType>::calcSurfaceValues
         Ts = td.cloud().constProps().TMin();
     }
 
-    // Assuming thermo props vary linearly with T for small dT
+    // Assuming thermo props vary linearly with T for small d(T)
     const scalar TRatio = Tc_/Ts;
 
     rhos = this->rhoc_*TRatio;
@@ -252,6 +252,16 @@ void Foam::ThermoParcel<ParcelType>::calc
 
         // Update sensible enthalpy coefficient
         td.cloud().hsCoeff()[cellI] += np0*Sph;
+
+        // Update radiation fields
+        if (td.cloud().radiation())
+        {
+            const scalar ap = this->areaP();
+            const scalar T4 = pow4(this->T_);
+            td.cloud().radAreaP()[cellI] += dt*np0*ap;
+            td.cloud().radT4()[cellI] += dt*np0*T4;
+            td.cloud().radAreaP()[cellI] += dt*np0*ap*T4;
+        }
     }
 }
 

src/lagrangian/intermediate/parcels/include/makeParcelInjectionModels.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -37,7 +37,7 @@ License
 #include "ManualInjection.H"
 #include "NoInjection.H"
 #include "PatchInjection.H"
-
+#include "PatchFlowRateInjection.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -53,7 +53,8 @@ License
     makeInjectionModelType(KinematicLookupTableInjection, CloudType);         \
     makeInjectionModelType(ManualInjection, CloudType);                       \
     makeInjectionModelType(NoInjection, CloudType);                           \
-    makeInjectionModelType(PatchInjection, CloudType);
+    makeInjectionModelType(PatchInjection, CloudType);                        \
+    makeInjectionModelType(PatchFlowRateInjection, CloudType);
 
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

src/lagrangian/intermediate/parcels/include/makeReactingMultiphaseParcelInjectionModels.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -35,6 +35,7 @@ License
 #include "ManualInjection.H"
 #include "NoInjection.H"
 #include "PatchInjection.H"
+#include "PatchFlowRateInjection.H"
 #include "ReactingMultiphaseLookupTableInjection.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -49,6 +50,7 @@ License
     makeInjectionModelType(ManualInjection, CloudType);                       \
     makeInjectionModelType(NoInjection, CloudType);                           \
     makeInjectionModelType(PatchInjection, CloudType);                        \
+    makeInjectionModelType(PatchFlowRateInjection, CloudType);                \
     makeInjectionModelType(ReactingMultiphaseLookupTableInjection, CloudType);
 
 

src/lagrangian/intermediate/parcels/include/makeReactingParcelInjectionModels.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -35,6 +35,7 @@ License
 #include "ManualInjection.H"
 #include "NoInjection.H"
 #include "PatchInjection.H"
+#include "PatchFlowRateInjection.H"
 #include "ReactingLookupTableInjection.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -49,6 +50,7 @@ License
     makeInjectionModelType(ManualInjection, CloudType);                       \
     makeInjectionModelType(NoInjection, CloudType);                           \
     makeInjectionModelType(PatchInjection, CloudType);                        \
+    makeInjectionModelType(PatchFlowRateInjection, CloudType);                \
     makeInjectionModelType(ReactingLookupTableInjection, CloudType);
 
 

src/lagrangian/intermediate/submodels/CloudFunctionObjects/CloudFunctionObject/CloudFunctionObject.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -107,8 +107,9 @@ template<class CloudType>
 void Foam::CloudFunctionObject<CloudType>::postPatch
 (
     const typename CloudType::parcelType&,
-    const label,
+    const polyPatch&,
-    const label
+    const scalar,
+    const tetIndices&
 )
 {
     // do nothing

src/lagrangian/intermediate/submodels/CloudFunctionObjects/CloudFunctionObject/CloudFunctionObject.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -141,8 +141,9 @@ public:
             virtual void postPatch
             (
                 const typename CloudType::parcelType& p,
-                const label patchI,
+                const polyPatch& pp,
-                const label patchFaceI
+                const scalar trackFraction,
+                const tetIndices& testIs
             );
 
             //- Post-face hook

src/lagrangian/intermediate/submodels/CloudFunctionObjects/CloudFunctionObjectList/CloudFunctionObjectList.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -146,13 +146,14 @@ template<class CloudType>
 void Foam::CloudFunctionObjectList<CloudType>::postPatch
 (
     const typename CloudType::parcelType& p,
-    const label patchI,
+    const polyPatch& pp,
-    const label patchFaceI
+    const scalar trackFraction,
+    const tetIndices& tetIs
 )
 {
     forAll(*this, i)
     {
-        this->operator[](i).postPatch(p, patchI, patchFaceI);
+        this->operator[](i).postPatch(p, pp, trackFraction, tetIs);
     }
 }
 

src/lagrangian/intermediate/submodels/CloudFunctionObjects/CloudFunctionObjectList/CloudFunctionObjectList.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -121,8 +121,9 @@ public:
             virtual void postPatch
             (
                 const typename CloudType::parcelType& p,
-                const label patchI,
+                const polyPatch& pp,
-                const label patchFaceI
+                const scalar trackFraction,
+                const tetIndices& tetIs
             );
 
             //- Post-face hook

src/lagrangian/intermediate/submodels/CloudFunctionObjects/ParticleErosion/ParticleErosion.C

@@ -166,29 +166,32 @@ template<class CloudType>
 void Foam::ParticleErosion<CloudType>::postPatch
 (
     const parcelType& p,
-    const label patchI,
+    const polyPatch& pp,
-    const label patchFaceI
+    const scalar trackFraction,
+    const tetIndices& tetIs
 )
 {
+    const label patchI = pp.index();
+
     const label localPatchI = applyToPatch(patchI);
 
     if (localPatchI != -1)
     {
-        const fvMesh& mesh = this->owner().mesh();
+        vector nw;
+        vector Up;
 
         // patch-normal direction
-        vector nw = p.currentTetIndices().faceTri(mesh).normal();
+        this->owner().patchData(p, pp, trackFraction, tetIs, nw, Up);
 
-        // particle direction of travel
+        // particle velocity reletive to patch
-        const vector& U = p.U();
+        const vector& U = p.U() - Up;
 
         // quick reject if particle travelling away from the patch
-        if ((-nw & U) < 0)
+        if ((nw & U) < 0)
         {
             return;
         }
 
-        nw /= mag(nw);
         const scalar magU = mag(U);
         const vector Udir = U/magU;
 
@@ -197,6 +200,7 @@ void Foam::ParticleErosion<CloudType>::postPatch
 
         const scalar coeff = p.nParticle()*p.mass()*sqr(magU)/(p_*psi_*K_);
 
+        const label patchFaceI = pp.whichFace(p.face());
         scalar& Q = QPtr_->boundaryField()[patchI][patchFaceI];
         if (tan(alpha) < K_/6.0)
         {

src/lagrangian/intermediate/submodels/CloudFunctionObjects/ParticleErosion/ParticleErosion.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -126,8 +126,9 @@ public:
             virtual void postPatch
             (
                 const parcelType& p,
-                const label patchI,
+                const polyPatch& pp,
-                const label patchFaceI
+                const scalar trackFraction,
+                const tetIndices& tetIs
             );
 };
 

src/lagrangian/intermediate/submodels/CloudFunctionObjects/PatchPostProcessing/PatchPostProcessing.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -218,11 +218,14 @@ template<class CloudType>
 void Foam::PatchPostProcessing<CloudType>::postPatch
 (
     const parcelType& p,
-    const label patchI,
+    const polyPatch& pp,
-    const label
+    const scalar,
+    const tetIndices& tetIs
 )
 {
+    const label patchI = pp.index();
     const label localPatchI = applyToPatch(patchI);
+
     if (localPatchI != -1 && patchData_[localPatchI].size() < maxStoredParcels_)
     {
         times_[localPatchI].append(this->owner().time().value());

src/lagrangian/intermediate/submodels/CloudFunctionObjects/PatchPostProcessing/PatchPostProcessing.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -127,8 +127,9 @@ public:
             virtual void postPatch
             (
                 const parcelType& p,
-                const label patchI,
+                const polyPatch& pp,
-                const label patchFaceI
+                const scalar trackFraction,
+                const tetIndices& tetIs
             );
 };
 

src/lagrangian/intermediate/submodels/Kinematic/CollisionModel/PairCollision/PairCollision.C

@@ -306,8 +306,9 @@ void Foam::PairCollision<CloudType>::wallInteraction()
                         this->owner().functions().postPatch
                         (
                             p,
-                            patchI,
+                            mesh.boundaryMesh()[patchI],
-                            patchFaceI
+                            1.0,
+                            p.currentTetIndices()
                         );
                      }
                 }