diff --git a/src/advance_field.F90 b/src/advance_field.F90
index 5e8c87d622973eb76990cfba6658d4a3711d49c4..071d1217018b1e1d6b730ab94c8bcff9bddf2592 100644
--- a/src/advance_field.F90
+++ b/src/advance_field.F90
@@ -57,7 +57,7 @@ CONTAINS
USE array
USE grid
use prec_const
- use initial_par, ONLY: WIPE_ZF
+ use initial_par, ONLY: ACT_ON_MODES
IMPLICIT NONE
COMPLEX(dp), DIMENSION ( ikxs:ikxe, ikys:ikye, izs:ize, ntimelevel ) :: f
diff --git a/src/array_mod.F90 b/src/array_mod.F90
index 5a9f101731cbef2f7e4a21999814240cf374e2f0..d77653125be8cccc7688630fe4748bb486a7c77c 100644
--- a/src/array_mod.F90
+++ b/src/array_mod.F90
@@ -11,10 +11,15 @@ 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)
+ ! Arrays to store special initial modes (semi linear simulation)
+ ! Zonal ones (ky=0)
COMPLEX(dp), DIMENSION(:,:,:,:), ALLOCATABLE :: moments_e_ZF
COMPLEX(dp), DIMENSION(:,:,:,:), ALLOCATABLE :: moments_i_ZF
COMPLEX(dp), DIMENSION(:,:), ALLOCATABLE :: phi_ZF
+ ! Entropy modes (kx=0)
+ COMPLEX(dp), DIMENSION(:,:,:,:), ALLOCATABLE :: moments_e_EM
+ COMPLEX(dp), DIMENSION(:,:,:,:), ALLOCATABLE :: moments_i_EM
+ COMPLEX(dp), DIMENSION(:,:), ALLOCATABLE :: phi_EM
! Non linear term array (ip,ij,ikx,iky,iz)
COMPLEX(dp), DIMENSION(:,:,:,:,:), ALLOCATABLE :: Sepj ! electron
diff --git a/src/auxval.F90 b/src/auxval.F90
index 241d59753ac379687228077267509d93295e49f9..1d978124759912634a114414a90f1fe64531fff6 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 .GT. 0) THEN
+ IF (LINEARITY .NE. 'linear') 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 .GT. 0 ) THEN;
+ IF ( LINEARITY .NE. 'linear' ) THEN;
CALL build_dnjs_table ! precompute the Laguerre nonlin product coeffs
ENDIF
diff --git a/src/compute_Sapj.F90 b/src/compute_Sapj.F90
index e7a4a96c0859246336ab727d6c6961ae7f210066..ecf8df5675a40815860124b52c6121ea378f439f 100644
--- a/src/compute_Sapj.F90
+++ b/src/compute_Sapj.F90
@@ -5,7 +5,7 @@ SUBROUTINE compute_Sapj
!! where # denotes the convolution.
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 initial_par, ONLY : ACT_ON_MODES
USE basic
USE fourier
USE fields, ONLY : phi, moments_e, moments_i
@@ -16,6 +16,7 @@ SUBROUTINE compute_Sapj
IMPLICIT NONE
INCLUDE 'fftw3-mpi.f03'
+ COMPLEX(dp), DIMENSION(ikxs:ikxe,ikys:ikye) :: F_cmpx, G_cmpx
COMPLEX(dp), DIMENSION(ikxs:ikxe,ikys:ikye) :: Fx_cmpx, Gy_cmpx
COMPLEX(dp), DIMENSION(ikxs:ikxe,ikys:ikye) :: Fy_cmpx, Gx_cmpx, F_conv_G
REAL(dp), DIMENSION(ixs:ixe,iys:iye) :: fr_real, gz_real
@@ -28,13 +29,19 @@ SUBROUTINE compute_Sapj
! 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
+ SELECT CASE(LINEARITY)
+ CASE ('nonlinear')
+ CALL compute_nonlinear
+ CASE ('ZF_semilin')
+ CALL compute_semi_linear_ZF
+ CASE ('NZ_semilin')
+ CALL compute_semi_linear_NZ
+ CASE ('linear')
+ Sepj = 0._dp; Sipj = 0._dp
+ CASE DEFAULT
+ IF(my_id.EQ.0) write(*,*) '/!\ Linearity not recognized /!\'
+ stop
+ END SELECT
! Execution time END
CALL cpu_time(t1_Sapj)
@@ -42,7 +49,7 @@ SUBROUTINE compute_Sapj
CONTAINS
-SUBROUTINE compute_non_linear
+SUBROUTINE compute_nonlinear
IMPLICIT NONE
!!!!!!!!!!!!!!!!!!!! ELECTRON non linear term computation (Sepj)!!!!!!!!!!
IF(KIN_E) THEN
@@ -51,7 +58,6 @@ SUBROUTINE compute_non_linear
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
@@ -60,69 +66,27 @@ SUBROUTINE compute_non_linear
ELSE
nmax = NL_CLOS
ENDIF
-
+ bracket_sum_r = 0._dp ! initialize sum over real nonlinear term
nloope: DO in = 1,nmax+1 ! Loop over laguerre for the sum
-
- 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)
-
- ! First convolution terms
- Fx_cmpx(ikx,iky) = imagu*kx* phi(ikx,iky,iz) * kerneln
- Fy_cmpx(ikx,iky) = imagu*ky* phi(ikx,iky,iz) * kerneln
- ! Second convolution terms
- Gy_cmpx(ikx,iky) = 0._dp ! initialization of the sum
- Gx_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)
- Gx_cmpx(ikx,iky) = Gx_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)
- Gx_cmpx(ikx,iky) = imagu*kx*Gx_cmpx(ikx,iky)
- ENDDO kyloope
- ENDDO kxloope
-
- ! 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) !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
+ ! First convolution terms
+ F_cmpx(:,:) = phi(:,:,iz) * kernel_e(in, :, :, iz, eo)
+ ! Second convolution terms
+ G_cmpx(:,:) = 0._dp ! initialization of the sum
+ smax = MIN( (in-1)+(ij-1), jmaxe );
+ DO is = 1, smax+1 ! sum truncation on number of moments
+ G_cmpx(:,:) = G_cmpx(:,:) + &
+ dnjs(in,ij,is) * moments_e(ip,is,:,:,iz,updatetlevel)
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) !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
-
+ !/!\ this function add its result to bracket_sum_r (hard to read sorry) /!\
+ CALL poisson_bracket_and_sum(F_cmpx,G_cmpx)
ENDDO nloope
! Put the real nonlinear product into k-space
- call fftw_mpi_execute_dft_r2c(planf, real_data_c, cmpx_data_c)
-
+ call fftw_mpi_execute_dft_r2c(planf, bracket_sum_r, bracket_sum_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
+ Sepj(ip,ij,ikx,iky,iz) = bracket_sum_c(iky,ikx-local_nkx_offset)*AA_x(ikx)*AA_y(iky) !Anti aliasing filter
ENDDO
ENDDO
ENDDO jloope
@@ -135,11 +99,8 @@ ENDIF
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
@@ -148,76 +109,33 @@ zloopi: DO iz = izs,ize
ELSE
nmax = NL_CLOS
ENDIF
-
+ bracket_sum_r = 0._dp ! initialize sum over real nonlinear term
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
- kx = kxarray(ikx)
- ky = kyarray(iky)
- kerneln = kernel_i(in, ikx, iky, iz, eo)
-
- ! First convolution terms
- Fx_cmpx(ikx,iky) = imagu*kx* phi(ikx,iky,iz) * kerneln
- Fy_cmpx(ikx,iky) = imagu*ky* phi(ikx,iky,iz) * kerneln
- ! Second convolution terms
- Gy_cmpx(ikx,iky) = 0._dp ! initialization of the sum
- Gx_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)
- Gx_cmpx(ikx,iky) = Gx_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)
- Gx_cmpx(ikx,iky) = imagu*kx*Gx_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
+ ! First convolution terms
+ F_cmpx(:,:) = phi(:,:,iz) * kernel_i(in, :, :, iz, eo)
+ ! Second convolution terms
+ G_cmpx(:,:) = 0._dp ! initialization of the sum
+ smax = MIN( (in-1)+(ij-1), jmaxi );
+ DO is = 1, smax+1 ! sum truncation on number of moments
+ G_cmpx(:,:) = G_cmpx(:,:) + &
+ dnjs(in,ij,is) * moments_i(ip,is,:,:,iz,updatetlevel)
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
-
+ !/!\ this function add its result to bracket_sum_r (hard to read sorry) /!\
+ CALL poisson_bracket_and_sum(F_cmpx,G_cmpx)
ENDDO nloopi
-
! Put the real nonlinear product into k-space
- call fftw_mpi_execute_dft_r2c(planf, real_data_c, cmpx_data_c)
-
+ call fftw_mpi_execute_dft_r2c(planf, bracket_sum_r, bracket_sum_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)
+ Sipj(ip,ij,ikx,iky,iz) = bracket_sum_c(iky,ikx-local_nkx_offset)*AA_x(ikx)*AA_y(iky)
ENDDO
ENDDO
ENDDO jloopi
ENDDO ploopi
ENDDO zloopi
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-END SUBROUTINE compute_non_linear
+END SUBROUTINE compute_nonlinear
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
@@ -225,15 +143,8 @@ END SUBROUTINE compute_non_linear
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Semi linear computation : Only NZ-ZF convolutions are kept
-SUBROUTINE compute_semi_linear
+SUBROUTINE compute_semi_linear_ZF
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
@@ -241,7 +152,6 @@ SUBROUTINE compute_semi_linear
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
@@ -250,7 +160,7 @@ SUBROUTINE compute_semi_linear
ELSE
nmax = NL_CLOS
ENDIF
-
+ bracket_sum_r = 0._dp ! initialize sum over real nonlinear term
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
@@ -258,20 +168,18 @@ SUBROUTINE compute_semi_linear
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
+ Fx_cmpx(ikx,iky) = imagu*kx* phi(ikx,iky,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)
+ dnjs(in,ij,is) * moments_e(ip,is,ikx,iky,iz,updatetlevel)
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
@@ -281,45 +189,19 @@ SUBROUTINE compute_semi_linear
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
-
+ ! First term df/dx x dg/dy
+ CALL convolve_and_add(Fx_cmpx,Gy_cmpx)
+ ! Second term -df/dy x dg/dx
+ CALL convolve_and_add(-Fy_cmpx,Gx_cmpx)
ENDDO nloope
-
! Put the real nonlinear product into k-space
- call fftw_mpi_execute_dft_r2c(planf, real_data_c, cmpx_data_c)
-
+ call fftw_mpi_execute_dft_r2c(planf, bracket_sum_r, bracket_sum_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
+ Sepj(ip,ij,ikx,iky,iz) = bracket_sum_c(iky,ikx-local_nkx_offset)*AA_x(ikx)*AA_y(iky) !Anti aliasing filter
ENDDO
ENDDO
ENDDO jloope
@@ -332,11 +214,8 @@ ENDIF
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
@@ -345,20 +224,19 @@ zloopi: DO iz = izs,ize
ELSE
nmax = NL_CLOS
ENDIF
-
+ bracket_sum_r = 0._dp ! initialize sum over real nonlinear term
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
+ Fx_cmpx(ikx,iky) = imagu*kx* phi(ikx,iky,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)
+ dnjs(in,ij,is) * moments_i(ip,is,ikx,iky,iz,updatetlevel)
ENDDO
Gx_cmpx(ikx,iky) = imagu*kx*Gx_cmpx(ikx,iky)
ENDIF
@@ -374,47 +252,154 @@ zloopi: DO iz = izs,ize
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
+ CALL convolve_and_add(Fy_cmpx,Gx_cmpx)
+ ! Second term -dzphi x drf
+ CALL convolve_and_add(Fy_cmpx,Gx_cmpx)
+ ENDDO nloopi
+ ! Put the real nonlinear product into k-space
+ call fftw_mpi_execute_dft_r2c(planf, bracket_sum_r, bracket_sum_c)
+ ! Retrieve convolution in input format
+ DO ikx = ikxs, ikxe
+ DO iky = ikys, ikye
+ Sipj(ip,ij,ikx,iky,iz) = bracket_sum_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_ZF
- 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
+! Semi linear computation : Only kx=0*all convolutions are kept
+SUBROUTINE compute_semi_linear_NZ
+ IMPLICIT NONE
+ !!!!!!!!!!!!!!!!!!!! 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)
+ ! 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
+ bracket_sum_r = 0._dp ! initialize sum over real nonlinear term
+ 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)
+ ! All terms
+ Fx_cmpx(ikx,iky) = imagu*kx* phi(ikx,iky,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(ip,is,ikx,iky,iz,updatetlevel)
+ ENDDO
+ Gx_cmpx(ikx,iky) = imagu*kx*Gx_cmpx(ikx,iky)
+ ! Kx = 0 terms
+ Fy_cmpx(ikx,iky) = 0._dp
+ Gy_cmpx(ikx,iky) = 0._dp
+ IF (ikx .EQ. ikx_0) THEN
+ 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)
+ ENDIF
+ ENDDO kyloope
+ ENDDO kxloope
+ ! First term df/dx x dg/dy
+ CALL convolve_and_add(Fx_cmpx,Gy_cmpx)
+ ! Second term -df/dy x dg/dx
+ CALL convolve_and_add(-Fy_cmpx,Gx_cmpx)
+ ENDDO nloope
+ ! Put the real nonlinear product into k-space
+ call fftw_mpi_execute_dft_r2c(planf, bracket_sum_r, bracket_sum_c)
+ ! Retrieve convolution in input format
+ DO ikx = ikxs, ikxe
+ DO iky = ikys, ikye
+ Sepj(ip,ij,ikx,iky,iz) = bracket_sum_c(iky,ikx-local_nkx_offset)*AA_x(ikx)*AA_y(iky) !Anti aliasing filter
ENDDO
+ ENDDO
+ ENDDO jloope
+ ENDDO ploope
+ENDDO zloope
+ENDIF
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
- 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
+ !!!!!!!!!!!!!!!!!!!! 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)
+ jloopi: DO ij = ijs_i, ije_i ! Loop over Laguerre moments
+ j_int=jarray_i(ij)
+ ! 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
+ bracket_sum_r = 0._dp ! initialize sum over real nonlinear term
+ 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) = imagu*kx* phi(ikx,iky,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(ip,is,ikx,iky,iz,updatetlevel)
+ ENDDO
+ Gx_cmpx(ikx,iky) = imagu*kx*Gx_cmpx(ikx,iky)
+ ! Kx = 0 terms
+ Fy_cmpx(ikx,iky) = 0._dp
+ Gy_cmpx(ikx,iky) = 0._dp
+ IF (ikx .EQ. ikx_0) THEN
+ 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)
+ ENDIF
+ ENDDO kyloopi
+ ENDDO kxloopi
+ ! First term drphi x dzf
+ CALL convolve_and_add(Fy_cmpx,Gx_cmpx)
+ ! Second term -dzphi x drf
+ CALL convolve_and_add(Fy_cmpx,Gx_cmpx)
ENDDO nloopi
-
! Put the real nonlinear product into k-space
- call fftw_mpi_execute_dft_r2c(planf, real_data_c, cmpx_data_c)
-
+ call fftw_mpi_execute_dft_r2c(planf, bracket_sum_r, bracket_sum_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)
+ Sipj(ip,ij,ikx,iky,iz) = bracket_sum_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
+END SUBROUTINE compute_semi_linear_NZ
END SUBROUTINE compute_Sapj
diff --git a/src/diagnose.F90 b/src/diagnose.F90
index 1ba38d81a0e333fa1b5757ea43fc87411a70157b..056661f9126b5d9e7a2ebd893d3fda91057fea54 100644
--- a/src/diagnose.F90
+++ b/src/diagnose.F90
@@ -11,6 +11,7 @@ SUBROUTINE diagnose(kstep)
USE time_integration
USE utility
USE prec_const
+ USE collision, ONLY: coll_outputinputs
IMPLICIT NONE
INCLUDE 'srcinfo.h'
@@ -197,6 +198,8 @@ SUBROUTINE diagnose(kstep)
CALL model_outputinputs(fidres, str)
+ CALL coll_outputinputs(fidres, str)
+
CALL initial_outputinputs(fidres, str)
CALL time_integration_outputinputs(fidres, str)
diff --git a/src/fourier_mod.F90 b/src/fourier_mod.F90
index 3c5a232b171c707201e0e4b7be6221ed717fd24e..5b947db46e95fbfcc79853df21afed50341d2d15 100644
--- a/src/fourier_mod.F90
+++ b/src/fourier_mod.F90
@@ -11,16 +11,16 @@ MODULE fourier
PRIVATE
- PUBLIC :: init_grid_distr_and_plans, convolve_2D_F2F, finalize_plans
+ PUBLIC :: init_grid_distr_and_plans, poisson_bracket_and_sum, convolve_and_add, finalize_plans
- real(C_DOUBLE), pointer, PUBLIC :: real_data_f(:,:), real_data_g(:,:), real_data_c(:,:)
- complex(C_DOUBLE_complex), pointer, PUBLIC :: cmpx_data_f(:,:), cmpx_data_g(:,:), cmpx_data_c(:,:)
+ real(C_DOUBLE), pointer, PUBLIC :: real_data_f(:,:), real_data_g(:,:), bracket_sum_r(:,:)
+ complex(C_DOUBLE_complex), pointer, PUBLIC :: cmpx_data_f(:,:), cmpx_data_g(:,:), bracket_sum_c(:,:)
type(C_PTR) :: cdatar_f, cdatar_g, cdatar_c
type(C_PTR) :: cdatac_f, cdatac_g, cdatac_c
type(C_PTR) , PUBLIC :: planf, planb
integer(C_INTPTR_T) :: i, ix, iy
integer(C_INTPTR_T), PUBLIC :: alloc_local_1, alloc_local_2
- integer(C_INTPTR_T) :: NR_, NZ_, NR_halved
+ integer(C_INTPTR_T) :: NX_, NY_, NX_halved
integer :: communicator
! many plan data variables
integer(C_INTPTR_T) :: howmany=9 ! numer of eleemnt of the tensor
@@ -33,39 +33,39 @@ MODULE fourier
IMPLICIT NONE
INTEGER, INTENT(IN) :: Nx,Ny
- NR_ = Nx; NZ_ = Ny
- NR_halved = NR_/2 + 1
+ NX_ = Nx; NY_ = Ny
+ NX_halved = NX_/2 + 1
! communicator = MPI_COMM_WORLD
communicator = comm_kx
!! Complex arrays F, G
! Compute the room to allocate
- alloc_local_1 = fftw_mpi_local_size_2d(NR_halved, NZ_, communicator, local_nkx, local_nkx_offset)
+ alloc_local_1 = fftw_mpi_local_size_2d(NX_halved, NY_, communicator, local_nkx, local_nkx_offset)
! Initalize pointers to this room
cdatac_f = fftw_alloc_complex(alloc_local_1)
cdatac_g = fftw_alloc_complex(alloc_local_1)
cdatac_c = fftw_alloc_complex(alloc_local_1)
! Initalize the arrays with the rooms pointed
- call c_f_pointer(cdatac_f, cmpx_data_f, [NZ_ ,local_nkx])
- call c_f_pointer(cdatac_g, cmpx_data_g, [NZ_ ,local_nkx])
- call c_f_pointer(cdatac_c, cmpx_data_c, [NZ_ ,local_nkx])
+ call c_f_pointer(cdatac_f, cmpx_data_f, [NY_ ,local_nkx])
+ call c_f_pointer(cdatac_g, cmpx_data_g, [NY_ ,local_nkx])
+ call c_f_pointer(cdatac_c, bracket_sum_c, [NY_ ,local_nkx])
!! Real arrays iFFT(F), iFFT(G)
! Compute the room to allocate
- alloc_local_2 = fftw_mpi_local_size_2d(NZ_, NR_halved, communicator, local_nky, local_nky_offset)
+ alloc_local_2 = fftw_mpi_local_size_2d(NY_, NX_halved, communicator, local_nky, local_nky_offset)
! Initalize pointers to this room
cdatar_f = fftw_alloc_real(2*alloc_local_2)
cdatar_g = fftw_alloc_real(2*alloc_local_2)
cdatar_c = fftw_alloc_real(2*alloc_local_2)
! Initalize the arrays with the rooms pointed
- call c_f_pointer(cdatar_f, real_data_f, [2*(NR_/2 + 1),local_nky])
- call c_f_pointer(cdatar_g, real_data_g, [2*(NR_/2 + 1),local_nky])
- call c_f_pointer(cdatar_c, real_data_c, [2*(NR_/2 + 1),local_nky])
+ call c_f_pointer(cdatar_f, real_data_f, [2*(NX_/2 + 1),local_nky])
+ call c_f_pointer(cdatar_g, real_data_g, [2*(NX_/2 + 1),local_nky])
+ call c_f_pointer(cdatar_c, bracket_sum_r, [2*(NX_/2 + 1),local_nky])
! Plan Creation (out-of-place forward and backward FFT)
- planf = fftw_mpi_plan_dft_r2c_2d(NZ_, NR_, real_data_f, cmpx_data_f, communicator, ior(FFTW_MEASURE, FFTW_MPI_TRANSPOSED_OUT))
- planb = fftw_mpi_plan_dft_c2r_2d(NZ_, NR_, cmpx_data_f, real_data_f, communicator, ior(FFTW_MEASURE, FFTW_MPI_TRANSPOSED_IN))
+ planf = fftw_mpi_plan_dft_r2c_2D(NY_, NX_, real_data_f, cmpx_data_f, communicator, ior(FFTW_MEASURE, FFTW_MPI_TRANSPOSED_OUT))
+ planb = fftw_mpi_plan_dft_c2r_2D(NY_, NX_, cmpx_data_f, real_data_f, communicator, ior(FFTW_MEASURE, FFTW_MPI_TRANSPOSED_IN))
if ((.not. c_associated(planf)) .OR. (.not. c_associated(planb))) then
write(*,*) "plan creation error!!"
@@ -75,37 +75,56 @@ MODULE fourier
END SUBROUTINE init_grid_distr_and_plans
- !!! Convolution 2D Fourier to Fourier
- ! - Compute the convolution using the convolution theorem and MKL
- SUBROUTINE convolve_2D_F2F( F_2D, G_2D, C_2D )
+ !!! Compute the poisson bracket of [F,G] to real space
+ ! - Compute the convolution using the convolution theorem
+ SUBROUTINE poisson_bracket_and_sum( F_, G_)
IMPLICIT NONE
-
- COMPLEX(C_DOUBLE_COMPLEX), DIMENSION(ikxs:ikxe,ikys:ikye), INTENT(IN) :: F_2D, G_2D ! input fields
- COMPLEX(C_DOUBLE_COMPLEX), DIMENSION(ikxs:ikxe,ikys:ikye), INTENT(OUT) :: C_2D ! output convolutioned field
-
- do ikx = ikxs, ikxe
- do iky = ikys, ikye
- cmpx_data_f(iky,ikx-local_nkx_offset) = F_2D(ikx,iky)*AA_x(ikx)*AA_y(iky)
- cmpx_data_g(iky,ikx-local_nkx_offset) = G_2D(ikx,iky)*AA_x(ikx)*AA_y(iky)
- end do
- end do
-
+ COMPLEX(C_DOUBLE_COMPLEX), DIMENSION(ikxs:ikxe,ikys:ikye),&
+ INTENT(IN) :: F_, G_ ! input fields
+ ! First term df/dx x dg/dy
+ DO ikx = ikxs, ikxe
+ DO iky = ikys, ikye
+ cmpx_data_f(iky,ikx-local_nkx_offset) = &
+ imagu*kxarray(ikx)*F_(ikx,iky)*AA_x(ikx)*AA_y(iky) !Anti aliasing filter
+ cmpx_data_g(iky,ikx-local_nkx_offset) = &
+ imagu*kyarray(iky)*G_(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)
+ bracket_sum_r = bracket_sum_r + real_data_f*inv_Ny*inv_Nx * real_data_g*inv_Ny*inv_Nx
+ ! Second term -df/dy x dg/dx
+ DO ikx = ikxs, ikxe
+ DO iky = ikys, ikye
+ cmpx_data_f(iky,ikx-local_nkx_offset) = &
+ imagu*kyarray(iky)*F_(ikx,iky)*AA_x(ikx)*AA_y(iky) !Anti aliasing filter
+ cmpx_data_g(iky,ikx-local_nkx_offset) = &
+ imagu*kxarray(ikx)*G_(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)
+ bracket_sum_r = bracket_sum_r - real_data_f*inv_Ny*inv_Nx * real_data_g*inv_Ny*inv_Nx
+END SUBROUTINE poisson_bracket_and_sum
- real_data_c = real_data_f/NZ_/NR_ * real_data_g/NZ_/NR_
-
- 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
- C_2D(ikx,iky) = cmpx_data_c(iky,ikx-local_nkx_offset)*AA_x(ikx)*AA_y(iky)
- end do
- end do
+!!! Compute the poisson bracket of [F,G] to real space
+! - Compute the convolution using the convolution theorem
+SUBROUTINE convolve_and_add( F_, G_)
+ IMPLICIT NONE
+ COMPLEX(C_DOUBLE_COMPLEX), DIMENSION(ikxs:ikxe,ikys:ikye),&
+ INTENT(IN) :: F_, G_ ! input fields
+ ! First term df/dx x dg/dy
+ DO ikx = ikxs, ikxe
+ DO iky = ikys, ikye
+ cmpx_data_f(iky,ikx-local_nkx_offset) = F_(ikx,iky)*AA_x(ikx)*AA_y(iky) !Anti aliasing filter
+ cmpx_data_g(iky,ikx-local_nkx_offset) = G_(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)
+ bracket_sum_r = bracket_sum_r + real_data_f*inv_Ny*inv_Nx * real_data_g*inv_Ny*inv_Nx
+END SUBROUTINE convolve_and_add
-END SUBROUTINE convolve_2D_F2F
SUBROUTINE finalize_plans
IMPLICIT NONE
diff --git a/src/grid_mod.F90 b/src/grid_mod.F90
index e382af16b9ff19d04662d65e54b2a94dca9a3319..5323138fb0a971546dd1c31467e9ce96287dcbb7 100644
--- a/src/grid_mod.F90
+++ b/src/grid_mod.F90
@@ -233,7 +233,7 @@ CONTAINS
SUBROUTINE set_kxgrid
USE prec_const
- USE model, ONLY: NON_LIN
+ USE model, ONLY: LINEARITY
IMPLICIT NONE
INTEGER :: i_
@@ -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. (NON_LIN .EQ. 0)) THEN
+ IF ( (kxarray(ikx) .LT. two_third_kxmax) .OR. (LINEARITY .EQ. 'linear')) THEN
AA_x(ikx) = 1._dp;
ELSE
AA_x(ikx) = 0._dp;
@@ -295,7 +295,7 @@ CONTAINS
SUBROUTINE set_kygrid
USE prec_const
- USE model, ONLY: NON_LIN
+ USE model, ONLY: LINEARITY
IMPLICIT NONE
INTEGER :: i_, counter
@@ -329,11 +329,11 @@ CONTAINS
contains_ky0 = .true.
ENDIF
! Finding local kymax
- IF (ABS(kyarray(ikx)) .GT. local_kymax) THEN
+ IF (ABS(kyarray(iky)) .GT. local_kymax) THEN
local_kymax = ABS(kyarray(iky))
ENDIF
! Finding kymax
- IF (kyarray(ikx) .EQ. ky_max) ikx_max = ikx
+ IF (kyarray(iky) .EQ. ky_max) ikx_max = ikx
END DO
! Build the full grids on process 0 to diagnose it without comm
! ky
@@ -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. (NON_LIN .EQ. 0)) THEN
+ (kyarray(iky) .LT. two_third_kymax)) .OR. (LINEARITY .EQ. 'linear')) THEN
AA_y(iky) = 1._dp;
ELSE
AA_y(iky) = 0._dp;
diff --git a/src/inital.F90 b/src/inital.F90
index 04e4896df26d4158876f41457b9003c92ad62edb..880143c8187543a9737995af62c3cb28a31f994e 100644
--- a/src/inital.F90
+++ b/src/inital.F90
@@ -4,7 +4,6 @@
SUBROUTINE inital
USE basic
- USE model, ONLY : CO
USE initial_par
USE prec_const
USE time_integration
@@ -14,7 +13,7 @@ SUBROUTINE inital
USE closure
USE ghosts
USE restarts
- USE numerics, ONLY: wipe_turbulence, wipe_zonalflow
+ USE numerics, ONLY: wipe_turbulence, play_with_modes, save_EM_ZF_modes
USE processing, ONLY: compute_nadiab_moments
IMPLICIT NONE
@@ -41,17 +40,11 @@ SUBROUTINE inital
ENDIF
ENDIF
- ! Option for wiping the ZF modes (ky==0)
- IF ( WIPE_ZF .GT. 0 ) THEN
- IF (my_id .EQ. 0) WRITE(*,*) '-Wiping ZF modes'
- CALL wipe_zonalflow
- ENDIF
+ ! Save zonal and entropy modes
+ CALL save_EM_ZF_modes
- ! Option for wiping the turbulence and check growth of secondary inst.
- IF ( WIPE_TURB .GT. 0 ) THEN
- IF (my_id .EQ. 0) WRITE(*,*) '-Wiping turbulence'
- CALL wipe_turbulence
- ENDIF
+ ! Freeze/Wipe some selected modes (entropy,zonal,turbulent)
+ CALL play_with_modes
! Option for initializing a gaussian blob on the zonal profile
IF ( INIT_BLOB ) THEN
@@ -68,21 +61,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 (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
- CALL load_COSOlver_mat
- ! Compute collision
- CALL compute_TColl ! compute C_0 = C(N_0)
- ENDIF
+ IF(cosolver_coll) CALL load_COSOlver_mat
+ ! Compute collision
+ CALL compute_TColl ! compute C_0 = C(N_0)
END SUBROUTINE inital
!******************************************************************************!
@@ -97,7 +83,7 @@ SUBROUTINE init_moments
USE prec_const
USE utility, ONLY: checkfield
USE initial_par
- USE model, ONLY : NON_LIN
+ USE model, ONLY : LINEARITY
IMPLICIT NONE
REAL(dp) :: noise
@@ -152,7 +138,7 @@ SUBROUTINE init_moments
END DO
! Putting to zero modes that are not in the 2/3 Orszag rule
- IF (NON_LIN .GT. 0) THEN
+ IF (LINEARITY .NE. 'linear') THEN
DO ikx=ikxs,ikxe
DO iky=ikys,ikye
DO iz=izs,ize
@@ -183,7 +169,7 @@ SUBROUTINE init_gyrodens
USE prec_const
USE utility, ONLY: checkfield
USE initial_par
- USE model, ONLY : NON_LIN
+ USE model, ONLY : LINEARITY
IMPLICIT NONE
REAL(dp) :: noise
@@ -246,7 +232,7 @@ SUBROUTINE init_gyrodens
END DO
! Putting to zero modes that are not in the 2/3 Orszag rule
- IF (NON_LIN .GT. 0) THEN
+ IF (LINEARITY .NE. 'linear') THEN
DO ikx=ikxs,ikxe
DO iky=ikys,ikye
DO iz=izs,ize
diff --git a/src/initial_par_mod.F90 b/src/initial_par_mod.F90
index ffb71a686177be43b77fa35d4191f5a73d6f0800..fc72d99e9a7afea02f9955391653819a100bfc8c 100644
--- a/src/initial_par_mod.F90
+++ b/src/initial_par_mod.F90
@@ -10,8 +10,8 @@ MODULE initial_par
! Initialization through a zonal flow phi
INTEGER, PUBLIC, PROTECTED :: INIT_ZF = 0
REAL(DP), PUBLIC, PROTECTED :: ZF_AMP = 1E+3_dp
- ! Wipe zonal flow in the restart (=1) or at each step (=2), maintain (=-1)
- INTEGER, PUBLIC, PROTECTED :: WIPE_ZF = 0
+ ! Act on modes artificially (keep/wipe, zonal, non zonal, entropy mode etc.)
+ CHARACTER(len=32), PUBLIC, PROTECTED :: ACT_ON_MODES = 'nothing'
! Wipe turbulence in the restart (=1) or at each step (=2)
INTEGER, PUBLIC, PROTECTED :: WIPE_TURB = 0
! Init a Gaussian blob density in the middle
@@ -23,8 +23,6 @@ MODULE initial_par
! Initialization for random number generator
INTEGER, PUBLIC, PROTECTED :: iseed=42
- CHARACTER(len=128), PUBLIC :: mat_file ! COSOlver matrix file names
-
PUBLIC :: initial_outputinputs, initial_readinputs
CONTAINS
@@ -39,13 +37,12 @@ CONTAINS
NAMELIST /INITIAL_CON/ INIT_NOISY_PHI
NAMELIST /INITIAL_CON/ INIT_ZF
- NAMELIST /INITIAL_CON/ WIPE_ZF
+ NAMELIST /INITIAL_CON/ ACT_ON_MODES
NAMELIST /INITIAL_CON/ WIPE_TURB
NAMELIST /INITIAL_CON/ INIT_BLOB
NAMELIST /INITIAL_CON/ init_background
NAMELIST /INITIAL_CON/ init_noiselvl
NAMELIST /INITIAL_CON/ iseed
- NAMELIST /INITIAL_CON/ mat_file
READ(lu_in,initial_con)
!WRITE(*,initial_con)
diff --git a/src/memory.F90 b/src/memory.F90
index f2c9c6c7bed8c9f0d3c063466165aa138daf69df..b1508b08e8ecaacb3945217bb620998f5f42315c 100644
--- a/src/memory.F90
+++ b/src/memory.F90
@@ -6,7 +6,8 @@ SUBROUTINE memory
USE fields
USE grid
USE time_integration
- USE model, ONLY: CO, NON_LIN, KIN_E
+ USE model, ONLY: LINEARITY, KIN_E
+ USE collision
USE prec_const
IMPLICIT NONE
@@ -14,6 +15,7 @@ SUBROUTINE memory
! Electrostatic potential
CALL allocate_array( phi, ikxs,ikxe, ikys,ikye, izs,ize)
CALL allocate_array( phi_ZF, ikxs,ikxe, izs,ize)
+ CALL allocate_array( phi_EM, ikys,ikye, izs,ize)
CALL allocate_array(inv_poisson_op, ikxs,ikxe, ikys,ikye, izs,ize)
!Electrons arrays
@@ -26,6 +28,7 @@ SUBROUTINE memory
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( moments_e_EM, ipsg_e,ipeg_e, ijsg_e,ijeg_e, ikys,ikye, 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)
@@ -53,6 +56,7 @@ SUBROUTINE memory
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( moments_i_EM, ipsg_i,ipeg_i, ijsg_i,ijeg_i, ikys,ikye, 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)
@@ -95,16 +99,7 @@ SUBROUTINE memory
!___________________ 2x5D ARRAYS __________________________
!! Collision matrices
- IF (CO .LT. -1) THEN !DK collision matrix (same for every k)
- IF (KIN_E) THEN
- CALL allocate_array( Ceepj, 1,(pmaxe+1)*(jmaxe+1), 1,(pmaxe+1)*(jmaxe+1), 1,1, 1,1)
- CALL allocate_array( CeipjT, 1,(pmaxe+1)*(jmaxe+1), 1,(pmaxe+1)*(jmaxe+1), 1,1, 1,1)
- CALL allocate_array( CeipjF, 1,(pmaxe+1)*(jmaxe+1), 1,(pmaxi+1)*(jmaxi+1), 1,1, 1,1)
- CALL allocate_array( CiepjT, 1,(pmaxi+1)*(jmaxi+1), 1,(pmaxi+1)*(jmaxi+1), 1,1, 1,1)
- CALL allocate_array( CiepjF, 1,(pmaxi+1)*(jmaxi+1), 1,(pmaxe+1)*(jmaxe+1), 1,1, 1,1)
- ENDIF
- CALL allocate_array( Ciipj, 1,(pmaxi+1)*(jmaxi+1), 1,(pmaxi+1)*(jmaxi+1), 1,1, 1,1)
- ELSEIF (CO .GT. 1) THEN !GK collision matrices (one for each kperp)
+ IF (gyrokin_CO) THEN !GK collision matrices (one for each kperp)
IF (KIN_E) THEN
CALL allocate_array( Ceepj, 1,(pmaxe+1)*(jmaxe+1), 1,(pmaxe+1)*(jmaxe+1), ikxs,ikxe, ikys,ikye)
CALL allocate_array( CeipjT, 1,(pmaxe+1)*(jmaxe+1), 1,(pmaxe+1)*(jmaxe+1), ikxs,ikxe, ikys,ikye)
@@ -113,6 +108,14 @@ SUBROUTINE memory
CALL allocate_array( CiepjF, 1,(pmaxi+1)*(jmaxi+1), 1,(pmaxe+1)*(jmaxe+1), ikxs,ikxe, ikys,ikye)
ENDIF
CALL allocate_array( Ciipj, 1,(pmaxi+1)*(jmaxi+1), 1,(pmaxi+1)*(jmaxi+1), ikxs,ikxe, ikys,ikye)
- ENDIF
-
+ ELSE !DK collision matrix (same for every k)
+ IF (KIN_E) THEN
+ CALL allocate_array( Ceepj, 1,(pmaxe+1)*(jmaxe+1), 1,(pmaxe+1)*(jmaxe+1), 1,1, 1,1)
+ CALL allocate_array( CeipjT, 1,(pmaxe+1)*(jmaxe+1), 1,(pmaxe+1)*(jmaxe+1), 1,1, 1,1)
+ CALL allocate_array( CeipjF, 1,(pmaxe+1)*(jmaxe+1), 1,(pmaxi+1)*(jmaxi+1), 1,1, 1,1)
+ CALL allocate_array( CiepjT, 1,(pmaxi+1)*(jmaxi+1), 1,(pmaxi+1)*(jmaxi+1), 1,1, 1,1)
+ CALL allocate_array( CiepjF, 1,(pmaxi+1)*(jmaxi+1), 1,(pmaxe+1)*(jmaxe+1), 1,1, 1,1)
+ ENDIF
+ CALL allocate_array( Ciipj, 1,(pmaxi+1)*(jmaxi+1), 1,(pmaxi+1)*(jmaxi+1), 1,1, 1,1)
+ ENDIF
END SUBROUTINE memory
diff --git a/src/model_mod.F90 b/src/model_mod.F90
index c6e3d37dfa55b3a3ac52d47e719ab35b2036b13a..2f56451f96761098c4a91da9c35d8558e963c160 100644
--- a/src/model_mod.F90
+++ b/src/model_mod.F90
@@ -5,11 +5,11 @@ MODULE model
IMPLICIT NONE
PRIVATE
- INTEGER, PUBLIC, PROTECTED :: CO = 0 ! Collision Operator
INTEGER, PUBLIC, PROTECTED :: CLOS = 0 ! linear truncation method
INTEGER, PUBLIC, PROTECTED :: NL_CLOS = 0 ! nonlinear truncation method
INTEGER, PUBLIC, PROTECTED :: KERN = 0 ! Kernel model
- INTEGER, PUBLIC, PROTECTED :: NON_LIN = 0 ! To turn on non linear bracket term
+ CHARACTER(len=16), &
+ PUBLIC, PROTECTED ::LINEARITY= 'linear' ! 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)
@@ -54,7 +54,7 @@ CONTAINS
USE prec_const
IMPLICIT NONE
- NAMELIST /MODEL_PAR/ CO, CLOS, NL_CLOS, KERN, NON_LIN, KIN_E, mu, mu_p, mu_j, nu, tau_e, tau_i, sigma_e, sigma_i, &
+ NAMELIST /MODEL_PAR/ CLOS, NL_CLOS, KERN, LINEARITY, KIN_E, mu, mu_p, mu_j, nu, tau_e, tau_i, sigma_e, sigma_i, &
q_e, q_i, K_n, K_T, K_E, GradB, CurvB, lambdaD
READ(lu_in,model_par)
@@ -99,10 +99,9 @@ CONTAINS
INTEGER, INTENT(in) :: fidres
CHARACTER(len=256), INTENT(in) :: str
- CALL attach(fidres, TRIM(str), "CO", CO)
CALL attach(fidres, TRIM(str), "CLOS", CLOS)
CALL attach(fidres, TRIM(str), "KERN", KERN)
- CALL attach(fidres, TRIM(str), "NON_LIN", NON_LIN)
+ CALL attach(fidres, TRIM(str), "LINEARITY", LINEARITY)
CALL attach(fidres, TRIM(str), "KIN_E", KIN_E)
CALL attach(fidres, TRIM(str), "nu", nu)
CALL attach(fidres, TRIM(str), "mu", mu)
diff --git a/src/moments_eq_rhs.F90 b/src/moments_eq_rhs.F90
index eed44047f87e731a06c29b645d0f115b6cb74020..29f34e8035c0403523c186be8d0f737964bc070f 100644
--- a/src/moments_eq_rhs.F90
+++ b/src/moments_eq_rhs.F90
@@ -23,7 +23,6 @@ SUBROUTINE moments_eq_rhs_e
COMPLEX(dp) :: Tnepj, Tnepp2j, Tnepm2j, Tnepjp1, Tnepjm1, Tpare, Tphi ! Terms from b x gradB and drives
COMPLEX(dp) :: Tmir, Tnepp1j, Tnepm1j, Tnepp1jm1, Tnepm1jm1 ! Terms from mirror force with non adiab moments
COMPLEX(dp) :: UNepm1j, UNepm1jp1, UNepm1jm1 ! Terms from mirror force with adiab moments
- COMPLEX(dp) :: TColl ! terms of the rhs
COMPLEX(dp) :: i_ky
INTEGER :: izm2, izm1, izp1, izp2 ! indices for centered FDF ddz
! To store derivatives and odd-even z grid interpolations
@@ -101,13 +100,13 @@ SUBROUTINE moments_eq_rhs_e
ENDIF
!! Collision
- IF (CO .EQ. 0) THEN ! Lenard Bernstein
- CALL LenardBernstein_e(ip,ij,ikx,iky,iz,TColl)
- ELSEIF (CO .EQ. 1) THEN ! GK Dougherty
- CALL DoughertyGK_e(ip,ij,ikx,iky,iz,TColl)
- ELSE ! COSOLver matrix
- TColl = TColl_e(ip,ij,ikx,iky,iz)
- ENDIF
+ ! IF (CO .EQ. 0) THEN ! Lenard Bernstein
+ ! CALL LenardBernstein_e(ip,ij,ikx,iky,iz,TColl)
+ ! ELSEIF (CO .EQ. 1) THEN ! GK Dougherty
+ ! CALL DoughertyGK_e(ip,ij,ikx,iky,iz,TColl)
+ ! ELSE ! COSOLver matrix
+ ! TColl = TColl_e(ip,ij,ikx,iky,iz)
+ ! ENDIF
!! Sum of all linear terms (the sign is inverted to match RHS)
moments_rhs_e(ip,ij,ikx,iky,iz,updatetlevel) = &
@@ -124,7 +123,7 @@ SUBROUTINE moments_eq_rhs_e
! Numerical hyperdiffusion (totally artificial, for stability purpose)
- mu*kperp2**2 * moments_e(ip,ij,ikx,iky,iz,updatetlevel) &
! Collision term
- + TColl
+ + TColl_e(ip,ij,ikx,iky,iz)
!! Adding non linearity
moments_rhs_e(ip,ij,ikx,iky,iz,updatetlevel) = &
@@ -165,7 +164,6 @@ SUBROUTINE moments_eq_rhs_i
COMPLEX(dp) :: Tnipj, Tnipp2j, Tnipm2j, Tnipjp1, Tnipjm1, Tpari, Tphi
COMPLEX(dp) :: Tmir, Tnipp1j, Tnipm1j, Tnipp1jm1, Tnipm1jm1 ! Terms from mirror force with non adiab moments
COMPLEX(dp) :: UNipm1j, UNipm1jp1, UNipm1jm1 ! Terms from mirror force with adiab moments
- COMPLEX(dp) :: TColl ! terms of the rhs
COMPLEX(dp) :: i_ky
INTEGER :: izm2, izm1, izp1, izp2 ! indices for centered FDF ddz
! To store derivatives and odd-even z grid interpolations
@@ -239,13 +237,13 @@ SUBROUTINE moments_eq_rhs_i
ENDIF
!! Collision
- IF (CO .EQ. 0) THEN ! Lenard Bernstein
- CALL LenardBernstein_i(ip,ij,ikx,iky,iz,TColl)
- ELSEIF (CO .EQ. 1) THEN ! GK Dougherty
- CALL DoughertyGK_i(ip,ij,ikx,iky,iz,TColl)
- ELSE! COSOLver matrix (Sugama, Coulomb)
- TColl = TColl_i(ip,ij,ikx,iky,iz)
- ENDIF
+ ! IF (CO .EQ. 0) THEN ! Lenard Bernstein
+ ! CALL LenardBernstein_i(ip,ij,ikx,iky,iz,TColl)
+ ! ELSEIF (CO .EQ. 1) THEN ! GK Dougherty
+ ! CALL DoughertyGK_i(ip,ij,ikx,iky,iz,TColl)
+ ! ELSE! COSOLver matrix (Sugama, Coulomb)
+ ! TColl = TColl_i(ip,ij,ikx,iky,iz)
+ ! ENDIF
!! Sum of all linear terms (the sign is inverted to match RHS)
moments_rhs_i(ip,ij,ikx,iky,iz,updatetlevel) = &
@@ -262,7 +260,7 @@ SUBROUTINE moments_eq_rhs_i
! Numerical hyperdiffusion (totally artificial, for stability purpose)
- mu*kperp2**2 * moments_i(ip,ij,ikx,iky,iz,updatetlevel) &
! Collision term
- + TColl
+ + TColl_i(ip,ij,ikx,iky,iz)
!! Adding non linearity
moments_rhs_i(ip,ij,ikx,iky,iz,updatetlevel) = &
diff --git a/src/numerics_mod.F90 b/src/numerics_mod.F90
index da59158ef16b5bd1d948b78ba0c7676fa11a1af4..1b3926173c113de07eaf802f64ae1eb35d8c3dbd 100644
--- a/src/numerics_mod.F90
+++ b/src/numerics_mod.F90
@@ -7,7 +7,7 @@ MODULE numerics
implicit none
PUBLIC :: build_dnjs_table, evaluate_kernels, evaluate_poisson_op, compute_lin_coeff
- PUBLIC :: wipe_turbulence, wipe_zonalflow
+ PUBLIC :: wipe_turbulence, play_with_modes, save_EM_ZF_modes
CONTAINS
@@ -282,41 +282,83 @@ SUBROUTINE wipe_turbulence
ENDDO
END SUBROUTINE
!******************************************************************************!
-!!!!!!! Remove all ky==0 modes to conserve only non zonal modes in a restart
+!!!!!!! Routine that can artificially increase or wipe modes
!******************************************************************************!
-SUBROUTINE wipe_zonalflow
+SUBROUTINE save_EM_ZF_modes
USE fields
+ USE array, ONLY : moments_e_ZF, moments_i_ZF, phi_ZF, moments_e_EM,moments_i_EM,phi_EM
USE grid
+ USE time_integration, ONLY: updatetlevel
+ USE initial_par, ONLY: ACT_ON_MODES
IMPLICIT NONE
- DO ikx=ikxs,ikxe
- DO iky=ikys,ikye
- DO iz=izs,ize
- DO ip=ips_e,ipe_e
- DO ij=ijs_e,ije_e
- IF( iky .EQ. iky_0) THEN
- moments_e( ip,ij,ikx,iky,iz, :) = 0e-3_dp*moments_e( ip,ij,ikx,iky,iz, :)
- ELSE
- moments_e( ip,ij,ikx,iky,iz, :) = 1e+0_dp*moments_e( ip,ij,ikx,iky,iz, :)
- ENDIF
- ENDDO
- ENDDO
- DO ip=ips_i,ipe_i
- DO ij=ijs_i,ije_i
- IF( iky .EQ. iky_0) THEN
- moments_i( ip,ij,ikx,iky,iz, :) = 0e-3_dp*moments_i( ip,ij,ikx,iky,iz, :)
- ELSE
- moments_i( ip,ij,ikx,iky,iz, :) = 1e+0_dp*moments_i( ip,ij,ikx,iky,iz, :)
- ENDIF
+ ! Store Zonal and entropy modes
+ 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)
+ IF(contains_kx0) THEN
+ moments_e_EM(ips_e:ipe_e,ijs_e:ije_e,ikys:ikye,izs:ize) = moments_e(ips_e:ipe_e,ijs_e:ije_e,ikx_0,ikys:ikye,izs:ize,updatetlevel)
+ moments_i_EM(ips_e:ipe_e,ijs_e:ije_e,ikys:ikye,izs:ize) = moments_i(ips_i:ipe_i,ijs_i:ije_i,ikx_0,ikys:ikye,izs:ize,updatetlevel)
+ phi_EM(ikys:ikye,izs:ize) = phi(ikx_0,ikys:ikye,izs:ize)
+ ELSE
+ moments_e_EM(ips_e:ipe_e,ijs_e:ije_e,ikys:ikye,izs:ize) = 0._dp
+ moments_i_EM(ips_e:ipe_e,ijs_e:ije_e,ikys:ikye,izs:ize) = 0._dp
+ phi_EM(ikys:ikye,izs:ize) = 0._dp
+ ENDIF
+END SUBROUTINE
+
+SUBROUTINE play_with_modes
+ USE fields
+ USE array, ONLY : moments_e_ZF, moments_i_ZF, phi_ZF, moments_e_EM,moments_i_EM,phi_EM
+ USE grid
+ USE time_integration, ONLY: updatetlevel
+ USE initial_par, ONLY: ACT_ON_MODES
+ IMPLICIT NONE
+ REAL(dp) :: AMP = 1.5_dp
+
+ SELECT CASE(ACT_ON_MODES)
+ CASE('wipe_zonal') ! Errase the zonal flow
+ moments_e(ips_e:ipe_e,ijs_e:ije_e,ikxs:ikxe,iky_0,izs:ize,updatetlevel) = 0._dp
+ moments_i(ips_i:ipe_i,ijs_i:ije_i,ikxs:ikxe,iky_0,izs:ize,updatetlevel) = 0._dp
+ phi(ikxs:ikxe,iky_0,izs:ize) = 0._dp
+ CASE('wipe_entropymode')
+ moments_e(ips_e:ipe_e,ijs_e:ije_e,ikx_0,ikys:ikye,izs:ize,updatetlevel) = 0._dp
+ moments_i(ips_i:ipe_i,ijs_i:ije_i,ikx_0,ikys:ikye,izs:ize,updatetlevel) = 0._dp
+ phi(ikx_0,ikys:ikye,izs:ize) = 0._dp
+ CASE('wipe_turbulence')
+ DO ikx = ikxs,ikxe
+ DO iky = ikys, ikye
+ IF ( (ikx .NE. ikx_0) .AND. (iky .NE. iky_0) ) THEN
+ moments_e(ips_e:ipe_e,ijs_e:ije_e,ikx,iky,izs:ize,updatetlevel) = 0._dp
+ moments_i(ips_i:ipe_i,ijs_i:ije_i,ikx,iky,izs:ize,updatetlevel) = 0._dp
+ phi(ikx,iky,izs:ize) = 0._dp
+ ENDIF
+ ENDDO
ENDDO
+ CASE('wipe_nonzonal')
+ DO ikx = ikxs,ikxe
+ DO iky = ikys, ikye
+ IF ( (ikx .NE. ikx_0) ) THEN
+ moments_e(ips_e:ipe_e,ijs_e:ije_e,ikx,iky,izs:ize,updatetlevel) = 0._dp
+ moments_i(ips_i:ipe_i,ijs_i:ije_i,ikx,iky,izs:ize,updatetlevel) = 0._dp
+ phi(ikx,iky,izs:ize) = 0._dp
+ ENDIF
+ ENDDO
ENDDO
- IF( iky .EQ. iky_0) THEN
- phi(ikx,iky,iz) = 0e-3_dp*phi(ikx,iky,iz)
- ELSE
- phi(ikx,iky,iz) = 1e+0_dp*phi(ikx,iky,iz)
+ CASE('freeze_zonal')
+ moments_e(ips_e:ipe_e,ijs_e:ije_e,ikxs:ikxe,iky_0,izs:ize,updatetlevel) = moments_e_ZF(ips_e:ipe_e,ijs_e:ije_e,ikxs:ikxe,izs:ize)
+ moments_i(ips_i:ipe_i,ijs_i:ije_i,ikxs:ikxe,iky_0,izs:ize,updatetlevel) = moments_i_ZF(ips_i:ipe_i,ijs_i:ije_i,ikxs:ikxe,izs:ize)
+ phi(ikxs:ikxe,iky_0,izs:ize) = phi_ZF(ikxs:ikxe,izs:ize)
+ CASE('freeze_entropymode')
+ IF(contains_kx0) THEN
+ moments_e(ips_e:ipe_e,ijs_e:ije_e,ikx_0,ikys:ikye,izs:ize,updatetlevel) = moments_e_EM(ips_e:ipe_e,ijs_e:ije_e,ikys:ikye,izs:ize)
+ moments_i(ips_i:ipe_i,ijs_i:ije_i,ikx_0,ikys:ikye,izs:ize,updatetlevel) = moments_i_EM(ips_i:ipe_i,ijs_i:ije_i,ikys:ikye,izs:ize)
+ phi(ikx_0,ikys:ikye,izs:ize) = phi_EM(ikys:ikye,izs:ize)
ENDIF
- ENDDO
- ENDDO
- ENDDO
+ CASE('amplify_zonal')
+ moments_e(ips_e:ipe_e,ijs_e:ije_e,ikxs:ikxe,iky_0,izs:ize,updatetlevel) = AMP*moments_e_ZF(ips_e:ipe_e,ijs_e:ije_e,ikxs:ikxe,izs:ize)
+ moments_i(ips_i:ipe_i,ijs_i:ije_i,ikxs:ikxe,iky_0,izs:ize,updatetlevel) = AMP*moments_i_ZF(ips_i:ipe_i,ijs_i:ije_i,ikxs:ikxe,izs:ize)
+ phi(ikxs:ikxe,iky_0,izs:ize) = AMP*phi_ZF(ikxs:ikxe,izs:ize)
+ END SELECT
END SUBROUTINE
END MODULE numerics
diff --git a/src/readinputs.F90 b/src/readinputs.F90
index e8c6a499c78b158c50783fe91ea9c0116124c077..1322ddeb871a2d633adfaeaa640ae5a7c9610d5f 100644
--- a/src/readinputs.F90
+++ b/src/readinputs.F90
@@ -3,6 +3,7 @@ SUBROUTINE readinputs
USE grid, ONLY: grid_readinputs
USE diagnostics_par, ONLY: diag_par_readinputs
+ USE collision, ONLY: collision_readinputs
USE model, ONLY: model_readinputs
USE initial_par, ONLY: initial_readinputs
USE time_integration, ONLY: time_integration_readinputs
@@ -23,6 +24,9 @@ SUBROUTINE readinputs
! Load model parameters from input file
CALL model_readinputs
+ ! Load collision parameters from input file
+ CALL collision_readinputs
+
! Load initial condition parameters from input file
CALL initial_readinputs
diff --git a/src/srcinfo.h b/src/srcinfo.h
index 5b3f294989538a1cc8233fc1241b30c097cf8000..7cc5e9ad5cd3212593dd6b361c4cf1a77640ab8f 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='0d685e3-dirty')
+parameter (VERSION='2df4eb7-dirty')
parameter (BRANCH='master')
parameter (AUTHOR='ahoffman')
-parameter (EXECDATE='Thu Nov 11 10:19:20 CET 2021')
+parameter (EXECDATE='Tue Nov 23 10:32:57 CET 2021')
parameter (HOST ='spcpc606')
diff --git a/src/srcinfo/srcinfo.h b/src/srcinfo/srcinfo.h
index 5b3f294989538a1cc8233fc1241b30c097cf8000..7cc5e9ad5cd3212593dd6b361c4cf1a77640ab8f 100644
--- a/src/srcinfo/srcinfo.h
+++ b/src/srcinfo/srcinfo.h
@@ -3,8 +3,8 @@ character(len=40) BRANCH
character(len=20) AUTHOR
character(len=40) EXECDATE
character(len=40) HOST
-parameter (VERSION='0d685e3-dirty')
+parameter (VERSION='2df4eb7-dirty')
parameter (BRANCH='master')
parameter (AUTHOR='ahoffman')
-parameter (EXECDATE='Thu Nov 11 10:19:20 CET 2021')
+parameter (EXECDATE='Tue Nov 23 10:32:57 CET 2021')
parameter (HOST ='spcpc606')
diff --git a/src/stepon.F90 b/src/stepon.F90
index a59d15fabbe0b0aa0d752fe4e3dcfa840767a258..a57b3e4bde6d2e6391f583d9b596dd80096d7d1a 100644
--- a/src/stepon.F90
+++ b/src/stepon.F90
@@ -6,13 +6,13 @@ SUBROUTINE stepon
USE closure
USE collision, ONLY : compute_TColl
USE fields, ONLY: moments_e, moments_i, phi
- USE initial_par, ONLY: WIPE_ZF, WIPE_TURB
+ USE initial_par, ONLY: ACT_ON_MODES, WIPE_TURB
USE ghosts
USE grid
- USE model, ONLY : NON_LIN, KIN_E
+ USE model, ONLY : LINEARITY, KIN_E
use prec_const
USE time_integration
- USE numerics, ONLY: wipe_zonalflow, wipe_turbulence
+ USE numerics, ONLY: play_with_modes, wipe_turbulence
USE processing, ONLY: compute_nadiab_moments
USE utility, ONLY: checkfield
@@ -45,8 +45,8 @@ SUBROUTINE stepon
CALL compute_nadiab_moments
! Update nonlinear term S_n -> S_n+1(phi_n+1,N_n+1)
CALL compute_Sapj
- ! Cancel zonal modes artificially
- IF ( WIPE_ZF .EQ. 2) CALL wipe_zonalflow
+ ! Store or cancel/maintain zonal modes artificially
+ CALL play_with_modes
! Cancel non zonal modes artificially
IF ( WIPE_TURB .EQ. 2) CALL wipe_turbulence
!- Check before next step
@@ -63,8 +63,8 @@ SUBROUTINE stepon
! Execution time start
CALL cpu_time(t0_checkfield)
- 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
+ IF(LINEARITY .NE. 'linear') CALL anti_aliasing ! ensure 0 mode for 2/3 rule
+ IF(LINEARITY .NE. 'linear') CALL enforce_symmetry ! Enforcing symmetry on kx = 0
mlend=.FALSE.
IF(.NOT.nlend) THEN
diff --git a/src/time_integration_mod.F90 b/src/time_integration_mod.F90
index 287e304e23fdcc6701d29ee3458bf6c5563d3bc3..7203bdf0538ec784de84dbf4cc3ad74f2f32d93f 100644
--- a/src/time_integration_mod.F90
+++ b/src/time_integration_mod.F90
@@ -32,7 +32,6 @@ CONTAINS
NAMELIST /TIME_INTEGRATION_PAR/ numerical_scheme
READ(lu_in,time_integration_par)
- !WRITE(*,time_integration_par)
CALL set_numerical_scheme