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);
85 while (fTest < fLower)
98 if (stepi) missedBins.
append(fTest/fRatio*fRatioL2C);
104 if (bandIDs.insert(i))
107 fc.append(fTest*fRatioL2C);
118 fBandIDs = bandIDs.sortedToc();
120 if (missedBins.size())
122 label nMiss = missedBins.size();
123 label nTotal = nMiss + fc.size() - 1;
125 <<
"Empty bands found: " << nMiss <<
" of " << nTotal
126 <<
" with centre frequencies " <<
flatOutput(missedBins)
135 fCentre.transfer(fc);
145 auto& result = tresult.ref();
172 option = option && Pstream::master();
190 scalar deltaT = -1.0;
192 if (times.size() > 1)
195 deltaT = (times.last() - times.first())/scalar(times.size() - 1);
197 for (label i = 1; i < times.size(); ++i)
199 scalar dT = times[i] - times[i-1];
201 if (
mag(dT/deltaT - 1) > 1
e-8)
204 <<
"Unable to process data with a variable time step"
212 <<
"Unable to create FFT with a single value"
227 <<
"Pressure data at position " << i
228 <<
" is outside of permitted bounds:" <<
nl
229 <<
" pressure: " <<
p[i] <<
nl
250 const instant& i = allTimes[timeI];
282 const auto& window = windowModelPtr_();
283 const label
N = window.nSamples();
288 const scalar deltaf = 1.0/(
N*deltaT);
295 if (
f[i] > fLower_ &&
f[i] < fUpper_)
301 if (
check && nFreq == 0)
304 <<
"No frequencies found in range "
305 << fLower_ <<
" to " << fUpper_
320 if (freqBandIDs.
size() < 2)
323 <<
"Octave frequency bands are not defined "
324 <<
"- skipping octaves calculation"
331 auto& octData = toctData.ref();
334 for (label bandI = 0; bandI < freqBandIDs.
size() - 1; ++bandI)
336 label fb0 = freqBandIDs[bandI];
337 label fb1 = freqBandIDs[bandI+1];
339 if (fb0 == fb1)
continue;
341 for (label freqI = fb0; freqI < fb1; ++freqI)
344 label f1 =
f[freqI + 1];
345 scalar dataAve = 0.5*(
data[freqI] +
data[freqI + 1]);
346 octData[bandI] += dataAve*(f1 - f0);
353 labelList bandUnused = bandUsed.sortedToc();
354 if (bandUnused.size())
357 <<
"Empty bands found: " << bandUnused.size() <<
" of "
358 << bandUsed.size() <<
endl;
367 if (planInfo_.active)
369 if (
p.size() != planInfo_.windowSize)
372 <<
"Expected pressure data to have " << planInfo_.windowSize
373 <<
" values, but received " <<
p.size() <<
" values"
384 ::fftw_execute(planInfo_.plan);
386 const label
n = planInfo_.windowSize;
387 const label nBy2 =
n/2;
389 auto& result = tresult.ref();
394 for (label i = 1; i <= nBy2; ++i)
396 const auto re = out[i];
397 const auto im = out[
n - i];
410 const auto& window = windowModelPtr_();
411 const label
N = window.nSamples();
412 const label nWindow = window.nWindow();
415 auto& meanPf = tmeanPf.ref();
417 for (label windowI = 0; windowI < nWindow; ++windowI)
419 meanPf += Pf(window.apply<scalar>(
p, windowI));
422 meanPf /= scalar(nWindow);
433 const auto& window = windowModelPtr_();
434 const label
N = window.nSamples();
435 const label nWindow = window.nWindow();
438 for (label windowI = 0; windowI < nWindow; ++windowI)
440 RMSMeanPf +=
sqr(Pf(window.apply<scalar>(
p, windowI)));
443 return sqrt(RMSMeanPf/scalar(nWindow))/scalar(
N);
453 const auto& window = windowModelPtr_();
454 const label
N = window.nSamples();
455 const label nWindow = window.nWindow();
458 auto& psd = tpsd.ref();
460 for (label windowI = 0; windowI < nWindow; ++windowI)
462 psd +=
sqr(Pf(window.apply<scalar>(
p, windowI)));
465 scalar fs = 1.0/deltaT;
466 psd /= scalar(nWindow)*fs*
N;
485 const scalar
c1 =
sqr(12194.0);
486 const scalar
c2 =
sqr(20.6);
487 const scalar c3 =
sqr(107.7);
488 const scalar c4 =
sqr(739.9);
490 const scalar f2 =
f*
f;
495 (f2 +
c2)*
sqrt((f2 + c3)*(f2 + c4))*(f2 +
c1)
513 const scalar
c1 =
sqr(12194.0);
514 const scalar
c2 =
sqr(20.6);
515 const scalar c3 =
sqr(158.5);
517 const scalar f2 =
f*
f;
540 const scalar
c1 =
sqr(12194.0);
541 const scalar
c2 =
sqr(20.6);
543 const scalar f2 =
f*
f;
545 return c1*f2/((f2 +
c2)*(f2 +
c1));
562 const scalar f2 =
f*
f;
565 (
sqr(1037918.48 - f2) + 1080768.16*f2)
566 /(
sqr(9837328 - f2) + 11723776*f2);
569 f/6.8966888496476e-5*
Foam::sqrt(hf/((f2 + 79919.29)*(f2 + 1345600)));
598 minPressure_(-0.5*VGREAT),
599 maxPressure_(0.5*VGREAT),
639 <<
"fl: lower frequency bound must be greater than zero"
646 <<
"fu: upper frequency bound must be greater than zero"
650 if (fUpper_ < fLower_)
653 <<
"fu: upper frequency bound must be greater than lower "
654 <<
"frequency bound (fl)"
659 Info<<
" Frequency bounds:" <<
nl
660 <<
" lower: " << fLower_ <<
nl
661 <<
" upper: " << fUpper_ <<
endl;
663 weightingTypeNames_.readIfPresent(
"SPLweighting",
dict, SPLweighting_);
665 Info<<
" Weighting: " << weightingTypeNames_[SPLweighting_] <<
endl;
669 Info<<
" Reference for dB calculation: " << dBRef_ <<
endl;
674 readWriteOption(optDict,
"writePrmsf", writePrmsf_);
675 readWriteOption(optDict,
"writeSPL", writeSPL_);
676 readWriteOption(optDict,
"writePSD", writePSD_);
677 readWriteOption(optDict,
"writePSDf", writePSDf_);
678 readWriteOption(optDict,
"writeOctaves", writeOctaves_);
685 const label windowSize = windowModelPtr_->nSamples();
689 planInfo_.active =
true;
690 planInfo_.windowSize = windowSize;
691 planInfo_.in.setSize(windowSize);
692 planInfo_.out.setSize(windowSize);
700 planInfo_.out.data(),
702 windowSize <= 8192 ? FFTW_MEASURE : FFTW_ESTIMATE
730 switch (SPLweighting_)
736 case weightingType::dBA:
741 case weightingType::dBB:
746 case weightingType::dBC:
751 case weightingType::dBD:
759 <<
"Unknown weighting " << weightingTypeNames_[SPLweighting_]
777 switch (SPLweighting_)
783 case weightingType::dBA:
787 spl[i] += gainA(
f[i]);
791 case weightingType::dBB:
795 spl[i] += gainB(
f[i]);
799 case weightingType::dBC:
803 spl[i] += gainC(
f[i]);
807 case weightingType::dBD:
811 spl[i] += gainD(
f[i]);
818 <<
"Unknown weighting " << weightingTypeNames_[SPLweighting_]
829 if (planInfo_.active)
831 planInfo_.active =
false;
832 fftw_destroy_plan(planInfo_.plan);
843 OFstream osA(
"noiseModel-weight-A");
844 OFstream osB(
"noiseModel-weight-B");
845 OFstream osC(
"noiseModel-weight-C");
846 OFstream osD(
"noiseModel-weight-D");
848 for (label
f = f0;
f <= f1; ++
f)
850 osA <<
f <<
" " << gainA(
f) <<
nl;
851 osB <<
f <<
" " << gainB(
f) <<
nl;
852 osC <<
f <<
" " << gainC(
f) <<
nl;
853 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.
A bitSet stores bits (elements with only two states) in packed internal format and supports a variety...
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.
tmp< scalarField > safeLog10(const scalarField &fld)
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].
static void check(const int retVal, const char *what)
messageStream Info
Information stream (stdout output on master, null elsewhere)
dimensionedScalar log10(const dimensionedScalar &ds)
DynamicList< T, SizeMin > & append(const T &val)
Append an element to the end of this list.
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.
gmvFile<< "tracers "<< particles.size()<< nl;for(const passiveParticle &p :particles){ gmvFile<< p.position().x()<< ' ';}gmvFile<< nl;for(const passiveParticle &p :particles){ gmvFile<< p.position().y()<< ' ';}gmvFile<< nl;for(const passiveParticle &p :particles){ gmvFile<< p.position().z()<< ' ';}gmvFile<< nl;for(const word &name :lagrangianScalarNames){ IOField< scalar > fld(IOobject(name, runTime.timeName(), cloud::prefix, mesh, IOobject::MUST_READ, IOobject::NO_WRITE))
scalar gainC(const scalar f) const
C weighting as gain in dB.
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.
FlatOutput::OutputAdaptor< Container, Delimiters > flatOutput(const Container &obj, Delimiters delim)
Global flatOutput() function with specified output delimiters.
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.
tmp< scalarField > uniformFrequencies(const scalar deltaT, const bool check) const
word name(const expressions::valueTypeCode typeCode)
A word representation of a valueTypeCode. Empty for INVALID.
#define FatalIOErrorInFunction(ios)
Report an error message using Foam::FatalIOError.
An instant of time. Contains the time value and name.
void size(const label n)
Older name for setAddressableSize.
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].