diff --git a/src/ExB_shear_flow_mod.F90 b/src/ExB_shear_flow_mod.F90
index 31d3051edbd7d3654eee66f25ba7a645b61b95fb..73377126fe8771d2e41f9c7fa0e3804036e8302c 100644
--- a/src/ExB_shear_flow_mod.F90
+++ b/src/ExB_shear_flow_mod.F90
@@ -9,6 +9,7 @@ MODULE ExB_shear_flow
REAL(xp), PUBLIC, PROTECTED :: gamma_E = 0._xp ! ExB background shearing rate \gamma_E
REAL(xp), PUBLIC, PROTECTED :: t0, inv_t0 = 0._xp ! charact. shear time
REAL(xp), DIMENSION(:), ALLOCATABLE, PUBLIC, PROTECTED :: sky_ExB ! shift of the kx modes, kx* = kx + s(ky)
+ REAL(xp), DIMENSION(:), ALLOCATABLE, PUBLIC, PROTECTED :: sky_ExB_full ! full ky version
REAL(xp), DIMENSION(:), ALLOCATABLE, PUBLIC, PROTECTED :: dkx_ExB ! correction to obtain the exact kx mode
INTEGER, DIMENSION(:), ALLOCATABLE, PUBLIC, PROTECTED :: jump_ExB ! jump to do to shift the kx grids
LOGICAL, DIMENSION(:), ALLOCATABLE, PUBLIC, PROTECTED :: shiftnow_ExB ! Indicates if there is a line to shift
@@ -16,19 +17,20 @@ MODULE ExB_shear_flow
COMPLEX(xp),DIMENSION(:,:), ALLOCATABLE, PUBLIC, PROTECTED :: inv_ExB_NL_factor
LOGICAL, PUBLIC, PROTECTED :: ExB = .false. ! presence of ExB background shearing rate
! Routines
- PUBLIC :: Setup_ExB_shear_flow, Apply_ExB_shear_flow, Update_ExB_shear_flow
+ PUBLIC :: Setup_ExB_shear_flow, Array_shift_ExB_shear_flow, Update_ExB_shear_flow
CONTAINS
! Setup the variables for the ExB shear
SUBROUTINE Setup_ExB_shear_flow(ExBrate)
- USE grid, ONLY : Nx, local_nky, local_nx, Ny, deltakx, deltaky
+ USE grid, ONLY: Nx, local_nky, total_nky, local_nx, Ny, deltakx, deltaky
+ USE geometry, ONLY: Cyq0_x0
IMPLICIT NONE
REAL(xp), INTENT(IN) :: ExBrate
! Setup the ExB shearing rate and aux var
- gamma_E = ExBrate
- IF(gamma_E .GT. 0._xp) THEN
+ gamma_E = -ExBrate*Cyq0_x0
+ IF(abs(gamma_E) .GT. EPSILON(gamma_E)) THEN
ExB = .TRUE.
t0 = deltakx/deltaky/gamma_E
inv_t0 = 1._xp/t0
@@ -41,6 +43,8 @@ CONTAINS
! Setup the ExB shift array
ALLOCATE(sky_ExB(local_nky))
sky_ExB = 0._xp
+ ALLOCATE(sky_ExB_full(total_nky))
+ sky_ExB_full = 0._xp
! Setup the kx correction array
ALLOCATE(dkx_ExB(local_nky))
@@ -63,16 +67,14 @@ CONTAINS
END SUBROUTINE Setup_ExB_shear_flow
! Update according to the current ExB shear value
- ! -the grids
+ ! -shift the grids
! -the spatial operators
! -the fields by imposing a shift on kx
- SUBROUTINE Apply_ExB_shear_flow
- USE basic, ONLY: time
- USE grid, ONLY: local_nky, kxarray, update_grids, &
- total_nkx, deltakx, kx_min, kx_max
+ SUBROUTINE Array_shift_ExB_shear_flow
+ USE grid, ONLY: local_nky, update_grids, &
+ total_nkx, deltakx, kx_min, kx_max, kxarray0
USE prec_const, ONLY: PI
USE fields, ONLY: moments, phi, psi
- USE geometry, ONLY: gxx,gxy,inv_hatB2, evaluate_magn_curv
USE numerics, ONLY: evaluate_EM_op, evaluate_kernels
IMPLICIT NONE
! local var
@@ -105,10 +107,17 @@ CONTAINS
ikx = i_ - total_nkx/2 + 1
ENDIF
ikx_s = ikx + jump_ExB(iky)
+ ! adjust the shift according
+ IF (ikx_s .LE. 0) &
+ ikx_s = ikx_s + total_nkx
+ IF (ikx_s .GT. total_nkx) &
+ ikx_s = ikx_s - total_nkx
! We test if the shifted modes are still in contained in our resolution
- IF ( ((ikx_s .GT. 0 ) .AND. (ikx_s .LE. total_nkx )) .AND. &
- (((ikx .LE. (total_nkx/2+1)) .AND. (ikx_s .LE. (total_nkx/2+1))) .OR. &
- ((ikx .GT. (total_nkx/2+1)) .AND. (ikx_s .GT. (total_nkx/2+1)))) ) THEN
+ ! IF ( ((ikx_s .GT. 0 ) .AND. (ikx_s .LE. total_nkx )) .AND. &
+ ! (((ikx .LE. (total_nkx/2+1)) .AND. (ikx_s .LE. (total_nkx/2+1))) .OR. &
+ ! ((ikx .GT. (total_nkx/2+1)) .AND. (ikx_s .GT. (total_nkx/2+1)))) ) THEN
+ IF ( (kxarray0(ikx)+jump_ExB(iky)*deltakx .LE. kx_max) .AND. &
+ (kxarray0(ikx)+jump_ExB(iky)*deltakx .GE. kx_min)) THEN
moments(:,:,:,iky,ikx,:,:) = moments(:,:,:,iky,ikx_s,:,:)
phi(iky,ikx,:) = phi(iky,ikx_s,:)
psi(iky,ikx,:) = psi(iky,ikx_s,:)
@@ -123,20 +132,20 @@ CONTAINS
ENDIF
ENDDO
ENDIF
- END SUBROUTINE Apply_ExB_shear_flow
+ END SUBROUTINE Array_shift_ExB_shear_flow
! update the ExB shear value for the next time step
SUBROUTINE Update_ExB_shear_flow
- USE basic, ONLY: dt, time, chrono_ExBs, start_chrono, stop_chrono
- USE grid, ONLY: local_nky, local_nky_offset, kyarray, kyarray_full, inv_dkx, xarray, Nx, Ny, &
- local_nx, local_nx_offset, deltax, kxarray, &
- ikyarray, inv_ikyarray, deltakx, deltaky, update_grids
+ USE basic, ONLY: dt, time
+ USE grid, ONLY: local_nky, local_nky_offset, kyarray, inv_dkx, xarray, Nx, Ny, &
+ local_nx, local_nx_offset, kyarray_full,&
+ ikyarray, inv_ikyarray, deltaky, update_grids
USE geometry, ONLY: gxx,gxy,gyy,inv_hatB2, evaluate_magn_curv
USE numerics, ONLY: evaluate_EM_op, evaluate_kernels
IMPLICIT NONE
! local var
INTEGER :: iky, ix
- REAL(xp):: dt_ExB, ky, kx, J_dp, inv_J, x, dtshift
+ REAL(xp):: dt_ExB, J_dp, inv_J, x
! do nothing if no ExB
IF(ExB) THEN
! reset the ExB shift, jumps and flags
@@ -149,7 +158,8 @@ CONTAINS
! zero-shiftings that may be majoritary.
shiftnow_ExB(iky) = (abs(jump_ExB(iky)) .GT. 0)
ENDDO
-
+ ! Update the full skyExB array too
+ sky_ExB_full = sky_ExB_full - kyarray_full*gamma_E*dt
! Update the ExB nonlinear factor...
DO iky = 1,local_Nky ! WARNING: Local indices ky loop
! for readability
@@ -161,7 +171,8 @@ CONTAINS
DO ix = 1,Nx
x = xarray(ix)
! assemble the ExB nonlin factor
- ExB_NL_factor(ix,iky) = EXP(imagu*x*dkx_ExB(iky))
+ !ExB_NL_factor(ix,iky) = EXP(imagu*x*dkx_ExB(iky))
+ ExB_NL_factor(ix,iky) = EXP(imagu*x*sky_ExB(iky)) !GENE does that???
ENDDO
ENDDO
! ... and the inverse
@@ -174,15 +185,14 @@ CONTAINS
DO ix = 1,local_nx
x = xarray(ix+local_nx_offset)
! assemble the inverse ExB nonlin factor
- inv_ExB_NL_factor(iky,ix) = EXP(-imagu*x*gamma_E*J_dp*deltaky*dt_ExB)
+ ! inv_ExB_NL_factor(iky,ix) = EXP(-imagu*x*gamma_E*J_dp*deltaky*dt_ExB)
+ inv_ExB_NL_factor(iky,ix) = EXP(-imagu*x*sky_ExB_full(iky)) !GENE does that???
ENDDO
ENDDO
ENDIF
! We update the operators and grids
! update the grids
CALL update_grids(dkx_ExB,gxx,gxy,gyy,inv_hatB2)
- !write(42,*), kxarray(5,3)
-
! update the EM op., the kernels and the curvature op.
CALL evaluate_kernels
CALL evaluate_EM_op
diff --git a/src/auxval.F90 b/src/auxval.F90
index 0cf3960c669324a667f5e424b8892ea649125a4d..f549564374ec62a35ebd9478b3385f94a423c2ae 100644
--- a/src/auxval.F90
+++ b/src/auxval.F90
@@ -2,9 +2,9 @@ subroutine auxval
! Set auxiliary values, at beginning of simulation
USE basic, ONLY: str, speak
- USE grid, ONLY: local_np, local_np_offset, total_np, local_nj, local_nj_offset, total_nj,&
+ USE grid, ONLY: local_np, local_np_offset, total_np, local_nj, local_nj_offset,&
local_nky, local_nky_offset, total_nky, local_nkx, local_nkx_offset, dmax,&
- local_nz, local_nz_offset, total_nz, init_grids_data, set_grids
+ local_nz, local_nz_offset, init_grids_data, set_grids
!USE array
USE model, ONLY: Na, EM, LINEARITY, N_HD, ExBrate
USE fourier, ONLY: init_grid_distr_and_plans
diff --git a/src/control.F90 b/src/control.F90
index bb743ee641a871f4361285e30945f621603a3b33..a537953125db97c76e8b5b8709e097402cc5f4ef 100644
--- a/src/control.F90
+++ b/src/control.F90
@@ -1,14 +1,16 @@
SUBROUTINE control
! Control the run
- use basic, ONLY: str,daytim,speak,basic_data,&
+ use basic, ONLY: str,daytim,speak,basic_data,&
nlend,step,increase_step,increase_time,increase_cstep,&
- chrono_runt,chrono_step, chrono_diag, start_chrono, stop_chrono
- use prec_const, ONLY: xp, stdout
- USE parallel, ONLY: ppinit
- USE initial, ONLY: initialize
- USE mpi, ONLY: MPI_COMM_WORLD
- USE diagnostics, ONLY: diagnose
+ chrono_runt,chrono_step, chrono_diag, chrono_ExBs,&
+ start_chrono, stop_chrono
+ use prec_const, ONLY: xp, stdout
+ USE parallel, ONLY: ppinit
+ USE initial, ONLY: initialize
+ USE mpi, ONLY: MPI_COMM_WORLD
+ USE diagnostics, ONLY: diagnose
+ USE ExB_shear_flow, ONLY: Array_shift_ExB_shear_flow, Update_ExB_shear_flow
IMPLICIT NONE
REAL(xp) :: t_init_diag_0, t_init_diag_1
INTEGER :: ierr
@@ -70,6 +72,18 @@ SUBROUTINE control
CALL tesend
IF( nlend ) EXIT ! exit do loop
+ ! 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 start_chrono(chrono_ExBs)
+ CALL Update_ExB_shear_flow
+ ! Shift the arrays if the shear value sky is too high
+ CALL Array_shift_ExB_shear_flow
+ CALL stop_chrono(chrono_ExBs)
+
! Increment steps and csteps (private in basic module)
CALL increase_step
CALL increase_cstep
diff --git a/src/diagnostics_mod.F90 b/src/diagnostics_mod.F90
index db03142ced064f5330848d69684b91895aa36053..8d4f9879281a950a165bc851dd950a7c8a51e0a1 100644
--- a/src/diagnostics_mod.F90
+++ b/src/diagnostics_mod.F90
@@ -120,7 +120,7 @@ CONTAINS
SUBROUTINE diagnose(kstep)
! Diagnostics, writing simulation state to disk
- USE basic, ONLY: lu_in, chrono_runt, cstep, dt, time, tmax, display_h_min_s
+ USE basic, ONLY: lu_in, chrono_runt, cstep, time, display_h_min_s
USE processing, ONLY: pflux_x, hflux_x
USE parallel, ONLY: my_id
IMPLICIT NONE
@@ -159,7 +159,7 @@ CONTAINS
USE basic, ONLY: speak,chrono_runt,&
cstep,iframe0d,iframe3d,iframe5d,crashed
USE grid, ONLY: &
- parray_full,pmax,jarray_full,jmax, kparray, &
+ parray_full,jarray_full, kparray, &
kyarray_full,kxarray_full,zarray_full, ngz, total_nz, local_nz, ieven,&
local_Nky, total_nky, local_nkx, total_nkx
USE geometry, ONLY: gxx, gxy, gxz, gyy, gyz, gzz, &
@@ -390,7 +390,6 @@ CONTAINS
USE futils, ONLY: append, attach, getatt, creatd
USE prec_const
USE processing
- USE model, ONLY: Na
IMPLICIT NONE
! Time measurement data
CALL append(fidres, "/profiler/Tc_rhs", REAL(chrono_mrhs%ttot,dp),ionode=0)
@@ -436,7 +435,6 @@ CONTAINS
USE CLA, ONLY: Sf
#endif
IMPLICIT NONE
- CHARACTER(50) :: dset_name
iframe2d=iframe2d+1
CALL append(fidres,"/data/var2d/time", REAL(time,dp), ionode=0)
CALL append(fidres,"/data/var2d/cstep",REAL(cstep,dp),ionode=0)
diff --git a/src/fourier_mod.F90 b/src/fourier_mod.F90
index d6d1af3173ef2d1d283d4c7b80e0eb6f7b396996..1f7951e3b6c398955fd53dec4414d5683d6bf0a6 100644
--- a/src/fourier_mod.F90
+++ b/src/fourier_mod.F90
@@ -44,11 +44,10 @@ MODULE fourier
!******************************************************************************!
!------------- Initialize the grid distribution and plans -------------
! If we use the 2D fftw routines, the fftw library decides the best data distribution
- SUBROUTINE init_grid_distr_and_plans(FFT2D,Nx,Ny,communicator,&
+ SUBROUTINE init_grid_distr_and_plans(Nx,Ny,communicator,&
local_nx_ptr,local_nx_ptr_offset,local_nky_ptr,local_nky_ptr_offset)
USE basic, ONLY: speak
IMPLICIT NONE
- LOGICAL, INTENT(IN) :: FFT2D
INTEGER, INTENT(IN) :: Nx,Ny, communicator
INTEGER(C_INTPTR_T), INTENT(OUT) :: local_nx_ptr,local_nx_ptr_offset
INTEGER(C_INTPTR_T), INTENT(OUT) :: local_nky_ptr,local_nky_ptr_offset
@@ -57,7 +56,7 @@ MODULE fourier
NY_halved = NY_/2 + 1
inv_Ny_ = 1._xp/REAL(2*NY_halved,xp)
! Call FFTW 2D mpi routines to distribute the data and init 2D MPI FFTW plans
- CALL fft2D_distr_and_plans(Nx,Ny,communicator,&
+ 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
@@ -67,10 +66,10 @@ MODULE fourier
END SUBROUTINE init_grid_distr_and_plans
!------------- 2D fft initialization and mpi distribution
- SUBROUTINE fft2D_distr_and_plans(Nx,Ny,communicator,&
+ SUBROUTINE fft2D_distr_and_plans(communicator,&
local_nx_ptr,local_nx_ptr_offset,local_nky_ptr,local_nky_ptr_offset)
IMPLICIT NONE
- INTEGER, INTENT(IN) :: Nx,Ny, communicator
+ INTEGER, INTENT(IN) :: communicator
! Distribution in the real space along x
INTEGER(C_INTPTR_T), INTENT(OUT) :: local_nx_ptr, local_nx_ptr_offset
! Distribution in the fourier space along ky
@@ -138,18 +137,14 @@ MODULE fourier
!------------- 1D initialization with balanced data distribution
SUBROUTINE fft1D_plans
USE utility, ONLY: decomp1D
- USE parallel, ONLY: num_procs_ky, rank_ky
IMPLICIT NONE
! local var
integer(C_INTPTR_T) :: rank ! rank of each 1D fourier transforms
integer(C_INTPTR_T) :: n ! size of the data to fft
integer(C_INTPTR_T) :: howmany ! howmany 1D fourier transforms
- COMPLEX, DIMENSION(:,:), ALLOCATABLE:: in, out
integer(C_INTPTR_T) :: inembed, onembed
integer(C_INTPTR_T) :: istride, ostride
integer(C_INTPTR_T) :: idist, odist
- integer(C_INTPTR_T) :: sign
- integer(C_INTPTR_T) :: flags
!! Plan of the 4 1D many transforms required
!----------- 1: FFTx and inv through local ky data
@@ -273,12 +268,11 @@ END SUBROUTINE fft1D_plans
! module variable (convolution theorem)
SUBROUTINE poisson_bracket_and_sum( ky_array, kx_array, inv_Ny, inv_Nx, AA_y, AA_x,&
local_nky, total_nkx, F_, G_,&
- ExB, ExB_NL_factor, dkx_ExB, sum_real_)
- USE parallel, ONLY: my_id, num_procs_ky, comm_ky, rank_ky
+ ExB, ExB_NL_factor, sum_real_)
IMPLICIT NONE
INTEGER, INTENT(IN) :: local_nky,total_nkx
REAL(xp), INTENT(IN) :: inv_Nx, inv_Ny
- REAL(xp), DIMENSION(local_nky), INTENT(IN) :: ky_array, AA_y, dkx_ExB
+ REAL(xp), DIMENSION(local_nky), INTENT(IN) :: ky_array, AA_y
REAL(xp), DIMENSION(total_nkx), INTENT(IN) :: AA_x
REAL(xp), DIMENSION(local_nky,total_nkx), INTENT(IN) :: kx_array
COMPLEX(c_xp_c), DIMENSION(local_nky,total_nkx), &
@@ -298,10 +292,10 @@ END SUBROUTINE fft1D_plans
DO ikx = 1,total_nkx
ky = ky_array(iky)
kxs = kx_array(iky,ikx)
- ikxF(ikx,iky) = imagu*kxs*F_(iky,ikx)*AA_y(iky)*AA_x(ikx)
- ikyG(ikx,iky) = imagu*ky *G_(iky,ikx)*AA_y(iky)*AA_x(ikx)
- ikyF(ikx,iky) = imagu*ky *F_(iky,ikx)*AA_y(iky)*AA_x(ikx)
- ikxG(ikx,iky) = imagu*kxs*G_(iky,ikx)*AA_y(iky)*AA_x(ikx)
+ ikxF(ikx,iky) = imagu*kxs*F_(iky,ikx)
+ ikyG(ikx,iky) = imagu*ky *G_(iky,ikx)
+ ikyF(ikx,iky) = imagu*ky *F_(iky,ikx)
+ ikxG(ikx,iky) = imagu*kxs*G_(iky,ikx)
ENDDO
ENDDO
IF(ExB) THEN
@@ -311,6 +305,13 @@ END SUBROUTINE fft1D_plans
CALL apply_ExB_NL_factor(ikyF,ExB_NL_factor)
CALL apply_ExB_NL_factor(ikxG,ExB_NL_factor)
ENDIF
+ ! Anti Aliasing
+ DO iky = 1,local_nky
+ ikxF(:,iky) = ikxF(:,iky)*AA_y(iky)*AA_x(:)
+ ikyG(:,iky) = ikyG(:,iky)*AA_y(iky)*AA_x(:)
+ ikyF(:,iky) = ikyF(:,iky)*AA_y(iky)*AA_x(:)
+ ikxG(:,iky) = ikxG(:,iky)*AA_y(iky)*AA_x(:)
+ ENDDO
!-------------------- First term df/dx x dg/dy --------------------
#ifdef SINGLE_PRECISION
call fftwf_mpi_execute_dft_c2r(planb, ikxF, real_data_f)
@@ -342,38 +343,63 @@ END SUBROUTINE fft1D_plans
COMPLEX(xp), DIMENSION(NX_,local_nky_), INTENT(INOUT) :: fkxky
COMPLEX(xp), DIMENSION(NX_,local_nky_), INTENT(IN) :: ExB_NL_factor
! local variables
- COMPLEX(xp), DIMENSION(NX_,local_nky_) :: fxky
- CALL iFFT_kxky_to_xky(fkxky,fxky)
- fxky = fxky*ExB_NL_factor*inv_Nx_
- CALL FFT_xky_to_kxky(fxky,fkxky)
+ COMPLEX(xp), DIMENSION(NX_,local_nky_) :: tmp_kxky, tmp_xky
+ INTEGER :: ix,ikx,iky
+ ! Fill the buffer
+ DO iky = 1,local_nky_
+ DO ikx = 1,NX_
+ tmp_kxky(ikx,iky) = fkxky(ikx,iky)
+ ENDDO
+ ENDDO
+ ! go to x,ky
+ CALL iFFT_kxky_to_xky(tmp_kxky,tmp_xky)
+ ! multiply result by NL factor
+ DO iky = 1,local_nky_
+ DO ix = 1,NX_
+ tmp_xky(ix,iky) = inv_Nx_*tmp_xky(ix,iky)*ExB_NL_factor(ix,iky)
+ ENDDO
+ ENDDO
+ ! back to kx,ky
+ CALL FFT_xky_to_kxky(tmp_xky,tmp_kxky)
+ ! fill output array
+ DO iky = 1,local_nky_
+ DO ikx = 1,NX_
+ fkxky(ikx,iky) = tmp_kxky(ikx,iky)
+ ENDDO
+ 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]
SUBROUTINE apply_inv_ExB_NL_factor(fyx,inv_ExB_NL_factor)
IMPLICIT NONE
- real(c_xp_r), pointer, INTENT(INOUT) :: fyx(:,:)
+ real(c_xp_r), pointer, INTENT(INOUT) :: fyx(:,:)
COMPLEX(xp), DIMENSION(NY_halved,local_nx_), INTENT(IN) :: inv_ExB_NL_factor
! local variables
- REAL(xp), DIMENSION(2*NY_halved,local_nx_) :: buffer_1, buffer_2
- COMPLEX(xp), DIMENSION(NY_halved+1,local_nx_) :: fkyx
- INTEGER :: ix, iy
- ! print*, SIZE(fyx), SIZE(buffer_1)
+ REAL(xp), DIMENSION(2*NY_halved,local_nx_) :: tmp_yx_1, tmp_yx_2
+ COMPLEX(xp), DIMENSION(NY_halved+1,local_nx_) :: tmp_kyx
+ INTEGER :: ix, iy, iky
+ ! Fill buffer
DO ix = 1,local_nx_
DO iy = 1,2*NY_halved
- buffer_1(iy,ix) = fyx(iy,ix)
+ tmp_yx_1(iy,ix) = fyx(iy,ix)
ENDDO
ENDDO
- ! print*, SUM(fyx), SUM(buffer_1)
- CALL FFT_yx_to_kyx(buffer_1,fkyx)
- fkyx = fkyx*inv_ExB_NL_factor
- CALL iFFT_kyx_to_yx(fkyx,buffer_2)
+ ! Fourier real to complex on the second buffer (first buffer is now unusable)
+ CALL FFT_yx_to_kyx(tmp_yx_1,tmp_kyx)
+ ! Treat the result with the ExB NL factor
+ DO iky = 1,NY_halved
+ DO ix = 1,local_nx_
+ 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)
+ ! Fill the result array with normalization of iFFT
DO ix = 1,local_nx_
DO iy = 1,2*NY_halved
- fyx(iy,ix) = buffer_2(iy,ix)*inv_Ny_
+ fyx(iy,ix) = tmp_yx_2(iy,ix)*inv_Ny_
ENDDO
ENDDO
- ! print*, SUM(fyx)
- ! stop
END SUBROUTINE apply_inv_ExB_NL_factor
!******************************************************************************!
diff --git a/src/grid_mod.F90 b/src/grid_mod.F90
index 823c3de6496025290d5f21742eb693af3580a6f7..c2ef3f0e84342f35f7b795b807d700f394081921 100644
--- a/src/grid_mod.F90
+++ b/src/grid_mod.F90
@@ -112,11 +112,6 @@ MODULE grid
! 1D Antialiasing arrays (2/3 rule)
REAL(xp), DIMENSION(:), ALLOCATABLE, PUBLIC,PROTECTED :: AA_x
REAL(xp), DIMENSION(:), ALLOCATABLE, PUBLIC,PROTECTED :: AA_y
- ! 2D fft routines or 1D methods (for ExBshear simulations, 1D is required)
- ! The 2D fft is using fftw mpi optimization
- ! The 1D is not using mpi and does a data transfer with redundant computations
- ! (each process computes the convolution)
- LOGICAL, PUBLIC, PROTECTED :: FFT2D = .TRUE.
! Public Functions
PUBLIC :: init_grids_data, set_grids, update_grids
@@ -134,7 +129,7 @@ CONTAINS
USE prec_const
IMPLICIT NONE
INTEGER :: lun = 90 ! File duplicated from STDIN
- NAMELIST /GRID/ pmax, jmax, Nx, Lx, Ny, Ly, Nz, SG, Nexc, FFT2D
+ NAMELIST /GRID/ pmax, jmax, Nx, Lx, Ny, Ly, Nz, SG, Nexc
READ(lun,grid)
IF(Nz .EQ. 1) & ! overwrite SG option if Nz = 1 for safety of use
SG = .FALSE.
@@ -202,7 +197,7 @@ CONTAINS
!! Parallel distribution of kx ky grid
IF (LINEARITY .NE. 'linear') THEN ! we let FFTW distribute if we use it
IF (my_id .EQ. 0) write(*,*) 'FFTW3 y-grid distribution'
- CALL init_grid_distr_and_plans(FFT2D,Nx,Ny,comm_ky,local_nx_ptr,local_nx_ptr_offset,local_nky_ptr,local_nky_ptr_offset)
+ CALL init_grid_distr_and_plans(Nx,Ny,comm_ky,local_nx_ptr,local_nx_ptr_offset,local_nky_ptr,local_nky_ptr_offset)
ELSE ! otherwise we distribute equally
IF (my_id .EQ. 0) write(*,*) 'Manual y-grid distribution'
! balanced distribution among the processes
@@ -453,7 +448,6 @@ CONTAINS
! indexation (|1 2 3||1 2 3|... local_nky|)
IF(Ny .EQ. 1) THEN
kyarray(iky) = deltaky
- kyarray(iky) = iky-1
kyarray_full(iky) = deltaky
SINGLE_KY = .TRUE.
ELSE
@@ -543,7 +537,7 @@ CONTAINS
ERROR STOP "Gyacomo is safer with an even Kx number"
ENDIF
! Orszag 2/3 filter
- two_third_kxmax = 2._xp/3._xp*kx_max;
+ two_third_kxmax = 2._xp/3._xp*(kx_max-deltakx);
! Antialiasing filter
DO iky = 1,local_nky
DO ikx = 1,local_nkx
@@ -673,11 +667,11 @@ CONTAINS
REAL(xp), DIMENSION(local_nky),INTENT(IN) :: dkx_ExB ! ExB correction dkx = gamma_E ky dtshift
REAL(xp), DIMENSION(local_nz+ngz,nzgrid), INTENT(IN) :: gxx,gxy,gyy,inv_hatB2
INTEGER :: eo,iz,iky,ikx
- REAL(xp) :: kx, ky, skp2
+ REAL(xp) :: kx, ky
! Update the kx grid
- DO iky = 1,local_nky
- DO ikx = 1,total_Nkx
- kxarray(iky,ikx) = kxarray0(ikx) + dkx_ExB(iky)
+ DO ikx = 1,total_Nkx
+ DO iky = 1,local_nky
+ kxarray(iky,ikx) = kxarray0(ikx) - dkx_ExB(iky)
ENDDO
ENDDO
! Update the kperp grid
diff --git a/src/initial_mod.F90 b/src/initial_mod.F90
index d02af4b8af855d7dc3938dc2d5447bd3f7380d3e..0a0c5c700102c90f3a3fdbedfafcca9df5a74402 100644
--- a/src/initial_mod.F90
+++ b/src/initial_mod.F90
@@ -46,7 +46,6 @@ CONTAINS
USE ghosts, ONLY: update_ghosts_moments, update_ghosts_EM
USE restarts, ONLY: load_moments, job2load
USE processing, ONLY: compute_fluid_moments, compute_radial_heatflux, compute_radial_transport
- USE model, ONLY: LINEARITY
USE nonlinear, ONLY: compute_Sapj, nonlinear_init
IMPLICIT NONE
CALL set_updatetlevel(1)
@@ -266,19 +265,10 @@ CONTAINS
!!!!!!! Initialize a single mode in the gyrocenter density
!******************************************************************************!
SUBROUTINE init_single_mode
- USE grid, ONLY: local_na, local_np, local_nj, total_nkx, local_nky, local_nz,&
- ngp, ngj, ngz, iky0, kxarray, kyarray, parray, jarray, &
- contains_ky0, deltakx, deltaky
USE fields, ONLY: moments
USE prec_const, ONLY: xp
- USE model, ONLY: LINEARITY
USE parallel, ONLY: my_id
IMPLICIT NONE
-
- REAL(xp) :: amplitude
- INTEGER :: ia,ip,ij,ikx,iky,iz, ikxi, ikyi
- INTEGER, DIMENSION(12) :: iseedarr
-
moments = 0._xp
IF (my_id .EQ. 0) THEN
moments(:,:,:,2,2,:,:) = init_amp
@@ -470,7 +460,7 @@ CONTAINS
USE fields, ONLY: moments
USE prec_const, ONLY: xp, pi
USE model, ONLY: LINEARITY
- USE geometry, ONLY: Jacobian, iInt_Jacobian, shear
+ USE geometry, ONLY: Jacobian, iInt_Jacobian
IMPLICIT NONE
REAL(xp) :: kx, ky, sigma_z, z
INTEGER :: ia,iky,ikx,iz,ip,ij
@@ -548,8 +538,8 @@ CONTAINS
USE basic, ONLY: maindir
USE grid, ONLY: local_na, local_np, local_nj, total_nkx, local_nky, local_nz,&
local_nkx_offset, local_nky_offset, kxarray, kyarray, &
- ngp, ngj, ngz, iky0, ikx0, parray, jarray,&
- deltakx, deltaky, contains_ky0, contains_kx0,&
+ ngp, ngj, ngz, parray, jarray,&
+ deltakx, deltaky,&
AA_x, AA_y
USE fields, ONLY: moments
USE prec_const, ONLY: xp, imagu
diff --git a/src/nonlinear_mod.F90 b/src/nonlinear_mod.F90
index 7eea23f6defa0ab8dd4ed3246dee89dda3c21405..946304ac4b328bedd15933932a3eb3eae7126f52 100644
--- a/src/nonlinear_mod.F90
+++ b/src/nonlinear_mod.F90
@@ -100,7 +100,7 @@ SUBROUTINE compute_nonlinear
! this function adds its result to bracket_sum_r
CALL poisson_bracket_and_sum( kyarray,kxarray,inv_Ny,inv_Nx,AA_y,AA_x,&
local_nky,total_nkx,F_cmpx,G_cmpx,&
- ExB, ExB_NL_factor, dkx_ExB, bracket_sum_r)
+ ExB, ExB_NL_factor, bracket_sum_r)
!-----------!! ELECTROMAGNETIC CONTRIBUTION -sqrt(tau)/sigma*{Sum_s dnjs [sqrt(p+1)Nap+1s + sqrt(p)Nap-1s], Kernel psi}
IF(EM) THEN
! First convolution terms
@@ -116,7 +116,7 @@ SUBROUTINE compute_nonlinear
! this function adds its result to bracket_sum_r
CALL poisson_bracket_and_sum( kyarray,kxarray,inv_Ny,inv_Nx,AA_y,AA_x,&
local_nky,total_nkx,F_cmpx,G_cmpx,&
- ExB, ExB_NL_factor, dkx_ExB, bracket_sum_r)
+ ExB, ExB_NL_factor, bracket_sum_r)
ENDIF
ENDDO n
! Apply the ExB shearing rate factor before going back to k-space
diff --git a/src/numerics_mod.F90 b/src/numerics_mod.F90
index 6ec4d1de2ea60ef855242cf3f7923ad5c05a009d..3a1e4a98874382453c358d3af26de72ac0dd7e7d 100644
--- a/src/numerics_mod.F90
+++ b/src/numerics_mod.F90
@@ -63,6 +63,7 @@ SUBROUTINE build_dv4Hp_table
!we scale it w.r.t. to the max degree since
!D_4^{v}\sim (\Delta v/2)^4 and \Delta v \sim 2pi/kvpar = pi/\sqrt{2P}
! dv4_Hp_coeff = dv4_Hp_coeff*(1._xp/2._xp/SQRT(REAL(pmax,xp)))**4
+ IF(pmax .GT. 0) &
dv4_Hp_coeff = dv4_Hp_coeff*(PI/2._xp/SQRT(2._xp*REAL(pmax,xp)))**4
END SUBROUTINE build_dv4Hp_table
!******************************************************************************!
diff --git a/src/stepon.F90 b/src/stepon.F90
index 232d7d68febadcd299e4d73351a7335c6b2da807..ee5fd1da84dbf9c9c5ab0300e3970993bdd279b0 100644
--- a/src/stepon.F90
+++ b/src/stepon.F90
@@ -2,14 +2,13 @@ SUBROUTINE stepon
! Advance one time step, (num_step=4 for Runge Kutta 4 scheme)
USE advance_field_routine, ONLY: advance_time_level, advance_moments
USE basic, ONLY: nlend, chrono_advf, chrono_pois,&
- chrono_chck, chrono_clos, chrono_ghst, chrono_ExBs,&
+ chrono_chck, chrono_clos, chrono_ghst,&
start_chrono, stop_chrono
USE closure, ONLY: apply_closure_model
USE ghosts, ONLY: update_ghosts_moments, update_ghosts_EM
use mpi, ONLY: MPI_COMM_WORLD
USE time_integration, ONLY: ntimelevel
USE prec_const, ONLY: xp
- USE ExB_shear_flow, ONLY: Apply_ExB_shear_flow, Update_ExB_shear_flow
#ifdef TEST_SVD
USE CLA, ONLY: test_svd,filter_sv_moments_ky_pj
#endif
@@ -18,17 +17,7 @@ SUBROUTINE stepon
INTEGER :: num_step, ierr
LOGICAL :: mlend
- ! Update the ExB backg. shear flow before the step
- ! This call includes :
- ! - the kx grid update
- ! - the kernel, poisson op. and ampere op update
- ! - kx-shift of the fields at required ky values
- ! - the ExB shear value (s(ky)) update for the next time step
- CALL start_chrono(chrono_ExBs)
- CALL Apply_ExB_shear_flow
- CALL stop_chrono(chrono_ExBs)
-
- DO num_step=1,ntimelevel ! eg RK4 compute successively k1, k2, k3, k4
+ SUBSTEPS:DO num_step=1,ntimelevel ! eg RK4 compute successively k1, k2, k3, k4
!----- BEFORE: All fields+ghosts are updated for step = n
! Compute right hand side from current fields
! N_rhs(N_n, nadia_n, phi_n, S_n, Tcoll_n)
@@ -80,12 +69,7 @@ SUBROUTINE stepon
CALL MPI_BARRIER(MPI_COMM_WORLD,ierr)
!----- AFTER: All fields are updated for step = n+1
- END DO
-
- ! Update the ExB shear flow for the next step
- CALL start_chrono(chrono_ExBs)
- CALL Update_ExB_shear_flow
- CALL stop_chrono(chrono_ExBs)
+ END DO SUBSTEPS
CONTAINS
!!!! Basic structure to simplify stepon