diff --git a/applications/utilities/postProcessing/miscellaneous/postChannel/Make/files b/applications/utilities/postProcessing/miscellaneous/postChannel/Make/files
deleted file mode 100644
index 90ae1fa7d6..0000000000
--- a/applications/utilities/postProcessing/miscellaneous/postChannel/Make/files
+++ /dev/null
@@ -1,4 +0,0 @@
-postChannel.C
-channelIndex.C
-
-EXE = $(FOAM_APPBIN)/postChannel
diff --git a/applications/utilities/postProcessing/miscellaneous/postChannel/Make/options b/applications/utilities/postProcessing/miscellaneous/postChannel/Make/options
deleted file mode 100644
index c5e4b6238b..0000000000
--- a/applications/utilities/postProcessing/miscellaneous/postChannel/Make/options
+++ /dev/null
@@ -1,10 +0,0 @@
-EXE_INC = \
-    -I$(LIB_SRC)/meshTools/lnInclude \
-    -I$(LIB_SRC)/finiteVolume/lnInclude \
-    -I$(LIB_SRC)/sampling/lnInclude
-
-EXE_LIBS = \
-    -lmeshTools \
-    -lfiniteVolume \
-    -lgenericPatchFields \
-    -lsampling
diff --git a/applications/utilities/postProcessing/miscellaneous/postChannel/calculateFields.H b/applications/utilities/postProcessing/miscellaneous/postChannel/calculateFields.H
deleted file mode 100644
index 026861d875..0000000000
--- a/applications/utilities/postProcessing/miscellaneous/postChannel/calculateFields.H
+++ /dev/null
@@ -1,73 +0,0 @@
-        /*
-        volTensorField gradU(fvc::grad(U));
-        volSymmTensorField D(symm(fvc::grad(U)));
-        volTensorField Dprim(symm(fvc::grad(U - UMean)));
-
-        volScalarField prod(-((U - UMean)*(U - UMean)) && D);
-        */
-
-        /*
-        volScalarField txx
-        (
-            IOobject
-            (
-                "txx",
-                runTime.timeName(),
-                mesh,
-                IOobject::NO_READ,
-                IOobject::NO_WRITE
-            ),
-            mesh,
-            dimensionSet(0, 1, -1, 0, 0)
-        );
-        txx =sqrt(Txx - (UMeanx*UMeanx));
-        txx.write();
-
-        volScalarField tyy
-        (
-            IOobject
-            (
-                "tyy",
-                runTime.timeName(),
-                mesh,
-                IOobject::NO_READ,
-                IOobject::NO_WRITE
-            ),
-            mesh,
-            dimensionSet(0, 1, -1, 0, 0)
-        );
-        tyy = sqrt(Tyy - (UMeany*UMeany));
-        tyy.write();
-
-        volScalarField tzz
-        (
-            IOobject
-            (
-                "tzz",
-                runTime.timeName(),
-                mesh,
-                IOobject::NO_READ,
-                IOobject::NO_WRITE
-            ),
-            mesh,
-            dimensionSet(0, 1, -1, 0, 0)
-        );
-        tzz = sqrt(Tzz - (UMeanz*UMeanz));
-        tzz.write();
-
-        volScalarField txy
-        (
-            IOobject
-            (
-                "txy",
-                runTime.timeName(),
-                mesh,
-                IOobject::NO_READ,
-                IOobject::NO_WRITE
-            ),
-            mesh,
-            dimensionSet(0, 2, -2, 0, 0)
-        );
-        txy = Txy - (UMeanx*UMeany);
-        txy.write();
-*/
diff --git a/applications/utilities/postProcessing/miscellaneous/postChannel/channelIndex.C b/applications/utilities/postProcessing/miscellaneous/postChannel/channelIndex.C
deleted file mode 100644
index be0aa492c3..0000000000
--- a/applications/utilities/postProcessing/miscellaneous/postChannel/channelIndex.C
+++ /dev/null
@@ -1,293 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2011-2018 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 "channelIndex.H"
-#include "boolList.H"
-#include "syncTools.H"
-#include "OFstream.H"
-#include "meshTools.H"
-#include "Time.H"
-#include "SortableList.H"
-
-// * * * * * * * * * * * * * Static Member Data  * * * * * * * * * * * * * * //
-
-namespace Foam
-{
-    template<>
-    const char* Foam::NamedEnum
-    <
-        Foam::vector::components,
-        3
-    >::names[] =
-    {
-        "x",
-        "y",
-        "z"
-    };
-}
-
-const Foam::NamedEnum<Foam::vector::components, 3>
-    Foam::channelIndex::vectorComponentsNames_;
-
-
-// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
-
-// Determines face blocking
-void Foam::channelIndex::walkOppositeFaces
-(
-    const polyMesh& mesh,
-    const labelList& startFaces,
-    boolList& blockedFace
-)
-{
-    const cellList& cells = mesh.cells();
-    const faceList& faces = mesh.faces();
-    label nBnd = mesh.nFaces() - mesh.nInternalFaces();
-
-    DynamicList<label> frontFaces(startFaces);
-    forAll(frontFaces, i)
-    {
-        label facei = frontFaces[i];
-        blockedFace[facei] = true;
-    }
-
-    while (returnReduce(frontFaces.size(), sumOp<label>()) > 0)
-    {
-        // Transfer across.
-        boolList isFrontBndFace(nBnd, false);
-        forAll(frontFaces, i)
-        {
-            label facei = frontFaces[i];
-
-            if (!mesh.isInternalFace(facei))
-            {
-                isFrontBndFace[facei-mesh.nInternalFaces()] = true;
-            }
-        }
-        syncTools::swapBoundaryFaceList(mesh, isFrontBndFace);
-
-        // Add
-        forAll(isFrontBndFace, i)
-        {
-            label facei = mesh.nInternalFaces()+i;
-            if (isFrontBndFace[i] && !blockedFace[facei])
-            {
-                blockedFace[facei] = true;
-                frontFaces.append(facei);
-            }
-        }
-
-        // Transfer across cells
-        DynamicList<label> newFrontFaces(frontFaces.size());
-
-        forAll(frontFaces, i)
-        {
-            label facei = frontFaces[i];
-
-            {
-                const cell& ownCell = cells[mesh.faceOwner()[facei]];
-
-                label oppositeFacei = ownCell.opposingFaceLabel(facei, faces);
-
-                if (oppositeFacei == -1)
-                {
-                    FatalErrorInFunction
-                        << "Face:" << facei << " owner cell:" << ownCell
-                        << " is not a hex?" << abort(FatalError);
-                }
-                else
-                {
-                    if (!blockedFace[oppositeFacei])
-                    {
-                        blockedFace[oppositeFacei] = true;
-                        newFrontFaces.append(oppositeFacei);
-                    }
-                }
-            }
-
-            if (mesh.isInternalFace(facei))
-            {
-                const cell& neiCell = mesh.cells()[mesh.faceNeighbour()[facei]];
-
-                label oppositeFacei = neiCell.opposingFaceLabel(facei, faces);
-
-                if (oppositeFacei == -1)
-                {
-                    FatalErrorInFunction
-                        << "Face:" << facei << " neighbour cell:" << neiCell
-                        << " is not a hex?" << abort(FatalError);
-                }
-                else
-                {
-                    if (!blockedFace[oppositeFacei])
-                    {
-                        blockedFace[oppositeFacei] = true;
-                        newFrontFaces.append(oppositeFacei);
-                    }
-                }
-            }
-        }
-
-        frontFaces.transfer(newFrontFaces);
-    }
-}
-
-
-// Calculate regions.
-void Foam::channelIndex::calcLayeredRegions
-(
-    const polyMesh& mesh,
-    const labelList& startFaces
-)
-{
-    boolList blockedFace(mesh.nFaces(), false);
-    walkOppositeFaces
-    (
-        mesh,
-        startFaces,
-        blockedFace
-    );
-
-
-    if (false)
-    {
-        OFstream str(mesh.time().path()/"blockedFaces.obj");
-        label vertI = 0;
-        forAll(blockedFace, facei)
-        {
-            if (blockedFace[facei])
-            {
-                const face& f = mesh.faces()[facei];
-                forAll(f, fp)
-                {
-                    meshTools::writeOBJ(str, mesh.points()[f[fp]]);
-                }
-                str<< 'f';
-                forAll(f, fp)
-                {
-                    str << ' ' << vertI+fp+1;
-                }
-                str << nl;
-                vertI += f.size();
-            }
-        }
-    }
-
-
-    // Do analysis for connected regions
-    cellRegion_.reset(new regionSplit(mesh, blockedFace));
-
-    Info<< "Detected " << cellRegion_().nRegions() << " layers." << nl << endl;
-
-    // Sum number of entries per region
-    regionCount_ = regionSum(scalarField(mesh.nCells(), 1.0));
-
-    // Average cell centres to determine ordering.
-    pointField regionCc
-    (
-        regionSum(mesh.cellCentres())
-      / regionCount_
-    );
-
-    SortableList<scalar> sortComponent(regionCc.component(dir_));
-
-    sortMap_ = sortComponent.indices();
-
-    y_ = sortComponent;
-
-    if (symmetric_)
-    {
-        y_.setSize(cellRegion_().nRegions()/2);
-    }
-}
-
-
-// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
-
-Foam::channelIndex::channelIndex
-(
-    const polyMesh& mesh,
-    const dictionary& dict
-)
-:
-    symmetric_(readBool(dict.lookup("symmetric"))),
-    dir_(vectorComponentsNames_.read(dict.lookup("component")))
-{
-    const polyBoundaryMesh& patches = mesh.boundaryMesh();
-
-    const wordList patchNames(dict.lookup("patches"));
-
-    label nFaces = 0;
-
-    forAll(patchNames, i)
-    {
-        const label patchi = patches.findPatchID(patchNames[i]);
-
-        if (patchi == -1)
-        {
-            FatalErrorInFunction
-                << "Illegal patch " << patchNames[i]
-                << ". Valid patches are " << patches.name()
-                << exit(FatalError);
-        }
-
-        nFaces += patches[patchi].size();
-    }
-
-    labelList startFaces(nFaces);
-    nFaces = 0;
-
-    forAll(patchNames, i)
-    {
-        const polyPatch& pp = patches[patchNames[i]];
-
-        forAll(pp, j)
-        {
-            startFaces[nFaces++] = pp.start()+j;
-        }
-    }
-
-    // Calculate regions.
-    calcLayeredRegions(mesh, startFaces);
-}
-
-
-Foam::channelIndex::channelIndex
-(
-    const polyMesh& mesh,
-    const labelList& startFaces,
-    const bool symmetric,
-    const direction dir
-)
-:
-    symmetric_(symmetric),
-    dir_(dir)
-{
-    // Calculate regions.
-    calcLayeredRegions(mesh, startFaces);
-}
-
-
-// ************************************************************************* //
diff --git a/applications/utilities/postProcessing/miscellaneous/postChannel/channelIndex.H b/applications/utilities/postProcessing/miscellaneous/postChannel/channelIndex.H
deleted file mode 100644
index 5559608b97..0000000000
--- a/applications/utilities/postProcessing/miscellaneous/postChannel/channelIndex.H
+++ /dev/null
@@ -1,160 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2011-2020 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/>.
-
-Class
-    Foam::channelIndex
-
-Description
-    Does averaging of fields over layers of cells. Assumes layered mesh.
-
-SourceFiles
-    channelIndex.C
-
-\*---------------------------------------------------------------------------*/
-
-#ifndef channelIndex_H
-#define channelIndex_H
-
-#include "regionSplit.H"
-#include "direction.H"
-#include "scalarField.H"
-#include "polyMesh.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-namespace Foam
-{
-
-/*---------------------------------------------------------------------------*\
-                        Class channelIndex Declaration
-\*---------------------------------------------------------------------------*/
-
-class channelIndex
-{
-
-    // Private Data
-
-        static const NamedEnum<vector::components, 3> vectorComponentsNames_;
-
-        //- Is mesh symmetric
-        const bool symmetric_;
-
-        //- Direction to sort
-        const direction dir_;
-
-        //- Per cell the global region
-        autoPtr<regionSplit> cellRegion_;
-
-        //- Per global region the number of cells (scalarField so we can use
-        //  field algebra)
-        scalarField regionCount_;
-
-        //- From sorted region back to unsorted global region
-        labelList sortMap_;
-
-        //- Sorted component of cell centres
-        scalarField y_;
-
-
-
-    // Private Member Functions
-
-        void walkOppositeFaces
-        (
-            const polyMesh& mesh,
-            const labelList& startFaces,
-            boolList& blockedFace
-        );
-
-        void calcLayeredRegions
-        (
-            const polyMesh& mesh,
-            const labelList& startFaces
-        );
-
-
-public:
-
-    // Constructors
-
-        //- Construct from dictionary
-        channelIndex(const polyMesh&, const dictionary&);
-
-        //- Construct from supplied starting faces
-        channelIndex
-        (
-            const polyMesh& mesh,
-            const labelList& startFaces,
-            const bool symmetric,
-            const direction dir
-        );
-
-        //- Disallow default bitwise copy construction
-        channelIndex(const channelIndex&) = delete;
-
-
-    // Member Functions
-
-        // Access
-
-            //- Sum field per region
-            template<class T>
-            Field<T> regionSum(const Field<T>& cellField) const;
-
-            //- Collapse a field to a line
-            template<class T>
-            Field<T> collapse
-            (
-                const Field<T>& vsf,
-                const bool asymmetric=false
-            ) const;
-
-            //- Return the field of Y locations from the cell centres
-            const scalarField& y() const
-            {
-                return y_;
-            }
-
-
-    // Member Operators
-
-        //- Disallow default bitwise assignment
-        void operator=(const channelIndex&);
-};
-
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-} // End namespace Foam
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-#ifdef NoRepository
-    #include "channelIndexTemplates.C"
-#endif
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-#endif
-
-// ************************************************************************* //
diff --git a/applications/utilities/postProcessing/miscellaneous/postChannel/channelIndexTemplates.C b/applications/utilities/postProcessing/miscellaneous/postChannel/channelIndexTemplates.C
deleted file mode 100644
index 157172015f..0000000000
--- a/applications/utilities/postProcessing/miscellaneous/postChannel/channelIndexTemplates.C
+++ /dev/null
@@ -1,102 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2011-2018 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 "channelIndex.H"
-
-// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
-
-template<class T>
-Foam::Field<T> Foam::channelIndex::regionSum(const Field<T>& cellField) const
-{
-    Field<T> regionField(cellRegion_().nRegions(), Zero);
-
-    forAll(cellRegion_(), celli)
-    {
-        regionField[cellRegion_()[celli]] += cellField[celli];
-    }
-
-    // Global sum
-    Pstream::listCombineGather(regionField, plusEqOp<T>());
-    Pstream::listCombineScatter(regionField);
-
-    return regionField;
-}
-
-
-template<class T>
-Foam::Field<T> Foam::channelIndex::collapse
-(
-    const Field<T>& cellField,
-    const bool asymmetric
-) const
-{
-    // Average and order
-    const Field<T> summedField(regionSum(cellField));
-
-    Field<T> regionField
-    (
-        summedField
-      / regionCount_,
-        sortMap_
-    );
-
-    // Symmetry?
-    if (symmetric_)
-    {
-        label nlb2 = cellRegion_().nRegions()/2;
-
-        if (asymmetric)
-        {
-            for (label j=0; j<nlb2; j++)
-            {
-                regionField[j] =
-                    0.5
-                  * (
-                        regionField[j]
-                      - regionField[cellRegion_().nRegions() - j - 1]
-                    );
-            }
-        }
-        else
-        {
-            for (label j=0; j<nlb2; j++)
-            {
-                regionField[j] =
-                    0.5
-                  * (
-                        regionField[j]
-                      + regionField[cellRegion_().nRegions() - j - 1]
-                    );
-            }
-        }
-
-        regionField.setSize(nlb2);
-    }
-
-    return regionField;
-}
-
-
-// ************************************************************************* //
diff --git a/applications/utilities/postProcessing/miscellaneous/postChannel/collapse.H b/applications/utilities/postProcessing/miscellaneous/postChannel/collapse.H
deleted file mode 100644
index d7afa04147..0000000000
--- a/applications/utilities/postProcessing/miscellaneous/postChannel/collapse.H
+++ /dev/null
@@ -1,69 +0,0 @@
-    scalarField UMeanXvalues
-    (
-        channelIndexing.collapse(UMean.component(vector::X)())
-    );
-
-    scalarField UMeanYvalues
-    (
-        channelIndexing.collapse(UMean.component(vector::Y)())
-    );
-
-    scalarField UMeanZvalues
-    (
-        channelIndexing.collapse(UMean.component(vector::Z)())
-    );
-
-    scalarField RxxValues(channelIndexing.collapse(Rxx));
-    scalarField RyyValues(channelIndexing.collapse(Ryy));
-    scalarField RzzValues(channelIndexing.collapse(Rzz));
-    scalarField RxyValues(channelIndexing.collapse(Rxy, true));
-
-    scalarField pPrime2MeanValues(channelIndexing.collapse(pPrime2Mean));
-
-    /*
-    scalarField epsilonValues(channelIndexing.collapse(epsilonMean));
-
-    scalarField nuMeanValues(channelIndexing.collapse(nuMean));
-    scalarField nuPrimeValues(channelIndexing.collapse(nuPrime));
-
-    scalarField gammaDotMeanValues(channelIndexing.collapse(gammaDotMean));
-    scalarField gammaDotPrimeValues(channelIndexing.collapse(gammaDotPrime));
-    */
-
-    scalarField urmsValues(sqrt(mag(RxxValues)));
-    scalarField vrmsValues(sqrt(mag(RyyValues)));
-    scalarField wrmsValues(sqrt(mag(RzzValues)));
-
-    scalarField kValues
-    (
-        0.5*(sqr(urmsValues) + sqr(vrmsValues) + sqr(wrmsValues))
-    );
-
-    setWriter::New(runTime.graphFormat())->write
-    (
-        runTime.globalPath()
-       /functionObjects::writeFile::outputPrefix
-       /args.executable()
-       /runTime.timeName(),
-
-       args.executable(),
-
-       coordSet(true, "y", channelIndexing.y()),
-
-       "Uf", UMeanXvalues,
-
-       "u", urmsValues,
-       "v", vrmsValues,
-       "w", wrmsValues,
-       "uv", RxyValues,
-
-       "k", kValues,
-
-       "pPrime2Mean", pPrime2MeanValues
-
-       // "epsilon", epsilonValues,
-       // "nu", nuMeanValues,
-       // "nuPrime", nuPrimeValues,
-       // "gammaDot", gammaDotMeanValues,
-       // "gammaDotPrime", gammaDotPrimeValues
-    );
diff --git a/applications/utilities/postProcessing/miscellaneous/postChannel/postChannel.C b/applications/utilities/postProcessing/miscellaneous/postChannel/postChannel.C
deleted file mode 100644
index a56a8a73b8..0000000000
--- a/applications/utilities/postProcessing/miscellaneous/postChannel/postChannel.C
+++ /dev/null
@@ -1,97 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2011-2021 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/>.
-
-Application
-    postChannel
-
-Description
-    Post-processes data from channel flow calculations.
-
-    For each time: calculate: txx, txy,tyy, txy,
-    eps, prod, vorticity, enstrophy and helicity. Assuming that the mesh
-    is periodic in the x and z directions, collapse Umeanx, Umeany, txx,
-    txy and tyy to a line and print them as standard output.
-
-\*---------------------------------------------------------------------------*/
-
-#include "argList.H"
-#include "timeSelector.H"
-#include "volFields.H"
-#include "channelIndex.H"
-#include "setWriter.H"
-#include "writeFile.H"
-
-using namespace Foam;
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-int main(int argc, char *argv[])
-{
-    argList::noParallel();
-    timeSelector::addOptions();
-
-    #include "setRootCase.H"
-    #include "createTime.H"
-
-    // Get times list
-    instantList timeDirs = timeSelector::select0(runTime, args);
-
-    #include "createMeshNoChangers.H"
-    #include "readPhysicalProperties.H"
-
-    // Setup channel indexing for averaging over channel down to a line
-
-    IOdictionary channelDict
-    (
-        IOobject
-        (
-            "postChannelDict",
-            mesh.time().system(),
-            mesh,
-            IOobject::MUST_READ_IF_MODIFIED,
-            IOobject::NO_WRITE
-        )
-    );
-    channelIndex channelIndexing(mesh, channelDict);
-
-
-    // For each time step read all fields
-    forAll(timeDirs, timeI)
-    {
-        runTime.setTime(timeDirs[timeI], timeI);
-        Info<< "Collapsing fields for time " << runTime.timeName() << endl;
-
-        #include "readFields.H"
-        #include "calculateFields.H"
-
-        // Average fields over channel down to a line
-        #include "collapse.H"
-    }
-
-    Info<< "\nEnd\n" << endl;
-
-    return 0;
-}
-
-
-// ************************************************************************* //
diff --git a/applications/utilities/postProcessing/miscellaneous/postChannel/readFields.H b/applications/utilities/postProcessing/miscellaneous/postChannel/readFields.H
deleted file mode 100644
index b257d73645..0000000000
--- a/applications/utilities/postProcessing/miscellaneous/postChannel/readFields.H
+++ /dev/null
@@ -1,112 +0,0 @@
-    typeIOobject<volVectorField> UMeanHeader
-    (
-        "UMean",
-        runTime.timeName(),
-        mesh,
-        IOobject::MUST_READ
-    );
-
-    if (!UMeanHeader.headerOk())
-    {
-        Info<< "    No UMean field" << endl;
-        continue;
-    }
-
-    volVectorField UMean
-    (
-        UMeanHeader,
-        mesh
-    );
-
-    volSymmTensorField UPrime2Mean
-    (
-        IOobject
-        (
-            "UPrime2Mean",
-            runTime.timeName(),
-            mesh,
-            IOobject::MUST_READ
-        ),
-        mesh
-    );
-    volScalarField Rxx(UPrime2Mean.component(symmTensor::XX));
-    volScalarField Ryy(UPrime2Mean.component(symmTensor::YY));
-    volScalarField Rzz(UPrime2Mean.component(symmTensor::ZZ));
-    volScalarField Rxy(UPrime2Mean.component(symmTensor::XY));
-
-    volScalarField pPrime2Mean
-    (
-        IOobject
-        (
-            "pPrime2Mean",
-            runTime.timeName(),
-            mesh,
-            IOobject::MUST_READ
-        ),
-        mesh
-    );
-
-    /*
-    volScalarField epsilonMean
-    (
-        IOobject
-        (
-            "epsilonMean",
-            runTime.timeName(),
-            mesh,
-            IOobject::MUST_READ
-        ),
-        mesh
-    );
-
-    volScalarField nuMean
-    (
-        IOobject
-        (
-            "nuMean",
-            runTime.timeName(),
-            mesh,
-            IOobject::MUST_READ
-        ),
-        mesh
-    );
-
-    volScalarField gammaDotMean
-    (
-        IOobject
-        (
-            "gammaDotMean",
-            runTime.timeName(),
-            mesh,
-            IOobject::MUST_READ
-        ),
-        mesh
-    );
-
-    volScalarField nuPrime2
-    (
-        IOobject
-        (
-            "nuPrime",
-            runTime.timeName(),
-            mesh,
-            IOobject::MUST_READ
-        ),
-        mesh
-    );
-    volScalarField nuPrime = sqrt(mag(nuPrime2 - sqr(nuMean)));
-
-
-    volScalarField gammaDotPrime2
-    (
-        IOobject
-        (
-            "gammaDotPrime",
-            runTime.timeName(),
-            mesh,
-            IOobject::MUST_READ
-        ),
-        mesh
-    );
-    volScalarField gammaDotPrime = sqrt(mag(gammaDotPrime2 -sqr(gammaDotMean)));
-    */
diff --git a/applications/utilities/postProcessing/miscellaneous/postChannel/readPhysicalProperties.H b/applications/utilities/postProcessing/miscellaneous/postChannel/readPhysicalProperties.H
deleted file mode 100644
index 0c7a609890..0000000000
--- a/applications/utilities/postProcessing/miscellaneous/postChannel/readPhysicalProperties.H
+++ /dev/null
@@ -1,13 +0,0 @@
-    Info<< nl << "Reading physicalProperties" << endl;
-
-    IOdictionary transportProperties
-    (
-        IOobject
-        (
-            "physicalProperties",
-            runTime.constant(),
-            mesh,
-            IOobject::MUST_READ_IF_MODIFIED,
-            IOobject::NO_WRITE
-        )
-    );
diff --git a/etc/caseDicts/postProcessing/graphs/graphLayerAverage b/etc/caseDicts/postProcessing/graphs/graphLayerAverage
new file mode 100644
index 0000000000..c114db9a79
--- /dev/null
+++ b/etc/caseDicts/postProcessing/graphs/graphLayerAverage
@@ -0,0 +1,26 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     | Website:  https://openfoam.org
+    \\  /    A nd           | Version:  dev
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+Description
+    Generates plots of fields averaged over the layers in the mesh
+
+\*---------------------------------------------------------------------------*/
+
+patches         (<patchNames>);  // Patches from which layers extrude
+
+zones           ();              // Zones from which layers extrude
+
+axis            <axisName>;      // The independent variable of the graph. Can
+                                 // be "x", "y" or "z".
+
+symmetric       false;           // Are the layers symmetric about the centre?
+
+fields          (<fieldsNames>); // Fields to plot
+
+#includeEtc "caseDicts/postProcessing/graphs/graphLayerAverage.cfg"
+
+// ************************************************************************* //
diff --git a/etc/caseDicts/postProcessing/graphs/graphLayerAverage.cfg b/etc/caseDicts/postProcessing/graphs/graphLayerAverage.cfg
new file mode 100644
index 0000000000..85e58beb64
--- /dev/null
+++ b/etc/caseDicts/postProcessing/graphs/graphLayerAverage.cfg
@@ -0,0 +1,17 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     | Website:  https://openfoam.org
+    \\  /    A nd           | Version:  dev
+     \\/     M anipulation  |
+\*---------------------------------------------------------------------------*/
+
+type            layerAverage;
+libs            ("libfieldFunctionObjects.so");
+
+writeControl    writeTime;
+writeInterval   1;
+
+setFormat       raw;
+
+// ************************************************************************* //
diff --git a/src/functionObjects/field/Make/files b/src/functionObjects/field/Make/files
index 8e10e05b2d..0b9a05330a 100644
--- a/src/functionObjects/field/Make/files
+++ b/src/functionObjects/field/Make/files
@@ -81,4 +81,6 @@ cylindrical/cylindricalFunctionObject.C
 
 uniform/uniform.C
 
+layerAverage/layerAverage.C
+
 LIB = $(FOAM_LIBBIN)/libfieldFunctionObjects
diff --git a/src/functionObjects/field/layerAverage/layerAverage.C b/src/functionObjects/field/layerAverage/layerAverage.C
new file mode 100644
index 0000000000..7304a53687
--- /dev/null
+++ b/src/functionObjects/field/layerAverage/layerAverage.C
@@ -0,0 +1,420 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     | Website:  https://openfoam.org
+    \\  /    A nd           | Copyright (C) 2011-2021 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 "layerAverage.H"
+#include "regionSplit.H"
+#include "syncTools.H"
+#include "volFields.H"
+#include "writeFile.H"
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    namespace functionObjects
+    {
+        defineTypeNameAndDebug(layerAverage, 0);
+        addToRunTimeSelectionTable(functionObject, layerAverage, dictionary);
+    }
+
+    template<>
+    const char*
+        Foam::NamedEnum<Foam::vector::components, 3>::names[] =
+        {"x", "y", "z"};
+}
+
+const Foam::NamedEnum<Foam::vector::components, 3>
+    Foam::functionObjects::layerAverage::axisNames_;
+
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+void Foam::functionObjects::layerAverage::walkOppositeFaces
+(
+    const labelList& startFaces,
+    const boolList& startFaceIntoOwners,
+    boolList& blockedFace,
+    boolList& cellIsLayer
+)
+{
+    // Initialise the front
+    DynamicList<label> frontFaces(startFaces);
+    DynamicList<bool> frontFaceIntoOwners(startFaceIntoOwners);
+    DynamicList<label> newFrontFaces(frontFaces.size());
+    DynamicList<bool> newFrontFaceIntoOwners(frontFaceIntoOwners.size());
+
+    // Set the start and end faces as blocked
+    UIndirectList<bool>(blockedFace, startFaces) = true;
+
+    // Iterate until the front is empty
+    while (returnReduce(frontFaces.size(), sumOp<label>()) > 0)
+    {
+        // Transfer front faces across couplings
+        boolList bndFaceIsFront(mesh_.nFaces() - mesh_.nInternalFaces(), false);
+        forAll(frontFaces, i)
+        {
+            const label facei = frontFaces[i];
+            const label intoOwner = frontFaceIntoOwners[i];
+
+            if (!mesh_.isInternalFace(facei) && !intoOwner)
+            {
+                bndFaceIsFront[facei - mesh_.nInternalFaces()] = true;
+            }
+        }
+        syncTools::swapBoundaryFaceList(mesh_, bndFaceIsFront);
+
+        // Set new front faces
+        forAll(bndFaceIsFront, i)
+        {
+            const label facei = mesh_.nInternalFaces() + i;
+
+            if (bndFaceIsFront[i] && !blockedFace[facei])
+            {
+                frontFaces.append(facei);
+                frontFaceIntoOwners.append(true);
+            }
+        }
+
+        // Transfer across cells
+        newFrontFaces.clear();
+        newFrontFaceIntoOwners.clear();
+
+        forAll(frontFaces, i)
+        {
+            const label facei = frontFaces[i];
+            const label intoOwner = frontFaceIntoOwners[i];
+
+            const label celli =
+                intoOwner
+              ? mesh_.faceOwner()[facei]
+              : mesh_.isInternalFace(facei) ? mesh_.faceNeighbour()[facei] : -1;
+
+            if (celli != -1)
+            {
+                cellIsLayer[celli] = true;
+
+                const label oppositeFacei =
+                    mesh_.cells()[celli]
+                   .opposingFaceLabel(facei, mesh_.faces());
+                const bool oppositeIntoOwner =
+                    mesh_.faceOwner()[oppositeFacei] != celli;
+
+                if (oppositeFacei == -1)
+                {
+                    FatalErrorInFunction
+                        << "Cannot find face on cell " << mesh_.cells()[celli]
+                        << " opposing face " << mesh_.faces()[facei]
+                        << ". Mesh is not layered." << exit(FatalError);
+                }
+                else if (!blockedFace[oppositeFacei])
+                {
+                    blockedFace[oppositeFacei] = true;
+                    newFrontFaces.append(oppositeFacei);
+                    newFrontFaceIntoOwners.append(oppositeIntoOwner);
+                }
+            }
+        }
+
+        frontFaces.transfer(newFrontFaces);
+        frontFaceIntoOwners.transfer(newFrontFaceIntoOwners);
+    }
+}
+
+
+void Foam::functionObjects::layerAverage::calcLayers()
+{
+    DynamicList<label> startFaces;
+    DynamicList<bool> startFaceIntoOwners;
+
+    forAll(patchIDs_, i)
+    {
+        const polyPatch& pp = mesh_.boundaryMesh()[patchIDs_[i]];
+        startFaces.append(identity(pp.size()) + pp.start());
+        startFaceIntoOwners.append(boolList(pp.size(), true));
+    }
+
+    forAll(zoneIDs_, i)
+    {
+        const faceZone& fz = mesh_.faceZones()[zoneIDs_[i]];
+        startFaces.append(fz);
+        startFaceIntoOwners.append(fz.flipMap());
+    }
+
+    // Identify faces that separate the layers
+    boolList blockedFace(mesh_.nFaces(), false);
+    boolList cellIsLayer(mesh_.nCells(), false);
+    walkOppositeFaces
+    (
+        startFaces,
+        startFaceIntoOwners,
+        blockedFace,
+        cellIsLayer
+    );
+
+    // Do analysis for connected layers
+    regionSplit rs(mesh_, blockedFace);
+    nLayers_ = rs.nRegions();
+    cellLayer_.transfer(rs);
+
+    // Get rid of regions that are not layers
+    label layeri0 = labelMax, layeri1 = -labelMax;
+    forAll(cellLayer_, celli)
+    {
+        if (cellIsLayer[celli])
+        {
+            layeri0 = min(layeri0, cellLayer_[celli]);
+            layeri1 = max(layeri1, cellLayer_[celli]);
+        }
+        else
+        {
+            cellLayer_[celli] = -1;
+        }
+    }
+    reduce(layeri0, minOp<label>());
+    reduce(layeri1, maxOp<label>());
+    nLayers_ = layeri0 != labelMax ? 1 + layeri1 - layeri0 : 0;
+    forAll(cellLayer_, celli)
+    {
+        if (cellLayer_[celli] != -1)
+        {
+            cellLayer_[celli] -= layeri0;
+        }
+    }
+
+    // Report
+    if (nLayers_ != 0)
+    {
+        Info<< "    Detected " << nLayers_ << " layers" << nl << endl;
+    }
+    else
+    {
+        WarningInFunction
+            << "No layers detected" << endl;
+    }
+
+    // Sum number of entries per layer
+    layerCount_ = sum(scalarField(mesh_.nCells(), 1));
+
+    // Average the cell centres
+    const pointField layerCentres(sum(mesh_.cellCentres())/layerCount_);
+
+    // Sort by the direction component
+    x_ = layerCentres.component(axis_);
+    sortMap_ = identity(layerCentres.size());
+    stableSort(sortMap_, UList<scalar>::less(x_));
+    inplaceReorder(sortMap_, x_);
+
+    // If symmetric, keep only half of the coordinates
+    if (symmetric_)
+    {
+        x_.setSize(nLayers_/2);
+    }
+}
+
+
+template<>
+Foam::vector
+Foam::functionObjects::layerAverage::symmetricCoeff<Foam::vector>() const
+{
+    vector result = vector::one;
+    result[axis_] = -1;
+    return result;
+}
+
+
+template<>
+Foam::symmTensor
+Foam::functionObjects::layerAverage::symmetricCoeff<Foam::symmTensor>() const
+{
+    return sqr(symmetricCoeff<vector>());
+}
+
+
+template<>
+Foam::tensor
+Foam::functionObjects::layerAverage::symmetricCoeff<Foam::tensor>() const
+{
+    return symmetricCoeff<symmTensor>();
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::functionObjects::layerAverage::layerAverage
+(
+    const word& name,
+    const Time& runTime,
+    const dictionary& dict
+)
+:
+    fvMeshFunctionObject(name, runTime, dict)
+{
+    read(dict);
+}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::functionObjects::layerAverage::~layerAverage()
+{}
+
+
+// * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
+
+bool Foam::functionObjects::layerAverage::read(const dictionary& dict)
+{
+    Info<< type() << " " << name() << ":" << nl;
+
+    patchIDs_ =
+        mesh_.boundaryMesh().patchSet
+        (
+            dict.lookupOrDefault<wordReList>("patches", wordReList())
+        ).toc();
+
+    zoneIDs_ =
+        findStrings
+        (
+            dict.lookupOrDefault<wordReList>("zones", wordReList()),
+            mesh_.faceZones().names()
+        );
+
+    if (patchIDs_.empty() && zoneIDs_.empty())
+    {
+        WarningInFunction
+            << "No patches or zones specified" << endl;
+    }
+
+    symmetric_ = dict.lookupOrDefault<bool>("symmetric", false);
+
+    axis_ = axisNames_.read(dict.lookup("axis"));
+
+    fields_ = dict.lookup<wordList>("fields");
+
+    formatter_ = setWriter::New(dict.lookup("setFormat"), dict);
+
+    calcLayers();
+
+    return true;
+}
+
+
+Foam::wordList Foam::functionObjects::layerAverage::fields() const
+{
+    return fields_;
+}
+
+
+bool Foam::functionObjects::layerAverage::execute()
+{
+    return true;
+}
+
+
+bool Foam::functionObjects::layerAverage::write()
+{
+    // Create list of available fields
+    wordList fieldNames;
+    forAll(fields_, fieldi)
+    {
+        if
+        (
+            false
+            #define FoundTypeField(Type, nullArg) \
+              || foundObject<VolField<Type>>(fields_[fieldi])
+            FOR_ALL_FIELD_TYPES(FoundTypeField)
+            #undef FoundTypeField
+        )
+        {
+            fieldNames.append(fields_[fieldi]);
+        }
+        else
+        {
+            cannotFindObject(fields_[fieldi]);
+        }
+    }
+
+    // Collapse the fields
+    #define DeclareTypeValueSets(Type, nullArg) \
+        PtrList<Field<Type>> Type##ValueSets(fieldNames.size());
+    FOR_ALL_FIELD_TYPES(DeclareTypeValueSets);
+    #undef DeclareTypeValueSets
+    forAll(fieldNames, fieldi)
+    {
+        #define CollapseTypeFields(Type, nullArg)                           \
+            if (mesh_.foundObject<VolField<Type>>(fieldNames[fieldi]))      \
+            {                                                               \
+                const VolField<Type>& field =                               \
+                    mesh_.lookupObject<VolField<Type>>(fieldNames[fieldi]); \
+                                                                            \
+                Type##ValueSets.set                                         \
+                (                                                           \
+                    fieldi,                                                 \
+                    average(field.primitiveField())                         \
+                );                                                          \
+            }
+        FOR_ALL_FIELD_TYPES(CollapseTypeFields);
+        #undef CollapseTypeFields
+    }
+
+    // Output directory
+    fileName outputPath =
+        mesh_.time().globalPath()/writeFile::outputPrefix/name();
+    if (mesh_.name() != fvMesh::defaultRegion)
+    {
+        outputPath = outputPath/mesh_.name();
+    }
+
+    // Write
+    formatter_->write
+    (
+        outputPath/mesh_.time().timeName(),
+        typeName,
+        coordSet(true, axisNames_[axis_], x_),
+        fieldNames
+        #define TypeValueSetsParameter(Type, nullArg) , Type##ValueSets
+        FOR_ALL_FIELD_TYPES(TypeValueSetsParameter)
+        #undef TypeValueSetsParameter
+    );
+
+    return true;
+}
+
+
+void Foam::functionObjects::layerAverage::updateMesh(const mapPolyMesh&)
+{
+    Info<< type() << " " << name() << ":" << nl;
+    calcLayers();
+}
+
+
+void Foam::functionObjects::layerAverage::movePoints(const polyMesh&)
+{
+    Info<< type() << " " << name() << ":" << nl;
+    calcLayers();
+}
+
+
+// ************************************************************************* //
diff --git a/src/functionObjects/field/layerAverage/layerAverage.H b/src/functionObjects/field/layerAverage/layerAverage.H
new file mode 100644
index 0000000000..8ca30fdf57
--- /dev/null
+++ b/src/functionObjects/field/layerAverage/layerAverage.H
@@ -0,0 +1,243 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     | Website:  https://openfoam.org
+    \\  /    A nd           | Copyright (C) 2011-2021 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/>.
+
+Class
+    Foam::layerAverage
+
+Description
+    Generates plots of fields averaged over the layers in the mesh
+
+    Example of function object specification:
+    \verbatim
+    layerAverage1
+    {
+        type            layerAverage;
+        libs            ("libfieldFunctionObjects.so");
+
+        writeControl    writeTime;
+
+        setFormat       raw;
+
+        patches         (bottom);
+        zones           (quarterPlane threeQuartersPlane);
+
+        component       y;
+
+        symmetric       true;
+
+        fields          (pMean pPrime2Mean UMean UPrime2Mean k);
+    }
+    \endverbatim
+
+Usage
+    \table
+        Property      | Description                | Required   | Default value
+        type          | Type name: layerAverage    | yes        |
+        setFormat     | Output plotting format     | yes        |
+        patches       | Patches that layers extrude from | no | ()
+        zones         | Face zones that the layers extrude from | no | ()
+        component     | Component of position to plot against | yes |
+        symmetric     | Is the geometry symmetric around the centre layer? \
+                                                   | no         | false
+        field         | Fields to average and plot | yes        |
+    \endtable
+
+SourceFiles
+    layerAverage.C
+    layerAverageTemplates.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef layerAverage_H
+#define layerAverage_H
+
+#include "fvMeshFunctionObject.H"
+#include "setWriter.H"
+#include "boolList.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace functionObjects
+{
+
+/*---------------------------------------------------------------------------*\
+                          Class layerAverage Declaration
+\*---------------------------------------------------------------------------*/
+
+class layerAverage
+:
+    public fvMeshFunctionObject
+{
+    // Private Static Data
+
+        //- Names of vector components
+        static const NamedEnum<vector::components, 3> axisNames_;
+
+
+    // Private Data
+
+        //- Patches which form the start of the layers
+        labelList patchIDs_;
+
+        //- Zones which form the start of the layers
+        labelList zoneIDs_;
+
+        //- Zones on which the layers are considered to end
+        labelList endZoneIDs_;
+
+        //- Is the case symmetric?
+        bool symmetric_;
+
+        //- The direction over which to plot the results
+        vector::components axis_;
+
+        //- Per cell the global layer
+        label nLayers_;
+
+        //- Per cell the global layer
+        labelList cellLayer_;
+
+        //- Per global layer the number of cells
+        scalarField layerCount_;
+
+        //- From sorted layer back to unsorted global layer
+        labelList sortMap_;
+
+        //- Sorted component of cell centres
+        scalarField x_;
+
+        //- Fields to sample
+        wordList fields_;
+
+        //- Set formatter
+        autoPtr<setWriter> formatter_;
+
+
+    // Private Member Functions
+
+        //- Walk through layers marking faces that separate layers and cells
+        //  that are within layers
+        void walkOppositeFaces
+        (
+            const labelList& startFaces,
+            const boolList& startFaceIntoOwners,
+            boolList& blockedFace,
+            boolList& cellIsLayer
+        );
+
+        //- Create the layer information, the sort map, and the scalar axis
+        void calcLayers();
+
+        //- Return the coefficient to multiply onto symmetric values
+        template<class T>
+        T symmetricCoeff() const;
+
+        //- Sum field per layer
+        template<class T>
+        Field<T> sum(const Field<T>& cellField) const;
+
+        //- Average a field per layer
+        template<class T>
+        Field<T> average(const Field<T>& cellField) const;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("layerAverage");
+
+
+    // Constructors
+
+        //- Construct from Time and dictionary
+        layerAverage
+        (
+            const word& name,
+            const Time& runTime,
+            const dictionary& dict
+        );
+
+        //- Disallow default bitwise copy construction
+        layerAverage(const layerAverage&) = delete;
+
+
+    //- Destructor
+    virtual ~layerAverage();
+
+
+    // Member Functions
+
+        //- Read the field average data
+        virtual bool read(const dictionary&);
+
+        //- Return the list of fields required
+        virtual wordList fields() const;
+
+        //- Do nothing
+        virtual bool execute();
+
+        //- Calculate and write the graphs
+        virtual bool write();
+
+        //- Update for changes of mesh
+        virtual void updateMesh(const mapPolyMesh&);
+
+        //- Update for mesh point-motion
+        virtual void movePoints(const polyMesh&);
+
+
+    // Member Operators
+
+        //- Disallow default bitwise assignment
+        void operator=(const layerAverage&) = delete;
+};
+
+
+template<>
+vector layerAverage::symmetricCoeff<vector>() const;
+
+template<>
+symmTensor layerAverage::symmetricCoeff<symmTensor>() const;
+
+template<>
+tensor layerAverage::symmetricCoeff<tensor>() const;
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace functionObjects
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+    #include "layerAverageTemplates.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/functionObjects/field/layerAverage/layerAverageTemplates.C b/src/functionObjects/field/layerAverage/layerAverageTemplates.C
new file mode 100644
index 0000000000..1c897c5977
--- /dev/null
+++ b/src/functionObjects/field/layerAverage/layerAverageTemplates.C
@@ -0,0 +1,89 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     | Website:  https://openfoam.org
+    \\  /    A nd           | Copyright (C) 2011-2021 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 "layerAverage.H"
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+template<class T>
+T Foam::functionObjects::layerAverage::symmetricCoeff() const
+{
+    return pTraits<T>::one;
+}
+
+
+template<class T>
+Foam::Field<T> Foam::functionObjects::layerAverage::sum
+(
+    const Field<T>& cellField
+) const
+{
+    Field<T> layerField(nLayers_, Zero);
+
+    forAll(cellLayer_, celli)
+    {
+        if (cellLayer_[celli] != -1)
+        {
+            layerField[cellLayer_[celli]] += cellField[celli];
+        }
+    }
+
+    Pstream::listCombineGather(layerField, plusEqOp<T>());
+    Pstream::listCombineScatter(layerField);
+
+    return layerField;
+}
+
+
+template<class T>
+Foam::Field<T> Foam::functionObjects::layerAverage::average
+(
+    const Field<T>& cellField
+) const
+{
+    // Sum, average and sort
+    Field<T> layerField(sum(cellField)/layerCount_, sortMap_);
+
+    // Handle symmetry
+    if (symmetric_)
+    {
+        const T coeff = symmetricCoeff<T>();
+
+        for (label i=0; i<nLayers_/2; i++)
+        {
+            const label j = nLayers_ - i - 1;
+
+            layerField[i] =
+                (layerField[i] + cmptMultiply(coeff, layerField[j]))/2;
+        }
+
+        layerField.setSize(nLayers_/2);
+    }
+
+    return layerField;
+}
+
+
+// ************************************************************************* //
diff --git a/test/postProcessing/channel/system/controlDict b/test/postProcessing/channel/system/controlDict
index a4bd4cd5c1..769810535a 100644
--- a/test/postProcessing/channel/system/controlDict
+++ b/test/postProcessing/channel/system/controlDict
@@ -248,6 +248,13 @@ graphsFunctions
         end=(0.500001 0.500001 0),
         fields=(p U)
     )
+
+    #includeFunc graphLayerAverage
+    (
+        patches=(inlet),
+        axis=x,
+        fields=(p U)
+    )
 }
 
 lagrangianFunctions
diff --git a/tutorials/incompressible/pimpleFoam/LES/channel395/Allrun b/tutorials/incompressible/pimpleFoam/LES/channel395/Allrun
index c34f238440..0ae462d1ad 100755
--- a/tutorials/incompressible/pimpleFoam/LES/channel395/Allrun
+++ b/tutorials/incompressible/pimpleFoam/LES/channel395/Allrun
@@ -17,6 +17,4 @@ runApplication decomposePar -cellDist
 runParallel $application
 runApplication reconstructPar
 
-runApplication postChannel
-
 #------------------------------------------------------------------------------
diff --git a/tutorials/incompressible/pimpleFoam/LES/channel395/system/controlDict b/tutorials/incompressible/pimpleFoam/LES/channel395/system/controlDict
index 8600a4e8b8..fe890d08ee 100644
--- a/tutorials/incompressible/pimpleFoam/LES/channel395/system/controlDict
+++ b/tutorials/incompressible/pimpleFoam/LES/channel395/system/controlDict
@@ -47,6 +47,31 @@ runTimeModifiable true;
 functions
 {
     #includeFunc fieldAverage(U, p, prime2Mean = yes)
+
+    #includeFunc uniform
+    (
+        fieldType = volScalarField,
+        name = half,
+        dimensions = [0 0 0 0 0 0 0],
+        value = 0.5
+    )
+
+    #includeFunc mag(UPrime2Mean)
+
+    #includeFunc multiply(half, mag(UPrime2Mean), result = k)
+
+    #includeFunc graphLayerAverage
+    (
+        funcName = layerAverage,
+        patches = (bottomWall),
+        axis = y,
+        symmetric = yes,
+        pMean,
+        pPrime2Mean,
+        UMean,
+        UPrime2Mean,
+        k
+    );
 }
 
 // ************************************************************************* //