Skip to content
Snippets Groups Projects
Commit 388bad0a authored by Antoine Cyril David Hoffmann's avatar Antoine Cyril David Hoffmann
Browse files

simplification of kernel computation and benchmark of parallel terms

parent 0bc911b9
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
Showing
with 57 additions and 117 deletions
src/numerics_mod.F90
+ 57
− 117
View file @ 388bad0a
...@@ -54,105 +54,34 @@ END SUBROUTINE build_dnjs_table ...@@ -54,105 +54,34 @@ END SUBROUTINE build_dnjs_table
!******************************************************************************! !******************************************************************************!
SUBROUTINE evaluate_kernels SUBROUTINE evaluate_kernels
USE basic USE basic
USE array, Only : kernel_e, kernel_i, gxy, gyy, gxx USE array, Only : kernel_e, kernel_i, kparray
USE grid USE grid
USE model, ONLY : tau_e, tau_i, sigma_e, sigma_i, q_e, q_i, lambdaD, CLOS, sigmae2_taue_o2, sigmai2_taui_o2 USE model, ONLY : tau_e, tau_i, sigma_e, sigma_i, q_e, q_i, lambdaD, CLOS, sigmae2_taue_o2, sigmai2_taui_o2
IMPLICIT NONE IMPLICIT NONE
INTEGER :: j_int
REAL(dp) :: j_dp, y_, kp2_, kx_, ky_
REAL(dp) :: factj, j_dp, j_int DO ikx = ikxs,ikxe
REAL(dp) :: be_2, bi_2, alphaD DO iky = ikys,ikye
REAL(dp) :: kx, ky, kperp2 DO iz = izs,ize
!!!!! Electron kernels !!!!! !!!!! Electron kernels !!!!!
!Precompute species dependant factors DO ij = ijsg_e, ijeg_e
factj = 1.0 ! Start of the recursive factorial
DO ij = 1, jmaxe+1
j_int = jarray_e(ij) j_int = jarray_e(ij)
j_dp = REAL(j_int,dp) ! REAL of degree j_dp = REAL(j_int,dp)
y_ = sigmae2_taue_o2 * kparray(ikx,iky,iz)**2
! Recursive factorial kernel_e(ij,ikx,iky,iz) = y_**j_int*EXP(-y_)/GAMMA(j_dp+1._dp)!factj
IF (j_dp .GT. 0) THEN
factj = factj * j_dp
ELSE
factj = 1._dp
ENDIF
DO ikx = ikxs,ikxe
kx = kxarray(ikx)
DO iky = ikys,ikye
ky = kyarray(iky)
DO iz = izs,ize
kperp2= gxx(iz)*kx**2 + 2._dp*gxy(iz)*kx*ky + gyy(iz)*ky**2
be_2 = kperp2 * sigmae2_taue_o2
kernel_e(ij, ikx, iky,iz) = be_2**j_int * exp(-be_2)/factj
ENDDO
ENDDO
ENDDO
ENDDO ENDDO
! Kernels closure
DO ikx = ikxs,ikxe
kx = kxarray(ikx)
DO iky = ikys,ikye
ky = kyarray(iky)
DO iz = izs,ize
kperp2= gxx(iz)*kx**2 + 2._dp*gxy(iz)*kx*ky + gyy(iz)*ky**2
be_2 = kperp2 * sigmae2_taue_o2
! Kernel ghost + 1 with Kj+1 = y/(j+1) Kj (/!\ ij = j+1)
kernel_e(ijeg_e,ikx,iky,iz) = be_2/(real(ijeg_e-1,dp))*kernel_e(ije_e,ikx,iky,iz)
! Kernel ghost - 1 with Kj-1 = j/y Kj(careful about the kperp=0)
IF ( be_2 .NE. 0 ) THEN
kernel_e(ijsg_e,ikx,iky,iz) = (real(ijsg_e,dp))/be_2*kernel_e(ijs_e,ikx,iky,iz)
ELSE
kernel_e(ijsg_e,ikx,iky,iz) = 0._dp
ENDIF
ENDDO
ENDDO
ENDDO
!!!!! Ion kernels !!!!! !!!!! Ion kernels !!!!!
factj = 1.0 ! Start of the recursive factorial DO ij = ijsg_i, ijeg_i
DO ij = 1, jmaxi+1 j_int = jarray_i(ij)
j_int = jarray_e(ij) j_dp = REAL(j_int,dp)
j_dp = REAL(j_int,dp) ! REAL of degree y_ = sigmai2_taui_o2 * kparray(ikx,iky,iz)**2
kernel_i(ij,ikx,iky,iz) = y_**j_int*EXP(-y_)/GAMMA(j_dp+1._dp)!factj
! Recursive factorial
IF (j_dp .GT. 0) THEN
factj = factj * j_dp
ELSE
factj = 1._dp
ENDIF
DO ikx = ikxs,ikxe
kx = kxarray(ikx)
DO iky = ikys,ikye
ky = kyarray(iky)
DO iz = izs,ize
kperp2= gxx(iz)*kx**2 + 2._dp*gxy(iz)*kx*ky + gyy(iz)*ky**2
bi_2 = kperp2 * sigmai2_taui_o2
kernel_i(ij, ikx, iky,iz) = bi_2**j_int * exp(-bi_2)/factj
ENDDO
ENDDO
ENDDO
ENDDO
! Kernels closure
DO ikx = ikxs,ikxe
kx = kxarray(ikx)
DO iky = ikys,ikye
ky = kyarray(iky)
DO iz = izs,ize
kperp2= gxx(iz)*kx**2 + 2._dp*gxy(iz)*kx*ky + gyy(iz)*ky**2
bi_2 = kperp2 * sigmai2_taui_o2
! Kernel ghost + 1 with Kj+1 = y/(j+1) Kj
kernel_i(ijeg_i,ikx,iky,iz) = bi_2/(real(ijeg_i-1,dp))*kernel_i(ije_i,ikx,iky,iz)
! Kernel ghost - 1 with Kj-1 = j/y Kj(careful about the kperp=0)
IF ( bi_2 .NE. 0 ) THEN
kernel_i(ijsg_i,ikx,iky,iz) = (real(ijsg_i,dp))/bi_2*kernel_i(ijs_i,ikx,iky,iz)
ELSE
kernel_i(ijsg_i,ikx,iky,iz) = 0._dp
ENDIF
ENDDO
ENDDO
ENDDO ENDDO
ENDDO
ENDDO
ENDDO
END SUBROUTINE evaluate_kernels END SUBROUTINE evaluate_kernels
!******************************************************************************! !******************************************************************************!
...@@ -174,30 +103,35 @@ SUBROUTINE compute_lin_coeff ...@@ -174,30 +103,35 @@ SUBROUTINE compute_lin_coeff
DO ij = ijs_e, ije_e DO ij = ijs_e, ije_e
j_int= jarray_e(ij) ! Laguerre degree j_int= jarray_e(ij) ! Laguerre degree
j_dp = REAL(j_int,dp) ! REAL of Laguerre degree j_dp = REAL(j_int,dp) ! REAL of Laguerre degree
! All Napj terms
xnepj(ip,ij) = taue_qe*(CurvB*(2._dp*p_dp + 1._dp) & xnepj(ip,ij) = taue_qe*(CurvB*(2._dp*p_dp + 1._dp) &
+GradB*(2._dp*j_dp + 1._dp)) +GradB*(2._dp*j_dp + 1._dp))
ynepp1j (ip,ij) = -SQRT(tau_e)/sigma_e * (j_dp+1)*SQRT(p_dp+1._dp) ! Mirror force terms
ynepm1j (ip,ij) = -SQRT(tau_e)/sigma_e * (j_dp+1)*SQRT(p_dp) ynepp1j (ip,ij) = -SQRT(tau_e)/sigma_e * (j_dp+1)*SQRT(p_dp+1._dp)
ynepp1jm1(ip,ij) = +SQRT(tau_e)/sigma_e * j_dp*SQRT(p_dp+1._dp) ynepm1j (ip,ij) = -SQRT(tau_e)/sigma_e * (j_dp+1)*SQRT(p_dp)
ynepm1jm1(ip,ij) = +SQRT(tau_e)/sigma_e * j_dp*SQRT(p_dp) ynepp1jm1(ip,ij) = +SQRT(tau_e)/sigma_e * j_dp*SQRT(p_dp+1._dp)
zNepm1j (ip,ij) = -SQRT(tau_e)/sigma_e * (2._dp*j_dp+1_dp)*SQRT(p_dp) ynepm1jm1(ip,ij) = +SQRT(tau_e)/sigma_e * j_dp*SQRT(p_dp)
zNepm1jp1(ip,ij) = +SQRT(tau_e)/sigma_e * (j_dp+1_dp)*SQRT(p_dp) zNepm1j (ip,ij) = +SQRT(tau_e)/sigma_e * (2._dp*j_dp+1_dp)*SQRT(p_dp)
zNepm1jm1(ip,ij) = +SQRT(tau_e)/sigma_e * j_dp*SQRT(p_dp) zNepm1jp1(ip,ij) = -SQRT(tau_e)/sigma_e * (j_dp+1_dp)*SQRT(p_dp)
zNepm1jm1(ip,ij) = -SQRT(tau_e)/sigma_e * j_dp*SQRT(p_dp)
ENDDO ENDDO
ENDDO ENDDO
DO ip = ips_e, ipe_e DO ip = ips_e, ipe_e
p_int= parray_e(ip) ! Hermite degree p_int= parray_e(ip) ! Hermite degree
p_dp = REAL(p_int,dp) ! REAL of Hermite degree p_dp = REAL(p_int,dp) ! REAL of Hermite degree
xnepp1j(ip) = sqrtTaue_qe * SQRT(p_dp + 1_dp) ! Landau damping coefficients (ddz napj term)
xnepm1j(ip) = sqrtTaue_qe * SQRT(p_dp) xnepp1j(ip) = SQRT(tau_e)/sigma_e * SQRT(p_dp + 1_dp)
xnepp2j(ip) = taue_qe * SQRT((p_dp + 1._dp) * (p_dp + 2._dp)) xnepm1j(ip) = SQRT(tau_e)/sigma_e * SQRT(p_dp)
xnepm2j(ip) = taue_qe * SQRT(p_dp * (p_dp - 1._dp)) ! Magnetic curvature term
xnepp2j(ip) = taue_qe * CurvB * SQRT((p_dp + 1._dp) * (p_dp + 2._dp))
xnepm2j(ip) = taue_qe * CurvB * SQRT(p_dp * (p_dp - 1._dp))
ENDDO ENDDO
DO ij = ijs_e, ije_e DO ij = ijs_e, ije_e
j_int= jarray_e(ij) ! Laguerre degree j_int= jarray_e(ij) ! Laguerre degree
j_dp = REAL(j_int,dp) ! REAL of Laguerre degree j_dp = REAL(j_int,dp) ! REAL of Laguerre degree
xnepjp1(ij) = -taue_qe * (j_dp + 1._dp) ! Magnetic gradient term
xnepjm1(ij) = -taue_qe * j_dp xnepjp1(ij) = -taue_qe * GradB * (j_dp + 1._dp)
xnepjm1(ij) = -taue_qe * GradB * j_dp
ENDDO ENDDO
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!! Ions linear coefficients for moment RHS !!!!!!!!!! !! Ions linear coefficients for moment RHS !!!!!!!!!!
...@@ -207,30 +141,36 @@ SUBROUTINE compute_lin_coeff ...@@ -207,30 +141,36 @@ SUBROUTINE compute_lin_coeff
DO ij = ijs_i, ije_i DO ij = ijs_i, ije_i
j_int= jarray_i(ij) ! Laguerre degree j_int= jarray_i(ij) ! Laguerre degree
j_dp = REAL(j_int,dp) ! REAL of Laguerre degree j_dp = REAL(j_int,dp) ! REAL of Laguerre degree
xnepj(ip,ij) = taui_qi*(CurvB*(2._dp*p_dp + 1._dp) & ! All Napj terms
xnipj(ip,ij) = taui_qi*(CurvB*(2._dp*p_dp + 1._dp) &
+GradB*(2._dp*j_dp + 1._dp)) +GradB*(2._dp*j_dp + 1._dp))
ynipp1j (ip,ij) = -SQRT(tau_i)/sigma_i * (j_dp+1)*SQRT(p_dp+1._dp) ! Mirror force terms
ynipm1j (ip,ij) = -SQRT(tau_i)/sigma_i * (j_dp+1)*SQRT(p_dp) ynipp1j (ip,ij) = -SQRT(tau_i)/sigma_i* (j_dp+1)*SQRT(p_dp+1._dp)
ynipp1jm1(ip,ij) = +SQRT(tau_i)/sigma_i * j_dp*SQRT(p_dp+1._dp) ynipm1j (ip,ij) = -SQRT(tau_i)/sigma_i* (j_dp+1)*SQRT(p_dp)
ynipm1jm1(ip,ij) = +SQRT(tau_i)/sigma_i * j_dp*SQRT(p_dp) ynipp1jm1(ip,ij) = +SQRT(tau_i)/sigma_i* j_dp*SQRT(p_dp+1._dp)
zNipm1j (ip,ij) = -SQRT(tau_i)/sigma_i * (2._dp*j_dp+1_dp)*SQRT(p_dp) ynipm1jm1(ip,ij) = +SQRT(tau_i)/sigma_i* j_dp*SQRT(p_dp)
zNipm1jp1(ip,ij) = +SQRT(tau_i)/sigma_i * (j_dp+1_dp)*SQRT(p_dp) ! Trapping terms
zNipm1jm1(ip,ij) = +SQRT(tau_i)/sigma_i * j_dp*SQRT(p_dp) zNipm1j (ip,ij) = +SQRT(tau_i)/sigma_i* (2._dp*j_dp+1_dp)*SQRT(p_dp)
zNipm1jp1(ip,ij) = -SQRT(tau_i)/sigma_i* (j_dp+1_dp)*SQRT(p_dp)
zNipm1jm1(ip,ij) = -SQRT(tau_i)/sigma_i* j_dp*SQRT(p_dp)
ENDDO ENDDO
ENDDO ENDDO
DO ip = ips_i, ipe_i DO ip = ips_i, ipe_i
p_int= parray_i(ip) ! Hermite degree p_int= parray_i(ip) ! Hermite degree
p_dp = REAL(p_int,dp) ! REAL of Hermite degree p_dp = REAL(p_int,dp) ! REAL of Hermite degree
xnipp1j(ip) = sqrtTaui_qi * SQRT(p_dp + 1._dp) ! Landau damping coefficients (ddz napj term)
xnipm1j(ip) = sqrtTaui_qi * SQRT(p_dp) xnipp1j(ip) = SQRT(tau_i)/sigma_i * SQRT(p_dp + 1._dp)
xnipp2j(ip) = taui_qi * SQRT((p_dp + 1._dp) * (p_dp + 2._dp)) xnipm1j(ip) = SQRT(tau_i)/sigma_i * SQRT(p_dp)
xnipm2j(ip) = taui_qi * SQRT(p_dp * (p_dp - 1._dp)) ! Magnetic curvature term
xnipp2j(ip) = taui_qi * CurvB * SQRT((p_dp + 1._dp) * (p_dp + 2._dp))
xnipm2j(ip) = taui_qi * CurvB * SQRT(p_dp * (p_dp - 1._dp))
ENDDO ENDDO
DO ij = ijs_i, ije_i DO ij = ijs_i, ije_i
j_int= jarray_i(ij) ! Laguerre degree j_int= jarray_i(ij) ! Laguerre degree
j_dp = REAL(j_int,dp) ! REAL of Laguerre degree j_dp = REAL(j_int,dp) ! REAL of Laguerre degree
xnipjp1(ij) = -taui_qi * (j_dp + 1._dp) ! Magnetic gradient term
xnipjm1(ij) = -taui_qi * j_dp xnipjp1(ij) = -taui_qi * GradB * (j_dp + 1._dp)
xnipjm1(ij) = -taui_qi * GradB * j_dp
ENDDO ENDDO
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!! ES linear coefficients for moment RHS !!!!!!!!!! !! ES linear coefficients for moment RHS !!!!!!!!!!
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment