For each primary RE, calculating probability to generate a secondary RE
Generating 1D and 2D ranges for secondary RE distribution
Saving maximum secondary RE source for use in rejection-acceptance sampling algorithm Trapezoidal integration of secondary RE source to find probabilty
If secondary RE generated, begin acceptance-rejection sampling algorithm
Write secondary RE degrees of freedom to particle derived type
Write changes to primary RE degrees of freedom to temporary arrays for passing back out to particle derived type
| Type | Intent | Optional | Attributes | Name | ||
|---|---|---|---|---|---|---|
| type(SPECIES), | intent(inout) | :: | spp | |||
| type(KORC_PARAMS), | intent(in) | :: | params | |||
| integer, | intent(in) | :: | achunk | |||
| type(FIELDS), | intent(in) | :: | F | |||
| real(kind=rp), | intent(in), | DIMENSION(achunk) | :: | Y_R | ||
| real(kind=rp), | intent(in), | DIMENSION(achunk) | :: | Y_PHI | ||
| real(kind=rp), | intent(in), | DIMENSION(achunk) | :: | Y_Z | ||
| real(kind=rp), | intent(inout), | DIMENSION(achunk) | :: | pm | ||
| real(kind=rp), | intent(inout), | DIMENSION(achunk) | :: | xi | ||
| real(kind=rp), | intent(in), | DIMENSION(achunk) | :: | ne | ||
| real(kind=rp), | intent(in), | DIMENSION(achunk) | :: | netot | ||
| real(kind=rp), | intent(in), | DIMENSION(achunk) | :: | Te | ||
| real(kind=rp), | intent(in), | DIMENSION(achunk) | :: | Bmag | ||
| real(kind=rp), | intent(in), | DIMENSION(achunk) | :: | E_PHI | ||
| real(kind=rp), | intent(in) | :: | me | |||
| integer(kind=is), | intent(in), | DIMENSION(achunk) | :: | flagCol | ||
| integer(kind=is), | intent(in), | DIMENSION(achunk) | :: | flagCon |
subroutine large_angle_source(spp,params,achunk,F,Y_R,Y_PHI,Y_Z, &
pm,xi,ne,netot,Te,Bmag,E_PHI,me,flagCol,flagCon)
TYPE(SPECIES), INTENT(INOUT) :: spp
TYPE(KORC_PARAMS), INTENT(IN) :: params
TYPE(FIELDS), INTENT(IN) :: F
INTEGER, INTENT(IN) :: achunk
REAL(rp), INTENT(INOUT), DIMENSION(achunk) :: pm,xi
REAL(rp), INTENT(IN), DIMENSION(achunk) :: Y_R,Y_PHI,Y_Z
REAL(rp), INTENT(IN), DIMENSION(achunk) :: ne,netot,Te
REAL(rp), INTENT(IN), DIMENSION(achunk) :: Bmag,E_PHI
INTEGER(is), INTENT(IN), DIMENSION(achunk) :: flagCol,flagCon
REAL(rp), INTENT(IN) :: me
REAL(rp), DIMENSION(achunk) :: gam,prob0,prob1
REAL(rp) :: gam_min,p_min,gammax,dt,gamsecmax,psecmax,ptrial
REAL(rp) :: gamtrial,cosgam1,xirad,xip,xim,xitrial,sinsq1,cossq1,pitchprob1
REAL(rp) :: dsigdgam1,S_LAmax,S_LA1,tmppm,gamvth,vmin,E_C,p_c,gam_c,pRE
INTEGER :: ngam1,neta1
INTEGER :: ii,jj,cc,seciter
REAL(rp), DIMENSION(cparams_ss%ngrid1) :: gam1,pm1,tmpgam1,tmpcosgam,tmpdsigdgam,tmpsecthreshgam,probtmp,intpitchprob
REAL(rp), DIMENSION(cparams_ss%ngrid1-1) :: dpm1
REAL(rp), DIMENSION(cparams_ss%ngrid1) :: eta1,tmpsinsq,tmpcossq
REAL(rp), DIMENSION(cparams_ss%ngrid1-1) :: deta1
REAL(rp), DIMENSION(cparams_ss%ngrid1,cparams_ss%ngrid1) :: cosgam,sinsq,cossq,tmpcossq1,pitchprob,dsigdgam
REAL(rp), DIMENSION(cparams_ss%ngrid1,cparams_ss%ngrid1) :: secthreshgam,pm11,gam11,eta11,S_LA,pitchrad
LOGICAL :: accepted
dt=cparams_ss%coll_per_dump_dt*params%cpp%time
ngam1=cparams_ss%ngrid1
neta1=cparams_ss%ngrid1
!$OMP SIMD
do cc=1_idef,achunk
gam(cc) = sqrt(1+pm(cc)*pm(cc))
#ifdef PARALLEL_RANDOM
prob0(cc) = get_random()
#else
call RANDOM_NUMBER(prob0)
#endif
end do
!$OMP END SIMD
vmin=1/sqrt(1+1/(cparams_ss%p_min*cparams_ss%pmin_scale)**2)
!! For each primary RE, calculating probability to generate a secondary RE
do cc=1_idef,achunk
E_C=Gammacee(vmin,netot(cc),Te(cc))
!write(6,*) 'E',E_PHI*params%cpp%Eo
!write(6,*) 'E_C',E_C*params%cpp%Eo
!write(6,*) 'E_c,min',cparams_ms%Ec_min*params%cpp%Eo
!write(6,*) 'ne',ne(cc)*params%cpp%density
!write(6,*) 'netot',netot(cc)*params%cpp%density
!write(6,*) 'Te',Te(cc)*params%cpp%temperature
!write(6,*) 'Clog',CLogee_wu(params,ne(cc)*params%cpp%density,Te(cc)*params%cpp%temperature)
if (E_C.gt.abs(E_PHI(cc))) cycle
p_c=cparams_ss%pmin_scale/sqrt(abs(E_PHI(cc))/E_C-1)
gam_c=sqrt(1+p_c**2)
if(cparams_ss%min_secRE.eq.'THERM') then
gam_min=(gam_c+1)/2
p_min=sqrt(gam_min**2-1)
else
gam_min=gam_c
p_min=p_c
end if
gammax=(gam(cc)+1._rp)/2._rp
if (gam_min.eq.1._rp) then
write(6,*) 'R',Y_R(cc)*params%cpp%length,'Z',Y_Z(cc)*params%cpp%length
write(6,*) 'vmin',vmin,'netot',netot(cc)*params%cpp%density,'Te',Te(cc)*params%cpp%temperature/C_E
write(6,*) 'E',E_PHI*params%cpp%Eo,'E_c',E_C*params%cpp%Eo
write(6,*) 'p_c',p_c,'gam_c',gam_c
write(6,*) 'p_min',p_min,'gam_min',gam_min
!write(6,*) 'LAC_gam_resolution: ',TRIM(cparams_ss%LAC_gam_resolution)
!write(6,*) 'gam_min,gammax',gam_min,gammax
end if
!! Generating 1D and 2D ranges for secondary RE distribution
if (TRIM(cparams_ss%LAC_gam_resolution).eq.'LIN') then
do ii=1,ngam1
gam1(ii)=gam_min+(gammax-gam_min)* &
REAL(ii-1)/REAL(ngam1-1)
end do
elseif (TRIM(cparams_ss%LAC_gam_resolution).eq.'EXP') then
do ii=1,ngam1
tmpgam1(ii)=log10(gam_min)+ &
(log10(gammax)-log10(gam_min))* &
REAL(ii-1)/REAL(ngam1-1)
end do
gam1=10**(tmpgam1)
elseif (TRIM(cparams_ss%LAC_gam_resolution).eq.'2EXP') then
do ii=1,ngam1
tmpgam1(ii)=log10(log10(gam_min))+ &
(log10(log10(gammax))-log10(log10(gam_min)))* &
REAL(ii-1)/REAL(ngam1-1)
end do
gam1=10**(10**(tmpgam1))
end if
!write(6,*) 'tmpgam1',tmpgam1
!write(6,*) 'gam1',gam1
pm1=sqrt(gam1**2-1)
do ii=1,ngam1-1
dpm1(ii)=pm1(ii+1)-pm1(ii)
end do
do ii=1,neta1
eta1(ii)=C_PI*(ii-1)/(neta1-1)
end do
!write(6,*) 'eta1',eta1
do ii=1,neta1-1
deta1(ii)=eta1(ii+1)-eta1(ii)
end do
do ii=1,neta1
pm11(:,ii)=pm1
gam11(:,ii)=gam1
end do
do ii=1,ngam1
eta11(ii,:)=eta1
end do
tmpcosgam=sqrt(((gam(cc)+1)*(gam1-1))/((gam(cc)-1)*(gam1+1)))
tmpdsigdgam=2*C_PI*C_RE**2/(gam(cc)**2-1)* &
(((gam(cc)-1)**2*gam(cc)**2)/((gam1-1)**2*(gam(cc)-gam1)**2)- &
(2*gam(cc)**2+2*gam(cc)-1)/((gam1-1)*(gam(cc)-gam1))+1)
tmpsecthreshgam=1._rp
where(gam1.gt.(gam(cc)+1)/2._rp) tmpsecthreshgam=0._rp
!write(6,*) 'tmpcosgam',tmpcosgam
!write(6,*) 'tmpdsigdgam',tmpdsigdgam
!write(6,*) 'tmpsecthreshgam',tmpsecthreshgam
do ii=1,neta1
cosgam(:,ii)=tmpcosgam
dsigdgam(:,ii)=tmpdsigdgam
secthreshgam(:,ii)=tmpsecthreshgam
end do
!if (cc.eq.1) then
!write(6,*) cosgam
!end if
tmpsinsq=(1-xi(cc)**2)*sin(eta1)**2
tmpcossq=xi(cc)*cos(eta1)
do ii=1,ngam1
sinsq(ii,:)=tmpsinsq
tmpcossq1(ii,:)=tmpcossq
end do
!if (cc.eq.1) then
!write(6,*) sinsq
!end if
cossq=(cosgam-tmpcossq1)**2
pitchrad=sinsq-cossq
where(pitchrad.lt.0) pitchrad=tiny(0._rp)
!if (cc.eq.1) then
!write(6,*) cossq
!write(6,*) 'sinsq-cossq',sqrt(sinsq-cossq)
!write(6,*) 'pitchrad',pitchrad
!end if
pitchprob=1/(C_PI*sqrt(pitchrad))
where(pitchprob.eq.1/(C_PI*sqrt(tiny(0._rp)))) pitchprob=0._rp
!if (cc.eq.1) then
!write(6,*) 'pitchprob',pitchprob
!end if
S_LA=netot(cc)*params%cpp%density*C_C/(2*C_PI)* &
(pm11/gam11)*(pm(cc)/gam(cc))* &
pitchprob*dsigdgam*secthreshgam
!! Saving maximum secondary RE source for use in rejection-acceptance
!! sampling algorithm
S_LAmax=maxval(S_LA)
S_LA=S_LA*sin(eta11)
!! Trapezoidal integration of secondary RE source to find probabilty
do ii=1,ngam1
probtmp(ii)=S_LA(ii,1)*deta1(1)/2+S_LA(ii,neta1)*deta1(neta1-1)/2
!intpitchprob(ii)=pitchprob(ii,1)*sin(eta1(1))*deta1(1)/2+ &
! pitchprob(ii,neta1)*sin(eta1(neta1))*deta1(neta1-1)/2
do jj=2,neta1-1
probtmp(ii)=probtmp(ii)+S_LA(ii,jj)*(deta1(jj)+deta1(jj-1))/2
! intpitchprob(ii)=intpitchprob(ii)+ &
! pitchprob(ii,jj)*sin(eta1(jj))*(deta1(jj)+deta1(jj-1))/2
end do
end do
prob1(cc)=probtmp(1)*dpm1(1)/2+probtmp(ngam1)*dpm1(ngam1-1)/2
do jj=2,ngam1-1
prob1(cc)=prob1(cc)+probtmp(jj)*(dpm1(jj)+dpm1(jj-1))/2
end do
!write(6,*) 'prob1pre',prob1(cc),'flagCol',flagCol(cc),'flagCon',flagCon(cc),'dt',dt
!write(6,*) 'intpitchprob',intpitchprob
prob1(cc)=prob1(cc)*dt*2*C_PI*REAL(flagCol(cc))*REAL(flagCon(cc))
if (ISNAN(prob1(cc))) then
write(6,*) 'NaN probability from secondary RE source'
write(6,*) 'p,xi',pm(cc),xi(cc)
write(6,*) 'gam_min,gammax',gam_min,gammax
write(6,*) 'E',E_PHI(cc)*params%cpp%Eo
write(6,*) 'E_C',E_C*params%cpp%Eo
!write(6,*) 'pitchprob',pitchprob
!write(6,*) 'S_LA',S_LA
call korc_abort(24)
end if
if (prob1(cc).gt.1._rp) then
write(6,*) 'Multiple secondary REs generated in a collisional time step'
write(6,*) 'p,xi',pm(cc),xi(cc)
write(6,*) 'gam_min,gammax',gam_min,gammax
write(6,*) 'E',E_PHI(cc)*params%cpp%Eo
write(6,*) 'E_C',E_C*params%cpp%Eo
call korc_abort(24)
end if
!write(6,*) 'gam',gam(cc),'xi',xi(cc)
!write(6,*) 'prob1',prob1(cc),'prob0',prob0(cc)
if (prob1(cc).gt.prob0(cc)) then
!! If secondary RE generated, begin acceptance-rejection sampling
!! algorithm
!write(6,*) 'secondary RE from ',prob1,prob0
accepted=.false.
gamsecmax=(gam(cc)+1)/2
psecmax=sqrt(gamsecmax**2-1)
seciter=0
do while (.not.accepted)
ptrial=p_min+(psecmax-p_min)*get_random()
gamtrial=sqrt(1+ptrial*ptrial)
cosgam1=sqrt(((gam(cc)+1)*(gamtrial-1))/((gam(cc)-1)*(gamtrial+1)))
!xirad=sqrt((cosgam1*xi(cc))**2-(xi(cc)**2+cosgam1**2-1))
!if (isnan(xirad)) then
! write(6,*) 'Sample not in allowable region of phase space'
! call korc_abort(24)
!end if
!xip=cosgam1*xi(cc)+xirad
!xim=cosgam1*xi(cc)-xirad
!xitrial=xim+(xip-xim)*get_random()
xitrial=-1+2*get_random()
sinsq1=(1-xi(cc)*xi(cc))*(1-xitrial*xitrial)
cossq1=(cosgam1-xi(cc)*xitrial)**2
pitchprob1=1/(C_PI*sqrt(sinsq1-cossq1))
!if (isnan(pitchprob1)) cycle
dsigdgam1=2*C_PI*C_RE**2/(gam(cc)**2-1)* &
(((gam(cc)-1)**2*gam(cc)**2)/ &
((gamtrial-1)**2*(gam(cc)-gamtrial)**2)- &
(2*gam(cc)**2+2*gam(cc)-1)/ &
((gamtrial-1)*(gam(cc)-gamtrial))+1)
S_LA1=netot(cc)*params%cpp%density*C_C/(2*C_PI)* &
(ptrial/gamtrial)*(pm(cc)/gam(cc))* &
pitchprob1*dsigdgam1
if (S_LA1/S_LAmax.gt.get_random()) accepted=.true.
seciter=seciter+1
!if (mod(seciter,100).eq.0) then
! write(6,*) 'iteration',seciter
!end if
end do
!! Write secondary RE degrees of freedom to particle derived type
!$OMP ATOMIC UPDATE
spp%pRE=spp%pRE+1
!$OMP ATOMIC WRITE
spp%vars%flagRE(spp%pRE)=1
!$OMP ATOMIC WRITE
spp%vars%flagCon(spp%pRE)=1
!$OMP ATOMIC WRITE
spp%vars%flagCol(spp%pRE)=1
!$OMP ATOMIC WRITE
spp%vars%V(spp%pRE,1)=ptrial*xitrial
!$OMP ATOMIC WRITE
spp%vars%V(spp%pRE,2)=ptrial*ptrial*(1-xitrial*xitrial)/ &
(2*me*Bmag(cc))
!$OMP ATOMIC WRITE
spp%vars%Y(spp%pRE,1)=Y_R(cc)
!$OMP ATOMIC WRITE
spp%vars%Y(spp%pRE,2)=Y_PHI(cc)
!$OMP ATOMIC WRITE
spp%vars%Y(spp%pRE,3)=Y_Z(cc)
!! Write changes to primary RE degrees of freedom to temporary
!! arrays for passing back out to particle derived type
#if DBG_CHECK
if (spp%vars%V(spp%pRE,2).lt.0._rp) then
write(6,*) 'mu is negative for secondary'
write(6,*) 'ppll,mu',spp%vars%V(spp%pRE,1),spp%vars%V(spp%pRE,2)
write(6,*) 'pm,xi',ptrial,xitrial
write(6,*) 'xim,xip',xim,xip
write(6,*) 'pmin,psecmax',p_min,psecmax
call korc_abort(24)
end if
# endif
if (cparams_ss%ConserveLA) then
tmppm=pm(cc)
gamvth=1/sqrt(1-2*Te(cc))
gam(cc)=gam(cc)-gamtrial+gamvth
pm(cc)=sqrt(gam(cc)*gam(cc)-1)
xi(cc)=(tmppm*xi(cc)-ptrial*xitrial)/pm(cc)
end if
!$OMP ATOMIC READ
pRE=spp%pRE
if (pRE.eq.spp%ppp) then
write(6,*) 'All REs allocated on proc',params%mpi_params%rank
call korc_abort(24)
end if
end if
end do
end subroutine large_angle_source