Merge branch 'feature-noise-fft' into 'develop'

Updated noise fft handling via fftw

See merge request Development/OpenFOAM-plus!186

src/randomProcesses/fft/fft.C

@@ -34,6 +34,60 @@ namespace Foam
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
+Foam::tmp<Foam::complexField> fft::realTransform1D(const scalarField& field)
+{
+    const label n = field.size();
+    const label nBy2 = n/2;
+
+    // Copy of input field for use by fftw
+    // - fftw requires non-const access to input and output...
+    scalar in[n], out[n];
+    forAll(field, i)
+    {
+        in[i] = field[i];
+    }
+
+    // Using real to half-complex fftw 'kind'
+    fftw_plan plan = fftw_plan_r2r_1d
+    (
+        n,
+        in,
+        out,
+        FFTW_R2HC,
+        FFTW_ESTIMATE
+    );
+
+    fftw_execute(plan);
+
+    // field[0] = DC component
+    tmp<complexField> tresult(new complexField(nBy2 + 1));
+    complexField& result = tresult.ref();
+
+    result[0].Re() = out[0];
+    result[nBy2].Re() = out[nBy2];
+    for (label i = 1; i < nBy2; ++i)
+    {
+        result[i].Re() = out[i];
+        result[i].Im() = out[n - i];
+    }
+
+    fftw_destroy_plan(plan);
+
+    return tresult;
+}
+
+
+Foam::tmp<Foam::complexField> fft::realTransform1D
+(
+    const tmp<scalarField>& tfield
+)
+{
+    tmp<complexField> tresult = realTransform1D(tfield());
+    tfield.clear();
+    return tresult;
+}
+
+
 void fft::transform
 (
     complexField& field,

src/randomProcesses/fft/fft.H

@@ -44,7 +44,6 @@ SourceFiles
 #define fft_H
 
 #include "complexFields.H"
-#include "UList.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -63,6 +62,19 @@ public:
     };
 
 
+    //- Transform real-value data
+    //  - uses the fftw real to half-complex method
+    //  - result size is field.size()/2 + 1
+    static tmp<complexField> realTransform1D(const scalarField& field);
+
+
+    //- Transform real-value data
+    //  - uses the fftw real to half-complex method
+    //  - result size is field.size()/2 + 1
+    static tmp<complexField> realTransform1D(const tmp<scalarField>& field);
+
+
+    //- Transform compex-value data
     static void transform
     (
         complexField& field,

src/randomProcesses/noise/noiseFFT/noiseFFT.C

@@ -3,7 +3,7 @@
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
     \\  /    A nd           | Copyright (C) 2011-2015 OpenFOAM Foundation
-     \\/     M anipulation  | Copyright (C) 2016 OpenCFD Ltd.
+     \\/     M anipulation  | Copyright (C) 2016-2017 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -26,7 +26,6 @@ License
 #include "noiseFFT.H"
 #include "IFstream.H"
 #include "DynamicList.H"
-#include "SubField.H"
 #include "mathematicalConstants.H"
 #include "HashSet.H"
 #include "fft.H"
@@ -139,24 +138,61 @@ void Foam::noiseFFT::octaveBandInfo
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
-Foam::noiseFFT::noiseFFT
+Foam::noiseFFT::noiseFFT(const scalar deltaT, const label windowSize)
-(
-    const scalar deltaT,
-    const scalarField& pressure
-)
 :
-    scalarField(pressure),
+    scalarField(),
     deltaT_(deltaT)
 {
+    if (windowSize > 1)
+    {
+        planInfo_.active = true;
+        planInfo_.windowSize = windowSize;
+        planInfo_.in.setSize(windowSize);
+        planInfo_.out.setSize(windowSize);
+
+        // Using real to half-complex fftw 'kind'
+        planInfo_.plan =
+            fftw_plan_r2r_1d
+            (
+                windowSize,
+                planInfo_.in.data(),
+                planInfo_.out.data(),
+                FFTW_R2HC,
+                windowSize <= 8192 ? FFTW_MEASURE : FFTW_ESTIMATE
+            );
+    }
+    else
+    {
+        planInfo_.active = false;
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::noiseFFT::~noiseFFT()
+{
+    if (planInfo_.active)
+    {
+        planInfo_.active = false;
+        fftw_destroy_plan(planInfo_.plan);
+        fftw_cleanup();
+    }
+}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+void Foam::noiseFFT::setData(scalarList& data)
+{
+    this->transfer(data);
+
     scalarField& p = *this;
     p -= average(p);
 }
 
 
-Foam::noiseFFT::noiseFFT(const fileName& pFileName, const label skip)
+void Foam::noiseFFT::setData(const fileName& pFileName, const label skip)
-:
-    scalarField(),
-    deltaT_(0.0)
 {
     // Construct pressure data file
     IFstream pFile(pFileName);
@@ -173,7 +209,7 @@ Foam::noiseFFT::noiseFFT(const fileName& pFileName, const label skip)
     {
         scalar dummyt, dummyp;
 
-        for (label i = 0; i < skip; i++)
+        for (label i = 0; i < skip; ++i)
         {
             pFile >> dummyt;
 
@@ -190,7 +226,7 @@ Foam::noiseFFT::noiseFFT(const fileName& pFileName, const label skip)
     }
 
     scalar t = 0, T0 = 0, T1 = 0;
-    DynamicList<scalar> pData(100000);
+    DynamicList<scalar> pData(1024);
     label i = 0;
 
     while (!(pFile >> t).eof())
@@ -214,8 +250,6 @@ Foam::noiseFFT::noiseFFT(const fileName& pFileName, const label skip)
 }
 
 
-// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
-
 Foam::graph Foam::noiseFFT::pt() const
 {
     scalarField t(size());
@@ -240,33 +274,48 @@ Foam::tmp<Foam::scalarField> Foam::noiseFFT::Pf
     const tmp<scalarField>& tpn
 ) const
 {
-    // Calculate the 2-sided fft
+    if (planInfo_.active)
-    // Note: result not scaled
+    {
-    tmp<scalarField> tPn2
+        const scalarField& pn = tpn();
-    (
+        if (pn.size() != planInfo_.windowSize)
-        mag
+        {
-        (
+            FatalErrorInFunction
-            fft::reverseTransform
+                << "Expected pressure data to have " << planInfo_.windowSize
-            (
+                << " values, but received " << pn.size() << " values"
-                ReComplexField(tpn),
+                << abort(FatalError);
-                List<int>(1, tpn().size())
+        }
-            )
-        )
-    );
 
+        List<scalar>& in = planInfo_.in;
+        const List<scalar>& out = planInfo_.out;
+        forAll(in, i)
+        {
+            in[i] = pn[i];
+        }
         tpn.clear();
 
-    // Only storing the positive half
+        ::fftw_execute(planInfo_.plan);
-    // Note: storing (N/2) values, DC component at position (0)
-    tmp<scalarField> tPn
-    (
-        new scalarField
-        (
-            scalarField::subField(tPn2(), tPn2().size()/2 + 1)
-        )
-    );
 
-    return tPn;
+        const label n = planInfo_.windowSize;
+        const label nBy2 = n/2;
+        tmp<scalarField> tresult(new scalarField(nBy2 + 1));
+        scalarField& result = tresult.ref();
+
+        // 0 th value = DC component
+        // nBy2 th value = real only if n is even
+        result[0] = out[0];
+        for (label i = 1; i <= nBy2; ++i)
+        {
+            scalar re = out[i];
+            scalar im = out[n - i];
+            result[i] = sqrt(re*re + im*im);
+        }
+
+        return tresult;
+    }
+    else
+    {
+        return mag(fft::realTransform1D(tpn));
+    }
 }
 
 
@@ -421,7 +470,7 @@ Foam::graph Foam::noiseFFT::octaves
     scalarField octData(freqBandIDs.size() - 1, 0.0);
     scalarField fm(freqBandIDs.size() -1, 0.0);
 
-    for (label bandI = 0; bandI < freqBandIDs.size() - 1; bandI++)
+    for (label bandI = 0; bandI < freqBandIDs.size() - 1; ++bandI)
     {
         label fb0 = freqBandIDs[bandI];
         label fb1 = freqBandIDs[bandI+1];
@@ -431,7 +480,7 @@ Foam::graph Foam::noiseFFT::octaves
 
         if (integrate)
         {
-            for (label freqI = fb0; freqI < fb1; freqI++)
+            for (label freqI = fb0; freqI < fb1; ++freqI)
             {
                 label f0 = f[freqI];
                 label f1 = f[freqI + 1];
@@ -441,7 +490,7 @@ Foam::graph Foam::noiseFFT::octaves
         }
         else
         {
-            for (label freqI = fb0; freqI < fb1; freqI++)
+            for (label freqI = fb0; freqI < fb1; ++freqI)
             {
                 octData[bandI] += data[freqI];
             }

src/randomProcesses/noise/noiseFFT/noiseFFT.H

@@ -3,7 +3,7 @@
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
     \\  /    A nd           | Copyright (C) 2011-2015 OpenFOAM Foundation
-     \\/     M anipulation  | Copyright (C) 2016 OpenCFD Ltd.
+     \\/     M anipulation  | Copyright (C) 2016-2017 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -53,6 +53,7 @@ SeeAlso
 #include "scalarField.H"
 #include "graph.H"
 #include "windowModel.H"
+#include <fftw3.h>
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -67,11 +68,17 @@ class noiseFFT
 :
     public scalarField
 {
-    // Private data
+    //- FFTW planner information
-
+    struct planInfo
-        //- Time spacing of the raw data
+    {
-        scalar deltaT_;
+        bool active;
+        label windowSize;
+        scalarList in;
+        scalarList out;
+        fftw_plan plan;
+    };
 
+    //- Octave band information
     struct octaveBandInfo
     {
         label octave;
@@ -84,26 +91,25 @@ class noiseFFT
     };
 
 
+    // Private data
+
+        //- Time spacing of the raw data (uniform)
+        scalar deltaT_;
+
+        //- Plan information for FFTW
+        mutable planInfo planInfo_;
+
+
 public:
 
     //- Reference pressure
     static scalar p0;
 
-
-    // Constructors
-
     //- Construct from pressure field
-        noiseFFT
+    noiseFFT(const scalar deltaT, const label windowSize = -1);
-        (
-            const scalar deltat,
-            const scalarField& pressure
-        );
 
-        //- Construct from Istream
+    //- Destructor
-        noiseFFT(Istream&);
+    ~noiseFFT();
-
-        //- Construct from pressure field file name
-        noiseFFT(const fileName& pFileName, const label skip = 0);
 
 
     // Member Functions
@@ -135,6 +141,14 @@ public:
             scalarField& fCentre
         );
 
+        //- Set the pressure data
+        //- \note transfers storage to *this
+        void setData(scalarList& data);
+
+        //- Set the pressure data by reading from file with an optional offset
+        void setData(const fileName& pFileName, const label skip = 0);
+
+
         //- Return the graph of pressure as a function of time
         graph pt() const;
 

src/randomProcesses/noise/noiseModels/pointNoise/pointNoise.C

@@ -101,7 +101,8 @@ void pointNoise::processData
     fileName outDir(baseFileDir(dataseti)/fNameBase);
 
     // Create the fft
-    noiseFFT nfft(deltaT, p);
+    noiseFFT nfft(deltaT, win.nSamples());
+    nfft.setData(p);
 
 
     // Narrow band data

src/randomProcesses/noise/noiseModels/surfaceNoise/surfaceNoise.C

@@ -142,9 +142,9 @@ void surfaceNoise::readSurfaceData
         // Complete pressure time history data for subset of faces
         pData.setSize(nLocalFace);
         const label nTimes = times_.size();
-        forAll(pData, faceI)
+        for (scalarField& pf : pData)
         {
-            pData[faceI].setSize(nTimes);
+            pf.setSize(nTimes);
         }
 
         // Read and send data
@@ -165,7 +165,7 @@ void surfaceNoise::readSurfaceData
                 p *= rhoRef_;
 
                 // Send subset of faces to each processor
-                for (label procI = 0; procI < Pstream::nProcs(); procI++)
+                for (label procI = 0; procI < Pstream::nProcs(); ++procI)
                 {
                     label face0 = procFaceOffset[procI];
                     label nLocalFace =
@@ -194,12 +194,9 @@ void surfaceNoise::readSurfaceData
         const label nLocalFace = procFaceOffset[0];
 
         pData.setSize(nLocalFace);
-        forAll(times_, timeI)
+        for (scalarField& pf : pData)
         {
-            forAll(pData, faceI)
+            pf.setSize(times_.size());
-            {
-                pData[faceI].setSize(times_.size());
-            }
         }
 
         forAll(times_, i)
@@ -266,7 +263,7 @@ Foam::scalar surfaceNoise::writeSurfaceData
                 allData[faceI] = data[faceI];
             }
 
-            for (label procI = 1; procI < Pstream::nProcs(); procI++)
+            for (label procI = 1; procI < Pstream::nProcs(); ++procI)
             {
                 UIPstream fromProc(procI, pBufs);
                 scalarList dataSlice(fromProc);
@@ -486,7 +483,7 @@ void surfaceNoise::calculate()
             const label nFacePerProc = floor(nFace_/nProcs) + 1;
 
             procFaceOffset.setSize(nProcs + 1, 0);
-            for (label i = 1; i < procFaceOffset.size(); i++)
+            for (label i = 1; i < procFaceOffset.size(); ++i)
             {
                 procFaceOffset[i] = min(i*nFacePerProc, nFace_);
             }
@@ -539,8 +536,8 @@ void surfaceNoise::calculate()
         if (octave13BandIDs.empty())
         {
             WarningInFunction
-                << "Ocatve band calculation failed (zero sized). "
+                << "Octave band calculation failed (zero sized). "
-                << "please check your input data"
+                << "Please check your input data"
                 << endl;
         }
         else
@@ -558,11 +555,15 @@ void surfaceNoise::calculate()
 
         const windowModel& win = windowModelPtr_();
 
+        {
+            noiseFFT nfft(deltaT_, win.nSamples());
+
             forAll(pData, faceI)
             {
-            const scalarField& p = pData[faceI];
+                // Note: passes storage from pData to noiseFFT
+                nfft.setData(pData[faceI]);
 
-            noiseFFT nfft(deltaT_, p);
+                // Generate the FFT-based data
                 graph Prmsf(nfft.RMSmeanPf(win));
                 graph PSDf(nfft.PSDf(win));
 
@@ -587,6 +588,7 @@ void surfaceNoise::calculate()
                     surfPrms13f2[freqI][faceI] = Prms13f2.y()[freqI];
                 }
             }
+        }
 
         const word fNameBase = fName.nameLessExt();
 
@@ -609,7 +611,7 @@ void surfaceNoise::calculate()
             fileName outDir(outDirBase/"fft");
 
             // Determine frequency range of interest
-            // Note: freqencies have fixed interval, and are in the range
+            // Note: frequencies have fixed interval, and are in the range
             //       0 to fftWriteInterval_*(n-1)*deltaf
             label f0 = ceil(fLower_/deltaf/scalar(fftWriteInterval_));
             label f1 = floor(fUpper_/deltaf/scalar(fftWriteInterval_));

src/randomProcesses/windowModels/windowModel/windowModelTemplates.C

@@ -23,6 +23,10 @@ License
 
 \*---------------------------------------------------------------------------*/
 
+#include "SubField.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
 template<class Type>
 Foam::tmp<Foam::Field<Type>> Foam::windowModel::apply
 (