search

icl-utk-edu / heffte

BSD 3-Clause "New" or "Revised" License
20 stars 15 forks source link

Method to calculate modes? #29

Closed pgrete closed 10 months ago

pgrete commented 1 year ago

Is there a method to calculate the modes/wavenumbers in Fourier space? The background is that I'd like to calculate power spectra and thus need to bin the data in Fourier space based on their wavenumber.

pgrete commented 1 year ago

ping @mkstoyanov (as you were just active on the repo yesterday)

mkstoyanov commented 1 year ago

The box assumption is that the index of the data and the wave number will match, i.e.,

input[ k * plane_stride  + j * line_stride + i] // box entry (i, j, k)
output[ k * plane_stride  + j * line_stride + i] // wave vector (i, j, k)

Note the orientation (order) of the box, it indicates which is the leading (fast changing), middle and trailing (slowest changing) dimension.

In the distributed case, you have to account for the fact that each MPI rank contains a local box and the box and index start and end offset.

There is no method in heffte to give you the wavenumber, but that should be enough to write your own method.

pgrete commented 1 year ago

Thanks for the pointer. I'm not sure I fully understand. Assuming that i,j,k.plane_stride,line_stride are all > 0, how would I get to a wave vector of, say, (-3, -2, 0)?

PS: I also assume that we might use different vocabulary

mkstoyanov commented 10 months ago

@pgrete For some reason I am not getting notifications on this one. Sorry for the delay.

I'm not sure what you're calling "wavenumber" and what would correspond to negative numbers.

pgrete commented 10 months ago

Thanks for the follow-up @mkstoyanov What I meant is similar to the numpy fftfreq output where the second half of the frequencies (wavenumbers) are negative.

mkstoyanov commented 10 months ago

heFFTe doesn't have such method, but it seems like a simple map from the index to the wavenumber.

Basically, what they call n in numpy would be (i, j, k) in the box index of heFFTe.

pgrete commented 10 months ago

Yes, so that was probably implicitly my question. How are indices mapped to wavenumber in heFFTe (as different convention exit in different libraries). Does it follow numpy convention, i.e., for even number of modes the first half of the wavemodes go from 0 to n/2 -1 and the second half goes from -n/2 to -1 ?

mkstoyanov commented 10 months ago

heFFTe doesn't decide those, it's a matter of the backend. However, every backend that we support will generate identical wave-numbers, which is the direct definition of DFT. Libraries that shuffle the wave-numbers by default would be in the minority.

See this about FFTW, which is the library that sets most of the conventions when it comes to FFTs

https://www.uni-muenster.de/Physik.TP/archive/fileadmin/lehre/NumMethoden/SoSe10/Skript/Ordening.pdf

pgrete commented 10 months ago

Thanks for clarifying. I was just double checking as I already learned the hard way in the past that different conventions are being used by different libs.

Feel free to close the issue (or keep it open if some callback function would be desirable to be available in heFFTe).

mkstoyanov commented 10 months ago

That's actually good to know, I would never expect that there can be a difference in the different libs. It's something to keep an eye out. Thanks for bringing it up.

I won't add a method to heFFTe (at least for now) but I should probably document that we are using the same conventions as FFTW.