diff --git a/applications/solvers/combustion/PDRFoam/Make/options b/applications/solvers/combustion/PDRFoam/Make/options
index 4768728619..a415be99b6 100644
--- a/applications/solvers/combustion/PDRFoam/Make/options
+++ b/applications/solvers/combustion/PDRFoam/Make/options
@@ -16,7 +16,8 @@ EXE_INC = \
     -I$(LIB_SRC)/dynamicFvMesh/lnInclude \
     -I$(LIB_SRC)/dynamicMesh/lnInclude \
     -I$(LIB_SRC)/meshTools/lnInclude \
-    -I$(LIB_SRC)/triSurface/lnInclude
+    -I$(LIB_SRC)/triSurface/lnInclude \
+    -I$(LIB_SRC)/fvOptions/lnInclude
 
 EXE_LIBS = \
     -lengine \
@@ -29,4 +30,5 @@ EXE_LIBS = \
     -lspecie \
     -llaminarFlameSpeedModels \
     -lfiniteVolume \
-    -ldynamicFvMesh
+    -ldynamicFvMesh \
+    -lfvOptions
diff --git a/applications/solvers/lagrangian/DPMFoam/DPMTurbulenceModels/Make/options b/applications/solvers/lagrangian/DPMFoam/DPMTurbulenceModels/Make/options
index 902afe59e4..50f1c4f124 100644
--- a/applications/solvers/lagrangian/DPMFoam/DPMTurbulenceModels/Make/options
+++ b/applications/solvers/lagrangian/DPMFoam/DPMTurbulenceModels/Make/options
@@ -7,4 +7,5 @@ EXE_INC = \
     -I$(LIB_SRC)/TurbulenceModels/incompressible/lnInclude \
     -I$(LIB_SRC)/TurbulenceModels/phaseIncompressible/lnInclude \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
-    -I$(LIB_SRC)/meshTools/lnInclude
+    -I$(LIB_SRC)/meshTools/lnInclude \
+    -I$(LIB_SRC)/fvOptions/lnInclude
diff --git a/src/TurbulenceModels/compressible/Make/options b/src/TurbulenceModels/compressible/Make/options
index fa99ca3a58..d0c99293d1 100644
--- a/src/TurbulenceModels/compressible/Make/options
+++ b/src/TurbulenceModels/compressible/Make/options
@@ -6,7 +6,8 @@ EXE_INC = \
     -I$(LIB_SRC)/thermophysicalModels/solidThermo/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/solidSpecie/lnInclude \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
-    -I$(LIB_SRC)/meshTools/lnInclude
+    -I$(LIB_SRC)/meshTools/lnInclude \
+    -I$(LIB_SRC)/fvOptions/lnInclude
 
 LIB_LIBS = \
     -lcompressibleTransportModels \
@@ -16,4 +17,5 @@ LIB_LIBS = \
     -lturbulenceModels \
     -lspecie \
     -lfiniteVolume \
-    -lmeshTools
+    -lmeshTools \
+    -lfvOptions
diff --git a/src/TurbulenceModels/incompressible/Make/options b/src/TurbulenceModels/incompressible/Make/options
index 8eceaf533f..3d11423ccd 100644
--- a/src/TurbulenceModels/incompressible/Make/options
+++ b/src/TurbulenceModels/incompressible/Make/options
@@ -2,10 +2,12 @@ EXE_INC = \
     -I../turbulenceModels/lnInclude \
     -I$(LIB_SRC)/transportModels \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
-    -I$(LIB_SRC)/meshTools/lnInclude
+    -I$(LIB_SRC)/meshTools/lnInclude \
+    -I$(LIB_SRC)/fvOptions/lnInclude
 
 LIB_LIBS = \
     -lincompressibleTransportModels \
     -lturbulenceModels \
     -lfiniteVolume \
-    -lmeshTools
+    -lmeshTools \
+    -lfvOptions
diff --git a/src/TurbulenceModels/turbulenceModels/RAS/kEpsilon/kEpsilon.C b/src/TurbulenceModels/turbulenceModels/RAS/kEpsilon/kEpsilon.C
index 7cf2b1fc10..1deb554d1c 100644
--- a/src/TurbulenceModels/turbulenceModels/RAS/kEpsilon/kEpsilon.C
+++ b/src/TurbulenceModels/turbulenceModels/RAS/kEpsilon/kEpsilon.C
@@ -24,6 +24,7 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "kEpsilon.H"
+#include "fvOptionList.H"
 #include "bound.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -238,6 +239,14 @@ void kEpsilon<BasicTurbulenceModel>::correct()
     const volVectorField& U = this->U_;
     volScalarField& nut = this->nut_;
 
+    // const_cast needed because the operators and functions of fvOptions
+    // are currently non-const.
+    fv::optionList& fvOptions = const_cast<fv::optionList&>
+    (
+        this->mesh_.objectRegistry::template
+            lookupObject<fv::optionList>("fvOptions")
+    );
+
     eddyViscosity<RASModel<BasicTurbulenceModel> >::correct();
 
     volScalarField divU(fvc::div(fvc::absolute(this->phi(), U)));
@@ -260,11 +269,14 @@ void kEpsilon<BasicTurbulenceModel>::correct()
       - fvm::SuSp(((2.0/3.0)*C1_ + C3_)*alpha*rho*divU, epsilon_)
       - fvm::Sp(C2_*alpha*rho*epsilon_/k_, epsilon_)
       + epsilonSource()
+      + fvOptions(alpha, rho, epsilon_)
     );
 
     epsEqn().relax();
+    fvOptions.constrain(epsEqn());
     epsEqn().boundaryManipulate(epsilon_.boundaryField());
     solve(epsEqn);
+    fvOptions.correct(epsilon_);
     bound(epsilon_, this->epsilonMin_);
 
     // Turbulent kinetic energy equation
@@ -278,13 +290,17 @@ void kEpsilon<BasicTurbulenceModel>::correct()
       - fvm::SuSp((2.0/3.0)*alpha*rho*divU, k_)
       - fvm::Sp(alpha*rho*epsilon_/k_, k_)
       + kSource()
+      + fvOptions(alpha, rho, k_)
     );
 
     kEqn().relax();
+    fvOptions.constrain(kEqn());
     solve(kEqn);
+    fvOptions.correct(k_);
     bound(k_, this->kMin_);
 
     correctNut();
+    fvOptions.correct(nut);
 }
 
 
diff --git a/src/fvOptions/fvOption/fvOptionList.H b/src/fvOptions/fvOption/fvOptionList.H
index f6e36d9a41..c881fb01b8 100644
--- a/src/fvOptions/fvOption/fvOptionList.H
+++ b/src/fvOptions/fvOption/fvOptionList.H
@@ -38,6 +38,7 @@ SourceFile
 #include "fvOption.H"
 #include "PtrList.H"
 #include "GeometricField.H"
+#include "geometricOneField.H"
 #include "fvPatchField.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -173,6 +174,33 @@ public:
                 const word& fieldName
             );
 
+            //- Return source for equation
+            template<class Type>
+            tmp<fvMatrix<Type> > operator()
+            (
+                const volScalarField& alpha,
+                const geometricOneField& rho,
+                GeometricField<Type, fvPatchField, volMesh>& field
+            );
+
+            //- Return source for equation
+            template<class Type>
+            tmp<fvMatrix<Type> > operator()
+            (
+                const geometricOneField& alpha,
+                const volScalarField& rho,
+                GeometricField<Type, fvPatchField, volMesh>& field
+            );
+
+            //- Return source for equation
+            template<class Type>
+            tmp<fvMatrix<Type> > operator()
+            (
+                const geometricOneField& alpha,
+                const geometricOneField& rho,
+                GeometricField<Type, fvPatchField, volMesh>& field
+            );
+
 
         // Constraints
 
diff --git a/src/fvOptions/fvOption/fvOptionListTemplates.C b/src/fvOptions/fvOption/fvOptionListTemplates.C
index 21805dd9ff..58e362df07 100644
--- a/src/fvOptions/fvOption/fvOptionListTemplates.C
+++ b/src/fvOptions/fvOption/fvOptionListTemplates.C
@@ -191,6 +191,57 @@ Foam::tmp<Foam::fvMatrix<Type> > Foam::fv::optionList::operator()
 }
 
 
+template<class Type>
+Foam::tmp<Foam::fvMatrix<Type> > Foam::fv::optionList::operator()
+(
+    const geometricOneField& alpha,
+    const geometricOneField& rho,
+    GeometricField<Type, fvPatchField, volMesh>& field
+)
+{
+    return this->operator()(field, field.name());
+}
+
+
+template<class Type>
+Foam::tmp<Foam::fvMatrix<Type> > Foam::fv::optionList::operator()
+(
+    const volScalarField& alpha,
+    const geometricOneField& rho,
+    GeometricField<Type, fvPatchField, volMesh>& field
+)
+{
+    volScalarField one
+    (
+        IOobject
+        (
+            "one",
+            this->mesh_.time().timeName(),
+            this->mesh_,
+            IOobject::NO_READ,
+            IOobject::NO_WRITE,
+            false
+        ),
+        this->mesh_,
+        dimensionedScalar("one", dimless, 1.0)
+    );
+
+    return this->operator()(alpha, one, field, field.name());
+}
+
+
+template<class Type>
+Foam::tmp<Foam::fvMatrix<Type> > Foam::fv::optionList::operator()
+(
+    const geometricOneField& alpha,
+    const volScalarField& rho,
+    GeometricField<Type, fvPatchField, volMesh>& field
+)
+{
+    return this->operator()(rho, field, field.name());
+}
+
+
 template<class Type>
 void Foam::fv::optionList::constrain(fvMatrix<Type>& eqn)
 {