diff --git a/applications/solvers/incompressible/pimpleFoam/SRFPimpleFoam/UrelEqn.H b/applications/solvers/incompressible/pimpleFoam/SRFPimpleFoam/UrelEqn.H
index a86b142276ba8cccd075949e562729427648ed96..8b496a6aed5f73765f5156be78af592d1c97e65d 100644
--- a/applications/solvers/incompressible/pimpleFoam/SRFPimpleFoam/UrelEqn.H
+++ b/applications/solvers/incompressible/pimpleFoam/SRFPimpleFoam/UrelEqn.H
@@ -9,4 +9,6 @@
 
     UrelEqn().relax();
 
-    solve(UrelEqn() == -fvc::grad(p));
+    sources.constrain(UrelEqn());
+
+    solve(UrelEqn() == -fvc::grad(p) + sources(Urel));
diff --git a/applications/solvers/incompressible/pimpleFoam/SRFPimpleFoam/pEqn.H b/applications/solvers/incompressible/pimpleFoam/SRFPimpleFoam/pEqn.H
index df3b09adc642da17702decfcd0e0d33f0296d840..72c0d66290957d3e37016cf535f2a34711b4fe51 100644
--- a/applications/solvers/incompressible/pimpleFoam/SRFPimpleFoam/pEqn.H
+++ b/applications/solvers/incompressible/pimpleFoam/SRFPimpleFoam/pEqn.H
@@ -1,5 +1,5 @@
 volScalarField rAUrel(1.0/UrelEqn().A());
-Urel = rAUrel*UrelEqn().H();
+Urel = rAUrel*(UrelEqn() == sources(Urel))().H();
 
 if (pimple.nCorrPISO() <= 1)
 {
@@ -37,3 +37,4 @@ p.relax();
 // Momentum corrector
 Urel -= rAUrel*fvc::grad(p);
 Urel.correctBoundaryConditions();
+sources.correct(Urel);
diff --git a/applications/solvers/incompressible/pimpleFoam/UEqn.H b/applications/solvers/incompressible/pimpleFoam/UEqn.H
index d7be010a9f91f0515516495b6abf0ea84becab9b..a6aef51dafb35a9fee2033b0a07a38f8404f2296 100644
--- a/applications/solvers/incompressible/pimpleFoam/UEqn.H
+++ b/applications/solvers/incompressible/pimpleFoam/UEqn.H
@@ -5,8 +5,6 @@ tmp<fvVectorMatrix> UEqn
     fvm::ddt(U)
   + fvm::div(phi, U)
   + turbulence->divDevReff(U)
-  ==
-    sources(U)
 );
 
 UEqn().relax();
@@ -17,5 +15,5 @@ volScalarField rAU(1.0/UEqn().A());
 
 if (pimple.momentumPredictor())
 {
-    solve(UEqn() == -fvc::grad(p));
+    solve(UEqn() == -fvc::grad(p) + sources(U));
 }
diff --git a/applications/solvers/incompressible/pimpleFoam/pEqn.H b/applications/solvers/incompressible/pimpleFoam/pEqn.H
index 61acd7d69788d444525800cbce5ccf778af0e461..55062637c215dd4cca03f4f0142505c3240f9f7a 100644
--- a/applications/solvers/incompressible/pimpleFoam/pEqn.H
+++ b/applications/solvers/incompressible/pimpleFoam/pEqn.H
@@ -1,4 +1,4 @@
-U = rAU*UEqn().H();
+U = rAU*(UEqn() == sources(U))().H();
 
 if (pimple.nCorrPISO() <= 1)
 {
@@ -36,3 +36,4 @@ p.relax();
 
 U -= rAU*fvc::grad(p);
 U.correctBoundaryConditions();
+sources.correct(U);
diff --git a/applications/solvers/incompressible/pimpleFoam/pimpleDyMFoam/UEqn.H b/applications/solvers/incompressible/pimpleFoam/pimpleDyMFoam/UEqn.H
index 5a211caa9c2f2e6bbf5614928d1eb6d35e307777..618816f5fe1e3a54d658632ab07d43f12920e9f0 100644
--- a/applications/solvers/incompressible/pimpleFoam/pimpleDyMFoam/UEqn.H
+++ b/applications/solvers/incompressible/pimpleFoam/pimpleDyMFoam/UEqn.H
@@ -9,9 +9,11 @@ tmp<fvVectorMatrix> UEqn
 
 UEqn().relax();
 
+sources.constrain(UEqn());
+
 rAU = 1.0/UEqn().A();
 
 if (pimple.momentumPredictor())
 {
-    solve(UEqn() == -fvc::grad(p));
+    solve(UEqn() == -fvc::grad(p) + sources(U));
 }
diff --git a/applications/solvers/incompressible/pimpleFoam/pimpleDyMFoam/pEqn.H b/applications/solvers/incompressible/pimpleFoam/pimpleDyMFoam/pEqn.H
index aeef26a83770ff7e235bf60f7ae15c0bd82417cb..dd90826e7d0f62b2e08f2d084dd30b99b788c2e6 100644
--- a/applications/solvers/incompressible/pimpleFoam/pimpleDyMFoam/pEqn.H
+++ b/applications/solvers/incompressible/pimpleFoam/pimpleDyMFoam/pEqn.H
@@ -1,4 +1,4 @@
-U = rAU*UEqn().H();
+U = rAU*(UEqn() == sources(U))().H();
 
 if (pimple.nCorrPISO() <= 1)
 {
@@ -46,3 +46,4 @@ fvc::makeRelative(phi, U);
 
 U -= rAU*fvc::grad(p);
 U.correctBoundaryConditions();
+sources.correct(U);
diff --git a/applications/solvers/incompressible/potentialFreeSurfaceFoam/UEqn.H b/applications/solvers/incompressible/potentialFreeSurfaceFoam/UEqn.H
index c869df39a235b7f8bdc94b7a8bd552be584defc7..74e84d9bfc9a9c41ce65b7f570ece27c735bcdb8 100644
--- a/applications/solvers/incompressible/potentialFreeSurfaceFoam/UEqn.H
+++ b/applications/solvers/incompressible/potentialFreeSurfaceFoam/UEqn.H
@@ -3,8 +3,6 @@ tmp<fvVectorMatrix> UEqn
     fvm::ddt(U)
   + fvm::div(phi, U)
   + turbulence->divDevReff(U)
-  ==
-    sources(U)
 );
 
 
@@ -14,5 +12,5 @@ sources.constrain(UEqn());
 
 if (pimple.momentumPredictor())
 {
-    solve(UEqn() == -fvc::grad(p_gh));
+    solve(UEqn() == -fvc::grad(p_gh) + sources(U));
 }
diff --git a/applications/solvers/incompressible/potentialFreeSurfaceFoam/pEqn.H b/applications/solvers/incompressible/potentialFreeSurfaceFoam/pEqn.H
index a5adeb9e9142f00802c3ba62db09f2834cc7f886..d7cf6cd5b17fe7dc27377bb65e5676768542738c 100644
--- a/applications/solvers/incompressible/potentialFreeSurfaceFoam/pEqn.H
+++ b/applications/solvers/incompressible/potentialFreeSurfaceFoam/pEqn.H
@@ -1,7 +1,7 @@
 volScalarField rAU(1.0/UEqn().A());
 surfaceScalarField rAUf(rAU.name() + 'f', fvc::interpolate(rAU));
 
-U = rAU*UEqn().H();
+U = rAU*(UEqn() == sources(U))().H();
 
 if (pimple.nCorrPISO() <= 1)
 {
@@ -41,3 +41,4 @@ p = p_gh + (g & (mesh.C() + zeta - refLevel));
 
 U -= rAU*fvc::grad(p_gh);
 U.correctBoundaryConditions();
+sources.correct(U);
diff --git a/applications/solvers/incompressible/simpleFoam/MRFSimpleFoam/pEqn.H b/applications/solvers/incompressible/simpleFoam/MRFSimpleFoam/pEqn.H
index 2d2540d3038a8e28c966ab847fb5fe6d22798749..93d6537b6a99dce1380fd3ecb50d319a2254c3c7 100644
--- a/applications/solvers/incompressible/simpleFoam/MRFSimpleFoam/pEqn.H
+++ b/applications/solvers/incompressible/simpleFoam/MRFSimpleFoam/pEqn.H
@@ -34,4 +34,5 @@
     // Momentum corrector
     U -= rAU*fvc::grad(p);
     U.correctBoundaryConditions();
+    sources.correct(U);
 }
diff --git a/applications/solvers/incompressible/simpleFoam/SRFSimpleFoam/pEqn.H b/applications/solvers/incompressible/simpleFoam/SRFSimpleFoam/pEqn.H
index 54eab9cba7b6c2aa9e110a534c9f6e00ee33ff41..b9a9079b92ca8f21ac0bf65fbd23c98c29a46508 100644
--- a/applications/solvers/incompressible/simpleFoam/SRFSimpleFoam/pEqn.H
+++ b/applications/solvers/incompressible/simpleFoam/SRFSimpleFoam/pEqn.H
@@ -34,4 +34,5 @@
     // Momentum corrector
     Urel -= rAUrel*fvc::grad(p);
     Urel.correctBoundaryConditions();
+    sources.correct(Urel);
 }
diff --git a/applications/solvers/incompressible/simpleFoam/pEqn.H b/applications/solvers/incompressible/simpleFoam/pEqn.H
index 966dedb61688d8c98d50cfd2b83163ce874903a4..78dd40500babb1cea0a56903032b5bb739f5cda7 100644
--- a/applications/solvers/incompressible/simpleFoam/pEqn.H
+++ b/applications/solvers/incompressible/simpleFoam/pEqn.H
@@ -34,4 +34,5 @@
     // Momentum corrector
     U -= rAU*fvc::grad(p);
     U.correctBoundaryConditions();
+    sources.correct(U);
 }
diff --git a/applications/solvers/lagrangian/LTSReactingParcelFoam/UEqn.H b/applications/solvers/lagrangian/LTSReactingParcelFoam/UEqn.H
index 70f752019068876aae9753ab101b6fd0a3fcb9ef..b381b2901d3776a69c403f623176f49d4761dfbd 100644
--- a/applications/solvers/lagrangian/LTSReactingParcelFoam/UEqn.H
+++ b/applications/solvers/lagrangian/LTSReactingParcelFoam/UEqn.H
@@ -7,7 +7,6 @@
      ==
         rho.dimensionedInternalField()*g
       + parcels.SU(U)
-      + sources(rho, U)
     );
 
     sources.constrain(UEqn);
@@ -16,5 +15,5 @@
 
     if (pimple.momentumPredictor())
     {
-        solve(UEqn == -fvc::grad(p));
+        solve(UEqn == -fvc::grad(p) + sources(rho, U));
     }
diff --git a/applications/solvers/lagrangian/LTSReactingParcelFoam/pEqn.H b/applications/solvers/lagrangian/LTSReactingParcelFoam/pEqn.H
index 141faa30bf2ea0fe7e738720fa5977a22d338bb2..8a77a2656e4e371f6fd0acb45c5662e9a746875a 100644
--- a/applications/solvers/lagrangian/LTSReactingParcelFoam/pEqn.H
+++ b/applications/solvers/lagrangian/LTSReactingParcelFoam/pEqn.H
@@ -6,7 +6,7 @@
     thermo.rho() -= psi*p;
 
     volScalarField rAU(1.0/UEqn.A());
-    U = rAU*UEqn.H();
+    U = rAU*(UEqn == sources(rho, U))().H();
 
     if (pZones.size() > 0)
     {
@@ -60,6 +60,7 @@
 
     U -= rAU*fvc::grad(p);
     U.correctBoundaryConditions();
+    sources.correct(U);
 
     rho = thermo.rho();
     rho = max(rho, rhoMin);
diff --git a/applications/solvers/lagrangian/coalChemistryFoam/UEqn.H b/applications/solvers/lagrangian/coalChemistryFoam/UEqn.H
index fe7b047d9a24e41d4baa7f54324f75e3cc08f160..24f56e83d23232a5a3e97018f420f73cd13f4abc 100644
--- a/applications/solvers/lagrangian/coalChemistryFoam/UEqn.H
+++ b/applications/solvers/lagrangian/coalChemistryFoam/UEqn.H
@@ -7,7 +7,6 @@
         rho.dimensionedInternalField()*g
       + coalParcels.SU(U)
       + limestoneParcels.SU(U)
-      + sources(rho, U)
     );
 
     UEqn.relax();
@@ -16,6 +15,6 @@
 
     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/lagrangian/coalChemistryFoam/pEqn.H b/applications/solvers/lagrangian/coalChemistryFoam/pEqn.H
index 56435d31ce12c0c67ef9b8730b49022048c824ef..cb8553bbfbb4615d6a0934faaaefca80a027f95e 100644
--- a/applications/solvers/lagrangian/coalChemistryFoam/pEqn.H
+++ b/applications/solvers/lagrangian/coalChemistryFoam/pEqn.H
@@ -1,7 +1,7 @@
 rho = thermo.rho();
 
 volScalarField rAU(1.0/UEqn.A());
-U = rAU*UEqn.H();
+U = rAU*(UEqn == sources(rho, U))().H();
 
 if (pimple.transonic())
 {
@@ -74,6 +74,8 @@ else
 
 U -= rAU*fvc::grad(p);
 U.correctBoundaryConditions();
+sources.correct(U);
+
 K = 0.5*magSqr(U);
 
 dpdt = fvc::ddt(p);
diff --git a/applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/UEqn.H b/applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/UEqn.H
index 28f6f81d9f0f2c1c8499d292d5d8473c0233bf84..ce25309532bd6f99d79bc79bb6052107c082540b 100644
--- a/applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/UEqn.H
+++ b/applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/UEqn.H
@@ -7,7 +7,6 @@
      ==
         rho.dimensionedInternalField()*g
       + parcels.SU(U)
-      + sources(rho, U)
     );
 
     UEqn.relax();
@@ -18,7 +17,7 @@
 
     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/lagrangian/porousExplicitSourceReactingParcelFoam/pEqn.H b/applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/pEqn.H
index 080faeccedd2c9cd2d4a85d762510d6ee796ea19..6ed7affe77c42f3280f66f4ec5df56520a2a7b74 100644
--- a/applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/pEqn.H
+++ b/applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/pEqn.H
@@ -6,7 +6,7 @@
     thermo.rho() -= psi*p;
 
     volScalarField rAU(1.0/UEqn.A());
-    U = rAU*UEqn.H();
+    U = rAU*(UEqn == sources(rho, U))().H();
 
     if (pZones.size() > 0)
     {
@@ -59,6 +59,8 @@
 
     U -= rAU*fvc::grad(p);
     U.correctBoundaryConditions();
+    sources.correct(U);
+
     K = 0.5*magSqr(U);
 
     dpdt = fvc::ddt(p);