Type | Intent | Optional | Attributes | Name | ||
---|---|---|---|---|---|---|
type(KORC_PARAMS), | intent(inout) | :: | params | Core KORC simulation parameters. |
||
type(SPECIES), | intent(inout), | DIMENSION(:), ALLOCATABLE | :: | spp | An instance of the derived type SPECIES containing all the parameters and simulation variables of the different species in the simulation. |
|
type(PROFILES), | intent(in) | :: | P | |||
type(FIELDS), | intent(inout) | :: | F |
subroutine adv_GCinterp_psiwE_top(params,spp,P,F)
USE omp_lib
IMPLICIT NONE
TYPE(KORC_PARAMS), INTENT(INOUT) :: params
!! Core KORC simulation parameters.
TYPE(PROFILES), INTENT(IN) :: P
TYPE(FIELDS), INTENT(INOUT) :: F
TYPE(SPECIES), DIMENSION(:), ALLOCATABLE, INTENT(INOUT) :: spp
!! An instance of the derived type SPECIES containing all the parameters
!! and simulation variables of the different species in the simulation.
REAL(rp), DIMENSION(params%pchunk) :: Bmag
REAL(rp),DIMENSION(params%pchunk) :: Y_R,Y_PHI,Y_Z
REAL(rp),DIMENSION(params%pchunk) :: B_R,B_PHI,B_Z
REAL(rp),DIMENSION(params%pchunk) :: E_R,E_PHI,E_Z
REAL(rp),DIMENSION(params%pchunk) :: ne,Te,Zeff
REAL(rp),DIMENSION(params%pchunk) :: V_PLL,V_MU
REAL(rp),DIMENSION(params%pchunk) :: PSIp
REAL(rp),DIMENSION(params%pchunk) :: curlb_R,curlb_PHI,curlb_Z
REAL(rp),DIMENSION(params%pchunk) :: gradB_R,gradB_PHI,gradB_Z
INTEGER(is),DIMENSION(params%pchunk) :: flagCon,flagCol
REAL(rp) :: m_cache,q_cache,B0,EF0,R0,q0,lam,ar
INTEGER :: ii
!! Species iterator.
INTEGER :: pp
!! Particles iterator.
INTEGER :: cc,pchunk,achunk
!! Chunk iterator.
INTEGER(ip) :: tt,ttt
!! time iterator.
INTEGER :: thread_num
do ii = 1_idef,params%num_species
pchunk=params%pchunk
q_cache=spp(ii)%q
m_cache=spp(ii)%m
if (.not.params%LargeCollisions) then
!$OMP PARALLEL DO default(none) &
!$OMP& FIRSTPRIVATE(q_cache,m_cache,pchunk) &
!$OMP& SHARED(params,ii,spp,P,F) &
!$OMP& PRIVATE(pp,tt,Bmag,cc,Y_R,Y_PHI,Y_Z,V_PLL,V_MU,B_R,B_PHI,B_Z, &
!$OMP& flagCon,flagCol,E_PHI,PSIp,curlb_R,curlb_PHI,curlb_Z, &
!$OMP& gradB_R,gradB_PHI,gradB_Z,ne,E_R,E_Z,thread_num)
do pp=1_idef,spp(ii)%ppp,pchunk
thread_num = OMP_GET_THREAD_NUM()
! write(output_unit_write,'("pp: ",I16)') pp
!$OMP SIMD
do cc=1_idef,pchunk
Y_R(cc)=spp(ii)%vars%Y(pp-1+cc,1)
Y_PHI(cc)=spp(ii)%vars%Y(pp-1+cc,2)
Y_Z(cc)=spp(ii)%vars%Y(pp-1+cc,3)
V_PLL(cc)=spp(ii)%vars%V(pp-1+cc,1)
V_MU(cc)=spp(ii)%vars%V(pp-1+cc,2)
PSIp(cc)=spp(ii)%vars%PSI_P(pp-1+cc)
flagCon(cc)=spp(ii)%vars%flagCon(pp-1+cc)
flagCol(cc)=spp(ii)%vars%flagCol(pp-1+cc)
end do
!$OMP END SIMD
if (.not.params%FokPlan) then
do tt=1_ip,params%t_skip
!if (params%t_skip.ge.10) then
! if(mod(tt,params%t_skip/10).eq.0) then
! if((params%mpi_params%rank.eq.0).and. &
! thread_num.eq.0) then
! write(6,'("tt iteration ",I8)') tt
! endif
! end if
!end if
call advance_GCinterp_psiwE_vars(spp(ii)%vars,pchunk,pp,tt, &
params, &
Y_R,Y_PHI,Y_Z,V_PLL,V_MU,q_cache,m_cache,flagCon,flagCol, &
F,P,B_R,B_PHI,B_Z,E_PHI,PSIp,curlb_R,curlb_PHI, &
curlb_Z,gradB_R,gradB_PHI,gradB_Z,ne)
if (params%collisions) then
call include_CoulombCollisions_GC_p(tt,params, &
Y_R,Y_PHI,Y_Z, V_PLL,V_MU,m_cache, &
flagCon,flagCol,F,P,E_PHI,ne,PSIp)
end if
end do !timestep iterator
!$OMP SIMD
do cc=1_idef,pchunk
spp(ii)%vars%Y(pp-1+cc,1)=Y_R(cc)
spp(ii)%vars%Y(pp-1+cc,2)=Y_PHI(cc)
spp(ii)%vars%Y(pp-1+cc,3)=Y_Z(cc)
spp(ii)%vars%V(pp-1+cc,1)=V_PLL(cc)
spp(ii)%vars%V(pp-1+cc,2)=V_MU(cc)
spp(ii)%vars%flagCon(pp-1+cc)=flagCon(cc)
spp(ii)%vars%flagCol(pp-1+cc)=flagCol(cc)
spp(ii)%vars%B(pp-1+cc,1) = B_R(cc)
spp(ii)%vars%B(pp-1+cc,2) = B_PHI(cc)
spp(ii)%vars%B(pp-1+cc,3) = B_Z(cc)
spp(ii)%vars%gradB(pp-1+cc,1) = gradB_R(cc)
spp(ii)%vars%gradB(pp-1+cc,2) = gradB_PHI(cc)
spp(ii)%vars%gradB(pp-1+cc,3) = gradB_Z(cc)
spp(ii)%vars%curlb(pp-1+cc,1) = curlb_R(cc)
spp(ii)%vars%curlb(pp-1+cc,2) = curlb_PHI(cc)
spp(ii)%vars%curlb(pp-1+cc,3) = curlb_Z(cc)
spp(ii)%vars%E(pp-1+cc,2) = E_PHI(cc)
spp(ii)%vars%PSI_P(pp-1+cc) = PSIp(cc)
spp(ii)%vars%ne(pp-1+cc) = ne(cc)
end do
!$OMP END SIMD
else if (params%FokPlan.and.params%collisions) then
do tt=1_ip,params%t_skip
call include_CoulombCollisions_GC_p(tt,params,Y_R,Y_PHI,Y_Z, &
V_PLL,V_MU,m_cache,flagCon,flagCol,F,P,E_PHI,ne,PSIp)
end do
!$OMP SIMD
do cc=1_idef,pchunk
spp(ii)%vars%V(pp-1+cc,1)=V_PLL(cc)
spp(ii)%vars%V(pp-1+cc,2)=V_MU(cc)
spp(ii)%vars%flagCol(pp-1+cc)=flagCol(cc)
spp(ii)%vars%E(pp-1+cc,2) = E_PHI(cc)
spp(ii)%vars%ne(pp-1+cc) = ne(cc)
end do
!$OMP END SIMD
else
do tt=1_ip,params%t_skip
call calculate_GCfields_p(pchunk,F,Y_R,Y_PHI,Y_Z,B_R,B_PHI,B_Z,&
E_R,E_PHI,E_Z,curlb_R,curlb_PHI,curlb_Z, &
gradB_R,gradB_PHI,gradB_Z,flagCon,PSIp)
end do
!$OMP SIMD
do cc=1_idef,pchunk
spp(ii)%vars%B(pp-1+cc,1) = B_R(cc)
spp(ii)%vars%B(pp-1+cc,2) = B_PHI(cc)
spp(ii)%vars%B(pp-1+cc,3) = B_Z(cc)
spp(ii)%vars%E(pp-1+cc,2) = E_PHI(cc)
spp(ii)%vars%PSI_P(pp-1+cc) = PSIp(cc)
end do
!$OMP END SIMD
end if
!$OMP SIMD
do cc=1_idef,pchunk
B_R(cc)=spp(ii)%vars%B(pp-1+cc,1)
B_PHI(cc)=spp(ii)%vars%B(pp-1+cc,2)
B_Z(cc)=spp(ii)%vars%B(pp-1+cc,3)
Bmag(cc)=sqrt(B_R(cc)*B_R(cc)+B_PHI(cc)*B_PHI(cc)+ &
B_Z(cc)*B_Z(cc))
spp(ii)%vars%g(pp-1+cc)=sqrt(1+V_PLL(cc)**2+ &
2*V_MU(cc)*Bmag(cc))
spp(ii)%vars%eta(pp-1+cc) = atan2(sqrt(2*m_cache*Bmag(cc)* &
spp(ii)%vars%V(pp-1+cc,2)),spp(ii)%vars%V(pp-1+cc,1))* &
180.0_rp/C_PI
end do
!$OMP END SIMD
end do !particle chunk iterator
!$OMP END PARALLEL DO
else
do tt=1_ip,params%coll_per_dump
#if DBG_CHECK
! if(params%mpi_params%rank.eq.6) then
! write(6,*) 'before loop load:ppll',spp(ii)%vars%V(1:8,1),'mu',spp(ii)%vars%V(1:8,2)
! write(6,*) 'before loop load:R',spp(ii)%vars%Y(1:8,1),'PHI',spp(ii)%vars%Y(1:8,2),'Z',spp(ii)%vars%Y(1:8,3)
! end if
# endif
!if (modulo(tt,params%coll_per_dump/10).eq.0) &
! write(6,*) 'mpi',params%mpi_params%rank,', Coll step',tt
!$OMP PARALLEL DO default(none) &
!$OMP& FIRSTPRIVATE(q_cache,m_cache,pchunk,tt) &
!$OMP& SHARED(params,ii,spp,P,F) &
!$OMP& PRIVATE(pp,ttt,Bmag,cc,Y_R,Y_PHI,Y_Z,V_PLL,V_MU, &
!$OMP& B_R,B_PHI,B_Z,achunk, &
!$OMP& flagCon,flagCol,E_PHI,PSIp,curlb_R,curlb_PHI,curlb_Z, &
!$OMP& gradB_R,gradB_PHI,gradB_Z,ne,Te,E_R,E_Z,thread_num)
do pp=1_idef,spp(ii)%pRE,pchunk
thread_num = OMP_GET_THREAD_NUM()
if ((spp(ii)%pRE-pp).lt.pchunk) then
achunk=spp(ii)%pRE-pp+1
else
achunk=pchunk
end if
! write(output_unit_write,'("pp: ",I16)') pp
!$OMP SIMD
do cc=1_idef,achunk
Y_R(cc)=spp(ii)%vars%Y(pp-1+cc,1)
Y_PHI(cc)=spp(ii)%vars%Y(pp-1+cc,2)
Y_Z(cc)=spp(ii)%vars%Y(pp-1+cc,3)
V_PLL(cc)=spp(ii)%vars%V(pp-1+cc,1)
V_MU(cc)=spp(ii)%vars%V(pp-1+cc,2)
PSIp(cc)=spp(ii)%vars%PSI_P(pp-1+cc)
flagCon(cc)=spp(ii)%vars%flagCon(pp-1+cc)
flagCol(cc)=spp(ii)%vars%flagCol(pp-1+cc)
end do
!$OMP END SIMD
#if DBG_CHECK
do cc=1_idef,achunk
if (V_MU(cc).le.0._rp) then
write(6,*) 'mu is negative before orbit loop'
write(6,*) 'coll_it',tt
write(6,*) 'V_PLL,V_MU',V_PLL(cc),V_MU(cc)
write(6,*) 'mpi,omp',params%mpi_params%rank,thread_num
write(6,*) 'pp,cc,part,pRE',pp,cc,pp-1+cc,spp%pRE
call korc_abort(25)
endif
end do
#endif
!if (params%t_skip.ge.10) then
! if(mod(tt,params%t_skip/10).eq.0) then
! if((params%mpi_params%rank.eq.0).and. &
! thread_num.eq.0) then
! write(6,'("tt iteration ",I8)') tt
! endif
! end if
!end if
do ttt=1_ip,params%orbits_per_coll
!if (modulo(ttt,params%orbits_per_coll/10).eq.0) &
! write(6,*) 'mpi',params%mpi_params%rank,' , OMP',thread_num,', Orbit step',ttt
#if DBG_CHECK
do cc=1_idef,achunk
if (V_MU(cc).le.0._rp) then
write(6,*) 'mu is negative in orbit loop'
write(6,*) 'coll_it,orb_it',tt,ttt
write(6,*) 'V_PLL,V_MU',V_PLL(cc),V_MU(cc)
write(6,*) 'mpi,omp',params%mpi_params%rank,thread_num
write(6,*) 'pp,cc,part,pRE',pp,cc,pp-1+cc,spp%pRE
call korc_abort(25)
endif
end do
# endif
call advance_GCinterp_psiwE_vars(spp(ii)%vars,achunk, &
pp,tt,params,Y_R,Y_PHI,Y_Z,V_PLL,V_MU, &
q_cache,m_cache,flagCon,flagCol, &
F,P,B_R,B_PHI,B_Z,E_PHI,PSIp,curlb_R,curlb_PHI, &
curlb_Z,gradB_R,gradB_PHI,gradB_Z,ne)
end do
#if DBG_CHECK
do cc=1_idef,achunk
if (V_MU(cc).le.0._rp) then
write(6,*) 'mu is negative after orbit loop'
write(6,*) 'coll_it,orb_it',tt,ttt
write(6,*) 'V_PLL,V_MU',V_PLL(cc),V_MU(cc)
write(6,*) 'mpi,omp',params%mpi_params%rank,thread_num
write(6,*) 'part,pRE',pp-1+cc,spp%pRE
call korc_abort(25)
endif
end do
# endif
call include_CoulombCollisionsLA_GC_p(spp(ii),achunk, &
tt,params,Y_R,Y_PHI,Y_Z,V_PLL,V_MU,m_cache, &
flagCon,flagCol,F,P,E_PHI,ne,Te,PSIp)
#if DBG_CHECK
do cc=1_idef,achunk
if (V_MU(cc).le.0._rp) then
write(6,*) 'mu is negative after coll'
write(6,*) 'coll_it,orb_it',tt,ttt
write(6,*) 'V_PLL,V_MU',V_PLL(cc),V_MU(cc)
write(6,*) 'mpi,omp',params%mpi_params%rank,thread_num
write(6,*) 'part,pRE',pp-1+cc,spp%pRE
call korc_abort(25)
endif
end do
if((params%mpi_params%rank.eq.6).and.(pp.eq.1)) then
write(6,*) 'before loop save:ppll',V_PLL,'mu',V_MU
end if
#endif
!$OMP SIMD
do cc=1_idef,achunk
spp(ii)%vars%Y(pp-1+cc,1)=Y_R(cc)
spp(ii)%vars%Y(pp-1+cc,2)=Y_PHI(cc)
spp(ii)%vars%Y(pp-1+cc,3)=Y_Z(cc)
spp(ii)%vars%V(pp-1+cc,1)=V_PLL(cc)
spp(ii)%vars%V(pp-1+cc,2)=V_MU(cc)
spp(ii)%vars%flagCon(pp-1+cc)=flagCon(cc)
spp(ii)%vars%flagCol(pp-1+cc)=flagCol(cc)
spp(ii)%vars%B(pp-1+cc,1) = B_R(cc)
spp(ii)%vars%B(pp-1+cc,2) = B_PHI(cc)
spp(ii)%vars%B(pp-1+cc,3) = B_Z(cc)
spp(ii)%vars%gradB(pp-1+cc,1) = gradB_R(cc)
spp(ii)%vars%gradB(pp-1+cc,2) = gradB_PHI(cc)
spp(ii)%vars%gradB(pp-1+cc,3) = gradB_Z(cc)
spp(ii)%vars%curlb(pp-1+cc,1) = curlb_R(cc)
spp(ii)%vars%curlb(pp-1+cc,2) = curlb_PHI(cc)
spp(ii)%vars%curlb(pp-1+cc,3) = curlb_Z(cc)
spp(ii)%vars%E(pp-1+cc,2) = E_PHI(cc)
spp(ii)%vars%PSI_P(pp-1+cc) = PSIp(cc)
spp(ii)%vars%ne(pp-1+cc) = ne(cc)
end do
!$OMP END SIMD
#if DBG_CHECK
! if(params%mpi_params%rank.eq.6.and.(pp.eq.1)) then
! write(6,*) 'after loop save 1:ppll',spp(ii)%vars%V(1:8,1),'mu',spp(ii)%vars%V(1:8,2)
! write(6,*) 'after loop save 1:R',spp(ii)%vars%Y(1:8,1),'PHI',spp(ii)%vars%Y(1:8,2),'Z',spp(ii)%vars%Y(1:8,3)
! end if
#endif
end do !particle chunk iterator
!$OMP END PARALLEL DO
#if DBG_CHECK
! if(params%mpi_params%rank.eq.6) then
! write(6,*) 'after loop save 2:ppll',spp(ii)%vars%V(1:8,1),'mu',spp(ii)%vars%V(1:8,2)
! write(6,*) 'after loop save 2:R',spp(ii)%vars%Y(1:8,1),'PHI',spp(ii)%vars%Y(1:8,2),'Z',spp(ii)%vars%Y(1:8,3)
! end if
#endif
end do !timestep iterator
!$OMP PARALLEL DO default(none) &
!$OMP& FIRSTPRIVATE(m_cache,pchunk) &
!$OMP& SHARED(ii,spp) &
!$OMP& PRIVATE(pp,Bmag,cc, &
!$OMP& B_R,B_PHI,B_Z,achunk)
do pp=1_idef,spp(ii)%pRE,pchunk
if ((spp(ii)%pRE-pp).lt.pchunk) then
achunk=spp(ii)%pRE-pp+1
else
achunk=pchunk
end if
!$OMP SIMD
do cc=1_idef,achunk
B_R(cc)=spp(ii)%vars%B(pp-1+cc,1)
B_PHI(cc)=spp(ii)%vars%B(pp-1+cc,2)
B_Z(cc)=spp(ii)%vars%B(pp-1+cc,3)
Bmag(cc)=sqrt(B_R(cc)*B_R(cc)+B_PHI(cc)*B_PHI(cc)+ &
B_Z(cc)*B_Z(cc))
spp(ii)%vars%g(pp-1+cc)=sqrt(1+spp(ii)%vars%V(pp-1+cc,1)**2 &
+2*spp(ii)%vars%V(pp-1+cc,2)*Bmag(cc))
spp(ii)%vars%eta(pp-1+cc) = atan2(sqrt(2*m_cache*Bmag(cc)* &
spp(ii)%vars%V(pp-1+cc,2)),spp(ii)%vars%V(pp-1+cc,1))* &
180.0_rp/C_PI
end do
!$OMP END SIMD
end do !particle chunk iterator
!$OMP END PARALLEL DO
endif
end do !species iterator
end subroutine adv_GCinterp_psiwE_top