diff --git a/src/randomProcesses/fft/fft.C b/src/randomProcesses/fft/fft.C index fd44c0109d..1f37590be9 100644 --- a/src/randomProcesses/fft/fft.C +++ b/src/randomProcesses/fft/fft.C @@ -34,6 +34,60 @@ namespace Foam // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // +Foam::tmp 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 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 fft::realTransform1D +( + const tmp& tfield +) +{ + tmp tresult = realTransform1D(tfield()); + tfield.clear(); + return tresult; +} + + void fft::transform ( complexField& field, diff --git a/src/randomProcesses/fft/fft.H b/src/randomProcesses/fft/fft.H index 11cc07575e..7b62e44757 100644 --- a/src/randomProcesses/fft/fft.H +++ b/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 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 realTransform1D(const tmp& field); + + + //- Transform compex-value data static void transform ( complexField& field, diff --git a/src/randomProcesses/noise/noiseFFT/noiseFFT.C b/src/randomProcesses/noise/noiseFFT/noiseFFT.C index e82d25b90f..2b3a3f5376 100644 --- a/src/randomProcesses/noise/noiseFFT/noiseFFT.C +++ b/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 -( - const scalar deltaT, - const scalarField& pressure -) +Foam::noiseFFT::noiseFFT(const scalar deltaT, const label windowSize) : - 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) -: - scalarField(), - deltaT_(0.0) +void Foam::noiseFFT::setData(const fileName& pFileName, const label skip) { // 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 pData(100000); + DynamicList 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::noiseFFT::Pf const tmp& tpn ) const { - // Calculate the 2-sided fft - // Note: result not scaled - tmp tPn2 - ( - mag - ( - fft::reverseTransform - ( - ReComplexField(tpn), - List(1, tpn().size()) - ) - ) - ); + if (planInfo_.active) + { + const scalarField& pn = tpn(); + if (pn.size() != planInfo_.windowSize) + { + FatalErrorInFunction + << "Expected pressure data to have " << planInfo_.windowSize + << " values, but received " << pn.size() << " values" + << abort(FatalError); + } - tpn.clear(); + List& in = planInfo_.in; + const List& out = planInfo_.out; + forAll(in, i) + { + in[i] = pn[i]; + } + tpn.clear(); - // Only storing the positive half - // Note: storing (N/2) values, DC component at position (0) - tmp tPn - ( - new scalarField - ( - scalarField::subField(tPn2(), tPn2().size()/2 + 1) - ) - ); + ::fftw_execute(planInfo_.plan); - return tPn; + const label n = planInfo_.windowSize; + const label nBy2 = n/2; + tmp 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]; } diff --git a/src/randomProcesses/noise/noiseFFT/noiseFFT.H b/src/randomProcesses/noise/noiseFFT/noiseFFT.H index 6439a2f429..cca1ee25d1 100644 --- a/src/randomProcesses/noise/noiseFFT/noiseFFT.H +++ b/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 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // @@ -67,11 +68,17 @@ class noiseFFT : public scalarField { - // Private data - - //- Time spacing of the raw data - scalar deltaT_; + //- FFTW planner information + struct planInfo + { + 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; + //- Construct from pressure field + noiseFFT(const scalar deltaT, const label windowSize = -1); - // Constructors - - //- Construct from pressure field - noiseFFT - ( - const scalar deltat, - const scalarField& pressure - ); - - //- Construct from Istream - noiseFFT(Istream&); - - //- Construct from pressure field file name - noiseFFT(const fileName& pFileName, const label skip = 0); + //- Destructor + ~noiseFFT(); // 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; diff --git a/src/randomProcesses/noise/noiseModels/pointNoise/pointNoise.C b/src/randomProcesses/noise/noiseModels/pointNoise/pointNoise.C index 209febd8c9..3a58db4cd2 100644 --- a/src/randomProcesses/noise/noiseModels/pointNoise/pointNoise.C +++ b/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 diff --git a/src/randomProcesses/noise/noiseModels/surfaceNoise/surfaceNoise.C b/src/randomProcesses/noise/noiseModels/surfaceNoise/surfaceNoise.C index 367459ce47..b76f4e8329 100644 --- a/src/randomProcesses/noise/noiseModels/surfaceNoise/surfaceNoise.C +++ b/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) - { - pData[faceI].setSize(times_.size()); - } + pf.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). " - << "please check your input data" + << "Octave band calculation failed (zero sized). " + << "Please check your input data" << endl; } else @@ -558,33 +555,38 @@ void surfaceNoise::calculate() const windowModel& win = windowModelPtr_(); - forAll(pData, faceI) { - const scalarField& p = pData[faceI]; + noiseFFT nfft(deltaT_, win.nSamples()); - noiseFFT nfft(deltaT_, p); - graph Prmsf(nfft.RMSmeanPf(win)); - graph PSDf(nfft.PSDf(win)); - - // Store the frequency results in slot for face of surface - forAll(surfPrmsf, i) + forAll(pData, faceI) { - label freqI = i*fftWriteInterval_; - surfPrmsf[i][faceI] = Prmsf.y()[freqI]; - surfPSDf[i][faceI] = PSDf.y()[freqI]; - } + // Note: passes storage from pData to noiseFFT + nfft.setData(pData[faceI]); - // PSD [Pa^2/Hz] - graph PSD13f(nfft.octaves(PSDf, octave13BandIDs, false)); + // Generate the FFT-based data + graph Prmsf(nfft.RMSmeanPf(win)); + graph PSDf(nfft.PSDf(win)); - // Integrated PSD = P(rms)^2 [Pa^2] - graph Prms13f2(nfft.octaves(PSDf, octave13BandIDs, true)); + // Store the frequency results in slot for face of surface + forAll(surfPrmsf, i) + { + label freqI = i*fftWriteInterval_; + surfPrmsf[i][faceI] = Prmsf.y()[freqI]; + surfPSDf[i][faceI] = PSDf.y()[freqI]; + } - // Store the 1/3 octave results in slot for face of surface - forAll(surfPSD13f, freqI) - { - surfPSD13f[freqI][faceI] = PSD13f.y()[freqI]; - surfPrms13f2[freqI][faceI] = Prms13f2.y()[freqI]; + // PSD [Pa^2/Hz] + graph PSD13f(nfft.octaves(PSDf, octave13BandIDs, false)); + + // Integrated PSD = P(rms)^2 [Pa^2] + graph Prms13f2(nfft.octaves(PSDf, octave13BandIDs, true)); + + // Store the 1/3 octave results in slot for face of surface + forAll(surfPSD13f, freqI) + { + surfPSD13f[freqI][faceI] = PSD13f.y()[freqI]; + surfPrms13f2[freqI][faceI] = Prms13f2.y()[freqI]; + } } } @@ -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_)); diff --git a/src/randomProcesses/windowModels/windowModel/windowModelTemplates.C b/src/randomProcesses/windowModels/windowModel/windowModelTemplates.C index 9d3e478c04..90f9971381 100644 --- a/src/randomProcesses/windowModels/windowModel/windowModelTemplates.C +++ b/src/randomProcesses/windowModels/windowModel/windowModelTemplates.C @@ -23,6 +23,10 @@ License \*---------------------------------------------------------------------------*/ +#include "SubField.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + template Foam::tmp> Foam::windowModel::apply (