# Spectral derivatives D2Z transform

Good morning everyone,

I have started using cuda fortran for two weeks now, and I must recognize I am stuck with a strange problem.

I want to compute the derivative of a 2D field using spectral derivatives according to the following program :

program example

use cudafor
use precision_m
use cufft_m

implicit none

integer,          parameter :: N = 8
integer,          parameter :: Nh = N/2 + 1
double precision, parameter :: Ns = real(N*N)
double precision, parameter :: pi = 3.141592653589793238d0

double precision, dimension(N,N)  :: signal, dx_signal, dy_signal
double precision, dimension(N)    :: kx
double precision, dimension(Nh)   :: ky

double precision, dimension(N,N),  device :: d_signal, d_dx_signal, d_dy_signal
double complex,   dimension(N,Nh), device :: d_fft,    d_dx_fft,    d_dy_fft
double precision, dimension(N),    device :: d_kx
double precision, dimension(Nh),   device :: d_ky

integer :: i, j, plan
double complex :: temp

!--- Initialisation

do i = 1, N
do j = 1, N
signal(i,j) = cos( 2.*pi * ( 0.*(i-1) + 1.*(j-1) ) / N)
end do
end do

kx(1:N/2)   = (/(dfloat(i-1)   * 2. * pi, i=1,N/2,1)/)
kx(N/2+1:N) = (/(dfloat(i-1-N) * 2. * pi, i=N/2+1,N,1)/)
ky(1:Nh)    = (/(dfloat(i-1)   * 2. * pi, i=1,Nh,1)/)

! Move to device
d_signal = signal
d_kx = kx
d_ky = ky

!--- Computations

! We compute the Fourier transform
call cufftPlan2D(plan, N, N, CUFFT_D2Z)
call cufftExecD2Z(plan, d_signal, d_fft)
call cufftDestroy(plan)

! We compute the x derivative
!$cuf kernel do(2) <<<*,*>>> do i = 1, N do j = 1, Nh temp = cmplx(0., 1.) * d_kx(i) temp = temp * d_fft(i,j) d_dx_fft(i,j) = temp / Ns end do end do ! We compute the y derivative !$ cuf kernel do(2) <<<*,*>>>
do i = 1, N
do j = 1, Nh
temp = cmplx(0., 1.) * d_ky(j)
temp = temp * d_fft(i,j)
d_dy_fft(i,j) = temp / Ns
end do
end do

! We compute the Fourier inverse transform
call cufftPlan2D(plan, N, N, CUFFT_Z2D)
call cufftExecZ2D(plan, d_dx_fft, d_dx_signal)
call cufftExecZ2D(plan, d_dy_fft, d_dy_signal)
call cufftDestroy(plan)

! Move to host
dx_signal = d_dx_signal
dy_signal = d_dy_signal

!--- Print dx_signal and dy_signal

write(*,*) "dx_signal"
do i = 1, N
write(*,*) ( dx_signal(i,j), j=1,N )
end do

write(*,*) "dy_signal"
do i = 1, N
write(*,*) ( dy_signal(i,j), j=1,N )
end do

end program example


What is strange is that the derivative according to x should be zero, but this is not what we observe. On the contrary, we can almost find the derivative according to y.

I compile using nvfortran and the following command:
nvfortran -fast -Mpreprocess -Mcuda -c precision_m.f90 -o precision_m.o -lcufft
nvfortran -fast -Mpreprocess -Mcuda -c cufft_m.f90 -o cufft_m.o -lcufft
nvfortran -fast -Mpreprocess -Mcuda -c example.f90 -o example.o -lcufft
nvfortran -fast -Mpreprocess -Mcuda -O3 -o example example.o cufft_m.o precision_m.o -lcufft

And I execute using:
./example

Thank you very much for your help

B.

fortran is column major for array storage. CUFFT expects row-major inputs (that is, the first dimension for CUFFT transforms will have elements that are adjacent to each other in memory). Have you accounted for that?

The complex domain (i.e. output) from a D2Z transform will exhibit hermitian symmetry. This has specific implications for the data that will actually be stored - it is not identical to the data storage for an equivalent Z2Z forward transform. Have you accounted for that?

The complex domain must have the same kind of hermitian symmetry for the inverse (Z2D) transform (to be valid). Are you certain that your complex domain arithmetic has that implication?