ENH: noiseModel updates

- Limit output to frequency range given by fLower and fUpper (if supplied)
- Enable noise models to be run outside of $FOAM_CASE directory
  - if relative paths are used, $FOAM_CASE is prepended to the noise
    dict and input file names
- Enable output to be customised, e.g.

    // Optional write options dictionary (all default to 'yes')
    writeOptions
    {
        writePrmsf  no;
        writeSPL    yes;
        writePSD    yes;
        writePSDf   no;
        writeOctaves yes;
    }
This commit is contained in:
Andrew Heather
2017-03-01 13:55:05 +00:00
parent bcde59e646
commit 34bc14a53d
7 changed files with 312 additions and 110 deletions

View File

@ -28,6 +28,7 @@ License
#include "IOmanip.H" #include "IOmanip.H"
#include "Pair.H" #include "Pair.H"
#include "OSspecific.H" #include "OSspecific.H"
#include "SubField.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
@ -160,6 +161,45 @@ Foam::scalarField& Foam::graph::y()
} }
void Foam::graph::setXRange(const scalar x0, const scalar x1)
{
if (x1 < x0)
{
FatalErrorInFunction
<< "When setting limits, x1 must be greater than x0" << nl
<< " x0: " << x0 << nl
<< " x1: " << x1 << nl
<< abort(FatalError);
}
label i0 = 0;
label i1 = 0;
forAll(x_, i)
{
if (x_[i] < x0)
{
i0 = i + 1;
}
if (x_[i] < x1)
{
i1 = i;
}
}
label nX = i1 - i0 + 1;
scalarField xNew(SubField<scalar>(x_, nX, i0));
x_.transfer(xNew);
forAllIter(HashPtrTable<curve>, *this, iter)
{
curve* c = iter();
scalarField cNew(SubField<scalar>(*c, nX, i0));
c->transfer(cNew);
}
}
Foam::autoPtr<Foam::graph::writer> Foam::graph::writer::New Foam::autoPtr<Foam::graph::writer> Foam::graph::writer::New
( (
const word& graphFormat const word& graphFormat

View File

@ -178,6 +178,9 @@ public:
scalarField& y(); scalarField& y();
// Limit the data for all axes based on x-limits
void setXRange(const scalar x0, const scalar x1);
// Write // Write

View File

@ -35,6 +35,26 @@ namespace Foam
// * * * * * * * * * * * * * Protected Member Functions * * * * * * * * * * // // * * * * * * * * * * * * * Protected Member Functions * * * * * * * * * * //
void Foam::noiseModel::readWriteOption
(
const dictionary& dict,
const word& lookup,
bool& option
) const
{
dict.readIfPresent(lookup, option);
if (option)
{
Info<< " " << lookup << ": " << "yes" << endl;
}
else
{
Info<< " " << lookup << ": " << "no" << endl;
}
}
Foam::scalar Foam::noiseModel::checkUniformTimeStep Foam::scalar Foam::noiseModel::checkUniformTimeStep
( (
const scalarList& times const scalarList& times
@ -92,6 +112,14 @@ Foam::label Foam::noiseModel::findStartTimeIndex
} }
Foam::fileName Foam::noiseModel::baseFileDir() const
{
fileName baseDir("$FOAM_CASE");
baseDir = baseDir.expand()/"postProcessing"/"noise";
return baseDir;
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::noiseModel::noiseModel(const dictionary& dict, const bool readFields) Foam::noiseModel::noiseModel(const dictionary& dict, const bool readFields)
@ -101,9 +129,15 @@ Foam::noiseModel::noiseModel(const dictionary& dict, const bool readFields)
nSamples_(65536), nSamples_(65536),
fLower_(25), fLower_(25),
fUpper_(10000), fUpper_(10000),
customBounds_(false),
startTime_(0), startTime_(0),
windowModelPtr_(), windowModelPtr_(),
graphFormat_("raw") graphFormat_("raw"),
writePrmsf_(true),
writeSPL_(true),
writePSD_(true),
writePSDf_(true),
writeOctaves_(true)
{ {
if (readFields) if (readFields)
{ {
@ -124,8 +158,11 @@ bool Foam::noiseModel::read(const dictionary& dict)
{ {
dict.readIfPresent("rhoRef", rhoRef_); dict.readIfPresent("rhoRef", rhoRef_);
dict.readIfPresent("N", nSamples_); dict.readIfPresent("N", nSamples_);
dict.readIfPresent("fl", fLower_); customBounds_ = false;
dict.readIfPresent("fu", fUpper_); if (dict.readIfPresent("fl", fLower_) || dict.readIfPresent("fu", fUpper_))
{
customBounds_ = true;
}
dict.readIfPresent("startTime", startTime_); dict.readIfPresent("startTime", startTime_);
dict.readIfPresent("graphFormat", graphFormat_); dict.readIfPresent("graphFormat", graphFormat_);
@ -164,6 +201,14 @@ bool Foam::noiseModel::read(const dictionary& dict)
} }
Info<< " Write options:" << endl;
dictionary optDict(dict.subOrEmptyDict("writeOptions"));
readWriteOption(optDict, "writePrmsf", writePrmsf_);
readWriteOption(optDict, "writeSPL", writeSPL_);
readWriteOption(optDict, "writePSD", writePSD_);
readWriteOption(optDict, "writeOctaves", writeOctaves_);
windowModelPtr_ = windowModel::New(dict, nSamples_); windowModelPtr_ = windowModel::New(dict, nSamples_);
return true; return true;

View File

@ -2,7 +2,7 @@
========= | ========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | \\ / O peration |
\\ / A nd | Copyright (C) 2015-2016 OpenCFD Ltd. \\ / A nd | Copyright (C) 2015-2017 OpenCFD Ltd.
\\/ M anipulation | \\/ M anipulation |
------------------------------------------------------------------------------- -------------------------------------------------------------------------------
License License
@ -35,6 +35,16 @@ Description
fl 25; fl 25;
fu 25; fu 25;
startTime 0; startTime 0;
// Optional write options dictionary
writeOptions
{
writePrmsf no;
writeSPL yes;
writePSD yes;
writePSDf no;
writeOctaves yes;
}
\endverbatim \endverbatim
where where
@ -46,6 +56,11 @@ Description
fu | Upper frequency bounds | no | 10000 fu | Upper frequency bounds | no | 10000
startTime | Start time | no | 0 startTime | Start time | no | 0
graphFormat | Graph format | no | raw graphFormat | Graph format | no | raw
writePrmsf | Write Prmsf data | no | yes
writeSPL | Write SPL data | no | yes
writePSD | Write PSD data | no | yes
writePSDf | Write PSDf data | no | yes
writeOctaves| Write octaves data | no | yes
\endtable \endtable
Note Note
@ -108,6 +123,9 @@ protected:
//- Upper frequency limit, default = 10kHz //- Upper frequency limit, default = 10kHz
scalar fUpper_; scalar fUpper_;
//- Flagto indicate that custom frequenct bounds are being used
bool customBounds_;
//- Start time, default = 0s //- Start time, default = 0s
scalar startTime_; scalar startTime_;
@ -118,8 +136,34 @@ protected:
word graphFormat_; word graphFormat_;
// Write options
//- Write Prmsf; default = yes
bool writePrmsf_;
//- Write SPL; default = yes
bool writeSPL_;
//- Write PSD; default = yes
bool writePSD_;
//- Write PSDf; default = yes
bool writePSDf_;
//- Write writeOctaves; default = yes
bool writeOctaves_;
// Protected Member Functions // Protected Member Functions
//- Helper function to read write options and provide info feedback
void readWriteOption
(
const dictionary& dict,
const word& lookup,
bool& option
) const;
//- Check and return uniform time step //- Check and return uniform time step
scalar checkUniformTimeStep scalar checkUniformTimeStep
( (
@ -133,6 +177,9 @@ protected:
const scalar startTime const scalar startTime
) const; ) const;
//- Return the base output directory
fileName baseFileDir() const;
public: public:

View File

@ -2,7 +2,7 @@
========= | ========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | \\ / O peration |
\\ / A nd | Copyright (C) 2015-2016 OpenCFD Ltd. \\ / A nd | Copyright (C) 2015-2017 OpenCFD Ltd.
\\/ M anipulation | \\/ M anipulation |
------------------------------------------------------------------------------- -------------------------------------------------------------------------------
License License
@ -86,7 +86,7 @@ void pointNoise::processData(const Function1Types::CSV<scalar>& data)
windowModelPtr_->validate(t.size()); windowModelPtr_->validate(t.size());
const windowModel& win = windowModelPtr_(); const windowModel& win = windowModelPtr_();
const scalar deltaf = 1.0/(deltaT*win.nSamples()); const scalar deltaf = 1.0/(deltaT*win.nSamples());
fileName outDir(fileName("postProcessing")/"noise"/typeName/fNameBase); fileName outDir(baseFileDir()/typeName/fNameBase);
// Create the fft // Create the fft
noiseFFT nfft(deltaT, p); noiseFFT nfft(deltaT, p);
@ -97,13 +97,27 @@ void pointNoise::processData(const Function1Types::CSV<scalar>& data)
// RMS pressure [Pa] // RMS pressure [Pa]
graph Prmsf(nfft.RMSmeanPf(win)); graph Prmsf(nfft.RMSmeanPf(win));
Info<< " Creating graph for " << Prmsf.title() << endl; if (customBounds_)
Prmsf.write(outDir, graph::wordify(Prmsf.title()), graphFormat_); {
Prmsf.setXRange(fLower_, fUpper_);
}
if (writePrmsf_)
{
Info<< " Creating graph for " << Prmsf.title() << endl;
Prmsf.write(outDir, graph::wordify(Prmsf.title()), graphFormat_);
}
// PSD [Pa^2/Hz] // PSD [Pa^2/Hz]
graph PSDf(nfft.PSDf(win)); graph PSDf(nfft.PSDf(win));
Info<< " Creating graph for " << PSDf.title() << endl; if (customBounds_)
PSDf.write(outDir, graph::wordify(PSDf.title()), graphFormat_); {
PSDf.setXRange(fLower_, fUpper_);
}
if (writePSDf_)
{
Info<< " Creating graph for " << PSDf.title() << endl;
PSDf.write(outDir, graph::wordify(PSDf.title()), graphFormat_);
}
// PSD [dB/Hz] // PSD [dB/Hz]
graph PSDg graph PSDg
@ -114,8 +128,12 @@ void pointNoise::processData(const Function1Types::CSV<scalar>& data)
Prmsf.x(), Prmsf.x(),
noiseFFT::PSD(PSDf.y()) noiseFFT::PSD(PSDf.y())
); );
Info<< " Creating graph for " << PSDg.title() << endl;
PSDg.write(outDir, graph::wordify(PSDg.title()), graphFormat_); if (writePSD_)
{
Info<< " Creating graph for " << PSDg.title() << endl;
PSDg.write(outDir, graph::wordify(PSDg.title()), graphFormat_);
}
// SPL [dB] // SPL [dB]
graph SPLg graph SPLg
@ -126,52 +144,59 @@ void pointNoise::processData(const Function1Types::CSV<scalar>& data)
Prmsf.x(), Prmsf.x(),
noiseFFT::SPL(PSDf.y()*deltaf) noiseFFT::SPL(PSDf.y()*deltaf)
); );
Info<< " Creating graph for " << SPLg.title() << endl;
SPLg.write(outDir, graph::wordify(SPLg.title()), graphFormat_);
labelList octave13BandIDs; if (writeSPL_)
scalarField octave13FreqCentre; {
noiseFFT::octaveBandInfo Info<< " Creating graph for " << SPLg.title() << endl;
( SPLg.write(outDir, graph::wordify(SPLg.title()), graphFormat_);
Prmsf.x(), }
fLower_,
fUpper_, if (writeOctaves_)
3, {
octave13BandIDs, labelList octave13BandIDs;
octave13FreqCentre scalarField octave13FreqCentre;
); noiseFFT::octaveBandInfo
(
Prmsf.x(),
fLower_,
fUpper_,
3,
octave13BandIDs,
octave13FreqCentre
);
// 1/3 octave data // 1/3 octave data
// --------------- // ---------------
// PSD [Pa^2/Hz] // PSD [Pa^2/Hz]
graph PSD13f(nfft.octaves(PSDf, octave13BandIDs, false)); graph PSD13f(nfft.octaves(PSDf, octave13BandIDs, false));
// Integrated PSD = P(rms)^2 [Pa^2] // Integrated PSD = P(rms)^2 [Pa^2]
graph Prms13f2(nfft.octaves(PSDf, octave13BandIDs, true)); graph Prms13f2(nfft.octaves(PSDf, octave13BandIDs, true));
graph PSD13g graph PSD13g
( (
"PSD13_dB_Hz(fm)", "PSD13_dB_Hz(fm)",
"fm [Hz]", "fm [Hz]",
"PSD(fm) [dB_Hz]", "PSD(fm) [dB_Hz]",
octave13FreqCentre, octave13FreqCentre,
noiseFFT::PSD(PSD13f.y()) noiseFFT::PSD(PSD13f.y())
); );
Info<< " Creating graph for " << PSD13g.title() << endl; Info<< " Creating graph for " << PSD13g.title() << endl;
PSD13g.write(outDir, graph::wordify(PSD13g.title()), graphFormat_); PSD13g.write(outDir, graph::wordify(PSD13g.title()), graphFormat_);
graph SPL13g graph SPL13g
( (
"SPL13_dB(fm)", "SPL13_dB(fm)",
"fm [Hz]", "fm [Hz]",
"SPL(fm) [dB]", "SPL(fm) [dB]",
octave13FreqCentre, octave13FreqCentre,
noiseFFT::SPL(Prms13f2.y()) noiseFFT::SPL(Prms13f2.y())
); );
Info<< " Creating graph for " << SPL13g.title() << endl; Info<< " Creating graph for " << SPL13g.title() << endl;
SPL13g.write(outDir, graph::wordify(SPL13g.title()), graphFormat_); SPL13g.write(outDir, graph::wordify(SPL13g.title()), graphFormat_);
}
} }
@ -186,18 +211,15 @@ void pointNoise::calculate()
} }
if (inputFileNames_.size()) forAll(inputFileNames_, i)
{ {
forAll(inputFileNames_, i) fileName fName = inputFileNames_[i];
if (!fName.isAbsolute())
{ {
const fileName fName = inputFileNames_[i].expand(); fName = "$FOAM_CASE"/fName;
Function1Types::CSV<scalar> data("pressure", dict_, "Data", fName);
processData(data);
} }
} fName.expand();
else Function1Types::CSV<scalar> data("pressure", dict_, "Data", fName);
{
Function1Types::CSV<scalar> data("pressure", dict_, "Data");
processData(data); processData(data);
} }
} }
@ -226,7 +248,15 @@ bool pointNoise::read(const dictionary& dict)
{ {
if (noiseModel::read(dict)) if (noiseModel::read(dict))
{ {
dict.readIfPresent("inputFiles", inputFileNames_); if (!dict.readIfPresent("inputFiles", inputFileNames_))
{
inputFileNames_.setSize(1);
// Note: unsafe lookup of file name from pressureData sub-dict!
dict.subDict("pressureData").lookup("fileName")
>> inputFileNames_[0];
}
return true; return true;
} }

View File

@ -2,7 +2,7 @@
========= | ========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | \\ / O peration |
\\ / A nd | Copyright (C) 2015-2016 OpenCFD Ltd. \\ / A nd | Copyright (C) 2015-2017 OpenCFD Ltd.
\\/ M anipulation | \\/ M anipulation |
------------------------------------------------------------------------------- -------------------------------------------------------------------------------
License License
@ -231,15 +231,13 @@ Foam::scalar surfaceNoise::writeSurfaceData
const word& title, const word& title,
const scalar freq, const scalar freq,
const scalarField& data, const scalarField& data,
const labelList& procFaceOffset const labelList& procFaceOffset,
const bool writeSurface
) const ) const
{ {
Info<< " processing " << title << " for frequency " << freq << endl; Info<< " processing " << title << " for frequency " << freq << endl;
fileName outDir fileName outDir(baseFileDir()/groupName/Foam::name(freq));
(
fileName("postProcessing")/"noise"/groupName/Foam::name(freq)
);
if (Pstream::parRun()) if (Pstream::parRun())
{ {
@ -279,20 +277,23 @@ Foam::scalar surfaceNoise::writeSurfaceData
} }
} }
// could also have meshedSurface implement meshedSurf // Could also have meshedSurface implement meshedSurf
fileName outFileName = writerPtr_->write if (writeSurface)
( {
outDir, fileName outFileName = writerPtr_->write
fName,
meshedSurfRef
( (
surf.points(), outDir,
surf.surfFaces() fName,
), meshedSurfRef
title, (
allData, surf.points(),
false surf.surfFaces()
); ),
title,
allData,
false
);
}
// TO BE VERIFIED: area-averaged values // TO BE VERIFIED: area-averaged values
// areaAverage = sum(allData*surf.magSf())/sum(surf.magSf()); // areaAverage = sum(allData*surf.magSf())/sum(surf.magSf());
@ -306,20 +307,23 @@ Foam::scalar surfaceNoise::writeSurfaceData
{ {
const meshedSurface& surf = readerPtr_->geometry(); const meshedSurface& surf = readerPtr_->geometry();
// could also have meshedSurface implement meshedSurf // Could also have meshedSurface implement meshedSurf
writerPtr_->write if (writeSurface)
( {
outDir, writerPtr_->write
fName,
meshedSurfRef
( (
surf.points(), outDir,
surf.surfFaces() fName,
), meshedSurfRef
title, (
data, surf.points(),
false surf.surfFaces()
); ),
title,
data,
false
);
}
// TO BE VERIFIED: area-averaged values // TO BE VERIFIED: area-averaged values
// return sum(data*surf.magSf())/sum(surf.magSf()); // return sum(data*surf.magSf())/sum(surf.magSf());
@ -463,6 +467,11 @@ void surfaceNoise::calculate()
{ {
fileName fName = inputFileNames_[i]; fileName fName = inputFileNames_[i];
if (!fName.isAbsolute())
{
fName = "$FOAM_CASE"/fName;
}
initialise(fName.expand()); initialise(fName.expand());
// Container for pressure time history data per face // Container for pressure time history data per face
@ -575,19 +584,23 @@ void surfaceNoise::calculate()
const word& fNameBase = fName.name(true); const word& fNameBase = fName.name(true);
// Output directory for graphs // Output directory for graphs
fileName outDir fileName outDir(baseFileDir()/typeName/fNameBase);
(
fileName("postProcessing")/"noise"/typeName/fNameBase
);
const scalar deltaf = 1.0/(deltaT_*win.nSamples()); const scalar deltaf = 1.0/(deltaT_*win.nSamples());
Info<< "Writing fft surface data" << endl; Info<< "Writing fft surface data" << endl;
{ {
scalarField PrmsfAve(surfPrmsf.size(), 0); // Determine frequency range of interest
scalarField PSDfAve(surfPrmsf.size(), 0); // Note: freqencies have fixed interval, and are in the range
scalarField fOut(surfPrmsf.size(), 0); // 0 to (n-1)*deltaf
label f0 = ceil(fLower_/deltaf);
label f1 = floor(fUpper_/deltaf);
label nFreq = f1 - f0 + 1;
forAll(surfPrmsf, i) scalarField PrmsfAve(nFreq, 0);
scalarField PSDfAve(nFreq, 0);
scalarField fOut(nFreq, 0);
for (label i = f0; i <= f1; ++i)
{ {
label freqI = (i + 1)*fftWriteInterval_ - 1; label freqI = (i + 1)*fftWriteInterval_ - 1;
fOut[i] = freq1[freqI]; fOut[i] = freq1[freqI];
@ -599,7 +612,8 @@ void surfaceNoise::calculate()
"Prmsf", "Prmsf",
freq1[freqI], freq1[freqI],
surfPrmsf[i], surfPrmsf[i],
procFaceOffset procFaceOffset,
writePrmsf_
); );
PSDfAve[i] = writeSurfaceData PSDfAve[i] = writeSurfaceData
@ -609,7 +623,8 @@ void surfaceNoise::calculate()
"PSDf", "PSDf",
freq1[freqI], freq1[freqI],
surfPSDf[i], surfPSDf[i],
procFaceOffset procFaceOffset,
writePSDf_
); );
writeSurfaceData writeSurfaceData
( (
@ -618,7 +633,8 @@ void surfaceNoise::calculate()
"PSD", "PSD",
freq1[freqI], freq1[freqI],
noiseFFT::PSD(surfPSDf[i]), noiseFFT::PSD(surfPSDf[i]),
procFaceOffset procFaceOffset,
writePSD_
); );
writeSurfaceData writeSurfaceData
( (
@ -627,7 +643,8 @@ void surfaceNoise::calculate()
"SPL", "SPL",
freq1[freqI], freq1[freqI],
noiseFFT::SPL(surfPSDf[i]*deltaf), noiseFFT::SPL(surfPSDf[i]*deltaf),
procFaceOffset procFaceOffset,
writeSPL_
); );
} }
@ -688,7 +705,8 @@ void surfaceNoise::calculate()
"PSD13f", "PSD13f",
octave13FreqCentre[i], octave13FreqCentre[i],
surfPSD13f[i], surfPSD13f[i],
procFaceOffset procFaceOffset,
writeOctaves_
); );
writeSurfaceData writeSurfaceData
( (
@ -697,7 +715,8 @@ void surfaceNoise::calculate()
"PSD13", "PSD13",
octave13FreqCentre[i], octave13FreqCentre[i],
noiseFFT::PSD(surfPSD13f[i]), noiseFFT::PSD(surfPSD13f[i]),
procFaceOffset procFaceOffset,
writeOctaves_
); );
writeSurfaceData writeSurfaceData
( (
@ -706,7 +725,8 @@ void surfaceNoise::calculate()
"SPL13", "SPL13",
octave13FreqCentre[i], octave13FreqCentre[i],
noiseFFT::SPL(surfPrms13f2[i]), noiseFFT::SPL(surfPrms13f2[i]),
procFaceOffset procFaceOffset,
writeOctaves_
); );
Prms13f2Ave[i] = Prms13f2Ave[i] =

View File

@ -2,7 +2,7 @@
========= | ========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | \\ / O peration |
\\ / A nd | Copyright (C) 2015-2016 OpenCFD Ltd. \\ / A nd | Copyright (C) 2015-2017 OpenCFD Ltd.
\\/ M anipulation | \\/ M anipulation |
------------------------------------------------------------------------------- -------------------------------------------------------------------------------
License License
@ -71,7 +71,23 @@ Description
{ {
collateTimes true; collateTimes true;
} }
// Write Prmsf; default = yes
writePrmsf no;
// Write SPL; default = yes
writeSPL yes;
// Write PSD; default = yes
writePSD yes;
// Write PSDf; default = yes
writePSDf no;
// Write writeOctaves; default = yes
writeOctaves yes;
} }
\endverbatim \endverbatim
SourceFiles SourceFiles
@ -172,7 +188,8 @@ protected:
const word& title, const word& title,
const scalar freq, const scalar freq,
const scalarField& data, const scalarField& data,
const labelList& procFaceOffset const labelList& procFaceOffset,
const bool writeSurface
) const; ) const;
//- Calculate the area average value //- Calculate the area average value