diff --git a/applications/solvers/incompressible/pimpleFoam/UEqn.H b/applications/solvers/incompressible/pimpleFoam/UEqn.H
index 4d5a99c800..6a424227f3 100644
--- a/applications/solvers/incompressible/pimpleFoam/UEqn.H
+++ b/applications/solvers/incompressible/pimpleFoam/UEqn.H
@@ -5,12 +5,14 @@ tmp<fvVectorMatrix> UEqn
     fvm::ddt(U)
   + fvm::div(phi, U)
   + turbulence->divDevReff(U)
+  ==
+    sources(U)
 );
 
-sources.apply(UEqn());
-
 UEqn().relax();
 
+sources.constrain(UEqn());
+
 volScalarField rAU(1.0/UEqn().A());
 
 if (pimple.momentumPredictor())
diff --git a/applications/solvers/incompressible/simpleFoam/MRFSimpleFoam/UEqn.H b/applications/solvers/incompressible/simpleFoam/MRFSimpleFoam/UEqn.H
index dd23e3d221..8c39a839f7 100644
--- a/applications/solvers/incompressible/simpleFoam/MRFSimpleFoam/UEqn.H
+++ b/applications/solvers/incompressible/simpleFoam/MRFSimpleFoam/UEqn.H
@@ -4,11 +4,13 @@
     (
         fvm::div(phi, U)
       + turbulence->divDevReff(U)
+      ==
+        sources(U)
     );
 
     mrfZones.addCoriolis(UEqn());
 
-    sources.apply(UEqn());
+    sources.constrain(UEqn());
 
     UEqn().relax();
 
diff --git a/applications/solvers/incompressible/simpleFoam/SRFSimpleFoam/UrelEqn.H b/applications/solvers/incompressible/simpleFoam/SRFSimpleFoam/UrelEqn.H
index 94380de28c..fb605fffd3 100644
--- a/applications/solvers/incompressible/simpleFoam/SRFSimpleFoam/UrelEqn.H
+++ b/applications/solvers/incompressible/simpleFoam/SRFSimpleFoam/UrelEqn.H
@@ -5,10 +5,12 @@
         fvm::div(phi, Urel)
       + turbulence->divDevReff(Urel)
       + SRF->Su()
+     ==
+        sources(Urel)
     );
 
-    sources.apply(UrelEqn());
-
     UrelEqn().relax();
 
+    sources.constrain(UrelEqn());
+
     solve(UrelEqn() == -fvc::grad(p));
diff --git a/applications/solvers/incompressible/simpleFoam/UEqn.H b/applications/solvers/incompressible/simpleFoam/UEqn.H
index df1d53a353..786c4f5919 100644
--- a/applications/solvers/incompressible/simpleFoam/UEqn.H
+++ b/applications/solvers/incompressible/simpleFoam/UEqn.H
@@ -4,10 +4,12 @@
     (
         fvm::div(phi, U)
       + turbulence->divDevReff(U)
+      ==
+        sources(U)
     );
 
-    sources.apply(UEqn());
-
     UEqn().relax();
 
+    sources.constrain(UEqn());
+
     solve(UEqn() == -fvc::grad(p));
diff --git a/applications/solvers/incompressible/simpleFoam/porousSimpleFoam/UEqn.H b/applications/solvers/incompressible/simpleFoam/porousSimpleFoam/UEqn.H
index 2dd998c484..632a160030 100644
--- a/applications/solvers/incompressible/simpleFoam/porousSimpleFoam/UEqn.H
+++ b/applications/solvers/incompressible/simpleFoam/porousSimpleFoam/UEqn.H
@@ -4,12 +4,14 @@
     (
         fvm::div(phi, U)
       + turbulence->divDevReff(U)
+      ==
+        sources(U)
     );
 
-    sources.apply(UEqn());
-
     UEqn().relax();
 
+    sources.constrain(UEqn());
+
     // Include the porous media resistance and solve the momentum equation
     // either implicit in the tensorial resistance or transport using by
     // including the spherical part of the resistance in the momentum diagonal
diff --git a/applications/solvers/lagrangian/LTSReactingParcelFoam/UEqn.H b/applications/solvers/lagrangian/LTSReactingParcelFoam/UEqn.H
index bf6f141875..70f7520190 100644
--- a/applications/solvers/lagrangian/LTSReactingParcelFoam/UEqn.H
+++ b/applications/solvers/lagrangian/LTSReactingParcelFoam/UEqn.H
@@ -7,9 +7,10 @@
      ==
         rho.dimensionedInternalField()*g
       + parcels.SU(U)
+      + sources(rho, U)
     );
 
-    sources.apply(UEqn);
+    sources.constrain(UEqn);
 
     pZones.addResistance(UEqn);
 
diff --git a/applications/solvers/lagrangian/LTSReactingParcelFoam/YEqn.H b/applications/solvers/lagrangian/LTSReactingParcelFoam/YEqn.H
index a5ee090d84..743b1d68ff 100644
--- a/applications/solvers/lagrangian/LTSReactingParcelFoam/YEqn.H
+++ b/applications/solvers/lagrangian/LTSReactingParcelFoam/YEqn.H
@@ -31,9 +31,10 @@ if (solveSpecies)
              ==
                 parcels.SYi(i, Yi)
               + combustion->R(Yi)
+              + sources(Yi)
             );
 
-            sources.apply(YEqn);
+            sources.constrain(YEqn);
 
             YEqn.solve(mesh.solver("Yi"));
 
diff --git a/applications/solvers/lagrangian/LTSReactingParcelFoam/hsEqn.H b/applications/solvers/lagrangian/LTSReactingParcelFoam/hsEqn.H
index c8c634bbd9..151e8e60bb 100644
--- a/applications/solvers/lagrangian/LTSReactingParcelFoam/hsEqn.H
+++ b/applications/solvers/lagrangian/LTSReactingParcelFoam/hsEqn.H
@@ -9,9 +9,10 @@
       + parcels.Sh(hs)
       + radiation->Shs(thermo)
       + combustion->Sh()
+      + sources(rho, hs)
     );
 
-    sources.apply(hsEqn);
+    sources.constrain(hsEqn);
 
     hsEqn.solve();
 
diff --git a/applications/solvers/lagrangian/LTSReactingParcelFoam/pEqn.H b/applications/solvers/lagrangian/LTSReactingParcelFoam/pEqn.H
index 5dc9058d9e..27eb76e42f 100644
--- a/applications/solvers/lagrangian/LTSReactingParcelFoam/pEqn.H
+++ b/applications/solvers/lagrangian/LTSReactingParcelFoam/pEqn.H
@@ -29,10 +29,9 @@
       + fvc::div(phi)
      ==
         parcels.Srho()
+      + sources(psi, p, rho.name())
     );
 
-    sources.apply(pDDtEqn);
-
     while (pimple.correctNonOrthogonal())
     {
         fvScalarMatrix pEqn
@@ -41,6 +40,8 @@
           - fvm::laplacian(rho*rAU, p)
         );
 
+        sources.constrain(pDDtEqn, rho.name());
+
         pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter())));
 
         if (pimple.finalNonOrthogonalIter())
diff --git a/applications/solvers/lagrangian/LTSReactingParcelFoam/rhoEqn.H b/applications/solvers/lagrangian/LTSReactingParcelFoam/rhoEqn.H
index 829ad231e7..88b7a83221 100644
--- a/applications/solvers/lagrangian/LTSReactingParcelFoam/rhoEqn.H
+++ b/applications/solvers/lagrangian/LTSReactingParcelFoam/rhoEqn.H
@@ -36,9 +36,10 @@ Description
       + fvc::div(phi)
       ==
         parcels.Srho(rho)
+      + sources(rho)
     );
 
-    sources.apply(rhoEqn);
+    sources.constrain(rhoEqn);
 
     rhoEqn.solve();
 
diff --git a/applications/solvers/lagrangian/LTSReactingParcelFoam/timeScales.H b/applications/solvers/lagrangian/LTSReactingParcelFoam/timeScales.H
index e61e6cdac7..43ebced834 100644
--- a/applications/solvers/lagrangian/LTSReactingParcelFoam/timeScales.H
+++ b/applications/solvers/lagrangian/LTSReactingParcelFoam/timeScales.H
@@ -76,7 +76,7 @@ Info<< "Time scales min/max:" << endl;
             (
                 DpDt
               + parcels.hsTrans()/(mesh.V()*runTime.deltaT())
-//              + energySource.Su()
+//              + sources(rho, hs)
               + combustion->Sh()()
             )
            /rho
diff --git a/applications/solvers/lagrangian/coalChemistryFoam/UEqn.H b/applications/solvers/lagrangian/coalChemistryFoam/UEqn.H
index 2956bf99d8..b7321872c3 100644
--- a/applications/solvers/lagrangian/coalChemistryFoam/UEqn.H
+++ b/applications/solvers/lagrangian/coalChemistryFoam/UEqn.H
@@ -7,12 +7,13 @@
         rho.dimensionedInternalField()*g
       + coalParcels.SU(U)
       + limestoneParcels.SU(U)
+      + sources(rho, U)
     );
 
-    sources.apply(UEqn);
-
     UEqn.relax();
 
+    sources.constrain(UEqn);
+
     if (pimple.momentumPredictor())
     {
         solve(UEqn == -fvc::grad(p));
diff --git a/applications/solvers/lagrangian/coalChemistryFoam/YEqn.H b/applications/solvers/lagrangian/coalChemistryFoam/YEqn.H
index c6abcb4db5..f07208909c 100644
--- a/applications/solvers/lagrangian/coalChemistryFoam/YEqn.H
+++ b/applications/solvers/lagrangian/coalChemistryFoam/YEqn.H
@@ -30,11 +30,13 @@ tmp<fv::convectionScheme<scalar> > mvConvection
               ==
                 coalParcels.SYi(i, Yi)
               + combustion->R(Yi)
+              + sources(rho, Yi)
             );
 
-            sources.apply(YiEqn);
-
             YiEqn.relax();
+
+            sources.constrain(YiEqn);
+
             YiEqn.solve(mesh.solver("Yi"));
 
             Yi.max(0.0);
diff --git a/applications/solvers/lagrangian/coalChemistryFoam/hsEqn.H b/applications/solvers/lagrangian/coalChemistryFoam/hsEqn.H
index fa9f3a0a44..6df66134c4 100644
--- a/applications/solvers/lagrangian/coalChemistryFoam/hsEqn.H
+++ b/applications/solvers/lagrangian/coalChemistryFoam/hsEqn.H
@@ -6,15 +6,17 @@
       - fvm::laplacian(turbulence->alphaEff(), hs)
      ==
         DpDt
-     +  combustion->Sh()
-     +  coalParcels.Sh(hs)
-     +  limestoneParcels.Sh(hs)
-     +  radiation->Shs(thermo)
+      + combustion->Sh()
+      + coalParcels.Sh(hs)
+      + limestoneParcels.Sh(hs)
+      + radiation->Shs(thermo)
+      + sources(rho, hs)
     );
 
-    sources.apply(hsEqn);
-
     hsEqn.relax();
+
+    sources.constrain(hsEqn);
+
     hsEqn.solve();
 
     thermo.correct();
diff --git a/applications/solvers/lagrangian/coalChemistryFoam/pEqn.H b/applications/solvers/lagrangian/coalChemistryFoam/pEqn.H
index 5ef8472210..539521ca9b 100644
--- a/applications/solvers/lagrangian/coalChemistryFoam/pEqn.H
+++ b/applications/solvers/lagrangian/coalChemistryFoam/pEqn.H
@@ -24,9 +24,10 @@ if (pimple.transonic())
           - fvm::laplacian(rho*rAU, p)
          ==
             coalParcels.Srho()
+          + sources(psi, p, rho.name())
         );
 
-        sources.apply(pEqn);
+        sources.constrain(pEqn, rho.name());
 
         pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter())));
 
@@ -39,8 +40,8 @@ if (pimple.transonic())
 else
 {
     phi =
-        fvc::interpolate(rho)*
-        (
+        fvc::interpolate(rho)
+       *(
             (fvc::interpolate(U) & mesh.Sf())
           + fvc::ddtPhiCorr(rAU, rho, U, phi)
         );
@@ -54,9 +55,10 @@ else
           - fvm::laplacian(rho*rAU, p)
          ==
             coalParcels.Srho()
+          + sources(psi, p, rho.name())
         );
 
-        sources.apply(pEqn);
+        sources.constrain(pEqn, rho.name());
 
         pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter())));
 
diff --git a/applications/solvers/lagrangian/coalChemistryFoam/rhoEqn.H b/applications/solvers/lagrangian/coalChemistryFoam/rhoEqn.H
index 8902cc95a8..c8e8193eb2 100644
--- a/applications/solvers/lagrangian/coalChemistryFoam/rhoEqn.H
+++ b/applications/solvers/lagrangian/coalChemistryFoam/rhoEqn.H
@@ -36,9 +36,10 @@ Description
       + fvc::div(phi)
       ==
         coalParcels.Srho(rho)
+      + sources(rho)
     );
 
-    sources.apply(rhoEqn);
+    sources.constrain(rhoEqn);
 
     rhoEqn.solve();
 }
diff --git a/applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/UEqn.H b/applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/UEqn.H
index a89058d2b5..9bb0de54b1 100644
--- a/applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/UEqn.H
+++ b/applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/UEqn.H
@@ -7,9 +7,10 @@
      ==
         rho.dimensionedInternalField()*g
       + parcels.SU(U)
+      + sources(rho, U)
     );
 
-    sources.apply(UEqn);
+    sources.constrain(UEqn);
 
     pZones.addResistance(UEqn);
 
diff --git a/applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/YEqn.H b/applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/YEqn.H
index 8811a5ffe3..e333a14712 100644
--- a/applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/YEqn.H
+++ b/applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/YEqn.H
@@ -32,9 +32,10 @@ if (solveSpecies)
               ==
                 parcels.SYi(i, Yi)
               + combustion->R(Yi)
+              + sources(rho, Yi)
             );
 
-            sources.apply(YiEqn);
+            sources.constrain(YiEqn);
 
             YiEqn.solve(mesh.solver("Yi"));
 
diff --git a/applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/hsEqn.H b/applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/hsEqn.H
index d15dd3367e..27a9f6c8b3 100644
--- a/applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/hsEqn.H
+++ b/applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/hsEqn.H
@@ -39,9 +39,10 @@
           + parcels.Sh(hs)
           + radiation->Shs(thermo)
           + combustion->Sh()
+          + sources(rho, hs)
         );
 
-        sources.apply(hsEqn);
+        sources.constrain(hsEqn);
 
         hsEqn.solve();
 
diff --git a/applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/pEqn.H b/applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/pEqn.H
index 581c03e104..23645a1b7b 100644
--- a/applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/pEqn.H
+++ b/applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/pEqn.H
@@ -29,9 +29,9 @@
       + fvc::div(phi)
      ==
         parcels.Srho()
+      + sources(psi, p, rho.name())
     );
 
-    sources.apply(pDDtEqn);
 
     while (pimple.correctNonOrthogonal())
     {
@@ -41,6 +41,8 @@
           - fvm::laplacian(rho*rAU, p)
         );
 
+        sources.constrain(pDDtEqn, rho.name());
+
         pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter())));
 
         if (pimple.finalNonOrthogonalIter())
diff --git a/applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/rhoEqn.H b/applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/rhoEqn.H
index d5d36f7c5f..9eacd97597 100644
--- a/applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/rhoEqn.H
+++ b/applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/rhoEqn.H
@@ -36,9 +36,10 @@ Description
       + fvc::div(phi)
       ==
         parcels.Srho(rho)
+      + sources(rho)
     );
 
-    sources.apply(rhoEqn);
+    sources.constrain(rhoEqn);
 
     rhoEqn.solve();
 }
diff --git a/src/finiteVolume/cfdTools/general/fieldSources/basicSource/actuationDiskSource/actuationDiskSourceTemplates.C b/src/finiteVolume/cfdTools/general/fieldSources/basicSource/actuationDiskSource/actuationDiskSourceTemplates.C
index 1b1436767e..dbaf95c4b9 100644
--- a/src/finiteVolume/cfdTools/general/fieldSources/basicSource/actuationDiskSource/actuationDiskSourceTemplates.C
+++ b/src/finiteVolume/cfdTools/general/fieldSources/basicSource/actuationDiskSource/actuationDiskSourceTemplates.C
@@ -55,7 +55,7 @@ void Foam::actuationDiskSource::addActuationDiskAxialInertialResistance
     }
     forAll(cells, i)
     {
-        Usource[cells[i]] += ((Vcells[cells[i]]/V())*T[i]*E) & U[cells[i]];
+        Usource[cells[i]] -= ((Vcells[cells[i]]/V())*T[i]*E) & U[cells[i]];
     }
 }
 
diff --git a/src/finiteVolume/cfdTools/general/fieldSources/basicSource/basicSource/basicSource.H b/src/finiteVolume/cfdTools/general/fieldSources/basicSource/basicSource/basicSource.H
index c70a2ba267..bf61cfe0ca 100644
--- a/src/finiteVolume/cfdTools/general/fieldSources/basicSource/basicSource/basicSource.H
+++ b/src/finiteVolume/cfdTools/general/fieldSources/basicSource/basicSource/basicSource.H
@@ -33,6 +33,8 @@ Description
         duration        1000.0;                 // duration
         selectionMode   cellSet;                // cellSet // points //cellZone
 
+    On evaluation, source expects to be added to the rhs of the equation
+
 
 SourceFiles
     basicSource.C
diff --git a/src/finiteVolume/cfdTools/general/fieldSources/basicSource/basicSource/basicSourceList.C b/src/finiteVolume/cfdTools/general/fieldSources/basicSource/basicSource/basicSourceList.C
index 1fe8bd12c2..819d7b8f6b 100644
--- a/src/finiteVolume/cfdTools/general/fieldSources/basicSource/basicSource/basicSourceList.C
+++ b/src/finiteVolume/cfdTools/general/fieldSources/basicSource/basicSource/basicSourceList.C
@@ -36,7 +36,7 @@ namespace Foam
 }
 
 
-// * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+// * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * //
 
 void Foam::basicSourceList::checkApplied() const
 {
diff --git a/src/finiteVolume/cfdTools/general/fieldSources/basicSource/basicSource/basicSourceList.H b/src/finiteVolume/cfdTools/general/fieldSources/basicSource/basicSource/basicSourceList.H
index 59134e64c2..5e5b5fda14 100644
--- a/src/finiteVolume/cfdTools/general/fieldSources/basicSource/basicSource/basicSourceList.H
+++ b/src/finiteVolume/cfdTools/general/fieldSources/basicSource/basicSource/basicSourceList.H
@@ -36,7 +36,8 @@ SourceFile
 #define basicSourceList_H
 
 #include "PtrList.H"
-#include "DimensionedField.H"
+#include "GeometricField.H"
+#include "fvPatchField.H"
 #include "basicSource.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -52,9 +53,9 @@ class basicSourceList
 :
     public PtrList<basicSource>
 {
-private:
+protected:
 
-    // Private data
+    // Protected data
 
         //- Reference to the mesh database
         const fvMesh& mesh_;
@@ -63,7 +64,7 @@ private:
         label checkTimeIndex_;
 
 
-    // Private Member Functions
+    // Protected Member Functions
 
         //- Check that all sources have been applied
         void checkApplied() const;
@@ -94,11 +95,50 @@ public:
 
     // Member Functions
 
-        // Evaluation
+        // Sources
 
-            //- Apply sources to equation
+            //- Return source for equation
             template<class Type>
-            void apply(fvMatrix<Type>& eqn);
+            tmp<fvMatrix<Type> > operator()
+            (
+                GeometricField<Type, fvPatchField, volMesh>& fld
+            );
+
+            //- Return source for equation with specified name
+            template<class Type>
+            tmp<fvMatrix<Type> > operator()
+            (
+                GeometricField<Type, fvPatchField, volMesh>& fld,
+                const word& fieldName
+            );
+
+            //- Return source for equation
+            template<class Type, class RhoType>
+            tmp<fvMatrix<Type> > operator()
+            (
+                const RhoType& rho,
+                GeometricField<Type, fvPatchField, volMesh>& fld
+            );
+
+            //- Return source for equation with specified name
+            template<class Type, class RhoType>
+            tmp<fvMatrix<Type> > operator()
+            (
+                const RhoType& rho,
+                GeometricField<Type, fvPatchField, volMesh>& fld,
+                const word& fieldName
+            );
+
+
+        // Constraints
+
+            //- Apply constraints to equation
+            template<class Type>
+            void constrain(fvMatrix<Type>& eqn);
+
+            //- Apply constraints to equation with specified name
+            template<class Type>
+            void constrain(fvMatrix<Type>& eqn, const word& fieldName);
 
 
         // I-O
diff --git a/src/finiteVolume/cfdTools/general/fieldSources/basicSource/basicSource/basicSourceListTemplates.C b/src/finiteVolume/cfdTools/general/fieldSources/basicSource/basicSource/basicSourceListTemplates.C
index fef93a9e95..aafa91d098 100644
--- a/src/finiteVolume/cfdTools/general/fieldSources/basicSource/basicSource/basicSourceListTemplates.C
+++ b/src/finiteVolume/cfdTools/general/fieldSources/basicSource/basicSource/basicSourceListTemplates.C
@@ -26,11 +26,30 @@ License
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
 template<class Type>
-void Foam::basicSourceList::apply(fvMatrix<Type>& eqn)
+Foam::tmp<Foam::fvMatrix<Type> > Foam::basicSourceList::operator()
+(
+    GeometricField<Type, fvPatchField, volMesh>& fld
+)
+{
+    return this->operator()(fld, fld.name());
+}
+
+
+template<class Type>
+Foam::tmp<Foam::fvMatrix<Type> > Foam::basicSourceList::operator()
+(
+    GeometricField<Type, fvPatchField, volMesh>& fld,
+    const word& fieldName
+)
 {
     checkApplied();
 
-    const word& fieldName = eqn.psi().name();
+    const dimensionSet ds = fld.dimensions()/dimTime*dimVolume;
+
+    tmp<fvMatrix<Type> > tmtx(new fvMatrix<Type>(fld, ds));
+
+    fvMatrix<Type>& mtx = tmtx();
+
 
     forAll(*this, i)
     {
@@ -50,7 +69,104 @@ void Foam::basicSourceList::apply(fvMatrix<Type>& eqn)
                         << fieldName << endl;
                 }
 
-                source.addSup(eqn, fieldI);
+                source.addSup(mtx, fieldI);
+            }
+        }
+    }
+
+    return tmtx;
+}
+
+
+template<class Type, class RhoType>
+Foam::tmp<Foam::fvMatrix<Type> > Foam::basicSourceList::operator()
+(
+    const RhoType& rho,
+    GeometricField<Type, fvPatchField, volMesh>& fld
+)
+{
+    return this->operator()(rho, fld, fld.name());
+}
+
+
+template<class Type, class RhoType>
+Foam::tmp<Foam::fvMatrix<Type> > Foam::basicSourceList::operator()
+(
+    const RhoType& rho,
+    GeometricField<Type, fvPatchField, volMesh>& fld,
+    const word& fieldName
+)
+{
+    checkApplied();
+
+    const dimensionSet ds = rho.dimensions()*fld.dimensions()/dimTime*dimVolume;
+
+    tmp<fvMatrix<Type> > tmtx(new fvMatrix<Type>(fld, ds));
+
+    fvMatrix<Type>& mtx = tmtx();
+
+
+    forAll(*this, i)
+    {
+        basicSource& source = this->operator[](i);
+
+        label fieldI = source.applyToField(fieldName);
+
+        if (fieldI != -1)
+        {
+            source.setApplied(fieldI);
+
+            if (source.isActive())
+            {
+                if (debug)
+                {
+                    Info<< "Applying source " << source.name() << " to field "
+                        << fieldName << endl;
+                }
+
+                source.addSup(mtx, fieldI);
+            }
+        }
+    }
+
+    return tmtx;
+}
+
+
+template<class Type>
+void Foam::basicSourceList::constrain(fvMatrix<Type>& eqn)
+{
+    constrain(eqn, eqn.psi().name());
+}
+
+
+template<class Type>
+void Foam::basicSourceList::constrain
+(
+    fvMatrix<Type>& eqn,
+    const word& fieldName
+)
+{
+    checkApplied();
+
+    forAll(*this, i)
+    {
+        basicSource& source = this->operator[](i);
+
+        label fieldI = source.applyToField(fieldName);
+
+        if (fieldI != -1)
+        {
+            source.setApplied(fieldI);
+
+            if (source.isActive())
+            {
+                if (debug)
+                {
+                    Info<< "Applying constraint " << source.name()
+                        << " to field " << fieldName << endl;
+                }
+
                 source.setValue(eqn, fieldI);
             }
         }
diff --git a/src/finiteVolume/cfdTools/general/fieldSources/basicSource/explicitSetValue/ExplicitSetValue.C b/src/finiteVolume/cfdTools/general/fieldSources/basicSource/explicitSetValue/ExplicitSetValue.C
index 6a14cb2bff..de957994d2 100644
--- a/src/finiteVolume/cfdTools/general/fieldSources/basicSource/explicitSetValue/ExplicitSetValue.C
+++ b/src/finiteVolume/cfdTools/general/fieldSources/basicSource/explicitSetValue/ExplicitSetValue.C
@@ -34,7 +34,7 @@ template<class Type>
 void Foam::ExplicitSetValue<Type>::setFieldData(const dictionary& dict)
 {
     fieldNames_.setSize(dict.toc().size());
-    fieldData_.setSize(fieldNames_.size());
+    injectionRate_.setSize(fieldNames_.size());
 
     applied_.setSize(fieldNames_.size(), false);
 
@@ -42,7 +42,7 @@ void Foam::ExplicitSetValue<Type>::setFieldData(const dictionary& dict)
     forAllConstIter(dictionary, dict, iter)
     {
         fieldNames_[i] = iter().keyword();
-        dict.lookup(iter().keyword()) >> fieldData_[i];
+        dict.lookup(iter().keyword()) >> injectionRate_[i];
         i++;
     }
 }
@@ -60,7 +60,7 @@ Foam::ExplicitSetValue<Type>::ExplicitSetValue
 )
 :
     basicSource(name, modelType, dict, mesh),
-    fieldData_()
+    injectionRate_()
 {
     read(dict);
 }
@@ -81,32 +81,9 @@ void Foam::ExplicitSetValue<Type>::setValue
             << ">::setValue for source " << name_ << endl;
     }
 
-    DimensionedField<Type, volMesh> rhs
-    (
-        IOobject
-        (
-            name_ + fieldNames_[fieldI] + "rhs",
-            eqn.psi().mesh().time().timeName(),
-            eqn.psi().mesh(),
-            IOobject::NO_READ,
-            IOobject::NO_WRITE,
-            false
-        ),
-        eqn.psi().mesh(),
-        dimensioned<Type>
-        (
-            "zero",
-            dimless,
-            pTraits<Type>::zero
-        )
-    );
-
     List<Type> values(cells_.size());
 
-    forAll(values, i)
-    {
-        values[i] = fieldData_[fieldI];
-    }
+    UIndirectList<Type>(values, cells_) = injectionRate_[fieldI];
 
     eqn.setValues(cells_, values);
 }
diff --git a/src/finiteVolume/cfdTools/general/fieldSources/basicSource/explicitSetValue/ExplicitSetValue.H b/src/finiteVolume/cfdTools/general/fieldSources/basicSource/explicitSetValue/ExplicitSetValue.H
index 58a50ea75c..ec0a9d2ce1 100644
--- a/src/finiteVolume/cfdTools/general/fieldSources/basicSource/explicitSetValue/ExplicitSetValue.H
+++ b/src/finiteVolume/cfdTools/general/fieldSources/basicSource/explicitSetValue/ExplicitSetValue.H
@@ -31,7 +31,7 @@ Description
 
         <Type>ExplicitSetValueCoeffs
         {
-            fieldData
+            injectionRate
             {
                 k           30.7;
                 epsilon     1.5;
@@ -70,7 +70,7 @@ protected:
     // Protected data
 
         //- Source value per field
-        List<Type> fieldData_;
+        List<Type> injectionRate_;
 
 
     // Protected functions
diff --git a/src/finiteVolume/cfdTools/general/fieldSources/basicSource/explicitSetValue/ExplicitSetValueIO.C b/src/finiteVolume/cfdTools/general/fieldSources/basicSource/explicitSetValue/ExplicitSetValueIO.C
index 7fd1d45592..a27190ea49 100644
--- a/src/finiteVolume/cfdTools/general/fieldSources/basicSource/explicitSetValue/ExplicitSetValueIO.C
+++ b/src/finiteVolume/cfdTools/general/fieldSources/basicSource/explicitSetValue/ExplicitSetValueIO.C
@@ -40,7 +40,7 @@ bool Foam::ExplicitSetValue<Type>::read(const dictionary& dict)
 {
     if (basicSource::read(dict))
     {
-        setFieldData(coeffs_.subDict("fieldData"));
+        setFieldData(coeffs_.subDict("injectionRate"));
         return true;
     }
     else
diff --git a/src/finiteVolume/cfdTools/general/fieldSources/basicSource/explicitSource/ExplicitSource.C b/src/finiteVolume/cfdTools/general/fieldSources/basicSource/explicitSource/ExplicitSource.C
index b95cd7ba0e..64ed783cb9 100644
--- a/src/finiteVolume/cfdTools/general/fieldSources/basicSource/explicitSource/ExplicitSource.C
+++ b/src/finiteVolume/cfdTools/general/fieldSources/basicSource/explicitSource/ExplicitSource.C
@@ -88,7 +88,7 @@ template<class Type>
 void Foam::ExplicitSource<Type>::setFieldData(const dictionary& dict)
 {
     fieldNames_.setSize(dict.toc().size());
-    fieldData_.setSize(fieldNames_.size());
+    injectionRate_.setSize(fieldNames_.size());
 
     applied_.setSize(fieldNames_.size(), false);
 
@@ -96,7 +96,7 @@ void Foam::ExplicitSource<Type>::setFieldData(const dictionary& dict)
     forAllConstIter(dictionary, dict, iter)
     {
         fieldNames_[i] = iter().keyword();
-        dict.lookup(iter().keyword()) >> fieldData_[i];
+        dict.lookup(iter().keyword()) >> injectionRate_[i];
         i++;
     }
 
@@ -122,7 +122,7 @@ Foam::ExplicitSource<Type>::ExplicitSource
     basicSource(name, modelType, dict, mesh),
     volumeMode_(vmAbsolute),
     VDash_(1.0),
-    fieldData_()
+    injectionRate_()
 {
     read(dict);
 }
@@ -137,7 +137,7 @@ void Foam::ExplicitSource<Type>::addSup
     const label fieldI
 )
 {
-//    if (debug)
+    if (debug)
     {
         Info<< "ExplicitSource<"<< pTraits<Type>::typeName
             << ">::addSup for source " << name_ << endl;
@@ -163,12 +163,9 @@ void Foam::ExplicitSource<Type>::addSup
         false
     );
 
-    forAll(cells_, i)
-    {
-        Su[cells_[i]] = fieldData_[fieldI]/VDash_;
-    }
+    UIndirectList<Type>(Su, cells_) = injectionRate_[fieldI]/VDash_;
 
-    eqn -= Su;
+    eqn += Su;
 }
 
 
diff --git a/src/finiteVolume/cfdTools/general/fieldSources/basicSource/explicitSource/ExplicitSource.H b/src/finiteVolume/cfdTools/general/fieldSources/basicSource/explicitSource/ExplicitSource.H
index 7248b5ef26..e53e6004be 100644
--- a/src/finiteVolume/cfdTools/general/fieldSources/basicSource/explicitSource/ExplicitSource.H
+++ b/src/finiteVolume/cfdTools/general/fieldSources/basicSource/explicitSource/ExplicitSource.H
@@ -32,7 +32,7 @@ Description
         <Type>ExplicitSourceCoeffs
         {
             volumeMode      absolute; // specific
-            fieldData
+            injectionRate
             {
                 k           30.7;
                 epsilon     1.5;
@@ -111,7 +111,7 @@ protected:
         scalar VDash_;
 
         //- Source field values
-        List<Type> fieldData_;
+        List<Type> injectionRate_;
 
 
     // Protected functions
@@ -152,7 +152,7 @@ public:
             inline const volumeModeType& volumeMode() const;
 
             //- Return const access to the source field values
-            inline const List<Type>& fieldData() const;
+            inline const List<Type>& injectionRate() const;
 
 
         // Edit
@@ -161,7 +161,7 @@ public:
             inline volumeModeType& volumeMode();
 
             //- Return access to the source field values
-            inline List<Type>& fieldData();
+            inline List<Type>& injectionRate();
 
 
         // Evaluation
diff --git a/src/finiteVolume/cfdTools/general/fieldSources/basicSource/explicitSource/ExplicitSourceI.H b/src/finiteVolume/cfdTools/general/fieldSources/basicSource/explicitSource/ExplicitSourceI.H
index 14276e0197..e3f4b8aa63 100644
--- a/src/finiteVolume/cfdTools/general/fieldSources/basicSource/explicitSource/ExplicitSourceI.H
+++ b/src/finiteVolume/cfdTools/general/fieldSources/basicSource/explicitSource/ExplicitSourceI.H
@@ -36,9 +36,9 @@ Foam::ExplicitSource<Type>::volumeMode() const
 
 
 template<class Type>
-inline const Foam::List<Type>& Foam::ExplicitSource<Type>::fieldData() const
+inline const Foam::List<Type>& Foam::ExplicitSource<Type>::injectionRate() const
 {
-    return fieldData_;
+    return injectionRate_;
 }
 
 
@@ -51,9 +51,9 @@ Foam::ExplicitSource<Type>::volumeMode()
 
 
 template<class Type>
-inline Foam::List<Type>& Foam::ExplicitSource<Type>::fieldData()
+inline Foam::List<Type>& Foam::ExplicitSource<Type>::injectionRate()
 {
-    return fieldData_;
+    return injectionRate_;
 }
 
 
diff --git a/src/finiteVolume/cfdTools/general/fieldSources/basicSource/explicitSource/ExplicitSourceIO.C b/src/finiteVolume/cfdTools/general/fieldSources/basicSource/explicitSource/ExplicitSourceIO.C
index d181f3f207..d59c3f35cf 100644
--- a/src/finiteVolume/cfdTools/general/fieldSources/basicSource/explicitSource/ExplicitSourceIO.C
+++ b/src/finiteVolume/cfdTools/general/fieldSources/basicSource/explicitSource/ExplicitSourceIO.C
@@ -40,7 +40,7 @@ bool Foam::ExplicitSource<Type>::read(const dictionary& dict)
 {
     if (basicSource::read(dict))
     {
-        setFieldData(coeffs_.subDict("fieldData"));
+        setFieldData(coeffs_.subDict("injectionRate"));
         volumeMode_ = wordToVolumeModeType(coeffs_.lookup("volumeMode"));
 
         return true;
diff --git a/src/finiteVolume/cfdTools/general/fieldSources/basicSource/radialActuationDiskSource/radialActuationDiskSourceTemplates.C b/src/finiteVolume/cfdTools/general/fieldSources/basicSource/radialActuationDiskSource/radialActuationDiskSourceTemplates.C
index fce5c0a170..127263c94f 100644
--- a/src/finiteVolume/cfdTools/general/fieldSources/basicSource/radialActuationDiskSource/radialActuationDiskSourceTemplates.C
+++ b/src/finiteVolume/cfdTools/general/fieldSources/basicSource/radialActuationDiskSource/radialActuationDiskSourceTemplates.C
@@ -79,7 +79,7 @@ addRadialActuationDiskAxialInertialResistance
 
     forAll(cells, i)
     {
-        Usource[cells[i]] += ((Vcells[cells[i]]/V_)*Tr[i]*E) & U[cells[i]];
+        Usource[cells[i]] -= ((Vcells[cells[i]]/V_)*Tr[i]*E) & U[cells[i]];
     }
 
     if (debug)
diff --git a/src/finiteVolume/cfdTools/general/fieldSources/basicSource/rotorDiskSource/rotorDiskSource.C b/src/finiteVolume/cfdTools/general/fieldSources/basicSource/rotorDiskSource/rotorDiskSource.C
index 02f29fcd77..00e9b59ff4 100644
--- a/src/finiteVolume/cfdTools/general/fieldSources/basicSource/rotorDiskSource/rotorDiskSource.C
+++ b/src/finiteVolume/cfdTools/general/fieldSources/basicSource/rotorDiskSource/rotorDiskSource.C
@@ -399,7 +399,7 @@ Foam::rotorDiskSource::~rotorDiskSource()
 
 void Foam::rotorDiskSource::addSup(fvMatrix<vector>& eqn, const label fieldI)
 {
-    // add source to lhs of eqn
+    // add source to rhs of eqn
 
     const volVectorField& U = eqn.psi();
 
@@ -410,7 +410,7 @@ void Foam::rotorDiskSource::addSup(fvMatrix<vector>& eqn, const label fieldI)
         const volScalarField& rho =
             mesh_.lookupObject<volScalarField>(rhoName_);
 
-        eqn += calculateForces
+        eqn -= calculateForces
             (
                 rho.internalField(),
                 inflowVelocity(U),
@@ -419,7 +419,7 @@ void Foam::rotorDiskSource::addSup(fvMatrix<vector>& eqn, const label fieldI)
     }
     else
     {
-        eqn += calculateForces
+        eqn -= calculateForces
             (
                 oneField(),
                 inflowVelocity(U),
diff --git a/tutorials/lagrangian/coalChemistryFoam/simplifiedSiwek/constant/sourcesProperties b/tutorials/lagrangian/coalChemistryFoam/simplifiedSiwek/constant/sourcesProperties
index 0c93adb632..06b09abe24 100644
--- a/tutorials/lagrangian/coalChemistryFoam/simplifiedSiwek/constant/sourcesProperties
+++ b/tutorials/lagrangian/coalChemistryFoam/simplifiedSiwek/constant/sourcesProperties
@@ -27,7 +27,7 @@ source1
     scalarExplicitSourceCoeffs
     {
         volumeMode      absolute;
-        fieldData
+        injectionRate
         {
             hs              20000;
         }
diff --git a/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/filter/constant/sourcesProperties b/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/filter/constant/sourcesProperties
index 0bcf0c982b..ab2601c34c 100644
--- a/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/filter/constant/sourcesProperties
+++ b/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/filter/constant/sourcesProperties
@@ -30,11 +30,10 @@ massSource1
     scalarExplicitSourceCoeffs
     {
         volumeMode      absolute;
-        fieldData
+        injectionRate
         {
             rho         1e-4; // kg/s
             H2O         1e-4; // kg/s
-            p           1e-4; // kg/s
         }
     }
 }
@@ -55,7 +54,7 @@ momentumSource1
     vectorExplicitSourceCoeffs
     {
         volumeMode      absolute;
-        fieldData
+        injectionRate
         {
             U           (0 0.005 0);
         }
@@ -78,7 +77,7 @@ energySource1
     scalarExplicitSourceCoeffs
     {
         volumeMode      absolute;
-        fieldData
+        injectionRate
         {
             hs          10;
         }