ENH: Updated solvers using fieldSources

applications/solvers/incompressible/pimpleFoam/UEqn.H

@@ -5,12 +5,13 @@ 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())

applications/solvers/incompressible/simpleFoam/MRFSimpleFoam/UEqn.H

@@ -4,11 +4,12 @@
     (
         fvm::div(phi, U)
       + turbulence->divDevReff(U)
+      + sources(U)
     );
 
     mrfZones.addCoriolis(UEqn());
 
-    sources.apply(UEqn());
+    sources.constrain(UEqn());
 
     UEqn().relax();
 

applications/solvers/incompressible/simpleFoam/SRFSimpleFoam/UrelEqn.H

@@ -5,10 +5,11 @@
         fvm::div(phi, Urel)
       + turbulence->divDevReff(Urel)
       + SRF->Su()
+      + sources(Urel)
     );
 
-    sources.apply(UrelEqn());
-
     UrelEqn().relax();
 
+    sources.constrain(UrelEqn());
+
     solve(UrelEqn() == -fvc::grad(p));

applications/solvers/incompressible/simpleFoam/UEqn.H

@@ -4,10 +4,11 @@
     (
         fvm::div(phi, U)
       + turbulence->divDevReff(U)
+      + sources(U)
     );
 
-    sources.apply(UEqn());
-
     UEqn().relax();
 
+    sources.constrain(UEqn());
+
     solve(UEqn() == -fvc::grad(p));

applications/solvers/incompressible/simpleFoam/porousSimpleFoam/UEqn.H

@@ -4,12 +4,13 @@
     (
         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

applications/solvers/lagrangian/LTSReactingParcelFoam/UEqn.H

@@ -4,12 +4,13 @@
         fvm::ddt(rho, U)
       + fvm::div(phi, U)
       + turbulence->divDevRhoReff(U)
+      + sources(rho, U)
      ==
         rho.dimensionedInternalField()*g
       + parcels.SU(U)
     );
 
-    sources.apply(UEqn);
+    sources.constrain(UEqn);
 
     pZones.addResistance(UEqn);
 

applications/solvers/lagrangian/LTSReactingParcelFoam/YEqn.H

@@ -28,12 +28,13 @@ 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.apply(YEqn);
+            sources.constrain(YEqn);
 
             YEqn.solve(mesh.solver("Yi"));
 

applications/solvers/lagrangian/LTSReactingParcelFoam/hsEqn.H

@@ -4,6 +4,7 @@
         fvm::ddt(rho, hs)
       + mvConvection->fvmDiv(phi, hs)
       - fvm::laplacian(turbulence->alphaEff(), hs)
+      + sources(rho, hs)
      ==
         DpDt
       + parcels.Sh(hs)
@@ -11,7 +12,7 @@
       + combustion->Sh()
     );
 
-    sources.apply(hsEqn);
+    sources.constrain(hsEqn);
 
     hsEqn.solve();
 

applications/solvers/lagrangian/LTSReactingParcelFoam/pEqn.H

@@ -27,12 +27,11 @@
     (
         fvc::ddt(rho) + psi*correction(fvm::ddt(p))
       + fvc::div(phi)
+      + sources(psi, p, rho.name())
      ==
         parcels.Srho()
     );
 
-    sources.apply(pDDtEqn, rho.name());
-
     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())

applications/solvers/lagrangian/LTSReactingParcelFoam/rhoEqn.H

@@ -34,11 +34,12 @@ Description
     (
         fvm::ddt(rho)
       + fvc::div(phi)
+      + sources(rho)
       ==
         parcels.Srho(rho)
     );
 
-    sources.apply(rhoEqn);
+    sources.constrain(rhoEqn);
 
     rhoEqn.solve();
 

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

applications/solvers/lagrangian/coalChemistryFoam/UEqn.H

@@ -3,16 +3,17 @@
         fvm::ddt(rho, U)
       + fvm::div(phi, U)
       + turbulence->divDevRhoReff(U)
+      + sources(rho, U)
      ==
         rho.dimensionedInternalField()*g
       + coalParcels.SU(U)
       + limestoneParcels.SU(U)
     );
 
-    sources.apply(UEqn);
-
     UEqn.relax();
 
+    sources.constrain(UEqn);
+
     if (pimple.momentumPredictor())
     {
         solve(UEqn == -fvc::grad(p));

applications/solvers/lagrangian/coalChemistryFoam/YEqn.H

@@ -27,14 +27,16 @@ 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.apply(YiEqn);
-
             YiEqn.relax();
+
+            sources.constrain(YiEqn);
+
             YiEqn.solve(mesh.solver("Yi"));
 
             Yi.max(0.0);

applications/solvers/lagrangian/coalChemistryFoam/hsEqn.H

@@ -4,6 +4,7 @@
         fvm::ddt(rho, hs)
       + mvConvection->fvmDiv(phi, hs)
       - fvm::laplacian(turbulence->alphaEff(), hs)
+      + sources(rho, hs)
      ==
         DpDt
      +  combustion->Sh()
@@ -12,9 +13,10 @@
      +  radiation->Shs(thermo)
     );
 
-    sources.apply(hsEqn);
-
     hsEqn.relax();
+
+    sources.constrain(hsEqn);
+
     hsEqn.solve();
 
     thermo.correct();

applications/solvers/lagrangian/coalChemistryFoam/pEqn.H

@@ -22,11 +22,12 @@ if (pimple.transonic())
             fvm::ddt(psi, p)
           + fvm::div(phid, p)
           - fvm::laplacian(rho*rAU, p)
+          + sources(psi, p, rho.name())
          ==
             coalParcels.Srho()
         );
 
-        sources.apply(pEqn, rho.name());
+        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)
         );
@@ -52,11 +53,12 @@ else
             fvm::ddt(psi, p)
           + fvc::div(phi)
           - fvm::laplacian(rho*rAU, p)
+          + sources(psi, p, rho.name())
          ==
             coalParcels.Srho()
         );
 
-        sources.apply(pEqn, rho.name());
+        sources.constrain(pEqn, rho.name());
 
         pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter())));
 

applications/solvers/lagrangian/coalChemistryFoam/rhoEqn.H

@@ -34,11 +34,12 @@ Description
     (
         fvm::ddt(rho)
       + fvc::div(phi)
+      + sources(rho)
       ==
         coalParcels.Srho(rho)
     );
 
-    sources.apply(rhoEqn);
+    sources.constrain(rhoEqn);
 
     rhoEqn.solve();
 }

applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/UEqn.H

@@ -4,12 +4,13 @@
         fvm::ddt(rho, U)
       + fvm::div(phi, U)
       + turbulence->divDevRhoReff(U)
+      + sources(rho, U)
      ==
         rho.dimensionedInternalField()*g
       + parcels.SU(U)
     );
 
-    sources.apply(UEqn);
+    sources.constrain(UEqn);
 
     pZones.addResistance(UEqn);
 

applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/YEqn.H

@@ -29,12 +29,13 @@ 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.apply(YiEqn);
+            sources.constrain(YiEqn);
 
             YiEqn.solve(mesh.solver("Yi"));
 

applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/hsEqn.H

@@ -34,6 +34,7 @@
             fvm::ddt(rho, hs)
           + mvConvection->fvmDiv(phi, hs)
           - fvm::laplacian(turbulence->alphaEff(), hs)
+          + sources(rho, hs)
          ==
             pWork()
           + parcels.Sh(hs)
@@ -41,7 +42,7 @@
           + combustion->Sh()
         );
 
-        sources.apply(hsEqn);
+        sources.constrain(hsEqn);
 
         hsEqn.solve();
 

applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/pEqn.H

@@ -27,11 +27,11 @@
     (
         fvc::ddt(rho) + psi*correction(fvm::ddt(p))
       + fvc::div(phi)
+      + sources(psi, p, rho.name())
      ==
         parcels.Srho()
     );
 
-    sources.apply(pDDtEqn, rho.name());
 
     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())