2D Plimpton FFTDriver when using 3D blocks

[kyriakosTapinou]

I'd like to get thoughts on the most efficient way to use the parallel FFT package in Athena++. I've seen the fft.cpp pgen but this doesn't answer all my questions.

My use case is summarised below but I am basically solving an ODE in Fourier space (this is a correction step for electric field driving of flux emergence).

1. Take a time step
2. Extract all block data on the lower z-boundary and store relevant data into a 2D array representing the lower z-boundary
3. Perform FFT, make some calculations, perform inverse FFT with MPI enabled
4. Store result in meshblock data to later distribute to each meshblock on the boundary

I think I have a rough method worked out but I don't think it is very efficient. Its below.

1. Define mesh class pointers to be used with the FFTDriver objects
2. Initialise FFTDriver objects in Mesh::Initialize
3. From within Mesh::UserWorkInLoop --> Run the data harvesting from meshblocks, FFT, and inverse FFT, and allocation to reuser_meshblocks

Eventually I will set this up as internal functions like the Turbulence FFT.

There seems to be a weird disconnect in setting up the FFTW plans and the data load/retrieval where the latter is based off the underlying block but not the dimensions of the plan. I will need to change the load/retrieval by the look of it.

1. Does anyone have efficiency advice for the FFT workflow?

2. Does anyone have any advice for the 2D FFT ambiguity when using Athena++ with 3D?

[changgoo]

Sounds like you are doing 2D FFT on boundaries, which look like a task required for the James method for a self-gravity solver (https://ui.adsabs.harvard.edu/abs/2019ApJS..241...24M/abstract). @sanghyukmoon has a version for that gravity solver, which may have a functionality you need. I think he didn't use the Athena's FFT wrapper for that task though.

[kyriakosTapinou]

@changgoo Thank you for your response. Yes, it is a 2D FFT on the physical boundary. The work from @sanghyukmoon sounds relevant.

That code or a snippet/explanation would be very helpful.

[sanghyukmoon]

Due to the problems you mentioned, I had to use FFTW and Plimpton's remap routines directly rather than utilizing AthenaFFT wrappers. To do the 2D FFT starting from block-decomposed data, I create communicators involving only those ranks at the boundary (https://github.com/sanghyukmoon/athena/blob/0d63120317b5bd95328b9501e1a017ae24cd7c2b/src/gravity/james_gravity.cpp#L346-L358). The rank and size within these boundary communicators (which is different from those associated with MPI_COMM_WORLD) are then used to assign the start/end indices for y- and z-pencil decomposition (for x-boundaries) needed to perform 1D FFTs along each dimension (https://github.com/sanghyukmoon/athena/blob/0d63120317b5bd95328b9501e1a017ae24cd7c2b/src/gravity/james_gravity.cpp#L379-L407). I hope this can provide a starting point for your application.

[sanghyukmoon]

For 2D FFT at the surface?

[kyriakosTapinou]

Yes, 2D FFT at a boundary.

[sanghyukmoon]

I haven't tried serial solve (which would require gathering all data in the surface to one process), so I don't have empirical answer. Very roughly, because (serial surface FFT) / (parallel volume FFT) ~ N^2 / n^3 ~ m^2 / n (neglecting log n factors; N is the total number of cells per dimension, n is the number of cells per MeshBlock per dimension, and m is the number of MeshBlocks per dimension), I guess the surface term can dominate in the usual problem sizes.

[kyriakosTapinou]

Ok, thanks.

I've had a dive into your code and it is very helpful. Is the main reason you have used the DCT and DST to enforce your boundary conditions? I have periodic boundaries so using a regular fft is fine (I think) and I think it won't break the pencil decomposition. What do you think?

[sanghyukmoon]

Yes, the use of DCT/DST is a very particular choice for the James algorithm to enforce boundary conditions. If your BC is periodic the regular complex FFT is entirely fine (in fact if your input is real, than use of r2c routines (https://www.fftw.org/fftw3_doc/Real_002ddata-DFTs.html) will cut the data size in half, thus in principle bring at least factor of two performance gain).