ENH: FFT - replaced Numerical Recipes-based FFT by the FFTW library

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
+
+#------------------------------------------------------------------------------

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
+    \\  /    A nd           | Copyright (C) 2011-2015 OpenFOAM Foundation
-     \\/     M anipulation  |
+     \\/     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 <fftw3.h>
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -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,127 +55,57 @@ void fft::transform
         }
     }
 
-    const label ndim = nn.size();
+    // Copy field into fftw containers
-
+    label N = field.size();
-    label i1, i2, i3, i2rev, i3rev, ip1, ip2, ip3, ifp1, ifp2;
+    fftw_complex in[N], out[N];
-    label ibit, k1, k2, n, nprev, nrem, idim;
-    scalar tempi, tempr;
-    scalar theta, wi, wpi, wpr, wr, wtemp;
-    scalar* data = reinterpret_cast<scalar*>(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));
 
     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
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 

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
+    Fast fourier transform using the fftw library.
-    Recipes in C routine.
 
     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,
+        FORWARD_TRANSFORM = -1,
-        REVERSE_TRANSFORM = -1
+        REVERSE_TRANSFORM = 1
     };
 
 

src/randomProcesses/fft/fftRenumber.C

@@ -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 <http://www.gnu.org/licenses/>.
-
-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<complex>& data,
-    List<complex>& 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<nn[ii]; i++)
-        {
-            // now evaluate the indices (both from array 1 and to
-            // array 2). These get multiplied by nnprod to (cumulatively)
-            // find the real position in the list corresponding to
-            // this set of indices.
-
-            if (i<nn[ii]/2)
-            {
-                i_1 = i + nn[ii]/2;
-            }
-            else
-            {
-                i_1 = i - nn[ii]/2;
-            }
-
-
-            // go to the next level of recursion.
-
-            fftRenumberRecurse
-            (
-                data,
-                renumData,
-                nn,
-                nnprod,
-                ii+1,
-                l1+i*nnprod,
-                l2+i_1*nnprod
-            );
-        }
-    }
-}
-
-
-void Foam::fftRenumber
-(
-    List<complex>& data,
-    const labelList& nn
-)
-{
-    List<complex> 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
-    );
-}
-
-
-// ************************************************************************* //

src/randomProcesses/fft/fftRenumber.H

@@ -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 <http://www.gnu.org/licenses/>.
-
-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<complex>& data,
-    List<complex>& 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<complex>& data,
-    const labelList& nn
-);
-
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-} // End namespace Foam
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-#endif
-
-// ************************************************************************* //