diff --git a/applications/solvers/compressible/rhoPimpleFoam/rhoPimpleFoam.C b/applications/solvers/compressible/rhoPimpleFoam/rhoPimpleFoam.C
index 046da72269..94f3f54de4 100644
--- a/applications/solvers/compressible/rhoPimpleFoam/rhoPimpleFoam.C
+++ b/applications/solvers/compressible/rhoPimpleFoam/rhoPimpleFoam.C
@@ -68,12 +68,6 @@ int main(int argc, char *argv[])
         // --- Pressure-velocity PIMPLE corrector loop
         for (pimple.start(); pimple.loop(); pimple++)
         {
-            if (pimple.nOuterCorr() != 1)
-            {
-                p.storePrevIter();
-                rho.storePrevIter();
-            }
-
             #include "UEqn.H"
             #include "hEqn.H"
 
diff --git a/applications/solvers/compressible/rhoPimpleFoam/rhoPorousMRFLTSPimpleFoam/rhoPorousMRFLTSPimpleFoam.C b/applications/solvers/compressible/rhoPimpleFoam/rhoPorousMRFLTSPimpleFoam/rhoPorousMRFLTSPimpleFoam.C
index c4cdf7df01..66ffb08275 100644
--- a/applications/solvers/compressible/rhoPimpleFoam/rhoPorousMRFLTSPimpleFoam/rhoPorousMRFLTSPimpleFoam.C
+++ b/applications/solvers/compressible/rhoPimpleFoam/rhoPorousMRFLTSPimpleFoam/rhoPorousMRFLTSPimpleFoam.C
@@ -79,12 +79,6 @@ int main(int argc, char *argv[])
         // --- Pressure-velocity PIMPLE corrector loop
         for (pimple.start(); pimple.loop(); pimple++)
         {
-            if (pimple.nOuterCorr() != 1)
-            {
-                p.storePrevIter();
-                rho.storePrevIter();
-            }
-
             turbulence->correct();
 
             #include "UEqn.H"
diff --git a/applications/solvers/compressible/rhoPimpleFoam/rhoPorousMRFPimpleFoam/rhoPorousMRFPimpleFoam.C b/applications/solvers/compressible/rhoPimpleFoam/rhoPorousMRFPimpleFoam/rhoPorousMRFPimpleFoam.C
index da5d3d886a..5bfbc84d56 100644
--- a/applications/solvers/compressible/rhoPimpleFoam/rhoPorousMRFPimpleFoam/rhoPorousMRFPimpleFoam.C
+++ b/applications/solvers/compressible/rhoPimpleFoam/rhoPorousMRFPimpleFoam/rhoPorousMRFPimpleFoam.C
@@ -73,12 +73,6 @@ int main(int argc, char *argv[])
         // --- Pressure-velocity PIMPLE corrector loop
         for (pimple.start(); pimple.loop(); pimple++)
         {
-            if (pimple.nOuterCorr() != 1)
-            {
-                p.storePrevIter();
-                rho.storePrevIter();
-            }
-
             #include "UEqn.H"
             #include "hEqn.H"
 
diff --git a/applications/solvers/compressible/rhoSimpleFoam/pEqn.H b/applications/solvers/compressible/rhoSimpleFoam/pEqn.H
index 54af64c003..df120f6dea 100644
--- a/applications/solvers/compressible/rhoSimpleFoam/pEqn.H
+++ b/applications/solvers/compressible/rhoSimpleFoam/pEqn.H
@@ -26,7 +26,7 @@ if (simple.transonic())
         );
 
         // Relax the pressure equation to ensure diagonal-dominance
-        pEqn.relax(mesh.relaxationFactor("pEqn"));
+        pEqn.relax(mesh.equationRelaxationFactor("pEqn"));
 
         pEqn.setReference(pRefCell, pRefValue);
 
diff --git a/applications/solvers/compressible/rhoSimpleFoam/rhoPorousMRFSimpleFoam/rhoPorousMRFSimpleFoam.C b/applications/solvers/compressible/rhoSimpleFoam/rhoPorousMRFSimpleFoam/rhoPorousMRFSimpleFoam.C
index 1bfc7d795f..5273f50eaa 100644
--- a/applications/solvers/compressible/rhoSimpleFoam/rhoPorousMRFSimpleFoam/rhoPorousMRFSimpleFoam.C
+++ b/applications/solvers/compressible/rhoSimpleFoam/rhoPorousMRFSimpleFoam/rhoPorousMRFSimpleFoam.C
@@ -59,9 +59,6 @@ int main(int argc, char *argv[])
     {
         Info<< "Time = " << runTime.timeName() << nl << endl;
 
-        p.storePrevIter();
-        rho.storePrevIter();
-
         // Pressure-velocity SIMPLE corrector
         {
             #include "UEqn.H"
diff --git a/applications/solvers/compressible/rhoSimpleFoam/rhoSimpleFoam.C b/applications/solvers/compressible/rhoSimpleFoam/rhoSimpleFoam.C
index 377734b7b4..454da31441 100644
--- a/applications/solvers/compressible/rhoSimpleFoam/rhoSimpleFoam.C
+++ b/applications/solvers/compressible/rhoSimpleFoam/rhoSimpleFoam.C
@@ -55,9 +55,6 @@ int main(int argc, char *argv[])
     {
         Info<< "Time = " << runTime.timeName() << nl << endl;
 
-        p.storePrevIter();
-        rho.storePrevIter();
-
         // Pressure-velocity SIMPLE corrector
         {
             #include "UEqn.H"
diff --git a/applications/solvers/compressible/rhoSimplecFoam/rhoSimplecFoam.C b/applications/solvers/compressible/rhoSimplecFoam/rhoSimplecFoam.C
index 111530979f..278e689bc2 100644
--- a/applications/solvers/compressible/rhoSimplecFoam/rhoSimplecFoam.C
+++ b/applications/solvers/compressible/rhoSimplecFoam/rhoSimplecFoam.C
@@ -57,14 +57,6 @@ int main(int argc, char *argv[])
     {
         Info<< "Time = " << runTime.timeName() << nl << endl;
 
-        p.storePrevIter();
-        rho.storePrevIter();
-
-        if (!simple.transonic())
-        {
-            rho.storePrevIter();
-        }
-
         // Velocity-pressure-enthalpy SIMPLEC corrector
         {
             #include "UEqn.H"
diff --git a/applications/solvers/heatTransfer/buoyantBaffleSimpleFoam/buoyantBaffleSimpleFoam.C b/applications/solvers/heatTransfer/buoyantBaffleSimpleFoam/buoyantBaffleSimpleFoam.C
index 67c1531e9c..1a5b66f0c7 100644
--- a/applications/solvers/heatTransfer/buoyantBaffleSimpleFoam/buoyantBaffleSimpleFoam.C
+++ b/applications/solvers/heatTransfer/buoyantBaffleSimpleFoam/buoyantBaffleSimpleFoam.C
@@ -58,9 +58,6 @@ int main(int argc, char *argv[])
     {
         Info<< "Time = " << runTime.timeName() << nl << endl;
 
-        p_rgh.storePrevIter();
-        rho.storePrevIter();
-
         // Pressure-velocity SIMPLE corrector
         {
             #include "UEqn.H"
diff --git a/applications/solvers/heatTransfer/buoyantBoussinesqPimpleFoam/buoyantBoussinesqPimpleFoam.C b/applications/solvers/heatTransfer/buoyantBoussinesqPimpleFoam/buoyantBoussinesqPimpleFoam.C
index 8d91f0f48b..c00c1fe8d8 100644
--- a/applications/solvers/heatTransfer/buoyantBoussinesqPimpleFoam/buoyantBoussinesqPimpleFoam.C
+++ b/applications/solvers/heatTransfer/buoyantBoussinesqPimpleFoam/buoyantBoussinesqPimpleFoam.C
@@ -81,11 +81,6 @@ int main(int argc, char *argv[])
         // --- Pressure-velocity PIMPLE corrector loop
         for (pimple.start(); pimple.loop(); pimple++)
         {
-            if (pimple.nOuterCorr() != 1)
-            {
-                p_rgh.storePrevIter();
-            }
-
             #include "UEqn.H"
             #include "TEqn.H"
 
diff --git a/applications/solvers/heatTransfer/buoyantBoussinesqSimpleFoam/buoyantBoussinesqSimpleFoam.C b/applications/solvers/heatTransfer/buoyantBoussinesqSimpleFoam/buoyantBoussinesqSimpleFoam.C
index aa8f0aa88e..a3d5c956fc 100644
--- a/applications/solvers/heatTransfer/buoyantBoussinesqSimpleFoam/buoyantBoussinesqSimpleFoam.C
+++ b/applications/solvers/heatTransfer/buoyantBoussinesqSimpleFoam/buoyantBoussinesqSimpleFoam.C
@@ -71,8 +71,6 @@ int main(int argc, char *argv[])
     {
         Info<< "Time = " << runTime.timeName() << nl << endl;
 
-        p_rgh.storePrevIter();
-
         // Pressure-velocity SIMPLE corrector
         {
             #include "UEqn.H"
diff --git a/applications/solvers/heatTransfer/buoyantPimpleFoam/buoyantPimpleFoam.C b/applications/solvers/heatTransfer/buoyantPimpleFoam/buoyantPimpleFoam.C
index 8b6dead722..4239771168 100644
--- a/applications/solvers/heatTransfer/buoyantPimpleFoam/buoyantPimpleFoam.C
+++ b/applications/solvers/heatTransfer/buoyantPimpleFoam/buoyantPimpleFoam.C
@@ -74,11 +74,6 @@ int main(int argc, char *argv[])
         // --- Pressure-velocity PIMPLE corrector loop
         for (pimple.start(); pimple.loop(); pimple++)
         {
-            if (pimple.nOuterCorr() != 1)
-            {
-                p_rgh.storePrevIter();
-            }
-
             #include "UEqn.H"
             #include "hEqn.H"
 
diff --git a/applications/solvers/heatTransfer/buoyantSimpleFoam/buoyantSimpleFoam.C b/applications/solvers/heatTransfer/buoyantSimpleFoam/buoyantSimpleFoam.C
index fdb29b8026..cebe1207d8 100644
--- a/applications/solvers/heatTransfer/buoyantSimpleFoam/buoyantSimpleFoam.C
+++ b/applications/solvers/heatTransfer/buoyantSimpleFoam/buoyantSimpleFoam.C
@@ -56,9 +56,6 @@ int main(int argc, char *argv[])
     {
         Info<< "Time = " << runTime.timeName() << nl << endl;
 
-        p_rgh.storePrevIter();
-        rho.storePrevIter();
-
         // Pressure-velocity SIMPLE corrector
         {
             #include "UEqn.H"
diff --git a/applications/solvers/heatTransfer/buoyantSimpleRadiationFoam/buoyantSimpleRadiationFoam.C b/applications/solvers/heatTransfer/buoyantSimpleRadiationFoam/buoyantSimpleRadiationFoam.C
index 47194e893f..fc6a37ea0f 100644
--- a/applications/solvers/heatTransfer/buoyantSimpleRadiationFoam/buoyantSimpleRadiationFoam.C
+++ b/applications/solvers/heatTransfer/buoyantSimpleRadiationFoam/buoyantSimpleRadiationFoam.C
@@ -59,9 +59,6 @@ int main(int argc, char *argv[])
     {
         Info<< "Time = " << runTime.timeName() << nl << endl;
 
-        p_rgh.storePrevIter();
-        rho.storePrevIter();
-
         // Pressure-velocity SIMPLE corrector
         {
             #include "UEqn.H"
diff --git a/applications/solvers/incompressible/adjointShapeOptimizationFoam/adjointShapeOptimizationFoam.C b/applications/solvers/incompressible/adjointShapeOptimizationFoam/adjointShapeOptimizationFoam.C
index 0f0c65355f..1ba675c4bf 100644
--- a/applications/solvers/incompressible/adjointShapeOptimizationFoam/adjointShapeOptimizationFoam.C
+++ b/applications/solvers/incompressible/adjointShapeOptimizationFoam/adjointShapeOptimizationFoam.C
@@ -85,15 +85,13 @@ int main(int argc, char *argv[])
     {
         Info<< "Time = " << runTime.timeName() << nl << endl;
 
-        p.storePrevIter();
-
         laminarTransport.lookup("lambda") >> lambda;
 
         //alpha +=
         //    mesh.relaxationFactor("alpha")
         //   *(lambda*max(Ua & U, zeroSensitivity) - alpha);
         alpha +=
-            mesh.relaxationFactor("alpha")
+            mesh.fieldRelaxationFactor("alpha")
            *(min(max(alpha + lambda*(Ua & U), zeroAlpha), alphaMax) - alpha);
 
         zeroCells(alpha, inletCells);
@@ -148,8 +146,6 @@ int main(int argc, char *argv[])
             U.correctBoundaryConditions();
         }
 
-        pa.storePrevIter();
-
         // Adjoint Pressure-velocity SIMPLE corrector
         {
             // Adjoint Momentum predictor
diff --git a/applications/solvers/incompressible/pimpleFoam/SRFPimpleFoam/SRFPimpleFoam.C b/applications/solvers/incompressible/pimpleFoam/SRFPimpleFoam/SRFPimpleFoam.C
index ff12ed3f33..30725cd070 100644
--- a/applications/solvers/incompressible/pimpleFoam/SRFPimpleFoam/SRFPimpleFoam.C
+++ b/applications/solvers/incompressible/pimpleFoam/SRFPimpleFoam/SRFPimpleFoam.C
@@ -67,11 +67,6 @@ int main(int argc, char *argv[])
         // --- Pressure-velocity PIMPLE corrector loop
         for (pimple.start(); pimple.loop(); pimple++)
         {
-            if (pimple.nOuterCorr() != 1)
-            {
-                p.storePrevIter();
-            }
-
             #include "UrelEqn.H"
 
             // --- PISO loop
diff --git a/applications/solvers/incompressible/pimpleFoam/pimpleDyMFoam/pimpleDyMFoam.C b/applications/solvers/incompressible/pimpleFoam/pimpleDyMFoam/pimpleDyMFoam.C
index 5e686a781a..55ee07dbc7 100644
--- a/applications/solvers/incompressible/pimpleFoam/pimpleDyMFoam/pimpleDyMFoam.C
+++ b/applications/solvers/incompressible/pimpleFoam/pimpleDyMFoam/pimpleDyMFoam.C
@@ -88,11 +88,6 @@ int main(int argc, char *argv[])
         // --- Pressure-velocity PIMPLE corrector loop
         for (pimple.start(); pimple.loop(); pimple++)
         {
-            if (pimple.nOuterCorr() != 1)
-            {
-                p.storePrevIter();
-            }
-
             #include "UEqn.H"
 
             // --- PISO loop
diff --git a/applications/solvers/incompressible/pimpleFoam/pimpleFoam.C b/applications/solvers/incompressible/pimpleFoam/pimpleFoam.C
index 69cb3aaf70..e8b58bbea5 100644
--- a/applications/solvers/incompressible/pimpleFoam/pimpleFoam.C
+++ b/applications/solvers/incompressible/pimpleFoam/pimpleFoam.C
@@ -66,11 +66,6 @@ int main(int argc, char *argv[])
         // --- Pressure-velocity PIMPLE corrector loop
         for (pimple.start(); pimple.loop(); pimple++)
         {
-            if (pimple.nOuterCorr() != 1)
-            {
-                p.storePrevIter();
-            }
-
             #include "UEqn.H"
 
             // --- PISO loop
diff --git a/applications/solvers/incompressible/simpleFoam/MRFSimpleFoam/MRFSimpleFoam.C b/applications/solvers/incompressible/simpleFoam/MRFSimpleFoam/MRFSimpleFoam.C
index 912b0a85de..c0cc824b49 100644
--- a/applications/solvers/incompressible/simpleFoam/MRFSimpleFoam/MRFSimpleFoam.C
+++ b/applications/solvers/incompressible/simpleFoam/MRFSimpleFoam/MRFSimpleFoam.C
@@ -60,8 +60,6 @@ int main(int argc, char *argv[])
     {
         Info<< "Time = " << runTime.timeName() << nl << endl;
 
-        p.storePrevIter();
-
         // --- Pressure-velocity SIMPLE corrector
         {
             #include "UEqn.H"
diff --git a/applications/solvers/incompressible/simpleFoam/SRFSimpleFoam/SRFSimpleFoam.C b/applications/solvers/incompressible/simpleFoam/SRFSimpleFoam/SRFSimpleFoam.C
index 8ecf1cdf1d..7f1fef3946 100644
--- a/applications/solvers/incompressible/simpleFoam/SRFSimpleFoam/SRFSimpleFoam.C
+++ b/applications/solvers/incompressible/simpleFoam/SRFSimpleFoam/SRFSimpleFoam.C
@@ -56,8 +56,6 @@ int main(int argc, char *argv[])
     {
         Info<< "Time = " << runTime.timeName() << nl << endl;
 
-        p.storePrevIter();
-
         // --- Pressure-velocity SIMPLE corrector
         {
             #include "UrelEqn.H"
diff --git a/applications/solvers/incompressible/simpleFoam/porousSimpleFoam/porousSimpleFoam.C b/applications/solvers/incompressible/simpleFoam/porousSimpleFoam/porousSimpleFoam.C
index bdeb2891c8..d4aaaf596a 100644
--- a/applications/solvers/incompressible/simpleFoam/porousSimpleFoam/porousSimpleFoam.C
+++ b/applications/solvers/incompressible/simpleFoam/porousSimpleFoam/porousSimpleFoam.C
@@ -58,8 +58,6 @@ int main(int argc, char *argv[])
     {
         Info<< "Time = " << runTime.timeName() << nl << endl;
 
-        p.storePrevIter();
-
         // Pressure-velocity SIMPLE corrector
         {
             #include "UEqn.H"
diff --git a/applications/solvers/incompressible/simpleFoam/simpleFoam.C b/applications/solvers/incompressible/simpleFoam/simpleFoam.C
index c578a27bf6..ebfd10bba2 100644
--- a/applications/solvers/incompressible/simpleFoam/simpleFoam.C
+++ b/applications/solvers/incompressible/simpleFoam/simpleFoam.C
@@ -54,8 +54,6 @@ int main(int argc, char *argv[])
     {
         Info<< "Time = " << runTime.timeName() << nl << endl;
 
-        p.storePrevIter();
-
         // --- Pressure-velocity SIMPLE corrector
         {
             #include "UEqn.H"
diff --git a/applications/solvers/incompressible/simpleFoam/windSimpleFoam/windSimpleFoam.C b/applications/solvers/incompressible/simpleFoam/windSimpleFoam/windSimpleFoam.C
index b06649b603..d3a8bcda19 100644
--- a/applications/solvers/incompressible/simpleFoam/windSimpleFoam/windSimpleFoam.C
+++ b/applications/solvers/incompressible/simpleFoam/windSimpleFoam/windSimpleFoam.C
@@ -58,8 +58,6 @@ int main(int argc, char *argv[])
     {
         Info<< "Time = " << runTime.timeName() << nl << endl;
 
-        p.storePrevIter();
-
         // Pressure-velocity SIMPLE corrector
         {
             #include "UEqn.H"
diff --git a/applications/solvers/lagrangian/LTSReactingParcelFoam/LTSReactingParcelFoam.C b/applications/solvers/lagrangian/LTSReactingParcelFoam/LTSReactingParcelFoam.C
index fc08c9e5ef..dc642e7c72 100644
--- a/applications/solvers/lagrangian/LTSReactingParcelFoam/LTSReactingParcelFoam.C
+++ b/applications/solvers/lagrangian/LTSReactingParcelFoam/LTSReactingParcelFoam.C
@@ -89,11 +89,6 @@ int main(int argc, char *argv[])
         // --- Pressure-velocity PIMPLE corrector loop
         for (pimple.start(); pimple.loop(); pimple++)
         {
-            if (pimple.nOuterCorr() != 1)
-            {
-                p.storePrevIter();
-            }
-
             turbulence->correct();
 
             #include "UEqn.H"
diff --git a/applications/solvers/multiphase/bubbleFoam/bubbleFoam.C b/applications/solvers/multiphase/bubbleFoam/bubbleFoam.C
index f6c3e9c3d7..0d095d3648 100644
--- a/applications/solvers/multiphase/bubbleFoam/bubbleFoam.C
+++ b/applications/solvers/multiphase/bubbleFoam/bubbleFoam.C
@@ -70,11 +70,6 @@ int main(int argc, char *argv[])
         // --- Pressure-velocity PIMPLE corrector loop
         for (pimple.start(); pimple.loop(); pimple++)
         {
-            if (pimple.nOuterCorr() != 1)
-            {
-                p.storePrevIter();
-            }
-
             #include "alphaEqn.H"
             #include "liftDragCoeffs.H"
             #include "UEqns.H"
diff --git a/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/compressibleTwoPhaseEulerFoam.C b/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/compressibleTwoPhaseEulerFoam.C
index cec5cbf313..c764472cb7 100644
--- a/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/compressibleTwoPhaseEulerFoam.C
+++ b/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/compressibleTwoPhaseEulerFoam.C
@@ -79,11 +79,6 @@ int main(int argc, char *argv[])
         // --- Pressure-velocity PIMPLE corrector loop
         for (pimple.start(); pimple.loop(); pimple++)
         {
-            if (pimple.nOuterCorr() != 1)
-            {
-                p.storePrevIter();
-            }
-
             #include "alphaEqn.H"
             #include "kEpsilon.H"
             #include "interfacialCoeffs.H"
diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/multiphaseEulerFoam.C b/applications/solvers/multiphase/multiphaseEulerFoam/multiphaseEulerFoam.C
index 814b4dcb6d..b4cc579835 100644
--- a/applications/solvers/multiphase/multiphaseEulerFoam/multiphaseEulerFoam.C
+++ b/applications/solvers/multiphase/multiphaseEulerFoam/multiphaseEulerFoam.C
@@ -77,11 +77,6 @@ int main(int argc, char *argv[])
         // --- Pressure-velocity PIMPLE corrector loop
         for (pimple.start(); pimple.loop(); pimple++)
         {
-            if (pimple.nOuterCorr() != 1)
-            {
-                p.storePrevIter();
-            }
-
             sgsModel->correct();
             fluid.solve();
             rho = fluid.rho();
diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseEulerFoam.C b/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseEulerFoam.C
index 5fb4b7f7f1..43f353cd1c 100644
--- a/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseEulerFoam.C
+++ b/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseEulerFoam.C
@@ -78,11 +78,6 @@ int main(int argc, char *argv[])
         // --- Pressure-velocity PIMPLE corrector loop
         for (pimple.start(); pimple.loop(); pimple++)
         {
-            if (pimple.nOuterCorr() != 1)
-            {
-                p.storePrevIter();
-            }
-
             #include "alphaEqn.H"
             #include "liftDragCoeffs.H"
             #include "UEqns.H"
diff --git a/src/OpenFOAM/fields/GeometricFields/GeometricField/GeometricField.C b/src/OpenFOAM/fields/GeometricFields/GeometricField/GeometricField.C
index 28d718b67f..8898b162f4 100644
--- a/src/OpenFOAM/fields/GeometricFields/GeometricField/GeometricField.C
+++ b/src/OpenFOAM/fields/GeometricFields/GeometricField/GeometricField.C
@@ -919,9 +919,9 @@ void Foam::GeometricField<Type, PatchField, GeoMesh>::relax()
         name += "Final";
     }
 
-    if (this->mesh().relax(name))
+    if (this->mesh().relaxField(name))
     {
-        relax(this->mesh().relaxationFactor(name));
+        relax(this->mesh().fieldRelaxationFactor(name));
     }
 }
 
diff --git a/src/OpenFOAM/matrices/solution/solution.C b/src/OpenFOAM/matrices/solution/solution.C
index 79090bdbd9..aa4676f7e5 100644
--- a/src/OpenFOAM/matrices/solution/solution.C
+++ b/src/OpenFOAM/matrices/solution/solution.C
@@ -56,10 +56,52 @@ void Foam::solution::read(const dictionary& dict)
 
     if (dict.found("relaxationFactors"))
     {
-        relaxationFactors_ = dict.subDict("relaxationFactors");
+        const dictionary& relaxDict(dict.subDict("relaxationFactors"));
+        if (relaxDict.found("variables"))
+        {
+            fieldRelaxDict_ = relaxDict.subDict("variables");
+            eqnRelaxDict_ = relaxDict.subDict("equations");
+        }
+        else
+        {
+            // backwards compatibility
+            const wordList entryNames(relaxDict.toc());
+            forAll(entryNames, i)
+            {
+                const word& e = entryNames[i];
+                scalar value = readScalar(relaxDict.lookup(e));
+
+                if (e(0, 1) == "p")
+                {
+                    fieldRelaxDict_.add(e, value);
+                }
+                else if (e.length() >= 3)
+                {
+                    if (e(0, 3) == "rho")
+                    {
+                        fieldRelaxDict_.add(e, value);
+                    }
+                }
+
+            }
+
+            eqnRelaxDict_ = relaxDict;
+        }
+
+        fieldRelaxDefault_ =
+            fieldRelaxDict_.lookupOrDefault<scalar>("default", 0.0);
+
+        eqnRelaxDefault_ =
+            eqnRelaxDict_.lookupOrDefault<scalar>("default", 0.0);
+
+        if (debug)
+        {
+            Info<< "relaxation factors:" << nl
+                << "fields: " << fieldRelaxDict_ << nl
+                << "equations: " << eqnRelaxDict_ << endl;
+        }
     }
 
-    relaxationFactors_.readIfPresent("default", defaultRelaxationFactor_);
 
     if (dict.found("solvers"))
     {
@@ -92,11 +134,13 @@ Foam::solution::solution
             IOobject::NO_WRITE
         )
     ),
-    cache_(ITstream("cache", tokenList())()),
+    cache_(dictionary::null),
     caching_(false),
-    relaxationFactors_(ITstream("relaxationFactors", tokenList())()),
-    defaultRelaxationFactor_(0),
-    solvers_(ITstream("solvers", tokenList())())
+    fieldRelaxDict_(dictionary::null),
+    eqnRelaxDict_(dictionary::null),
+    fieldRelaxDefault_(0),
+    eqnRelaxDefault_(0),
+    solvers_(dictionary::null)
 {
     if
     (
@@ -207,41 +251,80 @@ bool Foam::solution::cache(const word& name) const
 }
 
 
-bool Foam::solution::relax(const word& name) const
+bool Foam::solution::relaxField(const word& name) const
 {
     if (debug)
     {
-        Info<< "Find relax for " << name << endl;
+        Info<< "Find variable relaxation factor for " << name << endl;
     }
 
-    return
-        relaxationFactors_.found(name)
-     || relaxationFactors_.found("default");
+    return fieldRelaxDict_.found(name) || fieldRelaxDict_.found("default");
 }
 
 
-Foam::scalar Foam::solution::relaxationFactor(const word& name) const
+bool Foam::solution::relaxEquation(const word& name) const
 {
     if (debug)
     {
-        Info<< "Lookup relaxationFactor for " << name << endl;
+        Info<< "Find equation relaxation factor for " << name << endl;
     }
 
-    if (relaxationFactors_.found(name))
+    return eqnRelaxDict_.found(name) || eqnRelaxDict_.found("default");
+}
+
+
+Foam::scalar Foam::solution::fieldRelaxationFactor(const word& name) const
+{
+    if (debug)
     {
-        return readScalar(relaxationFactors_.lookup(name));
+        Info<< "Lookup variable relaxation factor for " << name << endl;
     }
-    else if (defaultRelaxationFactor_ > SMALL)
+
+    if (fieldRelaxDict_.found(name))
     {
-        return defaultRelaxationFactor_;
+        return readScalar(fieldRelaxDict_.lookup(name));
+    }
+    else if (fieldRelaxDefault_ > SMALL)
+    {
+        return fieldRelaxDefault_;
     }
     else
     {
         FatalIOErrorIn
         (
-            "Foam::solution::relaxationFactor(const word&)",
-            relaxationFactors_
-        )   << "Cannot find relaxationFactor for '" << name
+            "Foam::solution::fieldRelaxationFactor(const word&)",
+            fieldRelaxDict_
+        )   << "Cannot find variable relaxation factor for '" << name
+            << "' or a suitable default value."
+            << exit(FatalIOError);
+
+        return 0;
+    }
+}
+
+
+Foam::scalar Foam::solution::equationRelaxationFactor(const word& name) const
+{
+    if (debug)
+    {
+        Info<< "Lookup equation relaxation factor for " << name << endl;
+    }
+
+    if (eqnRelaxDict_.found(name))
+    {
+        return readScalar(eqnRelaxDict_.lookup(name));
+    }
+    else if (eqnRelaxDefault_ > SMALL)
+    {
+        return eqnRelaxDefault_;
+    }
+    else
+    {
+        FatalIOErrorIn
+        (
+            "Foam::solution::eqnRelaxationFactor(const word&)",
+            eqnRelaxDict_
+        )   << "Cannot find equation relaxation factor for '" << name
             << "' or a suitable default value."
             << exit(FatalIOError);
 
diff --git a/src/OpenFOAM/matrices/solution/solution.H b/src/OpenFOAM/matrices/solution/solution.H
index 627e83e49b..68bb0d8d7d 100644
--- a/src/OpenFOAM/matrices/solution/solution.H
+++ b/src/OpenFOAM/matrices/solution/solution.H
@@ -59,10 +59,16 @@ class solution
         bool caching_;
 
         //- Dictionary of relaxation factors for all the fields
-        dictionary relaxationFactors_;
+        dictionary fieldRelaxDict_;
+
+        //- Dictionary of relaxation factors for all the equations
+        dictionary eqnRelaxDict_;
 
         //- Optional default relaxation factor for all the fields
-        scalar defaultRelaxationFactor_;
+        scalar fieldRelaxDefault_;
+
+        //- Optional default relaxation factor for all the equations
+        scalar eqnRelaxDefault_;
 
         //- Dictionary of solver parameters for all the fields
         dictionary solvers_;
@@ -107,10 +113,16 @@ public:
             bool cache(const word& name) const;
 
             //- Return true if the relaxation factor is given for the field
-            bool relax(const word& name) const;
+            bool relaxField(const word& name) const;
+
+            //- Return true if the relaxation factor is given for the equation
+            bool relaxEquation(const word& name) const;
 
             //- Return the relaxation factor for the given field
-            scalar relaxationFactor(const word& name) const;
+            scalar fieldRelaxationFactor(const word& name) const;
+
+            //- Return the relaxation factor for the given eqation
+            scalar equationRelaxationFactor(const word& name) const;
 
             //- Return the selected sub-dictionary of solvers if the "select"
             //  keyword is given, otherwise return the complete dictionary
@@ -122,6 +134,7 @@ public:
             //- Return the solver controls dictionary for the given field
             const dictionary& solver(const word& name) const;
 
+
         // Read
 
             //- Read the solution dictionary
diff --git a/src/finiteVolume/cfdTools/general/solutionControl/pimpleControl/pimpleControlI.H b/src/finiteVolume/cfdTools/general/solutionControl/pimpleControl/pimpleControlI.H
index 471c97bf35..cc4062356f 100644
--- a/src/finiteVolume/cfdTools/general/solutionControl/pimpleControl/pimpleControlI.H
+++ b/src/finiteVolume/cfdTools/general/solutionControl/pimpleControl/pimpleControlI.H
@@ -75,6 +75,11 @@ inline bool Foam::pimpleControl::loop()
                 Info<< algorithmName_ << ": iteration " << corr_ + 1 << endl;
             }
 
+            if (corr_ != 0)
+            {
+                storePrevIterFields();
+            }
+
             return true;
         }
         else
diff --git a/src/finiteVolume/cfdTools/general/solutionControl/simpleControl/simpleControlI.H b/src/finiteVolume/cfdTools/general/solutionControl/simpleControl/simpleControlI.H
index db807b5eda..9080c9b245 100644
--- a/src/finiteVolume/cfdTools/general/solutionControl/simpleControl/simpleControlI.H
+++ b/src/finiteVolume/cfdTools/general/solutionControl/simpleControl/simpleControlI.H
@@ -43,12 +43,18 @@ inline bool Foam::simpleControl::loop()
             // Set to finalise calculation
             time.writeAndEnd();
         }
+        else
+        {
+            storePrevIterFields();
+        }
     }
     else
     {
         initialised_ = true;
+        storePrevIterFields();
     }
 
+
     return time.loop();
 }
 
diff --git a/src/finiteVolume/cfdTools/general/solutionControl/solutionControl/solutionControl.C b/src/finiteVolume/cfdTools/general/solutionControl/solutionControl/solutionControl.C
index ef5515f3aa..e1ea4007d5 100644
--- a/src/finiteVolume/cfdTools/general/solutionControl/solutionControl/solutionControl.C
+++ b/src/finiteVolume/cfdTools/general/solutionControl/solutionControl/solutionControl.C
@@ -102,6 +102,17 @@ Foam::label Foam::solutionControl::applyToField(const word& fieldName) const
 }
 
 
+void Foam::solutionControl::storePrevIterFields() const
+{
+//    storePrevIter<label>();
+    storePrevIter<scalar>();
+    storePrevIter<vector>();
+    storePrevIter<sphericalTensor>();
+    storePrevIter<symmTensor>();
+    storePrevIter<tensor>();
+}
+
+
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 Foam::solutionControl::solutionControl(fvMesh& mesh, const word& algorithmName)
diff --git a/src/finiteVolume/cfdTools/general/solutionControl/solutionControl/solutionControl.H b/src/finiteVolume/cfdTools/general/solutionControl/solutionControl/solutionControl.H
index 4fe934adbb..75565551a2 100644
--- a/src/finiteVolume/cfdTools/general/solutionControl/solutionControl/solutionControl.H
+++ b/src/finiteVolume/cfdTools/general/solutionControl/solutionControl/solutionControl.H
@@ -93,6 +93,13 @@ protected:
         //- Return true if all convergence checks are satified
         virtual bool criteriaSatisfied() = 0;
 
+        //- Store previous iteration fields
+        virtual void storePrevIterFields() const;
+
+        //- Store previous iteration field for vol<Type>Fields
+        template<class Type>
+        void storePrevIter() const;
+
         //- Disallow default bitwise copy construct
         solutionControl(const solutionControl&);
 
@@ -146,6 +153,12 @@ public:
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
+#ifdef NoRepository
+    #include "solutionControlTemplates.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
 #include "solutionControlI.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
diff --git a/src/finiteVolume/cfdTools/general/solutionControl/solutionControl/solutionControlTemplates.C b/src/finiteVolume/cfdTools/general/solutionControl/solutionControl/solutionControlTemplates.C
new file mode 100644
index 0000000000..a0f0a3579b
--- /dev/null
+++ b/src/finiteVolume/cfdTools/general/solutionControl/solutionControl/solutionControlTemplates.C
@@ -0,0 +1,58 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "GeometricField.H"
+#include "volMesh.H"
+#include "fvPatchField.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+template<class Type>
+void Foam::solutionControl::storePrevIter() const
+{
+    typedef GeometricField<Type, fvPatchField, volMesh> GeoField;
+
+    HashTable<const GeoField*>
+        flds(mesh_.objectRegistry::lookupClass<GeoField>());
+
+    forAllConstIter(typename HashTable<const GeoField*>, flds, iter)
+    {
+        GeoField& fld = const_cast<GeoField&>(*iter());
+
+        if (mesh_.relaxField(fld.name()))
+        {
+            if (debug)
+            {
+                Info<< algorithmName_ << ": storing previous iter for "
+                    << fld.name() << endl;
+            }
+
+            fld.storePrevIter();
+        }
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/src/finiteVolume/fvMatrices/fvMatrix/fvMatrix.C b/src/finiteVolume/fvMatrices/fvMatrix/fvMatrix.C
index ecc7350a84..33922b5019 100644
--- a/src/finiteVolume/fvMatrices/fvMatrix/fvMatrix.C
+++ b/src/finiteVolume/fvMatrices/fvMatrix/fvMatrix.C
@@ -627,9 +627,9 @@ void Foam::fvMatrix<Type>::relax()
         ("finalIteration", false)
     );
 
-    if (psi_.mesh().relax(name))
+    if (psi_.mesh().relaxEquation(name))
     {
-        relax(psi_.mesh().relaxationFactor(name));
+        relax(psi_.mesh().equationRelaxationFactor(name));
     }
 }