Reads in a MC particle fast-ion distribution and puts them in particles
| Type | Intent | Optional | Attributes | Name | ||
|---|---|---|---|---|---|---|
| integer(kind=HID_T), | intent(inout) | :: | fid | HDF5 file ID |
||
| integer, | intent(out) | :: | error | Error code |
subroutine read_mc(fid, error)
!+ Reads in a MC particle fast-ion distribution and puts them in [[libfida:particles]]
integer(HID_T), intent(inout) :: fid
!+ HDF5 file ID
integer, intent(out) :: error
!+ Error code
integer(HSIZE_T), dimension(1) :: dims
integer(Int32) :: i,j,ii,ir,iz,iphi,nphi
real(Float64) :: phi,beam_grid_phi_enter,beam_grid_phi_exit,delta_phi,xp,yp,zp
real(Float64), dimension(3) :: uvw,xyz,ri,vi,e1_xyz,e2_xyz,C_xyz,dum
real(Float64), dimension(:), allocatable :: weight
type(LocalEMFields) :: fields
integer :: cnt,num
logical :: inp,path_valid
character(len=50) :: dist_type_name = ''
if(inputs%verbose.ge.1) then
write(*,'(a)') '---- Fast-ion distribution settings ----'
endif
call h5ltread_dataset_int_scalar_f(fid, "/nparticle", particles%nparticle, error)
!!ALLOCATE SPACE
allocate(particles%fast_ion(particles%nparticle))
allocate(weight(particles%nparticle))
dims(1) = particles%nparticle
!call h5ltread_dataset_double_f(fid, "/A", particles%fast_ion%A, dims, error)
particles%fast_ion%A = beam_mass
call h5ltread_dataset_double_f(fid, "/r", particles%fast_ion%r, dims, error)
call h5ltread_dataset_double_f(fid, "/z", particles%fast_ion%z, dims, error)
call h5ltpath_valid_f(fid, "/phi", .True., path_valid, error)
if(path_valid) then
call h5ltread_dataset_double_f(fid, "/phi", particles%fast_ion%phi, dims, error)
particles%axisym = .False.
endif
call h5ltread_dataset_int_f(fid, "/class", particles%fast_ion%class, dims, error)
if(any(particles%fast_ion%class.le.0)) then
if(inputs%verbose.ge.0) then
write(*,'(a)') 'READ_MC: Orbit class ID must be greater than 0'
endif
stop
endif
if(inputs%split.ge.1) then
call h5ltread_dataset_int_scalar_f(fid, "/nclass", particles%nclass, error)
if(any(particles%fast_ion%class.gt.particles%nclass)) then
if(inputs%verbose.ge.0) then
write(*,'(a)') 'READ_MC: Orbit class ID greater than the number of classes'
endif
stop
endif
endif
if(inputs%dist_type.eq.2) then
if(particles%axisym) then
dist_type_name = "Axisymmetric Guiding Center Monte Carlo"
else
dist_type_name = "Non-axisymmetric Guiding Center Monte Carlo"
endif
call h5ltread_dataset_double_f(fid, "/energy", particles%fast_ion%energy, dims, error)
call h5ltread_dataset_double_f(fid, "/pitch", particles%fast_ion%pitch, dims, error)
particles%fast_ion%vabs = sqrt(particles%fast_ion%energy/(v2_to_E_per_amu*particles%fast_ion%A))
else
if(particles%axisym) then
dist_type_name = "Axisymmetric Full Orbit Monte Carlo"
else
dist_type_name = "Non-axisymmetric Full Orbit Monte Carlo"
endif
call h5ltread_dataset_double_f(fid, "/vr", particles%fast_ion%vr, dims, error)
call h5ltread_dataset_double_f(fid, "/vt", particles%fast_ion%vt, dims, error)
call h5ltread_dataset_double_f(fid, "/vz", particles%fast_ion%vz, dims, error)
particles%fast_ion%vabs = sqrt(particles%fast_ion%vr**2 + &
particles%fast_ion%vt**2 + &
particles%fast_ion%vz**2)
particles%fast_ion%energy = v2_to_E_per_amu*particles%fast_ion%A*particles%fast_ion%vabs**2
endif
call h5ltread_dataset_double_f(fid, "/weight", weight, dims, error)
cnt=0
e1_xyz = matmul(beam_grid%inv_basis,[1.0,0.0,0.0])
e2_xyz = matmul(beam_grid%inv_basis,[0.0,1.0,0.0])
!$OMP PARALLEL DO schedule(guided) private(i,ii,j,ir,iz,iphi,fields,uvw,phi,ri,vi, &
!$OMP& delta_phi,beam_grid_phi_enter,beam_grid_phi_exit,C_xyz,xyz,xp,yp,zp,dum,inp)
particle_loop: do i=1,particles%nparticle
if(inputs%verbose.ge.2) then
WRITE(*,'(f7.2,"% completed",a,$)') cnt/real(particles%nparticle)*100,char(13)
endif
if(particles%axisym) then
uvw = [particles%fast_ion(i)%r, 0.d0 , particles%fast_ion(i)%z]
else
xp = particles%fast_ion(i)%r * cos(particles%fast_ion(i)%phi)
yp = particles%fast_ion(i)%r * sin(particles%fast_ion(i)%phi)
zp = particles%fast_ion(i)%z
uvw = [xp,yp,zp]
endif
call in_plasma(uvw,inp,input_coords=1)
if(.not.inp) cycle particle_loop
if(particles%axisym) then
beam_grid_phi_enter = 0.0
beam_grid_phi_exit = 0.0
dum = [0.d0, 0.d0, particles%fast_ion(i)%z]
call uvw_to_xyz(dum, C_xyz)
call circle_grid_intersect(C_xyz,e1_xyz,e2_xyz,particles%fast_ion(i)%r,beam_grid_phi_enter,beam_grid_phi_exit)
delta_phi = beam_grid_phi_exit-beam_grid_phi_enter
if(delta_phi.gt.0) then
particles%fast_ion(i)%beam_grid_cross_grid = .True.
else
particles%fast_ion(i)%beam_grid_cross_grid = .False.
delta_phi = 2*pi
endif
particles%fast_ion(i)%beam_grid_phi_enter = beam_grid_phi_enter
else
delta_phi = 2*pi
call uvw_to_xyz(uvw,xyz)
particles%fast_ion(i)%beam_grid_cross_grid = in_grid(xyz)
particles%fast_ion(i)%beam_grid_phi_enter = particles%fast_ion(i)%phi
endif
particles%fast_ion(i)%delta_phi = delta_phi
particles%fast_ion(i)%weight = weight(i)
ir = minloc(abs(inter_grid%r - particles%fast_ion(i)%r),1)
iphi = minloc(abs(inter_grid%phi - particles%fast_ion(i)%phi),1)
iz = minloc(abs(inter_grid%z - particles%fast_ion(i)%z),1)
!$OMP ATOMIC UPDATE
equil%plasma(ir,iz,iphi)%denf = equil%plasma(ir,iz,iphi)%denf + weight(i) / &
(particles%fast_ion(i)%r*inter_grid%dv)
!$OMP END ATOMIC
cnt=cnt+1
enddo particle_loop
!$OMP END PARALLEL DO
num = count(particles%fast_ion%beam_grid_cross_grid)
if(num.le.0) then
if(inputs%verbose.ge.0) then
write(*,'(a)') 'READ_MC: No mc particles in beam grid'
endif
stop
endif
if(inputs%verbose.ge.1) then
write(*,'(T2,"Distribution type: ",a)') dist_type_name
write(*,'(T2,"Number of mc particles: ",i10)') particles%nparticle
write(*,'(T2,"Number of orbit classes: ",i6)') particles%nclass
write(*,*) ''
endif
end subroutine read_mc