diff --git a/CONTRIBUTORS.md b/CONTRIBUTORS.md
index bcf1033bcb338b674a52c118c5df79345ea84a10..30d191014adaddae0fba39e92472da54d57737a4 100644
--- a/CONTRIBUTORS.md
+++ b/CONTRIBUTORS.md
@@ -50,3 +50,7 @@ It is likely incomplete...
 - Norbert Weber
 - Henry Weller
 - Niklas Wikstrom
+- Thorsten Zirwes
+
+
+<!----------------------------------------------------------------------------->
diff --git a/src/thermophysicalModels/chemistryModel/Make/files b/src/thermophysicalModels/chemistryModel/Make/files
index 630137c0e6e9ef5a69c19b8974423ec22a2d8ec3..a77950e3bd4cd71fac6ded40949030088b951cbe 100644
--- a/src/thermophysicalModels/chemistryModel/Make/files
+++ b/src/thermophysicalModels/chemistryModel/Make/files
@@ -7,5 +7,6 @@ chemistryModel/TDACChemistryModel/tabulation/makeChemistryTabulationMethods.C
 chemistrySolver/chemistrySolver/makeChemistrySolvers.C
 
 functionObjects/specieReactionRates/specieReactionRates.C
+functionObjects/BilgerMixtureFraction/BilgerMixtureFraction.C
 
 LIB = $(FOAM_LIBBIN)/libchemistryModel
diff --git a/src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/TDACChemistryModel.C b/src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/TDACChemistryModel.C
index 187e30093243689fc7e297716051abc6aeb908ad..c4d502ebe0713823f48f663b16ebf65ea563980d 100644
--- a/src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/TDACChemistryModel.C
+++ b/src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/TDACChemistryModel.C
@@ -75,13 +75,15 @@ Foam::TDACChemistryModel<ReactionThermo, ThermoType>::TDACChemistryModel
     // Store the species composition according to the species index
     speciesTable speciesTab = composition.species();
 
-    const HashTable<List<specieElement>>& specComp =
+    autoPtr<HashTable<List<specieElement>>> specCompPtr
+    (
         dynamicCast<const reactingMixture<ThermoType>&>(this->thermo())
-       .specieComposition();
+       .specieComposition()
+    );
 
     forAll(specieComp_, i)
     {
-        specieComp_[i] = specComp[this->Y()[i].member()];
+        specieComp_[i] = (specCompPtr.ref())[this->Y()[i].member()];
     }
 
     mechRed_ = chemistryReductionMethod<ReactionThermo, ThermoType>::New
diff --git a/src/thermophysicalModels/chemistryModel/functionObjects/BilgerMixtureFraction/BilgerMixtureFraction.C b/src/thermophysicalModels/chemistryModel/functionObjects/BilgerMixtureFraction/BilgerMixtureFraction.C
new file mode 100644
index 0000000000000000000000000000000000000000..65a93a982f222e94991e440e8e2fb3a033c573bb
--- /dev/null
+++ b/src/thermophysicalModels/chemistryModel/functionObjects/BilgerMixtureFraction/BilgerMixtureFraction.C
@@ -0,0 +1,409 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2020 Thorsten Zirwes
+    Copyright (C) 2020 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+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 "BilgerMixtureFraction.H"
+#include "basicThermo.H"
+#include "reactingMixture.H"
+#include "thermoPhysicsTypes.H"
+#include "scalarRange.H"
+#include "basicChemistryModel.H"
+#include "psiReactionThermo.H"
+#include "rhoReactionThermo.H"
+#include "BasicChemistryModel.H"
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace functionObjects
+{
+    defineTypeNameAndDebug(BilgerMixtureFraction, 0);
+    addToRunTimeSelectionTable
+    (
+        functionObject,
+        BilgerMixtureFraction,
+        dictionary
+    );
+}
+}
+
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+void Foam::functionObjects::BilgerMixtureFraction::calcBilgerMixtureFraction()
+{
+    if (!mesh_.foundObject<volScalarField>(resultName_, false))
+    {
+        auto tCo = tmp<volScalarField>::New
+        (
+            IOobject
+            (
+                resultName_,
+                mesh_.time().timeName(),
+                mesh_,
+                IOobject::NO_READ,
+                IOobject::NO_WRITE
+            ),
+            mesh_,
+            dimensionedScalar(dimless, Zero)
+        );
+        mesh_.objectRegistry::store(tCo.ptr());
+    }
+
+    auto& f_Bilger = mesh_.lookupObjectRef<volScalarField>(resultName_);
+
+    auto& Y = thermo_.Y();
+
+    f_Bilger = -o2RequiredOx_;
+    forAll(Y, i)
+    {
+        f_Bilger +=
+            Y[i]
+           *(nAtomsC_[i] + nAtomsS_[i] + 0.25*nAtomsH_[i] - 0.5*nAtomsO_[i])
+           /thermo_.W(i);
+    }
+    f_Bilger /= o2RequiredFuelOx_;
+    f_Bilger.clip
+    (
+        dimensionedScalar(dimless, 0),
+        dimensionedScalar(dimless, 1)
+    );
+}
+
+
+bool Foam::functionObjects::BilgerMixtureFraction::readComposition
+(
+    const dictionary& subDict,
+    scalarField& comp
+)
+{
+    auto& Y = thermo_.Y();
+    const speciesTable& speciesTab = thermo_.species();
+
+    // Read mass fractions of all species for the oxidiser or fuel
+    forAll(Y, i)
+    {
+        comp[i] =
+            subDict.getCheckOrDefault<scalar>
+            (
+                speciesTab[i],
+                0,
+                scalarRange::ge0()
+            );
+    }
+
+    if (sum(comp) < SMALL)
+    {
+        FatalIOErrorInFunction(subDict)
+            << "No composition is given" << nl
+            << "Valid species are:" << nl
+            << speciesTab
+            << exit(FatalIOError);
+
+        return false;
+    }
+
+    const word fractionBasisType
+    (
+        subDict.getOrDefault<word>("fractionBasis", "mass")
+    );
+
+    if (fractionBasisType == "mass")
+    {
+        // Normalize fractionBasis to the unity
+        comp /= sum(comp);
+    }
+    else if (fractionBasisType == "mole")
+    {
+        // In case the fractionBasis is given in mole fractions,
+        // convert from mole fractions to normalized mass fractions
+        scalar W(0);
+        forAll(comp, i)
+        {
+            comp[i] *= thermo_.W(i);
+            W += comp[i];
+        }
+        comp /= W;
+    }
+    else
+    {
+        FatalIOErrorInFunction(subDict)
+            << "The given fractionBasis type is invalid" << nl
+            << "Valid fractionBasis types are" << nl
+            << "  \"mass\" (default)" << nl
+            << "  \"mole\""
+            << exit(FatalIOError);
+
+        return false;
+    }
+
+    return true;
+}
+
+
+Foam::scalar Foam::functionObjects::BilgerMixtureFraction::o2Required
+(
+    const scalarField& comp
+) const
+{
+    scalar o2req(0);
+    forAll(thermo_.Y(), i)
+    {
+        o2req +=
+            comp[i]/thermo_.W(i)*(nAtomsC_[i] + nAtomsS_[i] + 0.25*nAtomsH_[i]);
+    }
+
+    return o2req;
+}
+
+
+Foam::scalar Foam::functionObjects::BilgerMixtureFraction::o2Present
+(
+    const scalarField& comp
+) const
+{
+    scalar o2pres(0);
+    forAll(thermo_.Y(), i)
+    {
+        o2pres += comp[i]/thermo_.W(i)*nAtomsO_[i];
+    }
+
+    return 0.5*o2pres;
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::functionObjects::BilgerMixtureFraction::BilgerMixtureFraction
+(
+    const word& name,
+    const Time& runTime,
+    const dictionary& dict
+)
+:
+    fvMeshFunctionObject(name, runTime, dict),
+    phaseName_(dict.getOrDefault<word>("phase", word::null)),
+    resultName_
+    (
+        dict.getOrDefault<word>
+        (
+            "result",
+            IOobject::groupName("f_Bilger", phaseName_)
+        )
+    ),
+    thermo_
+    (
+        mesh_.lookupObject<basicSpecieMixture>
+        (
+            IOobject::groupName(basicThermo::dictName, phaseName_)
+        )
+    ),
+    nSpecies_(thermo_.Y().size()),
+    o2RequiredOx_(0),
+    o2RequiredFuelOx_(0),
+    nAtomsC_(nSpecies_, 0),
+    nAtomsS_(nSpecies_, 0),
+    nAtomsH_(nSpecies_, 0),
+    nAtomsO_(nSpecies_, 0),
+    Yoxidiser_(nSpecies_, 0),
+    Yfuel_(nSpecies_, 0)
+{
+    read(dict);
+
+    calcBilgerMixtureFraction();
+}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+bool Foam::functionObjects::BilgerMixtureFraction::read(const dictionary& dict)
+{
+    if (!fvMeshFunctionObject::read(dict))
+    {
+        return false;
+    }
+
+    Info<< nl << type() << " " << name() << ":" << nl;
+
+    phaseName_ = dict.getOrDefault<word>("phase", word::null);
+    resultName_ =
+        dict.getOrDefault<word>
+        (
+            "result",
+            IOobject::groupName("f_Bilger", phaseName_)
+        );
+
+    nSpecies_ = thermo_.Y().size();
+
+    if (nSpecies_ == 0)
+    {
+        FatalErrorInFunction
+            << "Number of input species is zero"
+            << exit(FatalError);
+    }
+
+    nAtomsC_.resize(nSpecies_, 0);
+    nAtomsS_.resize(nSpecies_, 0);
+    nAtomsH_.resize(nSpecies_, 0);
+    nAtomsO_.resize(nSpecies_, 0);
+
+    auto& Y = thermo_.Y();
+    const speciesTable& speciesTab = thermo_.species();
+
+    typedef BasicChemistryModel<psiReactionThermo> psiChemistryModelType;
+    typedef BasicChemistryModel<rhoReactionThermo> rhoChemistryModelType;
+
+    const auto* psiChemPtr =
+        mesh_.cfindObject<psiChemistryModelType>("chemistryProperties");
+
+    const auto* rhoChemPtr =
+        mesh_.cfindObject<rhoChemistryModelType>("chemistryProperties");
+
+    autoPtr<HashTable<List<specieElement>>> speciesCompPtr;
+
+    if (psiChemPtr)
+    {
+        speciesCompPtr.reset((*psiChemPtr).thermo().specieComposition());
+    }
+    else if (rhoChemPtr)
+    {
+        speciesCompPtr.reset((*rhoChemPtr).thermo().specieComposition());
+    }
+    else
+    {
+        FatalErrorInFunction
+            << "BasicChemistryModel not found"
+            << exit(FatalError);
+    }
+
+    forAll(Y, i)
+    {
+        const List<specieElement>& curSpecieComposition =
+            (speciesCompPtr.ref())[speciesTab[i]];
+
+        forAll(curSpecieComposition, j)
+        {
+            const word& e = curSpecieComposition[j].name();
+            const label nAtoms  = curSpecieComposition[j].nAtoms();
+
+            if (e == "C")
+            {
+                nAtomsC_[i] = nAtoms;
+            }
+            else if (e == "S")
+            {
+                nAtomsS_[i] = nAtoms;
+            }
+            else if (e == "H")
+            {
+                nAtomsH_[i] = nAtoms;
+            }
+            else if (e == "O")
+            {
+                nAtomsO_[i] = nAtoms;
+            }
+        }
+    }
+
+    if (sum(nAtomsO_) == 0)
+    {
+        FatalErrorInFunction
+            << "No specie contains oxygen"
+            << exit(FatalError);
+    }
+
+    Yoxidiser_.resize(nSpecies_, 0);
+    Yfuel_.resize(nSpecies_, 0);
+
+    if
+    (
+        !readComposition(dict.subDict("oxidiser"), Yoxidiser_)
+     || !readComposition(dict.subDict("fuel"), Yfuel_)
+    )
+    {
+        return false;
+    }
+
+    o2RequiredOx_ = o2Required(Yoxidiser_) - o2Present(Yoxidiser_);
+
+    if (o2RequiredOx_ > 0)
+    {
+        FatalErrorInFunction
+            << "Oxidiser composition contains not enough oxygen" << endl
+            << "Mixed up fuel and oxidiser compositions?"
+            << exit(FatalError);
+    }
+
+    const scalar o2RequiredFuel = o2Required(Yfuel_) - o2Present(Yfuel_);
+
+    if (o2RequiredFuel < 0)
+    {
+        FatalErrorInFunction
+            << "Fuel composition contains too much oxygen" << endl
+            << "Mixed up fuel and oxidiser compositions?"
+            << exit(FatalError);
+    }
+
+    o2RequiredFuelOx_ = o2RequiredFuel - o2RequiredOx_;
+
+    if (mag(o2RequiredFuelOx_) < SMALL)
+    {
+        FatalErrorInFunction
+            << "Fuel and oxidiser have the same composition"
+            << exit(FatalError);
+    }
+
+    return true;
+}
+
+
+bool Foam::functionObjects::BilgerMixtureFraction::execute()
+{
+    calcBilgerMixtureFraction();
+
+    return true;
+}
+
+
+bool Foam::functionObjects::BilgerMixtureFraction::clear()
+{
+    return clearObject(resultName_);
+}
+
+
+bool Foam::functionObjects::BilgerMixtureFraction::write()
+{
+    Log << type() << " " << name() << " write:" << nl
+        << "    writing field " << resultName_ << endl;
+
+    return writeObject(resultName_);
+}
+
+
+// ************************************************************************* //
diff --git a/src/thermophysicalModels/chemistryModel/functionObjects/BilgerMixtureFraction/BilgerMixtureFraction.H b/src/thermophysicalModels/chemistryModel/functionObjects/BilgerMixtureFraction/BilgerMixtureFraction.H
new file mode 100644
index 0000000000000000000000000000000000000000..66259b93fe279a5fb90e423abfc09e91b8a58ff5
--- /dev/null
+++ b/src/thermophysicalModels/chemistryModel/functionObjects/BilgerMixtureFraction/BilgerMixtureFraction.H
@@ -0,0 +1,301 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2020 Thorsten Zirwes
+    Copyright (C) 2020 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+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::functionObjects::BilgerMixtureFraction
+
+Group
+    grpThermophysicalFunctionObjects
+
+Description
+    Calculates the Bilger mixture-fraction field (i.e.
+    \f$f_{\mathrm{Bilger}}\f$) based on the elemental composition
+    of the mixture. Elements \c C, \c H, \c S and \c O are considered.
+
+    \f$f_{\mathrm{Bilger}}\f$ is the mass mixing ratio of fuel and oxidiser
+    in [kg fuel / kg mixture], and is invariant to the reaction progress of
+    the mixture (e.g. the same for unburnt and burnt mixtures).
+
+    \f$f_{\mathrm{Bilger}}\f$ equals to the unity
+    for pure fuel and to zero for pure oxidiser.
+
+    The Bilger mixture-fraction field is computed based on the following:
+
+    \f[
+        f_{\mathrm{Bilger}} =
+            \frac{\beta - \beta_{\mathrm{ox}}}{\beta_{\mathrm{fuel}}
+          - \beta_{\mathrm{ox}}}
+    \f]
+
+    with
+
+    \f[
+        \beta =
+            2\frac{Y_C}{W_C}
+          + 2\frac{Y_S}{W_S}
+          + \frac{1}{2}\frac{Y_H}{W_H}
+          - \frac{Y_O}{W_O}
+    \f]
+
+    where
+    \vartable
+      \beta               | Coupling function                        [kmol/kg]
+      Y_e                 | Elemental mass fraction of element, e
+      W_e                 | Atomic weight of element, e
+      {.}_{\mathrm{ox}}   | Subscript to denote oxidiser composition
+      {.}_{\mathrm{fuel}} | Subscript to denote fuel composition
+    \endvartable
+
+    Operands:
+    \table
+      Operand        | Type           | Location
+      input          | -              | -
+      output file    | -              | -
+      output field   | volScalarField | $FOAM_CASE/\<time\>/\<outField\>
+    \endtable
+
+    References:
+    \verbatim
+        Original method:
+            Bilger, R. W. (1979).
+            Turbulent jet diffusion flames.
+            Energy and Combustion Science, p. 109-131.
+            DOI:10.1016/B978-0-08-024780-9.50011-3
+
+        Implementation:
+            Zirwes, T., Zhang, F., Habisreuther, P., Hansinger, M.,
+            Bockhorn, H., Pfitzner, M., & Trimis, D. (2019).
+            Quasi-DNS dataset of a piloted flame
+            with inhomogeneous inlet conditions.
+            Flow, Turbulence and Combustion, vol 104, p. 997-1027.
+            DOI:10.1007/s10494-019-00081-5
+    \endverbatim
+
+Usage
+    Minimal example by using \c system/controlDict.functions:
+    \verbatim
+    BilgerMixtureFraction1
+    {
+        // Mandatory entries (unmodifiable)
+        type                 BilgerMixtureFraction;
+        libs                 (fieldFunctionObjects);
+
+        // Mandatory entries (runtime modifiable)
+        fuel
+        {
+            // Optional entries (runtime modifiable)
+            fractionBasis    mass;
+
+            // Conditional mandatory entries (runtime modifiable)
+            CH4              1;
+        }
+
+        oxidiser
+        {
+            // Optional entries (runtime modifiable)
+            fractionBasis    mole;
+
+            // Conditional mandatory entries (runtime modifiable)
+            O2               0.23;
+            N2               0.77;
+        }
+
+        // Optional entries (runtime modifiable)
+        phase                <phaseName>;
+        result               <resultName>;
+
+        // Optional (inherited) entries
+        ...
+    }
+    \endverbatim
+
+    where the entries mean:
+    \table
+      Property | Description                         | Type | Reqd | Dflt
+      type     | Type name: BilgerMixtureFraction    | word | yes  | -
+      libs     | Library name: fieldFunctionObjects  | word | yes  | -
+      fuel     | Dictionary for fuel composition     | dict | yes  | -
+      oxidiser | Dictionary for oxidiser composition | dict | yes  | -
+      phase    | Name of phase (e.g. "gas")          | word | no   | ""
+      result   | Name of resulting field             | word | no   | f_Bilger
+      fractionBasis | Species-fraction interpretation method | word | no | mass
+    \endtable
+
+    Options for the \c fractionBasis entry:
+    \verbatim
+      mass   | Interpret species fractions as mass fractions
+      mole   | Interpret species fractions as mole fractions
+    \endverbatim
+
+    The inherited entries are elaborated in:
+     - \link functionObject.H \endlink
+
+    Usage by the \c postProcess utility is not available.
+
+Note
+  - The mole or mass fractions are automatically normalized to the unity.
+
+See also
+  - Foam::functionObject
+  - Foam::functionObjects::fvMeshFunctionObject
+
+SourceFiles
+    BilgerMixtureFraction.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef BilgerMixtureFraction_H
+#define BilgerMixtureFraction_H
+
+#include "fvMeshFunctionObject.H"
+#include "specieElement.H"
+#include "basicSpecieMixture.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace functionObjects
+{
+/*---------------------------------------------------------------------------*\
+                    Class BilgerMixtureFraction Declaration
+\*---------------------------------------------------------------------------*/
+
+class BilgerMixtureFraction
+:
+    public fvMeshFunctionObject
+{
+    // Private Data
+
+        //- Name of the phase (e.g. "gas" for multiphase applications)
+        word phaseName_;
+
+        //- Name of the resulting mixture-fraction field
+        word resultName_;
+
+        //- Reference to thermo object
+        const basicSpecieMixture& thermo_;
+
+        //- Number of species
+        label nSpecies_;
+
+        // Amount of oxygen required to fully oxidise the oxidiser
+        scalar o2RequiredOx_;
+
+        // Amount of oxygen required to oxidise the fuel minus the oxidiser
+        scalar o2RequiredFuelOx_;
+
+        //- Number of carbon atoms for each species
+        labelField nAtomsC_;
+
+        //- Number of sulphur atoms for each species
+        labelField nAtomsS_;
+
+        //- Number of hydrogen atoms for each species
+        labelField nAtomsH_;
+
+        //- Number of oxygen atoms for each species
+        labelField nAtomsO_;
+
+        //- Mass fractions of species in the oxidiser
+        scalarField Yoxidiser_;
+
+        //- Mass fractions of species in the fuel
+        scalarField Yfuel_;
+
+
+    // Private Member Functions
+
+        //- Calculate the Bilger mixture-fraction
+        void calcBilgerMixtureFraction();
+
+        //- Read composition of fuel and oxidiser from subdictionary
+        bool readComposition
+        (
+            const dictionary& subDict,
+            scalarField& comp
+        );
+
+        //- Compute amount of oxygen required to oxidise a mixture
+        scalar o2Present(const scalarField&) const;
+
+        //- Compute amount of oxygen present in a mixture
+        scalar o2Required(const scalarField&) const;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("BilgerMixtureFraction");
+
+
+    // Constructors
+
+        //- Construct from Time and dictionary
+        BilgerMixtureFraction
+        (
+            const word& name,
+            const Time& runTime,
+            const dictionary& dict
+        );
+
+        //- No copy construct
+        BilgerMixtureFraction(const BilgerMixtureFraction&) = delete;
+
+        //- No copy assignment
+        void operator=(const BilgerMixtureFraction&) = delete;
+
+
+    //- Destructor
+    virtual ~BilgerMixtureFraction() = default;
+
+
+    // Member Functions
+
+        //- Read the BilgerMixtureFraction data
+        virtual bool read(const dictionary&);
+
+        //- Calculate the Bilger mixture-fraction field
+        virtual bool execute();
+
+        //- Clear the Bilger mixture-fraction field from registry
+        virtual bool clear();
+
+        //- Write Bilger mixture-fraction field
+        virtual bool write();
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace functionObjects
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/thermophysicalModels/reactionThermo/Make/files b/src/thermophysicalModels/reactionThermo/Make/files
index 8870103dfe05d72dc41520b0d43172d771f9103e..e31e2076ffe1c67de3cebc5276d7b84bbdcfb0c5 100644
--- a/src/thermophysicalModels/reactionThermo/Make/files
+++ b/src/thermophysicalModels/reactionThermo/Make/files
@@ -21,4 +21,5 @@ derivedFvPatchFields/mixedUnburntEnthalpy/mixedUnburntEnthalpyFvPatchScalarField
 
 functionObjects/moleFractions/moleFractionsFunctionObjects.C
 
+
 LIB = $(FOAM_LIBBIN)/libreactionThermophysicalModels
diff --git a/src/thermophysicalModels/reactionThermo/Make/options b/src/thermophysicalModels/reactionThermo/Make/options
index 26ee649fb42ad834cfcc7da4e92e521fc56d57e7..20754771fa29104edd861292d33525bcb38aa69f 100644
--- a/src/thermophysicalModels/reactionThermo/Make/options
+++ b/src/thermophysicalModels/reactionThermo/Make/options
@@ -4,6 +4,7 @@ EXE_INC = \
     -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/specie/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/solidSpecie/lnInclude \
+    -I$(LIB_SRC)/functionObjects/field/lnInclude
 
 LIB_LIBS = \
     -lfiniteVolume \
diff --git a/src/thermophysicalModels/reactionThermo/functionObjects/moleFractions/moleFractions.C b/src/thermophysicalModels/reactionThermo/functionObjects/moleFractions/moleFractions.C
index c32cf62734431a83358fb4bd1517fb3982fb057d..1787595dfd38ad8d39fe1545b5003946ecb88d70 100644
--- a/src/thermophysicalModels/reactionThermo/functionObjects/moleFractions/moleFractions.C
+++ b/src/thermophysicalModels/reactionThermo/functionObjects/moleFractions/moleFractions.C
@@ -5,7 +5,8 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2016-2017 OpenFOAM Foundation
+    Copyright (C) 2016-2020 OpenFOAM Foundation
+    Copyright (C) 2020 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -31,10 +32,13 @@ License
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
 template<class ThermoType>
-void Foam::moleFractions<ThermoType>::calculateMoleFractions()
+void Foam::moleFractions<ThermoType>::calcMoleFractions()
 {
-    const ThermoType& thermo =
-        mesh_.lookupObject<ThermoType>(basicThermo::dictName);
+    const auto& thermo =
+        mesh_.lookupObject<ThermoType>
+        (
+            IOobject::groupName(basicThermo::dictName, phaseName_)
+        );
 
     const PtrList<volScalarField>& Y = thermo.composition().Y();
 
@@ -44,7 +48,6 @@ void Foam::moleFractions<ThermoType>::calculateMoleFractions()
     {
         const dimensionedScalar Wi
         (
-            "W",
             dimMass/dimMoles,
             thermo.composition().W(i)
         );
@@ -64,14 +67,19 @@ Foam::moleFractions<ThermoType>::moleFractions
     const dictionary& dict
 )
 :
-    fvMeshFunctionObject(name, runTime, dict)
+    fvMeshFunctionObject(name, runTime, dict),
+    phaseName_(dict.getOrDefault<word>("phase", word::null))
 {
-    const ThermoType* thermo =
-        mesh_.findObject<ThermoType>(basicThermo::dictName);
+    const word dictName
+    (
+        IOobject::groupName(basicThermo::dictName, phaseName_)
+    );
 
-    if (thermo)
+    if (mesh_.foundObject<ThermoType>(dictName))
     {
-        const PtrList<volScalarField>& Y = thermo->composition().Y();
+        const auto& thermo = mesh_.lookupObject<ThermoType>(dictName);
+
+        const PtrList<volScalarField>& Y = thermo.composition().Y();
 
         X_.setSize(Y.size());
 
@@ -96,10 +104,19 @@ Foam::moleFractions<ThermoType>::moleFractions
             );
         }
 
-        calculateMoleFractions();
+        calcMoleFractions();
     }
     else
     {
+        if (phaseName_ != word::null)
+        {
+            FatalErrorInFunction
+                << "Cannot find thermodynamics model of type "
+                << ThermoType::typeName
+                << " for phase " << phaseName_
+                << exit(FatalError);
+        }
+
         FatalErrorInFunction
             << "Cannot find thermodynamics model of type "
             << ThermoType::typeName
@@ -108,13 +125,6 @@ Foam::moleFractions<ThermoType>::moleFractions
 }
 
 
-// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
-
-template<class ThermoType>
-Foam::moleFractions<ThermoType>::~moleFractions()
-{}
-
-
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
 template<class ThermoType>
@@ -123,21 +133,23 @@ bool Foam::moleFractions<ThermoType>::read
     const dictionary& dict
 )
 {
-    return true;
-}
+    if (functionObjects::fvMeshFunctionObject::read(dict))
+    {
+        phaseName_ = dict.getOrDefault<word>("phase", word::null);
 
+        return true;
+    }
 
-template<class ThermoType>
-bool Foam::moleFractions<ThermoType>::execute()
-{
-    calculateMoleFractions();
-    return true;
+    return false;
 }
 
 
+
 template<class ThermoType>
-bool Foam::moleFractions<ThermoType>::write()
+bool Foam::moleFractions<ThermoType>::execute()
 {
+    calcMoleFractions();
+
     return true;
 }
 
diff --git a/src/thermophysicalModels/reactionThermo/functionObjects/moleFractions/moleFractions.H b/src/thermophysicalModels/reactionThermo/functionObjects/moleFractions/moleFractions.H
index a85c7add199fc1bb52e9573966b9d1e7f6be6f56..3bd9f16d59acc6a725bd6c3db21b7bdd8928dd38 100644
--- a/src/thermophysicalModels/reactionThermo/functionObjects/moleFractions/moleFractions.H
+++ b/src/thermophysicalModels/reactionThermo/functionObjects/moleFractions/moleFractions.H
@@ -5,7 +5,8 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2016 OpenFOAM Foundation
+    Copyright (C) 2016-2020 OpenFOAM Foundation
+    Copyright (C) 2020 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -30,34 +31,67 @@ Group
     grpThermophysicalFunctionObjects
 
 Description
-    This function object calculates mole-fraction fields from the mass-fraction
-    fields of the psi/rhoReactionThermo and caches them for output and further
-    post-processing.
-
-    The names of the mole-fraction fields are obtained from the corresponding
-    mass-fraction fields prepended by "X_"
-
-    Example of function object specification:
+    Calculates mole-fraction fields from the mass-fraction
+    fields of the psi/rhoReactionThermo and caches them
+    for output and further post-processing.
+
+    The names of the mole-fraction fields are obtained from
+    the corresponding mass-fraction fields prepended by "X_".
+
+    Operands:
+    \table
+      Operand        | Type           | Location
+      input          | -              | -
+      output file    | -              | -
+      output field   | volScalarField | $FOAM_CASE/\<time\>/\<outField\>
+    \endtable
+
+Usage
+    Minimal example by using \c system/controlDict.functions:
     \verbatim
-    moleFractions
+    moleFractions1
     {
-        type psiReactionThermoMoleFractions;
-    }
-    \endverbatim
-    or
-    \verbatim
-    moleFractions
-    {
-        type rhoReactionThermoMoleFractions;
+        // Conditional mandatory entries (unmodifiable)
+        // Either of the below depending on
+        // the thermodynamics package used in the solver.
+
+            // Option-1
+            type                 psiReactionThermoMoleFractions;
+
+            // Option-2
+            type                 rhoReactionThermoMoleFractions;
+
+        // Mandatory entries (unmodifiable)
+        libs                 (fieldFunctionObjects);
+
+        // Optional entries (runtime modifiable)
+        phase                <phaseName>;
+
+        // Optional (inherited) entries
+        ...
     }
     \endverbatim
-    depending on the thermodynamics package used in the solver.
+
+    where the entries mean:
+    \table
+      Property | Description                         | Type | Reqd | Dflt
+      type     | Type name: psiReactionThermoMoleFractions or <!--
+                --> rhoReactionThermoMoleFractions   | word | yes  | -
+      libs     | Library name: fieldFunctionObjects  | word | yes  | -
+      phase    | Name of phase (e.g. "gas")          | word | no   | ""
+    \endtable
+
+    The inherited entries are elaborated in:
+     - \link functionObject.H \endlink
+
+    Usage by the \c postProcess utility is not available.
 
 See also
     Foam::functionObjects::fvMeshFunctionObject
 
 SourceFiles
     moleFractions.C
+    moleFractionsFunctionObjects.C
 
 \*---------------------------------------------------------------------------*/
 
@@ -81,22 +115,19 @@ class moleFractions
 :
     public functionObjects::fvMeshFunctionObject
 {
-    // Private data
+    // Private Data
 
         //- Species mole fractions
         PtrList<volScalarField> X_;
 
+        //- Name of phase
+        word phaseName_;
+
 
     // Private Member Functions
 
         //- Calculate the mole fraction fields
-        virtual void calculateMoleFractions();
-
-        //- No copy construct
-        moleFractions(const moleFractions&) = delete;
-
-        //- No copy assignment
-        void operator=(const moleFractions&) = delete;
+        virtual void calcMoleFractions();
 
 
 public:
@@ -115,21 +146,30 @@ public:
             const dictionary& dict
         );
 
+        //- No copy construct
+        moleFractions(const moleFractions&) = delete;
+
+        //- No copy assignment
+        void operator=(const moleFractions&) = delete;
+
 
     //- Destructor
-    virtual ~moleFractions();
+    virtual ~moleFractions() = default;
 
 
     // Member Functions
 
         //- Read the moleFractions data
-        virtual bool read(const dictionary&);
+        virtual bool read(const dictionary& dict);
 
         //- Calculate the mole-fraction fields
         virtual bool execute();
 
-        //- The mole-fraction fields auto-write
-        virtual bool write();
+        //- The mole-fraction fields auto-write - no-op
+        virtual bool write()
+        {
+            return true;
+        }
 };
 
 
diff --git a/src/thermophysicalModels/reactionThermo/mixtures/basicSpecieMixture/basicSpecieMixture.H b/src/thermophysicalModels/reactionThermo/mixtures/basicSpecieMixture/basicSpecieMixture.H
index d7db29cc424e7a034936c660ca9f1672757cb618..522c4be19e91a40a08bb99907ea89b0cc828552d 100644
--- a/src/thermophysicalModels/reactionThermo/mixtures/basicSpecieMixture/basicSpecieMixture.H
+++ b/src/thermophysicalModels/reactionThermo/mixtures/basicSpecieMixture/basicSpecieMixture.H
@@ -6,6 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2014-2017 OpenFOAM Foundation
+    Copyright (C) 2020 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -42,12 +43,15 @@ SourceFiles
 #define basicSpecieMixture_H
 
 #include "basicMultiComponentMixture.H"
+#include "specieElement.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 namespace Foam
 {
 
+typedef HashTable<List<specieElement>> speciesCompositionTable;
+
 /*---------------------------------------------------------------------------*\
                  Class basicSpecieMixture Declaration
 \*---------------------------------------------------------------------------*/
@@ -202,6 +206,16 @@ public:
                 const scalar p,
                 const scalar T
             ) const = 0;
+
+            //- Species composition
+            virtual autoPtr<speciesCompositionTable> specieComposition() const
+            {
+                return
+                    autoPtr<speciesCompositionTable>
+                    (
+                        new speciesCompositionTable()
+                    );
+            }
 };
 
 
diff --git a/src/thermophysicalModels/reactionThermo/mixtures/reactingMixture/reactingMixture.H b/src/thermophysicalModels/reactionThermo/mixtures/reactingMixture/reactingMixture.H
index 08763c8a2a513b7b4125f5ef745fe89239c5b52e..bb33b6cfbf41d7faa07c0a056d13693fbe97b72f 100644
--- a/src/thermophysicalModels/reactionThermo/mixtures/reactingMixture/reactingMixture.H
+++ b/src/thermophysicalModels/reactionThermo/mixtures/reactingMixture/reactingMixture.H
@@ -6,6 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2017 OpenFOAM Foundation
+    Copyright (C) 2020 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -115,9 +116,13 @@ public:
         using PtrList<Reaction<ThermoType>>::operator[];
 
         //- Table of species composition
-        const speciesCompositionTable& specieComposition() const
+        virtual autoPtr<speciesCompositionTable> specieComposition() const
         {
-            return speciesComposition_;
+            return
+                autoPtr<speciesCompositionTable>
+                (
+                    new speciesCompositionTable(speciesComposition_)
+                );
         }
 };
 
diff --git a/src/thermophysicalModels/reactionThermo/psiReactionThermo/psiReactionThermo.H b/src/thermophysicalModels/reactionThermo/psiReactionThermo/psiReactionThermo.H
index fe1035f45a9224ee42dc1e32d8927af427ef306c..080e26315e1755989fce903ad370cc357c8538db 100644
--- a/src/thermophysicalModels/reactionThermo/psiReactionThermo/psiReactionThermo.H
+++ b/src/thermophysicalModels/reactionThermo/psiReactionThermo/psiReactionThermo.H
@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2015 OpenFOAM Foundation
-    Copyright (C) 2017 OpenCFD Ltd.
+    Copyright (C) 2017-2020 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -45,12 +45,15 @@ SourceFiles
 #include "basicSpecieMixture.H"
 #include "autoPtr.H"
 #include "runTimeSelectionTables.H"
+#include "specieElement.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 namespace Foam
 {
 
+typedef HashTable<List<specieElement>> speciesCompositionTable;
+
 /*---------------------------------------------------------------------------*\
                      Class psiReactionThermo Declaration
 \*---------------------------------------------------------------------------*/
@@ -137,6 +140,16 @@ public:
 
         //- Return the composition of the multi-component mixture
         virtual const basicSpecieMixture& composition() const = 0;
+
+        //- Table of species composition
+        autoPtr<speciesCompositionTable> specieComposition() const
+        {
+            return
+                autoPtr<speciesCompositionTable>
+                (
+                    composition().specieComposition()
+                );
+        }
 };
 
 
diff --git a/src/thermophysicalModels/reactionThermo/rhoReactionThermo/rhoReactionThermo.H b/src/thermophysicalModels/reactionThermo/rhoReactionThermo/rhoReactionThermo.H
index e28274331a843ef8121751e7a2aa9454ad7379b3..72ce0177b6a998c975e8f80651c8b42d8194713a 100644
--- a/src/thermophysicalModels/reactionThermo/rhoReactionThermo/rhoReactionThermo.H
+++ b/src/thermophysicalModels/reactionThermo/rhoReactionThermo/rhoReactionThermo.H
@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2015 OpenFOAM Foundation
-    Copyright (C) 2017 OpenCFD Ltd.
+    Copyright (C) 2017-2020 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -45,12 +45,15 @@ SourceFiles
 #include "basicSpecieMixture.H"
 #include "autoPtr.H"
 #include "runTimeSelectionTables.H"
+#include "specieElement.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 namespace Foam
 {
 
+typedef HashTable<List<specieElement>> speciesCompositionTable;
+
 /*---------------------------------------------------------------------------*\
                      Class rhoReactionThermo Declaration
 \*---------------------------------------------------------------------------*/
@@ -137,6 +140,16 @@ public:
 
         //- Return the composition of the multi-component mixture
         virtual const basicSpecieMixture& composition() const = 0;
+
+        //- Table of species composition
+        autoPtr<speciesCompositionTable> specieComposition() const
+        {
+            return
+                autoPtr<speciesCompositionTable>
+                (
+                    composition().specieComposition()
+                );
+        }
 };
 
 
diff --git a/tutorials/combustion/reactingFoam/laminar/counterFlowFlame2D/0/CH4 b/tutorials/combustion/reactingFoam/laminar/counterFlowFlame2D/0/CH4
index aa5d46ff58bdddf70525aba2ec8a9df4bba35e27..5074f0d0e72f1df823c71e63cf6dba3629061a88 100644
--- a/tutorials/combustion/reactingFoam/laminar/counterFlowFlame2D/0/CH4
+++ b/tutorials/combustion/reactingFoam/laminar/counterFlowFlame2D/0/CH4
@@ -10,34 +10,35 @@ FoamFile
     version     2.0;
     format      ascii;
     class       volScalarField;
-    location    "0";
     object      CH4;
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 dimensions      [0 0 0 0 0 0 0];
 
-internalField   uniform 0.0;
+internalField   uniform 0;
 
 boundaryField
 {
     fuel
     {
         type            fixedValue;
-        value           uniform 1.0;
+        value           uniform 1;
     }
+
     air
     {
         type            fixedValue;
-        value           uniform 0.0;
+        value           uniform 0;
     }
+
     outlet
     {
         type            inletOutlet;
-        inletValue      uniform 0.0;
-        value           uniform 0.0;
-
+        inletValue      uniform 0;
+        value           uniform 0;
     }
+
     frontAndBack
     {
         type            empty;
diff --git a/tutorials/combustion/reactingFoam/laminar/counterFlowFlame2D/0/CO2 b/tutorials/combustion/reactingFoam/laminar/counterFlowFlame2D/0/CO2
index ab2f76c425e0725e05fcfab99ebe15fc1fdb525f..d3ea0e9680fd94ffc1a62d215a1e5f3d211e8260 100644
--- a/tutorials/combustion/reactingFoam/laminar/counterFlowFlame2D/0/CO2
+++ b/tutorials/combustion/reactingFoam/laminar/counterFlowFlame2D/0/CO2
@@ -10,7 +10,6 @@ FoamFile
     version     2.0;
     format      ascii;
     class       volScalarField;
-    location    "0";
     object      CO2;
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -21,23 +20,19 @@ internalField   uniform 0;
 
 boundaryField
 {
-    fuel
-    {
-        type            fixedValue;
-        value           uniform 0;
-    }
-    air
+    "(fuel|air)"
     {
         type            fixedValue;
         value           uniform 0;
     }
+
     outlet
     {
         type            inletOutlet;
         inletValue      uniform 0;
         value           uniform 0;
-
     }
+
     frontAndBack
     {
         type            empty;
diff --git a/tutorials/combustion/reactingFoam/laminar/counterFlowFlame2D/0/H2O b/tutorials/combustion/reactingFoam/laminar/counterFlowFlame2D/0/H2O
index 2ae517b27fa25bf2799eb70819570a2136094274..d0dbfef8872f1df65f7ba526a6d6be1c989ca769 100644
--- a/tutorials/combustion/reactingFoam/laminar/counterFlowFlame2D/0/H2O
+++ b/tutorials/combustion/reactingFoam/laminar/counterFlowFlame2D/0/H2O
@@ -10,7 +10,6 @@ FoamFile
     version     2.0;
     format      ascii;
     class       volScalarField;
-    location    "0";
     object      H2O;
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
diff --git a/tutorials/combustion/reactingFoam/laminar/counterFlowFlame2D/0/N2 b/tutorials/combustion/reactingFoam/laminar/counterFlowFlame2D/0/N2
index 62e1ef202e11536306a5a475e65a4d2c5162d8f0..fb12cb70cba18b90bdbe9b3e24eb39e54b7ab9ad 100644
--- a/tutorials/combustion/reactingFoam/laminar/counterFlowFlame2D/0/N2
+++ b/tutorials/combustion/reactingFoam/laminar/counterFlowFlame2D/0/N2
@@ -10,7 +10,6 @@ FoamFile
     version     2.0;
     format      ascii;
     class       volScalarField;
-    location    "0";
     object      N2;
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -26,18 +25,20 @@ boundaryField
         type            fixedValue;
         value           uniform 0.0;
     }
+
     air
     {
         type            fixedValue;
         value           uniform 0.77;
     }
+
     outlet
     {
         type            inletOutlet;
         inletValue      uniform 1;
         value           uniform 1;
-
     }
+
     frontAndBack
     {
         type            empty;
diff --git a/tutorials/combustion/reactingFoam/laminar/counterFlowFlame2D/0/O2 b/tutorials/combustion/reactingFoam/laminar/counterFlowFlame2D/0/O2
index 506f04ecb9809a8989cae2a1c7c27df2148a1d55..ee7433eaffcc07453cc06146e5dd6d03401a006f 100644
--- a/tutorials/combustion/reactingFoam/laminar/counterFlowFlame2D/0/O2
+++ b/tutorials/combustion/reactingFoam/laminar/counterFlowFlame2D/0/O2
@@ -10,7 +10,6 @@ FoamFile
     version     2.0;
     format      ascii;
     class       volScalarField;
-    location    "0";
     object      O2;
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -24,19 +23,22 @@ boundaryField
     fuel
     {
         type            fixedValue;
-        value           uniform 0.0;
+        value           uniform 0;
     }
+
     air
     {
         type            fixedValue;
         value           uniform 0.23;
     }
+
     outlet
     {
         type            inletOutlet;
         inletValue      uniform 0;
         value           uniform 0;
     }
+
     frontAndBack
     {
         type            empty;
diff --git a/tutorials/combustion/reactingFoam/laminar/counterFlowFlame2D/0/T b/tutorials/combustion/reactingFoam/laminar/counterFlowFlame2D/0/T
index 2a5121fb076a4ae62f15d5931cf321e71c8224c9..31ddd7db75bd2832b3dcfaea98754dae8c12cd5c 100644
--- a/tutorials/combustion/reactingFoam/laminar/counterFlowFlame2D/0/T
+++ b/tutorials/combustion/reactingFoam/laminar/counterFlowFlame2D/0/T
@@ -10,7 +10,6 @@ FoamFile
     version     2.0;
     format      ascii;
     class       volScalarField;
-    location    "0";
     object      T;
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -21,22 +20,19 @@ internalField   uniform 2000;
 
 boundaryField
 {
-    fuel
-    {
-        type            fixedValue;
-        value           uniform 293;
-    }
-    air
+    "(fuel|air)"
     {
         type            fixedValue;
         value           uniform 293;
     }
+
     outlet
     {
         type            inletOutlet;
         inletValue      uniform 293;
         value           uniform 293;
     }
+
     frontAndBack
     {
         type            empty;
diff --git a/tutorials/combustion/reactingFoam/laminar/counterFlowFlame2D/0/U b/tutorials/combustion/reactingFoam/laminar/counterFlowFlame2D/0/U
index 40f4727f524cfc1ba52dda14413f97ac73189792..c192391019f9cf69cca7dcb4278f8de37fe2c6c8 100644
--- a/tutorials/combustion/reactingFoam/laminar/counterFlowFlame2D/0/U
+++ b/tutorials/combustion/reactingFoam/laminar/counterFlowFlame2D/0/U
@@ -10,7 +10,6 @@ FoamFile
     version     2.0;
     format      ascii;
     class       volVectorField;
-    location    "0";
     object      U;
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -26,16 +25,19 @@ boundaryField
         type            fixedValue;
         value           uniform (0.1 0 0);
     }
+
     air
     {
         type            fixedValue;
         value           uniform (-0.1 0 0);
     }
+
     outlet
     {
         type            pressureInletOutletVelocity;
         value           $internalField;
     }
+
     frontAndBack
     {
         type            empty;
diff --git a/tutorials/combustion/reactingFoam/laminar/counterFlowFlame2D/0/Ydefault b/tutorials/combustion/reactingFoam/laminar/counterFlowFlame2D/0/Ydefault
index b13a5a7aa3cd407bff88b5a53ab3f1d6d970e5e8..dfc48588997ba50b646b6cf9f52d2429d08ba86c 100644
--- a/tutorials/combustion/reactingFoam/laminar/counterFlowFlame2D/0/Ydefault
+++ b/tutorials/combustion/reactingFoam/laminar/counterFlowFlame2D/0/Ydefault
@@ -10,34 +10,29 @@ FoamFile
     version     2.0;
     format      ascii;
     class       volScalarField;
-    location    "0";
     object      Ydefault;
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 dimensions      [0 0 0 0 0 0 0];
 
-internalField   uniform 0.0;
+internalField   uniform 0;
 
 boundaryField
 {
-    fuel
+    "(fuel|air)"
     {
         type            fixedValue;
-        value           uniform 0.0;
-    }
-    air
-    {
-        type            fixedValue;
-        value           uniform 0.0;
+        value           uniform 0;
     }
+
     outlet
     {
         type            inletOutlet;
-        inletValue      uniform 0.0;
-        value           uniform 0.0;
-
+        inletValue      uniform 0;
+        value           uniform 0;
     }
+
     frontAndBack
     {
         type            empty;
diff --git a/tutorials/combustion/reactingFoam/laminar/counterFlowFlame2D/0/alphat b/tutorials/combustion/reactingFoam/laminar/counterFlowFlame2D/0/alphat
index d13c439be3cf5bdbf2efaa3b981d677783a9a3aa..2f2a619e46813c94b19ee8089989e8b667b4de4d 100644
--- a/tutorials/combustion/reactingFoam/laminar/counterFlowFlame2D/0/alphat
+++ b/tutorials/combustion/reactingFoam/laminar/counterFlowFlame2D/0/alphat
@@ -10,7 +10,6 @@ FoamFile
     version     2.0;
     format      ascii;
     class       volScalarField;
-    location    "0";
     object      alphat;
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -21,20 +20,17 @@ internalField   uniform 0;
 
 boundaryField
 {
-    fuel
-    {
-        type            fixedValue;
-        value           uniform 0;
-    }
-    air
+    "(fuel|air)"
     {
         type            fixedValue;
         value           uniform 0;
     }
+
     outlet
     {
         type            zeroGradient;
     }
+
     frontAndBack
     {
         type            empty;
diff --git a/tutorials/combustion/reactingFoam/laminar/counterFlowFlame2D/0/p b/tutorials/combustion/reactingFoam/laminar/counterFlowFlame2D/0/p
index 470e37aad8c6df9bae198d63f6e15f76fd960b4c..2d13058f89dd62ca7040200faf53231988eabff0 100644
--- a/tutorials/combustion/reactingFoam/laminar/counterFlowFlame2D/0/p
+++ b/tutorials/combustion/reactingFoam/laminar/counterFlowFlame2D/0/p
@@ -10,7 +10,6 @@ FoamFile
     version     2.0;
     format      ascii;
     class       volScalarField;
-    location    "0";
     object      p;
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -21,19 +20,17 @@ internalField   uniform 1e5;
 
 boundaryField
 {
-    fuel
-    {
-        type            zeroGradient;
-    }
-    air
+    "(fuel|air)"
     {
         type            zeroGradient;
     }
+
     outlet
     {
         type            totalPressure;
         p0              $internalField;
     }
+
     frontAndBack
     {
         type            empty;
diff --git a/tutorials/combustion/reactingFoam/laminar/counterFlowFlame2D/constant/chemistryProperties b/tutorials/combustion/reactingFoam/laminar/counterFlowFlame2D/constant/chemistryProperties
index 39868c8f6335f210d5de890250e275cf809924fb..d48ed1ba9d94d40798783be14fe7d0ba0e83905a 100644
--- a/tutorials/combustion/reactingFoam/laminar/counterFlowFlame2D/constant/chemistryProperties
+++ b/tutorials/combustion/reactingFoam/laminar/counterFlowFlame2D/constant/chemistryProperties
@@ -10,14 +10,13 @@ FoamFile
     version     2.0;
     format      ascii;
     class       dictionary;
-    location    "constant";
     object      chemistryProperties;
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 chemistryType
 {
-    solver            EulerImplicit;
+    solver          EulerImplicit;
 }
 
 chemistry           on;
@@ -37,4 +36,5 @@ odeCoeffs
     relTol          0.01;
 }
 
+
 // ************************************************************************* //
diff --git a/tutorials/combustion/reactingFoam/laminar/counterFlowFlame2D/constant/combustionProperties b/tutorials/combustion/reactingFoam/laminar/counterFlowFlame2D/constant/combustionProperties
index 0d845cb4c5d5b9802689ba0e08e4f076dd0bfbec..8c202ef2e7133978ec382139456e4d468169c3bf 100644
--- a/tutorials/combustion/reactingFoam/laminar/counterFlowFlame2D/constant/combustionProperties
+++ b/tutorials/combustion/reactingFoam/laminar/counterFlowFlame2D/constant/combustionProperties
@@ -10,14 +10,13 @@ FoamFile
     version     2.0;
     format      ascii;
     class       dictionary;
-    location    "constant";
     object      combustionProperties;
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 combustionModel  laminar;
 
-active  true;
+active           true;
 
 laminarCoeffs
 {
diff --git a/tutorials/combustion/reactingFoam/laminar/counterFlowFlame2D/constant/thermo.compressibleGas b/tutorials/combustion/reactingFoam/laminar/counterFlowFlame2D/constant/thermo.compressibleGas
index 88f778a7c08664edf6efac5f46ecf46d8f7f282b..e009ca623e4763bc0333088f109ed4d54dbb3ef1 100644
--- a/tutorials/combustion/reactingFoam/laminar/counterFlowFlame2D/constant/thermo.compressibleGas
+++ b/tutorials/combustion/reactingFoam/laminar/counterFlowFlame2D/constant/thermo.compressibleGas
@@ -10,7 +10,6 @@ FoamFile
     version     2.0;
     format      ascii;
     class       dictionary;
-    location    "constant";
     object      thermo.compressibleGas;
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
diff --git a/tutorials/combustion/reactingFoam/laminar/counterFlowFlame2D/constant/thermophysicalProperties b/tutorials/combustion/reactingFoam/laminar/counterFlowFlame2D/constant/thermophysicalProperties
index 0266f6be0ff18d2bd91ecc0af720beb877a3e7f3..0933dc55d37f89636ab807d02633389f68e26c7a 100644
--- a/tutorials/combustion/reactingFoam/laminar/counterFlowFlame2D/constant/thermophysicalProperties
+++ b/tutorials/combustion/reactingFoam/laminar/counterFlowFlame2D/constant/thermophysicalProperties
@@ -10,7 +10,6 @@ FoamFile
     version     2.0;
     format      ascii;
     class       dictionary;
-    location    "constant";
     object      thermophysicalProperties;
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -26,11 +25,11 @@ thermoType
     specie          specie;
 }
 
-inertSpecie N2;
+inertSpecie         N2;
 
-chemistryReader foamChemistryReader;
+chemistryReader     foamChemistryReader;
 
-foamChemistryFile "<constant>/reactions";
+foamChemistryFile   "<constant>/reactions";
 
 foamChemistryThermoFile "<constant>/thermo.compressibleGas";
 
diff --git a/tutorials/combustion/reactingFoam/laminar/counterFlowFlame2D/constant/turbulenceProperties b/tutorials/combustion/reactingFoam/laminar/counterFlowFlame2D/constant/turbulenceProperties
index fbd1e3c277c13766d96341f6214f8c833ed4b691..c5b843164197cc6c0e83282226fd91201f577e15 100644
--- a/tutorials/combustion/reactingFoam/laminar/counterFlowFlame2D/constant/turbulenceProperties
+++ b/tutorials/combustion/reactingFoam/laminar/counterFlowFlame2D/constant/turbulenceProperties
@@ -10,7 +10,6 @@ FoamFile
     version     2.0;
     format      ascii;
     class       dictionary;
-    location    "constant";
     object      turbulenceProperties;
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
diff --git a/tutorials/combustion/reactingFoam/laminar/counterFlowFlame2D/system/FOBilgerMixtureFraction b/tutorials/combustion/reactingFoam/laminar/counterFlowFlame2D/system/FOBilgerMixtureFraction
new file mode 100644
index 0000000000000000000000000000000000000000..e7fd08212687a9e6f693589e262acf9a7bdba283
--- /dev/null
+++ b/tutorials/combustion/reactingFoam/laminar/counterFlowFlame2D/system/FOBilgerMixtureFraction
@@ -0,0 +1,46 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  v2006                                 |
+|   \\  /    A nd           | Website:  www.openfoam.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+
+BilgerMixtureFraction1
+{
+    // Mandatory entries
+    type            BilgerMixtureFraction;
+    libs            (fieldFunctionObjects);
+
+    fuel
+    {
+        composition moleFractions;
+        CH4         1;
+    }
+
+    oxidiser
+    {
+        composition massFractions;
+        O2          0.23;
+        N2          0.77;
+    }
+
+
+    // Optional entries
+    //phase           "";
+    result          f_Bilger;
+
+    // Optional (inherited) entries
+    region          region0;
+    enabled         true;
+    log             true;
+    timeStart       0;
+    timeEnd         1000;
+    executeControl  timeStep;
+    executeInterval 1;
+    writeControl    writeTime;
+    writeInterval   -1;
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/combustion/reactingFoam/laminar/counterFlowFlame2D/system/controlDict b/tutorials/combustion/reactingFoam/laminar/counterFlowFlame2D/system/controlDict
index 0b250276b04f20a7d58d7d12abd909a053695f72..c39369b9fc3e6270702f6ce5bc9f269cbc9e6f48 100644
--- a/tutorials/combustion/reactingFoam/laminar/counterFlowFlame2D/system/controlDict
+++ b/tutorials/combustion/reactingFoam/laminar/counterFlowFlame2D/system/controlDict
@@ -10,7 +10,6 @@ FoamFile
     version     2.0;
     format      ascii;
     class       dictionary;
-    location    "system";
     object      controlDict;
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -49,5 +48,10 @@ adjustTimeStep  yes;
 
 maxCo           0.4;
 
+functions
+{
+    #include "FOBilgerMixtureFraction"
+}
+
 
 // ************************************************************************* //
diff --git a/tutorials/combustion/reactingFoam/laminar/counterFlowFlame2D/system/fvSchemes b/tutorials/combustion/reactingFoam/laminar/counterFlowFlame2D/system/fvSchemes
index 289266088f7a6bd1ee14061b4d4eb1e9265e269c..d56a5372aba624d13ccc4f3a1ce5811273ec8514 100644
--- a/tutorials/combustion/reactingFoam/laminar/counterFlowFlame2D/system/fvSchemes
+++ b/tutorials/combustion/reactingFoam/laminar/counterFlowFlame2D/system/fvSchemes
@@ -10,7 +10,6 @@ FoamFile
     version     2.0;
     format      ascii;
     class       dictionary;
-    location    "system";
     object      fvSchemes;
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -34,7 +33,7 @@ divSchemes
     div(phi,K)      Gauss limitedLinear 1;
     div(phid,p)     Gauss limitedLinear 1;
     div(phi,epsilon) Gauss limitedLinear 1;
-    div(phi,k) Gauss limitedLinear 1;
+    div(phi,k)      Gauss limitedLinear 1;
     div(((rho*nuEff)*dev2(T(grad(U))))) Gauss linear;
 }
 
diff --git a/tutorials/combustion/reactingFoam/laminar/counterFlowFlame2D/system/fvSolution b/tutorials/combustion/reactingFoam/laminar/counterFlowFlame2D/system/fvSolution
index d5e23ddc42c94aac510773778424436ff5249588..b6da8b0e6b5ed92a0aebcc3b3a080ff965a8381b 100644
--- a/tutorials/combustion/reactingFoam/laminar/counterFlowFlame2D/system/fvSolution
+++ b/tutorials/combustion/reactingFoam/laminar/counterFlowFlame2D/system/fvSolution
@@ -10,7 +10,6 @@ FoamFile
     version     2.0;
     format      ascii;
     class       dictionary;
-    location    "system";
     object      fvSolution;
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -61,7 +60,7 @@ PIMPLE
 {
     momentumPredictor no;
     nOuterCorrectors  1;
-    nCorrectors     2;
+    nCorrectors       2;
     nNonOrthogonalCorrectors 0;
 }