diff --git a/src/ExB_shear_flow_mod.F90 b/src/ExB_shear_flow_mod.F90
index 6c43d8ccff4e7fc3d1f2146a0189a6eb98609199..b165cb73f5c0c0cff437b9113ac9031f872fbd2e 100644
--- a/src/ExB_shear_flow_mod.F90
+++ b/src/ExB_shear_flow_mod.F90
@@ -74,7 +74,7 @@ CONTAINS
! Setup nonlinear factor
ALLOCATE( ExB_NL_factor(Nx,local_nky))
- ALLOCATE(inv_ExB_NL_factor(Ny/2+2,local_nx))
+ ALLOCATE(inv_ExB_NL_factor(Ny/2+1,local_nx))
ExB_NL_factor = 1._xp
inv_ExB_NL_factor = 1._xp
@@ -231,7 +231,7 @@ CONTAINS
ENDDO
ENDDO
! ... and the inverse
- DO iky = 1,Ny/2+2 ! WARNING: Global indices ky loop
+ DO iky = 1,Ny/2+1 ! WARNING: Global indices ky loop
! for readability
J_xp = REAL(iky-1,xp)
IF(J_xp .GT. 0._xp) THEN
@@ -248,9 +248,6 @@ CONTAINS
! inv_ExB_NL_factor(iky,ix) = EXP(imagu*sky_ExB_full(iky)*xval)
ENDDO
ENDDO
- ! Cancel the additional point
- inv_ExB_NL_factor(Ny/2+1,:) = 0._xp
- inv_ExB_NL_factor(Ny/2+2,:) = 0._xp
END SUBROUTINE Update_nonlinear_ExB_factors
END MODULE ExB_shear_flow
\ No newline at end of file
diff --git a/src/fourier_mod.F90 b/src/fourier_mod.F90
index f8785a89748999b4592f211238ef45a106faa8f1..50b44d5a281269052060d07ae7544ca3808819f4 100644
--- a/src/fourier_mod.F90
+++ b/src/fourier_mod.F90
@@ -29,15 +29,27 @@ MODULE fourier
complex(c_xp_c), pointer, PUBLIC :: cmpx_data_f(:,:), cmpx_data_g(:,:), bracket_sum_c(:,:)
REAL(xp), PUBLIC :: inv_Nx_, inv_Ny_
!! 1D fft specific variables (full fortran interface)
+ ! Single DFT plans
type(C_PTR), PUBLIC :: plan_ky2y_c2r ! transform from ( x,ky) to ( x, y) (complex to real)
type(C_PTR), PUBLIC :: plan_y2ky_r2c ! transform from ( x, y) to ( x,ky) (real to complex)
type(C_PTR), PUBLIC :: plan_kx2x_c2c ! transform from (kx,ky) to ( x,ky) (complex to complex)
type(C_PTR), PUBLIC :: plan_x2kx_c2c ! transform from ( x,ky) to (kx,ky) (complex to complex)
+! COMPLEX(xp), DIMENSION(:), ALLOCATABLE :: f_kx_l,f_x_l,f_ky_g ! aux array for plans
+! REAL(xp), DIMENSION(:), ALLOCATABLE :: f_y_g
+ complex(c_xp_c), ALLOCATABLE, PUBLIC :: f_kx_l(:),f_x_l(:),f_ky_g(:)
+ real (c_xp_r), ALLOCATABLE, PUBLIC :: f_y_g(:)
+ ! Many DFT plans
+ type(C_PTR), PUBLIC :: plan_many_ky2y_c2r ! transform from ( x,ky) to ( x, y) (complex to real)
+ type(C_PTR), PUBLIC :: plan_many_y2ky_r2c ! transform from ( x, y) to ( x,ky) (real to complex)
+ type(C_PTR), PUBLIC :: plan_many_kx2x_c2c ! transform from (kx,ky) to ( x,ky) (complex to complex)
+ type(C_PTR), PUBLIC :: plan_many_x2kx_c2c ! transform from ( x,ky) to (kx,ky) (complex to complex)
COMPLEX(xp), DIMENSION(:,:), ALLOCATABLE :: f_kxky_l ! working arrays
COMPLEX(xp), DIMENSION(:,:), ALLOCATABLE :: f_xky_l ! temp. 1D ifft algo stage (local mpi)
- REAL(xp), DIMENSION(:,:), ALLOCATABLE :: bracket_sum_yx_g ! poisson bracket sum in real space
- COMPLEX(xp), DIMENSION(:,:), ALLOCATABLE :: bracket_sum_kyx_g ! final poisson bracket in complex space
+ REAL(xp), DIMENSION(:,:), ALLOCATABLE :: f_yx_g ! poisson bracket sum in real space
+ COMPLEX(xp), DIMENSION(:,:), ALLOCATABLE :: f_kyx_g ! final poisson bracket in complex space
+ ! to chose between a plan many or 1D plan loops
+ LOGICAL :: PLAN_MANY = .FALSE.
CONTAINS
!******************************************************************************!
!------------- Initialize the grid distribution and plans -------------
@@ -52,13 +64,14 @@ MODULE fourier
NX_ = Nx; NY_ = Ny
inv_Nx_ = 1._xp/REAL(NX_,xp)
NY_halved = NY_/2 + 1
- inv_Ny_ = 1._xp/REAL(2*NY_halved,xp)
+ ! Not clear which one should be the normalization factor
+ inv_Ny_ = 1._xp/REAL(NY_,xp)
! Call FFTW 2D mpi routines to distribute the data and init 2D MPI FFTW plans
CALL fft2D_distr_and_plans(communicator,&
local_nx_ptr,local_nx_ptr_offset,local_nky_ptr,local_nky_ptr_offset)
local_nx_ = local_nx_ptr ! store number of local x in the module
local_nky_ = local_nky_ptr ! store number of local ky in the module
- ! Init 1D MPI FFTW plans for ExB rate correction factor
+ ! Init 1D MPI FFTW plans for ExB rate correction factor
CALL fft1D_plans
! store data distr. in the module for the poisson_bracket function
END SUBROUTINE init_grid_distr_and_plans
@@ -143,8 +156,32 @@ MODULE fourier
integer(C_INTPTR_T) :: inembed, onembed
integer(C_INTPTR_T) :: istride, ostride
integer(C_INTPTR_T) :: idist, odist
-
- !! Plan of the 4 1D many transforms required
+ !! Plans of the 4 single 1D transforms
+ !----------- 1: FFTx and inv through local ky data
+ !(kx,ky) <-> (x,ky), C <-> C, transforms
+ ! in & out
+ ALLOCATE( f_kx_l(NX_))
+ ALLOCATE( f_x_l(NX_))
+#ifdef SINGLE_PRECISION
+ CALL sfftw_plan_dft_1d(plan_kx2x_c2c,NX_,f_kx_l,f_x_l,FFTW_BACKWARD,FFTW_PATIENT)
+ CALL sfftw_plan_dft_1d(plan_x2kx_c2c,NX_,f_x_l,f_kx_l,FFTW_FORWARD, FFTW_PATIENT)
+#else
+ CALL dfftw_plan_dft_1d(plan_kx2x_c2c,NX_,f_kx_l,f_x_l,FFTW_BACKWARD,FFTW_PATIENT)
+ CALL dfftw_plan_dft_1d(plan_x2kx_c2c,NX_,f_x_l,f_kx_l,FFTW_FORWARD, FFTW_PATIENT)
+#endif
+ !----------- 2: FFTy and inv through global ky data
+ !(y,x) <-> (ky,x), R <-> C, transforms (bplan_y in GENE)
+ ! in & out
+ ALLOCATE( f_y_g(NY_))
+ ALLOCATE(f_ky_g(NY_/2+1))
+#ifdef SINGLE_PRECISION
+ CALL sfftw_plan_dft_r2c_1d(plan_y2ky_r2c,NY_,f_y_g,f_ky_g,FFTW_PATIENT)
+ CALL sfftw_plan_dft_c2r_1d(plan_ky2y_c2r,NY_,f_ky_g,f_y_g,FFTW_PATIENT)
+#else
+ CALL dfftw_plan_dft_r2c_1d(plan_y2ky_r2c,NY_,f_y_g,f_ky_g,FFTW_PATIENT)
+ CALL dfftw_plan_dft_c2r_1d(plan_ky2y_c2r,NY_,f_ky_g,f_y_g,FFTW_PATIENT)
+#endif
+ !! Plans of the 4 1D many transforms required
!----------- 1: FFTx and inv through local ky data
!1.1 (kx,ky) -> (x,ky), C -> C, transforms
! in:
@@ -162,12 +199,12 @@ MODULE fourier
istride = 1 ! contiguous data
ostride = 1
#ifdef SINGLE_PRECISION
- CALL sfftw_plan_many_dft(plan_kx2x_c2c, rank, n, howmany,&
+ CALL sfftw_plan_many_dft(plan_many_kx2x_c2c, rank, n, howmany,&
f_kxky_l, inembed, istride, idist,&
f_xky_l, onembed, ostride, odist,&
FFTW_BACKWARD, FFTW_PATIENT)
#else
- CALL dfftw_plan_many_dft(plan_kx2x_c2c, rank, n, howmany,&
+ CALL dfftw_plan_many_dft(plan_many_kx2x_c2c, rank, n, howmany,&
f_kxky_l, inembed, istride, idist,&
f_xky_l, onembed, ostride, odist,&
FFTW_BACKWARD, FFTW_PATIENT)
@@ -186,78 +223,81 @@ MODULE fourier
istride = 1 ! contiguous data
ostride = 1
#ifdef SINGLE_PRECISION
- CALL sfftw_plan_many_dft(plan_x2kx_c2c, rank, n, howmany,&
+ CALL sfftw_plan_many_dft(plan_many_x2kx_c2c, rank, n, howmany,&
f_xky_l, inembed, istride, idist,&
f_kxky_l, onembed, ostride, odist,&
FFTW_FORWARD, FFTW_PATIENT)
#else
- CALL dfftw_plan_many_dft(plan_x2kx_c2c, rank, n, howmany,&
+ CALL dfftw_plan_many_dft(plan_many_x2kx_c2c, rank, n, howmany,&
f_xky_l, inembed, istride, idist,&
f_kxky_l, onembed, ostride, odist,&
FFTW_FORWARD, FFTW_PATIENT)
#endif
!----------- 2: FFTy and inv through global ky data
+ ! (can be improved by inputting directly the padded array of size 2(Ny/2+1))
! 2.1 (y,x) -> (ky,x), R -> C, transforms (bplan_y in GENE)
! in:
- ALLOCATE(bracket_sum_yx_g(2*NY_halved,local_nx_))
+ ALLOCATE(f_yx_g(NY_,local_nx_))
! out:
- ALLOCATE(bracket_sum_kyx_g(NY_halved+1,local_nx_))
+ ALLOCATE(f_kyx_g(NY_/2+1,local_nx_))
! transform parameters
rank = 1 ! 1D transform
- n = 2*NY_halved ! all y
+ n = NY_ ! all y
howmany = local_nx_ ! all x
- inembed = 2*NY_halved ! all y must be used
+ inembed = NY_ ! all y must be used
istride = 1 ! contiguous data
- idist = 2*NY_halved ! distance between two slice to transforms (y row)
- onembed = NY_halved+1
+ idist = NY_ ! distance between two slice to transforms (y row)
+ onembed = NY_/2+1
ostride = 1
- odist = NY_halved+1
+ odist = NY_/2+1
#ifdef SINGLE_PRECISION
- CALL sfftw_plan_many_dft_r2c(plan_y2ky_r2c, rank, n, howmany,&
- bracket_sum_yx_g, inembed, istride, idist,&
- bracket_sum_kyx_g, onembed, ostride, odist,&
+ CALL sfftw_plan_many_dft_r2c(plan_many_y2ky_r2c, rank, n, howmany,&
+ f_yx_g, inembed, istride, idist,&
+ f_kyx_g, onembed, ostride, odist,&
FFTW_PATIENT)
#else
- CALL dfftw_plan_many_dft_r2c(plan_y2ky_r2c, rank, n, howmany,&
- bracket_sum_yx_g, inembed, istride, idist,&
- bracket_sum_kyx_g, onembed, ostride, odist,&
+ CALL dfftw_plan_many_dft_r2c(plan_many_y2ky_r2c, rank, n, howmany,&
+ f_yx_g, inembed, istride, idist,&
+ f_kyx_g, onembed, ostride, odist,&
FFTW_PATIENT)
#endif
!-----------
! 2.2 (x,ky) -> (x,y), C -> R, transforms (fplan_y in GENE)
- ! in: bracket_sum_kyx_g
- ! out: bracket_sum_yx_g
+ ! in: f_kyx_g
+ ! out: f_yx_g
! transform parameters
rank = 1 ! 1D transform
- n = 2*NY_halved ! all ky
+ n = NY_ ! all ky
howmany = local_nx_ ! all x
- inembed = NY_halved+1 ! all ky must be used
+ inembed = NY_/2+1 ! all ky must be used
istride = 1 ! contiguous data
- idist = NY_halved+1 ! distance between two slice to transforms (y row)
- onembed = 2*NY_halved ! to all y
+ idist = NY_/2+1 ! distance between two slice to transforms (y row)
+ onembed = NY_ ! to all y
ostride = 1
- odist = 2*NY_halved
+ odist = NY_
#ifdef SINGLE_PRECISION
- CALL sfftw_plan_many_dft_c2r(plan_ky2y_c2r, rank, n, howmany,&
- bracket_sum_kyx_g, inembed, istride, idist,&
- bracket_sum_yx_g, onembed, ostride, odist,&
+ CALL sfftw_plan_many_dft_c2r(plan_many_ky2y_c2r, rank, n, howmany,&
+ f_kyx_g, inembed, istride, idist,&
+ f_yx_g, onembed, ostride, odist,&
FFTW_PATIENT)
#else
- CALL dfftw_plan_many_dft_c2r(plan_ky2y_c2r, rank, n, howmany,&
- bracket_sum_kyx_g, inembed, istride, idist,&
- bracket_sum_yx_g, onembed, ostride, odist,&
+ CALL dfftw_plan_many_dft_c2r(plan_many_ky2y_c2r, rank, n, howmany,&
+ f_kyx_g, inembed, istride, idist,&
+ f_yx_g, onembed, ostride, odist,&
FFTW_PATIENT)
#endif
- if ((.not. c_associated(plan_ky2y_c2r)) .OR. &
- (.not. c_associated(plan_y2ky_r2c)) .OR. &
- (.not. c_associated(plan_kx2x_c2c)) .OR. &
- (.not. c_associated(plan_x2kx_c2c))) then
+ if ((.not. c_associated(plan_kx2x_c2c)) .OR. &
+ (.not. c_associated(plan_x2kx_c2c)) .OR. &
+ (.not. c_associated(plan_ky2y_c2r)) .OR. &
+ (.not. c_associated(plan_ky2y_c2r)) .OR. &
+ (.not. c_associated(plan_many_x2kx_c2c)) .OR. &
+ (.not. c_associated(plan_many_kx2x_c2c)) .OR. &
+ (.not. c_associated(plan_many_ky2y_c2r)) .OR. &
+ (.not. c_associated(plan_many_y2ky_r2c)) ) then
ERROR STOP '>> ERROR << 1D plan creation error!!'
end if
- ! Free mem (optional)
- DEALLOCATE(f_xky_l,f_kxky_l)
- DEALLOCATE(bracket_sum_kyx_g,bracket_sum_yx_g)
+
END SUBROUTINE fft1D_plans
!******************************************************************************!
@@ -283,6 +323,7 @@ END SUBROUTINE fft1D_plans
! local variables
INTEGER :: ikx,iky
COMPLEX(xp), DIMENSION(total_nkx,local_nky) :: ikxF, ikyG, ikyF, ikxG
+ COMPLEX(xp), DIMENSION(total_nkx,local_nky) :: invfactor !TEST
REAL(xp):: kxs, ky
! Build the fields to convolve
@@ -303,6 +344,13 @@ END SUBROUTINE fft1D_plans
CALL apply_ExB_NL_factor(ikyG,ExB_NL_factor)
CALL apply_ExB_NL_factor(ikyF,ExB_NL_factor)
CALL apply_ExB_NL_factor(ikxG,ExB_NL_factor)
+ ! ! TEST apply inverse to test the identity
+ ! invfactor = 1._xp/ExB_NL_factor
+ ! CALL apply_ExB_NL_factor(ikxF,invfactor)
+ ! CALL apply_ExB_NL_factor(ikyG,invfactor)
+ ! CALL apply_ExB_NL_factor(ikyF,invfactor)
+ ! CALL apply_ExB_NL_factor(ikxG,invfactor)
+ ! ! END TEST
ENDIF
! Anti Aliasing
DO iky = 1,local_nky
@@ -334,6 +382,7 @@ END SUBROUTINE fft1D_plans
END SUBROUTINE poisson_bracket_and_sum
!******************************************************************************!
+ !******************************************************************************!
!******************************************************************************!
! Apply the exp(-i*x*ky*S*J*dt) factor to the Poisson bracket fields (ikxF, ikyG, ...)
@@ -369,86 +418,189 @@ END SUBROUTINE fft1D_plans
ENDDO
END SUBROUTINE apply_ExB_NL_factor
- ! Apply the exp(i*x*ky*S*J*dt) factor to the real Poisson Bracket sum [f,g]
+ ! High level FFT routines required
+ ! 1D inverse Fourier transform from (kx,ky) to (x,ky), complex to complex
+ SUBROUTINE iFFT_kxky_to_xky(in_kxky,out_xky)
+ IMPLICIT NONE
+ ! arguments
+ COMPLEX(xp), DIMENSION(NX_,local_nky_), INTENT(IN) :: in_kxky
+ COMPLEX(xp), DIMENSION(NX_,local_nky_), INTENT(OUT) :: out_xky
+ ! local variables
+ COMPLEX(xp), DIMENSION(NX_) :: in_kx, out_x
+ INTEGER :: iky
+#ifdef SINGLE_PRECISION
+ !plan many 1D transforms
+ CALL sfftw_execute_dft(plan_many_kx2x_c2c, in_kxky, out_xky)
+#else
+ IF (PLAN_MANY) THEN
+ ! IF (.TRUE.) THEN
+ !plan many 1D transforms
+ ! CALL dfftw_execute_dft(plan_many_kx2x_c2c, in_kxky, out_xky)
+ f_kxky_l = in_kxky
+ CALL dfftw_execute_dft(plan_many_kx2x_c2c, f_kxky_l, f_xky_l)
+ out_xky = f_xky_l
+ ELSE
+ !or loop over 1D transforms (less efficient)
+ DO iky = 1,local_nky_
+ ! in_kx = in_kxky(:,iky)
+ ! CALL dfftw_execute_dft(plan_kx2x_c2c,in_kx,out_x)
+ ! out_xky(:,iky) = out_x
+ f_kx_l = in_kxky(:,iky)
+ CALL dfftw_execute_dft(plan_kx2x_c2c,f_kx_l,f_x_l)
+ out_xky(:,iky) = f_x_l
+ ENDDO
+ ENDIF
+#endif
+ END SUBROUTINE iFFT_kxky_to_xky
+ ! 1D Fourier transform from kx,ky) to (kx,ky), complex to complex
+ SUBROUTINE FFT_xky_to_kxky(in_xky,out_kxky)
+ IMPLICIT NONE
+ ! arguments
+ COMPLEX(xp), DIMENSION(NX_,local_nky_), INTENT(IN) :: in_xky
+ COMPLEX(xp), DIMENSION(NX_,local_nky_), INTENT(OUT) :: out_kxky
+ ! local variables
+ COMPLEX(xp), DIMENSION(NX_) :: in_x, out_kx
+ INTEGER :: iky
+#ifdef SINGLE_PRECISION
+ CALL sfftw_execute_dft(plan_many_x2kx_c2c, in_xky, out_kxky)
+#else
+ IF (PLAN_MANY) THEN
+ ! IF (.TRUE.) THEN
+ !plan many 1D transforms
+ ! CALL dfftw_execute_dft(plan_many_x2kx_c2c, in_xky, out_kxky)
+ f_xky_l = in_xky
+ CALL dfftw_execute_dft(plan_many_x2kx_c2c, f_xky_l, f_kxky_l)
+ out_kxky = f_kxky_l
+ ELSE
+ !or loop over 1D transforms (less efficient)
+ DO iky = 1,local_nky_
+ ! in_x = in_xky(:,iky)
+ ! CALL dfftw_execute_dft(plan_x2kx_c2c,in_x,out_kx)
+ ! out_kxky(:,iky) = out_kx
+ f_x_l = in_xky(:,iky)
+ CALL dfftw_execute_dft(plan_x2kx_c2c,f_x_l,f_kx_l)
+ out_kxky(:,iky) = f_kx_l
+ ENDDO
+ ENDIF
+#endif
+ END SUBROUTINE FFT_xky_to_kxky
+
+ !*****************************************************************************!
+ ! Apply the inverse factor, exp(i*x*ky*S*J*dt), to the real Poisson Bracket sum [f,g]
SUBROUTINE apply_inv_ExB_NL_factor(fyx,inv_ExB_NL_factor)
IMPLICIT NONE
- real(c_xp_r), pointer, INTENT(INOUT) :: fyx(:,:)
- COMPLEX(xp), DIMENSION(NY_halved+1,local_nx_), INTENT(IN) :: inv_ExB_NL_factor
+ ! REAL(xp), DIMENSION(2*NY_halved,local_nx_), INTENT(INOUT) :: fyx
+ real(c_xp_r), pointer, INTENT(INOUT) :: fyx(:,:)
+ COMPLEX(xp), DIMENSION(NY_/2+1,local_nx_), INTENT(IN) :: inv_ExB_NL_factor
! local variables
- REAL(xp), DIMENSION(2*NY_halved,local_nx_) :: tmp_yx_1, tmp_yx_2
- COMPLEX(xp), DIMENSION(NY_halved+1,local_nx_) :: tmp_kyx
+ REAL(xp), DIMENSION(NY_, local_nx_) :: tmp_yx
+ COMPLEX(xp), DIMENSION(NY_/2+1,local_nx_) :: tmp_kyx
integer(C_INTPTR_T) :: ix, iy, iky
! Fill buffer
DO ix = 1,local_nx_
- DO iy = 1,2*NY_halved
- tmp_yx_1(iy,ix) = fyx(iy,ix)
+ DO iy = 1,NY_
+ tmp_yx(iy,ix) = fyx(iy,ix)
+ ! tmp_yx(iy,ix) = bracket_sum_r(iy,ix)
ENDDO
ENDDO
! Fourier real to complex on the second buffer (first buffer is now unusable)
- CALL FFT_yx_to_kyx(tmp_yx_1,tmp_kyx)
+ CALL FFT_yx_to_kyx(tmp_yx,tmp_kyx)
! Treat the result with the ExB NL factor
- DO iky = 1,NY_halved+1
+ DO iky = 1,NY_/2+1
DO ix = 1,local_nx_
- tmp_kyx(iky,ix) = tmp_kyx(iky,ix)*inv_ExB_NL_factor(iky,ix)
+ tmp_kyx(iky,ix) = tmp_kyx(iky,ix)!*inv_ExB_NL_factor(iky,ix)
ENDDO
ENDDO
! Back to Fourier space in the third buffer
- CALL iFFT_kyx_to_yx(tmp_kyx,tmp_yx_2)
+ CALL iFFT_kyx_to_yx(tmp_kyx,tmp_yx)
! Fill the result array with normalization of iFFT
DO ix = 1,local_nx_
- DO iy = 1,2*NY_halved
- fyx(iy,ix) = tmp_yx_2(iy,ix)*inv_Ny_
+ DO iy = 1,NY_
+ fyx(iy,ix) = tmp_yx(iy,ix)*inv_Ny_
+ ! bracket_sum_r(iy,ix) = tmp_yx(iy,ix)*inv_Ny_
ENDDO
ENDDO
END SUBROUTINE apply_inv_ExB_NL_factor
-
- !******************************************************************************!
- ! High level FFT routines
- ! 1D inverse Fourier transform from (kx,ky) to (x,ky), complex to complex
- SUBROUTINE iFFT_kxky_to_xky(in_kxky,out_xky)
- IMPLICIT NONE
- COMPLEX(xp), DIMENSION(NX_,local_nky_), INTENT(IN) :: in_kxky
- COMPLEX(xp), DIMENSION(NX_,local_nky_), INTENT(OUT) :: out_xky
-#ifdef SINGLE_PRECISION
- CALL sfftw_execute_dft(plan_kx2x_c2c, in_kxky, out_xky)
-#else
- CALL dfftw_execute_dft(plan_kx2x_c2c, in_kxky, out_xky)
-#endif
- END SUBROUTINE iFFT_kxky_to_xky
- ! 1D Fourier transform from kx,ky) to (kx,ky), complex to complex
- SUBROUTINE FFT_xky_to_kxky(in_xky,out_kxky)
- IMPLICIT NONE
- COMPLEX(xp), DIMENSION(NX_,local_nky_), INTENT(IN) :: in_xky
- COMPLEX(xp), DIMENSION(NX_,local_nky_), INTENT(OUT) :: out_kxky
-#ifdef SINGLE_PRECISION
- CALL sfftw_execute_dft(plan_x2kx_c2c, in_xky, out_kxky)
-#else
- CALL dfftw_execute_dft(plan_x2kx_c2c, in_xky, out_kxky)
-#endif
- END SUBROUTINE FFT_xky_to_kxky
! 1D Fourier transform from (y,x) to (ky,x), real to complex
SUBROUTINE FFT_yx_to_kyx(in_yx,out_kyx)
IMPLICIT NONE
- REAL(xp), DIMENSION(2*NY_halved,local_nx_), INTENT(IN) :: in_yx
- COMPLEX(xp), DIMENSION(NY_halved+1,local_nx_), INTENT(OUT) :: out_kyx
+ ! arguments
+ REAL(xp), DIMENSION(NY_, local_nx_), INTENT(IN) :: in_yx
+ COMPLEX(xp), DIMENSION(NY_/2+1,local_nx_), INTENT(OUT) :: out_kyx
+ ! local var.
+ REAL(xp), DIMENSION(NY_) :: in_y
+ COMPLEX(xp), DIMENSION(NY_/2+1) :: out_ky
+ ! COMPLEX(xp), DIMENSION(NY_halved+1,local_nx_) :: out_kyx_test
+ INTEGER :: ix, iy, iky
#ifdef SINGLE_PRECISION
- CALL sfftw_execute_dft_r2c(plan_y2ky_r2c, in_yx, out_kyx)
+ CALL sfftw_execute_dft_r2c(plan_many_y2ky_r2c, in_yx, out_kyx)
#else
- CALL dfftw_execute_dft_r2c(plan_y2ky_r2c, in_yx, out_kyx)
+ IF (PLAN_MANY) THEN
+ ! IF (.TRUE.) THEN
+ !plan many 1D transforms
+ ! CALL dfftw_execute_dft_r2c(plan_many_y2ky_r2c, in_yx, out_kyx)
+ f_yx_g = in_yx
+ CALL dfftw_execute_dft_r2c(plan_many_y2ky_r2c, f_yx_g, f_kyx_g)
+ out_kyx = f_kyx_g
+ ELSE
+ ! CALL dfftw_execute_dft_r2c(plan_many_y2ky_r2c, in_yx, out_kyx_test)
+ !or loop over 1D transforms (less efficient but checkable)
+ DO ix = 1,local_nx_
+ DO iy = 1,NY_
+ ! in_y(iy) = in_yx(iy,ix)
+ f_y_g(iy) = in_yx(iy,ix)
+ ENDDO
+ ! CALL dfftw_execute_dft_r2c(plan_y2ky_r2c,in_y,out_ky)
+ CALL dfftw_execute_dft_r2c(plan_y2ky_r2c,f_y_g,f_ky_g)
+ DO iky = 1,NY_/2+1
+ ! out_kyx(iky,ix) = out_ky(iky)
+ out_kyx(iky,ix) = f_ky_g(iky)
+ ENDDO
+ ENDDO
+ ENDIF
#endif
END SUBROUTINE FFT_yx_to_kyx
! 1D Fourier transform from (ky,x) to (y,x), complex to real
SUBROUTINE iFFT_kyx_to_yx(in_kyx,out_yx)
IMPLICIT NONE
- COMPLEX(xp), DIMENSION(NY_halved+1,local_nx_), INTENT(IN) :: in_kyx
- REAL(xp), DIMENSION(2*NY_halved,local_nx_), INTENT(OUT) :: out_yx
+ ! arguments
+ COMPLEX(xp), DIMENSION(NY_/2+1,local_nx_), INTENT(IN) :: in_kyx
+ REAL(xp), DIMENSION(NY_ ,local_nx_), INTENT(OUT) :: out_yx
+ ! local var.
+ COMPLEX(xp), DIMENSION(Ny_/2+1) :: in_ky
+ REAL(xp), DIMENSION(NY_) :: out_y
+ REAL(xp), DIMENSION(NY_,local_nx_) :: out_yx_test
+ INTEGER :: ix, iy, iky
#ifdef SINGLE_PRECISION
- CALL sfftw_execute_dft_c2r(plan_ky2y_c2r, in_kyx, out_yx)
+ CALL sfftw_execute_dft_c2r(plan_many_ky2y_c2r, in_kyx, out_yx)
#else
- CALL dfftw_execute_dft_c2r(plan_ky2y_c2r, in_kyx, out_yx)
+ IF (PLAN_MANY) THEN
+ ! IF (.TRUE.) THEN
+ !plan many 1D transforms
+ ! CALL dfftw_execute_dft_c2r(plan_many_ky2y_c2r, in_kyx, out_yx)
+ f_kyx_g = in_kyx
+ CALL dfftw_execute_dft_c2r(plan_many_ky2y_c2r, f_kyx_g, f_yx_g)
+ out_yx = f_yx_g
+ ELSE
+ ! CALL dfftw_execute_dft_c2r(plan_many_ky2y_c2r, in_kyx, out_yx_test)
+ !or loop over 1D transforms (less efficient)
+ DO ix = 1,local_nx_
+ DO iky = 1,NY_/2+1
+ ! in_ky(iky) = in_kyx(iky,ix)
+ f_ky_g(iky) = in_kyx(iky,ix)
+ ENDDO
+ ! CALL dfftw_execute_dft_c2r(plan_ky2y_c2r,in_ky,out_y)
+ CALL dfftw_execute_dft_c2r(plan_ky2y_c2r,f_ky_g,f_y_g)
+ DO iy = 1,NY_
+ ! out_yx(iy,ix) = out_y(iy)
+ out_yx(iy,ix) = f_y_g(iy)
+ ENDDO
+ ENDDO
+ ENDIF
#endif
END SUBROUTINE iFFT_kyx_to_yx
+ !******************************************************************************!
!******************************************************************************!
SUBROUTINE finalize_plans
@@ -464,6 +616,8 @@ END SUBROUTINE fft1D_plans
call fftw_free(cdatac_f)
call fftw_free(cdatac_g)
call fftw_free(cdatac_c)
+ DEALLOCATE(f_xky_l,f_kxky_l)
+ DEALLOCATE(f_kyx_g,f_yx_g)
END SUBROUTINE finalize_plans
END MODULE fourier
diff --git a/src/nonlinear_mod.F90 b/src/nonlinear_mod.F90
index 38dd61bd32eead22cb775ac6e36ae021aacd7732..a6e3c6faf0b5e7999475364976eb7ccb235b2163 100644
--- a/src/nonlinear_mod.F90
+++ b/src/nonlinear_mod.F90
@@ -7,7 +7,7 @@ MODULE nonlinear
local_np,ngp,parray,pmax,&
local_nj,ngj,jarray,jmax, local_nj_offset, dmax,&
kyarray, AA_y, local_nky, inv_Ny,&
- total_nkx,kxarray, AA_x, inv_Nx,&
+ total_nkx,kxarray, AA_x, inv_Nx,local_nx, Ny, &
local_nz,ngz,zarray,nzgrid, deltakx, iky0, contains_kx0, contains_ky0
USE model, ONLY : LINEARITY, EM, ikxZF, ZFamp, ExB_NL_CORRECTION
USE closure, ONLY : evolve_mom, nmaxarray
@@ -56,6 +56,7 @@ SUBROUTINE compute_nonlinear
IMPLICIT NONE
INTEGER :: iz,ij,ip,eo,ia,ikx,iky,izi,ipi,iji,ini,isi
INTEGER :: ikxExBp, ikxExBn ! Negative and positive ExB flow indices
+ COMPLEX(xp), DIMENSION(Ny/2+1,local_nx) :: invinvfactor ! TEST
z:DO iz = 1,local_nz
izi = iz + ngz/2
j:DO ij = 1,local_nj ! Loop over Laguerre moments
@@ -122,6 +123,8 @@ SUBROUTINE compute_nonlinear
! Apply the ExB shearing rate factor before going back to k-space
IF (ExB_NL_CORRECTION) THEN
CALL apply_inv_ExB_NL_factor(bracket_sum_r,inv_ExB_NL_factor)
+ ! invinvfactor = 1._xp/inv_ExB_NL_factor
+ ! CALL apply_inv_ExB_NL_factor(bracket_sum_r,invinvfactor)
ENDIF
! Put the real nonlinear product back into k-space
#ifdef SINGLE_PRECISION
diff --git a/src/stepon.F90 b/src/stepon.F90
index b9aee7e906a500cade18ba8bb6c0cd5faf3a95e0..77c761285b4af29c2ec280919bb3cdd7a3d101e6 100644
--- a/src/stepon.F90
+++ b/src/stepon.F90
@@ -20,13 +20,13 @@ SUBROUTINE stepon
SUBSTEPS:DO num_step=1,ntimelevel ! eg RK4 compute successively k1, k2, k3, k4
!----- TEST !-----
- ! Update the ExB shear flow for the next step
- ! This call includes :
- ! - the ExB shear value (s(ky)) update for the next time step
- ! - the kx grid update
- ! - the ExB NL correction factor update (exp(+/- ixkySdts))
- ! - (optional) the kernel, poisson op. and ampere op update
- CALL Update_ExB_shear_flow(num_step)
+ ! ! Update the ExB shear flow for the next step
+ ! ! This call includes :
+ ! ! - the ExB shear value (s(ky)) update for the next time step
+ ! ! - the kx grid update
+ ! ! - the ExB NL correction factor update (exp(+/- ixkySdts))
+ ! ! - (optional) the kernel, poisson op. and ampere op update
+ ! CALL Update_ExB_shear_flow(num_step)
!-----END TEST !-----
!----- BEFORE: All fields+ghosts are updated for step = n
! Compute right hand side from current fields