Merge branch 'master' of github.com:OpenFOAM/OpenFOAM-dev

applications/solvers/combustion/reactingFoam/reactingFoam.C

@@ -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-2021 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -101,15 +101,15 @@ int main(int argc, char *argv[])
 
         Info<< "Time = " << runTime.timeName() << nl << endl;
 
-        if (pimple.frozenFlow())
+        // --- Pressure-velocity PIMPLE corrector loop
+        while (pimple.loop())
         {
-            #include "YEqn.H"
+            if (pimple.frozenFlow())
-            #include "EEqn.H"
+            {
-        }
+                #include "YEqn.H"
-        else
+                #include "EEqn.H"
-        {
+            }
-            // --- Pressure-velocity PIMPLE corrector loop
+            else
-            while (pimple.loop())
             {
                 if (pimple.firstPimpleIter() || moveMeshOuterCorrectors)
                 {

applications/solvers/multiphase/multiphaseEulerFoam/Allwmake

@@ -9,6 +9,7 @@ wmake $targetType interfacialModels
 wmake $targetType interfacialCompositionModels
 wmake $targetType multiphaseCompressibleMomentumTransportModels
 wmake $targetType multiphaseThermophysicalTransportModels
+wmake $targetType multiphaseReactions
 multiphaseEulerFoam/Allwmake $targetType $*
 wmake $targetType functionObjects
 

applications/solvers/multiphase/multiphaseEulerFoam/multiphaseEulerFoam/Make/options

@@ -18,6 +18,7 @@ EXE_LIBS = \
     -leulerianInterfacialCompositionModels \
     -lmultiphaseMomentumTransportModels \
     -lmultiphaseThermophysicalTransportModels \
+    -lmultiphaseReactions \
     -lthermophysicalTransportModels \
     -lfiniteVolume \
     -lfvOptions \

applications/solvers/multiphase/multiphaseEulerFoam/multiphaseEulerFoam/multiphaseEulerFoam.C

@@ -101,25 +101,25 @@ int main(int argc, char *argv[])
         runTime++;
         Info<< "Time = " << runTime.timeName() << nl << endl;
 
-        if (pimple.frozenFlow())
+        // --- Pressure-velocity PIMPLE corrector loop
+        while (pimple.loop())
         {
-             fluid.solve(rAUs, rAUfs);
+            if (pimple.frozenFlow())
-             fluid.correct();
+            {
-             fluid.correctContinuityError();
+                fluid.solve(rAUs, rAUfs);
+                fluid.correct();
+                fluid.correctContinuityError();
 
-             #include "YEqns.H"
+                #include "YEqns.H"
-             #include "EEqns.H"
+                #include "EEqns.H"
-             #include "pEqnComps.H"
+                #include "pEqnComps.H"
 
-             forAll(phases, phasei)
+                forAll(phases, phasei)
-             {
+                {
-                 phases[phasei].divU(-pEqnComps[phasei] & p_rgh);
+                    phases[phasei].divU(-pEqnComps[phasei] & p_rgh);
-             }
+                }
-        }
+            }
-        else
+            else
-        {
-            // --- Pressure-velocity PIMPLE corrector loop
-            while (pimple.loop())
             {
                 if (pimple.firstPimpleIter() || moveMeshOuterCorrectors)
                 {

applications/solvers/multiphase/multiphaseEulerFoam/multiphaseReactions/Make/files

@@ -0,0 +1,3 @@
+phaseSurfaceArrheniusReactionRate/makephaseSurfaceArrheniusReactions.C
+
+LIB = $(FOAM_LIBBIN)/libmultiphaseReactions

applications/solvers/multiphase/multiphaseEulerFoam/multiphaseReactions/Make/options

@@ -0,0 +1,15 @@
+EXE_INC = \
+    -I$(LIB_SRC)/thermophysicalModels/specie/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/reactionThermo/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/chemistryModel/lnInclude \
+    -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I../phaseSystems/lnInclude
+
+LIB_LIBS = \
+    -lspecie \
+    -lfluidThermophysicalModels \
+    -lreactionThermophysicalModels \
+    -lchemistryModel \
+    -lfiniteVolume \
+    -lphaseSystem

applications/solvers/multiphase/multiphaseEulerFoam/multiphaseReactions/phaseSurfaceArrheniusReactionRate/makephaseSurfaceArrheniusReactions.C

@@ -0,0 +1,44 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     | Website:  https://openfoam.org
+    \\  /    A nd           | Copyright (C) 2021 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "makeReaction.H"
+
+#include "phaseSurfaceArrheniusReactionRate.H"
+
+#include "forCommonGases.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    forCommonGases
+    (
+        makeGeneralReaction,
+        IrreversibleReaction,
+        phaseSurfaceArrheniusReactionRate
+    )
+}
+
+// ************************************************************************* //

applications/solvers/multiphase/multiphaseEulerFoam/multiphaseReactions/phaseSurfaceArrheniusReactionRate/phaseSurfaceArrheniusReactionRate.H

@@ -0,0 +1,149 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     | Website:  https://openfoam.org
+    \\  /    A nd           | Copyright (C) 2021 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::phaseSurfaceArrheniusReactionRate
+
+Description
+    A modified Arrhenius reaction rate given by:
+
+        k = (A * T^beta * exp(-Ta/T))*a
+
+    Where a is the phase surface area per unit volume. The name of the phase is
+    specified by the user.
+
+SourceFiles
+    phaseSurfaceArrheniusReactionRateI.H
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef phaseSurfaceArrheniusReactionRate_H
+#define phaseSurfaceArrheniusReactionRate_H
+
+#include "ArrheniusReactionRate.H"
+#include "speciesTable.H"
+#include "volFields.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// Forward declaration of friend functions and operators
+
+class phaseSurfaceArrheniusReactionRate;
+
+Ostream& operator<<(Ostream&, const phaseSurfaceArrheniusReactionRate&);
+
+/*---------------------------------------------------------------------------*\
+                Class phaseSurfaceArrheniusReactionRate Declaration
+\*---------------------------------------------------------------------------*/
+
+class phaseSurfaceArrheniusReactionRate
+:
+    public ArrheniusReactionRate
+{
+    // Private Data
+
+        //- Name of the phase
+        const word phaseName_;
+
+        //- Reference to the object registry
+        const objectRegistry& ob_;
+
+        //- Surface area per unit volume field
+        mutable tmp<volScalarField> tA_;
+
+
+public:
+
+    // Constructors
+
+        //- Construct from dictionary
+        inline phaseSurfaceArrheniusReactionRate
+        (
+            const speciesTable& species,
+            const objectRegistry& ob,
+            const dictionary& dict
+        );
+
+
+    // Member Functions
+
+        //- Return the type name
+        static word type()
+        {
+            return "phaseSurfaceArrhenius";
+        }
+
+        //- Pre-evaluation hook
+        inline void preEvaluate() const;
+
+        //- Post-evaluation hook
+        inline void postEvaluate() const;
+
+        //- Evaluate the rate
+        inline scalar operator()
+        (
+            const scalar p,
+            const scalar T,
+            const scalarField& c,
+            const label li
+        ) const;
+
+        //- Evaluate the derivative
+        inline scalar ddT
+        (
+            const scalar p,
+            const scalar T,
+            const scalarField& c,
+            const label li
+        ) const;
+
+        //- Write to stream
+        inline void write(Ostream& os) const;
+
+
+    // Ostream Operator
+
+        inline friend Ostream& operator<<
+        (
+            Ostream&,
+            const phaseSurfaceArrheniusReactionRate&
+        );
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#include "phaseSurfaceArrheniusReactionRateI.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

applications/solvers/multiphase/multiphaseEulerFoam/multiphaseReactions/phaseSurfaceArrheniusReactionRate/phaseSurfaceArrheniusReactionRateI.H

@@ -0,0 +1,110 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     | Website:  https://openfoam.org
+    \\  /    A nd           | Copyright (C) 2021 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "phaseSurfaceArrheniusReactionRate.H"
+#include "phaseSystem.H"
+#include "diameterModel.H"
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+inline Foam::phaseSurfaceArrheniusReactionRate::
+phaseSurfaceArrheniusReactionRate
+(
+    const speciesTable& species,
+    const objectRegistry& ob,
+    const dictionary& dict
+)
+:
+    ArrheniusReactionRate(species, dict),
+    phaseName_(dict.lookup("phase")),
+    ob_(ob),
+    tA_(nullptr)
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+inline void Foam::phaseSurfaceArrheniusReactionRate::preEvaluate() const
+{
+    ArrheniusReactionRate::preEvaluate();
+
+    const phaseModel& phase =
+        ob_.lookupObject<phaseModel>(IOobject::groupName("alpha", phaseName_));
+
+    tA_ = phase.dPtr()->a();
+}
+
+
+inline void Foam::phaseSurfaceArrheniusReactionRate::postEvaluate() const
+{
+    ArrheniusReactionRate::postEvaluate();
+
+    tA_.clear();
+}
+
+
+inline Foam::scalar Foam::phaseSurfaceArrheniusReactionRate::operator()
+(
+    const scalar p,
+    const scalar T,
+    const scalarField& c,
+    const label li
+) const
+{
+    return ArrheniusReactionRate::operator()(p, T, c, li)*tA_()[li];
+}
+
+
+inline Foam::scalar Foam::phaseSurfaceArrheniusReactionRate::ddT
+(
+    const scalar p,
+    const scalar T,
+    const scalarField& c,
+    const label li
+) const
+{
+    return ArrheniusReactionRate::ddT(p, T, c, li)*tA_()[li];
+}
+
+
+inline void Foam::phaseSurfaceArrheniusReactionRate::write(Ostream& os) const
+{
+    ArrheniusReactionRate::write(os);
+    writeEntry(os, "phase", phaseName_);
+}
+
+
+inline Foam::Ostream& Foam::operator<<
+(
+    Ostream& os,
+    const phaseSurfaceArrheniusReactionRate& arr
+)
+{
+    arr.write(os);
+    return os;
+}
+
+
+// ************************************************************************* //

applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/diameterModels/IATE/IATE.C

@@ -48,17 +48,11 @@ namespace diameterModels
 }
 
 
-// * * * * * * * * * * * * Protected Member Functions * * * * * * * * * * * //
+// * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * * //
 
-Foam::tmp<Foam::volScalarField> Foam::diameterModels::IATE::calcD() const
+Foam::tmp<Foam::volScalarField> Foam::diameterModels::IATE::dsm() const
 {
-    return d_;
+    return max(6/max(kappai_, 6/dMax_), dMin_);
-}
-
-
-Foam::tmp<Foam::volScalarField> Foam::diameterModels::IATE::calcA() const
-{
-    return phase()*kappai_;
 }
 
 
@@ -86,11 +80,9 @@ Foam::diameterModels::IATE::IATE
     dMax_("dMax", dimLength, diameterProperties),
     dMin_("dMin", dimLength, diameterProperties),
     residualAlpha_("residualAlpha", dimless, diameterProperties),
-    d_(dRef()),
+    d_(IOobject::groupName("d", phase.name()), dsm()),
     sources_(diameterProperties.lookup("sources"), IATEsource::iNew(*this))
-{
+{}
-    d_ = dsm();
-}
 
 
 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
@@ -101,13 +93,19 @@ Foam::diameterModels::IATE::~IATE()
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-Foam::tmp<Foam::volScalarField> Foam::diameterModels::IATE::dsm() const
+Foam::tmp<Foam::volScalarField> Foam::diameterModels::IATE::d() const
 {
-    return max(6/max(kappai_, 6/dMax_), dMin_);
+    return d_;
 }
 
 
-void Foam::diameterModels::IATE::correctNoStore()
+Foam::tmp<Foam::volScalarField> Foam::diameterModels::IATE::a() const
+{
+    return phase()*kappai_;
+}
+
+
+void Foam::diameterModels::IATE::correct()
 {
     volScalarField alphaAv
     (

applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/diameterModels/IATE/IATE.H

@@ -83,7 +83,7 @@ class IATE
         dimensionedScalar residualAlpha_;
 
         //- The Sauter-mean diameter of the phase
-        volScalarField& d_;
+        volScalarField d_;
 
         //- IATE sources
         PtrList<IATEsource> sources_;
@@ -94,17 +94,6 @@ class IATE
         tmp<volScalarField> dsm() const;
 
 
-protected:
-
-    // Protected Member Functions
-
-        //- Get the diameter field
-        virtual tmp<volScalarField> calcD() const;
-
-        //- Get the surface area per unit volume field
-        virtual tmp<volScalarField> calcA() const;
-
-
 public:
 
     friend class IATEsource;
@@ -115,7 +104,7 @@ public:
 
     // Constructors
 
-        //- Construct from components
+        //- Construct from dictionary and phase
         IATE
         (
             const dictionary& diameterProperties,
@@ -135,8 +124,14 @@ public:
             return kappai_;
         }
 
-        //- Correct the diameter field
+        //- Get the diameter field
-        virtual void correctNoStore();
+        virtual tmp<volScalarField> d() const;
+
+        //- Get the surface area per unit volume field
+        virtual tmp<volScalarField> a() const;
+
+        //- Correct the model
+        virtual void correct();
 
         //- Read phaseProperties dictionary
         virtual bool read(const dictionary& phaseProperties);

applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/diameterModels/constantDiameter/constantDiameter.C

@@ -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-2021 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -38,19 +38,6 @@ namespace diameterModels
 }
 
 
-// * * * * * * * * * * * * Protected Member Functions * * * * * * * * * * * //
-
-Foam::tmp<Foam::volScalarField> Foam::diameterModels::constant::calcD() const
-{
-    return volScalarField::New
-    (
-        IOobject::groupName("d", phase().name()),
-        phase().mesh(),
-        d_
-    );
-}
-
-
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 Foam::diameterModels::constant::constant
@@ -72,6 +59,17 @@ Foam::diameterModels::constant::~constant()
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
+Foam::tmp<Foam::volScalarField> Foam::diameterModels::constant::d() const
+{
+    return volScalarField::New
+    (
+        IOobject::groupName("d", phase().name()),
+        phase().mesh(),
+        d_
+    );
+}
+
+
 bool Foam::diameterModels::constant::read(const dictionary& phaseProperties)
 {
     spherical::read(phaseProperties);

applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/diameterModels/constantDiameter/constantDiameter.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-2021 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -58,14 +58,6 @@ class constant
         dimensionedScalar d_;
 
 
-protected:
-
-    // Protected Member Functions
-
-        //- Get the diameter field
-        virtual tmp<volScalarField> calcD() const;
-
-
 public:
 
     //- Runtime type information
@@ -74,7 +66,7 @@ public:
 
     // Constructors
 
-        //- Construct from components
+        //- Construct from dictionary and phase
         constant
         (
             const dictionary& diameterProperties,
@@ -88,6 +80,9 @@ public:
 
     // Member Functions
 
+        //- Get the diameter field
+        virtual tmp<volScalarField> d() const;
+
         //- Read diameterProperties dictionary
         virtual bool read(const dictionary& diameterProperties);
 };

applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/diameterModels/diameterModel/diameterModel.C

@@ -34,56 +34,6 @@ namespace Foam
 }
 
 
-// * * * * * * * * * * * * Protected Member Functions * * * * * * * * * * * //
-
-Foam::volScalarField& Foam::diameterModel::dRef()
-{
-    if (!dPtr_.valid())
-    {
-        dPtr_.reset
-        (
-            new volScalarField
-            (
-                IOobject
-                (
-                    IOobject::groupName("d", phase_.name()),
-                    phase_.time().timeName(),
-                    phase_.mesh()
-                ),
-                phase_.mesh(),
-                dimensionedScalar(dimLength, 0)
-            )
-        );
-    }
-
-    return dPtr_();
-}
-
-
-Foam::volScalarField& Foam::diameterModel::aRef()
-{
-    if (!aPtr_.valid())
-    {
-        aPtr_.reset
-        (
-            new volScalarField
-            (
-                IOobject
-                (
-                    IOobject::groupName("a", phase_.name()),
-                    phase_.time().timeName(),
-                    phase_.mesh()
-                ),
-                phase_.mesh(),
-                dimensionedScalar(dimless/dimLength, 0)
-            )
-        );
-    }
-
-    return aPtr_();
-}
-
-
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 Foam::diameterModel::diameterModel
@@ -93,19 +43,8 @@ Foam::diameterModel::diameterModel
 )
 :
     diameterProperties_(diameterProperties),
-    phase_(phase),
+    phase_(phase)
-    dPtr_(nullptr),
+{}
-    aPtr_(nullptr)
-{
-    if (diameterProperties.lookupOrDefault("storeD", false))
-    {
-        dRef();
-    }
-    if (diameterProperties.lookupOrDefault("storeA", false))
-    {
-        aRef();
-    }
-}
 
 
 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
@@ -116,57 +55,8 @@ Foam::diameterModel::~diameterModel()
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-Foam::tmp<Foam::volScalarField> Foam::diameterModel::d() const
-{
-    if (dPtr_.valid())
-    {
-        return dPtr_();
-    }
-    else
-    {
-        return calcD();
-    }
-}
-
-
-Foam::tmp<Foam::volScalarField> Foam::diameterModel::a() const
-{
-    if (aPtr_.valid())
-    {
-        return aPtr_();
-    }
-    else
-    {
-        return calcA();
-    }
-}
-
-
-void Foam::diameterModel::correctNoStore()
-{}
-
-
 void Foam::diameterModel::correct()
-{
+{}
-    correctNoStore();
-
-    if (dPtr_.valid())
-    {
-        tmp<volScalarField> td = calcD();
-        if (td.isTmp())
-        {
-            dPtr_() = td;
-        }
-    }
-    if (aPtr_.valid())
-    {
-        tmp<volScalarField> tA = calcA();
-        if (tA.isTmp())
-        {
-            aPtr_() = tA;
-        }
-    }
-}
 
 
 bool Foam::diameterModel::read(const dictionary& phaseProperties)

applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/diameterModels/diameterModel/diameterModel.H

@@ -59,29 +59,6 @@ class diameterModel
         //- The phase that this model applies
         const phaseModel& phase_;
 
-        //- Optionally stored diameter field
-        autoPtr<volScalarField> dPtr_;
-
-        //- Optionally stored surface area per unit volume field
-        autoPtr<volScalarField> aPtr_;
-
-
-protected:
-
-    // Access
-
-        //- Get a reference to the stored diameter field
-        volScalarField& dRef();
-
-        //- Get a reference to the stored surface area per unit volume field
-        volScalarField& aRef();
-
-        //- Get the diameter field
-        virtual tmp<volScalarField> calcD() const = 0;
-
-        //- Get the surface area per unit volume field
-        virtual tmp<volScalarField> calcA() const = 0;
-
 
 public:
 
@@ -106,6 +83,7 @@ public:
 
     // Constructors
 
+        //- Construct from dictionary and phase
         diameterModel
         (
             const dictionary& diameterProperties,
@@ -119,6 +97,7 @@ public:
 
     // Selectors
 
+        //- Select from dictionary and phase
         static autoPtr<diameterModel> New
         (
             const dictionary& diameterProperties,
@@ -141,16 +120,13 @@ public:
         }
 
         //- Return the diameter
-        tmp<volScalarField> d() const;
+        virtual tmp<volScalarField> d() const = 0;
 
         //- Return the surface area per unit volume
-        tmp<volScalarField> a() const;
+        virtual tmp<volScalarField> a() const = 0;
 
-        //- Correct the diameter field
+        //- Correct the model
-        virtual void correctNoStore();
+        virtual void correct();
-
-        //- Correct the diameter field and optionally store the results
-        void correct();
 
         //- Read phaseProperties dictionary
         virtual bool read(const dictionary& phaseProperties);

applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/diameterModels/isothermalDiameter/isothermalDiameter.C

@@ -38,14 +38,6 @@ namespace diameterModels
 }
 
 
-// * * * * * * * * * * * * Protected Member Functions * * * * * * * * * * * //
-
-Foam::tmp<Foam::volScalarField> Foam::diameterModels::isothermal::calcD() const
-{
-    return d_;
-}
-
-
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 Foam::diameterModels::isothermal::isothermal
@@ -57,10 +49,18 @@ Foam::diameterModels::isothermal::isothermal
     spherical(diameterProperties, phase),
     d0_("d0", dimLength, diameterProperties),
     p0_("p0", dimPressure, diameterProperties),
-    d_(dRef())
+    d_
-{
+    (
-    d_ = d0_;
+        IOobject
-}
+        (
+            IOobject::groupName("d", phase.name()),
+            phase.time().timeName(),
+            phase.mesh()
+        ),
+        phase.mesh(),
+        d0_
+    )
+{}
 
 
 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
@@ -71,7 +71,13 @@ Foam::diameterModels::isothermal::~isothermal()
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-void Foam::diameterModels::isothermal::correctNoStore()
+Foam::tmp<Foam::volScalarField> Foam::diameterModels::isothermal::d() const
+{
+    return d_;
+}
+
+
+void Foam::diameterModels::isothermal::correct()
 {
     const volScalarField& p = phase().db().lookupObject<volScalarField>("p");
 

applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/diameterModels/isothermalDiameter/isothermalDiameter.H

@@ -61,15 +61,7 @@ class isothermal
         dimensionedScalar p0_;
 
         //- Diameter field
-        volScalarField& d_;
+        volScalarField d_;
-
-
-protected:
-
-    // Protected Member Functions
-
-        //- Get the diameter field
-        virtual tmp<volScalarField> calcD() const;
 
 
 public:
@@ -80,7 +72,7 @@ public:
 
     // Constructors
 
-        //- Construct from components
+        //- Construct from dictionary and phase
         isothermal
         (
             const dictionary& diameterProperties,
@@ -94,8 +86,11 @@ public:
 
     // Member Functions
 
-        //- Correct the diameter field
+        //- Get the diameter field
-        virtual void correctNoStore();
+        virtual tmp<volScalarField> d() const;
+
+        //- Correct the model
+        virtual void correct();
 
         //- Read phaseProperties dictionary
         virtual bool read(const dictionary& phaseProperties);

applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/diameterModels/linearTsubDiameter/linearTsubDiameter.C

@@ -40,14 +40,6 @@ namespace diameterModels
 }
 
 
-// * * * * * * * * * * * * Protected Member Functions * * * * * * * * * * * //
-
-Foam::tmp<Foam::volScalarField> Foam::diameterModels::linearTsub::calcD() const
-{
-    return d_;
-}
-
-
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 Foam::diameterModels::linearTsub::linearTsub
@@ -77,10 +69,18 @@ Foam::diameterModels::linearTsub::linearTsub
         dimTemperature,
         diameterProperties.lookupOrDefault("Tsub1", 13.5)
     ),
-    d_(dRef())
+    d_
-{
+    (
-    d_ = d1_;
+        IOobject
-}
+        (
+            IOobject::groupName("d", phase.name()),
+            phase.time().timeName(),
+            phase.mesh()
+        ),
+        phase.mesh(),
+        d1_
+    )
+{}
 
 
 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
@@ -91,7 +91,13 @@ Foam::diameterModels::linearTsub::~linearTsub()
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-void Foam::diameterModels::linearTsub::correctNoStore()
+Foam::tmp<Foam::volScalarField> Foam::diameterModels::linearTsub::d() const
+{
+    return d_;
+}
+
+
+void Foam::diameterModels::linearTsub::correct()
 {
     // Lookup the fluid model
     const phaseSystem& fluid =

applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/diameterModels/linearTsubDiameter/linearTsubDiameter.H

@@ -81,15 +81,7 @@ class linearTsub
         dimensionedScalar Tsub1_;
 
         //- Diameter field
-        volScalarField& d_;
+        volScalarField d_;
-
-
-protected:
-
-    // Protected Member Functions
-
-        //- Get the diameter field
-        virtual tmp<volScalarField> calcD() const;
 
 
 public:
@@ -100,7 +92,7 @@ public:
 
     // Constructors
 
-        //- Construct from components
+        //- Construct from dictionary and phase
         linearTsub
         (
             const dictionary& diameterProperties,
@@ -114,8 +106,11 @@ public:
 
     // Member Functions
 
-        //- Correct the diameter field
+        //- Get the diameter field
-        virtual void correctNoStore();
+        virtual tmp<volScalarField> d() const;
+
+        //- Correct the model
+        virtual void correct();
 
         //- Read phaseProperties dictionary
         virtual bool read(const dictionary& phaseProperties);

applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/diameterModels/sphericalDiameter/sphericalDiameter.C

@@ -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-2021 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -37,15 +37,6 @@ namespace diameterModels
 }
 
 
-// * * * * * * * * * * * * Protected Member Functions * * * * * * * * * * * //
-
-Foam::tmp<Foam::volScalarField>
-Foam::diameterModels::spherical::calcA() const
-{
-    return 6*phase()/d();
-}
-
-
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 Foam::diameterModels::spherical::spherical
@@ -58,6 +49,14 @@ Foam::diameterModels::spherical::spherical
 {}
 
 
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+Foam::tmp<Foam::volScalarField> Foam::diameterModels::spherical::a() const
+{
+    return 6*phase()/d();
+}
+
+
 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
 
 Foam::diameterModels::spherical::~spherical()

applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/diameterModels/sphericalDiameter/sphericalDiameter.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-2021 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -53,14 +53,6 @@ class spherical
 :
     public diameterModel
 {
-protected:
-
-    // Protected Member Functions
-
-        //- Get the diameter field
-        virtual tmp<volScalarField> calcA() const;
-
-
 public:
 
     //- Runtime type information
@@ -69,6 +61,7 @@ public:
 
     // Constructors
 
+        //- Construct from dictionary and phase
         spherical
         (
             const dictionary& diameterProperties,
@@ -78,6 +71,12 @@ public:
 
     //- Destructor
     virtual ~spherical();
+
+
+    // Member Functions
+
+        //- Get the surface area per unit volume field
+        virtual tmp<volScalarField> a() const;
 };
 
 

applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/diameterModels/velocityGroup/velocityGroup.C

@@ -135,41 +135,6 @@ void Foam::diameterModels::velocityGroup::scale()
 }
 
 
-// * * * * * * * * * * * * Protected Member Functions * * * * * * * * * * * //
-
-Foam::tmp<Foam::volScalarField>
-Foam::diameterModels::velocityGroup::calcD() const
-{
-    return d_;
-}
-
-
-Foam::tmp<Foam::volScalarField>
-Foam::diameterModels::velocityGroup::calcA() const
-{
-    tmp<volScalarField> tA
-    (
-        volScalarField::New
-        (
-            "a",
-            phase().mesh(),
-            dimensionedScalar(inv(dimLength), Zero)
-        )
-    );
-
-    volScalarField& a = tA.ref();
-
-    forAll(sizeGroups_, i)
-    {
-        const sizeGroup& fi = sizeGroups_[i];
-
-        a += fi.a()*fi/fi.x();
-    }
-
-    return phase()*a;
-}
-
-
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 Foam::diameterModels::velocityGroup::velocityGroup
@@ -205,10 +170,8 @@ Foam::diameterModels::velocityGroup::velocityGroup
         diameterProperties.lookup("sizeGroups"),
         sizeGroup::iNew(phase, *this)
     ),
-    d_(dRef())
+    d_(IOobject::groupName("d", phase.name()), dsm())
-{
+{}
-    d_ = dsm();
-}
 
 
 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
@@ -219,7 +182,38 @@ Foam::diameterModels::velocityGroup::~velocityGroup()
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-void Foam::diameterModels::velocityGroup::correctNoStore()
+Foam::tmp<Foam::volScalarField> Foam::diameterModels::velocityGroup::d() const
+{
+    return d_;
+}
+
+
+Foam::tmp<Foam::volScalarField> Foam::diameterModels::velocityGroup::a() const
+{
+    tmp<volScalarField> tA
+    (
+        volScalarField::New
+        (
+            "a",
+            phase().mesh(),
+            dimensionedScalar(inv(dimLength), Zero)
+        )
+    );
+
+    volScalarField& a = tA.ref();
+
+    forAll(sizeGroups_, i)
+    {
+        const sizeGroup& fi = sizeGroups_[i];
+
+        a += fi.a()*fi/fi.x();
+    }
+
+    return phase()*a;
+}
+
+
+void Foam::diameterModels::velocityGroup::correct()
 {
     forAll(sizeGroups_, i)
     {

applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/diameterModels/velocityGroup/velocityGroup.H

@@ -108,7 +108,7 @@ class velocityGroup
         PtrList<sizeGroup> sizeGroups_;
 
         //- Sauter-mean diameter of the phase
-        volScalarField& d_;
+        volScalarField d_;
 
 
     // Private Member Functions
@@ -122,17 +122,6 @@ class velocityGroup
         void scale();
 
 
-protected:
-
-    // Protected Member Functions
-
-        //- Get the diameter field
-        virtual tmp<volScalarField> calcD() const;
-
-        //- Get the surface area per unit volume field
-        virtual tmp<volScalarField> calcA() const;
-
-
 public:
 
     //- Runtime type information
@@ -141,7 +130,7 @@ public:
 
     // Constructors
 
-        //- Construct from components
+        //- Construct from dictionary and phase
         velocityGroup
         (
             const dictionary& diameterProperties,
@@ -164,8 +153,14 @@ public:
         //- Return sizeGroups belonging to this velocityGroup
         inline const PtrList<sizeGroup>& sizeGroups() const;
 
-        //- Correct the diameter field
+        //- Get the diameter field
-        virtual void correctNoStore();
+        virtual tmp<volScalarField> d() const;
+
+        //- Get the surface area per unit volume field
+        virtual tmp<volScalarField> a() const;
+
+        //- Correct the model
+        virtual void correct();
 
         //- Read diameterProperties dictionary
         virtual bool read(const dictionary& diameterProperties);

applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/populationBalanceModel/coalescenceModels/PrinceBlanch/PrinceBlanch.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2018-2020 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2018-2021 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -193,8 +193,8 @@ addToCoalescenceRate
     if (laminarShear_)
     {
         coalescenceRate +=
-            1.0/6.0*pow3(fi.d() + fj.d())*shearStrainRate_()
+            pow3(fi.dSph() + fj.dSph())/6
-           *collisionEfficiency;
+           *shearStrainRate_()*collisionEfficiency;
     }
 }
 

applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/populationBalanceModel/coalescenceModels/PrinceBlanch/PrinceBlanch.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2018-2020 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2018-2021 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -62,7 +62,7 @@ Description
 
     \f[
         \theta_{ij}^{LS} =
-            \frac{1}{6} \left(d_i + d_j\right)^{3} \gamma_c\;.
+            \frac{1}{6} \left(d_i + d_j\right)^{3} \dot{\gamma}_c\;.
     \f]
 
     The rise velocity of bubble i is calculated by
@@ -86,12 +86,12 @@ Description
     and the shear strain rate by
 
     \f[
-        \dot{\gamma_{b}} = \mathrm{mag}(\mathrm{symm}(\mathrm{grad}(U_c)))\;.
+        \dot{\gamma}_c = \mathrm{mag}(\mathrm{symm}(\mathrm{grad}(U_c)))\;.
     \f]
 
     Note that in equation 2, the bubble radius has been substituted by the
-    bubble diameter. Also the expression for the equivalent radius r_ij
+    bubble diameter, which leads to a different coefficient in equation 8.
-    (equation 19 in the paper of Prince and Blanch (1990)) was corrected.
+    Also the expression for the equivalent radius, equation 19 was corrected.
 
     \vartable
         \theta_{ij}^{T}   |  Turbulent collision rate [m3/s]
@@ -109,7 +109,7 @@ Description
         u_{ri}            |  Rise velocity of bubble i [m/s]
         S_{ij}            |  Collision cross sectional area [m2]
         g                 |  Gravitational constant [m/s2]
-        \gamma_c          |  Continuous phase shear strain rate [1/s]
+        \dot{\gamma}_c    |  Continuous phase shear strain rate [1/s]
         U_c               |  Continuous phase velocity field [m/s]
     \endvartable
 
@@ -123,7 +123,7 @@ Description
 Usage
     \table
         Property     | Description             | Required    | Default value
-        C1           | coefficient C1          | no          | 0.089
+        C1           | coefficient C1          | no          | 0.356
         h0           | initial film thickness  | no          | 1e-4m
         hf           | critical film thickness | no          | 1e-8m
         turbulence   | Switch for collisions due to turbulence    | yes | none
@@ -160,7 +160,7 @@ class PrinceBlanch
 {
     // Private Data
 
-        //- Optional coefficient C1, defaults to 0.089
+        //- Optional coefficient C1, defaults to 0.356
         dimensionedScalar C1_;
 
         //- Initial film thickness, defaults to 1e-4m

src/OpenFOAM/memory/tmp/tmp.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-2021 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -167,6 +167,9 @@ public:
 
         //- Assignment transferring the temporary T to this tmp
         inline void operator=(const tmp<T>&);
+
+        //- Move assignment transferring the temporary T to this tmp
+        inline void operator=(const tmp<T>&&);
 };
 
 

src/OpenFOAM/memory/tmp/tmpI.H

@@ -369,8 +369,6 @@ inline void Foam::tmp<T>::operator=(const tmp<T>& t)
 
     if (t.isAnyTmp())
     {
-        type_ = t.type_;
-
         // if (!t.ptr_)
         // {
         //     FatalErrorInFunction
@@ -378,7 +376,9 @@ inline void Foam::tmp<T>::operator=(const tmp<T>& t)
         //         << abort(FatalError);
         // }
 
+        type_ = t.type_;
         ptr_ = t.ptr_;
+
         t.ptr_ = 0;
     }
     else
@@ -391,4 +391,19 @@ inline void Foam::tmp<T>::operator=(const tmp<T>& t)
 }
 
 
+template<class T>
+inline void Foam::tmp<T>::operator=(const tmp<T>&& t)
+{
+    clear();
+
+    type_ = t.type_;
+    ptr_ = t.ptr_;
+
+    if (isAnyTmp())
+    {
+        t.ptr_ = 0;
+    }
+}
+
+
 // ************************************************************************* //

src/finiteVolume/cfdTools/general/fvOptions/fvOptionList.C

@@ -80,10 +80,10 @@ void Foam::fv::optionList::checkApplied() const
         {
             const option& source = this->operator[](i);
 
-            wordHashSet notAddedToFields(source.addSupFields());
+            wordHashSet notAddSupFields(source.addSupFields());
-            notAddedToFields -= addedToFields_[i];
+            notAddSupFields -= addSupFields_[i];
 
-            forAllConstIter(wordHashSet, notAddedToFields, iter)
+            forAllConstIter(wordHashSet, notAddSupFields, iter)
             {
                 WarningInFunction
                     << "Source " << source.name() << " defined for field "
@@ -123,7 +123,7 @@ Foam::fv::optionList::optionList(const fvMesh& mesh, const dictionary& dict)
     PtrListDictionary<option>(0),
     mesh_(mesh),
     checkTimeIndex_(mesh_.time().timeIndex() + 1),
-    addedToFields_(),
+    addSupFields_(),
     constrainedFields_(),
     correctedFields_()
 {
@@ -136,7 +136,7 @@ Foam::fv::optionList::optionList(const fvMesh& mesh)
     PtrListDictionary<option>(0),
     mesh_(mesh),
     checkTimeIndex_(mesh_.time().timeIndex() + 1),
-    addedToFields_(),
+    addSupFields_(),
     constrainedFields_(),
     correctedFields_()
 {}
@@ -160,7 +160,7 @@ void Foam::fv::optionList::reset(const dictionary& dict)
 
     this->setSize(count);
 
-    addedToFields_.setSize(count);
+    addSupFields_.setSize(count);
     constrainedFields_.setSize(count);
     correctedFields_.setSize(count);
 
@@ -179,7 +179,7 @@ void Foam::fv::optionList::reset(const dictionary& dict)
                 option::New(name, sourceDict, mesh_).ptr()
             );
 
-            addedToFields_.set(i, new wordHashSet());
+            addSupFields_.set(i, new wordHashSet());
             constrainedFields_.set(i, new wordHashSet());
             correctedFields_.set(i, new wordHashSet());
 

src/finiteVolume/cfdTools/general/fvOptions/fvOptionList.H

@@ -77,7 +77,7 @@ protected:
         mutable label checkTimeIndex_;
 
         //- Sets of the fields that have had sources added by the options
-        mutable PtrList<wordHashSet> addedToFields_;
+        mutable PtrList<wordHashSet> addSupFields_;
 
         //- Sets of the fields that have been constrained by the options
         mutable PtrList<wordHashSet> constrainedFields_;

src/finiteVolume/cfdTools/general/fvOptions/fvOptionListTemplates.C

@@ -52,7 +52,7 @@ Foam::tmp<Foam::fvMatrix<Type>> Foam::fv::optionList::source
 
         if (source.addsSupToField(fieldName))
         {
-            addedToFields_[i].insert(fieldName);
+            addSupFields_[i].insert(fieldName);
 
             if (debug)
             {

src/finiteVolume/cfdTools/general/porosityModel/DarcyForchheimer/DarcyForchheimer.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2012-2020 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2012-2021 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -244,22 +244,6 @@ void Foam::porosityModels::DarcyForchheimer::correct
 }
 
 
-void Foam::porosityModels::DarcyForchheimer::correct
-(
-    fvVectorMatrix& UEqn,
-    const volScalarField& rho,
-    const volScalarField& mu
-) const
-{
-    const vectorField& U = UEqn.psi();
-    const scalarField& V = mesh_.V();
-    scalarField& Udiag = UEqn.diag();
-    vectorField& Usource = UEqn.source();
-
-    apply(Udiag, Usource, V, rho, mu, U);
-}
-
-
 void Foam::porosityModels::DarcyForchheimer::correct
 (
     const fvVectorMatrix& UEqn,

src/finiteVolume/cfdTools/general/porosityModel/DarcyForchheimer/DarcyForchheimer.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2012-2019 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2012-2021 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -160,14 +160,6 @@ public:
         //- Add resistance
         virtual void correct(fvVectorMatrix& UEqn) const;
 
-        //- Add resistance
-        virtual void correct
-        (
-            fvVectorMatrix& UEqn,
-            const volScalarField& rho,
-            const volScalarField& mu
-        ) const;
-
         //- Add resistance
         virtual void correct
         (

src/finiteVolume/cfdTools/general/porosityModel/fixedCoeff/fixedCoeff.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2012-2020 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2012-2021 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -226,28 +226,6 @@ void Foam::porosityModels::fixedCoeff::correct
 }
 
 
-void Foam::porosityModels::fixedCoeff::correct
-(
-    fvVectorMatrix& UEqn,
-    const volScalarField&,
-    const volScalarField&
-) const
-{
-    const vectorField& U = UEqn.psi();
-    const scalarField& V = mesh_.V();
-    scalarField& Udiag = UEqn.diag();
-    vectorField& Usource = UEqn.source();
-
-    scalar rho = 1.0;
-    if (UEqn.dimensions() == dimForce)
-    {
-        coeffs_.lookup("rhoRef") >> rho;
-    }
-
-    apply(Udiag, Usource, V, U, rho);
-}
-
-
 void Foam::porosityModels::fixedCoeff::correct
 (
     const fvVectorMatrix& UEqn,

src/finiteVolume/cfdTools/general/porosityModel/fixedCoeff/fixedCoeff.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2012-2020 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2012-2021 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -137,14 +137,6 @@ public:
         //- Add resistance
         virtual void correct(fvVectorMatrix& UEqn) const;
 
-        //- Add resistance
-        virtual void correct
-        (
-            fvVectorMatrix& UEqn,
-            const volScalarField& rho,
-            const volScalarField& mu
-        ) const;
-
         //- Add resistance
         virtual void correct
         (

src/finiteVolume/cfdTools/general/porosityModel/porosityModel/porosityModel.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2012-2020 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2012-2021 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -150,9 +150,9 @@ Foam::tmp<Foam::vectorField> Foam::porosityModel::porosityModel::force
     const volVectorField& U,
     const volScalarField& rho,
     const volScalarField& mu
-)
+) const
 {
-    transformModelData();
+    const_cast<porosityModel&>(*this).transformModelData();
 
     tmp<vectorField> tforce(new vectorField(U.size(), Zero));
 
@@ -177,23 +177,6 @@ void Foam::porosityModel::addResistance(fvVectorMatrix& UEqn)
 }
 
 
-void Foam::porosityModel::addResistance
-(
-    fvVectorMatrix& UEqn,
-    const volScalarField& rho,
-    const volScalarField& mu
-)
-{
-    if (cellZoneIDs_.empty())
-    {
-        return;
-    }
-
-    transformModelData();
-    this->correct(UEqn, rho, mu);
-}
-
-
 void Foam::porosityModel::addResistance
 (
     const fvVectorMatrix& UEqn,

src/finiteVolume/cfdTools/general/porosityModel/porosityModel/porosityModel.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2012-2020 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2012-2021 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -37,12 +37,11 @@ SourceFiles
 #define porosityModel_H
 
 #include "fvMesh.H"
-#include "dictionary.H"
 #include "fvMatricesFwd.H"
-#include "runTimeSelectionTables.H"
 #include "coordinateSystem.H"
 #include "dimensionedVector.H"
 #include "wordRe.H"
+#include "runTimeSelectionTables.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -94,7 +93,6 @@ protected:
 
     // Protected Member Functions
 
-
         //- Transform the model data wrt mesh changes
         virtual void calcTransformModelData() = 0;
 
@@ -112,13 +110,6 @@ protected:
 
         virtual void correct(fvVectorMatrix& UEqn) const = 0;
 
-        virtual void correct
-        (
-            fvVectorMatrix& UEqn,
-            const volScalarField& rho,
-            const volScalarField& mu
-        ) const = 0;
-
         virtual void correct
         (
             const fvVectorMatrix& UEqn,
@@ -229,19 +220,11 @@ public:
             const volVectorField& U,
             const volScalarField& rho,
             const volScalarField& mu
-        );
+        ) const;
 
         //- Add resistance
         virtual void addResistance(fvVectorMatrix& UEqn);
 
-        //- Add resistance
-        virtual void addResistance
-        (
-            fvVectorMatrix& UEqn,
-            const volScalarField& rho,
-            const volScalarField& mu
-        );
-
         //- Add resistance
         virtual void addResistance
         (

src/finiteVolume/cfdTools/general/porosityModel/porosityModel/porosityModelList.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2012-2020 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2012-2021 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -116,20 +116,6 @@ void Foam::porosityModelList::addResistance
 }
 
 
-void Foam::porosityModelList::addResistance
-(
-    fvVectorMatrix& UEqn,
-    const volScalarField& rho,
-    const volScalarField& mu
-)
-{
-    forAll(*this, i)
-    {
-        this->operator[](i).addResistance(UEqn, rho, mu);
-    }
-}
-
-
 void Foam::porosityModelList::addResistance
 (
     const fvVectorMatrix& UEqn,

src/finiteVolume/cfdTools/general/porosityModel/porosityModel/porosityModelList.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2012-2020 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2012-2021 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -87,14 +87,6 @@ public:
         //- Add resistance
         void addResistance(fvVectorMatrix& UEqn);
 
-        //- Add resistance
-        void addResistance
-        (
-            fvVectorMatrix& UEqn,
-            const volScalarField& rho,
-            const volScalarField& mu
-        );
-
         //- Add resistance
         void addResistance
         (

src/finiteVolume/cfdTools/general/porosityModel/powerLaw/powerLaw.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2012-2019 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2012-2021 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -114,21 +114,6 @@ void Foam::porosityModels::powerLaw::correct
 }
 
 
-void Foam::porosityModels::powerLaw::correct
-(
-    fvVectorMatrix& UEqn,
-    const volScalarField& rho,
-    const volScalarField& mu
-) const
-{
-    const vectorField& U = UEqn.psi();
-    const scalarField& V = mesh_.V();
-    scalarField& Udiag = UEqn.diag();
-
-    apply(Udiag, V, rho, U);
-}
-
-
 void Foam::porosityModels::powerLaw::correct
 (
     const fvVectorMatrix& UEqn,

src/finiteVolume/cfdTools/general/porosityModel/powerLaw/powerLaw.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2012-2020 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2012-2021 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -138,14 +138,6 @@ public:
         //- Add resistance
         virtual void correct(fvVectorMatrix& UEqn) const;
 
-        //- Add resistance
-        virtual void correct
-        (
-            fvVectorMatrix& UEqn,
-            const volScalarField& rho,
-            const volScalarField& mu
-        ) const;
-
         //- Add resistance
         virtual void correct
         (

src/finiteVolume/cfdTools/general/porosityModel/solidification/solidification.C

@@ -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) 2017-2021 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -113,21 +113,6 @@ void Foam::porosityModels::solidification::correct
 }
 
 
-void Foam::porosityModels::solidification::correct
-(
-    fvVectorMatrix& UEqn,
-    const volScalarField& rho,
-    const volScalarField& mu
-) const
-{
-    const volVectorField& U = UEqn.psi();
-    const scalarField& V = mesh_.V();
-    scalarField& Udiag = UEqn.diag();
-
-    apply(Udiag, V, rho, U);
-}
-
-
 void Foam::porosityModels::solidification::correct
 (
     const fvVectorMatrix& UEqn,

src/finiteVolume/cfdTools/general/porosityModel/solidification/solidification.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) 2017-2021 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -201,14 +201,6 @@ public:
         //- Add resistance
         virtual void correct(fvVectorMatrix& UEqn) const;
 
-        //- Add resistance
-        virtual void correct
-        (
-            fvVectorMatrix& UEqn,
-            const volScalarField& rho,
-            const volScalarField& mu
-        ) const;
-
         //- Add resistance
         virtual void correct
         (

src/fvOptions/sources/interRegion/interRegionExplicitPorositySource/interRegionExplicitPorositySource.C

@@ -45,71 +45,12 @@ namespace fv
 }
 }
 
+
 // * * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * //
 
 void Foam::fv::interRegionExplicitPorositySource::readCoeffs()
 {
     UName_ = coeffs_.lookupOrDefault<word>("U", "U");
-
-    muName_ = coeffs_.lookupOrDefault<word>("mu", "thermo:mu");
-}
-
-
-Foam::porosityModel&
-Foam::fv::interRegionExplicitPorositySource::porosity() const
-{
-    if (!porosityPtr_.valid())
-    {
-        const word zoneName(name_ + ":porous");
-
-        const fvMesh& nbrMesh =
-            mesh_.time().lookupObject<fvMesh>(nbrRegionName());
-        const cellZoneMesh& cellZones = nbrMesh.cellZones();
-        label zoneID = cellZones.findZoneID(zoneName);
-
-        if (zoneID == -1)
-        {
-            cellZoneMesh& cz = const_cast<cellZoneMesh&>(cellZones);
-
-            zoneID = cz.size();
-
-            cz.setSize(zoneID + 1);
-
-            cz.set
-            (
-                zoneID,
-                new cellZone
-                (
-                    zoneName,
-                    nbrMesh.faceNeighbour(), // Neighbour internal cells
-                    zoneID,
-                    cellZones
-                )
-            );
-
-            cz.clearAddressing();
-        }
-        else
-        {
-            FatalErrorInFunction
-                << "Unable to create porous cellZone " << zoneName
-                << ": zone already exists"
-                << abort(FatalError);
-        }
-
-        porosityPtr_.reset
-        (
-            porosityModel::New
-            (
-                name_,
-                nbrMesh,
-                coeffs_,
-                zoneName
-            ).ptr()
-        );
-    }
-
-    return porosityPtr_();
 }
 
 
@@ -125,10 +66,83 @@ Foam::fv::interRegionExplicitPorositySource::interRegionExplicitPorositySource
 :
     interRegionOption(name, modelType, dict, mesh),
     UName_(word::null),
-    muName_(word::null),
+    filter_
+    (
+        volScalarField::Internal::New
+        (
+            "filter",
+            mesh,
+            dimensionedScalar(dimless, 0)
+        )
+    ),
     porosityPtr_(nullptr)
 {
     readCoeffs();
+
+    const fvMesh& nbrMesh = mesh_.time().lookupObject<fvMesh>(nbrRegionName());
+
+    meshInterp().mapTgtToSrc
+    (
+        scalarField(nbrMesh.nCells(), 1),
+        plusEqOp<scalar>(),
+        filter_
+    );
+
+    const word zoneName(name_ + ":porous");
+
+    const cellZoneMesh& cellZones = mesh_.cellZones();
+    label zoneID = cellZones.findZoneID(zoneName);
+
+    if (zoneID == -1)
+    {
+        cellZoneMesh& cz = const_cast<cellZoneMesh&>(cellZones);
+
+        zoneID = cz.size();
+
+        cz.setSize(zoneID + 1);
+
+        // Scan the porous region filter for all cells containing porosity
+        labelList porousCells(mesh_.nCells());
+
+        label i = 0;
+        forAll(filter_, celli)
+        {
+            if (filter_[celli] > small)
+            {
+                porousCells[i++] = celli;
+            }
+        }
+        porousCells.setSize(i);
+
+        cz.set
+        (
+            zoneID,
+            new cellZone
+            (
+                zoneName,
+                porousCells,
+                zoneID,
+                cellZones
+            )
+        );
+
+        cz.clearAddressing();
+    }
+    else
+    {
+        FatalErrorInFunction
+            << "Unable to create porous cellZone " << zoneName
+            << ": zone already exists"
+            << abort(FatalError);
+    }
+
+    porosityPtr_ = porosityModel::New
+    (
+        name_,
+        mesh_,
+        coeffs_,
+        zoneName
+    );
 }
 
 
@@ -138,58 +152,6 @@ Foam::wordList
 Foam::fv::interRegionExplicitPorositySource::addSupFields() const
 {
     return wordList(1, UName_);
-
-}
-
-
-void Foam::fv::interRegionExplicitPorositySource::addSup
-(
-    fvMatrix<vector>& eqn,
-    const word& fieldName
-) const
-{
-    const fvMesh& nbrMesh = mesh_.time().lookupObject<fvMesh>(nbrRegionName());
-
-    const volVectorField& U = eqn.psi();
-
-    volVectorField UNbr
-    (
-        IOobject
-        (
-            name_ + ":UNbr",
-            nbrMesh.time().timeName(),
-            nbrMesh,
-            IOobject::NO_READ,
-            IOobject::NO_WRITE
-        ),
-        nbrMesh,
-        dimensionedVector(U.dimensions(), Zero)
-    );
-
-    // Map local velocity onto neighbour region
-    meshInterp().mapSrcToTgt
-    (
-        U.primitiveField(),
-        plusEqOp<vector>(),
-        UNbr.primitiveFieldRef()
-    );
-
-    fvMatrix<vector> nbrEqn(UNbr, eqn.dimensions());
-
-    porosity().addResistance(nbrEqn);
-
-    // Convert source from neighbour to local region
-    fvMatrix<vector> porosityEqn(U, eqn.dimensions());
-    scalarField& Udiag = porosityEqn.diag();
-    vectorField& Usource = porosityEqn.source();
-
-    Udiag.setSize(eqn.diag().size(), 0.0);
-    Usource.setSize(eqn.source().size(), Zero);
-
-    meshInterp().mapTgtToSrc(nbrEqn.diag(), plusEqOp<scalar>(), Udiag);
-    meshInterp().mapTgtToSrc(nbrEqn.source(), plusEqOp<vector>(), Usource);
-
-    eqn -= porosityEqn;
 }
 
 
@@ -200,95 +162,9 @@ void Foam::fv::interRegionExplicitPorositySource::addSup
     const word& fieldName
 ) const
 {
-    const fvMesh& nbrMesh = mesh_.time().lookupObject<fvMesh>(nbrRegionName());
+    fvMatrix<vector> porosityEqn(eqn.psi(), eqn.dimensions());
-
+    porosityPtr_->addResistance(porosityEqn);
-    const volVectorField& U = eqn.psi();
+    eqn -= filter_*porosityEqn;
-
-    volVectorField UNbr
-    (
-        IOobject
-        (
-            name_ + ":UNbr",
-            nbrMesh.time().timeName(),
-            nbrMesh,
-            IOobject::NO_READ,
-            IOobject::NO_WRITE
-        ),
-        nbrMesh,
-        dimensionedVector(U.dimensions(), Zero)
-    );
-
-    // Map local velocity onto neighbour region
-    meshInterp().mapSrcToTgt
-    (
-        U.primitiveField(),
-        plusEqOp<vector>(),
-        UNbr.primitiveFieldRef()
-    );
-
-    fvMatrix<vector> nbrEqn(UNbr, eqn.dimensions());
-
-    volScalarField rhoNbr
-    (
-        IOobject
-        (
-            "rho:UNbr",
-            nbrMesh.time().timeName(),
-            nbrMesh,
-            IOobject::NO_READ,
-            IOobject::NO_WRITE
-        ),
-        nbrMesh,
-        dimensionedScalar(dimDensity, 0)
-    );
-
-    volScalarField muNbr
-    (
-        IOobject
-        (
-            "mu:UNbr",
-            nbrMesh.time().timeName(),
-            nbrMesh,
-            IOobject::NO_READ,
-            IOobject::NO_WRITE
-        ),
-        nbrMesh,
-        dimensionedScalar(dimViscosity, 0)
-    );
-
-    const volScalarField& mu =
-        mesh_.lookupObject<volScalarField>(muName_);
-
-    // Map local rho onto neighbour region
-    meshInterp().mapSrcToTgt
-    (
-        rho.primitiveField(),
-        plusEqOp<scalar>(),
-        rhoNbr.primitiveFieldRef()
-    );
-
-    // Map local mu onto neighbour region
-    meshInterp().mapSrcToTgt
-    (
-        mu.primitiveField(),
-        plusEqOp<scalar>(),
-        muNbr.primitiveFieldRef()
-    );
-
-    porosity().addResistance(nbrEqn, rhoNbr, muNbr);
-
-    // Convert source from neighbour to local region
-    fvMatrix<vector> porosityEqn(U, eqn.dimensions());
-    scalarField& Udiag = porosityEqn.diag();
-    vectorField& Usource = porosityEqn.source();
-
-    Udiag.setSize(eqn.diag().size(), 0.0);
-    Usource.setSize(eqn.source().size(), Zero);
-
-    meshInterp().mapTgtToSrc(nbrEqn.diag(), plusEqOp<scalar>(), Udiag);
-    meshInterp().mapTgtToSrc(nbrEqn.source(), plusEqOp<vector>(), Usource);
-
-    eqn -= porosityEqn;
 }
 
 

src/fvOptions/sources/interRegion/interRegionExplicitPorositySource/interRegionExplicitPorositySource.H

@@ -88,9 +88,8 @@ class interRegionExplicitPorositySource
         //- Velocity field name, default = U
         word UName_;
 
-        //- Dynamic viscosity field name (compressible case only)
+        //- Porosity mapping filter
-        //  default = thermo:mu
+        volScalarField::Internal filter_;
-        word muName_;
 
         //- Run-time selectable porosity model
         mutable autoPtr<porosityModel> porosityPtr_;
@@ -101,9 +100,6 @@ class interRegionExplicitPorositySource
         //- Non-virtual read
         void readCoeffs();
 
-        //- Get the porosity model
-        porosityModel& porosity() const;
-
 
 public:
 
@@ -143,15 +139,6 @@ public:
             virtual wordList addSupFields() const;
 
 
-        // Add explicit and implicit contributions
-
-            virtual void addSup
-            (
-                fvMatrix<vector>& eqn,
-                const word& fieldName
-            ) const;
-
-
         // Add explicit and implicit contributions to compressible equation
 
             virtual void addSup

src/lagrangian/parcel/submodels/Momentum/ParticleForces/Drag/DenseDrag/DenseDragForce.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2019-2020 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2019-2021 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -87,7 +87,7 @@ void Foam::DenseDragForce<CloudType>::cacheFields(const bool store)
 {
     if (store)
     {
-        if (!this->mesh().template foundObject<volVectorField>(alphacName_))
+        if (!this->mesh().template foundObject<volScalarField>(alphacName_))
         {
             alphacPtr_.reset
             (

src/thermophysicalModels/chemistryModel/chemistryModel/StandardChemistryModel/StandardChemistryModel.C

@@ -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-2021 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -107,7 +107,6 @@ void Foam::StandardChemistryModel<ThermoType>::omega
     scalarField& dcdt
 ) const
 {
-
     dcdt = Zero;
 
     forAll(reactions_, i)
@@ -306,6 +305,8 @@ Foam::StandardChemistryModel<ThermoType>::tc() const
 
     if (this->chemistry_)
     {
+        reactionEvaluationScope scope(*this);
+
         forAll(rho, celli)
         {
             const scalar rhoi = rho[celli];
@@ -358,6 +359,8 @@ Foam::StandardChemistryModel<ThermoType>::Qdot() const
 
     if (this->chemistry_)
     {
+        reactionEvaluationScope scope(*this);
+
         scalarField& Qdot = tQdot.ref();
 
         forAll(Y_, i)
@@ -400,6 +403,8 @@ Foam::StandardChemistryModel<ThermoType>::calculateRR
     const scalarField& T = this->thermo().T();
     const scalarField& p = this->thermo().p();
 
+    reactionEvaluationScope scope(*this);
+
     scalar pf, cf, pr, cr;
     label lRef, rRef;
 
@@ -458,6 +463,8 @@ void Foam::StandardChemistryModel<ThermoType>::calculate()
     const scalarField& T = this->thermo().T();
     const scalarField& p = this->thermo().p();
 
+    reactionEvaluationScope scope(*this);
+
     forAll(rho, celli)
     {
         const scalar rhoi = rho[celli];
@@ -502,6 +509,8 @@ Foam::scalar Foam::StandardChemistryModel<ThermoType>::solve
     const scalarField& T = this->thermo().T();
     const scalarField& p = this->thermo().p();
 
+    reactionEvaluationScope scope(*this);
+
     scalarField c0(nSpecie_);
 
     forAll(rho, celli)

src/thermophysicalModels/chemistryModel/chemistryModel/StandardChemistryModel/StandardChemistryModel.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-2021 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -71,6 +71,40 @@ class StandardChemistryModel
 
 protected:
 
+    // Protected classes
+
+        //- Class to define scope of reaction evaluation. Runs pre-evaluate
+        //  hook on all reactions on construction and post-evaluate on
+        //  destruction.
+        class reactionEvaluationScope
+        {
+            const StandardChemistryModel<ThermoType>& chemistry_;
+
+        public:
+
+            reactionEvaluationScope
+            (
+                const StandardChemistryModel<ThermoType>& chemistry
+            )
+            :
+                chemistry_(chemistry)
+            {
+                forAll(chemistry_.reactions_, i)
+                {
+                    chemistry_.reactions_[i].preEvaluate();
+                }
+            }
+
+            ~reactionEvaluationScope()
+            {
+                forAll(chemistry_.reactions_, i)
+                {
+                    chemistry_.reactions_[i].postEvaluate();
+                }
+            }
+        };
+
+
     // Protected data
 
         //- Reference to the field of specie mass fractions

src/thermophysicalModels/chemistryModel/reactions/fluxLimitedLangmuirHinshelwood/fluxLimitedLangmuirHinshelwoodReactionRate.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2019 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2019-2021 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -127,6 +127,12 @@ public:
             return "fluxLimitedLangmuirHinshelwood";
         }
 
+        //- Pre-evaluation hook
+        inline void preEvaluate() const;
+
+        //- Post-evaluation hook
+        inline void postEvaluate() const;
+
         inline scalar operator()
         (
             const scalar p,

src/thermophysicalModels/chemistryModel/reactions/fluxLimitedLangmuirHinshelwood/fluxLimitedLangmuirHinshelwoodReactionRateI.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2019-2020 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2019-2021 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -206,6 +206,16 @@ fluxLimitedLangmuirHinshelwoodReactionRate
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
+inline void
+Foam::fluxLimitedLangmuirHinshelwoodReactionRate::preEvaluate() const
+{}
+
+
+inline void
+Foam::fluxLimitedLangmuirHinshelwoodReactionRate::postEvaluate() const
+{}
+
+
 inline Foam::scalar
 Foam::fluxLimitedLangmuirHinshelwoodReactionRate::operator()
 (

src/thermophysicalModels/chemistryModel/reactions/surfaceArrheniusReactionRate/surfaceArrheniusReactionRate.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2019-2020 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2019-2021 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -92,6 +92,12 @@ public:
             return "surfaceArrhenius";
         }
 
+        //- Pre-evaluation hook
+        inline void preEvaluate() const;
+
+        //- Post-evaluation hook
+        inline void postEvaluate() const;
+
         //- Evaluate the rate
         inline scalar operator()
         (

src/thermophysicalModels/chemistryModel/reactions/surfaceArrheniusReactionRate/surfaceArrheniusReactionRateI.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2019 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2019-2021 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -43,6 +43,18 @@ inline Foam::surfaceArrheniusReactionRate::surfaceArrheniusReactionRate
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
+inline void Foam::surfaceArrheniusReactionRate::preEvaluate() const
+{
+    ArrheniusReactionRate::preEvaluate();
+}
+
+
+inline void Foam::surfaceArrheniusReactionRate::postEvaluate() const
+{
+    ArrheniusReactionRate::postEvaluate();
+}
+
+
 inline Foam::scalar Foam::surfaceArrheniusReactionRate::operator()
 (
     const scalar p,

src/thermophysicalModels/specie/reaction/Reactions/IrreversibleReaction/IrreversibleReaction.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2011-2019 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2021 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -84,6 +84,22 @@ IrreversibleReaction
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
+template<class ReactionThermo, class ReactionRate>
+void
+Foam::IrreversibleReaction<ReactionThermo, ReactionRate>::preEvaluate() const
+{
+    k_.preEvaluate();
+}
+
+
+template<class ReactionThermo, class ReactionRate>
+void
+Foam::IrreversibleReaction<ReactionThermo, ReactionRate>::postEvaluate() const
+{
+    k_.postEvaluate();
+}
+
+
 template<class ReactionThermo, class ReactionRate>
 Foam::scalar Foam::IrreversibleReaction<ReactionThermo, ReactionRate>::kf
 (

src/thermophysicalModels/specie/reaction/Reactions/IrreversibleReaction/IrreversibleReaction.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2011-2019 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2021 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -128,6 +128,15 @@ public:
 
     // Member Functions
 
+        // Hooks
+
+            //- Pre-evaluation hook
+            virtual void preEvaluate() const;
+
+            //- Post-evaluation hook
+            virtual void postEvaluate() const;
+
+
         // IrreversibleReaction rate coefficients
 
             //- Forward rate constant

src/thermophysicalModels/specie/reaction/Reactions/NonEquilibriumReversibleReaction/NonEquilibriumReversibleReaction.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2011-2019 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2021 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -89,6 +89,24 @@ NonEquilibriumReversibleReaction
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
+template<class ReactionThermo, class ReactionRate>
+void Foam::NonEquilibriumReversibleReaction<ReactionThermo, ReactionRate>::
+preEvaluate() const
+{
+    fk_.preEvaluate();
+    rk_.preEvaluate();
+}
+
+
+template<class ReactionThermo, class ReactionRate>
+void Foam::NonEquilibriumReversibleReaction<ReactionThermo, ReactionRate>::
+postEvaluate() const
+{
+    fk_.postEvaluate();
+    rk_.postEvaluate();
+}
+
+
 template<class ReactionThermo, class ReactionRate>
 Foam::scalar
 Foam::NonEquilibriumReversibleReaction<ReactionThermo, ReactionRate>::kf

src/thermophysicalModels/specie/reaction/Reactions/NonEquilibriumReversibleReaction/NonEquilibriumReversibleReaction.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-2021 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -143,6 +143,15 @@ public:
 
     // Member Functions
 
+        // Hooks
+
+            //- Pre-evaluation hook
+            virtual void preEvaluate() const;
+
+            //- Post-evaluation hook
+            virtual void postEvaluate() const;
+
+
         // NonEquilibriumReversibleReaction rate coefficients
 
             //- Forward rate constant

src/thermophysicalModels/specie/reaction/Reactions/Reaction/Reaction.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-2021 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -213,6 +213,15 @@ public:
             inline scalar Thigh() const;
 
 
+        // Hooks
+
+            //- Pre-evaluation hook
+            virtual void preEvaluate() const = 0;
+
+            //- Post-evaluation hook
+            virtual void postEvaluate() const = 0;
+
+
         // Reaction rate coefficients
 
             //- Forward reaction rate

src/thermophysicalModels/specie/reaction/Reactions/ReactionProxy/ReactionProxy.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2018-2019 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2018-2021 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -101,6 +101,16 @@ Foam::ReactionProxy<ReactionThermo>::clone
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
+template<class ReactionThermo>
+void Foam::ReactionProxy<ReactionThermo>::preEvaluate() const
+{}
+
+
+template<class ReactionThermo>
+void Foam::ReactionProxy<ReactionThermo>::postEvaluate() const
+{}
+
+
 template<class ReactionThermo>
 Foam::scalar Foam::ReactionProxy<ReactionThermo>::kf
 (

src/thermophysicalModels/specie/reaction/Reactions/ReactionProxy/ReactionProxy.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2018-2019 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2018-2021 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -100,6 +100,15 @@ public:
 
     // Member Functions
 
+        // Hooks
+
+            //- Pre-evaluation hook
+            virtual void preEvaluate() const;
+
+            //- Post-evaluation hook
+            virtual void postEvaluate() const;
+
+
         // Reaction rate coefficients
 
             //- Forward rate constant

src/thermophysicalModels/specie/reaction/Reactions/ReversibleReaction/ReversibleReaction.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2011-2019 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2021 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -84,6 +84,22 @@ ReversibleReaction
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
+template<class ReactionThermo, class ReactionRate>
+void
+Foam::ReversibleReaction<ReactionThermo, ReactionRate>::preEvaluate() const
+{
+    k_.preEvaluate();
+}
+
+
+template<class ReactionThermo, class ReactionRate>
+void
+Foam::ReversibleReaction<ReactionThermo, ReactionRate>::postEvaluate() const
+{
+    k_.postEvaluate();
+}
+
+
 template<class ReactionThermo, class ReactionRate>
 Foam::scalar Foam::ReversibleReaction<ReactionThermo, ReactionRate>::kf
 (

src/thermophysicalModels/specie/reaction/Reactions/ReversibleReaction/ReversibleReaction.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2011-2019 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2021 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -129,6 +129,15 @@ public:
 
     // Member Functions
 
+        // Hooks
+
+            //- Pre-evaluation hook
+            virtual void preEvaluate() const;
+
+            //- Post-evaluation hook
+            virtual void postEvaluate() const;
+
+
         // ReversibleReaction rate coefficients
 
             //- Forward rate constant

src/thermophysicalModels/specie/reaction/reactionRate/ArrheniusReactionRate/ArrheniusReactionRate.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2011-2019 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2021 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -93,6 +93,12 @@ public:
             return "Arrhenius";
         }
 
+        //- Pre-evaluation hook
+        inline void preEvaluate() const;
+
+        //- Post-evaluation hook
+        inline void postEvaluate() const;
+
         inline scalar operator()
         (
             const scalar p,

src/thermophysicalModels/specie/reaction/reactionRate/ArrheniusReactionRate/ArrheniusReactionRateI.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2011-2019 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2021 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -52,6 +52,14 @@ inline Foam::ArrheniusReactionRate::ArrheniusReactionRate
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
+inline void Foam::ArrheniusReactionRate::preEvaluate() const
+{}
+
+
+inline void Foam::ArrheniusReactionRate::postEvaluate() const
+{}
+
+
 inline Foam::scalar Foam::ArrheniusReactionRate::operator()
 (
     const scalar p,

src/thermophysicalModels/specie/reaction/reactionRate/ChemicallyActivatedReactionRate/ChemicallyActivatedReactionRate.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-2021 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -103,6 +103,12 @@ public:
                 + "ChemicallyActivated";
         }
 
+        //- Pre-evaluation hook
+        inline void preEvaluate() const;
+
+        //- Post-evaluation hook
+        inline void postEvaluate() const;
+
         inline scalar operator()
         (
             const scalar p,

src/thermophysicalModels/specie/reaction/reactionRate/ChemicallyActivatedReactionRate/ChemicallyActivatedReactionRateI.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2011-2019 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2021 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -82,6 +82,30 @@ inline Foam::ChemicallyActivatedReactionRate
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
+template<class ReactionRate, class ChemicallyActivationFunction>
+inline void Foam::ChemicallyActivatedReactionRate
+<
+    ReactionRate,
+    ChemicallyActivationFunction
+>::preEvaluate() const
+{
+    k0_.preEvaluate();
+    kInf_.preEvaluate();
+}
+
+
+template<class ReactionRate, class ChemicallyActivationFunction>
+inline void Foam::ChemicallyActivatedReactionRate
+<
+    ReactionRate,
+    ChemicallyActivationFunction
+>::postEvaluate() const
+{
+    k0_.postEvaluate();
+    kInf_.postEvaluate();
+}
+
+
 template<class ReactionRate, class ChemicallyActivationFunction>
 inline Foam::scalar Foam::ChemicallyActivatedReactionRate
 <

src/thermophysicalModels/specie/reaction/reactionRate/FallOffReactionRate/FallOffReactionRate.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-2021 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -100,6 +100,12 @@ public:
             return ReactionRate::type() + FallOffFunction::type() + "FallOff";
         }
 
+        //- Pre-evaluation hook
+        inline void preEvaluate() const;
+
+        //- Post-evaluation hook
+        inline void postEvaluate() const;
+
         inline scalar operator()
         (
             const scalar p,

src/thermophysicalModels/specie/reaction/reactionRate/FallOffReactionRate/FallOffReactionRateI.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2011-2019 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2021 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -76,6 +76,24 @@ FallOffReactionRate
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
+template<class ReactionRate, class FallOffFunction>
+inline void
+Foam::FallOffReactionRate<ReactionRate, FallOffFunction>::preEvaluate() const
+{
+    k0_.preEvaluate();
+    kInf_.preEvaluate();
+}
+
+
+template<class ReactionRate, class FallOffFunction>
+inline void
+Foam::FallOffReactionRate<ReactionRate, FallOffFunction>::postEvaluate() const
+{
+    k0_.postEvaluate();
+    kInf_.postEvaluate();
+}
+
+
 template<class ReactionRate, class FallOffFunction>
 inline Foam::scalar
 Foam::FallOffReactionRate<ReactionRate, FallOffFunction>::operator()

src/thermophysicalModels/specie/reaction/reactionRate/JanevReactionRate/JanevReactionRate.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-2021 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -96,6 +96,12 @@ public:
             return "Janev";
         }
 
+        //- Pre-evaluation hook
+        inline void preEvaluate() const;
+
+        //- Post-evaluation hook
+        inline void postEvaluate() const;
+
         inline scalar operator()
         (
             const scalar p,

src/thermophysicalModels/specie/reaction/reactionRate/JanevReactionRate/JanevReactionRateI.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2011-2019 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2021 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -55,6 +55,14 @@ inline Foam::JanevReactionRate::JanevReactionRate
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
+inline void Foam::JanevReactionRate::preEvaluate() const
+{}
+
+
+inline void Foam::JanevReactionRate::postEvaluate() const
+{}
+
+
 inline Foam::scalar Foam::JanevReactionRate::operator()
 (
     const scalar p,

src/thermophysicalModels/specie/reaction/reactionRate/LandauTellerReactionRate/LandauTellerReactionRate.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-2021 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -95,6 +95,12 @@ public:
             return "LandauTeller";
         }
 
+        //- Pre-evaluation hook
+        inline void preEvaluate() const;
+
+        //- Post-evaluation hook
+        inline void postEvaluate() const;
+
         inline scalar operator()
         (
             const scalar p,

src/thermophysicalModels/specie/reaction/reactionRate/LandauTellerReactionRate/LandauTellerReactionRateI.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2011-2019 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2021 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -58,6 +58,14 @@ inline Foam::LandauTellerReactionRate::LandauTellerReactionRate
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
+inline void Foam::LandauTellerReactionRate::preEvaluate() const
+{}
+
+
+inline void Foam::LandauTellerReactionRate::postEvaluate() const
+{}
+
+
 inline Foam::scalar Foam::LandauTellerReactionRate::operator()
 (
     const scalar p,

src/thermophysicalModels/specie/reaction/reactionRate/LangmuirHinshelwood/LangmuirHinshelwoodReactionRate.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2011-2019 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2021 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -98,6 +98,12 @@ public:
             return "LangmuirHinshelwood";
         }
 
+        //- Pre-evaluation hook
+        inline void preEvaluate() const;
+
+        //- Post-evaluation hook
+        inline void postEvaluate() const;
+
         inline scalar operator()
         (
             const scalar p,

src/thermophysicalModels/specie/reaction/reactionRate/LangmuirHinshelwood/LangmuirHinshelwoodReactionRateI.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-2021 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -64,6 +64,14 @@ inline Foam::LangmuirHinshelwoodReactionRate::LangmuirHinshelwoodReactionRate
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
+inline void Foam::LangmuirHinshelwoodReactionRate::preEvaluate() const
+{}
+
+
+inline void Foam::LangmuirHinshelwoodReactionRate::postEvaluate() const
+{}
+
+
 inline Foam::scalar Foam::LangmuirHinshelwoodReactionRate::operator()
 (
     const scalar p,

src/thermophysicalModels/specie/reaction/reactionRate/MichaelisMenten/MichaelisMentenReactionRate.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2018-2020 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2018-2021 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -100,6 +100,12 @@ public:
             return "MichaelisMenten";
         }
 
+        //- Pre-evaluation hook
+        inline void preEvaluate() const;
+
+        //- Post-evaluation hook
+        inline void postEvaluate() const;
+
         inline scalar operator()
         (
             const scalar p,

src/thermophysicalModels/specie/reaction/reactionRate/MichaelisMenten/MichaelisMentenReactionRateI.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2018-2019 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2018-2021 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -42,6 +42,14 @@ inline Foam::MichaelisMentenReactionRate::MichaelisMentenReactionRate
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
+inline void Foam::MichaelisMentenReactionRate::preEvaluate() const
+{}
+
+
+inline void Foam::MichaelisMentenReactionRate::postEvaluate() const
+{}
+
+
 inline Foam::scalar Foam::MichaelisMentenReactionRate::operator()
 (
     const scalar p,

src/thermophysicalModels/specie/reaction/reactionRate/powerSeries/powerSeriesReactionRate.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2011-2019 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2021 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -96,6 +96,12 @@ public:
             return "powerSeries";
         }
 
+        //- Pre-evaluation hook
+        inline void preEvaluate() const;
+
+        //- Post-evaluation hook
+        inline void postEvaluate() const;
+
         inline scalar operator()
         (
             const scalar p,

src/thermophysicalModels/specie/reaction/reactionRate/powerSeries/powerSeriesReactionRateI.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2011-2019 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2021 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -55,6 +55,14 @@ inline Foam::powerSeriesReactionRate::powerSeriesReactionRate
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
+inline void Foam::powerSeriesReactionRate::preEvaluate() const
+{}
+
+
+inline void Foam::powerSeriesReactionRate::postEvaluate() const
+{}
+
+
 inline Foam::scalar Foam::powerSeriesReactionRate::operator()
 (
     const scalar p,

src/thermophysicalModels/specie/reaction/reactionRate/thirdBodyArrheniusReactionRate/thirdBodyArrheniusReactionRate.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-2021 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -93,6 +93,12 @@ public:
             return "thirdBodyArrhenius";
         }
 
+        //- Pre-evaluation hook
+        inline void preEvaluate() const;
+
+        //- Post-evaluation hook
+        inline void postEvaluate() const;
+
         inline scalar operator()
         (
             const scalar p,

src/thermophysicalModels/specie/reaction/reactionRate/thirdBodyArrheniusReactionRate/thirdBodyArrheniusReactionRateI.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2011-2019 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2021 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -72,6 +72,18 @@ inline Foam::thirdBodyArrheniusReactionRate::thirdBodyArrheniusReactionRate
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
+inline void Foam::thirdBodyArrheniusReactionRate::preEvaluate() const
+{
+    ArrheniusReactionRate::preEvaluate();
+}
+
+
+inline void Foam::thirdBodyArrheniusReactionRate::postEvaluate() const
+{
+    ArrheniusReactionRate::postEvaluate();
+}
+
+
 inline Foam::scalar Foam::thirdBodyArrheniusReactionRate::operator()
 (
     const scalar p,

tutorials/multiphase/interFoam/laminar/containerDischarge2D/0/U

@@ -0,0 +1,36 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     | Website:  https://openfoam.org
+    \\  /    A nd           | Version:  dev
+     \\/     M anipulation  |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volVectorField;
+    object      U;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 1 -1 0 0 0 0];
+
+internalField   uniform (0 0 0);
+
+boundaryField
+{
+    "(left|right|atmosphere)"
+    {
+        type            pressureInletOutletVelocity;
+        value           uniform (0 0 0);
+    }
+
+    wall
+    {
+        type            noSlip;
+    }
+}
+
+
+// ************************************************************************* //

tutorials/multiphase/interFoam/laminar/containerDischarge2D/0/alpha.liquid.orig

@@ -0,0 +1,36 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     | Website:  https://openfoam.org
+    \\  /    A nd           | Version:  dev
+     \\/     M anipulation  |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    object      alpha.liquid;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 0 0 0 0 0 0];
+
+internalField   uniform 0;
+
+boundaryField
+{
+    "(left|right|atmosphere)"
+    {
+        type            inletOutlet;
+        inletValue      uniform 0;
+        value           uniform 0;
+    }
+
+    wall
+    {
+        type            zeroGradient;
+    }
+}
+
+// ************************************************************************* //

tutorials/multiphase/interFoam/laminar/containerDischarge2D/0/k

@@ -0,0 +1,37 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     | Website:  https://openfoam.org
+    \\  /    A nd           | Version:  dev
+     \\/     M anipulation  |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    object      k;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 2 -2 0 0 0 0];
+
+internalField   uniform 0.0001;
+
+boundaryField
+{
+    "(left|right|atmosphere)"
+    {
+        type            inletOutlet;
+        inletValue      $internalField;
+        value           $internalField;
+    }
+
+    wall
+    {
+        type            kqRWallFunction;
+        value           $internalField;
+    }
+}
+
+// ************************************************************************* //

tutorials/multiphase/interFoam/laminar/containerDischarge2D/0/nut

@@ -0,0 +1,38 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     | Website:  https://openfoam.org
+    \\  /    A nd           | Version:  dev
+     \\/     M anipulation  |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    location    "0";
+    object      nut;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 2 -1 0 0 0 0];
+
+internalField   uniform 0;
+
+boundaryField
+{
+    "(left|right|atmosphere)"
+    {
+        type            calculated;
+        value           uniform 0;
+    }
+
+    wall
+    {
+        type            nutkWallFunction;
+        value           uniform 0;
+    }
+}
+
+
+// ************************************************************************* //

tutorials/multiphase/interFoam/laminar/containerDischarge2D/0/omega

@@ -0,0 +1,37 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     | Website:  https://openfoam.org
+    \\  /    A nd           | Version:  dev
+     \\/     M anipulation  |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    object      omega;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 0 -1 0 0 0 0];
+
+internalField   uniform 0.003;
+
+boundaryField
+{
+    "(left|right|atmosphere)"
+    {
+        type            inletOutlet;
+        inletValue      $internalField;
+        value           $internalField;
+    }
+
+    wall
+    {
+        type            omegaWallFunction;
+        value           $internalField;
+    }
+}
+
+// ************************************************************************* //

tutorials/multiphase/interFoam/laminar/containerDischarge2D/0/p_rgh

@@ -0,0 +1,36 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     | Website:  https://openfoam.org
+    \\  /    A nd           | Version:  dev
+     \\/     M anipulation  |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    object      p_rgh;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [1 -1 -2 0 0 0 0];
+
+internalField   uniform 0;
+
+boundaryField
+{
+    "(left|right|atmosphere)"
+    {
+        type            prghEntrainmentPressure;
+        p0              uniform 0;
+    }
+
+    wall
+    {
+        type            fixedFluxPressure;
+        value           uniform 0;
+    }
+}
+
+// ************************************************************************* //

tutorials/multiphase/interFoam/laminar/containerDischarge2D/Allclean

@@ -0,0 +1,9 @@
+#!/bin/sh
+cd ${0%/*} || exit 1    # Run from this directory
+
+# Source tutorial clean functions
+. $WM_PROJECT_DIR/bin/tools/CleanFunctions
+
+rm 0/alpha.liquid.gz 0/alpha.liquid 2>/dev/null
+
+cleanCase

tutorials/multiphase/interFoam/laminar/containerDischarge2D/Allrun

@@ -0,0 +1,14 @@
+#!/bin/sh
+cd ${0%/*} || exit 1    # Run from this directory
+
+# Source tutorial run functions
+. $WM_PROJECT_DIR/bin/tools/RunFunctions
+
+application=$(getApplication)
+
+runApplication blockMesh
+
+runApplication setFields
+runApplication $application
+
+#------------------------------------------------------------------------------

tutorials/multiphase/interFoam/laminar/containerDischarge2D/constant/g

@@ -0,0 +1,22 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     | Website:  https://openfoam.org
+    \\  /    A nd           | Version:  dev
+     \\/     M anipulation  |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       uniformDimensionedVectorField;
+    location    "constant";
+    object      g;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 1 -2 0 0 0 0];
+value           (0 0 -9.81);
+
+
+// ************************************************************************* //

tutorials/multiphase/interFoam/laminar/containerDischarge2D/constant/momentumTransport

@@ -0,0 +1,33 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     | Website:  https://openfoam.org
+    \\  /    A nd           | Version:  dev
+     \\/     M anipulation  |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "constant";
+    object      momentumTransport;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+simulationType  laminar;
+
+/*
+simulationType  RAS;
+
+RAS
+{
+    model           kOmegaSST;
+
+    turbulence      on;
+
+    printCoeffs     on;
+}
+*/
+
+// ************************************************************************* //

tutorials/multiphase/interFoam/laminar/containerDischarge2D/constant/transportProperties

@@ -0,0 +1,37 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     | Website:  https://openfoam.org
+    \\  /    A nd           | Version:  dev
+     \\/     M anipulation  |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "constant";
+    object      transportProperties;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+phases (liquid gas);
+
+liquid
+{
+    transportModel  Newtonian;
+    nu              1e-4;
+    rho             950;
+}
+
+gas
+{
+    transportModel  Newtonian;
+    nu              1.48e-05;
+    rho             1;
+}
+
+sigma           0.07;
+
+
+// ************************************************************************* //

tutorials/multiphase/interFoam/laminar/containerDischarge2D/system/blockMeshDict

@@ -0,0 +1,192 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     | Website:  https://openfoam.org
+    \\  /    A nd           | Version:  dev
+     \\/     M anipulation  |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    object      blockMeshDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+convertToMeters 1;
+
+tankRad 0.05;
+holeRad 0.01;
+height  0.1;
+holeCellSize 0.0005;
+tankCellSize 0.001;
+halfDepth    0.01;
+
+nTankCells 20;
+nHoleCells 20;
+nVerCells  60;
+
+horizGrading    3;
+upperGrading    1;
+lowerGrading    1;
+
+tankNeg       #neg $tankRad;
+holeNeg       #neg $holeRad;
+heightNeg     #neg $height;
+halfDepthNeg  #neg $halfDepth;
+
+vertices
+(
+    ($tankNeg $halfDepthNeg $heightNeg)
+    ($holeNeg $halfDepthNeg $heightNeg)
+    ($holeRad $halfDepthNeg $heightNeg)
+    ($tankRad $halfDepthNeg $heightNeg)
+    ($tankNeg $halfDepthNeg 0)
+    ($holeNeg $halfDepthNeg 0)
+    ($holeRad $halfDepthNeg 0)
+    ($tankRad $halfDepthNeg 0)
+    ($tankNeg $halfDepthNeg 0)
+    ($tankRad $halfDepthNeg 0)
+    ($tankNeg $halfDepthNeg $height)
+    ($holeNeg $halfDepthNeg $height)
+    ($holeRad $halfDepthNeg $height)
+    ($tankRad $halfDepthNeg $height)
+
+    ($tankNeg $halfDepth $heightNeg)
+    ($holeNeg $halfDepth $heightNeg)
+    ($holeRad $halfDepth $heightNeg)
+    ($tankRad $halfDepth $heightNeg)
+    ($tankNeg $halfDepth 0)
+    ($holeNeg $halfDepth 0)
+    ($holeRad $halfDepth 0)
+    ($tankRad $halfDepth 0)
+    ($tankNeg $halfDepth 0)
+    ($tankRad $halfDepth 0)
+    ($tankNeg $halfDepth $height)
+    ($holeNeg $halfDepth $height)
+    ($holeRad $halfDepth $height)
+    ($tankRad $halfDepth $height)
+);
+
+blocks
+(
+    hex (15 14 0 1 19 18 4 5)
+    ($nTankCells 1 $nVerCells)
+    simpleGrading ($horizGrading 1 $lowerGrading)
+
+    hex (1 2 16 15 5 6 20 19)
+    ($nHoleCells 1 $nVerCells)
+    simpleGrading (1 1 $lowerGrading)
+
+    hex (2 3 17 16 6 7 21 20)
+    ($nTankCells 1 $nVerCells)
+    simpleGrading ($horizGrading 1 $lowerGrading)
+
+    hex (19 22 8 5 25 24 10 11)
+    ($nTankCells 1 $nVerCells)
+    simpleGrading ($horizGrading 1 $upperGrading)
+
+    hex (5 6 20 19 11 12 26 25)
+    ($nHoleCells 1 $nVerCells)
+    simpleGrading (1 1 $upperGrading)
+
+    hex (6 9 23 20 12 13 27 26)
+    ($nTankCells 1 $nVerCells)
+    simpleGrading ($horizGrading 1 $upperGrading)
+);
+
+edges
+(
+);
+
+boundary
+(
+    left
+    {
+        type patch;
+        faces
+        (
+            (0 4 18 14)
+        );
+    }
+
+    right
+    {
+        type patch;
+        faces
+        (
+            (3 7 21 17)
+        );
+    }
+
+    containerWalls
+    {
+        type wall;
+        faces
+        (
+            (8 10 24 22)
+            (9 13 27 23)
+        );
+    }
+
+    containerBase
+    {
+        type wall;
+        faces
+        (
+            (4 5 19 18)
+            (6 7 21 20)
+            (8 5 19 22)
+            (6 9 23 20)
+        );
+    }
+
+    ground
+    {
+        type wall;
+        faces
+        (
+            (0 1 15 14)
+            (1 2 16 15)
+            (2 3 17 16)
+        );
+    }
+
+    atmosphere
+    {
+        type patch;
+        faces
+        (
+            (10 11 25 24)
+            (11 12 26 25)
+            (12 13 27 26)
+        );
+    }
+
+    frontAndBack
+    {
+        type empty;
+        faces
+        (
+            (0 1 5 4)
+            (1 2 6 5)
+            (2 3 7 6)
+            (8 5 11 10)
+            (5 6 12 11)
+            (6 9 13 12)
+            (14 15 19 18)
+            (15 16 20 19)
+            (16 17 21 20)
+            (22 19 25 24)
+            (19 20 26 25)
+            (20 23 27 26)
+        );
+    }
+);
+
+mergePatchPairs
+(
+);
+
+// ************************************************************************* //

tutorials/multiphase/interFoam/laminar/containerDischarge2D/system/controlDict

@@ -0,0 +1,56 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     | Website:  https://openfoam.org
+    \\  /    A nd           | Version:  dev
+     \\/     M anipulation  |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "system";
+    object      controlDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+application     interFoam;
+
+startFrom       startTime;
+
+startTime       0;
+
+stopAt          endTime;
+
+endTime         1.5;
+
+deltaT          0.003;
+
+writeControl    adjustableRunTime;
+
+writeInterval   0.02;
+
+purgeWrite      0;
+
+writeFormat     ascii;
+
+writePrecision  6;
+
+writeCompression off;
+
+timeFormat      general;
+
+timePrecision   6;
+
+runTimeModifiable yes;
+
+adjustTimeStep  on;
+
+maxCo           5;
+maxAlphaCo      5;
+
+maxDeltaT       1;
+
+
+// ************************************************************************* //

tutorials/multiphase/interFoam/laminar/containerDischarge2D/system/fvSchemes

@@ -0,0 +1,59 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     | Website:  https://openfoam.org
+    \\  /    A nd           | Version:  dev
+     \\/     M anipulation  |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "system";
+    object      fvSchemes;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+ddtSchemes
+{
+    default         Euler;
+}
+
+gradSchemes
+{
+    default         Gauss linear;
+}
+
+divSchemes
+{
+    default             none;
+
+    div(rhoPhi,U)       Gauss linearUpwind grad(U);
+    div(phi,alpha)      Gauss interfaceCompression vanLeer 1;
+    "div\(phi,(k|omega)\)"              Gauss upwind;
+    div(((rho*nuEff)*dev2(T(grad(U))))) Gauss linear;
+}
+
+laplacianSchemes
+{
+    default         Gauss linear corrected;
+}
+
+interpolationSchemes
+{
+    default         linear;
+}
+
+snGradSchemes
+{
+    default         corrected;
+}
+
+wallDist
+{
+    method meshWave;
+}
+
+
+// ************************************************************************* //

tutorials/multiphase/interFoam/laminar/containerDischarge2D/system/fvSolution

@@ -0,0 +1,91 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     | Website:  https://openfoam.org
+    \\  /    A nd           | Version:  dev
+     \\/     M anipulation  |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "system";
+    object      fvSolution;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+solvers
+{
+    "alpha.*"
+    {
+        nAlphaCorr      1;
+        nAlphaSubCycles 1;
+
+        MULESCorr       yes;
+        nLimiterIter    3;
+
+        solver          smoothSolver;
+        smoother        symGaussSeidel;
+        tolerance       1e-8;
+        relTol          0;
+    }
+
+    "pcorr.*"
+    {
+        solver          PCG;
+        preconditioner
+        {
+            preconditioner  GAMG;
+            tolerance       1e-5;
+            relTol          0;
+            smoother        GaussSeidel;
+        }
+        tolerance       1e-5;
+        relTol          0;
+        maxIter         50;
+    }
+
+    p_rgh
+    {
+        solver           GAMG;
+        tolerance        5e-9;
+        relTol           0.01;
+        smoother         GaussSeidel;
+        maxIter          50;
+    };
+
+    p_rghFinal
+    {
+        $p_rgh;
+        tolerance       5e-9;
+        relTol          0;
+    }
+
+    "(U|k|omega).*"
+    {
+        solver          smoothSolver;
+        smoother        symGaussSeidel;
+        nSweeps         1;
+        tolerance       1e-6;
+        relTol          0.1;
+    };
+}
+
+PIMPLE
+{
+    momentumPredictor no;
+    nCorrectors     2;
+    nNonOrthogonalCorrectors 0;
+}
+
+relaxationFactors
+{
+    equations
+    {
+        ".*" 1;
+    }
+}
+
+
+// ************************************************************************* //

tutorials/multiphase/interFoam/laminar/containerDischarge2D/system/setFieldsDict

@@ -0,0 +1,36 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     | Website:  https://openfoam.org
+    \\  /    A nd           | Version:  dev
+     \\/     M anipulation  |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "system";
+    object      setFieldsDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+defaultFieldValues
+(
+    volScalarFieldValue alpha.liquid 0
+);
+
+regions
+(
+    boxToCell
+    {
+        box (-100 -100 0.0) (100 100 0.05);
+        fieldValues
+        (
+            volScalarFieldValue alpha.liquid 1
+        );
+    }
+);
+
+
+// ************************************************************************* //

tutorials/multiphase/multiphaseEulerFoam/laminar/titaniaSynthesis/0/T.particles

@@ -36,10 +36,9 @@ boundaryField
     wall
     {
         type            fixedProfile;
-        profile         table;
+        profile
-
-        profileCoeffs
         {
+            type                tableFile;
             format              csv;        // Input format
             nHeaderLine         0;          // Number of header lines
             refColumn           0;          // Reference column index

tutorials/multiphase/multiphaseEulerFoam/laminar/titaniaSynthesis/0/T.vapor

@@ -36,10 +36,9 @@ boundaryField
     wall
     {
         type            fixedProfile;
-        profile         table;
+        profile
-
-        profileCoeffs
         {
+            type                tableFile;
             format              csv;        // Input format
             nHeaderLine         0;          // Number of header lines
             refColumn           0;          // Reference column index

tutorials/multiphase/multiphaseEulerFoam/laminar/titaniaSynthesisSurface/0/T.particles

@@ -36,10 +36,9 @@ boundaryField
     wall
     {
         type            fixedProfile;
-        profile         table;
+        profile
-
-        profileCoeffs
         {
+            type                tableFile;
             format              csv;        // Input format
             nHeaderLine         0;          // Number of header lines
             refColumn           0;          // Reference column index

tutorials/multiphase/multiphaseEulerFoam/laminar/titaniaSynthesisSurface/0/T.vapor

@@ -36,10 +36,9 @@ boundaryField
     wall
     {
         type            fixedProfile;
-        profile         table;
+        profile
-
-        profileCoeffs
         {
+            type                tableFile;
             format              csv;        // Input format
             nHeaderLine         0;          // Number of header lines
             refColumn           0;          // Reference column index

tutorials/multiphase/multiphaseEulerFoam/laminar/titaniaSynthesisSurface/constant/phaseProperties

@@ -78,8 +78,6 @@ particles
             f28{dSph 7.873E-06; value 0.0; kappa 1.5E+10; Df 1.8; alphaC 1;}
             f29{dSph 1.136E-05; value 0.0; kappa 1.5E+10; Df 1.8; alphaC 1;}
         );
-
-        storeA              on;
     }
 
     Sct             1.0;

tutorials/multiphase/multiphaseEulerFoam/laminar/titaniaSynthesisSurface/constant/reactions.vapor

@@ -40,7 +40,7 @@ reactions
 
     oxidationAtSurface
     {
-        type        irreversibleSurfaceArrhenius;
+        type        irreversiblePhaseSurfaceArrhenius;
 
         reaction    "O2^0 + TiCl4 = TiO2_s + 2Cl2";
 
@@ -48,7 +48,7 @@ reactions
         beta        0.0;
         Ta          8993;
 
-        a           a.particles;
+        phase       particles;
     }
 }