diff --git a/src/collision_mod.F90 b/src/collision_mod.F90 index 17283101fb9b9bbc418eb5d63036eef729a0fc1b..076792b78d3cc6aeab3935366abda24346c3266a 100644 --- a/src/collision_mod.F90 +++ b/src/collision_mod.F90 @@ -17,133 +17,6 @@ PUBLIC :: apply_COSOlver_mat_e, apply_COSOlver_mat_i CONTAINS - ! !******************************************************************************! - ! !! Doughtery gyrokinetic collision operator for electrons - ! !******************************************************************************! - ! SUBROUTINE DoughertyGK(ip_,ij_,ikr_,ikz_,TColl_,specie_) - ! IMPLICIT NONE - ! INTEGER, INTENT(IN) :: ip_,ij_,ikr_,ikz_ - ! CHARACTER(len = 1), INTENT(IN) :: specie_ - ! COMPLEX(dp), INTENT(INOUT) :: TColl_ - ! - ! COMPLEX(dp) :: n_,upar_,uperp_,Tpar_, Tperp_ - ! COMPLEX(dp) :: Dpj, Ppj, T_ - ! COMPLEX(dp) :: nadiab_moment_0j - ! COMPLEX(dp), DIMENSION(:,:), ALLOCATABLE :: moments_ - ! REAL(dp), DIMENSION(:), ALLOCATABLE :: kernel_ - ! REAL(dp) :: Knp0, Knp1, Knm1 - ! INTEGER :: in_, jmax_ - ! REAL(dp) :: n_dp, j_dp, p_dp, b_, bo2_2_, q_tau_, nu_ - ! - ! !** Auxiliary variables ** - ! !! If electrons !! - ! IF ( specie_ .EQ. 'e' ) THEN - ! p_dp = REAL(parray_e(ip_),dp) - ! j_dp = REAL(jarray_e(ij_),dp) - ! jmax_ = jmaxe - ! bo2_2_ = (krarray(ikr_)**2 + kzarray(ikz_)**2) * sigmae2_taue_o2 ! this is (be/2)^2 - ! b_ = 2_dp*SQRT(bo2_2_) ! this is be - ! q_tau_ = q_e/tau_e - ! nu_ = nu_ee - ! CALL allocate_array(moments_, ips_e,ipe_e, ijs_e,ije_e) - ! moments_(ips_e:ipe_e,ijs_e:ije_e) = moments_e(ips_e:ipe_e,ijs_e:ije_e,ikr_,ikz_,updatetlevel) - ! CALL allocate_array(kernel_, ijs_e,ije_e) - ! kernel_(ijsg_e:ijeg_e) = kernel_e(ijsg_e:ijeg_e,ikr_,ikz_) - ! !! If ions !! - ! ELSEIF ( specie_ .EQ. 'i') THEN - ! p_dp = REAL(parray_i(ip_),dp) - ! j_dp = REAL(jarray_i(ij_),dp) - ! jmax_ = jmaxi - ! bo2_2_ = (krarray(ikr_)**2 + kzarray(ikz_)**2) * sigmai2_taui_o2 ! this is (bi/2)^2 - ! b_ = 2_dp*SQRT(bo2_2_) ! this is be - ! q_tau_ = q_i/tau_i - ! nu_ = nu_i - ! CALL allocate_array(moments_, ips_i,ipe_i, ijs_i,ije_i) - ! moments_(ips_i:ipe_i,ijs_i:ije_i) = moments_i(ips_i:ipe_i,ijs_i:ije_i,ikr_,ikz_,updatetlevel) - ! CALL allocate_array(kernel_, ijsg_i,ijeg_i) - ! kernel_(ijsg_i:ijeg_i) = kernel_i(ijsg_i:ijeg_i,ikr_,ikz_) - ! ENDIF - ! - ! !** Assembling collison operator ** - ! ! Velocity-space diffusion (similar to Lenhard Bernstein) - ! ! -nuee (p + 2j + b^2/2) Nepj - ! TColl_ = -(p_dp + 2._dp*j_dp + 2._dp*bo2_2_)*moments_(ip_,ij_) - ! - ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! Non zero term for p = 0 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - ! IF( p_dp .EQ. 0 ) THEN ! Kronecker p0 - ! ! Get adiabatic moment - ! TColl_ = TColl_ - (p_dp + 2._dp*j_dp + 2._dp*bo2_2_) * q_tau_ * kernel_(ij_)*phi(ikr_,ikz_) - ! !** build required fluid moments ** - ! n_ = 0._dp - ! upar_ = 0._dp; uperp_ = 0._dp - ! Tpar_ = 0._dp; Tperp_ = 0._dp - ! DO in_ = 1,jmaxe+1 - ! n_dp = REAL(in_-1,dp) - ! ! Store the kernels for sparing readings - ! Knp0 = kernel_(in_) - ! Knp1 = kernel_(in_+1) - ! Knm1 = kernel_(in_-1) - ! ! Nonadiabatic moments (only different from moments when p=0) - ! nadiab_moment_0j = moments_(1,in_) + q_tau_ * Knp0 *phi(ikr_,ikz_) - ! ! Density - ! n_ = n_ + Knp0 * nadiab_moment_0j - ! ! Perpendicular velocity - ! uperp_ = uperp_ + b_*0.5_dp*(Knp0 - Knm1) * nadiab_moment_0j - ! ! Parallel temperature - ! Tpar_ = Tpar_ + Knp0 * (SQRT2*moments_(3,in_) + nadiab_moment_0j) - ! ! Perpendicular temperature - ! Tperp_ = Tperp_ + ((2._dp*n_dp+1._dp)*Knp0 - (n_dp+1._dp)*Knp1 - n_dp*Knm1)*nadiab_moment_0j - ! ENDDO - ! T_ = (Tpar_ + 2._dp*Tperp_)/3._dp - n_ - ! ! Add energy restoring term - ! TColl_ = TColl_ + T_* 4._dp * j_dp * kernel_(ij_) - ! TColl_ = TColl_ - T_* 2._dp * (j_dp + 1._dp) * kernel_(ij_+1) - ! TColl_ = TColl_ - T_* 2._dp * j_dp * kernel_(ij_-1) - ! TColl_ = TColl_ + uperp_*b_* (kernel_(ij_) - kernel_(ij_-1)) - ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - ! - ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! Non zero term for p = 1 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - ! ELSEIF( p_dp .eq. 1 ) THEN ! kronecker p1 - ! !** build required fluid moments ** - ! upar_ = 0._dp - ! DO in_ = 1,jmax_+1 - ! ! Parallel velocity - ! upar_ = upar_ + Kernel_(in_) * moments_(2,in_) - ! ENDDO - ! TColl_ = TColl_ + upar_*Kernel_(ij_) - ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - ! - ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! Non zero term for p = 2 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - ! ELSEIF( p_dp .eq. 2 ) THEN ! kronecker p2 - ! !** build required fluid moments ** - ! n_ = 0._dp - ! upar_ = 0._dp; uperp_ = 0._dp - ! Tpar_ = 0._dp; Tperp_ = 0._dp - ! DO in_ = 1,jmaxe+1 - ! n_dp = REAL(in_-1,dp) - ! ! Store the kernels for sparing readings - ! Knp0 = kernel_(in_) - ! Knp1 = kernel_(in_+1) - ! Knm1 = kernel_(in_-1) - ! ! Nonadiabatic moments (only different from moments when p=0) - ! nadiab_moment_0j = moments_(1,in_) + q_tau_*Knp0*phi(ikr_,ikz_) - ! ! Density - ! n_ = n_ + Knp0 * nadiab_moment_0j - ! ! Parallel temperature - ! Tpar_ = Tpar_ + Knp0 * (SQRT2*moments_(3,in_) + nadiab_moment_0j) - ! ! Perpendicular temperature - ! Tperp_ = Tperp_ + ((2._dp*n_dp+1._dp)*Knp0 - (n_dp+1) * Knp1 - n_dp * Knm1)*nadiab_moment_0j - ! ENDDO - ! T_ = (Tpar_ + 2._dp*Tperp_)/3._dp - n_ - ! TColl_ = TColl_ + T_*SQRT2*kernel_(ij_) - ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - ! - ! ENDIF - ! ! Multiply by specieslike collision coefficient - ! TColl_ = nu_ * TColl_ - ! - ! END SUBROUTINE DoughertyGK - !******************************************************************************! !! Doughtery gyrokinetic collision operator for electrons !******************************************************************************!