diff --git a/applications/solvers/incompressible/pimpleFoam/UEqn.H b/applications/solvers/incompressible/pimpleFoam/UEqn.H
index 120c230e17..6a424227f3 100644
--- a/applications/solvers/incompressible/pimpleFoam/UEqn.H
+++ b/applications/solvers/incompressible/pimpleFoam/UEqn.H
@@ -5,7 +5,8 @@ tmp<fvVectorMatrix> UEqn
     fvm::ddt(U)
   + fvm::div(phi, U)
   + turbulence->divDevReff(U)
-  + sources(U)
+  ==
+    sources(U)
 );
 
 UEqn().relax();
diff --git a/applications/solvers/incompressible/simpleFoam/MRFSimpleFoam/UEqn.H b/applications/solvers/incompressible/simpleFoam/MRFSimpleFoam/UEqn.H
index 3c0d01e922..8c39a839f7 100644
--- a/applications/solvers/incompressible/simpleFoam/MRFSimpleFoam/UEqn.H
+++ b/applications/solvers/incompressible/simpleFoam/MRFSimpleFoam/UEqn.H
@@ -4,7 +4,8 @@
     (
         fvm::div(phi, U)
       + turbulence->divDevReff(U)
-      + sources(U)
+      ==
+        sources(U)
     );
 
     mrfZones.addCoriolis(UEqn());
diff --git a/applications/solvers/incompressible/simpleFoam/SRFSimpleFoam/UrelEqn.H b/applications/solvers/incompressible/simpleFoam/SRFSimpleFoam/UrelEqn.H
index ef88991aef..fb605fffd3 100644
--- a/applications/solvers/incompressible/simpleFoam/SRFSimpleFoam/UrelEqn.H
+++ b/applications/solvers/incompressible/simpleFoam/SRFSimpleFoam/UrelEqn.H
@@ -5,7 +5,8 @@
         fvm::div(phi, Urel)
       + turbulence->divDevReff(Urel)
       + SRF->Su()
-      + sources(Urel)
+     ==
+        sources(Urel)
     );
 
     UrelEqn().relax();
diff --git a/applications/solvers/incompressible/simpleFoam/UEqn.H b/applications/solvers/incompressible/simpleFoam/UEqn.H
index 854c86bec1..786c4f5919 100644
--- a/applications/solvers/incompressible/simpleFoam/UEqn.H
+++ b/applications/solvers/incompressible/simpleFoam/UEqn.H
@@ -4,7 +4,8 @@
     (
         fvm::div(phi, U)
       + turbulence->divDevReff(U)
-      + sources(U)
+      ==
+        sources(U)
     );
 
     UEqn().relax();
diff --git a/applications/solvers/incompressible/simpleFoam/porousSimpleFoam/UEqn.H b/applications/solvers/incompressible/simpleFoam/porousSimpleFoam/UEqn.H
index d9127ab754..632a160030 100644
--- a/applications/solvers/incompressible/simpleFoam/porousSimpleFoam/UEqn.H
+++ b/applications/solvers/incompressible/simpleFoam/porousSimpleFoam/UEqn.H
@@ -4,7 +4,8 @@
     (
         fvm::div(phi, U)
       + turbulence->divDevReff(U)
-      + sources(U)
+      ==
+        sources(U)
     );
 
     UEqn().relax();
diff --git a/applications/solvers/lagrangian/LTSReactingParcelFoam/UEqn.H b/applications/solvers/lagrangian/LTSReactingParcelFoam/UEqn.H
index c689f2667d..70f7520190 100644
--- a/applications/solvers/lagrangian/LTSReactingParcelFoam/UEqn.H
+++ b/applications/solvers/lagrangian/LTSReactingParcelFoam/UEqn.H
@@ -4,10 +4,10 @@
         fvm::ddt(rho, U)
       + fvm::div(phi, U)
       + turbulence->divDevRhoReff(U)
-      + sources(rho, U)
      ==
         rho.dimensionedInternalField()*g
       + parcels.SU(U)
+      + sources(rho, U)
     );
 
     sources.constrain(UEqn);
diff --git a/applications/solvers/lagrangian/LTSReactingParcelFoam/YEqn.H b/applications/solvers/lagrangian/LTSReactingParcelFoam/YEqn.H
index 2eeb4f7850..743b1d68ff 100644
--- a/applications/solvers/lagrangian/LTSReactingParcelFoam/YEqn.H
+++ b/applications/solvers/lagrangian/LTSReactingParcelFoam/YEqn.H
@@ -28,10 +28,10 @@ if (solveSpecies)
                 fvm::ddt(rho, Yi)
               + mvConvection->fvmDiv(phi, Yi)
               - fvm::laplacian(turbulence->muEff(), Yi)
-              + sources(Yi)
              ==
                 parcels.SYi(i, Yi)
               + combustion->R(Yi)
+              + sources(Yi)
             );
 
             sources.constrain(YEqn);
diff --git a/applications/solvers/lagrangian/LTSReactingParcelFoam/hsEqn.H b/applications/solvers/lagrangian/LTSReactingParcelFoam/hsEqn.H
index 9ed83206bc..151e8e60bb 100644
--- a/applications/solvers/lagrangian/LTSReactingParcelFoam/hsEqn.H
+++ b/applications/solvers/lagrangian/LTSReactingParcelFoam/hsEqn.H
@@ -4,12 +4,12 @@
         fvm::ddt(rho, hs)
       + mvConvection->fvmDiv(phi, hs)
       - fvm::laplacian(turbulence->alphaEff(), hs)
-      + sources(rho, hs)
      ==
         DpDt
       + parcels.Sh(hs)
       + radiation->Shs(thermo)
       + combustion->Sh()
+      + sources(rho, hs)
     );
 
     sources.constrain(hsEqn);
diff --git a/applications/solvers/lagrangian/LTSReactingParcelFoam/pEqn.H b/applications/solvers/lagrangian/LTSReactingParcelFoam/pEqn.H
index bb7b7be815..27eb76e42f 100644
--- a/applications/solvers/lagrangian/LTSReactingParcelFoam/pEqn.H
+++ b/applications/solvers/lagrangian/LTSReactingParcelFoam/pEqn.H
@@ -27,9 +27,9 @@
     (
         fvc::ddt(rho) + psi*correction(fvm::ddt(p))
       + fvc::div(phi)
-      + sources(psi, p, rho.name())
      ==
         parcels.Srho()
+      + sources(psi, p, rho.name())
     );
 
     while (pimple.correctNonOrthogonal())
diff --git a/applications/solvers/lagrangian/LTSReactingParcelFoam/rhoEqn.H b/applications/solvers/lagrangian/LTSReactingParcelFoam/rhoEqn.H
index f87fc5211d..88b7a83221 100644
--- a/applications/solvers/lagrangian/LTSReactingParcelFoam/rhoEqn.H
+++ b/applications/solvers/lagrangian/LTSReactingParcelFoam/rhoEqn.H
@@ -34,9 +34,9 @@ Description
     (
         fvm::ddt(rho)
       + fvc::div(phi)
-      + sources(rho)
       ==
         parcels.Srho(rho)
+      + sources(rho)
     );
 
     sources.constrain(rhoEqn);
diff --git a/applications/solvers/lagrangian/coalChemistryFoam/UEqn.H b/applications/solvers/lagrangian/coalChemistryFoam/UEqn.H
index 1211ce88ab..b7321872c3 100644
--- a/applications/solvers/lagrangian/coalChemistryFoam/UEqn.H
+++ b/applications/solvers/lagrangian/coalChemistryFoam/UEqn.H
@@ -3,11 +3,11 @@
         fvm::ddt(rho, U)
       + fvm::div(phi, U)
       + turbulence->divDevRhoReff(U)
-      + sources(rho, U)
      ==
         rho.dimensionedInternalField()*g
       + coalParcels.SU(U)
       + limestoneParcels.SU(U)
+      + sources(rho, U)
     );
 
     UEqn.relax();
diff --git a/applications/solvers/lagrangian/coalChemistryFoam/YEqn.H b/applications/solvers/lagrangian/coalChemistryFoam/YEqn.H
index c2e6d6944e..f07208909c 100644
--- a/applications/solvers/lagrangian/coalChemistryFoam/YEqn.H
+++ b/applications/solvers/lagrangian/coalChemistryFoam/YEqn.H
@@ -27,10 +27,10 @@ tmp<fv::convectionScheme<scalar> > mvConvection
                 fvm::ddt(rho, Yi)
               + mvConvection->fvmDiv(phi, Yi)
               - fvm::laplacian(turbulence->muEff(), Yi)
-              + sources(rho, Yi)
               ==
                 coalParcels.SYi(i, Yi)
               + combustion->R(Yi)
+              + sources(rho, Yi)
             );
 
             YiEqn.relax();
diff --git a/applications/solvers/lagrangian/coalChemistryFoam/hsEqn.H b/applications/solvers/lagrangian/coalChemistryFoam/hsEqn.H
index 1c4c20d18c..6df66134c4 100644
--- a/applications/solvers/lagrangian/coalChemistryFoam/hsEqn.H
+++ b/applications/solvers/lagrangian/coalChemistryFoam/hsEqn.H
@@ -4,13 +4,13 @@
         fvm::ddt(rho, hs)
       + mvConvection->fvmDiv(phi, hs)
       - fvm::laplacian(turbulence->alphaEff(), hs)
-      + sources(rho, 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)
     );
 
     hsEqn.relax();
diff --git a/applications/solvers/lagrangian/coalChemistryFoam/pEqn.H b/applications/solvers/lagrangian/coalChemistryFoam/pEqn.H
index 245f419bbf..539521ca9b 100644
--- a/applications/solvers/lagrangian/coalChemistryFoam/pEqn.H
+++ b/applications/solvers/lagrangian/coalChemistryFoam/pEqn.H
@@ -22,9 +22,9 @@ if (pimple.transonic())
             fvm::ddt(psi, p)
           + fvm::div(phid, p)
           - fvm::laplacian(rho*rAU, p)
-          + sources(psi, p, rho.name())
          ==
             coalParcels.Srho()
+          + sources(psi, p, rho.name())
         );
 
         sources.constrain(pEqn, rho.name());
@@ -53,9 +53,9 @@ else
             fvm::ddt(psi, p)
           + fvc::div(phi)
           - fvm::laplacian(rho*rAU, p)
-          + sources(psi, p, rho.name())
          ==
             coalParcels.Srho()
+          + sources(psi, p, rho.name())
         );
 
         sources.constrain(pEqn, rho.name());
diff --git a/applications/solvers/lagrangian/coalChemistryFoam/rhoEqn.H b/applications/solvers/lagrangian/coalChemistryFoam/rhoEqn.H
index 336aa8b48c..c8e8193eb2 100644
--- a/applications/solvers/lagrangian/coalChemistryFoam/rhoEqn.H
+++ b/applications/solvers/lagrangian/coalChemistryFoam/rhoEqn.H
@@ -34,9 +34,9 @@ Description
     (
         fvm::ddt(rho)
       + fvc::div(phi)
-      + sources(rho)
       ==
         coalParcels.Srho(rho)
+      + sources(rho)
     );
 
     sources.constrain(rhoEqn);
diff --git a/applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/UEqn.H b/applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/UEqn.H
index 134ec50950..9bb0de54b1 100644
--- a/applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/UEqn.H
+++ b/applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/UEqn.H
@@ -4,10 +4,10 @@
         fvm::ddt(rho, U)
       + fvm::div(phi, U)
       + turbulence->divDevRhoReff(U)
-      + sources(rho, U)
      ==
         rho.dimensionedInternalField()*g
       + parcels.SU(U)
+      + sources(rho, U)
     );
 
     sources.constrain(UEqn);
diff --git a/applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/YEqn.H b/applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/YEqn.H
index 71059f3e63..e333a14712 100644
--- a/applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/YEqn.H
+++ b/applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/YEqn.H
@@ -29,10 +29,10 @@ if (solveSpecies)
                 fvm::ddt(rho, Yi)
               + mvConvection->fvmDiv(phi, Yi)
               - fvm::laplacian(turbulence->muEff(), Yi)
-              + sources(rho, Yi)
               ==
                 parcels.SYi(i, Yi)
               + combustion->R(Yi)
+              + sources(rho, Yi)
             );
 
             sources.constrain(YiEqn);
diff --git a/applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/hsEqn.H b/applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/hsEqn.H
index 27628beb7c..27a9f6c8b3 100644
--- a/applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/hsEqn.H
+++ b/applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/hsEqn.H
@@ -34,12 +34,12 @@
             fvm::ddt(rho, hs)
           + mvConvection->fvmDiv(phi, hs)
           - fvm::laplacian(turbulence->alphaEff(), hs)
-          + sources(rho, hs)
          ==
             pWork()
           + parcels.Sh(hs)
           + radiation->Shs(thermo)
           + combustion->Sh()
+          + sources(rho, hs)
         );
 
         sources.constrain(hsEqn);
diff --git a/applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/pEqn.H b/applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/pEqn.H
index cd2f228659..23645a1b7b 100644
--- a/applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/pEqn.H
+++ b/applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/pEqn.H
@@ -27,9 +27,9 @@
     (
         fvc::ddt(rho) + psi*correction(fvm::ddt(p))
       + fvc::div(phi)
-      + sources(psi, p, rho.name())
      ==
         parcels.Srho()
+      + sources(psi, p, rho.name())
     );
 
 
diff --git a/applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/rhoEqn.H b/applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/rhoEqn.H
index e2816dd18c..9eacd97597 100644
--- a/applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/rhoEqn.H
+++ b/applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/rhoEqn.H
@@ -34,9 +34,9 @@ Description
     (
         fvm::ddt(rho)
       + fvc::div(phi)
-      + sources(rho)
       ==
         parcels.Srho(rho)
+      + sources(rho)
     );
 
     sources.constrain(rhoEqn);
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/explicitSetValue/ExplicitSetValue.C b/src/finiteVolume/cfdTools/general/fieldSources/basicSource/explicitSetValue/ExplicitSetValue.C
index 4ec825da34..de957994d2 100644
--- a/src/finiteVolume/cfdTools/general/fieldSources/basicSource/explicitSetValue/ExplicitSetValue.C
+++ b/src/finiteVolume/cfdTools/general/fieldSources/basicSource/explicitSetValue/ExplicitSetValue.C
@@ -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] = injectionRate_[fieldI];
-    }
+    UIndirectList<Type>(values, cells_) = injectionRate_[fieldI];
 
     eqn.setValues(cells_, values);
 }
diff --git a/src/finiteVolume/cfdTools/general/fieldSources/basicSource/explicitSource/ExplicitSource.C b/src/finiteVolume/cfdTools/general/fieldSources/basicSource/explicitSource/ExplicitSource.C
index b9f8204d87..64ed783cb9 100644
--- a/src/finiteVolume/cfdTools/general/fieldSources/basicSource/explicitSource/ExplicitSource.C
+++ b/src/finiteVolume/cfdTools/general/fieldSources/basicSource/explicitSource/ExplicitSource.C
@@ -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]] = injectionRate_[fieldI]/VDash_;
-    }
+    UIndirectList<Type>(Su, cells_) = injectionRate_[fieldI]/VDash_;
 
-    eqn -= Su;
+    eqn += Su;
 }
 
 
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),