diff --git a/applications/solvers/combustion/XiFoam/EaEqn.H b/applications/solvers/combustion/XiFoam/EaEqn.H
index 6840c4eb6d..5493722c93 100644
--- a/applications/solvers/combustion/XiFoam/EaEqn.H
+++ b/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();
diff --git a/applications/solvers/combustion/XiFoam/Make/options b/applications/solvers/combustion/XiFoam/Make/options
index 3a07a7de58..fdebdc8ce8 100644
--- a/applications/solvers/combustion/XiFoam/Make/options
+++ b/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)/meshTools/lnInclude \
-    -I$(LIB_SRC)/finiteVolume/lnInclude
+    -I$(LIB_SRC)/thermophysicalModels/laminarFlameSpeed/lnInclude
 
 EXE_LIBS = \
+    -lfiniteVolume \
+    -lfieldSources \
+    -lsampling \
+    -lmeshTools \
     -lengine \
     -lcompressibleTurbulenceModel \
     -lcompressibleRASModels \
@@ -16,6 +22,4 @@ EXE_LIBS = \
     -lfluidThermophysicalModels \
     -lreactionThermophysicalModels \
     -lspecie \
-    -llaminarFlameSpeedModels \
-    -lmeshTools \
-    -lfiniteVolume
+    -llaminarFlameSpeedModels
diff --git a/applications/solvers/combustion/XiFoam/UEqn.H b/applications/solvers/combustion/XiFoam/UEqn.H
index b9bc567aae..643a99f134 100644
--- a/applications/solvers/combustion/XiFoam/UEqn.H
+++ b/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));
diff --git a/applications/solvers/combustion/XiFoam/XiFoam.C b/applications/solvers/combustion/XiFoam/XiFoam.C
index 1573ede331..a3e90f7169 100644
--- a/applications/solvers/combustion/XiFoam/XiFoam.C
+++ b/applications/solvers/combustion/XiFoam/XiFoam.C
@@ -56,6 +56,7 @@ Description
 #include "ignition.H"
 #include "Switch.H"
 #include "pimpleControl.H"
+#include "IObasicSourceList.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
diff --git a/applications/solvers/combustion/XiFoam/bEqn.H b/applications/solvers/combustion/XiFoam/bEqn.H
index ff19598ba1..49c2c63ae7 100644
--- a/applications/solvers/combustion/XiFoam/bEqn.H
+++ b/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
diff --git a/applications/solvers/combustion/XiFoam/createFields.H b/applications/solvers/combustion/XiFoam/createFields.H
index 17103885e5..3b7cd9b7ae 100644
--- a/applications/solvers/combustion/XiFoam/createFields.H
+++ b/applications/solvers/combustion/XiFoam/createFields.H
@@ -138,3 +138,5 @@
     fields.add(b);
     fields.add(thermo.he());
     fields.add(thermo.heu());
+
+    IObasicSourceList sources(mesh);
diff --git a/applications/solvers/combustion/XiFoam/ftEqn.H b/applications/solvers/combustion/XiFoam/ftEqn.H
index 96dbc3f06b..da1d76fe69 100644
--- a/applications/solvers/combustion/XiFoam/ftEqn.H
+++ b/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();
 }
diff --git a/applications/solvers/combustion/XiFoam/pEqn.H b/applications/solvers/combustion/XiFoam/pEqn.H
index 1543ca539a..99d445e6e9 100644
--- a/applications/solvers/combustion/XiFoam/pEqn.H
+++ b/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, rho.name());
+
         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, rho.name());
+
         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())
diff --git a/applications/solvers/combustion/engineFoam/Make/options b/applications/solvers/combustion/engineFoam/Make/options
index ff38fdf051..c7b3219f2d 100644
--- a/applications/solvers/combustion/engineFoam/Make/options
+++ b/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)/finiteVolume/lnInclude
+    -I$(LIB_SRC)/thermophysicalModels/laminarFlameSpeed/lnInclude
 
 EXE_LIBS = \
+    -lfiniteVolume \
+    -lfieldSources \
+    -lsampling \
+    -lmeshTools \
     -lengine \
     -lcompressibleTurbulenceModel \
     -lcompressibleRASModels \
@@ -16,5 +23,4 @@ EXE_LIBS = \
     -lfluidThermophysicalModels \
     -lreactionThermophysicalModels \
     -lspecie \
-    -llaminarFlameSpeedModels \
-    -lfiniteVolume
+    -llaminarFlameSpeedModels
diff --git a/applications/solvers/combustion/engineFoam/UEqn.H b/applications/solvers/combustion/engineFoam/UEqn.H
index f110051946..de3dc757cc 100644
--- a/applications/solvers/combustion/engineFoam/UEqn.H
+++ b/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())
diff --git a/applications/solvers/combustion/engineFoam/engineFoam.C b/applications/solvers/combustion/engineFoam/engineFoam.C
index 0ee3df6b7b..65cd0356ed 100644
--- a/applications/solvers/combustion/engineFoam/engineFoam.C
+++ b/applications/solvers/combustion/engineFoam/engineFoam.C
@@ -59,6 +59,7 @@ Description
 #include "OFstream.H"
 #include "mathematicalConstants.H"
 #include "pimpleControl.H"
+#include "IObasicSourceList.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
diff --git a/applications/solvers/combustion/engineFoam/pEqn.H b/applications/solvers/combustion/engineFoam/pEqn.H
index 7b005cc311..2a45a1d7e0 100644
--- a/applications/solvers/combustion/engineFoam/pEqn.H
+++ b/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, rho.name());
+
         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, rho.name());
+
         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())
diff --git a/applications/solvers/combustion/fireFoam/Make/options b/applications/solvers/combustion/fireFoam/Make/options
index b7498e50f5..73804d8283 100644
--- a/applications/solvers/combustion/fireFoam/Make/options
+++ b/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 \
diff --git a/applications/solvers/combustion/fireFoam/UEqn.H b/applications/solvers/combustion/fireFoam/UEqn.H
index 31d01e4327..90e633fc36 100644
--- a/applications/solvers/combustion/fireFoam/UEqn.H
+++ b/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
diff --git a/applications/solvers/combustion/fireFoam/YEEqn.H b/applications/solvers/combustion/fireFoam/YEEqn.H
index f979d3bb7d..a23f99dab7 100644
--- a/applications/solvers/combustion/fireFoam/YEEqn.H
+++ b/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();
diff --git a/applications/solvers/combustion/fireFoam/createFields.H b/applications/solvers/combustion/fireFoam/createFields.H
index 8159ee73aa..9ebd05196c 100644
--- a/applications/solvers/combustion/fireFoam/createFields.H
+++ b/applications/solvers/combustion/fireFoam/createFields.H
@@ -146,3 +146,6 @@
     (
         additionalControlsDict.lookup("solvePrimaryRegion")
     );
+
+    IObasicSourceList sources(mesh);
+
diff --git a/applications/solvers/combustion/fireFoam/fireFoam.C b/applications/solvers/combustion/fireFoam/fireFoam.C
index 021ed3c6f0..3f30495a99 100644
--- a/applications/solvers/combustion/fireFoam/fireFoam.C
+++ b/applications/solvers/combustion/fireFoam/fireFoam.C
@@ -40,6 +40,7 @@ Description
 #include "solidChemistryModel.H"
 #include "psiCombustionModel.H"
 #include "pimpleControl.H"
+#include "IObasicSourceList.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
diff --git a/applications/solvers/combustion/fireFoam/pEqn.H b/applications/solvers/combustion/fireFoam/pEqn.H
index 4c33f5633e..4a0145a693 100644
--- a/applications/solvers/combustion/fireFoam/pEqn.H
+++ b/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, rho.name());
+
     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);
     }
 }
 
diff --git a/applications/solvers/combustion/fireFoam/rhoEqn.H b/applications/solvers/combustion/fireFoam/rhoEqn.H
index f939060181..c465309edb 100644
--- a/applications/solvers/combustion/fireFoam/rhoEqn.H
+++ b/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();
 }
 
 // ************************************************************************* //
diff --git a/applications/solvers/combustion/reactingFoam/EEqn.H b/applications/solvers/combustion/reactingFoam/EEqn.H
index e1bf055604..5568c7e5be 100644
--- a/applications/solvers/combustion/reactingFoam/EEqn.H
+++ b/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();
diff --git a/applications/solvers/combustion/reactingFoam/Make/options b/applications/solvers/combustion/reactingFoam/Make/options
index 46b69a63cd..d29dc0ed9a 100644
--- a/applications/solvers/combustion/reactingFoam/Make/options
+++ b/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
diff --git a/applications/solvers/combustion/reactingFoam/UEqn.H b/applications/solvers/combustion/reactingFoam/UEqn.H
index b9bc567aae..643a99f134 100644
--- a/applications/solvers/combustion/reactingFoam/UEqn.H
+++ b/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));
diff --git a/applications/solvers/combustion/reactingFoam/YEqn.H b/applications/solvers/combustion/reactingFoam/YEqn.H
index ccc4b135a5..2f0d29fd84 100644
--- a/applications/solvers/combustion/reactingFoam/YEqn.H
+++ b/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);
diff --git a/applications/solvers/combustion/reactingFoam/createFields.H b/applications/solvers/combustion/reactingFoam/createFields.H
index 32bdd372ee..21769875ea 100644
--- a/applications/solvers/combustion/reactingFoam/createFields.H
+++ b/applications/solvers/combustion/reactingFoam/createFields.H
@@ -99,3 +99,5 @@ volScalarField dQ
     mesh,
     dimensionedScalar("dQ", dimEnergy/dimTime, 0.0)
 );
+
+IObasicSourceList sources(mesh);
diff --git a/applications/solvers/combustion/reactingFoam/pEqn.H b/applications/solvers/combustion/reactingFoam/pEqn.H
index 1543ca539a..34c8bde310 100644
--- a/applications/solvers/combustion/reactingFoam/pEqn.H
+++ b/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, rho.name());
+
         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, rho.name());
+
         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())
diff --git a/applications/solvers/combustion/reactingFoam/reactingFoam.C b/applications/solvers/combustion/reactingFoam/reactingFoam.C
index cd7d371651..c68f94ebf7 100644
--- a/applications/solvers/combustion/reactingFoam/reactingFoam.C
+++ b/applications/solvers/combustion/reactingFoam/reactingFoam.C
@@ -34,6 +34,7 @@ Description
 #include "psiCombustionModel.H"
 #include "multivariateScheme.H"
 #include "pimpleControl.H"
+#include "IObasicSourceList.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
diff --git a/applications/solvers/combustion/rhoReactingFoam/Make/options b/applications/solvers/combustion/rhoReactingFoam/Make/options
index 6dbd401b63..0d1336e1ac 100644
--- a/applications/solvers/combustion/rhoReactingFoam/Make/options
+++ b/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
diff --git a/applications/solvers/combustion/rhoReactingFoam/createFields.H b/applications/solvers/combustion/rhoReactingFoam/createFields.H
index e99639e189..4c517fb646 100644
--- a/applications/solvers/combustion/rhoReactingFoam/createFields.H
+++ b/applications/solvers/combustion/rhoReactingFoam/createFields.H
@@ -101,3 +101,5 @@ volScalarField dQ
     mesh,
     dimensionedScalar("dQ", dimEnergy/dimTime, 0.0)
 );
+
+IObasicSourceList sources(mesh);
diff --git a/applications/solvers/combustion/rhoReactingFoam/pEqn.H b/applications/solvers/combustion/rhoReactingFoam/pEqn.H
index 85dbbf7f7e..55880747b2 100644
--- a/applications/solvers/combustion/rhoReactingFoam/pEqn.H
+++ b/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, rho.name());
+
             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, rho.name());
+
             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())
diff --git a/applications/solvers/combustion/rhoReactingFoam/rhoReactingFoam.C b/applications/solvers/combustion/rhoReactingFoam/rhoReactingFoam.C
index 3ac060612d..0a91aa4b03 100644
--- a/applications/solvers/combustion/rhoReactingFoam/rhoReactingFoam.C
+++ b/applications/solvers/combustion/rhoReactingFoam/rhoReactingFoam.C
@@ -35,6 +35,7 @@ Description
 #include "turbulenceModel.H"
 #include "multivariateScheme.H"
 #include "pimpleControl.H"
+#include "IObasicSourceList.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //