subroutine include_CoulombCollisions_GCfio_p(tt,params,Y_R,Y_PHI,Y_Z, &
Ppll,Pmu,me,flagCon,flagCol,F,P,E_PHI,ne,ni,Te,Zeff,nimp,PSIp,hint)
TYPE(PROFILES), INTENT(IN) :: P
TYPE(FIELDS), INTENT(IN) :: F
TYPE(KORC_PARAMS), INTENT(INOUT) :: params
REAL(rp), DIMENSION(params%pchunk), INTENT(INOUT) :: Ppll
REAL(rp), DIMENSION(params%pchunk), INTENT(INOUT) :: Pmu
REAL(rp), DIMENSION(params%pchunk) :: Bmag
REAL(rp), DIMENSION(params%pchunk) :: B_R,B_PHI,B_Z
REAL(rp), DIMENSION(params%pchunk) :: curlb_R,curlb_PHI,curlb_Z
REAL(rp), DIMENSION(params%pchunk) :: gradB_R,gradB_PHI,gradB_Z
REAL(rp), DIMENSION(params%pchunk) :: E_R,E_Z
REAL(rp), DIMENSION(params%pchunk) :: E_PHI,PSIp,E_PHI0
REAL(rp), DIMENSION(params%pchunk), INTENT(INOUT) :: ne,Te,ni
REAL(rp), DIMENSION(params%pchunk,params%num_impurity_species),&
& INTENT(INOUT) :: nimp
REAL(rp), DIMENSION(params%pchunk), INTENT(IN) :: Y_R,Y_PHI&
&,Y_Z
INTEGER(is), DIMENSION(params%pchunk), INTENT(INOUT) :: flagCol
INTEGER(is), DIMENSION(params%pchunk), INTENT(INOUT) :: flagCon
TYPE(C_PTR), DIMENSION(params%pchunk), INTENT(INOUT) :: hint
REAL(rp), INTENT(IN) :: me
REAL(rp), DIMENSION(params%pchunk), INTENT(OUT) :: Zeff
REAL(rp), DIMENSION(params%pchunk,2) :: dW
REAL(rp), DIMENSION(params%pchunk,2) :: rnd1
REAL(rp) :: dt,time
REAL(rp), DIMENSION(params%pchunk) :: pm
REAL(rp), DIMENSION(params%pchunk) :: dp
REAL(rp), DIMENSION(params%pchunk) :: xi
REAL(rp), DIMENSION(params%pchunk) :: dxi
REAL(rp), DIMENSION(params%pchunk) :: v,gam
!! speed of particle
REAL(rp), DIMENSION(params%pchunk) :: CAL
REAL(rp) , DIMENSION(params%pchunk) :: dCAL
REAL(rp), DIMENSION(params%pchunk) :: CFL
REAL(rp), DIMENSION(params%pchunk) :: CBL
REAL(rp), DIMENSION(params%pchunk) :: SC_p,SC_mu,BREM_p
REAL(rp) :: kappa
integer :: cc,pchunk,ii
integer(ip),INTENT(IN) :: tt
pchunk=params%pchunk
if (MODULO(params%it+tt,cparams_ss%subcycling_iterations) .EQ. 0_ip) then
dt = REAL(cparams_ss%subcycling_iterations,rp)*params%dt
call get_fio_GCmagnetic_fields_p(params,F,Y_R,Y_PHI,Y_Z, &
B_R,B_PHI,B_Z,gradB_R,gradB_PHI,gradB_Z, &
curlb_R,curlb_PHI,curlb_Z,flagCon,hint)
if (F%FIO_E .ge. 0) then
call get_fio_GCelectric_fields_p(params,F, &
Y_R,Y_PHI,Y_Z,E_R,E_PHI,E_Z,flagCon,hint)
end if
call get_fio_vector_potential_p(params,F,Y_R,Y_PHI,Y_Z, &
PSIp,flagCon,hint)
call get_fio_profile_p(params,P,Y_R,Y_PHI,Y_Z, &
ne,Te,flagCon,hint)
call get_fio_ion_p(params,P,Y_R,Y_PHI,Y_Z, &
ne,ni,nimp,Zeff,flagCon,hint)
!write(6,*) ne,Te,nimp,Zeff
!$OMP SIMD
! !$OMP& aligned (pm,xi,v,Ppll,Bmag,Pmu)
do cc=1_idef,pchunk
Bmag(cc)=sqrt(B_R(cc)*B_R(cc)+B_PHI(cc)*B_PHI(cc)+B_Z(cc)*B_Z(cc))
! Transform p_pll,mu to P,eta
pm(cc) = SQRT(Ppll(cc)*Ppll(cc)+2*me*Bmag(cc)*Pmu(cc))
xi(cc) = Ppll(cc)/pm(cc)
gam(cc) = sqrt(1+pm(cc)*pm(cc))
v(cc) = pm(cc)/gam(cc)
! normalized speed (v_K=v_P/c)
E_PHI0(cc)=E_PHI(cc)
end do
!$OMP END SIMD
!write(6,*) 'R',Y_R*params%cpp%length
!write(6,*) 'PHI',Y_PHI
!write(6,*) 'Z',Y_Z*params%cpp%length
!write(6,*) 'ne',ne*params%cpp%density
!write(6,*) 'ni',ni*params%cpp%density
!write(6,*) 'nimp',nimp(:,1:2)*params%cpp%density
!write(6,*) 'Zeff',Zeff
if (.not.params%FokPlan) E_PHI=0._rp
! write(output_unit_write,'("ne: "E17.10)') ne
! write(output_unit_write,'("Te: "E17.10)') Te
! write(output_unit_write,'("Bmag: "E17.10)') Bmag
! write(output_unit_write,'("v: ",E17.10)') v
! write(output_unit_write,'("xi: ",E17.10)') xi
! write(output_unit_write,'("size(E_PHI_GC): ",I16)') size(E_PHI)
!$OMP SIMD
! !$OMP& aligned(rnd1,dW,CAL,dCAL,CFL,CBL,v,ne,Te,Zeff,dp, &
! !$OMP& flagCon,flagCol,dxi,xi,pm,Ppll,Pmu,Bmag)
do cc=1_idef,pchunk
#ifdef PARALLEL_RANDOM
rnd1(cc,1) = get_random()
rnd1(cc,2) = get_random()
! rnd1(:,1) = get_random_mkl()
! rnd1(:,2) = get_random_mkl()
#else
call RANDOM_NUMBER(rnd1)
#endif
dW(cc,1) = SQRT(3*dt)*(-1+2*rnd1(cc,1))
dW(cc,2) = SQRT(3*dt)*(-1+2*rnd1(cc,2))
! write(output_unit_write,'("dW1: ",E17.10)') dW(cc,1)
! write(output_unit_write,'("dW2: ",E17.10)') dW(cc,2)
CAL(cc) = CA_SD(v(cc),ne(cc),Te(cc))
dCAL(cc)= dCA_SD(v(cc),me,ne(cc),Te(cc))
CFL(cc) = CF_SD_FIO(params,v(cc),ne(cc),Te(cc),nimp(cc,:))
CBL(cc) = (CB_ee_SD(v(cc),ne(cc),Te(cc),Zeff(cc))+ &
CB_ei_SD_FIO(params,v(cc),ne(cc),Te(cc),nimp(cc,:),Zeff(cc)))
dp(cc)=REAL(flagCol(cc))*REAL(flagCon(cc))* &
((-CFL(cc)+dCAL(cc)+E_PHI(cc)*xi(cc))*dt+ &
sqrt(2.0_rp*CAL(cc))*dW(cc,1))
dxi(cc)=REAL(flagCol(cc))*REAL(flagCon(cc))* &
((-2*xi(cc)*CBL(cc)/(pm(cc)*pm(cc))+ &
E_PHI(cc)*(1-xi(cc)*xi(cc))/pm(cc))*dt- &
sqrt(2.0_rp*CBL(cc)*(1-xi(cc)*xi(cc)))/pm(cc)*dW(cc,2))
! write(output_unit_write,'("dp: ",E17.10)') dp(cc)
! write(output_unit_write,'("dxi: ",E17.10)') dxi(cc)
end do
!$OMP END SIMD
if (params%FokPlan.and.params%radiation) then
if(params%GC_rad_model.eq.'SDE') then
!$OMP SIMD
do cc=1_idef,pchunk
SC_p(cc)=-gam(cc)*pm(cc)*(1-xi(cc)*xi(cc))/ &
(cparams_ss%taur/Bmag(cc)**2)
SC_mu(cc)=xi(cc)*(1-xi(cc)*xi(cc))/ &
((cparams_ss%taur/Bmag(cc)**2)*gam(cc))
kappa=2._rp*C_PI*C_RE**2._rp*C_ME*C_C**2._rp
BREM_p(cc)=-2._rp*ne(cc)*kappa*Zeff(cc)*(Zeff(cc)+1._rp)* &
C_a/C_PI*(gam(cc)-1._rp)*(log(2._rp*gam(cc))-1._rp/3._rp)
dp(cc)=dp(cc)+(SC_p(cc)+BREM_p(cc))*dt* &
REAL(flagCol(cc))*REAL(flagCon(cc))
dxi(cc)=dxi(cc)+(SC_mu(cc))*dt* &
REAL(flagCol(cc))*REAL(flagCon(cc))
end do
!$OMP END SIMD
end if
end if
!$OMP SIMD
do cc=1_idef,pchunk
pm(cc)=pm(cc)+dp(cc)
xi(cc)=xi(cc)+dxi(cc)
! if (pm(cc)<0) pm(cc)=-pm(cc)
! Keep xi between [-1,1]
if (xi(cc)>1) then
! write(output_unit_write,'("High xi at: ",E17.10," with dxi: ",E17.10)') &
! time*params%cpp%time, dxi(cc)
xi(cc)=1-mod(xi(cc),1._rp)
else if (xi(cc)<-1) then
xi(cc)=-1-mod(xi(cc),-1._rp)
! write(output_unit_write,'("Low xi at: ",E17.10," with dxi: ",E17.10)') &
! time*params%cpp%time, dxi(cc)
endif
! Transform P,xi to p_pll,mu
Ppll(cc)=pm(cc)*xi(cc)
Pmu(cc)=(pm(cc)*pm(cc)-Ppll(cc)*Ppll(cc))/(2*me*Bmag(cc))
end do
!$OMP END SIMD
#if DBG_CHECK
do cc=1_idef,pchunk
if((isnan(Ppll(cc)).or.isnan(Pmu(cc))).and. &
((flagCol(cc).eq.1_is).and.(flagCon(cc).eq.1_is))) then
write(6,*) 'End collision'
write(6,*) 'Ppll',Ppll(cc)
write(6,*) 'Pmu',Pmu(cc)
write(6,*) 'Bmag',Bmag(cc)
write(6,*) 'pm',pm(cc)
write(6,*) 'xi',xi(cc)
write(6,*) 'dp',dp(cc)
write(6,*) 'dxi',dxi(cc)
write(6,*) 'CFL',CFL(cc)
write(6,*) 'CBL',CBL(cc)
write(6,*) 'dCAL',dCAL(cc)
write(6,*) 'CAL',CAL(cc)
write(6,*) 'v',v(cc)
write(6,*) 'ne',ne(cc)
write(6,*) 'Te',Te(cc)
write(6,*) 'nimp',nimp(cc,:)
write(6,*) 'Zeff',Zeff(cc)
write(6,*) 'Y_R',Y_R(cc)
write(6,*) 'Y_PHI',Y_PHI(cc)
write(6,*) 'Y_Z',Y_Z(cc)
stop 'Ppll or Pmu is NaN'
endif
end do
#endif
! write(output_unit_write,'("rnd1: ",E17.10)') rnd1
! write(output_unit_write,'("flag: ",I16)') flag
! write(output_unit_write,'("CA: ",E17.10)') CAL
! write(output_unit_write,'("dCA: ",E17.10)') dCAL
! write(output_unit_write,'("CF ",E17.10)') CFL
! write(output_unit_write,'("CB: ",E17.10)') CBL
! write(output_unit_write,'("dp: ",E17.10)') dp
! write(output_unit_write,'("dxi: ",E17.10)') dxi
! write(output_unit_write,'("Ppll: ",E17.10)') Ppll
! write(output_unit_write,'("Pmu: ",E17.10)') Pmu
! write(output_unit_write,'("E_PHI_COL: ",E17.10)') E_PHI
do cc=1_idef,pchunk
if ((pm(cc).lt.min(cparams_ss%p_min*cparams_ss%pmin_scale, &
p_therm)).and.flagCol(cc).eq.1_ip) then
! write(output_unit_write,'("Momentum less than zero")')
! stop
! write(output_unit_write,'("Particle not tracked at: ",E17.10," &
! & with xi: ",E17.10)') time*params%cpp%time, xi(cc)
flagCol(cc)=0_ip
end if
end do
! if (tt .EQ. 1_ip) then
! write(output_unit_write,'("dp_rad: ",E17.10)') &
! -gam(1)*pm(1)*(1-xi(1)*xi(1))/ &
! (cparams_ss%taur/Bmag(1)**2)*dt
! write(output_unit_write,'("dxi_rad: ",E17.10)') &
! xi(1)*(1-xi(1)*xi(1))/ &
! ((cparams_ss%taur/Bmag(1)**2)*gam(1))*dt
! end if
! if (tt .EQ. 1_ip) then
! write(output_unit_write,'("CA: ",E17.10)') CAL(1)
! write(output_unit_write,'("dCA: ",E17.10)') dCAL(1)
! write(output_unit_write,'("CF ",E17.10)') CFL(1)
! write(output_unit_write,'("CB: ",E17.10)') CBL(1)
! end if
if (.not.params%FokPlan) E_PHI=E_PHI0
end if
end subroutine include_CoulombCollisions_GCfio_p