Calls get_fields in korc_fields.
| Type | Intent | Optional | Attributes | Name | ||
|---|---|---|---|---|---|---|
| type(KORC_PARAMS), | intent(in) | :: | params | Core KORC simulation parameters. |
||
| type(FIELDS), | intent(in) | :: | F | An instance of the KORC derived type FIELDS. |
||
| 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. |
|
| logical, | intent(in) | :: | output | |||
| logical, | intent(in) | :: | step |
subroutine FO_init(params,F,spp,output,step)
TYPE(KORC_PARAMS), INTENT(IN) :: params
!! Core KORC simulation parameters.
TYPE(FIELDS), INTENT(IN) :: F
!! An instance of the KORC derived type FIELDS.
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) :: dt
!! Time step used in the leapfrog step (\(\Delta t\)).
REAL(rp) :: Prad
!! Total radiated power of each particle.
REAL(rp) :: Bmag1
!! Magnitude of the magnetic field seen by each particle .
REAL(rp) :: v
!! Speed of each particle.
REAL(rp) :: vpar
!! Parallel velocity \(v_\parallel = \mathbf{v}\cdot \hat{b}\).
REAL(rp) :: vperp
!! Perpendicular velocity \(v_\parallel = |\mathbf{v} - (\mathbf{v}\cdot
!! \hat{b})\hat{b}|\).
REAL(rp) :: tmp
!! Temporary variable used for various computations.
REAL(rp) :: a
!! This variable is used to simplify notation in the code, and
!! is given by \(a=q\Delta t/m\),
REAL(rp), DIMENSION(3) :: Frad
!! Synchrotron radiation reaction force of each particle.
REAL(rp), DIMENSION(3) :: vec
!! Auxiliary vector used in various computations.
REAL(rp), DIMENSION(3) :: b_unit
!! Unitary vector pointing along the local magnetic field \(\hat{b}\).
INTEGER :: ii
!! Species iterator.
INTEGER :: pp
!! Particles iterator.
INTEGER :: cc
!! Chunk iterator.
LOGICAL,intent(in) :: output
LOGICAL,intent(in) :: step
REAL(rp),DIMENSION(8) :: X_X,X_Y,X_Z
REAL(rp),DIMENSION(8) :: Y_R,Y_PHI,Y_Z
REAL(rp),DIMENSION(8) :: B_X,B_Y,B_Z
REAL(rp),DIMENSION(8) :: E_X,E_Y,E_Z
REAL(rp),DIMENSION(8) :: PSIp
INTEGER(is) ,DIMENSION(8) :: flag_cache
do ii = 1_idef,params%num_species
if(output) then
!$OMP PARALLEL DO default(none) &
!$OMP& shared(params,ii,spp,F) &
!$OMP& PRIVATE(pp,cc,X_X,X_Y,X_Z,B_X,B_Y,B_Z, &
!$OMP& E_X,E_Y,E_Z,Y_R,Y_PHI,Y_Z,flag_cache,PSIp)
do pp=1_idef,spp(ii)%ppp,8
!$OMP SIMD
do cc=1_idef,8_idef
X_X(cc)=spp(ii)%vars%X(pp-1+cc,1)
X_Y(cc)=spp(ii)%vars%X(pp-1+cc,2)
X_Z(cc)=spp(ii)%vars%X(pp-1+cc,3)
flag_cache(cc)=spp(ii)%vars%flag(pp-1+cc)
end do
!$OMP END SIMD
call cart_to_cyl_p(X_X,X_Y,X_Z,Y_R,Y_PHI,Y_Z)
if (params%orbit_model(3:5).eq.'new') then
call interp_FOfields_p(F,Y_R,Y_PHI,Y_Z,B_X,B_Y,B_Z, &
E_X,E_Y,E_Z,PSIp,flag_cache)
else if (params%orbit_model(3:5).eq.'old') then
call interp_FOfields1_p(F,Y_R,Y_PHI,Y_Z,B_X,B_Y,B_Z, &
E_X,E_Y,E_Z,PSIp,flag_cache)
end if
!$OMP SIMD
do cc=1_idef,8_idef
spp(ii)%vars%B(pp-1+cc,1) = B_X(cc)
spp(ii)%vars%B(pp-1+cc,2) = B_Y(cc)
spp(ii)%vars%B(pp-1+cc,3) = B_Z(cc)
spp(ii)%vars%E(pp-1+cc,1) = E_X(cc)
spp(ii)%vars%E(pp-1+cc,2) = E_Y(cc)
spp(ii)%vars%E(pp-1+cc,3) = E_Z(cc)
spp(ii)%vars%PSI_P(pp-1+cc) = PSIp(cc)
end do
!$OMP END SIMD
end do
!$OMP END PARALLEL DO
!! Calls [[get_fields]] in [[korc_fields]].
! Interpolates fields at local particles' position and keeps in
! spp%vars. Fields in (R,\(\phi\),Z) coordinates.
! write(6,'("korc_ppusher")')
! write(6,'("B_X: ",E17.10)') spp(ii)%vars%B(:,1)
! write(6,'("B_Z: ",E17.10)') spp(ii)%vars%B(:,2)
! write(6,'("B_Y: ",E17.10)') spp(ii)%vars%B(:,3)
!$OMP PARALLEL DO DEFAULT(none) SHARED(ii,spp) &
!$OMP& FIRSTPRIVATE(E0) &
!$OMP& PRIVATE(pp,b_unit,Bmag1,vpar,v,vperp,vec,tmp)
do pp=1_idef,spp(ii)%ppp
Bmag1 = SQRT(DOT_PRODUCT(spp(ii)%vars%B(pp,:), &
spp(ii)%vars%B(pp,:)))
! Parallel unit vector
b_unit = spp(ii)%vars%B(pp,:)/Bmag1
v = SQRT(DOT_PRODUCT(spp(ii)%vars%V(pp,:),spp(ii)%vars%V(pp,:)))
if (v.GT.korc_zero) then
! Parallel and perpendicular components of velocity
vpar = DOT_PRODUCT(spp(ii)%vars%V(pp,:), b_unit)
vperp = DOT_PRODUCT(spp(ii)%vars%V(pp,:), &
spp(ii)%vars%V(pp,:)) &
- vpar**2
if ( vperp .GE. korc_zero ) then
vperp = SQRT( vperp )
else
vperp = 0.0_rp
end if
! Pitch angle
spp(ii)%vars%eta(pp) = 180.0_rp*MODULO(ATAN2(vperp,vpar), &
2.0_rp*C_PI)/C_PI
! Magnetic moment
spp(ii)%vars%mu(pp) = 0.5_rp*spp(ii)%m* &
spp(ii)%vars%g(pp)**2*vperp**2/Bmag1
! See Northrop's book (The adiabatic motion of charged
! particles)
! Radiated power
tmp = spp(ii)%q**4/(6.0_rp*C_PI*E0*spp(ii)%m**2)
vec = spp(ii)%vars%E(pp,:) + cross(spp(ii)%vars%V(pp,:), &
spp(ii)%vars%B(pp,:))
spp(ii)%vars%Prad(pp) = tmp*( DOT_PRODUCT(spp(ii)% &
vars%E(pp,:), &
spp(ii)%vars%E(pp,:)) + &
DOT_PRODUCT(cross(spp(ii)%vars%V(pp,:), &
spp(ii)%vars%B(pp,:)),spp(ii)%vars%E(pp,:))+ &
spp(ii)%vars%g(pp)**2* &
(DOT_PRODUCT(spp(ii)%vars%E(pp,:), &
spp(ii)%vars%V(pp,:))**2 - DOT_PRODUCT(vec,vec)) )
! Input power due to electric field
spp(ii)%vars%Pin(pp) = spp(ii)%q*DOT_PRODUCT( &
spp(ii)%vars%E(pp,:),spp(ii)%vars%V(pp,:))
else
spp(ii)%vars%eta(pp) = 0.0_rp
spp(ii)%vars%mu(pp) = 0.0_rp
spp(ii)%vars%Prad(pp) = 0.0_rp
spp(ii)%vars%Pin(pp) = 0.0_rp
end if
end do ! loop over particles on an mpi process
!$OMP END PARALLEL DO
end if !(if output)
if(step.and.(.not.params%FokPlan)) then
dt=0.5_rp*params%dt
!$OMP PARALLEL DO FIRSTPRIVATE(dt) PRIVATE(pp,cc) &
!$OMP& SHARED(ii,spp,params)
do pp=1_idef,spp(ii)%ppp,8
!$OMP SIMD
do cc=1_idef,8
spp(ii)%vars%X(pp-1+cc,1) = spp(ii)%vars%X(pp-1+cc,1) + &
dt*spp(ii)%vars%V(pp-1+cc,1)
spp(ii)%vars%X(pp-1+cc,2) = spp(ii)%vars%X(pp-1+cc,2) + &
dt*spp(ii)%vars%V(pp-1+cc,2)
spp(ii)%vars%X(pp-1+cc,3) = spp(ii)%vars%X(pp-1+cc,3) + &
dt*spp(ii)%vars%V(pp-1+cc,3)
end do
!$OMP END SIMD
end do
!$OMP END PARALLEL DO
end if !(if step)
end do ! over species
end subroutine FO_init