From b2f002ba457bc9849178198d8ae721755d6c5e33 Mon Sep 17 00:00:00 2001
From: Andrew Heather <a.heather@opencfd.co.uk>
Date: Thu, 16 Nov 2017 18:05:56 +0000
Subject: [PATCH] ENH: FFT - updated for fft of real data
---
src/randomProcesses/fft/fft.C | 50 +++++++++++++++++++
src/randomProcesses/fft/fft.H | 14 +++++-
src/randomProcesses/noise/noiseFFT/noiseFFT.C | 29 +----------
.../windowModel/windowModelTemplates.C | 4 ++
4 files changed, 68 insertions(+), 29 deletions(-)
diff --git a/src/randomProcesses/fft/fft.C b/src/randomProcesses/fft/fft.C
index fd44c0109df..608fd2fd3ab 100644
--- a/src/randomProcesses/fft/fft.C
+++ b/src/randomProcesses/fft/fft.C
@@ -34,6 +34,56 @@ namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+Foam::tmp<Foam::complexField> fft::realTransform1D(const scalarField& field)
+{
+ // Copy of input field for use by fftw
+ // - fftw requires non-const access to input and output...
+ scalarField fftInOut(field);
+
+ const label n = fftInOut.size();
+ const label nBy2 = n/2;
+
+ // Using real to half-complex fftw 'kind'
+ fftw_plan plan = fftw_plan_r2r_1d
+ (
+ n,
+ fftInOut.data(),
+ fftInOut.data(),
+ FFTW_R2HC,
+ FFTW_ESTIMATE
+ );
+
+ fftw_execute(plan);
+
+ fftw_destroy_plan(plan);
+
+ // field[0] = DC component
+ tmp<complexField> tresult(new complexField(nBy2 + 1));
+ complexField& result = tresult.ref();
+
+ result[0].Re() = fftInOut[0];
+ result[nBy2].Re() = fftInOut[nBy2];
+ for (label i = 1; i < nBy2; ++i)
+ {
+ result[i].Re() = fftInOut[i];
+ result[i].Im() = fftInOut[n - i];
+ }
+
+ return tresult;
+}
+
+
+Foam::tmp<Foam::complexField> fft::realTransform1D
+(
+ const tmp<scalarField>& tfield
+)
+{
+ tmp<complexField> tresult = realTransform1D(tfield());
+ tfield.clear();
+ return tresult;
+}
+
+
void fft::transform
(
complexField& field,
diff --git a/src/randomProcesses/fft/fft.H b/src/randomProcesses/fft/fft.H
index 11cc07575e9..7b62e44757e 100644
--- a/src/randomProcesses/fft/fft.H
+++ b/src/randomProcesses/fft/fft.H
@@ -44,7 +44,6 @@ SourceFiles
#define fft_H
#include "complexFields.H"
-#include "UList.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -63,6 +62,19 @@ public:
};
+ //- Transform real-value data
+ // - uses the fftw real to half-complex method
+ // - result size is field.size()/2 + 1
+ static tmp<complexField> realTransform1D(const scalarField& field);
+
+
+ //- Transform real-value data
+ // - uses the fftw real to half-complex method
+ // - result size is field.size()/2 + 1
+ static tmp<complexField> realTransform1D(const tmp<scalarField>& field);
+
+
+ //- Transform compex-value data
static void transform
(
complexField& field,
diff --git a/src/randomProcesses/noise/noiseFFT/noiseFFT.C b/src/randomProcesses/noise/noiseFFT/noiseFFT.C
index e82d25b90fb..57c46d7c0c7 100644
--- a/src/randomProcesses/noise/noiseFFT/noiseFFT.C
+++ b/src/randomProcesses/noise/noiseFFT/noiseFFT.C
@@ -26,7 +26,6 @@ License
#include "noiseFFT.H"
#include "IFstream.H"
#include "DynamicList.H"
-#include "SubField.H"
#include "mathematicalConstants.H"
#include "HashSet.H"
#include "fft.H"
@@ -240,33 +239,7 @@ Foam::tmp<Foam::scalarField> Foam::noiseFFT::Pf
const tmp<scalarField>& tpn
) const
{
- // Calculate the 2-sided fft
- // Note: result not scaled
- tmp<scalarField> tPn2
- (
- mag
- (
- fft::reverseTransform
- (
- ReComplexField(tpn),
- List<int>(1, tpn().size())
- )
- )
- );
-
- tpn.clear();
-
- // Only storing the positive half
- // Note: storing (N/2) values, DC component at position (0)
- tmp<scalarField> tPn
- (
- new scalarField
- (
- scalarField::subField(tPn2(), tPn2().size()/2 + 1)
- )
- );
-
- return tPn;
+ return mag(fft::realTransform1D(tpn));
}
diff --git a/src/randomProcesses/windowModels/windowModel/windowModelTemplates.C b/src/randomProcesses/windowModels/windowModel/windowModelTemplates.C
index 9d3e478c043..90f9971381f 100644
--- a/src/randomProcesses/windowModels/windowModel/windowModelTemplates.C
+++ b/src/randomProcesses/windowModels/windowModel/windowModelTemplates.C
@@ -23,6 +23,10 @@ License
\*---------------------------------------------------------------------------*/
+#include "SubField.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
template<class Type>
Foam::tmp<Foam::Field<Type>> Foam::windowModel::apply
(
--
GitLab