thermophysicalModels: Primitive mixture classes

Mixture classes (e.g., pureMixtrure, coefficientMulticomponentMixture),
now have no fvMesh or volScalarField dependence. They operate on
primitive values only. All the fvMesh-dependent functionality has been
moved into the base thermodynamic classes. The 'composition()' access
function has been removed from multi-component thermo models. Functions
that were once provided by composition base classes such as
basicSpecieMixture and basicCombustionMixture are now implemented
directly in the relevant multi-component thermo base class.

applications/legacy/combustion/PDRFoam/bEqn.H

@@ -94,9 +94,9 @@ if (ign.ignited())
 
     Info<< "min(b) = " << min(b).value() << endl;
 
-    if (composition.contains("ft"))
+    if (thermo.containsSpecie("ft"))
     {
-        volScalarField& ft = composition.Y("ft");
+        volScalarField& ft = thermo.Y("ft");
 
         Info<< "Combustion progress = "
             << 100*(1.0 - b)().weightedAverage(mesh.V()*ft).value() << "%"

applications/legacy/combustion/PDRFoam/createFields.H

@@ -7,8 +7,6 @@ autoPtr<psiuMulticomponentThermo> pThermo
 psiuMulticomponentThermo& thermo = pThermo();
 thermo.validate(args.executable(), "ha", "ea");
 
-basicCombustionMixture& composition = thermo.composition();
-
 volScalarField rho
 (
     IOobject
@@ -24,7 +22,7 @@ volScalarField rho
 
 volScalarField& p = thermo.p();
 
-volScalarField& b = composition.Y("b");
+volScalarField& b = thermo.Y("b");
 Info<< "min(b) = " << min(b).value() << endl;
 
 Info<< "\nReading field U\n" << endl;
@@ -229,9 +227,9 @@ volScalarField St
 
 multivariateSurfaceInterpolationScheme<scalar>::fieldTable fields;
 
-if (composition.contains("ft"))
+if (thermo.containsSpecie("ft"))
 {
-    fields.add(composition.Y("ft"));
+    fields.add(thermo.Y("ft"));
 }
 
 fields.add(b);

applications/legacy/combustion/PDRFoam/ftEqn.H

@@ -1,6 +1,6 @@
-if (composition.contains("ft"))
+if (thermo.containsSpecie("ft"))
 {
-    volScalarField& ft = composition.Y("ft");
+    volScalarField& ft = thermo.Y("ft");
 
     solve
     (

applications/legacy/combustion/PDRFoam/laminarFlameSpeed/SCOPE/SCOPELaminarFlameSpeed.C

@@ -368,10 +368,9 @@ Foam::tmp<Foam::volScalarField> Foam::laminarFlameSpeedModels::SCOPE::Ma
 Foam::tmp<Foam::volScalarField>
 Foam::laminarFlameSpeedModels::SCOPE::Ma() const
 {
-    if (psiuMulticomponentThermo_.composition().contains("ft"))
+    if (psiuMulticomponentThermo_.containsSpecie("ft"))
     {
-        const volScalarField& ft =
-            psiuMulticomponentThermo_.composition().Y("ft");
+        const volScalarField& ft = psiuMulticomponentThermo_.Y("ft");
 
         return Ma
         (
@@ -403,10 +402,9 @@ Foam::laminarFlameSpeedModels::SCOPE::Ma() const
 Foam::tmp<Foam::volScalarField>
 Foam::laminarFlameSpeedModels::SCOPE::operator()() const
 {
-    if (psiuMulticomponentThermo_.composition().contains("ft"))
+    if (psiuMulticomponentThermo_.containsSpecie("ft"))
     {
-        const volScalarField& ft =
-            psiuMulticomponentThermo_.composition().Y("ft");
+        const volScalarField& ft = psiuMulticomponentThermo_.Y("ft");
 
         return Su0pTphi
         (

applications/modules/XiFluid/XiFluid.C

@@ -51,9 +51,7 @@ Foam::solvers::XiFluid::XiFluid(fvMesh& mesh)
 
     thermo_(refCast<psiuMulticomponentThermo>(isothermalFluid::thermo_)),
 
-    composition(thermo_.composition()),
-
-    b_(composition.Y("b")),
+    b_(thermo_.Y("b")),
 
     unstrainedLaminarFlameSpeed(laminarFlameSpeed::New(thermo_)),
 
@@ -144,9 +142,9 @@ Foam::solvers::XiFluid::XiFluid(fvMesh& mesh)
 {
     thermo.validate(type(), "ha", "ea");
 
-    if (composition.contains("ft"))
+    if (thermo_.containsSpecie("ft"))
     {
-        fields.add(composition.Y("ft"));
+        fields.add(thermo_.Y("ft"));
     }
 
     fields.add(b);

applications/modules/XiFluid/XiFluid.H

@@ -106,11 +106,7 @@ protected:
 
     // Composition
 
-        //- Reference to the combustion mixture
-        basicCombustionMixture& composition;
-
         //- Reference to the combustion regress variable
-        //  obtained from the combustion mixture
         volScalarField& b_;
 
         //- Set of fields used for the multivariate convection scheme

applications/modules/XiFluid/thermophysicalPredictor.C

@@ -38,7 +38,7 @@ void Foam::solvers::XiFluid::ftSolve
     const fv::convectionScheme<scalar>& mvConvection
 )
 {
-    volScalarField& ft = composition.Y("ft");
+    volScalarField& ft = thermo_.Y("ft");
 
     fvScalarMatrix ftEqn
     (
@@ -569,7 +569,7 @@ void Foam::solvers::XiFluid::thermophysicalPredictor()
         )
     );
 
-    if (composition.contains("ft"))
+    if (thermo_.containsSpecie("ft"))
     {
         ftSolve(mvConvection());
     }

applications/modules/multicomponentFluid/multicomponentFluid.C

@@ -51,9 +51,7 @@ Foam::solvers::multicomponentFluid::multicomponentFluid(fvMesh& mesh)
 
     thermo_(refCast<fluidMulticomponentThermo>(isothermalFluid::thermo_)),
 
-    composition(thermo_.composition()),
-
-    Y_(composition.Y()),
+    Y_(thermo_.Y()),
 
     reaction(combustionModel::New(thermo_, momentumTransport())),
 

applications/modules/multicomponentFluid/multicomponentFluid.H

@@ -86,7 +86,6 @@ protected:
 
     // Composition
 
-        basicSpecieMixture& composition;
         PtrList<volScalarField>& Y_;
 
 

applications/modules/multicomponentFluid/setRDeltaT.C

@@ -106,7 +106,7 @@ void Foam::solvers::multicomponentFluid::setRDeltaT()
 
         forAll(Y, i)
         {
-            if (composition.solve(i))
+            if (thermo_.solveSpecie(i))
             {
                 volScalarField& Yi = Y_[i];
 

applications/modules/multicomponentFluid/thermophysicalPredictor.C

@@ -45,7 +45,7 @@ void Foam::solvers::multicomponentFluid::thermophysicalPredictor()
 
     forAll(Y, i)
     {
-        if (composition.solve(i))
+        if (thermo_.solveSpecie(i))
         {
             volScalarField& Yi = Y_[i];
 
@@ -69,7 +69,7 @@ void Foam::solvers::multicomponentFluid::thermophysicalPredictor()
         }
     }
 
-    composition.normalise();
+    thermo_.normaliseY();
 
 
     volScalarField& he = thermo_.he();

applications/modules/multiphaseEuler/interfacialCompositionModels/interfaceCompositionModels/Henry/Henry.C

@@ -100,12 +100,12 @@ Foam::tmp<Foam::volScalarField> Foam::interfaceCompositionModels::Henry::Yf
 
         return
             k_[index]
-           *otherComposition().Y(speciesName)
+           *otherMulticomponentThermo().Y(speciesName)
            *otherThermo().rho()/thermo().rho();
     }
     else
     {
-        return YSolvent_*composition().Y(speciesName);
+        return YSolvent_*thermo().Y(speciesName);
     }
 }
 

applications/modules/multiphaseEuler/interfacialCompositionModels/interfaceCompositionModels/Raoult/Raoult.C

@@ -110,11 +110,11 @@ void Foam::interfaceCompositionModels::Raoult::update(const volScalarField& Tf)
         iter()->update(Tf);
 
         YNonVapour_ -=
-            otherComposition().Y(iter.key())
+            otherMulticomponentThermo().Y(iter.key())
            *iter()->Yf(iter.key(), Tf);
 
         YNonVapourPrime_ -=
-            otherComposition().Y(iter.key())
+            otherMulticomponentThermo().Y(iter.key())
            *iter()->YfPrime(iter.key(), Tf);
     }
 }
@@ -129,12 +129,12 @@ Foam::tmp<Foam::volScalarField> Foam::interfaceCompositionModels::Raoult::Yf
     if (species().found(speciesName))
     {
         return
-            otherComposition().Y(speciesName)
+            otherMulticomponentThermo().Y(speciesName)
            *speciesModels_[speciesName]->Yf(speciesName, Tf);
     }
     else
     {
-        return composition().Y(speciesName)*YNonVapour_;
+        return thermo().Y(speciesName)*YNonVapour_;
     }
 }
 
@@ -149,12 +149,12 @@ Foam::interfaceCompositionModels::Raoult::YfPrime
     if (species().found(speciesName))
     {
         return
-            otherComposition().Y(speciesName)
+            otherMulticomponentThermo().Y(speciesName)
            *speciesModels_[speciesName]->YfPrime(speciesName, Tf);
     }
     else
     {
-        return composition().Y(speciesName)*YNonVapourPrime_;
+        return thermo().Y(speciesName)*YNonVapourPrime_;
     }
 }
 

applications/modules/multiphaseEuler/interfacialCompositionModels/interfaceCompositionModels/interfaceCompositionModel/interfaceCompositionModel.C

@@ -74,9 +74,9 @@ Foam::tmp<Foam::volScalarField> Foam::interfaceCompositionModel::dY
     const volScalarField& Tf
 ) const
 {
-    const label speciei = composition().species()[speciesName];
+    const label speciei = thermo().species()[speciesName];
 
-    return Yf(speciesName, Tf) - composition().Y()[speciei];
+    return Yf(speciesName, Tf) - thermo().Y()[speciei];
 }
 
 
@@ -95,16 +95,16 @@ Foam::tmp<Foam::volScalarField> Foam::interfaceCompositionModel::D
     const word& speciesName
 ) const
 {
-    const label speciei = composition().species()[speciesName];
+    const label speciei = thermo().species()[speciesName];
     const volScalarField& p(thermo_.p());
     const volScalarField& T(thermo_.T());
 
     return volScalarField::New
     (
         IOobject::groupName("D" + speciesName, interface_.name()),
-        composition().kappa(speciei, p, T)
-       /composition().Cp(speciei, p, T)
-       /composition().rho(speciei, p, T)
+        thermo().kappai(speciei, p, T)
+       /thermo().Cpi(speciei, p, T)
+       /thermo().rhoi(speciei, p, T)
        /Le_
     );
 }

applications/modules/multiphaseEuler/interfacialCompositionModels/interfaceCompositionModels/interfaceCompositionModel/interfaceCompositionModel.H

@@ -132,17 +132,12 @@ public:
             //- Return the thermo
             inline const rhoMulticomponentThermo& thermo() const;
 
-            //- Return the composition
-            inline const basicSpecieMixture& composition() const;
-
             //- Return the other thermo
             inline const rhoThermo& otherThermo() const;
 
-            //- Return whether the other side has a multi-specie composition
-            inline bool otherHasComposition() const;
-
-            //- Return the other composition
-            inline const basicSpecieMixture& otherComposition() const;
+            //- Return the other multicomponent thermo
+            inline const rhoMulticomponentThermo&
+                otherMulticomponentThermo() const;
 
 
         // Evaluation

applications/modules/multiphaseEuler/interfacialCompositionModels/interfaceCompositionModels/interfaceCompositionModel/interfaceCompositionModelI.H

@@ -24,6 +24,7 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "interfaceCompositionModel.H"
+#include "rhoMulticomponentThermo.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -47,29 +48,16 @@ Foam::interfaceCompositionModel::thermo() const
 }
 
 
-const Foam::basicSpecieMixture&
-Foam::interfaceCompositionModel::composition() const
-{
-    return thermo_.composition();
-}
-
-
 const Foam::rhoThermo& Foam::interfaceCompositionModel::otherThermo() const
 {
     return otherThermo_;
 }
 
 
-bool Foam::interfaceCompositionModel::otherHasComposition() const
+const Foam::rhoMulticomponentThermo&
+Foam::interfaceCompositionModel::otherMulticomponentThermo() const
 {
-    return isA<rhoMulticomponentThermo>(otherThermo_);
-}
-
-
-const Foam::basicSpecieMixture&
-Foam::interfaceCompositionModel::otherComposition() const
-{
-    return refCast<const rhoMulticomponentThermo>(otherThermo_).composition();
+    return refCast<const rhoMulticomponentThermo>(otherThermo_);
 }
 
 

applications/modules/multiphaseEuler/interfacialCompositionModels/interfaceCompositionModels/nonRandomTwoLiquid/nonRandomTwoLiquid.C

@@ -87,8 +87,8 @@ Foam::interfaceCompositionModels::nonRandomTwoLiquid::nonRandomTwoLiquid
     species1Name_ = species()[0];
     species2Name_ = species()[1];
 
-    species1Index_ = composition().species()[species1Name_];
-    species2Index_ = composition().species()[species2Name_];
+    species1Index_ = thermo().species()[species1Name_];
+    species2Index_ = thermo().species()[species2Name_];
 
     alpha12_ = dimensionedScalar
     (
@@ -171,25 +171,25 @@ void Foam::interfaceCompositionModels::nonRandomTwoLiquid::update
 
     const volScalarField X1
     (
-        composition().Y(species1Index_)
+        thermo().Y(species1Index_)
        *W
        /dimensionedScalar
         (
             "W",
             dimMass/dimMoles,
-            composition().Wi(species1Index_)
+            thermo().Wi(species1Index_)
         )
     );
 
     const volScalarField X2
     (
-        composition().Y(species2Index_)
+        thermo().Y(species2Index_)
        *W
        /dimensionedScalar
         (
             "W",
             dimMass/dimMoles,
-            composition().Wi(species2Index_)
+            thermo().Wi(species2Index_)
         )
     );
 
@@ -233,21 +233,21 @@ Foam::interfaceCompositionModels::nonRandomTwoLiquid::Yf
     if (speciesName == species1Name_)
     {
         return
-            otherComposition().Y(speciesName)
+            otherMulticomponentThermo().Y(speciesName)
            *speciesModel1_->Yf(speciesName, Tf)
            *gamma1_;
     }
     else if (speciesName == species2Name_)
     {
         return
-            otherComposition().Y(speciesName)
+            otherMulticomponentThermo().Y(speciesName)
            *speciesModel2_->Yf(speciesName, Tf)
            *gamma2_;
     }
     else
     {
         return
-            composition().Y(speciesName)
+            thermo().Y(speciesName)
            *(scalar(1) - Yf(species1Name_, Tf) - Yf(species2Name_, Tf));
     }
 }
@@ -263,21 +263,21 @@ Foam::interfaceCompositionModels::nonRandomTwoLiquid::YfPrime
     if (speciesName == species1Name_)
     {
         return
-            otherComposition().Y(speciesName)
+            otherMulticomponentThermo().Y(speciesName)
            *speciesModel1_->YfPrime(speciesName, Tf)
            *gamma1_;
     }
     else if (speciesName == species2Name_)
     {
         return
-            otherComposition().Y(speciesName)
+            otherMulticomponentThermo().Y(speciesName)
            *speciesModel2_->YfPrime(speciesName, Tf)
            *gamma2_;
     }
     else
     {
         return
-          - composition().Y(speciesName)
+          - thermo().Y(speciesName)
            *(YfPrime(species1Name_, Tf) + YfPrime(species2Name_, Tf));
     }
 }

applications/modules/multiphaseEuler/interfacialCompositionModels/interfaceCompositionModels/saturated/saturated.C

@@ -52,7 +52,7 @@ Foam::interfaceCompositionModels::saturated::wRatioByP() const
     (
         "W",
         dimMass/dimMoles,
-        composition().Wi(saturatedIndex_)
+        thermo().Wi(saturatedIndex_)
     );
 
     return Wi/thermo().W()/thermo().p();
@@ -69,7 +69,7 @@ Foam::interfaceCompositionModels::saturated::saturated
 :
     interfaceCompositionModel(dict, interface),
     saturatedName_(species()[0]),
-    saturatedIndex_(composition().species()[saturatedName_]),
+    saturatedIndex_(thermo().species()[saturatedName_]),
     saturationModel_(saturationPressureModel::New("pSat", dict))
 {
     if (species().size() != 1)
@@ -108,12 +108,12 @@ Foam::tmp<Foam::volScalarField> Foam::interfaceCompositionModels::saturated::Yf
     }
     else
     {
-        const label speciesIndex = composition().species()[speciesName];
+        const label speciesIndex = thermo().species()[speciesName];
 
         return
-            composition().Y()[speciesIndex]
+            thermo().Y()[speciesIndex]
            *(scalar(1) - wRatioByP()*saturationModel_->pSat(Tf))
-           /max(scalar(1) - composition().Y()[saturatedIndex_], small);
+           /max(scalar(1) - thermo().Y()[saturatedIndex_], small);
     }
 }
 
@@ -131,12 +131,12 @@ Foam::interfaceCompositionModels::saturated::YfPrime
     }
     else
     {
-        const label speciesIndex = composition().species()[speciesName];
+        const label speciesIndex = thermo().species()[speciesName];
 
         return
-          - composition().Y()[speciesIndex]
+          - thermo().Y()[speciesIndex]
            *wRatioByP()*saturationModel_->pSatPrime(Tf)
-           /max(scalar(1) - composition().Y()[saturatedIndex_], small);
+           /max(scalar(1) - thermo().Y()[saturatedIndex_], small);
     }
 }
 

applications/modules/multiphaseEuler/multiphaseThermophysicalTransportModels/HeatTransferPhaseSystem/HeatTransferPhaseSystem.C

@@ -91,14 +91,14 @@ void Foam::HeatTransferPhaseSystem<BasePhaseSystem>::addDmidtHefs
         const phaseModel& phase2 = interface.phase2();
         const rhoThermo& thermo1 = phase1.thermo();
         const rhoThermo& thermo2 = phase2.thermo();
-        const basicSpecieMixture* compositionPtr1 =
+        const rhoMulticomponentThermo* mcThermoPtr1 =
             isA<rhoMulticomponentThermo>(thermo1)
-          ? &refCast<const rhoMulticomponentThermo>(thermo1).composition()
-          : static_cast<const basicSpecieMixture*>(nullptr);
-        const basicSpecieMixture* compositionPtr2 =
+          ? &refCast<const rhoMulticomponentThermo>(thermo1)
+          : static_cast<const rhoMulticomponentThermo*>(nullptr);
+        const rhoMulticomponentThermo* mcThermoPtr2 =
             isA<rhoMulticomponentThermo>(thermo2)
-          ? &refCast<const rhoMulticomponentThermo>(thermo2).composition()
-          : static_cast<const basicSpecieMixture*>(nullptr);
+          ? &refCast<const rhoMulticomponentThermo>(thermo2)
+          : static_cast<const rhoMulticomponentThermo*>(nullptr);
         const volScalarField& he1 = thermo1.he();
         const volScalarField& he2 = thermo2.he();
         const volScalarField hs1(thermo1.hs());
@@ -118,21 +118,21 @@ void Foam::HeatTransferPhaseSystem<BasePhaseSystem>::addDmidtHefs
 
             // Specie indices
             const label speciei1 =
-                compositionPtr1 ? compositionPtr1->species()[specie] : -1;
+                mcThermoPtr1 ? mcThermoPtr1->species()[specie] : -1;
             const label speciei2 =
-                compositionPtr2 ? compositionPtr2->species()[specie] : -1;
+                mcThermoPtr2 ? mcThermoPtr2->species()[specie] : -1;
 
             // Enthalpies
             const volScalarField hsi1
             (
-                compositionPtr1
-              ? compositionPtr1->Hs(speciei1, thermo1.p(), thermo1.T())
+                mcThermoPtr1
+              ? mcThermoPtr1->hsi(speciei1, thermo1.p(), thermo1.T())
               : tmp<volScalarField>(hs1)
             );
             const volScalarField hsi2
             (
-                compositionPtr2
-              ? compositionPtr2->Hs(speciei2, thermo2.p(), thermo2.T())
+                mcThermoPtr2
+              ? mcThermoPtr2->hsi(speciei2, thermo2.p(), thermo2.T())
               : tmp<volScalarField>(hs2)
             );
 
@@ -141,12 +141,12 @@ void Foam::HeatTransferPhaseSystem<BasePhaseSystem>::addDmidtHefs
             if (residualY_ > 0)
             {
                 tYi1 =
-                    compositionPtr1
-                  ? max(compositionPtr1->Y(speciei1), residualY_)
+                    mcThermoPtr1
+                  ? max(mcThermoPtr1->Y(speciei1), residualY_)
                   : volScalarField::New("Yi1", this->mesh(), one);
                 tYi2 =
-                    compositionPtr2
-                  ? max(compositionPtr2->Y(speciei2), residualY_)
+                    mcThermoPtr2
+                  ? max(mcThermoPtr2->Y(speciei2), residualY_)
                   : volScalarField::New("Yi2", this->mesh(), one);
             }
 
@@ -307,14 +307,14 @@ void Foam::HeatTransferPhaseSystem<BasePhaseSystem>::addDmidtHefsWithoutL
         const phaseModel& phase2 = interface.phase2();
         const rhoThermo& thermo1 = phase1.thermo();
         const rhoThermo& thermo2 = phase2.thermo();
-        const basicSpecieMixture* compositionPtr1 =
+        const rhoMulticomponentThermo* mcThermoPtr1 =
             isA<rhoMulticomponentThermo>(thermo1)
-          ? &refCast<const rhoMulticomponentThermo>(thermo1).composition()
-          : static_cast<const basicSpecieMixture*>(nullptr);
-        const basicSpecieMixture* compositionPtr2 =
+          ? &refCast<const rhoMulticomponentThermo>(thermo1)
+          : static_cast<const rhoMulticomponentThermo*>(nullptr);
+        const rhoMulticomponentThermo* mcThermoPtr2 =
             isA<rhoMulticomponentThermo>(thermo2)
-          ? &refCast<const rhoMulticomponentThermo>(thermo2).composition()
-          : static_cast<const basicSpecieMixture*>(nullptr);
+          ? &refCast<const rhoMulticomponentThermo>(thermo2)
+          : static_cast<const rhoMulticomponentThermo*>(nullptr);
         const volScalarField& he1 = thermo1.he();
         const volScalarField& he2 = thermo2.he();
         const volScalarField K1(phase1.K());
@@ -336,21 +336,21 @@ void Foam::HeatTransferPhaseSystem<BasePhaseSystem>::addDmidtHefsWithoutL
 
             // Specie indices
             const label speciei1 =
-                compositionPtr1 ? compositionPtr1->species()[specie] : -1;
+                mcThermoPtr1 ? mcThermoPtr1->species()[specie] : -1;
             const label speciei2 =
-                compositionPtr2 ? compositionPtr2->species()[specie] : -1;
+                mcThermoPtr2 ? mcThermoPtr2->species()[specie] : -1;
 
             // Interface enthalpies
             const volScalarField hsfi1
             (
-                compositionPtr1
-              ? compositionPtr1->Hs(speciei1, thermo1.p(), Tf)
+                mcThermoPtr1
+              ? mcThermoPtr1->hsi(speciei1, thermo1.p(), Tf)
               : tmp<volScalarField>(hsf1)
             );
             const volScalarField hsfi2
             (
-                compositionPtr2
-              ? compositionPtr2->Hs(speciei2, thermo2.p(), Tf)
+                mcThermoPtr2
+              ? mcThermoPtr2->hsi(speciei2, thermo2.p(), Tf)
               : tmp<volScalarField>(hsf2)
             );
 
@@ -359,12 +359,12 @@ void Foam::HeatTransferPhaseSystem<BasePhaseSystem>::addDmidtHefsWithoutL
             if (this->residualY_ > 0)
             {
                 tYi1 =
-                    compositionPtr1
-                  ? max(compositionPtr1->Y(speciei1), this->residualY_)
+                    mcThermoPtr1
+                  ? max(mcThermoPtr1->Y(speciei1), this->residualY_)
                   : volScalarField::New("Yi1", this->mesh(), one);
                 tYi2 =
-                    compositionPtr2
-                  ? max(compositionPtr2->Y(speciei2), this->residualY_)
+                    mcThermoPtr2
+                  ? max(mcThermoPtr2->Y(speciei2), this->residualY_)
                   : volScalarField::New("Yi2", this->mesh(), one);
             }
 
@@ -383,14 +383,14 @@ void Foam::HeatTransferPhaseSystem<BasePhaseSystem>::addDmidtHefsWithoutL
                     // Bulk enthalpies
                     const volScalarField hsi1
                     (
-                        compositionPtr1
-                      ? compositionPtr1->Hs(speciei1, thermo1.p(), thermo1.T())
+                        mcThermoPtr1
+                      ? mcThermoPtr1->hsi(speciei1, thermo1.p(), thermo1.T())
                       : thermo1.hs()
                     );
                     const volScalarField hsi2
                     (
-                        compositionPtr2
-                      ? compositionPtr2->Hs(speciei2, thermo2.p(), thermo2.T())
+                        mcThermoPtr2
+                      ? mcThermoPtr2->hsi(speciei2, thermo2.p(), thermo2.T())
                       : thermo2.hs()
                     );
 
@@ -597,30 +597,30 @@ Foam::HeatTransferPhaseSystem<BasePhaseSystem>::Li
 {
     const rhoThermo& thermo1 = interface.phase1().thermo();
     const rhoThermo& thermo2 = interface.phase2().thermo();
-    const basicSpecieMixture* compositionPtr1 =
+    const rhoMulticomponentThermo* mcThermoPtr1 =
         isA<rhoMulticomponentThermo>(thermo1)
-      ? &refCast<const rhoMulticomponentThermo>(thermo1).composition()
-      : static_cast<const basicSpecieMixture*>(nullptr);
-    const basicSpecieMixture* compositionPtr2 =
+      ? &refCast<const rhoMulticomponentThermo>(thermo1)
+      : static_cast<const rhoMulticomponentThermo*>(nullptr);
+    const rhoMulticomponentThermo* mcThermoPtr2 =
         isA<rhoMulticomponentThermo>(thermo2)
-      ? &refCast<const rhoMulticomponentThermo>(thermo2).composition()
-      : static_cast<const basicSpecieMixture*>(nullptr);
+      ? &refCast<const rhoMulticomponentThermo>(thermo2)
+      : static_cast<const rhoMulticomponentThermo*>(nullptr);
     const label speciei1 =
-        compositionPtr1 ? compositionPtr1->species()[specie] : -1;
+        mcThermoPtr1 ? mcThermoPtr1->species()[specie] : -1;
     const label speciei2 =
-        compositionPtr2 ? compositionPtr2->species()[specie] : -1;
+        mcThermoPtr2 ? mcThermoPtr2->species()[specie] : -1;
 
     // Interface enthalpies
     const volScalarField hafi1
     (
-        compositionPtr1
-      ? compositionPtr1->Ha(speciei1, thermo1.p(), Tf)
+        mcThermoPtr1
+      ? mcThermoPtr1->hai(speciei1, thermo1.p(), Tf)
       : thermo1.ha(thermo1.p(), Tf)
     );
     const volScalarField hafi2
     (
-        compositionPtr2
-      ? compositionPtr2->Ha(speciei2, thermo2.p(), Tf)
+        mcThermoPtr2
+      ? mcThermoPtr2->hai(speciei2, thermo2.p(), Tf)
       : thermo2.ha(thermo1.p(), Tf)
     );
 
@@ -635,14 +635,14 @@ Foam::HeatTransferPhaseSystem<BasePhaseSystem>::Li
             // Bulk enthalpies
             const volScalarField hai1
             (
-                compositionPtr1
-              ? compositionPtr1->Ha(speciei1, thermo1.p(), thermo1.T())
+                mcThermoPtr1
+              ? mcThermoPtr1->hai(speciei1, thermo1.p(), thermo1.T())
               : thermo1.ha()
             );
             const volScalarField hai2
             (
-                compositionPtr2
-              ? compositionPtr2->Ha(speciei2, thermo2.p(), thermo2.T())
+                mcThermoPtr2
+              ? mcThermoPtr2->hai(speciei2, thermo2.p(), thermo2.T())
               : thermo2.ha()
             );
 
@@ -670,18 +670,18 @@ Foam::HeatTransferPhaseSystem<BasePhaseSystem>::Li
 {
     const rhoThermo& thermo1 = interface.phase1().thermo();
     const rhoThermo& thermo2 = interface.phase2().thermo();
-    const basicSpecieMixture* compositionPtr1 =
+    const rhoMulticomponentThermo* mcThermoPtr1 =
         isA<rhoMulticomponentThermo>(thermo1)
-      ? &refCast<const rhoMulticomponentThermo>(thermo1).composition()
-      : static_cast<const basicSpecieMixture*>(nullptr);
-    const basicSpecieMixture* compositionPtr2 =
+      ? &refCast<const rhoMulticomponentThermo>(thermo1)
+      : static_cast<const rhoMulticomponentThermo*>(nullptr);
+    const rhoMulticomponentThermo* mcThermoPtr2 =
         isA<rhoMulticomponentThermo>(thermo2)
-      ? &refCast<const rhoMulticomponentThermo>(thermo2).composition()
-      : static_cast<const basicSpecieMixture*>(nullptr);
+      ? &refCast<const rhoMulticomponentThermo>(thermo2)
+      : static_cast<const rhoMulticomponentThermo*>(nullptr);
     const label speciei1 =
-        compositionPtr1 ? compositionPtr1->species()[specie] : -1;
+        mcThermoPtr1 ? mcThermoPtr1->species()[specie] : -1;
     const label speciei2 =
-        compositionPtr2 ? compositionPtr2->species()[specie] : -1;
+        mcThermoPtr2 ? mcThermoPtr2->species()[specie] : -1;
 
     const scalarField p1(UIndirectList<scalar>(thermo1.p(), cells));
     const scalarField p2(UIndirectList<scalar>(thermo2.p(), cells));
@@ -689,14 +689,14 @@ Foam::HeatTransferPhaseSystem<BasePhaseSystem>::Li
     // Interface enthalpies
     const scalarField hafi1
     (
-        compositionPtr1
-      ? compositionPtr1->Ha(speciei1, p1, Tf)
+        mcThermoPtr1
+      ? mcThermoPtr1->hai(speciei1, p1, Tf)
       : thermo1.ha(Tf, cells)
     );
     const scalarField hafi2
     (
-        compositionPtr2
-      ? compositionPtr2->Ha(speciei2, p2, Tf)
+        mcThermoPtr2
+      ? mcThermoPtr2->hai(speciei2, p2, Tf)
       : thermo2.ha(Tf, cells)
     );
 
@@ -714,14 +714,14 @@ Foam::HeatTransferPhaseSystem<BasePhaseSystem>::Li
             // Bulk enthalpies
             const scalarField hai1
             (
-                compositionPtr1
-              ? compositionPtr1->Ha(speciei1, p1, T1)
+                mcThermoPtr1
+              ? mcThermoPtr1->hai(speciei1, p1, T1)
               : thermo1.ha(T1, cells)
             );
             const scalarField hai2
             (
-                compositionPtr2
-              ? compositionPtr2->Ha(speciei2, p2, T2)
+                mcThermoPtr2
+              ? mcThermoPtr2->hai(speciei2, p2, T2)
               : thermo2.ha(T2, cells)
             );
 

applications/modules/multiphaseEuler/phaseSystems/phaseModel/MulticomponentPhaseModel/MulticomponentPhaseModel.C

@@ -48,11 +48,11 @@ Foam::MulticomponentPhaseModel<BasePhaseModel>::MulticomponentPhaseModel
 :
     BasePhaseModel(fluid, phaseName, referencePhase, index)
 {
-    PtrList<volScalarField>& Y = this->thermo_->composition().Y();
+    PtrList<volScalarField>& Y = this->thermo_->Y();
 
     forAll(Y, i)
     {
-        if (this->thermo_->composition().solve(i))
+        if (this->thermo_->solveSpecie(i))
         {
             const label j = YActive_.size();
             YActive_.resize(j + 1);
@@ -74,7 +74,7 @@ Foam::MulticomponentPhaseModel<BasePhaseModel>::~MulticomponentPhaseModel()
 template<class BasePhaseModel>
 void Foam::MulticomponentPhaseModel<BasePhaseModel>::correctSpecies()
 {
-    this->thermo_->composition().normalise();
+    this->thermo_->normaliseY();
     BasePhaseModel::correctSpecies();
 }
 
@@ -121,7 +121,7 @@ template<class BasePhaseModel>
 const Foam::PtrList<Foam::volScalarField>&
 Foam::MulticomponentPhaseModel<BasePhaseModel>::Y() const
 {
-    return this->thermo_->composition().Y();
+    return this->thermo_->Y();
 }
 
 
@@ -129,7 +129,7 @@ template<class BasePhaseModel>
 const Foam::volScalarField&
 Foam::MulticomponentPhaseModel<BasePhaseModel>::Y(const word& name) const
 {
-    return this->thermo_->composition().Y(name);
+    return this->thermo_->Y(name);
 }
 
 
@@ -137,7 +137,7 @@ template<class BasePhaseModel>
 Foam::PtrList<Foam::volScalarField>&
 Foam::MulticomponentPhaseModel<BasePhaseModel>::YRef()
 {
-    return this->thermo_->composition().Y();
+    return this->thermo_->Y();
 }
 
 

applications/solvers/chemFoam/createFieldRefs.H

@@ -1,5 +1,4 @@
 basicChemistryModel& chemistry = pChemistry();
 scalar dtChem = min(chemistry.deltaTChem()[0], runTime.deltaT().value());
-basicSpecieMixture& composition = thermo.composition();
-PtrList<volScalarField>& Y = composition.Y();
+PtrList<volScalarField>& Y = thermo.Y();
 volScalarField& p = thermo.p();

applications/solvers/chemFoam/thermoTypeFunctions.H

@@ -28,13 +28,11 @@ namespace Foam
 
 scalarList W(const fluidMulticomponentThermo& thermo)
 {
-    const basicSpecieMixture& composition = thermo.composition();
-
-    scalarList W(composition.Y().size());
+    scalarList W(thermo.Y().size());
 
     forAll(W, i)
     {
-        W[i] = composition.Wi(i);
+        W[i] = thermo.Wi(i);
     }
 
     return W;
@@ -49,13 +47,11 @@ scalar h0
     const scalar T
 )
 {
-    const basicSpecieMixture& composition = thermo.composition();
-
     scalar h0 = 0;
 
     forAll(Y, i)
     {
-        h0 += Y[i]*composition.Hs(i, p, T);
+        h0 += Y[i]*thermo.hsi(i, p, T);
     }
 
     return h0;

etc/codeTemplates/dynamicCode/fluidMulticomponentThermoTemplate.C

@@ -23,13 +23,9 @@ License
 
 \*---------------------------------------------------------------------------*/
 
-#include "forThermo.H"
-#include "makeMulticomponentThermo.H"
-
+// Specie
 #include "${specie}.H"
 
-#include "thermo.H"
-
 // EoS
 #include "${equationOfState}.H"
 
@@ -47,6 +43,12 @@ License
 // Mixture
 #include "${mixture}.H"
 
+#include "thermo.H"
+#include "heMulticomponentThermo.H"
+#include "heFluidMulticomponentThermo.H"
+#include "typedefThermo.H"
+#include "makeThermo.H"
+
 
 // * * * * * * * * * * * * * * * Global Functions  * * * * * * * * * * * * * //
 
@@ -75,19 +77,37 @@ extern "C"
 
 namespace Foam
 {
-    forThermo
+    typedefThermo
     (
         ${transport}Transport,
         ${energy},
         ${thermo}Thermo,
         ${equationOfState},
-        ${specie},
-        makeMulticomponentThermo,
+        ${specie}
+    );
+
+    defineThermo
+    (
         ${typeBase},
-        ${type},
-        ${mixture}
+        ${mixture},
+        ${transport}Transport${energy}${thermo}Thermo${equationOfState}${specie}
+    );
+
+    addThermo
+    (
+        fluidMulticomponentThermo,
+        ${typeBase},
+        ${mixture},
+        ${transport}Transport${energy}${thermo}Thermo${equationOfState}${specie}
+    );
+
+    addThermo
+    (
+        ${typeBase},
+        ${typeBase},
+        ${mixture},
+        ${transport}Transport${energy}${thermo}Thermo${equationOfState}${specie}
     );
 }
 
-
 // ************************************************************************* //

etc/codeTemplates/dynamicCode/fluidThermoTemplate.C

@@ -23,13 +23,9 @@ License
 
 \*---------------------------------------------------------------------------*/
 
-#include "forThermo.H"
-#include "makeThermo.H"
-
+// Specie
 #include "${specie}.H"
 
-#include "thermo.H"
-
 // EoS
 #include "${equationOfState}.H"
 
@@ -47,6 +43,10 @@ License
 // Mixture
 #include "${mixture}.H"
 
+#include "thermo.H"
+#include "typedefThermo.H"
+#include "makeThermo.H"
+
 
 // * * * * * * * * * * * * * * * Global Functions  * * * * * * * * * * * * * //
 
@@ -75,17 +75,36 @@ extern "C"
 
 namespace Foam
 {
-    forThermo
+    typedefThermo
     (
         ${transport}Transport,
         ${energy},
         ${thermo}Thermo,
         ${equationOfState},
-        ${specie},
-        makeThermo,
+        ${specie}
+    );
+
+    defineThermo
+    (
         ${typeBase},
-        ${type},
-        ${mixture}
+        ${mixture},
+        ${transport}Transport${energy}${thermo}Thermo${equationOfState}${specie}
+    );
+
+    addThermo
+    (
+        fluidThermo,
+        ${typeBase},
+        ${mixture},
+        ${transport}Transport${energy}${thermo}Thermo${equationOfState}${specie}
+    );
+
+    addThermo
+    (
+        ${typeBase},
+        ${typeBase},
+        ${mixture},
+        ${transport}Transport${energy}${thermo}Thermo${equationOfState}${specie}
     );
 }
 

etc/codeTemplates/dynamicCode/psiuMulticomponentThermoTemplate.C

@@ -23,13 +23,9 @@ License
 
 \*---------------------------------------------------------------------------*/
 
-#include "forThermo.H"
-#include "makeMulticomponentThermo.H"
-
+// Specie
 #include "${specie}.H"
 
-#include "thermo.H"
-
 // EoS
 #include "${equationOfState}.H"
 
@@ -47,6 +43,10 @@ License
 // Mixture
 #include "${mixture}.H"
 
+#include "thermo.H"
+#include "typedefThermo.H"
+#include "makeThermo.H"
+
 
 // * * * * * * * * * * * * * * * Global Functions  * * * * * * * * * * * * * //
 
@@ -75,15 +75,28 @@ extern "C"
 
 namespace Foam
 {
-    forThermo
+    typedefThermo
     (
         ${transport}Transport,
         ${energy},
         ${thermo}Thermo,
         ${equationOfState},
-        ${specie},
-        makePsiuMulticomponentThermo,
-        ${mixture}
+        ${specie}
+    );
+
+    defineThermo
+    (
+        ${typeBase},
+        ${mixture},
+        ${transport}Transport${energy}${thermo}Thermo${equationOfState}${specie}
+    );
+
+    addThermo
+    (
+        ${typeBase},
+        ${typeBase},
+        ${mixture},
+        ${transport}Transport${energy}${thermo}Thermo${equationOfState}${specie}
     );
 }
 

etc/codeTemplates/dynamicCode/solidThermoTemplate.C

@@ -23,13 +23,9 @@ License
 
 \*---------------------------------------------------------------------------*/
 
-#include "forThermo.H"
-#include "makeSolidThermo.H"
-
+// Specie
 #include "${specie}.H"
 
-#include "thermo.H"
-
 // EoS
 #include "${equationOfState}.H"
 
@@ -47,6 +43,10 @@ License
 // Mixture
 #include "${mixture}.H"
 
+#include "thermo.H"
+#include "typedefThermo.H"
+#include "makeThermo.H"
+
 
 // * * * * * * * * * * * * * * * Global Functions  * * * * * * * * * * * * * //
 
@@ -75,17 +75,28 @@ extern "C"
 
 namespace Foam
 {
-    forThermo
+    typedefThermo
     (
         ${transport}Transport,
         ${energy},
         ${thermo}Thermo,
         ${equationOfState},
-        ${specie},
-        makeSolidThermo,
+        ${specie}
+    );
+
+    defineThermo
+    (
         ${typeBase},
-        ${type},
-        ${mixture}
+        ${mixture},
+        ${transport}Transport${energy}${thermo}Thermo${equationOfState}${specie}
+    );
+
+    addThermo
+    (
+        ${typeBase},
+        ${typeBase},
+        ${mixture},
+        ${transport}Transport${energy}${thermo}Thermo${equationOfState}${specie}
     );
 }
 

src/ThermophysicalTransportModels/fluid/laminar/Fickian/Fickian.C

@@ -41,8 +41,7 @@ namespace Foam
 template<class BasicThermophysicalTransportModel>
 void Fickian<BasicThermophysicalTransportModel>::updateDm() const
 {
-    const basicSpecieMixture& composition = this->thermo().composition();
-    const PtrList<volScalarField>& Y = composition.Y();
+    const PtrList<volScalarField>& Y = this->thermo().Y();
     const volScalarField& p = this->thermo().p();
     const volScalarField& T = this->thermo().T();
 
@@ -83,7 +82,7 @@ void Fickian<BasicThermophysicalTransportModel>::updateDm() const
                            (
                                "Wj",
                                Wm.dimensions(),
-                               composition.Wi(j)
+                               this->thermo().Wi(j)
                            )
                           *(
                                i < j
@@ -100,7 +99,12 @@ void Fickian<BasicThermophysicalTransportModel>::updateDm() const
                 (
                     1/Wm
                   - Y[i]
-                   /dimensionedScalar("Wi", Wm.dimensions(), composition.Wi(i))
+                   /dimensionedScalar
+                    (
+                        "Wi",
+                        Wm.dimensions(),
+                        this->thermo().Wi(i)
+                    )
                 )/max(sumXbyD, dimensionedScalar(sumXbyD.dimensions(), small))
             );
         }
@@ -142,14 +146,14 @@ Fickian<BasicThermophysicalTransportModel>::Fickian
 
     mixtureDiffusionCoefficients_(true),
 
-    DFuncs_(this->thermo().composition().species().size()),
+    DFuncs_(this->thermo().species().size()),
 
-    DmFuncs_(this->thermo().composition().species().size()),
+    DmFuncs_(this->thermo().species().size()),
 
     DTFuncs_
     (
         this->coeffDict_.found("DT")
-      ? this->thermo().composition().species().size()
+      ? this->thermo().species().size()
       : 0
     )
 {}
@@ -165,8 +169,7 @@ bool Fickian<BasicThermophysicalTransportModel>::read()
         BasicThermophysicalTransportModel::read()
     )
     {
-        const basicSpecieMixture& composition = this->thermo().composition();
-        const speciesTable& species = composition.species();
+        const speciesTable& species = this->thermo().species();
 
         this->coeffDict_.lookup("mixtureDiffusionCoefficients")
             >> mixtureDiffusionCoefficients_;
@@ -273,13 +276,11 @@ tmp<volScalarField> Fickian<BasicThermophysicalTransportModel>::DEff
     const volScalarField& Yi
 ) const
 {
-    const basicSpecieMixture& composition = this->thermo().composition();
-
     return volScalarField::New
     (
         "DEff",
         this->momentumTransport().rho()
-       *Dm()[composition.index(Yi)]
+       *Dm()[this->thermo().specieIndex(Yi)]
     );
 }
 
@@ -291,11 +292,9 @@ tmp<scalarField> Fickian<BasicThermophysicalTransportModel>::DEff
     const label patchi
 ) const
 {
-    const basicSpecieMixture& composition = this->thermo().composition();
-
     return
         this->momentumTransport().rho().boundaryField()[patchi]
-       *Dm()[composition.index(Yi)].boundaryField()[patchi];
+       *Dm()[this->thermo().specieIndex(Yi)].boundaryField()[patchi];
 }
 
 
@@ -316,8 +315,9 @@ tmp<surfaceScalarField> Fickian<BasicThermophysicalTransportModel>::q() const
         )
     );
 
-    const basicSpecieMixture& composition = this->thermo().composition();
-    const PtrList<volScalarField>& Y = composition.Y();
+    const PtrList<volScalarField>& Y = this->thermo().Y();
+    const volScalarField& p = this->thermo().p();
+    const volScalarField& T = this->thermo().T();
 
     if (Y.size())
     {
@@ -343,12 +343,9 @@ tmp<surfaceScalarField> Fickian<BasicThermophysicalTransportModel>::q() const
 
         forAll(Y, i)
         {
-            if (i != composition.defaultSpecie())
+            if (i != this->thermo().defaultSpecie())
             {
-                const volScalarField hi
-                (
-                    composition.Hs(i, this->thermo().p(), this->thermo().T())
-                );
+                const volScalarField hi(this->thermo().hsi(i, p, T));
 
                 const surfaceScalarField ji(this->j(Y[i]));
                 sumJ += ji;
@@ -358,12 +355,9 @@ tmp<surfaceScalarField> Fickian<BasicThermophysicalTransportModel>::q() const
         }
 
         {
-            const label i = composition.defaultSpecie();
+            const label i = this->thermo().defaultSpecie();
 
-            const volScalarField hi
-            (
-                composition.Hs(i, this->thermo().p(), this->thermo().T())
-            );
+            const volScalarField hi(this->thermo().hsi(i, p, T));
 
             sumJh -= sumJ*fvc::interpolate(hi);
         }
@@ -390,8 +384,9 @@ tmp<fvScalarMatrix> Fickian<BasicThermophysicalTransportModel>::divq
         )
     );
 
-    const basicSpecieMixture& composition = this->thermo().composition();
-    const PtrList<volScalarField>& Y = composition.Y();
+    const PtrList<volScalarField>& Y = this->thermo().Y();
+    const volScalarField& p = this->thermo().p();
+    const volScalarField& T = this->thermo().T();
 
     tmpDivq.ref() -=
         fvm::laplacianCorrection(this->alpha()*this->alphaEff(), he);
@@ -418,12 +413,9 @@ tmp<fvScalarMatrix> Fickian<BasicThermophysicalTransportModel>::divq
 
     forAll(Y, i)
     {
-        if (i != composition.defaultSpecie())
+        if (i != this->thermo().defaultSpecie())
         {
-            const volScalarField hi
-            (
-                composition.Hs(i, this->thermo().p(), this->thermo().T())
-            );
+            const volScalarField hi(this->thermo().hsi(i, p, T));
 
             const surfaceScalarField ji(this->j(Y[i]));
             sumJ += ji;
@@ -433,12 +425,9 @@ tmp<fvScalarMatrix> Fickian<BasicThermophysicalTransportModel>::divq
     }
 
     {
-        const label i = composition.defaultSpecie();
+        const label i = this->thermo().defaultSpecie();
 
-        const volScalarField hi
-        (
-            composition.Hs(i, this->thermo().p(), this->thermo().T())
-        );
+        const volScalarField hi(this->thermo().hsi(i, p, T));
 
         sumJh -= sumJ*fvc::interpolate(hi);
     }
@@ -457,7 +446,6 @@ tmp<surfaceScalarField> Fickian<BasicThermophysicalTransportModel>::j
 {
     if (DTFuncs_.size())
     {
-        const basicSpecieMixture& composition = this->thermo().composition();
         const volScalarField& p = this->thermo().p();
         const volScalarField& T = this->thermo().T();
 
@@ -467,7 +455,7 @@ tmp<surfaceScalarField> Fickian<BasicThermophysicalTransportModel>::j
             (
                 evaluate
                 (
-                    DTFuncs_[composition.index(Yi)],
+                    DTFuncs_[this->thermo().specieIndex(Yi)],
                     dimDynamicViscosity,
                     p,
                     T
@@ -490,7 +478,6 @@ tmp<fvScalarMatrix> Fickian<BasicThermophysicalTransportModel>::divj
 {
     if (DTFuncs_.size())
     {
-        const basicSpecieMixture& composition = this->thermo().composition();
         const volScalarField& p = this->thermo().p();
         const volScalarField& T = this->thermo().T();
 
@@ -502,7 +489,7 @@ tmp<fvScalarMatrix> Fickian<BasicThermophysicalTransportModel>::divj
                 (
                     evaluate
                     (
-                        DTFuncs_[composition.index(Yi)],
+                        DTFuncs_[this->thermo().specieIndex(Yi)],
                         dimDynamicViscosity,
                         p,
                         T

src/ThermophysicalTransportModels/fluid/laminar/MaxwellStefan/MaxwellStefan.C

@@ -42,8 +42,7 @@ template<class BasicThermophysicalTransportModel>
 void MaxwellStefan<BasicThermophysicalTransportModel>::
 transformDiffusionCoefficient() const
 {
-    const basicSpecieMixture& composition = this->thermo().composition();
-    const label d = composition.defaultSpecie();
+    const label d = this->thermo().defaultSpecie();
 
     // Calculate the molecular weight of the mixture and the mole fractions
     scalar Wm = 0;
@@ -110,8 +109,7 @@ template<class BasicThermophysicalTransportModel>
 void MaxwellStefan<BasicThermophysicalTransportModel>::
 transformDiffusionCoefficientFields() const
 {
-    const basicSpecieMixture& composition = this->thermo().composition();
-    const label d = composition.defaultSpecie();
+    const label d = this->thermo().defaultSpecie();
 
     // For each cell or patch face
     forAll(*(YPtrs[0]), pi)
@@ -166,13 +164,12 @@ void MaxwellStefan<BasicThermophysicalTransportModel>::transform
     List<PtrList<volScalarField>>& Dij
 ) const
 {
-    const basicSpecieMixture& composition = this->thermo().composition();
-    const PtrList<volScalarField>& Y = composition.Y();
+    const PtrList<volScalarField>& Y = this->thermo().Y();
     const volScalarField& Y0 = Y[0];
 
     forAll(W, i)
     {
-        // Map composition.Y() internal fields -> YPtrs
+        // Map this->thermo().Y() internal fields -> YPtrs
         YPtrs[i] = &Y[i].primitiveField();
 
         // Map Dii_ internal fields -> DijPtrs
@@ -196,7 +193,7 @@ void MaxwellStefan<BasicThermophysicalTransportModel>::transform
     {
         forAll(W, i)
         {
-            // Map composition.Y() patch fields -> YPtrs
+            // Map this->thermo().Y() patch fields -> YPtrs
             YPtrs[i] = &Y[i].boundaryField()[patchi];
 
             // Map Dii_ patch fields -> DijPtrs
@@ -222,10 +219,9 @@ void MaxwellStefan<BasicThermophysicalTransportModel>::transform
 template<class BasicThermophysicalTransportModel>
 void MaxwellStefan<BasicThermophysicalTransportModel>::updateDii() const
 {
-    const basicSpecieMixture& composition = this->thermo().composition();
-    const label d = composition.defaultSpecie();
+    const label d = this->thermo().defaultSpecie();
 
-    const PtrList<volScalarField>& Y = composition.Y();
+    const PtrList<volScalarField>& Y = this->thermo().Y();
     const volScalarField& p = this->thermo().p();
     const volScalarField& T = this->thermo().T();
     const volScalarField& rho = this->momentumTransport().rho();
@@ -363,16 +359,16 @@ MaxwellStefan<BasicThermophysicalTransportModel>::MaxwellStefan
 
     UpdateableMeshObject(*this, thermo.mesh()),
 
-    DFuncs_(this->thermo().composition().species().size()),
+    DFuncs_(this->thermo().species().size()),
 
     DTFuncs_
     (
         this->coeffDict_.found("DT")
-      ? this->thermo().composition().species().size()
+      ? this->thermo().species().size()
       : 0
     ),
 
-    W(this->thermo().composition().species().size()),
+    W(this->thermo().species().size()),
 
     YPtrs(W.size()),
     DijPtrs(W.size()),
@@ -385,12 +381,10 @@ MaxwellStefan<BasicThermophysicalTransportModel>::MaxwellStefan
     invA(A.m()),
     D(W.size())
 {
-    const basicSpecieMixture& composition = this->thermo().composition();
-
     // Set the molecular weights of the species
     forAll(W, i)
     {
-        W[i] = composition.Wi(i);
+        W[i] = this->thermo().Wi(i);
     }
 }
 
@@ -405,8 +399,7 @@ bool MaxwellStefan<BasicThermophysicalTransportModel>::read()
         BasicThermophysicalTransportModel::read()
     )
     {
-        const basicSpecieMixture& composition = this->thermo().composition();
-        const speciesTable& species = composition.species();
+        const speciesTable& species = this->thermo().species();
 
         const dictionary& Ddict = this->coeffDict_.subDict("D");
 
@@ -492,12 +485,10 @@ tmp<volScalarField> MaxwellStefan<BasicThermophysicalTransportModel>::DEff
     const volScalarField& Yi
 ) const
 {
-    const basicSpecieMixture& composition = this->thermo().composition();
-
     return volScalarField::New
     (
         "DEff",
-        this->momentumTransport().rho()*Dii()[composition.index(Yi)]
+        this->momentumTransport().rho()*Dii()[this->thermo().specieIndex(Yi)]
     );
 }
 
@@ -509,11 +500,9 @@ tmp<scalarField> MaxwellStefan<BasicThermophysicalTransportModel>::DEff
     const label patchi
 ) const
 {
-    const basicSpecieMixture& composition = this->thermo().composition();
-
     return
         this->momentumTransport().rho().boundaryField()[patchi]
-       *Dii()[composition.index(Yi)].boundaryField()[patchi];
+       *Dii()[this->thermo().specieIndex(Yi)].boundaryField()[patchi];
 }
 
 
@@ -535,10 +524,11 @@ MaxwellStefan<BasicThermophysicalTransportModel>::q() const
         )
     );
 
-    const basicSpecieMixture& composition = this->thermo().composition();
-    const label d = composition.defaultSpecie();
+    const label d = this->thermo().defaultSpecie();
 
-    const PtrList<volScalarField>& Y = composition.Y();
+    const PtrList<volScalarField>& Y = this->thermo().Y();
+    const volScalarField& p = this->thermo().p();
+    const volScalarField& T = this->thermo().T();
 
     if (Y.size())
     {
@@ -566,10 +556,7 @@ MaxwellStefan<BasicThermophysicalTransportModel>::q() const
         {
             if (i != d)
             {
-                const volScalarField hi
-                (
-                    composition.Hs(i, this->thermo().p(), this->thermo().T())
-                );
+                const volScalarField hi(this->thermo().hsi(i, p, T));
 
                 const surfaceScalarField ji(this->j(Y[i]));
                 sumJ += ji;
@@ -581,10 +568,7 @@ MaxwellStefan<BasicThermophysicalTransportModel>::q() const
         {
             const label i = d;
 
-            const volScalarField hi
-            (
-                composition.Hs(i, this->thermo().p(), this->thermo().T())
-            );
+            const volScalarField hi(this->thermo().hsi(i, p, T));
 
             sumJh -= sumJ*fvc::interpolate(hi);
         }
@@ -611,10 +595,11 @@ tmp<fvScalarMatrix> MaxwellStefan<BasicThermophysicalTransportModel>::divq
         )
     );
 
-    const basicSpecieMixture& composition = this->thermo().composition();
-    const label d = composition.defaultSpecie();
+    const label d = this->thermo().defaultSpecie();
 
-    const PtrList<volScalarField>& Y = composition.Y();
+    const PtrList<volScalarField>& Y = this->thermo().Y();
+    const volScalarField& p = this->thermo().p();
+    const volScalarField& T = this->thermo().T();
 
     tmpDivq.ref() -=
         fvm::laplacianCorrection(this->alpha()*this->alphaEff(), he);
@@ -643,10 +628,7 @@ tmp<fvScalarMatrix> MaxwellStefan<BasicThermophysicalTransportModel>::divq
     {
         if (i != d)
         {
-            const volScalarField hi
-            (
-                composition.Hs(i, this->thermo().p(), this->thermo().T())
-            );
+            const volScalarField hi(this->thermo().hsi(i, p, T));
 
             const surfaceScalarField ji(this->j(Y[i]));
             sumJ += ji;
@@ -658,10 +640,7 @@ tmp<fvScalarMatrix> MaxwellStefan<BasicThermophysicalTransportModel>::divq
     {
         const label i = d;
 
-        const volScalarField hi
-        (
-            composition.Hs(i, this->thermo().p(), this->thermo().T())
-        );
+        const volScalarField hi(this->thermo().hsi(i, p, T));
 
         sumJh -= sumJ*fvc::interpolate(hi);
     }
@@ -678,12 +657,11 @@ tmp<surfaceScalarField> MaxwellStefan<BasicThermophysicalTransportModel>::j
     const volScalarField& Yi
 ) const
 {
-    const basicSpecieMixture& composition = this->thermo().composition();
-    const label d = composition.defaultSpecie();
+    const label d = this->thermo().defaultSpecie();
 
-    if (composition.index(Yi) == d)
+    if (this->thermo().specieIndex(Yi) == d)
     {
-        const PtrList<volScalarField>& Y = composition.Y();
+        const PtrList<volScalarField>& Y = this->thermo().Y();
 
         tmp<surfaceScalarField> tjd
         (
@@ -715,7 +693,7 @@ tmp<surfaceScalarField> MaxwellStefan<BasicThermophysicalTransportModel>::j
     {
         return
             BasicThermophysicalTransportModel::j(Yi)
-          + jexp()[composition.index(Yi)];
+          + jexp()[this->thermo().specieIndex(Yi)];
     }
 }
 
@@ -726,10 +704,9 @@ tmp<fvScalarMatrix> MaxwellStefan<BasicThermophysicalTransportModel>::divj
     volScalarField& Yi
 ) const
 {
-    const basicSpecieMixture& composition = this->thermo().composition();
     return
         BasicThermophysicalTransportModel::divj(Yi)
-      + fvc::div(jexp()[composition.index(Yi)]*Yi.mesh().magSf());
+      + fvc::div(jexp()[this->thermo().specieIndex(Yi)]*Yi.mesh().magSf());
 }
 
 

src/ThermophysicalTransportModels/fluid/turbulence/nonUnityLewisEddyDiffusivity/nonUnityLewisEddyDiffusivity.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2020-2022 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2020-2023 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -102,8 +102,7 @@ nonUnityLewisEddyDiffusivity<TurbulenceThermophysicalTransportModel>::q() const
         )
     );
 
-    const basicSpecieMixture& composition = this->thermo().composition();
-    const PtrList<volScalarField>& Y = composition.Y();
+    const PtrList<volScalarField>& Y = this->thermo().Y();
 
     if (Y.size())
     {
@@ -121,7 +120,7 @@ nonUnityLewisEddyDiffusivity<TurbulenceThermophysicalTransportModel>::q() const
         {
             const volScalarField hi
             (
-                composition.Hs(i, this->thermo().p(), this->thermo().T())
+                this->thermo().hsi(i, this->thermo().p(), this->thermo().T())
             );
 
             hGradY += fvc::interpolate(hi)*fvc::snGrad(Y[i]);
@@ -158,8 +157,7 @@ nonUnityLewisEddyDiffusivity<TurbulenceThermophysicalTransportModel>::divq
         )
     );
 
-    const basicSpecieMixture& composition = this->thermo().composition();
-    const PtrList<volScalarField>& Y = composition.Y();
+    const PtrList<volScalarField>& Y = this->thermo().Y();
 
     tmpDivq.ref() -=
         fvm::laplacianCorrection(this->alpha()*this->alphaEff(), he);
@@ -178,7 +176,7 @@ nonUnityLewisEddyDiffusivity<TurbulenceThermophysicalTransportModel>::divq
     {
         const volScalarField hi
         (
-            composition.Hs(i, this->thermo().p(), this->thermo().T())
+            this->thermo().hsi(i, this->thermo().p(), this->thermo().T())
         );
 
         hGradY += fvc::interpolate(hi)*fvc::snGrad(Y[i]);

src/combustionModels/EDC/EDC.C

@@ -211,7 +211,7 @@ Foam::combustionModels::EDC::R(volScalarField& Y) const
     tmp<fvScalarMatrix> tSu(new fvScalarMatrix(Y, dimMass/dimTime));
     fvScalarMatrix& Su = tSu.ref();
 
-    const label speciei = this->thermo().composition().species()[Y.member()];
+    const label speciei = this->thermo().species()[Y.member()];
     Su += chemistryPtr_->RR()[speciei];
 
     return kappa_*tSu;

src/combustionModels/FSD/FSD.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2011-2022 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2023 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -107,9 +107,9 @@ void Foam::combustionModels::FSD::calculateSourceNorm()
 
     const label fuelI = this->fuelIndex();
 
-    const volScalarField& YFuel = this->thermo().composition().Y()[fuelI];
+    const volScalarField& YFuel = this->thermo().Y()[fuelI];
 
-    const volScalarField& YO2 = this->thermo().composition().Y("O2");
+    const volScalarField& YO2 = this->thermo().Y("O2");
 
     const dimensionedScalar s = this->s();
 
@@ -298,7 +298,7 @@ void Foam::combustionModels::FSD::calculateSourceNorm()
     forAll(productsIndex, j)
     {
         label specieI = productsIndex[j];
-        const volScalarField& Yp = this->thermo().composition().Y()[specieI];
+        const volScalarField& Yp = this->thermo().Y()[specieI];
         products += Yp;
     }
 

src/combustionModels/diffusion/diffusion.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2012-2022 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2012-2023 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -79,12 +79,11 @@ void Foam::combustionModels::diffusion::correct()
 
     const label fuelI = this->fuelIndex();
 
-    const volScalarField& YFuel = this->thermo().composition().Y()[fuelI];
+    const volScalarField& YFuel = this->thermo().Y()[fuelI];
 
-    if (this->thermo().composition().contains(oxidantName_))
+    if (this->thermo().containsSpecie(oxidantName_))
     {
-        const volScalarField& YO2 =
-            this->thermo().composition().Y(oxidantName_);
+        const volScalarField& YO2 = this->thermo().Y(oxidantName_);
 
         this->wFuel_ ==
             C_*this->thermo().rho()*this->turbulence().nuEff()

src/combustionModels/functionObjects/adjustTimeStepToCombustion/adjustTimeStepToCombustion.C

@@ -176,9 +176,9 @@ Foam::functionObjects::adjustTimeStepToCombustion::maxDeltaT() const
             dimensionedScalar(dimless/dimTime, 0)
         )
     );
-    forAll(thermo.composition().Y(), i)
+    forAll(thermo.Y(), i)
     {
-        if (thermo.composition().solve(i))
+        if (thermo.solveSpecie(i))
         {
             rhoDotByRho += mag(combustion.R(i))/2/thermo.rho()();
         }

src/combustionModels/infinitelyFastChemistry/infinitelyFastChemistry.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2011-2022 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2023 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -82,13 +82,13 @@ void Foam::combustionModels::infinitelyFastChemistry::correct()
 
     const label fuelI = this->fuelIndex();
 
-    const volScalarField& YFuel = this->thermo().composition().Y()[fuelI];
+    const volScalarField& YFuel = this->thermo().Y()[fuelI];
 
     const dimensionedScalar s = this->s();
 
-    if (this->thermo().composition().contains("O2"))
+    if (this->thermo().containsSpecie("O2"))
     {
-        const volScalarField& YO2 = this->thermo().composition().Y("O2");
+        const volScalarField& YO2 = this->thermo().Y("O2");
 
         this->wFuel_ ==
             this->rho()/(this->mesh().time().deltaT()*C_)

src/combustionModels/laminar/laminar.C

@@ -140,7 +140,7 @@ Foam::combustionModels::laminar::R(volScalarField& Y) const
     tmp<fvScalarMatrix> tSu(new fvScalarMatrix(Y, dimMass/dimTime));
     fvScalarMatrix& Su = tSu.ref();
 
-    const label specieI = this->thermo().composition().species()[Y.member()];
+    const label specieI = this->thermo().species()[Y.member()];
     Su += chemistryPtr_->RR()[specieI];
 
     return tSu;

src/combustionModels/noCombustion/noCombustion.C

@@ -71,7 +71,7 @@ Foam::combustionModels::noCombustion::R(const label speciei) const
     return
         volScalarField::Internal::New
         (
-            typedName("R_" + this->thermo().composition().Y()[speciei].name()),
+            typedName("R_" + this->thermo().Y()[speciei].name()),
             this->mesh(),
             dimensionedScalar(dimDensity/dimTime, 0)
         );

src/combustionModels/radiationModels/sootModels/mixtureFraction/mixtureFraction.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2013-2022 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2013-2023 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -87,8 +87,6 @@ Foam::radiationModels::sootModels::mixtureFraction::mixtureFraction
             combustionModel::combustionPropertiesName
         );
 
-    const basicSpecieMixture& mixture = combustion.mixture();
-
     const reaction& singleReaction = combustion.singleReaction();
 
     scalar totalMol = 0;
@@ -108,15 +106,15 @@ Foam::radiationModels::sootModels::mixtureFraction::mixtureFraction
         const label speciei = singleReaction.rhs()[i].index;
         const scalar stoichCoeff = singleReaction.rhs()[i].stoichCoeff;
         Xi[i] = mag(stoichCoeff)/totalMol;
-        Wm += Xi[i]*mixture.Wi(speciei);
+        Wm += Xi[i]*combustion.thermo().Wi(speciei);
     }
 
-    scalarList Yprod0(mixture.species().size(), 0.0);
+    scalarList Yprod0(combustion.thermo().species().size(), 0.0);
 
     forAll(singleReaction.rhs(), i)
     {
         const label speciei = singleReaction.rhs()[i].index;
-        Yprod0[speciei] = mixture.Wi(speciei)/Wm*Xi[i];
+        Yprod0[speciei] = combustion.thermo().Wi(speciei)/Wm*Xi[i];
     }
 
     const scalar XSoot = nuSoot_/totalMol;
@@ -129,10 +127,11 @@ Foam::radiationModels::sootModels::mixtureFraction::mixtureFraction
     if (mappingFieldName_ == "none")
     {
         const label index = singleReaction.rhs()[0].index;
-        mappingFieldName_ = mixture.Y(index).name();
+        mappingFieldName_ = combustion.thermo().Y(index).name();
     }
 
-    const label mapFieldIndex = mixture.species()[mappingFieldName_];
+    const label mapFieldIndex =
+        combustion.thermo().species()[mappingFieldName_];
 
     mapFieldMax_ = Yprod0[mapFieldIndex];
 }

src/combustionModels/singleStepCombustion/singleStepCombustion.C

@@ -32,7 +32,7 @@ License
 
 void Foam::combustionModels::singleStepCombustion::calculateqFuel()
 {
-    const scalar Wu = mixture_.Wi(fuelIndex_);
+    const scalar Wu = thermo_.Wi(fuelIndex_);
 
     forAll(reaction_.lhs(), i)
     {
@@ -40,7 +40,7 @@ void Foam::combustionModels::singleStepCombustion::calculateqFuel()
         const scalar stoichCoeff = reaction_.lhs()[i].stoichCoeff;
         specieStoichCoeffs_[speciei] = -stoichCoeff;
         qFuel_.value() +=
-            mixture_.Hf(speciei)*mixture_.Wi(speciei)*stoichCoeff/Wu;
+            thermo_.hfi(speciei)*thermo_.Wi(speciei)*stoichCoeff/Wu;
     }
 
     forAll(reaction_.rhs(), i)
@@ -49,7 +49,7 @@ void Foam::combustionModels::singleStepCombustion::calculateqFuel()
         const scalar stoichCoeff = reaction_.rhs()[i].stoichCoeff;
         specieStoichCoeffs_[speciei] = stoichCoeff;
         qFuel_.value() -=
-            mixture_.Hf(speciei)*mixture_.Wi(speciei)*stoichCoeff/Wu;
+            thermo_.hfi(speciei)*thermo_.Wi(speciei)*stoichCoeff/Wu;
         specieProd_[speciei] = -1;
     }
 
@@ -59,17 +59,17 @@ void Foam::combustionModels::singleStepCombustion::calculateqFuel()
 
 void Foam::combustionModels::singleStepCombustion::massAndAirStoichRatios()
 {
-    const label O2Index = mixture_.species()["O2"];
-    const scalar Wu = mixture_.Wi(fuelIndex_);
+    const label O2Index = thermo_.species()["O2"];
+    const scalar Wu = thermo_.Wi(fuelIndex_);
 
     stoicRatio_ =
         (
-            mixture_.Wi(mixture_.defaultSpecie())
-           *specieStoichCoeffs_[mixture_.defaultSpecie()]
-          + mixture_.Wi(O2Index)*mag(specieStoichCoeffs_[O2Index])
+            thermo_.Wi(thermo_.defaultSpecie())
+           *specieStoichCoeffs_[thermo_.defaultSpecie()]
+          + thermo_.Wi(O2Index)*mag(specieStoichCoeffs_[O2Index])
         )/(Wu*mag(specieStoichCoeffs_[fuelIndex_]));
 
-    s_ = mixture_.Wi(O2Index)*mag(specieStoichCoeffs_[O2Index])
+    s_ = thermo_.Wi(O2Index)*mag(specieStoichCoeffs_[O2Index])
         /(Wu*mag(specieStoichCoeffs_[fuelIndex_]));
 
     Info << "stoichiometric air-fuel ratio: " << stoicRatio_.value() << endl;
@@ -93,13 +93,13 @@ void Foam::combustionModels::singleStepCombustion::calculateMaxProducts()
     {
         const label speciei = reaction_.rhs()[i].index;
         Xi[i] = mag(specieStoichCoeffs_[speciei])/totalMol;
-        Wm += Xi[i]*mixture_.Wi(speciei);
+        Wm += Xi[i]*thermo_.Wi(speciei);
     }
 
     forAll(reaction_.rhs(), i)
     {
         const label speciei = reaction_.rhs()[i].index;
-        Yprod0_[speciei] =  mixture_.Wi(speciei)/Wm*Xi[i];
+        Yprod0_[speciei] =  thermo_.Wi(speciei)/Wm*Xi[i];
     }
 
     Info << "Maximum products mass concentrations: " << nl;
@@ -107,7 +107,7 @@ void Foam::combustionModels::singleStepCombustion::calculateMaxProducts()
     {
         if (Yprod0_[i] > 0)
         {
-            Info<< "    " << mixture_.species()[i] << ": " << Yprod0_[i] << nl;
+            Info<< "    " << thermo_.species()[i] << ": " << Yprod0_[i] << nl;
         }
     }
 
@@ -115,8 +115,8 @@ void Foam::combustionModels::singleStepCombustion::calculateMaxProducts()
     forAll(specieStoichCoeffs_, i)
     {
         specieStoichCoeffs_[i] =
-            specieStoichCoeffs_[i]*mixture_.Wi(i)
-           /(mixture_.Wi(fuelIndex_)*mag(specieStoichCoeffs_[fuelIndex_]));
+            specieStoichCoeffs_[i]*thermo_.Wi(i)
+           /(thermo_.Wi(fuelIndex_)*mag(specieStoichCoeffs_[fuelIndex_]));
     }
 }
 
@@ -132,15 +132,14 @@ Foam::combustionModels::singleStepCombustion::singleStepCombustion
 )
 :
     combustionModel(modelType, thermo, turb, combustionProperties),
-    mixture_(dynamic_cast<const basicSpecieMixture&>(this->thermo())),
-    reaction_(mixture_.species(), this->subDict("reaction")),
+    reaction_(thermo_.species(), this->subDict("reaction")),
     stoicRatio_(dimensionedScalar("stoicRatio", dimless, 0)),
     s_(dimensionedScalar("s", dimless, 0)),
     qFuel_(dimensionedScalar("qFuel", sqr(dimVelocity), 0)),
-    specieStoichCoeffs_(mixture_.species().size(), 0.0),
-    Yprod0_(mixture_.species().size(), 0.0),
+    specieStoichCoeffs_(thermo_.species().size(), 0.0),
+    Yprod0_(thermo_.species().size(), 0.0),
     fres_(Yprod0_.size()),
-    fuelIndex_(mixture_.species()[thermo.properties().lookup("fuel")]),
+    fuelIndex_(thermo_.species()[thermo.properties().lookup("fuel")]),
     specieProd_(Yprod0_.size(), 1),
     wFuel_
     (
@@ -161,7 +160,7 @@ Foam::combustionModels::singleStepCombustion::singleStepCombustion
     {
         IOobject header
         (
-            "fres_" + mixture_.species()[fresI],
+            "fres_" + thermo_.species()[fresI],
             this->mesh().time().name(),
             this->mesh()
         );
@@ -213,7 +212,7 @@ Foam::combustionModels::singleStepCombustion::R(const label speciei) const
 Foam::tmp<Foam::fvScalarMatrix>
 Foam::combustionModels::singleStepCombustion::R(volScalarField& Y) const
 {
-    const label specieI = mixture_.species()[Y.member()];
+    const label specieI = thermo_.species()[Y.member()];
 
     volScalarField wSpecie
     (
@@ -240,7 +239,7 @@ Foam::combustionModels::singleStepCombustion::Qdot() const
 {
     const label fuelI = fuelIndex();
     volScalarField& YFuel =
-        const_cast<volScalarField&>(this->thermo().composition().Y(fuelI));
+        const_cast<volScalarField&>(this->thermo().Y(fuelI));
 
     return -qFuel()*(R(YFuel) & YFuel);
 }
@@ -261,9 +260,9 @@ bool Foam::combustionModels::singleStepCombustion::read()
 
 void Foam::combustionModels::singleStepCombustion::fresCorrect()
 {
-    const label O2Index = mixture_.species()["O2"];
-    const volScalarField& YFuel = mixture_.Y()[fuelIndex_];
-    const volScalarField& YO2 = mixture_.Y()[O2Index];
+    const label O2Index = thermo_.species()["O2"];
+    const volScalarField& YFuel = thermo_.Y()[fuelIndex_];
+    const volScalarField& YO2 = thermo_.Y()[O2Index];
 
     // reactants
     forAll(reaction_.lhs(), i)
@@ -283,7 +282,7 @@ void Foam::combustionModels::singleStepCombustion::fresCorrect()
     forAll(reaction_.rhs(), i)
     {
         const label speciei = reaction_.rhs()[i].index;
-        if (speciei != mixture_.defaultSpecie())
+        if (speciei != thermo_.defaultSpecie())
         {
             forAll(fres_[speciei], celli)
             {

src/combustionModels/singleStepCombustion/singleStepCombustion.H

@@ -25,7 +25,7 @@ Class
     Foam::combustionModels::singleStepCombustion
 
 Description
-    Base class for combustion models using basicSpecieMixture.
+    Base class for single-step combustion models
 
 SourceFiles
     singleStepCombustion.C
@@ -36,7 +36,6 @@ SourceFiles
 #define singleStepCombustion_H
 
 #include "combustionModel.H"
-#include "basicSpecieMixture.H"
 #include "reaction.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -58,9 +57,6 @@ protected:
 
     // Protected data
 
-        //- Reference to the mixture
-        const basicSpecieMixture& mixture_;
-
         //- The single-step reaction
         reaction reaction_;
 
@@ -132,9 +128,6 @@ public:
 
         // Access functions
 
-            //- Return the mixture
-            const basicSpecieMixture& mixture() const;
-
             //- Return the single step reaction
             inline const reaction& singleReaction() const;
 

src/combustionModels/singleStepCombustion/singleStepCombustionI.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2011-2020 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2023 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -34,12 +34,6 @@ namespace combustionModels
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-inline const Foam::basicSpecieMixture& singleStepCombustion::mixture() const
-{
-    return mixture_;
-}
-
-
 inline const Foam::reaction& singleStepCombustion::singleReaction() const
 {
     return reaction_;

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

@@ -24,7 +24,7 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "clouds.H"
-#include "basicSpecieMixture.H"
+#include "multicomponentThermo.H"
 #include "fvMatrix.H"
 #include "addToRunTimeSelectionTable.H"
 
@@ -167,16 +167,16 @@ Foam::wordList Foam::fv::clouds::addSupFields() const
 
         fieldNames.append(carrierThermo.he().name());
 
-        if (isA<basicSpecieMixture>(carrierThermo))
+        if (isA<multicomponentThermo>(carrierThermo))
         {
-            const basicSpecieMixture& composition =
-                refCast<const basicSpecieMixture>(carrierThermo);
+            const multicomponentThermo& carrierMcThermo =
+                refCast<const multicomponentThermo>(carrierThermo);
 
-            const PtrList<volScalarField>& Y = composition.Y();
+            const PtrList<volScalarField>& Y = carrierMcThermo.Y();
 
             forAll(Y, i)
             {
-                if (composition.solve(i))
+                if (carrierMcThermo.solveSpecie(i))
                 {
                     fieldNames.append(Y[i].name());
                 }
@@ -262,20 +262,20 @@ void Foam::fv::clouds::addSup
     {
         eqn += cloudsPtr_().Sh(eqn.psi());
     }
-    else if
-    (
-        isA<basicSpecieMixture>(carrierThermo)
-     && refCast<const basicSpecieMixture>(carrierThermo).contains
-        (
-            eqn.psi().name()
-        )
-    )
+    else if (isA<multicomponentThermo>(carrierThermo))
     {
-        eqn += cloudsPtr_().SYi
-        (
-            refCast<const basicSpecieMixture>(carrierThermo).index(eqn.psi()),
-            eqn.psi()
-        );
+        const multicomponentThermo& carrierMcThermo =
+            refCast<const multicomponentThermo>(carrierThermo);
+
+        if (carrierMcThermo.containsSpecie(eqn.psi().name()))
+        {
+            eqn +=
+                cloudsPtr_().SYi
+                (
+                    carrierMcThermo.specieIndex(eqn.psi()),
+                    eqn.psi()
+                );
+        }
     }
 }
 

src/lagrangian/parcel/parcels/Templates/ReactingMultiphaseParcel/ReactingMultiphaseParcel.C

@@ -435,7 +435,7 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::calc
         {
             scalar dm = np0*dMassGas[i];
             label gid = composition.localToCarrierId(GAS, i);
-            scalar hs = composition.carrier().Hs(gid, pc, T0);
+            scalar hs = composition.carrier().hsi(gid, pc, T0);
             cloud.rhoTrans(gid)[this->cell()] += dm;
             cloud.UTransRef()[this->cell()] += dm*U0;
             cloud.hsTransRef()[this->cell()] += dm*hs;
@@ -444,7 +444,7 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::calc
         {
             scalar dm = np0*dMassLiquid[i];
             label gid = composition.localToCarrierId(LIQ, i);
-            scalar hs = composition.carrier().Hs(gid, pc, T0);
+            scalar hs = composition.carrier().hsi(gid, pc, T0);
             cloud.rhoTrans(gid)[this->cell()] += dm;
             cloud.UTransRef()[this->cell()] += dm*U0;
             cloud.hsTransRef()[this->cell()] += dm*hs;
@@ -456,7 +456,7 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::calc
         {
             scalar dm = np0*dMassSolid[i];
             label gid = composition.localToCarrierId(SLD, i);
-            scalar hs = composition.carrier().Hs(gid, pc, T0);
+            scalar hs = composition.carrier().hsi(gid, pc, T0);
             cloud.rhoTrans(gid)[this->cell()] += dm;
             cloud.UTransRef()[this->cell()] += dm*U0;
             cloud.hsTransRef()[this->cell()] += dm*hs;
@@ -466,7 +466,7 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::calc
         forAll(dMassSRCarrier, i)
         {
             scalar dm = np0*dMassSRCarrier[i];
-            scalar hs = composition.carrier().Hs(i, pc, T0);
+            scalar hs = composition.carrier().hsi(i, pc, T0);
             cloud.rhoTrans(i)[this->cell()] += dm;
             cloud.UTransRef()[this->cell()] += dm*U0;
             cloud.hsTransRef()[this->cell()] += dm*hs;
@@ -598,7 +598,7 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::calcDevolatilisation
         forAll(dMassDV, i)
         {
             const label id = composition.localToCarrierId(GAS, i);
-            const scalar Cp = composition.carrier().Cp(id, td.pc(), Ts);
+            const scalar Cp = composition.carrier().Cpi(id, td.pc(), Ts);
             const scalar W = composition.carrier().Wi(id);
             const scalar Ni = dMassDV[i]/(this->areaS(d)*dt*W);
 

src/lagrangian/parcel/parcels/Templates/ReactingParcel/ReactingParcel.C

@@ -132,7 +132,7 @@ void Foam::ReactingParcel<ParcelType>::calcPhaseChange
         {
             const label cid = composition.localToCarrierId(idPhase, i);
 
-            const scalar Cp = composition.carrier().Cp(cid, td.pc(), Tsdash);
+            const scalar Cp = composition.carrier().Cpi(cid, td.pc(), Tsdash);
             const scalar W = composition.carrier().Wi(cid);
             const scalar Ni = dMassPC[i]/(this->areaS(d)*dt*W);
 
@@ -254,7 +254,7 @@ void Foam::ReactingParcel<ParcelType>::cellValueSourceCorrection
     forAll(cloud.rhoTrans(), i)
     {
         scalar Y = cloud.rhoTrans(i)[this->cell()]/addedMass;
-        CpEff += Y*cloud.composition().carrier().Cp(i, td.pc(), td.Tc());
+        CpEff += Y*cloud.composition().carrier().Cpi(i, td.pc(), td.Tc());
     }
 
     const scalar Cpc = td.CpInterp().psi()[this->cell()];
@@ -296,14 +296,15 @@ void Foam::ReactingParcel<ParcelType>::correctSurfaceValues
         return;
     }
 
-    const basicSpecieMixture& carrier = cloud.composition().carrier();
+    const fluidMulticomponentThermo& carrierThermo =
+        cloud.composition().carrier();
 
     // Far field carrier  molar fractions
-    scalarField Xinf(carrier.species().size());
+    scalarField Xinf(carrierThermo.species().size());
 
     forAll(Xinf, i)
     {
-        Xinf[i] = carrier.Y(i)[this->cell()]/carrier.Wi(i);
+        Xinf[i] = carrierThermo.Y(i)[this->cell()]/carrierThermo.Wi(i);
     }
     Xinf /= sum(Xinf);
 
@@ -325,7 +326,7 @@ void Foam::ReactingParcel<ParcelType>::correctSurfaceValues
         const scalar Csi = Cs[i] + Xsff*Xinf[i]*CsTot;
 
         Xs[i] = (2.0*Csi + Xinf[i]*CsTot)/3.0;
-        Ys[i] = Xs[i]*carrier.Wi(i);
+        Ys[i] = Xs[i]*carrierThermo.Wi(i);
     }
     Xs /= sum(Xs);
     Ys /= sum(Ys);
@@ -340,14 +341,14 @@ void Foam::ReactingParcel<ParcelType>::correctSurfaceValues
 
     forAll(Ys, i)
     {
-        const scalar W = carrier.Wi(i);
+        const scalar W = carrierThermo.Wi(i);
         const scalar sqrtW = sqrt(W);
         const scalar cbrtW = cbrt(W);
 
         rhos += Xs[i]*W;
-        mus += Ys[i]*sqrtW*carrier.mu(i, td.pc(), T);
-        kappas += Ys[i]*cbrtW*carrier.kappa(i, td.pc(), T);
-        Cps += Xs[i]*carrier.Cp(i, td.pc(), T);
+        mus += Ys[i]*sqrtW*carrierThermo.mui(i, td.pc(), T);
+        kappas += Ys[i]*cbrtW*carrierThermo.kappai(i, td.pc(), T);
+        Cps += Xs[i]*carrierThermo.Cpi(i, td.pc(), T);
 
         sumYiSqrtW += Ys[i]*sqrtW;
         sumYiCbrtW += Ys[i]*cbrtW;
@@ -494,7 +495,7 @@ void Foam::ReactingParcel<ParcelType>::calc
             {
                 scalar dmi = dm*Y_[i];
                 label gid = composition.localToCarrierId(0, i);
-                scalar hs = composition.carrier().Hs(gid, td.pc(), T0);
+                scalar hs = composition.carrier().hsi(gid, td.pc(), T0);
 
                 cloud.rhoTrans(gid)[this->cell()] += dmi;
                 cloud.hsTransRef()[this->cell()] += dmi*hs;
@@ -555,7 +556,7 @@ void Foam::ReactingParcel<ParcelType>::calc
         {
             scalar dm = np0*dMass[i];
             label gid = composition.localToCarrierId(0, i);
-            scalar hs = composition.carrier().Hs(gid, td.pc(), T0);
+            scalar hs = composition.carrier().hsi(gid, td.pc(), T0);
 
             cloud.rhoTrans(gid)[this->cell()] += dm;
             cloud.UTransRef()[this->cell()] += dm*U0;

src/lagrangian/parcel/submodels/Reacting/CompositionModel/CompositionModel/CompositionModel.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2011-2021 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2023 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -24,6 +24,7 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "CompositionModel.H"
+#include "fluidMulticomponentThermo.H"
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
@@ -32,10 +33,10 @@ Foam::CompositionModel<CloudType>::CompositionModel(CloudType& owner)
 :
     CloudSubModelBase<CloudType>(owner),
     carrierThermo_(owner.carrierThermo()),
-    carrierMixture_
+    carrierMcThermoPtr_
     (
-        isA<basicSpecieMixture>(carrierThermo_)
-      ? &refCast<const basicSpecieMixture>(carrierThermo_)
+        isA<fluidMulticomponentThermo>(carrierThermo_)
+      ? &refCast<const fluidMulticomponentThermo>(carrierThermo_)
       : nullptr
     ),
     thermo_(owner.thermo()),
@@ -53,19 +54,19 @@ Foam::CompositionModel<CloudType>::CompositionModel
 :
     CloudSubModelBase<CloudType>(owner, dict, typeName, type),
     carrierThermo_(owner.carrierThermo()),
-    carrierMixture_
+    carrierMcThermoPtr_
     (
-        isA<basicSpecieMixture>(carrierThermo_)
-      ? &refCast<const basicSpecieMixture>(carrierThermo_)
+        isA<fluidMulticomponentThermo>(carrierThermo_)
+      ? &refCast<const fluidMulticomponentThermo>(carrierThermo_)
       : nullptr
     ),
     thermo_(owner.thermo()),
     phaseProps_
     (
         this->coeffDict().lookup("phases"),
-        carrierMixture_ == nullptr
+        carrierMcThermoPtr_ == nullptr
       ? hashedWordList::null()
-      : carrierMixture_->species(),
+      : carrierMcThermoPtr_->species(),
         thermo_.liquids().components(),
         thermo_.solids().components()
     )
@@ -80,7 +81,7 @@ Foam::CompositionModel<CloudType>::CompositionModel
 :
     CloudSubModelBase<CloudType>(cm),
     carrierThermo_(cm.carrierThermo_),
-    carrierMixture_(cm.carrierMixture_),
+    carrierMcThermoPtr_(cm.carrierMcThermoPtr_),
     thermo_(cm.thermo_),
     phaseProps_(cm.phaseProps_)
 {}
@@ -103,17 +104,17 @@ const Foam::parcelThermo& Foam::CompositionModel<CloudType>::thermo() const
 
 
 template<class CloudType>
-const Foam::basicSpecieMixture&
+const Foam::fluidMulticomponentThermo&
 Foam::CompositionModel<CloudType>::carrier() const
 {
-    if (carrierMixture_ == nullptr)
+    if (carrierMcThermoPtr_ == nullptr)
     {
         FatalErrorInFunction
             << "carrier requested, but object is not allocated"
             << abort(FatalError);
     }
 
-    return *carrierMixture_;
+    return *carrierMcThermoPtr_;
 }
 
 
@@ -187,9 +188,9 @@ Foam::label Foam::CompositionModel<CloudType>::carrierId
 {
     label id = -1;
 
-    forAll(carrierMixture_->species(), i)
+    forAll(carrierMcThermoPtr_->species(), i)
     {
-        if (cmptName == carrierMixture_->species()[i])
+        if (cmptName == carrierMcThermoPtr_->species()[i])
         {
             id = i;
         }
@@ -200,7 +201,7 @@ Foam::label Foam::CompositionModel<CloudType>::carrierId
         FatalErrorInFunction
             << "Unable to determine global id for requested component "
             << cmptName << ". Available components are " << nl
-            << carrierMixture_->species()
+            << carrierMcThermoPtr_->species()
             << abort(FatalError);
     }
 
@@ -278,7 +279,7 @@ Foam::scalarField Foam::CompositionModel<CloudType>::X
             forAll(Y, i)
             {
                 label cid = props.carrierId(i);
-                X[i] = Y[i]/carrierMixture_->Wi(cid);
+                X[i] = Y[i]/carrierMcThermoPtr_->Wi(cid);
                 WInv += X[i];
             }
             break;
@@ -324,7 +325,7 @@ Foam::scalar Foam::CompositionModel<CloudType>::H
             forAll(Y, i)
             {
                 label cid = props.carrierId(i);
-                HMixture += Y[i]*carrierMixture_->Ha(cid, p, T);
+                HMixture += Y[i]*carrierMcThermoPtr_->hai(cid, p, T);
             }
             break;
         }
@@ -373,7 +374,7 @@ Foam::scalar Foam::CompositionModel<CloudType>::Hs
             forAll(Y, i)
             {
                 label cid = props.carrierId(i);
-                HsMixture += Y[i]*carrierMixture_->Hs(cid, p, T);
+                HsMixture += Y[i]*carrierMcThermoPtr_->hsi(cid, p, T);
             }
             break;
         }
@@ -424,7 +425,7 @@ Foam::scalar Foam::CompositionModel<CloudType>::Hc
             forAll(Y, i)
             {
                 label cid = props.carrierId(i);
-                HcMixture += Y[i]*carrierMixture_->Hf(cid);
+                HcMixture += Y[i]*carrierMcThermoPtr_->hfi(cid);
             }
             break;
         }
@@ -474,7 +475,7 @@ Foam::scalar Foam::CompositionModel<CloudType>::Cp
             forAll(Y, i)
             {
                 label cid = props.carrierId(i);
-                CpMixture += Y[i]*carrierMixture_->Cp(cid, p, T);
+                CpMixture += Y[i]*carrierMcThermoPtr_->Cpi(cid, p, T);
             }
             break;
         }

src/lagrangian/parcel/submodels/Reacting/CompositionModel/CompositionModel/CompositionModel.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2011-2021 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2023 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -40,7 +40,7 @@ SourceFiles
 
 #include "CloudSubModelBase.H"
 #include "parcelThermo.H"
-#include "basicSpecieMixture.H"
+#include "fluidMulticomponentThermo.H"
 #include "phasePropertiesList.H"
 #include "runTimeSelectionTables.H"
 
@@ -60,11 +60,11 @@ class CompositionModel
 {
     // Private Data
 
-        //- Reference to thermophysical properties of the carrier fluid
+        //- Reference to thermo properties of the carrier fluid
         const fluidThermo& carrierThermo_;
 
-        //- Pointer to the multi-component carrier composition
-        const basicSpecieMixture* carrierMixture_;
+        //- Reference to multicomponent thermo properties of the carrier fluid
+        const fluidMulticomponentThermo* carrierMcThermoPtr_;
 
         //- Reference to the thermo database
         const parcelThermo& thermo_;
@@ -135,7 +135,7 @@ public:
             // Composition lists
 
                 //- Return the carrier components (wrapper function)
-                const basicSpecieMixture& carrier() const;
+                const fluidMulticomponentThermo& carrier() const;
 
                 //- Return the global (additional) liquids
                 const liquidMixtureProperties& liquids() const;

src/lagrangian/parcel/submodels/Reacting/PhaseChangeModel/LiquidEvaporation/LiquidEvaporation.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2011-2022 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2023 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -243,7 +243,7 @@ Foam::scalar Foam::LiquidEvaporation<CloudType>::dh
         }
         case (parent::etEnthalpyDifference):
         {
-            scalar hc = this->owner().composition().carrier().Ha(idc, p, T);
+            scalar hc = this->owner().composition().carrier().hai(idc, p, T);
             scalar hp = liquids_.properties()[idl].Ha(p, T);
 
             dh = hc - hp;

src/lagrangian/parcel/submodels/Reacting/PhaseChangeModel/LiquidEvaporationBoil/LiquidEvaporationBoil.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2011-2022 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2023 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -182,10 +182,10 @@ void Foam::LiquidEvaporationBoil<CloudType>::calculate
     forAll(this->owner().composition().carrier().Y(), i)
     {
         scalar Yc = this->owner().composition().carrier().Y()[i][p.cell()];
-        Hc += Yc*this->owner().composition().carrier().Ha(i, pc, Tc);
-        Hsc += Yc*this->owner().composition().carrier().Ha(i, ps, Ts);
-        Cpc += Yc*this->owner().composition().carrier().Cp(i, ps, Ts);
-        kappac += Yc*this->owner().composition().carrier().kappa(i, ps, Ts);
+        Hc += Yc*this->owner().composition().carrier().hai(i, pc, Tc);
+        Hsc += Yc*this->owner().composition().carrier().hai(i, ps, Ts);
+        Cpc += Yc*this->owner().composition().carrier().Cpi(i, ps, Ts);
+        kappac += Yc*this->owner().composition().carrier().kappai(i, ps, Ts);
     }
 
     // calculate mass transfer of each specie in liquid
@@ -326,7 +326,8 @@ Foam::scalar Foam::LiquidEvaporationBoil<CloudType>::dh
         }
         case (parent::etEnthalpyDifference):
         {
-            scalar hc = this->owner().composition().carrier().Ha(idc, p, TDash);
+            scalar hc =
+                this->owner().composition().carrier().hai(idc, p, TDash);
             scalar hp = liquids_.properties()[idl].Ha(p, TDash);
 
             dh = hc - hp;

src/lagrangian/parcel/submodels/ReactingMultiphase/SurfaceReactionModel/COxidationDiffusionLimitedRate/COxidationDiffusionLimitedRate.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2011-2021 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2023 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -56,7 +56,7 @@ Foam::COxidationDiffusionLimitedRate<CloudType>::COxidationDiffusionLimitedRate
     const scalar WCO2 = owner.composition().carrier().Wi(CO2GlobalId_);
     WC_ = WCO2 - WO2_;
 
-    HcCO2_ = owner.composition().carrier().Hf(CO2GlobalId_);
+    HcCO2_ = owner.composition().carrier().hfi(CO2GlobalId_);
 
     if (Sb_ < 0)
     {
@@ -132,10 +132,11 @@ Foam::scalar Foam::COxidationDiffusionLimitedRate<CloudType>::calculate
     }
 
     const parcelThermo& thermo = this->owner().thermo();
-    const basicSpecieMixture& carrier = this->owner().composition().carrier();
+    const fluidMulticomponentThermo& carrierThermo =
+        this->owner().composition().carrier();
 
     // Local mass fraction of O2 in the carrier phase
-    const scalar YO2 = carrier.Y(O2GlobalId_)[celli];
+    const scalar YO2 = carrierThermo.Y(O2GlobalId_)[celli];
 
     // Change in C mass [kg]
     scalar dmC = 4.0*mathematical::pi*d*D_*YO2*Tc*rhoc/(Sb_*(T + Tc))*dt;

src/lagrangian/parcel/submodels/ReactingMultiphase/SurfaceReactionModel/COxidationHurtMitchell/COxidationHurtMitchell.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2012-2021 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2012-2023 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -56,7 +56,7 @@ Foam::COxidationHurtMitchell<CloudType>::COxidationHurtMitchell
     const scalar WCO2 = owner.composition().carrier().Wi(CO2GlobalId_);
     WC_ = WCO2 - WO2_;
 
-    HcCO2_ = owner.composition().carrier().Hf(CO2GlobalId_);
+    HcCO2_ = owner.composition().carrier().hfi(CO2GlobalId_);
 
     const scalar YCloc = owner.composition().Y0(idSolid)[CsLocalId_];
     const scalar YSolidTot = owner.composition().YMixture0()[idSolid];
@@ -131,10 +131,11 @@ Foam::scalar Foam::COxidationHurtMitchell<CloudType>::calculate
     }
 
     const parcelThermo& thermo = this->owner().thermo();
-    const basicSpecieMixture& carrier = this->owner().composition().carrier();
+    const fluidMulticomponentThermo& carrierThermo =
+        this->owner().composition().carrier();
 
     // Local mass fraction of O2 in the carrier phase
-    const scalar YO2 = carrier.Y(O2GlobalId_)[celli];
+    const scalar YO2 = carrierThermo.Y(O2GlobalId_)[celli];
 
     // No combustion if no oxygen present
     if (YO2 < small)

src/lagrangian/parcel/submodels/ReactingMultiphase/SurfaceReactionModel/COxidationIntrinsicRate/COxidationIntrinsicRate.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2014-2021 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2014-2023 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -62,7 +62,7 @@ Foam::COxidationIntrinsicRate<CloudType>::COxidationIntrinsicRate
     const scalar WCO2 = owner.composition().carrier().Wi(CO2GlobalId_);
     WC_ = WCO2 - WO2_;
 
-    HcCO2_ = owner.composition().carrier().Hf(CO2GlobalId_);
+    HcCO2_ = owner.composition().carrier().hfi(CO2GlobalId_);
 
     if (Sb_ < 0)
     {
@@ -144,10 +144,11 @@ Foam::scalar Foam::COxidationIntrinsicRate<CloudType>::calculate
     }
 
     const parcelThermo& thermo = this->owner().thermo();
-    const basicSpecieMixture& carrier = this->owner().composition().carrier();
+    const fluidMulticomponentThermo& carrierThermo =
+        this->owner().composition().carrier();
 
     // Local mass fraction of O2 in the carrier phase []
-    const scalar YO2 = carrier.Y(O2GlobalId_)[celli];
+    const scalar YO2 = carrierThermo.Y(O2GlobalId_)[celli];
 
     // Quick exit if oxidant not present
     if (YO2 < rootVSmall)

src/lagrangian/parcel/submodels/ReactingMultiphase/SurfaceReactionModel/COxidationKineticDiffusionLimitedRate/COxidationKineticDiffusionLimitedRate.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2011-2021 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2023 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -57,7 +57,7 @@ COxidationKineticDiffusionLimitedRate
     const scalar WCO2 = owner.composition().carrier().Wi(CO2GlobalId_);
     WC_ = WCO2 - WO2_;
 
-    HcCO2_ = owner.composition().carrier().Hf(CO2GlobalId_);
+    HcCO2_ = owner.composition().carrier().hfi(CO2GlobalId_);
 
     const scalar YCloc = owner.composition().Y0(idSolid)[CsLocalId_];
     const scalar YSolidTot = owner.composition().YMixture0()[idSolid];
@@ -129,10 +129,11 @@ Foam::scalar Foam::COxidationKineticDiffusionLimitedRate<CloudType>::calculate
     }
 
     const parcelThermo& thermo = this->owner().thermo();
-    const basicSpecieMixture& carrier = this->owner().composition().carrier();
+    const fluidMulticomponentThermo& carrierThermo =
+        this->owner().composition().carrier();
 
     // Local mass fraction of O2 in the carrier phase
-    const scalar YO2 = carrier.Y(O2GlobalId_)[celli];
+    const scalar YO2 = carrierThermo.Y(O2GlobalId_)[celli];
 
     // Diffusion rate coefficient
     const scalar D0 = C1_/d*pow(0.5*(T + Tc), 0.75);

src/lagrangian/parcel/submodels/ReactingMultiphase/SurfaceReactionModel/COxidationMurphyShaddix/COxidationMurphyShaddix.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2011-2021 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2023 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -69,7 +69,7 @@ Foam::COxidationMurphyShaddix<CloudType>::COxidationMurphyShaddix
     const scalar WCO2 = owner.composition().carrier().Wi(CO2GlobalId_);
     WC_ = WCO2 - WO2_;
 
-    HcCO2_ = owner.composition().carrier().Hf(CO2GlobalId_);
+    HcCO2_ = owner.composition().carrier().hfi(CO2GlobalId_);
 
     const scalar YCloc = owner.composition().Y0(idSolid)[CsLocalId_];
     const scalar YSolidTot = owner.composition().YMixture0()[idSolid];
@@ -143,10 +143,11 @@ Foam::scalar Foam::COxidationMurphyShaddix<CloudType>::calculate
     }
 
     const parcelThermo& thermo = this->owner().thermo();
-    const basicSpecieMixture& carrier = this->owner().composition().carrier();
+    const fluidMulticomponentThermo& carrierThermo =
+        this->owner().composition().carrier();
 
     // Cell carrier phase O2 species density [kg/m^3]
-    const scalar rhoO2 = rhoc*carrier.Y(O2GlobalId_)[celli];
+    const scalar rhoO2 = rhoc*carrierThermo.Y(O2GlobalId_)[celli];
 
     if (rhoO2 < small)
     {

src/radiationModels/absorptionEmissionModels/greyMean/greyMean.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2011-2021 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2023 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -27,7 +27,7 @@ License
 #include "addToRunTimeSelectionTable.H"
 #include "unitConversion.H"
 #include "extrapolatedCalculatedFvPatchFields.H"
-#include "basicSpecieMixture.H"
+#include "fluidMulticomponentThermo.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
@@ -67,7 +67,7 @@ Foam::radiationModels::absorptionEmissionModels::greyMean::greyMean
     thermo_(mesh.lookupObject<fluidThermo>(physicalProperties::typeName)),
     Yj_(nSpecies_)
 {
-    if (!isA<basicSpecieMixture>(thermo_))
+    if (!isA<fluidMulticomponentThermo>(thermo_))
     {
         FatalErrorInFunction
             << "Model requires a multi-component thermo package"
@@ -180,8 +180,8 @@ Foam::radiationModels::absorptionEmissionModels::greyMean::aCont
     const label bandI
 ) const
 {
-    const basicSpecieMixture& mixture =
-        dynamic_cast<const basicSpecieMixture&>(thermo_);
+    const fluidMulticomponentThermo& mcThermo =
+        dynamic_cast<const fluidMulticomponentThermo&>(thermo_);
 
     const volScalarField& T = thermo_.T();
     const volScalarField& p = thermo_.p();
@@ -219,13 +219,13 @@ Foam::radiationModels::absorptionEmissionModels::greyMean::aCont
             else
             {
                 scalar invWt = 0.0;
-                forAll(mixture.Y(), s)
+                forAll(mcThermo.Y(), s)
                 {
-                    invWt += mixture.Y(s)[celli]/mixture.Wi(s);
+                    invWt += mcThermo.Y(s)[celli]/mcThermo.Wi(s);
                 }
 
-                label index = mixture.species()[iter.key()];
-                scalar Xk = mixture.Y(index)[celli]/(mixture.Wi(index)*invWt);
+                label index = mcThermo.species()[iter.key()];
+                scalar Xk = mcThermo.Y(index)[celli]/(mcThermo.Wi(index)*invWt);
 
                 Xipi = Xk*paToAtm(p[celli]);
             }

src/radiationModels/absorptionEmissionModels/wideBand/wideBand.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2011-2021 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2023 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -25,7 +25,7 @@ License
 
 #include "wideBand.H"
 #include "addToRunTimeSelectionTable.H"
-#include "basicSpecieMixture.H"
+#include "fluidMulticomponentThermo.H"
 #include "unitConversion.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
@@ -67,6 +67,13 @@ Foam::radiationModels::absorptionEmissionModels::wideBand::wideBand
     Yj_(nSpecies_),
     totalWaveLength_(0)
 {
+    if (!isA<fluidMulticomponentThermo>(thermo_))
+    {
+        FatalErrorInFunction
+            << "Model requires a multi-component thermo package"
+            << abort(FatalError);
+    }
+
     label nBand = 0;
     forAllConstIter(dictionary, coeffsDict_, iter)
     {
@@ -195,8 +202,8 @@ Foam::radiationModels::absorptionEmissionModels::wideBand::aCont
     const label bandi
 ) const
 {
-    const basicSpecieMixture& mixture =
-        dynamic_cast<const basicSpecieMixture&>(thermo_);
+    const fluidMulticomponentThermo& mcThermo =
+        dynamic_cast<const fluidMulticomponentThermo&>(thermo_);
 
     const volScalarField& T = thermo_.T();
     const volScalarField& p = thermo_.p();
@@ -232,15 +239,15 @@ Foam::radiationModels::absorptionEmissionModels::wideBand::aCont
             else
             {
                 scalar invWt = 0;
-                forAll(mixture.Y(), s)
+                forAll(mcThermo.Y(), s)
                 {
-                    invWt += mixture.Y(s)[celli]/mixture.Wi(s);
+                    invWt += mcThermo.Y(s)[celli]/mcThermo.Wi(s);
                 }
 
-                const label index = mixture.species()[iter.key()];
+                const label index = mcThermo.species()[iter.key()];
 
                 const scalar Xk =
-                    mixture.Y(index)[celli]/(mixture.Wi(index)*invWt);
+                    mcThermo.Y(index)[celli]/(mcThermo.Wi(index)*invWt);
 
                 Xipi = Xk*paToAtm(p[celli]);
             }

src/specieTransfer/derivedFvPatchFields/adsorptionMassFraction/adsorptionMassFractionFvPatchScalarField.C

@@ -27,7 +27,7 @@ License
 #include "volFields.H"
 #include "surfaceFields.H"
 #include "fluidThermophysicalTransportModel.H"
-#include "basicSpecieMixture.H"
+#include "fluidMulticomponentThermo.H"
 #include "addToRunTimeSelectionTable.H"
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
@@ -100,8 +100,8 @@ Foam::adsorptionMassFractionFvPatchScalarField::calcPhiYp() const
     scalar Wi = NaN;
     if (property_ != massFraction)
     {
-        const basicSpecieMixture& mixture = composition(db());
-        Wi = mixture.Wi(mixture.species()[YName]);
+        const fluidMulticomponentThermo& thermo = this->thermo(db());
+        Wi = thermo.Wi(thermo.species()[YName]);
     }
 
     // Get the mixture molecular weights, if needed

src/specieTransfer/derivedFvPatchFields/adsorptionMassFraction/adsorptionMassFractionFvPatchScalarField.H

@@ -91,8 +91,6 @@ SourceFiles
 namespace Foam
 {
 
-class basicSpecieMixture;
-
 /*---------------------------------------------------------------------------*\
           Class adsorptionMassFractionFvPatchScalarField Declaration
 \*---------------------------------------------------------------------------*/

src/specieTransfer/derivedFvPatchFields/semiPermeableBaffleMassFraction/semiPermeableBaffleMassFractionFvPatchScalarField.C

@@ -28,7 +28,7 @@ License
 #include "volFields.H"
 #include "surfaceFields.H"
 #include "fluidThermophysicalTransportModel.H"
-#include "basicSpecieMixture.H"
+#include "fluidMulticomponentThermo.H"
 #include "mappedPatchBase.H"
 #include "addToRunTimeSelectionTable.H"
 
@@ -134,8 +134,8 @@ Foam::semiPermeableBaffleMassFractionFvPatchScalarField::calcPhiYp() const
     scalar Wi = NaN;
     if (property_ != massFraction)
     {
-        const basicSpecieMixture& mixture = composition(db());
-        Wi = mixture.Wi(mixture.species()[YName]);
+        const fluidMulticomponentThermo& thermo = this->thermo(db());
+        Wi = thermo.Wi(thermo.species()[YName]);
     }
 
     // Get the mixture molecular weights, if needed

src/specieTransfer/derivedFvPatchFields/semiPermeableBaffleMassFraction/semiPermeableBaffleMassFractionFvPatchScalarField.H

@@ -92,8 +92,6 @@ SourceFiles
 namespace Foam
 {
 
-class basicSpecieMixture;
-
 /*---------------------------------------------------------------------------*\
       Class semiPermeableBaffleMassFractionFvPatchScalarField Declaration
 \*---------------------------------------------------------------------------*/

src/specieTransfer/derivedFvPatchFields/specieTransferMassFraction/specieTransferMassFractionFvPatchScalarField.C

@@ -57,15 +57,17 @@ const Foam::NamedEnum
 
 // * * * * * * * * * * * * * Static Member Functions * * * * * * * * * * * * //
 
-const Foam::basicSpecieMixture&
-Foam::specieTransferMassFractionFvPatchScalarField::composition
+const Foam::fluidMulticomponentThermo&
+Foam::specieTransferMassFractionFvPatchScalarField::thermo
 (
     const objectRegistry& db
 )
 {
-    const word& name = physicalProperties::typeName;
-
-    return db.lookupObject<fluidMulticomponentThermo>(name).composition();
+    return
+        db.lookupObject<fluidMulticomponentThermo>
+        (
+            physicalProperties::typeName
+        );
 }
 
 

src/specieTransfer/derivedFvPatchFields/specieTransferMassFraction/specieTransferMassFractionFvPatchScalarField.H

@@ -50,7 +50,7 @@ SourceFiles
 namespace Foam
 {
 
-class basicSpecieMixture;
+class fluidMulticomponentThermo;
 
 /*---------------------------------------------------------------------------*\
         Class specieTransferMassFractionFvPatchScalarField Declaration
@@ -107,8 +107,11 @@ public:
 
     // Static member functions
 
-        //- Access the composition for the given database
-        static const basicSpecieMixture& composition(const objectRegistry& db);
+        //- Get the thermo from the given database
+        static const fluidMulticomponentThermo& thermo
+        (
+            const objectRegistry& db
+        );
 
 
     // Constructors

src/specieTransfer/derivedFvPatchFields/specieTransferTemperature/specieTransferTemperatureFvPatchScalarField.C

@@ -29,7 +29,7 @@ License
 #include "volFields.H"
 #include "surfaceFields.H"
 #include "fluidThermophysicalTransportModel.H"
-#include "basicSpecieMixture.H"
+#include "fluidMulticomponentThermo.H"
 #include "addToRunTimeSelectionTable.H"
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
@@ -88,12 +88,9 @@ const Foam::tmp<Foam::scalarField>
 Foam::specieTransferTemperatureFvPatchScalarField::phiHep() const
 {
     typedef specieTransferMassFractionFvPatchScalarField YBCType;
-    const basicSpecieMixture& mixture = YBCType::composition(db());
-    const PtrList<volScalarField>& Y = mixture.Y();
+    const fluidMulticomponentThermo& thermo = YBCType::thermo(db());
+    const PtrList<volScalarField>& Y = thermo.Y();
 
-    // Get thermodynamic properties
-    const fluidThermo& thermo =
-        db().lookupObject<fluidThermo>(physicalProperties::typeName);
     const fvPatchScalarField& Tp = *this;
     const fvPatchScalarField& pp = thermo.p().boundaryField()[patch().index()];
 
@@ -112,7 +109,7 @@ Foam::specieTransferTemperatureFvPatchScalarField::phiHep() const
                 << exit(FatalError);
         }
 
-        phiHep += refCast<const YBCType>(Yp).phiYp()*mixture.HE(i, pp, Tp);
+        phiHep += refCast<const YBCType>(Yp).phiYp()*thermo.hei(i, pp, Tp);
     }
 
     return tPhiHep;

src/specieTransfer/derivedFvPatchFields/specieTransferVelocity/specieTransferVelocityFvPatchVectorField.C

@@ -28,7 +28,7 @@ License
 #include "addToRunTimeSelectionTable.H"
 #include "volFields.H"
 #include "surfaceFields.H"
-#include "basicSpecieMixture.H"
+#include "fluidMulticomponentThermo.H"
 #include "basicThermo.H"
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
@@ -78,7 +78,7 @@ const Foam::tmp<Foam::scalarField>
 Foam::specieTransferVelocityFvPatchVectorField::phip() const
 {
     typedef specieTransferMassFractionFvPatchScalarField YBCType;
-    const PtrList<volScalarField>& Y = YBCType::composition(db()).Y();
+    const PtrList<volScalarField>& Y = YBCType::thermo(db()).Y();
 
     // Sum up the phiYp-s from all the species
     tmp<scalarField> tPhip(new scalarField(this->size(), 0));

src/thermophysicalModels/basic/Make/files

@@ -1,5 +1,7 @@
 basicThermo/basicThermo.C
 
+pureThermo/pureThermo.C
+
 fluidThermo/fluidThermo.C
 fluidThermo/hydrostaticInitialisation.C
 
@@ -8,7 +10,9 @@ psiThermo/psiThermos.C
 
 rhoThermo/rhoThermo.C
 rhoThermo/rhoThermos.C
-rhoThermo/liquidThermo.C
+
+liquidThermo/liquidThermo.C
+liquidThermo/liquidThermos.C
 
 derivedFvPatchFields/energyJump/energyJumpFvPatchScalarField.C
 derivedFvPatchFields/fixedEnergy/fixedEnergyFvPatchScalarField.C
@@ -17,4 +21,6 @@ derivedFvPatchFields/gradientEnergy/gradientEnergyCalculatedTemperatureFvPatchSc
 derivedFvPatchFields/mixedEnergy/mixedEnergyFvPatchScalarField.C
 derivedFvPatchFields/mixedEnergy/mixedEnergyCalculatedTemperatureFvPatchScalarField.C
 
+heThermo/heThermoName.C
+
 LIB = $(FOAM_LIBBIN)/libfluidThermophysicalModels

src/thermophysicalModels/basic/PhysicalPropertiesThermo.C

@@ -0,0 +1,82 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     | Website:  https://openfoam.org
+    \\  /    A nd           | Copyright (C) 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 "PhysicalPropertiesThermo.H"
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+template<class BasicThermoType>
+Foam::PhysicalPropertiesThermo<BasicThermoType>::PhysicalPropertiesThermo
+(
+    const fvMesh& mesh,
+    const word& phaseName
+)
+:
+    physicalProperties(mesh, phaseName),
+    BasicThermoType(*this, mesh, phaseName)
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+template<class BasicThermoType>
+Foam::PhysicalPropertiesThermo<BasicThermoType>::~PhysicalPropertiesThermo()
+{}
+
+
+// * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
+
+template<class BasicThermoType>
+const Foam::IOdictionary&
+Foam::PhysicalPropertiesThermo<BasicThermoType>::properties() const
+{
+    return *this;
+}
+
+
+template<class BasicThermoType>
+Foam::IOdictionary&
+Foam::PhysicalPropertiesThermo<BasicThermoType>::properties()
+{
+    return *this;
+}
+
+
+template<class BasicThermoType>
+bool Foam::PhysicalPropertiesThermo<BasicThermoType>::read()
+{
+    if (physicalProperties::read())
+    {
+        BasicThermoType::read(*this);
+        return true;
+    }
+    else
+    {
+        return false;
+    }
+}
+
+
+// ************************************************************************* //

src/thermophysicalModels/basic/PhysicalPropertiesThermo.H

@@ -0,0 +1,96 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     | Website:  https://openfoam.org
+    \\  /    A nd           | Copyright (C) 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/>.
+
+Class
+    Foam::PhysicalPropertiesThermo
+
+Description
+    Wrapper around a thermo which also constructs the physical properties
+    dictionary
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef PhysicalPropertiesThermo_H
+#define PhysicalPropertiesThermo_H
+
+#include "physicalProperties.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                   Class PhysicalPropertiesThermo Declaration
+\*---------------------------------------------------------------------------*/
+
+template<class BasicThermoType>
+class PhysicalPropertiesThermo
+:
+    public physicalProperties,
+    public BasicThermoType
+{
+public:
+
+    // Constructors
+
+        //- Construct from mesh and phase name
+        PhysicalPropertiesThermo
+        (
+            const fvMesh& mesh,
+            const word& phaseName = word::null
+        );
+
+
+    //- Destructor
+    virtual ~PhysicalPropertiesThermo();
+
+
+    // Member Functions
+
+        //- Access the properties dictionary
+        virtual const IOdictionary& properties() const;
+
+        //- Non-const access the properties dictionary
+        virtual IOdictionary& properties();
+
+        //- Read the properties dictionary
+        virtual bool read();
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+    #include "PhysicalPropertiesThermo.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/thermophysicalModels/basic/basicThermo/basicThermo.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2011-2022 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2023 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -262,12 +262,11 @@ Foam::wordList Foam::basicThermo::heBoundaryTypes()
 
 Foam::basicThermo::implementation::implementation
 (
+    const dictionary& dict,
     const fvMesh& mesh,
     const word& phaseName
 )
 :
-    physicalProperties(mesh, phaseName),
-
     mesh_(mesh),
 
     phaseName_(phaseName),
@@ -299,7 +298,7 @@ Foam::basicThermo::implementation::implementation
         dimensionedScalar(dimEnergy/dimTime/dimLength/dimTemperature, Zero)
     ),
 
-    dpdt_(lookupOrDefault<Switch>("dpdt", true))
+    dpdt_(dict.lookupOrDefault<Switch>("dpdt", true))
 {}
 
 
@@ -384,10 +383,8 @@ const Foam::volScalarField& Foam::basicThermo::implementation::kappa() const
 }
 
 
-bool Foam::basicThermo::implementation::read()
-{
-    return regIOobject::read();
-}
+void Foam::basicThermo::implementation::read(const dictionary&)
+{}
 
 
 // ************************************************************************* //

src/thermophysicalModels/basic/basicThermo/basicThermo.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2011-2022 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2023 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -415,8 +415,7 @@ public:
 
 class basicThermo::implementation
 :
-    virtual public basicThermo,
-    public physicalProperties
+    virtual public basicThermo
 {
     // Private Member Data
 
@@ -447,16 +446,10 @@ protected:
 
 public:
 
-    // Static Member data
-
-        //- Name of the thermophysical properties dictionary
-        using physicalProperties::typeName;
-
-
     // Constructors
 
-        //- Construct from mesh and phase name
-        implementation(const fvMesh&, const word& phaseName);
+        //- Construct from dictionary, mesh and phase name
+        implementation(const dictionary&, const fvMesh&, const word&);
 
         //- Disallow default bitwise copy construction
         implementation(const implementation&) = delete;
@@ -469,18 +462,6 @@ public:
 
     // Member Functions
 
-        //- Properties dictionary
-        virtual const IOdictionary& properties() const
-        {
-            return *this;
-        }
-
-        //- Non-const access the properties dictionary
-        virtual IOdictionary& properties()
-        {
-            return *this;
-        }
-
         //- Return const access to the mesh
         virtual const fvMesh& mesh() const
         {
@@ -517,7 +498,7 @@ public:
 
 
         //- Read thermophysical properties dictionary
-        virtual bool read();
+        virtual void read(const dictionary&);
 
 
     // Member Operators
@@ -526,6 +507,7 @@ public:
         void operator=(const implementation&) = delete;
 };
 
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 } // End namespace Foam

src/thermophysicalModels/basic/basicThermo/basicThermoTemplates.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2012-2022 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2012-2023 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -55,7 +55,8 @@ typename Table::iterator Foam::basicThermo::lookupCstrIter
     {
         if
         (
-            dynamicCode::allowSystemOperations
+            nCmpt == 7
+         && dynamicCode::allowSystemOperations
          && !dynamicCode::resolveTemplate(Thermo::typeName).empty()
         )
         {

src/thermophysicalModels/basic/fluidThermo/fluidThermo.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2012-2022 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2012-2023 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -38,6 +38,7 @@ namespace Foam
 
 Foam::fluidThermo::implementation::implementation
 (
+    const dictionary& dict,
     const fvMesh& mesh,
     const word& phaseName
 )

src/thermophysicalModels/basic/fluidThermo/fluidThermo.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2012-2022 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2012-2023 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -180,8 +180,8 @@ public:
 
     // Constructors
 
-        //- Construct from mesh and phase name
-        implementation(const fvMesh&, const word& phaseName);
+        //- Construct from dictionary, mesh and phase name
+        implementation(const dictionary&, const fvMesh&, const word&);
 
         //- Disallow default bitwise copy construction
         implementation(const implementation&) = delete;

src/thermophysicalModels/basic/heThermo/heThermo.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2011-2022 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2023 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -29,21 +29,14 @@ License
 
 // * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * //
 
-template<class BasicThermo, class MixtureType>
-template
-<
-    class CellMixture,
-    class PatchFaceMixture,
-    class Method,
-    class ... Args
->
+template<class MixtureType, class BasicThermoType>
+template<class Mixture, class Method, class ... Args>
 Foam::tmp<Foam::volScalarField>
-Foam::heThermo<BasicThermo, MixtureType>::volScalarFieldProperty
+Foam::heThermo<MixtureType, BasicThermoType>::volScalarFieldProperty
 (
     const word& psiName,
     const dimensionSet& psiDim,
-    CellMixture cellMixture,
-    PatchFaceMixture patchFaceMixture,
+    Mixture mixture,
     Method psiMethod,
     const Args& ... args
 ) const
@@ -60,23 +53,26 @@ Foam::heThermo<BasicThermo, MixtureType>::volScalarFieldProperty
 
     volScalarField& psi = tPsi.ref();
 
-    forAll(this->T_, celli)
+    forAll(psi, celli)
     {
-        psi[celli] = ((this->*cellMixture)(celli).*psiMethod)(args[celli] ...);
+        auto composition = this->cellComposition(celli);
+
+        psi[celli] =
+            ((this->*mixture)(composition).*psiMethod)(args[celli] ...);
     }
 
     volScalarField::Boundary& psiBf = psi.boundaryFieldRef();
 
     forAll(psiBf, patchi)
     {
-        fvPatchScalarField& pPsi = psiBf[patchi];
-
-        forAll(this->T_.boundaryField()[patchi], facei)
+        forAll(psiBf[patchi], patchFacei)
         {
-            pPsi[facei] =
-                ((this->*patchFaceMixture)(patchi, facei).*psiMethod)
+            auto composition = this->patchFaceComposition(patchi, patchFacei);
+
+            psiBf[patchi][patchFacei] =
+                ((this->*mixture)(composition).*psiMethod)
                 (
-                    args.boundaryField()[patchi][facei] ...
+                    args.boundaryField()[patchi][patchFacei] ...
                 );
         }
     }
@@ -85,12 +81,12 @@ Foam::heThermo<BasicThermo, MixtureType>::volScalarFieldProperty
 }
 
 
-template<class BasicThermo, class MixtureType>
-template<class CellMixture, class Method, class ... Args>
+template<class MixtureType, class BasicThermoType>
+template<class Mixture, class Method, class ... Args>
 Foam::tmp<Foam::scalarField>
-Foam::heThermo<BasicThermo, MixtureType>::cellSetProperty
+Foam::heThermo<MixtureType, BasicThermoType>::cellSetProperty
 (
-    CellMixture cellMixture,
+    Mixture mixture,
     Method psiMethod,
     const labelList& cells,
     const Args& ... args
@@ -104,20 +100,22 @@ Foam::heThermo<BasicThermo, MixtureType>::cellSetProperty
 
     forAll(cells, celli)
     {
+        auto composition = this->cellComposition(cells[celli]);
+
         psi[celli] =
-            ((this->*cellMixture)(cells[celli]).*psiMethod)(args[celli] ...);
+            ((this->*mixture)(composition).*psiMethod)(args[celli] ...);
     }
 
     return tPsi;
 }
 
 
-template<class BasicThermo, class MixtureType>
-template<class PatchFaceMixture, class Method, class ... Args>
+template<class MixtureType, class BasicThermoType>
+template<class Mixture, class Method, class ... Args>
 Foam::tmp<Foam::scalarField>
-Foam::heThermo<BasicThermo, MixtureType>::patchFieldProperty
+Foam::heThermo<MixtureType, BasicThermoType>::patchFieldProperty
 (
-    PatchFaceMixture patchFaceMixture,
+    Mixture mixture,
     Method psiMethod,
     const label patchi,
     const Args& ... args
@@ -129,22 +127,21 @@ Foam::heThermo<BasicThermo, MixtureType>::patchFieldProperty
     );
     scalarField& psi = tPsi.ref();
 
-    forAll(this->T_.boundaryField()[patchi], facei)
+    forAll(psi, patchFacei)
     {
-        psi[facei] =
-            ((this->*patchFaceMixture)(patchi, facei).*psiMethod)
-            (
-                args[facei] ...
-            );
+        auto composition = this->patchFaceComposition(patchi, patchFacei);
+
+        psi[patchFacei] =
+            ((this->*mixture)(composition).*psiMethod)(args[patchFacei] ...);
     }
 
     return tPsi;
 }
 
 
-template<class BasicThermo, class MixtureType>
+template<class MixtureType, class BasicThermoType>
 Foam::UIndirectList<Foam::scalar>
-Foam::heThermo<BasicThermo, MixtureType>::cellSetScalarList
+Foam::heThermo<MixtureType, BasicThermoType>::cellSetScalarList
 (
     const volScalarField& psi,
     const labelList& cells
@@ -154,9 +151,9 @@ Foam::heThermo<BasicThermo, MixtureType>::cellSetScalarList
 }
 
 
-template<class BasicThermo, class MixtureType>
+template<class MixtureType, class BasicThermoType>
 Foam::UniformField<Foam::scalar>
-Foam::heThermo<BasicThermo, MixtureType>::cellSetScalarList
+Foam::heThermo<MixtureType, BasicThermoType>::cellSetScalarList
 (
     const uniformGeometricScalarField& psi,
     const labelList&
@@ -166,8 +163,8 @@ Foam::heThermo<BasicThermo, MixtureType>::cellSetScalarList
 }
 
 
-template<class BasicThermo, class MixtureType>
-void Foam::heThermo<BasicThermo, MixtureType>::
+template<class MixtureType, class BasicThermoType>
+void Foam::heThermo<MixtureType, BasicThermoType>::
 heBoundaryCorrection(volScalarField& h)
 {
     volScalarField::Boundary& hBf = h.boundaryFieldRef();
@@ -190,21 +187,22 @@ heBoundaryCorrection(volScalarField& h)
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
-template<class BasicThermo, class MixtureType>
-Foam::heThermo<BasicThermo, MixtureType>::heThermo
+template<class MixtureType, class BasicThermoType>
+Foam::heThermo<MixtureType, BasicThermoType>::heThermo
 (
     const fvMesh& mesh,
     const word& phaseName
 )
 :
-    BasicThermo(mesh, phaseName),
-    MixtureType(*this, mesh, phaseName),
+    physicalProperties(mesh, phaseName),
+    MixtureType(properties()),
+    BasicThermoType(properties(), mixture(), mesh, phaseName),
 
     he_
     (
         IOobject
         (
-            BasicThermo::phasePropertyName
+            BasicThermoType::phasePropertyName
             (
                 MixtureType::thermoType::heName(),
                 phaseName
@@ -218,8 +216,7 @@ Foam::heThermo<BasicThermo, MixtureType>::heThermo
         (
             "he",
             dimEnergy/dimMass,
-            &MixtureType::cellThermoMixture,
-            &MixtureType::patchFaceThermoMixture,
+            &MixtureType::thermoMixture,
             &MixtureType::thermoMixtureType::HE,
             this->p_,
             this->T_
@@ -232,7 +229,7 @@ Foam::heThermo<BasicThermo, MixtureType>::heThermo
     (
         IOobject
         (
-            BasicThermo::phasePropertyName("Cp", phaseName),
+            BasicThermoType::phasePropertyName("Cp", phaseName),
             mesh.time().name(),
             mesh
         ),
@@ -244,7 +241,7 @@ Foam::heThermo<BasicThermo, MixtureType>::heThermo
     (
         IOobject
         (
-            BasicThermo::phasePropertyName("Cv", phaseName),
+            BasicThermoType::phasePropertyName("Cv", phaseName),
             mesh.time().name(),
             mesh
         ),
@@ -258,16 +255,21 @@ Foam::heThermo<BasicThermo, MixtureType>::heThermo
 
 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
 
-template<class BasicThermo, class MixtureType>
-Foam::heThermo<BasicThermo, MixtureType>::~heThermo()
+template<class MixtureType, class BasicThermoType>
+Foam::heThermo<MixtureType, BasicThermoType>::~heThermo()
+{}
+
+
+template<class HeThermo>
+Foam::namedHeThermo<HeThermo>::~namedHeThermo()
 {}
 
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-template<class BasicThermo, class MixtureType>
+template<class MixtureType, class BasicThermoType>
 const Foam::volScalarField&
-Foam::heThermo<BasicThermo, MixtureType>::Cpv() const
+Foam::heThermo<MixtureType, BasicThermoType>::Cpv() const
 {
     if (MixtureType::thermoType::enthalpy())
     {
@@ -280,8 +282,8 @@ Foam::heThermo<BasicThermo, MixtureType>::Cpv() const
 }
 
 
-template<class BasicThermo, class MixtureType>
-Foam::tmp<Foam::volScalarField> Foam::heThermo<BasicThermo, MixtureType>::he
+template<class MixtureType, class BasicThermoType>
+Foam::tmp<Foam::volScalarField> Foam::heThermo<MixtureType, BasicThermoType>::he
 (
     const volScalarField& p,
     const volScalarField& T
@@ -291,8 +293,7 @@ Foam::tmp<Foam::volScalarField> Foam::heThermo<BasicThermo, MixtureType>::he
     (
         "he",
         dimEnergy/dimMass,
-        &MixtureType::cellThermoMixture,
-        &MixtureType::patchFaceThermoMixture,
+        &MixtureType::thermoMixture,
         &MixtureType::thermoMixtureType::HE,
         p,
         T
@@ -300,8 +301,8 @@ Foam::tmp<Foam::volScalarField> Foam::heThermo<BasicThermo, MixtureType>::he
 }
 
 
-template<class BasicThermo, class MixtureType>
-Foam::tmp<Foam::scalarField> Foam::heThermo<BasicThermo, MixtureType>::he
+template<class MixtureType, class BasicThermoType>
+Foam::tmp<Foam::scalarField> Foam::heThermo<MixtureType, BasicThermoType>::he
 (
     const scalarField& T,
     const labelList& cells
@@ -309,7 +310,7 @@ Foam::tmp<Foam::scalarField> Foam::heThermo<BasicThermo, MixtureType>::he
 {
     return cellSetProperty
     (
-        &MixtureType::cellThermoMixture,
+        &MixtureType::thermoMixture,
         &MixtureType::thermoMixtureType::HE,
         cells,
         cellSetScalarList(this->p_, cells),
@@ -318,8 +319,8 @@ Foam::tmp<Foam::scalarField> Foam::heThermo<BasicThermo, MixtureType>::he
 }
 
 
-template<class BasicThermo, class MixtureType>
-Foam::tmp<Foam::scalarField> Foam::heThermo<BasicThermo, MixtureType>::he
+template<class MixtureType, class BasicThermoType>
+Foam::tmp<Foam::scalarField> Foam::heThermo<MixtureType, BasicThermoType>::he
 (
     const scalarField& T,
     const label patchi
@@ -327,7 +328,7 @@ Foam::tmp<Foam::scalarField> Foam::heThermo<BasicThermo, MixtureType>::he
 {
     return patchFieldProperty
     (
-        &MixtureType::patchFaceThermoMixture,
+        &MixtureType::thermoMixture,
         &MixtureType::thermoMixtureType::HE,
         patchi,
         this->p_.boundaryField()[patchi],
@@ -336,16 +337,15 @@ Foam::tmp<Foam::scalarField> Foam::heThermo<BasicThermo, MixtureType>::he
 }
 
 
-template<class BasicThermo, class MixtureType>
+template<class MixtureType, class BasicThermoType>
 Foam::tmp<Foam::volScalarField>
-Foam::heThermo<BasicThermo, MixtureType>::hs() const
+Foam::heThermo<MixtureType, BasicThermoType>::hs() const
 {
     return volScalarFieldProperty
     (
         "hs",
         dimEnergy/dimMass,
-        &MixtureType::cellThermoMixture,
-        &MixtureType::patchFaceThermoMixture,
+        &MixtureType::thermoMixture,
         &MixtureType::thermoMixtureType::Hs,
         this->p_,
         this->T_
@@ -353,8 +353,8 @@ Foam::heThermo<BasicThermo, MixtureType>::hs() const
 }
 
 
-template<class BasicThermo, class MixtureType>
-Foam::tmp<Foam::volScalarField> Foam::heThermo<BasicThermo, MixtureType>::hs
+template<class MixtureType, class BasicThermoType>
+Foam::tmp<Foam::volScalarField> Foam::heThermo<MixtureType, BasicThermoType>::hs
 (
     const volScalarField& p,
     const volScalarField& T
@@ -364,8 +364,7 @@ Foam::tmp<Foam::volScalarField> Foam::heThermo<BasicThermo, MixtureType>::hs
     (
         "hs",
         dimEnergy/dimMass,
-        &MixtureType::cellThermoMixture,
-        &MixtureType::patchFaceThermoMixture,
+        &MixtureType::thermoMixture,
         &MixtureType::thermoMixtureType::Hs,
         p,
         T
@@ -373,8 +372,8 @@ Foam::tmp<Foam::volScalarField> Foam::heThermo<BasicThermo, MixtureType>::hs
 }
 
 
-template<class BasicThermo, class MixtureType>
-Foam::tmp<Foam::scalarField> Foam::heThermo<BasicThermo, MixtureType>::hs
+template<class MixtureType, class BasicThermoType>
+Foam::tmp<Foam::scalarField> Foam::heThermo<MixtureType, BasicThermoType>::hs
 (
     const scalarField& T,
     const labelList& cells
@@ -382,7 +381,7 @@ Foam::tmp<Foam::scalarField> Foam::heThermo<BasicThermo, MixtureType>::hs
 {
     return cellSetProperty
     (
-        &MixtureType::cellThermoMixture,
+        &MixtureType::thermoMixture,
         &MixtureType::thermoMixtureType::Hs,
         cells,
         cellSetScalarList(this->p_, cells),
@@ -391,8 +390,8 @@ Foam::tmp<Foam::scalarField> Foam::heThermo<BasicThermo, MixtureType>::hs
 }
 
 
-template<class BasicThermo, class MixtureType>
-Foam::tmp<Foam::scalarField> Foam::heThermo<BasicThermo, MixtureType>::hs
+template<class MixtureType, class BasicThermoType>
+Foam::tmp<Foam::scalarField> Foam::heThermo<MixtureType, BasicThermoType>::hs
 (
     const scalarField& T,
     const label patchi
@@ -400,7 +399,7 @@ Foam::tmp<Foam::scalarField> Foam::heThermo<BasicThermo, MixtureType>::hs
 {
     return patchFieldProperty
     (
-        &MixtureType::patchFaceThermoMixture,
+        &MixtureType::thermoMixture,
         &MixtureType::thermoMixtureType::Hs,
         patchi,
         this->p_.boundaryField()[patchi],
@@ -409,16 +408,15 @@ Foam::tmp<Foam::scalarField> Foam::heThermo<BasicThermo, MixtureType>::hs
 }
 
 
-template<class BasicThermo, class MixtureType>
+template<class MixtureType, class BasicThermoType>
 Foam::tmp<Foam::volScalarField>
-Foam::heThermo<BasicThermo, MixtureType>::ha() const
+Foam::heThermo<MixtureType, BasicThermoType>::ha() const
 {
     return volScalarFieldProperty
     (
         "ha",
         dimEnergy/dimMass,
-        &MixtureType::cellThermoMixture,
-        &MixtureType::patchFaceThermoMixture,
+        &MixtureType::thermoMixture,
         &MixtureType::thermoMixtureType::Ha,
         this->p_,
         this->T_
@@ -426,8 +424,8 @@ Foam::heThermo<BasicThermo, MixtureType>::ha() const
 }
 
 
-template<class BasicThermo, class MixtureType>
-Foam::tmp<Foam::volScalarField> Foam::heThermo<BasicThermo, MixtureType>::ha
+template<class MixtureType, class BasicThermoType>
+Foam::tmp<Foam::volScalarField> Foam::heThermo<MixtureType, BasicThermoType>::ha
 (
     const volScalarField& p,
     const volScalarField& T
@@ -437,8 +435,7 @@ Foam::tmp<Foam::volScalarField> Foam::heThermo<BasicThermo, MixtureType>::ha
     (
         "ha",
         dimEnergy/dimMass,
-        &MixtureType::cellThermoMixture,
-        &MixtureType::patchFaceThermoMixture,
+        &MixtureType::thermoMixture,
         &MixtureType::thermoMixtureType::Ha,
         p,
         T
@@ -446,8 +443,8 @@ Foam::tmp<Foam::volScalarField> Foam::heThermo<BasicThermo, MixtureType>::ha
 }
 
 
-template<class BasicThermo, class MixtureType>
-Foam::tmp<Foam::scalarField> Foam::heThermo<BasicThermo, MixtureType>::ha
+template<class MixtureType, class BasicThermoType>
+Foam::tmp<Foam::scalarField> Foam::heThermo<MixtureType, BasicThermoType>::ha
 (
     const scalarField& T,
     const labelList& cells
@@ -455,7 +452,7 @@ Foam::tmp<Foam::scalarField> Foam::heThermo<BasicThermo, MixtureType>::ha
 {
     return cellSetProperty
     (
-        &MixtureType::cellThermoMixture,
+        &MixtureType::thermoMixture,
         &MixtureType::thermoMixtureType::Ha,
         cells,
         cellSetScalarList(this->p_, cells),
@@ -464,8 +461,8 @@ Foam::tmp<Foam::scalarField> Foam::heThermo<BasicThermo, MixtureType>::ha
 }
 
 
-template<class BasicThermo, class MixtureType>
-Foam::tmp<Foam::scalarField> Foam::heThermo<BasicThermo, MixtureType>::ha
+template<class MixtureType, class BasicThermoType>
+Foam::tmp<Foam::scalarField> Foam::heThermo<MixtureType, BasicThermoType>::ha
 (
     const scalarField& T,
     const label patchi
@@ -473,7 +470,7 @@ Foam::tmp<Foam::scalarField> Foam::heThermo<BasicThermo, MixtureType>::ha
 {
     return patchFieldProperty
     (
-        &MixtureType::patchFaceThermoMixture,
+        &MixtureType::thermoMixture,
         &MixtureType::thermoMixtureType::Ha,
         patchi,
         this->p_.boundaryField()[patchi],
@@ -482,23 +479,22 @@ Foam::tmp<Foam::scalarField> Foam::heThermo<BasicThermo, MixtureType>::ha
 }
 
 
-template<class BasicThermo, class MixtureType>
+template<class MixtureType, class BasicThermoType>
 Foam::tmp<Foam::volScalarField>
-Foam::heThermo<BasicThermo, MixtureType>::hc() const
+Foam::heThermo<MixtureType, BasicThermoType>::hc() const
 {
     return volScalarFieldProperty
     (
         "hc",
         dimEnergy/dimMass,
-        &MixtureType::cellThermoMixture,
-        &MixtureType::patchFaceThermoMixture,
+        &MixtureType::thermoMixture,
         &MixtureType::thermoMixtureType::Hf
     );
 }
 
 
-template<class BasicThermo, class MixtureType>
-Foam::tmp<Foam::scalarField> Foam::heThermo<BasicThermo, MixtureType>::Cp
+template<class MixtureType, class BasicThermoType>
+Foam::tmp<Foam::scalarField> Foam::heThermo<MixtureType, BasicThermoType>::Cp
 (
     const scalarField& T,
     const label patchi
@@ -506,7 +502,7 @@ Foam::tmp<Foam::scalarField> Foam::heThermo<BasicThermo, MixtureType>::Cp
 {
     return patchFieldProperty
     (
-        &MixtureType::patchFaceThermoMixture,
+        &MixtureType::thermoMixture,
         &MixtureType::thermoMixtureType::Cp,
         patchi,
         this->p_.boundaryField()[patchi],
@@ -515,9 +511,8 @@ Foam::tmp<Foam::scalarField> Foam::heThermo<BasicThermo, MixtureType>::Cp
 }
 
 
-template<class BasicThermo, class MixtureType>
-Foam::tmp<Foam::scalarField>
-Foam::heThermo<BasicThermo, MixtureType>::Cv
+template<class MixtureType, class BasicThermoType>
+Foam::tmp<Foam::scalarField> Foam::heThermo<MixtureType, BasicThermoType>::Cv
 (
     const scalarField& T,
     const label patchi
@@ -525,7 +520,7 @@ Foam::heThermo<BasicThermo, MixtureType>::Cv
 {
     return patchFieldProperty
     (
-        &MixtureType::patchFaceThermoMixture,
+        &MixtureType::thermoMixture,
         &MixtureType::thermoMixtureType::Cv,
         patchi,
         this->p_.boundaryField()[patchi],
@@ -534,8 +529,8 @@ Foam::heThermo<BasicThermo, MixtureType>::Cv
 }
 
 
-template<class BasicThermo, class MixtureType>
-Foam::tmp<Foam::scalarField> Foam::heThermo<BasicThermo, MixtureType>::gamma
+template<class MixtureType, class BasicThermoType>
+Foam::tmp<Foam::scalarField> Foam::heThermo<MixtureType, BasicThermoType>::gamma
 (
     const scalarField& T,
     const label patchi
@@ -543,7 +538,7 @@ Foam::tmp<Foam::scalarField> Foam::heThermo<BasicThermo, MixtureType>::gamma
 {
     return patchFieldProperty
     (
-        &MixtureType::patchFaceThermoMixture,
+        &MixtureType::thermoMixture,
         &MixtureType::thermoMixtureType::gamma,
         patchi,
         this->p_.boundaryField()[patchi],
@@ -552,16 +547,16 @@ Foam::tmp<Foam::scalarField> Foam::heThermo<BasicThermo, MixtureType>::gamma
 }
 
 
-template<class BasicThermo, class MixtureType>
+template<class MixtureType, class BasicThermoType>
 Foam::tmp<Foam::volScalarField>
-Foam::heThermo<BasicThermo, MixtureType>::gamma() const
+Foam::heThermo<MixtureType, BasicThermoType>::gamma() const
 {
     return volScalarField::New("gamma", Cp_/Cv_);
 }
 
 
-template<class BasicThermo, class MixtureType>
-Foam::tmp<Foam::scalarField> Foam::heThermo<BasicThermo, MixtureType>::Cpv
+template<class MixtureType, class BasicThermoType>
+Foam::tmp<Foam::scalarField> Foam::heThermo<MixtureType, BasicThermoType>::Cpv
 (
     const scalarField& T,
     const label patchi
@@ -578,8 +573,9 @@ Foam::tmp<Foam::scalarField> Foam::heThermo<BasicThermo, MixtureType>::Cpv
 }
 
 
-template<class BasicThermo, class MixtureType>
-Foam::tmp<Foam::volScalarField> Foam::heThermo<BasicThermo, MixtureType>::THE
+template<class MixtureType, class BasicThermoType>
+Foam::tmp<Foam::volScalarField>
+Foam::heThermo<MixtureType, BasicThermoType>::THE
 (
     const volScalarField& h,
     const volScalarField& p,
@@ -590,8 +586,7 @@ Foam::tmp<Foam::volScalarField> Foam::heThermo<BasicThermo, MixtureType>::THE
     (
         "T",
         dimTemperature,
-        &MixtureType::cellThermoMixture,
-        &MixtureType::patchFaceThermoMixture,
+        &MixtureType::thermoMixture,
         &MixtureType::thermoMixtureType::THE,
         h,
         p,
@@ -600,8 +595,8 @@ Foam::tmp<Foam::volScalarField> Foam::heThermo<BasicThermo, MixtureType>::THE
 }
 
 
-template<class BasicThermo, class MixtureType>
-Foam::tmp<Foam::scalarField> Foam::heThermo<BasicThermo, MixtureType>::THE
+template<class MixtureType, class BasicThermoType>
+Foam::tmp<Foam::scalarField> Foam::heThermo<MixtureType, BasicThermoType>::THE
 (
     const scalarField& h,
     const scalarField& T0,
@@ -610,7 +605,7 @@ Foam::tmp<Foam::scalarField> Foam::heThermo<BasicThermo, MixtureType>::THE
 {
     return cellSetProperty
     (
-        &MixtureType::cellThermoMixture,
+        &MixtureType::thermoMixture,
         &MixtureType::thermoMixtureType::THE,
         cells,
         h,
@@ -620,8 +615,8 @@ Foam::tmp<Foam::scalarField> Foam::heThermo<BasicThermo, MixtureType>::THE
 }
 
 
-template<class BasicThermo, class MixtureType>
-Foam::tmp<Foam::scalarField> Foam::heThermo<BasicThermo, MixtureType>::THE
+template<class MixtureType, class BasicThermoType>
+Foam::tmp<Foam::scalarField> Foam::heThermo<MixtureType, BasicThermoType>::THE
 (
     const scalarField& h,
     const scalarField& T0,
@@ -630,7 +625,7 @@ Foam::tmp<Foam::scalarField> Foam::heThermo<BasicThermo, MixtureType>::THE
 {
     return patchFieldProperty
     (
-        &MixtureType::patchFaceThermoMixture,
+        &MixtureType::thermoMixture,
         &MixtureType::thermoMixtureType::THE,
         patchi,
         h,
@@ -640,42 +635,42 @@ Foam::tmp<Foam::scalarField> Foam::heThermo<BasicThermo, MixtureType>::THE
 }
 
 
-template<class BasicThermo, class MixtureType>
+template<class MixtureType, class BasicThermoType>
 Foam::tmp<Foam::volScalarField>
-Foam::heThermo<BasicThermo, MixtureType>::W() const
+Foam::heThermo<MixtureType, BasicThermoType>::W() const
 {
     return volScalarFieldProperty
     (
         "W",
         dimMass/dimMoles,
-        &MixtureType::cellThermoMixture,
-        &MixtureType::patchFaceThermoMixture,
+        &MixtureType::thermoMixture,
         &MixtureType::thermoMixtureType::W
     );
 }
 
 
-template<class BasicThermo, class MixtureType>
-Foam::tmp<Foam::scalarField> Foam::heThermo<BasicThermo, MixtureType>::W
+template<class MixtureType, class BasicThermoType>
+Foam::tmp<Foam::scalarField> Foam::heThermo<MixtureType, BasicThermoType>::W
 (
     const label patchi
 ) const
 {
     return patchFieldProperty
     (
-        &MixtureType::patchFaceThermoMixture,
+        &MixtureType::thermoMixture,
         &MixtureType::thermoMixtureType::W,
         patchi
     );
 }
 
 
-template<class BasicThermo, class MixtureType>
-bool Foam::heThermo<BasicThermo, MixtureType>::read()
+template<class MixtureType, class BasicThermoType>
+bool Foam::heThermo<MixtureType, BasicThermoType>::read()
 {
-    if (BasicThermo::read())
+    if (physicalProperties::read())
     {
         MixtureType::read(*this);
+        BasicThermoType::read(*this);
         return true;
     }
     else

src/thermophysicalModels/basic/heThermo/heThermo.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2011-2022 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2023 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -25,7 +25,9 @@ Class
     Foam::heThermo
 
 Description
-    Enthalpy/Internal energy for a mixture
+    Thermo implementation and storage of energy and heat capacities. Provides
+    overloads of the functions defined in the basic thermo type that depend on
+    the primitive thermo model.
 
 SourceFiles
     heThermo.C
@@ -35,8 +37,8 @@ SourceFiles
 #ifndef heThermo_H
 #define heThermo_H
 
-#include "basicMixture.H"
 #include "volFields.H"
+#include "physicalProperties.H"
 #include "uniformGeometricFields.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -44,15 +46,24 @@ SourceFiles
 namespace Foam
 {
 
+/*---------------------------------------------------------------------------*\
+                       Class heThermoName Declaration
+\*---------------------------------------------------------------------------*/
+
+TemplateName(heThermo);
+
+
 /*---------------------------------------------------------------------------*\
                           Class heThermo Declaration
 \*---------------------------------------------------------------------------*/
 
-template<class BasicThermo, class MixtureType>
+template<class MixtureType, class BasicThermoType>
 class heThermo
 :
-    public BasicThermo,
-    public MixtureType
+    public heThermoName,
+    public physicalProperties,
+    public MixtureType,
+    public BasicThermoType
 {
 protected:
 
@@ -71,38 +82,31 @@ protected:
     // Protected Member Functions
 
         //- Return a volScalarField of the given property
-        template
-        <
-            class CellMixture,
-            class PatchFaceMixture,
-            class Method,
-            class ... Args
-        >
+        template<class Mixture, class Method, class ... Args>
         tmp<volScalarField> volScalarFieldProperty
         (
             const word& psiName,
             const dimensionSet& psiDim,
-            CellMixture cellMixture,
-            PatchFaceMixture patchFaceMixture,
+            Mixture mixture,
             Method psiMethod,
             const Args& ... args
         ) const;
 
         //- Return a scalarField of the given property on a cell set
-        template<class CellMixture, class Method, class ... Args>
+        template<class Mixture, class Method, class ... Args>
         tmp<scalarField> cellSetProperty
         (
-            CellMixture cellMixture,
+            Mixture mixture,
             Method psiMethod,
             const labelList& cells,
             const Args& ... args
         ) const;
 
         //- Return a scalarField of the given property on a patch
-        template<class PatchFaceMixture, class Method, class ... Args>
+        template<class Mixture, class Method, class ... Args>
         tmp<scalarField> patchFieldProperty
         (
-            PatchFaceMixture patchFaceMixture,
+            Mixture mixture,
             Method psiMethod,
             const label patchi,
             const Args& ... args
@@ -128,17 +132,26 @@ protected:
 
 public:
 
+    // Typedefs
+
+        //- Mixture type
+        typedef MixtureType mixtureType;
+
+        //- Basic thermo
+        typedef BasicThermoType basicThermoType;
+
+
+    //- Disambiguate debug switch used by derivations
+    using heThermoName::debug;
+
+
     // Constructors
 
-        //- Construct from mesh
-        heThermo
-        (
-            const fvMesh&,
-            const word& phaseName
-        );
+        //- Construct from mesh and phase name
+        heThermo(const fvMesh&, const word& phaseName);
 
         //- Disallow default bitwise copy construction
-        heThermo(const heThermo<BasicThermo, MixtureType>&) = delete;
+        heThermo(const heThermo<MixtureType, BasicThermoType>&) = delete;
 
 
     //- Destructor
@@ -147,16 +160,20 @@ public:
 
     // Member Functions
 
-        //- Return the composition of the mixture
-        virtual typename MixtureType::basicMixtureType&
-        composition()
+        //- Return the properties dictionary
+        virtual IOdictionary& properties()
         {
             return *this;
         }
 
-        //- Return the composition of the mixture
-        virtual const typename MixtureType::basicMixtureType&
-        composition() const
+        //- Return the properties dictionary
+        virtual const IOdictionary& properties() const
+        {
+            return *this;
+        }
+
+        //- Return the mixture
+        const MixtureType& mixture() const
         {
             return *this;
         }
@@ -357,11 +374,37 @@ public:
 };
 
 
+/*---------------------------------------------------------------------------*\
+                       Class namedHeThermo Declaration
+\*---------------------------------------------------------------------------*/
+
+template<class HeThermo>
+class namedHeThermo
+:
+    public HeThermo
+{
+public:
+
+    //- Runtime type information
+    TypeName("heThermo");
+
+
+    // Constructors
+
+        //- Inherit constructors
+        using HeThermo::HeThermo;
+
+
+    //- Destructor
+    virtual ~namedHeThermo();
+};
+
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 } // End namespace Foam
 
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 #ifdef NoRepository
     #include "heThermo.C"

src/thermophysicalModels/basic/heThermo/heThermoName.C

@@ -0,0 +1,36 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     | Website:  https://openfoam.org
+    \\  /    A nd           | Copyright (C) 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 "heThermo.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    defineTypeNameAndDebug(heThermoName, 0);
+}
+
+
+// ************************************************************************* //

src/thermophysicalModels/basic/liquidThermo/heLiquidThermo.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2014-2022 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2023 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -23,49 +23,43 @@ License
 
 \*---------------------------------------------------------------------------*/
 
-#include "basicCombustionMixture.H"
-
-// * * * * * * * * * * * * * Static Member Functions * * * * * * * * * * * * //
-
-namespace Foam
-{
-    defineTypeNameAndDebug(basicCombustionMixture, 0);
-}
-
+#include "heLiquidThermo.H"
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
-Foam::basicCombustionMixture::basicCombustionMixture
+template<class HeThermo>
+Foam::heLiquidThermo<HeThermo>::heLiquidThermo
 (
-    const dictionary& thermoDict,
-    const wordList& specieNames,
     const fvMesh& mesh,
     const word& phaseName
 )
 :
-    basicMixture(thermoDict, mesh, phaseName),
-    species_(specieNames),
-    Y_(species_.size())
+    HeThermo(mesh, phaseName)
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+template<class HeThermo>
+Foam::heLiquidThermo<HeThermo>::~heLiquidThermo()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+template<class HeThermo>
+Foam::tmp<Foam::volScalarField>
+Foam::heLiquidThermo<HeThermo>::sigma() const
 {
-    forAll(species_, i)
-    {
-        Y_.set
-        (
-            i,
-            new volScalarField
-            (
-                IOobject
-                (
-                    IOobject::groupName(species_[i], phaseName),
-                    mesh.time().name(),
-                    mesh,
-                    IOobject::MUST_READ,
-                    IOobject::AUTO_WRITE
-                ),
-                mesh
-            )
-        );
-    }
+    return this->volScalarFieldProperty
+    (
+        "sigma",
+        dimForce/dimLength,
+        &HeThermo::mixtureType::thermoMixture,
+        &HeThermo::mixtureType::thermoMixtureType::sigma,
+        this->p_,
+        this->T_
+    );
 }
 
 

src/thermophysicalModels/basic/liquidThermo/heLiquidThermo.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2017-2018 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2023 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -21,78 +21,77 @@ License
     You should have received a copy of the GNU General Public License
     along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
 
+Class
+    Foam::heLiquidThermo
+
+Description
+    Liquid thermo implementation
+
+SourceFiles
+    heLiquidThermo.C
+
 \*---------------------------------------------------------------------------*/
 
+#ifndef heLiquidThermo_H
+#define heLiquidThermo_H
+
 #include "liquidThermo.H"
-#include "addToRunTimeSelectionTable.H"
+#include "heThermo.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 namespace Foam
 {
 
-/* * * * * * * * * * * * * * * private static data * * * * * * * * * * * * * */
+/*---------------------------------------------------------------------------*\
+                         Class heLiquidThermo Declaration
+\*---------------------------------------------------------------------------*/
 
-defineTemplateTypeNameAndDebugWithName
-(
-    heRhoThermopureMixtureliquidProperties,
-    "heRhoThermo<pureMixture<liquid,sensibleInternalEnergy>>",
-    0
-);
+template<class HeThermo>
+class heLiquidThermo
+:
+    public HeThermo
+{
+public:
 
-addToRunTimeSelectionTable
-(
-    basicThermo,
-    heRhoThermopureMixtureliquidProperties,
-    fvMesh
-);
+    // Constructors
 
-addToRunTimeSelectionTable
-(
-    fluidThermo,
-    heRhoThermopureMixtureliquidProperties,
-    fvMesh
-);
+        //- Construct from mesh and phase name
+        heLiquidThermo(const fvMesh&, const word& phaseName);
 
-addToRunTimeSelectionTable
-(
-    rhoThermo,
-    heRhoThermopureMixtureliquidProperties,
-    fvMesh
-);
+        //- Disallow default bitwise copy construction
+        heLiquidThermo(const heLiquidThermo<HeThermo>&) = delete;
 
 
-defineTemplateTypeNameAndDebugWithName
-(
-    heRhoThermopureMixtureEnthalpyliquidProperties,
-    "heRhoThermo<pureMixture<liquid,sensibleEnthalpy>>",
-    0
-);
+    //- Destructor
+    virtual ~heLiquidThermo();
 
-addToRunTimeSelectionTable
-(
-    basicThermo,
-    heRhoThermopureMixtureEnthalpyliquidProperties,
-    fvMesh
-);
 
-addToRunTimeSelectionTable
-(
-    fluidThermo,
-    heRhoThermopureMixtureEnthalpyliquidProperties,
-    fvMesh
-);
+    // Member Functions
 
-addToRunTimeSelectionTable
-(
-    rhoThermo,
-    heRhoThermopureMixtureEnthalpyliquidProperties,
-    fvMesh
-);
+        //- Surface tension [N/m]
+        virtual tmp<volScalarField> sigma() const;
+
+
+    // Member Operators
+
+        //- Disallow default bitwise assignment
+        void operator=(const heLiquidThermo<HeThermo>&) = delete;
+};
 
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 } // End namespace Foam
 
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+    #include "heLiquidThermo.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
 // ************************************************************************* //

src/thermophysicalModels/basic/liquidThermo/liquidThermo.C

@@ -0,0 +1,55 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     | Website:  https://openfoam.org
+    \\  /    A nd           | Copyright (C) 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 "liquidThermo.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    defineTypeNameAndDebug(liquidThermo, 0);
+    defineRunTimeSelectionTable(liquidThermo, fvMesh);
+}
+
+
+// * * * * * * * * * * * * * * * * Selectors * * * * * * * * * * * * * * * * //
+
+Foam::autoPtr<Foam::liquidThermo> Foam::liquidThermo::New
+(
+    const fvMesh& mesh,
+    const word& phaseName
+)
+{
+    return basicThermo::New<liquidThermo>(mesh, phaseName);
+}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::liquidThermo::~liquidThermo()
+{}
+
+
+// ************************************************************************* //

src/thermophysicalModels/basic/liquidThermo/liquidThermo.H

@@ -0,0 +1,153 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     | Website:  https://openfoam.org
+    \\  /    A nd           | Copyright (C) 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/>.
+
+Class
+    Foam::liquidThermo
+
+Description
+    Base-class for liquid thermodynamic properties.
+
+See also
+    Foam::basicThermo
+
+SourceFiles
+    liquidThermo.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef liquidThermo_H
+#define liquidThermo_H
+
+#include "rhoThermo.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+template<class, class> class heThermo;
+template<class> class heRhoThermo;
+template<class> class heLiquidThermo;
+
+/*---------------------------------------------------------------------------*\
+                          Class liquidThermo Declaration
+\*---------------------------------------------------------------------------*/
+
+class liquidThermo
+:
+    virtual public rhoThermo
+{
+public:
+
+    // Public Classes
+
+        //- Forward declare the composite class
+        class composite;
+
+
+    // Public Typedefs
+
+        //- The derived type
+        template<class MixtureType>
+        using heThermoType =
+            heLiquidThermo<heRhoThermo<heThermo<MixtureType, composite>>>;
+
+
+    //- Runtime type information
+    TypeName("liquidThermo");
+
+
+    //- Declare run-time constructor selection table
+    declareRunTimeSelectionTable
+    (
+        autoPtr,
+        liquidThermo,
+        fvMesh,
+        (const fvMesh& mesh, const word& phaseName),
+        (mesh, phaseName)
+    );
+
+
+    // Selectors
+
+        //- Standard selection based on fvMesh
+        static autoPtr<liquidThermo> New
+        (
+            const fvMesh&,
+            const word& phaseName=word::null
+        );
+
+
+    //- Destructor
+    virtual ~liquidThermo();
+
+
+    // Member Functions
+
+        //- Surface tension [N/m]
+        virtual tmp<volScalarField> sigma() const = 0;
+};
+
+
+/*---------------------------------------------------------------------------*\
+                    Class liquidThermo::composite Declaration
+\*---------------------------------------------------------------------------*/
+
+class liquidThermo::composite
+:
+    public basicThermo::implementation,
+    public pureThermo,
+    public fluidThermo::implementation,
+    public rhoThermo::implementation,
+    public liquidThermo
+{
+public:
+
+    // Constructors
+
+        //- Construct from dictionary, mesh and phase name
+        template<class MixtureType>
+        composite
+        (
+            const dictionary& dict,
+            const MixtureType& mixture,
+            const fvMesh& mesh,
+            const word& phaseName
+        )
+        :
+            basicThermo::implementation(dict, mesh, phaseName),
+            fluidThermo::implementation(dict, mesh, phaseName),
+            rhoThermo::implementation(dict, mesh, phaseName)
+        {}
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/thermophysicalModels/basic/liquidThermo/liquidThermos.C

@@ -0,0 +1,66 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     | Website:  https://openfoam.org
+    \\  /    A nd           | Copyright (C) 2017-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 "liquidPropertiesSelector.H"
+
+#include "sensibleInternalEnergy.H"
+#include "sensibleEnthalpy.H"
+
+#include "pureMixture.H"
+
+#include "thermo.H"
+
+#include "liquidThermo.H"
+#include "heRhoThermo.H"
+#include "heLiquidThermo.H"
+#include "makeThermo.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#define makeLiquidThermo(ThermoPhysics)                                        \
+                                                                               \
+    defineThermo(liquidThermo, pureMixture, ThermoPhysics);                    \
+                                                                               \
+    addThermo(basicThermo, liquidThermo, pureMixture, ThermoPhysics);          \
+    addThermo(fluidThermo, liquidThermo, pureMixture, ThermoPhysics);          \
+    addThermo(rhoThermo, liquidThermo, pureMixture, ThermoPhysics);            \
+    addThermo(liquidThermo, liquidThermo, pureMixture, ThermoPhysics)
+
+namespace Foam
+{
+    typedef
+        species::thermo<liquidPropertiesSelector, sensibleInternalEnergy>
+        liquidSensibleInternalEnergy;
+
+    makeLiquidThermo(liquidSensibleInternalEnergy);
+
+    typedef
+        species::thermo<liquidPropertiesSelector, sensibleEnthalpy>
+        liquidSensibleEnthalpy;
+
+    makeLiquidThermo(liquidSensibleEnthalpy);
+}
+
+// ************************************************************************* //

src/thermophysicalModels/basic/makeFluidThermo.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2017-2020 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2020-2023 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -21,54 +21,33 @@ License
     You should have received a copy of the GNU General Public License
     along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
 
+InClass
+    Foam::fluidThermo
+
+Description
+    Macros for creating basic fluid thermo packages
+
 \*---------------------------------------------------------------------------*/
 
-#include "rhoThermo.H"
-#include "heRhoThermo.H"
-#include "pureMixture.H"
-#include "thermo.H"
-#include "sensibleInternalEnergy.H"
-#include "sensibleEnthalpy.H"
-#include "thermophysicalPropertiesSelector.H"
-#include "liquidProperties.H"
+#ifndef makeFluidThermo_H
+#define makeFluidThermo_H
+
+#include "makeThermo.H"
+#include "fluidThermo.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-namespace Foam
-{
-
-/* * * * * * * * * * * * * * * private static data * * * * * * * * * * * * * */
-
-typedef heRhoThermo
-<
-    rhoThermo::composite,
-    pureMixture
-    <
-        species::thermo
-        <
-            thermophysicalPropertiesSelector<liquidProperties>,
-            sensibleInternalEnergy
-        >
-    >
-> heRhoThermopureMixtureliquidProperties;
-
-
-typedef heRhoThermo
-<
-    rhoThermo::composite,
-    pureMixture
-    <
-        species::thermo
-        <
-            thermophysicalPropertiesSelector<liquidProperties>,
-            sensibleEnthalpy
-        >
-    >
-> heRhoThermopureMixtureEnthalpyliquidProperties;
+#define makeFluidThermo(BaseThermo, Mixture, ThermoPhysics)                    \
+                                                                               \
+    defineThermo(BaseThermo, Mixture, ThermoPhysics);                          \
+                                                                               \
+    addThermo(basicThermo, BaseThermo, Mixture, ThermoPhysics);                \
+    addThermo(fluidThermo, BaseThermo, Mixture, ThermoPhysics);                \
+    addThermo(BaseThermo, BaseThermo, Mixture, ThermoPhysics)
 
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-} // End namespace Foam
+#endif
 
 // ************************************************************************* //

src/thermophysicalModels/basic/mixtures/pureMixture/pureMixture.C

@@ -24,29 +24,22 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "pureMixture.H"
-#include "fvMesh.H"
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 template<class ThermoType>
-Foam::pureMixture<ThermoType>::pureMixture
-(
-    const dictionary& thermoDict,
-    const fvMesh& mesh,
-    const word& phaseName
-)
+Foam::pureMixture<ThermoType>::pureMixture(const dictionary& dict)
 :
-    basicMixture(thermoDict, mesh, phaseName),
-    mixture_("mixture", thermoDict.subDict("mixture"))
+    mixture_("mixture", dict.subDict("mixture"))
 {}
 
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
 template<class ThermoType>
-void Foam::pureMixture<ThermoType>::read(const dictionary& thermoDict)
+void Foam::pureMixture<ThermoType>::read(const dictionary& dict)
 {
-    mixture_ = ThermoType("mixture", thermoDict.subDict("mixture"));
+    mixture_ = ThermoType("mixture", dict.subDict("mixture"));
 }
 
 

src/thermophysicalModels/basic/mixtures/pureMixture/pureMixture.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2011-2020 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2023 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -25,7 +25,8 @@ Class
     Foam::pureMixture
 
 Description
-    Foam::pureMixture
+    Pure mixture model. This does no mixing, it just returns the single
+    underlying thermo model.
 
 SourceFiles
     pureMixture.C
@@ -35,7 +36,7 @@ SourceFiles
 #ifndef pureMixture_H
 #define pureMixture_H
 
-#include "basicMixture.H"
+#include "dictionary.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -48,26 +49,26 @@ namespace Foam
 
 template<class ThermoType>
 class pureMixture
-:
-    public basicMixture
 {
-
 public:
 
-    //- The type of thermodynamics this mixture is instantiated for
-    typedef ThermoType thermoType;
+    // Public Typedefs
 
-    //- Mixing type for thermodynamic properties
-    typedef ThermoType thermoMixtureType;
+        //- The type of thermodynamics this mixture is instantiated for
+        typedef ThermoType thermoType;
 
-    //- Mixing type for transport properties
-    typedef ThermoType transportMixtureType;
+        //- Mixing type for thermodynamic properties
+        typedef ThermoType thermoMixtureType;
+
+        //- Mixing type for transport properties
+        typedef ThermoType transportMixtureType;
 
 
 private:
 
     // Private Data
 
+        //- Thermo model
         ThermoType mixture_;
 
 
@@ -75,8 +76,8 @@ public:
 
     // Constructors
 
-        //- Construct from dictionary, mesh and phase name
-        pureMixture(const dictionary&, const fvMesh&, const word&);
+        //- Construct from a dictionary
+        pureMixture(const dictionary&);
 
         //- Disallow default bitwise copy construction
         pureMixture(const pureMixture<ThermoType>&) = delete;
@@ -90,47 +91,22 @@ public:
             return "pureMixture<" + ThermoType::typeName() + '>';
         }
 
-        const thermoMixtureType& cellThermoMixture(const label) const
+        //- Return the mixture for thermodynamic properties
+        const thermoMixtureType& thermoMixture(const nil) const
         {
             return mixture_;
         }
 
-        const thermoMixtureType& patchFaceThermoMixture
+        //- Return the mixture for transport properties
+        const transportMixtureType& transportMixture(const nil) const
+        {
+            return mixture_;
+        }
+
+        //- Return the mixture for transport properties
+        const transportMixtureType& transportMixture
         (
-            const label,
-            const label
-        ) const
-        {
-            return mixture_;
-        }
-
-        const transportMixtureType& cellTransportMixture(const label) const
-        {
-            return mixture_;
-        }
-
-        const transportMixtureType& patchFaceTransportMixture
-        (
-            const label,
-            const label
-        ) const
-        {
-            return mixture_;
-        }
-
-        const transportMixtureType& cellTransportMixture
-        (
-            const label,
-            const thermoMixtureType&
-        ) const
-        {
-            return mixture_;
-        }
-
-        const transportMixtureType& patchFaceTransportMixture
-        (
-            const label,
-            const label,
+            const nil,
             const thermoMixtureType&
         ) const
         {

src/thermophysicalModels/basic/psiThermo/hePsiThermo.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2011-2022 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2023 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -27,8 +27,8 @@ License
 
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
-template<class BasicPsiThermo, class MixtureType>
-void Foam::hePsiThermo<BasicPsiThermo, MixtureType>::calculate()
+template<class HeThermo>
+void Foam::hePsiThermo<HeThermo>::calculate()
 {
     const scalarField& hCells = this->he_;
     const scalarField& pCells = this->p_;
@@ -42,11 +42,14 @@ void Foam::hePsiThermo<BasicPsiThermo, MixtureType>::calculate()
 
     forAll(TCells, celli)
     {
-        const typename MixtureType::thermoMixtureType& thermoMixture =
-            this->cellThermoMixture(celli);
+        auto composition = this->cellComposition(celli);
 
-        const typename MixtureType::transportMixtureType& transportMixture =
-            this->cellTransportMixture(celli, thermoMixture);
+        const typename HeThermo::mixtureType::thermoMixtureType&
+            thermoMixture = this->thermoMixture(composition);
+
+        const typename HeThermo::mixtureType::transportMixtureType&
+            transportMixture =
+            this->transportMixture(composition, thermoMixture);
 
         TCells[celli] = thermoMixture.THE
         (
@@ -103,13 +106,14 @@ void Foam::hePsiThermo<BasicPsiThermo, MixtureType>::calculate()
         {
             forAll(pT, facei)
             {
-                const typename MixtureType::thermoMixtureType&
-                    thermoMixture = this->patchFaceThermoMixture(patchi, facei);
+                auto composition = this->patchFaceComposition(patchi, facei);
 
-                const typename MixtureType::transportMixtureType&
+                const typename HeThermo::mixtureType::thermoMixtureType&
+                    thermoMixture = this->thermoMixture(composition);
+
+                const typename HeThermo::mixtureType::transportMixtureType&
                     transportMixture =
-                    this->patchFaceTransportMixture
-                    (patchi, facei, thermoMixture);
+                    this->transportMixture(composition, thermoMixture);
 
                 phe[facei] = thermoMixture.HE(pp[facei], pT[facei]);
 
@@ -125,13 +129,14 @@ void Foam::hePsiThermo<BasicPsiThermo, MixtureType>::calculate()
         {
             forAll(pT, facei)
             {
-                const typename MixtureType::thermoMixtureType& thermoMixture =
-                    this->patchFaceThermoMixture(patchi, facei);
+                auto composition = this->patchFaceComposition(patchi, facei);
 
-                const typename MixtureType::transportMixtureType&
+                const typename HeThermo::mixtureType::thermoMixtureType&
+                    thermoMixture = this->thermoMixture(composition);
+
+                const typename HeThermo::mixtureType::transportMixtureType&
                     transportMixture =
-                    this->patchFaceTransportMixture
-                    (patchi, facei, thermoMixture);
+                    this->transportMixture(composition, thermoMixture);
 
                 pT[facei] = thermoMixture.THE(phe[facei], pp[facei], pT[facei]);
 
@@ -149,14 +154,14 @@ void Foam::hePsiThermo<BasicPsiThermo, MixtureType>::calculate()
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
-template<class BasicPsiThermo, class MixtureType>
-Foam::hePsiThermo<BasicPsiThermo, MixtureType>::hePsiThermo
+template<class HeThermo>
+Foam::hePsiThermo<HeThermo>::hePsiThermo
 (
     const fvMesh& mesh,
     const word& phaseName
 )
 :
-    heThermo<BasicPsiThermo, MixtureType>(mesh, phaseName)
+    HeThermo(mesh, phaseName)
 {
     calculate();
 
@@ -167,17 +172,17 @@ Foam::hePsiThermo<BasicPsiThermo, MixtureType>::hePsiThermo
 
 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
 
-template<class BasicPsiThermo, class MixtureType>
-Foam::hePsiThermo<BasicPsiThermo, MixtureType>::~hePsiThermo()
+template<class HeThermo>
+Foam::hePsiThermo<HeThermo>::~hePsiThermo()
 {}
 
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-template<class BasicPsiThermo, class MixtureType>
-void Foam::hePsiThermo<BasicPsiThermo, MixtureType>::correct()
+template<class HeThermo>
+void Foam::hePsiThermo<HeThermo>::correct()
 {
-    if (debug)
+    if (HeThermo::debug)
     {
         InfoInFunction << endl;
     }
@@ -187,7 +192,7 @@ void Foam::hePsiThermo<BasicPsiThermo, MixtureType>::correct()
 
     calculate();
 
-    if (debug)
+    if (HeThermo::debug)
     {
         Info<< "    Finished" << endl;
     }

src/thermophysicalModels/basic/psiThermo/hePsiThermo.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2011-2020 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2023 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -25,7 +25,7 @@ Class
     Foam::hePsiThermo
 
 Description
-    Energy for a mixture based on compressibility
+    Thermo implementation based on compressibility
 
 SourceFiles
     hePsiThermo.C
@@ -47,10 +47,10 @@ namespace Foam
                          Class hePsiThermo Declaration
 \*---------------------------------------------------------------------------*/
 
-template<class BasicPsiThermo, class MixtureType>
+template<class HeThermo>
 class hePsiThermo
 :
-    public heThermo<BasicPsiThermo, MixtureType>
+    public HeThermo
 {
     // Private Member Functions
 
@@ -59,21 +59,13 @@ class hePsiThermo
 
 public:
 
-    //- Runtime type information
-    TypeName("hePsiThermo");
-
-
     // Constructors
 
         //- Construct from mesh and phase name
-        hePsiThermo
-        (
-            const fvMesh&,
-            const word& phaseName
-        );
+        hePsiThermo(const fvMesh&, const word& phaseName);
 
         //- Disallow default bitwise copy construction
-        hePsiThermo(const hePsiThermo<BasicPsiThermo, MixtureType>&) = delete;
+        hePsiThermo(const hePsiThermo<HeThermo>&) = delete;
 
 
     //- Destructor
@@ -89,10 +81,7 @@ public:
     // Member Operators
 
         //- Disallow default bitwise assignment
-        void operator=
-        (
-            const hePsiThermo<BasicPsiThermo, MixtureType>&
-        ) = delete;
+        void operator=(const hePsiThermo<HeThermo>&) = delete;
 };
 
 

src/thermophysicalModels/basic/psiThermo/psiThermo.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2011-2022 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2023 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -33,11 +33,14 @@ namespace Foam
     defineRunTimeSelectionTable(psiThermo, fvMesh);
 }
 
+const Foam::word Foam::psiThermo::heThermoName("hePsiThermo");
+
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 Foam::psiThermo::implementation::implementation
 (
+    const dictionary& dict,
     const fvMesh& mesh,
     const word& phaseName
 )

src/thermophysicalModels/basic/psiThermo/psiThermo.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2011-2022 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2023 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -39,12 +39,16 @@ SourceFiles
 #define psiThermo_H
 
 #include "fluidThermo.H"
+#include "pureThermo.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 namespace Foam
 {
 
+template<class, class> class heThermo;
+template<class> class hePsiThermo;
+
 /*---------------------------------------------------------------------------*\
                           Class psiThermo Declaration
 \*---------------------------------------------------------------------------*/
@@ -64,6 +68,16 @@ public:
         class composite;
 
 
+    // Public Typedefs
+
+        //- The derived type
+        template<class MixtureType>
+        using heThermoType = hePsiThermo<heThermo<MixtureType, composite>>;
+
+        //- The derived name
+        static const word heThermoName;
+
+
     //- Runtime type information
     TypeName("psiThermo");
 
@@ -127,8 +141,8 @@ public:
 
     // Constructors
 
-        //- Construct from mesh and phase name
-        implementation(const fvMesh&, const word& phaseName);
+        //- Construct from dictionary, mesh and phase name
+        implementation(const dictionary&, const fvMesh&, const word&);
 
         //- Disallow default bitwise copy construction
         implementation(const implementation&) = delete;
@@ -167,6 +181,7 @@ public:
 class psiThermo::composite
 :
     public basicThermo::implementation,
+    public pureThermo,
     public fluidThermo::implementation,
     public psiThermo::implementation
 {
@@ -174,16 +189,19 @@ public:
 
     // Constructors
 
-        //- Construct from mesh and phase name
+        //- Construct from dictionary, mesh and phase name
+        template<class MixtureType>
         composite
         (
+            const dictionary& dict,
+            const MixtureType& mixture,
             const fvMesh& mesh,
             const word& phaseName
         )
         :
-            basicThermo::implementation(mesh, phaseName),
-            fluidThermo::implementation(mesh, phaseName),
-            psiThermo::implementation(mesh, phaseName)
+            basicThermo::implementation(dict, mesh, phaseName),
+            fluidThermo::implementation(dict, mesh, phaseName),
+            psiThermo::implementation(dict, mesh, phaseName)
         {}
 };
 

src/thermophysicalModels/basic/psiThermo/psiThermos.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2011-2021 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2023 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -28,13 +28,13 @@ License
 #include "pureMixture.H"
 
 #include "forGases.H"
-#include "makeThermo.H"
+#include "makeFluidThermo.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 namespace Foam
 {
-    forGases(makeThermo, psiThermo, hePsiThermo, pureMixture);
+    forGases(makeFluidThermo, psiThermo, pureMixture);
 }
 
 // ************************************************************************* //

src/thermophysicalModels/basic/pureThermo/pureThermo.C

@@ -0,0 +1,34 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     | Website:  https://openfoam.org
+    \\  /    A nd           | Copyright (C) 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 "pureThermo.H"
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::pureThermo::~pureThermo()
+{}
+
+
+// ************************************************************************* //

src/thermophysicalModels/basic/pureThermo/pureThermo.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2011-2020 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2023 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -22,46 +22,54 @@ License
     along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
 
 Class
-    Foam::basicMixture
+    Foam::pureThermo
 
 Description
-    Foam::basicMixture
+    Base-class for multi-component thermodynamic properties.
+
+See also
+    Foam::basicThermo
 
 SourceFiles
-    basicMixture.C
+    pureThermo.C
 
 \*---------------------------------------------------------------------------*/
 
-#ifndef basicMixture_H
-#define basicMixture_H
+#ifndef pureThermo_H
+#define pureThermo_H
+
+#include "basicThermo.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 namespace Foam
 {
 
-class fvMesh;
-class dictionary;
-class word;
-
 /*---------------------------------------------------------------------------*\
-                        Class basicMixture Declaration
+                     Class pureThermo Declaration
 \*---------------------------------------------------------------------------*/
 
-class basicMixture
+class pureThermo
+:
+    virtual public basicThermo
 {
-
 public:
 
-    //- The base class of the mixture
-    typedef basicMixture basicMixtureType;
+    //- Destructor
+    virtual ~pureThermo();
 
 
-    // Constructors
+    // Member Functions
 
-        //- Construct from dictionary, mesh and phase name
-        basicMixture(const dictionary&, const fvMesh&, const word&)
-        {}
+        //- Get the composition of an internal cell
+        inline nil cellComposition(const label celli) const;
+
+        //- Get the composition of a boundary face
+        inline nil patchFaceComposition
+        (
+            const label patchi,
+            const label facei
+        ) const;
 };
 
 
@@ -71,6 +79,10 @@ public:
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
+#include "pureThermoI.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
 #endif
 
 // ************************************************************************* //

src/thermophysicalModels/basic/pureThermo/pureThermoI.H

@@ -0,0 +1,46 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     | Website:  https://openfoam.org
+    \\  /    A nd           | Copyright (C) 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 "pureThermo.H"
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+inline Foam::nil Foam::pureThermo::cellComposition(const label celli) const
+{
+    return nil();
+}
+
+
+inline Foam::nil Foam::pureThermo::patchFaceComposition
+(
+    const label patchi,
+    const label facei
+) const
+{
+    return nil();
+}
+
+
+// ************************************************************************* //

src/thermophysicalModels/basic/rhoThermo/heRhoThermo.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2011-2022 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2023 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -27,8 +27,8 @@ License
 
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
-template<class BasicRhoThermo, class MixtureType>
-void Foam::heRhoThermo<BasicRhoThermo, MixtureType>::calculate()
+template<class HeThermo>
+void Foam::heRhoThermo<HeThermo>::calculate()
 {
     const scalarField& hCells = this->he();
     const scalarField& pCells = this->p_;
@@ -43,11 +43,14 @@ void Foam::heRhoThermo<BasicRhoThermo, MixtureType>::calculate()
 
     forAll(TCells, celli)
     {
-        const typename MixtureType::thermoMixtureType& thermoMixture =
-            this->cellThermoMixture(celli);
+        auto composition = this->cellComposition(celli);
 
-        const typename MixtureType::transportMixtureType& transportMixture =
-            this->cellTransportMixture(celli, thermoMixture);
+        const typename HeThermo::mixtureType::thermoMixtureType&
+            thermoMixture = this->thermoMixture(composition);
+
+        const typename HeThermo::mixtureType::transportMixtureType&
+            transportMixture =
+            this->transportMixture(composition, thermoMixture);
 
         TCells[celli] = thermoMixture.THE
         (
@@ -109,13 +112,14 @@ void Foam::heRhoThermo<BasicRhoThermo, MixtureType>::calculate()
         {
             forAll(pT, facei)
             {
-                const typename MixtureType::thermoMixtureType& thermoMixture =
-                    this->patchFaceThermoMixture(patchi, facei);
+                auto composition = this->patchFaceComposition(patchi, facei);
 
-                const typename MixtureType::transportMixtureType&
+                const typename HeThermo::mixtureType::thermoMixtureType&
+                    thermoMixture = this->thermoMixture(composition);
+
+                const typename HeThermo::mixtureType::transportMixtureType&
                     transportMixture =
-                    this->patchFaceTransportMixture
-                    (patchi, facei, thermoMixture);
+                    this->transportMixture(composition, thermoMixture);
 
                 phe[facei] = thermoMixture.HE(pp[facei], pT[facei]);
 
@@ -132,13 +136,14 @@ void Foam::heRhoThermo<BasicRhoThermo, MixtureType>::calculate()
         {
             forAll(pT, facei)
             {
-                const typename MixtureType::thermoMixtureType& thermoMixture =
-                    this->patchFaceThermoMixture(patchi, facei);
+                auto composition = this->patchFaceComposition(patchi, facei);
 
-                const typename MixtureType::transportMixtureType&
+                const typename HeThermo::mixtureType::thermoMixtureType&
+                    thermoMixture = this->thermoMixture(composition);
+
+                const typename HeThermo::mixtureType::transportMixtureType&
                     transportMixture =
-                    this->patchFaceTransportMixture
-                    (patchi, facei, thermoMixture);
+                    this->transportMixture(composition, thermoMixture);
 
                 pT[facei] = thermoMixture.THE(phe[facei], pp[facei], pT[facei]);
 
@@ -157,14 +162,14 @@ void Foam::heRhoThermo<BasicRhoThermo, MixtureType>::calculate()
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
-template<class BasicRhoThermo, class MixtureType>
-Foam::heRhoThermo<BasicRhoThermo, MixtureType>::heRhoThermo
+template<class HeThermo>
+Foam::heRhoThermo<HeThermo>::heRhoThermo
 (
     const fvMesh& mesh,
     const word& phaseName
 )
 :
-    heThermo<BasicRhoThermo, MixtureType>(mesh, phaseName)
+    HeThermo(mesh, phaseName)
 {
     calculate();
 }
@@ -172,24 +177,24 @@ Foam::heRhoThermo<BasicRhoThermo, MixtureType>::heRhoThermo
 
 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
 
-template<class BasicRhoThermo, class MixtureType>
-Foam::heRhoThermo<BasicRhoThermo, MixtureType>::~heRhoThermo()
+template<class HeThermo>
+Foam::heRhoThermo<HeThermo>::~heRhoThermo()
 {}
 
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-template<class BasicRhoThermo, class MixtureType>
-void Foam::heRhoThermo<BasicRhoThermo, MixtureType>::correct()
+template<class HeThermo>
+void Foam::heRhoThermo<HeThermo>::correct()
 {
-    if (debug)
+    if (HeThermo::debug)
     {
         InfoInFunction << endl;
     }
 
     calculate();
 
-    if (debug)
+    if (HeThermo::debug)
     {
         Info<< "    Finished" << endl;
     }

src/thermophysicalModels/basic/rhoThermo/heRhoThermo.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2011-2020 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2023 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -25,7 +25,7 @@ Class
     Foam::heRhoThermo
 
 Description
-    Energy for a mixture based on density
+    Thermo implementation based on density
 
 SourceFiles
     heRhoThermo.C
@@ -47,10 +47,10 @@ namespace Foam
                          Class heRhoThermo Declaration
 \*---------------------------------------------------------------------------*/
 
-template<class BasicRhoThermo, class MixtureType>
+template<class HeThermo>
 class heRhoThermo
 :
-    public heThermo<BasicRhoThermo, MixtureType>
+    public HeThermo
 {
     // Private Member Functions
 
@@ -60,21 +60,13 @@ class heRhoThermo
 
 public:
 
-    //- Runtime type information
-    TypeName("heRhoThermo");
-
-
     // Constructors
 
         //- Construct from mesh and phase name
-        heRhoThermo
-        (
-            const fvMesh&,
-            const word& phaseName
-        );
+        heRhoThermo(const fvMesh&, const word& phaseName);
 
         //- Disallow default bitwise copy construction
-        heRhoThermo(const heRhoThermo<BasicRhoThermo, MixtureType>&) = delete;
+        heRhoThermo(const heRhoThermo<HeThermo>&) = delete;
 
 
     //- Destructor
@@ -90,10 +82,7 @@ public:
     // Member Operators
 
         //- Disallow default bitwise assignment
-        void operator=
-        (
-            const heRhoThermo<BasicRhoThermo, MixtureType>&
-        ) = delete;
+        void operator=(const heRhoThermo<HeThermo>&) = delete;
 };
 
 

src/thermophysicalModels/basic/rhoThermo/rhoThermo.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2011-2022 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2023 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -33,11 +33,14 @@ namespace Foam
     defineRunTimeSelectionTable(rhoThermo, fvMesh);
 }
 
+const Foam::word Foam::rhoThermo::heThermoName("heRhoThermo");
+
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 Foam::rhoThermo::implementation::implementation
 (
+    const dictionary& dict,
     const fvMesh& mesh,
     const word& phaseName
 )

src/thermophysicalModels/basic/rhoThermo/rhoThermo.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2011-2022 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2023 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -39,13 +39,16 @@ SourceFiles
 #define rhoThermo_H
 
 #include "fluidThermo.H"
-#include "runTimeSelectionTables.H"
+#include "pureThermo.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 namespace Foam
 {
 
+template<class, class> class heThermo;
+template<class> class heRhoThermo;
+
 /*---------------------------------------------------------------------------*\
                           Class rhoThermo Declaration
 \*---------------------------------------------------------------------------*/
@@ -65,6 +68,16 @@ public:
         class composite;
 
 
+    // Public Typedefs
+
+        //- The derived type
+        template<class MixtureType>
+        using heThermoType = heRhoThermo<heThermo<MixtureType, composite>>;
+
+        //- The derived name
+        static const word heThermoName;
+
+
     //- Runtime type information
     TypeName("rhoThermo");
 
@@ -138,8 +151,8 @@ public:
 
     // Constructors
 
-        //- Construct from mesh and phase name
-        implementation(const fvMesh&, const word& phaseName);
+        //- Construct from dictionary, mesh and phase name
+        implementation(const dictionary&, const fvMesh&, const word&);
 
         //- Disallow default bitwise copy construction
         implementation(const implementation&) = delete;
@@ -185,6 +198,7 @@ public:
 class rhoThermo::composite
 :
     public basicThermo::implementation,
+    public pureThermo,
     public fluidThermo::implementation,
     public rhoThermo::implementation
 {
@@ -192,16 +206,19 @@ public:
 
     // Constructors
 
-        //- Construct from mesh and phase name
+        //- Construct from dictionary, mesh and phase name
+        template<class MixtureType>
         composite
         (
+            const dictionary& dict,
+            const MixtureType& mixture,
             const fvMesh& mesh,
             const word& phaseName
         )
         :
-            basicThermo::implementation(mesh, phaseName),
-            fluidThermo::implementation(mesh, phaseName),
-            rhoThermo::implementation(mesh, phaseName)
+            basicThermo::implementation(dict, mesh, phaseName),
+            fluidThermo::implementation(dict, mesh, phaseName),
+            rhoThermo::implementation(dict, mesh, phaseName)
         {}
 };
 

src/thermophysicalModels/basic/rhoThermo/rhoThermos.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2011-2021 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2023 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -30,15 +30,15 @@ License
 #include "forGases.H"
 #include "forLiquids.H"
 #include "forTabulated.H"
-#include "makeThermo.H"
+#include "makeFluidThermo.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 namespace Foam
 {
-    forGases(makeThermo, rhoThermo, heRhoThermo, pureMixture);
-    forLiquids(makeThermo, rhoThermo, heRhoThermo, pureMixture);
-    forTabulated(makeThermo, rhoThermo, heRhoThermo, pureMixture);
+    forGases(makeFluidThermo, rhoThermo, pureMixture);
+    forLiquids(makeFluidThermo, rhoThermo, pureMixture);
+    forTabulated(makeFluidThermo, rhoThermo, pureMixture);
 }
 
 // ************************************************************************* //

src/thermophysicalModels/chemistryModel/chemistryModel/chemistryModel.C

@@ -45,9 +45,12 @@ Foam::chemistryModel<ThermoType>::chemistryModel
       ? jacobianTypeNames_.read(this->lookup("jacobian"))
       : jacobianType::fast
     ),
-    mixture_(refCast<const multicomponentMixture<ThermoType>>(this->thermo())),
+    mixture_
+    (
+        dynamicCast<const multicomponentMixture<ThermoType>>(this->thermo())
+    ),
     specieThermos_(mixture_.specieThermos()),
-    reactions_(mixture_.species(), specieThermos_, this->mesh(), *this),
+    reactions_(thermo.species(), specieThermos_, this->mesh(), *this),
     RR_(nSpecie_),
     Y_(nSpecie_),
     c_(nSpecie_),
@@ -94,8 +97,6 @@ Foam::chemistryModel<ThermoType>::chemistryModel
     // species should be initialised (by default 'active' is true)
     if (reduction_)
     {
-        const basicSpecieMixture& composition = this->thermo().composition();
-
         forAll(Yvf_, i)
         {
             typeIOobject<volScalarField> header
@@ -110,7 +111,7 @@ Foam::chemistryModel<ThermoType>::chemistryModel
             // and NO_WRITE
             if (!header.headerOk())
             {
-                composition.setInactive(i);
+                this->thermo().setSpecieInactive(i);
             }
         }
     }
@@ -778,21 +779,9 @@ Foam::scalar Foam::chemistryModel<ThermoType>::solve
     mechRed_.update();
     tabulation_.update();
 
-    if (reduction_ && Pstream::parRun())
+    if (reduction_)
     {
-        const basicSpecieMixture& composition = this->thermo().composition();
-
-        List<bool> active(composition.active());
-        Pstream::listCombineGather(active, orEqOp<bool>());
-        Pstream::listCombineScatter(active);
-
-        forAll(active, i)
-        {
-            if (active[i])
-            {
-                composition.setActive(i);
-            }
-        }
+        this->thermo().syncSpeciesActive();
     }
 
     return deltaTMin;

src/thermophysicalModels/chemistryModel/chemistryModel/chemistryModel.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2016-2022 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2016-2023 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -290,15 +290,6 @@ public:
             ) const = 0;
 
 
-        // Mechanism reduction functions
-
-            //- Return true if specie i is active
-            inline bool active(const label i) const;
-
-            //- Set specie i active
-            inline void setActive(const label i);
-
-
     // Member Operators
 
         //- Disallow default bitwise assignment

src/thermophysicalModels/chemistryModel/chemistryModel/chemistryModelI.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2016-2022 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2016-2023 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -64,18 +64,4 @@ Foam::chemistryModel<ThermoType>::RR() const
 }
 
 
-template<class ThermoType>
-inline void Foam::chemistryModel<ThermoType>::setActive(const label i)
-{
-    this->thermo().composition().setActive(i);
-}
-
-
-template<class ThermoType>
-inline bool Foam::chemistryModel<ThermoType>::active(const label i) const
-{
-    return this->thermo().composition().active(i);
-}
-
-
 // ************************************************************************* //

src/thermophysicalModels/chemistryModel/chemistryModel/reduction/DAC/DAC.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2016-2021 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2016-2023 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -90,7 +90,7 @@ Foam::chemistryReductionMethods::DAC<ThermoType>::DAC
     const wordHashSet initSet(this->coeffsDict_.lookup("initialSet"));
     forAllConstIter(wordHashSet, initSet, iter)
     {
-        searchInitSet_.append(chemistry.mixture().species()[iter.key()]);
+        searchInitSet_.append(chemistry.thermo().species()[iter.key()]);
     }
 
     if (this->coeffsDict_.found("automaticSIS"))
@@ -198,7 +198,7 @@ Foam::chemistryReductionMethods::DAC<ThermoType>::DAC
             fuelSpecies_[i] = fuelSpeciesEntry[i].first();
             fuelSpeciesProp_[i] = fuelSpeciesEntry[i].second();
             fuelSpeciesID_[i] =
-                this->chemistry_.mixture().species()[fuelSpecies_[i]];
+                this->chemistry_.thermo().species()[fuelSpecies_[i]];
             scalar curMm =
                 this->chemistry_.specieThermos()[fuelSpeciesID_[i]].W();
             Mmtot += fuelSpeciesProp_[i]/curMm;

src/thermophysicalModels/chemistryModel/chemistryModel/reduction/DRG/DRG.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2016-2021 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2016-2023 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -40,7 +40,7 @@ Foam::chemistryReductionMethods::DRG<ThermoType>::DRG
     const wordHashSet initSet(this->coeffsDict_.lookup("initialSet"));
     forAllConstIter(wordHashSet, initSet, iter)
     {
-        searchInitSet_.append(chemistry.mixture().species()[iter.key()]);
+        searchInitSet_.append(chemistry.thermo().species()[iter.key()]);
     }
 }
 

src/thermophysicalModels/chemistryModel/chemistryModel/reduction/DRGEP/DRGEP.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2016-2021 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2016-2023 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -48,7 +48,7 @@ Foam::chemistryReductionMethods::DRGEP<ThermoType>::DRGEP
     const wordHashSet initSet(this->coeffsDict_.lookup("initialSet"));
     forAllConstIter(wordHashSet, initSet, iter)
     {
-        searchInitSet_.append(chemistry.mixture().species()[iter.key()]);
+        searchInitSet_.append(chemistry.thermo().species()[iter.key()]);
     }
 
     if (this->coeffsDict_.found("NGroupBased"))