updates for reacting parcels

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

@@ -26,6 +26,18 @@ License
 
 #include "ReactingMultiphaseParcel.H"
 
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+template<class ParcelType>
+const Foam::label Foam::ReactingMultiphaseParcel<ParcelType>::GAS(0);
+
+template<class ParcelType>
+const Foam::label Foam::ReactingMultiphaseParcel<ParcelType>::LIQUID(1);
+
+template<class ParcelType>
+const Foam::label Foam::ReactingMultiphaseParcel<ParcelType>::SOLID(2);
+
+
 // * * * * * * * * * * * *  Private Member Functions * * * * * * * * * * * * //
 
 template<class ParcelType>
@@ -34,13 +46,16 @@ Foam::scalar Foam::ReactingMultiphaseParcel<ParcelType>::cpEff
 (
     TrackData& td,
     const scalar p,
-    const scalar T
+    const scalar T,
+    const label idG,
+    const label idL,
+    const label idS
 ) const
 {
     return
-        this->Y_[0]*td.cloud().composition().cp(0, YGas_, p, T)
+        this->Y_[GAS]*td.cloud().composition().cp(idG, YGas_, p, T)
-      + this->Y_[1]*td.cloud().composition().cp(0, YLiquid_, p, T)
+      + this->Y_[LIQUID]*td.cloud().composition().cp(idL, YLiquid_, p, T)
-      + this->Y_[2]*td.cloud().composition().cp(0, YSolid_, p, T);
+      + this->Y_[SOLID]*td.cloud().composition().cp(idS, YSolid_, p, T);
 }
 
 
@@ -50,13 +65,45 @@ Foam::scalar Foam::ReactingMultiphaseParcel<ParcelType>::HEff
 (
     TrackData& td,
     const scalar p,
-    const scalar T
+    const scalar T,
+    const label idG,
+    const label idL,
+    const label idS
 ) const
 {
     return
-        this->Y_[0]*td.cloud().composition().H(0, YGas_, p, T)
+        this->Y_[GAS]*td.cloud().composition().H(idG, YGas_, p, T)
-      + this->Y_[1]*td.cloud().composition().H(0, YLiquid_, p, T)
+      + this->Y_[LIQUID]*td.cloud().composition().H(idL, YLiquid_, p, T)
-      + this->Y_[2]*td.cloud().composition().H(0, YSolid_, p, T);
+      + this->Y_[SOLID]*td.cloud().composition().H(idS, YSolid_, p, T);
+}
+
+
+template<class ParcelType>
+Foam::scalar Foam::ReactingMultiphaseParcel<ParcelType>::updateMassFractions
+(
+    const scalar mass0,
+    const scalarField& dMassGas,
+    const scalarField& dMassLiquid,
+    const scalarField& dMassSolid
+)
+{
+    scalarList& YMix = this->Y_;
+
+    scalar dMassGasTot = sum(dMassGas);
+    scalar dMassLiquidTot = sum(dMassLiquid);
+    scalar dMassSolidTot = sum(dMassSolid);
+
+    this->updateMassFraction(mass0*YMix[GAS], dMassGas, YGas_);
+    this->updateMassFraction(mass0*YMix[LIQUID], dMassLiquid, YLiquid_);
+    this->updateMassFraction(mass0*YMix[SOLID], dMassSolid, YSolid_);
+
+    scalar massNew = mass0 - (dMassGasTot + dMassLiquidTot + dMassSolidTot);
+
+    YMix[GAS] = (mass0*YMix[GAS] - dMassGas)/massNew;
+    YMix[LIQUID] = (mass0*YMix[LIQUID] - dMassLiquid)/massNew;
+    YMix[SOLID] = 1.0 - YMix[GAS] - YMix[LIQUID];
+
+    return massNew;
 }
 
 
@@ -93,6 +140,7 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::calc
     scalarField& YMix = this->Y_;
     const label idG = td.cloud().composition().idGas();
     const label idL = td.cloud().composition().idLiquid();
+    const label idS = td.cloud().composition().idSolid();
 
 
     // Initialise transfer terms
@@ -105,13 +153,16 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::calc
     scalar dhTrans = 0.0;
 
     // Mass transfer due to phase change
-    scalarList dMassPC(td.cloud().gases().size(), 0.0);
+    scalarField dMassPC(YLiquid_.size(), 0.0);
 
     // Mass transfer due to devolatilisation
-    scalarList dMassDV(td.cloud().gases().size(), 0.0);
+    scalarField dMassDV(YGas_.size(), 0.0);
 
     // Change in carrier phase composition due to surface reactions
-    scalarList dMassSRc(td.cloud().gases().size(), 0.0);
+    scalarField dMassSRGas(YGas_.size(), 0.0);
+    scalarField dMassSRLiquid(YLiquid_.size(), 0.0);
+    scalarField dMassSRSolid(YSolid_.size(), 0.0);
+    scalarField dMassSRc(td.cloud().gases().size(), 0.0);
 
 
     // Phase change in liquid phase
@@ -121,26 +172,33 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::calc
     scalar ShPC =
         calcPhaseChange(td, dt, cellI, T0, idL, YMix[idL], YLiquid_, dMassPC);
 
-    // Update particle component mass fractions
-    this->updateMassFraction(mass0, dMassPC, YLiquid_);
-
 
     // Devolatilisation
     // ~~~~~~~~~~~~~~~~
 
     // Return enthalpy source and calc mass transfer due to devolatilisation
-    scalar ShDV = calcDevolatilisation(td, dt, T0, mass0, idG, YMix, dMassDV);
+    scalar ShDV =
+        calcDevolatilisation(td, dt, T0, mass0, idG, YMix, dMassDV);
 
 
     // Surface reactions
     // ~~~~~~~~~~~~~~~~~
 
-    // Mass transfer of volatile components from particle to carrier phase
-    const scalarList dMassMT = dMassPC + dMassDV;
-
     // Return enthalpy source and calc mass transfer(s) due to surface reaction
     scalar ShSR =
-        calcSurfaceReactions(td, dt, cellI, T0, dMassMT, dMassSRc, dhTrans);
+        calcSurfaceReactions
+        (
+            td,
+            dt,
+            cellI,
+            T0,
+            dMassPC + dMassDV, // total volatiles evolved
+            dMassSRGas,
+            dMassSRLiquid,
+            dMassSRSolid,
+            dMassSRc,
+            dhTrans
+        );
 
 
     // Heat transfer
@@ -154,12 +212,20 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::calc
     scalar T1 = calcHeatTransfer(td, dt, cellI, Sh, htc, dhTrans);
 
 
+    // Update component mass fractions
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+    scalarField dMassGas = dMassDV + dMassSRGas;
+    scalarField dMassLiquid = dMassPC + dMassSRLiquid;
+    scalarField dMassSolid = dMassSRSolid;
+
+    scalar mass1 =
+        updateMassFractions(mass0, dMassGas, dMassLiquid, dMassSolid);
+
+
     // Motion
     // ~~~~~~
 
-    // Update mass (not to include cMassSR)
-    scalar mass1 = mass0 - sum(dMassPC) - sum(dMassDV);
-
     // No additional forces
     vector Fx = vector::zero;
 
@@ -175,10 +241,25 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::calc
     if (td.cloud().coupled())
     {
         // Transfer mass lost from particle to carrier mass source
-        forAll(dMassMT, i)
+        forAll(YGas_, i)
         {
-            td.cloud().rhoTrans(i)[cellI] +=
+            label id = td.composition().localToGlobalGasId(GAS, i);
-                np0*(dMassPC[i] + dMassDV[i] + dMassSRc[i]);
+            td.cloud().rhoTrans(id)[cellI] += np0*dMassGas[i];
+        }
+        forAll(YLiquid_, i)
+        {
+            label id = td.composition().localToGlobalGasId(LIQUID, i);
+            td.cloud().rhoTrans(id)[cellI] += np0*dMassLiquid[i];
+        }
+        // No mapping between solid components and carrier phase
+//        forAll(YSolid_, i)
+//        {
+//            label id = td.composition().localToGlobalGasId(SOLID, i);
+//            td.cloud().rhoTrans(id)[cellI] += np0*dMassSolid[i];
+//        }
+        forAll(dMassSRc, i)
+        {
+            td.cloud().rhoTrans(i)[cellI] += np0*dMassSRc[i];
         }
 
         // Update momentum transfer
@@ -203,13 +284,29 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::calc
     {
         td.keepParticle = false;
 
-        // Absorb particle(s) into carrier phase
+        if (td.cloud().coupled())
-        forAll(dMassMT, i)
         {
-            td.cloud().rhoTrans(i)[cellI] += np0*dMassMT[i];
+            // Absorb particle(s) into carrier phase
+            forAll(YGas_, i)
+            {
+                label id = td.composition().localToGlobalGasId(GAS, i);
+                td.cloud().rhoTrans(id)[cellI] += np0*dMassGas[i];
+            }
+            forAll(YLiquid_, i)
+            {
+                label id = td.composition().localToGlobalGasId(LIQUID, i);
+                td.cloud().rhoTrans(id)[cellI] += np0*dMassLiquid[i];
+            }
+            // No mapping between solid components and carrier phase
+//            forAll(YSolid_, i)
+//            {
+//                label id = td.composition().localToGlobalGasId(SOLID, i);
+//                td.cloud().rhoTrans(id)[cellI] += np0*dMassSolid[i];
+//            }
+
+            td.cloud().hTrans()[cellI] += np0*mass1*HEff(td, pc, T1, idG, idL, idS);
+            td.cloud().UTrans()[cellI] += np0*mass1*U1;
         }
-        td.cloud().hTrans()[cellI] += np0*mass1*HEff(td, pc, T1);
-        td.cloud().UTrans()[cellI] += np0*mass1*U1;
     }
 
 
@@ -220,7 +317,7 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::calc
     {
         this->U_ = U1;
         this->T_ = T1;
-        this->cp_ = cpEff(td, pc, T1);
+        this->cp_ = cpEff(td, pc, T1, idG, idL, idS);
 
         // Update particle density or diameter
         if (td.constProps().constantVolume())
@@ -244,8 +341,8 @@ Foam::scalar Foam::ReactingMultiphaseParcel<ParcelType>::calcDevolatilisation
     const scalar T,
     const scalar mass,
     const label idVolatile,
-    scalarField& YMixture,
+    const scalarField& YMixture,
-    scalarList& dMassDV
+    scalarList& dMass
 )
 {
     // Check that model is active, and that the parcel temperature is
@@ -268,15 +365,10 @@ Foam::scalar Foam::ReactingMultiphaseParcel<ParcelType>::calcDevolatilisation
         mass,
         T,
         td.cloud().composition().YMixture0()[idVolatile],
-        YMixture[0],
+        YMixture[GAS],
         canCombust_
     );
 
-    // Update (total) mass fractions
-    YMixture[0] = (YMixture[0]*mass - dMassTot)/(mass - dMassTot);
-    YMixture[1] = YMixture[1]*mass/(mass - dMassTot);
-    YMixture[2] = 1.0 - YMixture[0] - YMixture[1];
-
     // Add to cummulative mass transfer
     forAll (YGas_, i)
     {
@@ -284,7 +376,7 @@ Foam::scalar Foam::ReactingMultiphaseParcel<ParcelType>::calcDevolatilisation
 
         // Volatiles mass transfer
         scalar volatileMass = YGas_[i]*dMassTot;
-        dMassDV[id] += volatileMass;
+        dMass[i] += volatileMass;
     }
 
     td.cloud().addToMassDevolatilisation(dMassTot);
@@ -301,7 +393,10 @@ Foam::scalar Foam::ReactingMultiphaseParcel<ParcelType>::calcSurfaceReactions
     const scalar dt,
     const label cellI,
     const scalar T,
-    const scalarList& dMassMT,
+    const scalarList& dMassVolatile,
+    scalarField& dMassSRGas,
+    scalarField& dMassSRLiquid,
+    scalarField& dMassSRSolid,
     scalarList& dMassSRc,
     scalar& dhTrans
 )
@@ -312,8 +407,7 @@ Foam::scalar Foam::ReactingMultiphaseParcel<ParcelType>::calcSurfaceReactions
         return;
     }
 
-    // Update mass transfer(s)
+    // Update surface reactions
-    // - Also updates Y()'s
     scalar HReaction = td.cloud().surfaceReaction().calculate
     (
         dt,
@@ -324,13 +418,15 @@ Foam::scalar Foam::ReactingMultiphaseParcel<ParcelType>::calcSurfaceReactions
         this->pc_,
         this->rhoc_,
         this->mass(),
-        dMassMT,
         YGas_,
         YLiquid_,
         YSolid_,
         this->Y_,
-        dMassSRc,
+        dMassVolatile,
-        HReaction
+        dMassSRGas,
+        dMassSRLiquid,
+        dMassSRSolid,
+        dMassSRc
     );
 
     // Heat of reaction divided between particle and carrier phase by the

src/lagrangian/intermediate/parcels/Templates/ReactingMultiphaseParcel/ReactingMultiphaseParcel.H

@@ -149,17 +149,49 @@ private:
 
         //- Return the mixture effective specific heat capacity
         template<class TrackData>
-        scalar cpEff(TrackData& td, const scalar p, const scalar T) const;
+        scalar cpEff
+        (
+            TrackData& td,
+            const scalar p,
+            const scalar T,
+            const label idG,
+            const label idL,
+            const label idS
+        ) const;
 
         //- Return the mixture effective enthalpy
         template<class TrackData>
-        scalar HEff(TrackData& td, const scalar p, const scalar T) const;
+        scalar HEff
+        (
+            TrackData& td,
+            const scalar p,
+            const scalar T,
+            const label idG,
+            const label idL,
+            const label idS
+        ) const;
+
+        //- Update the mass fractions (Y, YGas, YLiquid, YSolid)
+        scalar updateMassFractions
+        (
+            const scalar mass0,
+            const scalarField& dMassGas,
+            const scalarField& dMassLiquid,
+            const scalarField& dMassSolid
+        );
 
 
 protected:
 
     // Protected data
 
+        // IDs of phases in parent phase list
+
+            static const label GAS;
+            static const label LIQUID;
+            static const label SOLID;
+
+
         // Parcel properties
 
             //- Mass fractions of gases []
@@ -187,8 +219,8 @@ protected:
             const scalar T,
             const scalar mass,
             const label idVolatile,
-            scalarField& YMixture,
+            const scalarField& YMixture,
-            scalarList& dMassDV
+            scalarList& dMass
         );
 
         //- Calculate surface reactions
@@ -199,7 +231,10 @@ protected:
             const scalar dt,
             const label cellI,
             const scalar T,
-            const scalarList& dMassMT,
+            const scalarList& dMassVolatile,
+            scalarField& dMassSRGas,
+            scalarField& dMassSRLiquid,
+            scalarField& dMassSRSolid,
             scalarList& dMassSRc,
             scalar& dhTrans
         );

src/lagrangian/intermediate/parcels/Templates/ReactingMultiphaseParcel/ReactingMultiphaseParcelIO.C

@@ -69,9 +69,9 @@ Foam::ReactingMultiphaseParcel<ParcelType>::ReactingMultiphaseParcel
 
         // scale the mass fractions
         const scalarList& YMix = this->Y_;
-        YGas_ /= YMix[0] + ROOTVSMALL;
+        YGas_ /= YMix[GAS] + ROOTVSMALL;
-        YLiquid_ /= YMix[1] + ROOTVSMALL;
+        YLiquid_ /= YMix[LIQUID] + ROOTVSMALL;
-        YSolid_ /= YMix[2] + ROOTVSMALL;
+        YSolid_ /= YMix[SOLID] + ROOTVSMALL;
     }
 
     // Check state of Istream
@@ -134,7 +134,7 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::readFields
         forAllIter(typename Cloud<ParcelType>, c, iter)
         {
             ReactingMultiphaseParcel<ParcelType>& p = iter();
-            p.YGas_[j] = YGas[i++]/(p.Y()[0] + ROOTVSMALL);
+            p.YGas_[j] = YGas[i++]/(p.Y()[GAS] + ROOTVSMALL);
         }
     }
     // Populate YLiquid for each parcel
@@ -153,7 +153,7 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::readFields
         forAllIter(typename Cloud<ParcelType>, c, iter)
         {
             ReactingMultiphaseParcel<ParcelType>& p = iter();
-            p.YLiquid_[j] = YLiquid[i++]/(p.Y()[1] + ROOTVSMALL);
+            p.YLiquid_[j] = YLiquid[i++]/(p.Y()[LIQUID] + ROOTVSMALL);
         }
     }
     // Populate YSolid for each parcel
@@ -172,7 +172,7 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::readFields
         forAllIter(typename Cloud<ParcelType>, c, iter)
         {
             ReactingMultiphaseParcel<ParcelType>& p = iter();
-            p.YSolid_[j] = YSolid[i++]/(p.Y()[2] + ROOTVSMALL);
+            p.YSolid_[j] = YSolid[i++]/(p.Y()[SOLID] + ROOTVSMALL);
         }
     }
 }
@@ -211,7 +211,7 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::writeFields
             forAllConstIter(typename Cloud<ParcelType>, c, iter)
             {
                 const ReactingMultiphaseParcel<ParcelType>& p0 = iter();
-                YGas[i++] = p0.YGas()[j]*p0.Y()[0];
+                YGas[i++] = p0.YGas()[j]*p0.Y()[GAS];
             }
 
             YGas.write();
@@ -235,7 +235,7 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::writeFields
             forAllConstIter(typename Cloud<ParcelType>, c, iter)
             {
                 const ReactingMultiphaseParcel<ParcelType>& p0 = iter();
-                YLiquid[i++] = p0.YLiquid()[j]*p0.Y()[1];
+                YLiquid[i++] = p0.YLiquid()[j]*p0.Y()[LIQUID];
             }
 
             YLiquid.write();
@@ -259,7 +259,7 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::writeFields
             forAllConstIter(typename Cloud<ParcelType>, c, iter)
             {
                 const ReactingMultiphaseParcel<ParcelType>& p0 = iter();
-                YSolid[i++] = p0.YSolid()[j]*p0.Y()[2];
+                YSolid[i++] = p0.YSolid()[j]*p0.Y()[SOLID];
             }
 
             YSolid.write();

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

@@ -157,7 +157,8 @@ void Foam::ReactingParcel<ParcelType>::calc
             // Absorb particle(s) into carrier phase
             forAll(Y_, i)
             {
-                td.cloud().rhoTrans(i)[cellI] += np0*mass1*Y_[i];
+                label id = td.composition().localToGlobalGasId(0, i);
+                td.cloud().rhoTrans(id)[cellI] += np0*mass1*Y_[i];
             }
             td.cloud().UTrans()[cellI] += np0*mass1*U1;
             scalar HEff = td.cloud().composition().H(0, Y_, pc_, T1);

src/lagrangian/intermediate/phaseProperties/phaseProperties/phaseProperties.C

@@ -44,7 +44,7 @@ const Foam::NamedEnum<Foam::phaseProperties::phaseType, 4>
 
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
-void Foam::phaseProperties::setGlobalGasIds
+void Foam::phaseProperties::setGlobalIds
 (
     const PtrList<volScalarField>& YGas
 )
@@ -111,6 +111,27 @@ void Foam::phaseProperties::setGlobalIds(const wordList& globalNames)
 }
 
 
+void Foam::phaseProperties::setGlobalGasIds
+(
+    const PtrList<volScalarField>& YGas
+)
+{
+    forAll(names_, i)
+    {
+        forAll (YGas, j)
+        {
+            word specieName = YGas[j].name();
+
+            if (specieName == names_[i])
+            {
+                globalGasIds_[i] = j;
+                break;
+            }
+        }
+    }
+}
+
+
 void Foam::phaseProperties::checkTotalMassFraction() const
 {
     scalar total = 0.0;
@@ -209,17 +230,23 @@ void Foam::phaseProperties::initialiseGlobalIds
     {
         case GAS:
         {
-            setGlobalGasIds(YGas);
+            setGlobalIds(YGas);
+            forAll(globalGasIds_, i)
+            {
+                globalGasIds_[i] = globalIds_[i];
+            }
             break;
         }
         case LIQUID:
         {
             setGlobalIds(liquidNames);
+            setGlobalGasIds(YGas);
             break;
         }
         case SOLID:
         {
             setGlobalIds(solidNames);
+            setGlobalGasIds(YGas);
             break;
         }
         default:
@@ -330,6 +357,12 @@ const Foam::labelList& Foam::phaseProperties::globalIds() const
 }
 
 
+const Foam::labelList& Foam::phaseProperties::globalGasIds() const
+{
+    return globalGasIds_;
+}
+
+
 Foam::label Foam::phaseProperties::id(const word& cmptName) const
 {
     forAll(names_, cmptI)

src/lagrangian/intermediate/phaseProperties/phaseProperties/phaseProperties.H

@@ -88,15 +88,21 @@ private:
         //- Global ids
         labelList globalIds_;
 
+        //- Map to gas global id
+        labelList globalGasIds_;
+
 
     // Private member functions
 
         //- Set global ids - specialisation for carrier gas phases
-        void setGlobalGasIds(const PtrList<volScalarField>& YGas);
+        void setGlobalIds(const PtrList<volScalarField>& YGas);
 
         //- Set global ids - liquid and solid phases
         void setGlobalIds(const wordList& globalNames);
 
+        //- Set global gas ids - attempts to map component names to global gases
+        void setGlobalGasIds(const PtrList<volScalarField>& YGas);
+
         //- Check the total mass fraction
         void checkTotalMassFraction() const;
 
@@ -159,6 +165,9 @@ public:
             //- Return const acccess to the global ids
             const labelList& globalIds() const;
 
+            //- Return const acccess to the map to the gas global ids
+            const labelList& globalGasIds() const;
+
             //- Return the global id of a component in the local list by name
             //  Returns -1 if not found
             label globalId(const word& cmptName) const;

src/lagrangian/intermediate/phaseProperties/phaseProperties/phasePropertiesIO.C

@@ -35,7 +35,8 @@ Foam::phaseProperties::phaseProperties(Istream& is)
     stateLabel_("(unknown)"),
     names_(0),
     Y_(0),
-    globalIds_(0)
+    globalIds_(0),
+    globalGasIds_(0)
 {
     is.check("Foam::phaseProperties::phaseProperties(Istream& is)");
 
@@ -66,6 +67,7 @@ Foam::phaseProperties::phaseProperties(Istream& is)
 
     // initialise global ids to -1
     globalIds_.setSize(nComponents, -1);
+    globalGasIds_.setSize(nComponents, -1);
 
     checkTotalMassFraction();
 }
@@ -110,6 +112,7 @@ Foam::Istream& Foam::operator>>(Istream& is, phaseProperties& pp)
 
     // initialise global ids to -1
     pp.globalIds_.setSize(nComponents, -1);
+    pp.globalGasIds_.setSize(nComponents, -1);
 
     pp.checkTotalMassFraction();
 

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

@@ -240,6 +240,33 @@ Foam::label Foam::CompositionModel<CloudType>::localId
 }
 
 
+template<class CloudType>
+Foam::label Foam::CompositionModel<CloudType>::localToGlobalGaslId
+(
+    const label phaseI,
+    const label id
+) const
+{
+    label gid = phaseProps_[phaseI].globalGasIds()[id];
+
+    if (gid < 0)
+    {
+        FatalErrorIn
+        (
+            "Foam::label Foam::CompositionModel<CloudType>::localToGlobalGasId"
+            "("
+                "const label, "
+                "const label"
+            ") const"
+        )   << "Unable to determine global gas id for phase "
+            << phaseI << " with local id " << id
+            << nl << abort(FatalError);
+    }
+
+    return id;
+}
+
+
 template<class CloudType>
 const Foam::scalarField& Foam::CompositionModel<CloudType>::Y0
 (

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

@@ -187,6 +187,13 @@ public:
                 //- Return local id of component cmptName in phase phaseI
                 label localId(const label phaseI, const word& cmptName) const;
 
+                //- Return global gas id of component given local id
+                label localToGlobalGaslId
+                (
+                    const label phaseI,
+                    const label id
+                ) const;
+
                 //- Return the list of phase phaseI mass fractions
                 const scalarField& Y0(const label phaseI) const;
 

src/lagrangian/intermediate/submodels/ReactingMultiphase/SurfaceReactionModel/NoSurfaceReaction/NoSurfaceReaction.C

@@ -66,8 +66,11 @@ Foam::scalar Foam::NoSurfaceReaction<CloudType>::calculate
     const scalar,
     const scalar,
     const scalar,
+    const scalarField&,
+    const scalarField&,
+    const scalarField&,
+    const scalarList&,
     const scalarList&,
-    scalarField&,
     scalarField&,
     scalarField&,
     scalarField&,

src/lagrangian/intermediate/submodels/ReactingMultiphase/SurfaceReactionModel/NoSurfaceReaction/NoSurfaceReaction.H

@@ -84,11 +84,14 @@ public:
             const scalar pc,
             const scalar rhoc,
             const scalar mass,
-            const scalarList& dMassMT,
+            const scalarField& YGas,
-            scalarField& YGas,
+            const scalarField& YLiquid,
-            scalarField& YLiquid,
+            const scalarField& YSolid,
-            scalarField& YSolid,
+            const scalarList& YMixture,
-            scalarField& YMixture,
+            const scalarList& dMassVolatile,
+            scalarField& dMassGas,
+            scalarField& dMassLiquid,
+            scalarField& dMassSolid,
             scalarList& dMassSRc
         ) const;
 };

src/lagrangian/intermediate/submodels/ReactingMultiphase/SurfaceReactionModel/SurfaceReactionModel/SurfaceReactionModel.H

@@ -143,11 +143,14 @@ public:
             const scalar pc,
             const scalar rhoc,
             const scalar mass,
-            const scalarList& dMassMT,
+            const scalarField& YGas,
-            scalarField& YGas,
+            const scalarField& YLiquid,
-            scalarField& YLiquid,
+            const scalarField& YSolid,
-            scalarField& YSolid,
+            const scalarList& YMixture,
-            scalarField& YMixture,
+            const scalarList& dMassVolatile,
+            scalarField& dMassGas,
+            scalarField& dMassLiquid,
+            scalarField& dMassSolid,
             scalarList& dMassSRc
         ) const = 0;
 };