diff --git a/applications/solvers/compressible/sonicFoam/EEqn.H b/applications/solvers/compressible/sonicFoam/EEqn.H
index f10474c5..73ce8373 100644
--- a/applications/solvers/compressible/sonicFoam/EEqn.H
+++ b/applications/solvers/compressible/sonicFoam/EEqn.H
@@ -1,11 +1,21 @@
 {
-    solve
+    fvScalarMatrix EEqn
     (
         fvm::ddt(rho, e) + fvm::div(phi, e)
       + fvc::ddt(rho, K) + fvc::div(phi, K)
       + fvc::div(fvc::absolute(phi/fvc::interpolate(rho), U), p, "div(phiv,p)")
       - fvm::laplacian(turbulence->alphaEff(), e)
+     ==
+        fvOptions(rho, e)
     );
 
+    EEqn.relax();
+
+    fvOptions.constrain(EEqn);
+
+    EEqn.solve();
+
+    fvOptions.correct(e);
+
     thermo.correct();
 }
diff --git a/applications/solvers/compressible/sonicFoam/Make/options b/applications/solvers/compressible/sonicFoam/Make/options
index f9c097c8..b1663bb8 100644
--- a/applications/solvers/compressible/sonicFoam/Make/options
+++ b/applications/solvers/compressible/sonicFoam/Make/options
@@ -1,6 +1,10 @@
 EXE_INC = \
+    -I$(LIB_SRC)/transportModels/compressible/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
     -I$(LIB_SRC)/turbulenceModels/compressible/turbulenceModel \
+    -I$(LIB_SRC)/meshTools/lnInclude \
+    -I$(LIB_SRC)/sampling/lnInclude \
+    -I$(LIB_SRC)/fvOptions/lnInclude \
     -I$(LIB_SRC)/finiteVolume/lnInclude
 
 EXE_LIBS = \
@@ -9,5 +13,7 @@ EXE_LIBS = \
     -lcompressibleTurbulenceModel \
     -lcompressibleRASModels \
     -lcompressibleLESModels \
-    -lfiniteVolume \
-    -lmeshTools
+    -lmeshTools \
+    -lsampling \
+    -lfvOptions \
+    -lfiniteVolume
diff --git a/applications/solvers/compressible/sonicFoam/UEqn.H b/applications/solvers/compressible/sonicFoam/UEqn.H
index c002490b..c4065161 100644
--- a/applications/solvers/compressible/sonicFoam/UEqn.H
+++ b/applications/solvers/compressible/sonicFoam/UEqn.H
@@ -3,7 +3,18 @@ fvVectorMatrix UEqn
     fvm::ddt(rho, U)
   + fvm::div(phi, U)
   + turbulence->divDevRhoReff(U)
+ ==
+    fvOptions(rho, U)
 );
 
-solve(UEqn == -fvc::grad(p));
-K = 0.5*magSqr(U);
+UEqn.relax();
+
+fvOptions.constrain(UEqn);
+
+if (pimple.momentumPredictor())
+{
+    solve(UEqn == -fvc::grad(p));
+
+    fvOptions.correct(U);
+    K = 0.5*magSqr(U);
+}
diff --git a/applications/solvers/compressible/sonicFoam/pEqn.H b/applications/solvers/compressible/sonicFoam/pEqn.H
index 022abc46..693f9d60 100644
--- a/applications/solvers/compressible/sonicFoam/pEqn.H
+++ b/applications/solvers/compressible/sonicFoam/pEqn.H
@@ -16,6 +16,8 @@ surfaceScalarField phid
     )
 );
 
+fvOptions.makeRelative(fvc::interpolate(psi), phid);
+
 // Non-orthogonal pressure corrector loop
 while (pimple.correctNonOrthogonal())
 {
@@ -24,8 +26,12 @@ while (pimple.correctNonOrthogonal())
         fvm::ddt(psi, p)
       + fvm::div(phid, p)
       - fvm::laplacian(rhorAUf, p)
+     ==
+        fvOptions(psi, p, rho.name())
     );
 
+    fvOptions.constrain(pEqn);
+
     pEqn.solve();
 
     if (pimple.finalNonOrthogonalIter())
@@ -39,4 +45,5 @@ while (pimple.correctNonOrthogonal())
 
 U = HbyA - rAU*fvc::grad(p);
 U.correctBoundaryConditions();
+fvOptions.correct(U);
 K = 0.5*magSqr(U);
diff --git a/applications/solvers/compressible/sonicFoam/sonicDyMFoam/Make/options b/applications/solvers/compressible/sonicFoam/sonicDyMFoam/Make/options
index 263b9861..d374c750 100644
--- a/applications/solvers/compressible/sonicFoam/sonicDyMFoam/Make/options
+++ b/applications/solvers/compressible/sonicFoam/sonicDyMFoam/Make/options
@@ -1,10 +1,13 @@
 EXE_INC = \
     -I.. \
+    -I$(LIB_SRC)/transportModels/compressible/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
     -I$(LIB_SRC)/turbulenceModels/compressible/turbulenceModel \
+    -I$(LIB_SRC)/meshTools/lnInclude \
+    -I$(LIB_SRC)/sampling/lnInclude \
+    -I$(LIB_SRC)/fvOptions/lnInclude \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
-    -I$(LIB_SRC)/dynamicMesh/lnInclude \
-    -I$(LIB_SRC)/meshTools/lnInclude
+    -I$(LIB_SRC)/dynamicMesh/lnInclude
 
 EXE_LIBS = \
     -lfluidThermophysicalModels \
@@ -12,6 +15,8 @@ EXE_LIBS = \
     -lcompressibleTurbulenceModel \
     -lcompressibleRASModels \
     -lcompressibleLESModels \
+    -lmeshTools \
+    -lsampling \
+    -lfvOptions \
     -lfiniteVolume \
-    -ldynamicMesh \
-    -lmeshTools
+    -ldynamicMesh
diff --git a/applications/solvers/compressible/sonicFoam/sonicDyMFoam/pEqn.H b/applications/solvers/compressible/sonicFoam/sonicDyMFoam/pEqn.H
index b375f0e8..76e2a09a 100644
--- a/applications/solvers/compressible/sonicFoam/sonicDyMFoam/pEqn.H
+++ b/applications/solvers/compressible/sonicFoam/sonicDyMFoam/pEqn.H
@@ -19,6 +19,7 @@ surfaceScalarField phid
 );
 
 fvc::makeRelative(phid, psi, U);
+fvOptions.makeRelative(fvc::interpolate(psi), phid);
 
 for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
 {
@@ -27,8 +28,12 @@ for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
         fvm::ddt(psi, p)
       + fvm::div(phid, p)
       - fvm::laplacian(rhorAUf, p)
+     ==
+        fvOptions(psi, p, rho.name())
     );
 
+    fvOptions.constrain(pEqn);
+
     pEqn.solve();
 
     phi = pEqn.flux();
@@ -38,6 +43,8 @@ for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
 
 U = HbyA - rAU*fvc::grad(p);
 U.correctBoundaryConditions();
+fvOptions.correct(U);
+K = 0.5*magSqr(U);
 
 {
     rhoUf = fvc::interpolate(rho*U);
diff --git a/applications/solvers/compressible/sonicFoam/sonicDyMFoam/sonicDyMFoam.C b/applications/solvers/compressible/sonicFoam/sonicDyMFoam/sonicDyMFoam.C
index 62d4716f..746f95d9 100644
--- a/applications/solvers/compressible/sonicFoam/sonicDyMFoam/sonicDyMFoam.C
+++ b/applications/solvers/compressible/sonicFoam/sonicDyMFoam/sonicDyMFoam.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2014 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -35,6 +35,7 @@ Description
 #include "turbulenceModel.H"
 #include "motionSolver.H"
 #include "pimpleControl.H"
+#include "fvIOoptionList.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -44,6 +45,7 @@ int main(int argc, char *argv[])
     #include "createTime.H"
     #include "createMesh.H"
     #include "createFields.H"
+    #include "createFvOptions.H"
     #include "createRhoUf.H"
     #include "initContinuityErrs.H"
 
diff --git a/applications/solvers/compressible/sonicFoam/sonicFoam.C b/applications/solvers/compressible/sonicFoam/sonicFoam.C
index 71d032d4..04aab1d0 100644
--- a/applications/solvers/compressible/sonicFoam/sonicFoam.C
+++ b/applications/solvers/compressible/sonicFoam/sonicFoam.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2014 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -34,6 +34,7 @@ Description
 #include "psiThermo.H"
 #include "turbulenceModel.H"
 #include "pimpleControl.H"
+#include "fvIOoptionList.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -43,6 +44,7 @@ int main(int argc, char *argv[])
     #include "createTime.H"
     #include "createMesh.H"
     #include "createFields.H"
+    #include "createFvOptions.H"
     #include "initContinuityErrs.H"
 
     pimpleControl pimple(mesh);
diff --git a/applications/solvers/compressible/sonicFoam/sonicLiquidFoam/compressibleContinuityErrs.H b/applications/solvers/compressible/sonicFoam/sonicLiquidFoam/compressibleContinuityErrs.H
index b7a816fe..e264660c 100644
--- a/applications/solvers/compressible/sonicFoam/sonicLiquidFoam/compressibleContinuityErrs.H
+++ b/applications/solvers/compressible/sonicFoam/sonicLiquidFoam/compressibleContinuityErrs.H
@@ -1,5 +1,5 @@
 {
-#   include "rhoEqn.H"
+    solve(fvm::ddt(rho) + fvc::div(phi));
 }
 {
     scalar sumLocalContErr =
diff --git a/applications/solvers/compressible/sonicFoam/sonicLiquidFoam/sonicLiquidFoam.C b/applications/solvers/compressible/sonicFoam/sonicLiquidFoam/sonicLiquidFoam.C
index eb0a4a0d..acfb4d3b 100644
--- a/applications/solvers/compressible/sonicFoam/sonicLiquidFoam/sonicLiquidFoam.C
+++ b/applications/solvers/compressible/sonicFoam/sonicLiquidFoam/sonicLiquidFoam.C
@@ -58,7 +58,7 @@ int main(int argc, char *argv[])
         #include "readTimeControls.H"
         #include "compressibleCourantNo.H"
 
-        #include "rhoEqn.H"
+        solve(fvm::ddt(rho) + fvc::div(phi));
 
         // --- Pressure-velocity PIMPLE corrector loop
         while (pimple.loop())
diff --git a/src/finiteVolume/cfdTools/compressible/rhoEqn.H b/src/finiteVolume/cfdTools/compressible/rhoEqn.H
index 25ccb85a..353d027e 100644
--- a/src/finiteVolume/cfdTools/compressible/rhoEqn.H
+++ b/src/finiteVolume/cfdTools/compressible/rhoEqn.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2014 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -30,7 +30,19 @@ Description
 \*---------------------------------------------------------------------------*/
 
 {
-    solve(fvm::ddt(rho) + fvc::div(phi));
+    fvScalarMatrix rhoEqn
+    (
+        fvm::ddt(rho)
+      + fvc::div(phi)
+      ==
+        fvOptions(rho)
+    );
+
+    fvOptions.constrain(rhoEqn);
+
+    rhoEqn.solve();
+
+    fvOptions.correct(rho);
 }
 
 // ************************************************************************* //
diff --git a/src/fvOptions/fvOptions/fvIOoptionList.C b/src/fvOptions/fvOptions/fvIOoptionList.C
index 0a7ab4e6..45da9683 100644
--- a/src/fvOptions/fvOptions/fvIOoptionList.C
+++ b/src/fvOptions/fvOptions/fvIOoptionList.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2014 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -37,7 +37,7 @@ Foam::IOobject Foam::fv::IOoptionList::createIOobject
     IOobject io
     (
         "fvOptions",
-        mesh.time().system(),
+        mesh.time().constant(),
         mesh,
         IOobject::MUST_READ,
         IOobject::NO_WRITE
@@ -45,7 +45,8 @@ Foam::IOobject Foam::fv::IOoptionList::createIOobject
 
     if (io.headerOk())
     {
-        Info<< "Creating finite volume options from " << io.name() << nl
+        Info<< "Creating finite volume options from "
+            << io.instance()/io.name() << nl
             << endl;
 
         io.readOpt() = IOobject::MUST_READ_IF_MODIFIED;
@@ -53,10 +54,25 @@ Foam::IOobject Foam::fv::IOoptionList::createIOobject
     }
     else
     {
-        Info<< "No finite volume options present" << nl << endl;
+        // Check if the fvOptions file is in system
+        io.instance() = mesh.time().system();
 
-        io.readOpt() = IOobject::NO_READ;
-        return io;
+        if (io.headerOk())
+        {
+            Info<< "Creating finite volume options from "
+                << io.instance()/io.name() << nl
+                << endl;
+
+            io.readOpt() = IOobject::MUST_READ_IF_MODIFIED;
+            return io;
+        }
+        else
+        {
+            Info<< "No finite volume options present" << nl << endl;
+
+            io.readOpt() = IOobject::NO_READ;
+            return io;
+        }
     }
 }
 
diff --git a/src/fvOptions/fvOptions/fvOptionList.C b/src/fvOptions/fvOptions/fvOptionList.C
index b3ccc28f..e9a022a3 100644
--- a/src/fvOptions/fvOptions/fvOptionList.C
+++ b/src/fvOptions/fvOptions/fvOptionList.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2014 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -24,10 +24,7 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "fvOptionList.H"
-#include "addToRunTimeSelectionTable.H"
-#include "fvMesh.H"
 #include "surfaceFields.H"
-#include "Time.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
@@ -42,6 +39,37 @@ namespace fv
 
 // * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * //
 
+const Foam::dictionary& Foam::fv::optionList::optionsDict
+(
+    const dictionary& dict
+) const
+{
+    if (dict.found("options"))
+    {
+        return dict.subDict("options");
+    }
+    else
+    {
+        return dict;
+    }
+}
+
+
+bool Foam::fv::optionList::readOptions(const dictionary& dict)
+{
+    checkTimeIndex_ = mesh_.time().timeIndex() + 2;
+
+    bool allOk = true;
+    forAll(*this, i)
+    {
+        option& bs = this->operator[](i);
+        bool ok = bs.read(dict.subDict(bs.name()));
+        allOk = (allOk && ok);
+    }
+    return allOk;
+}
+
+
 void Foam::fv::optionList::checkApplied() const
 {
     if (mesh_.time().timeIndex() == checkTimeIndex_)
@@ -63,7 +91,7 @@ Foam::fv::optionList::optionList(const fvMesh& mesh, const dictionary& dict)
     mesh_(mesh),
     checkTimeIndex_(mesh_.time().startTimeIndex() + 2)
 {
-    reset(dict);
+    reset(optionsDict(dict));
 }
 
 
@@ -171,16 +199,7 @@ void Foam::fv::optionList::makeAbsolute
 
 bool Foam::fv::optionList::read(const dictionary& dict)
 {
-    checkTimeIndex_ = mesh_.time().timeIndex() + 2;
-
-    bool allOk = true;
-    forAll(*this, i)
-    {
-        option& bs = this->operator[](i);
-        bool ok = bs.read(dict.subDict(bs.name()));
-        allOk = (allOk && ok);
-    }
-    return allOk;
+    return readOptions(optionsDict(dict));
 }
 
 
diff --git a/src/fvOptions/fvOptions/fvOptionList.H b/src/fvOptions/fvOptions/fvOptionList.H
index 15e97331..d4d99df3 100644
--- a/src/fvOptions/fvOptions/fvOptionList.H
+++ b/src/fvOptions/fvOptions/fvOptionList.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2014 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -68,6 +68,12 @@ protected:
 
     // Protected Member Functions
 
+        //- Return the "options" sub-dictionary if present otherwise return dict
+        const dictionary& optionsDict(const dictionary& dict) const;
+
+        //- Read options dictionary
+        bool readOptions(const dictionary& dict);
+
         //- Check that all sources have been applied
         void checkApplied() const;
 
diff --git a/tutorials/multiphase/driftFluxFoam/ras/mixerVessel2D/system/fvOptions b/tutorials/multiphase/driftFluxFoam/ras/mixerVessel2D/constant/fvOptions
similarity index 100%
rename from tutorials/multiphase/driftFluxFoam/ras/mixerVessel2D/system/fvOptions
rename to tutorials/multiphase/driftFluxFoam/ras/mixerVessel2D/constant/fvOptions
diff --git a/tutorials/multiphase/interFoam/laminar/mixerVessel2D/system/fvOptions b/tutorials/multiphase/interFoam/laminar/mixerVessel2D/constant/fvOptions
similarity index 100%
rename from tutorials/multiphase/interFoam/laminar/mixerVessel2D/system/fvOptions
rename to tutorials/multiphase/interFoam/laminar/mixerVessel2D/constant/fvOptions
diff --git a/tutorials/multiphase/multiphaseInterFoam/laminar/mixerVessel2D/system/fvOptions b/tutorials/multiphase/multiphaseInterFoam/laminar/mixerVessel2D/constant/fvOptions
similarity index 100%
rename from tutorials/multiphase/multiphaseInterFoam/laminar/mixerVessel2D/system/fvOptions
rename to tutorials/multiphase/multiphaseInterFoam/laminar/mixerVessel2D/constant/fvOptions