Go to the documentation of this file.
47 {weightingType::dBA,
"dBA"},
48 {weightingType::dBB,
"dBB"},
49 {weightingType::dBC,
"dBC"},
50 {weightingType::dBD,
"dBD"},
71 scalar fTest = 15.625;
73 const scalar fRatio =
pow(2, 1.0/octave);
74 const scalar fRatioL2C =
pow(2, 0.5/octave);
96 if (bandIDs.insert(i))
99 fc.append(fTest*fRatioL2C);
110 fBandIDs = bandIDs.sortedToc();
117 fCentre.transfer(fc);
134 option = option && Pstream::master();
152 scalar deltaT = -1.0;
154 if (times.size() > 1)
157 deltaT = (times.last() - times.first())/scalar(times.size() - 1);
159 for (label i = 1; i < times.size(); ++i)
161 scalar dT = times[i] - times[i-1];
163 if (
mag(dT/deltaT - 1) > 1
e-8)
166 <<
"Unable to process data with a variable time step"
174 <<
"Unable to create FFT with a single value"
189 <<
"Pressure data at position " << i
190 <<
" is outside of permitted bounds:" <<
nl
191 <<
" pressure: " <<
p[i] <<
nl
212 const instant& i = allTimes[timeI];
243 const auto& window = windowModelPtr_();
244 const label
N = window.nSamples();
249 const scalar deltaf = 1.0/(
N*deltaT);
268 if (freqBandIDs.
size() < 2)
271 <<
"Octave frequency bands are not defined "
272 <<
"- skipping octaves calculation"
278 auto& octData = toctData.ref();
280 for (label bandI = 0; bandI < freqBandIDs.
size() - 1; ++bandI)
282 label fb0 = freqBandIDs[bandI];
283 label fb1 = freqBandIDs[bandI+1];
285 if (fb0 == fb1)
continue;
287 for (label freqI = fb0; freqI < fb1; ++freqI)
290 label f1 =
f[freqI + 1];
291 scalar dataAve = 0.5*(
data[freqI] +
data[freqI + 1]);
292 octData[bandI] += dataAve*(f1 - f0);
302 if (planInfo_.active)
304 if (
p.size() != planInfo_.windowSize)
307 <<
"Expected pressure data to have " << planInfo_.windowSize
308 <<
" values, but received " <<
p.size() <<
" values"
319 ::fftw_execute(planInfo_.plan);
321 const label
n = planInfo_.windowSize;
322 const label nBy2 =
n/2;
324 auto& result = tresult.ref();
329 for (label i = 1; i <= nBy2; ++i)
331 const auto re = out[i];
332 const auto im = out[
n - i];
345 const auto& window = windowModelPtr_();
346 const label
N = window.nSamples();
347 const label nWindow = window.nWindow();
350 auto& meanPf = tmeanPf.ref();
352 for (label windowI = 0; windowI < nWindow; ++windowI)
354 meanPf += Pf(window.apply<scalar>(
p, windowI));
357 meanPf /= scalar(nWindow);
368 const auto& window = windowModelPtr_();
369 const label
N = window.nSamples();
370 const label nWindow = window.nWindow();
373 for (label windowI = 0; windowI < nWindow; ++windowI)
375 RMSMeanPf +=
sqr(Pf(window.apply<scalar>(
p, windowI)));
378 return sqrt(RMSMeanPf/scalar(nWindow))/scalar(
N);
388 const auto& window = windowModelPtr_();
389 const label
N = window.nSamples();
390 const label nWindow = window.nWindow();
393 auto& psd = tpsd.ref();
395 for (label windowI = 0; windowI < nWindow; ++windowI)
397 psd +=
sqr(Pf(window.apply<scalar>(
p, windowI)));
400 scalar fs = 1.0/deltaT;
401 psd /= scalar(nWindow)*fs*
N;
420 const scalar
c1 =
sqr(12194.0);
421 const scalar
c2 =
sqr(20.6);
422 const scalar c3 =
sqr(107.7);
423 const scalar c4 =
sqr(739.9);
425 const scalar f2 =
f*
f;
430 (f2 +
c2)*
sqrt((f2 + c3)*(f2 + c4))*(f2 +
c1)
443 const scalar
c1 =
sqr(12194.0);
444 const scalar
c2 =
sqr(20.6);
445 const scalar c3 =
sqr(158.5);
447 const scalar f2 =
f*
f;
465 const scalar
c1 =
sqr(12194.0);
466 const scalar
c2 =
sqr(20.6);
468 const scalar f2 =
f*
f;
470 return c1*f2/((f2 +
c2)*(f2 +
c1));
482 const scalar f2 =
f*
f;
485 (
sqr(1037918.48 - f2) + 1080768.16*f2)
486 /(
sqr(9837328 - f2) + 11723776*f2);
489 f/6.8966888496476e-5*
Foam::sqrt(hf/((f2 + 79919.29)*(f2 + 1345600)));
512 minPressure_(-0.5*VGREAT),
513 maxPressure_(0.5*VGREAT),
553 <<
"fl: lower frequency bound must be greater than zero"
560 <<
"fu: upper frequency bound must be greater than zero"
564 if (fUpper_ < fLower_)
567 <<
"fu: upper frequency bound must be greater than lower "
568 <<
"frequency bound (fl)"
573 Info<<
" Frequency bounds:" <<
nl
574 <<
" lower: " << fLower_ <<
nl
575 <<
" upper: " << fUpper_ <<
endl;
577 weightingTypeNames_.readIfPresent(
"SPLweighting",
dict, SPLweighting_);
579 Info<<
" Weighting: " << weightingTypeNames_[SPLweighting_] <<
endl;
583 readWriteOption(optDict,
"writePrmsf", writePrmsf_);
584 readWriteOption(optDict,
"writeSPL", writeSPL_);
585 readWriteOption(optDict,
"writePSD", writePSD_);
586 readWriteOption(optDict,
"writePSDf", writePSDf_);
587 readWriteOption(optDict,
"writeOctaves", writeOctaves_);
594 const label windowSize = windowModelPtr_->nSamples();
598 planInfo_.active =
true;
599 planInfo_.windowSize = windowSize;
600 planInfo_.in.setSize(windowSize);
601 planInfo_.out.setSize(windowSize);
609 planInfo_.out.data(),
611 windowSize <= 8192 ? FFTW_MEASURE : FFTW_ESTIMATE
639 switch (SPLweighting_)
645 case weightingType::dBA:
650 case weightingType::dBB:
655 case weightingType::dBC:
660 case weightingType::dBD:
668 <<
"Unknown weighting " << weightingTypeNames_[SPLweighting_]
686 switch (SPLweighting_)
692 case weightingType::dBA:
696 spl[i] += gainA(
f[i]);
700 case weightingType::dBB:
704 spl[i] += gainB(
f[i]);
708 case weightingType::dBC:
712 spl[i] += gainC(
f[i]);
716 case weightingType::dBD:
720 spl[i] += gainD(
f[i]);
727 <<
"Unknown weighting " << weightingTypeNames_[SPLweighting_]
738 if (planInfo_.active)
740 planInfo_.active =
false;
741 fftw_destroy_plan(planInfo_.plan);
752 OFstream osA(
"noiseModel-weight-A");
753 OFstream osB(
"noiseModel-weight-B");
754 OFstream osC(
"noiseModel-weight-C");
755 OFstream osD(
"noiseModel-weight-D");
757 for (label
f = f0;
f <= f1; ++
f)
759 osA <<
f <<
" " << gainA(
f) <<
nl;
760 osB <<
f <<
" " << gainB(
f) <<
nl;
761 osC <<
f <<
" " << gainC(
f) <<
nl;
762 osD <<
f <<
" " << gainD(
f) <<
nl;
int debug
Static debugging option.
scalar checkUniformTimeStep(const scalarList ×) const
Check and return uniform time step.
Enum is a wrapper around a list of names/values that represent particular enumeration (or int) values...
virtual bool read(const dictionary &dict)
Read from dictionary.
void readWriteOption(const dictionary &dict, const word &lookup, bool &option) const
Helper function to read write options and provide info feedback.
A class for handling words, derived from Foam::string.
A class for handling file names.
static fileName envGlobalPath()
Global case (directory) from environment variable.
A class for managing temporary objects.
static constexpr const zero Zero
Global zero (0)
defineRunTimeSelectionTable(reactionRateFlameArea, dictionary)
scalar RBf(const scalar f) const
B weighting function.
tmp< scalarField > octaves(const scalarField &data, const scalarField &f, const labelUList &freqBandIDs) const
Generate octave data.
label findStartTimeIndex(const instantList &allTimes, const scalar startTime) const
Find and return start time index.
Ostream & endl(Ostream &os)
Add newline and flush stream.
prefixOSstream Pout
OSstream wrapped stdout (std::cout) with parallel prefix.
static autoPtr< windowModel > New(const dictionary &dict, const label nSamples)
Return a reference to the selected window model.
scalar RCf(const scalar f) const
C weighting function.
#define forAll(list, i)
Loop across all elements in list.
scalar gainB(const scalar f) const
B weighting as gain in dB.
scalar minPressure_
Min pressure value.
const dimensionedScalar c1
First radiation constant: default SI units: [W/m2].
messageStream Info
Information stream (uses stdout - output is on the master only)
word name(const complex &c)
Return string representation of complex.
dimensionedScalar log10(const dimensionedScalar &ds)
static void setOctaveBands(const scalarField &f, const scalar fLower, const scalar fUpper, const scalar octave, labelList &fBandIDs, scalarField &fCentre)
static word outputPrefix
Directory prefix.
void writeWeightings() const
Helper function to check weightings.
dimensionedScalar pow(const dimensionedScalar &ds, const dimensionedScalar &expt)
Lookup type of boundary radiation properties.
scalar gainC(const scalar f) const
C weighting as gain in dB.
tmp< scalarField > uniformFrequencies(const scalar deltaT) const
const dimensionedScalar re
Classical electron radius: default SI units: [m].
A list of keyword definitions, which are a keyword followed by a number of values (eg,...
bool validateBounds(const scalarList &p) const
Return true if all pressure data is within min/max bounds.
scalar maxPressure_
Min pressure value.
errorManip< error > abort(error &err)
const dimensionedScalar c2
Second radiation constant: default SI units: [m.K].
bool none(const UList< bool > &bools)
True if no entries are 'true'.
dictionary subOrEmptyDict(const word &keyword, enum keyType::option matchOpt=keyType::REGEX, const bool mandatory=false) const
scalar RDf(const scalar f) const
D weighting function.
static const Enum< weightingType > weightingTypeNames_
errorManipArg< error, int > exit(error &err, const int errNo=1)
Output to file stream, using an OSstream.
void readFields(const typename GeoFieldType::Mesh &mesh, const IOobjectList &objects, const wordHashSet &selectedFields, LIFOStack< regIOobject * > &storedObjects)
Read the selected GeometricFields of the templated type.
tmp< scalarField > RMSmeanPf(const scalarField &p) const
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
dimensionedSymmTensor sqr(const dimensionedVector &dv)
void cleanFFTW()
Clean up the FFTW.
tmp< scalarField > SPL(const scalarField &Prms2, const scalar f) const
SPL [dB].
dimensionedScalar sqrt(const dimensionedScalar &ds)
fileName baseFileDir(const label dataseti) const
Return the base output directory.
dimensioned< typename typeOfMag< Type >::type > mag(const dimensioned< Type > &dt)
scalar value() const
The value (const access)
fileName::Type type(const fileName &name, const bool followLink=true)
Return the file type: DIRECTORY or FILE, normally following symbolic links.
const dimensionedScalar e
Elementary charge.
scalar gainD(const scalar f) const
D weighting as gain in dB.
tmp< scalarField > PSDf(const scalarField &p, const scalar deltaT) const
static tmp< T > New(Args &&... args)
Construct tmp of T with forwarding arguments.
void size(const label n) noexcept
Override size to be inconsistent with allocated storage.
#define FatalIOErrorInFunction(ios)
Report an error message using Foam::FatalIOError.
An instant of time. Contains the time value and name.
const Vector< label > N(dict.get< Vector< label >>("N"))
scalar gainA(const scalar f) const
A weighting as gain in dB.
tmp< Foam::scalarField > PSD(const scalarField &PSDf) const
PSD [dB/Hz].
defineTypeNameAndDebug(combustionModel, 0)
#define WarningInFunction
Report a warning using Foam::Warning.
scalar RAf(const scalar f) const
A weighting function.
Database for solution data, solver performance and other reduced data.
planInfo planInfo_
Plan information for FFTW.
bool readIfPresent(const word &keyword, T &val, enum keyType::option matchOpt=keyType::REGEX) const
tmp< scalarField > Pf(const scalarField &p) const
Return the fft of the given pressure data.
static tmp< complexField > realTransform1D(const scalarField &field)
Transform real-value data.
dimensioned< Type > average(const DimensionedField< Type, GeoMesh > &df)
tmp< scalarField > meanPf(const scalarField &p) const
Return the multi-window mean fft of the complete pressure data [Pa].