diff --git a/applications/solvers/incompressible/pimpleFoam/Make/options b/applications/solvers/incompressible/pimpleFoam/Make/options
index bc0201f371..20feab872f 100644
--- a/applications/solvers/incompressible/pimpleFoam/Make/options
+++ b/applications/solvers/incompressible/pimpleFoam/Make/options
@@ -8,7 +8,6 @@ EXE_INC = \
     -I$(LIB_SRC)/fvOptions/lnInclude \
     -I$(LIB_SRC)/sampling/lnInclude
 
-
 EXE_LIBS = \
     -lturbulenceModels \
     -lincompressibleTurbulenceModels \
diff --git a/applications/solvers/incompressible/pimpleFoam/pimpleFoam.C b/applications/solvers/incompressible/pimpleFoam/pimpleFoam.C
index 281a9eb991..c37a8d76d9 100644
--- a/applications/solvers/incompressible/pimpleFoam/pimpleFoam.C
+++ b/applications/solvers/incompressible/pimpleFoam/pimpleFoam.C
@@ -30,7 +30,7 @@ Description
 
     Sub-models include:
     - turbulence modelling, i.e. laminar, RAS or LES
-    - run-time selectable finite volume options, e.g. MRF, explicit porosity
+    - run-time selectable MRF and finite volume options, e.g. explicit porosity
 
 \*---------------------------------------------------------------------------*/
 
diff --git a/applications/solvers/incompressible/pisoFoam/Make/options b/applications/solvers/incompressible/pisoFoam/Make/options
index 780009c876..20feab872f 100644
--- a/applications/solvers/incompressible/pisoFoam/Make/options
+++ b/applications/solvers/incompressible/pisoFoam/Make/options
@@ -4,11 +4,15 @@ EXE_INC = \
     -I$(LIB_SRC)/transportModels \
     -I$(LIB_SRC)/transportModels/incompressible/singlePhaseTransportModel \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
-    -I$(LIB_SRC)/meshTools/lnInclude
+    -I$(LIB_SRC)/meshTools/lnInclude \
+    -I$(LIB_SRC)/fvOptions/lnInclude \
+    -I$(LIB_SRC)/sampling/lnInclude
 
 EXE_LIBS = \
     -lturbulenceModels \
     -lincompressibleTurbulenceModels \
     -lincompressibleTransportModels \
     -lfiniteVolume \
-    -lmeshTools
+    -lmeshTools \
+    -lfvOptions \
+    -lsampling
diff --git a/applications/solvers/incompressible/pisoFoam/UEqn.H b/applications/solvers/incompressible/pisoFoam/UEqn.H
new file mode 100644
index 0000000000..50865de6bc
--- /dev/null
+++ b/applications/solvers/incompressible/pisoFoam/UEqn.H
@@ -0,0 +1,23 @@
+// Solve the Momentum equation
+
+MRF.correctBoundaryVelocity(U);
+
+fvVectorMatrix UEqn
+(
+    fvm::ddt(U) + fvm::div(phi, U)
+  + MRF.DDt(U)
+  + turbulence->divDevReff(U)
+ ==
+    fvOptions(U)
+);
+
+UEqn.relax();
+
+fvOptions.constrain(UEqn);
+
+if (piso.momentumPredictor())
+{
+    solve(UEqn == -fvc::grad(p));
+
+    fvOptions.correct(U);
+}
diff --git a/applications/solvers/incompressible/pisoFoam/pEqn.H b/applications/solvers/incompressible/pisoFoam/pEqn.H
new file mode 100644
index 0000000000..5e4d1560ff
--- /dev/null
+++ b/applications/solvers/incompressible/pisoFoam/pEqn.H
@@ -0,0 +1,40 @@
+volScalarField rAU(1.0/UEqn.A());
+
+volVectorField HbyA("HbyA", U);
+HbyA = rAU*UEqn.H();
+surfaceScalarField phiHbyA
+(
+    "phiHbyA",
+    (fvc::interpolate(HbyA) & mesh.Sf())
+  + fvc::interpolate(rAU)*fvc::ddtCorr(U, phi)
+);
+
+MRF.makeRelative(phiHbyA);
+
+adjustPhi(phiHbyA, U, p);
+
+// Non-orthogonal pressure corrector loop
+while (piso.correctNonOrthogonal())
+{
+    // Pressure corrector
+
+    fvScalarMatrix pEqn
+    (
+        fvm::laplacian(rAU, p) == fvc::div(phiHbyA)
+    );
+
+    pEqn.setReference(pRefCell, pRefValue);
+
+    pEqn.solve(mesh.solver(p.select(piso.finalInnerIter())));
+
+    if (piso.finalNonOrthogonalIter())
+    {
+        phi = phiHbyA - pEqn.flux();
+    }
+}
+
+#include "continuityErrs.H"
+
+U = HbyA - rAU*fvc::grad(p);
+U.correctBoundaryConditions();
+fvOptions.correct(U);
diff --git a/applications/solvers/incompressible/pisoFoam/pisoFoam.C b/applications/solvers/incompressible/pisoFoam/pisoFoam.C
index 7bf7abb3ac..482b695cb9 100644
--- a/applications/solvers/incompressible/pisoFoam/pisoFoam.C
+++ b/applications/solvers/incompressible/pisoFoam/pisoFoam.C
@@ -27,7 +27,9 @@ Application
 Description
     Transient solver for incompressible flow.
 
-    Turbulence modelling is generic, i.e. laminar, RAS or LES may be selected.
+    Sub-models include:
+    - turbulence modelling, i.e. laminar, RAS or LES
+    - run-time selectable MRF and finite volume options, e.g. explicit porosity
 
 \*---------------------------------------------------------------------------*/
 
@@ -35,6 +37,7 @@ Description
 #include "singlePhaseTransportModel.H"
 #include "turbulentTransportModel.H"
 #include "pisoControl.H"
+#include "fvIOoptionList.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -47,6 +50,8 @@ int main(int argc, char *argv[])
     pisoControl piso(mesh);
 
     #include "createFields.H"
+    #include "createMRF.H"
+    #include "createFvOptions.H"
     #include "initContinuityErrs.H"
 
     // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -61,61 +66,12 @@ int main(int argc, char *argv[])
 
         // Pressure-velocity PISO corrector
         {
-            // Momentum predictor
-            fvVectorMatrix UEqn
-            (
-                fvm::ddt(U)
-              + fvm::div(phi, U)
-              + turbulence->divDevReff(U)
-            );
-
-            UEqn.relax();
-
-            if (piso.momentumPredictor())
-            {
-                solve(UEqn == -fvc::grad(p));
-            }
+            #include "UEqn.H"
 
             // --- PISO loop
             while (piso.correct())
             {
-                volScalarField rAU(1.0/UEqn.A());
-
-                volVectorField HbyA("HbyA", U);
-                HbyA = rAU*UEqn.H();
-                surfaceScalarField phiHbyA
-                (
-                    "phiHbyA",
-                    (fvc::interpolate(HbyA) & mesh.Sf())
-                  + fvc::interpolate(rAU)*fvc::ddtCorr(U, phi)
-                );
-
-                adjustPhi(phiHbyA, U, p);
-
-                // Non-orthogonal pressure corrector loop
-                while (piso.correctNonOrthogonal())
-                {
-                    // Pressure corrector
-
-                    fvScalarMatrix pEqn
-                    (
-                        fvm::laplacian(rAU, p) == fvc::div(phiHbyA)
-                    );
-
-                    pEqn.setReference(pRefCell, pRefValue);
-
-                    pEqn.solve(mesh.solver(p.select(piso.finalInnerIter())));
-
-                    if (piso.finalNonOrthogonalIter())
-                    {
-                        phi = phiHbyA - pEqn.flux();
-                    }
-                }
-
-                #include "continuityErrs.H"
-
-                U = HbyA - rAU*fvc::grad(p);
-                U.correctBoundaryConditions();
+                #include "pEqn.H"
             }
         }