Commit cc36db19 authored by Andrew Heather's avatar Andrew Heather
Browse files

GIT: Resolved conflict

parents dd3fb622 cabd03d0
EXE_INC = \
-I$(LIB_SRC)/randomProcesses/lnInclude \
-I$(LIB_SRC)/sampling/lnInclude
-I$(LIB_SRC)/sampling/lnInclude \
-I$(LIB_SRC)/surfMesh/lnInclude
EXE_LIBS = \
-lrandomProcesses \
......
......@@ -2,8 +2,8 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2016 OpenFOAM Foundation
\\/ M anipulation |
\\ / A nd | Copyright (C) 2011-2015 OpenFOAM Foundation
\\/ M anipulation | Copyright (C) 2016 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
......@@ -28,157 +28,104 @@ Group
grpPostProcessingUtilities
Description
Utility to perform noise analysis of pressure data using the noiseFFT
library.
Utility to perform noise analysis of pressure data.
Control settings are read from the $FOAM_CASE/system/noiseDict dictionary,
or user-specified dictionary using the -dict option. Pressure data is
read using a CSV reader:
The utility provides a light wrapper around the run-time selectable
noise model. Current options include:
- point, and
- surface noise.
\heading Usage
\heading Example usage
\verbatim
pRef 101325;
N 65536;
nw 100;
f1 25;
fU 10000;
graphFormat raw;
pressureData
noiseModel surfaceNoise; // pointNoise
surfaceNoiseCoeffs
{
fileName "pressureData"
nHeaderLine 1; // number of header lines
refColumn 0; // reference column index
componentColumns (1); // component column indices
separator " "; // optional (defaults to ",")
mergeSeparators no; // merge multiple separators
outOfBounds clamp; // optional out-of-bounds handling
interpolationScheme linear; // optional interpolation scheme
windowModel Hanning;
HanningCoeffs
{
// Window overlap percentage
overlapPercent 50;
symmetric yes;
extended yes;
// Optional number of windows, default = all available
nWindow 5;
}
// Input file
inputFile "postProcessing/faceSource1/surface/patch/patch.case";
// Surface reader
reader ensight;
// Surface writer
writer ensight;
// Collate times for ensight output - ensures geometry is only written once
writeOptions
{
ensight
{
collateTimes 1;
}
}
// Number of samples in sampling window
// Must be a power of 2, default = 2^16 (=65536)
N 4096;
// Write interval for FFT data, default = 1
fftWriteInterval 100;
}
\endverbatim
where
\table
Property | Description | Required | Default value
pRef | Reference pressure | no | 0
N | Number of samples in sampling window | no | 65536
nw | Number of sampling windows | no | 100
fl | Lower frequency band | no | 25
fU | Upper frequency band | no | 10000
graphFormat | Output graph format | no | raw
\endtable
Current graph outputs include:
- FFT of the pressure data
- narrow-band PFL (pressure-fluctuation level) spectrum
- one-third-octave-band PFL spectrum
- one-third-octave-band pressure spectrum
SeeAlso
CSV.H
noiseFFT.H
noiseModel.H
windowModel.H
\*---------------------------------------------------------------------------*/
#include "noiseFFT.H"
#include "argList.H"
#include "Time.H"
#include "functionObjectFile.H"
#include "CSV.H"
#include "noiseModel.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
using namespace Foam;
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
Foam::scalar checkUniformTimeStep(const scalarField& t)
{
// check that a uniform time step has been applied
scalar deltaT = -1.0;
if (t.size() > 1)
{
for (label i = 1; i < t.size(); i++)
{
scalar dT = t[i] - t[i-1];
if (deltaT < 0)
{
deltaT = dT;
}
if (mag(deltaT - dT) > SMALL)
{
FatalErrorInFunction
<< "Unable to process data with a variable time step"
<< exit(FatalError);
}
}
}
else
{
FatalErrorInFunction
<< "Unable to create FFT with a single value"
<< exit(FatalError);
}
return deltaT;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
int main(int argc, char *argv[])
{
argList::noParallel();
#include "addDictOption.H"
#include "setRootCase.H"
#include "createTime.H"
#include "createFields.H"
Info<< "Reading data file" << endl;
Function1Types::CSV<scalar> pData("pressure", dict, "Data");
// time history data
const scalarField t(pData.x());
// pressure data
const scalarField p(pData.y());
if (t.size() < N)
word dictName("noiseDict");
if (args.optionFound("dict"))
{
FatalErrorInFunction
<< "Block size N = " << N
<< " is larger than number of data = " << t.size()
<< exit(FatalError);
dictName = args["dict"];
}
Info<< " read " << t.size() << " values" << nl << endl;
Info<< "Creating noise FFT" << endl;
noiseFFT nfft(checkUniformTimeStep(t), p);
nfft -= pRef;
fileName baseFileName(pData.fName().lessExt());
graph Pf(nfft.RMSmeanPf(N, min(nfft.size()/N, nw)));
Info<< " Creating graph for " << Pf.title() << endl;
Pf.write(baseFileName + graph::wordify(Pf.title()), graphFormat);
graph Lf(nfft.Lf(Pf));
Info<< " Creating graph for " << Lf.title() << endl;
Lf.write(baseFileName + graph::wordify(Lf.title()), graphFormat);
graph Ldelta(nfft.Ldelta(Lf, f1, fU));
Info<< " Creating graph for " << Ldelta.title() << endl;
Ldelta.write(baseFileName + graph::wordify(Ldelta.title()), graphFormat);
graph Pdelta(nfft.Pdelta(Pf, f1, fU));
Info<< " Creating graph for " << Pdelta.title() << endl;
Pdelta.write(baseFileName + graph::wordify(Pdelta.title()), graphFormat);
IOdictionary dict
(
IOobject
(
dictName,
runTime.system(),
runTime,
IOobject::MUST_READ
)
);
autoPtr<noiseModel> model(noiseModel::New(dict));
model->calculate();
Info<< nl << "End\n" << endl;
......
......@@ -31,7 +31,7 @@
#------------------------------------------------------------------------------
setenv WM_PROJECT OpenFOAM
setenv WM_PROJECT_VERSION plus
setenv WM_PROJECT_VERSION stage
################################################################################
# USER EDITABLE PART: Changes made here may be lost with the next upgrade
......
Kmesh = Kmesh
fft = fft
processes = processes
UOprocess = $(processes)/UOprocess
turbulence = turbulence
noise = noise
Kmesh = Kmesh
$(Kmesh)/Kmesh.C
fft = fft
$(fft)/fft.C
$(fft)/fftRenumber.C
$(fft)/calcEk.C
$(fft)/kShellIntegration.C
processes = processes
UOprocess = $(processes)/UOprocess
$(UOprocess)/UOprocess.C
turbulence = turbulence
$(turbulence)/turbGen.C
$(noise)/noiseFFT.C
noise = noise
$(noise)/noiseFFT/noiseFFT.C
$(noise)/noiseModels/noiseModel/noiseModel.C
$(noise)/noiseModels/noiseModel/noiseModelNew.C
$(noise)/noiseModels/pointNoise/pointNoise.C
$(noise)/noiseModels/surfaceNoise/surfaceNoise.C
windowModels = windowModels
$(windowModels)/windowModel/windowModel.C
$(windowModels)/windowModel/windowModelNew.C
$(windowModels)/Hanning/Hanning.C
LIB = $(FOAM_LIBBIN)/librandomProcesses
EXE_INC = \
-I$(LIB_SRC)/finiteVolume/lnInclude
-I$(LIB_SRC)/finiteVolume/lnInclude \
-I$(LIB_SRC)/sampling/lnInclude \
-I$(LIB_SRC)/surfMesh/lnInclude
LIB_LIBS = \
-lfiniteVolume
-lfiniteVolume \
-lsampling \
-lsurfMesh
......@@ -2,8 +2,8 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2016 OpenFOAM Foundation
\\/ M anipulation |
\\ / A nd | Copyright (C) 2011-2015 OpenFOAM Foundation
\\/ M anipulation | Copyright (C) 2016 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
......@@ -30,28 +30,116 @@ License
#include "SubField.H"
#include "mathematicalConstants.H"
using namespace Foam::constant;
// * * * * * * * * * * * * * * Static Member Data * * * * * * * * * * * * * //
Foam::scalar Foam::noiseFFT::p0 = 2e-5;
Foam::tmp<Foam::scalarField> Foam::noiseFFT::frequencies
(
const label N,
const scalar deltaT
)
{
tmp<scalarField> tf(new scalarField(N/2, 0));
scalarField& f = tf();
scalar deltaf = 1.0/(N*deltaT);
forAll(f, i)
{
f[i] = i*deltaf;
}
return tf;
}
void Foam::noiseFFT::octaveFrequenciesIDs
(
const scalarField& f,
const scalar fLower,
const scalar fUpper,
const scalar octave,
labelList& freqBandIDs
)
{
// Set the indices of to the lower frequency bands for the input frequency
// range. Ensure that the centre frequency passes though 1000 Hz
// Low frequency bound given by:
// fLow = f0*(2^(bandI/octave/2))
// Centre frequency given by:
// fCentre = f0*(2^(bandI/octave))
scalar f0 = 1000;
scalar minFrequency = max(fLower, min(f));
// Lower frequency band limit
label band0Low = ceil(2*octave*log(minFrequency/f0)/log(2.0));
// Centre frequency band limit
//label band0Centre = ceil(octave*log(fLower/f0)/log(2.0));
scalar fLowerBand = f0*pow(2, band0Low/octave/2);
scalar fRatio = pow(2, 1.0/octave);
bool complete = false;
DynamicList<label> bandIDs(f.size());
forAll(f, i)
{
while (f[i] >= fLowerBand)
{
bandIDs.append(i);
fLowerBand *= fRatio;
if (fLowerBand > fUpper)
{
complete = true;
break;
}
}
if (complete) break;
}
freqBandIDs.transfer(bandIDs);
}
Foam::tmp<Foam::scalarField> Foam::noiseFFT::octaveFrequencies
(
const scalarField& f,
const scalar fLower,
const scalar fUpper,
const scalar octave
)
{
labelList freqBandIDs;
octaveFrequenciesIDs(f, fLower, fUpper, octave, freqBandIDs);
tmp<scalarField> tf(new scalarField(f, freqBandIDs));
return tf;
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::noiseFFT::noiseFFT
(
const scalar deltat,
const scalar deltaT,
const scalarField& pressure
)
:
scalarField(pressure),
deltat_(deltat)
deltaT_(deltaT)
{}
Foam::noiseFFT::noiseFFT(const fileName& pFileName, const label skip)
:
scalarField(),
deltat_(0.0)
deltaT_(0.0)
{
// Construct pressure data file
IFstream pFile(pFileName);
......@@ -68,7 +156,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;
......@@ -84,17 +172,23 @@ Foam::noiseFFT::noiseFFT(const fileName& pFileName, const label skip)
}
}
scalar t = 0, T = 0;
scalar t = 0, T0 = 0, T1 = 0;
DynamicList<scalar> pData(100000);
label i = 0;
while (!(pFile >> t).eof())
{
T = t;
if (i == 0)
{
T0 = t;
}
T1 = t;
pFile >> pData(i++);
}
deltat_ = T/pData.size();
// Note: assumes fixed time step
deltaT_ = (T1 - T0)/pData.size();
this->transfer(pData);
}
......@@ -107,7 +201,7 @@ Foam::graph Foam::noiseFFT::pt() const
scalarField t(size());
forAll(t, i)
{
t[i] = i*deltat_;
t[i] = i*deltaT_;
}
return graph
......@@ -121,52 +215,6 @@ Foam::graph Foam::noiseFFT::pt() const
}
Foam::tmp<Foam::scalarField> Foam::noiseFFT::window
(
const label N,
const label ni
) const
{
label windowOffset = N;
if ((N + ni*windowOffset) > size())
{
FatalErrorInFunction
<< "Requested window is outside set of data" << endl
<< "number of data = " << size() << endl
<< "size of window = " << N << endl
<< "window = " << ni
<< exit(FatalError);
}
tmp<scalarField> tpw(new scalarField(N));
scalarField& pw = tpw.ref();
label offset = ni*windowOffset;
forAll(pw, i)
{
pw[i] = operator[](i + offset);
}
return tpw;
}
Foam::tmp<Foam::scalarField> Foam::noiseFFT::Hanning(const label N) const
{
scalarField t(N);
forAll(t, i)
{
t[i] = i*deltat_;
}
scalar T = N*deltat_;
return 2*(0.5 - 0.5*cos(constant::mathematical::twoPi*t/T));
}
Foam::tmp<Foam::scalarField> Foam::noiseFFT::Pf
(
const tmp<scalarField>& tpn
......@@ -193,7 +241,7 @@ Foam::tmp<Foam::scalarField> Foam::noiseFFT::Pf
scalarField::subField(tPn2(), tPn2().size()/2)
)
);
scalarField& Pn = tPn.ref();
scalarField& Pn = tPn();
Pn *= 2.0/sqrt(scalar(tPn2().size()));
Pn[0] /= 2.0;
......@@ -202,36 +250,22 @@ Foam::tmp<Foam::scalarField> Foam::noiseFFT::Pf
}
Foam::graph Foam::noiseFFT::meanPf
(
const label N,
const label nw
) const
Foam::graph Foam::noiseFFT::meanPf(const windowModel& window) const
{
if (N > size())
{
FatalErrorInFunction
<< "Requested window is outside set of data" << nl
<< "number of data = " << size() << nl
<< "size of window = " << N << nl
<< "window = " << nw
<< exit(FatalError);
}
const label N = window.nSamples();
const label nWindow = window.nWindow();
scalarField MeanPf(N/2, 0.0);
scalarField meanPf(N/2, 0.0);
scalarField Hwf(Hanning(N));
for (label wi=0; wi<nw; ++wi)
for (label windowI = 0; windowI < nWindow; ++windowI)
{
MeanPf += Pf(Hwf*window(N, wi));
meanPf += Pf(window.apply<scalar>(*this, windowI));
}
MeanPf /= nw;
scalarField f(MeanPf.size());
scalar deltaf = 1.0/(N*deltat_);
meanPf /= scalar(nWindow);
scalar deltaf = 1.0/(N*deltaT_);
scalarField f(meanPf.size());
forAll(f, i)
{
f[i] = i*deltaf;
......@@ -243,41 +277,27 @@ Foam::graph Foam::noiseFFT::meanPf
"f [Hz]",
"P(f) [Pa]",
f,
MeanPf
meanPf
);
}
Foam::graph Foam::noiseFFT::RMSmeanPf
(
const label N,
const label nw
) const
Foam::graph Foam::noiseFFT::RMSmeanPf(const windowModel& window) const
{
if (N > size())
{
FatalErrorInFunction
<< "Requested window is outside set of data" << endl
<< "number of data = " << size() << endl
<< "size of window = " << N << endl
<< "window = " << nw
<< exit(FatalError);
}
const label N = window.nSamples();