diff --git a/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/createMRF.H b/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/createMRF.H
deleted file mode 100644
index 16c7667be5..0000000000
--- a/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/createMRF.H
+++ /dev/null
@@ -1,3 +0,0 @@
-    IOMRFZoneList MRF(mesh);
-    MRF.correctBoundaryVelocity(U1);
-    MRF.correctBoundaryVelocity(U2);
diff --git a/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/generatePairsAndSubModels b/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/generatePairsAndSubModels
deleted file mode 100644
index e69de29bb2..0000000000
diff --git a/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/phaseCompressibleTurbulenceModels/Make/options b/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/phaseCompressibleTurbulenceModels/Make/options
index 96fecff16d..832cefa73d 100644
--- a/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/phaseCompressibleTurbulenceModels/Make/options
+++ b/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/phaseCompressibleTurbulenceModels/Make/options
@@ -8,4 +8,5 @@ EXE_INC = \
     -I$(LIB_SRC)/TurbulenceModels/compressible/lnInclude \
     -I$(LIB_SRC)/TurbulenceModels/phaseCompressible/lnInclude \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/fvOptions/lnInclude \
     -I$(LIB_SRC)/meshTools/lnInclude
diff --git a/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/phaseSystems/PhaseSystems/HeatAndMassTransferPhaseSystem/HeatAndMassTransferPhaseSystem.C b/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/phaseSystems/PhaseSystems/HeatAndMassTransferPhaseSystem/HeatAndMassTransferPhaseSystem.C
index f003d92b98..ba153f166a 100644
--- a/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/phaseSystems/PhaseSystems/HeatAndMassTransferPhaseSystem/HeatAndMassTransferPhaseSystem.C
+++ b/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/phaseSystems/PhaseSystems/HeatAndMassTransferPhaseSystem/HeatAndMassTransferPhaseSystem.C
@@ -171,13 +171,10 @@ Foam::HeatAndMassTransferPhaseSystem<BasePhaseSystem>::dmdt
 
 template<class BasePhaseSystem>
 Foam::autoPtr<Foam::phaseSystem::momentumTransferTable>
-Foam::HeatAndMassTransferPhaseSystem<BasePhaseSystem>::momentumTransfer
-(
-    IOMRFZoneList& MRF
-) const
+Foam::HeatAndMassTransferPhaseSystem<BasePhaseSystem>::momentumTransfer() const
 {
     autoPtr<phaseSystem::momentumTransferTable>
-        eqnsPtr(BasePhaseSystem::momentumTransfer(MRF));
+        eqnsPtr(BasePhaseSystem::momentumTransfer());
 
     phaseSystem::momentumTransferTable& eqns = eqnsPtr();
 
diff --git a/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/phaseSystems/PhaseSystems/HeatAndMassTransferPhaseSystem/HeatAndMassTransferPhaseSystem.H b/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/phaseSystems/PhaseSystems/HeatAndMassTransferPhaseSystem/HeatAndMassTransferPhaseSystem.H
index 02ec28e08d..27e4f32a77 100644
--- a/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/phaseSystems/PhaseSystems/HeatAndMassTransferPhaseSystem/HeatAndMassTransferPhaseSystem.H
+++ b/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/phaseSystems/PhaseSystems/HeatAndMassTransferPhaseSystem/HeatAndMassTransferPhaseSystem.H
@@ -142,10 +142,8 @@ public:
         ) const;
 
         //- Return the momentum transfer matrices
-        virtual autoPtr<phaseSystem::momentumTransferTable> momentumTransfer
-        (
-            IOMRFZoneList& mrfZones
-        ) const;
+        virtual autoPtr<phaseSystem::momentumTransferTable>
+            momentumTransfer() const;
 
         //- Return the heat transfer matrices
         virtual autoPtr<phaseSystem::heatTransferTable>
diff --git a/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/phaseSystems/PhaseSystems/MomentumTransferPhaseSystem/MomentumTransferPhaseSystem.C b/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/phaseSystems/PhaseSystems/MomentumTransferPhaseSystem/MomentumTransferPhaseSystem.C
index 083af1fed5..9646a6a80c 100644
--- a/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/phaseSystems/PhaseSystems/MomentumTransferPhaseSystem/MomentumTransferPhaseSystem.C
+++ b/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/phaseSystems/PhaseSystems/MomentumTransferPhaseSystem/MomentumTransferPhaseSystem.C
@@ -32,7 +32,6 @@ License
 #include "wallLubricationModel.H"
 #include "turbulentDispersionModel.H"
 
-#include "IOMRFZoneList.H"
 #include "HashPtrTable.H"
 
 #include "fvmDdt.H"
@@ -348,10 +347,7 @@ Foam::MomentumTransferPhaseSystem<BasePhaseSystem>::D
 
 template<class BasePhaseSystem>
 Foam::autoPtr<Foam::phaseSystem::momentumTransferTable>
-Foam::MomentumTransferPhaseSystem<BasePhaseSystem>::momentumTransfer
-(
-    IOMRFZoneList& MRF
-) const
+Foam::MomentumTransferPhaseSystem<BasePhaseSystem>::momentumTransfer() const
 {
     // Create a momentum transfer matrix for each phase
     autoPtr<phaseSystem::momentumTransferTable> eqnsPtr
@@ -452,7 +448,7 @@ Foam::MomentumTransferPhaseSystem<BasePhaseSystem>::momentumTransfer
                   - fvm::Sp(fvc::div(phi), U)
                   - otherPhase->DUDt()
                 )
-              - MRF.DDt(Vm, U - otherPhase->U());
+              - this->MRF_.DDt(Vm, U - otherPhase->U());
 
             Swap(phase, otherPhase);
         }
diff --git a/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/phaseSystems/PhaseSystems/MomentumTransferPhaseSystem/MomentumTransferPhaseSystem.H b/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/phaseSystems/PhaseSystems/MomentumTransferPhaseSystem/MomentumTransferPhaseSystem.H
index 87e180cb4a..5c9f057bed 100644
--- a/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/phaseSystems/PhaseSystems/MomentumTransferPhaseSystem/MomentumTransferPhaseSystem.H
+++ b/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/phaseSystems/PhaseSystems/MomentumTransferPhaseSystem/MomentumTransferPhaseSystem.H
@@ -171,10 +171,8 @@ public:
         virtual tmp<volScalarField> D(const phasePairKey& key) const;
 
         //- Return the momentum transfer matrices
-        virtual autoPtr<phaseSystem::momentumTransferTable> momentumTransfer
-        (
-            IOMRFZoneList& mrfZones
-        ) const;
+        virtual autoPtr<phaseSystem::momentumTransferTable>
+            momentumTransfer() const;
 
         //- Read base phaseProperties dictionary
         virtual bool read();
diff --git a/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/phaseSystems/phaseModel/MovingPhaseModel/MovingPhaseModel.C b/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/phaseSystems/phaseModel/MovingPhaseModel/MovingPhaseModel.C
index 61658e249c..8adb867401 100644
--- a/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/phaseSystems/phaseModel/MovingPhaseModel/MovingPhaseModel.C
+++ b/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/phaseSystems/phaseModel/MovingPhaseModel/MovingPhaseModel.C
@@ -42,12 +42,10 @@ License
 
 template<class BasePhaseModel>
 Foam::tmp<Foam::surfaceScalarField>
-Foam::MovingPhaseModel<BasePhaseModel>::generatePhi
-(
-    const word& phiName,
-    const volVectorField& U
-) const
+Foam::MovingPhaseModel<BasePhaseModel>::phi(const volVectorField& U) const
 {
+    word phiName(IOobject::groupName("phi", this->name()));
+
     IOobject phiHeader
     (
         phiName,
@@ -60,22 +58,21 @@ Foam::MovingPhaseModel<BasePhaseModel>::generatePhi
     {
         Info<< "Reading face flux field " << phiName << endl;
 
-        return
-            tmp<surfaceScalarField>
+        return tmp<surfaceScalarField>
+        (
+            new surfaceScalarField
             (
-                new surfaceScalarField
+                IOobject
                 (
-                    IOobject
-                    (
-                        phiName,
-                        U.mesh().time().timeName(),
-                        U.mesh(),
-                        IOobject::MUST_READ,
-                        IOobject::AUTO_WRITE
-                    ),
-                    U.mesh()
-                )
-            );
+                    phiName,
+                    U.mesh().time().timeName(),
+                    U.mesh(),
+                    IOobject::MUST_READ,
+                    IOobject::AUTO_WRITE
+                ),
+                U.mesh()
+            )
+        );
     }
     else
     {
@@ -142,14 +139,7 @@ Foam::MovingPhaseModel<BasePhaseModel>::MovingPhaseModel
         ),
         fluid.mesh()
     ),
-    phi_
-    (
-        generatePhi
-        (
-            IOobject::groupName("phi", this->name()),
-            U_
-        )
-    ),
+    phi_(phi(U_)),
     alphaPhi_
     (
         IOobject
@@ -225,6 +215,8 @@ void Foam::MovingPhaseModel<BasePhaseModel>::correct()
 {
     BasePhaseModel::correct();
 
+    this->fluid().MRF().correctBoundaryVelocity(U_);
+
     continuityError_ =
         fvc::ddt(*this, this->thermo().rho()) + fvc::div(alphaRhoPhi_);
 }
diff --git a/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/phaseSystems/phaseModel/MovingPhaseModel/MovingPhaseModel.H b/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/phaseSystems/phaseModel/MovingPhaseModel/MovingPhaseModel.H
index b92f259935..f8d32d46b3 100644
--- a/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/phaseSystems/phaseModel/MovingPhaseModel/MovingPhaseModel.H
+++ b/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/phaseSystems/phaseModel/MovingPhaseModel/MovingPhaseModel.H
@@ -91,12 +91,8 @@ class MovingPhaseModel
 
     // Private static member functions
 
-        //- Generate the flux field
-        tmp<surfaceScalarField> generatePhi
-        (
-            const word& phiName,
-            const volVectorField& U
-        ) const;
+        //- Calculate and return the flux field
+        tmp<surfaceScalarField> phi(const volVectorField& U) const;
 
 
 public:
diff --git a/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/phaseSystems/phaseSystem/phaseSystem.C b/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/phaseSystems/phaseSystem/phaseSystem.C
index c89d0eec3e..3bbb6e82de 100644
--- a/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/phaseSystems/phaseSystem/phaseSystem.C
+++ b/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/phaseSystems/phaseSystem/phaseSystem.C
@@ -183,7 +183,10 @@ Foam::phaseSystem::phaseSystem
         ),
         mesh,
         dimensionedScalar("dpdt", dimPressure/dimTime, 0)
-    )
+    ),
+
+    MRF_(mesh_),
+    fvOptions_(mesh_)
 {
     // Blending methods
     forAllConstIter(dictionary, subDict("blending"), iter)
diff --git a/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/phaseSystems/phaseSystem/phaseSystem.H b/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/phaseSystems/phaseSystem/phaseSystem.H
index 46b5ec20e0..e12f52df92 100644
--- a/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/phaseSystems/phaseSystem/phaseSystem.H
+++ b/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/phaseSystems/phaseSystem/phaseSystem.H
@@ -41,6 +41,8 @@ SourceFiles
 #include "phaseModel.H"
 #include "phasePair.H"
 #include "orderedPhasePair.H"
+#include "IOMRFZoneList.H"
+#include "fvIOoptionList.H"
 
 #include "volFields.H"
 #include "surfaceFields.H"
@@ -55,7 +57,6 @@ class blendingMethod;
 template <class modelType> class BlendedInterfacialModel;
 class surfaceTensionModel;
 class aspectRatioModel;
-class IOMRFZoneList;
 
 /*---------------------------------------------------------------------------*\
                          Class phaseSystem Declaration
@@ -174,6 +175,12 @@ protected:
         //- Rate of change of pressure
         volScalarField dpdt_;
 
+        //- Optional MRF zones
+        IOMRFZoneList MRF_;
+
+        //- Optional FV-options
+        fv::IOoptionList fvOptions_;
+
         //- Blending methods
         blendingMethodTable blendingMethods_;
 
@@ -316,10 +323,7 @@ public:
         virtual tmp<volScalarField> dmdt(const phasePairKey& key) const = 0;
 
         //- Return the momentum transfer matrices
-        virtual autoPtr<momentumTransferTable> momentumTransfer
-        (
-            IOMRFZoneList& MRF
-        ) const = 0;
+        virtual autoPtr<momentumTransferTable> momentumTransfer() const = 0;
 
         //- Return the heat transfer matrices
         virtual autoPtr<heatTransferTable> heatTransfer() const = 0;
@@ -363,6 +367,12 @@ public:
             //- Access the rate of change of the pressure
             inline volScalarField& dpdt();
 
+            //- Return MRF zones
+            inline const IOMRFZoneList& MRF() const;
+
+            //- Optional FV-options
+            inline fv::IOoptionList& fvOptions();
+
             //- Access a sub model between a phase pair
             template <class modelType>
             const modelType& lookupSubModel(const phasePair& key) const;
diff --git a/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/phaseSystems/phaseSystem/phaseSystemI.H b/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/phaseSystems/phaseSystem/phaseSystemI.H
index adb3813bbd..0613ba47ac 100644
--- a/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/phaseSystems/phaseSystem/phaseSystemI.H
+++ b/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/phaseSystems/phaseSystem/phaseSystemI.H
@@ -55,4 +55,16 @@ inline Foam::volScalarField& Foam::phaseSystem::dpdt()
 }
 
 
+inline const Foam::IOMRFZoneList& Foam::phaseSystem::MRF() const
+{
+    return MRF_;
+}
+
+
+inline Foam::fv::IOoptionList& Foam::phaseSystem::fvOptions()
+{
+    return fvOptions_;
+}
+
+
 // ************************************************************************* //
diff --git a/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/reactingTwoPhaseEulerFoam.C b/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/reactingTwoPhaseEulerFoam.C
index 5925fc4e24..087434ba5b 100644
--- a/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/reactingTwoPhaseEulerFoam.C
+++ b/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/reactingTwoPhaseEulerFoam.C
@@ -37,7 +37,6 @@ Description
 #include "twoPhaseSystem.H"
 #include "PhaseCompressibleTurbulenceModel.H"
 #include "pimpleControl.H"
-#include "fvIOoptionList.H"
 #include "fixedFluxPressureFvPatchScalarField.H"
 #include "HashPtrTable.H"
 
@@ -53,8 +52,6 @@ int main(int argc, char *argv[])
     pimpleControl pimple(mesh);
 
     #include "createFields.H"
-    #include "createMRF.H"
-    #include "createFvOptions.H"
     #include "initContinuityErrs.H"
     #include "readTimeControls.H"
     #include "CourantNos.H"
@@ -75,6 +72,9 @@ int main(int argc, char *argv[])
 
     #include "pUf/createDDtU.H"
 
+    const IOMRFZoneList& MRF = fluid.MRF();
+    fv::IOoptionList& fvOptions = fluid.fvOptions();
+
     // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
     Info<< "\nStarting time loop\n" << endl;
diff --git a/src/finiteVolume/cfdTools/general/MRF/MRFZoneList.C b/src/finiteVolume/cfdTools/general/MRF/MRFZoneList.C
index fc36a25f34..29c13a6fd2 100644
--- a/src/finiteVolume/cfdTools/general/MRF/MRFZoneList.C
+++ b/src/finiteVolume/cfdTools/general/MRF/MRFZoneList.C
@@ -164,7 +164,7 @@ void Foam::MRFZoneList::addAcceleration
 Foam::tmp<Foam::volVectorField> Foam::MRFZoneList::DDt
 (
     const volVectorField& U
-)
+) const
 {
     tmp<volVectorField> tacceleration
     (
@@ -195,7 +195,7 @@ Foam::tmp<Foam::volVectorField> Foam::MRFZoneList::DDt
 (
     const volScalarField& rho,
     const volVectorField& U
-)
+) const
 {
     return rho*DDt(U);
 }
diff --git a/src/finiteVolume/cfdTools/general/MRF/MRFZoneList.H b/src/finiteVolume/cfdTools/general/MRF/MRFZoneList.H
index 252f01c043..016bbded84 100644
--- a/src/finiteVolume/cfdTools/general/MRF/MRFZoneList.H
+++ b/src/finiteVolume/cfdTools/general/MRF/MRFZoneList.H
@@ -114,14 +114,14 @@ public:
         tmp<volVectorField> DDt
         (
             const volVectorField& U
-        );
+        ) const;
 
         //- Return the frame acceleration
         tmp<volVectorField> DDt
         (
             const volScalarField& rho,
             const volVectorField& U
-        );
+        ) const;
 
         //- Make the given absolute velocity relative within the MRF region
         void makeRelative(volVectorField& U) const;