ENH: Added run-time selectable sources to combustion solvers

2025-11-28 03:28:01 +00:00 · 2012-12-06 12:05:45 +00:00
parent 28c883564a
commit b2252f36b1
30 changed files with 172 additions and 21 deletions
--- a/applications/solvers/combustion/XiFoam/EaEqn.H
+++ b/applications/solvers/combustion/XiFoam/EaEqn.H
@ -16,9 +16,13 @@
          : -dpdt
        )
      - fvm::laplacian(turbulence->alphaEff(), hea)
+      + sources(rho, hea)
    );

    EaEqn.relax();
+
+    sources.constrain(EaEqn);
+
    EaEqn.solve();

    thermo.correct();
--- a/applications/solvers/combustion/XiFoam/Make/options
+++ b/applications/solvers/combustion/XiFoam/Make/options
@ -1,14 +1,20 @@
 EXE_INC = \
+    -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/fieldSources/lnInclude \
+    -I$(LIB_SRC)/sampling/lnInclude \
+    -I$(LIB_SRC)/meshTools/lnInclude\
    -I$(LIB_SRC)/engine/lnInclude \
    -I$(LIB_SRC)/thermophysicalModels/specie/lnInclude \
    -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
    -I$(LIB_SRC)/thermophysicalModels/reactionThermo/lnInclude \
    -I$(LIB_SRC)/turbulenceModels/compressible/turbulenceModel \
-    -I$(LIB_SRC)/thermophysicalModels/laminarFlameSpeed/lnInclude \
-    -I$(LIB_SRC)/meshTools/lnInclude \
-    -I$(LIB_SRC)/finiteVolume/lnInclude
+    -I$(LIB_SRC)/thermophysicalModels/laminarFlameSpeed/lnInclude

 EXE_LIBS = \
+    -lfiniteVolume \
+    -lfieldSources \
+    -lsampling \
+    -lmeshTools \
    -lengine \
    -lcompressibleTurbulenceModel \
    -lcompressibleRASModels \
@ -16,6 +22,4 @@ EXE_LIBS = \
    -lfluidThermophysicalModels \
    -lreactionThermophysicalModels \
    -lspecie \
-    -llaminarFlameSpeedModels \
-    -lmeshTools \
-    -lfiniteVolume
+    -llaminarFlameSpeedModels
--- a/applications/solvers/combustion/XiFoam/UEqn.H
+++ b/applications/solvers/combustion/XiFoam/UEqn.H
@ -5,10 +5,13 @@
      + turbulence->divDevRhoReff(U)
     ==
        rho*g
+      + sources(rho, U)
    );

    UEqn.relax();

+    sources.constrain(UEqn);
+
    if (pimple.momentumPredictor())
    {
        solve(UEqn == -fvc::grad(p));
--- a/applications/solvers/combustion/XiFoam/XiFoam.C
+++ b/applications/solvers/combustion/XiFoam/XiFoam.C
@ -56,6 +56,7 @@ Description
 #include "ignition.H"
 #include "Switch.H"
 #include "pimpleControl.H"
+#include "IObasicSourceList.H"

 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

--- a/applications/solvers/combustion/XiFoam/bEqn.H
+++ b/applications/solvers/combustion/XiFoam/bEqn.H
@ -34,7 +34,7 @@ if (ign.ignited())

    // Calculate turbulent flame speed flux
    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-    surfaceScalarField phiSt(fvc::interpolate(rhou*StCorr*Su*Xi)*nf);
+    surfaceScalarField phiSt("phiSt", fvc::interpolate(rhou*StCorr*Su*Xi)*nf);

    scalar StCoNum = max
    (
@ -50,9 +50,11 @@ if (ign.ignited())
    (
        fvm::ddt(rho, b)
      + mvConvection->fvmDiv(phi, b)
-      + fvm::div(phiSt, b, "div(phiSt,b)")
+      + fvm::div(phiSt, b)
      - fvm::Sp(fvc::div(phiSt), b)
      - fvm::laplacian(turbulence->alphaEff(), b)
+     ==
+        sources(rho, b)
    );


@ -64,6 +66,9 @@ if (ign.ignited())
    // Solve for b
    // ~~~~~~~~~~~
    bEqn.relax();
+
+    sources.constrain(bEqn);
+
    bEqn.solve();

    Info<< "min(b) = " << min(b).value() << endl;
@ -153,9 +158,13 @@ if (ign.ignited())
          ==
          - fvm::SuSp(-rho*Rc*Su0/Su, Su)
          - fvm::SuSp(rho*(sigmas + Rc), Su)
+          + sources(rho, Su)
        );

        SuEqn.relax();
+
+        sources.constrain(SuEqn);
+
        SuEqn.solve();

        // Limit the maximum Su
@ -234,9 +243,13 @@ if (ign.ignited())
                ),
                Xi
            )
+          + sources(rho, Xi)
        );

        XiEqn.relax();
+
+        sources.constrain(XiEqn);
+
        XiEqn.solve();

        // Correct boundedness of Xi
--- a/applications/solvers/combustion/XiFoam/createFields.H
+++ b/applications/solvers/combustion/XiFoam/createFields.H
@ -138,3 +138,5 @@
    fields.add(b);
    fields.add(thermo.he());
    fields.add(thermo.heu());
+
+    IObasicSourceList sources(mesh);
--- a/applications/solvers/combustion/XiFoam/ftEqn.H
+++ b/applications/solvers/combustion/XiFoam/ftEqn.H
@ -13,10 +13,16 @@ if (composition.contains("ft"))
 {
    volScalarField& ft = composition.Y("ft");

-    solve
+    fvScalarMatrix ftEqn
    (
        fvm::ddt(rho, ft)
      + mvConvection->fvmDiv(phi, ft)
      - fvm::laplacian(turbulence->alphaEff(), ft)
+     ==
+        sources(rho, ft)
    );
+
+    sources.constrain(ftEqn);
+
+    ftEqn.solve();
 }
--- a/applications/solvers/combustion/XiFoam/pEqn.H
+++ b/applications/solvers/combustion/XiFoam/pEqn.H
@ -16,6 +16,8 @@ if (pimple.transonic())
        )
    );

+    sources.relativeFlux(fvc::interpolate(psi), phid);
+
    while (pimple.correctNonOrthogonal())
    {
        fvScalarMatrix pEqn
@ -23,8 +25,12 @@ if (pimple.transonic())
            fvm::ddt(psi, p)
          + fvm::div(phid, p)
          - fvm::laplacian(rho*rAU, p)
+         ==
+            sources(psi, p, rho.name())
        );

+        sources.constrain(pEqn, rho.name());
+
        pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter())));

        if (pimple.finalNonOrthogonalIter())
@ -45,6 +51,8 @@ else
        )
    );

+    sources.relativeFlux(phiHbyA);
+
    while (pimple.correctNonOrthogonal())
    {
        fvScalarMatrix pEqn
@ -52,8 +60,12 @@ else
            fvm::ddt(psi, p)
          + fvc::div(phiHbyA)
          - fvm::laplacian(rho*rAU, p)
+         ==
+            sources(psi, p, rho.name())
        );

+        sources.constrain(pEqn, rho.name());
+
        pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter())));

        if (pimple.finalNonOrthogonalIter())
@ -68,6 +80,7 @@ else

 U = HbyA - rAU*fvc::grad(p);
 U.correctBoundaryConditions();
+sources.correct(U);
 K = 0.5*magSqr(U);

 if (thermo.dpdt())
--- a/applications/solvers/combustion/engineFoam/Make/options
+++ b/applications/solvers/combustion/engineFoam/Make/options
@ -1,14 +1,21 @@
 EXE_INC = \
-    -I../XiFoam \
+    -I$(FOAM_SOLVERS)/combustion/XiFoam \
+    -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/fieldSources/lnInclude \
+    -I$(LIB_SRC)/sampling/lnInclude \
+    -I$(LIB_SRC)/meshTools/lnInclude \
    -I$(LIB_SRC)/engine/lnInclude \
    -I$(LIB_SRC)/thermophysicalModels/specie/lnInclude \
    -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
    -I$(LIB_SRC)/thermophysicalModels/reactionThermo/lnInclude \
    -I$(LIB_SRC)/turbulenceModels/compressible/turbulenceModel \
-    -I$(LIB_SRC)/thermophysicalModels/laminarFlameSpeed/lnInclude \
-    -I$(LIB_SRC)/finiteVolume/lnInclude
+    -I$(LIB_SRC)/thermophysicalModels/laminarFlameSpeed/lnInclude

 EXE_LIBS = \
+    -lfiniteVolume \
+    -lfieldSources \
+    -lsampling \
+    -lmeshTools \
    -lengine \
    -lcompressibleTurbulenceModel \
    -lcompressibleRASModels \
@ -16,5 +23,4 @@ EXE_LIBS = \
    -lfluidThermophysicalModels \
    -lreactionThermophysicalModels \
    -lspecie \
-    -llaminarFlameSpeedModels \
-    -lfiniteVolume
+    -llaminarFlameSpeedModels
--- a/applications/solvers/combustion/engineFoam/UEqn.H
+++ b/applications/solvers/combustion/engineFoam/UEqn.H
@ -3,6 +3,8 @@
        fvm::ddt(rho, U)
      + fvm::div(phi, U)
      + turbulence->divDevRhoReff(U)
+     ==
+        sources(rho, U)
    );

    if (pimple.momentumPredictor())
--- a/applications/solvers/combustion/engineFoam/engineFoam.C
+++ b/applications/solvers/combustion/engineFoam/engineFoam.C
@ -59,6 +59,7 @@ Description
 #include "OFstream.H"
 #include "mathematicalConstants.H"
 #include "pimpleControl.H"
+#include "IObasicSourceList.H"

 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

--- a/applications/solvers/combustion/engineFoam/pEqn.H
+++ b/applications/solvers/combustion/engineFoam/pEqn.H
@ -13,6 +13,8 @@ if (pimple.transonic())
       *((fvc::interpolate(HbyA) & mesh.Sf()) - fvc::meshPhi(rho, U))
    );

+    sources.relativeFlux(fvc::interpolate(psi), phid);
+
    while (pimple.correctNonOrthogonal())
    {
        fvScalarMatrix pEqn
@ -20,8 +22,12 @@ if (pimple.transonic())
            fvm::ddt(psi, p)
          + fvm::div(phid, p)
          - fvm::laplacian(rho*rAU, p)
+         ==
+            sources(psi, p, rho.name())
        );

+        sources.constrain(pEqn, rho.name());
+
        pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter())));

        if (pimple.finalNonOrthogonalIter())
@ -39,6 +45,8 @@ else
       *((fvc::interpolate(HbyA) & mesh.Sf()) - fvc::meshPhi(rho, U))
    );

+    sources.relativeFlux(fvc::interpolate(rho), phiHbyA);
+
    while (pimple.correctNonOrthogonal())
    {
        fvScalarMatrix pEqn
@ -46,8 +54,12 @@ else
            fvm::ddt(psi, p)
          + fvc::div(phiHbyA)
          - fvm::laplacian(rho*rAU, p)
+         ==
+            sources(psi, p, rho.name())
        );

+        sources.constrain(pEqn, rho.name());
+
        pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter())));

        if (pimple.finalNonOrthogonalIter())
@ -62,6 +74,7 @@ else

 U = HbyA - rAU*fvc::grad(p);
 U.correctBoundaryConditions();
+sources.correct(U);
 K = 0.5*magSqr(U);

 if (thermo.dpdt())
--- a/applications/solvers/combustion/fireFoam/Make/options
+++ b/applications/solvers/combustion/fireFoam/Make/options
@ -1,6 +1,8 @@
 EXE_INC = \
    -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/fieldSources/lnInclude \
    -I${LIB_SRC}/meshTools/lnInclude \
+    -I${LIB_SRC}/sampling/lnInclude \
    -I$(LIB_SRC)/turbulenceModels/compressible/turbulenceModel \
    -I$(LIB_SRC)/lagrangian/distributionModels/lnInclude \
    -I$(LIB_SRC)/thermophysicalModels/specie/lnInclude \
@ -29,7 +31,9 @@ EXE_INC = \

 EXE_LIBS = \
    -lfiniteVolume \
+    -lfieldSources \
    -lmeshTools \
+    -lsampling \
    -lcompressibleRASModels \
    -lcompressibleLESModels \
    -lspecie \
--- a/applications/solvers/combustion/fireFoam/UEqn.H
+++ b/applications/solvers/combustion/fireFoam/UEqn.H
@ -5,10 +5,13 @@
      + turbulence->divDevRhoReff(U)
     ==
        parcels.SU(U)
+      + sources(rho, U)
    );

    UEqn.relax();

+    sources.constrain(UEqn);
+
    if (pimple.momentumPredictor())
    {
        solve
--- a/applications/solvers/combustion/fireFoam/YEEqn.H
+++ b/applications/solvers/combustion/fireFoam/YEEqn.H
@ -30,9 +30,13 @@ tmp<fv::convectionScheme<scalar> > mvConvection
                parcels.SYi(i, Yi)
              + surfaceFilm.Srho(i)
              + combustion->R(Yi)
+              + sources(rho, Yi)
            );

            YiEqn.relax();
+
+            sources.constrain(YiEqn);
+
            YiEqn.solve(mesh.solver("Yi"));

            Yi.max(0.0);
@ -69,9 +73,13 @@ tmp<fv::convectionScheme<scalar> > mvConvection
      + radiation->Sh(thermo)
      + parcels.Sh(he)
      + surfaceFilm.Sh()
+      + sources(rho, he)
    );

    EEqn.relax();
+
+    sources.constrain(EEqn);
+
    EEqn.solve();

    thermo.correct();
--- a/applications/solvers/combustion/fireFoam/createFields.H
+++ b/applications/solvers/combustion/fireFoam/createFields.H
@ -146,3 +146,6 @@
    (
        additionalControlsDict.lookup("solvePrimaryRegion")
    );
+
+    IObasicSourceList sources(mesh);
+
--- a/applications/solvers/combustion/fireFoam/fireFoam.C
+++ b/applications/solvers/combustion/fireFoam/fireFoam.C
@ -40,6 +40,7 @@ Description
 #include "solidChemistryModel.H"
 #include "psiCombustionModel.H"
 #include "pimpleControl.H"
+#include "IObasicSourceList.H"

 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

--- a/applications/solvers/combustion/fireFoam/pEqn.H
+++ b/applications/solvers/combustion/fireFoam/pEqn.H
@ -21,6 +21,7 @@ surfaceScalarField phiHbyA
  + phig
 );

+sources.relativeFlux(phiHbyA);

 while (pimple.correctNonOrthogonal())
 {
@ -33,8 +34,11 @@ while (pimple.correctNonOrthogonal())
     ==
        parcels.Srho()
      + surfaceFilm.Srho()
+      + sources(psi, p_rgh, rho.name())
    );

+    sources.constrain(p_rghEqn, rho.name());
+
    p_rghEqn.solve(mesh.solver(p_rgh.select(pimple.finalInnerIter())));

    if (pimple.finalNonOrthogonalIter())
@ -42,6 +46,7 @@ while (pimple.correctNonOrthogonal())
        phi = phiHbyA + p_rghEqn.flux();
        U = HbyA + rAU*fvc::reconstruct((p_rghEqn.flux() + phig)/rhorAUf);
        U.correctBoundaryConditions();
+        sources.correct(U);
    }
 }

--- a/applications/solvers/combustion/fireFoam/rhoEqn.H
+++ b/applications/solvers/combustion/fireFoam/rhoEqn.H
@ -30,14 +30,19 @@ Description
 \*---------------------------------------------------------------------------*/

 {
-    solve
+    fvScalarMatrix rhoEqn
    (
        fvm::ddt(rho)
      + fvc::div(phi)
      ==
        parcels.Srho(rho)
      + surfaceFilm.Srho()
+      + sources(rho)
    );
+
+    sources.constrain(rhoEqn);
+
+    rhoEqn.solve();
 }

 // ************************************************************************* //
--- a/applications/solvers/combustion/reactingFoam/EEqn.H
+++ b/applications/solvers/combustion/reactingFoam/EEqn.H
@ -19,9 +19,13 @@
 //    - fvm::laplacian(turbulence->muEff(), he)  // unit lewis no.
     ==
        reaction->Sh()
+      + sources(rho, he)
    );

    EEqn.relax();
+
+    sources.constrain(EEqn);
+
    EEqn.solve();

    thermo.correct();
--- a/applications/solvers/combustion/reactingFoam/Make/options
+++ b/applications/solvers/combustion/reactingFoam/Make/options
@ -1,14 +1,21 @@
 EXE_INC = \
+    -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/fieldSources/lnInclude \
+    -I$(LIB_SRC)/meshTools/lnInclude \
+    -I$(LIB_SRC)/sampling/lnInclude \
    -I$(LIB_SRC)/turbulenceModels/compressible/turbulenceModel \
    -I$(LIB_SRC)/thermophysicalModels/specie/lnInclude \
    -I$(LIB_SRC)/thermophysicalModels/reactionThermo/lnInclude \
    -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
    -I$(LIB_SRC)/thermophysicalModels/chemistryModel/lnInclude \
    -I$(LIB_SRC)/ODE/lnInclude \
-    -I$(LIB_SRC)/finiteVolume/lnInclude \
    -I$(LIB_SRC)/combustionModels/lnInclude

 EXE_LIBS = \
+    -lfiniteVolume \
+    -lfieldSources \
+    -lmeshTools \
+    -lsampling \
    -lcompressibleTurbulenceModel \
    -lcompressibleRASModels \
    -lcompressibleLESModels \
@ -17,5 +24,4 @@ EXE_LIBS = \
    -lfluidThermophysicalModels \
    -lchemistryModel \
    -lODE \
-    -lfiniteVolume \
    -lcombustionModels
--- a/applications/solvers/combustion/reactingFoam/UEqn.H
+++ b/applications/solvers/combustion/reactingFoam/UEqn.H
@ -5,10 +5,13 @@
      + turbulence->divDevRhoReff(U)
     ==
        rho*g
+      + sources(rho, U)
    );

    UEqn.relax();

+    sources.constrain(UEqn);
+
    if (pimple.momentumPredictor())
    {
        solve(UEqn == -fvc::grad(p));
--- a/applications/solvers/combustion/reactingFoam/YEqn.H
+++ b/applications/solvers/combustion/reactingFoam/YEqn.H
@ -28,9 +28,13 @@ tmp<fv::convectionScheme<scalar> > mvConvection
              - fvm::laplacian(turbulence->muEff(), Yi)
             ==
                reaction->R(Yi)
+              + sources(rho, Yi)
            );

            YiEqn.relax();
+
+            sources.constrain(YiEqn);
+
            YiEqn.solve(mesh.solver("Yi"));

            Yi.max(0.0);
--- a/applications/solvers/combustion/reactingFoam/createFields.H
+++ b/applications/solvers/combustion/reactingFoam/createFields.H
@ -99,3 +99,5 @@ volScalarField dQ
    mesh,
    dimensionedScalar("dQ", dimEnergy/dimTime, 0.0)
 );
+
+IObasicSourceList sources(mesh);
--- a/applications/solvers/combustion/reactingFoam/pEqn.H
+++ b/applications/solvers/combustion/reactingFoam/pEqn.H
@ -16,6 +16,8 @@ if (pimple.transonic())
        )
    );

+    sources.relativeFlux(fvc::interpolate(psi), phid);
+
    while (pimple.correctNonOrthogonal())
    {
        fvScalarMatrix pEqn
@ -23,8 +25,12 @@ if (pimple.transonic())
            fvm::ddt(psi, p)
          + fvm::div(phid, p)
          - fvm::laplacian(rho*rAU, p)
+         ==
+            sources(psi, p, rho.name())
        );

+        sources.constrain(pEqn, rho.name());
+
        pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter())));

        if (pimple.finalNonOrthogonalIter())
@ -45,6 +51,8 @@ else
        )
    );

+    sources.relativeFlux(fvc::interpolate(rho), phiHbyA);
+
    while (pimple.correctNonOrthogonal())
    {
        fvScalarMatrix pEqn
@ -52,8 +60,12 @@ else
            fvm::ddt(psi, p)
          + fvc::div(phiHbyA)
          - fvm::laplacian(rho*rAU, p)
+         ==
+            sources(psi, p, rho.name())
        );

+        sources.constrain(pEqn, rho.name());
+
        pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter())));

        if (pimple.finalNonOrthogonalIter())
@ -68,6 +80,7 @@ else

 U = HbyA - rAU*fvc::grad(p);
 U.correctBoundaryConditions();
+sources.correct(U);
 K = 0.5*magSqr(U);

 if (thermo.dpdt())
--- a/applications/solvers/combustion/reactingFoam/reactingFoam.C
+++ b/applications/solvers/combustion/reactingFoam/reactingFoam.C
@ -34,6 +34,7 @@ Description
 #include "psiCombustionModel.H"
 #include "multivariateScheme.H"
 #include "pimpleControl.H"
+#include "IObasicSourceList.H"

 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

--- a/applications/solvers/combustion/rhoReactingFoam/Make/options
+++ b/applications/solvers/combustion/rhoReactingFoam/Make/options
@ -1,17 +1,23 @@
 EXE_INC = \
-    -I../reactingFoam \
+    -I$(FOAM_SOLVERS)/combustion/reactingFoam \
+    -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/fieldSources/lnInclude \
+    -I$(LIB_SRC)/meshTools/lnInclude \
+    -I$(LIB_SRC)/sampling/lnInclude \
    -I$(LIB_SRC)/turbulenceModels/compressible/turbulenceModel \
    -I$(LIB_SRC)/thermophysicalModels/specie/lnInclude \
    -I$(LIB_SRC)/thermophysicalModels/reactionThermo/lnInclude \
    -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
    -I$(LIB_SRC)/thermophysicalModels/chemistryModel/lnInclude \
    -I$(LIB_SRC)/ODE/lnInclude \
-    -I$(LIB_SRC)/finiteVolume/lnInclude \
-    -I$(FOAM_SOLVERS)/combustion/reactingFoam \
    -I$(LIB_SRC)/combustionModels/lnInclude


 EXE_LIBS = \
+    -lfiniteVolume \
+    -lfieldSources \
+    -lmeshTools \
+    -lsampling \
    -lcompressibleTurbulenceModel \
    -lcompressibleRASModels \
    -lcompressibleLESModels \
@ -20,5 +26,4 @@ EXE_LIBS = \
    -lfluidThermophysicalModels \
    -lchemistryModel \
    -lODE \
-    -lfiniteVolume \
    -lcombustionModels
--- a/applications/solvers/combustion/rhoReactingFoam/createFields.H
+++ b/applications/solvers/combustion/rhoReactingFoam/createFields.H
@ -101,3 +101,5 @@ volScalarField dQ
    mesh,
    dimensionedScalar("dQ", dimEnergy/dimTime, 0.0)
 );
+
+IObasicSourceList sources(mesh);
--- a/applications/solvers/combustion/rhoReactingFoam/pEqn.H
+++ b/applications/solvers/combustion/rhoReactingFoam/pEqn.H
@ -18,6 +18,8 @@
          + fvc::ddtPhiCorr(rAU, rho, U, phi)
        );

+        sources.relativeFlux(phiHbyA);
+
        surfaceScalarField phid("phid", fvc::interpolate(thermo.psi())*phiHbyA);

        phiHbyA *= fvc::interpolate(rho);
@ -34,8 +36,12 @@
            (
                pDDtEqn
              - fvm::laplacian(rho*rAU, p)
+             ==
+                sources(psi, p, rho.name())
            );

+            sources.constrain(pEqn, rho.name());
+
            pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter())));

            if (pimple.finalNonOrthogonalIter())
@ -56,10 +62,14 @@
            )
        );

+        sources.relativeFlux(fvc::interpolate(rho), phiHbyA);
+
        fvScalarMatrix pDDtEqn
        (
            fvc::ddt(rho) + psi*correction(fvm::ddt(p))
          + fvc::div(phiHbyA)
+         ==
+            sources(psi, p, rho.name())
        );

        while (pimple.correctNonOrthogonal())
@ -70,6 +80,8 @@
              - fvm::laplacian(rho*rAU, p)
            );

+            sources.constrain(pEqn, rho.name());
+
            pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter())));

            if (pimple.finalNonOrthogonalIter())
@ -87,6 +99,7 @@

    U = HbyA - rAU*fvc::grad(p);
    U.correctBoundaryConditions();
+    sources.correct(U);
    K = 0.5*magSqr(U);

    if (thermo.dpdt())
--- a/applications/solvers/combustion/rhoReactingFoam/rhoReactingFoam.C
+++ b/applications/solvers/combustion/rhoReactingFoam/rhoReactingFoam.C
@ -35,6 +35,7 @@ Description
 #include "turbulenceModel.H"
 #include "multivariateScheme.H"
 #include "pimpleControl.H"
+#include "IObasicSourceList.H"

 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //