diff --git a/applications/solvers/compressible/rhoPimpleFoam/EEqn.H b/applications/solvers/compressible/rhoPimpleFoam/EEqn.H
index 6288b8a0af0842842ab061aa1148c804bfb9b7a4..1cc6ed584c6c3806864cc0a0f461182796b9799a 100644
--- a/applications/solvers/compressible/rhoPimpleFoam/EEqn.H
+++ b/applications/solvers/compressible/rhoPimpleFoam/EEqn.H
@@ -16,9 +16,12 @@
           : -dpdt
         )
       - fvm::laplacian(turbulence->alphaEff(), he)
+     ==
+        sources(rho, he)
     );
 
     EEqn.relax();
+    sources.constrain(EEqn);
     EEqn.solve();
 
     thermo.correct();
diff --git a/applications/solvers/compressible/rhoPimpleFoam/Make/options b/applications/solvers/compressible/rhoPimpleFoam/Make/options
index f21e7394e4db3c8d5cf3770ecebbe46e30e4ec50..c96a146eca26d9bfdb69c0526632707b4cea13e4 100644
--- a/applications/solvers/compressible/rhoPimpleFoam/Make/options
+++ b/applications/solvers/compressible/rhoPimpleFoam/Make/options
@@ -2,7 +2,10 @@ EXE_INC = \
     -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)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/meshTools/lnInclude \
+    -I$(LIB_SRC)/sampling/lnInclude \
+    -I$(LIB_SRC)/fieldSources/lnInclude
 
 EXE_LIBS = \
     -lfluidThermophysicalModels \
@@ -11,4 +14,6 @@ EXE_LIBS = \
     -lcompressibleRASModels \
     -lcompressibleLESModels \
     -lfiniteVolume \
-    -lmeshTools
+    -lmeshTools \
+    -lsampling \
+    -lfieldSources
diff --git a/applications/solvers/compressible/rhoPimpleFoam/UEqn.H b/applications/solvers/compressible/rhoPimpleFoam/UEqn.H
index 0f5ec2487b30854d2982e39196ea194bd9d32c4e..397e8930357dfd548a96f62b9f63acfc943f2423 100644
--- a/applications/solvers/compressible/rhoPimpleFoam/UEqn.H
+++ b/applications/solvers/compressible/rhoPimpleFoam/UEqn.H
@@ -9,8 +9,10 @@ tmp<fvVectorMatrix> UEqn
 
 UEqn().relax();
 
+sources.constrain(UEqn());
+
 if (pimple.momentumPredictor())
 {
-    solve(UEqn() == -fvc::grad(p));
+    solve(UEqn() == -fvc::grad(p) + sources(rho, U));
     K = 0.5*magSqr(U);
 }
diff --git a/applications/solvers/compressible/rhoPimpleFoam/createFields.H b/applications/solvers/compressible/rhoPimpleFoam/createFields.H
index 67cc0c3e456573e71fe9f4937a5fd5d3e9253692..8b3cfd10b9a1180d574c1610ccec487b3ab5ae68 100644
--- a/applications/solvers/compressible/rhoPimpleFoam/createFields.H
+++ b/applications/solvers/compressible/rhoPimpleFoam/createFields.H
@@ -69,3 +69,6 @@
 
     Info<< "Creating field kinetic energy K\n" << endl;
     volScalarField K("K", 0.5*magSqr(U));
+
+    Info<< "Creating sources\n" << endl;
+    IObasicSourceList sources(mesh);
diff --git a/applications/solvers/compressible/rhoPimpleFoam/pEqn.H b/applications/solvers/compressible/rhoPimpleFoam/pEqn.H
index d40c66d40252239b172f7d493da57b51c92e7d93..5f28a3456d82c8d6c1b8f05d90b1519e9a5f5931 100644
--- a/applications/solvers/compressible/rhoPimpleFoam/pEqn.H
+++ b/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().H();
+HbyA = rAU*(UEqn() == sources(rho, U))().H();
 
 if (pimple.nCorrPISO() <= 1)
 {
@@ -33,8 +33,12 @@ if (pimple.transonic())
             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())
@@ -65,8 +69,12 @@ else
             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())
@@ -92,6 +100,7 @@ Info<< "rho max/min : " << max(rho).value()
 
 U = HbyA - rAU*fvc::grad(p);
 U.correctBoundaryConditions();
+sources.correct(U);
 K = 0.5*magSqr(U);
 
 if (thermo.dpdt())
diff --git a/applications/solvers/compressible/rhoPimpleFoam/rhoPimpleFoam.C b/applications/solvers/compressible/rhoPimpleFoam/rhoPimpleFoam.C
index adb6882268a7315be173b4b8b5c175937fc7ada4..f9b8a901f721265fe10cd8c478651f8334ba8497 100644
--- a/applications/solvers/compressible/rhoPimpleFoam/rhoPimpleFoam.C
+++ b/applications/solvers/compressible/rhoPimpleFoam/rhoPimpleFoam.C
@@ -38,6 +38,7 @@ Description
 #include "turbulenceModel.H"
 #include "bound.H"
 #include "pimpleControl.H"
+#include "IObasicSourceList.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
diff --git a/applications/solvers/compressible/rhoPimpleFoam/rhoPimplecFoam/Make/options b/applications/solvers/compressible/rhoPimpleFoam/rhoPimplecFoam/Make/options
index 6d3c41f0dfa6bde2b2e91e36a9458c172ca6cf13..c78d474b4fb3d4649d7afad5957132ee2da37253 100644
--- a/applications/solvers/compressible/rhoPimpleFoam/rhoPimplecFoam/Make/options
+++ b/applications/solvers/compressible/rhoPimpleFoam/rhoPimplecFoam/Make/options
@@ -3,7 +3,10 @@ EXE_INC = \
     -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)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/meshTools/lnInclude \
+    -I$(LIB_SRC)/sampling/lnInclude \
+    -I$(LIB_SRC)/fieldSources/lnInclude
 
 EXE_LIBS = \
     -lfluidThermophysicalModels \
@@ -12,4 +15,6 @@ EXE_LIBS = \
     -lcompressibleRASModels \
     -lcompressibleLESModels \
     -lfiniteVolume \
-    -lmeshTools
+    -lmeshTools \
+    -lsampling \
+    -lfieldSources
diff --git a/applications/solvers/compressible/rhoPimpleFoam/rhoPimplecFoam/pEqn.H b/applications/solvers/compressible/rhoPimpleFoam/rhoPimplecFoam/pEqn.H
index 75ed20e271981cb2d350171947c2bc680593b389..628eb71f2961392d5c1173fbbf30aeb07cd5f81b 100644
--- a/applications/solvers/compressible/rhoPimpleFoam/rhoPimplecFoam/pEqn.H
+++ b/applications/solvers/compressible/rhoPimpleFoam/rhoPimplecFoam/pEqn.H
@@ -7,7 +7,7 @@ volScalarField rAU(1.0/UEqn().A());
 volScalarField rAtU(1.0/(1.0/rAU - UEqn().H1()));
 
 volVectorField HbyA("HbyA", U);
-HbyA = rAU*UEqn().H();
+HbyA = rAU*(UEqn() == sources(rho, U))().H();
 
 if (pimple.nCorrPIMPLE() <= 1)
 {
@@ -44,11 +44,15 @@ if (pimple.transonic())
           + fvm::div(phid, p)
           + fvc::div(phic)
           - fvm::laplacian(Dp, p)
+         ==
+            sources(psi, p, rho.name())
         );
 
         // Relax the pressure equation to maintain diagonal dominance
         pEqn.relax();
 
+        sources.constrain(pEqn, rho.name());
+
         pEqn.solve();
 
         if (pimple.finalNonOrthogonalIter())
@@ -81,8 +85,12 @@ else
             fvm::ddt(psi, p)
           + fvc::div(phiHbyA)
           - fvm::laplacian(Dp, p)
+         ==
+            sources(psi, p, rho.name())
         );
 
+        sources.constrain(pEqn, rho.name());
+
         pEqn.solve();
 
         if (pimple.finalNonOrthogonalIter())
@@ -100,6 +108,7 @@ p.relax();
 
 U = HbyA - rAtU*fvc::grad(p);
 U.correctBoundaryConditions();
+sources.correct(U);
 K = 0.5*magSqr(U);
 
 if (thermo.dpdt())
diff --git a/applications/solvers/compressible/rhoPimpleFoam/rhoPimplecFoam/rhoPimplecFoam.C b/applications/solvers/compressible/rhoPimpleFoam/rhoPimplecFoam/rhoPimplecFoam.C
index cf8aed0bcad3870312599f76237f1ae6cd6f8be2..7e1664dbeac745cf8a24e634cfd72b87e4f41213 100644
--- a/applications/solvers/compressible/rhoPimpleFoam/rhoPimplecFoam/rhoPimplecFoam.C
+++ b/applications/solvers/compressible/rhoPimpleFoam/rhoPimplecFoam/rhoPimplecFoam.C
@@ -38,6 +38,7 @@ Description
 #include "turbulenceModel.H"
 #include "bound.H"
 #include "pimpleControl.H"
+#include "IObasicSourceList.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
diff --git a/applications/solvers/compressible/rhoPimpleFoam/rhoPorousMRFLTSPimpleFoam/Make/options b/applications/solvers/compressible/rhoPimpleFoam/rhoPorousMRFLTSPimpleFoam/Make/options
index 5c83c910d95a4aa6b5de30075e4f764309a4bec3..2976366c08ac366a8c0abe6c3291714ed86f1908 100644
--- a/applications/solvers/compressible/rhoPimpleFoam/rhoPorousMRFLTSPimpleFoam/Make/options
+++ b/applications/solvers/compressible/rhoPimpleFoam/rhoPorousMRFLTSPimpleFoam/Make/options
@@ -5,7 +5,9 @@ EXE_INC = \
     -I$(LIB_SRC)/turbulenceModels/compressible/turbulenceModel \
     -I$(LIB_SRC)/finiteVolume/cfdTools \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
-    -I$(LIB_SRC)/meshTools/lnInclude
+    -I$(LIB_SRC)/meshTools/lnInclude \
+    -I$(LIB_SRC)/sampling/lnInclude \
+    -I$(LIB_SRC)/fieldSources/lnInclude
 
 EXE_LIBS = \
     -lfluidThermophysicalModels \
@@ -14,4 +16,6 @@ EXE_LIBS = \
     -lcompressibleRASModels \
     -lcompressibleLESModels \
     -lfiniteVolume \
-    -lmeshTools
+    -lmeshTools \
+    -lsampling \
+    -lfieldSources
diff --git a/applications/solvers/compressible/rhoPimpleFoam/rhoPorousMRFLTSPimpleFoam/rhoPorousMRFLTSPimpleFoam.C b/applications/solvers/compressible/rhoPimpleFoam/rhoPorousMRFLTSPimpleFoam/rhoPorousMRFLTSPimpleFoam.C
index d729b46c37cd5faa6b55aab813845ed397b69b75..82ba5cbac17371b9147b4057d5f657426cff1afc 100644
--- a/applications/solvers/compressible/rhoPimpleFoam/rhoPorousMRFLTSPimpleFoam/rhoPorousMRFLTSPimpleFoam.C
+++ b/applications/solvers/compressible/rhoPimpleFoam/rhoPorousMRFLTSPimpleFoam/rhoPorousMRFLTSPimpleFoam.C
@@ -39,6 +39,7 @@ Description
 #include "turbulenceModel.H"
 #include "MRFZones.H"
 #include "IOporosityModelList.H"
+#include "IObasicSourceList.H"
 #include "fvcSmooth.H"
 #include "pimpleControl.H"
 #include "bound.H"
diff --git a/applications/solvers/compressible/rhoPimpleFoam/rhoPorousMRFPimpleFoam/Make/options b/applications/solvers/compressible/rhoPimpleFoam/rhoPorousMRFPimpleFoam/Make/options
index f01ebda533b3a30f359683e4b72fde503798f830..c78d474b4fb3d4649d7afad5957132ee2da37253 100644
--- a/applications/solvers/compressible/rhoPimpleFoam/rhoPorousMRFPimpleFoam/Make/options
+++ b/applications/solvers/compressible/rhoPimpleFoam/rhoPorousMRFPimpleFoam/Make/options
@@ -4,7 +4,9 @@ EXE_INC = \
     -I$(LIB_SRC)/turbulenceModels/compressible/turbulenceModel \
     -I$(LIB_SRC)/finiteVolume/cfdTools \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
-    -I$(LIB_SRC)/meshTools/lnInclude
+    -I$(LIB_SRC)/meshTools/lnInclude \
+    -I$(LIB_SRC)/sampling/lnInclude \
+    -I$(LIB_SRC)/fieldSources/lnInclude
 
 EXE_LIBS = \
     -lfluidThermophysicalModels \
@@ -13,4 +15,6 @@ EXE_LIBS = \
     -lcompressibleRASModels \
     -lcompressibleLESModels \
     -lfiniteVolume \
-    -lmeshTools
+    -lmeshTools \
+    -lsampling \
+    -lfieldSources
diff --git a/applications/solvers/compressible/rhoPimpleFoam/rhoPorousMRFPimpleFoam/UEqn.H b/applications/solvers/compressible/rhoPimpleFoam/rhoPorousMRFPimpleFoam/UEqn.H
index c24b1f587a59bac6cd7b305b5dc6ccdf55986613..9651610157c483596253639cb1380be1cc8e2f9b 100644
--- a/applications/solvers/compressible/rhoPimpleFoam/rhoPorousMRFPimpleFoam/UEqn.H
+++ b/applications/solvers/compressible/rhoPimpleFoam/rhoPorousMRFPimpleFoam/UEqn.H
@@ -11,9 +11,12 @@ tmp<fvVectorMatrix> UEqn
 UEqn().relax();
 
 mrfZones.addCoriolis(rho, UEqn());
+
 pZones.addResistance(UEqn());
 
+sources.constrain(UEqn());
+
 if (pimple.momentumPredictor())
 {
-    solve(UEqn() == -fvc::grad(p));
+    solve(UEqn() == -fvc::grad(p) + sources(rho, U));
 }
diff --git a/applications/solvers/compressible/rhoPimpleFoam/rhoPorousMRFPimpleFoam/pEqn.H b/applications/solvers/compressible/rhoPimpleFoam/rhoPorousMRFPimpleFoam/pEqn.H
index b68823c6e3d1beaa68743ba7a973e5a17b2d44eb..304798b9662d6c7bb8dd3bfc0bded47e552db19e 100644
--- a/applications/solvers/compressible/rhoPimpleFoam/rhoPorousMRFPimpleFoam/pEqn.H
+++ b/applications/solvers/compressible/rhoPimpleFoam/rhoPorousMRFPimpleFoam/pEqn.H
@@ -5,7 +5,7 @@ rho.relax();
 
 volScalarField rAU(1.0/UEqn().A());
 volVectorField HbyA("HbyA", U);
-HbyA = rAU*UEqn().H();
+HbyA = rAU*(UEqn() == sources(rho, U))().H();
 
 if (pimple.nCorrPISO() <= 1)
 {
@@ -34,8 +34,12 @@ if (pimple.transonic())
             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())
@@ -68,8 +72,12 @@ else
             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())
@@ -94,6 +102,7 @@ Info<< "rho max/min : " << max(rho).value()
 
 U = HbyA - rAU*fvc::grad(p);
 U.correctBoundaryConditions();
+sources.correct(U);
 K = 0.5*magSqr(U);
 
 if (thermo.dpdt())
diff --git a/applications/solvers/compressible/rhoPimpleFoam/rhoPorousMRFPimpleFoam/rhoPorousMRFPimpleFoam.C b/applications/solvers/compressible/rhoPimpleFoam/rhoPorousMRFPimpleFoam/rhoPorousMRFPimpleFoam.C
index 47b3d6c35a6fa76738e74047d29790b7691a1c96..ac5dc1eadcb3d6c9108160928e93fa67d6e4569d 100644
--- a/applications/solvers/compressible/rhoPimpleFoam/rhoPorousMRFPimpleFoam/rhoPorousMRFPimpleFoam.C
+++ b/applications/solvers/compressible/rhoPimpleFoam/rhoPorousMRFPimpleFoam/rhoPorousMRFPimpleFoam.C
@@ -39,6 +39,7 @@ Description
 #include "bound.H"
 #include "MRFZones.H"
 #include "IOporosityModelList.H"
+#include "IObasicSourceList.H"
 #include "pimpleControl.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
diff --git a/applications/solvers/compressible/rhoSimpleFoam/EEqn.H b/applications/solvers/compressible/rhoSimpleFoam/EEqn.H
index e4c79b40cc5685ddbeeab18c32e8611d745656b8..e496906d082b96a7e2dc8a9aa731dba5372608aa 100644
--- a/applications/solvers/compressible/rhoSimpleFoam/EEqn.H
+++ b/applications/solvers/compressible/rhoSimpleFoam/EEqn.H
@@ -10,9 +10,12 @@
           : fvc::div(phi, volScalarField("K", 0.5*magSqr(U)))
         )
       - fvm::laplacian(turbulence->alphaEff(), he)
+     ==
+        sources(rho, he)
     );
 
     EEqn.relax();
+    sources.constrain(EEqn);
     EEqn.solve();
 
     thermo.correct();
diff --git a/applications/solvers/compressible/rhoSimpleFoam/Make/options b/applications/solvers/compressible/rhoSimpleFoam/Make/options
index 339cc53bd9e9e2e045b859cec7a26b8076168565..741dc4f822bc8343080f8afa237be831ba110e00 100644
--- a/applications/solvers/compressible/rhoSimpleFoam/Make/options
+++ b/applications/solvers/compressible/rhoSimpleFoam/Make/options
@@ -3,7 +3,10 @@ EXE_INC = \
     -I$(LIB_SRC)/turbulenceModels \
     -I$(LIB_SRC)/turbulenceModels/compressible/RAS/RASModel \
     -I$(LIB_SRC)/finiteVolume/cfdTools \
-    -I$(LIB_SRC)/finiteVolume/lnInclude
+    -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/sampling/lnInclude \
+    -I$(LIB_SRC)/meshTools/lnInclude \
+    -I$(LIB_SRC)/fieldSources/lnInclude
 
 EXE_LIBS = \
     -lfluidThermophysicalModels \
@@ -11,4 +14,6 @@ EXE_LIBS = \
     -lcompressibleTurbulenceModel \
     -lcompressibleRASModels \
     -lfiniteVolume \
-    -lmeshTools
+    -lsampling \
+    -lmeshTools \
+    -lfieldSources
diff --git a/applications/solvers/compressible/rhoSimpleFoam/createFields.H b/applications/solvers/compressible/rhoSimpleFoam/createFields.H
index dea35b76575477e953c521d77f6bb12cf28c3dc0..3f4b2fd0217cb800ae9c25793630f8266ca4d8d4 100644
--- a/applications/solvers/compressible/rhoSimpleFoam/createFields.H
+++ b/applications/solvers/compressible/rhoSimpleFoam/createFields.H
@@ -60,3 +60,6 @@
     );
 
     dimensionedScalar initialMass = fvc::domainIntegrate(rho);
+    
+    Info<< "Creating sources\n" << endl;
+    IObasicSourceList sources(mesh);
diff --git a/applications/solvers/compressible/rhoSimpleFoam/pEqn.H b/applications/solvers/compressible/rhoSimpleFoam/pEqn.H
index 4d7334c5f8640b44f31a06802f845a7d7c730407..2ea228f172d80bc4d4a330662723e9754840c0a2 100644
--- a/applications/solvers/compressible/rhoSimpleFoam/pEqn.H
+++ b/applications/solvers/compressible/rhoSimpleFoam/pEqn.H
@@ -6,7 +6,7 @@
 
     volScalarField rAU(1.0/UEqn().A());
     volVectorField HbyA("HbyA", U);
-    HbyA = rAU*UEqn().H();
+    HbyA = rAU*(UEqn() == sources(rho, U))().H();
 
     UEqn.clear();
 
@@ -25,12 +25,16 @@
             fvScalarMatrix pEqn
             (
                 fvm::div(phid, p)
-                - fvm::laplacian(rho*rAU, p)
+              - fvm::laplacian(rho*rAU, p)
+              ==
+                sources(psi, p, rho.name())
             );
 
             // Relax the pressure equation to ensure diagonal-dominance
             pEqn.relax();
 
+            sources.constrain(pEqn, rho.name());
+
             pEqn.setReference(pRefCell, pRefValue);
 
             pEqn.solve();
@@ -56,11 +60,15 @@
             fvScalarMatrix pEqn
             (
                 fvc::div(phiHbyA)
-                - fvm::laplacian(rho*rAU, p)
+              - fvm::laplacian(rho*rAU, p)
+              ==
+                sources(psi, p, rho.name())
             );
 
             pEqn.setReference(pRefCell, pRefValue);
 
+            sources.constrain(pEqn, rho.name());
+
             pEqn.solve();
 
             if (simple.finalNonOrthogonalIter())
@@ -78,6 +86,7 @@
 
     U = HbyA - rAU*fvc::grad(p);
     U.correctBoundaryConditions();
+    sources.correct(U);
 
     // For closed-volume cases adjust the pressure and density levels
     // to obey overall mass continuity
diff --git a/applications/solvers/compressible/rhoSimpleFoam/rhoPorousMRFSimpleFoam/EEqn.H b/applications/solvers/compressible/rhoSimpleFoam/rhoPorousMRFSimpleFoam/EEqn.H
index ff467c0382cfcb444f66d64cbbcb549c0dac1736..e496906d082b96a7e2dc8a9aa731dba5372608aa 100644
--- a/applications/solvers/compressible/rhoSimpleFoam/rhoPorousMRFSimpleFoam/EEqn.H
+++ b/applications/solvers/compressible/rhoSimpleFoam/rhoPorousMRFSimpleFoam/EEqn.H
@@ -10,11 +10,12 @@
           : fvc::div(phi, volScalarField("K", 0.5*magSqr(U)))
         )
       - fvm::laplacian(turbulence->alphaEff(), he)
+     ==
+        sources(rho, he)
     );
 
-    pZones.addEnergySource(thermo, rho, EEqn);
-
     EEqn.relax();
+    sources.constrain(EEqn);
     EEqn.solve();
 
     thermo.correct();
diff --git a/applications/solvers/compressible/rhoSimpleFoam/rhoPorousMRFSimpleFoam/Make/options b/applications/solvers/compressible/rhoSimpleFoam/rhoPorousMRFSimpleFoam/Make/options
index 2d328b3ed1f86ecd981c6c3fa3bc0225bd67aa4a..4ba7ed1f7d10bc843dcc0c9d55cd310b08bf98b0 100644
--- a/applications/solvers/compressible/rhoSimpleFoam/rhoPorousMRFSimpleFoam/Make/options
+++ b/applications/solvers/compressible/rhoSimpleFoam/rhoPorousMRFSimpleFoam/Make/options
@@ -1,18 +1,20 @@
 EXE_INC = \
     -I.. \
     -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
-    -I$(LIB_SRC)/thermophysicalModels/thermalPorousZone/lnInclude \
     -I$(LIB_SRC)/turbulenceModels \
     -I$(LIB_SRC)/turbulenceModels/compressible/RAS/RASModel \
     -I$(LIB_SRC)/finiteVolume/cfdTools \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
-    -I$(LIB_SRC)/meshTools/lnInclude
+    -I$(LIB_SRC)/sampling/lnInclude \
+    -I$(LIB_SRC)/meshTools/lnInclude \
+    -I$(LIB_SRC)/fieldSources/lnInclude
 
 EXE_LIBS = \
     -lfluidThermophysicalModels \
-    -lthermalPorousZone \
     -lspecie \
     -lcompressibleTurbulenceModel \
     -lcompressibleRASModels \
     -lfiniteVolume \
-    -lmeshTools
+    -lsampling \
+    -lmeshTools \
+    -lfieldSources
diff --git a/applications/solvers/compressible/rhoSimpleFoam/rhoPorousMRFSimpleFoam/UEqn.H b/applications/solvers/compressible/rhoSimpleFoam/rhoPorousMRFSimpleFoam/UEqn.H
index c0e087b536c3ffdffb01214e7204b13abfba8a46..52e15e07b5143f0e078cec89c9d1e57ebf04108a 100644
--- a/applications/solvers/compressible/rhoSimpleFoam/rhoPorousMRFSimpleFoam/UEqn.H
+++ b/applications/solvers/compressible/rhoSimpleFoam/rhoPorousMRFSimpleFoam/UEqn.H
@@ -24,11 +24,13 @@
         trTU = inv(tTU());
         trTU().rename("rAU");
 
+        sources.constrain(UEqn());
+
         volVectorField gradp(fvc::grad(p));
 
         for (int UCorr=0; UCorr<nUCorr; UCorr++)
         {
-            U = trTU() & (UEqn().H() - gradp);
+            U = trTU() & ((UEqn() == sources(rho, U))().H() - gradp);
         }
         U.correctBoundaryConditions();
     }
@@ -36,7 +38,9 @@
     {
         pZones.addResistance(UEqn());
 
-        solve(UEqn() == -fvc::grad(p));
+        sources.constrain(UEqn());
+
+        solve(UEqn() == -fvc::grad(p) + sources(rho, U));
 
         trAU = 1.0/UEqn().A();
         trAU().rename("rAU");
diff --git a/applications/solvers/compressible/rhoSimpleFoam/rhoPorousMRFSimpleFoam/createFields.H b/applications/solvers/compressible/rhoSimpleFoam/rhoPorousMRFSimpleFoam/createFields.H
index 4fff74d224cadb920dec4f6c42b936dff55f0452..ad924283eef67e97e6fcb193e0a12a6ce834ce9d 100644
--- a/applications/solvers/compressible/rhoSimpleFoam/rhoPorousMRFSimpleFoam/createFields.H
+++ b/applications/solvers/compressible/rhoSimpleFoam/rhoPorousMRFSimpleFoam/createFields.H
@@ -59,3 +59,6 @@
     );
 
     dimensionedScalar initialMass = fvc::domainIntegrate(rho);
+
+    Info<< "Creating sources\n" << endl;
+    IObasicSourceList sources(mesh);
diff --git a/applications/solvers/compressible/rhoSimpleFoam/rhoPorousMRFSimpleFoam/createZones.H b/applications/solvers/compressible/rhoSimpleFoam/rhoPorousMRFSimpleFoam/createZones.H
index e3cfd61f43e9875e6d43037f87e619e9131c8167..4eb2cb193cc43ea551ec1a0e03548bd8f498bdf2 100644
--- a/applications/solvers/compressible/rhoSimpleFoam/rhoPorousMRFSimpleFoam/createZones.H
+++ b/applications/solvers/compressible/rhoSimpleFoam/rhoPorousMRFSimpleFoam/createZones.H
@@ -1,7 +1,7 @@
     MRFZones mrfZones(mesh);
     mrfZones.correctBoundaryVelocity(U);
 
-    thermalPorousZones pZones(mesh);
+    IOporosityModelList pZones(mesh);
     Switch pressureImplicitPorosity(false);
 
     // nUCorrectors used for pressureImplicitPorosity
diff --git a/applications/solvers/compressible/rhoSimpleFoam/rhoPorousMRFSimpleFoam/pEqn.H b/applications/solvers/compressible/rhoSimpleFoam/rhoPorousMRFSimpleFoam/pEqn.H
index 81b194823cb3c22012db5117ce7b7f7f5382b604..708449eb3d872a522f3a6c234937cef21aafe782 100644
--- a/applications/solvers/compressible/rhoSimpleFoam/rhoPorousMRFSimpleFoam/pEqn.H
+++ b/applications/solvers/compressible/rhoSimpleFoam/rhoPorousMRFSimpleFoam/pEqn.H
@@ -4,15 +4,17 @@
     rho = min(rho, rhoMax);
     rho.relax();
 
+    const volScalarField& psi = thermo.psi();
+
     volVectorField HbyA("HbyA", U);
 
     if (pressureImplicitPorosity)
     {
-        HbyA = trTU() & UEqn().H();
+        HbyA = trTU() & (UEqn() == sources(rho, U))().H();
     }
     else
     {
-        HbyA = trAU()*UEqn().H();
+        HbyA = trAU()*(UEqn() == sources(rho, U))().H();
     }
 
     UEqn.clear();
@@ -35,15 +37,29 @@
 
         if (pressureImplicitPorosity)
         {
-            tpEqn = (fvm::laplacian(rho*trTU(), p) == fvc::div(phiHbyA));
+            tpEqn =
+            (
+                fvm::laplacian(rho*trTU(), p)
+              + sources(psi, p, rho.name())
+             ==
+                fvc::div(phiHbyA)
+            );
         }
         else
         {
-            tpEqn = (fvm::laplacian(rho*trAU(), p) == fvc::div(phiHbyA));
+            tpEqn =
+            (
+                fvm::laplacian(rho*trAU(), p)
+              + sources(psi, p, rho.name())
+             ==
+                fvc::div(phiHbyA)
+            );
         }
 
         tpEqn().setReference(pRefCell, pRefValue);
 
+        sources.constrain(tpEqn(), rho.name());
+
         tpEqn().solve();
 
         if (simple.finalNonOrthogonalIter())
@@ -67,12 +83,12 @@
     }
 
     U.correctBoundaryConditions();
+    sources.correct(U);
 
     // For closed-volume cases adjust the pressure and density levels
     // to obey overall mass continuity
     if (closedVolume)
     {
-        const volScalarField& psi = thermo.psi();
         p += (initialMass - fvc::domainIntegrate(psi*p))
             /fvc::domainIntegrate(psi);
     }
diff --git a/applications/solvers/compressible/rhoSimpleFoam/rhoPorousMRFSimpleFoam/rhoPorousMRFSimpleFoam.C b/applications/solvers/compressible/rhoSimpleFoam/rhoPorousMRFSimpleFoam/rhoPorousMRFSimpleFoam.C
index 3898295b915061669daa94fa773218a5877ef357..b552cadfcef253000d0e5352d941ae912d65348b 100644
--- a/applications/solvers/compressible/rhoSimpleFoam/rhoPorousMRFSimpleFoam/rhoPorousMRFSimpleFoam.C
+++ b/applications/solvers/compressible/rhoSimpleFoam/rhoPorousMRFSimpleFoam/rhoPorousMRFSimpleFoam.C
@@ -35,7 +35,8 @@ Description
 #include "rhoThermo.H"
 #include "RASModel.H"
 #include "MRFZones.H"
-#include "thermalPorousZones.H"
+#include "IObasicSourceList.H"
+#include "IOporosityModelList.H"
 #include "simpleControl.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
diff --git a/applications/solvers/compressible/rhoSimpleFoam/rhoSimpleFoam.C b/applications/solvers/compressible/rhoSimpleFoam/rhoSimpleFoam.C
index 312196583eaae9a90ce43b22350f2ecb5d4ac383..efc1067db233ca92b1a286fad6f799a6c77ea6c7 100644
--- a/applications/solvers/compressible/rhoSimpleFoam/rhoSimpleFoam.C
+++ b/applications/solvers/compressible/rhoSimpleFoam/rhoSimpleFoam.C
@@ -34,6 +34,7 @@ Description
 #include "psiThermo.H"
 #include "RASModel.H"
 #include "simpleControl.H"
+#include "IObasicSourceList.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
diff --git a/applications/solvers/compressible/rhoSimpleFoam/rhoSimplecFoam/Make/options b/applications/solvers/compressible/rhoSimpleFoam/rhoSimplecFoam/Make/options
index 8a8553525f39f3e90f00669cec5ee8eb4633f149..facff4ac649182a75207b28def3408506e51bc90 100644
--- a/applications/solvers/compressible/rhoSimpleFoam/rhoSimplecFoam/Make/options
+++ b/applications/solvers/compressible/rhoSimpleFoam/rhoSimplecFoam/Make/options
@@ -4,11 +4,16 @@ EXE_INC = \
     -I$(LIB_SRC)/turbulenceModels \
     -I$(LIB_SRC)/turbulenceModels/compressible/RAS/RASModel \
     -I$(LIB_SRC)/finiteVolume/cfdTools \
-    -I$(LIB_SRC)/finiteVolume/lnInclude
+    -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/sampling/lnInclude \
+    -I$(LIB_SRC)/meshTools/lnInclude \
+    -I$(LIB_SRC)/fieldSources/lnInclude
 
 EXE_LIBS = \
     -lfluidThermophysicalModels \
     -lspecie \
     -lcompressibleRASModels \
     -lfiniteVolume \
-    -lmeshTools
+    -lsampling \
+    -lmeshTools \
+    -lfieldSources
diff --git a/applications/solvers/compressible/rhoSimpleFoam/rhoSimplecFoam/pEqn.H b/applications/solvers/compressible/rhoSimpleFoam/rhoSimplecFoam/pEqn.H
index 8c7405c346ffdebce86d052fe60be505cb52d584..1ffacc1afaf9c5a2495e8e321c84f65a7ea592ac 100644
--- a/applications/solvers/compressible/rhoSimpleFoam/rhoSimplecFoam/pEqn.H
+++ b/applications/solvers/compressible/rhoSimpleFoam/rhoSimplecFoam/pEqn.H
@@ -7,7 +7,7 @@ volScalarField rAU(1.0/UEqn().A());
 volScalarField rAtU(1.0/(1.0/rAU - UEqn().H1()));
 
 volVectorField HbyA("HbyA", U);
-HbyA = rAU*UEqn().H();
+HbyA = rAU*(UEqn() == sources(rho, U))().H();
 
 UEqn.clear();
 
@@ -38,11 +38,15 @@ if (simple.transonic())
             fvm::div(phid, p)
           + fvc::div(phic)
           - fvm::laplacian(Dp, p)
+         ==
+            sources(psi, p, rho.name())
         );
 
         // Relax the pressure equation to maintain diagonal dominance
         pEqn.relax();
 
+        sources.constrain(pEqn, rho.name());
+
         pEqn.setReference(pRefCell, pRefValue);
 
         pEqn.solve();
@@ -74,8 +78,12 @@ else
         (
             fvc::div(phiHbyA)
           - fvm::laplacian(Dp, p)
+          ==
+            sources(psi, p, rho.name())
         );
 
+        sources.constrain(pEqn, rho.name());
+
         pEqn.setReference(pRefCell, pRefValue);
 
         pEqn.solve();
@@ -96,6 +104,7 @@ p.relax();
 
 U = HbyA - rAtU*fvc::grad(p);
 U.correctBoundaryConditions();
+sources.correct(U);
 
 // For closed-volume cases adjust the pressure and density levels
 // to obey overall mass continuity
diff --git a/applications/solvers/compressible/rhoSimpleFoam/rhoSimplecFoam/rhoSimplecFoam.C b/applications/solvers/compressible/rhoSimpleFoam/rhoSimplecFoam/rhoSimplecFoam.C
index f707570f60cd0e0a73162bd7ebabb4857b2b325e..41e15dbfa842e5c8b75228e95c37edf5cc3d2234 100644
--- a/applications/solvers/compressible/rhoSimpleFoam/rhoSimplecFoam/rhoSimplecFoam.C
+++ b/applications/solvers/compressible/rhoSimpleFoam/rhoSimplecFoam/rhoSimplecFoam.C
@@ -36,6 +36,7 @@ Description
 #include "mixedFvPatchFields.H"
 #include "bound.H"
 #include "simpleControl.H"
+#include "IObasicSourceList.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //