From 8d85caee1d52a18cbf0b0467b22a0c59a1cd0449 Mon Sep 17 00:00:00 2001
From: andy <andy>
Date: Fri, 19 Oct 2012 13:52:41 +0100
Subject: [PATCH] ENH: Added run-time selectable sources to compressible
 solvers

---
 .../solvers/compressible/rhoPimpleFoam/EEqn.H |  3 +++
 .../compressible/rhoPimpleFoam/Make/options   |  9 +++++--
 .../solvers/compressible/rhoPimpleFoam/UEqn.H |  4 ++-
 .../compressible/rhoPimpleFoam/createFields.H |  3 +++
 .../solvers/compressible/rhoPimpleFoam/pEqn.H | 11 +++++++-
 .../rhoPimpleFoam/rhoPimpleFoam.C             |  1 +
 .../rhoPimpleFoam/rhoPimplecFoam/Make/options |  9 +++++--
 .../rhoPimpleFoam/rhoPimplecFoam/pEqn.H       | 11 +++++++-
 .../rhoPimplecFoam/rhoPimplecFoam.C           |  1 +
 .../rhoPorousMRFLTSPimpleFoam/Make/options    |  8 ++++--
 .../rhoPorousMRFLTSPimpleFoam.C               |  1 +
 .../rhoPorousMRFPimpleFoam/Make/options       |  8 ++++--
 .../rhoPorousMRFPimpleFoam/UEqn.H             |  5 +++-
 .../rhoPorousMRFPimpleFoam/pEqn.H             | 11 +++++++-
 .../rhoPorousMRFPimpleFoam.C                  |  1 +
 .../solvers/compressible/rhoSimpleFoam/EEqn.H |  3 +++
 .../compressible/rhoSimpleFoam/Make/options   |  9 +++++--
 .../compressible/rhoSimpleFoam/createFields.H |  3 +++
 .../solvers/compressible/rhoSimpleFoam/pEqn.H | 15 ++++++++---
 .../rhoPorousMRFSimpleFoam/EEqn.H             |  5 ++--
 .../rhoPorousMRFSimpleFoam/Make/options       | 10 ++++---
 .../rhoPorousMRFSimpleFoam/UEqn.H             |  8 ++++--
 .../rhoPorousMRFSimpleFoam/createFields.H     |  3 +++
 .../rhoPorousMRFSimpleFoam/createZones.H      |  2 +-
 .../rhoPorousMRFSimpleFoam/pEqn.H             | 26 +++++++++++++++----
 .../rhoPorousMRFSimpleFoam.C                  |  3 ++-
 .../rhoSimpleFoam/rhoSimpleFoam.C             |  1 +
 .../rhoSimpleFoam/rhoSimplecFoam/Make/options |  9 +++++--
 .../rhoSimpleFoam/rhoSimplecFoam/pEqn.H       | 11 +++++++-
 .../rhoSimplecFoam/rhoSimplecFoam.C           |  1 +
 30 files changed, 159 insertions(+), 36 deletions(-)

diff --git a/applications/solvers/compressible/rhoPimpleFoam/EEqn.H b/applications/solvers/compressible/rhoPimpleFoam/EEqn.H
index 6288b8a0af0..1cc6ed584c6 100644
--- a/applications/solvers/compressible/rhoPimpleFoam/EEqn.H
+++ b/applications/solvers/compressible/rhoPimpleFoam/EEqn.H
@@ -16,9 +16,12 @@
           : -dpdt
         )
       - fvm::laplacian(turbulence->alphaEff(), he)
+     ==
+        sources(rho, he)
     );
 
     EEqn.relax();
+    sources.constrain(EEqn);
     EEqn.solve();
 
     thermo.correct();
diff --git a/applications/solvers/compressible/rhoPimpleFoam/Make/options b/applications/solvers/compressible/rhoPimpleFoam/Make/options
index f21e7394e4d..c96a146eca2 100644
--- a/applications/solvers/compressible/rhoPimpleFoam/Make/options
+++ b/applications/solvers/compressible/rhoPimpleFoam/Make/options
@@ -2,7 +2,10 @@ EXE_INC = \
     -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
     -I$(LIB_SRC)/turbulenceModels/compressible/turbulenceModel \
     -I$(LIB_SRC)/finiteVolume/cfdTools \
-    -I$(LIB_SRC)/finiteVolume/lnInclude
+    -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/meshTools/lnInclude \
+    -I$(LIB_SRC)/sampling/lnInclude \
+    -I$(LIB_SRC)/fieldSources/lnInclude
 
 EXE_LIBS = \
     -lfluidThermophysicalModels \
@@ -11,4 +14,6 @@ EXE_LIBS = \
     -lcompressibleRASModels \
     -lcompressibleLESModels \
     -lfiniteVolume \
-    -lmeshTools
+    -lmeshTools \
+    -lsampling \
+    -lfieldSources
diff --git a/applications/solvers/compressible/rhoPimpleFoam/UEqn.H b/applications/solvers/compressible/rhoPimpleFoam/UEqn.H
index 0f5ec2487b3..397e8930357 100644
--- a/applications/solvers/compressible/rhoPimpleFoam/UEqn.H
+++ b/applications/solvers/compressible/rhoPimpleFoam/UEqn.H
@@ -9,8 +9,10 @@ tmp<fvVectorMatrix> UEqn
 
 UEqn().relax();
 
+sources.constrain(UEqn());
+
 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/compressible/rhoPimpleFoam/createFields.H b/applications/solvers/compressible/rhoPimpleFoam/createFields.H
index 67cc0c3e456..8b3cfd10b9a 100644
--- a/applications/solvers/compressible/rhoPimpleFoam/createFields.H
+++ b/applications/solvers/compressible/rhoPimpleFoam/createFields.H
@@ -69,3 +69,6 @@
 
     Info<< "Creating field kinetic energy K\n" << endl;
     volScalarField K("K", 0.5*magSqr(U));
+
+    Info<< "Creating sources\n" << endl;
+    IObasicSourceList sources(mesh);
diff --git a/applications/solvers/compressible/rhoPimpleFoam/pEqn.H b/applications/solvers/compressible/rhoPimpleFoam/pEqn.H
index d40c66d4025..5f28a3456d8 100644
--- a/applications/solvers/compressible/rhoPimpleFoam/pEqn.H
+++ b/applications/solvers/compressible/rhoPimpleFoam/pEqn.H
@@ -5,7 +5,7 @@ rho.relax();
 
 volScalarField rAU(1.0/UEqn().A());
 volVectorField HbyA("HbyA", U);
-HbyA = rAU*UEqn().H();
+HbyA = rAU*(UEqn() == sources(rho, U))().H();
 
 if (pimple.nCorrPISO() <= 1)
 {
@@ -33,8 +33,12 @@ if (pimple.transonic())
             fvm::ddt(psi, p)
           + fvm::div(phid, p)
           - fvm::laplacian(Dp, p)
+          ==
+            sources(psi, p, rho.name())
         );
 
+        sources.constrain(pEqn, rho.name());
+
         pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter())));
 
         if (pimple.finalNonOrthogonalIter())
@@ -65,8 +69,12 @@ else
             fvm::ddt(psi, p)
           + fvc::div(phiHbyA)
           - fvm::laplacian(Dp, p)
+          ==
+            sources(psi, p, rho.name())
         );
 
+        sources.constrain(pEqn, rho.name());
+
         pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter())));
 
         if (pimple.finalNonOrthogonalIter())
@@ -92,6 +100,7 @@ Info<< "rho max/min : " << max(rho).value()
 
 U = HbyA - rAU*fvc::grad(p);
 U.correctBoundaryConditions();
+sources.correct(U);
 K = 0.5*magSqr(U);
 
 if (thermo.dpdt())
diff --git a/applications/solvers/compressible/rhoPimpleFoam/rhoPimpleFoam.C b/applications/solvers/compressible/rhoPimpleFoam/rhoPimpleFoam.C
index adb6882268a..f9b8a901f72 100644
--- a/applications/solvers/compressible/rhoPimpleFoam/rhoPimpleFoam.C
+++ b/applications/solvers/compressible/rhoPimpleFoam/rhoPimpleFoam.C
@@ -38,6 +38,7 @@ Description
 #include "turbulenceModel.H"
 #include "bound.H"
 #include "pimpleControl.H"
+#include "IObasicSourceList.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
diff --git a/applications/solvers/compressible/rhoPimpleFoam/rhoPimplecFoam/Make/options b/applications/solvers/compressible/rhoPimpleFoam/rhoPimplecFoam/Make/options
index 6d3c41f0dfa..c78d474b4fb 100644
--- a/applications/solvers/compressible/rhoPimpleFoam/rhoPimplecFoam/Make/options
+++ b/applications/solvers/compressible/rhoPimpleFoam/rhoPimplecFoam/Make/options
@@ -3,7 +3,10 @@ EXE_INC = \
     -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
     -I$(LIB_SRC)/turbulenceModels/compressible/turbulenceModel \
     -I$(LIB_SRC)/finiteVolume/cfdTools \
-    -I$(LIB_SRC)/finiteVolume/lnInclude
+    -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/meshTools/lnInclude \
+    -I$(LIB_SRC)/sampling/lnInclude \
+    -I$(LIB_SRC)/fieldSources/lnInclude
 
 EXE_LIBS = \
     -lfluidThermophysicalModels \
@@ -12,4 +15,6 @@ EXE_LIBS = \
     -lcompressibleRASModels \
     -lcompressibleLESModels \
     -lfiniteVolume \
-    -lmeshTools
+    -lmeshTools \
+    -lsampling \
+    -lfieldSources
diff --git a/applications/solvers/compressible/rhoPimpleFoam/rhoPimplecFoam/pEqn.H b/applications/solvers/compressible/rhoPimpleFoam/rhoPimplecFoam/pEqn.H
index 75ed20e2719..628eb71f296 100644
--- a/applications/solvers/compressible/rhoPimpleFoam/rhoPimplecFoam/pEqn.H
+++ b/applications/solvers/compressible/rhoPimpleFoam/rhoPimplecFoam/pEqn.H
@@ -7,7 +7,7 @@ volScalarField rAU(1.0/UEqn().A());
 volScalarField rAtU(1.0/(1.0/rAU - UEqn().H1()));
 
 volVectorField HbyA("HbyA", U);
-HbyA = rAU*UEqn().H();
+HbyA = rAU*(UEqn() == sources(rho, U))().H();
 
 if (pimple.nCorrPIMPLE() <= 1)
 {
@@ -44,11 +44,15 @@ if (pimple.transonic())
           + fvm::div(phid, p)
           + fvc::div(phic)
           - fvm::laplacian(Dp, p)
+         ==
+            sources(psi, p, rho.name())
         );
 
         // Relax the pressure equation to maintain diagonal dominance
         pEqn.relax();
 
+        sources.constrain(pEqn, rho.name());
+
         pEqn.solve();
 
         if (pimple.finalNonOrthogonalIter())
@@ -81,8 +85,12 @@ else
             fvm::ddt(psi, p)
           + fvc::div(phiHbyA)
           - fvm::laplacian(Dp, p)
+         ==
+            sources(psi, p, rho.name())
         );
 
+        sources.constrain(pEqn, rho.name());
+
         pEqn.solve();
 
         if (pimple.finalNonOrthogonalIter())
@@ -100,6 +108,7 @@ p.relax();
 
 U = HbyA - rAtU*fvc::grad(p);
 U.correctBoundaryConditions();
+sources.correct(U);
 K = 0.5*magSqr(U);
 
 if (thermo.dpdt())
diff --git a/applications/solvers/compressible/rhoPimpleFoam/rhoPimplecFoam/rhoPimplecFoam.C b/applications/solvers/compressible/rhoPimpleFoam/rhoPimplecFoam/rhoPimplecFoam.C
index cf8aed0bcad..7e1664dbeac 100644
--- a/applications/solvers/compressible/rhoPimpleFoam/rhoPimplecFoam/rhoPimplecFoam.C
+++ b/applications/solvers/compressible/rhoPimpleFoam/rhoPimplecFoam/rhoPimplecFoam.C
@@ -38,6 +38,7 @@ Description
 #include "turbulenceModel.H"
 #include "bound.H"
 #include "pimpleControl.H"
+#include "IObasicSourceList.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
diff --git a/applications/solvers/compressible/rhoPimpleFoam/rhoPorousMRFLTSPimpleFoam/Make/options b/applications/solvers/compressible/rhoPimpleFoam/rhoPorousMRFLTSPimpleFoam/Make/options
index 5c83c910d95..2976366c08a 100644
--- a/applications/solvers/compressible/rhoPimpleFoam/rhoPorousMRFLTSPimpleFoam/Make/options
+++ b/applications/solvers/compressible/rhoPimpleFoam/rhoPorousMRFLTSPimpleFoam/Make/options
@@ -5,7 +5,9 @@ EXE_INC = \
     -I$(LIB_SRC)/turbulenceModels/compressible/turbulenceModel \
     -I$(LIB_SRC)/finiteVolume/cfdTools \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
-    -I$(LIB_SRC)/meshTools/lnInclude
+    -I$(LIB_SRC)/meshTools/lnInclude \
+    -I$(LIB_SRC)/sampling/lnInclude \
+    -I$(LIB_SRC)/fieldSources/lnInclude
 
 EXE_LIBS = \
     -lfluidThermophysicalModels \
@@ -14,4 +16,6 @@ EXE_LIBS = \
     -lcompressibleRASModels \
     -lcompressibleLESModels \
     -lfiniteVolume \
-    -lmeshTools
+    -lmeshTools \
+    -lsampling \
+    -lfieldSources
diff --git a/applications/solvers/compressible/rhoPimpleFoam/rhoPorousMRFLTSPimpleFoam/rhoPorousMRFLTSPimpleFoam.C b/applications/solvers/compressible/rhoPimpleFoam/rhoPorousMRFLTSPimpleFoam/rhoPorousMRFLTSPimpleFoam.C
index d729b46c37c..82ba5cbac17 100644
--- a/applications/solvers/compressible/rhoPimpleFoam/rhoPorousMRFLTSPimpleFoam/rhoPorousMRFLTSPimpleFoam.C
+++ b/applications/solvers/compressible/rhoPimpleFoam/rhoPorousMRFLTSPimpleFoam/rhoPorousMRFLTSPimpleFoam.C
@@ -39,6 +39,7 @@ Description
 #include "turbulenceModel.H"
 #include "MRFZones.H"
 #include "IOporosityModelList.H"
+#include "IObasicSourceList.H"
 #include "fvcSmooth.H"
 #include "pimpleControl.H"
 #include "bound.H"
diff --git a/applications/solvers/compressible/rhoPimpleFoam/rhoPorousMRFPimpleFoam/Make/options b/applications/solvers/compressible/rhoPimpleFoam/rhoPorousMRFPimpleFoam/Make/options
index f01ebda533b..c78d474b4fb 100644
--- a/applications/solvers/compressible/rhoPimpleFoam/rhoPorousMRFPimpleFoam/Make/options
+++ b/applications/solvers/compressible/rhoPimpleFoam/rhoPorousMRFPimpleFoam/Make/options
@@ -4,7 +4,9 @@ EXE_INC = \
     -I$(LIB_SRC)/turbulenceModels/compressible/turbulenceModel \
     -I$(LIB_SRC)/finiteVolume/cfdTools \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
-    -I$(LIB_SRC)/meshTools/lnInclude
+    -I$(LIB_SRC)/meshTools/lnInclude \
+    -I$(LIB_SRC)/sampling/lnInclude \
+    -I$(LIB_SRC)/fieldSources/lnInclude
 
 EXE_LIBS = \
     -lfluidThermophysicalModels \
@@ -13,4 +15,6 @@ EXE_LIBS = \
     -lcompressibleRASModels \
     -lcompressibleLESModels \
     -lfiniteVolume \
-    -lmeshTools
+    -lmeshTools \
+    -lsampling \
+    -lfieldSources
diff --git a/applications/solvers/compressible/rhoPimpleFoam/rhoPorousMRFPimpleFoam/UEqn.H b/applications/solvers/compressible/rhoPimpleFoam/rhoPorousMRFPimpleFoam/UEqn.H
index c24b1f587a5..9651610157c 100644
--- a/applications/solvers/compressible/rhoPimpleFoam/rhoPorousMRFPimpleFoam/UEqn.H
+++ b/applications/solvers/compressible/rhoPimpleFoam/rhoPorousMRFPimpleFoam/UEqn.H
@@ -11,9 +11,12 @@ tmp<fvVectorMatrix> UEqn
 UEqn().relax();
 
 mrfZones.addCoriolis(rho, UEqn());
+
 pZones.addResistance(UEqn());
 
+sources.constrain(UEqn());
+
 if (pimple.momentumPredictor())
 {
-    solve(UEqn() == -fvc::grad(p));
+    solve(UEqn() == -fvc::grad(p) + sources(rho, U));
 }
diff --git a/applications/solvers/compressible/rhoPimpleFoam/rhoPorousMRFPimpleFoam/pEqn.H b/applications/solvers/compressible/rhoPimpleFoam/rhoPorousMRFPimpleFoam/pEqn.H
index b68823c6e3d..304798b9662 100644
--- a/applications/solvers/compressible/rhoPimpleFoam/rhoPorousMRFPimpleFoam/pEqn.H
+++ b/applications/solvers/compressible/rhoPimpleFoam/rhoPorousMRFPimpleFoam/pEqn.H
@@ -5,7 +5,7 @@ rho.relax();
 
 volScalarField rAU(1.0/UEqn().A());
 volVectorField HbyA("HbyA", U);
-HbyA = rAU*UEqn().H();
+HbyA = rAU*(UEqn() == sources(rho, U))().H();
 
 if (pimple.nCorrPISO() <= 1)
 {
@@ -34,8 +34,12 @@ if (pimple.transonic())
             fvm::ddt(psi, p)
           + fvm::div(phid, p)
           - fvm::laplacian(Dp, p)
+         ==
+            sources(psi, p, rho.name())
         );
 
+        sources.constrain(pEqn, rho.name());
+
         pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter())));
 
         if (pimple.finalNonOrthogonalIter())
@@ -68,8 +72,12 @@ else
             fvm::ddt(psi, p)
           + fvc::div(phiHbyA)
           - fvm::laplacian(Dp, p)
+         ==
+            sources(psi, p, rho.name())
         );
 
+        sources.constrain(pEqn, rho.name());
+
         pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter())));
 
         if (pimple.finalNonOrthogonalIter())
@@ -94,6 +102,7 @@ Info<< "rho max/min : " << max(rho).value()
 
 U = HbyA - rAU*fvc::grad(p);
 U.correctBoundaryConditions();
+sources.correct(U);
 K = 0.5*magSqr(U);
 
 if (thermo.dpdt())
diff --git a/applications/solvers/compressible/rhoPimpleFoam/rhoPorousMRFPimpleFoam/rhoPorousMRFPimpleFoam.C b/applications/solvers/compressible/rhoPimpleFoam/rhoPorousMRFPimpleFoam/rhoPorousMRFPimpleFoam.C
index 47b3d6c35a6..ac5dc1eadcb 100644
--- a/applications/solvers/compressible/rhoPimpleFoam/rhoPorousMRFPimpleFoam/rhoPorousMRFPimpleFoam.C
+++ b/applications/solvers/compressible/rhoPimpleFoam/rhoPorousMRFPimpleFoam/rhoPorousMRFPimpleFoam.C
@@ -39,6 +39,7 @@ Description
 #include "bound.H"
 #include "MRFZones.H"
 #include "IOporosityModelList.H"
+#include "IObasicSourceList.H"
 #include "pimpleControl.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
diff --git a/applications/solvers/compressible/rhoSimpleFoam/EEqn.H b/applications/solvers/compressible/rhoSimpleFoam/EEqn.H
index e4c79b40cc5..e496906d082 100644
--- a/applications/solvers/compressible/rhoSimpleFoam/EEqn.H
+++ b/applications/solvers/compressible/rhoSimpleFoam/EEqn.H
@@ -10,9 +10,12 @@
           : fvc::div(phi, volScalarField("K", 0.5*magSqr(U)))
         )
       - fvm::laplacian(turbulence->alphaEff(), he)
+     ==
+        sources(rho, he)
     );
 
     EEqn.relax();
+    sources.constrain(EEqn);
     EEqn.solve();
 
     thermo.correct();
diff --git a/applications/solvers/compressible/rhoSimpleFoam/Make/options b/applications/solvers/compressible/rhoSimpleFoam/Make/options
index 339cc53bd9e..741dc4f822b 100644
--- a/applications/solvers/compressible/rhoSimpleFoam/Make/options
+++ b/applications/solvers/compressible/rhoSimpleFoam/Make/options
@@ -3,7 +3,10 @@ EXE_INC = \
     -I$(LIB_SRC)/turbulenceModels \
     -I$(LIB_SRC)/turbulenceModels/compressible/RAS/RASModel \
     -I$(LIB_SRC)/finiteVolume/cfdTools \
-    -I$(LIB_SRC)/finiteVolume/lnInclude
+    -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/sampling/lnInclude \
+    -I$(LIB_SRC)/meshTools/lnInclude \
+    -I$(LIB_SRC)/fieldSources/lnInclude
 
 EXE_LIBS = \
     -lfluidThermophysicalModels \
@@ -11,4 +14,6 @@ EXE_LIBS = \
     -lcompressibleTurbulenceModel \
     -lcompressibleRASModels \
     -lfiniteVolume \
-    -lmeshTools
+    -lsampling \
+    -lmeshTools \
+    -lfieldSources
diff --git a/applications/solvers/compressible/rhoSimpleFoam/createFields.H b/applications/solvers/compressible/rhoSimpleFoam/createFields.H
index dea35b76575..3f4b2fd0217 100644
--- a/applications/solvers/compressible/rhoSimpleFoam/createFields.H
+++ b/applications/solvers/compressible/rhoSimpleFoam/createFields.H
@@ -60,3 +60,6 @@
     );
 
     dimensionedScalar initialMass = fvc::domainIntegrate(rho);
+    
+    Info<< "Creating sources\n" << endl;
+    IObasicSourceList sources(mesh);
diff --git a/applications/solvers/compressible/rhoSimpleFoam/pEqn.H b/applications/solvers/compressible/rhoSimpleFoam/pEqn.H
index 4d7334c5f86..2ea228f172d 100644
--- a/applications/solvers/compressible/rhoSimpleFoam/pEqn.H
+++ b/applications/solvers/compressible/rhoSimpleFoam/pEqn.H
@@ -6,7 +6,7 @@
 
     volScalarField rAU(1.0/UEqn().A());
     volVectorField HbyA("HbyA", U);
-    HbyA = rAU*UEqn().H();
+    HbyA = rAU*(UEqn() == sources(rho, U))().H();
 
     UEqn.clear();
 
@@ -25,12 +25,16 @@
             fvScalarMatrix pEqn
             (
                 fvm::div(phid, p)
-                - fvm::laplacian(rho*rAU, p)
+              - fvm::laplacian(rho*rAU, p)
+              ==
+                sources(psi, p, rho.name())
             );
 
             // Relax the pressure equation to ensure diagonal-dominance
             pEqn.relax();
 
+            sources.constrain(pEqn, rho.name());
+
             pEqn.setReference(pRefCell, pRefValue);
 
             pEqn.solve();
@@ -56,11 +60,15 @@
             fvScalarMatrix pEqn
             (
                 fvc::div(phiHbyA)
-                - fvm::laplacian(rho*rAU, p)
+              - fvm::laplacian(rho*rAU, p)
+              ==
+                sources(psi, p, rho.name())
             );
 
             pEqn.setReference(pRefCell, pRefValue);
 
+            sources.constrain(pEqn, rho.name());
+
             pEqn.solve();
 
             if (simple.finalNonOrthogonalIter())
@@ -78,6 +86,7 @@
 
     U = HbyA - rAU*fvc::grad(p);
     U.correctBoundaryConditions();
+    sources.correct(U);
 
     // For closed-volume cases adjust the pressure and density levels
     // to obey overall mass continuity
diff --git a/applications/solvers/compressible/rhoSimpleFoam/rhoPorousMRFSimpleFoam/EEqn.H b/applications/solvers/compressible/rhoSimpleFoam/rhoPorousMRFSimpleFoam/EEqn.H
index ff467c0382c..e496906d082 100644
--- a/applications/solvers/compressible/rhoSimpleFoam/rhoPorousMRFSimpleFoam/EEqn.H
+++ b/applications/solvers/compressible/rhoSimpleFoam/rhoPorousMRFSimpleFoam/EEqn.H
@@ -10,11 +10,12 @@
           : fvc::div(phi, volScalarField("K", 0.5*magSqr(U)))
         )
       - fvm::laplacian(turbulence->alphaEff(), he)
+     ==
+        sources(rho, he)
     );
 
-    pZones.addEnergySource(thermo, rho, EEqn);
-
     EEqn.relax();
+    sources.constrain(EEqn);
     EEqn.solve();
 
     thermo.correct();
diff --git a/applications/solvers/compressible/rhoSimpleFoam/rhoPorousMRFSimpleFoam/Make/options b/applications/solvers/compressible/rhoSimpleFoam/rhoPorousMRFSimpleFoam/Make/options
index 2d328b3ed1f..4ba7ed1f7d1 100644
--- a/applications/solvers/compressible/rhoSimpleFoam/rhoPorousMRFSimpleFoam/Make/options
+++ b/applications/solvers/compressible/rhoSimpleFoam/rhoPorousMRFSimpleFoam/Make/options
@@ -1,18 +1,20 @@
 EXE_INC = \
     -I.. \
     -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
-    -I$(LIB_SRC)/thermophysicalModels/thermalPorousZone/lnInclude \
     -I$(LIB_SRC)/turbulenceModels \
     -I$(LIB_SRC)/turbulenceModels/compressible/RAS/RASModel \
     -I$(LIB_SRC)/finiteVolume/cfdTools \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
-    -I$(LIB_SRC)/meshTools/lnInclude
+    -I$(LIB_SRC)/sampling/lnInclude \
+    -I$(LIB_SRC)/meshTools/lnInclude \
+    -I$(LIB_SRC)/fieldSources/lnInclude
 
 EXE_LIBS = \
     -lfluidThermophysicalModels \
-    -lthermalPorousZone \
     -lspecie \
     -lcompressibleTurbulenceModel \
     -lcompressibleRASModels \
     -lfiniteVolume \
-    -lmeshTools
+    -lsampling \
+    -lmeshTools \
+    -lfieldSources
diff --git a/applications/solvers/compressible/rhoSimpleFoam/rhoPorousMRFSimpleFoam/UEqn.H b/applications/solvers/compressible/rhoSimpleFoam/rhoPorousMRFSimpleFoam/UEqn.H
index c0e087b536c..52e15e07b51 100644
--- a/applications/solvers/compressible/rhoSimpleFoam/rhoPorousMRFSimpleFoam/UEqn.H
+++ b/applications/solvers/compressible/rhoSimpleFoam/rhoPorousMRFSimpleFoam/UEqn.H
@@ -24,11 +24,13 @@
         trTU = inv(tTU());
         trTU().rename("rAU");
 
+        sources.constrain(UEqn());
+
         volVectorField gradp(fvc::grad(p));
 
         for (int UCorr=0; UCorr<nUCorr; UCorr++)
         {
-            U = trTU() & (UEqn().H() - gradp);
+            U = trTU() & ((UEqn() == sources(rho, U))().H() - gradp);
         }
         U.correctBoundaryConditions();
     }
@@ -36,7 +38,9 @@
     {
         pZones.addResistance(UEqn());
 
-        solve(UEqn() == -fvc::grad(p));
+        sources.constrain(UEqn());
+
+        solve(UEqn() == -fvc::grad(p) + sources(rho, U));
 
         trAU = 1.0/UEqn().A();
         trAU().rename("rAU");
diff --git a/applications/solvers/compressible/rhoSimpleFoam/rhoPorousMRFSimpleFoam/createFields.H b/applications/solvers/compressible/rhoSimpleFoam/rhoPorousMRFSimpleFoam/createFields.H
index 4fff74d224c..ad924283eef 100644
--- a/applications/solvers/compressible/rhoSimpleFoam/rhoPorousMRFSimpleFoam/createFields.H
+++ b/applications/solvers/compressible/rhoSimpleFoam/rhoPorousMRFSimpleFoam/createFields.H
@@ -59,3 +59,6 @@
     );
 
     dimensionedScalar initialMass = fvc::domainIntegrate(rho);
+
+    Info<< "Creating sources\n" << endl;
+    IObasicSourceList sources(mesh);
diff --git a/applications/solvers/compressible/rhoSimpleFoam/rhoPorousMRFSimpleFoam/createZones.H b/applications/solvers/compressible/rhoSimpleFoam/rhoPorousMRFSimpleFoam/createZones.H
index e3cfd61f43e..4eb2cb193cc 100644
--- a/applications/solvers/compressible/rhoSimpleFoam/rhoPorousMRFSimpleFoam/createZones.H
+++ b/applications/solvers/compressible/rhoSimpleFoam/rhoPorousMRFSimpleFoam/createZones.H
@@ -1,7 +1,7 @@
     MRFZones mrfZones(mesh);
     mrfZones.correctBoundaryVelocity(U);
 
-    thermalPorousZones pZones(mesh);
+    IOporosityModelList pZones(mesh);
     Switch pressureImplicitPorosity(false);
 
     // nUCorrectors used for pressureImplicitPorosity
diff --git a/applications/solvers/compressible/rhoSimpleFoam/rhoPorousMRFSimpleFoam/pEqn.H b/applications/solvers/compressible/rhoSimpleFoam/rhoPorousMRFSimpleFoam/pEqn.H
index 81b194823cb..708449eb3d8 100644
--- a/applications/solvers/compressible/rhoSimpleFoam/rhoPorousMRFSimpleFoam/pEqn.H
+++ b/applications/solvers/compressible/rhoSimpleFoam/rhoPorousMRFSimpleFoam/pEqn.H
@@ -4,15 +4,17 @@
     rho = min(rho, rhoMax);
     rho.relax();
 
+    const volScalarField& psi = thermo.psi();
+
     volVectorField HbyA("HbyA", U);
 
     if (pressureImplicitPorosity)
     {
-        HbyA = trTU() & UEqn().H();
+        HbyA = trTU() & (UEqn() == sources(rho, U))().H();
     }
     else
     {
-        HbyA = trAU()*UEqn().H();
+        HbyA = trAU()*(UEqn() == sources(rho, U))().H();
     }
 
     UEqn.clear();
@@ -35,15 +37,29 @@
 
         if (pressureImplicitPorosity)
         {
-            tpEqn = (fvm::laplacian(rho*trTU(), p) == fvc::div(phiHbyA));
+            tpEqn =
+            (
+                fvm::laplacian(rho*trTU(), p)
+              + sources(psi, p, rho.name())
+             ==
+                fvc::div(phiHbyA)
+            );
         }
         else
         {
-            tpEqn = (fvm::laplacian(rho*trAU(), p) == fvc::div(phiHbyA));
+            tpEqn =
+            (
+                fvm::laplacian(rho*trAU(), p)
+              + sources(psi, p, rho.name())
+             ==
+                fvc::div(phiHbyA)
+            );
         }
 
         tpEqn().setReference(pRefCell, pRefValue);
 
+        sources.constrain(tpEqn(), rho.name());
+
         tpEqn().solve();
 
         if (simple.finalNonOrthogonalIter())
@@ -67,12 +83,12 @@
     }
 
     U.correctBoundaryConditions();
+    sources.correct(U);
 
     // For closed-volume cases adjust the pressure and density levels
     // to obey overall mass continuity
     if (closedVolume)
     {
-        const volScalarField& psi = thermo.psi();
         p += (initialMass - fvc::domainIntegrate(psi*p))
             /fvc::domainIntegrate(psi);
     }
diff --git a/applications/solvers/compressible/rhoSimpleFoam/rhoPorousMRFSimpleFoam/rhoPorousMRFSimpleFoam.C b/applications/solvers/compressible/rhoSimpleFoam/rhoPorousMRFSimpleFoam/rhoPorousMRFSimpleFoam.C
index 3898295b915..b552cadfcef 100644
--- a/applications/solvers/compressible/rhoSimpleFoam/rhoPorousMRFSimpleFoam/rhoPorousMRFSimpleFoam.C
+++ b/applications/solvers/compressible/rhoSimpleFoam/rhoPorousMRFSimpleFoam/rhoPorousMRFSimpleFoam.C
@@ -35,7 +35,8 @@ Description
 #include "rhoThermo.H"
 #include "RASModel.H"
 #include "MRFZones.H"
-#include "thermalPorousZones.H"
+#include "IObasicSourceList.H"
+#include "IOporosityModelList.H"
 #include "simpleControl.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
diff --git a/applications/solvers/compressible/rhoSimpleFoam/rhoSimpleFoam.C b/applications/solvers/compressible/rhoSimpleFoam/rhoSimpleFoam.C
index 312196583ea..efc1067db23 100644
--- a/applications/solvers/compressible/rhoSimpleFoam/rhoSimpleFoam.C
+++ b/applications/solvers/compressible/rhoSimpleFoam/rhoSimpleFoam.C
@@ -34,6 +34,7 @@ Description
 #include "psiThermo.H"
 #include "RASModel.H"
 #include "simpleControl.H"
+#include "IObasicSourceList.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
diff --git a/applications/solvers/compressible/rhoSimpleFoam/rhoSimplecFoam/Make/options b/applications/solvers/compressible/rhoSimpleFoam/rhoSimplecFoam/Make/options
index 8a8553525f3..facff4ac649 100644
--- a/applications/solvers/compressible/rhoSimpleFoam/rhoSimplecFoam/Make/options
+++ b/applications/solvers/compressible/rhoSimpleFoam/rhoSimplecFoam/Make/options
@@ -4,11 +4,16 @@ EXE_INC = \
     -I$(LIB_SRC)/turbulenceModels \
     -I$(LIB_SRC)/turbulenceModels/compressible/RAS/RASModel \
     -I$(LIB_SRC)/finiteVolume/cfdTools \
-    -I$(LIB_SRC)/finiteVolume/lnInclude
+    -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/sampling/lnInclude \
+    -I$(LIB_SRC)/meshTools/lnInclude \
+    -I$(LIB_SRC)/fieldSources/lnInclude
 
 EXE_LIBS = \
     -lfluidThermophysicalModels \
     -lspecie \
     -lcompressibleRASModels \
     -lfiniteVolume \
-    -lmeshTools
+    -lsampling \
+    -lmeshTools \
+    -lfieldSources
diff --git a/applications/solvers/compressible/rhoSimpleFoam/rhoSimplecFoam/pEqn.H b/applications/solvers/compressible/rhoSimpleFoam/rhoSimplecFoam/pEqn.H
index 8c7405c346f..1ffacc1afaf 100644
--- a/applications/solvers/compressible/rhoSimpleFoam/rhoSimplecFoam/pEqn.H
+++ b/applications/solvers/compressible/rhoSimpleFoam/rhoSimplecFoam/pEqn.H
@@ -7,7 +7,7 @@ volScalarField rAU(1.0/UEqn().A());
 volScalarField rAtU(1.0/(1.0/rAU - UEqn().H1()));
 
 volVectorField HbyA("HbyA", U);
-HbyA = rAU*UEqn().H();
+HbyA = rAU*(UEqn() == sources(rho, U))().H();
 
 UEqn.clear();
 
@@ -38,11 +38,15 @@ if (simple.transonic())
             fvm::div(phid, p)
           + fvc::div(phic)
           - fvm::laplacian(Dp, p)
+         ==
+            sources(psi, p, rho.name())
         );
 
         // Relax the pressure equation to maintain diagonal dominance
         pEqn.relax();
 
+        sources.constrain(pEqn, rho.name());
+
         pEqn.setReference(pRefCell, pRefValue);
 
         pEqn.solve();
@@ -74,8 +78,12 @@ else
         (
             fvc::div(phiHbyA)
           - fvm::laplacian(Dp, p)
+          ==
+            sources(psi, p, rho.name())
         );
 
+        sources.constrain(pEqn, rho.name());
+
         pEqn.setReference(pRefCell, pRefValue);
 
         pEqn.solve();
@@ -96,6 +104,7 @@ p.relax();
 
 U = HbyA - rAtU*fvc::grad(p);
 U.correctBoundaryConditions();
+sources.correct(U);
 
 // For closed-volume cases adjust the pressure and density levels
 // to obey overall mass continuity
diff --git a/applications/solvers/compressible/rhoSimpleFoam/rhoSimplecFoam/rhoSimplecFoam.C b/applications/solvers/compressible/rhoSimpleFoam/rhoSimplecFoam/rhoSimplecFoam.C
index f707570f60c..41e15dbfa84 100644
--- a/applications/solvers/compressible/rhoSimpleFoam/rhoSimplecFoam/rhoSimplecFoam.C
+++ b/applications/solvers/compressible/rhoSimpleFoam/rhoSimplecFoam/rhoSimplecFoam.C
@@ -36,6 +36,7 @@ Description
 #include "mixedFvPatchFields.H"
 #include "bound.H"
 #include "simpleControl.H"
+#include "IObasicSourceList.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-- 
GitLab