46 {weightingType::none,
"dB"},
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);
107 fc.
append(fTest*fRatioL2C);
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)
145 auto& result = tresult.ref();
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"
224 if ((
p[i] < minPressure_) || (
p[i] > maxPressure_))
227 <<
"Pressure data at position " << i
228 <<
" is outside of permitted bounds:" <<
nl
229 <<
" pressure: " <<
p[i] <<
nl
230 <<
" minimum pressure: " << minPressure_ <<
nl
231 <<
" maximum pressure: " << maxPressure_ <<
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);
354 if (bandUnused.
size())
357 <<
"Empty bands found: " << bandUnused.
size() <<
" of "
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];
398 result[i] =
sqrt(re*re + im*im);
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;
522 (f2 + c2)*
sqrt(f2 + c3)*(f2 + c1)
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_)
732 case weightingType::none:
736 case weightingType::dBA:
741 case weightingType::dBB:
746 case weightingType::dBC:
751 case weightingType::dBD:
759 <<
"Unknown weighting " << weightingTypeNames_[SPLweighting_]
777 switch (SPLweighting_)
779 case weightingType::none:
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;
Info<< nl<< "Wrote faMesh in vtk format: "<< writer.output().name()<< nl;}{ vtk::lineWriter writer(aMesh.points(), aMesh.edges(), fileName(aMesh.mesh().time().globalPath()/"finiteArea-edges"));writer.writeGeometry();writer.beginCellData(4);writer.writeProcIDs();{ Field< scalar > fld(faMeshTools::flattenEdgeField(aMesh.magLe(), true))
A 1D vector of objects of type <T> that resizes itself as necessary to accept the new objects.
T remove()
Remove and return the last element. Fatal on an empty list.
void append(const T &val)
Copy append an element to the end of this list.
Enum is a wrapper around a list of names/values that represent particular enumeration (or int) values...
bool insert(const Key &key)
Insert a new entry, not overwriting existing entries.
List< Key > sortedToc() const
The table of contents (the keys) in sorted order.
scalar value() const noexcept
The value (const access)
void transfer(List< T > &list)
Output to file stream, using an OSstream.
label size() const noexcept
Number of entries.
virtual bool read()
Re-read model coefficients if they have changed.
static autoPtr< Time > New()
Construct (dummy) Time - no functionObjects or libraries.
T & first()
Return the first element of the list.
void size(const label n)
Older name for setAddressableSize.
T & last()
Return the last element of the list.
static fileName envGlobalPath()
Global case (directory) from environment variable.
A bitSet stores bits (elements with only two states) in packed internal format and supports a variety...
void flip()
Invert all bits in the addressable region.
labelList sortedToc() const
The indices of the on bits as a sorted labelList.
void set(const bitSet &bitset)
Set specified bits from another bitset.
Database for solution data, solver performance and other reduced data.
A list of keyword definitions, which are a keyword followed by a number of values (eg,...
dictionary subOrEmptyDict(const word &keyword, enum keyType::option matchOpt=keyType::REGEX, const bool mandatory=false) const
bool readIfPresent(const word &keyword, T &val, enum keyType::option matchOpt=keyType::REGEX) const
static tmp< complexField > realTransform1D(const scalarField &field)
Transform real-value data.
A class for handling file names.
static word outputPrefix
Directory prefix.
fileName baseFileDir() const
Return the base directory for output.
An instant of time. Contains the time value and name. Uses Foam::Time when formatting the name.
Base class for noise models.
tmp< scalarField > uniformFrequencies(const scalar deltaT, const bool check) const
label findStartTimeIndex(const instantList &allTimes, const scalar startTime) const
Find and return start time index.
scalar RAf(const scalar f) const
A weighting function.
tmp< scalarField > octaves(const scalarField &data, const scalarField &f, const labelUList &freqBandIDs) const
Generate octave data.
tmp< scalarField > RMSmeanPf(const scalarField &p) const
scalar gainC(const scalar f) const
C weighting as gain in dB.
tmp< Foam::scalarField > PSD(const scalarField &PSDf) const
PSD [dB/Hz].
tmp< scalarField > meanPf(const scalarField &p) const
Return the multi-window mean fft of the complete pressure data [Pa].
scalar gainA(const scalar f) const
A weighting as gain in dB.
static void setOctaveBands(const scalarField &f, const scalar fLower, const scalar fUpper, const scalar octave, labelList &fBandIDs, scalarField &fCentre)
virtual bool read(const dictionary &dict)
Read from dictionary.
scalar RDf(const scalar f) const
D weighting function.
tmp< scalarField > Pf(const scalarField &p) const
Return the fft of the given pressure data.
scalar RCf(const scalar f) const
C weighting function.
tmp< scalarField > SPL(const scalarField &Prms2, const scalar f) const
SPL [dB].
planInfo planInfo_
Plan information for FFTW.
static const Enum< weightingType > weightingTypeNames_
scalar checkUniformTimeStep(const scalarList ×) const
Check and return uniform time step.
void readWriteOption(const dictionary &dict, const word &lookup, bool &option) const
Helper function to read write options and provide info feedback.
scalar gainD(const scalar f) const
D weighting as gain in dB.
tmp< scalarField > PSDf(const scalarField &p, const scalar deltaT) const
void writeWeightings() const
Helper function to check weightings.
void cleanFFTW()
Clean up the FFTW.
bool validateBounds(const scalarList &p) const
Return true if all pressure data is within min/max bounds.
scalar gainB(const scalar f) const
B weighting as gain in dB.
scalar RBf(const scalar f) const
B weighting function.
Lookup type of boundary radiation properties.
splitCell * master() const
A class for managing temporary objects.
A class for handling words, derived from Foam::string.
#define defineTypeNameAndDebug(Type, DebugSwitch)
Define the typeName and debug information.
#define FatalIOErrorInFunction(ios)
Report an error message using Foam::FatalIOError.
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
#define WarningInFunction
Report a warning using Foam::Warning.
tmp< scalarField > safeLog10(const scalarField &fld)
dimensionedSymmTensor sqr(const dimensionedVector &dv)
messageStream Info
Information stream (stdout output on master, null elsewhere)
static void check(const int retVal, const char *what)
dimensionedScalar log10(const dimensionedScalar &ds)
dimensioned< Type > average(const DimensionedField< Type, GeoMesh > &df)
void readFields(const typename GeoFieldType::Mesh &mesh, const IOobjectList &objects, const wordHashSet &selectedFields, LIFOStack< regIOobject * > &storedObjects)
Read the selected GeometricFields of the templated type.
fileName::Type type(const fileName &name, const bool followLink=true)
Return the file type: DIRECTORY or FILE, normally following symbolic links.
dimensionedScalar pow(const dimensionedScalar &ds, const dimensionedScalar &expt)
Ostream & endl(Ostream &os)
Add newline and flush stream.
dimensionedScalar sqrt(const dimensionedScalar &ds)
dimensioned< typename typeOfMag< Type >::type > mag(const dimensioned< Type > &dt)
FlatOutput::OutputAdaptor< Container, Delimiters > flatOutput(const Container &obj, Delimiters delim)
Global flatOutput() function with specified output delimiters.
errorManip< error > abort(error &err)
static constexpr const zero Zero
Global zero (0)
word name(const expressions::valueTypeCode typeCode)
A word representation of a valueTypeCode. Empty for INVALID.
prefixOSstream Pout
OSstream wrapped stdout (std::cout) with parallel prefix.
errorManipArg< error, int > exit(error &err, const int errNo=1)
constexpr char nl
The newline '\n' character (0x0a)
#define defineRunTimeSelectionTable(baseType, argNames)
Define run-time selection table.
#define forAll(list, i)
Loop across all elements in list.
const Vector< label > N(dict.get< Vector< label > >("N"))