Skip to content
GitLab
Commit 1428310b authored by Andrew Heather's avatar Andrew Heather
Browse files

ENH: FFT - only create the fftw plan once

parent b2f002ba
Hide whitespace changes
Inline Side-by-side
src/randomProcesses/fft/fft.C
View file @ 1428310b
......@@ -36,39 +36,43 @@ 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...
scalarField fftInOut(field);
const label n = fftInOut.size();
const label nBy2 = n/2;
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,
fftInOut.data(),
fftInOut.data(),
in,
out,
FFTW_R2HC,
FFTW_ESTIMATE
);
fftw_execute(plan);
fftw_destroy_plan(plan);
// field[0] = DC component
tmp<complexField> tresult(new complexField(nBy2 + 1));
complexField& result = tresult.ref();
result[0].Re() = fftInOut[0];
result[nBy2].Re() = fftInOut[nBy2];
result[0].Re() = out[0];
result[nBy2].Re() = out[nBy2];
for (label i = 1; i < nBy2; ++i)
{
result[i].Re() = fftInOut[i];
result[i].Im() = fftInOut[n - i];
result[i].Re() = out[i];
result[i].Im() = out[n - i];
}
fftw_destroy_plan(plan);
return tresult;
}
......
src/randomProcesses/noise/noiseFFT/noiseFFT.C
View file @ 1428310b
......@@ -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.
......@@ -138,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);
......@@ -172,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;
......@@ -189,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())
......@@ -213,8 +250,6 @@ Foam::noiseFFT::noiseFFT(const fileName& pFileName, const label skip)
}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
Foam::graph Foam::noiseFFT::pt() const
{
scalarField t(size());
......@@ -239,7 +274,48 @@ Foam::tmp<Foam::scalarField> Foam::noiseFFT::Pf
const tmp<scalarField>& tpn
) const
{
return mag(fft::realTransform1D(tpn));
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);
}
List<scalar>& in = planInfo_.in;
const List<scalar>& out = planInfo_.out;
forAll(in, i)
{
in[i] = pn[i];
}
tpn.clear();
::fftw_execute(planInfo_.plan);
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));
}
}
......@@ -394,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];
......@@ -404,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];
......@@ -414,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
View file @ 1428310b
......@@ -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
//- 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
};
public:
// Private data
//- Reference pressure
static scalar p0;
//- Time spacing of the raw data (uniform)
scalar deltaT_;
//- Plan information for FFTW
mutable planInfo planInfo_;
// Constructors
//- Construct from pressure field
noiseFFT
(
const scalar deltat,
const scalarField& pressure
);
public:
//- Construct from Istream
noiseFFT(Istream&);
//- Reference pressure
static scalar p0;
//- Construct from pressure field file name
noiseFFT(const fileName& pFileName, const label skip = 0);
//- Construct from pressure field
noiseFFT(const scalar deltaT, const label windowSize = -1);
//- 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;
......
src/randomProcesses/noise/noiseModels/pointNoise/pointNoise.C
View file @ 1428310b
......@@ -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
View file @ 1428310b
......@@ -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_, p);
graph Prmsf(nfft.RMSmeanPf(win));
graph PSDf(nfft.PSDf(win));
noiseFFT nfft(deltaT_, win.nSamples());
// 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_));
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment