ENH: Updated use of field sources in simple/pimple solvers

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));

applications/solvers/incompressible/pimpleFoam/SRFPimpleFoam/pEqn.H

 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);

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));
 }

applications/solvers/incompressible/pimpleFoam/pEqn.H

-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);

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));
 }

applications/solvers/incompressible/pimpleFoam/pimpleDyMFoam/pEqn.H

-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);

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));
 }

applications/solvers/incompressible/potentialFreeSurfaceFoam/pEqn.H

 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);

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

@@ -34,4 +34,5 @@
     // Momentum corrector
     U -= rAU*fvc::grad(p);
     U.correctBoundaryConditions();
+    sources.correct(U);
 }

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

@@ -34,4 +34,5 @@
     // Momentum corrector
     Urel -= rAUrel*fvc::grad(p);
     Urel.correctBoundaryConditions();
+    sources.correct(Urel);
 }

applications/solvers/incompressible/simpleFoam/pEqn.H

@@ -34,4 +34,5 @@
     // Momentum corrector
     U -= rAU*fvc::grad(p);
     U.correctBoundaryConditions();
+    sources.correct(U);
 }

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));
     }

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);

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);
     }

applications/solvers/lagrangian/coalChemistryFoam/pEqn.H

 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);

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);
     }
 

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);