diff --git a/src/OpenFOAM/matrices/lduMatrix/lduMatrix/lduMatrix.H b/src/OpenFOAM/matrices/lduMatrix/lduMatrix/lduMatrix.H
index 881f88f..6c2db2c 100644
--- a/src/OpenFOAM/matrices/lduMatrix/lduMatrix/lduMatrix.H
+++ b/src/OpenFOAM/matrices/lduMatrix/lduMatrix/lduMatrix.H
@@ -96,6 +96,8 @@ public:
         label  noIterations_;
         bool   converged_;
         bool   singular_;
+        //- Normalization factor for residuals
+        scalar normFactor_;
 
 
     public:
@@ -109,7 +111,8 @@ public:
                 finalResidual_(0),
                 noIterations_(0),
                 converged_(false),
-                singular_(false)
+                singular_(false),
+                normFactor_(Foam::lduMatrix::small_)
             {}
 
             //- Construct from components
@@ -121,7 +124,8 @@ public:
                 const scalar fRes = 0,
                 const label  nIter = 0,
                 const bool   converged = false,
-                const bool   singular = false
+                const bool   singular = false,
+                const scalar normFactor = Foam::lduMatrix::small_
             )
             :
                 solverName_(solverName),
@@ -130,7 +134,8 @@ public:
                 finalResidual_(fRes),
                 noIterations_(nIter),
                 converged_(converged),
-                singular_(singular)
+                singular_(singular),
+                normFactor_(normFactor)
             {}
 
             //- Construct from Istream
@@ -210,6 +215,12 @@ public:
                 return singular_;
             }
 
+            //- Return final residual
+            scalar& normFactor()
+            {
+                return normFactor_;
+            }
+
             //- Convergence test
             bool checkConvergence
             (
@@ -413,9 +424,9 @@ public:
                 const scalarField& Apsi,
                 scalarField& tmpField
             ) const;
-    };
-
 
+    };
+    
     //- Abstract base-class for lduMatrix smoothers
     class smoother
     {
diff --git a/src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/GAMGSolverSolve.C b/src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/GAMGSolverSolve.C
index 3e347e0..d93440e 100644
--- a/src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/GAMGSolverSolve.C
+++ b/src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/GAMGSolverSolve.C
@@ -62,6 +62,7 @@ Foam::lduMatrix::solverPerformance Foam::GAMGSolver::solve
     // Calculate normalised residual for convergence test
     solverPerf.initialResidual() = gSumMag(finestResidual)/normFactor;
     solverPerf.finalResidual() = solverPerf.initialResidual();
+    solverPerf.normFactor() = normFactor;
 
 
     // Check convergence, solve if not converged
diff --git a/src/OpenFOAM/matrices/lduMatrix/solvers/PBiCG/PBiCG.C b/src/OpenFOAM/matrices/lduMatrix/solvers/PBiCG/PBiCG.C
index c433a8d..821c455 100644
--- a/src/OpenFOAM/matrices/lduMatrix/solvers/PBiCG/PBiCG.C
+++ b/src/OpenFOAM/matrices/lduMatrix/solvers/PBiCG/PBiCG.C
@@ -116,6 +116,7 @@ Foam::lduMatrix::solverPerformance Foam::PBiCG::solve
     // --- Calculate normalised residual norm
     solverPerf.initialResidual() = gSumMag(rA)/normFactor;
     solverPerf.finalResidual() = solverPerf.initialResidual();
+    solverPerf.normFactor() = normFactor;
 
     // --- Check convergence, solve if not converged
     if (!solverPerf.checkConvergence(tolerance_, relTol_))
diff --git a/src/OpenFOAM/matrices/lduMatrix/solvers/PCG/PCG.C b/src/OpenFOAM/matrices/lduMatrix/solvers/PCG/PCG.C
index 965c072..6b57aea 100644
--- a/src/OpenFOAM/matrices/lduMatrix/solvers/PCG/PCG.C
+++ b/src/OpenFOAM/matrices/lduMatrix/solvers/PCG/PCG.C
@@ -107,6 +107,7 @@ Foam::lduMatrix::solverPerformance Foam::PCG::solve
     // --- Calculate normalised residual norm
     solverPerf.initialResidual() = gSumMag(rA)/normFactor;
     solverPerf.finalResidual() = solverPerf.initialResidual();
+    solverPerf.normFactor() = normFactor;
 
     // --- Check convergence, solve if not converged
     if (!solverPerf.checkConvergence(tolerance_, relTol_))
diff --git a/src/OpenFOAM/matrices/lduMatrix/solvers/smoothSolver/smoothSolver.C b/src/OpenFOAM/matrices/lduMatrix/solvers/smoothSolver/smoothSolver.C
index a3d9a8f..1412d8e 100644
--- a/src/OpenFOAM/matrices/lduMatrix/solvers/smoothSolver/smoothSolver.C
+++ b/src/OpenFOAM/matrices/lduMatrix/solvers/smoothSolver/smoothSolver.C
@@ -124,6 +124,7 @@ Foam::lduMatrix::solverPerformance Foam::smoothSolver::solve
             // Calculate residual magnitude
             solverPerf.initialResidual() = gSumMag(source - Apsi)/normFactor;
             solverPerf.finalResidual() = solverPerf.initialResidual();
+            solverPerf.normFactor() = normFactor;
         }
 
         if (lduMatrix::debug >= 2)
diff --git a/src/finiteVolume/fvMatrices/fvMatrix/fvMatrix.C b/src/finiteVolume/fvMatrices/fvMatrix/fvMatrix.C
index d5086f5..017bdb7 100644
--- a/src/finiteVolume/fvMatrices/fvMatrix/fvMatrix.C
+++ b/src/finiteVolume/fvMatrices/fvMatrix/fvMatrix.C
@@ -1345,6 +1345,17 @@ Foam::lduMatrix::solverPerformance Foam::solve
 template<class Type>
 Foam::lduMatrix::solverPerformance Foam::solve
 (
+    fvMatrix<Type>& fvm,
+    const dictionary& solverControls,
+    Type& normFactor
+)
+{
+    return fvm.solve(solverControls, normFactor);
+}
+
+template<class Type>
+Foam::lduMatrix::solverPerformance Foam::solve
+(
     const tmp<fvMatrix<Type> >& tfvm,
     const dictionary& solverControls
 )
@@ -1357,6 +1368,21 @@ Foam::lduMatrix::solverPerformance Foam::solve
     return solverPerf;
 }
 
+template<class Type>
+Foam::lduMatrix::solverPerformance Foam::solve
+(
+    const tmp<fvMatrix<Type> >& tfvm,
+    const dictionary& solverControls,
+    Type& normFactor
+)
+{
+    lduMatrix::solverPerformance solverPerf =
+        const_cast<fvMatrix<Type>&>(tfvm()).solve(solverControls, normFactor);
+
+    tfvm.clear();
+
+    return solverPerf;
+}
 
 template<class Type>
 Foam::lduMatrix::solverPerformance Foam::solve(fvMatrix<Type>& fvm)
@@ -1365,6 +1391,12 @@ Foam::lduMatrix::solverPerformance Foam::solve(fvMatrix<Type>& fvm)
 }
 
 template<class Type>
+Foam::lduMatrix::solverPerformance Foam::solve(fvMatrix<Type>& fvm, Type& normFactor)
+{
+    return fvm.solve(normFactor);
+}
+
+template<class Type>
 Foam::lduMatrix::solverPerformance Foam::solve(const tmp<fvMatrix<Type> >& tfvm)
 {
     lduMatrix::solverPerformance solverPerf =
@@ -1375,6 +1407,16 @@ Foam::lduMatrix::solverPerformance Foam::solve(const tmp<fvMatrix<Type> >& tfvm)
     return solverPerf;
 }
 
+template<class Type>
+Foam::lduMatrix::solverPerformance Foam::solve(const tmp<fvMatrix<Type> >& tfvm, Type& normFactor)
+{
+    lduMatrix::solverPerformance solverPerf =
+        const_cast<fvMatrix<Type>&>(tfvm()).solve(normFactor);
+
+    tfvm.clear();
+
+    return solverPerf;
+}
 
 template<class Type>
 Foam::tmp<Foam::fvMatrix<Type> > Foam::correction
diff --git a/src/finiteVolume/fvMatrices/fvMatrix/fvMatrix.H b/src/finiteVolume/fvMatrices/fvMatrix/fvMatrix.H
index 4763714..0946f09 100644
--- a/src/finiteVolume/fvMatrices/fvMatrix/fvMatrix.H
+++ b/src/finiteVolume/fvMatrices/fvMatrix/fvMatrix.H
@@ -142,6 +142,9 @@ class fvMatrix
         mutable GeometricField<Type, fvsPatchField, surfaceMesh>
             *faceFluxCorrectionPtr_;
 
+        // Normalization factors for residuals
+        Type normFactor_;
+
 protected:
 
     //- Declare friendship with the fvSolver class
@@ -240,9 +243,13 @@ public:
             //  Use the given solver controls
             lduMatrix::solverPerformance solve(const dictionary&);
 
+            lduMatrix::solverPerformance solve(const dictionary&, Type& normFactor);
+
             //- Solve returning the solution statistics.
             //  Solver controls read from fvSolution
             lduMatrix::solverPerformance solve();
+
+            lduMatrix::solverPerformance solve(Type& normFactor);
     };
 
 
@@ -323,6 +330,11 @@ public:
                 return faceFluxCorrectionPtr_;
             }
 
+            Type normFactor() const
+            {
+                return normFactor_;
+            }
+
 
         // Operations
 
@@ -385,14 +397,20 @@ public:
             //  Solver controls read from fvSolution
             autoPtr<fvSolver> solver();
 
+            autoPtr<fvSolver> solver(Type& normFactor);
+
             //- Solve returning the solution statistics.
             //  Use the given solver controls
             lduMatrix::solverPerformance solve(const dictionary&);
 
+            lduMatrix::solverPerformance solve(const dictionary&, Type& normFactor);
+
             //- Solve returning the solution statistics.
             //  Solver controls read from fvSolution
             lduMatrix::solverPerformance solve();
 
+            lduMatrix::solverPerformance solve(Type& normFactor);
+
             //- Return the matrix residual
             tmp<Field<Type> > residual() const;
 
@@ -559,6 +577,8 @@ lduMatrix::solverPerformance solve
 template<class Type>
 lduMatrix::solverPerformance solve(fvMatrix<Type>&);
 
+template<class Type>
+lduMatrix::solverPerformance solve(fvMatrix<Type>&, Type& normFactor);
 
 //- Solve returning the solution statistics given convergence tolerance,
 //  deleting temporary matrix after solution.
@@ -572,7 +592,6 @@ lduMatrix::solverPerformance solve(const tmp<fvMatrix<Type> >&);
 template<class Type>
 tmp<fvMatrix<Type> > correction(const fvMatrix<Type>&);
 
-
 //- Return the correction form of the given temporary matrix
 //  by subtracting the matrix multiplied by the current field
 template<class Type>
diff --git a/src/finiteVolume/fvMatrices/fvMatrix/fvMatrixSolve.C b/src/finiteVolume/fvMatrices/fvMatrix/fvMatrixSolve.C
index 0b1f510..80df261 100644
--- a/src/finiteVolume/fvMatrices/fvMatrix/fvMatrixSolve.C
+++ b/src/finiteVolume/fvMatrices/fvMatrix/fvMatrixSolve.C
@@ -152,6 +152,9 @@ Foam::lduMatrix::solverPerformance Foam::fvMatrix<Type>::solve
         solverPerfVec = max(solverPerfVec, solverPerf);
         solverPerfVec.solverName() = solverPerf.solverName();
 
+        // Save normalization factor for residuals
+        normFactor_.component(cmpt) = solverPerf.normFactor();
+
         psi.internalField().replace(cmpt, psiCmpt);
         diag() = saveDiag;
     }
@@ -163,6 +166,18 @@ Foam::lduMatrix::solverPerformance Foam::fvMatrix<Type>::solve
     return solverPerfVec;
 }
 
+template<class Type>
+Foam::lduMatrix::solverPerformance Foam::fvMatrix<Type>::solve
+(
+    const dictionary& solverControls,
+    Type& normFactor
+)
+{
+    lduMatrix::solverPerformance solverPerfVec = solve(solverControls);
+    normFactor = normFactor_;
+    return solverPerfVec;
+}
+
 
 template<class Type>
 Foam::autoPtr<typename Foam::fvMatrix<Type>::fvSolver>
@@ -181,6 +196,23 @@ Foam::fvMatrix<Type>::solver()
     );
 }
 
+template<class Type>
+Foam::autoPtr<typename Foam::fvMatrix<Type>::fvSolver>
+Foam::fvMatrix<Type>::solver(Type& normFactor)
+{
+    return solver
+    (
+        psi_.mesh().solverDict
+        (
+            psi_.select
+            (
+                psi_.mesh().data::template lookupOrDefault<bool>
+                ("finalIteration", false)
+            )
+        ),
+        normFactor
+    );
+}
 
 template<class Type>
 Foam::lduMatrix::solverPerformance Foam::fvMatrix<Type>::fvSolver::solve()
@@ -198,6 +230,22 @@ Foam::lduMatrix::solverPerformance Foam::fvMatrix<Type>::fvSolver::solve()
     );
 }
 
+template<class Type>
+Foam::lduMatrix::solverPerformance Foam::fvMatrix<Type>::fvSolver::solve(Type& normFactor)
+{
+    return solve
+    (
+        fvMat_.psi_.mesh().solverDict
+        (
+            fvMat_.psi_.select
+            (
+                fvMat_.psi_.mesh().data::template lookupOrDefault<bool>
+                ("finalIteration", false)
+            )
+        ),
+        normFactor
+    );
+}
 
 template<class Type>
 Foam::lduMatrix::solverPerformance Foam::fvMatrix<Type>::solve()
@@ -215,6 +263,22 @@ Foam::lduMatrix::solverPerformance Foam::fvMatrix<Type>::solve()
     );
 }
 
+template<class Type>
+Foam::lduMatrix::solverPerformance Foam::fvMatrix<Type>::solve(Type& normFactor)
+{
+    return solve
+    (
+        psi_.mesh().solverDict
+        (
+            psi_.select
+            (
+                psi_.mesh().data::template lookupOrDefault<bool>
+                ("finalIteration", false)
+            )
+        ),
+        normFactor
+    );
+}
 
 template<class Type>
 Foam::tmp<Foam::Field<Type> > Foam::fvMatrix<Type>::residual() const
diff --git a/src/finiteVolume/fvMatrices/fvScalarMatrix/fvScalarMatrix.C b/src/finiteVolume/fvMatrices/fvScalarMatrix/fvScalarMatrix.C
index c17cf14..91c556b 100644
--- a/src/finiteVolume/fvMatrices/fvScalarMatrix/fvScalarMatrix.C
+++ b/src/finiteVolume/fvMatrices/fvScalarMatrix/fvScalarMatrix.C
@@ -170,6 +170,9 @@ Foam::lduMatrix::solverPerformance Foam::fvMatrix<Foam::scalar>::solve
 
     psi.mesh().setSolverPerformance(psi.name(), solverPerf);
 
+    // Calculate normalization factor for residuals
+    normFactor_ = solverPerf.normFactor();
+
     return solverPerf;
 }