diff --git a/applications/solvers/compressible/rhoSimpleFoam/pEqn.H b/applications/solvers/compressible/rhoSimpleFoam/pEqn.H
index 54af64c003e36af6fceb68487c0cbcc287a428f2..df120f6dea2720316e6e0e49ae5b4c3b1c5add15 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/incompressible/adjointShapeOptimizationFoam/adjointShapeOptimizationFoam.C b/applications/solvers/incompressible/adjointShapeOptimizationFoam/adjointShapeOptimizationFoam.C
index 0f0c65355fe14f17c9fd0b61d576dd6c28485d76..9f0171a02f269cb99297a6aad30b150b7c6588ca 100644
--- a/applications/solvers/incompressible/adjointShapeOptimizationFoam/adjointShapeOptimizationFoam.C
+++ b/applications/solvers/incompressible/adjointShapeOptimizationFoam/adjointShapeOptimizationFoam.C
@@ -93,7 +93,7 @@ int main(int argc, char *argv[])
         //    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);
diff --git a/src/OpenFOAM/fields/GeometricFields/GeometricField/GeometricField.C b/src/OpenFOAM/fields/GeometricFields/GeometricField/GeometricField.C
index 28d718b67fe6c33bb29029ede46ac8a01e166b8a..8898b162f44dbe436538ed9dddb951cf7d42d3f7 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 79090bdbd9fda97ddb50ff6446b20e4d88ee24cb..aa4676f7e59c8625ca016ec3ac9243e3b02ee569 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 variable relaxation factor for " << name << endl;
+    }
+
+    return fieldRelaxDict_.found(name) || fieldRelaxDict_.found("default");
+}
+
+
+bool Foam::solution::relaxEquation(const word& name) const
 {
     if (debug)
     {
-        Info<< "Find relax for " << name << endl;
+        Info<< "Find equation relaxation factor for " << name << endl;
     }
 
-    return
-        relaxationFactors_.found(name)
-     || relaxationFactors_.found("default");
+    return eqnRelaxDict_.found(name) || eqnRelaxDict_.found("default");
+}
+
+
+Foam::scalar Foam::solution::fieldRelaxationFactor(const word& name) const
+{
+    if (debug)
+    {
+        Info<< "Lookup variable relaxation factor for " << name << endl;
+    }
+
+    if (fieldRelaxDict_.found(name))
+    {
+        return readScalar(fieldRelaxDict_.lookup(name));
+    }
+    else if (fieldRelaxDefault_ > SMALL)
+    {
+        return fieldRelaxDefault_;
+    }
+    else
+    {
+        FatalIOErrorIn
+        (
+            "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::relaxationFactor(const word& name) const
+Foam::scalar Foam::solution::equationRelaxationFactor(const word& name) const
 {
     if (debug)
     {
-        Info<< "Lookup relaxationFactor for " << name << endl;
+        Info<< "Lookup equation relaxation factor for " << name << endl;
     }
 
-    if (relaxationFactors_.found(name))
+    if (eqnRelaxDict_.found(name))
     {
-        return readScalar(relaxationFactors_.lookup(name));
+        return readScalar(eqnRelaxDict_.lookup(name));
     }
-    else if (defaultRelaxationFactor_ > SMALL)
+    else if (eqnRelaxDefault_ > SMALL)
     {
-        return defaultRelaxationFactor_;
+        return eqnRelaxDefault_;
     }
     else
     {
         FatalIOErrorIn
         (
-            "Foam::solution::relaxationFactor(const word&)",
-            relaxationFactors_
-        )   << "Cannot find relaxationFactor for '" << name
+            "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 627e83e49b776d139e2a0e840876b8dc5c8b1f1d..68bb0d8d7db037c9e0565fa9d20eadc3a80ef3cd 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 471c97bf350070ab4f7b06c094aec180eb681080..cc4062356f2867ab3a9f811333dd438f35b38f01 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 db807b5edad2331ea05c01f162592a7b46a44d2a..9080c9b24523aebb4470ec5f770f9c321e7e4544 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 ef5515f3aa049dd0d42fc996278d80a1fdde099c..e1ea4007d539ad2f340543294d3a75b4b19cbba5 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 4fe934adbb6f95bce693d9aaf453b2bbe1e2fc2e..75565551a290e5b4d9a6dd736f73ed8b188457bc 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 0000000000000000000000000000000000000000..a0f0a3579bb2d4709ad80f3dd79c6c88e4e0b88a
--- /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 ecc7350a8427dbaa4eed5f9ca0d91acda26c3cf0..33922b5019b241971720d0452aff117ce2f03253 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));
     }
 }