diff --git a/src/parallel/decompose/decompositionMethods/hierarchGeomDecomp/hierarchGeomDecomp.C b/src/parallel/decompose/decompositionMethods/hierarchGeomDecomp/hierarchGeomDecomp.C
index fd51978a4eeab9ac3b6e29e509c8ada50770582b..766adb25af8c779c3444e6d0dfefe2407746bc52 100644
--- a/src/parallel/decompose/decompositionMethods/hierarchGeomDecomp/hierarchGeomDecomp.C
+++ b/src/parallel/decompose/decompositionMethods/hierarchGeomDecomp/hierarchGeomDecomp.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2015-2018 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2015-2019 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
                             | Copyright (C) 2011-2016 OpenFOAM Foundation
@@ -54,36 +54,37 @@ namespace Foam
 
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
-void Foam::hierarchGeomDecomp::setDecompOrder()
+void Foam::hierarchGeomDecomp::setOrder()
 {
-    word order;
+    const word order(coeffsDict_.lookupOrDefault<word>("order", ""));
 
-    if (coeffsDict_.readIfPresent("order", order))
+    if (order.empty())
     {
-        if (order.size() != 3)
-        {
-            FatalIOErrorInFunction(decompDict_)
-                << "number of characters in order (" << order << ") != 3"
-                << exit(FatalIOError);
-        }
+        return;
+    }
+    else if (order.size() != 3)
+    {
+        FatalIOErrorInFunction(decompDict_)
+            << "Number of characters in order (" << order << ") != 3"
+            << exit(FatalIOError);
+    }
 
-        for (int i = 0; i < 3; ++i)
-        {
-            // Change [x-z] -> [0-2]
+    for (int i = 0; i < 3; ++i)
+    {
+        // Change [x-z] -> [0-2]
 
-            switch (order[i])
-            {
-                case 'x': decompOrder_[i] = 0; break;
-                case 'y': decompOrder_[i] = 1; break;
-                case 'z': decompOrder_[i] = 2; break;
-
-                default:
-                    FatalIOErrorInFunction(decompDict_)
-                        << "Illegal decomposition order " << order << nl
-                        << "It should only contain x, y or z"
-                        << exit(FatalError);
-                    break;
-            }
+        switch (order[i])
+        {
+            case 'x': order_[i] = 0; break;
+            case 'y': order_[i] = 1; break;
+            case 'z': order_[i] = 2; break;
+
+            default:
+                FatalIOErrorInFunction(decompDict_)
+                    << "Illegal decomposition order " << order << nl
+                    << "It should only contain x, y or z"
+                    << exit(FatalError);
+                break;
         }
     }
 }
@@ -91,25 +92,25 @@ void Foam::hierarchGeomDecomp::setDecompOrder()
 
 Foam::label Foam::hierarchGeomDecomp::findLower
 (
-    const List<scalar>& l,
-    const scalar t,
-    const label initLow,
-    const label initHigh
+    const UList<scalar>& list,
+    const scalar val,
+    const label first,
+    const label last
 )
 {
-    if (initHigh <= initLow)
+    label low = first;
+    label high = last;
+
+    if (high <= low)
     {
-        return initLow;
+        return low;
     }
 
-    label low = initLow;
-    label high = initHigh;
-
     while ((high - low) > 1)
     {
-        label mid = (low + high)/2;
+        const label mid = (low + high)/2;
 
-        if (l[mid] < t)
+        if (list[mid] < val)
         {
             low = mid;
         }
@@ -121,18 +122,14 @@ Foam::label Foam::hierarchGeomDecomp::findLower
 
     // high and low can still differ by one. Choose best.
 
-    label tIndex = -1;
-
-    if (l[high-1] < t)
+    if (list[high-1] < val)
     {
-        tIndex = high;
+        return high;
     }
     else
     {
-        tIndex = low;
+        return low;
     }
-
-    return tIndex;
 }
 
 
@@ -184,12 +181,11 @@ bool Foam::hierarchGeomDecomp::findBinary
     scalar& midValue            // value at mid
 )
 {
-    label low = minIndex;
     scalar lowValue = minValue;
-
     scalar highValue = maxValue;
-    // (one beyond) index of highValue
-    label high = values.size();
+
+    label low = minIndex;
+    label high = values.size();  // (one beyond) index of highValue
 
     // Safeguards to avoid infinite loop.
     scalar midValuePrev = VGREAT;
@@ -340,13 +336,13 @@ Foam::label Foam::hierarchGeomDecomp::sortComponent
     const label sizeTol,
     const pointField& points,
     const labelList& current,       // slice of points to decompose
-    const direction componentIndex, // index in decompOrder_
+    const direction componentIndex, // index in order_
     const label mult,               // multiplication factor for finalDecomp
     labelList& finalDecomp
 ) const
 {
     // Current component
-    const label compI = decompOrder_[componentIndex];
+    const label compI = order_[componentIndex];
 
     // Track the number of times that findBinary() did not converge
     label nWarnings = 0;
@@ -532,13 +528,13 @@ Foam::label Foam::hierarchGeomDecomp::sortComponent
     const scalarField& weights,
     const pointField& points,
     const labelList& current,       // slice of points to decompose
-    const direction componentIndex, // index in decompOrder_
+    const direction componentIndex, // index in order_
     const label mult,               // multiplication factor for finalDecomp
     labelList& finalDecomp
 ) const
 {
     // Current component
-    const label compI = decompOrder_[componentIndex];
+    const label compI = order_[componentIndex];
 
     // Track the number of times that findBinary() did not converge
     label nWarnings = 0;
@@ -727,9 +723,9 @@ Foam::hierarchGeomDecomp::hierarchGeomDecomp
 )
 :
     geomDecomp(typeName, decompDict),
-    decompOrder_({0,1,2})
+    order_({0,1,2})
 {
-    setDecompOrder();
+    setOrder();
 }
 
 
@@ -740,9 +736,9 @@ Foam::hierarchGeomDecomp::hierarchGeomDecomp
 )
 :
     geomDecomp(typeName, decompDict, regionName),
-    decompOrder_({0,1,2})
+    order_({0,1,2})
 {
-    setDecompOrder();
+    setOrder();
 }
 
 
@@ -773,7 +769,7 @@ Foam::labelList Foam::hierarchGeomDecomp::decompose
         sizeTol,
         rotatedPoints,
         slice,
-        0,              // Sort first component in decompOrder.
+        0,              // Sort first component in order_
         1,              // Offset for different x bins.
         finalDecomp
     );
@@ -816,7 +812,7 @@ Foam::labelList Foam::hierarchGeomDecomp::decompose
         weights,
         rotatedPoints,
         slice,
-        0,              // Sort first component in decompOrder.
+        0,              // Sort first component in order_
         1,              // Offset for different x bins.
         finalDecomp
     );
diff --git a/src/parallel/decompose/decompositionMethods/hierarchGeomDecomp/hierarchGeomDecomp.H b/src/parallel/decompose/decompositionMethods/hierarchGeomDecomp/hierarchGeomDecomp.H
index c042ca2e316833cb42a26855cd0de1af54cf85e3..b91a6ff7a73fb480d6fad4b740168495a8caa2df 100644
--- a/src/parallel/decompose/decompositionMethods/hierarchGeomDecomp/hierarchGeomDecomp.H
+++ b/src/parallel/decompose/decompositionMethods/hierarchGeomDecomp/hierarchGeomDecomp.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2017-2018 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2017-2019 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
                             | Copyright (C) 2011-2015 OpenFOAM Foundation
@@ -49,10 +49,10 @@ Description
 
     Method coefficients:
     \table
-        Property  | Description                 | Required | Default
-        n         | (nx ny nz)                  | yes
-        delta     | delta for rotation matrix   | no | 0.001
-        order     | order of operation          | no | xyz
+        Property  | Description                             | Required | Default
+        n         | (nx ny nz)                              | yes |
+        delta     | delta for rotation matrix               | no  | 0.001
+        order     | order of operation                      | no  | xyz
     \endtable
 
 SourceFiles
@@ -78,17 +78,27 @@ class hierarchGeomDecomp
 :
     public geomDecomp
 {
-    // Private data
+    // Private Data
 
         //- Decomposition order in terms of components.
-        FixedList<direction, 3> decompOrder_;
+        FixedList<direction, 3> order_;
 
 
     // Private Member Functions
 
         //- Convert ordering string ("xyz") into list of components.
         //  Checks for bad entries, but no check for duplicate entries.
-        void setDecompOrder();
+        void setOrder();
+
+        //- Find index of value in list between
+        //- first (inclusive) and last (exclusive)
+        static label findLower
+        (
+            const UList<scalar>& list,
+            const scalar val,
+            const label first,
+            const label last
+        );
 
         //- Evaluates the weighted sizes for each sorted point.
         static void calculateSortedWeightedSizes
@@ -101,15 +111,6 @@ class hierarchGeomDecomp
             scalarField& sortedWeightedSizes
         );
 
-        //- Find index of t in list inbetween indices left and right
-        static label findLower
-        (
-            const List<scalar>&,
-            const scalar t,
-            const label left,
-            const label right
-        );
-
         //- Find midValue (at local index mid) such that the number of
         //  elements between mid and leftIndex are (globally summed) the
         //  wantedSize. Binary search.