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

add full coulomb operator and completion with dougherty

parent a5564278
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
Showing
with 48 additions and 19 deletions
src/collision_mod.F90
+ 48
− 19
View file @ a1ac1cd5
...@@ -464,7 +464,7 @@ CONTAINS ...@@ -464,7 +464,7 @@ CONTAINS
USE basic USE basic
USE time_integration, ONLY : updatetlevel USE time_integration, ONLY : updatetlevel
USE utility USE utility
USE model, ONLY: CO, NON_LIN USE model, ONLY: CO, NON_LIN, sigmae2_taue_o2, sigmai2_taui_o2
IMPLICIT NONE IMPLICIT NONE
! Indices for row and columns of the COSOlver matrix (4D compressed 2D matrices) ! Indices for row and columns of the COSOlver matrix (4D compressed 2D matrices)
INTEGER :: irow_sub, irow_full, icol_sub, icol_full INTEGER :: irow_sub, irow_full, icol_sub, icol_full
...@@ -488,7 +488,7 @@ CONTAINS ...@@ -488,7 +488,7 @@ CONTAINS
REAL(dp), DIMENSION(:), ALLOCATABLE :: kp_grid_mat ! kperp grid of the matrices REAL(dp), DIMENSION(:), ALLOCATABLE :: kp_grid_mat ! kperp grid of the matrices
INTEGER :: ikp_next, ikp_prev, nkp_mat, ikp_mat INTEGER :: ikp_next, ikp_prev, nkp_mat, ikp_mat
REAL(dp) :: kp_next, kp_prev, kperp_sim, kperp_mat, zerotoone REAL(dp) :: kp_next, kp_prev, kperp_sim, kperp_mat, zerotoone, be_2, bi_2
CHARACTER(len=128) :: var_name, kperp_string, ikp_string CHARACTER(len=128) :: var_name, kperp_string, ikp_string
...@@ -499,10 +499,14 @@ CONTAINS ...@@ -499,10 +499,14 @@ CONTAINS
IF (my_id .EQ. 0) WRITE(*,*) '=== Load GK Sugama matrix ===' IF (my_id .EQ. 0) WRITE(*,*) '=== Load GK Sugama matrix ==='
ELSEIF(CO .EQ. 3) THEN ELSEIF(CO .EQ. 3) THEN
IF (my_id .EQ. 0) WRITE(*,*) '=== Load GK pitch angle matrix ===' IF (my_id .EQ. 0) WRITE(*,*) '=== Load GK pitch angle matrix ==='
ELSEIF(CO .EQ. 4) THEN
IF (my_id .EQ. 0) WRITE(*,*) '=== Load GK Coulomb matrix ==='
ELSEIF(CO .EQ. -2) THEN ELSEIF(CO .EQ. -2) THEN
IF (my_id .EQ. 0) WRITE(*,*) '=== Load DK Sugama matrix ===' IF (my_id .EQ. 0) WRITE(*,*) '=== Load DK Sugama matrix ==='
ELSEIF(CO .EQ. -3) THEN ELSEIF(CO .EQ. -3) THEN
IF (my_id .EQ. 0) WRITE(*,*) '=== Load DK pitch angle matrix ===' IF (my_id .EQ. 0) WRITE(*,*) '=== Load DK pitch angle matrix ==='
ELSEIF(CO .EQ. -4) THEN
IF (my_id .EQ. 0) WRITE(*,*) '=== Load DK Coulomb matrix ==='
ENDIF ENDIF
! Opening the compiled cosolver matrices results ! Opening the compiled cosolver matrices results
...@@ -540,7 +544,7 @@ CONTAINS ...@@ -540,7 +544,7 @@ CONTAINS
CALL allocate_array( CiepjF_kp, 1,(pmaxe+1)*(jmaxe+1), 1,(pmaxe+1)*(jmaxe+1), ikps_C,ikpe_C) CALL allocate_array( CiepjF_kp, 1,(pmaxe+1)*(jmaxe+1), 1,(pmaxe+1)*(jmaxe+1), ikps_C,ikpe_C)
DO ikp = ikps_C,ikpe_C ! Loop over everz kperp values DO ikp = ikps_C,ikpe_C ! Loop over everz kperp values
! we put zeros if kp>2/3kpmax because thoses frequenvies are filtered through AA ! 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) THEN
CiepjT_kp(:,:,ikp) = 0._dp CiepjT_kp(:,:,ikp) = 0._dp
CiepjF_kp(:,:,ikp) = 0._dp CiepjF_kp(:,:,ikp) = 0._dp
...@@ -682,23 +686,48 @@ CONTAINS ...@@ -682,23 +686,48 @@ CONTAINS
DO ikp=1,nkp_mat DO ikp=1,nkp_mat
ikp_mat = ikp ! the first indice of the interval (k0) ikp_mat = ikp ! the first indice of the interval (k0)
kperp_mat = kp_grid_mat(ikp) kperp_mat = kp_grid_mat(ikp)
IF(kperp_mat .GT. kperp_sim) EXIT IF(kperp_mat .GT. kperp_sim) EXIT ! a matrix with kperp2 > current kx2 + ky2 has been found
ENDDO ENDDO
! interval boundaries IF((abs(CO) .eq. 4) .AND. (kperp_sim .GT. 3.0)) THEN ! hybrid CO for Coulomb collisions
ikp_prev = ikp_mat - 1 !index of k0 ! write(*,*) kp_grid_mat(ikp_mat), '<', kperp_sim
ikp_next = ikp_mat !index of k1 !! Fill non existing matrices with diagonal Dougherty damping
! write(*,*) kp_grid_mat(ikp_prev), '<', kperp_sim, '<', kp_grid_mat(ikp_next) DO ip_e = 0,Pmaxe ! Loop over rows
DO ij_e = 0,Jmaxe
! 0->1 variable for linear interp, i.e. zero2one = (k-k0)/(k1-k0) irow_sub = (Jmaxe +1)*ip_e + ij_e +1
zerotoone = (kperp_sim - kp_grid_mat(ikp_prev))/(kp_grid_mat(ikp_next) - kp_grid_mat(ikp_prev)) be_2 = (kxarray(ikx)**2 + kyarray(iky)**2) * sigmae2_taue_o2
Ceepj (irow_sub,irow_sub,ikx,iky) = -(real(ip_e,dp) + 2._dp*real(ij_e,dp) + 2._dp*be_2)!Ceepj__kp(:,:,nkp_mat)
! Linear interpolation between previous and next kperp matrix values ENDDO
Ceepj (:,:,ikx,iky) = (Ceepj__kp(:,:,ikp_next) - Ceepj__kp(:,:,ikp_prev))*zerotoone + Ceepj__kp(:,:,ikp_prev) ENDDO
CeipjT(:,:,ikx,iky) = (CeipjT_kp(:,:,ikp_next) - CeipjT_kp(:,:,ikp_prev))*zerotoone + CeipjT_kp(:,:,ikp_prev) DO ip_i = 0,Pmaxi ! Loop over rows
CeipjF(:,:,ikx,iky) = (CeipjF_kp(:,:,ikp_next) - CeipjF_kp(:,:,ikp_prev))*zerotoone + CeipjF_kp(:,:,ikp_prev) DO ij_i = 0,Jmaxi
Ciipj (:,:,ikx,iky) = (Ciipj__kp(:,:,ikp_next) - Ciipj__kp(:,:,ikp_prev))*zerotoone + Ciipj__kp(:,:,ikp_prev) irow_sub = (Jmaxi +1)*ip_i + ij_i +1
CiepjT(:,:,ikx,iky) = (CiepjT_kp(:,:,ikp_next) - CiepjT_kp(:,:,ikp_prev))*zerotoone + CiepjT_kp(:,:,ikp_prev) bi_2 = (kxarray(ikx)**2 + kyarray(iky)**2) * sigmai2_taui_o2
CiepjF(:,:,ikx,iky) = (CiepjF_kp(:,:,ikp_next) - CiepjF_kp(:,:,ikp_prev))*zerotoone + CiepjF_kp(:,:,ikp_prev) Ciipj (irow_sub,irow_sub,ikx,iky) = -(real(ip_i,dp) + 2._dp*real(ij_i,dp) + 2._dp*bi_2)!0*Ciipj__kp(:,:,nkp_mat)
ENDDO
ENDDO
! Ceepj (:,:,ikx,iky) = 0*Ceepj__kp(:,:,nkp_mat)
! Ciipj (:,:,ikx,iky) = 0*Ciipj__kp(:,:,nkp_mat)
CeipjT(:,:,ikx,iky) = 0*CeipjT_kp(:,:,nkp_mat)
CeipjF(:,:,ikx,iky) = 0*CeipjF_kp(:,:,nkp_mat)
CiepjT(:,:,ikx,iky) = 0*CiepjT_kp(:,:,nkp_mat)
CiepjF(:,:,ikx,iky) = 0*CiepjF_kp(:,:,nkp_mat)
ELSE ! Interpolation
! interval boundaries
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)
! 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))
! Linear interpolation between previous and next kperp matrix values
Ceepj (:,:,ikx,iky) = (Ceepj__kp(:,:,ikp_next) - Ceepj__kp(:,:,ikp_prev))*zerotoone + Ceepj__kp(:,:,ikp_prev)
CeipjT(:,:,ikx,iky) = (CeipjT_kp(:,:,ikp_next) - CeipjT_kp(:,:,ikp_prev))*zerotoone + CeipjT_kp(:,:,ikp_prev)
CeipjF(:,:,ikx,iky) = (CeipjF_kp(:,:,ikp_next) - CeipjF_kp(:,:,ikp_prev))*zerotoone + CeipjF_kp(:,:,ikp_prev)
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)
ENDIF
ENDDO ENDDO
ENDDO ENDDO
ELSE ! DK -> No kperp dep, copy simply to final collision matrices ELSE ! DK -> No kperp dep, copy simply to final collision matrices
......
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