fvModels: Improved interface for mass/volume sources

The interface for fvModels has been modified to improve its application
to "proxy" equations. That is, equations that are not straightforward
statements of conservation laws in OpenFOAM's usual convention.

A standard conservation law typically takes the following form:

    fvMatrix<scalar> psiEqn
    (
        fvm::ddt(alpha, rho, psi)
      + <fluxes>
     ==
        <sources>
    );

A proxy equation, on the other hand, may be a derivation or
rearrangement of a law like this, and may be linearised in terms of a
different variable.

The pressure equation is the most common example of a proxy equation. It
represents a statement of the conservation of volume or mass, but it is
a rearrangement of the original continuity equation, and it has been
linearised in terms of a different variable; the pressure. Another
example is that in the pre-predictor of a VoF solver the
phase-continuity equation is constructed, but it is linearised in terms
of volume fraction rather than density.

In these situations, fvModels sources are now applied by calling:

    fvModels().sourceProxy(<conserved-fields ...>, <equation-field>)

Where <conserved-fields ...> are (alpha, rho, psi), (rho, psi), just
(psi), or are omitted entirely (for volume continuity), and the
<equation-field> is the field associated with the proxy equation. This
produces a source term identical in value to the following call:

    fvModels().source(<conserved-fields ...>)

It is only the linearisation in terms of <equation-field> that differs
between these two calls.

This change permits much greater flexibility in the handling of mass and
volume sources than the previous name-based system did. All the relevant
fields are available, dimensions can be used in the logic to determine
what sources are being constructed, and sources relating to a given
conservation law all share the same function.

This commit adds the functionality for injection-type sources in the
compressibleVoF solver. A following commit will add a volume source
model for use in incompressible solvers.

applications/legacy/combustion/PDRFoam/pEqn.H

@@ -24,7 +24,7 @@ if (pimple.transonic())
           + fvm::div(phid, p)
           - fvm::laplacian(rho*invA, p)
          ==
-            betav*fvModels.source(psi, p, rho.name())
+            betav*fvModels.sourceProxy(rho, p)
         );
 
         pEqn.solve();
@@ -54,7 +54,7 @@ else
           + fvc::div(phiHbyA)
           - fvm::laplacian(rho*invA, p)
          ==
-            betav*fvModels.source(psi, p, rho.name())
+            betav*fvModels.sourceProxy(rho, p)
         );
 
         pEqn.solve();

applications/legacy/compressible/rhoPorousSimpleFoam/pEqn.H

@@ -34,7 +34,7 @@
             tpEqn =
             (
                 fvm::laplacian(rho*trTU(), p)
-              + fvModels.source(psi, p, rho.name())
+              + fvModels.sourceProxy(rho, p)
              ==
                 fvc::div(phiHbyA)
             );
@@ -44,7 +44,7 @@
             tpEqn =
             (
                 fvm::laplacian(rho*trAU(), p)
-              + fvModels.source(psi, p, rho.name())
+              + fvModels.sourceProxy(rho, p)
              ==
                 fvc::div(phiHbyA)
             );

applications/modules/compressibleMultiphaseVoF/pressureCorrector.C

@@ -88,7 +88,7 @@ void Foam::solvers::compressibleMultiphaseVoF::pressureCorrector()
                 (
                     fvc::ddt(rho) + thermo.psi()*correction(fvm::ddt(p_rgh))
                   + fvc::div(phi, rho) - fvc::Sp(fvc::div(phi), rho)
-                  - (fvModels().source(phase, rho)&rho)
+                  - fvModels().sourceProxy(phase, rho, p_rgh)
                 ).ptr()
             );
         }

applications/modules/compressibleVoF/alphaSuSp.C

@@ -30,46 +30,52 @@ License
 
 void Foam::solvers::compressibleVoF::alphaSuSp
 (
-    tmp<volScalarField::Internal>& Su,
+    tmp<volScalarField::Internal>& tSu,
-    tmp<volScalarField::Internal>& Sp
+    tmp<volScalarField::Internal>& tSp
 )
 {
-    Sp = volScalarField::Internal::New
+    const dimensionedScalar Szero(dgdt.dimensions(), 0);
-    (
-        "Sp",
-        mesh,
-        dimensionedScalar(dgdt.dimensions(), 0)
-    );
 
-    Su = volScalarField::Internal::New
+    tSp = volScalarField::Internal::New("Sp", mesh, Szero);
-    (
+    tSu = volScalarField::Internal::New("Su", mesh, Szero);
-        "Su",
-        mesh,
-        dimensionedScalar(dgdt.dimensions(), 0)
-    );
 
-    if (fvModels().addsSupToField(alpha1.name()))
+    volScalarField::Internal& Sp = tSp.ref();
+    volScalarField::Internal& Su = tSu.ref();
+
+    if (fvModels().addsSupToField(mixture.rho1().name()))
     {
-        // Phase change alpha1 source
+        const volScalarField::Internal alpha2ByRho1(alpha2()/mixture.rho1()());
-        const fvScalarMatrix alphaSup(fvModels().source(alpha1));
+        const fvScalarMatrix alphaRho1Sup
+        (
+            fvModels().sourceProxy(alpha1, mixture.rho1(), alpha1)
+        );
 
-        Su = alphaSup.Su();
+        Su += alpha2ByRho1*alphaRho1Sup.Su();
-        Sp = alphaSup.Sp();
+        Sp += alpha2ByRho1*alphaRho1Sup.Sp();
     }
 
-    volScalarField::Internal& SpRef = Sp.ref();
+    if (fvModels().addsSupToField(mixture.rho2().name()))
-    volScalarField::Internal& SuRef = Su.ref();
+    {
+        const volScalarField::Internal alpha1ByRho2(alpha1()/mixture.rho2()());
+        const fvScalarMatrix alphaRho2Sup
+        (
+            fvModels().sourceProxy(alpha2, mixture.rho2(), alpha2)
+        );
+
+        Su -= alpha1ByRho2*(alphaRho2Sup.Su() + alphaRho2Sup.Sp());
+        Sp += alpha1ByRho2*alphaRho2Sup.Sp();
+    }
 
     forAll(dgdt, celli)
     {
         if (dgdt[celli] > 0.0)
         {
-            SpRef[celli] -= dgdt[celli]/max(1.0 - alpha1[celli], 1e-4);
+            Sp[celli] -= dgdt[celli]/max(1.0 - alpha1[celli], 1e-4);
-            SuRef[celli] += dgdt[celli]/max(1.0 - alpha1[celli], 1e-4);
+            Su[celli] += dgdt[celli]/max(1.0 - alpha1[celli], 1e-4);
         }
         else if (dgdt[celli] < 0.0)
         {
-            SpRef[celli] += dgdt[celli]/max(alpha1[celli], 1e-4);
+            Sp[celli] += dgdt[celli]/max(alpha1[celli], 1e-4);
         }
     }
 }

applications/modules/compressibleVoF/compressibleVoF.H

@@ -156,7 +156,8 @@ protected:
         {
             return
                !incompressible()
-             || fvModels().addsSupToField(alpha1.name());
+             || fvModels().addsSupToField(alpha1.name())
+             || fvModels().addsSupToField(alpha2.name());
         }
 
         //- Return the mixture compressibility/density

applications/modules/compressibleVoF/fvModels/VoFCavitation/VoFCavitation.C

@@ -69,9 +69,7 @@ Foam::fv::compressible::VoFCavitation::VoFCavitation
         )
     ),
 
-    cavitation_(Foam::compressible::cavitationModel::New(dict, mixture_)),
+    cavitation_(Foam::compressible::cavitationModel::New(dict, mixture_))
-
-    alphaName_(mixture_.alpha1().name())
 {}
 
 
@@ -79,14 +77,19 @@ Foam::fv::compressible::VoFCavitation::VoFCavitation
 
 Foam::wordList Foam::fv::compressible::VoFCavitation::addSupFields() const
 {
-    return {alphaName_, "p_rgh"};
+    return
+    {
+        mixture_.rho1().name(),
+        mixture_.rho2().name()
+    };
 }
 
 
 void Foam::fv::compressible::VoFCavitation::addSup
 (
-    fvMatrix<scalar>& eqn,
+    const volScalarField& alpha,
-    const word& fieldName
+    const volScalarField& rho,
+    fvMatrix<scalar>& eqn
 ) const
 {
     if (debug)
@@ -94,64 +97,68 @@ void Foam::fv::compressible::VoFCavitation::addSup
         Info<< type() << ": applying source to " << eqn.psi().name() << endl;
     }
 
-    if (fieldName == alphaName_)
+    if (&rho == &mixture_.rho1() || &rho == &mixture_.rho2())
     {
-        const volScalarField::Internal alpha1Coeff
+        const volScalarField& alpha1 = mixture_.alpha1();
-        (
+        const volScalarField& alpha2 = mixture_.alpha2();
-            1.0/mixture_.rho1()()
-          - mixture_.alpha1()()*(1.0/mixture_.rho1()() - 1.0/mixture_.rho2()())
-        );
 
-        const Pair<tmp<volScalarField::Internal>> mDot12Alpha
+        const scalar s = &alpha == &alpha1 ? +1 : -1;
-        (
-            cavitation_->mDot12Alpha()
-        );
 
-        const volScalarField::Internal vDot1Alpha(alpha1Coeff*mDot12Alpha[0]());
+        // Volume-fraction linearisation
-        const volScalarField::Internal vDot2Alpha(alpha1Coeff*mDot12Alpha[1]());
+        if (&alpha == &eqn.psi())
+        {
+            const Pair<tmp<volScalarField::Internal>> mDot12Alpha
+            (
+                cavitation_->mDot12Alpha()
+            );
+            const volScalarField::Internal& mDot1Alpha2 = mDot12Alpha[0]();
+            const volScalarField::Internal& mDot2Alpha1 = mDot12Alpha[1]();
 
-        eqn += fvm::Sp(-vDot2Alpha - vDot1Alpha, eqn.psi()) + vDot1Alpha;
+            eqn +=
+                (&alpha == &alpha1 ? mDot1Alpha2 : mDot2Alpha1)
+              - fvm::Sp(mDot1Alpha2 + mDot2Alpha1, eqn.psi());
+        }
+
+        // Pressure linearisation
+        else if (eqn.psi().name() == "p_rgh")
+        {
+            const Pair<tmp<volScalarField::Internal>> mDot12P
+            (
+                cavitation_->mDot12P()
+            );
+            const volScalarField::Internal& mDot1P = mDot12P[0];
+            const volScalarField::Internal& mDot2P = mDot12P[1];
+
+            const volScalarField::Internal& rho =
+                mesh().lookupObject<volScalarField>("rho");
+            const volScalarField::Internal& gh =
+                mesh().lookupObject<volScalarField>("gh");
+
+            eqn +=
+                fvm::Sp(s*(mDot1P - mDot2P), eqn.psi())
+              + s*(mDot1P - mDot2P)*rho*gh
+              - s*(mDot1P*cavitation_->pSat1() - mDot2P*cavitation_->pSat2());
+        }
+
+        // Explicit non-linearised value. Used in density predictors and
+        // continuity error terms.
+        else
+        {
+            const Pair<tmp<volScalarField::Internal>> mDot12Alpha
+            (
+                cavitation_->mDot12Alpha()
+            );
+            const volScalarField::Internal mDot1(mDot12Alpha[0]*alpha2);
+            const volScalarField::Internal mDot2(mDot12Alpha[1]*alpha1);
+
+            eqn += s*(mDot1 - mDot2);
+        }
     }
-}
+    else
-
-
-void Foam::fv::compressible::VoFCavitation::addSup
-(
-    const volScalarField&,
-    fvMatrix<scalar>& eqn,
-    const word& fieldName
-) const
-{
-    if (debug)
     {
-        Info<< type() << ": applying source to " << eqn.psi().name() << endl;
+        FatalErrorInFunction
-    }
+            << "Support for field " << alpha.name() << " is not implemented"
-
+            << exit(FatalError);
-    if (fieldName == "p_rgh")
-    {
-        const volScalarField::Internal& rho =
-            mesh().lookupObject<volScalarField>("rho");
-
-        const volScalarField::Internal& gh =
-            mesh().lookupObject<volScalarField>("gh");
-
-        const volScalarField::Internal pCoeff
-        (
-            1.0/mixture_.rho1()() - 1.0/mixture_.rho2()()
-        );
-
-        const Pair<tmp<volScalarField::Internal>> mDot12P
-        (
-            cavitation_->mDot12P()
-        );
-
-        const volScalarField::Internal vDot1P(pCoeff*mDot12P[0]);
-        const volScalarField::Internal vDot2P(pCoeff*mDot12P[1]);
-
-        eqn +=
-            vDot2P*cavitation_->pSat1() - vDot1P*cavitation_->pSat2()
-          - (vDot2P - vDot1P)*rho*gh
-          - fvm::Sp(vDot2P - vDot1P, eqn.psi());
     }
 }
 

applications/modules/compressibleVoF/fvModels/VoFCavitation/VoFCavitation.H

@@ -104,11 +104,9 @@ class VoFCavitation
         //- Reference to the mixture properties
         const compressibleTwoPhaseVoFMixture& mixture_;
 
+        //- The cavitation model
         autoPtr<Foam::compressible::cavitationModel> cavitation_;
 
-        //- The name of the VoF phase-fraction
-        word alphaName_;
-
 
 public:
 
@@ -140,19 +138,12 @@ public:
         //  to the transport equation
         virtual wordList addSupFields() const;
 
-        //- Add implicit/explicit contributions to VoF phase-fraction equation
+        //- Add a source to the phase continuity equation
         virtual void addSup
         (
-            fvMatrix<scalar>& eqn,
+            const volScalarField& alpha,
-            const word& fieldName
+            const volScalarField& rho,
-        ) const;
+            fvMatrix<scalar>& eqn
-
-        //- Add implicit/explicit contributions to p_rgh equation
-        virtual void addSup
-        (
-            const volScalarField& psi,
-            fvMatrix<scalar>& eqn,
-            const word& fieldName
         ) const;
 
 

applications/modules/compressibleVoF/fvModels/VoFClouds/VoFClouds.C

@@ -109,8 +109,8 @@ void Foam::fv::VoFClouds::correct()
 void Foam::fv::VoFClouds::addSup
 (
     const volScalarField& alpha,
-    fvMatrix<scalar>& eqn,
+    const volScalarField& rho,
-    const word& fieldName
+    fvMatrix<scalar>& eqn
 ) const
 {
     if (debug)
@@ -118,14 +118,14 @@ void Foam::fv::VoFClouds::addSup
         Info<< type() << ": applying source to " << eqn.psi().name() << endl;
     }
 
-    if (fieldName == thermo_.rho()().name())
+    if (&rho == &thermo_.rho()())
     {
         eqn += clouds_.Srho();
     }
     else
     {
         FatalErrorInFunction
-            << "Support for field " << fieldName << " is not implemented"
+            << "Support for field " << rho.name() << " is not implemented"
             << exit(FatalError);
     }
 }
@@ -135,8 +135,8 @@ void Foam::fv::VoFClouds::addSup
 (
     const volScalarField& alpha,
     const volScalarField& rho,
-    fvMatrix<scalar>& eqn,
+    const volScalarField& he,
-    const word& fieldName
+    fvMatrix<scalar>& eqn
 ) const
 {
     if (debug)
@@ -144,14 +144,14 @@ void Foam::fv::VoFClouds::addSup
         Info<< type() << ": applying source to " << eqn.psi().name() << endl;
     }
 
-    if (fieldName == thermo_.he().name())
+    if (&he == &thermo_.he())
     {
         eqn += clouds_.Sh(eqn.psi());
     }
     else
     {
         FatalErrorInFunction
-            << "Support for field " << fieldName << " is not implemented"
+            << "Support for field " << he.name() << " is not implemented"
             << exit(FatalError);
     }
 }
@@ -160,8 +160,8 @@ void Foam::fv::VoFClouds::addSup
 void Foam::fv::VoFClouds::addSup
 (
     const volScalarField& rho,
-    fvMatrix<vector>& eqn,
+    const volVectorField& U,
-    const word& fieldName
+    fvMatrix<vector>& eqn
 ) const
 {
     if (debug)
@@ -169,7 +169,16 @@ void Foam::fv::VoFClouds::addSup
         Info<< type() << ": applying source to " << eqn.psi().name() << endl;
     }
 
-    eqn += clouds_.SU(eqn.psi());
+    if (U.name() == "U")
+    {
+        eqn += clouds_.SU(eqn.psi());
+    }
+    else
+    {
+        FatalErrorInFunction
+            << "Support for field " << U.name() << " is not implemented"
+            << exit(FatalError);
+    }
 }
 
 

applications/modules/compressibleVoF/fvModels/VoFClouds/VoFClouds.H

@@ -130,8 +130,8 @@ public:
             virtual void addSup
             (
                 const volScalarField& alpha,
-                fvMatrix<scalar>& eqn,
+                const volScalarField& rho,
-                const word& fieldName
+                fvMatrix<scalar>& eqn
             ) const;
 
             //- Add explicit contribution to phase energy equation
@@ -139,16 +139,16 @@ public:
             (
                 const volScalarField& alpha,
                 const volScalarField& rho,
-                fvMatrix<scalar>& eqn,
+                const volScalarField& he,
-                const word& fieldName
+                fvMatrix<scalar>& eqn
             ) const;
 
             //- Add implicit contribution to mixture momentum equation
             virtual void addSup
             (
                 const volScalarField& rho,
-                fvMatrix<vector>& eqn,
+                const volVectorField& U,
-                const word& fieldName
+                fvMatrix<vector>& eqn
             ) const;
 
 

applications/modules/compressibleVoF/fvModels/VoFSolidificationMeltingSource/VoFSolidificationMeltingSource.C

@@ -132,8 +132,8 @@ void Foam::fv::VoFSolidificationMeltingSource::addSup
 (
     const volScalarField& alpha,
     const volScalarField& rho,
-    fvMatrix<scalar>& eqn,
+    const volScalarField& he,
-    const word& fieldName
+    fvMatrix<scalar>& eqn
 ) const
 {
     if (debug)
@@ -148,8 +148,8 @@ void Foam::fv::VoFSolidificationMeltingSource::addSup
 void Foam::fv::VoFSolidificationMeltingSource::addSup
 (
     const volScalarField& rho,
-    fvMatrix<vector>& eqn,
+    const volVectorField& U,
-    const word& fieldName
+    fvMatrix<vector>& eqn
 ) const
 {
     if (debug)

applications/modules/compressibleVoF/fvModels/VoFSolidificationMeltingSource/VoFSolidificationMeltingSource.H

@@ -182,16 +182,16 @@ public:
             (
                 const volScalarField& alpha,
                 const volScalarField& rho,
-                fvMatrix<scalar>& eqn,
+                const volScalarField& he,
-                const word& fieldName
+                fvMatrix<scalar>& eqn
             ) const;
 
             //- Add implicit contribution to mixture momentum equation
             virtual void addSup
             (
                 const volScalarField& rho,
-                fvMatrix<vector>& eqn,
+                const volVectorField& U,
-                const word& fieldName
+                fvMatrix<vector>& eqn
             ) const;
 
 

applications/modules/compressibleVoF/fvModels/VoFTurbulenceDamping/VoFTurbulenceDamping.C

@@ -127,8 +127,8 @@ Foam::fv::compressible::VoFTurbulenceDamping::addSupFields() const
 void Foam::fv::compressible::VoFTurbulenceDamping::addSup
 (
     const volScalarField& rho,
-    fvMatrix<scalar>& eqn,
+    const volScalarField& field,
-    const word& fieldName
+    fvMatrix<scalar>& eqn
 ) const
 {
     if (debug)
@@ -142,12 +142,12 @@ void Foam::fv::compressible::VoFTurbulenceDamping::addSup
       + mixture_.alpha2()()*mixture_.rho2()()*sqr(mixture_.thermo2().nu()()())
     );
 
-    if (fieldName == "epsilon")
+    if (field.name() == "epsilon")
     {
         eqn += mixture_.interfaceFraction()
            *C2_*aRhoSqrnu*momentumTransport_.k()()/pow4(delta_);
     }
-    else if (fieldName == "omega")
+    else if (field.name() == "omega")
     {
         eqn += mixture_.interfaceFraction()
            *beta_*aRhoSqrnu/(sqr(betaStar_)*pow4(delta_));
@@ -155,7 +155,7 @@ void Foam::fv::compressible::VoFTurbulenceDamping::addSup
     else
     {
         FatalErrorInFunction
-            << "Support for field " << fieldName << " is not implemented"
+            << "Support for field " << field.name() << " is not implemented"
             << exit(FatalError);
     }
 }
@@ -164,9 +164,9 @@ void Foam::fv::compressible::VoFTurbulenceDamping::addSup
 void Foam::fv::compressible::VoFTurbulenceDamping::addSup
 (
     const volScalarField& alpha,
-    const volScalarField&,
+    const volScalarField& rho,
-    fvMatrix<scalar>& eqn,
+    const volScalarField& field,
-    const word& fieldName
+    fvMatrix<scalar>& eqn
 ) const
 {
     if (debug)
@@ -193,12 +193,12 @@ void Foam::fv::compressible::VoFTurbulenceDamping::addSup
             << exit(FatalError);
     }
 
-    if (fieldName == IOobject::groupName("epsilon", phaseName_))
+    if (field.name() == IOobject::groupName("epsilon", phaseName_))
     {
         eqn += mixture_.interfaceFraction()
            *C2_*taRhoSqrnu*momentumTransport_.k()()/pow4(delta_);
     }
-    else if (fieldName == IOobject::groupName("omega", phaseName_))
+    else if (field.name() == IOobject::groupName("omega", phaseName_))
     {
         eqn += mixture_.interfaceFraction()
            *beta_*taRhoSqrnu/(sqr(betaStar_)*pow4(delta_));
@@ -206,7 +206,7 @@ void Foam::fv::compressible::VoFTurbulenceDamping::addSup
     else
     {
         FatalErrorInFunction
-            << "Support for field " << fieldName << " is not implemented"
+            << "Support for field " << field.name() << " is not implemented"
             << exit(FatalError);
     }
 }

applications/modules/compressibleVoF/fvModels/VoFTurbulenceDamping/VoFTurbulenceDamping.H

@@ -146,26 +146,31 @@ public:
 
     // Member Functions
 
-        //- Return the list of fields for which the option adds source term
+        // Checks
-        //  to the transport equation
-        virtual wordList addSupFields() const;
 
-        //- Add explicit contribution to epsilon or omega equation
+            //- Return the list of fields for which the option adds source term
-        virtual void addSup
+            //  to the transport equation
-        (
+            virtual wordList addSupFields() const;
-            const volScalarField& rho,
-            fvMatrix<scalar>& eqn,
-            const word& fieldName
-        ) const;
 
-        //- Add explicit contribution to phase epsilon or omega equation
+
-        virtual void addSup
+        // Sources
-        (
+
-            const volScalarField& alpha,
+            //- Add explicit contribution to epsilon or omega equation
-            const volScalarField& rho,
+            virtual void addSup
-            fvMatrix<scalar>& eqn,
+            (
-            const word& fieldName
+                const volScalarField& rho,
-        ) const;
+                const volScalarField& field,
+                fvMatrix<scalar>& eqn
+            ) const;
+
+            //- Add explicit contribution to phase epsilon or omega equation
+            virtual void addSup
+            (
+                const volScalarField& alpha,
+                const volScalarField& rho,
+                const volScalarField& field,
+                fvMatrix<scalar>& eqn
+            ) const;
 
 
         // Mesh changes

applications/modules/compressibleVoF/pressureCorrector.C

@@ -81,19 +81,11 @@ void Foam::solvers::compressibleVoF::pressureCorrector()
         // Update the pressure BCs to ensure flux consistency
         constrainPressure(p_rgh, U, phiHbyA, rAUf, MRF);
 
-        // Cache the phase change pressure source
+        // Cache any sources
-        fvScalarMatrix Sp_rgh
+        fvScalarMatrix p_rghEqnSource
         (
-            fvModels().source
+            fvModels().sourceProxy(alpha1, rho1, p_rgh)/rho1
-            (
+          + fvModels().sourceProxy(alpha2, rho2, p_rgh)/rho2
-                volScalarField::New
-                (
-                    "1",
-                    mesh,
-                    dimensionedScalar(dimless/dimPressure, 1)
-                ),
-                p_rgh
-            )
         );
 
         // Make the fluxes relative to the mesh motion
@@ -163,11 +155,6 @@ void Foam::solvers::compressibleVoF::pressureCorrector()
         p_rghEqnComp1.ref() *= pos(alpha1);
         p_rghEqnComp2.ref() *= pos(alpha2);
 
-        p_rghEqnComp1.ref() -=
-            (fvModels().source(alpha1, mixture_.thermo1().rho())&rho1)/rho1;
-        p_rghEqnComp2.ref() -=
-            (fvModels().source(alpha2, mixture_.thermo2().rho())&rho2)/rho2;
-
         if (pimple.transonic())
         {
             p_rghEqnComp1.ref().relax();
@@ -182,7 +169,7 @@ void Foam::solvers::compressibleVoF::pressureCorrector()
             fvScalarMatrix p_rghEqnIncomp
             (
                 fvc::div(phiHbyA) - fvm::laplacian(rAUf, p_rgh)
-             == Sp_rgh
+             == p_rghEqnSource
             );
 
             {

applications/modules/incompressibleDenseParticleFluid/correctPressure.C

@@ -89,6 +89,7 @@ void Foam::solvers::incompressibleDenseParticleFluid::correctPressure()
          ==
             fvc::ddt(alphac)
           + fvc::div(alphacf*phiHbyAD)
+          + fvModels().sourceProxy(alphac, p)
         );
 
         pEqn.setReference

applications/modules/incompressibleDriftFlux/alphaSuSp.C

@@ -45,19 +45,37 @@ Foam::solvers::incompressibleDriftFlux::alphaPhi
     );
 }
 
+
 void Foam::solvers::incompressibleDriftFlux::alphaSuSp
 (
-    tmp<volScalarField::Internal>& Su,
+    tmp<volScalarField::Internal>& tSu,
-    tmp<volScalarField::Internal>& Sp
+    tmp<volScalarField::Internal>& tSp
 )
 {
-    if (divergent())
+    if (!divergent()) return;
-    {
-        // Phase change alpha1 source
-        const fvScalarMatrix alphaSup(fvModels().source(alpha1));
 
-        Su = alphaSup.Su();
+    const dimensionedScalar Szero(dimless/dimTime, 0);
-        Sp = alphaSup.Sp();
+
+    tSp = volScalarField::Internal::New("Sp", mesh, Szero);
+    tSu = volScalarField::Internal::New("Su", mesh, Szero);
+
+    volScalarField::Internal& Sp = tSp.ref();
+    volScalarField::Internal& Su = tSu.ref();
+
+    if (fvModels().addsSupToField(alpha1.name()))
+    {
+        const fvScalarMatrix alpha1Sup(fvModels().source(alpha1));
+
+        Su += alpha2()*alpha1Sup.Su();
+        Sp += alpha2()*alpha1Sup.Sp();
+    }
+
+    if (fvModels().addsSupToField(alpha2.name()))
+    {
+        const fvScalarMatrix alpha2Sup(fvModels().source(alpha2));
+
+        Su -= alpha1()*(alpha2Sup.Su() + alpha2Sup.Sp());
+        Sp += alpha1()*alpha2Sup.Sp();
     }
 }
 

applications/modules/incompressibleDriftFlux/incompressibleDriftFlux.H

@@ -133,7 +133,9 @@ protected:
         //  i.e. includes phase-fraction sources
         virtual bool divergent() const
         {
-            return fvModels().addsSupToField(alpha1.name());
+            return
+                fvModels().addsSupToField(alpha1.name())
+             || fvModels().addsSupToField(alpha2.name());
         }
 
         //- Return the mixture compressibility/density

applications/modules/incompressibleFluid/correctPressure.C

@@ -80,12 +80,18 @@ void Foam::solvers::incompressibleFluid::correctPressure()
     // Update the pressure BCs to ensure flux consistency
     constrainPressure(p, U, phiHbyA, rAtU(), MRF);
 
+    // Evaluate any volume sources
+    fvScalarMatrix p_rghEqnSource(fvModels().sourceProxy(p));
+
     // Non-orthogonal pressure corrector loop
     while (pimple.correctNonOrthogonal())
     {
         fvScalarMatrix pEqn
         (
-            fvm::laplacian(rAtU(), p) == fvc::div(phiHbyA)
+            fvm::laplacian(rAtU(), p)
+         ==
+            fvc::div(phiHbyA)
+          - p_rghEqnSource
         );
 
         pEqn.setReference

applications/modules/incompressibleMultiphaseVoF/pressureCorrector.C

@@ -80,27 +80,20 @@ void Foam::solvers::incompressibleMultiphaseVoF::pressureCorrector()
         // Update the pressure BCs to ensure flux consistency
         constrainPressure(p_rgh, U, phiHbyA, rAUf, MRF);
 
-        // Cache the phase change pressure source
+        // Evaluate any phase sources
-        fvScalarMatrix Sp_rgh
+        fvScalarMatrix p_rghEqnSource(p_rgh, dimVolume/dimTime);
-        (
+        forAll(phases, phasei)
-            fvModels().source
+        {
-            (
+            p_rghEqnSource +=
-                volScalarField::New
+                fvModels().sourceProxy(phases[phasei], p_rgh);
-                (
+        }
-                    "1",
-                    mesh,
-                    dimensionedScalar(dimless/dimPressure, 1)
-                ),
-                p_rgh
-            )
-        );
 
         while (pimple.correctNonOrthogonal())
         {
             fvScalarMatrix p_rghEqn
             (
                 fvc::div(phiHbyA) - fvm::laplacian(rAUf, p_rgh)
-             == Sp_rgh
+             == p_rghEqnSource
             );
 
             p_rghEqn.setReference

applications/modules/incompressibleVoF/alphaSuSp.C

@@ -30,17 +30,34 @@ License
 
 void Foam::solvers::incompressibleVoF::alphaSuSp
 (
-    tmp<volScalarField::Internal>& Su,
+    tmp<volScalarField::Internal>& tSu,
-    tmp<volScalarField::Internal>& Sp
+    tmp<volScalarField::Internal>& tSp
 )
 {
-    if (divergent())
+    if (!divergent()) return;
-    {
-        // Phase change alpha1 source
-        const fvScalarMatrix alphaSup(fvModels().source(alpha1));
 
-        Su = alphaSup.Su();
+    const dimensionedScalar Szero(dimless/dimTime, 0);
-        Sp = alphaSup.Sp();
+
+    tSp = volScalarField::Internal::New("Sp", mesh, Szero);
+    tSu = volScalarField::Internal::New("Su", mesh, Szero);
+
+    volScalarField::Internal& Sp = tSp.ref();
+    volScalarField::Internal& Su = tSu.ref();
+
+    if (fvModels().addsSupToField(alpha1.name()))
+    {
+        const fvScalarMatrix alpha1Sup(fvModels().source(alpha1));
+
+        Su += alpha2()*alpha1Sup.Su();
+        Sp += alpha2()*alpha1Sup.Sp();
+    }
+
+    if (fvModels().addsSupToField(alpha2.name()))
+    {
+        const fvScalarMatrix alpha2Sup(fvModels().source(alpha2));
+
+        Su -= alpha1()*(alpha2Sup.Su() + alpha2Sup.Sp());
+        Sp += alpha1()*alpha2Sup.Sp();
     }
 }
 

applications/modules/incompressibleVoF/fvModels/VoFCavitation/VoFCavitation.C

@@ -68,9 +68,7 @@ Foam::fv::VoFCavitation::VoFCavitation
         )
     ),
 
-    cavitation_(cavitationModel::New(dict, mixture_)),
+    cavitation_(cavitationModel::New(dict, mixture_))
-
-    alphaName_(mixture_.alpha1().name())
 {}
 
 
@@ -78,14 +76,19 @@ Foam::fv::VoFCavitation::VoFCavitation
 
 Foam::wordList Foam::fv::VoFCavitation::addSupFields() const
 {
-    return {alphaName_, "p_rgh", "U"};
+    return
+    {
+        mixture_.alpha1().name(),
+        mixture_.alpha2().name(),
+        "U"
+    };
 }
 
 
 void Foam::fv::VoFCavitation::addSup
 (
-    fvMatrix<scalar>& eqn,
+    const volScalarField& alpha,
-    const word& fieldName
+    fvMatrix<scalar>& eqn
 ) const
 {
     if (debug)
@@ -93,63 +96,70 @@ void Foam::fv::VoFCavitation::addSup
         Info<< type() << ": applying source to " << eqn.psi().name() << endl;
     }
 
-    if (fieldName == alphaName_)
+    if (&alpha == &mixture_.alpha1() || &alpha == &mixture_.alpha2())
     {
-        const volScalarField::Internal alpha1Coeff
+        const volScalarField& alpha1 = mixture_.alpha1();
-        (
+        const volScalarField& alpha2 = mixture_.alpha2();
-            1.0/mixture_.rho1()
-          - mixture_.alpha1()()*(1.0/mixture_.rho1() - 1.0/mixture_.rho2())
-        );
 
-        const Pair<tmp<volScalarField::Internal>> mDot12Alpha
+        const dimensionedScalar& rho =
-        (
+            &alpha == &alpha1 ? mixture_.rho1() : mixture_.rho2();
-            cavitation_->mDot12Alpha()
-        );
 
-        const volScalarField::Internal vDot1Alpha(alpha1Coeff*mDot12Alpha[0]());
+        const scalar s = &alpha == &alpha1 ? +1 : -1;
-        const volScalarField::Internal vDot2Alpha(alpha1Coeff*mDot12Alpha[1]());
 
-        eqn += fvm::Sp(-vDot2Alpha - vDot1Alpha, eqn.psi()) + vDot1Alpha;
+        // Volume-fraction linearisation
+        if (&alpha == &eqn.psi())
+        {
+            const Pair<tmp<volScalarField::Internal>> mDot12Alpha
+            (
+                cavitation_->mDot12Alpha()
+            );
+            const volScalarField::Internal vDot1Alpha2(mDot12Alpha[0]/rho);
+            const volScalarField::Internal vDot2Alpha1(mDot12Alpha[1]/rho);
+
+            eqn +=
+                (&alpha == &alpha1 ? vDot1Alpha2 : vDot2Alpha1)
+              - fvm::Sp(vDot1Alpha2 + vDot2Alpha1, eqn.psi());
+        }
+
+        // Pressure linearisation
+        else if (eqn.psi().name() == "p_rgh")
+        {
+            const Pair<tmp<volScalarField::Internal>> mDot12P
+            (
+                cavitation_->mDot12P()
+            );
+            const volScalarField::Internal vDot1P(mDot12P[0]/rho);
+            const volScalarField::Internal vDot2P(mDot12P[1]/rho);
+
+            const volScalarField::Internal& rho =
+                mesh().lookupObject<volScalarField>("rho");
+            const volScalarField::Internal& gh =
+                mesh().lookupObject<volScalarField>("gh");
+
+            eqn +=
+                fvm::Sp(s*(vDot1P - vDot2P), eqn.psi())
+              + s*(vDot1P - vDot2P)*rho*gh
+              - s*(vDot1P - vDot2P)*cavitation_->pSat();
+        }
+
+        // Explicit non-linearised value. Probably not used.
+        else
+        {
+            const Pair<tmp<volScalarField::Internal>> mDot12Alpha
+            (
+                cavitation_->mDot12Alpha()
+            );
+            const volScalarField::Internal vDot1(mDot12Alpha[0]*alpha2/rho);
+            const volScalarField::Internal vDot2(mDot12Alpha[1]*alpha1/rho);
+
+            eqn += s*(vDot1 - vDot2);
+        }
     }
-}
+    else
-
-
-void Foam::fv::VoFCavitation::addSup
-(
-    const volScalarField&,
-    fvMatrix<scalar>& eqn,
-    const word& fieldName
-) const
-{
-    if (debug)
     {
-        Info<< type() << ": applying source to " << eqn.psi().name() << endl;
+        FatalErrorInFunction
-    }
+            << "Support for field " << alpha.name() << " is not implemented"
-
+            << exit(FatalError);
-    if (fieldName == "p_rgh")
-    {
-        const volScalarField::Internal& rho =
-            mesh().lookupObject<volScalarField>("rho");
-
-        const volScalarField::Internal& gh =
-            mesh().lookupObject<volScalarField>("gh");
-
-        const dimensionedScalar pCoeff
-        (
-            1.0/mixture_.rho1() - 1.0/mixture_.rho2()
-        );
-
-        const Pair<tmp<volScalarField::Internal>> mDot12P
-        (
-            cavitation_->mDot12P()
-        );
-
-        const volScalarField::Internal vDot1P(pCoeff*mDot12P[0]);
-        const volScalarField::Internal vDot2P(pCoeff*mDot12P[1]);
-
-        eqn +=
-            (vDot2P - vDot1P)*(cavitation_->pSat() - rho*gh)
-          - fvm::Sp(vDot2P - vDot1P, eqn.psi());
     }
 }
 
@@ -157,8 +167,8 @@ void Foam::fv::VoFCavitation::addSup
 void Foam::fv::VoFCavitation::addSup
 (
     const volScalarField& rho,
-    fvMatrix<vector>& eqn,
+    const volVectorField& U,
-    const word& fieldName
+    fvMatrix<vector>& eqn
 ) const
 {
     if (debug)
@@ -166,12 +176,25 @@ void Foam::fv::VoFCavitation::addSup
         Info<< type() << ": applying source to " << eqn.psi().name() << endl;
     }
 
-    if (fieldName == "U")
+    if (U.name() == "U")
     {
         const surfaceScalarField& rhoPhi =
             mesh().lookupObject<surfaceScalarField>("rhoPhi");
 
-        eqn += fvm::Sp(fvc::ddt(rho) + fvc::div(rhoPhi), eqn.psi());
+        if (&U == &eqn.psi())
+        {
+            eqn += fvm::Sp(fvc::ddt(rho) + fvc::div(rhoPhi), eqn.psi());
+        }
+        else
+        {
+            eqn += (fvc::ddt(rho) + fvc::div(rhoPhi))*U;
+        }
+    }
+    else
+    {
+        FatalErrorInFunction
+            << "Support for field " << U.name() << " is not implemented"
+            << exit(FatalError);
     }
 }
 

applications/modules/incompressibleVoF/fvModels/VoFCavitation/VoFCavitation.H

@@ -102,11 +102,9 @@ class VoFCavitation
         //- Reference to the mixture properties
         const incompressibleTwoPhaseVoFMixture& mixture_;
 
+        //- The cavitation model
         autoPtr<cavitationModel> cavitation_;
 
-        //- The name of the VoF phase-fraction
-        word alphaName_;
-
 
 public:
 
@@ -138,27 +136,19 @@ public:
         //  to the transport equation
         virtual wordList addSupFields() const;
 
-        //- Add implicit/explicit contributions to VoF phase-fraction equation
+        //- Add a source to the phase continuity equation
         virtual void addSup
         (
-            fvMatrix<scalar>& eqn,
+            const volScalarField& alpha,
-            const word& fieldName
+            fvMatrix<scalar>& eqn
         ) const;
 
-        //- Add implicit/explicit contributions to p_rgh equation
+        //- Add a source to the mixture momentum equation
-        virtual void addSup
-        (
-            const volScalarField& psi,
-            fvMatrix<scalar>& eqn,
-            const word& fieldName
-        ) const;
-
-        //- Add contribution to U equation
         virtual void addSup
         (
             const volScalarField& rho,
-            fvMatrix<vector>& eqn,
+            const volVectorField& U,
-            const word& fieldName
+            fvMatrix<vector>& eqn
         ) const;
 
 

applications/modules/incompressibleVoF/fvModels/VoFTurbulenceDamping/VoFTurbulenceDamping.C

@@ -118,8 +118,8 @@ Foam::wordList Foam::fv::VoFTurbulenceDamping::addSupFields() const
 
 void Foam::fv::VoFTurbulenceDamping::addSup
 (
-    fvMatrix<scalar>& eqn,
+    const volScalarField& field,
-    const word& fieldName
+    fvMatrix<scalar>& eqn
 ) const
 {
     if (debug)
@@ -133,12 +133,12 @@ void Foam::fv::VoFTurbulenceDamping::addSup
       + mixture_.alpha2()()*sqr(mixture_.nuModel2().nu()()())
     );
 
-    if (fieldName == "epsilon")
+    if (field.name() == "epsilon")
     {
         eqn += mixture_.interfaceFraction()*C2_*aSqrnu*turbulence_.k()()
            /pow4(delta_);
     }
-    else if (fieldName == "omega")
+    else if (field.name() == "omega")
     {
         eqn += mixture_.interfaceFraction()*beta_*aSqrnu/(sqr(betaStar_)
            *pow4(delta_));
@@ -146,7 +146,7 @@ void Foam::fv::VoFTurbulenceDamping::addSup
     else
     {
         FatalErrorInFunction
-            << "Support for field " << fieldName << " is not implemented"
+            << "Support for field " << field.name() << " is not implemented"
             << exit(FatalError);
     }
 }
@@ -155,9 +155,9 @@ void Foam::fv::VoFTurbulenceDamping::addSup
 void Foam::fv::VoFTurbulenceDamping::addSup
 (
     const volScalarField& alpha,
-    const volScalarField&,
+    const volScalarField& rho,
-    fvMatrix<scalar>& eqn,
+    const volScalarField& field,
-    const word& fieldName
+    fvMatrix<scalar>& eqn
 ) const
 {
     if (debug)
@@ -182,12 +182,12 @@ void Foam::fv::VoFTurbulenceDamping::addSup
             << exit(FatalError);
     }
 
-    if (fieldName == IOobject::groupName("epsilon", phaseName_))
+    if (field.name() == IOobject::groupName("epsilon", phaseName_))
     {
         eqn += mixture_.interfaceFraction()*C2_*taSqrnu*turbulence_.k()()
            /pow4(delta_);
     }
-    else if (fieldName == IOobject::groupName("omega", phaseName_))
+    else if (field.name() == IOobject::groupName("omega", phaseName_))
     {
         eqn += mixture_.interfaceFraction()*beta_*taSqrnu
            /(sqr(betaStar_)*pow4(delta_));
@@ -195,7 +195,7 @@ void Foam::fv::VoFTurbulenceDamping::addSup
     else
     {
         FatalErrorInFunction
-            << "Support for field " << fieldName << " is not implemented"
+            << "Support for field " << field.name() << " is not implemented"
             << exit(FatalError);
     }
 }

applications/modules/incompressibleVoF/fvModels/VoFTurbulenceDamping/VoFTurbulenceDamping.H

@@ -144,25 +144,30 @@ public:
 
     // Member Functions
 
-        //- Return the list of fields for which the option adds source term
+        // Checks
-        //  to the transport equation
-        virtual wordList addSupFields() const;
 
-        //- Add explicit contribution to epsilon or omega equation
+            //- Return the list of fields for which the option adds source term
-        virtual void addSup
+            //  to the transport equation
-        (
+            virtual wordList addSupFields() const;
-            fvMatrix<scalar>& eqn,
-            const word& fieldName
-        ) const;
 
-        //- Add explicit contribution to phase epsilon or omega equation
+
-        virtual void addSup
+        // Sources
-        (
+
-            const volScalarField& alpha,
+            //- Add explicit contribution to epsilon or omega equation
-            const volScalarField& rho,
+            virtual void addSup
-            fvMatrix<scalar>& eqn,
+            (
-            const word& fieldName
+                const volScalarField& field,
-        ) const;
+                fvMatrix<scalar>& eqn
+            ) const;
+
+            //- Add explicit contribution to phase epsilon or omega equation
+            virtual void addSup
+            (
+                const volScalarField& alpha,
+                const volScalarField& rho,
+                const volScalarField& field,
+                fvMatrix<scalar>& eqn
+            ) const;
 
 
         // Mesh changes

applications/modules/incompressibleVoF/incompressibleVoF.H

@@ -122,7 +122,9 @@ protected:
         //  i.e. includes phase-fraction sources
         virtual bool divergent() const
         {
-            return fvModels().addsSupToField(alpha1.name());
+            return
+                fvModels().addsSupToField(alpha1.name())
+             || fvModels().addsSupToField(alpha2.name());
         }
 
         //- Return the mixture compressibility/density

applications/modules/isothermalFilm/fvModels/filmCloudTransfer/filmCloudTransfer.C

@@ -138,8 +138,8 @@ inline Foam::fv::filmCloudTransfer::CloudToFilmTransferRate
 void Foam::fv::filmCloudTransfer::addSup
 (
     const volScalarField& rho,
-    fvMatrix<scalar>& eqn,
+    const volScalarField& alpha,
-    const word& fieldName
+    fvMatrix<scalar>& eqn
 ) const
 {
     if (debug)
@@ -147,7 +147,7 @@ void Foam::fv::filmCloudTransfer::addSup
         Info<< type() << ": applying source to " << eqn.psi().name() << endl;
     }
 
-    if (fieldName == film_.alpha.name())
+    if (&alpha == &film_.alpha && &eqn.psi() == &film_.alpha)
     {
         eqn += CloudToFilmTransferRate<scalar>(massFromCloud_, dimMass);
 
@@ -159,7 +159,7 @@ void Foam::fv::filmCloudTransfer::addSup
     else
     {
         FatalErrorInFunction
-            << "Support for field " << fieldName << " is not implemented"
+            << "Support for field " << alpha.name() << " is not implemented"
             << exit(FatalError);
     }
 }
@@ -169,8 +169,8 @@ void Foam::fv::filmCloudTransfer::addSup
 (
     const volScalarField& alpha,
     const volScalarField& rho,
-    fvMatrix<scalar>& eqn,
+    const volScalarField& he,
-    const word& fieldName
+    fvMatrix<scalar>& eqn
 ) const
 {
     if (debug)
@@ -178,7 +178,7 @@ void Foam::fv::filmCloudTransfer::addSup
         Info<< type() << ": applying source to " << eqn.psi().name() << endl;
     }
 
-    if (fieldName == film_.thermo.he().name())
+    if (&he == &film_.thermo.he() && &eqn.psi() == &film_.thermo.he())
     {
         eqn += CloudToFilmTransferRate<scalar>(energyFromCloud_, dimEnergy);
 
@@ -190,7 +190,7 @@ void Foam::fv::filmCloudTransfer::addSup
     else
     {
         FatalErrorInFunction
-            << "Support for field " << fieldName << " is not implemented"
+            << "Support for field " << he.name() << " is not implemented"
             << exit(FatalError);
     }
 }
@@ -200,8 +200,8 @@ void Foam::fv::filmCloudTransfer::addSup
 (
     const volScalarField& alpha,
     const volScalarField& rho,
-    fvMatrix<vector>& eqn,
+    const volVectorField& U,
-    const word& fieldName
+    fvMatrix<vector>& eqn
 ) const
 {
     if (debug)
@@ -209,7 +209,7 @@ void Foam::fv::filmCloudTransfer::addSup
         Info<< type() << ": applying source to " << eqn.psi().name() << endl;
     }
 
-    if (fieldName == film_.U.name())
+    if (&U == &film_.U && &U == &film_.U)
     {
         eqn += CloudToFilmTransferRate<vector>(momentumFromCloud_, dimMomentum);
 
@@ -221,7 +221,7 @@ void Foam::fv::filmCloudTransfer::addSup
     else
     {
         FatalErrorInFunction
-            << "Support for field " << fieldName << " is not implemented"
+            << "Support for field " << U.name() << " is not implemented"
             << exit(FatalError);
     }
 }

applications/modules/isothermalFilm/fvModels/filmCloudTransfer/filmCloudTransfer.H

@@ -144,32 +144,32 @@ public:
             virtual void correct();
 
 
-        // Add explicit and implicit contributions to compressible equation
+        // Sources
 
-            //- Add explicit contribution to phase continuity
+            //- Add source to phase continuity equation
             virtual void addSup
             (
+                const volScalarField& rho,
                 const volScalarField& alpha,
-                fvMatrix<scalar>& eqn,
+                fvMatrix<scalar>& eqn
-                const word& fieldName
             ) const;
 
-            //- Add explicit contribution to phase energy equation
+            //- Add source to phase energy equation
             virtual void addSup
             (
                 const volScalarField& alpha,
                 const volScalarField& rho,
-                fvMatrix<scalar>& eqn,
+                const volScalarField& he,
-                const word& fieldName
+                fvMatrix<scalar>& eqn
             ) const;
 
-            //- Add implicit contribution to mixture momentum equation
+            //- Add source to mixture momentum equation
             virtual void addSup
             (
                 const volScalarField& alpha,
                 const volScalarField& rho,
-                fvMatrix<vector>& eqn,
+                const volVectorField& U,
-                const word& fieldName
+                fvMatrix<vector>& eqn
             ) const;
 
 

applications/modules/isothermalFilm/fvModels/filmVoFTransfer/VoFFilmTransfer.C

@@ -105,7 +105,6 @@ Foam::wordList Foam::fv::VoFFilmTransfer::addSupFields() const
 {
     return wordList
     {
-        alpha_.name(),
         thermo_.rho()().name(),
         thermo_.he().name(),
         VoF_.U.name()
@@ -241,41 +240,11 @@ inline Foam::fv::VoFFilmTransfer::filmVoFTransferRate
 }
 
 
-void Foam::fv::VoFFilmTransfer::addSup
-(
-    fvMatrix<scalar>& eqn,
-    const word& fieldName
-) const
-{
-    if (debug)
-    {
-        Info<< type() << ": applying source to " << eqn.psi().name() << endl;
-    }
-
-    if (fieldName == alpha_.name())
-    {
-        eqn +=
-            filmVoFTransferRate<scalar>
-            (
-                &filmVoFTransfer::transferRate,
-                dimVolume
-            )
-          - fvm::Sp(transferRate_, eqn.psi());
-    }
-    else
-    {
-        FatalErrorInFunction
-            << "Support for field " << fieldName << " is not implemented"
-            << exit(FatalError);
-    }
-}
-
-
 void Foam::fv::VoFFilmTransfer::addSup
 (
     const volScalarField& alpha,
-    fvMatrix<scalar>& eqn,
+    const volScalarField& rho,
-    const word& fieldName
+    fvMatrix<scalar>& eqn
 ) const
 {
     if (debug)
@@ -283,20 +252,34 @@ void Foam::fv::VoFFilmTransfer::addSup
         Info<< type() << ": applying source to " << eqn.psi().name() << endl;
     }
 
-    if (fieldName == thermo_.rho()().name())
+    if (&rho == &thermo_.rho()())
     {
+        // Explicit transfer into the VoF
         eqn +=
             filmVoFTransferRate<scalar>
             (
                 &filmVoFTransfer::rhoTransferRate,
                 dimMass
-            )
+            );
-          - fvm::Sp(alpha()*transferRate_, eqn.psi());
+
+        // Potentially implicit transfer out of the VoF
+        if (&rho == &eqn.psi())
+        {
+            eqn -= fvm::Sp(alpha()*transferRate_, eqn.psi());
+        }
+        else if (&alpha == &eqn.psi())
+        {
+            eqn -= fvm::Sp(rho()*transferRate_, eqn.psi());
+        }
+        else
+        {
+            eqn -= alpha()*rho()*transferRate_;
+        }
     }
     else
     {
         FatalErrorInFunction
-            << "Support for field " << fieldName << " is not implemented"
+            << "Support for field " << rho.name() << " is not implemented"
             << exit(FatalError);
     }
 }
@@ -306,8 +289,8 @@ void Foam::fv::VoFFilmTransfer::addSup
 (
     const volScalarField& alpha,
     const volScalarField& rho,
-    fvMatrix<scalar>& eqn,
+    const volScalarField& he,
-    const word& fieldName
+    fvMatrix<scalar>& eqn
 ) const
 {
     if (debug)
@@ -315,20 +298,30 @@ void Foam::fv::VoFFilmTransfer::addSup
         Info<< type() << ": applying source to " << eqn.psi().name() << endl;
     }
 
-    if (fieldName == thermo_.he().name())
+    if (&he == &thermo_.he())
     {
+        // Explicit transfer into the VoF
         eqn +=
             filmVoFTransferRate<scalar>
             (
                 &filmVoFTransfer::heTransferRate,
                 dimEnergy
-            )
+            );
-          - fvm::Sp(alpha()*rho()*transferRate_, eqn.psi());
+
+        // Potentially implicit transfer out of the VoF
+        if (&he == &eqn.psi())
+        {
+            eqn -= fvm::Sp(alpha()*rho()*transferRate_, eqn.psi());
+        }
+        else
+        {
+            eqn -= alpha()*rho()*he*transferRate_;
+        }
     }
     else
     {
         FatalErrorInFunction
-            << "Support for field " << fieldName << " is not implemented"
+            << "Support for field " << he.name() << " is not implemented"
             << exit(FatalError);
     }
 }
@@ -337,8 +330,8 @@ void Foam::fv::VoFFilmTransfer::addSup
 void Foam::fv::VoFFilmTransfer::addSup
 (
     const volScalarField& rho,
-    fvMatrix<vector>& eqn,
+    const volVectorField& U,
-    const word& fieldName
+    fvMatrix<vector>& eqn
 ) const
 {
     if (debug)
@@ -346,13 +339,32 @@ void Foam::fv::VoFFilmTransfer::addSup
         Info<< type() << ": applying source to " << eqn.psi().name() << endl;
     }
 
-    eqn +=
+    if (&U == &VoF_.U)
-        filmVoFTransferRate<vector>
+    {
-        (
+        // Explicit transfer into the VoF
-            &filmVoFTransfer::UTransferRate,
+        eqn +=
-            dimMass*dimVelocity
+            filmVoFTransferRate<vector>
-        )
+            (
-      - fvm::Sp(alpha_()*thermo_.rho()()*transferRate_, eqn.psi());
+                &filmVoFTransfer::UTransferRate,
+                dimMass*dimVelocity
+            );
+
+        // Potentially implicit transfer out of the VoF
+        if (&U == &eqn.psi())
+        {
+            eqn -= fvm::Sp(alpha_()*rho()*transferRate_, eqn.psi());
+        }
+        else
+        {
+            eqn -= alpha_()*rho()*U*transferRate_;
+        }
+    }
+    else
+    {
+        FatalErrorInFunction
+            << "Support for field " << U.name() << " is not implemented"
+            << exit(FatalError);
+    }
 }
 
 

applications/modules/isothermalFilm/fvModels/filmVoFTransfer/VoFFilmTransfer.H

@@ -172,21 +172,14 @@ public:
             virtual void correct();
 
 
-        // Add explicit and implicit contributions to compressible equation
+        // Sources
 
-            //- Add implicit contribution to phase-fraction equation
+            //- Add implicit contribution to phase continuity equation
-            virtual void addSup
-            (
-                fvMatrix<scalar>& eqn,
-                const word& fieldName
-            ) const;
-
-            //- Add implicit contribution to phase density equation
             virtual void addSup
             (
                 const volScalarField& alpha,
-                fvMatrix<scalar>& eqn,
+                const volScalarField& rho,
-                const word& fieldName
+                fvMatrix<scalar>& eqn
             ) const;
 
             //- Add implicit contribution to phase energy equation
@@ -194,16 +187,16 @@ public:
             (
                 const volScalarField& alpha,
                 const volScalarField& rho,
-                fvMatrix<scalar>& eqn,
+                const volScalarField& he,
-                const word& fieldName
+                fvMatrix<scalar>& eqn
             ) const;
 
             //- Add implicit contribution to mixture momentum equation
             virtual void addSup
             (
                 const volScalarField& rho,
-                fvMatrix<vector>& eqn,
+                const volVectorField& U,
-                const word& fieldName
+                fvMatrix<vector>& eqn
             ) const;
 
 

applications/modules/isothermalFilm/fvModels/filmVoFTransfer/filmVoFTransfer.C

@@ -243,8 +243,8 @@ inline Foam::fv::filmVoFTransfer::VoFToFilmTransferRate
 void Foam::fv::filmVoFTransfer::addSup
 (
     const volScalarField& rho,
-    fvMatrix<scalar>& eqn,
+    const volScalarField& alpha,
-    const word& fieldName
+    fvMatrix<scalar>& eqn
 ) const
 {
     if (debug)
@@ -252,20 +252,30 @@ void Foam::fv::filmVoFTransfer::addSup
         Info<< type() << ": applying source to " << eqn.psi().name() << endl;
     }
 
-    if (fieldName == film_.alpha.name())
+    if (&alpha == &film_.alpha)
     {
+        // Explicit transfer into the film
         eqn +=
             VoFToFilmTransferRate<scalar>
             (
                 &VoFFilmTransfer::rhoTransferRate,
                 dimMass
-            )
+            );
-          - fvm::Sp(transferRate_*rho(), eqn.psi());
+
+        // Potentially implicit transfer out of the film
+        if (&alpha == &eqn.psi())
+        {
+            eqn -= fvm::Sp(rho()*transferRate_, eqn.psi());
+        }
+        else
+        {
+            eqn -= alpha()*rho()*transferRate_;
+        }
     }
     else
     {
         FatalErrorInFunction
-            << "Support for field " << fieldName << " is not implemented"
+            << "Support for field " << alpha.name() << " is not implemented"
             << exit(FatalError);
     }
 }
@@ -275,8 +285,8 @@ void Foam::fv::filmVoFTransfer::addSup
 (
     const volScalarField& alpha,
     const volScalarField& rho,
-    fvMatrix<scalar>& eqn,
+    const volScalarField& he,
-    const word& fieldName
+    fvMatrix<scalar>& eqn
 ) const
 {
     if (debug)
@@ -284,20 +294,30 @@ void Foam::fv::filmVoFTransfer::addSup
         Info<< type() << ": applying source to " << eqn.psi().name() << endl;
     }
 
-    if (fieldName == film_.thermo.he().name())
+    if (&he == &film_.thermo.he())
     {
+        // Explicit transfer into the film
         eqn +=
             VoFToFilmTransferRate<scalar>
             (
                 &VoFFilmTransfer::heTransferRate,
                 dimEnergy
-            )
+            );
-          - fvm::Sp(alpha()*rho()*transferRate_, eqn.psi());
+
+        // Potentially implicit transfer out of the film
+        if (&he == &eqn.psi())
+        {
+            eqn -= fvm::Sp(alpha()*rho()*transferRate_, eqn.psi());
+        }
+        else
+        {
+            eqn -= alpha()*rho()*he*transferRate_;
+        }
     }
     else
     {
         FatalErrorInFunction
-            << "Support for field " << fieldName << " is not implemented"
+            << "Support for field " << he.name() << " is not implemented"
             << exit(FatalError);
     }
 }
@@ -307,8 +327,8 @@ void Foam::fv::filmVoFTransfer::addSup
 (
     const volScalarField& alpha,
     const volScalarField& rho,
-    fvMatrix<vector>& eqn,
+    const volVectorField& U,
-    const word& fieldName
+    fvMatrix<vector>& eqn
 ) const
 {
     if (debug)
@@ -316,13 +336,32 @@ void Foam::fv::filmVoFTransfer::addSup
         Info<< type() << ": applying source to " << eqn.psi().name() << endl;
     }
 
-    eqn +=
+    if (&U == &film_.U)
-        VoFToFilmTransferRate<vector>
+    {
-        (
+        // Explicit transfer into the film
-            &VoFFilmTransfer::UTransferRate,
+        eqn +=
-            dimMomentum
+            VoFToFilmTransferRate<vector>
-        )
+            (
-      - fvm::Sp(alpha()*rho()*transferRate_, eqn.psi());
+                &VoFFilmTransfer::UTransferRate,
+                dimMomentum
+            );
+
+        // Potentially implicit transfer out of the film
+        if (&U == &eqn.psi())
+        {
+            eqn -= fvm::Sp(alpha()*rho()*transferRate_, eqn.psi());
+        }
+        else
+        {
+            eqn -= alpha()*rho()*U()*transferRate_;
+        }
+    }
+    else
+    {
+        FatalErrorInFunction
+            << "Support for field " << U.name() << " is not implemented"
+            << exit(FatalError);
+    }
 }
 
 

applications/modules/isothermalFilm/fvModels/filmVoFTransfer/filmVoFTransfer.H

@@ -147,32 +147,32 @@ public:
             virtual void correct();
 
 
-        // Add explicit and implicit contributions to compressible equation
+        // Sources
 
-            //- Add explicit contribution to phase continuity
+            //- Add source to phase continuity equation
             virtual void addSup
             (
+                const volScalarField& rho,
                 const volScalarField& alpha,
-                fvMatrix<scalar>& eqn,
+                fvMatrix<scalar>& eqn
-                const word& fieldName
             ) const;
 
-            //- Add explicit contribution to phase energy equation
+            //- Add source to phase energy equation
             virtual void addSup
             (
                 const volScalarField& alpha,
                 const volScalarField& rho,
-                fvMatrix<scalar>& eqn,
+                const volScalarField& he,
-                const word& fieldName
+                fvMatrix<scalar>& eqn
             ) const;
 
-            //- Add implicit contribution to mixture momentum equation
+            //- Add source to mixture momentum equation
             virtual void addSup
             (
                 const volScalarField& alpha,
                 const volScalarField& rho,
-                fvMatrix<vector>& eqn,
+                const volVectorField& U,
-                const word& fieldName
+                fvMatrix<vector>& eqn
             ) const;
 
 

applications/modules/isothermalFluid/correctBuoyantPressure.C

@@ -142,7 +142,7 @@ void Foam::solvers::isothermalFluid::correctBuoyantPressure()
             fvc::ddt(rho) + psi*correction(fvm::ddt(p_rgh))
           + fvc::div(phiHbyA) + fvm::div(phid, p_rgh)
          ==
-            fvModels().source(psi, p_rgh, rho.name())
+            fvModels().sourceProxy(rho, p_rgh)
         );
 
         while (pimple.correctNonOrthogonal())
@@ -183,7 +183,7 @@ void Foam::solvers::isothermalFluid::correctBuoyantPressure()
             fvc::ddt(rho) + psi*correction(fvm::ddt(p_rgh))
           + fvc::div(phiHbyA)
          ==
-            fvModels().source(psi, p_rgh, rho.name())
+            fvModels().sourceProxy(rho, p_rgh)
         );
 
         while (pimple.correctNonOrthogonal())

applications/modules/isothermalFluid/correctPressure.C

@@ -130,7 +130,7 @@ void Foam::solvers::isothermalFluid::correctPressure()
             fvc::ddt(rho) + psi*correction(fvm::ddt(p))
           + fvc::div(phiHbyA) + fvm::div(phid, p)
          ==
-            fvModels().source(psi, p, rho.name())
+            fvModels().sourceProxy(rho, p)
         );
 
         while (pimple.correctNonOrthogonal())
@@ -179,7 +179,7 @@ void Foam::solvers::isothermalFluid::correctPressure()
             fvc::ddt(rho) + psi*correction(fvm::ddt(p))
           + fvc::div(phiHbyA)
          ==
-            fvModels().source(psi, p, rho.name())
+            fvModels().sourceProxy(rho, p)
         );
 
         while (pimple.correctNonOrthogonal())

applications/modules/multiphaseEuler/compressibilityEqns.C

@@ -102,7 +102,7 @@ Foam::solvers::multiphaseEuler::compressibilityEqns
         // Option sources
         if (fvModels().addsSupToField(rho.name()))
         {
-            pEqnComp -= (fvModels().source(alpha, rho) & rho)/rho;
+            pEqnComp -= fvModels().sourceProxy(alpha, rho, p_rgh)/rho;
         }
 
         // Mass transfer

applications/modules/multiphaseEuler/fvModels/interfaceTurbulenceDamping/interfaceTurbulenceDamping.C

@@ -131,8 +131,8 @@ template<class RhoType>
 void Foam::fv::interfaceTurbulenceDamping::addRhoSup
 (
     const RhoType& rho,
-    fvMatrix<scalar>& eqn,
+    const volScalarField& field,
-    const word& fieldName
+    fvMatrix<scalar>& eqn
 ) const
 {
     if (debug)
@@ -154,12 +154,12 @@ void Foam::fv::interfaceTurbulenceDamping::addRhoSup
             movingPhases[phasei]*sqr(movingPhases[phasei].thermo().nu()()());
     }
 
-    if (fieldName == "epsilon")
+    if (field.name() == "epsilon")
     {
         eqn += rho*interfaceFraction(phase_)*C2_*aSqrnu*turbulence_.k()()
             /pow4(delta_);
     }
-    else if (fieldName == "omega")
+    else if (field.name() == "omega")
     {
         eqn += rho*interfaceFraction(phase_)*beta_*aSqrnu
             /(sqr(betaStar_)*pow4(delta_));
@@ -167,7 +167,7 @@ void Foam::fv::interfaceTurbulenceDamping::addRhoSup
     else
     {
         FatalErrorInFunction
-            << "Support for field " << fieldName << " is not implemented"
+            << "Support for field " << field.name() << " is not implemented"
             << exit(FatalError);
     }
 }
@@ -241,22 +241,22 @@ Foam::wordList Foam::fv::interfaceTurbulenceDamping::addSupFields() const
 
 void Foam::fv::interfaceTurbulenceDamping::addSup
 (
-    fvMatrix<scalar>& eqn,
+    const volScalarField& field,
-    const word& fieldName
+    fvMatrix<scalar>& eqn
 ) const
 {
-    addRhoSup(one(), eqn, fieldName);
+    addRhoSup(one(), field, eqn);
 }
 
 
 void Foam::fv::interfaceTurbulenceDamping::addSup
 (
     const volScalarField& rho,
-    fvMatrix<scalar>& eqn,
+    const volScalarField& field,
-    const word& fieldName
+    fvMatrix<scalar>& eqn
 ) const
 {
-    addRhoSup(rho(), eqn, fieldName);
+    addRhoSup(rho(), field, eqn);
 }
 
 
@@ -264,8 +264,8 @@ void Foam::fv::interfaceTurbulenceDamping::addSup
 (
     const volScalarField& alpha,
     const volScalarField& rho,
-    fvMatrix<scalar>& eqn,
+    const volScalarField& field,
-    const word& fieldName
+    fvMatrix<scalar>& eqn
 ) const
 {
     if (debug)
@@ -273,17 +273,17 @@ void Foam::fv::interfaceTurbulenceDamping::addSup
         Info<< type() << ": applying source to " << eqn.psi().name() << endl;
     }
 
-    volScalarField::Internal aSqrnu
+    const volScalarField::Internal aSqrnu
     (
         alpha*sqr(phase_.thermo().nu()()())
     );
 
-    if (fieldName == IOobject::groupName("epsilon", phaseName_))
+    if (field.name() == IOobject::groupName("epsilon", phaseName_))
     {
         eqn += rho()*interfaceFraction(alpha)
             *C2_*aSqrnu*turbulence_.k()()/pow4(delta_);
     }
-    else if (fieldName == IOobject::groupName("omega", phaseName_))
+    else if (field.name() == IOobject::groupName("omega", phaseName_))
     {
         eqn += rho()*interfaceFraction(alpha)
             *beta_*aSqrnu/(sqr(betaStar_)*pow4(delta_));
@@ -291,7 +291,7 @@ void Foam::fv::interfaceTurbulenceDamping::addSup
     else
     {
         FatalErrorInFunction
-            << "Support for field " << fieldName << " is not implemented"
+            << "Support for field " << field.name() << " is not implemented"
             << exit(FatalError);
     }
 }

applications/modules/multiphaseEuler/fvModels/interfaceTurbulenceDamping/interfaceTurbulenceDamping.H

@@ -131,8 +131,8 @@ class interfaceTurbulenceDamping
         void addRhoSup
         (
             const RhoType& rho,
-            fvMatrix<scalar>& eqn,
+            const volScalarField& field,
-            const word& fieldName
+            fvMatrix<scalar>& eqn
         ) const;
 
 
@@ -162,34 +162,38 @@ public:
 
     // Member Functions
 
-        //- Return the list of fields for which the option adds source term
+        // Checks
-        //  to the transport equation
-        virtual wordList addSupFields() const;
 
-        //- Add explicit contribution to mixture epsilon or omega equation
+            //- Return the list of fields for which the option adds source term
-        virtual void addSup
+            //  to the transport equation
-        (
+            virtual wordList addSupFields() const;
-            fvMatrix<scalar>& eqn,
-            const word& fieldName
-        ) const;
 
-        //- Add explicit contribution to compressible
-        //  mixture epsilon or omega equation
-        virtual void addSup
-        (
-            const volScalarField& rho,
-            fvMatrix<scalar>& eqn,
-            const word& fieldName
-        ) const;
 
-        //- Add explicit contribution to phase epsilon or omega equation
+        // Sources
-        virtual void addSup
+
-        (
+            //- Add source to mixture epsilon or omega equation
-            const volScalarField& alpha,
+            virtual void addSup
-            const volScalarField& rho,
+            (
-            fvMatrix<scalar>& eqn,
+                const volScalarField& field,
-            const word& fieldName
+                fvMatrix<scalar>& eqn
-        ) const;
+            ) const;
+
+            //- Add source to compressible mixture epsilon or omega equation
+            virtual void addSup
+            (
+                const volScalarField& rho,
+                const volScalarField& field,
+                fvMatrix<scalar>& eqn
+            ) const;
+
+            //- Add source to phase epsilon or omega equation
+            virtual void addSup
+            (
+                const volScalarField& alpha,
+                const volScalarField& rho,
+                const volScalarField& field,
+                fvMatrix<scalar>& eqn
+            ) const;
 
 
         // Mesh changes

applications/modules/multiphaseEuler/fvModels/phaseTurbulenceStabilisation/phaseTurbulenceStabilisation.C

@@ -50,10 +50,11 @@ namespace Foam
 
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
-void Foam::fv::phaseTurbulenceStabilisation::addSup
+void Foam::fv::phaseTurbulenceStabilisation::addAlphaRhoSup
 (
     const volScalarField& alpha,
     const volScalarField& rho,
+    const volScalarField& field,
     fvMatrix<scalar>& eqn,
     tmp<volScalarField>
     (phaseCompressible::momentumTransportModel::*psi)() const
@@ -184,36 +185,39 @@ void Foam::fv::phaseTurbulenceStabilisation::addSup
 (
     const volScalarField& alpha,
     const volScalarField& rho,
-    fvMatrix<scalar>& eqn,
+    const volScalarField& field,
-    const word& fieldName
+    fvMatrix<scalar>& eqn
 ) const
 {
-    if (fieldName == IOobject::groupName("k", phaseName_))
+    if (field.name() == IOobject::groupName("k", phaseName_))
     {
-        addSup
+        addAlphaRhoSup
         (
             alpha,
             rho,
+            field,
             eqn,
             &phaseCompressible::momentumTransportModel::k
         );
     }
-    else if (fieldName == IOobject::groupName("epsilon", phaseName_))
+    else if (field.name() == IOobject::groupName("epsilon", phaseName_))
     {
-        addSup
+        addAlphaRhoSup
         (
             alpha,
             rho,
+            field,
             eqn,
             &phaseCompressible::momentumTransportModel::epsilon
         );
     }
-    else if (fieldName == IOobject::groupName("omega", phaseName_))
+    else if (field.name() == IOobject::groupName("omega", phaseName_))
     {
-        addSup
+        addAlphaRhoSup
         (
             alpha,
             rho,
+            field,
             eqn,
             &phaseCompressible::momentumTransportModel::omega
         );
@@ -221,7 +225,7 @@ void Foam::fv::phaseTurbulenceStabilisation::addSup
     else
     {
         FatalErrorInFunction
-            << "Support for field " << fieldName << " is not implemented"
+            << "Support for field " << field.name() << " is not implemented"
             << exit(FatalError);
     }
 }

applications/modules/multiphaseEuler/fvModels/phaseTurbulenceStabilisation/phaseTurbulenceStabilisation.H

@@ -112,10 +112,11 @@ class phaseTurbulenceStabilisation
     // Private Member Functions
 
         //- Add contribution to phase psi equation
-        void addSup
+        void addAlphaRhoSup
         (
             const volScalarField& alpha,
             const volScalarField& rho,
+            const volScalarField& field,
             fvMatrix<scalar>& eqn,
             tmp<volScalarField>
             (phaseCompressible::momentumTransportModel::*psi)() const
@@ -148,20 +149,23 @@ public:
 
     // Member Functions
 
-        //- Return the list of fields for which the option adds source term
+        // Checks
-        //  to the transport equation
-        virtual wordList addSupFields() const;
 
-        using fvModel::addSup;
+            //- Return the list of fields for which the option adds source term
+            //  to the transport equation
+            virtual wordList addSupFields() const;
 
-        //- Add contribution to phase k, epsilon or omega equation
+
-        virtual void addSup
+        // Sources
-        (
+
-            const volScalarField& alpha,
+            //- Add contribution to phase k, epsilon or omega equation
-            const volScalarField& rho,
+            virtual void addSup
-            fvMatrix<scalar>& eqn,
+            (
-            const word& fieldName
+                const volScalarField& alpha,
-        ) const;
+                const volScalarField& rho,
+                const volScalarField& field,
+                fvMatrix<scalar>& eqn
+            ) const;
 
 
         // Mesh changes

applications/modules/twoPhaseSolver/pressureCorrector.C

@@ -83,19 +83,11 @@ void Foam::solvers::twoPhaseSolver::incompressiblePressureCorrector
         // Update the pressure BCs to ensure flux consistency
         constrainPressure(p_rgh, U, phiHbyA, rAUf, MRF);
 
-        // Cache the phase change pressure source
+        // Cache any sources
-        fvScalarMatrix Sp_rgh
+        fvScalarMatrix p_rghEqnSource
         (
-            fvModels().source
+            fvModels().sourceProxy(alpha1, p_rgh)
-            (
+          + fvModels().sourceProxy(alpha2, p_rgh)
-                volScalarField::New
-                (
-                    "1",
-                    mesh,
-                    dimensionedScalar(dimless/dimPressure, 1)
-                ),
-                p_rgh
-            )
         );
 
         while (pimple.correctNonOrthogonal())
@@ -103,7 +95,7 @@ void Foam::solvers::twoPhaseSolver::incompressiblePressureCorrector
             fvScalarMatrix p_rghEqn
             (
                 fvc::div(phiHbyA) - fvm::laplacian(rAUf, p_rgh)
-             == Sp_rgh
+             == p_rghEqnSource
             );
 
             p_rghEqn.setReference

etc/codeTemplates/dynamicCode/codedFvModelTemplate.C

@@ -127,8 +127,8 @@ ${typeName}FvModel${SourceType}::
 
 void ${typeName}FvModel${SourceType}::addSup
 (
-    fvMatrix<${TemplateType}>& eqn,
+    const VolField<${TemplateType}>& field,
-    const word& fieldName
+    fvMatrix<${TemplateType}>& eqn
 ) const
 {
     if (${verbose})
@@ -145,8 +145,8 @@ void ${typeName}FvModel${SourceType}::addSup
 void ${typeName}FvModel${SourceType}::addSup
 (
     const volScalarField& rho,
-    fvMatrix<${TemplateType}>& eqn,
+    const VolField<${TemplateType}>& field,
-    const word& fieldName
+    fvMatrix<${TemplateType}>& eqn
 ) const
 {
     if (${verbose})
@@ -164,8 +164,8 @@ void ${typeName}FvModel${SourceType}::addSup
 (
     const volScalarField& alpha,
     const volScalarField& rho,
-    fvMatrix<${TemplateType}>& eqn,
+    const VolField<${TemplateType}>& field,
-    const word& fieldName
+    fvMatrix<${TemplateType}>& eqn
 ) const
 {
     if (${verbose})

etc/codeTemplates/dynamicCode/codedFvModelTemplate.H

@@ -89,8 +89,8 @@ public:
             //- Explicit and implicit matrix contributions
             virtual void addSup
             (
-                fvMatrix<${TemplateType}>& eqn,
+                const VolField<${TemplateType}>& field,
-                const word& fieldName
+                fvMatrix<${TemplateType}>& eqn
             ) const;
 
             //- Explicit and implicit matrix contributions for compressible
@@ -98,8 +98,8 @@ public:
             virtual void addSup
             (
                 const volScalarField& rho,
-                fvMatrix<${TemplateType}>& eqn,
+                const VolField<${TemplateType}>& field,
-                const word& fieldName
+                fvMatrix<${TemplateType}>& eqn
             ) const;
 
             //- Explicit and implicit matrix contributions for phase equations
@@ -107,8 +107,8 @@ public:
             (
                 const volScalarField& alpha,
                 const volScalarField& rho,
-                fvMatrix<${TemplateType}>& eqn,
+                const VolField<${TemplateType}>& field,
-                const word& fieldName
+                fvMatrix<${TemplateType}>& eqn
             ) const;
 
 

src/finiteVolume/cfdTools/general/fvModels/fvModel.C

@@ -37,11 +37,16 @@ namespace Foam
 
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
+template<class Type>
+void Foam::fvModel::addSupType(fvMatrix<Type>& eqn) const
+{}
+
+
 template<class Type>
 void Foam::fvModel::addSupType
 (
-    fvMatrix<Type>& eqn,
+    const VolField<Type>& field,
-    const word& fieldName
+    fvMatrix<Type>& eqn
 ) const
 {}
 
@@ -50,8 +55,8 @@ template<class Type>
 void Foam::fvModel::addSupType
 (
     const volScalarField& rho,
-    fvMatrix<Type>& eqn,
+    const VolField<Type>& field,
-    const word& fieldName
+    fvMatrix<Type>& eqn
 ) const
 {}
 
@@ -61,8 +66,8 @@ void Foam::fvModel::addSupType
 (
     const volScalarField& alpha,
     const volScalarField& rho,
-    fvMatrix<Type>& eqn,
+    const VolField<Type>& field,
-    const word& fieldName
+    fvMatrix<Type>& eqn
 ) const
 {}
 
@@ -175,13 +180,16 @@ Foam::scalar Foam::fvModel::maxDeltaT() const
 }
 
 
-FOR_ALL_FIELD_TYPES(IMPLEMENT_FV_MODEL_ADD_SUP, fvModel);
+FOR_ALL_FIELD_TYPES(IMPLEMENT_FV_MODEL_ADD_SUP, fvModel)
 
 
-FOR_ALL_FIELD_TYPES(IMPLEMENT_FV_MODEL_ADD_RHO_SUP, fvModel);
+FOR_ALL_FIELD_TYPES(IMPLEMENT_FV_MODEL_ADD_FIELD_SUP, fvModel)
 
 
-FOR_ALL_FIELD_TYPES(IMPLEMENT_FV_MODEL_ADD_ALPHA_RHO_SUP, fvModel);
+FOR_ALL_FIELD_TYPES(IMPLEMENT_FV_MODEL_ADD_RHO_FIELD_SUP, fvModel)
+
+
+FOR_ALL_FIELD_TYPES(IMPLEMENT_FV_MODEL_ADD_ALPHA_RHO_FIELD_SUP, fvModel)
 
 
 void Foam::fvModel::preUpdateMesh()

src/finiteVolume/cfdTools/general/fvModels/fvModel.H

@@ -80,12 +80,16 @@ protected:
 
     // Protected Member Functions
 
+        //- Add a source term to an equation
+        template<class Type>
+        void addSupType(fvMatrix<Type>& eqn) const;
+
         //- Add a source term to an equation
         template<class Type>
         void addSupType
         (
-            fvMatrix<Type>& eqn,
+            const VolField<Type>& field,
-            const word& fieldName
+            fvMatrix<Type>& eqn
         ) const;
 
         //- Add a source term to a compressible equation
@@ -93,8 +97,8 @@ protected:
         void addSupType
         (
             const volScalarField& rho,
-            fvMatrix<Type>& eqn,
+            const VolField<Type>& field,
-            const word& fieldName
+            fvMatrix<Type>& eqn
         ) const;
 
         //- Add a source term to a phase equation
@@ -103,19 +107,17 @@ protected:
         (
             const volScalarField& alpha,
             const volScalarField& rho,
-            fvMatrix<Type>& eqn,
+            const VolField<Type>& field,
-            const word& fieldName
+            fvMatrix<Type>& eqn
         ) const;
 
-
+        //- Return a source for an equation
-        //- Return source for equation with specified name and dimensions
         template<class Type, class ... AlphaRhoFieldTypes>
-        tmp<fvMatrix<Type>> source
+        tmp<fvMatrix<Type>> sourceTerm
         (
-            const VolField<Type>& field,
+            const VolField<Type>& eqnField,
-            const word& fieldName,
             const dimensionSet& ds,
-            const AlphaRhoFieldTypes& ... alphaRhos
+            const AlphaRhoFieldTypes& ... alphaRhoFields
         ) const;
 
 
@@ -145,29 +147,35 @@ public:
     // Static Member Functions
 
         //- Return the dimensions of the matrix of a source term
-        template
+        template<class AlphaRhoFieldType, class ... AlphaRhoFieldTypes>
-        <
+        static dimensionSet sourceDims
-            class Type,
-            class AlphaRhoFieldType,
-            class ... AlphaRhoFieldTypes
-        >
-        inline static dimensionSet sourceDims
         (
-            const VolField<Type>& field,
             const dimensionSet& ds,
-            const AlphaRhoFieldType& alphaRho,
+            const AlphaRhoFieldType& alphaRhoField,
-            const AlphaRhoFieldTypes& ... alphaRhos
+            const AlphaRhoFieldTypes& ... alphaRhoFields
         );
 
         //- Return the dimensions of the matrix of a source term (base
         //  condition for the above)
-        template<class Type>
+        inline static const dimensionSet& sourceDims(const dimensionSet& ds);
-        inline static dimensionSet sourceDims
+
+        //- Return the name of the field associated with a source term
+        template<class AlphaRhoFieldType, class ... AlphaRhoFieldTypes>
+        static const word& fieldName
         (
-            const VolField<Type>& field,
+            const AlphaRhoFieldType& alphaRhoField,
-            const dimensionSet& ds
+            const AlphaRhoFieldTypes& ... alphaRhoFields
         );
 
+        //- Return the name of the field associated with a source term (base
+        //  condition for the above)
+        template<class AlphaRhoFieldType>
+        static const word& fieldName(const AlphaRhoFieldType& alphaRhoField);
+
+        //- Return the name of the field associated with a source term. Special
+        //  overload for volume sources with no associated field.
+        static const word& fieldName();
+
 
     // Constructors
 
@@ -267,13 +275,23 @@ public:
         // Sources
 
             //- Add a source term to an equation
-            FOR_ALL_FIELD_TYPES(DEFINE_FV_MODEL_ADD_SUP);
+            FOR_ALL_FIELD_TYPES(DEFINE_FV_MODEL_ADD_SUP)
+
+            //- Add a source term to an equation
+            FOR_ALL_FIELD_TYPES(DEFINE_FV_MODEL_ADD_FIELD_SUP)
 
             //- Add a source term to a compressible equation
-            FOR_ALL_FIELD_TYPES(DEFINE_FV_MODEL_ADD_RHO_SUP);
+            FOR_ALL_FIELD_TYPES(DEFINE_FV_MODEL_ADD_RHO_FIELD_SUP)
 
             //- Add a source term to a phase equation
-            FOR_ALL_FIELD_TYPES(DEFINE_FV_MODEL_ADD_ALPHA_RHO_SUP);
+            FOR_ALL_FIELD_TYPES(DEFINE_FV_MODEL_ADD_ALPHA_RHO_FIELD_SUP)
+
+            //- Return source for an equation
+            template<class Type>
+            tmp<fvMatrix<Type>> sourceProxy
+            (
+                const VolField<Type>& eqnField
+            ) const;
 
             //- Return source for an equation
             template<class Type>
@@ -282,12 +300,12 @@ public:
                 const VolField<Type>& field
             ) const;
 
-            //- Return source for an equation with a specified name
+            //- Return source for an equation
             template<class Type>
-            tmp<fvMatrix<Type>> source
+            tmp<fvMatrix<Type>> sourceProxy
             (
                 const VolField<Type>& field,
-                const word& fieldName
+                const VolField<Type>& eqnField
             ) const;
 
             //- Return source for a compressible equation
@@ -298,13 +316,13 @@ public:
                 const VolField<Type>& field
             ) const;
 
-            //- Return source for a compressible equation with a specified name
+            //- Return source for a compressible equation
             template<class Type>
-            tmp<fvMatrix<Type>> source
+            tmp<fvMatrix<Type>> sourceProxy
             (
                 const volScalarField& rho,
                 const VolField<Type>& field,
-                const word& fieldName
+                const VolField<Type>& eqnField
             ) const;
 
             //- Return source for a phase equation
@@ -316,14 +334,14 @@ public:
                 const VolField<Type>& field
             ) const;
 
-            //- Return source for a phase equation with a specified name
+            //- Return source for a phase equation
             template<class Type>
-            tmp<fvMatrix<Type>> source
+            tmp<fvMatrix<Type>> sourceProxy
             (
                 const volScalarField& alpha,
                 const volScalarField& rho,
                 const VolField<Type>& field,
-                const word& fieldName
+                const VolField<Type>& eqnField
             ) const;
 
             //- Return source for a phase equation
@@ -360,14 +378,6 @@ public:
                 const VolField<Type>& field
             ) const;
 
-            //- Return source for an equation with a second time derivative
-            template<class Type>
-            tmp<fvMatrix<Type>> d2dt2
-            (
-                const VolField<Type>& field,
-                const word& fieldName
-            ) const;
-
 
         // Mesh changes
 

src/finiteVolume/cfdTools/general/fvModels/fvModelI.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2021 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2021-2023 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -23,6 +23,25 @@ License
 
 \*---------------------------------------------------------------------------*/
 
+#include "fvModel.H"
+
+// * * * * * * * * * * * * Static Member Functions * * * * * * * * * * * * * //
+
+inline const Foam::dimensionSet& Foam::fvModel::sourceDims
+(
+    const dimensionSet& ds
+)
+{
+    return ds;
+}
+
+
+inline const Foam::word& Foam::fvModel::fieldName()
+{
+    return word::null;
+}
+
+
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
 inline const Foam::word& Foam::fvModel::name() const

src/finiteVolume/cfdTools/general/fvModels/fvModelM.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2021 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2021-2023 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -24,60 +24,69 @@ License
 \*---------------------------------------------------------------------------*/
 
 #define DEFINE_FV_MODEL_ADD_SUP(Type, nullArg)                                 \
-    virtual void addSup                                                        \
+    virtual void addSup(fvMatrix<Type>& eqn) const;
-    (                                                                          \
-        fvMatrix<Type>& eqn,                                                   \
-        const word& fieldName                                                  \
-    ) const;
 
 #define IMPLEMENT_FV_MODEL_ADD_SUP(Type, modelType)                            \
-    void Foam::modelType::addSup                                               \
+    void Foam::modelType::addSup(fvMatrix<Type>& eqn) const                    \
-    (                                                                          \
-        fvMatrix<Type>& eqn,                                                   \
-        const word& fieldName                                                  \
-    ) const                                                                    \
     {                                                                          \
-        addSupType(eqn, fieldName);                                            \
+        addSupType(eqn);                                                       \
     }
 
-#define DEFINE_FV_MODEL_ADD_RHO_SUP(Type, nullArg)                             \
+#define DEFINE_FV_MODEL_ADD_FIELD_SUP(Type, nullArg)                           \
+    virtual void addSup                                                        \
+    (                                                                          \
+        const VolField<Type>& field,                                           \
+        fvMatrix<Type>& eqn                                                    \
+    ) const;
+
+#define IMPLEMENT_FV_MODEL_ADD_FIELD_SUP(Type, modelType)                      \
+    void Foam::modelType::addSup                                               \
+    (                                                                          \
+        const VolField<Type>& field,                                           \
+        fvMatrix<Type>& eqn                                                    \
+    ) const                                                                    \
+    {                                                                          \
+        addSupType(field, eqn);                                                \
+    }
+
+#define DEFINE_FV_MODEL_ADD_RHO_FIELD_SUP(Type, nullArg)                       \
     virtual void addSup                                                        \
     (                                                                          \
         const volScalarField& rho,                                             \
-        fvMatrix<Type>& eqn,                                                   \
+        const VolField<Type>& field,                                           \
-        const word& fieldName                                                  \
+        fvMatrix<Type>& eqn                                                    \
     ) const;
 
-#define IMPLEMENT_FV_MODEL_ADD_RHO_SUP(Type, modelType)                        \
+#define IMPLEMENT_FV_MODEL_ADD_RHO_FIELD_SUP(Type, modelType)                  \
     void Foam::modelType::addSup                                               \
     (                                                                          \
         const volScalarField& rho,                                             \
-        fvMatrix<Type>& eqn,                                                   \
+        const VolField<Type>& field,                                           \
-        const word& fieldName                                                  \
+        fvMatrix<Type>& eqn                                                    \
     ) const                                                                    \
     {                                                                          \
-        addSupType(rho, eqn, fieldName);                                       \
+        addSupType(rho, field, eqn);                                           \
     }
 
-#define DEFINE_FV_MODEL_ADD_ALPHA_RHO_SUP(Type, nullArg)                       \
+#define DEFINE_FV_MODEL_ADD_ALPHA_RHO_FIELD_SUP(Type, nullArg)                 \
     virtual void addSup                                                        \
     (                                                                          \
         const volScalarField& alpha,                                           \
         const volScalarField& rho,                                             \
-        fvMatrix<Type>& eqn,                                                   \
+        const VolField<Type>& field,                                           \
-        const word& fieldName                                                  \
+        fvMatrix<Type>& eqn                                                    \
     ) const;
 
-#define IMPLEMENT_FV_MODEL_ADD_ALPHA_RHO_SUP(Type, modelType)                  \
+#define IMPLEMENT_FV_MODEL_ADD_ALPHA_RHO_FIELD_SUP(Type, modelType)            \
     void Foam::modelType::addSup                                               \
     (                                                                          \
         const volScalarField& alpha,                                           \
         const volScalarField& rho,                                             \
-        fvMatrix<Type>& eqn,                                                   \
+        const VolField<Type>& field,                                           \
-        const word& fieldName                                                  \
+        fvMatrix<Type>& eqn                                                    \
     ) const                                                                    \
     {                                                                          \
-        addSupType(alpha, rho, eqn, fieldName);                                \
+        addSupType(alpha, rho, field, eqn);                                    \
     }
 
 

src/finiteVolume/cfdTools/general/fvModels/fvModelTemplates.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2021-2022 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2021-2023 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -23,56 +23,66 @@ License
 
 \*---------------------------------------------------------------------------*/
 
-// * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * * //
+#include "fvModel.H"
 
-template<class Type, class AlphaRhoFieldType, class ... AlphaRhoFieldTypes>
+// * * * * * * * * * * * * Static Member Functions * * * * * * * * * * * * * //
+
+template<class AlphaRhoFieldType, class ... AlphaRhoFieldTypes>
 Foam::dimensionSet Foam::fvModel::sourceDims
 (
-    const VolField<Type>& field,
     const dimensionSet& ds,
-    const AlphaRhoFieldType& alphaRho,
+    const AlphaRhoFieldType& alphaRhoField,
-    const AlphaRhoFieldTypes& ... alphaRhos
+    const AlphaRhoFieldTypes& ... alphaRhoFields
 )
 {
-    return alphaRho.dimensions()*sourceDims(field, ds, alphaRhos ...);
+    return alphaRhoField.dimensions()*sourceDims(ds, alphaRhoFields ...);
 }
 
 
-template<class Type>
+template<class AlphaRhoFieldType, class ... AlphaRhoFieldTypes>
-Foam::dimensionSet Foam::fvModel::sourceDims
+const Foam::word& Foam::fvModel::fieldName
 (
-    const VolField<Type>& field,
+    const AlphaRhoFieldType& alphaRhoField,
-    const dimensionSet& ds
+    const AlphaRhoFieldTypes& ... alphaRhoFields
 )
 {
-    return field.dimensions()*ds;
+    return fieldName(alphaRhoFields ...);
+}
+
+
+template<class AlphaRhoFieldType>
+const Foam::word& Foam::fvModel::fieldName
+(
+    const AlphaRhoFieldType& alphaRhoField
+)
+{
+    return alphaRhoField.name();
 }
 
 
 // * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * //
 
 template<class Type, class ... AlphaRhoFieldTypes>
-Foam::tmp<Foam::fvMatrix<Type>> Foam::fvModel::source
+Foam::tmp<Foam::fvMatrix<Type>> Foam::fvModel::sourceTerm
 (
-    const VolField<Type>& field,
+    const VolField<Type>& eqnField,
-    const word& fieldName,
     const dimensionSet& ds,
-    const AlphaRhoFieldTypes& ... alphaRhos
+    const AlphaRhoFieldTypes& ... alphaRhoFields
 ) const
 {
     tmp<fvMatrix<Type>> tmtx
     (
         new fvMatrix<Type>
         (
-            field,
+            eqnField,
-            sourceDims(field, ds, alphaRhos ...)
+            sourceDims(ds, alphaRhoFields ...)
         )
     );
     fvMatrix<Type>& mtx = tmtx.ref();
 
-    if (addsSupToField(fieldName))
+    if (addsSupToField(fieldName(alphaRhoFields ...)))
     {
-        addSup(alphaRhos ..., mtx, fieldName);
+        addSup(alphaRhoFields ..., mtx);
     }
 
     return tmtx;
@@ -82,23 +92,33 @@ Foam::tmp<Foam::fvMatrix<Type>> Foam::fvModel::source
 // * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * * //
 
 template<class Type>
-Foam::tmp<Foam::fvMatrix<Type>> Foam::fvModel::source
+Foam::tmp<Foam::fvMatrix<Type>> Foam::fvModel::sourceProxy
 (
-    const VolField<Type>& field
+    const VolField<Type>& eqnField
 ) const
 {
-    return this->source(field, field.name());
+    return sourceTerm(eqnField, dimVolume/dimTime);
 }
 
 
 template<class Type>
 Foam::tmp<Foam::fvMatrix<Type>> Foam::fvModel::source
 (
-    const VolField<Type>& field,
+    const VolField<Type>& field
-    const word& fieldName
 ) const
 {
-    return source(field, fieldName, dimVolume/dimTime);
+    return sourceTerm(field, dimVolume/dimTime, field);
+}
+
+
+template<class Type>
+Foam::tmp<Foam::fvMatrix<Type>> Foam::fvModel::sourceProxy
+(
+    const VolField<Type>& field,
+    const VolField<Type>& eqnField
+) const
+{
+    return sourceTerm(eqnField, dimVolume/dimTime, field);
 }
 
 
@@ -109,19 +129,19 @@ Foam::tmp<Foam::fvMatrix<Type>> Foam::fvModel::source
     const VolField<Type>& field
 ) const
 {
-    return this->source(rho, field, field.name());
+    return sourceTerm(field, dimVolume/dimTime, rho, field);
 }
 
 
 template<class Type>
-Foam::tmp<Foam::fvMatrix<Type>> Foam::fvModel::source
+Foam::tmp<Foam::fvMatrix<Type>> Foam::fvModel::sourceProxy
 (
     const volScalarField& rho,
     const VolField<Type>& field,
-    const word& fieldName
+    const VolField<Type>& eqnField
 ) const
 {
-    return source(field, fieldName, dimVolume/dimTime, rho);
+    return sourceTerm(eqnField, dimVolume/dimTime, rho, field);
 }
 
 
@@ -133,20 +153,20 @@ Foam::tmp<Foam::fvMatrix<Type>> Foam::fvModel::source
     const VolField<Type>& field
 ) const
 {
-    return this->source(alpha, rho, field, field.name());
+    return sourceTerm(field, dimVolume/dimTime, alpha, rho, field);
 }
 
 
 template<class Type>
-Foam::tmp<Foam::fvMatrix<Type>> Foam::fvModel::source
+Foam::tmp<Foam::fvMatrix<Type>> Foam::fvModel::sourceProxy
 (
     const volScalarField& alpha,
     const volScalarField& rho,
     const VolField<Type>& field,
-    const word& fieldName
+    const VolField<Type>& eqnField
 ) const
 {
-    return source(field, fieldName, dimVolume/dimTime, alpha, rho);
+    return sourceTerm(eqnField, dimVolume/dimTime, alpha, rho, field);
 }
 
 
@@ -158,7 +178,7 @@ Foam::tmp<Foam::fvMatrix<Type>> Foam::fvModel::source
     const VolField<Type>& field
 ) const
 {
-    return this->source(field, field.name());
+    return this->source(field);
 }
 
 
@@ -185,7 +205,7 @@ Foam::tmp<Foam::fvMatrix<Type>> Foam::fvModel::source
         dimensionedScalar(dimless, 1.0)
     );
 
-    return this->source(alpha, one, field, field.name());
+    return this->source(alpha, one, field);
 }
 
 
@@ -197,30 +217,18 @@ Foam::tmp<Foam::fvMatrix<Type>> Foam::fvModel::source
     const VolField<Type>& field
 ) const
 {
-    return this->source(rho, field, field.name());
+    return this->source(rho, field);
 }
 
 
-// * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
-
 template<class Type>
 Foam::tmp<Foam::fvMatrix<Type>> Foam::fvModel::d2dt2
 (
     const VolField<Type>& field
 ) const
 {
-    return this->d2dt2(field, field.name());
+    return sourceTerm(field, dimVolume/sqr(dimTime), field);
 }
 
 
-template<class Type>
-Foam::tmp<Foam::fvMatrix<Type>> Foam::fvModel::d2dt2
-(
-    const VolField<Type>& field,
-    const word& fieldName
-) const
-{
-    return source(field, fieldName, dimVolume/sqr(dimTime));
-}
-
 // ************************************************************************* //

src/finiteVolume/cfdTools/general/fvModels/fvModels.H

@@ -81,14 +81,13 @@ class fvModels
         //- Check that all fvModels have been applied
         void checkApplied() const;
 
-        //- Return source for equation with specified name and dimensions
+        //- Return a source for an equation
         template<class Type, class ... AlphaRhoFieldTypes>
-        tmp<fvMatrix<Type>> source
+        tmp<fvMatrix<Type>> sourceTerm
         (
-            const VolField<Type>& field,
+            const VolField<Type>& eqnField,
-            const word& fieldName,
             const dimensionSet& ds,
-            const AlphaRhoFieldTypes& ... alphaRhos
+            const AlphaRhoFieldTypes& ... alphaRhoFields
         ) const;
 
 
@@ -153,6 +152,13 @@ public:
             //  e.g. solve equations, update model, for film, Lagrangian etc.
             virtual void correct();
 
+            //- Return source for an equation
+            template<class Type>
+            tmp<fvMatrix<Type>> sourceProxy
+            (
+                const VolField<Type>& eqnField
+            ) const;
+
             //- Return source for an equation
             template<class Type>
             tmp<fvMatrix<Type>> source
@@ -160,12 +166,12 @@ public:
                 const VolField<Type>& field
             ) const;
 
-            //- Return source for an equation with a specified name
+            //- Return source for an equation
             template<class Type>
-            tmp<fvMatrix<Type>> source
+            tmp<fvMatrix<Type>> sourceProxy
             (
                 const VolField<Type>& field,
-                const word& fieldName
+                const VolField<Type>& eqnField
             ) const;
 
             //- Return source for a compressible equation
@@ -176,13 +182,13 @@ public:
                 const VolField<Type>& field
             ) const;
 
-            //- Return source for a compressible equation with a specified name
+            //- Return source for a compressible equation
             template<class Type>
-            tmp<fvMatrix<Type>> source
+            tmp<fvMatrix<Type>> sourceProxy
             (
                 const volScalarField& rho,
                 const VolField<Type>& field,
-                const word& fieldName
+                const VolField<Type>& eqnField
             ) const;
 
             //- Return source for a phase equation
@@ -194,14 +200,14 @@ public:
                 const VolField<Type>& field
             ) const;
 
-            //- Return source for a phase equation with a specified name
+            //- Return source for a phase equation
             template<class Type>
-            tmp<fvMatrix<Type>> source
+            tmp<fvMatrix<Type>> sourceProxy
             (
                 const volScalarField& alpha,
                 const volScalarField& rho,
                 const VolField<Type>& field,
-                const word& fieldName
+                const VolField<Type>& eqnField
             ) const;
 
             //- Return source for a phase equation
@@ -238,14 +244,6 @@ public:
                 const VolField<Type>& field
             ) const;
 
-            //- Return source for an equation with a second time derivative
-            template<class Type>
-            tmp<fvMatrix<Type>> d2dt2
-            (
-                const VolField<Type>& field,
-                const word& fieldName
-            ) const;
-
 
         // Mesh changes
 

src/finiteVolume/cfdTools/general/fvModels/fvModelsTemplates.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2021-2022 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2021-2023 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -23,15 +23,16 @@ License
 
 \*---------------------------------------------------------------------------*/
 
+#include "fvModels.H"
+
 // * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * //
 
 template<class Type, class ... AlphaRhoFieldTypes>
-Foam::tmp<Foam::fvMatrix<Type>> Foam::fvModels::source
+Foam::tmp<Foam::fvMatrix<Type>> Foam::fvModels::sourceTerm
 (
-    const VolField<Type>& field,
+    const VolField<Type>& eqnField,
-    const word& fieldName,
     const dimensionSet& ds,
-    const AlphaRhoFieldTypes& ... alphaRhos
+    const AlphaRhoFieldTypes& ... alphaRhoFields
 ) const
 {
     checkApplied();
@@ -40,14 +41,16 @@ Foam::tmp<Foam::fvMatrix<Type>> Foam::fvModels::source
     (
         new fvMatrix<Type>
         (
-            field,
+            eqnField,
-            fvModel::sourceDims(field, ds, alphaRhos ...)
+            fvModel::sourceDims(ds, alphaRhoFields ...)
         )
     );
     fvMatrix<Type>& mtx = tmtx.ref();
 
     const PtrListDictionary<fvModel>& modelList(*this);
 
+    const word& fieldName = fvModel::fieldName(alphaRhoFields ...);
+
     forAll(modelList, i)
     {
         const fvModel& model = modelList[i];
@@ -62,7 +65,7 @@ Foam::tmp<Foam::fvMatrix<Type>> Foam::fvModels::source
                     << fieldName << endl;
             }
 
-            model.addSup(alphaRhos ..., mtx, fieldName);
+            model.addSup(alphaRhoFields ..., mtx);
         }
     }
 
@@ -73,23 +76,33 @@ Foam::tmp<Foam::fvMatrix<Type>> Foam::fvModels::source
 // * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * * //
 
 template<class Type>
-Foam::tmp<Foam::fvMatrix<Type>> Foam::fvModels::source
+Foam::tmp<Foam::fvMatrix<Type>> Foam::fvModels::sourceProxy
 (
-    const VolField<Type>& field
+    const VolField<Type>& eqnField
 ) const
 {
-    return this->source(field, field.name());
+    return sourceTerm(eqnField, dimVolume/dimTime);
 }
 
 
 template<class Type>
 Foam::tmp<Foam::fvMatrix<Type>> Foam::fvModels::source
 (
-    const VolField<Type>& field,
+    const VolField<Type>& field
-    const word& fieldName
 ) const
 {
-    return source(field, fieldName, dimVolume/dimTime);
+    return sourceTerm(field, dimVolume/dimTime, field);
+}
+
+
+template<class Type>
+Foam::tmp<Foam::fvMatrix<Type>> Foam::fvModels::sourceProxy
+(
+    const VolField<Type>& field,
+    const VolField<Type>& eqnField
+) const
+{
+    return sourceTerm(eqnField, dimVolume/dimTime, field);
 }
 
 
@@ -100,19 +113,19 @@ Foam::tmp<Foam::fvMatrix<Type>> Foam::fvModels::source
     const VolField<Type>& field
 ) const
 {
-    return this->source(rho, field, field.name());
+    return sourceTerm(field, dimVolume/dimTime, rho, field);
 }
 
 
 template<class Type>
-Foam::tmp<Foam::fvMatrix<Type>> Foam::fvModels::source
+Foam::tmp<Foam::fvMatrix<Type>> Foam::fvModels::sourceProxy
 (
     const volScalarField& rho,
     const VolField<Type>& field,
-    const word& fieldName
+    const VolField<Type>& eqnField
 ) const
 {
-    return source(field, fieldName, dimVolume/dimTime, rho);
+    return sourceTerm(eqnField, dimVolume/dimTime, rho, field);
 }
 
 
@@ -124,20 +137,20 @@ Foam::tmp<Foam::fvMatrix<Type>> Foam::fvModels::source
     const VolField<Type>& field
 ) const
 {
-    return this->source(alpha, rho, field, field.name());
+    return sourceTerm(field, dimVolume/dimTime, alpha, rho, field);
 }
 
 
 template<class Type>
-Foam::tmp<Foam::fvMatrix<Type>> Foam::fvModels::source
+Foam::tmp<Foam::fvMatrix<Type>> Foam::fvModels::sourceProxy
 (
     const volScalarField& alpha,
     const volScalarField& rho,
     const VolField<Type>& field,
-    const word& fieldName
+    const VolField<Type>& eqnField
 ) const
 {
-    return source(field, fieldName, dimVolume/dimTime, alpha, rho);
+    return sourceTerm(eqnField, dimVolume/dimTime, alpha, rho, field);
 }
 
 
@@ -149,7 +162,7 @@ Foam::tmp<Foam::fvMatrix<Type>> Foam::fvModels::source
     const VolField<Type>& field
 ) const
 {
-    return this->source(field, field.name());
+    return this->source(field);
 }
 
 
@@ -176,7 +189,7 @@ Foam::tmp<Foam::fvMatrix<Type>> Foam::fvModels::source
         dimensionedScalar(dimless, 1.0)
     );
 
-    return this->source(alpha, one, field, field.name());
+    return this->source(alpha, one, field);
 }
 
 
@@ -188,30 +201,17 @@ Foam::tmp<Foam::fvMatrix<Type>> Foam::fvModels::source
     const VolField<Type>& field
 ) const
 {
-    return this->source(rho, field, field.name());
+    return this->source(rho, field);
 }
 
 
-// * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
-
 template<class Type>
 Foam::tmp<Foam::fvMatrix<Type>> Foam::fvModels::d2dt2
 (
     const VolField<Type>& field
 ) const
 {
-    return this->d2dt2(field, field.name());
+    return sourceTerm(field, dimVolume/sqr(dimTime), field);
-}
-
-
-template<class Type>
-Foam::tmp<Foam::fvMatrix<Type>> Foam::fvModels::d2dt2
-(
-    const VolField<Type>& field,
-    const word& fieldName
-) const
-{
-    return source(field, fieldName, dimVolume/sqr(dimTime));
 }
 
 

src/finiteVolume/fvMatrices/fvMatrix/fvMatrix.C

@@ -854,7 +854,9 @@ Foam::tmp<Foam::volScalarField::Internal> Foam::fvMatrix<Type>::Sp() const
             "Sp(" + psi_.name() + ')',
             psi_.mesh(),
             dimensions_/psi_.dimensions()/dimVol,
-            diag()/psi_.mesh().V()
+            hasDiag()
+          ? diag()/psi_.mesh().V()
+          : tmp<scalarField>(new scalarField(lduAddr().size(), scalar(0)))
         )
     );
 

src/fvConstraints/zeroDimensionalFixedPressure/zeroDimensionalFixedPressureConstraint.C

@@ -168,37 +168,46 @@ Foam::fv::zeroDimensionalFixedPressureConstraint::constrainedFields() const
 }
 
 
-const Foam::volScalarField::Internal&
+Foam::tmp<Foam::volScalarField::Internal>
 Foam::fv::zeroDimensionalFixedPressureConstraint::pEqnSource
 (
+    const volScalarField& rho,
     fvMatrix<scalar>& pEqn
 ) const
 {
     // Ensure the corresponding fvModel exits
     model();
 
-    // Construct the source if it does not yet exist
+    // Return zero if the source does not yet exist
     if (!sourcePtr_.valid())
     {
-        sourcePtr_.set
+        return
-        (
+            volScalarField::Internal::New
-            new volScalarField::Internal
             (
-                IOobject
+                typedName("source"),
-                (
-                    typedName("source"),
-                    mesh().time().timeName(),
-                    mesh(),
-                    IOobject::READ_IF_PRESENT,
-                    IOobject::AUTO_WRITE
-                ),
                 mesh(),
                 dimensionedScalar(pEqn.dimensions()/dimVolume, 0)
-            )
+            );
-        );
     }
 
-    return sourcePtr_();
+    // Return the source, multiplying by density if needed
+    if (sourcePtr_->dimensions() == pEqn.dimensions()/dimVolume)
+    {
+        return sourcePtr_();
+    }
+    else if (sourcePtr_->dimensions() == pEqn.dimensions()/dimMass)
+    {
+        return rho()*sourcePtr_();
+    }
+    else
+    {
+        FatalErrorInFunction
+            << "Dimensions of equation for pressure "
+            << pEqn.psi().name() << " not recognised"
+            << exit(FatalError);
+
+        return tmp<volScalarField::Internal>(nullptr);
+    }
 }
 
 
@@ -229,11 +238,26 @@ bool Foam::fv::zeroDimensionalFixedPressureConstraint::constrain
     const word& fieldName
 ) const
 {
-    // Construct the source if it does not yet exist
+    // Create the source field if it does not already exist
-    pEqnSource(pEqn);
+    if (!sourcePtr_.valid())
-
+    {
-    // Check the dimensions have not changed
+        sourcePtr_.set
-    sourcePtr_->dimensions() = pEqn.dimensions()/dimVolume;
+        (
+            new volScalarField::Internal
+            (
+                IOobject
+                (
+                    typedName("source"),
+                    mesh().time().timeName(),
+                    mesh(),
+                    IOobject::READ_IF_PRESENT,
+                    IOobject::AUTO_WRITE
+                ),
+                mesh(),
+                dimensionedScalar(pEqn.dimensions()/dimVolume, 0)
+            )
+        );
+    }
 
     // Remove the previous iteration's source from the pressure equation
     pEqn += sourcePtr_();

src/fvConstraints/zeroDimensionalFixedPressure/zeroDimensionalFixedPressureConstraint.H

@@ -160,8 +160,9 @@ public:
         // Constraints
 
             //- Return the mass or volume source for the pressure equation
-            const volScalarField::Internal& pEqnSource
+            tmp<volScalarField::Internal> pEqnSource
             (
+                const volScalarField& rho,
                 fvMatrix<scalar>& pEqn
             ) const;
 

src/fvConstraints/zeroDimensionalFixedPressure/zeroDimensionalFixedPressureModel.C

@@ -76,12 +76,12 @@ Foam::fv::zeroDimensionalFixedPressureModel::constraint() const
 template<class Type>
 void Foam::fv::zeroDimensionalFixedPressureModel::addSupType
 (
-    fvMatrix<Type>& eqn,
+    const VolField<Type>& field,
-    const word& fieldName
+    fvMatrix<Type>& eqn
 ) const
 {
     FatalErrorInFunction
-        << "Cannot add a fixed pressure source to field " << fieldName
+        << "Cannot add a fixed pressure source of field " << field.name()
         << " because this field's equation is not in mass-conservative form"
         << exit(FatalError);
 }
@@ -89,17 +89,26 @@ void Foam::fv::zeroDimensionalFixedPressureModel::addSupType
 
 void Foam::fv::zeroDimensionalFixedPressureModel::addSupType
 (
-    fvMatrix<scalar>& eqn,
+    const volScalarField& rho,
-    const word& fieldName
+    fvMatrix<scalar>& eqn
 ) const
 {
-    if (IOobject::member(fieldName) == constraint().rhoName())
+    if (IOobject::member(rho.name()) == constraint().rhoName())
     {
-        eqn += constraint().massSource(eqn.psi()());
+        if (IOobject::member(eqn.psi().name()) == constraint().pName())
+        {
+            eqn += constraint().pEqnSource(rho, eqn);
+        }
+        else
+        {
+            eqn += constraint().massSource(rho());
+        }
     }
     else
     {
-        addSupType<scalar>(eqn, fieldName); // error above
+        // This is actually an incompressible single-phase equation. Rho is
+        // actually a property field. Fall back (and error).
+        addSupType<scalar>(rho, eqn);
     }
 }
 
@@ -108,46 +117,45 @@ template<class Type>
 void Foam::fv::zeroDimensionalFixedPressureModel::addSupType
 (
     const volScalarField& rho,
-    fvMatrix<Type>& eqn,
+    const VolField<Type>& field,
-    const word& fieldName
+    fvMatrix<Type>& eqn
 ) const
 {
+    if (&field != &eqn.psi())
+    {
+        FatalErrorInFunction
+            << "Cannot add a fixed pressure source of field " << field.name()
+            << " into an equation for field " << eqn.psi().name()
+            << exit(FatalError);
+    }
+
     eqn -= fvm::SuSp(-constraint().massSource(rho()), eqn.psi());
 }
 
 
 void Foam::fv::zeroDimensionalFixedPressureModel::addSupType
 (
+    const volScalarField& alpha,
     const volScalarField& rho,
-    fvMatrix<scalar>& eqn,
+    fvMatrix<scalar>& eqn
-    const word& fieldName
 ) const
 {
-    if (IOobject::member(fieldName) == constraint().rhoName())
+    if (IOobject::member(rho.name()) == constraint().rhoName())
     {
         if (IOobject::member(eqn.psi().name()) == constraint().pName())
         {
-            eqn += constraint().pEqnSource(eqn);
+            eqn += alpha()*constraint().pEqnSource(rho, eqn);
-        }
-        else if (IOobject::member(eqn.psi().name()) == constraint().rhoName())
-        {
-            // Phase density equation. Argument names are misleading.
-            const volScalarField& alpha = rho;
-            const volScalarField& rho = eqn.psi();
-
-            eqn += constraint().massSource(alpha(), rho());
         }
         else
         {
-            FatalErrorInFunction
+            eqn += constraint().massSource(alpha(), rho());
-                << "Cannot add source for density field " << fieldName
-                << " into an equation for " << eqn.psi().name()
-                << exit(FatalError);
         }
     }
     else
     {
-        addSupType<scalar>(rho, eqn, fieldName);
+        // This is actually a compressible single-phase equation. Alpha is
+        // actually rho, and rho is actually a property field. Fall back.
+        addSupType<scalar>(alpha, rho, eqn);
     }
 }
 
@@ -157,50 +165,22 @@ void Foam::fv::zeroDimensionalFixedPressureModel::addSupType
 (
     const volScalarField& alpha,
     const volScalarField& rho,
-    fvMatrix<Type>& eqn,
+    const VolField<Type>& field,
-    const word& fieldName
+    fvMatrix<Type>& eqn
 ) const
 {
+    if (&field != &eqn.psi())
+    {
+        FatalErrorInFunction
+            << "Cannot add a fixed pressure source of field " << field.name()
+            << " into an equation for field " << eqn.psi().name()
+            << exit(FatalError);
+    }
+
     eqn -= fvm::SuSp(-constraint().massSource(alpha(), rho()), eqn.psi());
 }
 
 
-void Foam::fv::zeroDimensionalFixedPressureModel::addSupType
-(
-    const volScalarField& alpha,
-    const volScalarField& rho,
-    fvMatrix<scalar>& eqn,
-    const word& fieldName
-) const
-{
-    if (IOobject::member(fieldName) == constraint().rhoName())
-    {
-        if (IOobject::member(eqn.psi().name()) == constraint().pName())
-        {
-            eqn += alpha*constraint().pEqnSource(eqn);
-        }
-        else if (IOobject::member(eqn.psi().name()) == constraint().rhoName())
-        {
-            FatalErrorInFunction
-                << "Cannot add source for density field " << fieldName
-                << " into a phase-conservative equation for "
-                << eqn.psi().name() << exit(FatalError);
-        }
-        else
-        {
-            FatalErrorInFunction
-                << "Cannot add source for density field " << fieldName
-                << " into an equation for " << eqn.psi().name()
-                << exit(FatalError);
-        }
-    }
-    else
-    {
-        addSupType<scalar>(alpha, rho, eqn, fieldName);
-    }
-}
-
-
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 Foam::fv::zeroDimensionalFixedPressureModel::zeroDimensionalFixedPressureModel
@@ -243,23 +223,23 @@ bool Foam::fv::zeroDimensionalFixedPressureModel::addsSupToField
 
 FOR_ALL_FIELD_TYPES
 (
-    IMPLEMENT_FV_MODEL_ADD_SUP,
+    IMPLEMENT_FV_MODEL_ADD_FIELD_SUP,
     fv::zeroDimensionalFixedPressureModel
-);
+)
 
 
 FOR_ALL_FIELD_TYPES
 (
-    IMPLEMENT_FV_MODEL_ADD_RHO_SUP,
+    IMPLEMENT_FV_MODEL_ADD_RHO_FIELD_SUP,
     fv::zeroDimensionalFixedPressureModel
-);
+)
 
 
 FOR_ALL_FIELD_TYPES
 (
-    IMPLEMENT_FV_MODEL_ADD_ALPHA_RHO_SUP,
+    IMPLEMENT_FV_MODEL_ADD_ALPHA_RHO_FIELD_SUP,
     fv::zeroDimensionalFixedPressureModel
-);
+)
 
 
 bool Foam::fv::zeroDimensionalFixedPressureModel::movePoints()

src/fvConstraints/zeroDimensionalFixedPressure/zeroDimensionalFixedPressureModel.H

@@ -79,26 +79,34 @@ class zeroDimensionalFixedPressureModel
 
             //- Add a source term to an equation
             template<class Type>
-            void addSupType(fvMatrix<Type>& eqn, const word& fieldName) const;
+            void addSupType
+            (
+                const VolField<Type>& field,
+                fvMatrix<Type>& eqn
+            ) const;
 
-            //- Add a source term to a scalar equation
+            //- Add a source term to a compressible continuity equation
-            void addSupType(fvMatrix<scalar>& eqn, const word& fieldName) const;
+            void addSupType
+            (
+                const volScalarField& rho,
+                fvMatrix<scalar>& eqn
+            ) const;
 
             //- Add a source term to a compressible equation
             template<class Type>
             void addSupType
             (
                 const volScalarField& rho,
-                fvMatrix<Type>& eqn,
+                const VolField<Type>& field,
-                const word& fieldName
+                fvMatrix<Type>& eqn
             ) const;
 
-            //- Add a source term to a scalar compressible equation
+            //- Add a source term to a phase continuity equation
             void addSupType
             (
+                const volScalarField& alpha,
                 const volScalarField& rho,
-                fvMatrix<scalar>& eqn,
+                fvMatrix<scalar>& eqn
-                const word& fieldName
             ) const;
 
             //- Add a source term to a phase equation
@@ -107,17 +115,8 @@ class zeroDimensionalFixedPressureModel
             (
                 const volScalarField& alpha,
                 const volScalarField& rho,
-                fvMatrix<Type>& eqn,
+                const VolField<Type>& field,
-                const word& fieldName
+                fvMatrix<Type>& eqn
-            ) const;
-
-            //- Add a source term to a scalar phase equation
-            void addSupType
-            (
-                const volScalarField& alpha,
-                const volScalarField& rho,
-                fvMatrix<scalar>& eqn,
-                const word& fieldName
             ) const;
 
 
@@ -155,13 +154,13 @@ public:
         // Sources
 
             //- Add a source term to an equation
-            FOR_ALL_FIELD_TYPES(DEFINE_FV_MODEL_ADD_SUP);
+            FOR_ALL_FIELD_TYPES(DEFINE_FV_MODEL_ADD_FIELD_SUP)
 
             //- Add a source term to a compressible equation
-            FOR_ALL_FIELD_TYPES(DEFINE_FV_MODEL_ADD_RHO_SUP);
+            FOR_ALL_FIELD_TYPES(DEFINE_FV_MODEL_ADD_RHO_FIELD_SUP)
 
             //- Add a source term to a phase equation
-            FOR_ALL_FIELD_TYPES(DEFINE_FV_MODEL_ADD_ALPHA_RHO_SUP);
+            FOR_ALL_FIELD_TYPES(DEFINE_FV_MODEL_ADD_ALPHA_RHO_FIELD_SUP)
 
 
         // Mesh changes

src/fvModels/Make/files

@@ -23,6 +23,7 @@ derived/phaseLimitStabilisation/phaseLimitStabilisation.C
 derived/accelerationSource/accelerationSource.C
 derived/volumeFractionSource/volumeFractionSource.C
 derived/solidEquilibriumEnergySource/solidEquilibriumEnergySource.C
+derived/volumeSource/volumeSource.C
 derived/massSource/massSource.C
 derived/heatSource/heatSource.C
 derived/heatTransfer/heatTransfer.C

src/fvModels/derived/accelerationSource/accelerationSource.C

@@ -51,6 +51,32 @@ void Foam::fv::accelerationSource::readCoeffs()
 }
 
 
+template<class AlphaRhoFieldType>
+void Foam::fv::accelerationSource::add
+(
+    const AlphaRhoFieldType& alphaRho,
+    fvMatrix<vector>& eqn
+) const
+{
+    const DimensionedField<scalar, volMesh>& V = mesh().V();
+
+    const scalar t = mesh().time().value();
+    const scalar dt = mesh().time().deltaTValue();
+    const vector dU =
+        velocity_->value(mesh().time().timeToUserTime(t))
+      - velocity_->value(mesh().time().timeToUserTime(t - dt));
+    const vector a = dU/mesh().time().deltaTValue();
+
+    const labelUList cells = set_.cells();
+
+    forAll(cells, i)
+    {
+        const label celli = cells[i];
+        eqn.source()[celli] -= V[celli]*alphaRho[celli]*a;
+    }
+}
+
+
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 Foam::fv::accelerationSource::accelerationSource
@@ -80,22 +106,22 @@ Foam::wordList Foam::fv::accelerationSource::addSupFields() const
 
 void Foam::fv::accelerationSource::addSup
 (
-    fvMatrix<vector>& eqn,
+    const volVectorField& U,
-    const word& fieldName
+    fvMatrix<vector>& eqn
 ) const
 {
-    add(geometricOneField(), eqn, fieldName);
+    add(geometricOneField(), eqn);
 }
 
 
 void Foam::fv::accelerationSource::addSup
 (
     const volScalarField& rho,
-    fvMatrix<vector>& eqn,
+    const volVectorField& U,
-    const word& fieldName
+    fvMatrix<vector>& eqn
 ) const
 {
-    add(rho, eqn, fieldName);
+    add(rho, eqn);
 }
 
 
@@ -103,11 +129,11 @@ void Foam::fv::accelerationSource::addSup
 (
     const volScalarField& alpha,
     const volScalarField& rho,
-    fvMatrix<vector>& eqn,
+    const volVectorField& U,
-    const word& fieldName
+    fvMatrix<vector>& eqn
 ) const
 {
-    add((alpha*rho)(), eqn, fieldName);
+    add((alpha*rho)(), eqn);
 }
 
 

src/fvModels/derived/accelerationSource/accelerationSource.H

@@ -96,11 +96,10 @@ class accelerationSource
 
         //- Source term to momentum equation
         template<class AlphaRhoFieldType>
-        void add
+        inline void add
         (
             const AlphaRhoFieldType& rho,
-            fvMatrix<vector>& eqn,
+            fvMatrix<vector>& eqn
-            const word& fieldName
         ) const;
 
 
@@ -141,16 +140,16 @@ public:
             //- Source term to momentum equation
             virtual void addSup
             (
-                fvMatrix<vector>& eqn,
+                const volVectorField& U,
-                const word& fieldName
+                fvMatrix<vector>& eqn
             ) const;
 
             //- Source term to compressible momentum equation
             virtual void addSup
             (
                 const volScalarField& rho,
-                fvMatrix<vector>& eqn,
+                const volVectorField& U,
-                const word& fieldName
+                fvMatrix<vector>& eqn
             ) const;
 
             //- Source term to phase momentum equation
@@ -158,8 +157,8 @@ public:
             (
                 const volScalarField& alpha,
                 const volScalarField& rho,
-                fvMatrix<vector>& eqn,
+                const volVectorField& U,
-                const word& fieldName
+                fvMatrix<vector>& eqn
             ) const;
 
 
@@ -192,12 +191,6 @@ public:
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-#ifdef NoRepository
-    #include "accelerationSourceTemplates.C"
-#endif
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
 #endif
 
 // ************************************************************************* //

src/fvModels/derived/accelerationSource/accelerationSourceTemplates.C

@@ -1,55 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2018-2023 OpenFOAM Foundation
-     \\/     M anipulation  |
--------------------------------------------------------------------------------
-License
-    This file is part of OpenFOAM.
-
-    OpenFOAM is free software: you can redistribute it and/or modify it
-    under the terms of the GNU General Public License as published by
-    the Free Software Foundation, either version 3 of the License, or
-    (at your option) any later version.
-
-    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
-    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
-    for more details.
-
-    You should have received a copy of the GNU General Public License
-    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
-
-\*---------------------------------------------------------------------------*/
-
-// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
-
-template<class AlphaRhoFieldType>
-void Foam::fv::accelerationSource::add
-(
-    const AlphaRhoFieldType& alphaRho,
-    fvMatrix<vector>& eqn,
-    const word& fieldName
-) const
-{
-    const DimensionedField<scalar, volMesh>& V = mesh().V();
-
-    const scalar t = mesh().time().value();
-    const scalar dt = mesh().time().deltaTValue();
-    const vector dU =
-        velocity_->value(mesh().time().timeToUserTime(t))
-      - velocity_->value(mesh().time().timeToUserTime(t - dt));
-    const vector a = dU/mesh().time().deltaTValue();
-
-    const labelUList cells = set_.cells();
-
-    forAll(cells, i)
-    {
-        const label celli = cells[i];
-        eqn.source()[celli] -= V[celli]*alphaRho[celli]*a;
-    }
-}
-
-
-// ************************************************************************* //

src/fvModels/derived/actuationDiskSource/actuationDiskSource.C

@@ -95,6 +95,37 @@ void Foam::fv::actuationDiskSource::readCoeffs()
 }
 
 
+template<class AlphaFieldType, class RhoFieldType>
+void Foam::fv::actuationDiskSource::addActuationDiskAxialInertialResistance
+(
+    vectorField& Usource,
+    const labelList& cells,
+    const scalarField& Vcells,
+    const AlphaFieldType& alpha,
+    const RhoFieldType& rho,
+    const vectorField& U
+) const
+{
+    const scalar a = 1 - Cp_/Ct_;
+    const vector dHat(diskDir_/mag(diskDir_));
+
+    scalar dHatUo(vGreat);
+    if (upstreamCellId_ != -1)
+    {
+        dHatUo = dHat & U[upstreamCellId_];
+    }
+    reduce(dHatUo, minOp<scalar>());
+
+    const vector T = 2*diskArea_*sqr(dHatUo)*a*(1 - a)*dHat;
+
+    forAll(cells, i)
+    {
+        Usource[cells[i]] +=
+            (alpha[cells[i]]*rho[cells[i]]*(Vcells[cells[i]]/set_.V()))*T;
+    }
+}
+
+
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 Foam::fv::actuationDiskSource::actuationDiskSource
@@ -130,19 +161,15 @@ Foam::wordList Foam::fv::actuationDiskSource::addSupFields() const
 
 void Foam::fv::actuationDiskSource::addSup
 (
-    fvMatrix<vector>& eqn,
+    const volVectorField& U,
-    const word& fieldName
+    fvMatrix<vector>& eqn
 ) const
 {
-    const scalarField& cellsV = mesh().V();
-    vectorField& Usource = eqn.source();
-    const vectorField& U = eqn.psi();
-
     addActuationDiskAxialInertialResistance
     (
-        Usource,
+        eqn.source(),
         set_.cells(),
-        cellsV,
+        mesh().V(),
         geometricOneField(),
         geometricOneField(),
         U
@@ -153,19 +180,15 @@ void Foam::fv::actuationDiskSource::addSup
 void Foam::fv::actuationDiskSource::addSup
 (
     const volScalarField& rho,
-    fvMatrix<vector>& eqn,
+    const volVectorField& U,
-    const word& fieldName
+    fvMatrix<vector>& eqn
 ) const
 {
-    const scalarField& cellsV = mesh().V();
-    vectorField& Usource = eqn.source();
-    const vectorField& U = eqn.psi();
-
     addActuationDiskAxialInertialResistance
     (
-        Usource,
+        eqn.source(),
         set_.cells(),
-        cellsV,
+        mesh().V(),
         geometricOneField(),
         rho,
         U
@@ -177,19 +200,15 @@ void Foam::fv::actuationDiskSource::addSup
 (
     const volScalarField& alpha,
     const volScalarField& rho,
-    fvMatrix<vector>& eqn,
+    const volVectorField& U,
-    const word& fieldName
+    fvMatrix<vector>& eqn
 ) const
 {
-    const scalarField& cellsV = mesh().V();
-    vectorField& Usource = eqn.source();
-    const vectorField& U = eqn.psi();
-
     addActuationDiskAxialInertialResistance
     (
-        Usource,
+        eqn.source(),
         set_.cells(),
-        cellsV,
+        mesh().V(),
         alpha,
         rho,
         U

src/fvModels/derived/actuationDiskSource/actuationDiskSource.H

@@ -129,7 +129,7 @@ private:
 
         //- Add resistance to the UEqn
         template<class AlphaFieldType, class RhoFieldType>
-        void addActuationDiskAxialInertialResistance
+        inline void addActuationDiskAxialInertialResistance
         (
             vectorField& Usource,
             const labelList& cells,
@@ -180,16 +180,16 @@ public:
             //- Source term to momentum equation
             virtual void addSup
             (
-                fvMatrix<vector>& eqn,
+                const volVectorField& U,
-                const word& fieldName
+                fvMatrix<vector>& eqn
             ) const;
 
             //- Source term to compressible momentum equation
             virtual void addSup
             (
                 const volScalarField& rho,
-                fvMatrix<vector>& eqn,
+                const volVectorField& U,
-                const word& fieldName
+                fvMatrix<vector>& eqn
             ) const;
 
             //- Explicit and implicit sources for phase equations
@@ -197,8 +197,8 @@ public:
             (
                 const volScalarField& alpha,
                 const volScalarField& rho,
-                fvMatrix<vector>& eqn,
+                const volVectorField& U,
-                const word& fieldName
+                fvMatrix<vector>& eqn
             ) const;
 
 
@@ -237,12 +237,6 @@ public:
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-#ifdef NoRepository
-    #include "actuationDiskSourceTemplates.C"
-#endif
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
 #endif
 
 // ************************************************************************* //

src/fvModels/derived/actuationDiskSource/actuationDiskSourceTemplates.C

@@ -1,62 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2011-2022 OpenFOAM Foundation
-     \\/     M anipulation  |
--------------------------------------------------------------------------------
-License
-    This file is part of OpenFOAM.
-
-    OpenFOAM is free software: you can redistribute it and/or modify it
-    under the terms of the GNU General Public License as published by
-    the Free Software Foundation, either version 3 of the License, or
-    (at your option) any later version.
-
-    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
-    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
-    for more details.
-
-    You should have received a copy of the GNU General Public License
-    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
-
-\*---------------------------------------------------------------------------*/
-
-#include "actuationDiskSource.H"
-#include "volFields.H"
-
-// * * * * * * * * * * * * * * *  Member Functions * * * * * * * * * * * * * //
-
-template<class AlphaFieldType, class RhoFieldType>
-void Foam::fv::actuationDiskSource::addActuationDiskAxialInertialResistance
-(
-    vectorField& Usource,
-    const labelList& cells,
-    const scalarField& Vcells,
-    const AlphaFieldType& alpha,
-    const RhoFieldType& rho,
-    const vectorField& U
-) const
-{
-    const scalar a = 1 - Cp_/Ct_;
-    const vector dHat(diskDir_/mag(diskDir_));
-
-    scalar dHatUo(vGreat);
-    if (upstreamCellId_ != -1)
-    {
-        dHatUo = dHat & U[upstreamCellId_];
-    }
-    reduce(dHatUo, minOp<scalar>());
-
-    const vector T = 2*diskArea_*sqr(dHatUo)*a*(1 - a)*dHat;
-
-    forAll(cells, i)
-    {
-        Usource[cells[i]] +=
-            (alpha[cells[i]]*rho[cells[i]]*(Vcells[cells[i]]/set_.V()))*T;
-    }
-}
-
-
-// ************************************************************************* //

src/fvModels/derived/buoyancyEnergy/buoyancyEnergy.C

@@ -96,8 +96,8 @@ Foam::wordList Foam::fv::buoyancyEnergy::addSupFields() const
 void Foam::fv::buoyancyEnergy::addSup
 (
     const volScalarField& rho,
-    fvMatrix<scalar>& eqn,
+    const volScalarField& he,
-    const word& fieldName
+    fvMatrix<scalar>& eqn
 ) const
 {
     const uniformDimensionedVectorField& g =
@@ -113,8 +113,8 @@ void Foam::fv::buoyancyEnergy::addSup
 (
     const volScalarField& alpha,
     const volScalarField& rho,
-    fvMatrix<scalar>& eqn,
+    const volScalarField& he,
-    const word& fieldName
+    fvMatrix<scalar>& eqn
 ) const
 {
     const uniformDimensionedVectorField& g =

src/fvModels/derived/buoyancyEnergy/buoyancyEnergy.H

@@ -116,8 +116,8 @@ public:
             virtual void addSup
             (
                 const volScalarField& rho,
-                fvMatrix<scalar>& eqn,
+                const volScalarField& he,
-                const word& fieldName
+                fvMatrix<scalar>& eqn
             ) const;
 
             //- Add explicit contribution to phase energy equation
@@ -125,8 +125,8 @@ public:
             (
                 const volScalarField& alpha,
                 const volScalarField& rho,
-                fvMatrix<scalar>& eqn,
+                const volScalarField& he,
-                const word& fieldName
+                fvMatrix<scalar>& eqn
             ) const;
 
 

src/fvModels/derived/buoyancyForce/buoyancyForce.C

@@ -99,8 +99,8 @@ Foam::wordList Foam::fv::buoyancyForce::addSupFields() const
 
 void Foam::fv::buoyancyForce::addSup
 (
-    fvMatrix<vector>& eqn,
+    const volVectorField& U,
-    const word& fieldName
+    fvMatrix<vector>& eqn
 ) const
 {
     eqn += g_;
@@ -110,8 +110,8 @@ void Foam::fv::buoyancyForce::addSup
 void Foam::fv::buoyancyForce::addSup
 (
     const volScalarField& rho,
-    fvMatrix<vector>& eqn,
+    const volVectorField& U,
-    const word& fieldName
+    fvMatrix<vector>& eqn
 ) const
 {
     eqn += rho*g_;
@@ -122,8 +122,8 @@ void Foam::fv::buoyancyForce::addSup
 (
     const volScalarField& alpha,
     const volScalarField& rho,
-    fvMatrix<vector>& eqn,
+    const volVectorField& U,
-    const word& fieldName
+    fvMatrix<vector>& eqn
 ) const
 {
     eqn += alpha*rho*g_;

src/fvModels/derived/buoyancyForce/buoyancyForce.H

@@ -118,16 +118,16 @@ public:
             //- Add explicit contribution to incompressible momentum equation
             virtual void addSup
             (
-                fvMatrix<vector>& eqn,
+                const volVectorField& U,
-                const word& fieldName
+                fvMatrix<vector>& eqn
             ) const;
 
             //- Add explicit contribution to compressible momentum equation
             virtual void addSup
             (
                 const volScalarField& rho,
-                fvMatrix<vector>& eqn,
+                const volVectorField& U,
-                const word& fieldName
+                fvMatrix<vector>& eqn
             ) const;
 
             //- Add explicit contribution to phase momentum equation
@@ -135,8 +135,8 @@ public:
             (
                 const volScalarField& alpha,
                 const volScalarField& rho,
-                fvMatrix<vector>& eqn,
+                const volVectorField& U,
-                const word& fieldName
+                fvMatrix<vector>& eqn
             ) const;
 
 

src/fvModels/derived/effectivenessHeatExchangerSource/effectivenessHeatExchangerSource.C

@@ -210,8 +210,8 @@ Foam::wordList Foam::fv::effectivenessHeatExchangerSource::addSupFields() const
 void Foam::fv::effectivenessHeatExchangerSource::addSup
 (
     const volScalarField& rho,
-    fvMatrix<scalar>& eqn,
+    const volScalarField& he,
-    const word& fieldName
+    fvMatrix<scalar>& eqn
 ) const
 {
     const basicThermo& thermo =

src/fvModels/derived/effectivenessHeatExchangerSource/effectivenessHeatExchangerSource.H

@@ -213,8 +213,8 @@ public:
         virtual void addSup
         (
             const volScalarField& rho,
-            fvMatrix<scalar>& eqn,
+            const volScalarField& he,
-            const word& fieldName
+            fvMatrix<scalar>& eqn
         ) const;
 
 

src/fvModels/derived/explicitPorositySource/explicitPorositySource.C

@@ -102,8 +102,8 @@ Foam::wordList Foam::fv::explicitPorositySource::addSupFields() const
 
 void Foam::fv::explicitPorositySource::addSup
 (
-    fvMatrix<vector>& eqn,
+    const volVectorField& U,
-    const word& fieldName
+    fvMatrix<vector>& eqn
 ) const
 {
     fvMatrix<vector> porosityEqn(eqn.psi(), eqn.dimensions());
@@ -115,8 +115,8 @@ void Foam::fv::explicitPorositySource::addSup
 void Foam::fv::explicitPorositySource::addSup
 (
     const volScalarField& rho,
-    fvMatrix<vector>& eqn,
+    const volVectorField& U,
-    const word& fieldName
+    fvMatrix<vector>& eqn
 ) const
 {
     fvMatrix<vector> porosityEqn(eqn.psi(), eqn.dimensions());
@@ -129,8 +129,8 @@ void Foam::fv::explicitPorositySource::addSup
 (
     const volScalarField& alpha,
     const volScalarField& rho,
-    fvMatrix<vector>& eqn,
+    const volVectorField& U,
-    const word& fieldName
+    fvMatrix<vector>& eqn
 ) const
 {
     fvMatrix<vector> porosityEqn(eqn.psi(), eqn.dimensions());

src/fvModels/derived/explicitPorositySource/explicitPorositySource.H

@@ -154,16 +154,16 @@ public:
             //- Add implicit contribution to momentum equation
             virtual void addSup
             (
-                fvMatrix<vector>& eqn,
+                const volVectorField& U,
-                const word& fieldName
+                fvMatrix<vector>& eqn
             ) const;
 
             //- Add implicit contribution to compressible momentum equation
             virtual void addSup
             (
                 const volScalarField& rho,
-                fvMatrix<vector>& eqn,
+                const volVectorField& U,
-                const word& fieldName
+                fvMatrix<vector>& eqn
             ) const;
 
             //- Add implicit contribution to phase momentum equation
@@ -171,8 +171,8 @@ public:
             (
                 const volScalarField& alpha,
                 const volScalarField& rho,
-                fvMatrix<vector>& eqn,
+                const volVectorField& U,
-                const word& fieldName
+                fvMatrix<vector>& eqn
             ) const;
 
 

src/fvModels/derived/heatSource/heatSource.C

@@ -123,8 +123,8 @@ Foam::wordList Foam::fv::heatSource::addSupFields() const
 
 void Foam::fv::heatSource::addSup
 (
-    fvMatrix<scalar>& eqn,
+    const volScalarField& he,
-    const word& fieldName
+    fvMatrix<scalar>& eqn
 ) const
 {
     const labelUList cells = set_.cells();
@@ -142,11 +142,11 @@ void Foam::fv::heatSource::addSup
 void Foam::fv::heatSource::addSup
 (
     const volScalarField& rho,
-    fvMatrix<scalar>& eqn,
+    const volScalarField& he,
-    const word& fieldName
+    fvMatrix<scalar>& eqn
 ) const
 {
-    addSup(eqn, fieldName);
+    addSup(he, eqn);
 }
 
 

src/fvModels/derived/heatSource/heatSource.H

@@ -117,16 +117,16 @@ public:
             //- Source term to energy equation
             virtual void addSup
             (
-                fvMatrix<scalar>& eqn,
+                const volScalarField& he,
-                const word& fieldName
+                fvMatrix<scalar>& eqn
             ) const;
 
             //- Source term to compressible energy equation
             virtual void addSup
             (
                 const volScalarField& rho,
-                fvMatrix<scalar>& eqn,
+                const volScalarField& he,
-                const word& fieldName
+                fvMatrix<scalar>& eqn
             ) const;
 
 

src/fvModels/derived/heatTransfer/heatTransfer.C

@@ -71,8 +71,7 @@ template<class AlphaFieldType>
 void Foam::fv::heatTransfer::add
 (
     const AlphaFieldType& alpha,
-    fvMatrix<scalar>& eqn,
+    fvMatrix<scalar>& eqn
-    const word& fieldName
 ) const
 {
     const volScalarField& he = eqn.psi();
@@ -162,22 +161,22 @@ Foam::wordList Foam::fv::heatTransfer::addSupFields() const
 
 void Foam::fv::heatTransfer::addSup
 (
-    fvMatrix<scalar>& eqn,
+    const volScalarField& he,
-    const word& fieldName
+    fvMatrix<scalar>& eqn
 ) const
 {
-    add(geometricOneField(), eqn, fieldName);
+    add(geometricOneField(), eqn);
 }
 
 
 void Foam::fv::heatTransfer::addSup
 (
     const volScalarField& rho,
-    fvMatrix<scalar>& eqn,
+    const volScalarField& he,
-    const word& fieldName
+    fvMatrix<scalar>& eqn
 ) const
 {
-    add(geometricOneField(), eqn, fieldName);
+    add(geometricOneField(), eqn);
 }
 
 
@@ -185,11 +184,11 @@ void Foam::fv::heatTransfer::addSup
 (
     const volScalarField& alpha,
     const volScalarField& rho,
-    fvMatrix<scalar>& eqn,
+    const volScalarField& he,
-    const word& fieldName
+    fvMatrix<scalar>& eqn
 ) const
 {
-    add(alpha, eqn, fieldName);
+    add(alpha, eqn);
 }
 
 

src/fvModels/derived/heatTransfer/heatTransfer.H

@@ -124,11 +124,10 @@ class heatTransfer
 
         //- Source term to energy equation
         template<class AlphaFieldType>
-        void add
+        inline void add
         (
             const AlphaFieldType& alpha,
-            fvMatrix<scalar>& eqn,
+            fvMatrix<scalar>& eqn
-            const word& fieldName
         ) const;
 
 
@@ -168,16 +167,16 @@ public:
             //- Source term to energy equation
             virtual void addSup
             (
-                fvMatrix<scalar>& eqn,
+                const volScalarField& he,
-                const word& fieldName
+                fvMatrix<scalar>& eqn
             ) const;
 
             //- Source term to compressible energy equation
             virtual void addSup
             (
                 const volScalarField& rho,
-                fvMatrix<scalar>& eqn,
+                const volScalarField& he,
-                const word& fieldName
+                fvMatrix<scalar>& eqn
             ) const;
 
             //- Source term to phase energy equation
@@ -185,8 +184,8 @@ public:
             (
                 const volScalarField& alpha,
                 const volScalarField& rho,
-                fvMatrix<scalar>& eqn,
+                const volScalarField& he,
-                const word& fieldName
+                fvMatrix<scalar>& eqn
             ) const;
 
 

src/fvModels/derived/massSource/massSource.C

@@ -74,20 +74,68 @@ void Foam::fv::massSourceBase::readCoeffs()
 
 
 template<class Type>
-void Foam::fv::massSourceBase::addGeneralSupType
+void Foam::fv::massSourceBase::addSupType
 (
-    fvMatrix<Type>& eqn,
+    const VolField<Type>& field,
-    const word& fieldName
+    fvMatrix<Type>& eqn
+) const
+{
+    FatalErrorInFunction
+        << "Cannot add a mass source for field " << field.name()
+        << " because this field's equation is not in mass-conservative form"
+        << exit(FatalError);
+}
+
+
+void Foam::fv::massSourceBase::addSupType
+(
+    const volScalarField& field,
+    fvMatrix<scalar>& eqn
+) const
+{
+    // Continuity equation. Add the mass flow rate.
+    if (field.name() == rhoName_)
+    {
+        const labelUList cells = set_.cells();
+
+        const scalar massFlowRate = this->massFlowRate();
+
+        forAll(cells, i)
+        {
+            eqn.source()[cells[i]] -=
+                mesh().V()[cells[i]]/set_.V()*massFlowRate;
+        }
+
+        return;
+    }
+
+    // Non-mass conservative property equation. Fail.
+    addSupType<scalar>(field, eqn);
+}
+
+
+template<class Type>
+void Foam::fv::massSourceBase::addSupType
+(
+    const volScalarField& rho,
+    const VolField<Type>& field,
+    fvMatrix<Type>& eqn
 ) const
 {
     const labelUList cells = set_.cells();
 
     const scalar massFlowRate = this->massFlowRate();
 
+    // Property equation. If the source is positive, introduce the value
+    // specified by the user. If negative, then sink the current internal value
+    // using an implicit term.
     if (massFlowRate > 0)
     {
         const Type value =
-            fieldValues_[fieldName]->value<Type>(mesh().time().userTimeValue());
+            fieldValues_[field.name()]->template value<Type>
+            (
+                mesh().time().userTimeValue()
+            );
 
         forAll(cells, i)
         {
@@ -106,37 +154,25 @@ void Foam::fv::massSourceBase::addGeneralSupType
 }
 
 
-template<class Type>
 void Foam::fv::massSourceBase::addSupType
 (
-    fvMatrix<Type>& eqn,
+    const volScalarField& rho,
-    const word& fieldName
+    const volScalarField& field,
+    fvMatrix<scalar>& eqn
 ) const
 {
-    addGeneralSupType(eqn, fieldName);
+    // Multiphase continuity equation. Same source as single-phase case.
-}
+    if (field.name() == rhoName_)
-
-
-void Foam::fv::massSourceBase::addSupType
-(
-    fvMatrix<scalar>& eqn,
-    const word& fieldName
-) const
-{
-    const labelUList cells = set_.cells();
-
-    if (fieldName == rhoName_)
     {
-        const scalar massFlowRate = this->massFlowRate();
+        addSupType(field, eqn);
-
+        return;
-        forAll(cells, i)
-        {
-            eqn.source()[cells[i]] -=
-                mesh().V()[cells[i]]/set_.V()*massFlowRate;
-        }
     }
-    else if (fieldName == heName_ && fieldValues_.found(TName_))
+
+    // Energy equation. Special handling for if temperature is specified.
+    if (field.name() == heName_ && fieldValues_.found(TName_))
     {
+        const labelUList cells = set_.cells();
+
         const scalar massFlowRate = this->massFlowRate();
 
         if (massFlowRate > 0)
@@ -184,23 +220,12 @@ void Foam::fv::massSourceBase::addSupType
                     mesh().V()[cells[i]]/set_.V()*massFlowRate;
             }
         }
-    }
-    else
-    {
-        addGeneralSupType(eqn, fieldName);
-    }
-}
 
+        return;
+    }
 
-template<class Type>
+    // Property equation
-void Foam::fv::massSourceBase::addSupType
+    addSupType<scalar>(rho, field, eqn);
-(
-    const volScalarField& rho,
-    fvMatrix<Type>& eqn,
-    const word& fieldName
-) const
-{
-    addSupType(eqn, fieldName);
 }
 
 
@@ -209,11 +234,12 @@ void Foam::fv::massSourceBase::addSupType
 (
     const volScalarField& alpha,
     const volScalarField& rho,
-    fvMatrix<Type>& eqn,
+    const VolField<Type>& field,
-    const word& fieldName
+    fvMatrix<Type>& eqn
 ) const
 {
-    addSupType(eqn, fieldName);
+    // Multiphase property equation. Same source as the single phase case.
+    addSupType(rho, field, eqn);
 }
 
 
@@ -300,11 +326,12 @@ Foam::fv::massSource::massSource
 
 bool Foam::fv::massSourceBase::addsSupToField(const word& fieldName) const
 {
+    const bool isMixture = IOobject::group(fieldName) == word::null;
     const bool isThisPhase = IOobject::group(fieldName) == phaseName_;
 
     if
     (
-        isThisPhase
+        (isMixture || isThisPhase)
      && massFlowRate() > 0
      && !(fieldName == rhoName_)
      && !(fieldName == heName_ && fieldValues_.found(TName_))
@@ -318,7 +345,7 @@ bool Foam::fv::massSourceBase::addsSupToField(const word& fieldName) const
         return false;
     }
 
-    return isThisPhase;
+    return isMixture || isThisPhase;
 }
 
 
@@ -335,13 +362,25 @@ Foam::wordList Foam::fv::massSourceBase::addSupFields() const
 }
 
 
-FOR_ALL_FIELD_TYPES(IMPLEMENT_FV_MODEL_ADD_SUP, fv::massSourceBase);
+FOR_ALL_FIELD_TYPES
+(
+    IMPLEMENT_FV_MODEL_ADD_FIELD_SUP,
+    fv::massSourceBase
+)
 
 
-FOR_ALL_FIELD_TYPES(IMPLEMENT_FV_MODEL_ADD_RHO_SUP, fv::massSourceBase);
+FOR_ALL_FIELD_TYPES
+(
+    IMPLEMENT_FV_MODEL_ADD_RHO_FIELD_SUP,
+    fv::massSourceBase
+)
 
 
-FOR_ALL_FIELD_TYPES(IMPLEMENT_FV_MODEL_ADD_ALPHA_RHO_SUP, fv::massSourceBase);
+FOR_ALL_FIELD_TYPES
+(
+    IMPLEMENT_FV_MODEL_ADD_ALPHA_RHO_FIELD_SUP,
+    fv::massSourceBase
+)
 
 
 bool Foam::fv::massSourceBase::movePoints()

src/fvModels/derived/massSource/massSource.H

@@ -123,26 +123,34 @@ private:
 
             //- Add a source term to an equation
             template<class Type>
-            void addGeneralSupType
+            void addSupType
             (
-                fvMatrix<Type>& eqn,
+                const VolField<Type>& field,
-                const word& fieldName
+                fvMatrix<Type>& eqn
             ) const;
 
-            //- Add a source term to an equation
-            template<class Type>
-            void addSupType(fvMatrix<Type>& eqn, const word& fieldName) const;
-
             //- Add a source term to a scalar equation
-            void addSupType(fvMatrix<scalar>& eqn, const word& fieldName) const;
+            void addSupType
+            (
+                const volScalarField& field,
+                fvMatrix<scalar>& eqn
+            ) const;
 
             //- Add a source term to a compressible equation
             template<class Type>
             void addSupType
             (
                 const volScalarField& rho,
-                fvMatrix<Type>& eqn,
+                const VolField<Type>& field,
-                const word& fieldName
+                fvMatrix<Type>& eqn
+            ) const;
+
+            //- Add a source term to a compressible scalar equation
+            void addSupType
+            (
+                const volScalarField& rho,
+                const volScalarField& field,
+                fvMatrix<scalar>& eqn
             ) const;
 
             //- Add a source term to a phase equation
@@ -151,8 +159,8 @@ private:
             (
                 const volScalarField& alpha,
                 const volScalarField& rho,
-                fvMatrix<Type>& eqn,
+                const VolField<Type>& field,
-                const word& fieldName
+                fvMatrix<Type>& eqn
             ) const;
 
 
@@ -204,13 +212,13 @@ public:
         // Sources
 
             //- Add a source term to an equation
-            FOR_ALL_FIELD_TYPES(DEFINE_FV_MODEL_ADD_SUP);
+            FOR_ALL_FIELD_TYPES(DEFINE_FV_MODEL_ADD_FIELD_SUP)
 
             //- Add a source term to a compressible equation
-            FOR_ALL_FIELD_TYPES(DEFINE_FV_MODEL_ADD_RHO_SUP);
+            FOR_ALL_FIELD_TYPES(DEFINE_FV_MODEL_ADD_RHO_FIELD_SUP)
 
             //- Add a source term to a phase equation
-            FOR_ALL_FIELD_TYPES(DEFINE_FV_MODEL_ADD_ALPHA_RHO_SUP);
+            FOR_ALL_FIELD_TYPES(DEFINE_FV_MODEL_ADD_ALPHA_RHO_FIELD_SUP)
 
 
         // Mesh changes

src/fvModels/derived/phaseLimitStabilisation/phaseLimitStabilisation.C

@@ -62,16 +62,14 @@ void Foam::fv::phaseLimitStabilisation::addSupType
 (
     const volScalarField& alpha,
     const volScalarField& rho,
-    fvMatrix<Type>& eqn,
+    const VolField<Type>& field,
-    const word& fieldName
+    fvMatrix<Type>& eqn
 ) const
 {
-    const VolField<Type>& psi = eqn.psi();
+    const uniformDimensionedScalarField& rate =
-
-    uniformDimensionedScalarField& rate =
         mesh().lookupObjectRef<uniformDimensionedScalarField>(rateName_);
 
-    eqn -= fvm::Sp(max(residualAlpha_ - alpha, scalar(0))*rho*rate, psi);
+    eqn -= fvm::Sp(max(residualAlpha_ - alpha, scalar(0))*rho*rate, eqn.psi());
 }
 
 
@@ -104,9 +102,9 @@ Foam::wordList Foam::fv::phaseLimitStabilisation::addSupFields() const
 
 FOR_ALL_FIELD_TYPES
 (
-    IMPLEMENT_FV_MODEL_ADD_ALPHA_RHO_SUP,
+    IMPLEMENT_FV_MODEL_ADD_ALPHA_RHO_FIELD_SUP,
     fv::phaseLimitStabilisation
-);
+)
 
 
 bool Foam::fv::phaseLimitStabilisation::movePoints()

src/fvModels/derived/phaseLimitStabilisation/phaseLimitStabilisation.H

@@ -98,8 +98,8 @@ class phaseLimitStabilisation
         (
             const volScalarField& alpha,
             const volScalarField& rho,
-            fvMatrix<Type>& eqn,
+            const VolField<Type>& field,
-            const word& fieldName
+            fvMatrix<Type>& eqn
         ) const;
 
 
@@ -141,7 +141,7 @@ public:
         // Sources
 
             //- Add a source term to a phase equation
-            FOR_ALL_FIELD_TYPES(DEFINE_FV_MODEL_ADD_ALPHA_RHO_SUP);
+            FOR_ALL_FIELD_TYPES(DEFINE_FV_MODEL_ADD_ALPHA_RHO_FIELD_SUP)
 
 
         // Mesh changes

src/fvModels/derived/radialActuationDiskSource/radialActuationDiskSource.C

@@ -24,6 +24,8 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "radialActuationDiskSource.H"
+#include "volFields.H"
+#include "fvMatrix.H"
 #include "geometricOneField.H"
 #include "addToRunTimeSelectionTable.H"
 
@@ -51,6 +53,70 @@ void Foam::fv::radialActuationDiskSource::readCoeffs()
 }
 
 
+
+template<class RhoFieldType>
+void Foam::fv::radialActuationDiskSource::
+addRadialActuationDiskAxialInertialResistance
+(
+    vectorField& Usource,
+    const labelList& cells,
+    const scalarField& Vcells,
+    const RhoFieldType& rho,
+    const vectorField& U
+) const
+{
+    scalar a = 1.0 - Cp_/Ct_;
+    scalarField Tr(cells.size());
+    const vector uniDiskDir = diskDir_/mag(diskDir_);
+
+    tensor E(Zero);
+    E.xx() = uniDiskDir.x();
+    E.yy() = uniDiskDir.y();
+    E.zz() = uniDiskDir.z();
+
+    const Field<vector> zoneCellCentres(mesh().cellCentres(), cells);
+    const Field<scalar> zoneCellVolumes(mesh().cellVolumes(), cells);
+
+    const vector avgCentre = gSum(zoneCellVolumes*zoneCellCentres)/set_.V();
+    const scalar maxR = gMax(mag(zoneCellCentres - avgCentre));
+
+    scalar intCoeffs =
+        radialCoeffs_[0]
+      + radialCoeffs_[1]*sqr(maxR)/2.0
+      + radialCoeffs_[2]*pow4(maxR)/3.0;
+
+    vector upU = vector(vGreat, vGreat, vGreat);
+    scalar upRho = vGreat;
+    if (upstreamCellId_ != -1)
+    {
+        upU =  U[upstreamCellId_];
+        upRho = rho[upstreamCellId_];
+    }
+    reduce(upU, minOp<vector>());
+    reduce(upRho, minOp<scalar>());
+
+    scalar T = 2.0*upRho*diskArea_*mag(upU)*a*(1.0 - a);
+    forAll(cells, i)
+    {
+        scalar r2 = magSqr(mesh().cellCentres()[cells[i]] - avgCentre);
+
+        Tr[i] =
+            T
+           *(radialCoeffs_[0] + radialCoeffs_[1]*r2 + radialCoeffs_[2]*sqr(r2))
+           /intCoeffs;
+
+        Usource[cells[i]] += ((Vcells[cells[i]]/set_.V())*Tr[i]*E) & upU;
+    }
+
+    if (debug)
+    {
+        Info<< "Source name: " << name() << nl
+            << "Average centre: " << avgCentre << nl
+            << "Maximum radius: " << maxR << endl;
+    }
+}
+
+
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 Foam::fv::radialActuationDiskSource::radialActuationDiskSource
@@ -72,19 +138,15 @@ Foam::fv::radialActuationDiskSource::radialActuationDiskSource
 
 void Foam::fv::radialActuationDiskSource::addSup
 (
-    fvMatrix<vector>& eqn,
+    const volVectorField& U,
-    const word& fieldName
+    fvMatrix<vector>& eqn
 ) const
 {
-    const scalarField& cellsV = mesh().V();
-    vectorField& Usource = eqn.source();
-    const vectorField& U = eqn.psi();
-
     addRadialActuationDiskAxialInertialResistance
     (
-        Usource,
+        eqn.source(),
         set_.cells(),
-        cellsV,
+        mesh().V(),
         geometricOneField(),
         U
     );
@@ -94,19 +156,15 @@ void Foam::fv::radialActuationDiskSource::addSup
 void Foam::fv::radialActuationDiskSource::addSup
 (
     const volScalarField& rho,
-    fvMatrix<vector>& eqn,
+    const volVectorField& U,
-    const word& fieldName
+    fvMatrix<vector>& eqn
 ) const
 {
-    const scalarField& cellsV = mesh().V();
-    vectorField& Usource = eqn.source();
-    const vectorField& U = eqn.psi();
-
     addRadialActuationDiskAxialInertialResistance
     (
-        Usource,
+        eqn.source(),
         set_.cells(),
-        cellsV,
+        mesh().V(),
         rho,
         U
     );

src/fvModels/derived/radialActuationDiskSource/radialActuationDiskSource.H

@@ -115,7 +115,7 @@ class radialActuationDiskSource
 
         //- Add resistance to the UEqn
         template<class RhoFieldType>
-        void addRadialActuationDiskAxialInertialResistance
+        inline void addRadialActuationDiskAxialInertialResistance
         (
             vectorField& Usource,
             const labelList& cells,
@@ -156,16 +156,16 @@ public:
         //- Source term to momentum equation
         virtual void addSup
         (
-            fvMatrix<vector>& eqn,
+            const volVectorField& U,
-            const word& fieldName
+            fvMatrix<vector>& eqn
         ) const;
 
         //- Source term to compressible momentum equation
         virtual void addSup
         (
             const volScalarField& rho,
-            fvMatrix<vector>& eqn,
+            const volVectorField& U,
-            const word& fieldName
+            fvMatrix<vector>& eqn
         ) const;
 
 
@@ -189,12 +189,6 @@ public:
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-#ifdef NoRepository
-    #include "radialActuationDiskSourceTemplates.C"
-#endif
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
 #endif
 
 // ************************************************************************* //

src/fvModels/derived/rotorDiskSource/rotorDiskSource.C

@@ -554,8 +554,8 @@ Foam::wordList Foam::fv::rotorDiskSource::addSupFields() const
 
 void Foam::fv::rotorDiskSource::addSup
 (
-    fvMatrix<vector>& eqn,
+    const volVectorField& U,
-    const word& fieldName
+    fvMatrix<vector>& eqn
 ) const
 {
     volVectorField::Internal force
@@ -578,7 +578,7 @@ void Foam::fv::rotorDiskSource::addSup
     // Read the reference density for incompressible flow
     coeffs().lookup("rhoRef") >> rhoRef_;
 
-    const vectorField Uin(inflowVelocity(eqn.psi()));
+    const vectorField Uin(inflowVelocity(U));
     trim_->correct(Uin, force);
     calculate(geometricOneField(), Uin, trim_->thetag(), force);
 
@@ -595,8 +595,8 @@ void Foam::fv::rotorDiskSource::addSup
 void Foam::fv::rotorDiskSource::addSup
 (
     const volScalarField& rho,
-    fvMatrix<vector>& eqn,
+    const volVectorField& U,
-    const word& fieldName
+    fvMatrix<vector>& eqn
 ) const
 {
     volVectorField::Internal force
@@ -616,7 +616,7 @@ void Foam::fv::rotorDiskSource::addSup
         )
     );
 
-    const vectorField Uin(inflowVelocity(eqn.psi()));
+    const vectorField Uin(inflowVelocity(U));
     trim_->correct(rho, Uin, force);
     calculate(rho, Uin, trim_->thetag(), force);
 

src/fvModels/derived/rotorDiskSource/rotorDiskSource.H

@@ -318,16 +318,16 @@ public:
             //- Source term to momentum equation
             virtual void addSup
             (
-                fvMatrix<vector>& eqn,
+                const volVectorField& U,
-                const word& fieldName
+                fvMatrix<vector>& eqn
             ) const;
 
             //- Source term to compressible momentum equation
             virtual void addSup
             (
                 const volScalarField& rho,
-                fvMatrix<vector>& eqn,
+                const volVectorField& U,
-                const word& fieldName
+                fvMatrix<vector>& eqn
             ) const;
 
 

src/fvModels/derived/sixDoFAccelerationSource/sixDoFAccelerationSource.C

@@ -26,6 +26,7 @@ License
 #include "sixDoFAccelerationSource.H"
 #include "fvMatrices.H"
 #include "geometricOneField.H"
+#include "uniformDimensionedFields.H"
 #include "makeFunction1s.H"
 #include "makeTableReaders.H"
 #include "addToRunTimeSelectionTable.H"
@@ -105,6 +106,68 @@ void Foam::fv::sixDoFAccelerationSource::readCoeffs()
 }
 
 
+template<class AlphaFieldType, class RhoFieldType>
+void Foam::fv::sixDoFAccelerationSource::addForce
+(
+    const AlphaFieldType& alpha,
+    const RhoFieldType& rho,
+    const volVectorField& U,
+    fvMatrix<vector>& eqn
+) const
+{
+    const Vector<vector> accelerations
+    (
+        accelerations_->value(mesh().time().userTimeValue())
+    );
+
+    // If gravitational force is present combine with the linear acceleration
+    if (mesh().foundObject<uniformDimensionedVectorField>("g"))
+    {
+        uniformDimensionedVectorField& g =
+            mesh().lookupObjectRef<uniformDimensionedVectorField>("g");
+
+        const uniformDimensionedScalarField& hRef =
+            mesh().lookupObject<uniformDimensionedScalarField>("hRef");
+
+        g = g_ - dimensionedVector("a", dimAcceleration, accelerations.x());
+
+        dimensionedScalar ghRef(- mag(g)*hRef);
+
+        mesh().lookupObjectRef<volScalarField>("gh") = (g & mesh().C()) - ghRef;
+
+        mesh().lookupObjectRef<surfaceScalarField>("ghf") =
+            (g & mesh().Cf()) - ghRef;
+    }
+    // ... otherwise include explicitly in the momentum equation
+    else
+    {
+        const dimensionedVector a("a", dimAcceleration, accelerations.x());
+        eqn -= alpha*rho*a;
+    }
+
+    dimensionedVector Omega
+    (
+        "Omega",
+        dimensionSet(0, 0, -1, 0, 0),
+        accelerations.y()
+    );
+
+    dimensionedVector dOmegaDT
+    (
+        "dOmegaDT",
+        dimensionSet(0, 0, -2, 0, 0),
+        accelerations.z()
+    );
+
+    eqn -=
+    (
+        alpha*rho*(2*Omega ^ U)                  // Coriolis force
+      + alpha*rho*(Omega ^ (Omega ^ mesh().C())) // Centrifugal force
+      + alpha*rho*(dOmegaDT ^ mesh().C())        // Angular acceleration force
+    );
+}
+
+
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 Foam::fv::sixDoFAccelerationSource::sixDoFAccelerationSource
@@ -139,22 +202,22 @@ Foam::wordList Foam::fv::sixDoFAccelerationSource::addSupFields() const
 
 void Foam::fv::sixDoFAccelerationSource::addSup
 (
-    fvMatrix<vector>& eqn,
+    const volVectorField& U,
-    const word& fieldName
+    fvMatrix<vector>& eqn
 ) const
 {
-    addForce(geometricOneField(), geometricOneField(), eqn, fieldName);
+    addForce(geometricOneField(), geometricOneField(), U, eqn);
 }
 
 
 void Foam::fv::sixDoFAccelerationSource::addSup
 (
     const volScalarField& rho,
-    fvMatrix<vector>& eqn,
+    const volVectorField& U,
-    const word& fieldName
+    fvMatrix<vector>& eqn
 ) const
 {
-    addForce(geometricOneField(), rho, eqn, fieldName);
+    addForce(geometricOneField(), rho, U, eqn);
 }
 
 
@@ -162,11 +225,11 @@ void Foam::fv::sixDoFAccelerationSource::addSup
 (
     const volScalarField& alpha,
     const volScalarField& rho,
-    fvMatrix<vector>& eqn,
+    const volVectorField& U,
-    const word& fieldName
+    fvMatrix<vector>& eqn
 ) const
 {
-    addForce(alpha, rho, eqn, fieldName);
+    addForce(alpha, rho, U, eqn);
 }
 
 

src/fvModels/derived/sixDoFAccelerationSource/sixDoFAccelerationSource.H

@@ -98,12 +98,12 @@ private:
 
         //- Add force to a momentum equation
         template<class AlphaFieldType, class RhoFieldType>
-        void addForce
+        inline void addForce
         (
             const AlphaFieldType& alpha,
             const RhoFieldType& rho,
-            fvMatrix<vector>& eqn,
+            const volVectorField& U,
-            const word& fieldName
+            fvMatrix<vector>& eqn
         ) const;
 
 
@@ -145,16 +145,16 @@ public:
         //- Source term to momentum equation
         virtual void addSup
         (
-            fvMatrix<vector>& eqn,
+            const volVectorField& U,
-            const word& fieldName
+            fvMatrix<vector>& eqn
         ) const;
 
         //- Source term to compressible momentum equation
         virtual void addSup
         (
             const volScalarField& rho,
-            fvMatrix<vector>& eqn,
+            const volVectorField& U,
-            const word& fieldName
+            fvMatrix<vector>& eqn
         ) const;
 
         //- Source term to phase momentum equation
@@ -162,8 +162,8 @@ public:
         (
             const volScalarField& alpha,
             const volScalarField& rho,
-            fvMatrix<vector>& eqn,
+            const volVectorField& U,
-            const word& fieldName
+            fvMatrix<vector>& eqn
         ) const;
 
 
@@ -200,12 +200,6 @@ public:
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-#ifdef NoRepository
-    #include "sixDoFAccelerationSourceTemplates.C"
-#endif
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
 #endif
 
 // ************************************************************************* //

src/fvModels/derived/sixDoFAccelerationSource/sixDoFAccelerationSourceTemplates.C

@@ -1,95 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2015-2023 OpenFOAM Foundation
-     \\/     M anipulation  |
--------------------------------------------------------------------------------
-License
-    This file is part of OpenFOAM.
-
-    OpenFOAM is free software: you can redistribute it and/or modify it
-    under the terms of the GNU General Public License as published by
-    the Free Software Foundation, either version 3 of the License, or
-    (at your option) any later version.
-
-    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
-    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
-    for more details.
-
-    You should have received a copy of the GNU General Public License
-    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
-
-\*---------------------------------------------------------------------------*/
-
-#include "sixDoFAccelerationSource.H"
-#include "fvMesh.H"
-#include "fvMatrices.H"
-#include "uniformDimensionedFields.H"
-
-// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
-
-template<class AlphaFieldType, class RhoFieldType>
-void Foam::fv::sixDoFAccelerationSource::addForce
-(
-    const AlphaFieldType& alpha,
-    const RhoFieldType& rho,
-    fvMatrix<vector>& eqn,
-    const word& fieldName
-) const
-{
-    const Vector<vector> accelerations
-    (
-        accelerations_->value(mesh().time().userTimeValue())
-    );
-
-    // If gravitational force is present combine with the linear acceleration
-    if (mesh().foundObject<uniformDimensionedVectorField>("g"))
-    {
-        uniformDimensionedVectorField& g =
-            mesh().lookupObjectRef<uniformDimensionedVectorField>("g");
-
-        const uniformDimensionedScalarField& hRef =
-            mesh().lookupObject<uniformDimensionedScalarField>("hRef");
-
-        g = g_ - dimensionedVector("a", dimAcceleration, accelerations.x());
-
-        dimensionedScalar ghRef(- mag(g)*hRef);
-
-        mesh().lookupObjectRef<volScalarField>("gh") = (g & mesh().C()) - ghRef;
-
-        mesh().lookupObjectRef<surfaceScalarField>("ghf") =
-            (g & mesh().Cf()) - ghRef;
-    }
-    // ... otherwise include explicitly in the momentum equation
-    else
-    {
-        const dimensionedVector a("a", dimAcceleration, accelerations.x());
-        eqn -= alpha*rho*a;
-    }
-
-    dimensionedVector Omega
-    (
-        "Omega",
-        dimensionSet(0, 0, -1, 0, 0),
-        accelerations.y()
-    );
-
-    dimensionedVector dOmegaDT
-    (
-        "dOmegaDT",
-        dimensionSet(0, 0, -2, 0, 0),
-        accelerations.z()
-    );
-
-    eqn -=
-    (
-        alpha*rho*(2*Omega ^ eqn.psi())          // Coriolis force
-      + alpha*rho*(Omega ^ (Omega ^ mesh().C())) // Centrifugal force
-      + alpha*rho*(dOmegaDT ^ mesh().C())        // Angular acceleration force
-    );
-}
-
-
-// ************************************************************************* //

src/fvModels/derived/solidEquilibriumEnergySource/solidEquilibriumEnergySource.C

@@ -144,8 +144,8 @@ Foam::wordList Foam::fv::solidEquilibriumEnergySource::addSupFields() const
 void Foam::fv::solidEquilibriumEnergySource::addSup
 (
     const volScalarField& rho,
-    fvMatrix<scalar>& eqn,
+    const volScalarField& he,
-    const word& fieldName
+    fvMatrix<scalar>& eqn
 ) const
 {
     const volScalarField alphahe
@@ -172,8 +172,8 @@ void Foam::fv::solidEquilibriumEnergySource::addSup
 (
     const volScalarField& alpha,
     const volScalarField& rho,
-    fvMatrix<scalar>& eqn,
+    const volScalarField& he,
-    const word& fieldName
+    fvMatrix<scalar>& eqn
 ) const
 {
     const volScalarField alphahe

src/fvModels/derived/solidEquilibriumEnergySource/solidEquilibriumEnergySource.H

@@ -138,18 +138,18 @@ public:
             //- Explicit and implicit sources for compressible equations
             virtual void addSup
             (
-                const volScalarField&,
+                const volScalarField& rho,
-                fvMatrix<scalar>&,
+                const volScalarField& he,
-                const word& fieldName
+                fvMatrix<scalar>& eqn
             ) const;
 
             //- Explicit and implicit sources for phase equations
             virtual void addSup
             (
-                const volScalarField&,
+                const volScalarField& alpha,
-                const volScalarField&,
+                const volScalarField& rho,
-                fvMatrix<scalar>&,
+                const volScalarField& he,
-                const word& fieldName
+                fvMatrix<scalar>& eqn
             ) const;
 
 

src/fvModels/derived/solidificationMeltingSource/solidificationMeltingSource.C

@@ -24,6 +24,7 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "solidificationMeltingSource.H"
+#include "fvcDdt.H"
 #include "fvMatrices.H"
 #include "basicThermo.H"
 #include "uniformDimensionedFields.H"
@@ -214,6 +215,36 @@ void Foam::fv::solidificationMeltingSource::update
 }
 
 
+template<class RhoFieldType>
+void Foam::fv::solidificationMeltingSource::apply
+(
+    const RhoFieldType& rho,
+    fvMatrix<scalar>& eqn
+) const
+{
+    if (debug)
+    {
+        Info<< type() << ": applying source to " << eqn.psi().name() << endl;
+    }
+
+    const volScalarField Cp(this->Cp());
+
+    update(Cp);
+
+    dimensionedScalar L("L", dimEnergy/dimMass, L_);
+
+    // Contributions added to rhs of solver equation
+    if (eqn.psi().dimensions() == dimTemperature)
+    {
+        eqn -= L/Cp*(fvc::ddt(rho, alpha1_));
+    }
+    else
+    {
+        eqn -= L*(fvc::ddt(rho, alpha1_));
+    }
+}
+
+
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 Foam::fv::solidificationMeltingSource::solidificationMeltingSource
@@ -286,8 +317,8 @@ Foam::wordList Foam::fv::solidificationMeltingSource::addSupFields() const
 
 void Foam::fv::solidificationMeltingSource::addSup
 (
-    fvMatrix<scalar>& eqn,
+    const volScalarField& he,
-    const word& fieldName
+    fvMatrix<scalar>& eqn
 ) const
 {
     apply(geometricOneField(), eqn);
@@ -297,8 +328,8 @@ void Foam::fv::solidificationMeltingSource::addSup
 void Foam::fv::solidificationMeltingSource::addSup
 (
     const volScalarField& rho,
-    fvMatrix<scalar>& eqn,
+    const volScalarField& he,
-    const word& fieldName
+    fvMatrix<scalar>& eqn
 ) const
 {
     apply(rho, eqn);
@@ -307,8 +338,8 @@ void Foam::fv::solidificationMeltingSource::addSup
 
 void Foam::fv::solidificationMeltingSource::addSup
 (
-    fvMatrix<vector>& eqn,
+    const volVectorField& U,
-    const word& fieldName
+    fvMatrix<vector>& eqn
 ) const
 {
     if (debug)
@@ -347,12 +378,11 @@ void Foam::fv::solidificationMeltingSource::addSup
 void Foam::fv::solidificationMeltingSource::addSup
 (
     const volScalarField& rho,
-    fvMatrix<vector>& eqn,
+    const volVectorField& U,
-    const word& fieldName
+    fvMatrix<vector>& eqn
 ) const
 {
-    // Momentum source uses a Boussinesq approximation - redirect
+    addSup(U, eqn);
-    addSup(eqn, fieldName);
 }
 
 

src/fvModels/derived/solidificationMeltingSource/solidificationMeltingSource.H

@@ -230,7 +230,7 @@ private:
 
         //- Helper function to apply to the energy equation
         template<class RhoFieldType>
-        void apply(const RhoFieldType& rho, fvMatrix<scalar>& eqn) const;
+        inline void apply(const RhoFieldType& rho, fvMatrix<scalar>& eqn) const;
 
 
 public:
@@ -266,39 +266,36 @@ public:
             virtual wordList addSupFields() const;
 
 
-        // Add explicit and implicit contributions
+        // Evaluation
 
             //- Add explicit contribution to enthalpy equation
             virtual void addSup
             (
-                fvMatrix<scalar>& eqn,
+                const volScalarField& he,
-                const word& fieldName
+                fvMatrix<scalar>& eqn
             ) const;
 
             //- Add implicit contribution to momentum equation
             virtual void addSup
             (
-                fvMatrix<vector>& eqn,
+                const volVectorField& U,
-                const word& fieldName
+                fvMatrix<vector>& eqn
             ) const;
 
-
-        // Add explicit and implicit contributions to compressible equation
-
             //- Add explicit contribution to compressible enthalpy equation
             virtual void addSup
             (
                 const volScalarField& rho,
-                fvMatrix<scalar>& eqn,
+                const volScalarField& he,
-                const word& fieldName
+                fvMatrix<scalar>& eqn
             ) const;
 
             //- Add implicit contribution to compressible momentum equation
             virtual void addSup
             (
                 const volScalarField& rho,
-                fvMatrix<vector>& eqn,
+                const volVectorField& U,
-                const word& fieldName
+                fvMatrix<vector>& eqn
             ) const;
 
 
@@ -337,12 +334,6 @@ public:
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-#ifdef NoRepository
-    #include "solidificationMeltingSourceTemplates.C"
-#endif
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
 #endif
 
 // ************************************************************************* //

src/fvModels/derived/solidificationMeltingSource/solidificationMeltingSourceTemplates.C

@@ -1,61 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2014-2021 OpenFOAM Foundation
-     \\/     M anipulation  |
--------------------------------------------------------------------------------
-License
-    This file is part of OpenFOAM.
-
-    OpenFOAM is free software: you can redistribute it and/or modify it
-    under the terms of the GNU General Public License as published by
-    the Free Software Foundation, either version 3 of the License, or
-    (at your option) any later version.
-
-    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
-    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
-    for more details.
-
-    You should have received a copy of the GNU General Public License
-    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
-
-\*---------------------------------------------------------------------------*/
-
-#include "fvMatrices.H"
-#include "fvcDdt.H"
-
-// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
-
-template<class RhoFieldType>
-void Foam::fv::solidificationMeltingSource::apply
-(
-    const RhoFieldType& rho,
-    fvMatrix<scalar>& eqn
-) const
-{
-    if (debug)
-    {
-        Info<< type() << ": applying source to " << eqn.psi().name() << endl;
-    }
-
-    const volScalarField Cp(this->Cp());
-
-    update(Cp);
-
-    dimensionedScalar L("L", dimEnergy/dimMass, L_);
-
-    // Contributions added to rhs of solver equation
-    if (eqn.psi().dimensions() == dimTemperature)
-    {
-        eqn -= L/Cp*(fvc::ddt(rho, alpha1_));
-    }
-    else
-    {
-        eqn -= L*(fvc::ddt(rho, alpha1_));
-    }
-}
-
-
-// ************************************************************************* //

src/fvModels/derived/volumeFractionSource/volumeFractionSource.C

@@ -128,21 +128,20 @@ Foam::tmp<Foam::volScalarField> Foam::fv::volumeFractionSource::D
 
 
 template <class Type, class AlphaFieldType>
-void Foam::fv::volumeFractionSource::addGeneralSup
+void Foam::fv::volumeFractionSource::addGeneralSupType
 (
     const AlphaFieldType& alpha,
-    fvMatrix<Type>& eqn,
+    fvMatrix<Type>& eqn
-    const word& fieldName
 ) const
 {
     const word phiName =
-        IOobject::groupName(phiName_, IOobject::group(fieldName));
+        IOobject::groupName(phiName_, IOobject::group(eqn.psi().name()));
     const surfaceScalarField& phi =
         mesh().lookupObject<surfaceScalarField>(phiName);
 
     const volScalarField B(1 - volumeAlpha());
     const volScalarField AByB(volumeAlpha()/B);
-    const volScalarField D(this->D(fieldName));
+    const volScalarField D(this->D(eqn.psi().name()));
 
     // Divergence term
     const word divScheme = "div(" + phiName + "," + eqn.psi().name() + ")";
@@ -161,11 +160,11 @@ template<class Type, class AlphaFieldType>
 void Foam::fv::volumeFractionSource::addAlphaSupType
 (
     const AlphaFieldType& alpha,
-    fvMatrix<Type>& eqn,
+    const VolField<Type>& field,
-    const word& fieldName
+    fvMatrix<Type>& eqn
 ) const
 {
-    addGeneralSup(alpha, eqn, fieldName);
+    addGeneralSupType(alpha, eqn);
 }
 
 
@@ -173,14 +172,14 @@ template<class AlphaFieldType>
 void Foam::fv::volumeFractionSource::addAlphaSupType
 (
     const AlphaFieldType& alpha,
-    fvMatrix<scalar>& eqn,
+    const volScalarField& field,
-    const word& fieldName
+    fvMatrix<scalar>& eqn
 ) const
 {
-    if (IOobject::member(fieldName) == rhoName_)
+    if (IOobject::member(field.name()) == rhoName_)
     {
         const word phiName =
-            IOobject::groupName(phiName_, IOobject::group(fieldName));
+            IOobject::groupName(phiName_, IOobject::group(field.name()));
         const surfaceScalarField& phi =
             mesh().lookupObject<surfaceScalarField>(phiName);
 
@@ -190,7 +189,7 @@ void Foam::fv::volumeFractionSource::addAlphaSupType
     }
     else
     {
-        addGeneralSup(alpha, eqn, fieldName);
+        addGeneralSupType(alpha, eqn);
     }
 }
 
@@ -199,14 +198,14 @@ template<class AlphaFieldType>
 void Foam::fv::volumeFractionSource::addAlphaSupType
 (
     const AlphaFieldType& alpha,
-    fvMatrix<vector>& eqn,
+    const volVectorField& field,
-    const word& fieldName
+    fvMatrix<vector>& eqn
 ) const
 {
-    if (IOobject::member(fieldName) == UName_)
+    if (IOobject::member(field.name()) == UName_)
     {
         const word phiName =
-            IOobject::groupName(phiName_, IOobject::group(fieldName));
+            IOobject::groupName(phiName_, IOobject::group(field.name()));
         const surfaceScalarField& phi =
             mesh().lookupObject<surfaceScalarField>(phiName);
 
@@ -218,7 +217,7 @@ void Foam::fv::volumeFractionSource::addAlphaSupType
     }
     else
     {
-        addGeneralSup(alpha, eqn, fieldName);
+        addGeneralSupType(alpha, eqn);
     }
 }
 
@@ -226,11 +225,11 @@ void Foam::fv::volumeFractionSource::addAlphaSupType
 template<class Type>
 void Foam::fv::volumeFractionSource::addSupType
 (
-    fvMatrix<Type>& eqn,
+    const VolField<Type>& field,
-    const word& fieldName
+    fvMatrix<Type>& eqn
 ) const
 {
-    addAlphaSupType(geometricOneField(), eqn, fieldName);
+    addAlphaSupType(geometricOneField(), field, eqn);
 }
 
 
@@ -238,11 +237,11 @@ template<class Type>
 void Foam::fv::volumeFractionSource::addSupType
 (
     const volScalarField& rho,
-    fvMatrix<Type>& eqn,
+    const VolField<Type>& field,
-    const word& fieldName
+    fvMatrix<Type>& eqn
 ) const
 {
-    addAlphaSupType(geometricOneField(), eqn, fieldName);
+    addAlphaSupType(geometricOneField(), field, eqn);
 }
 
 
@@ -251,11 +250,11 @@ void Foam::fv::volumeFractionSource::addSupType
 (
     const volScalarField& alpha,
     const volScalarField& rho,
-    fvMatrix<Type>& eqn,
+    const VolField<Type>& field,
-    const word& fieldName
+    fvMatrix<Type>& eqn
 ) const
 {
-    addAlphaSupType(alpha, eqn, fieldName);
+    addAlphaSupType(alpha, field, eqn);
 }
 
 
@@ -300,17 +299,25 @@ Foam::wordList Foam::fv::volumeFractionSource::addSupFields() const
 }
 
 
-FOR_ALL_FIELD_TYPES(IMPLEMENT_FV_MODEL_ADD_SUP, fv::volumeFractionSource);
+FOR_ALL_FIELD_TYPES
-
+(
-
+    IMPLEMENT_FV_MODEL_ADD_FIELD_SUP,
-FOR_ALL_FIELD_TYPES(IMPLEMENT_FV_MODEL_ADD_RHO_SUP, fv::volumeFractionSource);
+    fv::volumeFractionSource
+)
 
 
 FOR_ALL_FIELD_TYPES
 (
-    IMPLEMENT_FV_MODEL_ADD_ALPHA_RHO_SUP,
+    IMPLEMENT_FV_MODEL_ADD_RHO_FIELD_SUP,
     fv::volumeFractionSource
-);
+)
+
+
+FOR_ALL_FIELD_TYPES
+(
+    IMPLEMENT_FV_MODEL_ADD_ALPHA_RHO_FIELD_SUP,
+    fv::volumeFractionSource
+)
 
 
 bool Foam::fv::volumeFractionSource::movePoints()

src/fvModels/derived/volumeFractionSource/volumeFractionSource.H

@@ -120,11 +120,10 @@ class volumeFractionSource
 
             //- Add source terms to an equation
             template<class Type, class AlphaFieldType>
-            void addGeneralSup
+            void addGeneralSupType
             (
                 const AlphaFieldType& alpha,
-                fvMatrix<Type>&,
+                fvMatrix<Type>&
-                const word& fieldName
             ) const;
 
             //- Add a source term to an equation
@@ -132,8 +131,8 @@ class volumeFractionSource
             void addAlphaSupType
             (
                 const AlphaFieldType& alpha,
-                fvMatrix<Type>& eqn,
+                const VolField<Type>& field,
-                const word& fieldName
+                fvMatrix<Type>& eqn
             ) const;
 
             //- Add a source term to a scalar equation
@@ -141,8 +140,8 @@ class volumeFractionSource
             void addAlphaSupType
             (
                 const AlphaFieldType& alpha,
-                fvMatrix<scalar>& eqn,
+                const volScalarField& field,
-                const word& fieldName
+                fvMatrix<scalar>& eqn
             ) const;
 
             //- Add a source term to a vector equation
@@ -150,21 +149,25 @@ class volumeFractionSource
             void addAlphaSupType
             (
                 const AlphaFieldType& alpha,
-                fvMatrix<vector>& eqn,
+                const volVectorField& field,
-                const word& fieldName
+                fvMatrix<vector>& eqn
             ) const;
 
             //- Add a source term to an equation
             template<class Type>
-            void addSupType(fvMatrix<Type>& eqn, const word& fieldName) const;
+            void addSupType
+            (
+                const VolField<Type>& field,
+                fvMatrix<Type>& eqn
+            ) const;
 
             //- Add a source term to a compressible equation
             template<class Type>
             void addSupType
             (
                 const volScalarField& rho,
-                fvMatrix<Type>& eqn,
+                const VolField<Type>& field,
-                const word& fieldName
+                fvMatrix<Type>& eqn
             ) const;
 
             //- Add a source term to a phase equation
@@ -173,13 +176,11 @@ class volumeFractionSource
             (
                 const volScalarField& alpha,
                 const volScalarField& rho,
-                fvMatrix<Type>& eqn,
+                const VolField<Type>& field,
-                const word& fieldName
+                fvMatrix<Type>& eqn
             ) const;
 
 
-
-
 public:
 
     //- Runtime type information
@@ -221,13 +222,13 @@ public:
         // Sources
 
             //- Add a source term to an equation
-            FOR_ALL_FIELD_TYPES(DEFINE_FV_MODEL_ADD_SUP);
+            FOR_ALL_FIELD_TYPES(DEFINE_FV_MODEL_ADD_FIELD_SUP)
 
             //- Add a source term to a compressible equation
-            FOR_ALL_FIELD_TYPES(DEFINE_FV_MODEL_ADD_RHO_SUP);
+            FOR_ALL_FIELD_TYPES(DEFINE_FV_MODEL_ADD_RHO_FIELD_SUP)
 
             //- Add a source term to a phase equation
-            FOR_ALL_FIELD_TYPES(DEFINE_FV_MODEL_ADD_ALPHA_RHO_SUP);
+            FOR_ALL_FIELD_TYPES(DEFINE_FV_MODEL_ADD_ALPHA_RHO_FIELD_SUP)
 
 
         // Mesh changes

src/fvModels/general/codedFvModel/codedFvModel.C

@@ -133,7 +133,6 @@ const Foam::dictionary& Foam::fv::codedFvModel::codeDict() const
 
 Foam::wordList Foam::fv::codedFvModel::codeKeys() const
 {
-
     return
     {
         "codeAddSup",
@@ -165,8 +164,8 @@ Foam::fvModel& Foam::fv::codedFvModel::redirectFvModel() const
 template<class Type>
 void Foam::fv::codedFvModel::addSupType
 (
-    fvMatrix<Type>& eqn,
+    const VolField<Type>& field,
-    const word& fieldName
+    fvMatrix<Type>& eqn
 ) const
 {
     if (fieldPrimitiveTypeName() != word::null)
@@ -177,7 +176,7 @@ void Foam::fv::codedFvModel::addSupType
         }
 
         updateLibrary();
-        redirectFvModel().addSup(eqn, fieldName);
+        redirectFvModel().addSup(field, eqn);
     }
 }
 
@@ -186,8 +185,8 @@ template<class Type>
 void Foam::fv::codedFvModel::addSupType
 (
     const volScalarField& rho,
-    fvMatrix<Type>& eqn,
+    const VolField<Type>& field,
-    const word& fieldName
+    fvMatrix<Type>& eqn
 ) const
 {
     if (fieldPrimitiveTypeName() != word::null)
@@ -198,7 +197,7 @@ void Foam::fv::codedFvModel::addSupType
         }
 
         updateLibrary();
-        redirectFvModel().addSup(rho, eqn, fieldName);
+        redirectFvModel().addSup(rho, field, eqn);
     }
 }
 
@@ -208,8 +207,8 @@ void Foam::fv::codedFvModel::addSupType
 (
     const volScalarField& alpha,
     const volScalarField& rho,
-    fvMatrix<Type>& eqn,
+    const VolField<Type>& field,
-    const word& fieldName
+    fvMatrix<Type>& eqn
 ) const
 {
     if (fieldPrimitiveTypeName() != word::null)
@@ -220,7 +219,7 @@ void Foam::fv::codedFvModel::addSupType
         }
 
         updateLibrary();
-        redirectFvModel().addSup(alpha, rho, eqn, fieldName);
+        redirectFvModel().addSup(alpha, rho, field, eqn);
     }
 }
 
@@ -250,13 +249,17 @@ Foam::wordList Foam::fv::codedFvModel::addSupFields() const
 }
 
 
-FOR_ALL_FIELD_TYPES(IMPLEMENT_FV_MODEL_ADD_SUP, fv::codedFvModel);
+FOR_ALL_FIELD_TYPES(IMPLEMENT_FV_MODEL_ADD_FIELD_SUP, fv::codedFvModel)
 
 
-FOR_ALL_FIELD_TYPES(IMPLEMENT_FV_MODEL_ADD_RHO_SUP, fv::codedFvModel);
+FOR_ALL_FIELD_TYPES(IMPLEMENT_FV_MODEL_ADD_RHO_FIELD_SUP, fv::codedFvModel)
 
 
-FOR_ALL_FIELD_TYPES(IMPLEMENT_FV_MODEL_ADD_ALPHA_RHO_SUP, fv::codedFvModel);
+FOR_ALL_FIELD_TYPES
+(
+    IMPLEMENT_FV_MODEL_ADD_ALPHA_RHO_FIELD_SUP,
+    fv::codedFvModel
+)
 
 
 bool Foam::fv::codedFvModel::movePoints()

src/fvModels/general/codedFvModel/codedFvModel.H

@@ -136,8 +136,8 @@ class codedFvModel
             template<class Type>
             void addSupType
             (
-                fvMatrix<Type>& eqn,
+                const VolField<Type>& field,
-                const word& fieldName
+                fvMatrix<Type>& eqn
             ) const;
 
             //- Add a source term to a compressible equation
@@ -145,8 +145,8 @@ class codedFvModel
             void addSupType
             (
                 const volScalarField& rho,
-                fvMatrix<Type>& eqn,
+                const VolField<Type>& field,
-                const word& fieldName
+                fvMatrix<Type>& eqn
             ) const;
 
             //- Add a source term to a phase equation
@@ -155,8 +155,8 @@ class codedFvModel
             (
                 const volScalarField& alpha,
                 const volScalarField& rho,
-                fvMatrix<Type>& eqn,
+                const VolField<Type>& field,
-                const word& fieldName
+                fvMatrix<Type>& eqn
             ) const;
 
 
@@ -190,13 +190,13 @@ public:
         // Sources
 
             //- Add a source term to an equation
-            FOR_ALL_FIELD_TYPES(DEFINE_FV_MODEL_ADD_SUP);
+            FOR_ALL_FIELD_TYPES(DEFINE_FV_MODEL_ADD_FIELD_SUP)
 
             //- Add a source term to a compressible equation
-            FOR_ALL_FIELD_TYPES(DEFINE_FV_MODEL_ADD_RHO_SUP);
+            FOR_ALL_FIELD_TYPES(DEFINE_FV_MODEL_ADD_RHO_FIELD_SUP)
 
             //- Add a source term to a phase equation
-            FOR_ALL_FIELD_TYPES(DEFINE_FV_MODEL_ADD_ALPHA_RHO_SUP);
+            FOR_ALL_FIELD_TYPES(DEFINE_FV_MODEL_ADD_ALPHA_RHO_FIELD_SUP)
 
 
         // Mesh changes

src/fvModels/general/semiImplicitSource/semiImplicitSource.C

@@ -86,16 +86,10 @@ void Foam::fv::semiImplicitSource::readCoeffs()
 template<class Type>
 void Foam::fv::semiImplicitSource::addSupType
 (
-    fvMatrix<Type>& eqn,
+    const VolField<Type>& field,
-    const word& fieldName
+    fvMatrix<Type>& eqn
 ) const
 {
-    if (debug)
-    {
-        Info<< "semiImplicitSource<" << pTraits<Type>::typeName
-            << ">::addSup for source " << name() << endl;
-    }
-
     const scalar t = mesh().time().userTimeValue();
 
     const VolField<Type>& psi = eqn.psi();
@@ -104,7 +98,7 @@ void Foam::fv::semiImplicitSource::addSupType
     (
         IOobject
         (
-            name() + fieldName + "Su",
+            name() + field.name() + "Su",
             mesh().time().name(),
             mesh(),
             IOobject::NO_READ,
@@ -134,13 +128,13 @@ void Foam::fv::semiImplicitSource::addSupType
 
     // Explicit source function for the field
     UIndirectList<Type>(Su, set_.cells()) =
-        fieldSu_[fieldName]->value<Type>(t)/VDash;
+        fieldSu_[field.name()]->template value<Type>(t)/VDash;
 
     volScalarField::Internal Sp
     (
         IOobject
         (
-            name() + fieldName + "Sp",
+            name() + field.name() + "Sp",
             mesh().time().name(),
             mesh(),
             IOobject::NO_READ,
@@ -158,7 +152,7 @@ void Foam::fv::semiImplicitSource::addSupType
 
     // Implicit source function for the field
     UIndirectList<scalar>(Sp, set_.cells()) =
-        fieldSp_[fieldName]->value(t)/VDash;
+        fieldSp_[field.name()]->value(t)/VDash;
 
     eqn += Su - fvm::SuSp(-Sp, psi);
 }
@@ -168,11 +162,11 @@ template<class Type>
 void Foam::fv::semiImplicitSource::addSupType
 (
     const volScalarField& rho,
-    fvMatrix<Type>& eqn,
+    const VolField<Type>& field,
-    const word& fieldName
+    fvMatrix<Type>& eqn
 ) const
 {
-    return this->addSup(eqn, fieldName);
+    return addSup(field, eqn);
 }
 
 
@@ -181,11 +175,11 @@ void Foam::fv::semiImplicitSource::addSupType
 (
     const volScalarField& alpha,
     const volScalarField& rho,
-    fvMatrix<Type>& eqn,
+    const VolField<Type>& field,
-    const word& fieldName
+    fvMatrix<Type>& eqn
 ) const
 {
-    return this->addSup(eqn, fieldName);
+    return addSup(field, eqn);
 }
 
 
@@ -221,17 +215,25 @@ Foam::wordList Foam::fv::semiImplicitSource::addSupFields() const
 }
 
 
-FOR_ALL_FIELD_TYPES(IMPLEMENT_FV_MODEL_ADD_SUP, fv::semiImplicitSource);
+FOR_ALL_FIELD_TYPES
-
+(
-
+    IMPLEMENT_FV_MODEL_ADD_FIELD_SUP,
-FOR_ALL_FIELD_TYPES(IMPLEMENT_FV_MODEL_ADD_RHO_SUP, fv::semiImplicitSource);
+    fv::semiImplicitSource
+)
 
 
 FOR_ALL_FIELD_TYPES
 (
-    IMPLEMENT_FV_MODEL_ADD_ALPHA_RHO_SUP,
+    IMPLEMENT_FV_MODEL_ADD_RHO_FIELD_SUP,
     fv::semiImplicitSource
-);
+)
+
+
+FOR_ALL_FIELD_TYPES
+(
+    IMPLEMENT_FV_MODEL_ADD_ALPHA_RHO_FIELD_SUP,
+    fv::semiImplicitSource
+)
 
 
 bool Foam::fv::semiImplicitSource::movePoints()

src/fvModels/general/semiImplicitSource/semiImplicitSource.H

@@ -154,15 +154,19 @@ private:
 
             //- Add a source term to an equation
             template <class Type>
-            void addSupType(fvMatrix<Type>& eqn, const word& fieldName) const;
+            void addSupType
+            (
+                const VolField<Type>& field,
+                fvMatrix<Type>& eqn
+            ) const;
 
             //- Add a source term to a compressible equation
             template <class Type>
             void addSupType
             (
                 const volScalarField& rho,
-                fvMatrix<Type>& eqn,
+                const VolField<Type>& field,
-                const word& fieldName
+                fvMatrix<Type>& eqn
             ) const;
 
             //- Add a source term to a phase equation
@@ -171,8 +175,8 @@ private:
             (
                 const volScalarField& alpha,
                 const volScalarField& rho,
-                fvMatrix<Type>& eqn,
+                const VolField<Type>& field,
-                const word& fieldName
+                fvMatrix<Type>& eqn
             ) const;
 
 
@@ -209,13 +213,13 @@ public:
         // Sources
 
             //- Add a source term to an equation
-            FOR_ALL_FIELD_TYPES(DEFINE_FV_MODEL_ADD_SUP);
+            FOR_ALL_FIELD_TYPES(DEFINE_FV_MODEL_ADD_FIELD_SUP)
 
             //- Add a source term to a compressible equation
-            FOR_ALL_FIELD_TYPES(DEFINE_FV_MODEL_ADD_RHO_SUP);
+            FOR_ALL_FIELD_TYPES(DEFINE_FV_MODEL_ADD_RHO_FIELD_SUP)
 
             //- Add a source term to a phase equation
-            FOR_ALL_FIELD_TYPES(DEFINE_FV_MODEL_ADD_ALPHA_RHO_SUP);
+            FOR_ALL_FIELD_TYPES(DEFINE_FV_MODEL_ADD_ALPHA_RHO_FIELD_SUP)
 
 
         // Mesh changes

src/fvModels/interRegion/interRegionExplicitPorositySource/interRegionExplicitPorositySource.C

@@ -152,9 +152,8 @@ Foam::fv::interRegionExplicitPorositySource::addSupFields() const
 
 void Foam::fv::interRegionExplicitPorositySource::addSup
 (
-    const volScalarField& rho,
+    const volVectorField& U,
-    fvMatrix<vector>& eqn,
+    fvMatrix<vector>& eqn
-    const word& fieldName
 ) const
 {
     fvMatrix<vector> porosityEqn(eqn.psi(), eqn.dimensions());
@@ -163,6 +162,33 @@ void Foam::fv::interRegionExplicitPorositySource::addSup
 }
 
 
+void Foam::fv::interRegionExplicitPorositySource::addSup
+(
+    const volScalarField& rho,
+    const volVectorField& U,
+    fvMatrix<vector>& eqn
+) const
+{
+    fvMatrix<vector> porosityEqn(eqn.psi(), eqn.dimensions());
+    porosityPtr_->addResistance(porosityEqn);
+    eqn -= filter_*porosityEqn;
+}
+
+
+void Foam::fv::interRegionExplicitPorositySource::addSup
+(
+    const volScalarField& alpha,
+    const volScalarField& rho,
+    const volVectorField& U,
+    fvMatrix<vector>& eqn
+) const
+{
+    fvMatrix<vector> porosityEqn(eqn.psi(), eqn.dimensions());
+    porosityPtr_->addResistance(porosityEqn);
+    eqn -= alpha*filter_*porosityEqn;
+}
+
+
 bool Foam::fv::interRegionExplicitPorositySource::movePoints()
 {
     NotImplemented;

src/fvModels/interRegion/interRegionExplicitPorositySource/interRegionExplicitPorositySource.H

@@ -139,13 +139,30 @@ public:
             virtual wordList addSupFields() const;
 
 
-        // Add explicit and implicit contributions to compressible equation
+        // Evaluation
 
+            //- Add source to momentum equation
+            virtual void addSup
+            (
+                const volVectorField& U,
+                fvMatrix<vector>& eqn
+            ) const;
+
+            //- Add source to compressible momentum equation
             virtual void addSup
             (
                 const volScalarField& rho,
-                fvMatrix<vector>& eqn,
+                const volVectorField& U,
-                const word& fieldName
+                fvMatrix<vector>& eqn
+            ) const;
+
+            //- Add source to phase momentum equation
+            virtual void addSup
+            (
+                const volScalarField& alpha,
+                const volScalarField& rho,
+                const volVectorField& U,
+                fvMatrix<vector>& eqn
             ) const;
 
 

src/fvModels/interRegion/interRegionHeatTransfer/interRegionHeatTransfer.C

@@ -114,12 +114,10 @@ Foam::wordList Foam::fv::interRegionHeatTransfer::addSupFields() const
 
 void Foam::fv::interRegionHeatTransfer::addSup
 (
-    fvMatrix<scalar>& eqn,
+    const volScalarField& he,
-    const word& fieldName
+    fvMatrix<scalar>& eqn
 ) const
 {
-    const volScalarField& he = eqn.psi();
-
     const volScalarField& T =
         mesh().lookupObject<volScalarField>(TName_);
 
@@ -199,11 +197,11 @@ void Foam::fv::interRegionHeatTransfer::addSup
 void Foam::fv::interRegionHeatTransfer::addSup
 (
     const volScalarField& rho,
-    fvMatrix<scalar>& eqn,
+    const volScalarField& he,
-    const word& fieldName
+    fvMatrix<scalar>& eqn
 ) const
 {
-    addSup(eqn, fieldName);
+    addSup(he, eqn);
 }
 
 

src/fvModels/interRegion/interRegionHeatTransfer/interRegionHeatTransfer.H

@@ -111,7 +111,7 @@ class interRegionHeatTransfer
 
     // Private member functions
 
-        //- Non-virtual read
+        //- Non-virtual readalpha*
         void readCoeffs();
 
         //- Get the neighbour heat transfer
@@ -160,16 +160,16 @@ public:
             //- Source term to energy equation
             virtual void addSup
             (
-                fvMatrix<scalar>& eqn,
+                const volScalarField& he,
-                const word& fieldName
+                fvMatrix<scalar>& eqn
             ) const;
 
             //- Source term to compressible energy equation
             virtual void addSup
             (
                 const volScalarField& rho,
-                fvMatrix<scalar>& eqn,
+                const volScalarField& he,
-                const word& fieldName
+                fvMatrix<scalar>& eqn
             ) const;
 
 

src/lagrangian/parcel/fvModels/clouds/clouds.C

@@ -172,13 +172,11 @@ Foam::wordList Foam::fv::clouds::addSupFields() const
             const multicomponentThermo& carrierMcThermo =
                 refCast<const multicomponentThermo>(carrierThermo);
 
-            const PtrList<volScalarField>& Y = carrierMcThermo.Y();
+            forAll(carrierMcThermo.Y(), i)
-
-            forAll(Y, i)
             {
                 if (carrierMcThermo.solveSpecie(i))
                 {
-                    fieldNames.append(Y[i].name());
+                    fieldNames.append(carrierMcThermo.Y()[i].name());
                 }
             }
         }
@@ -208,8 +206,8 @@ void Foam::fv::clouds::correct()
 
 void Foam::fv::clouds::addSup
 (
-    fvMatrix<scalar>& eqn,
+    const volScalarField& rho,
-    const word& fieldName
+    fvMatrix<scalar>& eqn
 ) const
 {
     if (debug)
@@ -225,18 +223,24 @@ void Foam::fv::clouds::addSup
             << exit(FatalError);
     }
 
-    if (fieldName == rhoName_)
+    if (rho.name() == rhoName_)
     {
         eqn += cloudsPtr_().Srho(eqn.psi());
     }
+    else
+    {
+        FatalErrorInFunction
+            << "Support for field " << rho.name() << " is not implemented"
+            << exit(FatalError);
+    }
 }
 
 
 void Foam::fv::clouds::addSup
 (
     const volScalarField& rho,
-    fvMatrix<scalar>& eqn,
+    const volScalarField& field,
-    const word& fieldName
+    fvMatrix<scalar>& eqn
 ) const
 {
     if (debug)
@@ -254,36 +258,41 @@ void Foam::fv::clouds::addSup
 
     const fluidThermo& carrierThermo = tCarrierThermo_();
 
-    if (fieldName == rhoName_)
+    if (&field == &carrierThermo.he())
-    {
-        eqn += cloudsPtr_().Srho(eqn.psi());
-    }
-    else if (fieldName == carrierThermo.he().name())
     {
         eqn += cloudsPtr_().Sh(eqn.psi());
+        return;
     }
-    else if (isA<multicomponentThermo>(carrierThermo))
+
+    if (isA<multicomponentThermo>(carrierThermo))
     {
         const multicomponentThermo& carrierMcThermo =
             refCast<const multicomponentThermo>(carrierThermo);
 
-        if (carrierMcThermo.containsSpecie(eqn.psi().name()))
+        if (carrierMcThermo.containsSpecie(field.name()))
         {
             eqn +=
                 cloudsPtr_().SYi
                 (
-                    carrierMcThermo.specieIndex(eqn.psi()),
+                    carrierMcThermo.specieIndex(field),
-                    eqn.psi()
+                    field
                 );
+            return;
         }
     }
+
+    {
+        FatalErrorInFunction
+            << "Support for field " << field.name() << " is not implemented"
+            << exit(FatalError);
+    }
 }
 
 
 void Foam::fv::clouds::addSup
 (
-    fvMatrix<vector>& eqn,
+    const volVectorField& U,
-    const word& fieldName
+    fvMatrix<vector>& eqn
 ) const
 {
     if (debug)
@@ -299,18 +308,24 @@ void Foam::fv::clouds::addSup
             << exit(FatalError);
     }
 
-    if (fieldName == UName_)
+    if (U.name() == UName_)
     {
         eqn += cloudsPtr_().SU(eqn.psi())/tRho_();
     }
+    else
+    {
+        FatalErrorInFunction
+            << "Support for field " << U.name() << " is not implemented"
+            << exit(FatalError);
+    }
 }
 
 
 void Foam::fv::clouds::addSup
 (
     const volScalarField& rho,
-    fvMatrix<vector>& eqn,
+    const volVectorField& U,
-    const word& fieldName
+    fvMatrix<vector>& eqn
 ) const
 {
     if (debug)
@@ -326,10 +341,16 @@ void Foam::fv::clouds::addSup
             << exit(FatalError);
     }
 
-    if (fieldName == UName_)
+    if (U.name() == UName_)
     {
         eqn += cloudsPtr_().SU(eqn.psi());
     }
+    else
+    {
+        FatalErrorInFunction
+            << "Support for field " << U.name() << " is not implemented"
+            << exit(FatalError);
+    }
 }
 
 

src/lagrangian/parcel/fvModels/clouds/clouds.H

@@ -177,31 +177,31 @@ public:
             //- Add source to continuity equation
             virtual void addSup
             (
-                fvMatrix<scalar>& eqn,
+                const volScalarField& rho,
-                const word& fieldName
+                fvMatrix<scalar>& eqn
             ) const;
 
             //- Add source to enthalpy or species equation
             virtual void addSup
             (
                 const volScalarField& rho,
-                fvMatrix<scalar>& eqn,
+                const volScalarField& field,
-                const word& fieldName
+                fvMatrix<scalar>& eqn
             ) const;
 
             //- Add source to incompressible momentum equation
             virtual void addSup
             (
-                fvMatrix<vector>& eqn,
+                const volVectorField& U,
-                const word& fieldName
+                fvMatrix<vector>& eqn
             ) const;
 
             //- Add source to compressible momentum equation
             virtual void addSup
             (
                 const volScalarField& rho,
-                fvMatrix<vector>& eqn,
+                const volVectorField& U,
-                const word& fieldName
+                fvMatrix<vector>& eqn
             ) const;