diff --git a/src/OpenFOAM/matrices/LduMatrix/LduMatrix/LduMatrix.H b/src/OpenFOAM/matrices/LduMatrix/LduMatrix/LduMatrix.H
index 6ed2091a77..db730d7a97 100644
--- a/src/OpenFOAM/matrices/LduMatrix/LduMatrix/LduMatrix.H
+++ b/src/OpenFOAM/matrices/LduMatrix/LduMatrix/LduMatrix.H
@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2017 OpenFOAM Foundation
-    Copyright (C) 2019-2021 OpenCFD Ltd.
+    Copyright (C) 2019-2022 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -54,6 +54,7 @@ SourceFiles
 #define Foam_LduMatrix_H
 
 #include "lduMesh.H"
+#include "lduMatrix.H"
 #include "Field.H"
 #include "FieldField.H"
 #include "LduInterfaceFieldPtrsList.H"
@@ -69,8 +70,7 @@ namespace Foam
 
 // Forward Declarations
 
-template<class Type, class DType, class LUType>
-class LduMatrix;
+template<class Type, class DType, class LUType> class LduMatrix;
 
 template<class Type, class DType, class LUType>
 Ostream& operator<<
@@ -119,16 +119,13 @@ public:
 
         // Protected Data
 
-            //- Default maximum number of iterations in the solver
-            static const label defaultMaxIter_ = 1000;
-
             word fieldName_;
             const LduMatrix<Type, DType, LUType>& matrix_;
 
-            //- Dictionary of controls
+            //- Dictionary of solution controls
             dictionary controlDict_;
 
-            //- Level of verbosity in the solver output statements
+            //- Verbosity level for solver output statements
             int log_;
 
             //- Minimum number of iterations in the solver
@@ -137,6 +134,9 @@ public:
             //- Maximum number of iterations in the solver
             label maxIter_;
 
+            //- The matrix normalisation type
+            lduMatrix::normTypes normType_;
+
             //- Final convergence tolerance
             Type tolerance_;
 
@@ -146,7 +146,7 @@ public:
 
         // Protected Member Functions
 
-            //- Read the control parameters from the controlDict_
+            //- Read the control parameters from controlDict_
             virtual void readControls();
 
 
@@ -206,6 +206,7 @@ public:
 
         // Constructors
 
+            //- Construct for given field name, matrix and controls
             solver
             (
                 const word& fieldName,
@@ -225,9 +226,8 @@ public:
             );
 
 
-        // Destructor
-
-            virtual ~solver() = default;
+        //- Destructor
+        virtual ~solver() = default;
 
 
         // Member Functions
@@ -244,21 +244,33 @@ public:
 
 
             //- Read and reset the solver parameters from the given dictionary
-            virtual void read(const dictionary& solverDict);
+            virtual void read(const dictionary&);
 
             virtual SolverPerformance<Type> solve
             (
                 Field<Type>& psi
             ) const = 0;
 
+            //- Return the matrix norm using the specified norm method
+            Type normFactor
+            (
+                const Field<Type>& psi,
+                const Field<Type>& Apsi,
+                Field<Type>& tmpField,
+                const lduMatrix::normTypes normType
+            ) const;
+
             //- Return the matrix norm used to normalise the residual for the
-            //  stopping criterion
+            //- stopping criterion
             Type normFactor
             (
                 const Field<Type>& psi,
                 const Field<Type>& Apsi,
                 Field<Type>& tmpField
-            ) const;
+            ) const
+            {
+                return this->normFactor(psi, Apsi, tmpField, normType_);
+            }
     };
 
 
@@ -314,6 +326,7 @@ public:
 
         // Constructors
 
+            //- Construct for given field name and matrix
             smoother
             (
                 const word& fieldName,
@@ -332,9 +345,8 @@ public:
             );
 
 
-        // Destructor
-
-            virtual ~smoother() = default;
+        //- Destructor
+        virtual ~smoother() = default;
 
 
         // Member Functions
@@ -405,10 +417,8 @@ public:
 
         // Constructors
 
-            preconditioner
-            (
-                const solver& sol
-            )
+            //- Construct for given solver
+            preconditioner(const solver& sol)
             :
                 solver_(sol)
             {}
@@ -424,16 +434,15 @@ public:
             );
 
 
-        // Destructor
-
-            virtual ~preconditioner() = default;
+        //- Destructor
+        virtual ~preconditioner() = default;
 
 
         // Member functions
 
             //- Read and reset the preconditioner parameters
-            //  from the given dictionary
-            virtual void read(const dictionary& preconditionerDict)
+            //- from the given dictionary
+            virtual void read(const dictionary&)
             {}
 
             //- Return wA the preconditioned form of residual rA
@@ -444,7 +453,7 @@ public:
             ) const = 0;
 
             //- Return wT the transpose-matrix preconditioned form of
-            //  residual rT.
+            //- residual rT.
             //  This is only required for preconditioning asymmetric matrices.
             virtual void preconditionT
             (
@@ -480,17 +489,16 @@ public:
         LduMatrix(const lduMesh&, Istream&);
 
 
-    // Destructor
-
-        ~LduMatrix();
+    //- Destructor
+    ~LduMatrix();
 
 
-    // Member functions
+    // Member Functions
 
         // Access to addressing
 
             //- Return the LDU mesh from which the addressing is obtained
-            const lduMesh& mesh() const
+            const lduMesh& mesh() const noexcept
             {
                 return lduMesh_;
             }
@@ -554,43 +562,43 @@ public:
             }
 
 
-            bool hasDiag() const
+            bool hasDiag() const noexcept
             {
                 return (diagPtr_);
             }
 
-            bool hasUpper() const
+            bool hasUpper() const noexcept
             {
                 return (upperPtr_);
             }
 
-            bool hasLower() const
+            bool hasLower() const noexcept
             {
                 return (lowerPtr_);
             }
 
-            bool hasSource() const
+            bool hasSource() const noexcept
             {
                 return (sourcePtr_);
             }
 
-            bool diagonal() const
+            bool diagonal() const noexcept
             {
                 return (diagPtr_ && !lowerPtr_ && !upperPtr_);
             }
 
-            bool symmetric() const
+            bool symmetric() const noexcept
             {
                 return (diagPtr_ && (!lowerPtr_ && upperPtr_));
             }
 
-            bool asymmetric() const
+            bool asymmetric() const noexcept
             {
                 return (diagPtr_ && lowerPtr_ && upperPtr_);
             }
 
 
-        // operations
+        // Operations
 
             void sumDiag();
             void negSumDiag();
@@ -640,7 +648,7 @@ public:
             tmp<Field<Type>> faceH(const tmp<Field<Type>>&) const;
 
 
-    // Member operators
+    // Member Operators
 
         void operator=(const LduMatrix<Type, DType, LUType>&);
 
@@ -653,7 +661,7 @@ public:
         void operator*=(scalar);
 
 
-    // Ostream operator
+    // Ostream Operator
 
         friend Ostream& operator<< <Type, DType, LUType>
         (
diff --git a/src/OpenFOAM/matrices/LduMatrix/LduMatrix/LduMatrixSolver.C b/src/OpenFOAM/matrices/LduMatrix/LduMatrix/LduMatrixSolver.C
index 59d117381e..a9598b03a6 100644
--- a/src/OpenFOAM/matrices/LduMatrix/LduMatrix/LduMatrixSolver.C
+++ b/src/OpenFOAM/matrices/LduMatrix/LduMatrix/LduMatrixSolver.C
@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2017 OpenFOAM Foundation
-    Copyright (C) 2019-2021 OpenCFD Ltd.
+    Copyright (C) 2019-2022 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -131,8 +131,9 @@ Foam::LduMatrix<Type, DType, LUType>::solver::solver
 
     log_(1),
     minIter_(0),
-    maxIter_(defaultMaxIter_),
-    tolerance_(1e-6*pTraits<Type>::one),
+    maxIter_(lduMatrix::defaultMaxIter),
+    normType_(lduMatrix::normTypes::DEFAULT_NORM),
+    tolerance_(lduMatrix::defaultTolerance*pTraits<Type>::one),
     relTol_(Zero)
 {
     readControls();
@@ -145,6 +146,8 @@ template<class Type, class DType, class LUType>
 void Foam::LduMatrix<Type, DType, LUType>::solver::readControls()
 {
     controlDict_.readIfPresent("log", log_);
+    normType_ = lduMatrix::normTypes::DEFAULT_NORM;
+    lduMatrix::normTypesNames_.readIfPresent("norm", controlDict_, normType_);
     controlDict_.readIfPresent("minIter", minIter_);
     controlDict_.readIfPresent("maxIter", maxIter_);
     controlDict_.readIfPresent("tolerance", tolerance_);
@@ -168,21 +171,45 @@ Type Foam::LduMatrix<Type, DType, LUType>::solver::normFactor
 (
     const Field<Type>& psi,
     const Field<Type>& Apsi,
-    Field<Type>& tmpField
+    Field<Type>& tmpField,
+    const lduMatrix::normTypes normType
 ) const
 {
-    // --- Calculate A dot reference value of psi
-    matrix_.sumA(tmpField);
-    cmptMultiply(tmpField, tmpField, gAverage(psi));
+    switch (normType)
+    {
+        case lduMatrix::normTypes::NO_NORM :
+        {
+            break;
+        }
 
-    return stabilise
-    (
-        gSum(cmptMag(Apsi - tmpField) + cmptMag(matrix_.source() - tmpField)),
-        SolverPerformance<Type>::small_
-    );
+        case lduMatrix::normTypes::DEFAULT_NORM :
+        case lduMatrix::normTypes::L1_SCALED_NORM :
+        {
+            // --- Calculate A dot reference value of psi
+            matrix_.sumA(tmpField);
+            cmptMultiply(tmpField, tmpField, gAverage(psi));
 
-    // At convergence this simpler method is equivalent to the above
-    // return stabilise(2*gSumCmptMag(matrix_.source()), matrix_.small_);
+            return stabilise
+            (
+                gSum
+                (
+                    cmptMag(Apsi - tmpField)
+                  + cmptMag(matrix_.source() - tmpField)
+                ),
+                SolverPerformance<Type>::small_
+            );
+
+            // Equivalent at convergence:
+            // return stabilise
+            // (
+            //     2*gSumCmptMag(matrix_.source()), matrix_.small_
+            // );
+            break;
+        }
+    }
+
+    // Fall-through: no norm
+    return pTraits<Type>::one;
 }
 
 
diff --git a/src/OpenFOAM/matrices/LduMatrix/Solvers/DiagonalSolver/DiagonalSolver.H b/src/OpenFOAM/matrices/LduMatrix/Solvers/DiagonalSolver/DiagonalSolver.H
index 47055b7808..8918142637 100644
--- a/src/OpenFOAM/matrices/LduMatrix/Solvers/DiagonalSolver/DiagonalSolver.H
+++ b/src/OpenFOAM/matrices/LduMatrix/Solvers/DiagonalSolver/DiagonalSolver.H
@@ -34,8 +34,8 @@ SourceFiles
 
 \*---------------------------------------------------------------------------*/
 
-#ifndef DiagonalSolver_H
-#define DiagonalSolver_H
+#ifndef Foam_DiagonalSolver_H
+#define Foam_DiagonalSolver_H
 
 #include "LduMatrix.H"
 
@@ -53,7 +53,9 @@ class DiagonalSolver
 :
     public LduMatrix<Type, DType, LUType>::solver
 {
-    // Private Member Functions
+public:
+
+    // Generated Methods
 
         //- No copy construct
         DiagonalSolver(const DiagonalSolver&) = delete;
@@ -62,8 +64,6 @@ class DiagonalSolver
         void operator=(const DiagonalSolver&) = delete;
 
 
-public:
-
     //- Runtime type information
     TypeName("diagonal");
 
diff --git a/src/OpenFOAM/matrices/LduMatrix/Solvers/SmoothSolver/SmoothSolver.H b/src/OpenFOAM/matrices/LduMatrix/Solvers/SmoothSolver/SmoothSolver.H
index 278d29e776..9909a5db6d 100644
--- a/src/OpenFOAM/matrices/LduMatrix/Solvers/SmoothSolver/SmoothSolver.H
+++ b/src/OpenFOAM/matrices/LduMatrix/Solvers/SmoothSolver/SmoothSolver.H
@@ -37,8 +37,8 @@ SourceFiles
 
 \*---------------------------------------------------------------------------*/
 
-#ifndef SmoothSolver_H
-#define SmoothSolver_H
+#ifndef Foam_SmoothSolver_H
+#define Foam_SmoothSolver_H
 
 #include "lduMatrix.H"
 
@@ -56,10 +56,9 @@ class SmoothSolver
 :
     public LduMatrix<Type, DType, LUType>::solver
 {
-
 protected:
 
-    // Protected data
+    // Protected Data
 
         //- Number of sweeps before the evaluation of residual
         label nSweeps_;
diff --git a/src/OpenFOAM/matrices/lduMatrix/lduMatrix/lduMatrix.C b/src/OpenFOAM/matrices/lduMatrix/lduMatrix/lduMatrix.C
index 546d5db0c7..0ecf015fd4 100644
--- a/src/OpenFOAM/matrices/lduMatrix/lduMatrix/lduMatrix.C
+++ b/src/OpenFOAM/matrices/lduMatrix/lduMatrix/lduMatrix.C
@@ -41,7 +41,16 @@ namespace Foam
 }
 
 
-const Foam::label Foam::lduMatrix::solver::defaultMaxIter_ = 1000;
+const Foam::Enum
+<
+    Foam::lduMatrix::normTypes
+>
+Foam::lduMatrix::normTypesNames_
+({
+    { normTypes::NO_NORM, "none" },
+    { normTypes::DEFAULT_NORM, "default" },
+    { normTypes::L1_SCALED_NORM, "L1_scaled" },
+});
 
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
diff --git a/src/OpenFOAM/matrices/lduMatrix/lduMatrix/lduMatrix.H b/src/OpenFOAM/matrices/lduMatrix/lduMatrix/lduMatrix.H
index d19365b63b..6ebb028e56 100644
--- a/src/OpenFOAM/matrices/lduMatrix/lduMatrix/lduMatrix.H
+++ b/src/OpenFOAM/matrices/lduMatrix/lduMatrix/lduMatrix.H
@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2017 OpenFOAM Foundation
-    Copyright (C) 2016-2021 OpenCFD Ltd.
+    Copyright (C) 2016-2022 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -50,8 +50,8 @@ SourceFiles
 
 \*---------------------------------------------------------------------------*/
 
-#ifndef lduMatrix_H
-#define lduMatrix_H
+#ifndef Foam_lduMatrix_H
+#define Foam_lduMatrix_H
 
 #include "lduMesh.H"
 #include "primitiveFieldsFwd.H"
@@ -62,6 +62,7 @@ SourceFiles
 #include "runTimeSelectionTables.H"
 #include "solverPerformance.H"
 #include "InfoProxy.H"
+#include "Enum.H"
 #include "profilingTrigger.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -70,7 +71,6 @@ namespace Foam
 {
 
 // Forward Declarations
-
 class lduMatrix;
 
 Ostream& operator<<(Ostream&, const lduMatrix&);
@@ -95,6 +95,26 @@ class lduMatrix
 
 public:
 
+    // Public Types
+
+        //- Enumerated matrix normalisation types
+        enum class normTypes : char
+        {
+            NO_NORM,            //!< "none" norm (returns 1)
+            DEFAULT_NORM,       //!< "default" norm (== L1_scaled)
+            L1_SCALED_NORM,     //!< "L1_scaled" norm
+        };
+
+        //- Names for the normTypes
+        static const Enum<normTypes> normTypesNames_;
+
+        //- Default maximum number of iterations for solvers
+        static constexpr label defaultMaxIter = 1000;
+
+        //- Default (absolute) tolerance
+        static constexpr scalar defaultTolerance = 1e-6;
+
+
     //- Abstract base-class for lduMatrix solvers
     class solver
     {
@@ -102,19 +122,16 @@ public:
 
         // Protected Data
 
-            //- Default maximum number of iterations in the solver
-            static const label defaultMaxIter_;
-
             word fieldName_;
             const lduMatrix& matrix_;
             const FieldField<Field, scalar>& interfaceBouCoeffs_;
             const FieldField<Field, scalar>& interfaceIntCoeffs_;
             lduInterfaceFieldPtrsList interfaces_;
 
-            //- Dictionary of controls
+            //- Dictionary of solution controls
             dictionary controlDict_;
 
-            //- Level of verbosity in the solver output statements
+            //- Verbosity level for solver output statements
             int log_;
 
             //- Minimum number of iterations in the solver
@@ -123,18 +140,22 @@ public:
             //- Maximum number of iterations in the solver
             label maxIter_;
 
+            //- The normalisation type
+            lduMatrix::normTypes normType_;
+
             //- Final convergence tolerance
             scalar tolerance_;
 
             //- Convergence tolerance relative to the initial
             scalar relTol_;
 
+            //- Profiling instrumentation
             profilingTrigger profiling_;
 
 
         // Protected Member Functions
 
-            //- Read the control parameters from the controlDict_
+            //- Read the control parameters from controlDict_
             virtual void readControls();
 
 
@@ -195,6 +216,7 @@ public:
 
         // Constructors
 
+            //- Construct solver for given field name, matrix etc
             solver
             (
                 const word& fieldName,
@@ -272,6 +294,16 @@ public:
                 const direction cmpt=0
             ) const;
 
+            //- Return the matrix norm using the specified norm method
+            solveScalarField::cmptType normFactor
+            (
+                const solveScalarField& psi,
+                const solveScalarField& source,
+                const solveScalarField& Apsi,
+                solveScalarField& tmpField,
+                const lduMatrix::normTypes normType
+            ) const;
+
             //- Return the matrix norm used to normalise the residual for the
             //- stopping criterion
             solveScalarField::cmptType normFactor
@@ -280,7 +312,10 @@ public:
                 const solveScalarField& source,
                 const solveScalarField& Apsi,
                 solveScalarField& tmpField
-            ) const;
+            ) const
+            {
+                return this->normFactor(psi, source, Apsi, tmpField, normType_);
+            }
     };
 
 
@@ -354,6 +389,7 @@ public:
 
         // Constructors
 
+            //- Construct for given field name, matrix etc
             smoother
             (
                 const word& fieldName,
@@ -479,10 +515,8 @@ public:
 
         // Constructors
 
-            preconditioner
-            (
-                const solver& sol
-            )
+            //- Construct for given solver
+            explicit preconditioner(const solver& sol)
             :
                 solver_(sol)
             {}
@@ -564,7 +598,7 @@ public:
         // Access to addressing
 
             //- Return the LDU mesh from which the addressing is obtained
-            const lduMesh& mesh() const
+            const lduMesh& mesh() const noexcept
             {
                 return lduMesh_;
             }
@@ -606,38 +640,38 @@ public:
             const scalarField& diag() const;
             const scalarField& upper() const;
 
-            bool hasDiag() const
+            bool hasDiag() const noexcept
             {
                 return (diagPtr_);
             }
 
-            bool hasUpper() const
+            bool hasUpper() const noexcept
             {
                 return (upperPtr_);
             }
 
-            bool hasLower() const
+            bool hasLower() const noexcept
             {
                 return (lowerPtr_);
             }
 
-            bool diagonal() const
+            bool diagonal() const noexcept
             {
                 return (diagPtr_ && !lowerPtr_ && !upperPtr_);
             }
 
-            bool symmetric() const
+            bool symmetric() const noexcept
             {
                 return (diagPtr_ && (!lowerPtr_ && upperPtr_));
             }
 
-            bool asymmetric() const
+            bool asymmetric() const noexcept
             {
                 return (diagPtr_ && lowerPtr_ && upperPtr_);
             }
 
 
-        // operations
+        // Operations
 
             void sumDiag();
             void negSumDiag();
diff --git a/src/OpenFOAM/matrices/lduMatrix/lduMatrix/lduMatrixSolver.C b/src/OpenFOAM/matrices/lduMatrix/lduMatrix/lduMatrixSolver.C
index bae0f3d9e7..1ffcf3bdb0 100644
--- a/src/OpenFOAM/matrices/lduMatrix/lduMatrix/lduMatrixSolver.C
+++ b/src/OpenFOAM/matrices/lduMatrix/lduMatrix/lduMatrixSolver.C
@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2017 OpenFOAM Foundation
-    Copyright (C) 2016-2021 OpenCFD Ltd.
+    Copyright (C) 2016-2022 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -152,6 +152,14 @@ Foam::lduMatrix::solver::solver
     interfaceIntCoeffs_(interfaceIntCoeffs),
     interfaces_(interfaces),
     controlDict_(solverControls),
+
+    log_(1),
+    minIter_(0),
+    maxIter_(lduMatrix::defaultMaxIter),
+    normType_(lduMatrix::normTypes::DEFAULT_NORM),
+    tolerance_(lduMatrix::defaultTolerance),
+    relTol_(Zero),
+
     profiling_("lduMatrix::solver." + fieldName)
 {
     readControls();
@@ -162,11 +170,19 @@ Foam::lduMatrix::solver::solver
 
 void Foam::lduMatrix::solver::readControls()
 {
-    log_ = controlDict_.getOrDefault<int>("log", 1);
-    minIter_ = controlDict_.getOrDefault<label>("minIter", 0);
-    maxIter_ = controlDict_.getOrDefault<label>("maxIter", defaultMaxIter_);
-    tolerance_ = controlDict_.getOrDefault<scalar>("tolerance", 1e-6);
-    relTol_ = controlDict_.getOrDefault<scalar>("relTol", 0);
+    log_ = 1;
+    minIter_ = 0;
+    maxIter_ = lduMatrix::defaultMaxIter;
+    normType_ = lduMatrix::normTypes::DEFAULT_NORM;
+    tolerance_ = lduMatrix::defaultTolerance;
+    relTol_ = 0;
+
+    controlDict_.readIfPresent("log", log_);
+    lduMatrix::normTypesNames_.readIfPresent("norm", controlDict_, normType_);
+    controlDict_.readIfPresent("minIter", minIter_);
+    controlDict_.readIfPresent("maxIter", maxIter_);
+    controlDict_.readIfPresent("tolerance", tolerance_);
+    controlDict_.readIfPresent("relTol", relTol_);
 }
 
 
@@ -199,24 +215,40 @@ Foam::solveScalarField::cmptType Foam::lduMatrix::solver::normFactor
     const solveScalarField& psi,
     const solveScalarField& source,
     const solveScalarField& Apsi,
-    solveScalarField& tmpField
+    solveScalarField& tmpField,
+    const lduMatrix::normTypes normType
 ) const
 {
-    // --- Calculate A dot reference value of psi
-    matrix_.sumA(tmpField, interfaceBouCoeffs_, interfaces_);
+    switch (normType)
+    {
+        case lduMatrix::normTypes::NO_NORM :
+        {
+            break;
+        }
 
-    tmpField *= gAverage(psi, matrix_.mesh().comm());
+        case lduMatrix::normTypes::DEFAULT_NORM :
+        case lduMatrix::normTypes::L1_SCALED_NORM :
+        {
+            // --- Calculate A dot reference value of psi
+            matrix_.sumA(tmpField, interfaceBouCoeffs_, interfaces_);
 
-    return
-        gSum
-        (
-            (mag(Apsi - tmpField) + mag(source - tmpField))(),
-            matrix_.mesh().comm()
-        )
-      + solverPerformance::small_;
+            tmpField *= gAverage(psi, matrix_.mesh().comm());
 
-    // At convergence this simpler method is equivalent to the above
-    // return 2*gSumMag(source) + solverPerformance::small_;
+            return
+                gSum
+                (
+                    (mag(Apsi - tmpField) + mag(source - tmpField))(),
+                    matrix_.mesh().comm()
+                ) + solverPerformance::small_;
+
+            // Equivalent at convergence:
+            // return 2*gSumMag(source) + solverPerformance::small_;
+            break;
+        }
+    }
+
+    // Fall-through: no norm
+    return solveScalarField::cmptType(1);
 }
 
 
diff --git a/src/OpenFOAM/matrices/lduMatrix/preconditioners/GAMGPreconditioner/GAMGPreconditioner.C b/src/OpenFOAM/matrices/lduMatrix/preconditioners/GAMGPreconditioner/GAMGPreconditioner.C
index 4619815c81..daf6b62ad2 100644
--- a/src/OpenFOAM/matrices/lduMatrix/preconditioners/GAMGPreconditioner/GAMGPreconditioner.C
+++ b/src/OpenFOAM/matrices/lduMatrix/preconditioners/GAMGPreconditioner/GAMGPreconditioner.C
@@ -67,12 +67,6 @@ Foam::GAMGPreconditioner::GAMGPreconditioner
 }
 
 
-// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
-
-Foam::GAMGPreconditioner::~GAMGPreconditioner()
-{}
-
-
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
 void Foam::GAMGPreconditioner::readControls()
@@ -89,7 +83,7 @@ void Foam::GAMGPreconditioner::precondition
     const direction cmpt
 ) const
 {
-    wA = 0.0;
+    wA = Zero;
     solveScalarField AwA(wA.size());
     solveScalarField finestCorrection(wA.size());
     solveScalarField finestResidual(rA_ss);
diff --git a/src/OpenFOAM/matrices/lduMatrix/preconditioners/GAMGPreconditioner/GAMGPreconditioner.H b/src/OpenFOAM/matrices/lduMatrix/preconditioners/GAMGPreconditioner/GAMGPreconditioner.H
index b70af39402..8b3427c089 100644
--- a/src/OpenFOAM/matrices/lduMatrix/preconditioners/GAMGPreconditioner/GAMGPreconditioner.H
+++ b/src/OpenFOAM/matrices/lduMatrix/preconditioners/GAMGPreconditioner/GAMGPreconditioner.H
@@ -61,7 +61,8 @@ class GAMGPreconditioner
     public lduMatrix::preconditioner
 {
 protected:
-    // Protected data
+
+    // Protected Data
 
         //- Number of V-cycles to perform
         label nVcycles_;
@@ -86,7 +87,7 @@ public:
 
 
     //- Destructor
-    virtual ~GAMGPreconditioner();
+    virtual ~GAMGPreconditioner() = default;
 
 
     // Member Functions
diff --git a/src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/GAMGProcAgglomerations/procFacesGAMGProcAgglomeration/procFacesGAMGProcAgglomeration.H b/src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/GAMGProcAgglomerations/procFacesGAMGProcAgglomeration/procFacesGAMGProcAgglomeration.H
index 3ed7a5f2c8..7847a5c70d 100644
--- a/src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/GAMGProcAgglomerations/procFacesGAMGProcAgglomeration/procFacesGAMGProcAgglomeration.H
+++ b/src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/GAMGProcAgglomerations/procFacesGAMGProcAgglomeration/procFacesGAMGProcAgglomeration.H
@@ -51,6 +51,7 @@ SourceFiles
 namespace Foam
 {
 
+// Forward Declarations
 class GAMGAgglomeration;
 class lduMesh;
 class lduPrimitiveMesh;
@@ -63,7 +64,7 @@ class procFacesGAMGProcAgglomeration
 :
     public GAMGProcAgglomeration
 {
-    // Private data
+    // Private Data
 
         //- When to processor agglomerate
         const label nAgglomeratingCells_;
@@ -121,9 +122,8 @@ public:
 
     // Member Functions
 
-       //- Modify agglomeration. Return true if modified
+        //- Modify agglomeration. Return true if modified
         virtual bool agglomerate();
-
 };
 
 
diff --git a/src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/GAMGSolver.C b/src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/GAMGSolver.C
index bd79b8c07f..45fd5b9c0e 100644
--- a/src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/GAMGSolver.C
+++ b/src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/GAMGSolver.C
@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2017 OpenFOAM Foundation
-    Copyright (C) 2021 OpenCFD Ltd.
+    Copyright (C) 2021-2022 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -69,7 +69,6 @@ Foam::GAMGSolver::GAMGSolver
 
     // Default values for all controls
     // which may be overridden by those in controlDict
-    cacheAgglomeration_(true),
     nPreSweeps_(0),
     preSweepsLevelMultiplier_(1),
     maxPreSweeps_(4),
@@ -77,9 +76,12 @@ Foam::GAMGSolver::GAMGSolver
     postSweepsLevelMultiplier_(1),
     maxPostSweeps_(4),
     nFinestSweeps_(2),
+
+    cacheAgglomeration_(true),
     interpolateCorrection_(false),
     scaleCorrection_(matrix.symmetric()),
     directSolveCoarsest_(false),
+
     agglomeration_(GAMGAgglomeration::New(matrix_, controlDict_)),
 
     matrixLevels_(agglomeration_.size()),
@@ -389,14 +391,7 @@ void Foam::GAMGSolver::readControls()
 
 const Foam::lduMatrix& Foam::GAMGSolver::matrixLevel(const label i) const
 {
-    if (i == 0)
-    {
-        return matrix_;
-    }
-    else
-    {
-        return matrixLevels_[i - 1];
-    }
+    return i ? matrixLevels_[i-1] : matrix_;
 }
 
 
@@ -405,14 +400,7 @@ const Foam::lduInterfaceFieldPtrsList& Foam::GAMGSolver::interfaceLevel
     const label i
 ) const
 {
-    if (i == 0)
-    {
-        return interfaces_;
-    }
-    else
-    {
-        return interfaceLevels_[i - 1];
-    }
+    return i ? interfaceLevels_[i-1] : interfaces_;
 }
 
 
@@ -422,14 +410,7 @@ Foam::GAMGSolver::interfaceBouCoeffsLevel
     const label i
 ) const
 {
-    if (i == 0)
-    {
-        return interfaceBouCoeffs_;
-    }
-    else
-    {
-        return interfaceLevelsBouCoeffs_[i - 1];
-    }
+    return i ? interfaceLevelsBouCoeffs_[i-1] : interfaceBouCoeffs_;
 }
 
 
@@ -439,14 +420,7 @@ Foam::GAMGSolver::interfaceIntCoeffsLevel
     const label i
 ) const
 {
-    if (i == 0)
-    {
-        return interfaceIntCoeffs_;
-    }
-    else
-    {
-        return interfaceLevelsIntCoeffs_[i - 1];
-    }
+    return i ? interfaceLevelsIntCoeffs_[i-1] : interfaceIntCoeffs_;
 }
 
 
diff --git a/src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/GAMGSolver.H b/src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/GAMGSolver.H
index d3a9405e2a..ace02304d8 100644
--- a/src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/GAMGSolver.H
+++ b/src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/GAMGSolver.H
@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2016 OpenFOAM Foundation
-    Copyright (C) 2019 OpenCFD Ltd.
+    Copyright (C) 2019-2022 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -56,12 +56,11 @@ SourceFiles
 
 \*---------------------------------------------------------------------------*/
 
-#ifndef GAMGSolver_H
-#define GAMGSolver_H
+#ifndef Foam_GAMGSolver_H
+#define Foam_GAMGSolver_H
 
 #include "GAMGAgglomeration.H"
 #include "lduMatrix.H"
-#include "labelField.H"
 #include "primitiveFields.H"
 #include "LUscalarMatrix.H"
 
@@ -78,9 +77,7 @@ class GAMGSolver
 :
     public lduMatrix::solver
 {
-    // Private data
-
-        bool cacheAgglomeration_;
+    // Private Data
 
         //- Number of pre-smoothing sweeps
         label nPreSweeps_;
@@ -103,6 +100,9 @@ class GAMGSolver
         //- Number of smoothing sweeps on finest mesh
         label nFinestSweeps_;
 
+        //- Cache the agglomeration (default: true)
+        bool cacheAgglomeration_;
+
         //- Choose if the corrections should be interpolated after injection.
         //  By default corrections are not interpolated.
         bool interpolateCorrection_;
diff --git a/src/OpenFOAM/matrices/lduMatrix/solvers/PBiCG/PBiCG.C b/src/OpenFOAM/matrices/lduMatrix/solvers/PBiCG/PBiCG.C
index d033d292fc..e4c9b7a75f 100644
--- a/src/OpenFOAM/matrices/lduMatrix/solvers/PBiCG/PBiCG.C
+++ b/src/OpenFOAM/matrices/lduMatrix/solvers/PBiCG/PBiCG.C
@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2017 OpenFOAM Foundation
-    Copyright (C) 2019-2021 OpenCFD Ltd.
+    Copyright (C) 2019-2022 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -224,11 +224,13 @@ Foam::solverPerformance Foam::PBiCG::solve
     }
 
     // Recommend PBiCGStab if PBiCG fails to converge
-    if (solverPerf.nIterations() > max(defaultMaxIter_, maxIter_))
+    const label upperMaxIters = max(maxIter_, lduMatrix::defaultMaxIter);
+
+    if (solverPerf.nIterations() > upperMaxIters)
     {
         FatalErrorInFunction
-            << "PBiCG has failed to converge within the maximum number"
-               " of iterations " << max(defaultMaxIter_, maxIter_) << nl
+            << "PBiCG has failed to converge within the maximum iterations: "
+            << upperMaxIters << nl
             << "    Please try the more robust PBiCGStab solver."
             << exit(FatalError);
     }
diff --git a/src/OpenFOAM/matrices/lduMatrix/solvers/diagonalSolver/diagonalSolver.H b/src/OpenFOAM/matrices/lduMatrix/solvers/diagonalSolver/diagonalSolver.H
index 54f62c7d60..9249c763f2 100644
--- a/src/OpenFOAM/matrices/lduMatrix/solvers/diagonalSolver/diagonalSolver.H
+++ b/src/OpenFOAM/matrices/lduMatrix/solvers/diagonalSolver/diagonalSolver.H
@@ -37,8 +37,8 @@ SourceFiles
 
 \*---------------------------------------------------------------------------*/
 
-#ifndef diagonalSolver_H
-#define diagonalSolver_H
+#ifndef Foam_diagonalSolver_H
+#define Foam_diagonalSolver_H
 
 #include "lduMatrix.H"
 
@@ -55,7 +55,9 @@ class diagonalSolver
 :
     public lduMatrix::solver
 {
-    // Private Member Functions
+public:
+
+    // Generated Methods
 
         //- No copy construct
         diagonalSolver(const diagonalSolver&) = delete;
@@ -64,8 +66,6 @@ class diagonalSolver
         void operator=(const diagonalSolver&) = delete;
 
 
-public:
-
     //- Runtime type information
     TypeName("diagonal");
 
diff --git a/src/OpenFOAM/matrices/lduMatrix/solvers/smoothSolver/smoothSolver.H b/src/OpenFOAM/matrices/lduMatrix/solvers/smoothSolver/smoothSolver.H
index 4c9522539f..0c33854c2c 100644
--- a/src/OpenFOAM/matrices/lduMatrix/solvers/smoothSolver/smoothSolver.H
+++ b/src/OpenFOAM/matrices/lduMatrix/solvers/smoothSolver/smoothSolver.H
@@ -43,8 +43,8 @@ SourceFiles
 
 \*---------------------------------------------------------------------------*/
 
-#ifndef smoothSolver_H
-#define smoothSolver_H
+#ifndef Foam_smoothSolver_H
+#define Foam_smoothSolver_H
 
 #include "lduMatrix.H"
 
@@ -63,7 +63,7 @@ class smoothSolver
 {
 protected:
 
-    // Protected data
+    // Protected Data
 
         //- Number of sweeps before the evaluation of residual
         label nSweeps_;
@@ -95,6 +95,7 @@ public:
     //- Destructor
     virtual ~smoothSolver() = default;
 
+
     // Member Functions
 
         //- Solve the matrix with this solver
diff --git a/src/OpenFOAM/meshes/lduMesh/lduMesh.H b/src/OpenFOAM/meshes/lduMesh/lduMesh.H
index 84edbffb39..85e33ed605 100644
--- a/src/OpenFOAM/meshes/lduMesh/lduMesh.H
+++ b/src/OpenFOAM/meshes/lduMesh/lduMesh.H
@@ -33,8 +33,8 @@ Description
 
 \*---------------------------------------------------------------------------*/
 
-#ifndef lduMesh_H
-#define lduMesh_H
+#ifndef Foam_lduMesh_H
+#define Foam_lduMesh_H
 
 #include "lduAddressing.H"
 #include "lduInterfacePtrsList.H"
@@ -46,11 +46,8 @@ Description
 namespace Foam
 {
 
+// Forward Declarations
 class objectRegistry;
-
-
-// Forward declaration of friend functions and operators
-
 class lduMesh;
 
 Ostream& operator<<(Ostream&, const InfoProxy<lduMesh>&);
@@ -62,15 +59,12 @@ Ostream& operator<<(Ostream&, const InfoProxy<lduMesh>&);
 
 class lduMesh
 {
-
 public:
 
     //- Runtime type information
     TypeName("lduMesh");
 
 
-    // Constructors
-
     //- Destructor
     virtual ~lduMesh() = default;
 
@@ -114,7 +108,7 @@ public:
             }
 
 
-    // Ostream operator
+    // Ostream Operator
 
         friend Ostream& operator<<(Ostream&, const InfoProxy<lduMesh>&);
 };
diff --git a/src/OpenFOAM/meshes/lduMesh/lduPrimitiveMesh.C b/src/OpenFOAM/meshes/lduMesh/lduPrimitiveMesh.C
index 0c2e4c087b..656c58128f 100644
--- a/src/OpenFOAM/meshes/lduMesh/lduPrimitiveMesh.C
+++ b/src/OpenFOAM/meshes/lduMesh/lduPrimitiveMesh.C
@@ -37,25 +37,6 @@ License
 namespace Foam
 {
     defineTypeNameAndDebug(lduPrimitiveMesh, 0);
-
-    //- Less operator for pairs of \<processor\>\<index\>
-    class procLess
-    {
-        const labelPairList& list_;
-
-    public:
-
-        procLess(const labelPairList& list)
-        :
-            list_(list)
-        {}
-
-        bool operator()(const label a, const label b)
-        {
-            return list_[a].first() < list_[b].first();
-        }
-    };
-
 }
 
 
diff --git a/src/OpenFOAM/meshes/lduMesh/lduPrimitiveMesh.H b/src/OpenFOAM/meshes/lduMesh/lduPrimitiveMesh.H
index 5f1a8a6d68..3f6fcc78ee 100644
--- a/src/OpenFOAM/meshes/lduMesh/lduPrimitiveMesh.H
+++ b/src/OpenFOAM/meshes/lduMesh/lduPrimitiveMesh.H
@@ -35,8 +35,8 @@ SourceFiles
 
 \*---------------------------------------------------------------------------*/
 
-#ifndef lduPrimitiveMesh_H
-#define lduPrimitiveMesh_H
+#ifndef Foam_lduPrimitiveMesh_H
+#define Foam_lduPrimitiveMesh_H
 
 #include "lduMesh.H"
 #include "labelList.H"
@@ -55,7 +55,7 @@ class lduPrimitiveMesh
     public lduMesh,
     public lduAddressing
 {
-    // Private data
+    // Private Data
 
         //- Lower addressing
         labelList lowerAddr_;
diff --git a/tutorials/basic/scalarTransportFoam/Allclean b/tutorials/basic/scalarTransportFoam/Allclean
new file mode 100755
index 0000000000..ad79547de5
--- /dev/null
+++ b/tutorials/basic/scalarTransportFoam/Allclean
@@ -0,0 +1,19 @@
+#!/bin/sh
+cd "${0%/*}" || exit                                # Run from this directory
+. ${WM_PROJECT_DIR:?}/bin/tools/CleanFunctions      # Tutorial clean functions
+#------------------------------------------------------------------------------
+
+keepCases="pitzDaily"
+loseCases="pitzDaily-stepFunction"
+
+for caseName in $keepCases
+do
+    foamCleanTutorials -case="$caseName"
+done
+
+for caseName in $loseCases
+do
+    removeCase $caseName
+done
+
+#------------------------------------------------------------------------------
diff --git a/tutorials/basic/scalarTransportFoam/Allrun b/tutorials/basic/scalarTransportFoam/Allrun
new file mode 100755
index 0000000000..b392cb4e80
--- /dev/null
+++ b/tutorials/basic/scalarTransportFoam/Allrun
@@ -0,0 +1,36 @@
+#!/bin/sh
+cd "${0%/*}" || exit                                # Run from this directory
+. ${WM_PROJECT_DIR:?}/bin/tools/RunFunctions        # Tutorial run functions
+#------------------------------------------------------------------------------
+
+# Do pitzDaily
+( cd pitzDaily && foamRunTutorials )
+
+if true  ## if notTest "$@"
+then
+    # Clone case for additional tests
+    cloneCase pitzDaily pitzDaily-stepFunction
+
+    # Modify and execute
+    (
+        cd pitzDaily-stepFunction || exit
+
+        # Run a bit longer
+        foamDictionary system/controlDict -entry endTime -set 0.2
+
+        # Use table input to start scalar at 0.1s
+        ##runApplication changeDictionary -time 0
+        foamDictionary 0/T \
+            -entry boundaryField/inlet/uniformValue/type \
+            -set table
+
+        # Use 'none' for matrix norm
+        foamDictionary system/fvSolution \
+            -entry solvers/T/norm \
+            -set none
+
+        foamRunTutorials
+    )
+fi
+
+#------------------------------------------------------------------------------
diff --git a/tutorials/basic/scalarTransportFoam/pitzDaily/0/T b/tutorials/basic/scalarTransportFoam/pitzDaily/0/T
index 97d4d2f65d..160533c3a3 100644
--- a/tutorials/basic/scalarTransportFoam/pitzDaily/0/T
+++ b/tutorials/basic/scalarTransportFoam/pitzDaily/0/T
@@ -22,28 +22,41 @@ boundaryField
 {
     inlet
     {
-        type            fixedValue;
-        value           uniform 1;
+        type    uniformFixedValue;
+
+        uniformValue
+        {
+            type    constant;
+            value   1.0;
+
+            // Table entries (for modified version)
+            values
+            (
+                (0    1e-12)
+                (0.1  1e-12)
+                (0.10001 1.0)
+            );
+        }
     }
 
     outlet
     {
-        type            zeroGradient;
+        type    zeroGradient;
     }
 
     upperWall
     {
-        type            zeroGradient;
+        type    zeroGradient;
     }
 
     lowerWall
     {
-        type            zeroGradient;
+        type    zeroGradient;
     }
 
     frontAndBack
     {
-        type            empty;
+        type    empty;
     }
 }
 
diff --git a/tutorials/basic/scalarTransportFoam/pitzDaily/system/changeDictionaryDict b/tutorials/basic/scalarTransportFoam/pitzDaily/system/changeDictionaryDict
new file mode 100644
index 0000000000..d29010ec86
--- /dev/null
+++ b/tutorials/basic/scalarTransportFoam/pitzDaily/system/changeDictionaryDict
@@ -0,0 +1,41 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  v2206                                 |
+|   \\  /    A nd           | Website:  www.openfoam.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    object      changeDictionaryDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+T
+{
+    boundaryField
+    {
+        inlet
+        {
+            type    uniformFixedValue;
+
+            uniformValue
+            {
+                type    table;
+
+                values
+                (
+                    (0 1e-12)
+                    (0.1 1e-12)
+                    (0.10001 20)
+                );
+            }
+        }
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/basic/scalarTransportFoam/pitzDaily/system/fvSolution b/tutorials/basic/scalarTransportFoam/pitzDaily/system/fvSolution
index 5e9ff72ccc..0f308aa34a 100644
--- a/tutorials/basic/scalarTransportFoam/pitzDaily/system/fvSolution
+++ b/tutorials/basic/scalarTransportFoam/pitzDaily/system/fvSolution
@@ -1,7 +1,7 @@
 /*--------------------------------*- C++ -*----------------------------------*\
 | =========                 |                                                 |
 | \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
-|  \\    /   O peration     | Version:  v2206                                 |
+|  \\    /   O peration     | Version:  v2212                                 |
 |   \\  /    A nd           | Website:  www.openfoam.com                      |
 |    \\/     M anipulation  |                                                 |
 \*---------------------------------------------------------------------------*/
@@ -20,8 +20,9 @@ solvers
     {
         solver          PBiCGStab;
         preconditioner  DILU;
-        tolerance       1e-06;
+        tolerance       1e-6;
         relTol          0;
+        norm            default;
     }
 }