diff --git a/applications/solvers/combustion/PDRFoam/EaEqn.H b/applications/solvers/combustion/PDRFoam/EaEqn.H
index 1baaa7180f5d270d7a907f71c9b9c2e3b776d22f..8b844a8159198fcbde6ab1589e6bba45bd92e15e 100644
--- a/applications/solvers/combustion/PDRFoam/EaEqn.H
+++ b/applications/solvers/combustion/PDRFoam/EaEqn.H
@@ -1,7 +1,7 @@
 {
     volScalarField& hea = thermo.he();
 
-    solve
+    fvScalarMatrix EaEqn
     (
         betav*fvm::ddt(rho, hea) + mvConvection->fvmDiv(phi, hea)
       + betav*fvc::ddt(rho, K) + fvc::div(phi, K)
@@ -16,7 +16,16 @@
           : -betav*dpdt
         )
       - fvm::laplacian(Db, hea)
+      + betav*fvOptions(rho, hea)
     );
 
+    EaEqn.relax();
+
+    fvOptions.constrain(EaEqn);
+
+    EaEqn.solve();
+
+    fvOptions.correct(hea);
+
     thermo.correct();
 }
diff --git a/applications/solvers/combustion/PDRFoam/EauEqn.H b/applications/solvers/combustion/PDRFoam/EauEqn.H
index f031d907323602212d4331312b7c8bb62e65e1e2..91bb49b9136f659d4a3b28e07d1ef4fa8179974b 100644
--- a/applications/solvers/combustion/PDRFoam/EauEqn.H
+++ b/applications/solvers/combustion/PDRFoam/EauEqn.H
@@ -2,7 +2,7 @@ if (ign.ignited())
 {
     volScalarField& heau = thermo.heu();
 
-    solve
+    fvScalarMatrix heauEqn
     (
         betav*fvm::ddt(rho, heau) + mvConvection->fvmDiv(phi, heau)
       + (betav*fvc::ddt(rho, K) + fvc::div(phi, K))*rho/thermo.rhou()
@@ -23,5 +23,13 @@ if (ign.ignited())
         // A possible solution would be to solve for ftu as well as ft.
         //- fvm::div(muEff*fvc::grad(b)/(b + 0.001), heau)
         //+ fvm::Sp(fvc::div(muEff*fvc::grad(b)/(b + 0.001)), heau)
+     ==
+        betav*fvOptions(rho, heau)
     );
+
+    fvOptions.constrain(heauEqn);
+
+    heauEqn.solve();
+
+    fvOptions.correct(heau);
 }
diff --git a/applications/solvers/combustion/PDRFoam/PDRFoam.C b/applications/solvers/combustion/PDRFoam/PDRFoam.C
index fa8518ba495f105c3159be9d482b22b892778a46..65f52e2ee4cb4861175b44fca83724d593fc265d 100644
--- a/applications/solvers/combustion/PDRFoam/PDRFoam.C
+++ b/applications/solvers/combustion/PDRFoam/PDRFoam.C
@@ -77,6 +77,7 @@ Description
 #include "Switch.H"
 #include "bound.H"
 #include "pimpleControl.H"
+#include "fvIOoptionList.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -93,11 +94,13 @@ int main(int argc, char *argv[])
     #include "readGravitationalAcceleration.H"
     #include "createFields.H"
     #include "createMRF.H"
+    #include "createFvOptions.H"
     #include "initContinuityErrs.H"
     #include "createTimeControls.H"
     #include "compressibleCourantNo.H"
     #include "setInitialDeltaT.H"
 
+    turbulence->validate();
     scalar StCoNum = 0.0;
 
     // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
diff --git a/applications/solvers/combustion/PDRFoam/PDRFoamAutoRefine.C b/applications/solvers/combustion/PDRFoam/PDRFoamAutoRefine.C
index 057286f1e945dbc0f5edbbac8eabb202d211f3bb..4e561052fd0206b7ec9a3ce82a2c3df2d0f9d7ce 100644
--- a/applications/solvers/combustion/PDRFoam/PDRFoamAutoRefine.C
+++ b/applications/solvers/combustion/PDRFoam/PDRFoamAutoRefine.C
@@ -87,6 +87,7 @@ int main(int argc, char *argv[])
     #include "compressibleCourantNo.H"
     #include "setInitialDeltaT.H"
 
+    turbulence->validate();
     scalar StCoNum = 0.0;
 
     // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
diff --git a/applications/solvers/combustion/PDRFoam/UEqn.H b/applications/solvers/combustion/PDRFoam/UEqn.H
index 1d949d17b7c5e787c6366421cda84d819b4eec72..983e020cab7659458fd951f71a86be60e6d1a441 100644
--- a/applications/solvers/combustion/PDRFoam/UEqn.H
+++ b/applications/solvers/combustion/PDRFoam/UEqn.H
@@ -7,13 +7,17 @@
       + turbulence->divDevRhoReff(U)
      ==
         betav*rho*g
+      + betav*fvOptions(rho, U)
     );
 
+    fvOptions.constrain(UEqn);
+
     volSymmTensorField invA(inv(I*UEqn.A() + drag->Dcu()));
 
     if (pimple.momentumPredictor())
     {
         U = invA & (UEqn.H() - betav*fvc::grad(p));
         U.correctBoundaryConditions();
+        fvOptions.correct(U);
         K = 0.5*magSqr(U);
     }
diff --git a/applications/solvers/combustion/PDRFoam/bEqn.H b/applications/solvers/combustion/PDRFoam/bEqn.H
index abec1938c32ce5a09d0ed2c6591f6026dbdd4c0d..8ad1161871305276d7a68e40a04fdc9d44e12ada 100644
--- a/applications/solvers/combustion/PDRFoam/bEqn.H
+++ b/applications/solvers/combustion/PDRFoam/bEqn.H
@@ -73,6 +73,8 @@ if (ign.ignited())
       + fvm::div(phiSt, b)
       - fvm::Sp(fvc::div(phiSt), b)
       - fvm::laplacian(Db, b)
+     ==
+        betav*fvOptions(rho, b)
     );
 
 
@@ -82,8 +84,14 @@ if (ign.ignited())
 
     // Solve for b
     // ~~~~~~~~~~~
+    bEqn.relax();
+
+    fvOptions.constrain(bEqn);
+
     bEqn.solve();
 
+    fvOptions.correct(b);
+
     Info<< "min(b) = " << min(b).value() << endl;
 
     if (composition.contains("ft"))
diff --git a/applications/solvers/combustion/PDRFoam/pEqn.H b/applications/solvers/combustion/PDRFoam/pEqn.H
index 714a3e84e8b2ad35a26ff5c0e079f6a2b8e52de7..499885875bd859e5c4cd1ff12099c29f11832cf1 100644
--- a/applications/solvers/combustion/PDRFoam/pEqn.H
+++ b/applications/solvers/combustion/PDRFoam/pEqn.H
@@ -25,6 +25,8 @@ if (pimple.transonic())
             betav*fvm::ddt(psi, p)
           + fvm::div(phid, p)
           - fvm::laplacian(rho*invA, p)
+         ==
+            betav*fvOptions(psi, p, rho.name())
         );
 
         pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter())));
@@ -53,6 +55,8 @@ else
             betav*fvm::ddt(psi, p)
           + fvc::div(phiHbyA)
           - fvm::laplacian(rho*invA, p)
+         ==
+            betav*fvOptions(psi, p, rho.name())
         );
 
         pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter())));
@@ -69,6 +73,7 @@ else
 
 U = HbyA - (invA & (betav*fvc::grad(p)));
 U.correctBoundaryConditions();
+fvOptions.correct(U);
 K = 0.5*magSqr(U);
 
 if (thermo.dpdt())