This subroutine performs a Stochastic collision process consistent with the Fokker-Planck model for relativitic electron colliding with a thermal (Maxwellian) plasma. The collision operator is in spherical coordinates of the form found in Papp et al., NF (2011). CA corresponds to the parallel (speed diffusion) process, CF corresponds to a slowing down (momentum loss) process, and CB corresponds to a perpendicular diffusion process. Ordering of the processes are and only the dominant terms are kept.
| Type | Intent | Optional | Attributes | Name | ||
|---|---|---|---|---|---|---|
| integer(kind=ip), | intent(in) | :: | tt | |||
| type(KORC_PARAMS), | intent(in) | :: | params | |||
| real(kind=rp), | intent(in), | DIMENSION(8) | :: | X_X | ||
| real(kind=rp), | intent(in), | DIMENSION(8) | :: | X_Y | ||
| real(kind=rp), | intent(in), | DIMENSION(8) | :: | X_Z | ||
| real(kind=rp), | intent(inout), | DIMENSION(8) | :: | U_X | ||
| real(kind=rp), | intent(inout), | DIMENSION(8) | :: | U_Y | ||
| real(kind=rp), | intent(inout), | DIMENSION(8) | :: | U_Z | ||
| real(kind=rp), | intent(in), | DIMENSION(8) | :: | B_X | ||
| real(kind=rp), | intent(in), | DIMENSION(8) | :: | B_Y | ||
| real(kind=rp), | intent(in), | DIMENSION(8) | :: | B_Z | ||
| real(kind=rp), | intent(in) | :: | me | |||
| type(PROFILES), | intent(in) | :: | P | |||
| type(FIELDS), | intent(in) | :: | F | |||
| integer(kind=is), | intent(inout), | DIMENSION(8) | :: | flag | ||
| real(kind=rp), | intent(in), | DIMENSION(8) | :: | PSIp |
subroutine include_CoulombCollisions_FO_p(tt,params,X_X,X_Y,X_Z, &
U_X,U_Y,U_Z,B_X,B_Y,B_Z,me,P,F,flag,PSIp)
!! This subroutine performs a Stochastic collision process consistent
!! with the Fokker-Planck model for relativitic electron colliding with
!! a thermal (Maxwellian) plasma. The collision operator is in spherical
!! coordinates of the form found in Papp et al., NF (2011). CA
!! corresponds to the parallel (speed diffusion) process, CF corresponds
!! to a slowing down (momentum loss) process, and CB corresponds to a
!! perpendicular diffusion process. Ordering of the processes are
!! $$ \sqrt{CB}\gg CB \gg CF \sim \sqrt{CA} \gg CA,$$
!! and only the dominant terms are kept.
TYPE(PROFILES), INTENT(IN) :: P
TYPE(FIELDS), INTENT(IN) :: F
TYPE(KORC_PARAMS), INTENT(IN) :: params
REAL(rp), DIMENSION(8), INTENT(IN) :: X_X,X_Y,X_Z,PSIp
REAL(rp), DIMENSION(8) :: Y_R,Y_PHI,Y_Z
REAL(rp), DIMENSION(8), INTENT(INOUT) :: U_X,U_Y,U_Z
REAL(rp), DIMENSION(8) :: ne,Te,Zeff
INTEGER(is), DIMENSION(8), INTENT(INOUT) :: flag
REAL(rp), INTENT(IN) :: me
INTEGER(ip), INTENT(IN) :: tt
REAL(rp), DIMENSION(8), INTENT(IN) :: B_X,B_Y,B_Z
REAL(rp), DIMENSION(8) :: b_unit_X,b_unit_Y,b_unit_Z
REAL(rp), DIMENSION(8) :: b1_X,b1_Y,b1_Z
REAL(rp), DIMENSION(8) :: b2_X,b2_Y,b2_Z
REAL(rp), DIMENSION(8) :: b3_X,b3_Y,b3_Z
REAL(rp), DIMENSION(8) :: Bmag
REAL(rp), DIMENSION(8,3) :: dW
!! 3D Weiner process
REAL(rp), DIMENSION(8,3) :: rnd1
REAL(rp) :: dt,time
REAL(rp), DIMENSION(8) :: um
REAL(rp), DIMENSION(8) :: dpm
REAL(rp), DIMENSION(8) :: vm
REAL(rp), DIMENSION(8) :: pm
REAL(rp),DIMENSION(8) :: Ub_X,Ub_Y,Ub_Z
REAL(rp), DIMENSION(8) :: xi
REAL(rp), DIMENSION(8) :: dxi
REAL(rp), DIMENSION(8) :: phi
REAL(rp), DIMENSION(8) :: dphi
!! speed of particle
REAL(rp),DIMENSION(8) :: CAL
REAL(rp),DIMENSION(8) :: dCAL
REAL(rp),DIMENSION(8) :: CFL
REAL(rp),DIMENSION(8) :: CBL
integer(ip) :: cc
if (MODULO(params%it+tt,cparams_ss%subcycling_iterations) .EQ. 0_ip) then
dt = REAL(cparams_ss%subcycling_iterations,rp)*params%dt
time=params%init_time+(params%it-1+tt)*params%dt
! subcylcling iterations a fraction of fastest collision frequency,
! where fraction set by dTau in namelist &CollisionParamsSingleSpecies
call cart_to_cyl_p(X_X,X_Y,X_Z,Y_R,Y_PHI,Y_Z)
if (params%profile_model(1:10).eq.'ANALYTICAL') then
call analytical_profiles_p(time,params,Y_R,Y_Z,P,F,ne,Te,Zeff,PSIp)
else if (params%profile_model(1:8).eq.'EXTERNAL') then
call interp_FOcollision_p(Y_R,Y_PHI,Y_Z,ne,Te,Zeff,flag)
end if
!$OMP SIMD
! !$OMP& aligned(um,pm,vm,U_X,U_Y,U_Z,Bmag,B_X,B_Y,B_Z, &
! !$OMP& b_unit_X,b_unit_Y,b_unit_Z,xi)
do cc=1_idef,8_idef
um(cc) = SQRT(U_X(cc)*U_X(cc)+U_Y(cc)*U_Y(cc)+U_Z(cc)*U_Z(cc))
pm(cc)=me*um(cc)
vm(cc) = um(cc)/SQRT(1.0_rp + um(cc)*um(cc))
! um is gamma times v, this solves for v
Bmag(cc)= SQRT(B_X(cc)*B_X(cc)+B_Y(cc)*B_Y(cc)+B_Z(cc)*B_Z(cc))
b_unit_X(cc)=B_X(cc)/Bmag(cc)
b_unit_Y(cc)=B_Y(cc)/Bmag(cc)
b_unit_Z(cc)=B_Z(cc)/Bmag(cc)
xi(cc)=(U_X(cc)*b_unit_X(cc)+U_Y(cc)*b_unit_Y(cc)+ &
U_Z(cc)*b_unit_Z(cc))/um(cc)
! pitch angle in b_unit reference frame
end do
!$OMP END SIMD
! write(6,'("vm: ",E17.10)') vm
! write(6,'("xi: ",E17.10)') xi
call unitVectors_p(b_unit_X,b_unit_Y,b_unit_Z,b1_X,b1_Y,b1_Z, &
b2_X,b2_Y,b2_Z,b3_X,b3_Y,b3_Z)
! b1=b_unit, (b1,b2,b3) is right-handed
!$OMP SIMD
! !$OMP& aligned(phi,U_X,U_Y,U_Z,b3_X,b3_Y,b3_Z,b2_X,b2_Y,b2_Z)
do cc=1_idef,8_idef
phi(cc) = atan2((U_X(cc)*b3_X(cc)+U_Y(cc)*b3_Y(cc)+ &
U_Z(cc)*b3_Z(cc)), &
(U_X(cc)*b2_X(cc)+U_Y(cc)*b2_Y(cc)+U_Z(cc)*b2_Z(cc)))
! azimuthal angle in b_unit refernce frame
end do
!$OMP END SIMD
! write(6,'("phi: ",E17.10)') phi
!$OMP SIMD
! !$OMP& aligned(rnd1,dW,CAL,dCAL,CFL,CBL,vm,ne,Te,Zeff,dpm, &
! !$OMP& flag,dxi,xi,pm,dphi,um,Ub_X,Ub_Y,Ub_Z,U_X,U_Y,U_Z, &
! !$OMP& b1_X,b1_Y,b1_Z,b2_X,b2_Y,b2_Z,b3_X,b3_Y,b3_Z)
do cc=1_idef,8_idef
#ifdef PARALLEL_RANDOM
! uses C library to generate normal_distribution random variables,
! preserving parallelization where Fortran random number generator
! does not
rnd1(cc,1) = get_random()
rnd1(cc,2) = get_random()
rnd1(cc,3) = get_random()
#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))
dW(cc,3) = SQRT(3*dt)*(-1+2*rnd1(cc,3))
! 3D Weiner process
CAL(cc) = CA_SD(vm(cc),ne(cc),Te(cc))
dCAL(cc)= dCA_SD(vm(cc),me,ne(cc),Te(cc))
CFL(cc) = CF_SD(params,vm(cc),ne(cc),Te(cc))
CBL(cc) = (CB_ee_SD(vm(cc),ne(cc),Te(cc),Zeff(cc))+ &
CB_ei_SD(params,vm(cc),ne(cc),Te(cc),Zeff(cc)))
dpm(cc)=REAL(flag(cc))*((-CFL(cc)+dCAL(cc))*dt+ &
sqrt(2.0_rp*CAL(cc))*dW(cc,1))
dxi(cc)=REAL(flag(cc))*(-2*xi(cc)*CBL(cc)/(pm(cc)*pm(cc))*dt- &
sqrt(2.0_rp*CBL(cc)*(1-xi(cc)*xi(cc)))/pm(cc)*dW(cc,2))
dphi(cc)=REAL(flag(cc))*(sqrt(2*CBL(cc))/(pm(cc)* &
sqrt(1-xi(cc)*xi(cc)))*dW(cc,3))
pm(cc)=pm(cc)+dpm(cc)
xi(cc)=xi(cc)+dxi(cc)
phi(cc)=phi(cc)+dphi(cc)
! if (pm(cc)<0) pm(cc)=-pm(cc)
! Keep xi between [-1,1]
if (xi(cc)>1) then
xi(cc)=1-mod(xi(cc),1._rp)
else if (xi(cc)<-1) then
xi(cc)=-1-mod(xi(cc),-1._rp)
endif
! Keep phi between [0,pi]
! if (phi(cc)>C_PI) then
! phi(cc)=C_PI-mod(phi(cc),C_PI)
! else if (phi(cc)<0) then
! phi(cc)=mod(-phi(cc),C_PI)
! endif
um(cc)=pm(cc)/me
Ub_X(cc)=um(cc)*xi(cc)
Ub_Y(cc)=um(cc)*sqrt(1-xi(cc)*xi(cc))*cos(phi(cc))
Ub_Z(cc)=um(cc)*sqrt(1-xi(cc)*xi(cc))*sin(phi(cc))
U_X(cc) = Ub_X(cc)*b1_X(cc)+Ub_Y(cc)*b2_X(cc)+Ub_Z(cc)*b3_X(cc)
U_Y(cc) = Ub_X(cc)*b1_Y(cc)+Ub_Y(cc)*b2_Y(cc)+Ub_Z(cc)*b3_Y(cc)
U_Z(cc) = Ub_X(cc)*b1_Z(cc)+Ub_Y(cc)*b2_Z(cc)+Ub_Z(cc)*b3_Z(cc)
end do
!$OMP END SIMD
! if (tt .EQ. 1_ip) then
! write(6,'("CA: ",E17.10)') CAL(1)
! write(6,'("dCA: ",E17.10)') dCAL(1)
! write(6,'("CF ",E17.10)') CFL(1)
! write(6,'("CB: ",E17.10)') CBL(1)
! end if
do cc=1_idef,8_idef
if (pm(cc).lt.0) then
write(6,'("Momentum less than zero")')
stop
end if
end do
end if
end subroutine include_CoulombCollisions_FO_p