ENH: Added new function object to calculate the energy spectrum

Description
    Calculates the energy spectrum for a structured IJK mesh

Usage
    Example of function object specification:
    energySpectrum1
    {
        type        energySpectrum;
        libs        ("libfieldFunctionObjects.so");
    }

    Where the entries comprise:
    \table
        Property     | Description               | Required    | Default value
        type         | type name: energySpectrum | yes         |
        log          | write info to standard output | no      | yes
    \endtable

    Output data is written to the file \<timeDir\>/energySpectrum.dat

etc/caseDicts/postProcessing/fields/energySpectrum

@@ -0,0 +1,19 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Web:      www.OpenFOAM.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+Description
+    Calculates the energy spectrum for a box of turbulence
+
+\*---------------------------------------------------------------------------*/
+
+type            energySpectrum;
+libs            ("librandomProcessesFunctionObjects.so");
+
+executeControl  writeTime;
+writeControl    writeTime;
+
+// ************************************************************************* //

src/functionObjects/randomProcesses/Allwmake

@@ -0,0 +1,16 @@
+#!/bin/sh
+cd ${0%/*} || exit 1                        # Run from this directory
+targetType=libso                            # Preferred library type
+. $WM_PROJECT_DIR/wmake/scripts/AllwmakeParseArguments
+. $WM_PROJECT_DIR/wmake/scripts/have_fftw
+
+#------------------------------------------------------------------------------
+
+if have_fftw
+then
+    wmake $targetType
+else
+    echo "==> skip randomProcesses library (no FFTW)"
+fi
+
+#------------------------------------------------------------------------------

src/functionObjects/randomProcesses/Make/files

@@ -0,0 +1,3 @@
+energySpectrum/energySpectrum.C
+
+LIB = $(FOAM_LIBBIN)/librandomProcessesFunctionObjects

src/functionObjects/randomProcesses/Make/options

@@ -0,0 +1,9 @@
+EXE_INC = \
+    -I$(LIB_SRC)/randomProcesses/lnInclude \
+    -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/meshTools/lnInclude
+
+LIB_LIBS = \
+    -lrandomProcesses \
+    -lfiniteVolume \
+    -lmeshTools

src/functionObjects/randomProcesses/energySpectrum/energySpectrum.C

@@ -0,0 +1,263 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2018 OpenCFD Ltd.
+     \\/     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 "energySpectrum.H"
+#include "fft.H"
+#include "fvMesh.H"
+#include "boundBox.H"
+#include "OFstream.H"
+#include "mathematicalConstants.H"
+#include "volFields.H"
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace functionObjects
+{
+    defineTypeNameAndDebug(energySpectrum, 0);
+    addToRunTimeSelectionTable(functionObject, energySpectrum, dictionary);
+}
+}
+
+
+// * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * //
+
+void Foam::functionObjects::energySpectrum::writeFileHeader(Ostream& os)
+{
+    writeHeader(os, "Turbulence energy spectra");
+
+    writeCommented(os, "kappa E(kappa)");
+
+    os  << endl;
+}
+
+
+void Foam::functionObjects::energySpectrum::calcAndWriteSpectrum
+(
+    const vectorField& U,
+    const vectorField& C,
+    const vector& c0,
+    const vector& deltaC,
+    const Vector<label>& N,
+    const scalar kappaNorm
+)
+{
+    Log << "Computing FFT" << endl;
+    const complexVectorField Uf
+    (
+        fft::forwardTransform
+        (
+            ReComplexField(U),
+            List<label>({N.x(), N.y(), N.z()})
+        )
+       /scalar(cmptProduct(N))
+    );
+
+
+    Log << "Computing wave numbers" << endl;
+    const label nMax = cmptMax(N);
+    scalarField kappa(nMax);
+    forAll(kappa, kappai)
+    {
+        kappa[kappai] = kappai*kappaNorm;
+    }
+
+    Log << "Computing energy spectrum" << endl;
+    scalarField E(nMax, 0);
+    const scalarField Ec(0.5*magSqr(Uf));
+    forAll(C, celli)
+    {
+        point Cc(C[celli] - c0);
+        label i = round((Cc.x() - 0.5*deltaC.x())/(deltaC.x())*(N.x() - 1));
+        label j = round((Cc.y() - 0.5*deltaC.y())/(deltaC.y())*(N.y() - 1));
+        label k = round((Cc.z() - 0.5*deltaC.z())/(deltaC.z())*(N.z() - 1));
+        label kappai = round(Foam::sqrt(scalar(i*i + j*j + k*k)));
+
+        E[kappai] += Ec[celli];
+    }
+
+    E /= kappaNorm;
+
+    Log << "Writing spectrum" << endl;
+    autoPtr<OFstream> osPtr = createFile(name(), time_.value());
+    OFstream& os = osPtr.ref();
+    writeFileHeader(os);
+
+    forAll(kappa, kappai)
+    {
+        os  << kappa[kappai] << tab << E[kappai] << nl;
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::functionObjects::energySpectrum::energySpectrum
+(
+    const word& name,
+    const Time& runTime,
+    const dictionary& dict
+)
+:
+    fvMeshFunctionObject(name, runTime, dict),
+    writeFile(mesh_, name),
+    cellAddr_(mesh_.nCells()),
+    UName_("U"),
+    N_(Zero),
+    c0_(Zero),
+    deltaC_(Zero),
+    kappaNorm_(0)
+{
+    read(dict);
+}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+bool Foam::functionObjects::energySpectrum::read(const dictionary& dict)
+{
+    fvMeshFunctionObject::read(dict);
+    writeFile::read(dict);
+
+    dict.readIfPresent("U", UName_);
+
+    const boundBox meshBb(mesh_.bounds());
+
+    // Assume all cells are the same size...
+    const cell& c = mesh_.cells()[0];
+    boundBox cellBb(boundBox::invertedBox);
+    forAll(c, facei)
+    {
+        const face& f = mesh_.faces()[c[facei]];
+        cellBb.add(mesh_.points(), f);
+    }
+
+    const vector L(meshBb.max() - meshBb.min());
+    const vector nCellXYZ(cmptDivide(L, cellBb.max() - cellBb.min()));
+
+    N_ = Vector<label>
+    (
+        round(nCellXYZ.x()),
+        round(nCellXYZ.z()),
+        round(nCellXYZ.z())
+    );
+
+    // Check that the mesh is a structured box
+    vector cellDx(cellBb.max() - cellBb.min());
+    vector expectedMax(N_.x()*cellDx.x(), N_.y()*cellDx.y(), N_.z()*cellDx.z());
+    vector relativeSize(cmptDivide(L, expectedMax));
+    for (direction i = 0; i < 3; ++i)
+    {
+        if (mag(relativeSize[i] - 1) > 1e-3)
+        {
+            FatalErrorInFunction
+                << name() << " function object is only appropriate for "
+                << "isotropic structured IJK meshes. Mesh extents: " << L
+                << ", computed IJK mesh extents: " << expectedMax
+                << exit(FatalError);
+        }
+    }
+    Log << "Mesh extents (deltax,deltay,deltaz): " << L << endl;
+    Log << "Number of cells (Nx,Ny,Nz): " << N_ << endl;
+
+    // Map into i-j-k co-ordinates
+    const vectorField& C = mesh_.C();
+    c0_ = returnReduce(min(C), minOp<vector>());
+    const vector cMax = returnReduce(max(C), maxOp<vector>());
+    deltaC_ = cMax - c0_;
+
+    forAll(C, celli)
+    {
+        label i = round((C[celli].x() - c0_.x())/(deltaC_.x())*(N_.x() - 1));
+        label j = round((C[celli].y() - c0_.y())/(deltaC_.y())*(N_.y() - 1));
+        label k = round((C[celli].z() - c0_.z())/(deltaC_.z())*(N_.z() - 1));
+
+        cellAddr_[celli] = k + j*N_.y() + i*N_.y()*N_.z();
+    }
+
+    kappaNorm_ = constant::mathematical::twoPi/cmptMax(L);
+
+    return true;
+}
+
+
+bool Foam::functionObjects::energySpectrum::execute()
+{
+    return true;
+}
+
+
+bool Foam::functionObjects::energySpectrum::write()
+{
+    const auto& U = mesh_.lookupObject<volVectorField>(UName_);
+    const vectorField& Uc = U.primitiveField();
+    const vectorField& C = mesh_.C();
+
+    if (Pstream::parRun())
+    {
+        PstreamBuffers pBufs(Pstream::commsTypes::nonBlocking);
+
+        UOPstream toProc(0, pBufs);
+
+        toProc << Uc << C << cellAddr_;
+
+        pBufs.finishedSends();
+
+        if (Pstream::master())
+        {
+            const label nGlobalCells(cmptProduct(N_));
+            vectorField Uijk(nGlobalCells);
+            vectorField Cijk(nGlobalCells);
+
+            for (label proci = 0; proci < Pstream::nProcs(); ++proci)
+            {
+                UIPstream fromProc(proci, pBufs);
+                vectorField Up;
+                vectorField Cp;
+                labelList cellAddrp;
+                fromProc >> Up >> Cp >> cellAddrp;
+                UIndirectList<vector>(Uijk, cellAddrp) = Up;
+                UIndirectList<vector>(Cijk, cellAddrp) = Cp;
+            }
+
+            calcAndWriteSpectrum(Uijk, Cijk, c0_, deltaC_, N_, kappaNorm_);
+        }
+
+    }
+    else
+    {
+        vectorField Uijk(Uc, cellAddr_);
+        vectorField Cijk(C, cellAddr_);
+
+        calcAndWriteSpectrum(Uijk, Cijk, c0_, deltaC_, N_, kappaNorm_);
+    }
+
+    return true;
+}
+
+
+// ************************************************************************* //

src/functionObjects/randomProcesses/energySpectrum/energySpectrum.H

@@ -0,0 +1,178 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2018 OpenCFD Ltd.
+     \\/     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::functionObjects::energySpectrum
+
+Group
+    grpFieldFunctionObjects
+
+Description
+    Calculates the energy spectrum for a structured IJK mesh
+
+Usage
+    Example of function object specification:
+    \verbatim
+    energySpectrum1
+    {
+        type        energySpectrum;
+        libs        ("libfieldFunctionObjects.so");
+        ...
+        log         yes;
+    }
+    \endverbatim
+
+    Where the entries comprise:
+    \table
+        Property     | Description               | Required    | Default value
+        type         | type name: energySpectrum | yes         |
+        log          | write info to standard output | no      | yes
+    \endtable
+
+    Output data is written to the file \<timeDir\>/energySpectrum.dat
+
+See also
+    Foam::functionObjects::fvMeshFunctionObject
+    Foam::functionObjects::writeFile
+
+SourceFiles
+    energySpectrum.C
+    energySpectrumTemplates.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef functionObjects_energySpectrum_H
+#define functionObjects_energySpectrum_H
+
+#include "fvMeshFunctionObject.H"
+#include "writeFile.H"
+#include "Vector.H"
+#include "vectorField.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace functionObjects
+{
+
+/*---------------------------------------------------------------------------*\
+                       Class energySpectrum Declaration
+\*---------------------------------------------------------------------------*/
+
+class energySpectrum
+:
+    public fvMeshFunctionObject,
+    public writeFile
+{
+protected:
+
+    // Protected data
+
+        //- I-J-K mesh addressing
+        labelList cellAddr_;
+
+        //- Name of velocity field, default = U
+        word UName_;
+
+        //- Number of cells in I-J-K directions
+        Vector<label> N_;
+
+        //- Reference point
+        vector c0_;
+
+        //- Cell length scale
+        vector deltaC_;
+
+        //- Wave number
+        scalar kappaNorm_;
+
+
+    // Protected Member Functions
+
+
+        //- Output file header information
+        virtual void writeFileHeader(Ostream& os);
+
+        //- Calculate and write the spectrum
+        void calcAndWriteSpectrum
+        (
+            const vectorField& U,
+            const vectorField& C,
+            const vector& c0,
+            const vector& deltaC,
+            const Vector<label>& N,
+            const scalar kappaNorm
+        );
+
+        //- No copy construct
+        energySpectrum(const energySpectrum&) = delete;
+
+        //- No copy assignment
+        void operator=(const energySpectrum&) = delete;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("energySpectrum");
+
+
+    // Constructors
+
+        //- Construct from Time and dictionary
+        energySpectrum
+        (
+            const word& name,
+            const Time& runTime,
+            const dictionary& dict
+        );
+
+
+    //- Destructor
+    virtual ~energySpectrum() = default;
+
+
+    // Member Functions
+
+        //- Read the field min/max data
+        virtual bool read(const dictionary&);
+
+        //- Execute, currently does nothing
+        virtual bool execute();
+
+        //- Write the energySpectrum
+        virtual bool write();
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace functionObjects
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //