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

applications/solvers/DNS/dnsFoam/Make/options

@@ -1,10 +1,7 @@
 EXE_INC = \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
-    -I$(LIB_SRC)/randomProcesses/lnInclude \
+    -I$(LIB_SRC)/randomProcesses/lnInclude
-    -I$(LIB_SRC)/sampling/lnInclude
 
 EXE_LIBS = \
     -lfiniteVolume \
-    -lrandomProcesses \
+    -lrandomProcesses
-    -lsampling \
-    -lmeshTools

applications/solvers/basic/scalarTransportFoam/Make/options

@@ -1,4 +1,11 @@
 EXE_INC = \
-    -I$(LIB_SRC)/finiteVolume/lnInclude
+    -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/fieldSources/lnInclude \
+    -I$(LIB_SRC)/meshTools/lnInclude \
+    -I$(LIB_SRC)/sampling/lnInclude
 
-EXE_LIBS = -lfiniteVolume
+EXE_LIBS = \
+    -lfiniteVolume \
+    -lfieldSources \
+    -lmeshTools \
+    -lsampling

applications/solvers/basic/scalarTransportFoam/createFields.H

@@ -45,7 +45,7 @@
     );
 
 
-    Info<< "Reading diffusivity D\n" << endl;
+    Info<< "Reading diffusivity DT\n" << endl;
 
     dimensionedScalar DT
     (

applications/solvers/basic/scalarTransportFoam/scalarTransportFoam.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
@@ -30,6 +30,7 @@ Description
 \*---------------------------------------------------------------------------*/
 
 #include "fvCFD.H"
+#include "IObasicSourceList.H"
 #include "simpleControl.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -40,6 +41,7 @@ int main(int argc, char *argv[])
     #include "createTime.H"
     #include "createMesh.H"
     #include "createFields.H"
+    #include "createSources.H"
 
     simpleControl simple(mesh);
 
@@ -60,6 +62,8 @@ int main(int argc, char *argv[])
                 fvm::ddt(T)
               + fvm::div(phi, T)
               - fvm::laplacian(DT, T)
+             ==
+                sources(T)
             );
         }
 

applications/solvers/combustion/PDRFoam/Make/options

@@ -5,7 +5,6 @@ EXE_INC = \
     -IPDRModels/dragModels/PDRDragModel \
     -IlaminarFlameSpeed/SCOPE \
     -I$(LIB_SRC)/engine/lnInclude \
-    -I$(LIB_SRC)/sampling/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/specie/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/reactionThermo/lnInclude \

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

@@ -116,7 +116,7 @@ void PDRkEpsilon::correct()
     }
 
     tmp<volTensorField> tgradU = fvc::grad(U_);
-    volScalarField G("RASModel.G", mut_*(tgradU() && dev(twoSymm(tgradU()))));
+    volScalarField G(type() + ".G", mut_*(tgradU() && dev(twoSymm(tgradU()))));
     tgradU.clear();
 
     // Update espsilon and G at the wall

applications/solvers/combustion/XiFoam/EaEqn.H

@@ -16,10 +16,16 @@
           : -dpdt
         )
       - fvm::laplacian(turbulence->alphaEff(), hea)
+      + sources(rho, hea)
     );
 
     EaEqn.relax();
+
+    sources.constrain(EaEqn);
+
     EaEqn.solve();
 
+    sources.correct(hea);
+
     thermo.correct();
 }

applications/solvers/combustion/XiFoam/EauEqn.H

@@ -2,7 +2,7 @@ if (ign.ignited())
 {
     volScalarField& heau = thermo.heu();
 
-    solve
+    fvScalarMatrix heauEqn
     (
         fvm::ddt(rho, heau) + mvConvection->fvmDiv(phi, heau)
       + (fvc::ddt(rho, K) + fvc::div(phi, K))*rho/thermo.rhou()
@@ -23,5 +23,14 @@ if (ign.ignited())
         // A possible solution would be to solve for ftu as well as ft.
         //- fvm::div(muEff*fvc::grad(b)/(b + 0.001), heau)
         //+ fvm::Sp(fvc::div(muEff*fvc::grad(b)/(b + 0.001)), heau)
+
+     ==
+        sources(rho, heau)
     );
+
+    sources.constrain(heauEqn);
+
+    heauEqn.solve();
+
+    sources.correct(heau);
 }

applications/solvers/combustion/XiFoam/Make/options

@@ -1,14 +1,20 @@
 EXE_INC = \
+    -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/fieldSources/lnInclude \
+    -I$(LIB_SRC)/sampling/lnInclude \
+    -I$(LIB_SRC)/meshTools/lnInclude\
     -I$(LIB_SRC)/engine/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/specie/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/reactionThermo/lnInclude \
     -I$(LIB_SRC)/turbulenceModels/compressible/turbulenceModel \
-    -I$(LIB_SRC)/thermophysicalModels/laminarFlameSpeed/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/laminarFlameSpeed/lnInclude
-    -I$(LIB_SRC)/sampling/lnInclude \
-    -I$(LIB_SRC)/finiteVolume/lnInclude
 
 EXE_LIBS = \
+    -lfiniteVolume \
+    -lfieldSources \
+    -lsampling \
+    -lmeshTools \
     -lengine \
     -lcompressibleTurbulenceModel \
     -lcompressibleRASModels \
@@ -16,6 +22,4 @@ EXE_LIBS = \
     -lfluidThermophysicalModels \
     -lreactionThermophysicalModels \
     -lspecie \
-    -llaminarFlameSpeedModels \
+    -llaminarFlameSpeedModels
-    -lfiniteVolume \
-    -lmeshTools

applications/solvers/combustion/XiFoam/UEqn.H

@@ -5,12 +5,17 @@
       + turbulence->divDevRhoReff(U)
      ==
         rho*g
+      + sources(rho, U)
     );
 
     UEqn.relax();
 
+    sources.constrain(UEqn);
+
     if (pimple.momentumPredictor())
     {
         solve(UEqn == -fvc::grad(p));
+
+        sources.correct(U);
         K = 0.5*magSqr(U);
     }

applications/solvers/combustion/XiFoam/XiFoam.C

@@ -56,6 +56,7 @@ Description
 #include "ignition.H"
 #include "Switch.H"
 #include "pimpleControl.H"
+#include "IObasicSourceList.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -68,6 +69,7 @@ int main(int argc, char *argv[])
     #include "readCombustionProperties.H"
     #include "readGravitationalAcceleration.H"
     #include "createFields.H"
+    #include "createSources.H"
     #include "initContinuityErrs.H"
     #include "readTimeControls.H"
     #include "compressibleCourantNo.H"

applications/solvers/combustion/XiFoam/bEqn.H

@@ -34,7 +34,7 @@ if (ign.ignited())
 
     // Calculate turbulent flame speed flux
     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-    surfaceScalarField phiSt(fvc::interpolate(rhou*StCorr*Su*Xi)*nf);
+    surfaceScalarField phiSt("phiSt", fvc::interpolate(rhou*StCorr*Su*Xi)*nf);
 
     scalar StCoNum = max
     (
@@ -50,9 +50,11 @@ if (ign.ignited())
     (
         fvm::ddt(rho, b)
       + mvConvection->fvmDiv(phi, b)
-      + fvm::div(phiSt, b, "div(phiSt,b)")
+      + fvm::div(phiSt, b)
       - fvm::Sp(fvc::div(phiSt), b)
       - fvm::laplacian(turbulence->alphaEff(), b)
+     ==
+        sources(rho, b)
     );
 
 
@@ -64,8 +66,13 @@ if (ign.ignited())
     // Solve for b
     // ~~~~~~~~~~~
     bEqn.relax();
+
+    sources.constrain(bEqn);
+
     bEqn.solve();
 
+    sources.correct(b);
+
     Info<< "min(b) = " << min(b).value() << endl;
 
 
@@ -153,11 +160,17 @@ if (ign.ignited())
           ==
           - fvm::SuSp(-rho*Rc*Su0/Su, Su)
           - fvm::SuSp(rho*(sigmas + Rc), Su)
+          + sources(rho, Su)
         );
 
         SuEqn.relax();
+
+        sources.constrain(SuEqn);
+
         SuEqn.solve();
 
+        sources.correct(Su);
+
         // Limit the maximum Su
         // ~~~~~~~~~~~~~~~~~~~~
         Su.min(SuMax);
@@ -234,11 +247,17 @@ if (ign.ignited())
                 ),
                 Xi
             )
+          + sources(rho, Xi)
         );
 
         XiEqn.relax();
+
+        sources.constrain(XiEqn);
+
         XiEqn.solve();
 
+        sources.correct(Xi);
+
         // Correct boundedness of Xi
         // ~~~~~~~~~~~~~~~~~~~~~~~~~
         Xi.max(1.0);

applications/solvers/combustion/XiFoam/ftEqn.H

@@ -13,10 +13,18 @@ if (composition.contains("ft"))
 {
     volScalarField& ft = composition.Y("ft");
 
-    solve
+    fvScalarMatrix ftEqn
     (
         fvm::ddt(rho, ft)
       + mvConvection->fvmDiv(phi, ft)
       - fvm::laplacian(turbulence->alphaEff(), ft)
+     ==
+        sources(rho, ft)
     );
+
+    sources.constrain(ftEqn);
+
+    ftEqn.solve();
+
+    sources.correct(ft);
 }

applications/solvers/combustion/XiFoam/pEqn.H

@@ -16,6 +16,8 @@ if (pimple.transonic())
         )
     );
 
+    sources.relativeFlux(fvc::interpolate(psi), phid);
+
     while (pimple.correctNonOrthogonal())
     {
         fvScalarMatrix pEqn
@@ -23,8 +25,12 @@ if (pimple.transonic())
             fvm::ddt(psi, p)
           + fvm::div(phid, p)
           - fvm::laplacian(rho*rAU, p)
+         ==
+            sources(psi, p, rho.name())
         );
 
+        sources.constrain(pEqn);
+
         pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter())));
 
         if (pimple.finalNonOrthogonalIter())
@@ -45,6 +51,8 @@ else
         )
     );
 
+    sources.relativeFlux(phiHbyA);
+
     while (pimple.correctNonOrthogonal())
     {
         fvScalarMatrix pEqn
@@ -52,8 +60,12 @@ else
             fvm::ddt(psi, p)
           + fvc::div(phiHbyA)
           - fvm::laplacian(rho*rAU, p)
+         ==
+            sources(psi, p, rho.name())
         );
 
+        sources.constrain(pEqn);
+
         pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter())));
 
         if (pimple.finalNonOrthogonalIter())
@@ -68,6 +80,7 @@ else
 
 U = HbyA - rAU*fvc::grad(p);
 U.correctBoundaryConditions();
+sources.correct(U);
 K = 0.5*magSqr(U);
 
 if (thermo.dpdt())

applications/solvers/combustion/chemFoam/Make/options

@@ -12,6 +12,7 @@ EXE_INC = \
 
 EXE_LIBS = \
     -lfiniteVolume \
+    -lcompressibleTurbulenceModel \
     -lcompressibleRASModels \
     -lreactionThermophysicalModels \
     -lfluidThermophysicalModels \

applications/solvers/combustion/coldEngineFoam/Make/options

@@ -20,6 +20,5 @@ EXE_LIBS = \
     -lfluidThermophysicalModels \
     -lspecie \
     -lfiniteVolume \
-    -lsampling \
     -lmeshTools \
     -lfieldSources

applications/solvers/combustion/coldEngineFoam/coldEngineFoam.C

@@ -47,6 +47,7 @@ int main(int argc, char *argv[])
     #include "createEngineTime.H"
     #include "createEngineMesh.H"
     #include "createFields.H"
+    #include "createSources.H"
     #include "initContinuityErrs.H"
     #include "readEngineTimeControls.H"
     #include "compressibleCourantNo.H"

applications/solvers/combustion/coldEngineFoam/createFields.H

@@ -69,6 +69,3 @@
 
     Info<< "Creating field kinetic energy K\n" << endl;
     volScalarField K("K", 0.5*magSqr(U));
-
-    Info<< "Creating sources\n" << endl;
-    IObasicSourceList sources(mesh);

applications/solvers/combustion/engineFoam/Make/options

@@ -1,14 +1,21 @@
 EXE_INC = \
-    -I../XiFoam \
+    -I$(FOAM_SOLVERS)/combustion/XiFoam \
+    -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/fieldSources/lnInclude \
+    -I$(LIB_SRC)/sampling/lnInclude \
+    -I$(LIB_SRC)/meshTools/lnInclude \
     -I$(LIB_SRC)/engine/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/specie/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/reactionThermo/lnInclude \
     -I$(LIB_SRC)/turbulenceModels/compressible/turbulenceModel \
-    -I$(LIB_SRC)/thermophysicalModels/laminarFlameSpeed/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/laminarFlameSpeed/lnInclude
-    -I$(LIB_SRC)/finiteVolume/lnInclude
 
 EXE_LIBS = \
+    -lfiniteVolume \
+    -lfieldSources \
+    -lsampling \
+    -lmeshTools \
     -lengine \
     -lcompressibleTurbulenceModel \
     -lcompressibleRASModels \
@@ -16,5 +23,4 @@ EXE_LIBS = \
     -lfluidThermophysicalModels \
     -lreactionThermophysicalModels \
     -lspecie \
-    -llaminarFlameSpeedModels \
+    -llaminarFlameSpeedModels
-    -lfiniteVolume

applications/solvers/combustion/engineFoam/UEqn.H

@@ -3,10 +3,14 @@
         fvm::ddt(rho, U)
       + fvm::div(phi, U)
       + turbulence->divDevRhoReff(U)
+     ==
+        sources(rho, U)
     );
 
     if (pimple.momentumPredictor())
     {
         solve(UEqn == -fvc::grad(p));
+
+        sources.correct(U);
         K = 0.5*magSqr(U);
     }

applications/solvers/combustion/engineFoam/engineFoam.C

@@ -59,6 +59,7 @@ Description
 #include "OFstream.H"
 #include "mathematicalConstants.H"
 #include "pimpleControl.H"
+#include "IObasicSourceList.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -70,6 +71,7 @@ int main(int argc, char *argv[])
     #include "createEngineMesh.H"
     #include "readCombustionProperties.H"
     #include "createFields.H"
+    #include "createSources.H"
     #include "initContinuityErrs.H"
     #include "readEngineTimeControls.H"
     #include "compressibleCourantNo.H"

applications/solvers/combustion/engineFoam/pEqn.H

@@ -13,6 +13,8 @@ if (pimple.transonic())
        *((fvc::interpolate(HbyA) & mesh.Sf()) - fvc::meshPhi(rho, U))
     );
 
+    sources.relativeFlux(fvc::interpolate(psi), phid);
+
     while (pimple.correctNonOrthogonal())
     {
         fvScalarMatrix pEqn
@@ -20,8 +22,12 @@ if (pimple.transonic())
             fvm::ddt(psi, p)
           + fvm::div(phid, p)
           - fvm::laplacian(rho*rAU, p)
+         ==
+            sources(psi, p, rho.name())
         );
 
+        sources.constrain(pEqn);
+
         pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter())));
 
         if (pimple.finalNonOrthogonalIter())
@@ -39,6 +45,8 @@ else
        *((fvc::interpolate(HbyA) & mesh.Sf()) - fvc::meshPhi(rho, U))
     );
 
+    sources.relativeFlux(fvc::interpolate(rho), phiHbyA);
+
     while (pimple.correctNonOrthogonal())
     {
         fvScalarMatrix pEqn
@@ -46,8 +54,12 @@ else
             fvm::ddt(psi, p)
           + fvc::div(phiHbyA)
           - fvm::laplacian(rho*rAU, p)
+         ==
+            sources(psi, p, rho.name())
         );
 
+        sources.constrain(pEqn);
+
         pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter())));
 
         if (pimple.finalNonOrthogonalIter())
@@ -62,6 +74,7 @@ else
 
 U = HbyA - rAU*fvc::grad(p);
 U.correctBoundaryConditions();
+sources.correct(U);
 K = 0.5*magSqr(U);
 
 if (thermo.dpdt())

applications/solvers/combustion/fireFoam/Make/options

@@ -1,6 +1,8 @@
 EXE_INC = \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/fieldSources/lnInclude \
     -I${LIB_SRC}/meshTools/lnInclude \
+    -I${LIB_SRC}/sampling/lnInclude \
     -I$(LIB_SRC)/turbulenceModels/compressible/turbulenceModel \
     -I$(LIB_SRC)/lagrangian/distributionModels/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/specie/lnInclude \
@@ -29,7 +31,9 @@ EXE_INC = \
 
 EXE_LIBS = \
     -lfiniteVolume \
+    -lfieldSources \
     -lmeshTools \
+    -lsampling \
     -lcompressibleRASModels \
     -lcompressibleLESModels \
     -lspecie \
@@ -46,5 +50,6 @@ EXE_LIBS = \
     -lradiationModels \
     -lsurfaceFilmModels \
     -lpyrolysisModels \
+    -lregionCoupling \
     -llagrangianIntermediate \
     -lODE

applications/solvers/combustion/fireFoam/UEqn.H

@@ -5,10 +5,13 @@
       + turbulence->divDevRhoReff(U)
      ==
         parcels.SU(U)
+      + sources(rho, U)
     );
 
     UEqn.relax();
 
+    sources.constrain(UEqn);
+
     if (pimple.momentumPredictor())
     {
         solve
@@ -24,5 +27,6 @@
             )
         );
 
+        sources.correct(U);
         K = 0.5*magSqr(U);
     }

applications/solvers/combustion/fireFoam/YEEqn.H

@@ -30,11 +30,17 @@ tmp<fv::convectionScheme<scalar> > mvConvection
                 parcels.SYi(i, Yi)
               + surfaceFilm.Srho(i)
               + combustion->R(Yi)
+              + sources(rho, Yi)
             );
 
             YiEqn.relax();
+
+            sources.constrain(YiEqn);
+
             YiEqn.solve(mesh.solver("Yi"));
 
+            sources.correct(Yi);
+
             Yi.max(0.0);
             Yt += Yi;
         }
@@ -69,11 +75,17 @@ tmp<fv::convectionScheme<scalar> > mvConvection
       + radiation->Sh(thermo)
       + parcels.Sh(he)
       + surfaceFilm.Sh()
+      + sources(rho, he)
     );
 
     EEqn.relax();
+
+    sources.constrain(EEqn);
+
     EEqn.solve();
 
+    sources.correct(he);
+
     thermo.correct();
 
     Info<< "min/max(T) = "

applications/solvers/combustion/fireFoam/fireFoam.C

@@ -30,7 +30,6 @@ Description
 
 \*---------------------------------------------------------------------------*/
 
-#include "mapDistribute.H"
 #include "fvCFD.H"
 #include "turbulenceModel.H"
 #include "basicReactingCloud.H"
@@ -41,6 +40,7 @@ Description
 #include "solidChemistryModel.H"
 #include "psiCombustionModel.H"
 #include "pimpleControl.H"
+#include "IObasicSourceList.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -52,6 +52,7 @@ int main(int argc, char *argv[])
     #include "createMesh.H"
     #include "readGravitationalAcceleration.H"
     #include "createFields.H"
+    #include "createSources.H"
     #include "createClouds.H"
     #include "createSurfaceFilmModel.H"
     #include "createPyrolysisModel.H"

applications/solvers/combustion/fireFoam/pEqn.H

@@ -21,6 +21,7 @@ surfaceScalarField phiHbyA
   + phig
 );
 
+sources.relativeFlux(phiHbyA);
 
 while (pimple.correctNonOrthogonal())
 {
@@ -33,8 +34,11 @@ while (pimple.correctNonOrthogonal())
      ==
         parcels.Srho()
       + surfaceFilm.Srho()
+      + sources(psi, p_rgh, rho.name())
     );
 
+    sources.constrain(p_rghEqn);
+
     p_rghEqn.solve(mesh.solver(p_rgh.select(pimple.finalInnerIter())));
 
     if (pimple.finalNonOrthogonalIter())
@@ -42,6 +46,7 @@ while (pimple.correctNonOrthogonal())
         phi = phiHbyA + p_rghEqn.flux();
         U = HbyA + rAU*fvc::reconstruct((p_rghEqn.flux() + phig)/rhorAUf);
         U.correctBoundaryConditions();
+        sources.correct(U);
     }
 }
 

applications/solvers/combustion/fireFoam/rhoEqn.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
@@ -30,14 +30,21 @@ Description
 \*---------------------------------------------------------------------------*/
 
 {
-    solve
+    fvScalarMatrix rhoEqn
     (
         fvm::ddt(rho)
       + fvc::div(phi)
       ==
         parcels.Srho(rho)
       + surfaceFilm.Srho()
+      + sources(rho)
     );
+
+    sources.constrain(rhoEqn);
+
+    rhoEqn.solve();
+
+    sources.correct(rho);
 }
 
 // ************************************************************************* //

applications/solvers/combustion/reactingFoam/Allwmake

@@ -3,6 +3,7 @@ cd ${0%/*} || exit 1    # run from this directory
 set -x
 
 wmake
-wmake buoyantSimpleRadiationFoam
+wmake rhoReactingFoam
+wmake rhoReactingBuoyantFoam
 
 # ----------------------------------------------------------------- end-of-file

applications/solvers/combustion/reactingFoam/EEqn.H

@@ -19,11 +19,17 @@
 //    - fvm::laplacian(turbulence->muEff(), he)  // unit lewis no.
      ==
         reaction->Sh()
+      + sources(rho, he)
     );
 
     EEqn.relax();
+
+    sources.constrain(EEqn);
+
     EEqn.solve();
 
+    sources.correct(he);
+
     thermo.correct();
 
     Info<< "min/max(T) = "

applications/solvers/combustion/reactingFoam/Make/options

@@ -1,14 +1,21 @@
 EXE_INC = \
+    -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/fieldSources/lnInclude \
+    -I$(LIB_SRC)/meshTools/lnInclude \
+    -I$(LIB_SRC)/sampling/lnInclude \
     -I$(LIB_SRC)/turbulenceModels/compressible/turbulenceModel \
     -I$(LIB_SRC)/thermophysicalModels/specie/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/reactionThermo/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/chemistryModel/lnInclude \
     -I$(LIB_SRC)/ODE/lnInclude \
-    -I$(LIB_SRC)/finiteVolume/lnInclude \
     -I$(LIB_SRC)/combustionModels/lnInclude
 
 EXE_LIBS = \
+    -lfiniteVolume \
+    -lfieldSources \
+    -lmeshTools \
+    -lsampling \
     -lcompressibleTurbulenceModel \
     -lcompressibleRASModels \
     -lcompressibleLESModels \
@@ -17,5 +24,4 @@ EXE_LIBS = \
     -lfluidThermophysicalModels \
     -lchemistryModel \
     -lODE \
-    -lfiniteVolume \
     -lcombustionModels

applications/solvers/combustion/reactingFoam/UEqn.H

@@ -5,12 +5,17 @@
       + turbulence->divDevRhoReff(U)
      ==
         rho*g
+      + sources(rho, U)
     );
 
     UEqn.relax();
 
+    sources.constrain(UEqn);
+
     if (pimple.momentumPredictor())
     {
         solve(UEqn == -fvc::grad(p));
+
+        sources.correct(U);
         K = 0.5*magSqr(U);
     }

applications/solvers/combustion/reactingFoam/YEqn.H

@@ -28,11 +28,17 @@ tmp<fv::convectionScheme<scalar> > mvConvection
               - fvm::laplacian(turbulence->muEff(), Yi)
              ==
                 reaction->R(Yi)
+              + sources(rho, Yi)
             );
 
             YiEqn.relax();
+
+            sources.constrain(YiEqn);
+
             YiEqn.solve(mesh.solver("Yi"));
 
+            sources.correct(Yi);
+
             Yi.max(0.0);
             Yt += Yi;
         }

applications/solvers/combustion/reactingFoam/pEqn.H

@@ -16,6 +16,8 @@ if (pimple.transonic())
         )
     );
 
+    sources.relativeFlux(fvc::interpolate(psi), phid);
+
     while (pimple.correctNonOrthogonal())
     {
         fvScalarMatrix pEqn
@@ -23,8 +25,12 @@ if (pimple.transonic())
             fvm::ddt(psi, p)
           + fvm::div(phid, p)
           - fvm::laplacian(rho*rAU, p)
+         ==
+            sources(psi, p, rho.name())
         );
 
+        sources.constrain(pEqn);
+
         pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter())));
 
         if (pimple.finalNonOrthogonalIter())
@@ -45,6 +51,8 @@ else
         )
     );
 
+    sources.relativeFlux(fvc::interpolate(rho), phiHbyA);
+
     while (pimple.correctNonOrthogonal())
     {
         fvScalarMatrix pEqn
@@ -52,8 +60,12 @@ else
             fvm::ddt(psi, p)
           + fvc::div(phiHbyA)
           - fvm::laplacian(rho*rAU, p)
+         ==
+            sources(psi, p, rho.name())
         );
 
+        sources.constrain(pEqn);
+
         pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter())));
 
         if (pimple.finalNonOrthogonalIter())
@@ -68,6 +80,7 @@ else
 
 U = HbyA - rAU*fvc::grad(p);
 U.correctBoundaryConditions();
+sources.correct(U);
 K = 0.5*magSqr(U);
 
 if (thermo.dpdt())

applications/solvers/combustion/reactingFoam/reactingFoam.C

@@ -34,6 +34,7 @@ Description
 #include "psiCombustionModel.H"
 #include "multivariateScheme.H"
 #include "pimpleControl.H"
+#include "IObasicSourceList.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -44,6 +45,7 @@ int main(int argc, char *argv[])
     #include "createMesh.H"
     #include "readGravitationalAcceleration.H"
     #include "createFields.H"
+    #include "createSources.H"
     #include "initContinuityErrs.H"
     #include "readTimeControls.H"
     #include "compressibleCourantNo.H"

applications/solvers/combustion/reactingFoam/rhoReactingBuoyantFoam/Make/files

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

applications/solvers/combustion/reactingFoam/rhoReactingBuoyantFoam/Make/options

@@ -1,17 +1,23 @@
 EXE_INC = \
-    -I../reactingFoam \
+    -I$(FOAM_SOLVERS)/combustion/reactingFoam \
+    -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/fieldSources/lnInclude \
+    -I$(LIB_SRC)/meshTools/lnInclude \
+    -I$(LIB_SRC)/sampling/lnInclude \
     -I$(LIB_SRC)/turbulenceModels/compressible/turbulenceModel \
     -I$(LIB_SRC)/thermophysicalModels/specie/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/reactionThermo/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/chemistryModel/lnInclude \
     -I$(LIB_SRC)/ODE/lnInclude \
-    -I$(LIB_SRC)/finiteVolume/lnInclude \
-    -I$(FOAM_SOLVERS)/combustion/reactingFoam \
     -I$(LIB_SRC)/combustionModels/lnInclude
 
 
 EXE_LIBS = \
+    -lfiniteVolume \
+    -lfieldSources \
+    -lmeshTools \
+    -lsampling \
     -lcompressibleTurbulenceModel \
     -lcompressibleRASModels \
     -lcompressibleLESModels \
@@ -20,5 +26,4 @@ EXE_LIBS = \
     -lfluidThermophysicalModels \
     -lchemistryModel \
     -lODE \
-    -lfiniteVolume \
     -lcombustionModels

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

@@ -0,0 +1,31 @@
+    fvVectorMatrix UEqn
+    (
+        fvm::ddt(rho, U)
+      + fvm::div(phi, U)
+      + turbulence->divDevRhoReff(U)
+     ==
+        sources(rho, U)
+    );
+
+    UEqn.relax();
+
+    sources.constrain(UEqn);
+
+    if (pimple.momentumPredictor())
+    {
+        solve
+        (
+            UEqn
+         ==
+            fvc::reconstruct
+            (
+                (
+                  - ghf*fvc::snGrad(rho)
+                  - fvc::snGrad(p_rgh)
+                )*mesh.magSf()
+            )
+        );
+
+        sources.correct(U);
+        K = 0.5*magSqr(U);
+    }

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

@@ -0,0 +1,124 @@
+Info<< "Creating reaction model\n" << endl;
+
+autoPtr<combustionModels::rhoCombustionModel> reaction
+(
+    combustionModels::rhoCombustionModel::New(mesh)
+);
+
+rhoReactionThermo& thermo = reaction->thermo();
+thermo.validate(args.executable(), "h", "e");
+
+basicMultiComponentMixture& composition = thermo.composition();
+PtrList<volScalarField>& Y = composition.Y();
+
+word inertSpecie(thermo.lookup("inertSpecie"));
+
+volScalarField rho
+(
+    IOobject
+    (
+        "rho",
+        runTime.timeName(),
+        mesh
+    ),
+    thermo.rho()
+);
+
+Info<< "Reading field U\n" << endl;
+volVectorField U
+(
+    IOobject
+    (
+        "U",
+        runTime.timeName(),
+        mesh,
+        IOobject::MUST_READ,
+        IOobject::AUTO_WRITE
+    ),
+    mesh
+);
+
+
+volScalarField& p = thermo.p();
+const volScalarField& psi = thermo.psi();
+const volScalarField& T = thermo.T();
+
+
+#include "compressibleCreatePhi.H"
+
+
+Info << "Creating turbulence model.\n" << nl;
+autoPtr<compressible::turbulenceModel> turbulence
+(
+    compressible::turbulenceModel::New
+    (
+        rho,
+        U,
+        phi,
+        thermo
+    )
+);
+
+// Set the turbulence into the reaction model
+reaction->setTurbulence(turbulence());
+
+
+Info<< "Calculating field g.h\n" << endl;
+volScalarField gh("gh", g & mesh.C());
+surfaceScalarField ghf("ghf", g & mesh.Cf());
+
+Info<< "Reading field p_rgh\n" << endl;
+volScalarField p_rgh
+(
+    IOobject
+    (
+        "p_rgh",
+        runTime.timeName(),
+        mesh,
+        IOobject::MUST_READ,
+        IOobject::AUTO_WRITE
+    ),
+    mesh
+);
+
+// Force p_rgh to be consistent with p
+p_rgh = p - rho*gh;
+
+Info<< "Creating field dpdt\n" << endl;
+volScalarField dpdt
+(
+    IOobject
+    (
+        "dpdt",
+        runTime.timeName(),
+        mesh
+    ),
+    mesh,
+    dimensionedScalar("dpdt", p.dimensions()/dimTime, 0)
+);
+
+Info<< "Creating field kinetic energy K\n" << endl;
+volScalarField K("K", 0.5*magSqr(U));
+
+
+multivariateSurfaceInterpolationScheme<scalar>::fieldTable fields;
+
+forAll(Y, i)
+{
+    fields.add(Y[i]);
+}
+fields.add(thermo.he());
+
+volScalarField dQ
+(
+    IOobject
+    (
+        "dQ",
+        runTime.timeName(),
+        mesh,
+        IOobject::NO_READ,
+        IOobject::AUTO_WRITE
+    ),
+    mesh,
+    dimensionedScalar("dQ", dimEnergy/dimTime, 0.0)
+);

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

@@ -0,0 +1,78 @@
+{
+    rho = thermo.rho();
+
+    // Thermodynamic density needs to be updated by psi*d(p) after the
+    // pressure solution - done in 2 parts. Part 1:
+    thermo.rho() -= psi*p;
+
+    volScalarField rAU(1.0/UEqn.A());
+    surfaceScalarField rhorAUf("Dp", fvc::interpolate(rho*rAU));
+
+    volVectorField HbyA("HbyA", U);
+    HbyA = rAU*UEqn.H();
+
+    surfaceScalarField phig(-rhorAUf*ghf*fvc::snGrad(rho)*mesh.magSf());
+
+    surfaceScalarField phiHbyA
+    (
+        "phiHbyA",
+        fvc::interpolate(rho)
+       *(
+            (fvc::interpolate(U) & mesh.Sf())
+          + fvc::ddtPhiCorr(rAU, rho, U, phi)
+        )
+      + phig
+    );
+
+    sources.relativeFlux(fvc::interpolate(rho), phiHbyA);
+
+    fvScalarMatrix p_rghDDtEqn
+    (
+        fvc::ddt(rho) + psi*correction(fvm::ddt(p_rgh))
+      + fvc::div(phiHbyA)
+      ==
+        sources(psi, p_rgh, rho.name())
+    );
+
+    while (pimple.correctNonOrthogonal())
+    {
+        fvScalarMatrix p_rghEqn
+        (
+            p_rghDDtEqn
+          - fvm::laplacian(rhorAUf, p_rgh)
+        );
+
+        sources.constrain(p_rghEqn);
+
+        p_rghEqn.solve(mesh.solver(p_rgh.select(pimple.finalInnerIter())));
+
+        if (pimple.finalNonOrthogonalIter())
+        {
+            // Calculate the conservative fluxes
+            phi = phiHbyA + p_rghEqn.flux();
+
+            // Explicitly relax pressure for momentum corrector
+            p_rgh.relax();
+
+            // Correct the momentum source with the pressure gradient flux
+            // calculated from the relaxed pressure
+            U = HbyA + rAU*fvc::reconstruct((phig + p_rghEqn.flux())/rhorAUf);
+            U.correctBoundaryConditions();
+            sources.correct(U);
+            K = 0.5*magSqr(U);
+        }
+    }
+
+    p = p_rgh + rho*gh;
+
+    // Second part of thermodynamic density update
+    thermo.rho() += psi*p;
+
+    if (thermo.dpdt())
+    {
+        dpdt = fvc::ddt(p);
+    }
+
+    #include "rhoEqn.H"
+    #include "compressibleContinuityErrs.H"
+}

applications/solvers/combustion/reactingFoam/rhoReactingBuoyantFoam/rhoReactingBuoyantFoam.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -22,25 +22,20 @@ License
     along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
 
 Application
-    rhoPorousMRFPimpleFoam
+    rhoReactingBuoyantFoam
 
 Description
-    Transient solver for laminar or turbulent flow of compressible fluids
+    Solver for combustion with chemical reactions using density based
-    with support for porous media and MRF for HVAC and similar applications.
+    thermodynamics package, using enahanced buoyancy treatment.
-
-    Uses the flexible PIMPLE (PISO-SIMPLE) solution for time-resolved and
-    pseudo-transient simulations.
 
 \*---------------------------------------------------------------------------*/
 
 #include "fvCFD.H"
-#include "psiThermo.H"
+#include "rhoCombustionModel.H"
 #include "turbulenceModel.H"
-#include "bound.H"
+#include "multivariateScheme.H"
-#include "MRFZones.H"
-#include "IOporosityModelList.H"
-#include "IObasicSourceList.H"
 #include "pimpleControl.H"
+#include "IObasicSourceList.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -49,13 +44,16 @@ int main(int argc, char *argv[])
     #include "setRootCase.H"
     #include "createTime.H"
     #include "createMesh.H"
+    #include "readGravitationalAcceleration.H"
+    #include "createFields.H"
+    #include "createSources.H"
+    #include "initContinuityErrs.H"
+    #include "readTimeControls.H"
+    #include "compressibleCourantNo.H"
+    #include "setInitialDeltaT.H"
 
     pimpleControl pimple(mesh);
 
-    #include "createFields.H"
-    #include "createZones.H"
-    #include "initContinuityErrs.H"
-
     // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
     Info<< "\nStarting time loop\n" << endl;
@@ -67,18 +65,15 @@ int main(int argc, char *argv[])
         #include "setDeltaT.H"
 
         runTime++;
-
         Info<< "Time = " << runTime.timeName() << nl << endl;
 
-        if (pimple.nCorrPIMPLE() <= 1)
-        {
         #include "rhoEqn.H"
-        }
 
         // --- Pressure-velocity PIMPLE corrector loop
         while (pimple.loop())
         {
             #include "UEqn.H"
+            #include "YEqn.H"
             #include "EEqn.H"
 
             // --- Pressure corrector loop
@@ -93,6 +88,8 @@ int main(int argc, char *argv[])
             }
         }
 
+        rho = thermo.rho();
+
         runTime.write();
 
         Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"

applications/solvers/combustion/reactingFoam/rhoReactingFoam/Make/options

@@ -0,0 +1,29 @@
+EXE_INC = \
+    -I$(FOAM_SOLVERS)/combustion/reactingFoam \
+    -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/fieldSources/lnInclude \
+    -I$(LIB_SRC)/meshTools/lnInclude \
+    -I$(LIB_SRC)/sampling/lnInclude \
+    -I$(LIB_SRC)/turbulenceModels/compressible/turbulenceModel \
+    -I$(LIB_SRC)/thermophysicalModels/specie/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/reactionThermo/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/chemistryModel/lnInclude \
+    -I$(LIB_SRC)/ODE/lnInclude \
+    -I$(LIB_SRC)/combustionModels/lnInclude
+
+
+EXE_LIBS = \
+    -lfiniteVolume \
+    -lfieldSources \
+    -lmeshTools \
+    -lsampling \
+    -lcompressibleTurbulenceModel \
+    -lcompressibleRASModels \
+    -lcompressibleLESModels \
+    -lreactionThermophysicalModels \
+    -lspecie \
+    -lfluidThermophysicalModels \
+    -lchemistryModel \
+    -lODE \
+    -lcombustionModels

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

@@ -18,6 +18,8 @@
           + fvc::ddtPhiCorr(rAU, rho, U, phi)
         );
 
+        sources.relativeFlux(phiHbyA);
+
         surfaceScalarField phid("phid", fvc::interpolate(thermo.psi())*phiHbyA);
 
         phiHbyA *= fvc::interpolate(rho);
@@ -34,8 +36,12 @@
             (
                 pDDtEqn
               - fvm::laplacian(rho*rAU, p)
+             ==
+                sources(psi, p, rho.name())
             );
 
+            sources.constrain(pEqn);
+
             pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter())));
 
             if (pimple.finalNonOrthogonalIter())
@@ -56,10 +62,14 @@
             )
         );
 
+        sources.relativeFlux(fvc::interpolate(rho), phiHbyA);
+
         fvScalarMatrix pDDtEqn
         (
             fvc::ddt(rho) + psi*correction(fvm::ddt(p))
           + fvc::div(phiHbyA)
+         ==
+            sources(psi, p, rho.name())
         );
 
         while (pimple.correctNonOrthogonal())
@@ -70,6 +80,8 @@
               - fvm::laplacian(rho*rAU, p)
             );
 
+            sources.constrain(pEqn);
+
             pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter())));
 
             if (pimple.finalNonOrthogonalIter())
@@ -87,6 +99,7 @@
 
     U = HbyA - rAU*fvc::grad(p);
     U.correctBoundaryConditions();
+    sources.correct(U);
     K = 0.5*magSqr(U);
 
     if (thermo.dpdt())

applications/solvers/combustion/reactingFoam/rhoReactingFoam/rhoReactingFoam.C

@@ -35,6 +35,7 @@ Description
 #include "turbulenceModel.H"
 #include "multivariateScheme.H"
 #include "pimpleControl.H"
+#include "IObasicSourceList.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -45,6 +46,7 @@ int main(int argc, char *argv[])
     #include "createMesh.H"
     #include "readGravitationalAcceleration.H"
     #include "createFields.H"
+    #include "createSources.H"
     #include "initContinuityErrs.H"
     #include "readTimeControls.H"
     #include "compressibleCourantNo.H"

applications/solvers/compressible/rhoCentralFoam/BCs/T/smoluchowskiJumpTFvPatchScalarField.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
@@ -38,6 +38,10 @@ Foam::smoluchowskiJumpTFvPatchScalarField::smoluchowskiJumpTFvPatchScalarField
 )
 :
     mixedFvPatchScalarField(p, iF),
+    UName_("U"),
+    rhoName_("rho"),
+    psiName_("thermo:psi"),
+    muName_("thermo:mu"),
     accommodationCoeff_(1.0),
     Twall_(p.size(), 0.0),
     gamma_(1.4)
@@ -57,6 +61,10 @@ Foam::smoluchowskiJumpTFvPatchScalarField::smoluchowskiJumpTFvPatchScalarField
 )
 :
     mixedFvPatchScalarField(ptf, p, iF, mapper),
+    UName_(ptf.UName_),
+    rhoName_(ptf.rhoName_),
+    psiName_(ptf.psiName_),
+    muName_(ptf.muName_),
     accommodationCoeff_(ptf.accommodationCoeff_),
     Twall_(ptf.Twall_),
     gamma_(ptf.gamma_)
@@ -71,6 +79,10 @@ Foam::smoluchowskiJumpTFvPatchScalarField::smoluchowskiJumpTFvPatchScalarField
 )
 :
     mixedFvPatchScalarField(p, iF),
+    UName_(dict.lookupOrDefault<word>("U", "U")),
+    rhoName_(dict.lookupOrDefault<word>("rho", "rho")),
+    psiName_(dict.lookupOrDefault<word>("psi", "thermo:psi")),
+    muName_(dict.lookupOrDefault<word>("mu", "thermo:mu")),
     accommodationCoeff_(readScalar(dict.lookup("accommodationCoeff"))),
     Twall_("Twall", dict, p.size()),
     gamma_(dict.lookupOrDefault<scalar>("gamma", 1.4))
@@ -93,7 +105,7 @@ Foam::smoluchowskiJumpTFvPatchScalarField::smoluchowskiJumpTFvPatchScalarField
             dict
         )   << "unphysical accommodationCoeff specified"
             << "(0 < accommodationCoeff <= 1)" << endl
-            << exit(FatalError);
+            << exit(FatalIOError);
     }
 
     if (dict.found("value"))
@@ -159,13 +171,13 @@ void Foam::smoluchowskiJumpTFvPatchScalarField::updateCoeffs()
     }
 
     const fvPatchScalarField& pmu =
-        patch().lookupPatchField<volScalarField, scalar>("mu");
+        patch().lookupPatchField<volScalarField, scalar>(muName_);
     const fvPatchScalarField& prho =
-        patch().lookupPatchField<volScalarField, scalar>("rho");
+        patch().lookupPatchField<volScalarField, scalar>(rhoName_);
     const fvPatchField<scalar>& ppsi =
-        patch().lookupPatchField<volScalarField, scalar>("psi");
+        patch().lookupPatchField<volScalarField, scalar>(psiName_);
     const fvPatchVectorField& pU =
-        patch().lookupPatchField<volVectorField, vector>("U");
+        patch().lookupPatchField<volVectorField, vector>(UName_);
 
     // Prandtl number reading consistent with rhoCentralFoam
     const dictionary& thermophysicalProperties =
@@ -204,6 +216,12 @@ void Foam::smoluchowskiJumpTFvPatchScalarField::updateCoeffs()
 void Foam::smoluchowskiJumpTFvPatchScalarField::write(Ostream& os) const
 {
     fvPatchScalarField::write(os);
+
+    writeEntryIfDifferent<word>(os, "U", "U", UName_);
+    writeEntryIfDifferent<word>(os, "rho", "rho", rhoName_);
+    writeEntryIfDifferent<word>(os, "psi", "thermo:psi", psiName_);
+    writeEntryIfDifferent<word>(os, "mu", "thermo:mu", muName_);
+
     os.writeKeyword("accommodationCoeff")
         << accommodationCoeff_ << token::END_STATEMENT << nl;
     Twall_.writeEntry("Twall", os);

applications/solvers/compressible/rhoCentralFoam/BCs/T/smoluchowskiJumpTFvPatchScalarField.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
@@ -53,13 +53,25 @@ class smoluchowskiJumpTFvPatchScalarField
 
     // Private data
 
-        // Accommodation coefficient
+        //- Velocity field name, default = "U"
+        word UName_;
+
+        //- Density field name, default = "rho"
+        word rhoName_;
+
+        //- Compressibility field name, default = "thermo:psi"
+        word psiName_;
+
+        //- Dynamic viscosity field name, default = "thermo:mu"
+        word muName_;
+
+        //- Accommodation coefficient
         scalar accommodationCoeff_;
 
-        // Wall surface temperature
+        //- Wall surface temperature
         scalarField Twall_;
 
-        // Heat capacity ratio (default 1.4)
+        //- Heat capacity ratio (default 1.4)
         scalar gamma_;
 
 public:

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

@@ -30,20 +30,20 @@ License
 #include "volFields.H"
 #include "fvcGrad.H"
 
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-namespace Foam
-{
-
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
-maxwellSlipUFvPatchVectorField::maxwellSlipUFvPatchVectorField
+Foam::maxwellSlipUFvPatchVectorField::maxwellSlipUFvPatchVectorField
 (
     const fvPatch& p,
     const DimensionedField<vector, volMesh>& iF
 )
 :
     mixedFixedValueSlipFvPatchVectorField(p, iF),
+    TName_("T"),
+    rhoName_("rho"),
+    psiName_("thermo:psi"),
+    muName_("thermo:mu"),
+    tauMCName_("tauMC"),
     accommodationCoeff_(1.0),
     Uwall_(p.size(), vector(0.0, 0.0, 0.0)),
     thermalCreep_(true),
@@ -51,23 +51,28 @@ maxwellSlipUFvPatchVectorField::maxwellSlipUFvPatchVectorField
 {}
 
 
-maxwellSlipUFvPatchVectorField::maxwellSlipUFvPatchVectorField
+Foam::maxwellSlipUFvPatchVectorField::maxwellSlipUFvPatchVectorField
 (
-    const maxwellSlipUFvPatchVectorField& tdpvf,
+    const maxwellSlipUFvPatchVectorField& mspvf,
     const fvPatch& p,
     const DimensionedField<vector, volMesh>& iF,
     const fvPatchFieldMapper& mapper
 )
 :
-    mixedFixedValueSlipFvPatchVectorField(tdpvf, p, iF, mapper),
+    mixedFixedValueSlipFvPatchVectorField(mspvf, p, iF, mapper),
-    accommodationCoeff_(tdpvf.accommodationCoeff_),
+    TName_(mspvf.TName_),
-    Uwall_(tdpvf.Uwall_),
+    rhoName_(mspvf.rhoName_),
-    thermalCreep_(tdpvf.thermalCreep_),
+    psiName_(mspvf.psiName_),
-    curvature_(tdpvf.curvature_)
+    muName_(mspvf.muName_),
+    tauMCName_(mspvf.tauMCName_),
+    accommodationCoeff_(mspvf.accommodationCoeff_),
+    Uwall_(mspvf.Uwall_),
+    thermalCreep_(mspvf.thermalCreep_),
+    curvature_(mspvf.curvature_)
 {}
 
 
-maxwellSlipUFvPatchVectorField::maxwellSlipUFvPatchVectorField
+Foam::maxwellSlipUFvPatchVectorField::maxwellSlipUFvPatchVectorField
 (
     const fvPatch& p,
     const DimensionedField<vector, volMesh>& iF,
@@ -75,6 +80,11 @@ maxwellSlipUFvPatchVectorField::maxwellSlipUFvPatchVectorField
 )
 :
     mixedFixedValueSlipFvPatchVectorField(p, iF),
+    TName_(dict.lookupOrDefault<word>("T", "T")),
+    rhoName_(dict.lookupOrDefault<word>("rho", "rho")),
+    psiName_(dict.lookupOrDefault<word>("psi", "thermo:psi")),
+    muName_(dict.lookupOrDefault<word>("mu", "thermo:mu")),
+    tauMCName_(dict.lookupOrDefault<word>("tauMC", "tauMC")),
     accommodationCoeff_(readScalar(dict.lookup("accommodationCoeff"))),
     Uwall_("Uwall", dict, p.size()),
     thermalCreep_(dict.lookupOrDefault("thermalCreep", true)),
@@ -88,13 +98,16 @@ maxwellSlipUFvPatchVectorField::maxwellSlipUFvPatchVectorField
     {
         FatalIOErrorIn
         (
-            "maxwellSlipUFvPatchScalarField::"
+            "maxwellSlipUFvPatchScalarField::maxwellSlipUFvPatchScalarField"
-            "maxwellSlipUFvPatchScalarField"
+            "("
-            "(const fvPatch&, const scalarField&, const dictionary&)",
+                "const fvPatch&, "
+                "const DimensionedField<vector, volMesh>&, "
+                "const dictionary&"
+            ")",
             dict
         )   << "unphysical accommodationCoeff_ specified"
             << "(0 < accommodationCoeff_ <= 1)" << endl
-            << exit(FatalError);
+            << exit(FatalIOError);
     }
 
     if (dict.found("value"))
@@ -119,23 +132,28 @@ maxwellSlipUFvPatchVectorField::maxwellSlipUFvPatchVectorField
 }
 
 
-maxwellSlipUFvPatchVectorField::maxwellSlipUFvPatchVectorField
+Foam::maxwellSlipUFvPatchVectorField::maxwellSlipUFvPatchVectorField
 (
-    const maxwellSlipUFvPatchVectorField& tdpvf,
+    const maxwellSlipUFvPatchVectorField& mspvf,
     const DimensionedField<vector, volMesh>& iF
 )
 :
-    mixedFixedValueSlipFvPatchVectorField(tdpvf, iF),
+    mixedFixedValueSlipFvPatchVectorField(mspvf, iF),
-    accommodationCoeff_(tdpvf.accommodationCoeff_),
+    TName_(mspvf.TName_),
-    Uwall_(tdpvf.Uwall_),
+    rhoName_(mspvf.rhoName_),
-    thermalCreep_(tdpvf.thermalCreep_),
+    psiName_(mspvf.psiName_),
-    curvature_(tdpvf.curvature_)
+    muName_(mspvf.muName_),
+    tauMCName_(mspvf.tauMCName_),
+    accommodationCoeff_(mspvf.accommodationCoeff_),
+    Uwall_(mspvf.Uwall_),
+    thermalCreep_(mspvf.thermalCreep_),
+    curvature_(mspvf.curvature_)
 {}
 
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-void maxwellSlipUFvPatchVectorField::updateCoeffs()
+void Foam::maxwellSlipUFvPatchVectorField::updateCoeffs()
 {
     if (updated())
     {
@@ -143,11 +161,11 @@ void maxwellSlipUFvPatchVectorField::updateCoeffs()
     }
 
     const fvPatchScalarField& pmu =
-        patch().lookupPatchField<volScalarField, scalar>("mu");
+        patch().lookupPatchField<volScalarField, scalar>(muName_);
     const fvPatchScalarField& prho =
-        patch().lookupPatchField<volScalarField, scalar>("rho");
+        patch().lookupPatchField<volScalarField, scalar>(rhoName_);
     const fvPatchField<scalar>& ppsi =
-        patch().lookupPatchField<volScalarField, scalar>("psi");
+        patch().lookupPatchField<volScalarField, scalar>(psiName_);
 
     Field<scalar> C1
     (
@@ -163,7 +181,7 @@ void maxwellSlipUFvPatchVectorField::updateCoeffs()
     if (thermalCreep_)
     {
         const volScalarField& vsfT =
-            this->db().objectRegistry::lookupObject<volScalarField>("T");
+            this->db().objectRegistry::lookupObject<volScalarField>(TName_);
         label patchi = this->patch().index();
         const fvPatchScalarField& pT = vsfT.boundaryField()[patchi];
         Field<vector> gradpT(fvc::grad(vsfT)().boundaryField()[patchi]);
@@ -175,7 +193,7 @@ void maxwellSlipUFvPatchVectorField::updateCoeffs()
     if (curvature_)
     {
         const fvPatchTensorField& ptauMC =
-            patch().lookupPatchField<volTensorField, tensor>("tauMC");
+            patch().lookupPatchField<volTensorField, tensor>(tauMCName_);
         vectorField n(patch().nf());
 
         refValue() -= C1/prho*transform(I - n*n, (n & ptauMC));
@@ -185,9 +203,15 @@ void maxwellSlipUFvPatchVectorField::updateCoeffs()
 }
 
 
-void maxwellSlipUFvPatchVectorField::write(Ostream& os) const
+void Foam::maxwellSlipUFvPatchVectorField::write(Ostream& os) const
 {
     fvPatchVectorField::write(os);
+    writeEntryIfDifferent<word>(os, "T", "T", TName_);
+    writeEntryIfDifferent<word>(os, "rho", "rho", rhoName_);
+    writeEntryIfDifferent<word>(os, "psi", "thermo:psi", psiName_);
+    writeEntryIfDifferent<word>(os, "mu", "thermo:mu", muName_);
+    writeEntryIfDifferent<word>(os, "tauMC", "tauMC", tauMCName_);
+
     os.writeKeyword("accommodationCoeff")
         << accommodationCoeff_ << token::END_STATEMENT << nl;
     Uwall_.writeEntry("Uwall", os);
@@ -204,14 +228,13 @@ void maxwellSlipUFvPatchVectorField::write(Ostream& os) const
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-makePatchTypeField
+namespace Foam
-(
+{
+    makePatchTypeField
+    (
         fvPatchVectorField,
         maxwellSlipUFvPatchVectorField
-);
+    );
-
+}
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-} // End namespace Foam
 
 // ************************************************************************* //

applications/solvers/compressible/rhoCentralFoam/BCs/U/maxwellSlipUFvPatchVectorField.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
@@ -54,6 +54,21 @@ class maxwellSlipUFvPatchVectorField
 {
     // Private data
 
+        //- Temperature field name, default = "T"
+        word TName_;
+
+        //- Density field name, default = "rho"
+        word rhoName_;
+
+        //- Compressibility field name, default = "thermo:psi"
+        word psiName_;
+
+        //- Dynamic viscosity field name, default = "thermo:mu"
+        word muName_;
+
+        //- tauMC field name, default = "tauMC"
+        word tauMCName_;
+
         // Accommodation coefficient
         scalar accommodationCoeff_;
 

applications/solvers/compressible/rhoCentralFoam/BCs/mixedFixedValueSlip/mixedFixedValueSlipFvPatchField.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,6 +98,7 @@ mixedFixedValueSlipFvPatchField<Type>::mixedFixedValueSlipFvPatchField
     valueFraction_(ptf.valueFraction_)
 {}
 
+
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
 // Map from self

applications/solvers/compressible/rhoCentralFoam/BCs/rho/fixedRhoFvPatchScalarField.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
@@ -27,27 +27,22 @@ License
 #include "addToRunTimeSelectionTable.H"
 #include "fvPatchFieldMapper.H"
 #include "volFields.H"
-#include "surfaceFields.H"
-#include "fvCFD.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-namespace Foam
-{
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
-fixedRhoFvPatchScalarField::fixedRhoFvPatchScalarField
+Foam::fixedRhoFvPatchScalarField::fixedRhoFvPatchScalarField
 (
     const fvPatch& p,
     const DimensionedField<scalar, volMesh>& iF
 )
 :
-    fixedValueFvPatchScalarField(p, iF)
+    fixedValueFvPatchScalarField(p, iF),
+    pName_("p"),
+    psiName_("thermo:psi")
 {}
 
 
-fixedRhoFvPatchScalarField::fixedRhoFvPatchScalarField
+Foam::fixedRhoFvPatchScalarField::fixedRhoFvPatchScalarField
 (
     const fixedRhoFvPatchScalarField& ptf,
     const fvPatch& p,
@@ -55,43 +50,51 @@ fixedRhoFvPatchScalarField::fixedRhoFvPatchScalarField
     const fvPatchFieldMapper& mapper
 )
 :
-    fixedValueFvPatchScalarField(ptf, p, iF, mapper)
+    fixedValueFvPatchScalarField(ptf, p, iF, mapper),
+    pName_(ptf.pName_),
+    psiName_(ptf.psiName_)
 {}
 
 
-fixedRhoFvPatchScalarField::fixedRhoFvPatchScalarField
+Foam::fixedRhoFvPatchScalarField::fixedRhoFvPatchScalarField
 (
     const fvPatch& p,
     const DimensionedField<scalar, volMesh>& iF,
     const dictionary& dict
 )
 :
-    fixedValueFvPatchScalarField(p, iF, dict)
+    fixedValueFvPatchScalarField(p, iF, dict),
+    pName_(dict.lookupOrDefault<word>("p", "p")),
+    psiName_(dict.lookupOrDefault<word>("psi", "thermo:psi"))
 {}
 
 
-fixedRhoFvPatchScalarField::fixedRhoFvPatchScalarField
+Foam::fixedRhoFvPatchScalarField::fixedRhoFvPatchScalarField
 (
-    const fixedRhoFvPatchScalarField& tppsf
+    const fixedRhoFvPatchScalarField& frpsf
 )
 :
-    fixedValueFvPatchScalarField(tppsf)
+    fixedValueFvPatchScalarField(frpsf),
+    pName_(frpsf.pName_),
+    psiName_(frpsf.psiName_)
 {}
 
 
-fixedRhoFvPatchScalarField::fixedRhoFvPatchScalarField
+Foam::fixedRhoFvPatchScalarField::fixedRhoFvPatchScalarField
 (
-    const fixedRhoFvPatchScalarField& tppsf,
+    const fixedRhoFvPatchScalarField& frpsf,
     const DimensionedField<scalar, volMesh>& iF
 )
 :
-    fixedValueFvPatchScalarField(tppsf, iF)
+    fixedValueFvPatchScalarField(frpsf, iF),
+    pName_(frpsf.pName_),
+    psiName_(frpsf.psiName_)
 {}
 
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-void fixedRhoFvPatchScalarField::updateCoeffs()
+void Foam::fixedRhoFvPatchScalarField::updateCoeffs()
 {
     if (updated())
     {
@@ -99,10 +102,10 @@ void fixedRhoFvPatchScalarField::updateCoeffs()
     }
 
     const fvPatchField<scalar>& psip =
-        patch().lookupPatchField<volScalarField, scalar>("psi");
+        patch().lookupPatchField<volScalarField, scalar>(psiName_);
 
     const fvPatchField<scalar>& pp =
-        patch().lookupPatchField<volScalarField, scalar>("p");
+        patch().lookupPatchField<volScalarField, scalar>(pName_);
 
     operator==(psip*pp);
 
@@ -110,16 +113,24 @@ void fixedRhoFvPatchScalarField::updateCoeffs()
 }
 
 
+void Foam::fixedRhoFvPatchScalarField::write(Ostream& os) const
+{
+    fvPatchScalarField::write(os);
+
+    writeEntryIfDifferent<word>(os, "p", "p", this->pName_);
+    writeEntryIfDifferent<word>(os, "psi", "thermo:psi", psiName_);
+    writeEntry("value", os);
+}
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-makePatchTypeField
+namespace Foam
-(
+{
+    makePatchTypeField
+    (
         fvPatchScalarField,
         fixedRhoFvPatchScalarField
-);
+    );
-
+}
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-} // End namespace Foam
 
 // ************************************************************************* //

applications/solvers/compressible/rhoCentralFoam/BCs/rho/fixedRhoFvPatchScalarField.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
@@ -24,9 +24,42 @@ License
 Class
     Foam::fixedRhoFvPatchScalarField
 
+Group
+    grpInletBoundaryConditions
+
 Description
     Foam::fixedRhoFvPatchScalarField
 
+    This boundary condition provides a fixed density inlet condition for
+    compressible solvers, where the density of calculated using:
+
+    \f{
+        \rho = \psi p
+    \f]
+
+    where
+    \vartable
+        p     | pressure [Pa]
+        \rho  | density [kg/m3]
+    \endvartable
+
+
+    \heading Patch usage
+
+    \table
+        Property     | Description             | Required    | Default value
+        pName        | Pressure field name     | no          | p
+        psiName      | Compressibility field name | no       | thermo:psi
+    \endtable
+
+    Example of the boundary condition specification:
+    \verbatim
+    myPatch
+    {
+        type            fixedRho;
+    }
+    \endverbatim
+
 SourceFiles
     fixedRhoFvPatchScalarField.C
 
@@ -51,6 +84,17 @@ class fixedRhoFvPatchScalarField
     public fixedValueFvPatchScalarField
 {
 
+private:
+
+    // Private data
+
+        //- Pressure field name, default = "p"
+        word pName_;
+
+        //- Compressibility field name, default = "thermo:psi"
+        word psiName_;
+
+
 public:
 
     //- Runtime type information
@@ -125,6 +169,10 @@ public:
 
             //- Update the coefficients associated with the patch field
             virtual void updateCoeffs();
+
+
+        //- Write
+        virtual void write(Ostream&) const;
 };
 
 

applications/solvers/compressible/rhoPimpleFoam/Allwmake

@@ -4,7 +4,6 @@ set -x
 
 wmake
 wmake rhoPimplecFoam
-wmake rhoPorousMRFPimpleFoam
+wmake rhoLTSPimpleFoam
-wmake rhoPorousMRFLTSPimpleFoam
 
 # ----------------------------------------------------------------- end-of-file

applications/solvers/compressible/rhoPimpleFoam/EEqn.H

@@ -21,8 +21,12 @@
     );
 
     EEqn.relax();
+
     sources.constrain(EEqn);
+
     EEqn.solve();
 
+    sources.correct(he);
+
     thermo.correct();
 }

applications/solvers/compressible/rhoPimpleFoam/UEqn.H

@@ -5,6 +5,8 @@ tmp<fvVectorMatrix> UEqn
     fvm::ddt(rho, U)
   + fvm::div(phi, U)
   + turbulence->divDevRhoReff(U)
+ ==
+    sources(rho, U)
 );
 
 UEqn().relax();
@@ -13,6 +15,8 @@ sources.constrain(UEqn());
 
 if (pimple.momentumPredictor())
 {
-    solve(UEqn() == -fvc::grad(p) + sources(rho, U));
+    solve(UEqn() == -fvc::grad(p));
+
+    sources.correct(U);
     K = 0.5*magSqr(U);
 }

applications/solvers/compressible/rhoPimpleFoam/createFields.H

@@ -69,6 +69,3 @@
 
     Info<< "Creating field kinetic energy K\n" << endl;
     volScalarField K("K", 0.5*magSqr(U));
-
-    Info<< "Creating sources\n" << endl;
-    IObasicSourceList sources(mesh);

applications/solvers/compressible/rhoPimpleFoam/pEqn.H

@@ -5,7 +5,7 @@ rho.relax();
 
 volScalarField rAU(1.0/UEqn().A());
 volVectorField HbyA("HbyA", U);
-HbyA = rAU*(UEqn() == sources(rho, U))().H();
+HbyA = rAU*UEqn().H();
 
 if (pimple.nCorrPISO() <= 1)
 {
@@ -24,6 +24,8 @@ if (pimple.transonic())
         )
     );
 
+    sources.relativeFlux(fvc::interpolate(psi), phid);
+
     volScalarField Dp("Dp", rho*rAU);
 
     while (pimple.correctNonOrthogonal())
@@ -37,7 +39,7 @@ if (pimple.transonic())
             sources(psi, p, rho.name())
         );
 
-        sources.constrain(pEqn, rho.name());
+        sources.constrain(pEqn);
 
         pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter())));
 
@@ -59,6 +61,8 @@ else
         )
     );
 
+    sources.relativeFlux(fvc::interpolate(rho), phiHbyA);
+
     volScalarField Dp("Dp", rho*rAU);
 
     while (pimple.correctNonOrthogonal())
@@ -73,7 +77,7 @@ else
             sources(psi, p, rho.name())
         );
 
-        sources.constrain(pEqn, rho.name());
+        sources.constrain(pEqn);
 
         pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter())));
 

applications/solvers/compressible/rhoPimpleFoam/rhoLTSPimpleFoam/Make/files

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

applications/solvers/compressible/rhoPimpleFoam/rhoLTSPimpleFoam/rhoLTSPimpleFoam.C

@@ -26,7 +26,7 @@ Application
 
 Description
     Transient solver for laminar or turbulent flow of compressible fluids
-    with support for porous media and MRF for HVAC and similar applications.
+    with support for run-time selectable sources, e.g. MRF, explicit porosity.
 
     Uses the flexible PIMPLE (PISO-SIMPLE) solution for time-resolved and
     pseudo-transient simulations with support for local time-stepping for
@@ -37,8 +37,6 @@ Description
 #include "fvCFD.H"
 #include "psiThermo.H"
 #include "turbulenceModel.H"
-#include "MRFZones.H"
-#include "IOporosityModelList.H"
 #include "IObasicSourceList.H"
 #include "fvcSmooth.H"
 #include "pimpleControl.H"
@@ -56,7 +54,7 @@ int main(int argc, char *argv[])
 
     #include "setInitialrDeltaT.H"
     #include "createFields.H"
-    #include "createZones.H"
+    #include "createSources.H"
     #include "initContinuityErrs.H"
 
     // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

applications/solvers/compressible/rhoPimpleFoam/rhoPimpleFoam.C

@@ -51,6 +51,7 @@ int main(int argc, char *argv[])
     pimpleControl pimple(mesh);
 
     #include "createFields.H"
+    #include "createSources.H"
     #include "initContinuityErrs.H"
 
     // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

applications/solvers/compressible/rhoPimpleFoam/rhoPimplecFoam/pEqn.H

@@ -51,7 +51,7 @@ if (pimple.transonic())
         // Relax the pressure equation to maintain diagonal dominance
         pEqn.relax();
 
-        sources.constrain(pEqn, rho.name());
+        sources.constrain(pEqn);
 
         pEqn.solve();
 
@@ -89,7 +89,7 @@ else
             sources(psi, p, rho.name())
         );
 
-        sources.constrain(pEqn, rho.name());
+        sources.constrain(pEqn);
 
         pEqn.solve();
 

applications/solvers/compressible/rhoPimpleFoam/rhoPimplecFoam/rhoPimplecFoam.C

@@ -51,6 +51,7 @@ int main(int argc, char *argv[])
     pimpleControl pimple(mesh);
 
     #include "createFields.H"
+    #include "createSources.H"
     #include "initContinuityErrs.H"
 
     // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

applications/solvers/compressible/rhoPimpleFoam/rhoPorousMRFLTSPimpleFoam/Make/files

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

applications/solvers/compressible/rhoPimpleFoam/rhoPorousMRFPimpleFoam/Make/files

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

applications/solvers/compressible/rhoPimpleFoam/rhoPorousMRFPimpleFoam/Make/options

@@ -1,20 +0,0 @@
-EXE_INC = \
-    -I.. \
-    -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
-    -I$(LIB_SRC)/turbulenceModels/compressible/turbulenceModel \
-    -I$(LIB_SRC)/finiteVolume/cfdTools \
-    -I$(LIB_SRC)/finiteVolume/lnInclude \
-    -I$(LIB_SRC)/meshTools/lnInclude \
-    -I$(LIB_SRC)/sampling/lnInclude \
-    -I$(LIB_SRC)/fieldSources/lnInclude
-
-EXE_LIBS = \
-    -lfluidThermophysicalModels \
-    -lspecie \
-    -lcompressibleTurbulenceModel \
-    -lcompressibleRASModels \
-    -lcompressibleLESModels \
-    -lfiniteVolume \
-    -lmeshTools \
-    -lsampling \
-    -lfieldSources

applications/solvers/compressible/rhoPimpleFoam/rhoPorousMRFPimpleFoam/UEqn.H

@@ -1,22 +0,0 @@
-// Solve the Momentum equation
-
-tmp<fvVectorMatrix> UEqn
-(
-    //pZones.ddt(rho, U)
-    fvm::ddt(rho, U)
-  + fvm::div(phi, U)
-  + turbulence->divDevRhoReff(U)
-);
-
-UEqn().relax();
-
-mrfZones.addCoriolis(rho, UEqn());
-
-pZones.addResistance(UEqn());
-
-sources.constrain(UEqn());
-
-if (pimple.momentumPredictor())
-{
-    solve(UEqn() == -fvc::grad(p) + sources(rho, U));
-}

applications/solvers/compressible/rhoPimpleFoam/rhoPorousMRFPimpleFoam/createZones.H

@@ -1,5 +0,0 @@
-    MRFZones mrfZones(mesh);
-    mrfZones.correctBoundaryVelocity(U);
-
-    IOporosityModelList pZones(mesh);
-    Switch pressureImplicitPorosity(false);

applications/solvers/compressible/rhoPimpleFoam/rhoPorousMRFPimpleFoam/pEqn.H

@@ -1,111 +0,0 @@
-rho = thermo.rho();
-rho = max(rho, rhoMin);
-rho = min(rho, rhoMax);
-rho.relax();
-
-volScalarField rAU(1.0/UEqn().A());
-volVectorField HbyA("HbyA", U);
-HbyA = rAU*(UEqn() == sources(rho, U))().H();
-
-if (pimple.nCorrPISO() <= 1)
-{
-    UEqn.clear();
-}
-
-if (pimple.transonic())
-{
-    surfaceScalarField phid
-    (
-        "phid",
-        fvc::interpolate(psi)
-       *(
-            (fvc::interpolate(HbyA) & mesh.Sf())
-          + fvc::ddtPhiCorr(rAU, rho, U, phi)
-        )
-    );
-    mrfZones.relativeFlux(fvc::interpolate(psi), phid);
-
-    volScalarField Dp("Dp", rho*rAU);
-
-    while (pimple.correctNonOrthogonal())
-    {
-        fvScalarMatrix pEqn
-        (
-            fvm::ddt(psi, p)
-          + fvm::div(phid, p)
-          - fvm::laplacian(Dp, p)
-         ==
-            sources(psi, p, rho.name())
-        );
-
-        sources.constrain(pEqn, rho.name());
-
-        pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter())));
-
-        if (pimple.finalNonOrthogonalIter())
-        {
-            phi == pEqn.flux();
-        }
-    }
-}
-else
-{
-    surfaceScalarField phiHbyA
-    (
-        "phiHbyA",
-        fvc::interpolate(rho)
-       *(
-            (fvc::interpolate(HbyA) & mesh.Sf())
-          + fvc::ddtPhiCorr(rAU, rho, U, phi)
-        )
-    );
-
-    mrfZones.relativeFlux(fvc::interpolate(rho), phiHbyA);
-
-    volScalarField Dp("Dp", rho*rAU);
-
-    while (pimple.correctNonOrthogonal())
-    {
-        // Pressure corrector
-        fvScalarMatrix pEqn
-        (
-            fvm::ddt(psi, p)
-          + fvc::div(phiHbyA)
-          - fvm::laplacian(Dp, p)
-         ==
-            sources(psi, p, rho.name())
-        );
-
-        sources.constrain(pEqn, rho.name());
-
-        pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter())));
-
-        if (pimple.finalNonOrthogonalIter())
-        {
-            phi = phiHbyA + pEqn.flux();
-        }
-    }
-}
-
-#include "rhoEqn.H"
-#include "compressibleContinuityErrs.H"
-
-// Explicitly relax pressure for momentum corrector
-p.relax();
-
-rho = thermo.rho();
-rho = max(rho, rhoMin);
-rho = min(rho, rhoMax);
-rho.relax();
-Info<< "rho max/min : " << max(rho).value()
-    << " " << min(rho).value() << endl;
-
-U = HbyA - rAU*fvc::grad(p);
-U.correctBoundaryConditions();
-sources.correct(U);
-K = 0.5*magSqr(U);
-
-if (thermo.dpdt())
-{
-    dpdt = fvc::ddt(p);
-}

applications/solvers/compressible/rhoSimpleFoam/Allwclean

@@ -3,7 +3,7 @@ cd ${0%/*} || exit 1    # run from this directory
 set -x
 
 wclean
-wclean rhoPorousMRFSimpleFoam
+wclean rhoPorousSimpleFoam
 wclean rhoSimplecFoam
 
 # ----------------------------------------------------------------- end-of-file

applications/solvers/compressible/rhoSimpleFoam/Allwmake

@@ -3,7 +3,7 @@ cd ${0%/*} || exit 1    # run from this directory
 set -x
 
 wmake
-wmake rhoPorousMRFSimpleFoam
+wmake rhoPorousSimpleFoam
 wmake rhoSimplecFoam
 
 # ----------------------------------------------------------------- end-of-file

applications/solvers/compressible/rhoSimpleFoam/EEqn.H

@@ -15,8 +15,12 @@
     );
 
     EEqn.relax();
+
     sources.constrain(EEqn);
+
     EEqn.solve();
 
+    sources.correct(he);
+
     thermo.correct();
 }

applications/solvers/compressible/rhoSimpleFoam/UEqn.H

@@ -4,8 +4,14 @@
     (
         fvm::div(phi, U)
       + turbulence->divDevRhoReff(U)
+     ==
+        sources(rho, U)
     );
 
     UEqn().relax();
 
+    sources.constrain(UEqn());
+
     solve(UEqn() == -fvc::grad(p));
+
+    sources.correct(U);

applications/solvers/compressible/rhoSimpleFoam/createFields.H

@@ -60,6 +60,3 @@
     );
 
     dimensionedScalar initialMass = fvc::domainIntegrate(rho);
-    
-    Info<< "Creating sources\n" << endl;
-    IObasicSourceList sources(mesh);

applications/solvers/compressible/rhoSimpleFoam/pEqn.H

@@ -6,7 +6,7 @@
 
     volScalarField rAU(1.0/UEqn().A());
     volVectorField HbyA("HbyA", U);
-    HbyA = rAU*(UEqn() == sources(rho, U))().H();
+    HbyA = rAU*UEqn().H();
 
     UEqn.clear();
 
@@ -20,6 +20,8 @@
             fvc::interpolate(psi)*(fvc::interpolate(HbyA) & mesh.Sf())
         );
 
+        sources.relativeFlux(fvc::interpolate(psi), phid);
+
         while (simple.correctNonOrthogonal())
         {
             fvScalarMatrix pEqn
@@ -33,7 +35,7 @@
             // Relax the pressure equation to ensure diagonal-dominance
             pEqn.relax();
 
-            sources.constrain(pEqn, rho.name());
+            sources.constrain(pEqn);
 
             pEqn.setReference(pRefCell, pRefValue);
 
@@ -53,6 +55,8 @@
             fvc::interpolate(rho)*(fvc::interpolate(HbyA) & mesh.Sf())
         );
 
+        sources.relativeFlux(fvc::interpolate(rho), phiHbyA);
+
         closedVolume = adjustPhi(phiHbyA, U, p);
 
         while (simple.correctNonOrthogonal())
@@ -67,7 +71,7 @@
 
             pEqn.setReference(pRefCell, pRefValue);
 
-            sources.constrain(pEqn, rho.name());
+            sources.constrain(pEqn);
 
             pEqn.solve();
 

applications/solvers/compressible/rhoSimpleFoam/rhoPorousMRFSimpleFoam/Make/files

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

applications/solvers/compressible/rhoSimpleFoam/rhoPorousSimpleFoam/EEqn.H

@@ -15,8 +15,12 @@
     );
 
     EEqn.relax();
+
     sources.constrain(EEqn);
+
     EEqn.solve();
 
+    sources.correct(he);
+
     thermo.correct();
 }

applications/solvers/compressible/rhoSimpleFoam/rhoPorousSimpleFoam/Make/files

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

applications/solvers/compressible/rhoSimpleFoam/rhoPorousSimpleFoam/UEqn.H

@@ -33,6 +33,8 @@
             U = trTU() & ((UEqn() == sources(rho, U))().H() - gradp);
         }
         U.correctBoundaryConditions();
+
+        sources.correct(U);
     }
     else
     {
@@ -42,6 +44,8 @@
 
         solve(UEqn() == -fvc::grad(p) + sources(rho, U));
 
+        sources.correct(U);
+
         trAU = 1.0/UEqn().A();
         trAU().rename("rAU");
     }

applications/solvers/compressible/rhoSimpleFoam/rhoPorousSimpleFoam/createFields.H

@@ -59,6 +59,3 @@
     );
 
     dimensionedScalar initialMass = fvc::domainIntegrate(rho);
-
-    Info<< "Creating sources\n" << endl;
-    IObasicSourceList sources(mesh);

applications/solvers/compressible/rhoSimpleFoam/rhoPorousSimpleFoam/createZones.H

@@ -1,4 +1,4 @@
-    MRFZones mrfZones(mesh);
+    IOMRFZoneList mrfZones(mesh);
     mrfZones.correctBoundaryVelocity(U);
 
     IOporosityModelList pZones(mesh);

applications/solvers/compressible/rhoSimpleFoam/rhoPorousSimpleFoam/rhoPorousSimpleFoam.C

@@ -34,7 +34,7 @@ Description
 #include "fvCFD.H"
 #include "rhoThermo.H"
 #include "RASModel.H"
-#include "MRFZones.H"
+#include "IOMRFZoneList.H"
 #include "IObasicSourceList.H"
 #include "IOporosityModelList.H"
 #include "simpleControl.H"
@@ -50,6 +50,7 @@ int main(int argc, char *argv[])
     simpleControl simple(mesh);
 
     #include "createFields.H"
+    #include "createSources.H"
     #include "createZones.H"
     #include "initContinuityErrs.H"
 

applications/solvers/compressible/rhoSimpleFoam/rhoSimpleFoam.C

@@ -47,6 +47,7 @@ int main(int argc, char *argv[])
     simpleControl simple(mesh);
 
     #include "createFields.H"
+    #include "createSources.H"
     #include "initContinuityErrs.H"
 
     // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

applications/solvers/compressible/rhoSimpleFoam/rhoSimplecFoam/pEqn.H

@@ -45,7 +45,7 @@ if (simple.transonic())
         // Relax the pressure equation to maintain diagonal dominance
         pEqn.relax();
 
-        sources.constrain(pEqn, rho.name());
+        sources.constrain(pEqn);
 
         pEqn.setReference(pRefCell, pRefValue);
 
@@ -82,7 +82,7 @@ else
             sources(psi, p, rho.name())
         );
 
-        sources.constrain(pEqn, rho.name());
+        sources.constrain(pEqn);
 
         pEqn.setReference(pRefCell, pRefValue);
 

applications/solvers/compressible/rhoSimpleFoam/rhoSimplecFoam/rhoSimplecFoam.C

@@ -49,6 +49,7 @@ int main(int argc, char *argv[])
     simpleControl simple(mesh);
 
     #include "createFields.H"
+    #include "createSources.H"
     #include "initContinuityErrs.H"
 
     // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

applications/solvers/heatTransfer/buoyantPimpleFoam/EEqn.H

@@ -0,0 +1,32 @@
+{
+    volScalarField& he = thermo.he();
+
+    fvScalarMatrix EEqn
+    (
+        fvm::ddt(rho, he) + fvm::div(phi, he)
+      + fvc::ddt(rho, K) + fvc::div(phi, K)
+      + (
+            he.name() == "e"
+          ? fvc::div
+            (
+                fvc::absolute(phi/fvc::interpolate(rho), U),
+                p,
+                "div(phiv,p)"
+            )
+          : -dpdt
+        )
+      - fvm::laplacian(turbulence->alphaEff(), he)
+     ==
+        sources(rho, he)
+    );
+
+    EEqn.relax();
+
+    sources.constrain(EEqn);
+
+    EEqn.solve();
+
+    sources.correct(he);
+
+    thermo.correct();
+}

applications/solvers/heatTransfer/buoyantPimpleFoam/Make/options

@@ -1,21 +1,21 @@
 EXE_INC = \
-    -I../../compressible/rhoPimpleFoam \
-    -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
-    -I$(LIB_SRC)/turbulenceModels/compressible/turbulenceModel \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
     -I$(LIB_SRC)/sampling/lnInclude \
     -I$(LIB_SRC)/meshTools/lnInclude \
     -I$(LIB_SRC)/fieldSources/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/radiationModels/lnInclude \
+    -I$(LIB_SRC)/turbulenceModels/compressible/turbulenceModel
 
 EXE_LIBS = \
-    -lmeshTools \
-    -lfluidThermophysicalModels \
-    -lspecie \
-    -lcompressibleTurbulenceModel \
-    -lcompressibleRASModels \
-    -lcompressibleLESModels \
     -lfiniteVolume \
     -lsampling \
     -lmeshTools \
-    -lfieldSources
+    -lfieldSources \
+    -lfluidThermophysicalModels \
+    -lradiationModels \
+    -lspecie \
+    -lcompressibleTurbulenceModel \
+    -lcompressibleRASModels \
+    -lcompressibleLESModels
 

applications/solvers/heatTransfer/buoyantPimpleFoam/UEqn.H

@@ -5,10 +5,14 @@
         fvm::ddt(rho, U)
       + fvm::div(phi, U)
       + turbulence->divDevRhoReff(U)
+     ==
+        sources(rho, U)
     );
 
     UEqn.relax();
 
+    sources.constrain(UEqn);
+
     if (pimple.momentumPredictor())
     {
         solve
@@ -22,7 +26,8 @@
                   - fvc::snGrad(p_rgh)
                 )*mesh.magSf()
             )
-          + sources(rho, U)
         );
+
+        sources.correct(U);
         K = 0.5*magSqr(U);
     }

applications/solvers/heatTransfer/buoyantPimpleFoam/buoyantPimpleFoam.C

@@ -36,7 +36,7 @@ Description
 #include "fvCFD.H"
 #include "rhoThermo.H"
 #include "turbulenceModel.H"
-#include "fixedGradientFvPatchFields.H"
+#include "radiationModel.H"
 #include "IObasicSourceList.H"
 #include "pimpleControl.H"
 
@@ -49,6 +49,8 @@ int main(int argc, char *argv[])
     #include "createMesh.H"
     #include "readGravitationalAcceleration.H"
     #include "createFields.H"
+    #include "createSources.H"
+    #include "createRadiationModel.H"
     #include "initContinuityErrs.H"
     #include "readTimeControls.H"
     #include "compressibleCourantNo.H"

applications/solvers/heatTransfer/buoyantPimpleFoam/createFields.H

@@ -1,9 +1,6 @@
     Info<< "Reading thermophysical properties\n" << endl;
 
-    autoPtr<rhoThermo> pThermo
+    autoPtr<rhoThermo> pThermo(rhoThermo::New(mesh));
-    (
-        rhoThermo::New(mesh)
-    );
     rhoThermo& thermo = pThermo();
     thermo.validate(args.executable(), "h", "e");
 
@@ -89,6 +86,3 @@
 
     Info<< "Creating field kinetic energy K\n" << endl;
     volScalarField K("K", 0.5*magSqr(U));
-
-    Info<< "Creating sources\n" << endl;
-    IObasicSourceList sources(mesh);

applications/solvers/heatTransfer/buoyantPimpleFoam/pEqn.H

@@ -9,7 +9,7 @@
     surfaceScalarField rhorAUf("Dp", fvc::interpolate(rho*rAU));
 
     volVectorField HbyA("HbyA", U);
-    HbyA = rAU*(UEqn == sources(rho, U))().H();
+    HbyA = rAU*UEqn.H();
 
     surfaceScalarField phig(-rhorAUf*ghf*fvc::snGrad(rho)*mesh.magSf());
 
@@ -24,6 +24,8 @@
       + phig
     );
 
+    sources.relativeFlux(fvc::interpolate(rho), phiHbyA);
+
     fvScalarMatrix p_rghDDtEqn
     (
         fvc::ddt(rho) + psi*correction(fvm::ddt(p_rgh))
@@ -40,7 +42,7 @@
           - fvm::laplacian(rhorAUf, p_rgh)
         );
 
-        sources.constrain(p_rghEqn, rho.name());
+        sources.constrain(p_rghEqn);
 
         p_rghEqn.solve(mesh.solver(p_rgh.select(pimple.finalInnerIter())));
 

applications/solvers/heatTransfer/buoyantSimpleFoam/EEqn.H

@@ -10,10 +10,19 @@
           : fvc::div(phi, volScalarField("K", 0.5*magSqr(U)))
         )
       - fvm::laplacian(turbulence->alphaEff(), he)
+     ==
+        radiation->Sh(thermo)
+      + sources(rho, he)
     );
 
     EEqn.relax();
+
+    sources.constrain(EEqn);
+
     EEqn.solve();
 
+    sources.correct(he);
+
     thermo.correct();
+    radiation->correct();
 }

applications/solvers/heatTransfer/buoyantSimpleFoam/Make/options

@@ -1,14 +1,22 @@
 EXE_INC = \
+    -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/fieldSources/lnInclude \
+    -I$(LIB_SRC)/sampling/lnInclude \
+    -I$(LIB_SRC)/meshTools/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/radiationModels/lnInclude \
     -I$(LIB_SRC)/turbulenceModels \
     -I$(LIB_SRC)/turbulenceModels/compressible/RAS/lnInclude \
-    -I$(LIB_SRC)/finiteVolume/cfdTools \
+    -I$(LIB_SRC)/turbulenceModels/RAS \
-    -I$(LIB_SRC)/finiteVolume/lnInclude
 
 EXE_LIBS = \
+    -lfiniteVolume \
+    -lfieldSources \
+    -lsampling \
     -lmeshTools \
     -lfluidThermophysicalModels \
     -lspecie \
+    -lradiationModels \
     -lcompressibleTurbulenceModel \
     -lcompressibleRASModels \
-    -lfiniteVolume
+    -lmeshTools

applications/solvers/heatTransfer/buoyantSimpleFoam/UEqn.H

@@ -4,10 +4,14 @@
     (
         fvm::div(phi, U)
       + turbulence->divDevRhoReff(U)
+     ==
+        sources(rho, U)
     );
 
     UEqn().relax();
 
+    sources.constrain(UEqn());
+
     if (simple.momentumPredictor())
     {
         solve
@@ -22,4 +26,6 @@
                 )*mesh.magSf()
             )
         );
+
+        sources.correct(U);
     }

applications/solvers/heatTransfer/buoyantSimpleFoam/buoyantSimpleFoam.C

@@ -22,18 +22,20 @@ License
     along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
 
 Application
-    buoyantSimpleFoam
+    buoyantSimpleRadiationFoam
 
 Description
-    Steady-state solver for buoyant, turbulent flow of compressible fluids
+    Steady-state solver for buoyant, turbulent flow of compressible fluids,
+    including radiation, for ventilation and heat-transfer.
 
 \*---------------------------------------------------------------------------*/
 
 #include "fvCFD.H"
 #include "psiThermo.H"
 #include "RASModel.H"
-#include "fixedGradientFvPatchFields.H"
+#include "radiationModel.H"
 #include "simpleControl.H"
+#include "IObasicSourceList.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -44,6 +46,8 @@ int main(int argc, char *argv[])
     #include "createMesh.H"
     #include "readGravitationalAcceleration.H"
     #include "createFields.H"
+    #include "createSources.H"
+    #include "createRadiationModel.H"
     #include "initContinuityErrs.H"
 
     simpleControl simple(mesh);