Skip to content
Snippets Groups Projects
Commit 869994cb authored by Antoine Cyril David Hoffmann's avatar Antoine Cyril David Hoffmann :seedling:
Browse files

+ nonlinear electromagnetic effects

parent bc031375
No related branches found
No related tags found
No related merge requests found
......@@ -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
......
......@@ -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
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment