diff --git a/src/grid_mod.F90 b/src/grid_mod.F90 index 83f5fec52390ae36df6039a7d0c89dab7995bdcb..4f39aa206fcbbdc652aaa76ea86a5fd3c7ac4d54 100644 --- a/src/grid_mod.F90 +++ b/src/grid_mod.F90 @@ -49,8 +49,10 @@ MODULE grid REAL(dp), DIMENSION(:), ALLOCATABLE, PUBLIC :: kparray ! kperp array REAL(dp), PUBLIC, PROTECTED :: deltakr, deltakz INTEGER, PUBLIC, PROTECTED :: ikrs, ikre, ikzs, ikze, ikps, ikpe - INTEGER, PUBLIC, PROTECTED :: ikr_0, ikz_0 ! Indices of k-grid origin + INTEGER, PUBLIC, PROTECTED :: ikr_0, ikz_0, ikr_max, ikz_max ! Indices of k-grid origin and max INTEGER, PUBLIC :: ikr, ikz, ip, ij, ikp ! counters + LOGICAL, PUBLIC, PROTECTED :: contains_kr0 = .false. ! rank of the proc containing kr=0 indices + LOGICAL, PUBLIC, PROTECTED :: contains_krmax = .false. ! rank of the proc containing kr=max indices ! Grid containing the polynomials degrees INTEGER, DIMENSION(:), ALLOCATABLE, PUBLIC :: parray_e @@ -164,10 +166,19 @@ CONTAINS deltakr = 2._dp*PI/Lr ENDIF - ! Discretized kr positions ordered as dk*(0 1 2 3) + ! Creating a grid ordered as dk*(0 1 2 3) DO ikr = ikrs,ikre krarray(ikr) = REAL(ikr-1,dp) * deltakr - IF (krarray(ikr) .EQ. 0) ikr_0 = ikr + ! Finding kr=0 + IF (krarray(ikr) .EQ. 0) THEN + ikr_0 = ikr + contains_kr0 = .true. + ENDIF + ! Finding krmax + IF (krarray(ikr) .EQ. (Nr/2+1)*deltakr) THEN + ikr_max = ikr + contains_krmax = .true. + ENDIF END DO ! Orszag 2/3 filter @@ -193,18 +204,23 @@ CONTAINS ikze = Nkz ALLOCATE(kzarray(ikzs:ikze)) - ! Grid spacings and discretized kz positions ordered as dk*(0 1 2 3 -2 -1) IF (Lz .EQ. 0) THEN ! 1D linear case deltakz = 1._dp kzarray(1) = 0 ikz_0 = 1 + ikz_max = 1 ELSE deltakz = 2._dp*PI/Lz + ! Creating a grid ordered as dk*(0 1 2 3 -2 -1) DO ikz = ikzs,ikze kzarray(ikz) = deltakz*(MODULO(ikz-1,Nkz/2)-Nkz/2*FLOOR(2.*real(ikz-1)/real(Nkz))) if (ikz .EQ. Nz/2+1) kzarray(ikz) = -kzarray(ikz) + ! Finding kz=0 IF (kzarray(ikz) .EQ. 0) ikz_0 = ikz + ! Finding kzmax + IF (kzarray(ikr) .EQ. (Nz/2)*deltakr) ikr_max = ikr END DO + IF(my_id .EQ. 0) write(*,*) kzarray ENDIF ! Orszag 2/3 filter diff --git a/src/inital.F90 b/src/inital.F90 index 164a084be2aaaa433b7e851d41bd4cc776cbab12..a4502cb54a3b0665faab9737f826a34947be1404 100644 --- a/src/inital.F90 +++ b/src/inital.F90 @@ -100,9 +100,9 @@ SUBROUTINE init_moments END DO END DO - IF ( ikrs .EQ. 1 ) THEN + IF ( contains_kr0 ) THEN DO ikz=2,Nkz/2 !symmetry at kr = 0 - moments_e( ip,ij,1,ikz, :) = moments_e( ip,ij,1,Nkz+2-ikz, :) + moments_e( ip,ij,ikr_0,ikz, :) = moments_e( ip,ij,ikr_0,Nkz+2-ikz, :) END DO ENDIF @@ -119,9 +119,9 @@ SUBROUTINE init_moments END DO END DO - IF ( ikrs .EQ. 1 ) THEN + IF ( contains_kr0 ) THEN DO ikz=2,Nkz/2 !symmetry at kr = 0 - moments_i( ip,ij,1,ikz, :) = moments_i( ip,ij,1,Nkz+2-ikz, :) + moments_i( ip,ij,ikr_0,ikz, :) = moments_i( ip,ij,ikr_0,Nkz+2-ikz, :) END DO ENDIF @@ -184,10 +184,11 @@ SUBROUTINE init_phi END DO !symmetry at kr = 0 to keep real inverse transform - IF ( ikrs .EQ. 1 ) THEN + IF ( contains_kr0 ) THEN DO ikz=2,Nkz/2 - phi(1,ikz) = phi(1,Nkz+2-ikz) + phi(ikr_0,ikz) = phi(ikr_0,Nkz+2-ikz) END DO + phi(ikr_0,Nz/2) = REAL(phi(ikr_0,Nz/2)) !origin must be real ENDIF !**** Cancel previous moments initialization diff --git a/src/stepon.F90 b/src/stepon.F90 index b784a429fce26fd1e569cf291ea7ebbc0e0d6415..190faf56326550aacb86541b14e3e1be4beb61a9 100644 --- a/src/stepon.F90 +++ b/src/stepon.F90 @@ -129,34 +129,34 @@ SUBROUTINE stepon END SUBROUTINE anti_aliasing SUBROUTINE enforce_symetry ! Force X(k) = X(N-k)* complex conjugate symmetry - ! Electron moments - DO ip=ips_e,ipe_e - DO ij=ijs_e,ije_e - DO ikz=2,Nkz/2 !symmetry at kr = 0 - moments_e( ip,ij,1,ikz, :) = CONJG(moments_e( ip,ij,1,Nkz+2-ikz, :)) + IF ( contains_kr0 ) THEN + ! Electron moments + DO ip=ips_e,ipe_e + DO ij=ijs_e,ije_e + DO ikz=2,Nkz/2 !symmetry at kr = 0 + moments_e( ip,ij,ikr_0,ikz, :) = CONJG(moments_e( ip,ij,ikr_0,Nkz+2-ikz, :)) + END DO + ! must be real at origin + moments_e(ip,ij, ikr_0,ikz_0, :) = REAL(moments_e(ip,ij, ikr_0,ikz_0, :)) END DO - ! must be real at origin and top right - moments_e(ip,ij, ikr_0,ikz_0, :) = REAL(moments_e(ip,ij, ikr_0,ikz_0, :)) - moments_e(ip,ij, Nr/2,Nz/2 , :) = REAL(moments_e(ip,ij, Nr/2,Nz/2 , :)) END DO - END DO - ! Ion moments - DO ip=ips_i,ipe_i - DO ij=ijs_i,ije_i - DO ikz=2,Nkz/2 !symmetry at kr = 0 - moments_i( ip,ij,1,ikz, :) = CONJG(moments_i( ip,ij,1,Nkz+2-ikz, :)) + ! Ion moments + DO ip=ips_i,ipe_i + DO ij=ijs_i,ije_i + DO ikz=2,Nkz/2 !symmetry at kr = 0 + moments_i( ip,ij,ikr_0,ikz, :) = CONJG(moments_i( ip,ij,ikr_0,Nkz+2-ikz, :)) + END DO + ! must be real at origin and top right + moments_i(ip,ij, ikr_0,ikz_0, :) = REAL(moments_i(ip,ij, ikr_0,ikz_0, :)) END DO - ! must be real at origin and top right - moments_i(ip,ij, ikr_0,ikz_0, :) = REAL(moments_i(ip,ij, ikr_0,ikz_0, :)) - moments_i(ip,ij, Nr/2,Nz/2 , :) = REAL(moments_i(ip,ij, Nr/2,Nz/2 , :)) END DO - END DO - ! Phi - DO ikz=2,Nkz/2 !symmetry at kr = 0 - phi(1,ikz) = phi(1,Nkz+2-ikz) - END DO - ! must be real at origin and top right - phi(ikr_0,ikz_0) = REAL(phi(ikr_0,ikz_0)) + ! Phi + DO ikz=2,Nkz/2 !symmetry at kr = 0 + phi(ikr_0,ikz) = phi(ikr_0,Nkz+2-ikz) + END DO + ! must be real at origin + phi(ikr_0,ikz_0) = REAL(phi(ikr_0,ikz_0)) + ENDIF END SUBROUTINE enforce_symetry END SUBROUTINE stepon