Subroutine to advance GC variables
Comment this section further with evolution equations, numerical methods, and descriptions of both.
| Type | Intent | Optional | Attributes | Name | ||
|---|---|---|---|---|---|---|
| type(PARTICLES), | intent(inout) | :: | vars | |||
| integer, | intent(in) | :: | pp | |||
| integer(kind=ip), | intent(in) | :: | tt | time iterator. |
||
| type(KORC_PARAMS), | intent(inout) | :: | params | Core KORC simulation parameters. |
||
| real(kind=rp), | intent(inout), | DIMENSION(params%pchunk) | :: | Y_R | ||
| real(kind=rp), | intent(inout), | DIMENSION(params%pchunk) | :: | Y_PHI | ||
| real(kind=rp), | intent(inout), | DIMENSION(params%pchunk) | :: | Y_Z | ||
| real(kind=rp), | intent(inout), | DIMENSION(params%pchunk) | :: | V_PLL | ||
| real(kind=rp), | intent(inout), | DIMENSION(params%pchunk) | :: | V_MU | ||
| integer(kind=is), | intent(inout), | dimension(params%pchunk) | :: | flagCon | ||
| integer(kind=is), | intent(inout), | dimension(params%pchunk) | :: | flagCol | ||
| real(kind=rp), | intent(in) | :: | q_cache | |||
| real(kind=rp), | intent(in) | :: | m_cache | |||
| real(kind=rp), | intent(out), | DIMENSION(params%pchunk) | :: | B_R | ||
| real(kind=rp), | intent(out), | DIMENSION(params%pchunk) | :: | B_PHI | ||
| real(kind=rp), | intent(out), | DIMENSION(params%pchunk) | :: | B_Z | ||
| type(FIELDS), | intent(in) | :: | F | An instance of the KORC derived type PROFILES. |
||
| type(PROFILES), | intent(in) | :: | P | |||
| real(kind=rp), | intent(out), | DIMENSION(params%pchunk) | :: | PSIp | ||
| real(kind=rp), | intent(out), | DIMENSION(params%pchunk) | :: | E_PHI |
subroutine advance_GCeqn_vars(vars,pp,tt,params,Y_R,Y_PHI,Y_Z,V_PLL,V_MU, &
flagCon,flagCol,q_cache,m_cache,B_R,B_PHI,B_Z,F,P,PSIp,E_PHI)
!! @note Subroutine to advance GC variables \(({\bf X},p_\parallel)\)
!! @endnote
!! Comment this section further with evolution equations, numerical
!! methods, and descriptions of both.
TYPE(KORC_PARAMS), INTENT(INOUT) :: params
!! Core KORC simulation parameters.
TYPE(PARTICLES), INTENT(INOUT) :: vars
TYPE(PROFILES), INTENT(IN) :: P
TYPE(FIELDS), INTENT(IN) :: F
!! An instance of the KORC derived type PROFILES.
REAL(rp) :: dt
!! Time step used in the leapfrog step (\(\Delta t\)).
INTEGER :: cc,pchunk
!! Chunk iterator.
INTEGER(ip),INTENT(IN) :: tt
!! time iterator.
INTEGER,INTENT(IN) :: pp
REAL(rp) :: a1 = 1./5._rp
REAL(rp) :: a21 = 3./40._rp,a22=9./40._rp
REAL(rp) :: a31 = 3./10._rp,a32=-9./10._rp,a33=6./5._rp
REAL(rp) :: a41 = -11./54._rp,a42=5./2._rp,a43=-70./27._rp,a44=35./27._rp
REAL(rp) :: a51 = 1631./55296._rp,a52=175./512._rp,a53=575./13824._rp,a54=44275./110592._rp,a55=253./4096._rp
REAL(rp) :: b1=37./378._rp,b2=0._rp,b3=250./621._rp,b4=125./594._rp,b5=0._rp,b6=512./1771._rp
REAL(rp),DIMENSION(params%pchunk) :: k1_R,k1_PHI,k1_Z,k1_PLL
REAL(rp),DIMENSION(params%pchunk) :: k2_R,k2_PHI,k2_Z,k2_PLL
REAL(rp),DIMENSION(params%pchunk) :: k3_R,k3_PHI,k3_Z,k3_PLL
REAL(rp),DIMENSION(params%pchunk) :: k4_R,k4_PHI,k4_Z,k4_PLL
REAL(rp),DIMENSION(params%pchunk) :: k5_R,k5_PHI,k5_Z,k5_PLL
REAL(rp),DIMENSION(params%pchunk) :: k6_R,k6_PHI,k6_Z,k6_PLL
REAL(rp),DIMENSION(params%pchunk) :: Y0_R,Y0_PHI,Y0_Z
REAL(rp),DIMENSION(params%pchunk),INTENT(INOUT) :: Y_R,Y_PHI,Y_Z
REAL(rp),DIMENSION(params%pchunk),INTENT(OUT) :: B_R,B_PHI,B_Z
REAL(rp),DIMENSION(params%pchunk),INTENT(OUT) :: PSIp
REAL(rp),DIMENSION(params%pchunk) :: curlb_R,curlb_PHI,curlb_Z,gradB_R,gradB_PHI,gradB_Z
REAL(rp),DIMENSION(params%pchunk),INTENT(INOUT) :: V_PLL,V_MU
REAL(rp),DIMENSION(params%pchunk) :: RHS_R,RHS_PHI,RHS_Z,RHS_PLL,V0,E_Z,E_R
REAL(rp),DIMENSION(params%pchunk),INTENT(OUT) :: E_PHI
REAL(rp),DIMENSION(params%pchunk) :: Bmag,ne,Te,Zeff
INTEGER(is),dimension(params%pchunk), intent(inout) :: flagCon,flagCol
REAL(rp) :: ar,R0
REAL(rp),intent(IN) :: q_cache,m_cache
ar=F%AB%a
R0=F%AB%Ro
pchunk=params%pchunk
dt=params%dt
! write(output_unit_write,'("Y_R 0: ",E17.10)') Y_R(1)
! write(output_unit_write,'("Y_PHI 0: ",E17.10)') Y_PHI(1)
! write(output_unit_write,'("Y_Z 0: ",E17.10)') Y_Z(1)
!$OMP SIMD
! !$OMP& aligned(Y0_R,Y0_PHI,Y0_Z,V0,Y_R,Y_PHI,Y_Z,V_PLL)
do cc=1_idef,pchunk
Y0_R(cc)=Y_R(cc)
Y0_PHI(cc)=Y_PHI(cc)
Y0_Z(cc)=Y_Z(cc)
V0(cc)=V_PLL(cc)
end do
!$OMP END SIMD
call analytical_fields_GC_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,PSIp)
if (params%SC_E_add) then
#ifdef PSPLINE
call add_interp_SCE_p(params,F,Y_R,Y_PHI,Y_Z,E_PHI)
#endif
end if
! write(output_unit_write,'("ER:",E17.10)') E_R
! write(output_unit_write,'("EPHI:",E17.10)') E_PHI
! write(output_unit_write,'("EZ:",E17.10)') E_Z
call GCEoM_p(params,RHS_R,RHS_PHI,RHS_Z,RHS_PLL,B_R,B_PHI, &
B_Z,E_R,E_PHI,E_Z,curlb_R,curlb_PHI,curlb_Z,gradB_R, &
gradB_PHI,gradB_Z,V_PLL,V_MU,Y_R,q_cache,m_cache)
!$OMP SIMD
! !$OMP& aligned(Y0_R,Y0_PHI,Y0_Z,V0,Y_R,Y_PHI,Y_Z,V_PLL, &
! !$OMP& RHS_R,RHS_PHI,RHS_Z,RHS_PLL,k1_R,k1_PHI,k1_Z,k1_PLL)
do cc=1_idef,pchunk
k1_R(cc)=dt*RHS_R(cc)
k1_PHI(cc)=dt*RHS_PHI(cc)
k1_Z(cc)=dt*RHS_Z(cc)
k1_PLL(cc)=dt*RHS_PLL(cc)
! vars%RHS(pp-1+cc,1)=RHS_R(cc)
! vars%RHS(pp-1+cc,2)=RHS_PHI(cc)
! vars%RHS(pp-1+cc,3)=RHS_Z(cc)
! vars%RHS(pp-1+cc,4)=RHS_PLL(cc)
Y_R(cc)=Y0_R(cc)+a1*k1_R(cc)
Y_PHI(cc)=Y0_PHI(cc)+a1*k1_PHI(cc)
Y_Z(cc)=Y0_Z(cc)+a1*k1_Z(cc)
V_PLL(cc)=V0(cc) +a1*k1_PLL(cc)
end do
!$OMP END SIMD
! write(output_unit_write,'("Y_R 1: ",E17.10)') Y_R(1)
! write(output_unit_write,'("Y_PHI 1: ",E17.10)') Y_PHI(1)
! write(output_unit_write,'("Y_Z 1: ",E17.10)') Y_Z(1)
! write(output_unit_write,'("k1R: ",E17.10)') k1_R(1)
! write(output_unit_write,'("k1PHI: ",E17.10)') k1_PHI(1)
! write(output_unit_write,'("k1Z: ",E17.10)') k1_Z(1)
! write(output_unit_write,'("k1PLL: ",E17.10)') k1_PLL(1)
call analytical_fields_GC_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,PSIp)
if (params%SC_E_add) then
#ifdef PSPLINE
call add_interp_SCE_p(params,F,Y_R,Y_PHI,Y_Z,E_PHI)
#endif
end if
call GCEoM_p(params,RHS_R,RHS_PHI,RHS_Z,RHS_PLL,B_R,B_PHI, &
B_Z,E_R,E_PHI,E_Z,curlb_R,curlb_PHI,curlb_Z,gradB_R, &
gradB_PHI,gradB_Z,V_PLL,V_MU,Y_R,q_cache,m_cache)
!$OMP SIMD
! !$OMP& aligned(Y0_R,Y0_PHI,Y0_Z,V0,Y_R,Y_PHI,Y_Z,V_PLL, &
! !$OMP& RHS_R,RHS_PHI,RHS_Z,RHS_PLL,k2_R,k2_PHI,k2_Z,k2_PLL)
do cc=1_idef,pchunk
k2_R(cc)=dt*RHS_R(cc)
k2_PHI(cc)=dt*RHS_PHI (cc)
k2_Z(cc)=dt*RHS_Z(cc)
k2_PLL(cc)=dt*RHS_PLL(cc)
Y_R(cc)=Y0_R(cc)+a21*k1_R(cc)+a22*k2_R(cc)
Y_PHI(cc)=Y0_PHI(cc)+a21*k1_PHI(cc)+a22*k2_PHI(cc)
Y_Z(cc)=Y0_Z(cc)+a21*k1_Z(cc)+a22*k2_Z(cc)
V_PLL(cc)=V0(cc) +a21*k1_PLL(cc)+a22*k2_PLL(cc)
end do
!$OMP END SIMD
! write(output_unit_write,'("Y_R 2: ",E17.10)') Y_R(1)
! write(output_unit_write,'("Y_PHI 2: ",E17.10)') Y_PHI(1)
! write(output_unit_write,'("Y_Z 2: ",E17.10)') Y_Z(1)
call analytical_fields_GC_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,PSIp)
if (params%SC_E_add) then
#ifdef PSPLINE
call add_interp_SCE_p(params,F,Y_R,Y_PHI,Y_Z,E_PHI)
#endif
end if
call GCEoM_p(params,RHS_R,RHS_PHI,RHS_Z,RHS_PLL,B_R,B_PHI, &
B_Z,E_R,E_PHI,E_Z,curlb_R,curlb_PHI,curlb_Z,gradB_R, &
gradB_PHI,gradB_Z,V_PLL,V_MU,Y_R,q_cache,m_cache)
!$OMP SIMD
! !$OMP& aligned(Y0_R,Y0_PHI,Y0_Z,V0,Y_R,Y_PHI,Y_Z,V_PLL, &
! !$OMP& RHS_R,RHS_PHI,RHS_Z,RHS_PLL,k3_R,k3_PHI,k3_Z,k3_PLL)
do cc=1_idef,pchunk
k3_R(cc)=dt*RHS_R(cc)
k3_PHI(cc)=dt*RHS_PHI(cc)
k3_Z(cc)=dt*RHS_Z(cc)
k3_PLL(cc)=dt*RHS_PLL(cc)
Y_R(cc)=Y0_R(cc)+a31*k1_R(cc)+a32*k2_R(cc)+a33*k3_R(cc)
Y_PHI(cc)=Y0_PHI(cc)+a31*k1_PHI(cc)+a32*k2_PHI(cc)+ &
a33*k3_PHI(cc)
Y_Z(cc)=Y0_Z(cc)+a31*k1_Z(cc)+a32*k2_Z(cc)+a33*k3_Z(cc)
V_PLL(cc)=V0(cc) +a31*k1_PLL(cc)+a32*k2_PLL(cc)+a33*k3_PLL(cc)
end do
!$OMP END SIMD
! write(output_unit_write,'("Y_R 3: ",E17.10)') Y_R(1)
! write(output_unit_write,'("Y_PHI 3: ",E17.10)') Y_PHI(1)
! write(output_unit_write,'("Y_Z 3: ",E17.10)') Y_Z(1)
call analytical_fields_GC_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,PSIp)
if (params%SC_E_add) then
#ifdef PSPLINE
call add_interp_SCE_p(params,F,Y_R,Y_PHI,Y_Z,E_PHI)
#endif
end if
call GCEoM_p(params,RHS_R,RHS_PHI,RHS_Z,RHS_PLL,B_R,B_PHI, &
B_Z,E_R,E_PHI,E_Z,curlb_R,curlb_PHI,curlb_Z,gradB_R, &
gradB_PHI,gradB_Z,V_PLL,V_MU,Y_R,q_cache,m_cache)
!$OMP SIMD
! !$OMP& aligned(Y0_R,Y0_PHI,Y0_Z,V0,Y_R,Y_PHI,Y_Z,V_PLL, &
! !$OMP& RHS_R,RHS_PHI,RHS_Z,RHS_PLL,k4_R,k4_PHI,k4_Z,k4_PLL)
do cc=1_idef,pchunk
k4_R(cc)=dt*RHS_R(cc)
k4_PHI(cc)=dt*RHS_PHI(cc)
k4_Z(cc)=dt*RHS_Z(cc)
k4_PLL(cc)=dt*RHS_PLL(cc)
Y_R(cc)=Y0_R(cc)+a41*k1_R(cc)+a42*k2_R(cc)+a43*k3_R(cc)+ &
a44*k4_R(cc)
Y_PHI(cc)=Y0_PHI(cc)+a41*k1_PHI(cc)+a42*k2_PHI(cc)+ &
a43*k3_PHI(cc)+a44*k4_PHI(cc)
Y_Z(cc)=Y0_Z(cc)+a41*k1_Z(cc)+a42*k2_Z(cc)+a43*k3_Z(cc)+ &
a44*k4_Z(cc)
V_PLL(cc)=V0(cc) +a41*k1_PLL(cc)+a42*k2_PLL(cc)+ &
a43*k3_PLL(cc)+a44*k4_PLL(cc)
end do
!$OMP END SIMD
! write(output_unit_write,'("Y_R 4: ",E17.10)') Y_R(1)
! write(output_unit_write,'("Y_PHI 4: ",E17.10)') Y_PHI(1)
! write(output_unit_write,'("Y_Z 4: ",E17.10)') Y_Z(1)
call analytical_fields_GC_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,PSIp)
if (params%SC_E_add) then
#ifdef PSPLINE
call add_interp_SCE_p(params,F,Y_R,Y_PHI,Y_Z,E_PHI)
#endif
end if
call GCEoM_p(params,RHS_R,RHS_PHI,RHS_Z,RHS_PLL,B_R,B_PHI, &
B_Z,E_R,E_PHI,E_Z,curlb_R,curlb_PHI,curlb_Z,gradB_R, &
gradB_PHI,gradB_Z,V_PLL,V_MU,Y_R,q_cache,m_cache)
!$OMP SIMD
! !$OMP& aligned(Y0_R,Y0_PHI,Y0_Z,V0,Y_R,Y_PHI,Y_Z,V_PLL, &
! !$OMP& RHS_R,RHS_PHI,RHS_Z,RHS_PLL,k5_R,k5_PHI,k5_Z,k5_PLL)
do cc=1_idef,pchunk
k5_R(cc)=dt*RHS_R(cc)
k5_PHI(cc)=dt*RHS_PHI(cc)
k5_Z(cc)=dt*RHS_Z(cc)
k5_PLL(cc)=dt*RHS_PLL(cc)
Y_R(cc)=Y0_R(cc)+a51*k1_R(cc)+a52*k2_R(cc)+a53*k3_R(cc)+ &
a54*k4_R(cc)+a55*k5_R(cc)
Y_PHI(cc)=Y0_PHI(cc)+a51*k1_PHI(cc)+a52*k2_PHI(cc)+ &
a53*k3_PHI(cc)+a54*k4_PHI(cc)+a55*k5_PHI(cc)
Y_Z(cc)=Y0_Z(cc)+a51*k1_Z(cc)+a52*k2_Z(cc)+a53*k3_Z(cc)+ &
a54*k4_Z(cc)+a55*k5_Z(cc)
V_PLL(cc)=V0(cc) +a51*k1_PLL(cc)+a52*k2_PLL(cc)+ &
a53*k3_PLL(cc)+a54*k4_PLL(cc)+a55*k5_PLL(cc)
end do
!$OMP END SIMD
! write(output_unit_write,'("Y_R 5: ",E17.10)') Y_R(1)
! write(output_unit_write,'("Y_PHI 5: ",E17.10)') Y_PHI(1)
! write(output_unit_write,'("Y_Z 5: ",E17.10)') Y_Z(1)
call analytical_fields_GC_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,PSIp)
if (params%SC_E_add) then
#ifdef PSPLINE
call add_interp_SCE_p(params,F,Y_R,Y_PHI,Y_Z,E_PHI)
#endif
end if
call GCEoM_p(params,RHS_R,RHS_PHI,RHS_Z,RHS_PLL,B_R,B_PHI, &
B_Z,E_R,E_PHI,E_Z,curlb_R,curlb_PHI,curlb_Z,gradB_R, &
gradB_PHI,gradB_Z,V_PLL,V_MU,Y_R,q_cache,m_cache)
!$OMP SIMD
! !$OMP& aligned(Y0_R,Y0_PHI,Y0_Z,V0,Y_R,Y_PHI,Y_Z,V_PLL, &
! !$OMP& RHS_R,RHS_PHI,RHS_Z,RHS_PLL,k6_R,k6_PHI,k6_Z,k6_PLL)
do cc=1_idef,pchunk
k6_R(cc)=dt*RHS_R(cc)
k6_PHI(cc)=dt*RHS_PHI(cc)
k6_Z(cc)=dt*RHS_Z(cc)
k6_PLL(cc)=dt*RHS_PLL(cc)
Y_R(cc)=Y0_R(cc)+b1*k1_R(cc)+b2*k2_R(cc)+ &
b3*k3_R(cc)+b4*k4_R(cc)+b5*k5_R(cc)+b6*k6_R(cc)
Y_PHI(cc)=Y0_PHI(cc)+b1*k1_PHI(cc)+b2*k2_PHI(cc)+ &
b3*k3_PHI(cc)+b4*k4_PHI(cc)+b5*k5_PHI(cc)+b6*k6_PHI(cc)
Y_Z(cc)=Y0_Z(cc)+b1*k1_Z(cc)+b2*k2_Z(cc)+ &
b3*k3_Z(cc)+b4*k4_Z(cc)+b5*k5_Z(cc)+b6*k6_Z(cc)
V_PLL(cc)=V0(cc)+b1*k1_PLL(cc)+b2*k2_PLL(cc)+ &
b3*k3_PLL(cc)+b4*k4_PLL(cc)+b5*k5_PLL(cc)+b6*k6_PLL(cc)
end do
!$OMP END SIMD
! write(output_unit_write,'("Y_R 6: ",E17.10)') Y_R(1)
! write(output_unit_write,'("Y_PHI 6: ",E17.10)') Y_PHI(1)
! write(output_unit_write,'("Y_Z 6: ",E17.10)') Y_Z(1)
call cyl_check_if_confined_p(pchunk,ar,R0,Y_R,Y_Z,flagCon)
!$OMP SIMD
! !$OMP& aligned(Y_R,Y_PHI,Y_Z,V_PLL,Y0_R,Y0_PHI,Y0_Z,V0)
do cc=1_idef,pchunk
if ((flagCon(cc).eq.0_is).or.(flagCol(cc).eq.0_is)) then
Y_R(cc)=Y0_R(cc)
Y_PHI(cc)=Y0_PHI(cc)
Y_Z(cc)=Y0_Z(cc)
V_PLL(cc)=V0(cc)
end if
end do
!$OMP END SIMD
call analytical_fields_GC_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,PSIp)
if (params%SC_E_add) then
#ifdef PSPLINE
call add_interp_SCE_p(params,F,Y_R,Y_PHI,Y_Z,E_PHI)
#endif
end if
end subroutine advance_GCeqn_vars