diff --git a/applications/legacy/combustion/PDRFoam/pEqn.H b/applications/legacy/combustion/PDRFoam/pEqn.H
index 6239b7f8b4..d1db2cb1b2 100644
--- a/applications/legacy/combustion/PDRFoam/pEqn.H
+++ b/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();
diff --git a/applications/legacy/compressible/rhoPorousSimpleFoam/pEqn.H b/applications/legacy/compressible/rhoPorousSimpleFoam/pEqn.H
index 71e27f12ba..0679f096de 100644
--- a/applications/legacy/compressible/rhoPorousSimpleFoam/pEqn.H
+++ b/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)
             );
diff --git a/applications/modules/compressibleMultiphaseVoF/pressureCorrector.C b/applications/modules/compressibleMultiphaseVoF/pressureCorrector.C
index 1247a8815c..febb53bc0e 100644
--- a/applications/modules/compressibleMultiphaseVoF/pressureCorrector.C
+++ b/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()
             );
         }
diff --git a/applications/modules/compressibleVoF/alphaSuSp.C b/applications/modules/compressibleVoF/alphaSuSp.C
index 98b63a7c83..2ea27cd564 100644
--- a/applications/modules/compressibleVoF/alphaSuSp.C
+++ b/applications/modules/compressibleVoF/alphaSuSp.C
@@ -30,46 +30,52 @@ License
 
 void Foam::solvers::compressibleVoF::alphaSuSp
 (
-    tmp<volScalarField::Internal>& Su,
-    tmp<volScalarField::Internal>& Sp
+    tmp<volScalarField::Internal>& tSu,
+    tmp<volScalarField::Internal>& tSp
 )
 {
-    Sp = volScalarField::Internal::New
-    (
-        "Sp",
-        mesh,
-        dimensionedScalar(dgdt.dimensions(), 0)
-    );
+    const dimensionedScalar Szero(dgdt.dimensions(), 0);
 
-    Su = volScalarField::Internal::New
-    (
-        "Su",
-        mesh,
-        dimensionedScalar(dgdt.dimensions(), 0)
-    );
+    tSp = volScalarField::Internal::New("Sp", mesh, Szero);
+    tSu = volScalarField::Internal::New("Su", mesh, Szero);
 
-    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 fvScalarMatrix alphaSup(fvModels().source(alpha1));
+        const volScalarField::Internal alpha2ByRho1(alpha2()/mixture.rho1()());
+        const fvScalarMatrix alphaRho1Sup
+        (
+            fvModels().sourceProxy(alpha1, mixture.rho1(), alpha1)
+        );
 
-        Su = alphaSup.Su();
-        Sp = alphaSup.Sp();
+        Su += alpha2ByRho1*alphaRho1Sup.Su();
+        Sp += alpha2ByRho1*alphaRho1Sup.Sp();
     }
 
-    volScalarField::Internal& SpRef = Sp.ref();
-    volScalarField::Internal& SuRef = Su.ref();
+    if (fvModels().addsSupToField(mixture.rho2().name()))
+    {
+        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);
-            SuRef[celli] += dgdt[celli]/max(1.0 - alpha1[celli], 1e-4);
+            Sp[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);
         }
     }
 }
diff --git a/applications/modules/compressibleVoF/compressibleVoF.H b/applications/modules/compressibleVoF/compressibleVoF.H
index 0a73013281..b56773ada6 100644
--- a/applications/modules/compressibleVoF/compressibleVoF.H
+++ b/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
diff --git a/applications/modules/compressibleVoF/fvModels/VoFCavitation/VoFCavitation.C b/applications/modules/compressibleVoF/fvModels/VoFCavitation/VoFCavitation.C
index 6f292c3532..1b232133f2 100644
--- a/applications/modules/compressibleVoF/fvModels/VoFCavitation/VoFCavitation.C
+++ b/applications/modules/compressibleVoF/fvModels/VoFCavitation/VoFCavitation.C
@@ -69,9 +69,7 @@ Foam::fv::compressible::VoFCavitation::VoFCavitation
         )
     ),
 
-    cavitation_(Foam::compressible::cavitationModel::New(dict, mixture_)),
-
-    alphaName_(mixture_.alpha1().name())
+    cavitation_(Foam::compressible::cavitationModel::New(dict, mixture_))
 {}
 
 
@@ -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 word& fieldName
+    const volScalarField& alpha,
+    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
-        (
-            1.0/mixture_.rho1()()
-          - mixture_.alpha1()()*(1.0/mixture_.rho1()() - 1.0/mixture_.rho2()())
-        );
+        const volScalarField& alpha1 = mixture_.alpha1();
+        const volScalarField& alpha2 = mixture_.alpha2();
 
-        const Pair<tmp<volScalarField::Internal>> mDot12Alpha
-        (
-            cavitation_->mDot12Alpha()
-        );
+        const scalar s = &alpha == &alpha1 ? +1 : -1;
 
-        const volScalarField::Internal vDot1Alpha(alpha1Coeff*mDot12Alpha[0]());
-        const volScalarField::Internal vDot2Alpha(alpha1Coeff*mDot12Alpha[1]());
+        // Volume-fraction linearisation
+        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);
+        }
     }
-}
-
-
-void Foam::fv::compressible::VoFCavitation::addSup
-(
-    const volScalarField&,
-    fvMatrix<scalar>& eqn,
-    const word& fieldName
-) const
-{
-    if (debug)
+    else
     {
-        Info<< type() << ": applying source to " << eqn.psi().name() << endl;
-    }
-
-    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());
+        FatalErrorInFunction
+            << "Support for field " << alpha.name() << " is not implemented"
+            << exit(FatalError);
     }
 }
 
diff --git a/applications/modules/compressibleVoF/fvModels/VoFCavitation/VoFCavitation.H b/applications/modules/compressibleVoF/fvModels/VoFCavitation/VoFCavitation.H
index 55843f00f6..4d941f6bbf 100644
--- a/applications/modules/compressibleVoF/fvModels/VoFCavitation/VoFCavitation.H
+++ b/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 word& fieldName
-        ) const;
-
-        //- Add implicit/explicit contributions to p_rgh equation
-        virtual void addSup
-        (
-            const volScalarField& psi,
-            fvMatrix<scalar>& eqn,
-            const word& fieldName
+            const volScalarField& alpha,
+            const volScalarField& rho,
+            fvMatrix<scalar>& eqn
         ) const;
 
 
diff --git a/applications/modules/compressibleVoF/fvModels/VoFClouds/VoFClouds.C b/applications/modules/compressibleVoF/fvModels/VoFClouds/VoFClouds.C
index 74d376f4d2..5d51fc36e4 100644
--- a/applications/modules/compressibleVoF/fvModels/VoFClouds/VoFClouds.C
+++ b/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 word& fieldName
+    const volScalarField& rho,
+    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 word& fieldName
+    const volScalarField& he,
+    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 word& fieldName
+    const volVectorField& U,
+    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);
+    }
 }
 
 
diff --git a/applications/modules/compressibleVoF/fvModels/VoFClouds/VoFClouds.H b/applications/modules/compressibleVoF/fvModels/VoFClouds/VoFClouds.H
index bfebbd2247..26956b6a25 100644
--- a/applications/modules/compressibleVoF/fvModels/VoFClouds/VoFClouds.H
+++ b/applications/modules/compressibleVoF/fvModels/VoFClouds/VoFClouds.H
@@ -130,8 +130,8 @@ public:
             virtual void addSup
             (
                 const volScalarField& alpha,
-                fvMatrix<scalar>& eqn,
-                const word& fieldName
+                const volScalarField& rho,
+                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 word& fieldName
+                const volScalarField& he,
+                fvMatrix<scalar>& eqn
             ) const;
 
             //- Add implicit contribution to mixture momentum equation
             virtual void addSup
             (
                 const volScalarField& rho,
-                fvMatrix<vector>& eqn,
-                const word& fieldName
+                const volVectorField& U,
+                fvMatrix<vector>& eqn
             ) const;
 
 
diff --git a/applications/modules/compressibleVoF/fvModels/VoFSolidificationMeltingSource/VoFSolidificationMeltingSource.C b/applications/modules/compressibleVoF/fvModels/VoFSolidificationMeltingSource/VoFSolidificationMeltingSource.C
index dc01a7ebe5..95732accb5 100644
--- a/applications/modules/compressibleVoF/fvModels/VoFSolidificationMeltingSource/VoFSolidificationMeltingSource.C
+++ b/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 word& fieldName
+    const volScalarField& he,
+    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 word& fieldName
+    const volVectorField& U,
+    fvMatrix<vector>& eqn
 ) const
 {
     if (debug)
diff --git a/applications/modules/compressibleVoF/fvModels/VoFSolidificationMeltingSource/VoFSolidificationMeltingSource.H b/applications/modules/compressibleVoF/fvModels/VoFSolidificationMeltingSource/VoFSolidificationMeltingSource.H
index 05295ba963..a880d456a1 100644
--- a/applications/modules/compressibleVoF/fvModels/VoFSolidificationMeltingSource/VoFSolidificationMeltingSource.H
+++ b/applications/modules/compressibleVoF/fvModels/VoFSolidificationMeltingSource/VoFSolidificationMeltingSource.H
@@ -182,16 +182,16 @@ public:
             (
                 const volScalarField& alpha,
                 const volScalarField& rho,
-                fvMatrix<scalar>& eqn,
-                const word& fieldName
+                const volScalarField& he,
+                fvMatrix<scalar>& eqn
             ) const;
 
             //- Add implicit contribution to mixture momentum equation
             virtual void addSup
             (
                 const volScalarField& rho,
-                fvMatrix<vector>& eqn,
-                const word& fieldName
+                const volVectorField& U,
+                fvMatrix<vector>& eqn
             ) const;
 
 
diff --git a/applications/modules/compressibleVoF/fvModels/VoFTurbulenceDamping/VoFTurbulenceDamping.C b/applications/modules/compressibleVoF/fvModels/VoFTurbulenceDamping/VoFTurbulenceDamping.C
index d2099bb1cf..ef2a7eae3a 100644
--- a/applications/modules/compressibleVoF/fvModels/VoFTurbulenceDamping/VoFTurbulenceDamping.C
+++ b/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 word& fieldName
+    const volScalarField& field,
+    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&,
-    fvMatrix<scalar>& eqn,
-    const word& fieldName
+    const volScalarField& rho,
+    const volScalarField& field,
+    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);
     }
 }
diff --git a/applications/modules/compressibleVoF/fvModels/VoFTurbulenceDamping/VoFTurbulenceDamping.H b/applications/modules/compressibleVoF/fvModels/VoFTurbulenceDamping/VoFTurbulenceDamping.H
index 78028b08ff..a53f83bf11 100644
--- a/applications/modules/compressibleVoF/fvModels/VoFTurbulenceDamping/VoFTurbulenceDamping.H
+++ b/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
-        //  to the transport equation
-        virtual wordList addSupFields() const;
+        // Checks
 
-        //- Add explicit contribution to epsilon or omega equation
-        virtual void addSup
-        (
-            const volScalarField& rho,
-            fvMatrix<scalar>& eqn,
-            const word& fieldName
-        ) const;
+            //- Return the list of fields for which the option adds source term
+            //  to the transport equation
+            virtual wordList addSupFields() const;
 
-        //- Add explicit contribution to phase epsilon or omega equation
-        virtual void addSup
-        (
-            const volScalarField& alpha,
-            const volScalarField& rho,
-            fvMatrix<scalar>& eqn,
-            const word& fieldName
-        ) const;
+
+        // Sources
+
+            //- Add explicit contribution to epsilon or omega equation
+            virtual void addSup
+            (
+                const volScalarField& rho,
+                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
diff --git a/applications/modules/compressibleVoF/pressureCorrector.C b/applications/modules/compressibleVoF/pressureCorrector.C
index 4987549c14..c953233bb0 100644
--- a/applications/modules/compressibleVoF/pressureCorrector.C
+++ b/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
-        fvScalarMatrix Sp_rgh
+        // Cache any sources
+        fvScalarMatrix p_rghEqnSource
         (
-            fvModels().source
-            (
-                volScalarField::New
-                (
-                    "1",
-                    mesh,
-                    dimensionedScalar(dimless/dimPressure, 1)
-                ),
-                p_rgh
-            )
+            fvModels().sourceProxy(alpha1, rho1, p_rgh)/rho1
+          + fvModels().sourceProxy(alpha2, rho2, p_rgh)/rho2
         );
 
         // 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
             );
 
             {
diff --git a/applications/modules/incompressibleDenseParticleFluid/correctPressure.C b/applications/modules/incompressibleDenseParticleFluid/correctPressure.C
index 367b87ccea..d9cc2caed1 100644
--- a/applications/modules/incompressibleDenseParticleFluid/correctPressure.C
+++ b/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
diff --git a/applications/modules/incompressibleDriftFlux/alphaSuSp.C b/applications/modules/incompressibleDriftFlux/alphaSuSp.C
index f715ac0885..ceac90431b 100644
--- a/applications/modules/incompressibleDriftFlux/alphaSuSp.C
+++ b/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>& Sp
+    tmp<volScalarField::Internal>& tSu,
+    tmp<volScalarField::Internal>& tSp
 )
 {
-    if (divergent())
-    {
-        // Phase change alpha1 source
-        const fvScalarMatrix alphaSup(fvModels().source(alpha1));
+    if (!divergent()) return;
 
-        Su = alphaSup.Su();
-        Sp = alphaSup.Sp();
+    const dimensionedScalar Szero(dimless/dimTime, 0);
+
+    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();
     }
 }
 
diff --git a/applications/modules/incompressibleDriftFlux/incompressibleDriftFlux.H b/applications/modules/incompressibleDriftFlux/incompressibleDriftFlux.H
index a302e1fa0d..e83dae3f4e 100644
--- a/applications/modules/incompressibleDriftFlux/incompressibleDriftFlux.H
+++ b/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
diff --git a/applications/modules/incompressibleFluid/correctPressure.C b/applications/modules/incompressibleFluid/correctPressure.C
index 7c6a01e305..a2779d908f 100644
--- a/applications/modules/incompressibleFluid/correctPressure.C
+++ b/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
diff --git a/applications/modules/incompressibleMultiphaseVoF/pressureCorrector.C b/applications/modules/incompressibleMultiphaseVoF/pressureCorrector.C
index 2487207c70..ece2b135bc 100644
--- a/applications/modules/incompressibleMultiphaseVoF/pressureCorrector.C
+++ b/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
-        fvScalarMatrix Sp_rgh
-        (
-            fvModels().source
-            (
-                volScalarField::New
-                (
-                    "1",
-                    mesh,
-                    dimensionedScalar(dimless/dimPressure, 1)
-                ),
-                p_rgh
-            )
-        );
+        // Evaluate any phase sources
+        fvScalarMatrix p_rghEqnSource(p_rgh, dimVolume/dimTime);
+        forAll(phases, phasei)
+        {
+            p_rghEqnSource +=
+                fvModels().sourceProxy(phases[phasei], p_rgh);
+        }
 
         while (pimple.correctNonOrthogonal())
         {
             fvScalarMatrix p_rghEqn
             (
                 fvc::div(phiHbyA) - fvm::laplacian(rAUf, p_rgh)
-             == Sp_rgh
+             == p_rghEqnSource
             );
 
             p_rghEqn.setReference
diff --git a/applications/modules/incompressibleVoF/alphaSuSp.C b/applications/modules/incompressibleVoF/alphaSuSp.C
index ec8836b770..1099ff91e2 100644
--- a/applications/modules/incompressibleVoF/alphaSuSp.C
+++ b/applications/modules/incompressibleVoF/alphaSuSp.C
@@ -30,17 +30,34 @@ License
 
 void Foam::solvers::incompressibleVoF::alphaSuSp
 (
-    tmp<volScalarField::Internal>& Su,
-    tmp<volScalarField::Internal>& Sp
+    tmp<volScalarField::Internal>& tSu,
+    tmp<volScalarField::Internal>& tSp
 )
 {
-    if (divergent())
-    {
-        // Phase change alpha1 source
-        const fvScalarMatrix alphaSup(fvModels().source(alpha1));
+    if (!divergent()) return;
 
-        Su = alphaSup.Su();
-        Sp = alphaSup.Sp();
+    const dimensionedScalar Szero(dimless/dimTime, 0);
+
+    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();
     }
 }
 
diff --git a/applications/modules/incompressibleVoF/fvModels/VoFCavitation/VoFCavitation.C b/applications/modules/incompressibleVoF/fvModels/VoFCavitation/VoFCavitation.C
index 8d113b8868..d4a8491280 100644
--- a/applications/modules/incompressibleVoF/fvModels/VoFCavitation/VoFCavitation.C
+++ b/applications/modules/incompressibleVoF/fvModels/VoFCavitation/VoFCavitation.C
@@ -68,9 +68,7 @@ Foam::fv::VoFCavitation::VoFCavitation
         )
     ),
 
-    cavitation_(cavitationModel::New(dict, mixture_)),
-
-    alphaName_(mixture_.alpha1().name())
+    cavitation_(cavitationModel::New(dict, mixture_))
 {}
 
 
@@ -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 word& fieldName
+    const volScalarField& alpha,
+    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
-        (
-            1.0/mixture_.rho1()
-          - mixture_.alpha1()()*(1.0/mixture_.rho1() - 1.0/mixture_.rho2())
-        );
+        const volScalarField& alpha1 = mixture_.alpha1();
+        const volScalarField& alpha2 = mixture_.alpha2();
 
-        const Pair<tmp<volScalarField::Internal>> mDot12Alpha
-        (
-            cavitation_->mDot12Alpha()
-        );
+        const dimensionedScalar& rho =
+            &alpha == &alpha1 ? mixture_.rho1() : mixture_.rho2();
 
-        const volScalarField::Internal vDot1Alpha(alpha1Coeff*mDot12Alpha[0]());
-        const volScalarField::Internal vDot2Alpha(alpha1Coeff*mDot12Alpha[1]());
+        const scalar s = &alpha == &alpha1 ? +1 : -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);
+        }
     }
-}
-
-
-void Foam::fv::VoFCavitation::addSup
-(
-    const volScalarField&,
-    fvMatrix<scalar>& eqn,
-    const word& fieldName
-) const
-{
-    if (debug)
+    else
     {
-        Info<< type() << ": applying source to " << eqn.psi().name() << endl;
-    }
-
-    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());
+        FatalErrorInFunction
+            << "Support for field " << alpha.name() << " is not implemented"
+            << exit(FatalError);
     }
 }
 
@@ -157,8 +167,8 @@ void Foam::fv::VoFCavitation::addSup
 void Foam::fv::VoFCavitation::addSup
 (
     const volScalarField& rho,
-    fvMatrix<vector>& eqn,
-    const word& fieldName
+    const volVectorField& U,
+    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);
     }
 }
 
diff --git a/applications/modules/incompressibleVoF/fvModels/VoFCavitation/VoFCavitation.H b/applications/modules/incompressibleVoF/fvModels/VoFCavitation/VoFCavitation.H
index 1635753b3a..7dcc547769 100644
--- a/applications/modules/incompressibleVoF/fvModels/VoFCavitation/VoFCavitation.H
+++ b/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 word& fieldName
+            const volScalarField& alpha,
+            fvMatrix<scalar>& eqn
         ) const;
 
-        //- Add implicit/explicit contributions to p_rgh equation
-        virtual void addSup
-        (
-            const volScalarField& psi,
-            fvMatrix<scalar>& eqn,
-            const word& fieldName
-        ) const;
-
-        //- Add contribution to U equation
+        //- Add a source to the mixture momentum equation
         virtual void addSup
         (
             const volScalarField& rho,
-            fvMatrix<vector>& eqn,
-            const word& fieldName
+            const volVectorField& U,
+            fvMatrix<vector>& eqn
         ) const;
 
 
diff --git a/applications/modules/incompressibleVoF/fvModels/VoFTurbulenceDamping/VoFTurbulenceDamping.C b/applications/modules/incompressibleVoF/fvModels/VoFTurbulenceDamping/VoFTurbulenceDamping.C
index 5b638357f2..0cb22276db 100644
--- a/applications/modules/incompressibleVoF/fvModels/VoFTurbulenceDamping/VoFTurbulenceDamping.C
+++ b/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 word& fieldName
+    const volScalarField& field,
+    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&,
-    fvMatrix<scalar>& eqn,
-    const word& fieldName
+    const volScalarField& rho,
+    const volScalarField& field,
+    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);
     }
 }
diff --git a/applications/modules/incompressibleVoF/fvModels/VoFTurbulenceDamping/VoFTurbulenceDamping.H b/applications/modules/incompressibleVoF/fvModels/VoFTurbulenceDamping/VoFTurbulenceDamping.H
index 49b721dd65..97d057c5dc 100644
--- a/applications/modules/incompressibleVoF/fvModels/VoFTurbulenceDamping/VoFTurbulenceDamping.H
+++ b/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
-        //  to the transport equation
-        virtual wordList addSupFields() const;
+        // Checks
 
-        //- Add explicit contribution to epsilon or omega equation
-        virtual void addSup
-        (
-            fvMatrix<scalar>& eqn,
-            const word& fieldName
-        ) const;
+            //- Return the list of fields for which the option adds source term
+            //  to the transport equation
+            virtual wordList addSupFields() const;
 
-        //- Add explicit contribution to phase epsilon or omega equation
-        virtual void addSup
-        (
-            const volScalarField& alpha,
-            const volScalarField& rho,
-            fvMatrix<scalar>& eqn,
-            const word& fieldName
-        ) const;
+
+        // Sources
+
+            //- Add explicit contribution to epsilon or omega equation
+            virtual void addSup
+            (
+                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
diff --git a/applications/modules/incompressibleVoF/incompressibleVoF.H b/applications/modules/incompressibleVoF/incompressibleVoF.H
index 6c645b77dd..faa213689e 100644
--- a/applications/modules/incompressibleVoF/incompressibleVoF.H
+++ b/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
diff --git a/applications/modules/isothermalFilm/fvModels/filmCloudTransfer/filmCloudTransfer.C b/applications/modules/isothermalFilm/fvModels/filmCloudTransfer/filmCloudTransfer.C
index 7af2188488..1fd2d839b2 100644
--- a/applications/modules/isothermalFilm/fvModels/filmCloudTransfer/filmCloudTransfer.C
+++ b/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 word& fieldName
+    const volScalarField& alpha,
+    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 word& fieldName
+    const volScalarField& he,
+    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 word& fieldName
+    const volVectorField& U,
+    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);
     }
 }
diff --git a/applications/modules/isothermalFilm/fvModels/filmCloudTransfer/filmCloudTransfer.H b/applications/modules/isothermalFilm/fvModels/filmCloudTransfer/filmCloudTransfer.H
index a1e58e61b8..9e4fd8f624 100644
--- a/applications/modules/isothermalFilm/fvModels/filmCloudTransfer/filmCloudTransfer.H
+++ b/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,
-                const word& fieldName
+                fvMatrix<scalar>& eqn
             ) 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 word& fieldName
+                const volScalarField& he,
+                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 word& fieldName
+                const volVectorField& U,
+                fvMatrix<vector>& eqn
             ) const;
 
 
diff --git a/applications/modules/isothermalFilm/fvModels/filmVoFTransfer/VoFFilmTransfer.C b/applications/modules/isothermalFilm/fvModels/filmVoFTransfer/VoFFilmTransfer.C
index cac396cdf3..8f3f6a9bbd 100644
--- a/applications/modules/isothermalFilm/fvModels/filmVoFTransfer/VoFFilmTransfer.C
+++ b/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 word& fieldName
+    const volScalarField& rho,
+    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 word& fieldName
+    const volScalarField& he,
+    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 word& fieldName
+    const volVectorField& U,
+    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 +=
-        filmVoFTransferRate<vector>
-        (
-            &filmVoFTransfer::UTransferRate,
-            dimMass*dimVelocity
-        )
-      - fvm::Sp(alpha_()*thermo_.rho()()*transferRate_, eqn.psi());
+    if (&U == &VoF_.U)
+    {
+        // Explicit transfer into the VoF
+        eqn +=
+            filmVoFTransferRate<vector>
+            (
+                &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);
+    }
 }
 
 
diff --git a/applications/modules/isothermalFilm/fvModels/filmVoFTransfer/VoFFilmTransfer.H b/applications/modules/isothermalFilm/fvModels/filmVoFTransfer/VoFFilmTransfer.H
index d985161881..97165fa31c 100644
--- a/applications/modules/isothermalFilm/fvModels/filmVoFTransfer/VoFFilmTransfer.H
+++ b/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
-            virtual void addSup
-            (
-                fvMatrix<scalar>& eqn,
-                const word& fieldName
-            ) const;
-
-            //- Add implicit contribution to phase density equation
+            //- Add implicit contribution to phase continuity equation
             virtual void addSup
             (
                 const volScalarField& alpha,
-                fvMatrix<scalar>& eqn,
-                const word& fieldName
+                const volScalarField& rho,
+                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 word& fieldName
+                const volScalarField& he,
+                fvMatrix<scalar>& eqn
             ) const;
 
             //- Add implicit contribution to mixture momentum equation
             virtual void addSup
             (
                 const volScalarField& rho,
-                fvMatrix<vector>& eqn,
-                const word& fieldName
+                const volVectorField& U,
+                fvMatrix<vector>& eqn
             ) const;
 
 
diff --git a/applications/modules/isothermalFilm/fvModels/filmVoFTransfer/filmVoFTransfer.C b/applications/modules/isothermalFilm/fvModels/filmVoFTransfer/filmVoFTransfer.C
index 0c42ab0506..c73c45e897 100644
--- a/applications/modules/isothermalFilm/fvModels/filmVoFTransfer/filmVoFTransfer.C
+++ b/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 word& fieldName
+    const volScalarField& alpha,
+    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 word& fieldName
+    const volScalarField& he,
+    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 word& fieldName
+    const volVectorField& U,
+    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 +=
-        VoFToFilmTransferRate<vector>
-        (
-            &VoFFilmTransfer::UTransferRate,
-            dimMomentum
-        )
-      - fvm::Sp(alpha()*rho()*transferRate_, eqn.psi());
+    if (&U == &film_.U)
+    {
+        // Explicit transfer into the film
+        eqn +=
+            VoFToFilmTransferRate<vector>
+            (
+                &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);
+    }
 }
 
 
diff --git a/applications/modules/isothermalFilm/fvModels/filmVoFTransfer/filmVoFTransfer.H b/applications/modules/isothermalFilm/fvModels/filmVoFTransfer/filmVoFTransfer.H
index f9c8b9d375..f75bcd58ab 100644
--- a/applications/modules/isothermalFilm/fvModels/filmVoFTransfer/filmVoFTransfer.H
+++ b/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,
-                const word& fieldName
+                fvMatrix<scalar>& eqn
             ) 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 word& fieldName
+                const volScalarField& he,
+                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 word& fieldName
+                const volVectorField& U,
+                fvMatrix<vector>& eqn
             ) const;
 
 
diff --git a/applications/modules/isothermalFluid/correctBuoyantPressure.C b/applications/modules/isothermalFluid/correctBuoyantPressure.C
index 536cadb784..53fd2002d7 100644
--- a/applications/modules/isothermalFluid/correctBuoyantPressure.C
+++ b/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())
diff --git a/applications/modules/isothermalFluid/correctPressure.C b/applications/modules/isothermalFluid/correctPressure.C
index 558ea4823e..20c9190bbe 100644
--- a/applications/modules/isothermalFluid/correctPressure.C
+++ b/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())
diff --git a/applications/modules/multiphaseEuler/compressibilityEqns.C b/applications/modules/multiphaseEuler/compressibilityEqns.C
index d1352b75e3..806d573987 100644
--- a/applications/modules/multiphaseEuler/compressibilityEqns.C
+++ b/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
diff --git a/applications/modules/multiphaseEuler/fvModels/interfaceTurbulenceDamping/interfaceTurbulenceDamping.C b/applications/modules/multiphaseEuler/fvModels/interfaceTurbulenceDamping/interfaceTurbulenceDamping.C
index 1854990d68..cd496b5a2d 100644
--- a/applications/modules/multiphaseEuler/fvModels/interfaceTurbulenceDamping/interfaceTurbulenceDamping.C
+++ b/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 word& fieldName
+    const volScalarField& field,
+    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 word& fieldName
+    const volScalarField& field,
+    fvMatrix<scalar>& eqn
 ) const
 {
-    addRhoSup(one(), eqn, fieldName);
+    addRhoSup(one(), field, eqn);
 }
 
 
 void Foam::fv::interfaceTurbulenceDamping::addSup
 (
     const volScalarField& rho,
-    fvMatrix<scalar>& eqn,
-    const word& fieldName
+    const volScalarField& field,
+    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 word& fieldName
+    const volScalarField& field,
+    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);
     }
 }
diff --git a/applications/modules/multiphaseEuler/fvModels/interfaceTurbulenceDamping/interfaceTurbulenceDamping.H b/applications/modules/multiphaseEuler/fvModels/interfaceTurbulenceDamping/interfaceTurbulenceDamping.H
index 005e4b8898..4fecb3ff40 100644
--- a/applications/modules/multiphaseEuler/fvModels/interfaceTurbulenceDamping/interfaceTurbulenceDamping.H
+++ b/applications/modules/multiphaseEuler/fvModels/interfaceTurbulenceDamping/interfaceTurbulenceDamping.H
@@ -131,8 +131,8 @@ class interfaceTurbulenceDamping
         void addRhoSup
         (
             const RhoType& rho,
-            fvMatrix<scalar>& eqn,
-            const word& fieldName
+            const volScalarField& field,
+            fvMatrix<scalar>& eqn
         ) const;
 
 
@@ -162,34 +162,38 @@ public:
 
     // Member Functions
 
-        //- Return the list of fields for which the option adds source term
-        //  to the transport equation
-        virtual wordList addSupFields() const;
+        // Checks
 
-        //- Add explicit contribution to mixture epsilon or omega equation
-        virtual void addSup
-        (
-            fvMatrix<scalar>& eqn,
-            const word& fieldName
-        ) const;
+            //- Return the list of fields for which the option adds source term
+            //  to the transport equation
+            virtual wordList addSupFields() 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
-        virtual void addSup
-        (
-            const volScalarField& alpha,
-            const volScalarField& rho,
-            fvMatrix<scalar>& eqn,
-            const word& fieldName
-        ) const;
+        // Sources
+
+            //- Add source to mixture epsilon or omega equation
+            virtual void addSup
+            (
+                const volScalarField& field,
+                fvMatrix<scalar>& eqn
+            ) 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
diff --git a/applications/modules/multiphaseEuler/fvModels/phaseTurbulenceStabilisation/phaseTurbulenceStabilisation.C b/applications/modules/multiphaseEuler/fvModels/phaseTurbulenceStabilisation/phaseTurbulenceStabilisation.C
index 390df3f23e..fae9c1054d 100644
--- a/applications/modules/multiphaseEuler/fvModels/phaseTurbulenceStabilisation/phaseTurbulenceStabilisation.C
+++ b/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 word& fieldName
+    const volScalarField& field,
+    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);
     }
 }
diff --git a/applications/modules/multiphaseEuler/fvModels/phaseTurbulenceStabilisation/phaseTurbulenceStabilisation.H b/applications/modules/multiphaseEuler/fvModels/phaseTurbulenceStabilisation/phaseTurbulenceStabilisation.H
index 9c2468ad4c..4ddafcc3cc 100644
--- a/applications/modules/multiphaseEuler/fvModels/phaseTurbulenceStabilisation/phaseTurbulenceStabilisation.H
+++ b/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
-        //  to the transport equation
-        virtual wordList addSupFields() const;
+        // Checks
 
-        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
-        (
-            const volScalarField& alpha,
-            const volScalarField& rho,
-            fvMatrix<scalar>& eqn,
-            const word& fieldName
-        ) const;
+
+        // Sources
+
+            //- Add contribution to phase k, epsilon or omega equation
+            virtual void addSup
+            (
+                const volScalarField& alpha,
+                const volScalarField& rho,
+                const volScalarField& field,
+                fvMatrix<scalar>& eqn
+            ) const;
 
 
         // Mesh changes
diff --git a/applications/modules/twoPhaseSolver/pressureCorrector.C b/applications/modules/twoPhaseSolver/pressureCorrector.C
index b2ac2d7de7..10d0c0895d 100644
--- a/applications/modules/twoPhaseSolver/pressureCorrector.C
+++ b/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
-        fvScalarMatrix Sp_rgh
+        // Cache any sources
+        fvScalarMatrix p_rghEqnSource
         (
-            fvModels().source
-            (
-                volScalarField::New
-                (
-                    "1",
-                    mesh,
-                    dimensionedScalar(dimless/dimPressure, 1)
-                ),
-                p_rgh
-            )
+            fvModels().sourceProxy(alpha1, p_rgh)
+          + fvModels().sourceProxy(alpha2, 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
diff --git a/etc/codeTemplates/dynamicCode/codedFvModelTemplate.C b/etc/codeTemplates/dynamicCode/codedFvModelTemplate.C
index 0f6137a355..cea03d2c56 100644
--- a/etc/codeTemplates/dynamicCode/codedFvModelTemplate.C
+++ b/etc/codeTemplates/dynamicCode/codedFvModelTemplate.C
@@ -127,8 +127,8 @@ ${typeName}FvModel${SourceType}::
 
 void ${typeName}FvModel${SourceType}::addSup
 (
-    fvMatrix<${TemplateType}>& eqn,
-    const word& fieldName
+    const VolField<${TemplateType}>& field,
+    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 word& fieldName
+    const VolField<${TemplateType}>& field,
+    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 word& fieldName
+    const VolField<${TemplateType}>& field,
+    fvMatrix<${TemplateType}>& eqn
 ) const
 {
     if (${verbose})
diff --git a/etc/codeTemplates/dynamicCode/codedFvModelTemplate.H b/etc/codeTemplates/dynamicCode/codedFvModelTemplate.H
index 07be5e798a..267c3a8acf 100644
--- a/etc/codeTemplates/dynamicCode/codedFvModelTemplate.H
+++ b/etc/codeTemplates/dynamicCode/codedFvModelTemplate.H
@@ -89,8 +89,8 @@ public:
             //- Explicit and implicit matrix contributions
             virtual void addSup
             (
-                fvMatrix<${TemplateType}>& eqn,
-                const word& fieldName
+                const VolField<${TemplateType}>& field,
+                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 word& fieldName
+                const VolField<${TemplateType}>& field,
+                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 word& fieldName
+                const VolField<${TemplateType}>& field,
+                fvMatrix<${TemplateType}>& eqn
             ) const;
 
 
diff --git a/src/finiteVolume/cfdTools/general/fvModels/fvModel.C b/src/finiteVolume/cfdTools/general/fvModels/fvModel.C
index cc6178999d..db8a1fe3a7 100644
--- a/src/finiteVolume/cfdTools/general/fvModels/fvModel.C
+++ b/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 word& fieldName
+    const VolField<Type>& field,
+    fvMatrix<Type>& eqn
 ) const
 {}
 
@@ -50,8 +55,8 @@ template<class Type>
 void Foam::fvModel::addSupType
 (
     const volScalarField& rho,
-    fvMatrix<Type>& eqn,
-    const word& fieldName
+    const VolField<Type>& field,
+    fvMatrix<Type>& eqn
 ) const
 {}
 
@@ -61,8 +66,8 @@ void Foam::fvModel::addSupType
 (
     const volScalarField& alpha,
     const volScalarField& rho,
-    fvMatrix<Type>& eqn,
-    const word& fieldName
+    const VolField<Type>& field,
+    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()
diff --git a/src/finiteVolume/cfdTools/general/fvModels/fvModel.H b/src/finiteVolume/cfdTools/general/fvModels/fvModel.H
index fbaa0bd027..4def6a81a5 100644
--- a/src/finiteVolume/cfdTools/general/fvModels/fvModel.H
+++ b/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 word& fieldName
+            const VolField<Type>& field,
+            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 word& fieldName
+            const VolField<Type>& field,
+            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 word& fieldName
+            const VolField<Type>& field,
+            fvMatrix<Type>& eqn
         ) 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 word& fieldName,
+            const VolField<Type>& eqnField,
             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
-        <
-            class Type,
-            class AlphaRhoFieldType,
-            class ... AlphaRhoFieldTypes
-        >
-        inline static dimensionSet sourceDims
+        template<class AlphaRhoFieldType, class ... AlphaRhoFieldTypes>
+        static dimensionSet sourceDims
         (
-            const VolField<Type>& field,
             const dimensionSet& ds,
-            const AlphaRhoFieldType& alphaRho,
-            const AlphaRhoFieldTypes& ... alphaRhos
+            const AlphaRhoFieldType& alphaRhoField,
+            const AlphaRhoFieldTypes& ... alphaRhoFields
         );
 
         //- Return the dimensions of the matrix of a source term (base
         //  condition for the above)
-        template<class Type>
-        inline static dimensionSet sourceDims
+        inline static const dimensionSet& sourceDims(const dimensionSet& ds);
+
+        //- 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 dimensionSet& ds
+            const AlphaRhoFieldType& alphaRhoField,
+            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
 
diff --git a/src/finiteVolume/cfdTools/general/fvModels/fvModelI.H b/src/finiteVolume/cfdTools/general/fvModels/fvModelI.H
index 782f4f23c3..664d311e20 100644
--- a/src/finiteVolume/cfdTools/general/fvModels/fvModelI.H
+++ b/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
diff --git a/src/finiteVolume/cfdTools/general/fvModels/fvModelM.H b/src/finiteVolume/cfdTools/general/fvModels/fvModelM.H
index 30f5ea7477..cbd45e4a8f 100644
--- a/src/finiteVolume/cfdTools/general/fvModels/fvModelM.H
+++ b/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                                                        \
-    (                                                                          \
-        fvMatrix<Type>& eqn,                                                   \
-        const word& fieldName                                                  \
-    ) const;
+    virtual void addSup(fvMatrix<Type>& eqn) const;
 
 #define IMPLEMENT_FV_MODEL_ADD_SUP(Type, modelType)                            \
-    void Foam::modelType::addSup                                               \
-    (                                                                          \
-        fvMatrix<Type>& eqn,                                                   \
-        const word& fieldName                                                  \
-    ) const                                                                    \
+    void Foam::modelType::addSup(fvMatrix<Type>& eqn) 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 word& fieldName                                                  \
+        const VolField<Type>& field,                                           \
+        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 word& fieldName                                                  \
+        const VolField<Type>& field,                                           \
+        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 word& fieldName                                                  \
+        const VolField<Type>& field,                                           \
+        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 word& fieldName                                                  \
+        const VolField<Type>& field,                                           \
+        fvMatrix<Type>& eqn                                                    \
     ) const                                                                    \
     {                                                                          \
-        addSupType(alpha, rho, eqn, fieldName);                                \
+        addSupType(alpha, rho, field, eqn);                                    \
     }
 
 
diff --git a/src/finiteVolume/cfdTools/general/fvModels/fvModelTemplates.C b/src/finiteVolume/cfdTools/general/fvModels/fvModelTemplates.C
index 92e6795dad..0a3987b7f6 100644
--- a/src/finiteVolume/cfdTools/general/fvModels/fvModelTemplates.C
+++ b/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 AlphaRhoFieldTypes& ... alphaRhos
+    const AlphaRhoFieldType& alphaRhoField,
+    const AlphaRhoFieldTypes& ... alphaRhoFields
 )
 {
-    return alphaRho.dimensions()*sourceDims(field, ds, alphaRhos ...);
+    return alphaRhoField.dimensions()*sourceDims(ds, alphaRhoFields ...);
 }
 
 
-template<class Type>
-Foam::dimensionSet Foam::fvModel::sourceDims
+template<class AlphaRhoFieldType, class ... AlphaRhoFieldTypes>
+const Foam::word& Foam::fvModel::fieldName
 (
-    const VolField<Type>& field,
-    const dimensionSet& ds
+    const AlphaRhoFieldType& alphaRhoField,
+    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 word& fieldName,
+    const VolField<Type>& eqnField,
     const dimensionSet& ds,
-    const AlphaRhoFieldTypes& ... alphaRhos
+    const AlphaRhoFieldTypes& ... alphaRhoFields
 ) const
 {
     tmp<fvMatrix<Type>> tmtx
     (
         new fvMatrix<Type>
         (
-            field,
-            sourceDims(field, ds, alphaRhos ...)
+            eqnField,
+            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 word& fieldName
+    const VolField<Type>& field
 ) 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));
-}
-
 // ************************************************************************* //
diff --git a/src/finiteVolume/cfdTools/general/fvModels/fvModels.H b/src/finiteVolume/cfdTools/general/fvModels/fvModels.H
index da13ca49ce..836d517c7d 100644
--- a/src/finiteVolume/cfdTools/general/fvModels/fvModels.H
+++ b/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 word& fieldName,
+            const VolField<Type>& eqnField,
             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
 
diff --git a/src/finiteVolume/cfdTools/general/fvModels/fvModelsTemplates.C b/src/finiteVolume/cfdTools/general/fvModels/fvModelsTemplates.C
index 933e003dd6..f00d577200 100644
--- a/src/finiteVolume/cfdTools/general/fvModels/fvModelsTemplates.C
+++ b/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 word& fieldName,
+    const VolField<Type>& eqnField,
     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,
-            fvModel::sourceDims(field, ds, alphaRhos ...)
+            eqnField,
+            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 word& fieldName
+    const VolField<Type>& field
 ) 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());
-}
-
-
-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));
+    return sourceTerm(field, dimVolume/sqr(dimTime), field);
 }
 
 
diff --git a/src/finiteVolume/fvMatrices/fvMatrix/fvMatrix.C b/src/finiteVolume/fvMatrices/fvMatrix/fvMatrix.C
index e8302523dd..869793da60 100644
--- a/src/finiteVolume/fvMatrices/fvMatrix/fvMatrix.C
+++ b/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)))
         )
     );
 
diff --git a/src/fvConstraints/zeroDimensionalFixedPressure/zeroDimensionalFixedPressureConstraint.C b/src/fvConstraints/zeroDimensionalFixedPressure/zeroDimensionalFixedPressureConstraint.C
index 29c6b3cbe1..968b8f1a23 100644
--- a/src/fvConstraints/zeroDimensionalFixedPressure/zeroDimensionalFixedPressureConstraint.C
+++ b/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
-        (
-            new volScalarField::Internal
+        return
+            volScalarField::Internal::New
             (
-                IOobject
-                (
-                    typedName("source"),
-                    mesh().time().timeName(),
-                    mesh(),
-                    IOobject::READ_IF_PRESENT,
-                    IOobject::AUTO_WRITE
-                ),
+                typedName("source"),
                 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
-    pEqnSource(pEqn);
-
-    // Check the dimensions have not changed
-    sourcePtr_->dimensions() = pEqn.dimensions()/dimVolume;
+    // Create the source field if it does not already exist
+    if (!sourcePtr_.valid())
+    {
+        sourcePtr_.set
+        (
+            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_();
diff --git a/src/fvConstraints/zeroDimensionalFixedPressure/zeroDimensionalFixedPressureConstraint.H b/src/fvConstraints/zeroDimensionalFixedPressure/zeroDimensionalFixedPressureConstraint.H
index 5f6906b3da..6220214383 100644
--- a/src/fvConstraints/zeroDimensionalFixedPressure/zeroDimensionalFixedPressureConstraint.H
+++ b/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;
 
diff --git a/src/fvConstraints/zeroDimensionalFixedPressure/zeroDimensionalFixedPressureModel.C b/src/fvConstraints/zeroDimensionalFixedPressure/zeroDimensionalFixedPressureModel.C
index 47c94429ce..0f40273f4d 100644
--- a/src/fvConstraints/zeroDimensionalFixedPressure/zeroDimensionalFixedPressureModel.C
+++ b/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 word& fieldName
+    const VolField<Type>& field,
+    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 word& fieldName
+    const volScalarField& rho,
+    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 word& fieldName
+    const VolField<Type>& field,
+    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,
-    const word& fieldName
+    fvMatrix<scalar>& eqn
 ) 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);
-        }
-        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());
+            eqn += alpha()*constraint().pEqnSource(rho, eqn);
         }
         else
         {
-            FatalErrorInFunction
-                << "Cannot add source for density field " << fieldName
-                << " into an equation for " << eqn.psi().name()
-                << exit(FatalError);
+            eqn += constraint().massSource(alpha(), rho());
         }
     }
     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 word& fieldName
+    const VolField<Type>& field,
+    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()
diff --git a/src/fvConstraints/zeroDimensionalFixedPressure/zeroDimensionalFixedPressureModel.H b/src/fvConstraints/zeroDimensionalFixedPressure/zeroDimensionalFixedPressureModel.H
index c4ca2c7b24..c17da126d5 100644
--- a/src/fvConstraints/zeroDimensionalFixedPressure/zeroDimensionalFixedPressureModel.H
+++ b/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
-            void addSupType(fvMatrix<scalar>& eqn, const word& fieldName) const;
+            //- Add a source term to a compressible continuity equation
+            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 word& fieldName
+                const VolField<Type>& field,
+                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,
-                const word& fieldName
+                fvMatrix<scalar>& eqn
             ) 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 word& fieldName
-            ) 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 VolField<Type>& field,
+                fvMatrix<Type>& eqn
             ) 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
diff --git a/src/fvModels/Make/files b/src/fvModels/Make/files
index c26d420284..daf3fca336 100644
--- a/src/fvModels/Make/files
+++ b/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
diff --git a/src/fvModels/derived/accelerationSource/accelerationSource.C b/src/fvModels/derived/accelerationSource/accelerationSource.C
index a57acbf4f3..6b14f020ab 100644
--- a/src/fvModels/derived/accelerationSource/accelerationSource.C
+++ b/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 word& fieldName
+    const volVectorField& U,
+    fvMatrix<vector>& eqn
 ) const
 {
-    add(geometricOneField(), eqn, fieldName);
+    add(geometricOneField(), eqn);
 }
 
 
 void Foam::fv::accelerationSource::addSup
 (
     const volScalarField& rho,
-    fvMatrix<vector>& eqn,
-    const word& fieldName
+    const volVectorField& U,
+    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 word& fieldName
+    const volVectorField& U,
+    fvMatrix<vector>& eqn
 ) const
 {
-    add((alpha*rho)(), eqn, fieldName);
+    add((alpha*rho)(), eqn);
 }
 
 
diff --git a/src/fvModels/derived/accelerationSource/accelerationSource.H b/src/fvModels/derived/accelerationSource/accelerationSource.H
index 335aad9f98..e1ad4c73d3 100644
--- a/src/fvModels/derived/accelerationSource/accelerationSource.H
+++ b/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,
-            const word& fieldName
+            fvMatrix<vector>& eqn
         ) const;
 
 
@@ -141,16 +140,16 @@ public:
             //- Source term to momentum equation
             virtual void addSup
             (
-                fvMatrix<vector>& eqn,
-                const word& fieldName
+                const volVectorField& U,
+                fvMatrix<vector>& eqn
             ) const;
 
             //- Source term to compressible momentum equation
             virtual void addSup
             (
                 const volScalarField& rho,
-                fvMatrix<vector>& eqn,
-                const word& fieldName
+                const volVectorField& U,
+                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 word& fieldName
+                const volVectorField& U,
+                fvMatrix<vector>& eqn
             ) const;
 
 
@@ -192,12 +191,6 @@ public:
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-#ifdef NoRepository
-    #include "accelerationSourceTemplates.C"
-#endif
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
 #endif
 
 // ************************************************************************* //
diff --git a/src/fvModels/derived/accelerationSource/accelerationSourceTemplates.C b/src/fvModels/derived/accelerationSource/accelerationSourceTemplates.C
deleted file mode 100644
index edd5d65096..0000000000
--- a/src/fvModels/derived/accelerationSource/accelerationSourceTemplates.C
+++ /dev/null
@@ -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;
-    }
-}
-
-
-// ************************************************************************* //
diff --git a/src/fvModels/derived/actuationDiskSource/actuationDiskSource.C b/src/fvModels/derived/actuationDiskSource/actuationDiskSource.C
index 01fe18687f..4884f1a979 100644
--- a/src/fvModels/derived/actuationDiskSource/actuationDiskSource.C
+++ b/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 word& fieldName
+    const volVectorField& U,
+    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 word& fieldName
+    const volVectorField& U,
+    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 word& fieldName
+    const volVectorField& U,
+    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
diff --git a/src/fvModels/derived/actuationDiskSource/actuationDiskSource.H b/src/fvModels/derived/actuationDiskSource/actuationDiskSource.H
index 9749a109c1..1039be4736 100644
--- a/src/fvModels/derived/actuationDiskSource/actuationDiskSource.H
+++ b/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 word& fieldName
+                const volVectorField& U,
+                fvMatrix<vector>& eqn
             ) const;
 
             //- Source term to compressible momentum equation
             virtual void addSup
             (
                 const volScalarField& rho,
-                fvMatrix<vector>& eqn,
-                const word& fieldName
+                const volVectorField& U,
+                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 word& fieldName
+                const volVectorField& U,
+                fvMatrix<vector>& eqn
             ) const;
 
 
@@ -237,12 +237,6 @@ public:
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-#ifdef NoRepository
-    #include "actuationDiskSourceTemplates.C"
-#endif
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
 #endif
 
 // ************************************************************************* //
diff --git a/src/fvModels/derived/actuationDiskSource/actuationDiskSourceTemplates.C b/src/fvModels/derived/actuationDiskSource/actuationDiskSourceTemplates.C
deleted file mode 100644
index a14d98549f..0000000000
--- a/src/fvModels/derived/actuationDiskSource/actuationDiskSourceTemplates.C
+++ /dev/null
@@ -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;
-    }
-}
-
-
-// ************************************************************************* //
diff --git a/src/fvModels/derived/buoyancyEnergy/buoyancyEnergy.C b/src/fvModels/derived/buoyancyEnergy/buoyancyEnergy.C
index b5d9126fd1..ed2b5d2f35 100644
--- a/src/fvModels/derived/buoyancyEnergy/buoyancyEnergy.C
+++ b/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 word& fieldName
+    const volScalarField& he,
+    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 word& fieldName
+    const volScalarField& he,
+    fvMatrix<scalar>& eqn
 ) const
 {
     const uniformDimensionedVectorField& g =
diff --git a/src/fvModels/derived/buoyancyEnergy/buoyancyEnergy.H b/src/fvModels/derived/buoyancyEnergy/buoyancyEnergy.H
index 4cf9c86bf3..75a4169bc6 100644
--- a/src/fvModels/derived/buoyancyEnergy/buoyancyEnergy.H
+++ b/src/fvModels/derived/buoyancyEnergy/buoyancyEnergy.H
@@ -116,8 +116,8 @@ public:
             virtual void addSup
             (
                 const volScalarField& rho,
-                fvMatrix<scalar>& eqn,
-                const word& fieldName
+                const volScalarField& he,
+                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 word& fieldName
+                const volScalarField& he,
+                fvMatrix<scalar>& eqn
             ) const;
 
 
diff --git a/src/fvModels/derived/buoyancyForce/buoyancyForce.C b/src/fvModels/derived/buoyancyForce/buoyancyForce.C
index 38052d2634..ddcba19ba3 100644
--- a/src/fvModels/derived/buoyancyForce/buoyancyForce.C
+++ b/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 word& fieldName
+    const volVectorField& U,
+    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 word& fieldName
+    const volVectorField& U,
+    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 word& fieldName
+    const volVectorField& U,
+    fvMatrix<vector>& eqn
 ) const
 {
     eqn += alpha*rho*g_;
diff --git a/src/fvModels/derived/buoyancyForce/buoyancyForce.H b/src/fvModels/derived/buoyancyForce/buoyancyForce.H
index c8f730a673..8fd260dfb9 100644
--- a/src/fvModels/derived/buoyancyForce/buoyancyForce.H
+++ b/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 word& fieldName
+                const volVectorField& U,
+                fvMatrix<vector>& eqn
             ) const;
 
             //- Add explicit contribution to compressible momentum equation
             virtual void addSup
             (
                 const volScalarField& rho,
-                fvMatrix<vector>& eqn,
-                const word& fieldName
+                const volVectorField& U,
+                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 word& fieldName
+                const volVectorField& U,
+                fvMatrix<vector>& eqn
             ) const;
 
 
diff --git a/src/fvModels/derived/effectivenessHeatExchangerSource/effectivenessHeatExchangerSource.C b/src/fvModels/derived/effectivenessHeatExchangerSource/effectivenessHeatExchangerSource.C
index 285a8856d9..5d4a65cfb3 100644
--- a/src/fvModels/derived/effectivenessHeatExchangerSource/effectivenessHeatExchangerSource.C
+++ b/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 word& fieldName
+    const volScalarField& he,
+    fvMatrix<scalar>& eqn
 ) const
 {
     const basicThermo& thermo =
diff --git a/src/fvModels/derived/effectivenessHeatExchangerSource/effectivenessHeatExchangerSource.H b/src/fvModels/derived/effectivenessHeatExchangerSource/effectivenessHeatExchangerSource.H
index bb9c07dea6..1323b2a615 100644
--- a/src/fvModels/derived/effectivenessHeatExchangerSource/effectivenessHeatExchangerSource.H
+++ b/src/fvModels/derived/effectivenessHeatExchangerSource/effectivenessHeatExchangerSource.H
@@ -213,8 +213,8 @@ public:
         virtual void addSup
         (
             const volScalarField& rho,
-            fvMatrix<scalar>& eqn,
-            const word& fieldName
+            const volScalarField& he,
+            fvMatrix<scalar>& eqn
         ) const;
 
 
diff --git a/src/fvModels/derived/explicitPorositySource/explicitPorositySource.C b/src/fvModels/derived/explicitPorositySource/explicitPorositySource.C
index 5b7a2c818d..97f680245f 100644
--- a/src/fvModels/derived/explicitPorositySource/explicitPorositySource.C
+++ b/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 word& fieldName
+    const volVectorField& U,
+    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 word& fieldName
+    const volVectorField& U,
+    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 word& fieldName
+    const volVectorField& U,
+    fvMatrix<vector>& eqn
 ) const
 {
     fvMatrix<vector> porosityEqn(eqn.psi(), eqn.dimensions());
diff --git a/src/fvModels/derived/explicitPorositySource/explicitPorositySource.H b/src/fvModels/derived/explicitPorositySource/explicitPorositySource.H
index d2f53a64f6..9db51eeae5 100644
--- a/src/fvModels/derived/explicitPorositySource/explicitPorositySource.H
+++ b/src/fvModels/derived/explicitPorositySource/explicitPorositySource.H
@@ -154,16 +154,16 @@ public:
             //- Add implicit contribution to momentum equation
             virtual void addSup
             (
-                fvMatrix<vector>& eqn,
-                const word& fieldName
+                const volVectorField& U,
+                fvMatrix<vector>& eqn
             ) const;
 
             //- Add implicit contribution to compressible momentum equation
             virtual void addSup
             (
                 const volScalarField& rho,
-                fvMatrix<vector>& eqn,
-                const word& fieldName
+                const volVectorField& U,
+                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 word& fieldName
+                const volVectorField& U,
+                fvMatrix<vector>& eqn
             ) const;
 
 
diff --git a/src/fvModels/derived/heatSource/heatSource.C b/src/fvModels/derived/heatSource/heatSource.C
index b00c306463..f9cc69acec 100644
--- a/src/fvModels/derived/heatSource/heatSource.C
+++ b/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 word& fieldName
+    const volScalarField& he,
+    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 word& fieldName
+    const volScalarField& he,
+    fvMatrix<scalar>& eqn
 ) const
 {
-    addSup(eqn, fieldName);
+    addSup(he, eqn);
 }
 
 
diff --git a/src/fvModels/derived/heatSource/heatSource.H b/src/fvModels/derived/heatSource/heatSource.H
index d78a8a9f90..8e33e74d6f 100644
--- a/src/fvModels/derived/heatSource/heatSource.H
+++ b/src/fvModels/derived/heatSource/heatSource.H
@@ -117,16 +117,16 @@ public:
             //- Source term to energy equation
             virtual void addSup
             (
-                fvMatrix<scalar>& eqn,
-                const word& fieldName
+                const volScalarField& he,
+                fvMatrix<scalar>& eqn
             ) const;
 
             //- Source term to compressible energy equation
             virtual void addSup
             (
                 const volScalarField& rho,
-                fvMatrix<scalar>& eqn,
-                const word& fieldName
+                const volScalarField& he,
+                fvMatrix<scalar>& eqn
             ) const;
 
 
diff --git a/src/fvModels/derived/heatTransfer/heatTransfer.C b/src/fvModels/derived/heatTransfer/heatTransfer.C
index 282c703bcf..8958b2bbac 100644
--- a/src/fvModels/derived/heatTransfer/heatTransfer.C
+++ b/src/fvModels/derived/heatTransfer/heatTransfer.C
@@ -71,8 +71,7 @@ template<class AlphaFieldType>
 void Foam::fv::heatTransfer::add
 (
     const AlphaFieldType& alpha,
-    fvMatrix<scalar>& eqn,
-    const word& fieldName
+    fvMatrix<scalar>& eqn
 ) 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 word& fieldName
+    const volScalarField& he,
+    fvMatrix<scalar>& eqn
 ) const
 {
-    add(geometricOneField(), eqn, fieldName);
+    add(geometricOneField(), eqn);
 }
 
 
 void Foam::fv::heatTransfer::addSup
 (
     const volScalarField& rho,
-    fvMatrix<scalar>& eqn,
-    const word& fieldName
+    const volScalarField& he,
+    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 word& fieldName
+    const volScalarField& he,
+    fvMatrix<scalar>& eqn
 ) const
 {
-    add(alpha, eqn, fieldName);
+    add(alpha, eqn);
 }
 
 
diff --git a/src/fvModels/derived/heatTransfer/heatTransfer.H b/src/fvModels/derived/heatTransfer/heatTransfer.H
index c78272a356..ab4537f8a0 100644
--- a/src/fvModels/derived/heatTransfer/heatTransfer.H
+++ b/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,
-            const word& fieldName
+            fvMatrix<scalar>& eqn
         ) const;
 
 
@@ -168,16 +167,16 @@ public:
             //- Source term to energy equation
             virtual void addSup
             (
-                fvMatrix<scalar>& eqn,
-                const word& fieldName
+                const volScalarField& he,
+                fvMatrix<scalar>& eqn
             ) const;
 
             //- Source term to compressible energy equation
             virtual void addSup
             (
                 const volScalarField& rho,
-                fvMatrix<scalar>& eqn,
-                const word& fieldName
+                const volScalarField& he,
+                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 word& fieldName
+                const volScalarField& he,
+                fvMatrix<scalar>& eqn
             ) const;
 
 
diff --git a/src/fvModels/derived/massSource/massSource.C b/src/fvModels/derived/massSource/massSource.C
index 3d2f9b9a0c..e879b1a2bd 100644
--- a/src/fvModels/derived/massSource/massSource.C
+++ b/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 word& fieldName
+    const VolField<Type>& field,
+    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 word& fieldName
+    const volScalarField& rho,
+    const volScalarField& field,
+    fvMatrix<scalar>& eqn
 ) const
 {
-    addGeneralSupType(eqn, fieldName);
-}
-
-
-void Foam::fv::massSourceBase::addSupType
-(
-    fvMatrix<scalar>& eqn,
-    const word& fieldName
-) const
-{
-    const labelUList cells = set_.cells();
-
-    if (fieldName == rhoName_)
+    // Multiphase continuity equation. Same source as single-phase case.
+    if (field.name() == rhoName_)
     {
-        const scalar massFlowRate = this->massFlowRate();
-
-        forAll(cells, i)
-        {
-            eqn.source()[cells[i]] -=
-                mesh().V()[cells[i]]/set_.V()*massFlowRate;
-        }
+        addSupType(field, eqn);
+        return;
     }
-    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>
-void Foam::fv::massSourceBase::addSupType
-(
-    const volScalarField& rho,
-    fvMatrix<Type>& eqn,
-    const word& fieldName
-) const
-{
-    addSupType(eqn, fieldName);
+    // Property equation
+    addSupType<scalar>(rho, field, eqn);
 }
 
 
@@ -209,11 +234,12 @@ void Foam::fv::massSourceBase::addSupType
 (
     const volScalarField& alpha,
     const volScalarField& rho,
-    fvMatrix<Type>& eqn,
-    const word& fieldName
+    const VolField<Type>& field,
+    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()
diff --git a/src/fvModels/derived/massSource/massSource.H b/src/fvModels/derived/massSource/massSource.H
index 65cc1fb37d..ba2da30662 100644
--- a/src/fvModels/derived/massSource/massSource.H
+++ b/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 word& fieldName
+                const VolField<Type>& field,
+                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 word& fieldName
+                const VolField<Type>& field,
+                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 word& fieldName
+                const VolField<Type>& field,
+                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
diff --git a/src/fvModels/derived/phaseLimitStabilisation/phaseLimitStabilisation.C b/src/fvModels/derived/phaseLimitStabilisation/phaseLimitStabilisation.C
index 111c443c8c..afcd191acc 100644
--- a/src/fvModels/derived/phaseLimitStabilisation/phaseLimitStabilisation.C
+++ b/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 word& fieldName
+    const VolField<Type>& field,
+    fvMatrix<Type>& eqn
 ) const
 {
-    const VolField<Type>& psi = eqn.psi();
-
-    uniformDimensionedScalarField& rate =
+    const 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()
diff --git a/src/fvModels/derived/phaseLimitStabilisation/phaseLimitStabilisation.H b/src/fvModels/derived/phaseLimitStabilisation/phaseLimitStabilisation.H
index c06ed563f8..ba2abe5b78 100644
--- a/src/fvModels/derived/phaseLimitStabilisation/phaseLimitStabilisation.H
+++ b/src/fvModels/derived/phaseLimitStabilisation/phaseLimitStabilisation.H
@@ -98,8 +98,8 @@ class phaseLimitStabilisation
         (
             const volScalarField& alpha,
             const volScalarField& rho,
-            fvMatrix<Type>& eqn,
-            const word& fieldName
+            const VolField<Type>& field,
+            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
diff --git a/src/fvModels/derived/radialActuationDiskSource/radialActuationDiskSource.C b/src/fvModels/derived/radialActuationDiskSource/radialActuationDiskSource.C
index 52caa1c2d6..7c4a64919f 100644
--- a/src/fvModels/derived/radialActuationDiskSource/radialActuationDiskSource.C
+++ b/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 word& fieldName
+    const volVectorField& U,
+    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 word& fieldName
+    const volVectorField& U,
+    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
     );
diff --git a/src/fvModels/derived/radialActuationDiskSource/radialActuationDiskSource.H b/src/fvModels/derived/radialActuationDiskSource/radialActuationDiskSource.H
index 9e90065a58..6812f05220 100644
--- a/src/fvModels/derived/radialActuationDiskSource/radialActuationDiskSource.H
+++ b/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 word& fieldName
+            const volVectorField& U,
+            fvMatrix<vector>& eqn
         ) const;
 
         //- Source term to compressible momentum equation
         virtual void addSup
         (
             const volScalarField& rho,
-            fvMatrix<vector>& eqn,
-            const word& fieldName
+            const volVectorField& U,
+            fvMatrix<vector>& eqn
         ) const;
 
 
@@ -189,12 +189,6 @@ public:
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-#ifdef NoRepository
-    #include "radialActuationDiskSourceTemplates.C"
-#endif
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
 #endif
 
 // ************************************************************************* //
diff --git a/src/fvModels/derived/rotorDiskSource/rotorDiskSource.C b/src/fvModels/derived/rotorDiskSource/rotorDiskSource.C
index 11167c6d93..969ccf8bed 100644
--- a/src/fvModels/derived/rotorDiskSource/rotorDiskSource.C
+++ b/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 word& fieldName
+    const volVectorField& U,
+    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 word& fieldName
+    const volVectorField& U,
+    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);
 
diff --git a/src/fvModels/derived/rotorDiskSource/rotorDiskSource.H b/src/fvModels/derived/rotorDiskSource/rotorDiskSource.H
index 29422c89fa..cd9a36ebe5 100644
--- a/src/fvModels/derived/rotorDiskSource/rotorDiskSource.H
+++ b/src/fvModels/derived/rotorDiskSource/rotorDiskSource.H
@@ -318,16 +318,16 @@ public:
             //- Source term to momentum equation
             virtual void addSup
             (
-                fvMatrix<vector>& eqn,
-                const word& fieldName
+                const volVectorField& U,
+                fvMatrix<vector>& eqn
             ) const;
 
             //- Source term to compressible momentum equation
             virtual void addSup
             (
                 const volScalarField& rho,
-                fvMatrix<vector>& eqn,
-                const word& fieldName
+                const volVectorField& U,
+                fvMatrix<vector>& eqn
             ) const;
 
 
diff --git a/src/fvModels/derived/sixDoFAccelerationSource/sixDoFAccelerationSource.C b/src/fvModels/derived/sixDoFAccelerationSource/sixDoFAccelerationSource.C
index e93926bffd..d5c3da2afa 100644
--- a/src/fvModels/derived/sixDoFAccelerationSource/sixDoFAccelerationSource.C
+++ b/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 word& fieldName
+    const volVectorField& U,
+    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 word& fieldName
+    const volVectorField& U,
+    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 word& fieldName
+    const volVectorField& U,
+    fvMatrix<vector>& eqn
 ) const
 {
-    addForce(alpha, rho, eqn, fieldName);
+    addForce(alpha, rho, U, eqn);
 }
 
 
diff --git a/src/fvModels/derived/sixDoFAccelerationSource/sixDoFAccelerationSource.H b/src/fvModels/derived/sixDoFAccelerationSource/sixDoFAccelerationSource.H
index 3be269a3db..1234b57a26 100644
--- a/src/fvModels/derived/sixDoFAccelerationSource/sixDoFAccelerationSource.H
+++ b/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 word& fieldName
+            const volVectorField& U,
+            fvMatrix<vector>& eqn
         ) const;
 
 
@@ -145,16 +145,16 @@ public:
         //- Source term to momentum equation
         virtual void addSup
         (
-            fvMatrix<vector>& eqn,
-            const word& fieldName
+            const volVectorField& U,
+            fvMatrix<vector>& eqn
         ) const;
 
         //- Source term to compressible momentum equation
         virtual void addSup
         (
             const volScalarField& rho,
-            fvMatrix<vector>& eqn,
-            const word& fieldName
+            const volVectorField& U,
+            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 word& fieldName
+            const volVectorField& U,
+            fvMatrix<vector>& eqn
         ) const;
 
 
@@ -200,12 +200,6 @@ public:
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-#ifdef NoRepository
-    #include "sixDoFAccelerationSourceTemplates.C"
-#endif
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
 #endif
 
 // ************************************************************************* //
diff --git a/src/fvModels/derived/sixDoFAccelerationSource/sixDoFAccelerationSourceTemplates.C b/src/fvModels/derived/sixDoFAccelerationSource/sixDoFAccelerationSourceTemplates.C
deleted file mode 100644
index 52181f82fa..0000000000
--- a/src/fvModels/derived/sixDoFAccelerationSource/sixDoFAccelerationSourceTemplates.C
+++ /dev/null
@@ -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
-    );
-}
-
-
-// ************************************************************************* //
diff --git a/src/fvModels/derived/solidEquilibriumEnergySource/solidEquilibriumEnergySource.C b/src/fvModels/derived/solidEquilibriumEnergySource/solidEquilibriumEnergySource.C
index 21530c00ab..5e9a9a4f3b 100644
--- a/src/fvModels/derived/solidEquilibriumEnergySource/solidEquilibriumEnergySource.C
+++ b/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 word& fieldName
+    const volScalarField& he,
+    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 word& fieldName
+    const volScalarField& he,
+    fvMatrix<scalar>& eqn
 ) const
 {
     const volScalarField alphahe
diff --git a/src/fvModels/derived/solidEquilibriumEnergySource/solidEquilibriumEnergySource.H b/src/fvModels/derived/solidEquilibriumEnergySource/solidEquilibriumEnergySource.H
index 244d369c88..c8cd7c08bf 100644
--- a/src/fvModels/derived/solidEquilibriumEnergySource/solidEquilibriumEnergySource.H
+++ b/src/fvModels/derived/solidEquilibriumEnergySource/solidEquilibriumEnergySource.H
@@ -138,18 +138,18 @@ public:
             //- Explicit and implicit sources for compressible equations
             virtual void addSup
             (
-                const volScalarField&,
-                fvMatrix<scalar>&,
-                const word& fieldName
+                const volScalarField& rho,
+                const volScalarField& he,
+                fvMatrix<scalar>& eqn
             ) const;
 
             //- Explicit and implicit sources for phase equations
             virtual void addSup
             (
-                const volScalarField&,
-                const volScalarField&,
-                fvMatrix<scalar>&,
-                const word& fieldName
+                const volScalarField& alpha,
+                const volScalarField& rho,
+                const volScalarField& he,
+                fvMatrix<scalar>& eqn
             ) const;
 
 
diff --git a/src/fvModels/derived/solidificationMeltingSource/solidificationMeltingSource.C b/src/fvModels/derived/solidificationMeltingSource/solidificationMeltingSource.C
index 9605312395..7bf160696e 100644
--- a/src/fvModels/derived/solidificationMeltingSource/solidificationMeltingSource.C
+++ b/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 word& fieldName
+    const volScalarField& he,
+    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 word& fieldName
+    const volScalarField& he,
+    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 word& fieldName
+    const volVectorField& U,
+    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 word& fieldName
+    const volVectorField& U,
+    fvMatrix<vector>& eqn
 ) const
 {
-    // Momentum source uses a Boussinesq approximation - redirect
-    addSup(eqn, fieldName);
+    addSup(U, eqn);
 }
 
 
diff --git a/src/fvModels/derived/solidificationMeltingSource/solidificationMeltingSource.H b/src/fvModels/derived/solidificationMeltingSource/solidificationMeltingSource.H
index 8b7b320269..94767def6e 100644
--- a/src/fvModels/derived/solidificationMeltingSource/solidificationMeltingSource.H
+++ b/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 word& fieldName
+                const volScalarField& he,
+                fvMatrix<scalar>& eqn
             ) const;
 
             //- Add implicit contribution to momentum equation
             virtual void addSup
             (
-                fvMatrix<vector>& eqn,
-                const word& fieldName
+                const volVectorField& U,
+                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 word& fieldName
+                const volScalarField& he,
+                fvMatrix<scalar>& eqn
             ) const;
 
             //- Add implicit contribution to compressible momentum equation
             virtual void addSup
             (
                 const volScalarField& rho,
-                fvMatrix<vector>& eqn,
-                const word& fieldName
+                const volVectorField& U,
+                fvMatrix<vector>& eqn
             ) const;
 
 
@@ -337,12 +334,6 @@ public:
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-#ifdef NoRepository
-    #include "solidificationMeltingSourceTemplates.C"
-#endif
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
 #endif
 
 // ************************************************************************* //
diff --git a/src/fvModels/derived/solidificationMeltingSource/solidificationMeltingSourceTemplates.C b/src/fvModels/derived/solidificationMeltingSource/solidificationMeltingSourceTemplates.C
deleted file mode 100644
index 7781614d01..0000000000
--- a/src/fvModels/derived/solidificationMeltingSource/solidificationMeltingSourceTemplates.C
+++ /dev/null
@@ -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_));
-    }
-}
-
-
-// ************************************************************************* //
diff --git a/src/fvModels/derived/volumeFractionSource/volumeFractionSource.C b/src/fvModels/derived/volumeFractionSource/volumeFractionSource.C
index 8969ce5c16..5b2693f800 100644
--- a/src/fvModels/derived/volumeFractionSource/volumeFractionSource.C
+++ b/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,
-    const word& fieldName
+    fvMatrix<Type>& eqn
 ) 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 word& fieldName
+    const VolField<Type>& field,
+    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 word& fieldName
+    const volScalarField& field,
+    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 word& fieldName
+    const volVectorField& field,
+    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 word& fieldName
+    const VolField<Type>& field,
+    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 word& fieldName
+    const VolField<Type>& field,
+    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 word& fieldName
+    const VolField<Type>& field,
+    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_RHO_SUP, fv::volumeFractionSource);
+FOR_ALL_FIELD_TYPES
+(
+    IMPLEMENT_FV_MODEL_ADD_FIELD_SUP,
+    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()
diff --git a/src/fvModels/derived/volumeFractionSource/volumeFractionSource.H b/src/fvModels/derived/volumeFractionSource/volumeFractionSource.H
index e3149f3c9c..37d14b06ca 100644
--- a/src/fvModels/derived/volumeFractionSource/volumeFractionSource.H
+++ b/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>&,
-                const word& fieldName
+                fvMatrix<Type>&
             ) const;
 
             //- Add a source term to an equation
@@ -132,8 +131,8 @@ class volumeFractionSource
             void addAlphaSupType
             (
                 const AlphaFieldType& alpha,
-                fvMatrix<Type>& eqn,
-                const word& fieldName
+                const VolField<Type>& field,
+                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 word& fieldName
+                const volScalarField& field,
+                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 word& fieldName
+                const volVectorField& field,
+                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 word& fieldName
+                const VolField<Type>& field,
+                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 word& fieldName
+                const VolField<Type>& field,
+                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
diff --git a/src/fvModels/general/codedFvModel/codedFvModel.C b/src/fvModels/general/codedFvModel/codedFvModel.C
index e84bd06774..e2998a2bfd 100644
--- a/src/fvModels/general/codedFvModel/codedFvModel.C
+++ b/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 word& fieldName
+    const VolField<Type>& field,
+    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 word& fieldName
+    const VolField<Type>& field,
+    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 word& fieldName
+    const VolField<Type>& field,
+    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()
diff --git a/src/fvModels/general/codedFvModel/codedFvModel.H b/src/fvModels/general/codedFvModel/codedFvModel.H
index 571f4a12e0..5938a1f754 100644
--- a/src/fvModels/general/codedFvModel/codedFvModel.H
+++ b/src/fvModels/general/codedFvModel/codedFvModel.H
@@ -136,8 +136,8 @@ class codedFvModel
             template<class Type>
             void addSupType
             (
-                fvMatrix<Type>& eqn,
-                const word& fieldName
+                const VolField<Type>& field,
+                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 word& fieldName
+                const VolField<Type>& field,
+                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 word& fieldName
+                const VolField<Type>& field,
+                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
diff --git a/src/fvModels/general/semiImplicitSource/semiImplicitSource.C b/src/fvModels/general/semiImplicitSource/semiImplicitSource.C
index 622094c74b..e76e1c1c51 100644
--- a/src/fvModels/general/semiImplicitSource/semiImplicitSource.C
+++ b/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 word& fieldName
+    const VolField<Type>& field,
+    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 word& fieldName
+    const VolField<Type>& field,
+    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 word& fieldName
+    const VolField<Type>& field,
+    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_RHO_SUP, fv::semiImplicitSource);
+FOR_ALL_FIELD_TYPES
+(
+    IMPLEMENT_FV_MODEL_ADD_FIELD_SUP,
+    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()
diff --git a/src/fvModels/general/semiImplicitSource/semiImplicitSource.H b/src/fvModels/general/semiImplicitSource/semiImplicitSource.H
index 79aedbb726..c014bbb0bf 100644
--- a/src/fvModels/general/semiImplicitSource/semiImplicitSource.H
+++ b/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 word& fieldName
+                const VolField<Type>& field,
+                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 word& fieldName
+                const VolField<Type>& field,
+                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
diff --git a/src/fvModels/interRegion/interRegionExplicitPorositySource/interRegionExplicitPorositySource.C b/src/fvModels/interRegion/interRegionExplicitPorositySource/interRegionExplicitPorositySource.C
index 34d402747b..634af86e2d 100644
--- a/src/fvModels/interRegion/interRegionExplicitPorositySource/interRegionExplicitPorositySource.C
+++ b/src/fvModels/interRegion/interRegionExplicitPorositySource/interRegionExplicitPorositySource.C
@@ -152,9 +152,8 @@ Foam::fv::interRegionExplicitPorositySource::addSupFields() const
 
 void Foam::fv::interRegionExplicitPorositySource::addSup
 (
-    const volScalarField& rho,
-    fvMatrix<vector>& eqn,
-    const word& fieldName
+    const volVectorField& U,
+    fvMatrix<vector>& eqn
 ) 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;
diff --git a/src/fvModels/interRegion/interRegionExplicitPorositySource/interRegionExplicitPorositySource.H b/src/fvModels/interRegion/interRegionExplicitPorositySource/interRegionExplicitPorositySource.H
index 1e8af43fdf..195acc729d 100644
--- a/src/fvModels/interRegion/interRegionExplicitPorositySource/interRegionExplicitPorositySource.H
+++ b/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 word& fieldName
+                const volVectorField& U,
+                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;
 
 
diff --git a/src/fvModels/interRegion/interRegionHeatTransfer/interRegionHeatTransfer.C b/src/fvModels/interRegion/interRegionHeatTransfer/interRegionHeatTransfer.C
index 300599c499..de7c1359d0 100644
--- a/src/fvModels/interRegion/interRegionHeatTransfer/interRegionHeatTransfer.C
+++ b/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 word& fieldName
+    const volScalarField& he,
+    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 word& fieldName
+    const volScalarField& he,
+    fvMatrix<scalar>& eqn
 ) const
 {
-    addSup(eqn, fieldName);
+    addSup(he, eqn);
 }
 
 
diff --git a/src/fvModels/interRegion/interRegionHeatTransfer/interRegionHeatTransfer.H b/src/fvModels/interRegion/interRegionHeatTransfer/interRegionHeatTransfer.H
index cf6409655b..29fd1584a1 100644
--- a/src/fvModels/interRegion/interRegionHeatTransfer/interRegionHeatTransfer.H
+++ b/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 word& fieldName
+                const volScalarField& he,
+                fvMatrix<scalar>& eqn
             ) const;
 
             //- Source term to compressible energy equation
             virtual void addSup
             (
                 const volScalarField& rho,
-                fvMatrix<scalar>& eqn,
-                const word& fieldName
+                const volScalarField& he,
+                fvMatrix<scalar>& eqn
             ) const;
 
 
diff --git a/src/lagrangian/parcel/fvModels/clouds/clouds.C b/src/lagrangian/parcel/fvModels/clouds/clouds.C
index 7a59fddfca..5a90284821 100644
--- a/src/lagrangian/parcel/fvModels/clouds/clouds.C
+++ b/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(Y, i)
+            forAll(carrierMcThermo.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 word& fieldName
+    const volScalarField& rho,
+    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 word& fieldName
+    const volScalarField& field,
+    fvMatrix<scalar>& eqn
 ) const
 {
     if (debug)
@@ -254,36 +258,41 @@ void Foam::fv::clouds::addSup
 
     const fluidThermo& carrierThermo = tCarrierThermo_();
 
-    if (fieldName == rhoName_)
-    {
-        eqn += cloudsPtr_().Srho(eqn.psi());
-    }
-    else if (fieldName == carrierThermo.he().name())
+    if (&field == &carrierThermo.he())
     {
         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()),
-                    eqn.psi()
+                    carrierMcThermo.specieIndex(field),
+                    field
                 );
+            return;
         }
     }
+
+    {
+        FatalErrorInFunction
+            << "Support for field " << field.name() << " is not implemented"
+            << exit(FatalError);
+    }
 }
 
 
 void Foam::fv::clouds::addSup
 (
-    fvMatrix<vector>& eqn,
-    const word& fieldName
+    const volVectorField& U,
+    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 word& fieldName
+    const volVectorField& U,
+    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);
+    }
 }
 
 
diff --git a/src/lagrangian/parcel/fvModels/clouds/clouds.H b/src/lagrangian/parcel/fvModels/clouds/clouds.H
index 617c98d5af..582fe94ec7 100644
--- a/src/lagrangian/parcel/fvModels/clouds/clouds.H
+++ b/src/lagrangian/parcel/fvModels/clouds/clouds.H
@@ -177,31 +177,31 @@ public:
             //- Add source to continuity equation
             virtual void addSup
             (
-                fvMatrix<scalar>& eqn,
-                const word& fieldName
+                const volScalarField& rho,
+                fvMatrix<scalar>& eqn
             ) const;
 
             //- Add source to enthalpy or species equation
             virtual void addSup
             (
                 const volScalarField& rho,
-                fvMatrix<scalar>& eqn,
-                const word& fieldName
+                const volScalarField& field,
+                fvMatrix<scalar>& eqn
             ) const;
 
             //- Add source to incompressible momentum equation
             virtual void addSup
             (
-                fvMatrix<vector>& eqn,
-                const word& fieldName
+                const volVectorField& U,
+                fvMatrix<vector>& eqn
             ) const;
 
             //- Add source to compressible momentum equation
             virtual void addSup
             (
                 const volScalarField& rho,
-                fvMatrix<vector>& eqn,
-                const word& fieldName
+                const volVectorField& U,
+                fvMatrix<vector>& eqn
             ) const;
 
 
diff --git a/src/radiationModels/fvModels/radiation/radiation.C b/src/radiationModels/fvModels/radiation/radiation.C
index 5e6abc22e8..a8b92c26fc 100644
--- a/src/radiationModels/fvModels/radiation/radiation.C
+++ b/src/radiationModels/fvModels/radiation/radiation.C
@@ -81,8 +81,8 @@ Foam::wordList Foam::fv::radiation::addSupFields() const
 void Foam::fv::radiation::addSup
 (
     const volScalarField& rho,
-    fvMatrix<scalar>& eqn,
-    const word& fieldName
+    const volScalarField& he,
+    fvMatrix<scalar>& eqn
 ) const
 {
     const basicThermo& thermo =
diff --git a/src/radiationModels/fvModels/radiation/radiation.H b/src/radiationModels/fvModels/radiation/radiation.H
index 8ccf33068b..d62b2f2828 100644
--- a/src/radiationModels/fvModels/radiation/radiation.H
+++ b/src/radiationModels/fvModels/radiation/radiation.H
@@ -106,8 +106,8 @@ public:
             virtual void addSup
             (
                 const volScalarField& rho,
-                fvMatrix<scalar>& eqn,
-                const word& fieldName
+                const volScalarField& he,
+                fvMatrix<scalar>& eqn
             ) const;
 
 
diff --git a/src/randomProcesses/OUForce/OUForce.C b/src/randomProcesses/OUForce/OUForce.C
index 2614553b26..eb17f95593 100644
--- a/src/randomProcesses/OUForce/OUForce.C
+++ b/src/randomProcesses/OUForce/OUForce.C
@@ -84,8 +84,8 @@ Foam::wordList Foam::fv::OUForce::addSupFields() const
 
 void Foam::fv::OUForce::addSup
 (
-    fvMatrix<vector>& eqn,
-    const word& fieldName
+    const volVectorField& field,
+    fvMatrix<vector>& eqn
 ) const
 {
     eqn += volVectorField::Internal
diff --git a/src/randomProcesses/OUForce/OUForce.H b/src/randomProcesses/OUForce/OUForce.H
index e7d39d6ca0..5b740487cd 100644
--- a/src/randomProcesses/OUForce/OUForce.H
+++ b/src/randomProcesses/OUForce/OUForce.H
@@ -133,8 +133,8 @@ public:
             //- Add explicit contribution to incompressible momentum equation
             virtual void addSup
             (
-                fvMatrix<vector>& eqn,
-                const word& fieldName
+                const volVectorField& field,
+                fvMatrix<vector>& eqn
             ) const;
 
 
diff --git a/src/waves/fvModels/isotropicDamping/isotropicDamping.C b/src/waves/fvModels/isotropicDamping/isotropicDamping.C
index 80e33e9ba0..78621893e0 100644
--- a/src/waves/fvModels/isotropicDamping/isotropicDamping.C
+++ b/src/waves/fvModels/isotropicDamping/isotropicDamping.C
@@ -96,8 +96,8 @@ Foam::wordList Foam::fv::isotropicDamping::addSupFields() const
 
 void Foam::fv::isotropicDamping::addSup
 (
-    fvMatrix<vector>& eqn,
-    const word& fieldName
+    const volVectorField& U,
+    fvMatrix<vector>& eqn
 ) const
 {
     add(this->forceCoeff(), eqn);
@@ -107,11 +107,11 @@ void Foam::fv::isotropicDamping::addSup
 void Foam::fv::isotropicDamping::addSup
 (
     const volScalarField& rho,
-    fvMatrix<vector>& eqn,
-    const word& fieldName
+    const volVectorField& U,
+    fvMatrix<vector>& eqn
 ) const
 {
-    add(rho*forceCoeff(), eqn);
+    add(rho*this->forceCoeff(), eqn);
 }
 
 
@@ -119,8 +119,8 @@ void Foam::fv::isotropicDamping::addSup
 (
     const volScalarField& alpha,
     const volScalarField& rho,
-    fvMatrix<vector>& eqn,
-    const word& fieldName
+    const volVectorField& U,
+    fvMatrix<vector>& eqn
 ) const
 {
     add(alpha()*rho()*this->forceCoeff(), eqn);
diff --git a/src/waves/fvModels/isotropicDamping/isotropicDamping.H b/src/waves/fvModels/isotropicDamping/isotropicDamping.H
index f43dc1ac21..8419dac915 100644
--- a/src/waves/fvModels/isotropicDamping/isotropicDamping.H
+++ b/src/waves/fvModels/isotropicDamping/isotropicDamping.H
@@ -162,16 +162,16 @@ public:
             //- Source term to momentum equation
             virtual void addSup
             (
-                fvMatrix<vector>& eqn,
-                const word& fieldName
+                const volVectorField& U,
+                fvMatrix<vector>& eqn
             ) const;
 
             //- Source term to compressible momentum equation
             virtual void addSup
             (
                 const volScalarField& rho,
-                fvMatrix<vector>& eqn,
-                const word& fieldName
+                const volVectorField& U,
+                fvMatrix<vector>& eqn
             ) const;
 
             //- Source term to phase momentum equation
@@ -179,8 +179,8 @@ public:
             (
                 const volScalarField& alpha,
                 const volScalarField& rho,
-                fvMatrix<vector>& eqn,
-                const word& fieldName
+                const volVectorField& U,
+                fvMatrix<vector>& eqn
             ) const;
 
 
diff --git a/src/waves/fvModels/verticalDamping/verticalDamping.C b/src/waves/fvModels/verticalDamping/verticalDamping.C
index 4f4cdb650d..5d04867efd 100644
--- a/src/waves/fvModels/verticalDamping/verticalDamping.C
+++ b/src/waves/fvModels/verticalDamping/verticalDamping.C
@@ -90,22 +90,22 @@ Foam::wordList Foam::fv::verticalDamping::addSupFields() const
 
 void Foam::fv::verticalDamping::addSup
 (
-    fvMatrix<vector>& eqn,
-    const word& fieldName
+    const volVectorField& U,
+    fvMatrix<vector>& eqn
 ) const
 {
-    add(eqn.psi(), eqn);
+    add(U, eqn);
 }
 
 
 void Foam::fv::verticalDamping::addSup
 (
     const volScalarField& rho,
-    fvMatrix<vector>& eqn,
-    const word& fieldName
+    const volVectorField& U,
+    fvMatrix<vector>& eqn
 ) const
 {
-    add(rho*eqn.psi(), eqn);
+    add(rho*U, eqn);
 }
 
 
@@ -113,11 +113,11 @@ void Foam::fv::verticalDamping::addSup
 (
     const volScalarField& alpha,
     const volScalarField& rho,
-    fvMatrix<vector>& eqn,
-    const word& fieldName
+    const volVectorField& U,
+    fvMatrix<vector>& eqn
 ) const
 {
-    add(alpha*rho*eqn.psi(), eqn);
+    add(alpha*rho*U, eqn);
 }
 
 
diff --git a/src/waves/fvModels/verticalDamping/verticalDamping.H b/src/waves/fvModels/verticalDamping/verticalDamping.H
index b73927ac1c..2e26debfc2 100644
--- a/src/waves/fvModels/verticalDamping/verticalDamping.H
+++ b/src/waves/fvModels/verticalDamping/verticalDamping.H
@@ -173,16 +173,16 @@ public:
             //- Source term to momentum equation
             virtual void addSup
             (
-                fvMatrix<vector>& eqn,
-                const word& fieldName
+                const volVectorField& U,
+                fvMatrix<vector>& eqn
             ) const;
 
             //- Source term to compressible momentum equation
             virtual void addSup
             (
                 const volScalarField& rho,
-                fvMatrix<vector>& eqn,
-                const word& fieldName
+                const volVectorField& U,
+                fvMatrix<vector>& eqn
             ) const;
 
             //- Source term to phase momentum equation
@@ -190,8 +190,8 @@ public:
             (
                 const volScalarField& alpha,
                 const volScalarField& rho,
-                fvMatrix<vector>& eqn,
-                const word& fieldName
+                const volVectorField& U,
+                fvMatrix<vector>& eqn
             ) const;
 
 
diff --git a/src/waves/fvModels/waveForcing/waveForcing.C b/src/waves/fvModels/waveForcing/waveForcing.C
index ce480e7174..9be3be83fc 100644
--- a/src/waves/fvModels/waveForcing/waveForcing.C
+++ b/src/waves/fvModels/waveForcing/waveForcing.C
@@ -127,11 +127,11 @@ Foam::wordList Foam::fv::waveForcing::addSupFields() const
 
 void Foam::fv::waveForcing::addSup
 (
-    fvMatrix<scalar>& eqn,
-    const word& fieldName
+    const volScalarField& alpha,
+    fvMatrix<scalar>& eqn
 ) const
 {
-    if (fieldName == alphaName_)
+    if (alpha.name() == alphaName_ && &alpha == &eqn.psi())
     {
         const volScalarField::Internal forceCoeff(this->forceCoeff());
 
@@ -144,11 +144,11 @@ void Foam::fv::waveForcing::addSup
 void Foam::fv::waveForcing::addSup
 (
     const volScalarField& rho,
-    fvMatrix<vector>& eqn,
-    const word& fieldName
+    const volVectorField& U,
+    fvMatrix<vector>& eqn
 ) const
 {
-    if (fieldName == UName_)
+    if (U.name() == UName_ && &U == &eqn.psi())
     {
         const volScalarField::Internal forceCoeff(rho*this->forceCoeff());
 
diff --git a/src/waves/fvModels/waveForcing/waveForcing.H b/src/waves/fvModels/waveForcing/waveForcing.H
index 8c8cc5bc3f..2fde935542 100644
--- a/src/waves/fvModels/waveForcing/waveForcing.H
+++ b/src/waves/fvModels/waveForcing/waveForcing.H
@@ -193,16 +193,16 @@ public:
             //- Source term to VoF phase-fraction equation
             virtual void addSup
             (
-                fvMatrix<scalar>& eqn,
-                const word& fieldName
+                const volScalarField& alpha,
+                fvMatrix<scalar>& eqn
             ) const;
 
             //- Source term to momentum equation
             virtual void addSup
             (
                 const volScalarField& rho,
-                fvMatrix<vector>& eqn,
-                const word& fieldName
+                const volVectorField& U,
+                fvMatrix<vector>& eqn
             ) const;
 
 
diff --git a/tutorials/compressibleVoF/damBreakInjection/Allclean b/tutorials/compressibleVoF/damBreakInjection/Allclean
new file mode 100755
index 0000000000..74dacf475b
--- /dev/null
+++ b/tutorials/compressibleVoF/damBreakInjection/Allclean
@@ -0,0 +1,11 @@
+#!/bin/sh
+cd ${0%/*} || exit 1    # Run from this directory
+
+# Source tutorial clean functions
+. $WM_PROJECT_DIR/bin/tools/CleanFunctions
+
+cleanVoFCase && rm -rf 0 system
+
+find constant -type f -not -name fvModels.injection -delete
+
+#------------------------------------------------------------------------------
diff --git a/tutorials/compressibleVoF/damBreakInjection/Allrun b/tutorials/compressibleVoF/damBreakInjection/Allrun
new file mode 100755
index 0000000000..845acdbeb0
--- /dev/null
+++ b/tutorials/compressibleVoF/damBreakInjection/Allrun
@@ -0,0 +1,18 @@
+#!/bin/sh
+cd ${0%/*} || exit 1    # Run from this directory
+
+# Source tutorial run functions
+. $WM_PROJECT_DIR/bin/tools/RunFunctions
+
+# Copy the source case and add the mass source model
+isTest "$@" && path=.. || path=$FOAM_TUTORIALS/compressibleVoF
+cp -rn $path/damBreak/0 $path/damBreak/constant $path/damBreak/system .
+runApplication foamDictionary constant/fvModels \
+    -entry injection -dict -merge constant/fvModels.injection
+
+# Run
+runApplication blockMesh
+runApplication setFields
+runApplication $(getApplication)
+
+#------------------------------------------------------------------------------
diff --git a/tutorials/compressibleVoF/damBreakInjection/constant/fvModels.injection b/tutorials/compressibleVoF/damBreakInjection/constant/fvModels.injection
new file mode 100644
index 0000000000..bbff0b8204
--- /dev/null
+++ b/tutorials/compressibleVoF/damBreakInjection/constant/fvModels.injection
@@ -0,0 +1,37 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     | Website:  https://openfoam.org
+    \\  /    A nd           | Version:  dev
+     \\/     M anipulation  |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    format      ascii;
+    class       dictionary;
+    location    "constant";
+    object      fvModels.injection;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+injection
+{
+    type            massSource;
+
+    points          ((0.438 0.438 0.0073));
+
+    massFlowRate    0.3;
+
+    phase           water;
+
+    fieldValues
+    {
+        T.water         350;
+        U               (-1 0 0);
+        k               0.1;
+        epsilon         0.1;
+    }
+}
+
+
+//************************************************************************* //