From 0efe8b8ffa80d7043007f96368eef403548c355c Mon Sep 17 00:00:00 2001 From: Andrew Heather Date: Wed, 29 Jun 2016 20:37:39 +0100 Subject: [PATCH] ENH: FFT - replaced Numerical Recipes-based FFT by the FFTW library --- src/randomProcesses/Allwmake | 18 +++ src/randomProcesses/fft/fft.C | 174 ++++++++------------------ src/randomProcesses/fft/fft.H | 11 +- src/randomProcesses/fft/fftRenumber.C | 124 ------------------ src/randomProcesses/fft/fftRenumber.H | 81 ------------ 5 files changed, 73 insertions(+), 335 deletions(-) create mode 100755 src/randomProcesses/Allwmake delete mode 100644 src/randomProcesses/fft/fftRenumber.C delete mode 100644 src/randomProcesses/fft/fftRenumber.H 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/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/fft/fftRenumber.H b/src/randomProcesses/fft/fftRenumber.H deleted file mode 100644 index fa2e50c580..0000000000 --- a/src/randomProcesses/fft/fftRenumber.H +++ /dev/null @@ -1,81 +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 . - -Global - fftRenumber - -Description - Multi-dimensional renumbering used in the Numerical Recipes - fft routine. - -SourceFiles - fftRenumber.C - -\*---------------------------------------------------------------------------*/ - -#ifndef fftRenumber_H -#define fftRenumber_H - -#include "complex.H" -#include "List.H" -#include "labelList.H" - -// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // - -namespace Foam -{ - -// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // - -// Recursively evaluate the indexing necessary to do the folding of the fft -// data. We recurse until we have the indexing ncessary for the folding in all -// directions. -void fftRenumberRecurse -( - List& data, - List& renumData, - const labelList& nn, - label nnprod, - label ii, - label l1, - label l2 -); - - -// Fold the n-d data array to get the fft components in the right places. -void fftRenumber -( - List& data, - const labelList& nn -); - - -// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // - -} // End namespace Foam - -// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // - -#endif - -// ************************************************************************* //