diff --git a/matlab/setup.m b/matlab/setup.m
index aa8683ce69ffce482afbf7f1c38030f292c3da9b..b972a0c2deb4ed0bfd53a160c8dd28137efaafc4 100644
--- a/matlab/setup.m
+++ b/matlab/setup.m
@@ -18,8 +18,8 @@ if SG; GRID.SG = '.true.'; else; GRID.SG = '.false.';end;
MODEL.CO = CO; % Collision operator (0 : L.Bernstein, -1 : Full Coulomb, -2 : Dougherty)
MODEL.CLOS = CLOS;
MODEL.NL_CLOS = NL_CLOS;
-if NON_LIN; MODEL.NON_LIN = '.true.'; else; MODEL.NON_LIN = '.false.';end;
-MODEL.KIN_E = KIN_E;
+MODEL.NON_LIN = NON_LIN;
+MODEL.KIN_E = KIN_E;
if KIN_E; MODEL.KIN_E = '.true.'; else; MODEL.KIN_E = '.false.';end;
MODEL.mu = MU;
MODEL.mu_p = MU_P;
diff --git a/src/array_mod.F90 b/src/array_mod.F90
index b63b4ec08825faea0b5d9787b628171dcc7477b1..5a9f101731cbef2f7e4a21999814240cf374e2f0 100644
--- a/src/array_mod.F90
+++ b/src/array_mod.F90
@@ -11,6 +11,11 @@ MODULE array
COMPLEX(dp), DIMENSION(:,:,:,:,:), ALLOCATABLE :: nadiab_moments_e
COMPLEX(dp), DIMENSION(:,:,:,:,:), ALLOCATABLE :: nadiab_moments_i
+ ! Arrays to store the initial zonal modes (semi linear simulation)
+ COMPLEX(dp), DIMENSION(:,:,:,:), ALLOCATABLE :: moments_e_ZF
+ COMPLEX(dp), DIMENSION(:,:,:,:), ALLOCATABLE :: moments_i_ZF
+ COMPLEX(dp), DIMENSION(:,:), ALLOCATABLE :: phi_ZF
+
! Non linear term array (ip,ij,ikx,iky,iz)
COMPLEX(dp), DIMENSION(:,:,:,:,:), ALLOCATABLE :: Sepj ! electron
COMPLEX(dp), DIMENSION(:,:,:,:,:), ALLOCATABLE :: Sipj ! ion
diff --git a/src/auxval.F90 b/src/auxval.F90
index c177f1f85456bb448b3b6f818562a289854c6c3e..241d59753ac379687228077267509d93295e49f9 100644
--- a/src/auxval.F90
+++ b/src/auxval.F90
@@ -14,7 +14,7 @@ subroutine auxval
INTEGER :: irows,irowe, irow, icol, i_
IF (my_id .EQ. 0) WRITE(*,*) '=== Set auxiliary values ==='
- IF (NON_LIN) THEN
+ IF (NON_LIN .GT. 0) THEN
CALL init_grid_distr_and_plans(Nx,Ny)
ELSE
CALL init_1Dgrid_distr
@@ -41,7 +41,7 @@ subroutine auxval
CALL evaluate_poisson_op ! precompute the kernels
- IF ( NON_LIN ) THEN;
+ IF ( NON_LIN .GT. 0 ) THEN;
CALL build_dnjs_table ! precompute the Laguerre nonlin product coeffs
ENDIF
diff --git a/src/calculus_mod.F90 b/src/calculus_mod.F90
index d4ae848add5ff015324d0fa0af940a1ee3cc5686..dd2eea85c0a13c960916318e60c138d2d5a3a8d6 100644
--- a/src/calculus_mod.F90
+++ b/src/calculus_mod.F90
@@ -12,7 +12,7 @@ MODULE calculus
(/ onetwentyfourth,-nineeighths, nineeighths,-onetwentyfourth /) ! fd4 odd to even stencil
REAL(dp), dimension(-1:2) :: dz_e2o = &
(/ onetwentyfourth,-nineeighths, nineeighths,-onetwentyfourth /) ! fd4 odd to even stencil
-
+
PUBLIC :: simpson_rule_z, interp_z, grad_z
CONTAINS
@@ -26,26 +26,30 @@ SUBROUTINE grad_z(target,f,ddzf)
INTEGER, INTENT(IN) :: target
COMPLEX(dp),dimension( izs:ize ), intent(in) :: f
COMPLEX(dp),dimension( izs:ize ), intent(out) :: ddzf
- IF(SG) THEN
- IF(TARGET .EQ. 0) THEN
- CALL grad_z_o2e(f,ddzf)
- ELSE
- CALL grad_z_e2o(f,ddzf)
+ IF(Nz .GT. 3) THEN ! Cannot apply four points stencil on less than four points grid
+ IF(SG) THEN
+ IF(TARGET .EQ. 0) THEN
+ CALL grad_z_o2e(f,ddzf)
+ ELSE
+ CALL grad_z_e2o(f,ddzf)
+ ENDIF
+ ELSE ! No staggered grid -> usual centered finite differences
+ DO iz = 3,Nz-2
+ ddzf(iz) = dz_usu(-2)*f(iz-2) + dz_usu(-1)*f(iz-1) &
+ +dz_usu( 1)*f(iz+1) + dz_usu( 2)*f(iz+2)
+ ENDDO
+ ! Periodic boundary conditions
+ ddzf(Nz-1) = dz_usu(-2)*f(Nz-3) + dz_usu(-1)*f(Nz-2) &
+ +dz_usu(+1)*f(Nz ) + dz_usu(+2)*f(1 )
+ ddzf(Nz ) = dz_usu(-2)*f(Nz-2) + dz_usu(-1)*f(Nz-1) &
+ +dz_usu(+1)*f(1 ) + dz_usu(+2)*f(2 )
+ ddzf(1 ) = dz_usu(-2)*f(Nz-1) + dz_usu(-1)*f(Nz ) &
+ +dz_usu(+1)*f(2 ) + dz_usu(+2)*f(3 )
+ ddzf(2 ) = dz_usu(-2)*f(Nz ) + dz_usu(-1)*f(1 ) &
+ +dz_usu(+1)*f(3 ) + dz_usu(+2)*f(4)
ENDIF
- ELSE ! No staggered grid -> usual centered finite differences
- DO iz = 3,Nz-2
- ddzf(iz) = dz_usu(-2)*f(iz-2) + dz_usu(-1)*f(iz-1) &
- +dz_usu( 1)*f(iz+1) + dz_usu( 2)*f(iz+2)
- ENDDO
- ! Periodic boundary conditions
- ddzf(Nz-1) = dz_usu(-2)*f(Nz-3) + dz_usu(-1)*f(Nz-2) &
- +dz_usu(+1)*f(Nz ) + dz_usu(+2)*f(1 )
- ddzf(Nz ) = dz_usu(-2)*f(Nz-2) + dz_usu(-1)*f(Nz-1) &
- +dz_usu(+1)*f(1 ) + dz_usu(+2)*f(2 )
- ddzf(1 ) = dz_usu(-2)*f(Nz-1) + dz_usu(-1)*f(Nz ) &
- +dz_usu(+1)*f(2 ) + dz_usu(+2)*f(3 )
- ddzf(2 ) = dz_usu(-2)*f(Nz ) + dz_usu(-1)*f(1 ) &
- +dz_usu(+1)*f(3 ) + dz_usu(+2)*f(4)
+ ELSE
+ ddzf = 0._dp
ENDIF
END SUBROUTINE grad_z
diff --git a/src/collision_mod.F90 b/src/collision_mod.F90
index e81b6ed3d9cf085b3fccebc6e36b709ca43ae760..c82024512f1ef1a9e5741e006f6a3725870b8b60 100644
--- a/src/collision_mod.F90
+++ b/src/collision_mod.F90
@@ -541,8 +541,8 @@ CONTAINS
CALL getarr(fid, '/coordkperp', kp_grid_mat)
! check that we have enough kperps mat
- kp_max = SQRT(MAXVAL(kxarray)**2+MAXVAL(kyarray)**2)
- IF (NON_LIN) THEN
+ WRITE(*,*) kp_max
+ IF (NON_LIN .GT. 0) THEN
IF ( (kp_grid_mat(nkp_mat) .LT. 2./3.*kp_max) .AND. (my_id .EQ. 0)) WRITE(*,*) '!! Matrix kperp grid too small !!'
ELSE
IF ( (kp_grid_mat(nkp_mat) .LT. kp_max) .AND. (my_id .EQ. 0)) WRITE(*,*) '!! Matrix kperp grid too small !!'
@@ -563,7 +563,7 @@ CONTAINS
DO ikp = ikps_C,ikpe_C ! Loop over everz kperp values
! we put zeros if kp>2/3kpmax because thoses frequencies are filtered through AA
- IF( (kp_grid_mat(ikp) .GT. two_third_kpmax) .AND. NON_LIN) THEN
+ IF( (kp_grid_mat(ikp) .GT. two_third_kpmax) .AND. (NON_LIN .GT. 0)) THEN
CiepjT_kp(:,:,ikp) = 0._dp
CiepjF_kp(:,:,ikp) = 0._dp
CeipjT_kp(:,:,ikp) = 0._dp
@@ -696,8 +696,9 @@ CONTAINS
IF (my_id .EQ. 0 ) WRITE(*,*) '...Interpolation from matrices kperp to simulation kx,ky...'
DO ikx = ikxs,ikxe
DO iky = ikys,ikye
+ ! Check for nonlinear case if we are in the anti aliased domain or the filtered one
kperp_sim = kparray(ikx,iky,1,0) ! current simulation kperp
-
+ IF( (NON_LIN .EQ. 0) .OR. (kperp_sim .LT. two_third_kpmax) ) THEN
! Find the interval in kp grid mat where kperp_sim is contained
! Loop over the whole kp mat grid to find the smallest kperp that is
! larger than the current kperp_sim (brute force...)
@@ -711,9 +712,9 @@ CONTAINS
ikp_next = ikp_mat !index of k1 (larger than kperp_sim thanks to previous loop)
ikp_prev = ikp_mat - 1 !index of k0 (smaller neighbour to interpolate inbetween)
! write(*,*) kp_grid_mat(ikp_prev), '<', kperp_sim, '<', kp_grid_mat(ikp_next)
- if ( (kp_grid_mat(ikp_prev) .GT. kperp_sim) .OR. (kp_grid_mat(ikp_next) .LT. kperp_sim) )&
- write(*,*) 'Warning, linear interp of collision matrix failed!!'
-
+ if ( (kp_grid_mat(ikp_prev) .GT. kperp_sim) .OR. (kp_grid_mat(ikp_next) .LT. kperp_sim) ) THEN
+ write(*,*) 'Warning, linear interp of collision matrix failed!! '
+ ENDIF
! 0->1 variable for linear interp, i.e. zero2one = (k-k0)/(k1-k0)
zerotoone = (kperp_sim - kp_grid_mat(ikp_prev))/(kp_grid_mat(ikp_next) - kp_grid_mat(ikp_prev))
@@ -724,6 +725,14 @@ CONTAINS
Ciipj (:,:,ikx,iky) = (Ciipj__kp(:,:,ikp_next) - Ciipj__kp(:,:,ikp_prev))*zerotoone + Ciipj__kp(:,:,ikp_prev)
CiepjT(:,:,ikx,iky) = (CiepjT_kp(:,:,ikp_next) - CiepjT_kp(:,:,ikp_prev))*zerotoone + CiepjT_kp(:,:,ikp_prev)
CiepjF(:,:,ikx,iky) = (CiepjF_kp(:,:,ikp_next) - CiepjF_kp(:,:,ikp_prev))*zerotoone + CiepjF_kp(:,:,ikp_prev)
+ ELSE ! Out of anti anti aliased domain
+ Ceepj (:,:,ikx,iky) = 0._dp
+ CeipjT(:,:,ikx,iky) = 0._dp
+ CeipjF(:,:,ikx,iky) = 0._dp
+ Ciipj (:,:,ikx,iky) = 0._dp
+ CiepjT(:,:,ikx,iky) = 0._dp
+ CiepjF(:,:,ikx,iky) = 0._dp
+ ENDIF
ENDDO
ENDDO
ELSE ! DK -> No kperp dep, copy simply to final collision matrices
diff --git a/src/compute_Sapj.F90 b/src/compute_Sapj.F90
index 77c12b8a19b1dceb8ea5422d55b2e27950d345b8..e7a4a96c0859246336ab727d6c6961ae7f210066 100644
--- a/src/compute_Sapj.F90
+++ b/src/compute_Sapj.F90
@@ -3,14 +3,16 @@ SUBROUTINE compute_Sapj
!! In real space Sapj ~ b*(grad(phi) x grad(g)) which in moments in fourier becomes
!! Sapj = Sum_n (ikx Kn phi)#(iky Sum_s d_njs Naps) - (iky Kn phi)#(ikx Sum_s d_njs Naps)
!! where # denotes the convolution.
- USE array, ONLY : dnjs, Sepj, Sipj, kernel_i, kernel_e
+ USE array, ONLY : dnjs, Sepj, Sipj, kernel_i, kernel_e,&
+ moments_e_ZF, moments_i_ZF, phi_ZF
+ USE initial_par, ONLY : WIPE_ZF
USE basic
USE fourier
- USE fields!, ONLY : phi, moments_e, moments_i
+ USE fields, ONLY : phi, moments_e, moments_i
USE grid
USE model
USE prec_const
- USE time_integration!, ONLY : updatetlevel
+ USE time_integration, ONLY : updatetlevel
IMPLICIT NONE
INCLUDE 'fftw3-mpi.f03'
@@ -22,9 +24,26 @@ SUBROUTINE compute_Sapj
INTEGER :: in, is, p_int, j_int
INTEGER :: nmax, smax ! Upper bound of the sums
REAL(dp):: kx, ky, kerneln
+
! Execution time start
CALL cpu_time(t0_Sapj)
+ IF(NON_LIN .EQ. 1) THEN
+ CALL compute_non_linear
+ ELSEIF(NON_LIN .EQ. 0.5) THEN
+ CALL compute_semi_linear
+ ELSE
+ Sepj = 0._dp; Sipj = 0._dp
+ ENDIF
+
+ ! Execution time END
+ CALL cpu_time(t1_Sapj)
+ tc_Sapj = tc_Sapj + (t1_Sapj - t0_Sapj)
+
+CONTAINS
+
+SUBROUTINE compute_non_linear
+ IMPLICIT NONE
!!!!!!!!!!!!!!!!!!!! ELECTRON non linear term computation (Sepj)!!!!!!!!!!
IF(KIN_E) THEN
zloope: DO iz = izs,ize
@@ -198,8 +217,204 @@ zloopi: DO iz = izs,ize
ENDDO ploopi
ENDDO zloopi
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+END SUBROUTINE compute_non_linear
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
+! Semi linear computation : Only NZ-ZF convolutions are kept
+SUBROUTINE compute_semi_linear
+ IMPLICIT NONE
+ ! Update the ZF modes if not frozen
+ IF (WIPE_ZF .NE. -2) THEN
+ moments_e_ZF(ips_e:ipe_e,ijs_e:ije_e,ikxs:ikxe,izs:ize) = moments_e(ips_e:ipe_e,ijs_e:ije_e,ikxs:ikxe,iky_0,izs:ize,updatetlevel)
+ moments_i_ZF(ips_i:ipe_i,ijs_i:ije_i,ikxs:ikxe,izs:ize) = moments_i(ips_i:ipe_i,ijs_i:ije_i,ikxs:ikxe,iky_0,izs:ize,updatetlevel)
+ phi_ZF(ikxs:ikxe,izs:ize) = phi(ikxs:ikxe,iky_0,izs:ize)
+ ENDIF
+
+ !!!!!!!!!!!!!!!!!!!! ELECTRON semi linear term computation (Sepj)!!!!!!!!!!
+ IF(KIN_E) THEN
+ zloope: DO iz = izs,ize
+ ploope: DO ip = ips_e,ipe_e ! Loop over Hermite moments
+ eo = MODULO(parray_e(ip),2)
+ jloope: DO ij = ijs_e, ije_e ! Loop over Laguerre moments
+ j_int=jarray_e(ij)
+ real_data_c = 0._dp ! initialize sum over real nonlinear term
+ ! Set non linear sum truncation
+ IF (NL_CLOS .EQ. -2) THEN
+ nmax = Jmaxe
+ ELSEIF (NL_CLOS .EQ. -1) THEN
+ nmax = Jmaxe-(ij-1)
+ ELSE
+ nmax = NL_CLOS
+ ENDIF
+
+ nloope: DO in = 1,nmax+1 ! Loop over laguerre for the sum
+ ! Build the terms to convolve
+ kxloope: DO ikx = ikxs,ikxe ! Loop over kx
+ kyloope: DO iky = ikys,ikye ! Loop over ky
+ kx = kxarray(ikx)
+ ky = kyarray(iky)
+ kerneln = kernel_e(in, ikx, iky, iz, eo)
+
+ ! Zonal terms (=0 for all ky not 0)
+ Fx_cmpx(ikx,iky) = 0._dp
+ Gx_cmpx(ikx,iky) = 0._dp
+ IF(iky .EQ. iky_0) THEN
+ Fx_cmpx(ikx,iky) = imagu*kx* phi_ZF(ikx,iz) * kerneln
+ smax = MIN( (in-1)+(ij-1), jmaxe );
+ DO is = 1, smax+1 ! sum truncation on number of moments
+ Gx_cmpx(ikx,iky) = Gx_cmpx(ikx,iky) + &
+ dnjs(in,ij,is) * moments_e_ZF(ip,is,ikx,iz)
+ ENDDO
+ Gx_cmpx(ikx,iky) = imagu*kx*Gx_cmpx(ikx,iky)
+ ENDIF
+
+ ! NZ terms
+ Fy_cmpx(ikx,iky) = imagu*ky* phi(ikx,iky,iz) * kerneln
+ Gy_cmpx(ikx,iky) = 0._dp ! initialization of the sum
+ smax = MIN( (in-1)+(ij-1), jmaxe );
+ DO is = 1, smax+1 ! sum truncation on number of moments
+ Gy_cmpx(ikx,iky) = Gy_cmpx(ikx,iky) + &
+ dnjs(in,ij,is) * moments_e(ip,is,ikx,iky,iz,updatetlevel)
+ ENDDO
+ Gy_cmpx(ikx,iky) = imagu*ky*Gy_cmpx(ikx,iky)
+
+ ENDDO kyloope
+ ENDDO kxloope
+
+ ! First term ddx(phi_ZF) x ddy(f)
+ DO ikx = ikxs, ikxe
+ DO iky = ikys, ikye
+ cmpx_data_f(iky,ikx-local_nkx_offset) = Fx_cmpx(ikx,iky)*AA_x(ikx)*AA_y(iky) !Anti aliasing filter
+ cmpx_data_g(iky,ikx-local_nkx_offset) = Gy_cmpx(ikx,iky)*AA_x(ikx)*AA_y(iky) !Anti aliasing filter
+ ENDDO
+ ENDDO
+
+ call fftw_mpi_execute_dft_c2r(planb, cmpx_data_f, real_data_f)
+ call fftw_mpi_execute_dft_c2r(planb, cmpx_data_g, real_data_g)
+
+ real_data_c = real_data_c + real_data_f/Ny/Nx * real_data_g/Ny/Nx
+
+ ! Second term -dyphi x dxf_ZF
+ DO ikx = ikxs, ikxe
+ DO iky = ikys, ikye
+ cmpx_data_f(iky,ikx-local_nkx_offset) = Fy_cmpx(ikx,iky)*AA_x(ikx)*AA_y(iky) !Anti aliasing filter
+ cmpx_data_g(iky,ikx-local_nkx_offset) = Gx_cmpx(ikx,iky)*AA_x(ikx)*AA_y(iky) !Anti aliasing filter
+ ENDDO
+ ENDDO
+
+ call fftw_mpi_execute_dft_c2r(planb, cmpx_data_f, real_data_f)
+ call fftw_mpi_execute_dft_c2r(planb, cmpx_data_g, real_data_g)
+
+ real_data_c = real_data_c - real_data_f/Ny/Nx * real_data_g/Ny/Nx
+
+ ENDDO nloope
+
+ ! Put the real nonlinear product into k-space
+ call fftw_mpi_execute_dft_r2c(planf, real_data_c, cmpx_data_c)
+
+ ! Retrieve convolution in input format
+ DO ikx = ikxs, ikxe
+ DO iky = ikys, ikye
+ Sepj(ip,ij,ikx,iky,iz) = cmpx_data_c(iky,ikx-local_nkx_offset)*AA_x(ikx)*AA_y(iky) !Anti aliasing filter
+ ENDDO
+ ENDDO
+ ENDDO jloope
+ ENDDO ploope
+ENDDO zloope
+ENDIF
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
+ !!!!!!!!!!!!!!!!!!!! ION non linear term computation (Sipj)!!!!!!!!!!
+zloopi: DO iz = izs,ize
+ ploopi: DO ip = ips_i,ipe_i ! Loop over Hermite moments
+ eo = MODULO(parray_i(ip),2)
+ ! we check if poly degree is even (eq to index is odd) to spare computation
+ !EVEN_P_i: IF (.TRUE. .OR. (MODULO(ip,2) .EQ. 1) .OR. (.NOT. COMPUTE_ONLY_EVEN_P)) THEN
+ jloopi: DO ij = ijs_i, ije_i ! Loop over Laguerre moments
+ j_int=jarray_i(ij)
+ real_data_c = 0._dp ! initialize sum over real nonlinear term
+ ! Set non linear sum truncation
+ IF (NL_CLOS .EQ. -2) THEN
+ nmax = Jmaxi
+ ELSEIF (NL_CLOS .EQ. -1) THEN
+ nmax = Jmaxi-(ij-1)
+ ELSE
+ nmax = NL_CLOS
+ ENDIF
+
+ nloopi: DO in = 1,nmax+1 ! Loop over laguerre for the sum
+
+ kxloopi: DO ikx = ikxs,ikxe ! Loop over kx
+ kyloopi: DO iky = ikys,ikye ! Loop over ky
+ ! Zonal terms (=0 for all ky not 0)
+ Fx_cmpx(ikx,iky) = 0._dp
+ Gx_cmpx(ikx,iky) = 0._dp
+ IF(iky .EQ. iky_0) THEN
+ Fx_cmpx(ikx,iky) = imagu*kx* phi_ZF(ikx,iz) * kerneln
+ smax = MIN( (in-1)+(ij-1), jmaxe );
+ DO is = 1, smax+1 ! sum truncation on number of moments
+ Gx_cmpx(ikx,iky) = Gx_cmpx(ikx,iky) + &
+ dnjs(in,ij,is) * moments_i_ZF(ip,is,ikx,iz)
+ ENDDO
+ Gx_cmpx(ikx,iky) = imagu*kx*Gx_cmpx(ikx,iky)
+ ENDIF
+
+ ! NZ terms
+ Fy_cmpx(ikx,iky) = imagu*ky* phi(ikx,iky,iz) * kerneln
+ Gy_cmpx(ikx,iky) = 0._dp ! initialization of the sum
+ smax = MIN( (in-1)+(ij-1), jmaxi );
+ DO is = 1, smax+1 ! sum truncation on number of moments
+ Gy_cmpx(ikx,iky) = Gy_cmpx(ikx,iky) + &
+ dnjs(in,ij,is) * moments_i(ip,is,ikx,iky,iz,updatetlevel)
+ ENDDO
+ Gy_cmpx(ikx,iky) = imagu*ky*Gy_cmpx(ikx,iky)
+ ENDDO kyloopi
+ ENDDO kxloopi
+
+ ! First term drphi x dzf
+ DO ikx = ikxs, ikxe
+ DO iky = ikys, ikye
+ cmpx_data_f(iky,ikx-local_nkx_offset) = Fx_cmpx(ikx,iky)*AA_x(ikx)*AA_y(iky)
+ cmpx_data_g(iky,ikx-local_nkx_offset) = Gy_cmpx(ikx,iky)*AA_x(ikx)*AA_y(iky)
+ ENDDO
+ ENDDO
+
+ call fftw_mpi_execute_dft_c2r(planb, cmpx_data_f, real_data_f)
+ call fftw_mpi_execute_dft_c2r(planb, cmpx_data_g, real_data_g)
+
+ real_data_c = real_data_c + real_data_f/Ny/Nx * real_data_g/Ny/Nx
+
+ ! Second term -dzphi x drf
+ DO ikx = ikxs, ikxe
+ DO iky = ikys, ikye
+ cmpx_data_f(iky,ikx-local_nkx_offset) = Fy_cmpx(ikx,iky)*AA_x(ikx)*AA_y(iky)
+ cmpx_data_g(iky,ikx-local_nkx_offset) = Gx_cmpx(ikx,iky)*AA_x(ikx)*AA_y(iky)
+ ENDDO
+ ENDDO
+
+ call fftw_mpi_execute_dft_c2r(planb, cmpx_data_f, real_data_f)
+ call fftw_mpi_execute_dft_c2r(planb, cmpx_data_g, real_data_g)
+
+ real_data_c = real_data_c - real_data_f/Ny/Nx * real_data_g/Ny/Nx
+
+ ENDDO nloopi
+
+ ! Put the real nonlinear product into k-space
+ call fftw_mpi_execute_dft_r2c(planf, real_data_c, cmpx_data_c)
+
+ ! Retrieve convolution in input format
+ DO ikx = ikxs, ikxe
+ DO iky = ikys, ikye
+ Sipj(ip,ij,ikx,iky,iz) = cmpx_data_c(iky,ikx-local_nkx_offset)*AA_x(ikx)*AA_y(iky)
+ ENDDO
+ ENDDO
+ ENDDO jloopi
+ ENDDO ploopi
+ENDDO zloopi
+END SUBROUTINE compute_semi_linear
- ! Execution time END
- CALL cpu_time(t1_Sapj)
- tc_Sapj = tc_Sapj + (t1_Sapj - t0_Sapj)
END SUBROUTINE compute_Sapj
diff --git a/src/geometry_mod.F90 b/src/geometry_mod.F90
index b874797220b86f2e83a59596c503e683473bf03c..d8bc22b7501b92cde94267c4598383a8494e023d 100644
--- a/src/geometry_mod.F90
+++ b/src/geometry_mod.F90
@@ -17,7 +17,7 @@ implicit none
contains
subroutine eval_magnetic_geometry
- ! evalute metrix, elements, jacobian and gradient
+ ! evalute metrix, elementwo_third_kpmaxts, jacobian and gradient
implicit none
REAL(dp) :: kx,ky
COMPLEX(dp), DIMENSION(izs:ize) :: integrant
@@ -46,6 +46,7 @@ contains
ENDDO
ENDDO
ENDDO
+ two_third_kpmax = 2._dp/3._dp * MAXVAL(kparray)
!
! Compute the inverse z integrated Jacobian (useful for flux averaging)
integrant = Jacobian(:,0) ! Convert into complex array
diff --git a/src/grid_mod.F90 b/src/grid_mod.F90
index 0cf0f39471eda127acdb41f3492dcf03e7566060..e382af16b9ff19d04662d65e54b2a94dca9a3319 100644
--- a/src/grid_mod.F90
+++ b/src/grid_mod.F90
@@ -30,7 +30,7 @@ MODULE grid
! For Orszag filter
REAL(dp), PUBLIC, PROTECTED :: two_third_kxmax
REAL(dp), PUBLIC, PROTECTED :: two_third_kymax
- REAL(dp), PUBLIC, PROTECTED :: two_third_kpmax
+ REAL(dp), PUBLIC :: two_third_kpmax
! 1D Antialiasing arrays (2/3 rule)
REAL(dp), DIMENSION(:), ALLOCATABLE, PUBLIC :: AA_x
@@ -285,7 +285,7 @@ CONTAINS
two_third_kxmax = 2._dp/3._dp*deltakx*(Nkx-1)
ALLOCATE(AA_x(ikxs:ikxe))
DO ikx = ikxs,ikxe
- IF ( (kxarray(ikx) .LT. two_third_kxmax) .OR. (.NOT. NON_LIN)) THEN
+ IF ( (kxarray(ikx) .LT. two_third_kxmax) .OR. (NON_LIN .EQ. 0)) THEN
AA_x(ikx) = 1._dp;
ELSE
AA_x(ikx) = 0._dp;
@@ -347,7 +347,7 @@ CONTAINS
ALLOCATE(AA_y(ikys:ikye))
DO iky = ikys,ikye
IF ( ((kyarray(iky) .GT. -two_third_kymax) .AND. &
- (kyarray(iky) .LT. two_third_kymax)) .OR. (.NOT. NON_LIN)) THEN
+ (kyarray(iky) .LT. two_third_kymax)) .OR. (NON_LIN .EQ. 0)) THEN
AA_y(iky) = 1._dp;
ELSE
AA_y(iky) = 0._dp;
diff --git a/src/inital.F90 b/src/inital.F90
index c7e1db9de7d71e996282e7ac7758bb56880486bd..04e4896df26d4158876f41457b9003c92ad62edb 100644
--- a/src/inital.F90
+++ b/src/inital.F90
@@ -4,11 +4,12 @@
SUBROUTINE inital
USE basic
- USE model, ONLY : CO, NON_LIN
+ USE model, ONLY : CO
USE initial_par
USE prec_const
USE time_integration
- USE array, ONLY : Sepj,Sipj
+ USE array, ONLY : moments_e_ZF, moments_i_ZF, phi_ZF
+ USE fields
USE collision
USE closure
USE ghosts
@@ -67,11 +68,14 @@ SUBROUTINE inital
IF (my_id .EQ. 0) WRITE(*,*) 'Computing non adiab moments'
CALL compute_nadiab_moments
+ ! Store Zonal moments for semi linear run or for freezing ZF
+ moments_e_ZF(ips_e:ipe_e,ijs_e:ije_e,ikxs:ikxe,izs:ize) = moments_e(ips_e:ipe_e,ijs_e:ije_e,ikxs:ikxe,iky_0,izs:ize,updatetlevel)
+ moments_i_ZF(ips_i:ipe_i,ijs_i:ije_i,ikxs:ikxe,izs:ize) = moments_i(ips_i:ipe_i,ijs_i:ije_i,ikxs:ikxe,iky_0,izs:ize,updatetlevel)
+ phi_ZF(ikxs:ikxe,izs:ize) = phi(ikxs:ikxe,iky_0,izs:ize)
+
!!!!!! Set Sepj, Sipj and dnjs coeff table !!!!!!
- IF ( NON_LIN ) THEN;
- IF (my_id .EQ. 0) WRITE(*,*) 'Init Sapj'
- CALL compute_Sapj ! compute S_0 = S(phi_0,N_0)
- ENDIF
+ IF (my_id .EQ. 0) WRITE(*,*) 'Init Sapj'
+ CALL compute_Sapj ! compute S_0 = S(phi_0,N_0)
!!!!!! Load the COSOlver collision operator coefficients !!!!!!
IF (ABS(CO) .GT. 1) THEN
@@ -79,6 +83,7 @@ SUBROUTINE inital
! Compute collision
CALL compute_TColl ! compute C_0 = C(N_0)
ENDIF
+
END SUBROUTINE inital
!******************************************************************************!
@@ -147,7 +152,7 @@ SUBROUTINE init_moments
END DO
! Putting to zero modes that are not in the 2/3 Orszag rule
- IF (NON_LIN) THEN
+ IF (NON_LIN .GT. 0) THEN
DO ikx=ikxs,ikxe
DO iky=ikys,ikye
DO iz=izs,ize
@@ -241,7 +246,7 @@ SUBROUTINE init_gyrodens
END DO
! Putting to zero modes that are not in the 2/3 Orszag rule
- IF (NON_LIN) THEN
+ IF (NON_LIN .GT. 0) THEN
DO ikx=ikxs,ikxe
DO iky=ikys,ikye
DO iz=izs,ize
diff --git a/src/memory.F90 b/src/memory.F90
index 79249e7572f80faf763f5fe76210a569a27836a7..f2c9c6c7bed8c9f0d3c063466165aa138daf69df 100644
--- a/src/memory.F90
+++ b/src/memory.F90
@@ -12,7 +12,8 @@ SUBROUTINE memory
IMPLICIT NONE
! Electrostatic potential
- CALL allocate_array(phi, ikxs,ikxe, ikys,ikye, izs,ize)
+ CALL allocate_array( phi, ikxs,ikxe, ikys,ikye, izs,ize)
+ CALL allocate_array( phi_ZF, ikxs,ikxe, izs,ize)
CALL allocate_array(inv_poisson_op, ikxs,ikxe, ikys,ikye, izs,ize)
!Electrons arrays
@@ -24,6 +25,7 @@ SUBROUTINE memory
CALL allocate_array( moments_e, ipsg_e,ipeg_e, ijsg_e,ijeg_e, ikxs,ikxe, ikys,ikye, izs,ize, 1,ntimelevel )
CALL allocate_array( moments_rhs_e, ips_e,ipe_e, ijs_e,ije_e, ikxs,ikxe, ikys,ikye, izs,ize, 1,ntimelevel )
CALL allocate_array( nadiab_moments_e, ipsg_e,ipeg_e, ijsg_e,ijeg_e, ikxs,ikxe, ikys,ikye, izs,ize)
+ CALL allocate_array( moments_e_ZF, ipsg_e,ipeg_e, ijsg_e,ijeg_e, ikxs,ikxe, izs,ize)
CALL allocate_array( TColl_e, ips_e,ipe_e, ijs_e,ije_e , ikxs,ikxe, ikys,ikye, izs,ize)
CALL allocate_array( Sepj, ips_e,ipe_e, ijs_e,ije_e, ikxs,ikxe, ikys,ikye, izs,ize)
CALL allocate_array( xnepj, ips_e,ipe_e, ijs_e,ije_e)
@@ -50,6 +52,7 @@ SUBROUTINE memory
CALL allocate_array( moments_i, ipsg_i,ipeg_i, ijsg_i,ijeg_i, ikxs,ikxe, ikys,ikye, izs,ize, 1,ntimelevel )
CALL allocate_array( moments_rhs_i, ips_i,ipe_i, ijs_i,ije_i, ikxs,ikxe, ikys,ikye, izs,ize, 1,ntimelevel )
CALL allocate_array( nadiab_moments_i, ipsg_i,ipeg_i, ijsg_i,ijeg_i, ikxs,ikxe, ikys,ikye, izs,ize)
+ CALL allocate_array( moments_i_ZF, ipsg_i,ipeg_i, ijsg_i,ijeg_i, ikxs,ikxe, izs,ize)
CALL allocate_array( TColl_i, ips_i,ipe_i, ijs_i,ije_i, ikxs,ikxe, ikys,ikye, izs,ize)
CALL allocate_array( Sipj, ips_i,ipe_i, ijs_i,ije_i, ikxs,ikxe, ikys,ikye, izs,ize)
CALL allocate_array( xnipj, ips_i,ipe_i, ijs_i,ije_i)
diff --git a/src/model_mod.F90 b/src/model_mod.F90
index dac8b664894f6833d129a574e4401b9f91389abd..c6e3d37dfa55b3a3ac52d47e719ab35b2036b13a 100644
--- a/src/model_mod.F90
+++ b/src/model_mod.F90
@@ -9,7 +9,7 @@ MODULE model
INTEGER, PUBLIC, PROTECTED :: CLOS = 0 ! linear truncation method
INTEGER, PUBLIC, PROTECTED :: NL_CLOS = 0 ! nonlinear truncation method
INTEGER, PUBLIC, PROTECTED :: KERN = 0 ! Kernel model
- LOGICAL, PUBLIC, PROTECTED :: NON_LIN = .true. ! To turn on non linear bracket term
+ INTEGER, PUBLIC, PROTECTED :: NON_LIN = 0 ! To turn on non linear bracket term
LOGICAL, PUBLIC, PROTECTED :: KIN_E = .true. ! Simulate kinetic electron (adiabatic otherwise)
REAL(dp), PUBLIC, PROTECTED :: mu = 0._dp ! spatial Hyperdiffusivity coefficient (for num. stability)
REAL(dp), PUBLIC, PROTECTED :: mu_p = 0._dp ! kinetic para hyperdiffusivity coefficient (for num. stability)
diff --git a/src/moments_eq_rhs.F90 b/src/moments_eq_rhs.F90
index ba1fc93e68dd211da1af6981d4d31be12edd1f75..eed44047f87e731a06c29b645d0f115b6cb74020 100644
--- a/src/moments_eq_rhs.F90
+++ b/src/moments_eq_rhs.F90
@@ -127,10 +127,9 @@ SUBROUTINE moments_eq_rhs_e
+ TColl
!! Adding non linearity
- IF ( NON_LIN ) THEN
moments_rhs_e(ip,ij,ikx,iky,iz,updatetlevel) = &
moments_rhs_e(ip,ij,ikx,iky,iz,updatetlevel) - Sepj(ip,ij,ikx,iky,iz)
- ENDIF
+
END DO zloope
END DO kyloope
END DO kxloope
@@ -266,10 +265,8 @@ SUBROUTINE moments_eq_rhs_i
+ TColl
!! Adding non linearity
- IF ( NON_LIN ) THEN
moments_rhs_i(ip,ij,ikx,iky,iz,updatetlevel) = &
moments_rhs_i(ip,ij,ikx,iky,iz,updatetlevel) - Sipj(ip,ij,ikx,iky,iz)
- ENDIF
END DO zloopi
END DO kyloopi
END DO kxloopi
diff --git a/src/srcinfo.h b/src/srcinfo.h
index 9b4b2428e4e58f77736196858f5e960af03cdd4d..5b3f294989538a1cc8233fc1241b30c097cf8000 100644
--- a/src/srcinfo.h
+++ b/src/srcinfo.h
@@ -3,8 +3,8 @@ character(len=40) BRANCH
character(len=20) AUTHOR
character(len=40) EXECDATE
character(len=40) HOST
-parameter (VERSION='95b37e1-dirty')
+parameter (VERSION='0d685e3-dirty')
parameter (BRANCH='master')
parameter (AUTHOR='ahoffman')
-parameter (EXECDATE='Mon Nov 8 11:53:45 CET 2021')
+parameter (EXECDATE='Thu Nov 11 10:19:20 CET 2021')
parameter (HOST ='spcpc606')
diff --git a/src/stepon.F90 b/src/stepon.F90
index a5b65bb5efa875033fd36651faf561057e5cba21..a59d15fabbe0b0aa0d752fe4e3dcfa840767a258 100644
--- a/src/stepon.F90
+++ b/src/stepon.F90
@@ -44,7 +44,7 @@ SUBROUTINE stepon
! Update non adiabatic moments n -> n+1
CALL compute_nadiab_moments
! Update nonlinear term S_n -> S_n+1(phi_n+1,N_n+1)
- IF ( NON_LIN ) CALL compute_Sapj
+ CALL compute_Sapj
! Cancel zonal modes artificially
IF ( WIPE_ZF .EQ. 2) CALL wipe_zonalflow
! Cancel non zonal modes artificially
@@ -63,8 +63,8 @@ SUBROUTINE stepon
! Execution time start
CALL cpu_time(t0_checkfield)
- IF(NON_LIN) CALL anti_aliasing ! ensure 0 mode for 2/3 rule
- IF(NON_LIN) CALL enforce_symmetry ! Enforcing symmetry on kx = 0
+ IF(NON_LIN .GT. 0) CALL anti_aliasing ! ensure 0 mode for 2/3 rule
+ IF(NON_LIN .GT. 0) CALL enforce_symmetry ! Enforcing symmetry on kx = 0
mlend=.FALSE.
IF(.NOT.nlend) THEN