diff --git a/applications/solvers/incompressible/pimpleDyMFoam/UEqn.H b/applications/solvers/incompressible/pimpleDyMFoam/UEqn.H
index 1fd514634..5f4a38c9f 100644
--- a/applications/solvers/incompressible/pimpleDyMFoam/UEqn.H
+++ b/applications/solvers/incompressible/pimpleDyMFoam/UEqn.H
@@ -10,7 +10,7 @@
 
     // Get under-relaxation factor
     scalar UUrf =
-        mesh.solutionDict().relaxationFactor(U.select(pimple.finalIter()));
+        mesh.solutionDict().equationRelaxationFactor(U.select(pimple.finalIter()));
 
     if (pimple.momentumPredictor())
     {
@@ -23,4 +23,4 @@
           - fvc::grad(p),
             mesh.solutionDict().solver((U.select(pimple.finalIter())))
         );
-    }
\ No newline at end of file
+    }
diff --git a/applications/solvers/incompressible/pimpleFoam/UEqn.H b/applications/solvers/incompressible/pimpleFoam/UEqn.H
index ae7f8ce4d..df018319c 100644
--- a/applications/solvers/incompressible/pimpleFoam/UEqn.H
+++ b/applications/solvers/incompressible/pimpleFoam/UEqn.H
@@ -10,7 +10,7 @@ fvVectorMatrix ddtUEqn(fvm::ddt(U));
 
 // Get under-relaxation factor
 scalar UUrf =
-    mesh.solutionDict().relaxationFactor(U.select(pimple.finalIter()));
+    mesh.solutionDict().equationRelaxationFactor(U.select(pimple.finalIter()));
 
 if (pimple.momentumPredictor())
 {
diff --git a/applications/solvers/incompressible/simpleFoam/UEqn.H b/applications/solvers/incompressible/simpleFoam/UEqn.H
index 6bb9b59bd..76c83cdf0 100644
--- a/applications/solvers/incompressible/simpleFoam/UEqn.H
+++ b/applications/solvers/incompressible/simpleFoam/UEqn.H
@@ -7,7 +7,7 @@
     );
 
     // Get under-relaxation factor
-    const scalar UUrf = mesh.solutionDict().relaxationFactor(U.name());
+    const scalar UUrf = mesh.solutionDict().equationRelaxationFactor(U.name());
 
     // Momentum solution
     solve
diff --git a/applications/solvers/incompressible/simpleSRFFoam/UEqn.H b/applications/solvers/incompressible/simpleSRFFoam/UEqn.H
index 070e768d8..fd8521e0c 100644
--- a/applications/solvers/incompressible/simpleSRFFoam/UEqn.H
+++ b/applications/solvers/incompressible/simpleSRFFoam/UEqn.H
@@ -8,7 +8,7 @@
     );
 
     // Get under-relaxation factor
-    const scalar UUrf = mesh.solutionDict().relaxationFactor(Urel.name());
+    const scalar UUrf = mesh.solutionDict().equationRelaxationFactor(Urel.name());
 
     // Momentum solution
     solve(relax(HUrelEqn(), UUrf) == -fvc::grad(p));
diff --git a/src/finiteArea/faMatrices/faMatrix/faMatrix.C b/src/finiteArea/faMatrices/faMatrix/faMatrix.C
index 1f5c9a5f5..cf27c2ef3 100644
--- a/src/finiteArea/faMatrices/faMatrix/faMatrix.C
+++ b/src/finiteArea/faMatrices/faMatrix/faMatrix.C
@@ -548,16 +548,18 @@ void faMatrix<Type>::relax(const scalar alpha)
 template<class Type>
 void faMatrix<Type>::relax()
 {
-    scalar alpha = 0;
-
     if (psi_.mesh().solutionDict().relax(psi_.name()))
     {
-        alpha = psi_.mesh().solutionDict().relaxationFactor(psi_.name());
+        relax(psi_.mesh().solutionDict().equationRelaxationFactor(psi_.name()));
     }
-
-    if (alpha > 0)
+    else
     {
-        relax(alpha);
+        if (debug)
+        {
+            InfoIn("void faMatrix<Type>::relax()")
+                << "Relaxation factor for field " << psi_.name()
+                << " not found.  Relaxation will not be used." << endl;
+        }
     }
 }
 
diff --git a/src/finiteVolume/cfdTools/general/solutionControl/solutionControl/solutionControlTemplates.C b/src/finiteVolume/cfdTools/general/solutionControl/solutionControl/solutionControlTemplates.C
index 011c8af5f..152562216 100644
--- a/src/finiteVolume/cfdTools/general/solutionControl/solutionControl/solutionControlTemplates.C
+++ b/src/finiteVolume/cfdTools/general/solutionControl/solutionControl/solutionControlTemplates.C
@@ -45,7 +45,11 @@ void Foam::solutionControl::storePrevIter() const
 
         size_t prevIterField = fName.find("PrevIter");
 
-        if ((prevIterField == word::npos) && mesh_.solutionDict().relax(fName))
+        if
+        (
+            (prevIterField == word::npos)
+         && mesh_.solutionDict().relaxField(fName)
+        )
         {
             if (debug)
             {
diff --git a/src/finiteVolume/fvMatrices/fvMatrix/fvMatrix.C b/src/finiteVolume/fvMatrices/fvMatrix/fvMatrix.C
index e52be5a91..79802028a 100644
--- a/src/finiteVolume/fvMatrices/fvMatrix/fvMatrix.C
+++ b/src/finiteVolume/fvMatrices/fvMatrix/fvMatrix.C
@@ -696,9 +696,9 @@ void Foam::fvMatrix<Type>::relax(const scalar alpha)
 template<class Type>
 void Foam::fvMatrix<Type>::relax()
 {
-    if (psi_.mesh().solutionDict().relax(psi_.name()))
+    if (psi_.mesh().solutionDict().relaxEquation(psi_.name()))
     {
-        relax(psi_.mesh().solutionDict().relaxationFactor(psi_.name()));
+        relax(psi_.mesh().solutionDict().equationRelaxationFactor(psi_.name()));
     }
     else
     {
diff --git a/src/foam/fields/GeometricFields/GeometricField/GeometricField.C b/src/foam/fields/GeometricFields/GeometricField/GeometricField.C
index d8bcf09d7..1dfc0c2c8 100644
--- a/src/foam/fields/GeometricFields/GeometricField/GeometricField.C
+++ b/src/foam/fields/GeometricFields/GeometricField/GeometricField.C
@@ -666,7 +666,7 @@ Foam::GeometricField<Type, PatchField, GeoMesh>::GeometricField
     if (debug)
     {
         Info<< "GeometricField<Type, PatchField, GeoMesh>::GeometricField : "
-               "constructing as copy resetting IO params"
+               "constructing as copy resetting IO params and patch types"
             << endl << this->info() << endl;
     }
 
@@ -911,6 +911,20 @@ bool Foam::GeometricField<Type, PatchField, GeoMesh>::needReference() const
 template<class Type, template<class> class PatchField, class GeoMesh>
 void Foam::GeometricField<Type, PatchField, GeoMesh>::relax(const scalar alpha)
 {
+    if(alpha <= 0)
+    {
+        return;
+    }
+
+    //if (debug)
+    {
+        InfoIn
+        (
+            "GeometricField<Type, PatchField, GeoMesh>::relax"
+            "(const scalar alpha)"
+        )  << "Relaxing" << endl << this->info() << " by " << alpha << endl;
+    }
+
     operator==(prevIter() + alpha*(*this - prevIter()));
 }
 
@@ -918,16 +932,11 @@ void Foam::GeometricField<Type, PatchField, GeoMesh>::relax(const scalar alpha)
 template<class Type, template<class> class PatchField, class GeoMesh>
 void Foam::GeometricField<Type, PatchField, GeoMesh>::relax()
 {
-    scalar alpha = 0;
+    word name = this->name();
 
-    if (this->mesh().solutionDict().relax(this->name()))
+    if (this->mesh().solutionDict().relaxField(name))
     {
-        alpha = this->mesh().solutionDict().relaxationFactor(this->name());
-    }
-
-    if (alpha > 0)
-    {
-        relax(alpha);
+        relax(this->mesh().solutionDict().fieldRelaxationFactor(name));
     }
 }
 
@@ -949,6 +958,19 @@ Foam::word Foam::GeometricField<Type, PatchField, GeoMesh>::select
 }
 
 
+template<class Type, template<class> class PatchField, class GeoMesh>
+void Foam::GeometricField<Type, PatchField, GeoMesh>::writeMinMax
+(
+    Ostream& os
+) const
+{
+    os  << "min/max(" << this->name() << ") = "
+        << Foam::min(*this).value() << ", "
+        << Foam::max(*this).value()
+        << endl;
+}
+
+
 // writeData member function required by regIOobject
 template<class Type, template<class> class PatchField, class GeoMesh>
 bool Foam::GeometricField<Type, PatchField, GeoMesh>::
diff --git a/src/foam/fields/GeometricFields/GeometricField/GeometricField.H b/src/foam/fields/GeometricFields/GeometricField/GeometricField.H
index 6bfd6e31a..36a34e8f8 100644
--- a/src/foam/fields/GeometricFields/GeometricField/GeometricField.H
+++ b/src/foam/fields/GeometricFields/GeometricField/GeometricField.H
@@ -491,14 +491,17 @@ public:
         //  alpha = 0 : do nothing
         void relax(const scalar alpha);
 
+        //- Relax field (for steady-state solution).
+        //  alpha is read from controlDict
+        void relax();
+
         //- Select the final iteration parameters if 'final' is true
         //  by returning the field name + "Final"
         //  otherwise the standard parameters by returning the field name
         word select(bool final) const;
 
-        //- Relax field (for steady-state solution).
-        //  alpha is read from controlDict
-        void relax();
+        //- Helper function to write the min and max to an Ostream
+        void writeMinMax(Ostream& os) const;
 
 
     // Member function *this operators
diff --git a/src/foam/matrices/solution/solution.C b/src/foam/matrices/solution/solution.C
index 758d3507c..5a91b6d2c 100644
--- a/src/foam/matrices/solution/solution.C
+++ b/src/foam/matrices/solution/solution.C
@@ -49,6 +49,82 @@ static const Foam::List<Foam::word> subDictNames
 //! @endcond localScope
 
 
+// * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * * //
+
+void Foam::solution::read(const dictionary& dict)
+{
+    if (dict.found("cache"))
+    {
+        cache_ = dict.subDict("cache");
+        caching_ = cache_.lookupOrDefault("active", true);
+    }
+
+    if (dict.found("relaxationFactors"))
+    {
+        const dictionary& relaxDict(dict.subDict("relaxationFactors"));
+        if (relaxDict.found("fields") || relaxDict.found("equations"))
+        {
+            if (relaxDict.found("fields"))
+            {
+                fieldRelaxDict_ = relaxDict.subDict("fields");
+            }
+
+            if (relaxDict.found("equations"))
+            {
+                eqnRelaxDict_ = relaxDict.subDict("equations");
+            }
+        }
+        else
+        {
+            // backwards compatibility
+            fieldRelaxDict_.clear();
+
+            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;
+        }
+    }
+
+
+    if (dict.found("solvers"))
+    {
+        solvers_ = dict.subDict("solvers");
+        upgradeSolverDict(solvers_);
+    }
+}
+
+
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 Foam::solution::solution(const objectRegistry& obr, const fileName& dictName)
@@ -66,12 +142,11 @@ Foam::solution::solution(const objectRegistry& obr, const fileName& dictName)
     ),
     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 (!headerOk())
     {
@@ -86,7 +161,7 @@ Foam::solution::solution(const objectRegistry& obr, const fileName& dictName)
         }
     }
 
-    read();
+    read(solutionDict());
 }
 
 
@@ -170,43 +245,6 @@ Foam::label Foam::solution::upgradeSolverDict
 }
 
 
-bool Foam::solution::read()
-{
-    bool readOk = false;
-
-    if (headerOk())
-    {
-        readOk = regIOobject::read();
-    }
-
-    if (readOk)
-    {
-        if (found("cache"))
-        {
-            cache_ = subDict("cache");
-            caching_ = cache_.lookupOrDefault("active", true);
-        }
-
-        if (found("relaxationFactors"))
-        {
-            relaxationFactors_ = subDict("relaxationFactors");
-        }
-
-        relaxationFactors_.readIfPresent("default", defaultRelaxationFactor_);
-
-        if (found("solvers"))
-        {
-            solvers_ = subDict("solvers");
-            upgradeSolverDict(solvers_);
-        }
-
-        return true;
-    }
-
-    return readOk;
-}
-
-
 bool Foam::solution::cache(const word& name) const
 {
     if (caching_)
@@ -226,41 +264,52 @@ 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<< "Field relaxation factor for " << name
+            << " is " << (fieldRelaxDict_.found(name) ? "set" : "unset")
+            << 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);
 
@@ -269,6 +318,49 @@ Foam::scalar Foam::solution::relaxationFactor(const word& name) const
 }
 
 
+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);
+
+        return 0;
+    }
+}
+
+
+const Foam::dictionary& Foam::solution::solutionDict() const
+{
+    if (found("select"))
+    {
+        return subDict(word(lookup("select")));
+    }
+    else
+    {
+        return *this;
+    }
+}
+
+
 const Foam::dictionary& Foam::solution::solverDict(const word& name) const
 {
     if (debug)
@@ -293,15 +385,22 @@ const Foam::dictionary& Foam::solution::solver(const word& name) const
 }
 
 
+bool Foam::solution::read()
+{
+    if (regIOobject::read())
+    {
+        read(solutionDict());
+
+        return true;
+    }
+    else
+    {
+        return false;
+    }
+}
+
 bool Foam::solution::writeData(Ostream& os) const
 {
-    // Write dictionaries
-    os << nl << "solvers";
-    solvers_.write(os, true);
-
-    os << nl << "relaxationFactors";
-    relaxationFactors_.write(os, true);
-
     // Write direct entries of the solution dictionary
     // HJ, 16/Feb/2010
     dictionary::write(os, false);
diff --git a/src/foam/matrices/solution/solution.H b/src/foam/matrices/solution/solution.H
index 00e8e1fb3..310e288e0 100644
--- a/src/foam/matrices/solution/solution.H
+++ b/src/foam/matrices/solution/solution.H
@@ -60,10 +60,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_;
@@ -117,10 +123,20 @@ public:
             );
 
             //- 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
+            const dictionary& solutionDict() const;
 
             //- Return the solver controls dictionary for the given field
             const dictionary& solverDict(const word& name) const;
@@ -131,10 +147,16 @@ public:
 
         // Edit
 
-            //- Return access to relaxation factors dictionary
-            dictionary& relaxationFactors()
+            //- Return access to field relaxation factors dictionary
+            dictionary& fieldRelaxationFactors()
             {
-                return relaxationFactors_;
+                return fieldRelaxDict_;
+            }
+
+            //- Return access to equation relaxation factors dictionary
+            dictionary& equationRelaxationFactors()
+            {
+                return eqnRelaxDict_;
             }
 
             //- Return access to solvers dictionary
diff --git a/src/foam/matrices/tolerances/tolerances.C b/src/foam/matrices/tolerances/tolerances.C
index 214659f4a..e5518f8a8 100644
--- a/src/foam/matrices/tolerances/tolerances.C
+++ b/src/foam/matrices/tolerances/tolerances.C
@@ -62,7 +62,7 @@ bool tolerances::read()
 {
     if (regIOobject::read())
     {
-        word toleranceSetName(lookup("toleranceSet"));
+        const word toleranceSetName(lookup("toleranceSet"));
         const dictionary& toleranceSet(subDict(toleranceSetName));
 
         if (toleranceSet.found("relaxationFactors"))
diff --git a/src/tetFiniteElement/tetFemMatrix/tetFemMatrix.C b/src/tetFiniteElement/tetFemMatrix/tetFemMatrix.C
index ddac7e9c1..76429f44c 100644
--- a/src/tetFiniteElement/tetFemMatrix/tetFemMatrix.C
+++ b/src/tetFiniteElement/tetFemMatrix/tetFemMatrix.C
@@ -215,16 +215,18 @@ void tetFemMatrix<Type>::relax(const scalar alpha)
 template<class Type>
 void tetFemMatrix<Type>::relax()
 {
-    scalar alpha = 0;
-
     if (psi_.mesh().solutionDict().relax(psi_.name()))
     {
-        alpha = psi_.mesh().solutionDict().relaxationFactor(psi_.name());
+        relax(psi_.mesh().solutionDict().equationRelaxationFactor(psi_.name()));
     }
-
-    if (alpha > 0)
+    else
     {
-        relax(alpha);
+        if (debug)
+        {
+            InfoIn("void tetFemMatrix<Type>::relax()")
+                << "Relaxation factor for field " << psi_.name()
+                << " not found.  Relaxation will not be used." << endl;
+        }
     }
 }
 
diff --git a/tutorials/incompressible/simpleFoam/pitzDaily/system/fvSolution b/tutorials/incompressible/simpleFoam/pitzDaily/system/fvSolution
index 94616d261..4a7b2732e 100644
--- a/tutorials/incompressible/simpleFoam/pitzDaily/system/fvSolution
+++ b/tutorials/incompressible/simpleFoam/pitzDaily/system/fvSolution
@@ -74,12 +74,19 @@ SIMPLE
 
 relaxationFactors
 {
-    p               0.3;
-    U               0.7;
-    k               0.7;
-    epsilon         0.7;
-    R               0.7;
-    nuTilda         0.7;
+    fields
+    {
+        p               0.3;
+    }
+
+    equations
+    {
+        U               0.7;
+        k               0.7;
+        epsilon         0.7;
+        R               0.7;
+        nuTilda         0.7;
+    }
 }
 
 // ************************************************************************* //