diff --git a/src/fourier_mod.F90 b/src/fourier_mod.F90
index 53dfde72f2c86ef4807c3320e743f8101ae4a6d4..6a51df0abb0990ed5487ae6b4443d50af526ab56 100644
--- a/src/fourier_mod.F90
+++ b/src/fourier_mod.F90
@@ -92,7 +92,6 @@ MODULE fourier
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
bracket_sum_r = bracket_sum_r + real_data_f * real_data_g*inv_Ny*inv_Nx
! Second term -df/dy x dg/dx
DO ikx = ikxs, ikxe
@@ -105,7 +104,6 @@ MODULE fourier
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
bracket_sum_r = bracket_sum_r - real_data_f * real_data_g*inv_Ny*inv_Nx
END SUBROUTINE poisson_bracket_and_sum
@@ -124,7 +122,6 @@ SUBROUTINE convolve_and_add( F_, G_)
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
bracket_sum_r = bracket_sum_r + real_data_f * real_data_g*inv_Ny*inv_Nx
END SUBROUTINE convolve_and_add
diff --git a/src/nonlinear_mod.F90 b/src/nonlinear_mod.F90
index 90fe637b9979ed6341bc9c697d4c6cad8820cd45..8e9a2ce5cb3d0374acfbf7846a192cc37f56011b 100644
--- a/src/nonlinear_mod.F90
+++ b/src/nonlinear_mod.F90
@@ -4,7 +4,7 @@ MODULE nonlinear
USE initial_par, ONLY : ACT_ON_MODES
USE basic
USE fourier
- USE fields, ONLY : phi, moments_e, moments_i
+ USE fields, ONLY : phi, psi, moments_e, moments_i
USE grid
USE model
USE prec_const
@@ -18,7 +18,7 @@ MODULE nonlinear
INTEGER :: in, is, p_int, j_int
INTEGER :: nmax, smax ! Upper bound of the sums
- REAL(dp):: kx, ky, kerneln
+ REAL(dp):: kx, ky, kerneln, sqrt_p, sqrt_pp1
PUBLIC :: compute_Sapj, nonlinear_init
CONTAINS
@@ -69,92 +69,129 @@ SUBROUTINE compute_nonlinear
!!!!!!!!!!!!!!!!!!!! ELECTRON non 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)
- p_int = parray_e(ip)
- jloope: DO ij = ijs_e, ije_e ! Loop over Laguerre moments
- j_int=jarray_e(ij)
- IF((CLOS .NE. 1) .OR. (p_int+2*j_int .LE. dmaxe)) THEN !compute
- ! Set non linear sum truncation
- IF (NL_CLOS .EQ. -2) THEN
- nmax = Jmaxe
- ELSEIF (NL_CLOS .EQ. -1) THEN
- nmax = Jmaxe-j_int
+
+ ploope: DO ip = ips_e,ipe_e ! Loop over Hermite moments
+ eo = MODULO(parray_e(ip),2)
+ p_int = parray_e(ip)
+ sqrt_p = SQRT(REAL(p_int,dp)); sqrt_pp1 = SQRT(REAL(p_int,dp)+1._dp);
+
+ jloope: DO ij = ijs_e, ije_e ! Loop over Laguerre moments
+ j_int=jarray_e(ij)
+ IF((CLOS .NE. 1) .OR. (p_int+2*j_int .LE. dmaxe)) THEN !compute
+ ! Set non linear sum truncation
+ IF (NL_CLOS .EQ. -2) THEN
+ nmax = Jmaxe
+ ELSEIF (NL_CLOS .EQ. -1) THEN
+ nmax = Jmaxe-j_int
+ ELSE
+ nmax = min(NL_CLOS,Jmaxe-j_int)
+ ENDIF
+ bracket_sum_r = 0._dp ! initialize sum over real nonlinear term
+
+ nloope: DO in = 1,nmax+1 ! Loop over laguerre for the sum
+!-----------!! ELECTROSTATIC CONTRIBUTION {Sum_s dnjs Naps, Kernel phi}
+ ! First convolution terms
+ F_cmpx(ikys:ikye,ikxs:ikxe) = phi(ikys:ikye,ikxs:ikxe,iz) * kernel_e(in, ikys:ikye,ikxs:ikxe, iz, eo)
+ ! Second convolution terms
+ G_cmpx(ikys:ikye,ikxs:ikxe) = 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(ikys:ikye,ikxs:ikxe) = G_cmpx(ikys:ikye,ikxs:ikxe) + &
+ dnjs(in,ij,is) * moments_e(ip,is,ikys:ikye,ikxs:ikxe,iz,updatetlevel)
+ ENDDO
+ !/!\ this function add its result to bracket_sum_r (hard to read sorry) /!\
+ CALL poisson_bracket_and_sum(F_cmpx,G_cmpx)
+
+!-----------!! 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
+ F_cmpx(ikys:ikye,ikxs:ikxe) = -sqrt_tau_o_sigma_e * psi(ikys:ikye,ikxs:ikxe,iz) * kernel_e(in, ikys:ikye,ikxs:ikxe, iz, eo)
+ ! Second convolution terms
+ G_cmpx(ikys:ikye,ikxs:ikxe) = 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(ikys:ikye,ikxs:ikxe) = G_cmpx(ikys:ikye,ikxs:ikxe) + &
+ dnjs(in,ij,is) * (sqrt_pp1*moments_e(ip+1,is,ikys:ikye,ikxs:ikxe,iz,updatetlevel)&
+ +sqrt_p *moments_e(ip-1,is,ikys:ikye,ikxs:ikxe,iz,updatetlevel))
+ ENDDO
+ !/!\ this function add its result to bracket_sum_r (hard to read sorry) /!\
+ CALL poisson_bracket_and_sum(F_cmpx,G_cmpx)
+ ENDIF
+ ENDDO nloope
+
+!---------! Put back 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 iky = ikys, ikye
+ Sepj(ip,ij,iky,ikxs:ikxe,iz) = bracket_sum_c(ikxs:ikxe,iky-local_nky_offset)*AA_x(ikxs:ikxe)*AA_y(iky) !Anti aliasing filter
+ ENDDO
ELSE
- nmax = min(NL_CLOS,Jmaxe-j_int)
+ Sepj(ip,ij,:,:,iz) = 0._dp
ENDIF
- bracket_sum_r = 0._dp ! initialize sum over real nonlinear term
- nloope: DO in = 1,nmax+1 ! Loop over laguerre for the sum
- ! First convolution terms
- F_cmpx(ikys:ikye,ikxs:ikxe) = phi(ikys:ikye,ikxs:ikxe,iz) * kernel_e(in, ikys:ikye,ikxs:ikxe, iz, eo)
- ! Second convolution terms
- G_cmpx(ikys:ikye,ikxs:ikxe) = 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(ikys:ikye,ikxs:ikxe) = G_cmpx(ikys:ikye,ikxs:ikxe) + &
- dnjs(in,ij,is) * moments_e(ip,is,ikys:ikye,ikxs:ikxe,iz,updatetlevel)
- ENDDO
- !/!\ 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, bracket_sum_r, bracket_sum_c)
- ! Retrieve convolution in input format
- DO iky = ikys, ikye
- Sepj(ip,ij,iky,ikxs:ikxe,iz) = bracket_sum_c(ikxs:ikxe,iky-local_nky_offset)*AA_x(ikxs:ikxe)*AA_y(iky) !Anti aliasing filter
- ENDDO
- ELSE
- Sepj(ip,ij,:,:,iz) = 0._dp
- ENDIF
- ENDDO jloope
- ENDDO ploope
-ENDDO zloope
+ 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
- p_int = parray_i(ip)
- eo = MODULO(parray_i(ip),2)
- jloopi: DO ij = ijs_i, ije_i ! Loop over Laguerre moments
- j_int=jarray_i(ij)
- IF((CLOS .NE. 1) .OR. (p_int+2*j_int .LE. dmaxi)) THEN !compute
- ! Set non linear sum truncation
- IF (NL_CLOS .EQ. -2) THEN
- nmax = Jmaxi
- ELSEIF (NL_CLOS .EQ. -1) THEN
- nmax = Jmaxi-j_int
+ zloopi: DO iz = izs,ize
+ ploopi: DO ip = ips_i,ipe_i ! Loop over Hermite moments
+ p_int = parray_i(ip)
+ eo = MODULO(parray_i(ip),2)
+ jloopi: DO ij = ijs_i, ije_i ! Loop over Laguerre moments
+ j_int=jarray_i(ij)
+ IF((CLOS .NE. 1) .OR. (p_int+2*j_int .LE. dmaxi)) THEN !compute
+ ! Set non linear sum truncation
+ IF (NL_CLOS .EQ. -2) THEN
+ nmax = Jmaxi
+ ELSEIF (NL_CLOS .EQ. -1) THEN
+ nmax = Jmaxi-j_int
+ ELSE
+ nmax = min(NL_CLOS,Jmaxi-j_int)
+ ENDIF
+ bracket_sum_r = 0._dp ! initialize sum over real nonlinear term
+ nloopi: DO in = 1,nmax+1 ! Loop over laguerre for the sum
+!-----------!! ELECTROMAGNETIC CONTRIBUTION -sqrt(tau)/sigma*{Sum_s dnjs [sqrt(p+1)Nap+1s + sqrt(p)Nap-1s], Kernel psi}
+ ! First convolution terms
+ F_cmpx(ikys:ikye,ikxs:ikxe) = phi(ikys:ikye,ikxs:ikxe,iz) * kernel_i(in, ikys:ikye,ikxs:ikxe, iz, eo)
+ ! Second convolution terms
+ G_cmpx(ikys:ikye,ikxs:ikxe) = 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(ikys:ikye,ikxs:ikxe) = G_cmpx(ikys:ikye,ikxs:ikxe) + &
+ dnjs(in,ij,is) * moments_i(ip,is,ikys:ikye,ikxs:ikxe,iz,updatetlevel)
+ ENDDO
+ !/!\ this function add its result to bracket_sum_r (hard to read sorry) /!\
+ CALL poisson_bracket_and_sum(F_cmpx,G_cmpx)
+!-----------!! 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
+ F_cmpx(ikys:ikye,ikxs:ikxe) = -sqrt_tau_o_sigma_i * psi(ikys:ikye,ikxs:ikxe,iz) * kernel_i(in, ikys:ikye,ikxs:ikxe, iz, eo)
+ ! Second convolution terms
+ G_cmpx(ikys:ikye,ikxs:ikxe) = 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(ikys:ikye,ikxs:ikxe) = G_cmpx(ikys:ikye,ikxs:ikxe) + &
+ dnjs(in,ij,is) * (sqrt_pp1*moments_i(ip+1,is,ikys:ikye,ikxs:ikxe,iz,updatetlevel)&
+ +sqrt_p *moments_i(ip-1,is,ikys:ikye,ikxs:ikxe,iz,updatetlevel))
+ ENDDO
+ !/!\ this function add its result to bracket_sum_r (hard to read sorry) /!\
+ CALL poisson_bracket_and_sum(F_cmpx,G_cmpx)
+ ENDIF
+ 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 iky = ikys, ikye
+ Sipj(ip,ij,iky,ikxs:ikxe,iz) = bracket_sum_c(ikxs:ikxe,iky-local_nky_offset)*AA_x(ikxs:ikxe)*AA_y(iky)
+ ENDDO
ELSE
- nmax = min(NL_CLOS,Jmaxi-j_int)
+ Sipj(ip,ij,:,:,iz) = 0._dp
ENDIF
- bracket_sum_r = 0._dp ! initialize sum over real nonlinear term
- nloopi: DO in = 1,nmax+1 ! Loop over laguerre for the sum
- ! First convolution terms
- F_cmpx(ikys:ikye,ikxs:ikxe) = phi(ikys:ikye,ikxs:ikxe,iz) * kernel_i(in, ikys:ikye,ikxs:ikxe, iz, eo)
- ! Second convolution terms
- G_cmpx(ikys:ikye,ikxs:ikxe) = 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(ikys:ikye,ikxs:ikxe) = G_cmpx(ikys:ikye,ikxs:ikxe) + &
- dnjs(in,ij,is) * moments_i(ip,is,ikys:ikye,ikxs:ikxe,iz,updatetlevel)
- ENDDO
- !/!\ 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, bracket_sum_r, bracket_sum_c)
- ! Retrieve convolution in input format
- DO iky = ikys, ikye
- Sipj(ip,ij,iky,ikxs:ikxe,iz) = bracket_sum_c(ikxs:ikxe,iky-local_nky_offset)*AA_x(ikxs:ikxe)*AA_y(iky)
- ENDDO
- ELSE
- Sipj(ip,ij,:,:,iz) = 0._dp
- ENDIF
- ENDDO jloopi
- ENDDO ploopi
-ENDDO zloopi
+ ENDDO jloopi
+ ENDDO ploopi
+ ENDDO zloopi
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
END SUBROUTINE compute_nonlinear
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!