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

applications/solvers/combustion/XiFoam/EaEqn.H

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

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)/meshTools/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
-    -lmeshTools \
-    -lfiniteVolume

applications/solvers/combustion/XiFoam/UEqn.H

@@ -5,10 +5,13 @@
       + turbulence->divDevRhoReff(U)
      ==
         rho*g
+      + sources(rho, U)
     );
 
     UEqn.relax();
 
+    sources.constrain(UEqn);
+
     if (pimple.momentumPredictor())
     {
         solve(UEqn == -fvc::grad(p));

applications/solvers/combustion/XiFoam/XiFoam.C

@@ -56,6 +56,7 @@ Description
 #include "ignition.H"
 #include "Switch.H"
 #include "pimpleControl.H"
+#include "IObasicSourceList.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,6 +66,9 @@ if (ign.ignited())
     // Solve for b
     // ~~~~~~~~~~~
     bEqn.relax();
+
+    sources.constrain(bEqn);
+
     bEqn.solve();
 
     Info<< "min(b) = " << min(b).value() << endl;
@@ -153,9 +158,13 @@ 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();
 
         // Limit the maximum Su
@@ -234,9 +243,13 @@ if (ign.ignited())
                 ),
                 Xi
             )
+          + sources(rho, Xi)
         );
 
         XiEqn.relax();
+
+        sources.constrain(XiEqn);
+
         XiEqn.solve();
 
         // Correct boundedness of Xi

applications/solvers/combustion/XiFoam/createFields.H

@@ -138,3 +138,5 @@
     fields.add(b);
     fields.add(thermo.he());
     fields.add(thermo.heu());
+
+    IObasicSourceList sources(mesh);

applications/solvers/combustion/XiFoam/ftEqn.H

@@ -13,10 +13,16 @@ 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();
 }

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/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,6 +3,8 @@
         fvm::ddt(rho, U)
       + fvm::div(phi, U)
       + turbulence->divDevRhoReff(U)
+     ==
+        sources(rho, U)
     );
 
     if (pimple.momentumPredictor())

applications/solvers/combustion/engineFoam/engineFoam.C

@@ -59,6 +59,7 @@ Description
 #include "OFstream.H"
 #include "mathematicalConstants.H"
 #include "pimpleControl.H"
+#include "IObasicSourceList.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 \

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

applications/solvers/combustion/fireFoam/YEEqn.H

@@ -30,9 +30,13 @@ 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"));
 
             Yi.max(0.0);
@@ -69,9 +73,13 @@ tmp<fv::convectionScheme<scalar> > mvConvection
       + radiation->Sh(thermo)
       + parcels.Sh(he)
       + surfaceFilm.Sh()
+      + sources(rho, he)
     );
 
     EEqn.relax();
+
+    sources.constrain(EEqn);
+
     EEqn.solve();
 
     thermo.correct();

applications/solvers/combustion/fireFoam/createFields.H

@@ -146,3 +146,6 @@
     (
         additionalControlsDict.lookup("solvePrimaryRegion")
     );
+
+    IObasicSourceList sources(mesh);
+

applications/solvers/combustion/fireFoam/fireFoam.C

@@ -40,6 +40,7 @@ Description
 #include "solidChemistryModel.H"
 #include "psiCombustionModel.H"
 #include "pimpleControl.H"
+#include "IObasicSourceList.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

@@ -30,14 +30,19 @@ Description
 \*---------------------------------------------------------------------------*/
 
 {
-    solve
+    fvScalarMatrix rhoEqn
     (
         fvm::ddt(rho)
       + fvc::div(phi)
       ==
         parcels.Srho(rho)
       + surfaceFilm.Srho()
+      + sources(rho)
     );
+
+    sources.constrain(rhoEqn);
+
+    rhoEqn.solve();
 }
 
 // ************************************************************************* //

applications/solvers/combustion/reactingFoam/EEqn.H

@@ -19,9 +19,13 @@
 //    - fvm::laplacian(turbulence->muEff(), he)  // unit lewis no.
      ==
         reaction->Sh()
+      + sources(rho, he)
     );
 
     EEqn.relax();
+
+    sources.constrain(EEqn);
+
     EEqn.solve();
 
     thermo.correct();

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,10 +5,13 @@
       + turbulence->divDevRhoReff(U)
      ==
         rho*g
+      + sources(rho, U)
     );
 
     UEqn.relax();
 
+    sources.constrain(UEqn);
+
     if (pimple.momentumPredictor())
     {
         solve(UEqn == -fvc::grad(p));

applications/solvers/combustion/reactingFoam/YEqn.H

@@ -28,9 +28,13 @@ 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"));
 
             Yi.max(0.0);

applications/solvers/combustion/reactingFoam/createFields.H

@@ -99,3 +99,5 @@ volScalarField dQ
     mesh,
     dimensionedScalar("dQ", dimEnergy/dimTime, 0.0)
 );
+
+IObasicSourceList sources(mesh);

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"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 

applications/solvers/combustion/rhoReactingFoam/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/rhoReactingFoam/createFields.H

@@ -101,3 +101,5 @@ volScalarField dQ
     mesh,
     dimensionedScalar("dQ", dimEnergy/dimTime, 0.0)
 );
+
+IObasicSourceList sources(mesh);

applications/solvers/combustion/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/rhoReactingFoam/rhoReactingFoam.C

@@ -35,6 +35,7 @@ Description
 #include "turbulenceModel.H"
 #include "multivariateScheme.H"
 #include "pimpleControl.H"
+#include "IObasicSourceList.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 

applications/solvers/compressible/rhoCentralFoam/BCs/T/smoluchowskiJumpTFvPatchScalarField.C

@@ -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))
@@ -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

@@ -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,9 +98,12 @@ 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
@@ -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

@@ -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/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/pEqn.H

@@ -39,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())));
 
@@ -77,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/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/rhoSimpleFoam/pEqn.H

@@ -35,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);
 
@@ -71,7 +71,7 @@
 
             pEqn.setReference(pRefCell, pRefValue);
 
-            sources.constrain(pEqn, rho.name());
+            sources.constrain(pEqn);
 
             pEqn.solve();
 

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/heatTransfer/buoyantPimpleFoam/EEqn.H

@@ -0,0 +1,28 @@
+{
+    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();
+
+    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,6 @@
                   - fvc::snGrad(p_rgh)
                 )*mesh.magSf()
             )
-          + sources(rho, 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,7 @@ int main(int argc, char *argv[])
     #include "createMesh.H"
     #include "readGravitationalAcceleration.H"
     #include "createFields.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");
 

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/Allwmake

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

applications/solvers/heatTransfer/buoyantSimpleFoam/EEqn.H

@@ -10,10 +10,17 @@
           : 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();
 
     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,11 +4,13 @@
     (
         fvm::div(phi, U)
       + turbulence->divDevRhoReff(U)
+     ==
+        sources(rho, U)
     );
 
     UEqn().relax();
 
-    mrfZones.addCoriolis(rho, UEqn());
+    sources.constrain(UEqn());
 
     if (simple.momentumPredictor())
     {

applications/solvers/heatTransfer/buoyantSimpleFoam/buoyantSimpleFoam.C

@@ -22,19 +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 "IOMRFZoneList.H"
+#include "IObasicSourceList.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -45,7 +46,7 @@ int main(int argc, char *argv[])
     #include "createMesh.H"
     #include "readGravitationalAcceleration.H"
     #include "createFields.H"
-    #include "createZones.H"
+    #include "createRadiationModel.H"
     #include "initContinuityErrs.H"
 
     simpleControl simple(mesh);

applications/solvers/heatTransfer/buoyantSimpleFoam/buoyantSimpleRadiationFoam/EEqn.H

@@ -1,22 +0,0 @@
-{
-    volScalarField& he = thermo.he();
-
-    fvScalarMatrix EEqn
-    (
-        fvm::div(phi, he)
-      + (
-            he.name() == "e"
-          ? fvc::div(phi, volScalarField("Ekp", 0.5*magSqr(U) + p/rho))
-          : fvc::div(phi, volScalarField("K", 0.5*magSqr(U)))
-        )
-      - fvm::laplacian(turbulence->alphaEff(), he)
-     ==
-        radiation->Sh(thermo)
-    );
-
-    EEqn.relax();
-    EEqn.solve();
-
-    thermo.correct();
-    radiation->correct();
-}

applications/solvers/heatTransfer/buoyantSimpleFoam/buoyantSimpleRadiationFoam/Make/files

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

applications/solvers/heatTransfer/buoyantSimpleFoam/buoyantSimpleRadiationFoam/Make/options

@@ -1,18 +0,0 @@
-EXE_INC = \
-    -I.. \
-    -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)/turbulenceModels/RAS \
-    -I$(LIB_SRC)/finiteVolume/cfdTools \
-    -I$(LIB_SRC)/finiteVolume/lnInclude
-
-EXE_LIBS = \
-    -lfluidThermophysicalModels \
-    -lspecie \
-    -lradiationModels \
-    -lcompressibleTurbulenceModel \
-    -lcompressibleRASModels \
-    -lfiniteVolume \
-    -lmeshTools

applications/solvers/heatTransfer/buoyantSimpleFoam/buoyantSimpleRadiationFoam/buoyantSimpleRadiationFoam.C

@@ -1,86 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-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/>.
-
-Application
-    buoyantSimpleRadiationFoam
-
-Description
-    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 "IOMRFZoneList.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-int main(int argc, char *argv[])
-{
-    #include "setRootCase.H"
-    #include "createTime.H"
-    #include "createMesh.H"
-    #include "readGravitationalAcceleration.H"
-    #include "createFields.H"
-    #include "createZones.H"
-    #include "createRadiationModel.H"
-    #include "initContinuityErrs.H"
-
-    simpleControl simple(mesh);
-
-    // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-    Info<< "\nStarting time loop\n" << endl;
-
-    while (simple.loop())
-    {
-        Info<< "Time = " << runTime.timeName() << nl << endl;
-
-        // Pressure-velocity SIMPLE corrector
-        {
-            #include "UEqn.H"
-            #include "EEqn.H"
-            #include "pEqn.H"
-        }
-
-        turbulence->correct();
-
-        runTime.write();
-
-        Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
-            << "  ClockTime = " << runTime.elapsedClockTime() << " s"
-            << nl << endl;
-    }
-
-    Info<< "End\n" << endl;
-
-    return 0;
-}
-
-
-// ************************************************************************* //

applications/solvers/heatTransfer/buoyantSimpleFoam/createFields.H

@@ -87,3 +87,5 @@
 
     dimensionedScalar initialMass = fvc::domainIntegrate(rho);
     dimensionedScalar totalVolume = sum(mesh.V());
+
+    IObasicSourceList sources(mesh);

applications/solvers/heatTransfer/buoyantSimpleFoam/createZones.H

@@ -1,3 +0,0 @@
-    IOMRFZoneList mrfZones(mesh);
-    mrfZones.correctBoundaryVelocity(U);
-

applications/solvers/heatTransfer/buoyantSimpleFoam/pEqn.H

@@ -17,7 +17,7 @@
         fvc::interpolate(rho)*(fvc::interpolate(HbyA) & mesh.Sf())
     );
 
-    mrfZones.relativeFlux(fvc::interpolate(rho), phiHbyA);
+    sources.relativeFlux(fvc::interpolate(rho), phiHbyA);
 
     bool closedVolume = adjustPhi(phiHbyA, U, p_rgh);
 
@@ -45,6 +45,7 @@
             // calculated from the relaxed pressure
             U = HbyA + rAU*fvc::reconstruct((phig - p_rghEqn.flux())/rhorAUf);
             U.correctBoundaryConditions();
+            sources.correct(U);
         }
     }
 

applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/chtMultiRegionSimpleFoam.C

@@ -36,7 +36,6 @@ Description
 #include "regionProperties.H"
 #include "solidThermo.H"
 #include "radiationModel.H"
-#include "IOporosityModelList.H"
 #include "IObasicSourceList.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/fluid/EEqn.H

@@ -16,6 +16,9 @@
     );
 
     EEqn.relax();
+
+    sources.constrain(EEqn);
+
     EEqn.solve();
 
     thermo.correct();

applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/fluid/UEqn.H

@@ -3,13 +3,13 @@
     (
         fvm::div(phi, U)
       + turb.divDevRhoReff(U)
-      + sources(rho, U)
+     ==
+        sources(rho, U)
     );
 
     UEqn().relax();
 
-    // Add porous resistance
+    sources.constrain(UEqn());
-    pZones.addResistance(UEqn());
 
     solve
     (

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

@@ -18,7 +18,6 @@
     PtrList<dimensionedScalar> rhoMin(fluidRegions.size());
 
     PtrList<IObasicSourceList> heatSources(fluidRegions.size());
-    PtrList<IOporosityModelList> porousZonesFluid(fluidRegions.size());
 
     // Populate fluid field pointer lists
     forAll(fluidRegions, i)
@@ -200,13 +199,6 @@
             i,
             new IObasicSourceList(fluidRegions[i])
         );
-
-        Info<< "    Adding porous zones\n" << endl;
-        porousZonesFluid.set
-        (
-            i,
-            new IOporosityModelList(fluidRegions[i])
-        );
     }
 
 

applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/fluid/pEqn.H

@@ -23,6 +23,8 @@
 
     phiHbyA += phig;
 
+    sources.relativeFlux(fvc::interpolate(rho), phiHbyA);
+
     dimensionedScalar compressibility = fvc::domainIntegrate(psi);
     bool compressible = (compressibility.value() > SMALL);
 
@@ -54,6 +56,7 @@
             // calculated from the relaxed pressure
             U = HbyA + rAU*fvc::reconstruct((phig - p_rghEqn.flux())/rhorAUf);
             U.correctBoundaryConditions();
+            sources.correct(U);
         }
     }
 

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

@@ -14,8 +14,6 @@
 
     IObasicSourceList& sources = heatSources[i];
 
-    const IOporosityModelList& pZones = porousZonesFluid[i];
-
     const dimensionedScalar initialMass
     (
         "initialMass",

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

@@ -24,6 +24,8 @@
       + phig
     );
 
+    sources.relativeFlux(fvc::interpolate(rho), phiHbyA);
+
     {
         fvScalarMatrix p_rghDDtEqn
         (
@@ -64,6 +66,7 @@
                 U = HbyA
                   + rAU*fvc::reconstruct((phig + p_rghEqn.flux())/rhorAUf);
                 U.correctBoundaryConditions();
+                sources.correct(U);
                 K = 0.5*magSqr(U);
             }
         }

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

@@ -10,9 +10,14 @@ if (finalIter)
         (
             fvm::ddt(betav*rho, h)
           - fvm::laplacian(betav*alpha, h, "laplacian(alpha,h)")
-          + sources(rho, h)
+          ==
+            sources(rho, h)
         );
+
         hEqn().relax();
+
+        sources.constrain(hEqn());
+
         hEqn().solve(mesh.solver(h.select(finalIter)));
     }
 }

applications/solvers/incompressible/simpleFoam/Allwmake

@@ -5,7 +5,6 @@ set -x
 wmake
 
 wmake SRFSimpleFoam
-wmake MRFSimpleFoam
 wmake porousSimpleFoam
 
 # ----------------------------------------------------------------- end-of-file

applications/solvers/lagrangian/coalChemistryFoam/UEqn.H

@@ -7,6 +7,7 @@
         rho.dimensionedInternalField()*g
       + coalParcels.SU(U)
       + limestoneParcels.SU(U)
+      + sources(rho, U)
     );
 
     UEqn.relax();
@@ -15,6 +16,6 @@
 
     if (pimple.momentumPredictor())
     {
-        solve(UEqn == -fvc::grad(p) + sources(rho, U));
+        solve(UEqn == -fvc::grad(p));
         K = 0.5*magSqr(U);
     }

applications/solvers/lagrangian/coalChemistryFoam/pEqn.H

@@ -2,7 +2,7 @@ rho = thermo.rho();
 
 volScalarField rAU(1.0/UEqn.A());
 volVectorField HbyA("HbyA", U);
-HbyA = rAU*(UEqn == sources(rho, U))().H();
+HbyA = rAU*UEqn.H();
 
 if (pimple.transonic())
 {
@@ -16,6 +16,8 @@ if (pimple.transonic())
         )
     );
 
+    sources.relativeFlux(fvc::interpolate(psi), phid);
+
     while (pimple.correctNonOrthogonal())
     {
         fvScalarMatrix pEqn
@@ -28,7 +30,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())));
 
@@ -50,6 +52,8 @@ else
         )
     );
 
+    sources.relativeFlux(fvc::interpolate(rho), phiHbyA);
+
     while (pimple.correctNonOrthogonal())
     {
         fvScalarMatrix pEqn
@@ -62,7 +66,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/lagrangian/reactingParcelFilmFoam/EEqn.H

@@ -21,9 +21,11 @@
       + surfaceFilm.Sh()
       + radiation->Sh(thermo)
       + combustion->Sh()
+      + sources(rho, he)
     );
 
     EEqn.relax();
+    sources.constrain(EEqn);
     EEqn.solve();
 
     thermo.correct();

applications/solvers/lagrangian/reactingParcelFilmFoam/Make/options

@@ -1,5 +1,7 @@
 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)/turbulenceModels/compressible/turbulenceModel \
     -I$(LIB_SRC)/lagrangian/distributionModels/lnInclude \
@@ -25,6 +27,8 @@ EXE_INC = \
 
 EXE_LIBS = \
     -lfiniteVolume \
+    -lfieldSources \
+    -lsampling \
     -lmeshTools \
     -lcompressibleRASModels \
     -lcompressibleLESModels \

applications/solvers/lagrangian/reactingParcelFilmFoam/UEqn.H

@@ -5,10 +5,13 @@
       + turbulence->divDevRhoReff(U)
      ==
         parcels.SU(U)
+      + sources(rho, U)
     );
 
     UEqn.relax();
 
+    sources.constrain(UEqn);
+
     if (pimple.momentumPredictor())
     {
         solve

applications/solvers/lagrangian/reactingParcelFilmFoam/YEqn.H

@@ -22,18 +22,24 @@ tmp<fv::convectionScheme<scalar> > mvConvection
         {
             volScalarField& Yi = Y[i];
 
-            solve
+            fvScalarMatrix YEqn
             (
                 fvm::ddt(rho, Yi)
               + mvConvection->fvmDiv(phi, Yi)
               - fvm::laplacian(turbulence->alphaEff(), Yi)
               ==
                 parcels.SYi(i, Yi)
+              + sources(rho, Yi)
+              + combustion->R(Yi)
               + surfaceFilm.Srho(i)
-              + combustion->R(Yi),
-                mesh.solver("Yi")
             );
 
+            YEqn.relax();
+
+            sources.constrain(YEqn);
+
+            YEqn.solve(mesh.solver("Yi"));
+
             Yi.max(0.0);
             Yt += Yi;
         }

applications/solvers/lagrangian/reactingParcelFilmFoam/createFields.H

@@ -139,3 +139,5 @@
         mesh,
         dimensionedScalar("dQ", dimEnergy/dimTime, 0.0)
     );
+
+    IObasicSourceList sources(mesh);

applications/solvers/lagrangian/reactingParcelFilmFoam/pEqn.H

@@ -18,6 +18,7 @@ surfaceScalarField phiHbyA
   + phig
 );
 
+sources.relativeFlux(fvc::interpolate(rho), phiHbyA);
 
 while (pimple.correctNonOrthogonal())
 {
@@ -30,8 +31,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())
@@ -39,6 +43,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/lagrangian/reactingParcelFilmFoam/reactingParcelFilmFoam.C

@@ -37,6 +37,7 @@ Description
 #include "psiCombustionModel.H"
 #include "radiationModel.H"
 #include "SLGThermo.H"
+#include "IObasicSourceList.H"
 #include "pimpleControl.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

applications/solvers/lagrangian/reactingParcelFilmFoam/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
@@ -37,6 +37,7 @@ Description
       ==
         parcels.Srho(rho)
       + surfaceFilm.Srho()
+      + sources(rho)
     );
 }
 

applications/solvers/lagrangian/reactingParcelFoam/UEqn.H

@@ -7,6 +7,7 @@
      ==
         rho.dimensionedInternalField()*g
       + parcels.SU(U)
+      + sources(rho, U)
     );
 
     UEqn.relax();
@@ -17,6 +18,6 @@
 
     if (pimple.momentumPredictor())
     {
-        solve(UEqn == -fvc::grad(p) + sources(rho, U));
+        solve(UEqn == -fvc::grad(p));
         K = 0.5*magSqr(U);
     }

applications/solvers/lagrangian/reactingParcelFoam/pEqn.H

@@ -7,17 +7,19 @@
 
     volScalarField rAU(1.0/UEqn.A());
     volVectorField HbyA("HbyA", U);
-    HbyA = rAU*(UEqn == sources(rho, U))().H();
+    HbyA = rAU*UEqn.H();
 
-    surfaceScalarField phiHbyA("phiHbyA", fvc::interpolate(HbyA) & mesh.Sf());
+    surfaceScalarField phiHbyA
-    if (!pZones.active())
+    (
-    {
+        "phiHbyA",
-        // ddtPhiCorr only used without porosity
+        fvc::interpolate(rho)
-        phiHbyA += fvc::ddtPhiCorr(rAU, rho, U, phi);
+       *(
-    }
+            (fvc::interpolate(HbyA) & mesh.Sf())
-
+          + fvc::ddtPhiCorr(rAU, rho, U, phi)
-    phiHbyA *= fvc::interpolate(rho);
+        )
+    );
 
+    sources.relativeFlux(fvc::interpolate(rho), phiHbyA);
 
     fvScalarMatrix pDDtEqn
     (
@@ -36,7 +38,7 @@
           - fvm::laplacian(rho*rAU, p)
         );
 
-        sources.constrain(pDDtEqn, rho.name());
+        sources.constrain(pEqn);
 
         pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter())));
 

applications/solvers/lagrangian/sprayFoam/UEqn.H

@@ -6,6 +6,7 @@
      ==
         rho.dimensionedInternalField()*g
       + parcels.SU(U)
+      + sources(rho, U)
     );
 
     UEqn.relax();
@@ -14,6 +15,6 @@
 
     if (pimple.momentumPredictor())
     {
-        solve(UEqn == -fvc::grad(p) + sources(rho, U));
+        solve(UEqn == -fvc::grad(p));
         K = 0.5*magSqr(U);
     }

applications/solvers/lagrangian/sprayFoam/pEqn.H

@@ -2,7 +2,7 @@ rho = thermo.rho();
 
 volScalarField rAU(1.0/UEqn.A());
 volVectorField HbyA("HbyA", U);
-HbyA = rAU*(UEqn == sources(rho, U))().H();
+HbyA = rAU*UEqn.H();
 
 if (pimple.transonic())
 {
@@ -16,6 +16,8 @@ if (pimple.transonic())
         )
     );
 
+    sources.relativeFlux(fvc::interpolate(psi), phid);
+
     while (pimple.correctNonOrthogonal())
     {
         fvScalarMatrix pEqn
@@ -28,6 +30,8 @@ if (pimple.transonic())
           + sources(psi, p, rho.name())
         );
 
+        sources.constrain(pEqn);
+
         pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter())));
 
         if (pimple.finalNonOrthogonalIter())
@@ -48,6 +52,8 @@ else
         )
     );
 
+    sources.relativeFlux(fvc::interpolate(psi), phiHbyA);
+
     while (pimple.correctNonOrthogonal())
     {
         fvScalarMatrix pEqn
@@ -60,6 +66,8 @@ else
           + sources(psi, p, rho.name())
         );
 
+        sources.constrain(pEqn);
+
         pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter())));
 
         if (pimple.finalNonOrthogonalIter())

applications/solvers/lagrangian/sprayFoam/sprayEngineFoam/pEqn.H

@@ -2,7 +2,7 @@ rho = thermo.rho();
 
 volScalarField rAU(1.0/UEqn.A());
 volVectorField HbyA("HbyA", U);
-HbyA = rAU*(UEqn == sources(rho, U))().H();
+HbyA = rAU*UEqn.H();
 
 if (pimple.transonic())
 {
@@ -16,6 +16,8 @@ if (pimple.transonic())
         )
     );
 
+    sources.relativeFlux(fvc::interpolate(psi), phid);
+
     while (pimple.correctNonOrthogonal())
     {
         fvScalarMatrix pEqn
@@ -28,6 +30,8 @@ if (pimple.transonic())
           + sources(psi, p, rho.name())
         );
 
+        sources.constrain(pEqn);
+
         pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter())));
 
         if (pimple.finalNonOrthogonalIter())
@@ -48,6 +52,8 @@ else
         )
     );
 
+    sources.relativeFlux(fvc::interpolate(rho), phiHbyA);
+
     while (pimple.correctNonOrthogonal())
     {
         fvScalarMatrix pEqn
@@ -60,6 +66,8 @@ else
           + sources(psi, p, rho.name())
         );
 
+        sources.constrain(pEqn);
+
         pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter())));
 
         if (pimple.finalNonOrthogonalIter())

src/finiteVolume/cfdTools/general/porosityModel/DarcyForchheimer/DarcyForchheimer.C

@@ -56,7 +56,7 @@ Foam::porosityModels::DarcyForchheimer::DarcyForchheimer
     D_("D", dimless/sqr(dimLength), tensor::zero),
     F_("F", dimless/dimLength, tensor::zero),
     rhoName_(coeffs_.lookupOrDefault<word>("rho", "rho")),
-    muName_(coeffs_.lookupOrDefault<word>("mu", "mu")),
+    muName_(coeffs_.lookupOrDefault<word>("mu", "thermo:mu")),
     nuName_(coeffs_.lookupOrDefault<word>("nu", "nu"))
 {
     // local-to-global transformation tensor

src/finiteVolume/cfdTools/general/porosityModel/porosityModel/porosityModel.C

@@ -87,7 +87,7 @@ Foam::porosityModel::porosityModel
 {
     if (zoneName_ == word::null)
     {
-        dict.lookup("actuve") >> active_;
+        dict.lookup("active") >> active_;
         dict_.lookup("cellZone") >> zoneName_;
     }
 

src/finiteVolume/fields/fvPatchFields/derived/supersonicFreestream/supersonicFreestreamFvPatchVectorField.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,9 @@ supersonicFreestreamFvPatchVectorField
 )
 :
     mixedFvPatchVectorField(p, iF),
+    TName_("T"),
+    pName_("p"),
+    psiName_("thermo:psi"),
     UInf_(vector::zero),
     pInf_(0),
     TInf_(0),
@@ -59,6 +62,9 @@ supersonicFreestreamFvPatchVectorField
 )
 :
     mixedFvPatchVectorField(ptf, p, iF, mapper),
+    TName_(ptf.TName_),
+    pName_(ptf.pName_),
+    psiName_(ptf.psiName_),
     UInf_(ptf.UInf_),
     pInf_(ptf.pInf_),
     TInf_(ptf.TInf_),
@@ -75,6 +81,9 @@ supersonicFreestreamFvPatchVectorField
 )
 :
     mixedFvPatchVectorField(p, iF),
+    TName_(dict.lookupOrDefault<word>("T", "T")),
+    pName_(dict.lookupOrDefault<word>("p", "p")),
+    psiName_(dict.lookupOrDefault<word>("psi", "thermo:psi")),
     UInf_(dict.lookup("UInf")),
     pInf_(readScalar(dict.lookup("pInf"))),
     TInf_(readScalar(dict.lookup("TInf"))),
@@ -120,6 +129,9 @@ supersonicFreestreamFvPatchVectorField
 )
 :
     mixedFvPatchVectorField(sfspvf),
+    TName_(sfspvf.TName_),
+    pName_(sfspvf.pName_),
+    psiName_(sfspvf.psiName_),
     UInf_(sfspvf.UInf_),
     pInf_(sfspvf.pInf_),
     TInf_(sfspvf.TInf_),
@@ -135,6 +147,9 @@ supersonicFreestreamFvPatchVectorField
 )
 :
     mixedFvPatchVectorField(sfspvf, iF),
+    TName_(sfspvf.TName_),
+    pName_(sfspvf.pName_),
+    psiName_(sfspvf.psiName_),
     UInf_(sfspvf.UInf_),
     pInf_(sfspvf.pInf_),
     TInf_(sfspvf.TInf_),
@@ -152,13 +167,13 @@ void Foam::supersonicFreestreamFvPatchVectorField::updateCoeffs()
     }
 
     const fvPatchField<scalar>& pT =
-        patch().lookupPatchField<volScalarField, scalar>("T");
+        patch().lookupPatchField<volScalarField, scalar>(TName_);
 
     const fvPatchField<scalar>& pp =
-        patch().lookupPatchField<volScalarField, scalar>("p");
+        patch().lookupPatchField<volScalarField, scalar>(pName_);
 
     const fvPatchField<scalar>& ppsi =
-        patch().lookupPatchField<volScalarField, scalar>("psi");
+        patch().lookupPatchField<volScalarField, scalar>(psiName_);
 
     // Need R of the free-stream flow.  Assume R is independent of location
     // along patch so use face 0
@@ -288,6 +303,9 @@ void Foam::supersonicFreestreamFvPatchVectorField::updateCoeffs()
 void Foam::supersonicFreestreamFvPatchVectorField::write(Ostream& os) const
 {
     fvPatchVectorField::write(os);
+    writeEntryIfDifferent<word>(os, "T", "T", TName_);
+    writeEntryIfDifferent<word>(os, "p", "p", pName_);
+    writeEntryIfDifferent<word>(os, "psi", "thermo:psi", psiName_);
     os.writeKeyword("UInf") << UInf_ << token::END_STATEMENT << nl;
     os.writeKeyword("pInf") << pInf_ << token::END_STATEMENT << nl;
     os.writeKeyword("TInf") << TInf_ << token::END_STATEMENT << nl;

src/finiteVolume/fields/fvPatchFields/derived/supersonicFreestream/supersonicFreestreamFvPatchVectorField.H

@@ -40,6 +40,9 @@ Description
 
     \table
         Property     | Description             | Required    | Default value
+        TName        | Temperature field name  | no          | T
+        pName        | Pressure field name     | no          | p
+        psiName      | Compressibility field name | no       | thermo:psi
         UInf         | free-stream velocity    | yes         |
         pInf         | free-stream pressure    | yes         |
         TInf         | free-stream temperature | yes         |
@@ -88,6 +91,15 @@ class supersonicFreestreamFvPatchVectorField
 {
     // Private data
 
+        //- Name of temperature field, default = "T"
+        word TName_;
+
+        //- Name of pressure field, default = "p"
+        word pName_;
+
+        //- Name of compressibility field field, default = "thermo:psi"
+        word psiName_;
+
         //- Velocity of the free stream
         vector UInf_;
 

src/finiteVolume/fields/fvPatchFields/derived/totalTemperature/totalTemperatureFvPatchScalarField.C

@@ -40,7 +40,7 @@ Foam::totalTemperatureFvPatchScalarField::totalTemperatureFvPatchScalarField
     fixedValueFvPatchScalarField(p, iF),
     UName_("U"),
     phiName_("phi"),
-    psiName_("psi"),
+    psiName_("thermo:psi"),
     gamma_(0.0),
     T0_(p.size(), 0.0)
 {}
@@ -73,7 +73,7 @@ Foam::totalTemperatureFvPatchScalarField::totalTemperatureFvPatchScalarField
     fixedValueFvPatchScalarField(p, iF),
     UName_(dict.lookupOrDefault<word>("U", "U")),
     phiName_(dict.lookupOrDefault<word>("phi", "phi")),
-    psiName_(dict.lookupOrDefault<word>("psi", "psi")),
+    psiName_(dict.lookupOrDefault<word>("psi", "thermo:psi")),
     gamma_(readScalar(dict.lookup("gamma"))),
     T0_("T0", dict, p.size())
 {
@@ -179,7 +179,7 @@ void Foam::totalTemperatureFvPatchScalarField::write(Ostream& os) const
     fvPatchScalarField::write(os);
     writeEntryIfDifferent<word>(os, "U", "U", UName_);
     writeEntryIfDifferent<word>(os, "phi", "phi", phiName_);
-    writeEntryIfDifferent<word>(os, "psi", "psi", psiName_);
+    writeEntryIfDifferent<word>(os, "psi", "thermo:psi", psiName_);
     os.writeKeyword("gamma") << gamma_ << token::END_STATEMENT << nl;
     T0_.writeEntry("T0", os);
     writeEntry("value", os);

src/finiteVolume/fields/fvPatchFields/derived/totalTemperature/totalTemperatureFvPatchScalarField.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
@@ -36,7 +36,7 @@ Description
         Property     | Description             | Required    | Default value
         U            | Velocity field name     | no          | U
         phi          | Flux field name         | no          | phi
-        psi          | Compressibility field name | no       | psi
+        psi          | Compressibility field name | no       | thermo:psi
         gamma        | ratio of specific heats (Cp/Cv) | yes |
         T0           | reference temperature   | yes         |
     \endtable

src/finiteVolume/fields/fvPatchFields/derived/waveTransmissive/waveTransmissiveFvPatchField.C

@@ -41,7 +41,7 @@ Foam::waveTransmissiveFvPatchField<Type>::waveTransmissiveFvPatchField
 )
 :
     advectiveFvPatchField<Type>(p, iF),
-    psiName_("psi"),
+    psiName_("thermo:psi"),
     gamma_(0.0)
 {}
 
@@ -70,7 +70,7 @@ Foam::waveTransmissiveFvPatchField<Type>::waveTransmissiveFvPatchField
 )
 :
     advectiveFvPatchField<Type>(p, iF, dict),
-    psiName_(dict.lookupOrDefault<word>("psi", "psi")),
+    psiName_(dict.lookupOrDefault<word>("psi", "thermo:psi")),
     gamma_(readScalar(dict.lookup("gamma")))
 {}
 
@@ -108,27 +108,21 @@ Foam::waveTransmissiveFvPatchField<Type>::advectionSpeed() const
 {
     // Lookup the velocity and compressibility of the patch
     const fvPatchField<scalar>& psip =
-        this->patch().template lookupPatchField<volScalarField, scalar>
+        this->patch().template
-        (
+            lookupPatchField<volScalarField, scalar>(psiName_);
-            psiName_
-        );
 
     const surfaceScalarField& phi =
         this->db().template lookupObject<surfaceScalarField>(this->phiName_);
 
     fvsPatchField<scalar> phip =
-        this->patch().template lookupPatchField<surfaceScalarField, scalar>
+        this->patch().template
-        (
+            lookupPatchField<surfaceScalarField, scalar>(this->phiName_);
-            this->phiName_
-        );
 
     if (phi.dimensions() == dimDensity*dimVelocity*dimArea)
     {
         const fvPatchScalarField& rhop =
-            this->patch().template lookupPatchField<volScalarField, scalar>
+            this->patch().template
-            (
+                lookupPatchField<volScalarField, scalar>(this->rhoName_);
-                this->rhoName_
-            );
 
         phip /= rhop;
     }
@@ -145,18 +139,12 @@ void Foam::waveTransmissiveFvPatchField<Type>::write(Ostream& os) const
 {
     fvPatchField<Type>::write(os);
 
-    if (this->phiName_ != "phi")
+    this->template
-    {
+        writeEntryIfDifferent<word>(os, "phi", "phi", this->phiName_);
-        os.writeKeyword("phi") << this->phiName_ << token::END_STATEMENT << nl;
+    this->template
-    }
+        writeEntryIfDifferent<word>(os, "rho", "rho", this->rhoName_);
-    if (this->rhoName_ != "rho")
+    this->template
-    {
+        writeEntryIfDifferent<word>(os, "psi", "thermo:psi", psiName_);
-        os.writeKeyword("rho") << this->rhoName_ << token::END_STATEMENT << nl;
-    }
-    if (psiName_ != "psi")
-    {
-        os.writeKeyword("psi") << psiName_ << token::END_STATEMENT << nl;
-    }
 
     os.writeKeyword("gamma") << gamma_ << token::END_STATEMENT << nl;
 

src/regionModels/surfaceFilmModels/kinematicSingleLayer/kinematicSingleLayer.C

@@ -765,7 +765,7 @@ kinematicSingleLayer::kinematicSingleLayer
     (
         IOobject
         (
-            "mu", // must have same name as mu to enable mapping
+            "thermo:mu", // must have same name as mu to enable mapping
             time().timeName(),
             regionMesh(),
             IOobject::NO_READ,

src/thermophysicalModels/radiationModels/radiationModel/noRadiation/noRadiation.C

@@ -52,7 +52,7 @@ Foam::radiation::noRadiation::noRadiation
     const volScalarField& T
 )
 :
-    radiationModel(dict, T)
+    radiationModel(T)
 {}
 
 

src/thermophysicalModels/radiationModels/radiationModel/radiationModel/radiationModel.C

@@ -43,6 +43,33 @@ namespace Foam
 
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
+Foam::IOobject Foam::radiation::radiationModel::createIOobject
+(
+    const fvMesh& mesh
+) const
+{
+    IOobject io
+    (
+        "radiationProperties",
+        mesh.time().constant(),
+        mesh,
+        IOobject::MUST_READ,
+        IOobject::NO_WRITE
+    );
+
+    if (io.headerOk())
+    {
+        io.readOpt() = IOobject::MUST_READ_IF_MODIFIED;
+        return io;
+    }
+    else
+    {
+        io.readOpt() = IOobject::NO_READ;
+        return io;
+    }
+}
+
+
 void Foam::radiation::radiationModel::initialise()
 {
     if (radiation_)
@@ -62,35 +89,6 @@ void Foam::radiation::radiationModel::initialise()
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 Foam::radiation::radiationModel::radiationModel(const volScalarField& T)
-:
-    IOdictionary
-    (
-        IOobject
-        (
-            "radiationProperties",
-            T.time().constant(),
-            T.mesh(),
-            IOobject::MUST_READ_IF_MODIFIED,
-            IOobject::NO_WRITE
-        )
-    ),
-    mesh_(T.mesh()),
-    time_(T.time()),
-    T_(T),
-    radiation_(false),
-    coeffs_(dictionary::null),
-    solverFreq_(0),
-    firstIter_(true),
-    absorptionEmission_(NULL),
-    scatter_(NULL)
-{}
-
-
-Foam::radiation::radiationModel::radiationModel
-(
-    const dictionary& dict,
-    const volScalarField& T
-)
 :
     IOdictionary
     (
@@ -101,8 +99,7 @@ Foam::radiation::radiationModel::radiationModel
             T.mesh(),
             IOobject::NO_READ,
             IOobject::NO_WRITE
-        ),
+        )
-        dict
     ),
     mesh_(T.mesh()),
     time_(T.time()),
@@ -122,17 +119,7 @@ Foam::radiation::radiationModel::radiationModel
     const volScalarField& T
 )
 :
-    IOdictionary
+    IOdictionary(createIOobject(T.mesh())),
-    (
-        IOobject
-        (
-            "radiationProperties",
-            T.time().constant(),
-            T.mesh(),
-            IOobject::MUST_READ_IF_MODIFIED,
-            IOobject::NO_WRITE
-        )
-    ),
     mesh_(T.mesh()),
     time_(T.time()),
     T_(T),
@@ -143,6 +130,11 @@ Foam::radiation::radiationModel::radiationModel
     absorptionEmission_(NULL),
     scatter_(NULL)
 {
+    if (readOpt() == IOobject::NO_READ)
+    {
+        radiation_ = false;
+    }
+
     initialise();
 }
 

src/thermophysicalModels/radiationModels/radiationModel/radiationModel/radiationModel.H

@@ -109,6 +109,9 @@ private:
 
     // Private Member Functions
 
+        //- Create IO object if dictionary is present
+        IOobject createIOobject(const fvMesh& mesh) const;
+
         //- Initialise
         void initialise();
 
@@ -156,9 +159,6 @@ public:
         //- Null constructor
         radiationModel(const volScalarField& T);
 
-        //- Construct with dictionary
-        radiationModel(const dictionary& dict, const volScalarField& T);
-
         //- Construct from components
         radiationModel(const word& type, const volScalarField& T);
 

src/thermophysicalModels/radiationModels/radiationModel/radiationModel/radiationModelNew.C

@@ -33,12 +33,7 @@ Foam::radiation::radiationModel::New
     const volScalarField& T
 )
 {
-    // get model name, but do not register the dictionary
+    IOobject radIO
-    const word modelType
-    (
-        IOdictionary
-        (
-            IOobject
     (
         "radiationProperties",
         T.time().constant(),
@@ -46,10 +41,19 @@ Foam::radiation::radiationModel::New
         IOobject::MUST_READ_IF_MODIFIED,
         IOobject::NO_WRITE,
         false
-            )
-        ).lookup("radiationModel")
     );
 
+    word modelType("none");
+    if (radIO.headerOk())
+    {
+        IOdictionary(radIO).lookup("radiationModel") >> modelType;
+    }
+    else
+    {
+        Info<< "Radiation model not active: radiationProperties not found"
+            << endl;
+    }
+    
     Info<< "Selecting radiationModel " << modelType << endl;
 
     TConstructorTable::iterator cstrIter =

tutorials/combustion/fireFoam/les/oppositeBurningPanels/system/fvSchemes

@@ -50,7 +50,7 @@ laplacianSchemes
     default         Gauss linear corrected;
     laplacian(muEff,U) Gauss linear corrected;
     laplacian(DkEff,k) Gauss linear corrected;
-    laplacian(alphaEff,h) Gauss linear corrected;
+    laplacian(thermo:alpha,h) Gauss linear corrected;
     laplacian((((rho*(1|A(U)))*rho)*gh)) Gauss linear corrected;
     laplacian(interpolate((rho*(1|A(U)))),p) Gauss linear corrected;
     laplacian(gammaRad,G) Gauss linear corrected;

tutorials/combustion/fireFoam/les/oppositeBurningPanels/system/panelRegion/fvSchemes

@@ -32,7 +32,7 @@ divSchemes
 laplacianSchemes
 {
     default         none;
-    laplacian(alpha,h)  Gauss linear uncorrected;
+    laplacian(thermo:alpha,h) Gauss linear uncorrected;
 }
 
 interpolationSchemes

tutorials/compressible/rhoCentralFoam/LadenburgJet60psi/0.org/T

@@ -39,7 +39,7 @@ boundaryField
         U               U;
         phi             phi;
         rho             none;
-        psi             psi;
+        psi             thermo:psi;
         gamma           1.4;
     }
 

tutorials/compressible/rhoCentralFoam/LadenburgJet60psi/0/T

@@ -1243,7 +1243,7 @@ boundaryField
         type            totalTemperature;
         U               U;
         phi             phi;
-        psi             psi;
+        psi             thermo:psi;
         gamma           1.4;
         T0              uniform 297;
         value           nonuniform List<scalar>

tutorials/compressible/rhoCentralFoam/LadenburgJet60psi/0/p

@@ -1235,7 +1235,7 @@ boundaryField
         field           p;
         phi             phi;
         rho             rho;
-        psi             psi;
+        psi             thermo:psi;
         gamma           1.4;
         fieldInf        101325;
         lInf            0.025;
@@ -1275,7 +1275,7 @@ boundaryField
         U               U;
         phi             phi;
         rho             none;
-        psi             psi;
+        psi             thermo:psi;
         gamma           1.4;
         p0              uniform 101325;
         value           nonuniform List<scalar>