diff --git a/applications/solvers/combustion/coldEngineFoam/Make/options b/applications/solvers/combustion/coldEngineFoam/Make/options
index a59f3b8053149415dcc2d7972afaf005f0e18971..25217ab23f79c72d9f81c2aba8eb51c0ae138c30 100644
--- a/applications/solvers/combustion/coldEngineFoam/Make/options
+++ b/applications/solvers/combustion/coldEngineFoam/Make/options
@@ -7,7 +7,10 @@ EXE_INC = \
     -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/reactionThermo/lnInclude \
     -I$(LIB_SRC)/turbulenceModels/compressible/turbulenceModel \
-    -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 = \
     -lengine \
@@ -16,4 +19,7 @@ EXE_LIBS = \
     -lcompressibleLESModels \
     -lfluidThermophysicalModels \
     -lspecie \
-    -lfiniteVolume
+    -lfiniteVolume \
+    -lsampling \
+    -lmeshTools \
+    -lfieldSources
diff --git a/applications/solvers/combustion/coldEngineFoam/coldEngineFoam.C b/applications/solvers/combustion/coldEngineFoam/coldEngineFoam.C
index 87474e778ffe25124a56264d8a95bcaae382d6aa..401690c3e49a263443da462ebca49c6a11863329 100644
--- a/applications/solvers/combustion/coldEngineFoam/coldEngineFoam.C
+++ b/applications/solvers/combustion/coldEngineFoam/coldEngineFoam.C
@@ -35,6 +35,7 @@ Description
 #include "psiThermo.H"
 #include "turbulenceModel.H"
 #include "OFstream.H"
+#include "IObasicSourceList.H"
 #include "pimpleControl.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
diff --git a/applications/solvers/combustion/coldEngineFoam/createFields.H b/applications/solvers/combustion/coldEngineFoam/createFields.H
index f78e5bde2969cd6935e3bccec61887da6f9a2905..72e88435aceacbf2d6fa546b068157dc7caf0de3 100644
--- a/applications/solvers/combustion/coldEngineFoam/createFields.H
+++ b/applications/solvers/combustion/coldEngineFoam/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/heatTransfer/buoyantPimpleFoam/Make/options b/applications/solvers/heatTransfer/buoyantPimpleFoam/Make/options
index f6e12a3b7ce6dcd564b0eacedd671e7f2d1b75c1..5d86b094de32587a79d2bade587ef6f47ccbeb90 100644
--- a/applications/solvers/heatTransfer/buoyantPimpleFoam/Make/options
+++ b/applications/solvers/heatTransfer/buoyantPimpleFoam/Make/options
@@ -2,7 +2,10 @@ 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)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/sampling/lnInclude \
+    -I$(LIB_SRC)/meshTools/lnInclude \
+    -I$(LIB_SRC)/fieldSources/lnInclude \
 
 EXE_LIBS = \
     -lmeshTools \
@@ -11,4 +14,8 @@ EXE_LIBS = \
     -lcompressibleTurbulenceModel \
     -lcompressibleRASModels \
     -lcompressibleLESModels \
-    -lfiniteVolume
+    -lfiniteVolume \
+    -lsampling \
+    -lmeshTools \
+    -lfieldSources
+
diff --git a/applications/solvers/heatTransfer/buoyantPimpleFoam/UEqn.H b/applications/solvers/heatTransfer/buoyantPimpleFoam/UEqn.H
index 9597f6ac08ddcac42bba70642fe61c013f1c87ed..d8f3bb56ce13a0455c52f7ccbe598f353ea21b2e 100644
--- a/applications/solvers/heatTransfer/buoyantPimpleFoam/UEqn.H
+++ b/applications/solvers/heatTransfer/buoyantPimpleFoam/UEqn.H
@@ -22,6 +22,7 @@
                   - fvc::snGrad(p_rgh)
                 )*mesh.magSf()
             )
+          + sources(rho, U)
         );
         K = 0.5*magSqr(U);
     }
diff --git a/applications/solvers/heatTransfer/buoyantPimpleFoam/buoyantPimpleFoam.C b/applications/solvers/heatTransfer/buoyantPimpleFoam/buoyantPimpleFoam.C
index dfa88b59a9aebb32486807f80705b1381925946c..79d3fb7f2f7668b4dbf1696c25ac7f38a6a0c25b 100644
--- a/applications/solvers/heatTransfer/buoyantPimpleFoam/buoyantPimpleFoam.C
+++ b/applications/solvers/heatTransfer/buoyantPimpleFoam/buoyantPimpleFoam.C
@@ -37,6 +37,7 @@ Description
 #include "rhoThermo.H"
 #include "turbulenceModel.H"
 #include "fixedGradientFvPatchFields.H"
+#include "IObasicSourceList.H"
 #include "pimpleControl.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
diff --git a/applications/solvers/heatTransfer/buoyantPimpleFoam/createFields.H b/applications/solvers/heatTransfer/buoyantPimpleFoam/createFields.H
index f6b183e3fcbf95c113c1a81d7da789c70478ecaa..9838978028edd00a3f9fbce7563b59bce0dcc9cb 100644
--- a/applications/solvers/heatTransfer/buoyantPimpleFoam/createFields.H
+++ b/applications/solvers/heatTransfer/buoyantPimpleFoam/createFields.H
@@ -89,3 +89,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/heatTransfer/buoyantPimpleFoam/pEqn.H b/applications/solvers/heatTransfer/buoyantPimpleFoam/pEqn.H
index d67011dba107dfaa529a214691903cfb81db6dc7..6aa559a763a3df9b1f160a73418b313021fc86f2 100644
--- a/applications/solvers/heatTransfer/buoyantPimpleFoam/pEqn.H
+++ b/applications/solvers/heatTransfer/buoyantPimpleFoam/pEqn.H
@@ -9,7 +9,7 @@
     surfaceScalarField rhorAUf("Dp", fvc::interpolate(rho*rAU));
 
     volVectorField HbyA("HbyA", U);
-    HbyA = rAU*UEqn.H();
+    HbyA = rAU*(UEqn == sources(rho, U))().H();
 
     surfaceScalarField phig(-rhorAUf*ghf*fvc::snGrad(rho)*mesh.magSf());
 
@@ -28,6 +28,8 @@
     (
         fvc::ddt(rho) + psi*correction(fvm::ddt(p_rgh))
       + fvc::div(phiHbyA)
+      ==
+        sources(psi, p, rho.name())
     );
 
     while (pimple.correctNonOrthogonal())
@@ -38,6 +40,8 @@
           - fvm::laplacian(rhorAUf, p_rgh)
         );
 
+        sources.constrain(p_rghEqn, rho.name());
+
         p_rghEqn.solve(mesh.solver(p_rgh.select(pimple.finalInnerIter())));
 
         if (pimple.finalNonOrthogonalIter())
@@ -52,6 +56,7 @@
             // calculated from the relaxed pressure
             U = HbyA + rAU*fvc::reconstruct((phig + p_rghEqn.flux())/rhorAUf);
             U.correctBoundaryConditions();
+            sources.correct(U);
             K = 0.5*magSqr(U);
         }
     }