zhyang-dev/openfoam
1
0
Fork 0
mirror of https://develop.openfoam.com/Development/openfoam.git synced 2025-11-28 03:28:01 +00:00
Code Issues Actions Packages Projects Releases Wiki Activity

ENH: fft - use fftw_malloc to create in/out arrays. Fixes #1095

Browse Source
This commit is contained in:
Andrew Heather
2018-11-09 09:20:50 +00:00
parent 584d6d066a
commit 2a3e2b2cfa
1 changed files with 13 additions and 8 deletions
Download Patch File Download Diff File Expand all files Collapse all files

21
src/randomProcesses/fft/fft.C
View File

@ -180,10 +180,12 @@ void Foam::fft::transform
) )
{ {
// Copy field into fftw containers // Copy field into fftw containers
// NB: this is not fully compliant, since array sizing is non-constexpr.
// However, cannot instantiate List with fftw_complex
const label N = field.size(); const label N = field.size();
fftw_complex in[N], out[N];
fftw_complex* inPtr =
static_cast<fftw_complex*>(fftw_malloc(sizeof(fftw_complex)*N));
fftw_complex* outPtr =
static_cast<fftw_complex*>(fftw_malloc(sizeof(fftw_complex)*N));
// If reverse transform : renumber before transform // If reverse transform : renumber before transform
if (dir == REVERSE_TRANSFORM) if (dir == REVERSE_TRANSFORM)
@ -193,8 +195,8 @@ void Foam::fft::transform
forAll(field, i) forAll(field, i)
{ {
in[i][0] = field[i].Re(); inPtr[i][0] = field[i].Re();
in[i][1] = field[i].Im(); inPtr[i][1] = field[i].Im();
} }
// Create the plan // Create the plan
@ -208,19 +210,22 @@ void Foam::fft::transform
// Generic 1..3-D plan // Generic 1..3-D plan
const label rank = nn.size(); const label rank = nn.size();
fftw_plan plan = fftw_plan plan =
fftw_plan_dft(rank, nn.begin(), in, out, dir, FFTW_ESTIMATE); fftw_plan_dft(rank, nn.begin(), inPtr, outPtr, dir, FFTW_ESTIMATE);
// Compute the FFT // Compute the FFT
fftw_execute(plan); fftw_execute(plan);
forAll(field, i) forAll(field, i)
{ {
field[i].Re() = out[i][0]; field[i].Re() = outPtr[i][0];
field[i].Im() = out[i][1]; field[i].Im() = outPtr[i][1];
} }
fftw_destroy_plan(plan); fftw_destroy_plan(plan);
fftw_free(inPtr);
fftw_free(outPtr);
// If forward transform : renumber after transform // If forward transform : renumber after transform
if (dir == FORWARD_TRANSFORM) if (dir == FORWARD_TRANSFORM)
{ {