diff --git a/applications/solvers/DNS/dnsFoam/Allwmake b/applications/solvers/DNS/dnsFoam/Allwmake new file mode 100755 index 0000000000..759041a131 --- /dev/null +++ b/applications/solvers/DNS/dnsFoam/Allwmake @@ -0,0 +1,14 @@ +#!/bin/sh +cd ${0%/*} || exit 1 # run from this directory + +if [ -f "$FFTW_ARCH_PATH/include/fftw3.h" ] || \ + [ "${FFTW_ARCH_PATH##*-}" = system -a -f "/usr/include/fftw3.h" ] +then + wmake +else + echo + echo "Skipping dnsFoam solver (no FFTW)" + echo +fi + +#------------------------------------------------------------------------------ diff --git a/applications/utilities/postProcessing/noise/Allwmake b/applications/utilities/postProcessing/noise/Allwmake new file mode 100755 index 0000000000..84c6a2ac82 --- /dev/null +++ b/applications/utilities/postProcessing/noise/Allwmake @@ -0,0 +1,14 @@ +#!/bin/sh +cd ${0%/*} || exit 1 # run from this directory + +if [ -f "$FFTW_ARCH_PATH/include/fftw3.h" ] || \ + [ "${FFTW_ARCH_PATH##*-}" = system -a -f "/usr/include/fftw3.h" ] +then + wmake +else + echo + echo "Skipping noise utility (no FFTW)" + echo +fi + +#------------------------------------------------------------------------------ diff --git a/applications/utilities/postProcessing/noise/Make/options b/applications/utilities/postProcessing/noise/Make/options index b29b939580..2fc476e741 100644 --- a/applications/utilities/postProcessing/noise/Make/options +++ b/applications/utilities/postProcessing/noise/Make/options @@ -1,6 +1,7 @@ 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 \ diff --git a/applications/utilities/postProcessing/noise/noise.C b/applications/utilities/postProcessing/noise/noise.C index 1b0ca07fec..cf5794075f 100644 --- a/applications/utilities/postProcessing/noise/noise.C +++ b/applications/utilities/postProcessing/noise/noise.C @@ -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,103 @@ 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; + noiseModel surfaceNoise; // pointNoise - pressureData + 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 pData("pressure", dict, "Data"); - - // time history data - const scalarField t(pData.x()); - - // pressure data - const scalarField p(pData.y()); - - if (t.size() < N) + fileName dictName(runTime.system()/"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; + IOdictionary dict + ( + IOobject + ( + dictName, + runTime, + IOobject::MUST_READ + ) + ); - - 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); + autoPtr model(noiseModel::New(dict)); + model->calculate(); Info<< nl << "End\n" << endl; diff --git a/applications/utilities/postProcessing/noise/noiseDict b/applications/utilities/postProcessing/noise/noiseDict new file mode 100644 index 0000000000..e983cabc12 --- /dev/null +++ b/applications/utilities/postProcessing/noise/noiseDict @@ -0,0 +1,136 @@ +/*--------------------------------*- C++ -*----------------------------------*\ +| ========= | | +| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox | +| \\ / O peration | Version: plus | +| \\ / A nd | Web: www.OpenFOAM.com | +| \\/ M anipulation | | +\*---------------------------------------------------------------------------*/ +FoamFile +{ + version 2.0; + format ascii; + class dictionary; + location "system"; + object noiseDict; +} +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +noiseModel surfaceNoise; + +surfaceNoiseCoeffs +{ + windowModel Hanning; + + HanningCoeffs + { + // Window overlap percentage + overlapPercent 50; + symmetric yes; + extended yes; + + // Optional number of windows, default = all available + // nWindow 1; + } + +/* + windowModel uniform; + + uniformCoeffs + { + // Window overlap percentage + overlapPercent 50; + + value 1; + + // Optional number of windows, default = all available + // nWindow 1; + } +*/ + + + // Input file + inputFile "postProcessing/faceSource1/surface/patch_motorBike_rider-helmet%65/patch_motorBike_rider-helmet%65.case"; + + // Surface reader + reader ensight; + + // Surface writer + writer ensight; + + // Collate times for ensight output - ensures geometry is only written once + writeOptions + { + ensight + { + collateTimes 1; + } + } + + // Reference density (to convert from kinematic to static pressure) + rhoRef 1.205; + + // Number of samples in sampling window + // Must be a power of 2, default = 2^16 (=65536) + N 4096; // 8192; // 4096; + + // Lower frequency limit, default = 25Hz + //fl 25; + + // Upper frequency limit, default = 10kHz + fu 15000; + + // Start time, default = 0s + //startTime 0; + + // Write interval for FFT data, default = 1 +// fftWriteInterval 100; +} + +pointNoiseCoeffs +{ + csvFileData + { + fileName "pressureData"; + nHeaderLine 1; + refColumn 0; + componentColumns (1); + separator " "; + mergeSeparators yes; + } + + HanningCoeffs + { + // Window overlap percentage + overlapPercent 50; + symmetric yes; + extended yes; + + // Optional number of windows, default = all available + //nWindow 5; + } + + // Graph format, default = raw + graphFormat raw; + + // Reference density (to convert from kinematic to static pressure) + rhoRef 1.2; + + // Number of samples in sampling window + // Must be a power of 2, default = 2^16 (=65536) + N 4096; + + // Lower frequency limit, default = 25Hz + //fl 25; + + // Upper frequency limit, default = 10kHz + //fu 10000; + + // Start time, default = 0s + //startTime 0; + + // Write interval for FFT data, default = 1 + fftWriteInterval 100; +} + + +// ************************************************************************* // diff --git a/applications/utilities/preProcessing/boxTurb/Allwmake b/applications/utilities/preProcessing/boxTurb/Allwmake new file mode 100755 index 0000000000..46d04eebf5 --- /dev/null +++ b/applications/utilities/preProcessing/boxTurb/Allwmake @@ -0,0 +1,14 @@ +#!/bin/sh +cd ${0%/*} || exit 1 # run from this directory + +if [ -f "$FFTW_ARCH_PATH/include/fftw3.h" ] || \ + [ "${FFTW_ARCH_PATH##*-}" = system -a -f "/usr/include/fftw3.h" ] +then + wmake +else + echo + echo "Skipping boxTurb utility (no FFTW)" + echo +fi + +#------------------------------------------------------------------------------ diff --git a/etc/cshrc b/etc/cshrc index 099662b296..24154f8b86 100644 --- a/etc/cshrc +++ b/etc/cshrc @@ -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 diff --git a/src/conversion/Make/files b/src/conversion/Make/files index 7c5a0aa91c..b4a38476a1 100644 --- a/src/conversion/Make/files +++ b/src/conversion/Make/files @@ -1,5 +1,6 @@ ensight/file/ensightFile.C ensight/file/ensightGeoFile.C +ensight/readFile/ensightReadFile.C ensight/part/ensightPart.C ensight/part/ensightPartIO.C ensight/part/ensightPartCells.C diff --git a/src/conversion/ensight/readFile/ensightReadFile.C b/src/conversion/ensight/readFile/ensightReadFile.C new file mode 100644 index 0000000000..f73d95199e --- /dev/null +++ b/src/conversion/ensight/readFile/ensightReadFile.C @@ -0,0 +1,157 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 2016 OpenCFD Ltd. + \\/ M anipulation | +------------------------------------------------------------------------------- +License + This file is part of OpenFOAM. + + OpenFOAM is free software: you can redistribute it and/or modify it + under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + OpenFOAM is distributed in the hope that it will be useful, but WITHOUT + ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or + FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License + for more details. + + You should have received a copy of the GNU General Public License + along with OpenFOAM. If not, see . + +\*---------------------------------------------------------------------------*/ + +#include "ensightReadFile.H" +#include + +// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // + +Foam::ensightReadFile::ensightReadFile +( + const fileName& pathname, + IOstream::streamFormat format +) +: + IFstream(pathname, format) +{} + + +// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * // + +Foam::ensightReadFile::~ensightReadFile() +{} + + +// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // + +Foam::Istream& Foam::ensightReadFile::read +( + char* buf, + std::streamsize count +) +{ + stdStream().read(buf, count); + return *this; +} + + +Foam::Istream& Foam::ensightReadFile::read(string& value) +{ + if (format() == IOstream::BINARY) + { + char buf[80]; + + read(reinterpret_cast(buf), sizeof(buf)); + + string strBuf(value); + + const size_t iEnd = strBuf.find('\0', 0); + if (iEnd == string::npos) + { + value = buf; + } + else + { + value = strBuf.substr(0, iEnd - 1); + } + } + else + { + value = ""; + while (value.empty() && !eof()) + { + getLine(value); + } + } + + return *this; +} + + +Foam::Istream& Foam::ensightReadFile::read(label& value) +{ + int ivalue; + + if (format() == IOstream::BINARY) + { + read + ( + reinterpret_cast(&ivalue), + sizeof(ivalue) + ); + } + else + { + stdStream() >> ivalue; + } + + value = ivalue; + return *this; +} + + +Foam::Istream& Foam::ensightReadFile::read(scalar& value) +{ + float fvalue; + + if (format() == IOstream::BINARY) + { + read + ( + reinterpret_cast(&fvalue), + sizeof(fvalue) + ); + + value = fvalue; + } + else + { + stdStream() >> value; + } + + return *this; +} + + +Foam::Istream& Foam::ensightReadFile::readKeyword(string& key) +{ + read(key); + return *this; +} + + +Foam::Istream& Foam::ensightReadFile::readBinaryHeader() +{ + if (format() == IOstream::BINARY) + { + string buffer; + read(buffer); + } + + return *this; +} + + +// ************************************************************************* // diff --git a/src/conversion/ensight/readFile/ensightReadFile.H b/src/conversion/ensight/readFile/ensightReadFile.H new file mode 100644 index 0000000000..aa6723352b --- /dev/null +++ b/src/conversion/ensight/readFile/ensightReadFile.H @@ -0,0 +1,110 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 2016 OpenCFD Ltd. + \\/ M anipulation | +------------------------------------------------------------------------------- +License + This file is part of OpenFOAM. + + OpenFOAM is free software: you can redistribute it and/or modify it + under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + OpenFOAM is distributed in the hope that it will be useful, but WITHOUT + ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or + FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License + for more details. + + You should have received a copy of the GNU General Public License + along with OpenFOAM. If not, see . + +Class + Foam::ensightReadFile + +Description + Ensight output with specialized read() for strings, integers and floats. + Correctly handles binary read as well. + +\*---------------------------------------------------------------------------*/ + +#ifndef ensightReadFile_H +#define ensightReadFile_H + +#include "IFstream.H" +#include "IOstream.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + +/*---------------------------------------------------------------------------*\ + Class ensightReadFile Declaration +\*---------------------------------------------------------------------------*/ + +class ensightReadFile +: + public IFstream +{ + // Private Member Functions + + //- Disallow default bitwise assignment + void operator=(const ensightReadFile&); + + //- Disallow default copy constructor + ensightReadFile(const ensightReadFile&); + + +public: + + // Constructors + + //- Construct from pathname + ensightReadFile + ( + const fileName& pathname, + IOstream::streamFormat format=IOstream::BINARY + ); + + + //- Destructor + ~ensightReadFile(); + + + // Output + + //- Inherit read from Istream + using Istream::read; + + //- Binary read + virtual Istream& read(char* buf, std::streamsize count); + + //- Read string as "%80s" or as binary + Istream& read(string& value); + + //- Read integer as "%10d" or as binary + Istream& read(label& value); + + //- Read float as "%12.5e" or as binary + Istream& read(scalar& value); + + //- Read element keyword + virtual Istream& readKeyword(string& key); + + //- Read "C Binary" for binary files (eg, geometry/measured) + Istream& readBinaryHeader(); +}; + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace Foam + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#endif + +// ************************************************************************* // diff --git a/src/randomProcesses/Allwmake b/src/randomProcesses/Allwmake new file mode 100755 index 0000000000..e1ce66ba6a --- /dev/null +++ b/src/randomProcesses/Allwmake @@ -0,0 +1,18 @@ +#!/bin/sh +cd ${0%/*} || exit 1 # run from this directory + +# Parse arguments for library compilation +targetType=libso +. $WM_PROJECT_DIR/wmake/scripts/AllwmakeParseArguments + +if [ -f "$FFTW_ARCH_PATH/include/fftw3.h" ] || \ + [ "${FFTW_ARCH_PATH##*-}" = system -a -f "/usr/include/fftw3.h" ] +then + wmake $targetType +else + echo + echo "Skipping randomProcesses library (no FFTW)" + echo +fi + +#------------------------------------------------------------------------------ diff --git a/src/randomProcesses/Make/files b/src/randomProcesses/Make/files index 90cb360357..11138bd932 100644 --- a/src/randomProcesses/Make/files +++ b/src/randomProcesses/Make/files @@ -1,25 +1,31 @@ -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 +$(windowModels)/uniform/uniform.C + LIB = $(FOAM_LIBBIN)/librandomProcesses diff --git a/src/randomProcesses/Make/options b/src/randomProcesses/Make/options index 71b7873964..e1eddd51d2 100644 --- a/src/randomProcesses/Make/options +++ b/src/randomProcesses/Make/options @@ -1,5 +1,11 @@ EXE_INC = \ - -I$(LIB_SRC)/finiteVolume/lnInclude + -I$(FFTW_ARCH_PATH)/include \ + -I$(LIB_SRC)/finiteVolume/lnInclude \ + -I$(LIB_SRC)/sampling/lnInclude \ + -I$(LIB_SRC)/surfMesh/lnInclude LIB_LIBS = \ - -lfiniteVolume + -L$(FFTW_ARCH_PATH)/lib$(WM_COMPILER_LIB_ARCH) -lfftw3 \ + -lfiniteVolume \ + -lsampling \ + -lsurfMesh diff --git a/src/randomProcesses/fft/fft.C b/src/randomProcesses/fft/fft.C index f859b37f52..e8194da8dd 100644 --- a/src/randomProcesses/fft/fft.C +++ b/src/randomProcesses/fft/fft.C @@ -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. @@ -24,7 +24,8 @@ License \*---------------------------------------------------------------------------*/ #include "fft.H" -#include "fftRenumber.H" + +#include // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // @@ -33,16 +34,11 @@ namespace Foam // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // -#define SWAP(a,b) tempr=(a); (a)=(b); (b)=tempr -#define TWOPI 6.28318530717959 - -// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // - void fft::transform ( complexField& field, const labelList& nn, - transformDirection isign + transformDirection dir ) { forAll(nn, idim) @@ -59,128 +55,58 @@ void fft::transform } } - const label ndim = nn.size(); - - label i1, i2, i3, i2rev, i3rev, ip1, ip2, ip3, ifp1, ifp2; - label ibit, k1, k2, n, nprev, nrem, idim; - scalar tempi, tempr; - scalar theta, wi, wpi, wpr, wr, wtemp; - scalar* data = reinterpret_cast(field.begin()) - 1; - - - // if inverse transform : renumber before transform - - if (isign == REVERSE_TRANSFORM) - { - fftRenumber(field, nn); - } - - - label ntot = 1; - forAll(nn, idim) - { - ntot *= nn[idim]; - } - - - nprev = 1; - - for (idim=ndim; idim>=1; idim--) - { - n = nn[idim-1]; - nrem = ntot/(n*nprev); - ip1 = nprev << 1; - ip2 = ip1*n; - ip3 = ip2*nrem; - i2rev = 1; - - for (i2=1; i2<=ip2; i2+=ip1) - { - if (i2 < i2rev) - { - for (i1=i2; i1<=i2 + ip1 - 2; i1+=2) - { - for (i3=i1; i3<=ip3; i3+=ip2) - { - i3rev = i2rev + i3 - i2; - SWAP(data[i3], data[i3rev]); - SWAP(data[i3 + 1], data[i3rev + 1]); - } - } - } - - ibit = ip2 >> 1; - while (ibit >= ip1 && i2rev > ibit) - { - i2rev -= ibit; - ibit >>= 1; - } - - i2rev += ibit; - } - - ifp1 = ip1; - - while (ifp1 < ip2) - { - ifp2 = ifp1 << 1; - theta = isign*TWOPI/(ifp2/ip1); - wtemp = sin(0.5*theta); - wpr = -2.0*wtemp*wtemp; - wpi = sin(theta); - wr = 1.0; - wi = 0.0; - - for (i3 = 1; i3 <= ifp1; i3 += ip1) - { - for (i1 = i3; i1 <= i3 + ip1 - 2; i1 += 2) - { - for (i2 = i1; i2 <= ip3; i2 += ifp2) - { - k1 = i2; - k2 = k1 + ifp1; - tempr = scalar(wr*data[k2]) - scalar(wi*data[k2 + 1]); - tempi = scalar(wr*data[k2 + 1]) + scalar(wi*data[k2]); - data[k2] = data[k1] - tempr; - data[k2 + 1] = data[k1 + 1] - tempi; - data[k1] += tempr; - data[k1 + 1] += tempi; - } - } - - wr = (wtemp = wr)*wpr - wi*wpi + wr; - wi = wi*wpr + wtemp*wpi + wi; - } - ifp1 = ifp2; - } - nprev *= n; - } - - - // if forward transform : renumber after transform - - if (isign == FORWARD_TRANSFORM) - { - fftRenumber(field, nn); - } - - - // scale result (symmetric scale both forward and inverse transform) - - scalar recRootN = 1.0/sqrt(scalar(ntot)); + // Copy field into fftw containers + label N = field.size(); + fftw_complex in[N], out[N]; forAll(field, i) { - field[i] *= recRootN; + in[i][0] = field[i].Re(); + in[i][1] = field[i].Im(); } + + // Create the plan + // FFTW_FORWARD = -1 + // FFTW_BACKWARD = 1 + + // 1-D plan + // fftw_plan plan = + // fftw_plan_dft_1d(N, in, out, FFTW_FORWARD, FFTW_ESTIMATE); + + // Generic 1..3-D plan + label rank = nn.size(); + fftw_plan plan = + fftw_plan_dft(rank, nn.begin(), in, out, dir, FFTW_ESTIMATE); + + // Compute the FFT + fftw_execute(plan); + + forAll(field, i) + { + field[i].Re() = out[i][0]; + field[i].Im() = out[i][1]; + } + + fftw_destroy_plan(plan); + +/* + fftw_real in[N], out[N]; + // Create a plan for real data input + // - generates 1-sided FFT + // - direction given by FFTW_R2HC or FFTW_HC2R. + // - in[N] is the array of real input val ues + // - out[N] is the array of N/2 real valuea followed by N/2 complex values + // - 0 component = DC component + fftw_plan plan = fftw_plan_r2r_2d(N, in, out, FFTW_R2HC, FFTW_ESTIMATE) + + // Compute the FFT + fftw_execute(plan); + + fftw_destroy_plan(plan); +*/ } -// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // - -#undef SWAP -#undef TWOPI - // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // tmp fft::forwardTransform diff --git a/src/randomProcesses/fft/fft.H b/src/randomProcesses/fft/fft.H index f176261eb5..dc5a0c2c69 100644 --- a/src/randomProcesses/fft/fft.H +++ b/src/randomProcesses/fft/fft.H @@ -3,7 +3,7 @@ \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | \\ / A nd | Copyright (C) 2011 OpenFOAM Foundation - \\/ M anipulation | + \\/ M anipulation | Copyright (C) 2016 OpenCFD Ltd. ------------------------------------------------------------------------------- License This file is part of OpenFOAM. @@ -25,11 +25,10 @@ Class Foam::fft Description - Fast fourier transform derived from the Numerical - Recipes in C routine. + Fast fourier transform using the fftw library. The complex transform field is returned in the field supplied. The - direction of transform is supplied as an argument (1 = forward, -1 = + direction of transform is supplied as an argument (-1 = forward, 1 = reverse). The dimensionality and organisation of the array of values in space is supplied in the nn indexing array. @@ -56,8 +55,8 @@ public: enum transformDirection { - FORWARD_TRANSFORM = 1, - REVERSE_TRANSFORM = -1 + FORWARD_TRANSFORM = -1, + REVERSE_TRANSFORM = 1 }; diff --git a/src/randomProcesses/fft/fftRenumber.C b/src/randomProcesses/fft/fftRenumber.C deleted file mode 100644 index 07f8b238be..0000000000 --- a/src/randomProcesses/fft/fftRenumber.C +++ /dev/null @@ -1,124 +0,0 @@ -/*---------------------------------------------------------------------------*\ - ========= | - \\ / F ield | OpenFOAM: The Open Source CFD Toolbox - \\ / O peration | - \\ / A nd | Copyright (C) 2011-2016 OpenFOAM Foundation - \\/ M anipulation | -------------------------------------------------------------------------------- -License - This file is part of OpenFOAM. - - OpenFOAM is free software: you can redistribute it and/or modify it - under the terms of the GNU General Public License as published by - the Free Software Foundation, either version 3 of the License, or - (at your option) any later version. - - OpenFOAM is distributed in the hope that it will be useful, but WITHOUT - ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or - FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License - for more details. - - You should have received a copy of the GNU General Public License - along with OpenFOAM. If not, see . - -Description - Multi-dimensional renumbering used in the Numerical Recipes - fft routine. This version is recursive, so works in n-d : - determined by the length of array nn - -\*---------------------------------------------------------------------------*/ - -#include "fftRenumber.H" - -// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // - -void Foam::fftRenumberRecurse -( - List& data, - List& renumData, - const labelList& nn, - label nnprod, - label ii, - label l1, - label l2 -) -{ - if (ii == nn.size()) - { - // we've worked out the renumbering scheme. Now copy - // the components across - - data[l1] = complex(renumData[l2].Re(),renumData[l2].Im()); - } - else - { - // do another level of folding. First work out the - // multiplicative value of the index - - nnprod /= nn[ii]; - label i_1(0); - - for (label i=0; i& data, - const labelList& nn -) -{ - List renumData(data); - - label nnprod(1); - forAll(nn, i) - { - nnprod *= nn[i]; - } - - label ii(0), l1(0), l2(0); - - fftRenumberRecurse - ( - data, - renumData, - nn, - nnprod, - ii, - l1, - l2 - ); -} - - -// ************************************************************************* // diff --git a/src/randomProcesses/noise/noiseFFT.C b/src/randomProcesses/noise/noiseFFT.C deleted file mode 100644 index 0f09cdf817..0000000000 --- a/src/randomProcesses/noise/noiseFFT.C +++ /dev/null @@ -1,454 +0,0 @@ -/*---------------------------------------------------------------------------*\ - ========= | - \\ / F ield | OpenFOAM: The Open Source CFD Toolbox - \\ / O peration | - \\ / A nd | Copyright (C) 2011-2016 OpenFOAM Foundation - \\/ M anipulation | -------------------------------------------------------------------------------- -License - This file is part of OpenFOAM. - - OpenFOAM is free software: you can redistribute it and/or modify it - under the terms of the GNU General Public License as published by - the Free Software Foundation, either version 3 of the License, or - (at your option) any later version. - - OpenFOAM is distributed in the hope that it will be useful, but WITHOUT - ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or - FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License - for more details. - - You should have received a copy of the GNU General Public License - along with OpenFOAM. If not, see . - -\*---------------------------------------------------------------------------*/ - -#include "noiseFFT.H" -#include "IFstream.H" -#include "DynamicList.H" -#include "fft.H" -#include "SubField.H" -#include "mathematicalConstants.H" - -// * * * * * * * * * * * * * * Static Member Data * * * * * * * * * * * * * // - -Foam::scalar Foam::noiseFFT::p0 = 2e-5; - - -// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // - -Foam::noiseFFT::noiseFFT -( - const scalar deltat, - const scalarField& pressure -) -: - scalarField(pressure), - deltat_(deltat) -{} - - -Foam::noiseFFT::noiseFFT(const fileName& pFileName, const label skip) -: - scalarField(), - deltat_(0.0) -{ - // Construct pressure data file - IFstream pFile(pFileName); - - // Check pFile stream is OK - if (!pFile.good()) - { - FatalErrorInFunction - << "Cannot read file " << pFileName - << exit(FatalError); - } - - if (skip) - { - scalar dummyt, dummyp; - - for (label i=0; i> dummyt; - - if (!pFile.good() || pFile.eof()) - { - FatalErrorInFunction - << "Number of points in file " << pFileName - << " is less than the number to be skipped = " << skip - << exit(FatalError); - } - - pFile >> dummyp; - } - } - - scalar t = 0, T = 0; - DynamicList pData(100000); - label i = 0; - - while (!(pFile >> t).eof()) - { - T = t; - pFile >> pData(i++); - } - - deltat_ = T/pData.size(); - - this->transfer(pData); -} - - -// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // - -Foam::graph Foam::noiseFFT::pt() const -{ - scalarField t(size()); - forAll(t, i) - { - t[i] = i*deltat_; - } - - return graph - ( - "p(t)", - "t [s]", - "p(t) [Pa]", - t, - *this - ); -} - - -Foam::tmp 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 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::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::noiseFFT::Pf -( - const tmp& tpn -) const -{ - tmp tPn2 - ( - mag - ( - fft::reverseTransform - ( - ReComplexField(tpn), - labelList(1, tpn().size()) - ) - ) - ); - - tpn.clear(); - - tmp tPn - ( - new scalarField - ( - scalarField::subField(tPn2(), tPn2().size()/2) - ) - ); - scalarField& Pn = tPn.ref(); - - Pn *= 2.0/sqrt(scalar(tPn2().size())); - Pn[0] /= 2.0; - - return tPn; -} - - -Foam::graph Foam::noiseFFT::meanPf -( - const label N, - const label nw -) 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); - } - - scalarField MeanPf(N/2, 0.0); - - scalarField Hwf(Hanning(N)); - - for (label wi=0; wi size()) - { - FatalErrorInFunction - << "Requested window is outside set of data" << endl - << "number of data = " << size() << endl - << "size of window = " << N << endl - << "window = " << nw - << exit(FatalError); - } - - scalarField RMSMeanPf(N/2, 0.0); - - scalarField Hwf(Hanning(N)); - - for (label wi=0; wi fU + 1) break; - - if (f[i] >= fu) - { - fm[j] = fmi; - ldelta[j] = 10*log10(ldelta[j]); - - j++; - - fl = fu; - fu *= fratio; - } - - ldelta[j] += pow(10, Lf[i]/10.0); - } - - fm.setSize(j); - ldelta.setSize(j); - - return graph - ( - "Ldelta", - "fm [Hz]", - "Ldelta [dB]", - fm, - ldelta - ); -} - - -Foam::graph Foam::noiseFFT::Pdelta -( - const graph& gPf, - const scalar f1, - const scalar fU -) const -{ - const scalarField& f = gPf.x(); - const scalarField& Pf = gPf.y(); - - scalarField pdelta(Pf.size(), 0.0); - scalarField fm(pdelta.size()); - - scalar fratio = cbrt(2.0); - scalar deltaf = 1.0/(2*Pf.size()*deltat_); - - scalar fl = f1/sqrt(fratio); - scalar fu = fratio*fl; - - label istart = label(fl/deltaf + 1.0 - SMALL); - label j = 0; - - for (label i = istart; i fU + 1) break; - - if (f[i] >= fu) - { - fm[j] = fmi; - pdelta[j] = sqrt((2.0/3.0)*pdelta[j]); - - j++; - - fl = fu; - fu *= fratio; - } - - pdelta[j] += sqr(Pf[i]); - } - - fm.setSize(j); - pdelta.setSize(j); - - return graph - ( - "Pdelta", - "fm [Hz]", - "Pdelta [dB]", - fm, - pdelta - ); -} - - -Foam::scalar Foam::noiseFFT::Lsum(const graph& gLf) const -{ - const scalarField& Lf = gLf.y(); - - scalar lsum = 0.0; - - forAll(Lf, i) - { - lsum += pow(10, Lf[i]/10.0); - } - - lsum = 10*log10(lsum); - - return lsum; -} - - -Foam::scalar Foam::noiseFFT::dbToPa(const scalar db) const -{ - return p0*pow(10.0, db/20.0); -} - - -Foam::tmp Foam::noiseFFT::dbToPa -( - const tmp& db -) const -{ - return p0*pow(10.0, db/20.0); -} - - -// ************************************************************************* // diff --git a/src/randomProcesses/noise/noiseFFT/noiseFFT.C b/src/randomProcesses/noise/noiseFFT/noiseFFT.C new file mode 100644 index 0000000000..abaf199514 --- /dev/null +++ b/src/randomProcesses/noise/noiseFFT/noiseFFT.C @@ -0,0 +1,474 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 2011-2015 OpenFOAM Foundation + \\/ M anipulation | Copyright (C) 2016 OpenCFD Ltd. +------------------------------------------------------------------------------- +License + This file is part of OpenFOAM. + + OpenFOAM is free software: you can redistribute it and/or modify it + under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + OpenFOAM is distributed in the hope that it will be useful, but WITHOUT + ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or + FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License + for more details. + + You should have received a copy of the GNU General Public License + along with OpenFOAM. If not, see . + +\*---------------------------------------------------------------------------*/ + +#include "noiseFFT.H" +#include "IFstream.H" +#include "DynamicList.H" +#include "SubField.H" +#include "mathematicalConstants.H" +#include "HashSet.H" +#include "fft.H" + +using namespace Foam::constant; + +// * * * * * * * * * * * * * * Static Member Data * * * * * * * * * * * * * // + +Foam::scalar Foam::noiseFFT::p0 = 2e-5; + + +Foam::tmp Foam::noiseFFT::frequencies +( + const label N, + const scalar deltaT +) +{ + tmp tf(new scalarField(N/2, 0)); + scalarField& f = tf.ref(); + + scalar deltaf = 1.0/(N*deltaT); + forAll(f, i) + { + f[i] = i*deltaf; + } + + return tf; +} + + +Foam::tmp Foam::noiseFFT::PSD(const scalarField& PSDf) +{ + return 10*log10(PSDf/sqr(p0)); +} + + +Foam::tmp Foam::noiseFFT::SPL(const scalarField& Prms2) +{ + return 10*log10(Prms2/sqr(p0)); +} + + +void Foam::noiseFFT::octaveBandInfo +( + const scalarField& f, + const scalar fLower, + const scalar fUpper, + const scalar octave, + labelList& fBandIDs, + scalarField& fCentre +) +{ + // 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^(0.5*bandI/octave)) + + // Initial (lowest centre frequency) + scalar fTest = 15.625; + + const scalar fRatio = pow(2, 1.0/octave); + const scalar fRatioL2C = pow(2, 0.5/octave); + + // IDs of band IDs + labelHashSet bandIDs(f.size()); + + // Centre frequencies + DynamicList fc; + + // Convert to lower band limit + fTest /= fRatioL2C; + + forAll(f, i) + { + if (f[i] >= fTest) + { + // Advance band if appropriate + while (f[i] > fTest) + { + fTest *= fRatio; + } + fTest /= fRatio; + + if (bandIDs.insert(i)) + { + // Also store (next) centre frequency + fc.append(fTest*fRatioL2C); + } + fTest *= fRatio; + + if (fTest > fUpper) + { + break; + } + } + } + + fBandIDs = bandIDs.toc(); + + // Remove the last centre frequency (beyond upper frequency limit) + fc.remove(); + + fCentre.transfer(fc); +} + + +// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // + +Foam::noiseFFT::noiseFFT +( + const scalar deltaT, + const scalarField& pressure +) +: + scalarField(pressure), + deltaT_(deltaT) +{ + scalarField& p = *this; + p -= average(p); +} + + +Foam::noiseFFT::noiseFFT(const fileName& pFileName, const label skip) +: + scalarField(), + deltaT_(0.0) +{ + // Construct pressure data file + IFstream pFile(pFileName); + + // Check pFile stream is OK + if (!pFile.good()) + { + FatalErrorInFunction + << "Cannot read file " << pFileName + << exit(FatalError); + } + + if (skip) + { + scalar dummyt, dummyp; + + for (label i = 0; i < skip; i++) + { + pFile >> dummyt; + + if (!pFile.good() || pFile.eof()) + { + FatalErrorInFunction + << "Number of points in file " << pFileName + << " is less than the number to be skipped = " << skip + << exit(FatalError); + } + + pFile >> dummyp; + } + } + + scalar t = 0, T0 = 0, T1 = 0; + DynamicList pData(100000); + label i = 0; + + while (!(pFile >> t).eof()) + { + if (i == 0) + { + T0 = t; + } + + T1 = t; + pFile >> pData(i++); + } + + // Note: assumes fixed time step + deltaT_ = (T1 - T0)/pData.size(); + + this->transfer(pData); + + scalarField& p = *this; + p -= average(p); +} + + +// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // + +Foam::graph Foam::noiseFFT::pt() const +{ + scalarField t(size()); + forAll(t, i) + { + t[i] = i*deltaT_; + } + + return graph + ( + "p(t)", + "t [s]", + "p(t) [Pa]", + t, + *this + ); +} + + +Foam::tmp Foam::noiseFFT::Pf +( + const tmp& tpn +) const +{ + // Calculate the 2-sided fft + // Note: result not scaled + tmp tPn2 + ( + mag + ( + fft::reverseTransform + ( + ReComplexField(tpn), + labelList(1, tpn().size()) + ) + ) + ); + + tpn.clear(); + + // Only storing the positive half + // Note: storing (N/2) values, DC component at position (0) + tmp tPn + ( + new scalarField + ( + scalarField::subField(tPn2(), tPn2().size()/2 + 1) + ) + ); + + return tPn; +} + + +Foam::graph Foam::noiseFFT::meanPf(const windowModel& window) const +{ + const label N = window.nSamples(); + const label nWindow = window.nWindow(); + + scalarField meanPf(N/2 + 1, 0.0); + + for (label windowI = 0; windowI < nWindow; ++windowI) + { + meanPf += Pf(window.apply(*this, windowI)); + } + + meanPf /= scalar(nWindow); + + scalar deltaf = 1.0/(N*deltaT_); + scalarField f(meanPf.size()); + forAll(f, i) + { + f[i] = i*deltaf; + } + + return graph + ( + "P(f)", + "f [Hz]", + "P(f) [Pa]", + f, + meanPf + ); +} + + +Foam::graph Foam::noiseFFT::RMSmeanPf(const windowModel& window) const +{ + const label N = window.nSamples(); + const label nWindow = window.nWindow(); + + scalarField RMSMeanPf(N/2 + 1, 0.0); + for (label windowI = 0; windowI < nWindow; ++windowI) + { + RMSMeanPf += sqr(Pf(window.apply(*this, windowI))); + } + + RMSMeanPf = sqrt(RMSMeanPf/scalar(nWindow))/scalar(N); + + scalar deltaf = 1.0/(N*deltaT_); + scalarField f(RMSMeanPf.size()); + forAll(f, i) + { + f[i] = i*deltaf; + } + + return graph + ( + "Prms(f)", + "f [Hz]", + "Prms(f) [Pa]", + f, + RMSMeanPf + ); +} + + +Foam::graph Foam::noiseFFT::PSDf(const windowModel& window) const +{ + const label N = window.nSamples(); + const label nWindow = window.nWindow(); + + scalarField psd(N/2 + 1, 0.0); + + for (label windowI = 0; windowI < nWindow; ++windowI) + { + psd += sqr(Pf(window.apply(*this, windowI))); + } + + scalar deltaf = 1.0/(N*deltaT_); + scalar fs = 1.0/deltaT_; + psd /= scalar(nWindow)*fs*N; + + // Scaling due to use of 1-sided FFT + // Note: do not scale DC component + psd *= 2; + psd.first() /= 2; + psd.last() /= 2; + + scalarField f(psd.size()); + + if (0) // if (debug) + { + Pout<< "Average PSD: " << average(psd) << endl; + } + + forAll(f, i) + { + f[i] = i*deltaf; + } + + return graph + ( + "PSD(f)", + "f [Hz]", + "PSD(f) [PaPa_Hz]", + f, + psd + ); +} + + +Foam::graph Foam::noiseFFT::PSD(const graph& gPSDf) const +{ + return graph + ( + "PSD(f)", + "f [Hz]", + "PSD_dB(f) [dB_Hz]", + gPSDf.x(), + 10*log10(gPSDf.y()/sqr(p0)) + ); +} + + +Foam::graph Foam::noiseFFT::octaves +( + const graph& g, + const labelList& freqBandIDs, + bool integrate +) const +{ + if (freqBandIDs.size() < 2) + { + WarningInFunction + << "Octave frequency bands are not defined " + << "- skipping octaves calculation" + << endl; + + return graph + ( + "octave", + "x", + "y", + scalarField(), + scalarField() + ); + } + + const scalarField& f = g.x(); + const scalarField& data = g.y(); + + scalarField octData(freqBandIDs.size() - 1, 0.0); + scalarField fm(freqBandIDs.size() -1, 0.0); + + for (label bandI = 0; bandI < freqBandIDs.size() - 1; bandI++) + { + label fb0 = freqBandIDs[bandI]; + label fb1 = freqBandIDs[bandI+1]; + fm[bandI] = f[fb0]; + + if (fb0 == fb1) continue; + + if (integrate) + { + for (label freqI = fb0; freqI < fb1; freqI++) + { + label f0 = f[freqI]; + label f1 = f[freqI + 1]; + scalar dataAve = 0.5*(data[freqI] + data[freqI + 1]); + octData[bandI] += dataAve*(f1 - f0); + } + } + else + { + for (label freqI = fb0; freqI < fb1; freqI++) + { + octData[bandI] += data[freqI]; + } + } + } + + return graph + ( + "octaves(f)", + "fm [Hz]", + "octave data", + fm, + octData + ); +} + + +Foam::scalar Foam::noiseFFT::dbToPa(const scalar db) const +{ + return p0*pow(10.0, db/20.0); +} + + +Foam::tmp Foam::noiseFFT::dbToPa +( + const tmp& db +) const +{ + return p0*pow(10.0, db/20.0); +} + + +// ************************************************************************* // diff --git a/src/randomProcesses/noise/noiseFFT.H b/src/randomProcesses/noise/noiseFFT/noiseFFT.H similarity index 53% rename from src/randomProcesses/noise/noiseFFT.H rename to src/randomProcesses/noise/noiseFFT/noiseFFT.H index 085702964c..6439a2f429 100644 --- a/src/randomProcesses/noise/noiseFFT.H +++ b/src/randomProcesses/noise/noiseFFT/noiseFFT.H @@ -2,8 +2,8 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2011 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. @@ -25,11 +25,26 @@ Class Foam::noiseFFT Description - FFT of the pressure field + Performs FFT of pressure field to generate noise data. + + General functionality: + - Pf: fft of the pressure data + - meanPf: multi-window mean fft + - RMSmeanPf: multi-window RMS mean fft + - PSDf: multi-window power spectral density (PSD) in frequency domain + - PSD: power spectral density in dB/Hz + - SPL: sound pressure level in dB + + Octave-based data: + - PSD spectrum + - SPL spectrum SourceFiles noiseFFT.C +SeeAlso + windowModel.H + \*---------------------------------------------------------------------------*/ #ifndef noiseFFT_H @@ -37,6 +52,7 @@ SourceFiles #include "scalarField.H" #include "graph.H" +#include "windowModel.H" // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // @@ -44,7 +60,7 @@ namespace Foam { /*---------------------------------------------------------------------------*\ - Class noiseFFT Declaration + Class noiseFFT Declaration \*---------------------------------------------------------------------------*/ class noiseFFT @@ -54,7 +70,18 @@ class noiseFFT // Private data //- Time spacing of the raw data - scalar deltat_; + scalar deltaT_; + + struct octaveBandInfo + { + label octave; + + // IDs of bin boundaries in pressure data + labelList binIDs; + + // Centre frequencies for each bin + scalarField centreFreq; + }; public: @@ -81,37 +108,61 @@ public: // Member Functions - //- Return the graph of p(t) + //- Return the FFT frequencies + static tmp frequencies + ( + const label N, + const scalar deltaT + ); + + //- Return the PSD [dB/Hz] + // Input PSD in [Pa^2/Hz] + static tmp PSD(const scalarField& PSDf); + + //- Return the SPL [dB] + // Input P(rms)^2 in [Pa^2] + static tmp SPL(const scalarField& Prms2); + + //- Return a list of the frequency indices wrt f field that + // correspond to the bands limits for a given octave + static void octaveBandInfo + ( + const scalarField& f, + const scalar fLower, + const scalar fUpper, + const scalar octave, + labelList& fBandIDs, + scalarField& fCentre + ); + + //- Return the graph of pressure as a function of time graph pt() const; - //- Return the nth window - tmp window(const label N, const label n) const; - - //- Return the Hanning window function - tmp Hanning(const label N) const; - //- Return the fft of the given pressure data tmp Pf(const tmp& pn) const; - //- Return the multi-window mean fft of the complete pressure data - graph meanPf(const label N, const label nw) const; + //- Return the multi-window mean fft of the complete pressure data [Pa] + graph meanPf(const windowModel& window) const; - //- Return the multi-window RMS mean fft of the complete pressure data - graph RMSmeanPf(const label N, const label nw) const; + //- Return the multi-window RMS mean fft of the complete pressure + // data [Pa] + graph RMSmeanPf(const windowModel& window) const; - //- Return the narrow-band PFL (pressure-fluctuation level) spectrum - graph Lf(const graph& gPf) const; + //- Return the multi-window PSD (power spectral density) of the complete + // pressure data [Pa^2/Hz] + graph PSDf(const windowModel& window) const; - //- Return the one-third-octave-band PFL spectrum - // starting at octave with mean frequency f1 - graph Ldelta(const graph& gLf, const scalar f1, const scalar fU) const; + //- Return the PSD [dB/Hz] + // Takes PSD in [Pa^2/Hz] + graph PSD(const graph& gPSDf) const; - //- Return the one-third-octave-band pressure spectrum - // starting at octave with mean frequency f1 - graph Pdelta(const graph& gLf, const scalar f1, const scalar fU) const; - - //- Return the total PFL as the sum of Lf over all frequencies - scalar Lsum(const graph& gLf) const; + //- Generate octave data + graph octaves + ( + const graph& g, + const labelList& freqBandIDs, + bool integrate + ) const; //- Convert the db into Pa scalar dbToPa(const scalar db) const; diff --git a/src/randomProcesses/noise/noiseModels/noiseModel/noiseModel.C b/src/randomProcesses/noise/noiseModels/noiseModel/noiseModel.C new file mode 100644 index 0000000000..145f992725 --- /dev/null +++ b/src/randomProcesses/noise/noiseModels/noiseModel/noiseModel.C @@ -0,0 +1,151 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 2015-2016 OpenCFD Ltd. + \\/ M anipulation | +------------------------------------------------------------------------------- +License + This file is part of OpenFOAM. + + OpenFOAM is free software: you can redistribute it and/or modify it + under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + OpenFOAM is distributed in the hope that it will be useful, but WITHOUT + ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or + FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License + for more details. + + You should have received a copy of the GNU General Public License + along with OpenFOAM. If not, see . + +\*---------------------------------------------------------------------------*/ + +#include "noiseModel.H" + +// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // + +namespace Foam +{ + defineTypeNameAndDebug(noiseModel, 0); + defineRunTimeSelectionTable(noiseModel, dictionary); +} + +// * * * * * * * * * * * * * Protected Member Functions * * * * * * * * * * // + +Foam::scalar Foam::noiseModel::checkUniformTimeStep +( + const scalarList& times +) const +{ + scalar deltaT = -1.0; + + if (times.size() > 1) + { + for (label i = 1; i < times.size(); i++) + { + scalar dT = times[i] - times[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; +} + + +Foam::label Foam::noiseModel::findStartTimeIndex +( + const instantList& allTimes, + const scalar startTime +) const +{ + forAll(allTimes, timeI) + { + const instant& i = allTimes[timeI]; + + if (i.value() >= startTime) + { + return timeI; + } + } + + return 0; +} + + +// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // + +Foam::noiseModel::noiseModel(const dictionary& dict) +: + dict_(dict), + rhoRef_(dict.lookupOrDefault("rhoRef", 1)), + nSamples_(dict.lookupOrDefault