pfida_mc Subroutine

  1. fidasim.f90
  2. libfida
  3. pfida_mc

public subroutine pfida_mc()

Calculate Passive FIDA emission using a Monte Carlo Fast-ion distribution

Arguments

None

Calls

proc~~pfida_mc~~CallsGraph proc~pfida_mc pfida_mc proc~bt_cx_rates bt_cx_rates proc~pfida_mc->proc~bt_cx_rates proc~store_fida_photons store_fida_photons proc~pfida_mc->proc~store_fida_photons proc~get_plasma get_plasma proc~pfida_mc->proc~get_plasma proc~parallel_merge_reservoirs parallel_merge_reservoirs proc~pfida_mc->proc~parallel_merge_reservoirs proc~track_cylindrical track_cylindrical proc~pfida_mc->proc~track_cylindrical proc~store_photon_birth store_photon_birth proc~pfida_mc->proc~store_photon_birth interface~randu randu proc~pfida_mc->interface~randu proc~get_fields get_fields proc~pfida_mc->proc~get_fields proc~colrad colrad proc~pfida_mc->proc~colrad interface~parallel_sum parallel_sum proc~pfida_mc->interface~parallel_sum proc~gyro_correction gyro_correction proc~pfida_mc->proc~gyro_correction interface~interpol_coeff interpol_coeff proc~bt_cx_rates->interface~interpol_coeff proc~store_photons store_photons proc~store_fida_photons->proc~store_photons proc~xyz_to_uvw xyz_to_uvw proc~get_plasma->proc~xyz_to_uvw proc~uvw_to_xyz uvw_to_xyz proc~get_plasma->proc~uvw_to_xyz proc~in_plasma in_plasma proc~get_plasma->proc~in_plasma proc~get_position get_position proc~get_plasma->proc~get_position proc~parallel_merge_reservoirs->interface~parallel_sum proc~init_reservoir init_reservoir proc~parallel_merge_reservoirs->proc~init_reservoir proc~my_rank my_rank proc~parallel_merge_reservoirs->proc~my_rank proc~merge_reservoirs merge_reservoirs proc~parallel_merge_reservoirs->proc~merge_reservoirs proc~num_ranks num_ranks proc~parallel_merge_reservoirs->proc~num_ranks proc~track_cylindrical->proc~get_plasma proc~track_cylindrical->proc~get_fields proc~get_passive_grid_indices get_passive_grid_indices proc~track_cylindrical->proc~get_passive_grid_indices proc~cyl_to_uvw cyl_to_uvw proc~track_cylindrical->proc~cyl_to_uvw proc~line_cylinder_intersect line_cylinder_intersect proc~track_cylindrical->proc~line_cylinder_intersect proc~doppler_stark doppler_stark proc~track_cylindrical->proc~doppler_stark proc~track_cylindrical->proc~in_plasma proc~plane_basis plane_basis proc~track_cylindrical->proc~plane_basis proc~line_plane_intersect line_plane_intersect proc~track_cylindrical->proc~line_plane_intersect proc~store_photon_birth->proc~get_passive_grid_indices proc~update_reservoir update_reservoir proc~store_photon_birth->proc~update_reservoir proc~get_indices get_indices proc~store_photon_birth->proc~get_indices proc~calc_perp_vectors calc_perp_vectors proc~get_fields->proc~calc_perp_vectors proc~get_fields->proc~xyz_to_uvw proc~get_fields->proc~uvw_to_xyz proc~get_fields->proc~in_plasma proc~get_fields->proc~get_position proc~get_rate_matrix get_rate_matrix proc~colrad->proc~get_rate_matrix proc~eigen eigen proc~colrad->proc~eigen proc~linsolve linsolve proc~colrad->proc~linsolve proc~parallel_sum_i3 parallel_sum_i3 interface~parallel_sum->proc~parallel_sum_i3 proc~parallel_sum_d0 parallel_sum_d0 interface~parallel_sum->proc~parallel_sum_d0 proc~parallel_sum_d1 parallel_sum_d1 interface~parallel_sum->proc~parallel_sum_d1 proc~parallel_sum_d5 parallel_sum_d5 interface~parallel_sum->proc~parallel_sum_d5 proc~parallel_sum_i1 parallel_sum_i1 interface~parallel_sum->proc~parallel_sum_i1 proc~parallel_sum_i0 parallel_sum_i0 interface~parallel_sum->proc~parallel_sum_i0 proc~parallel_sum_i2 parallel_sum_i2 interface~parallel_sum->proc~parallel_sum_i2 proc~parallel_sum_d4 parallel_sum_d4 interface~parallel_sum->proc~parallel_sum_d4 proc~parallel_sum_d3 parallel_sum_d3 interface~parallel_sum->proc~parallel_sum_d3 proc~parallel_sum_d2 parallel_sum_d2 interface~parallel_sum->proc~parallel_sum_d2 proc~gyro_correction->interface~randu proc~gyro_step gyro_step proc~gyro_correction->proc~gyro_step proc~pitch_to_vec pitch_to_vec proc~gyro_correction->proc~pitch_to_vec proc~xyz_to_cyl xyz_to_cyl proc~get_passive_grid_indices->proc~xyz_to_cyl proc~parallel_sum_i3->proc~parallel_sum_i1 mpi_allreduce mpi_allreduce proc~parallel_sum_i3->mpi_allreduce proc~init_reservoir->interface~randu proc~parallel_sum_d0->mpi_allreduce proc~parallel_sum_d1->proc~parallel_sum_d1 proc~parallel_sum_d1->mpi_allreduce proc~parallel_sum_d5->proc~parallel_sum_d1 proc~parallel_sum_d5->mpi_allreduce proc~parallel_sum_i1->proc~parallel_sum_i1 proc~parallel_sum_i1->mpi_allreduce proc~parallel_sum_i0->mpi_allreduce proc~merge_reservoirs->interface~randu proc~merge_reservoirs->proc~init_reservoir proc~merge_reservoirs->proc~update_reservoir float float proc~merge_reservoirs->float proc~parallel_sum_i2->proc~parallel_sum_i1 proc~parallel_sum_i2->mpi_allreduce proc~update_reservoir->interface~randu proc~update_reservoir->proc~init_reservoir interface~randind randind proc~update_reservoir->interface~randind proc~get_rate_matrix->interface~interpol_coeff proc~elmhes elmhes proc~eigen->proc~elmhes proc~balance balance proc~eigen->proc~balance proc~hqr2 hqr2 proc~eigen->proc~hqr2 proc~balback balback proc~eigen->proc~balback proc~elmtrans elmtrans proc~eigen->proc~elmtrans proc~cross_product cross_product proc~gyro_step->proc~cross_product proc~in_plasma->proc~xyz_to_uvw proc~in_plasma->proc~cyl_to_uvw proc~in_plasma->interface~interpol_coeff proc~parallel_sum_d4->proc~parallel_sum_d1 proc~parallel_sum_d4->mpi_allreduce proc~plane_basis->proc~cross_product dgetrs dgetrs proc~linsolve->dgetrs proc~matinv matinv proc~linsolve->proc~matinv dgetrf dgetrf proc~linsolve->dgetrf proc~parallel_sum_d3->proc~parallel_sum_d1 proc~parallel_sum_d3->mpi_allreduce proc~cyl_to_xyz cyl_to_xyz proc~get_position->proc~cyl_to_xyz proc~store_photons->proc~get_fields proc~store_photons->proc~get_passive_grid_indices proc~store_photons->proc~uvw_to_xyz proc~store_photons->proc~get_indices proc~spectrum spectrum proc~store_photons->proc~spectrum proc~parallel_sum_d2->proc~parallel_sum_d1 proc~parallel_sum_d2->mpi_allreduce proc~rswap RSWAP proc~elmhes->proc~rswap proc~xyz_to_cyl->proc~xyz_to_uvw proc~uvw_to_cyl uvw_to_cyl proc~xyz_to_cyl->proc~uvw_to_cyl proc~balance->proc~rswap proc~hqrvec hqrvec proc~hqr2->proc~hqrvec proc~lubksb lubksb proc~matinv->proc~lubksb proc~ludcmp ludcmp proc~matinv->proc~ludcmp proc~balback->proc~rswap proc~cyl_to_xyz->proc~cyl_to_uvw proc~cyl_to_xyz->proc~uvw_to_xyz proc~comdiv Comdiv proc~hqrvec->proc~comdiv proc~outerprod outerprod proc~ludcmp->proc~outerprod proc~swap swap proc~ludcmp->proc~swap

Called by

proc~~pfida_mc~~CalledByGraph proc~pfida_mc pfida_mc program~fidasim fidasim program~fidasim->proc~pfida_mc

Contents

Source Code

pfida_mc

Source Code

subroutine pfida_mc
    !+ Calculate Passive FIDA emission using a Monte Carlo Fast-ion distribution
    integer :: iion,igamma,ngamma
    type(FastIon) :: fast_ion
    type(LocalEMFields) :: fields
    type(LocalProfiles) :: plasma
    real(Float64) :: phi
    real(Float64), dimension(3) :: ri      !! start position
    real(Float64), dimension(3) :: vi      !! velocity of fast ions
    real(Float64), dimension(3) :: xyz_vi
    !! Determination of the CX probability
    real(Float64), dimension(nlevs) :: denn    !!  neutral dens (n=1-4)
    real(Float64), dimension(nlevs) :: rates, rates_is    !! CX rates
    !! Collisiional radiative model along track
    real(Float64), dimension(nlevs) :: states  ! Density of n-states
    type(ParticleTrack), dimension(pass_grid%ntrack) :: tracks
    integer :: ntrack
    logical :: los_intersect
    integer :: jj, is
    real(Float64) :: photons
    real(Float64)  :: s, c
    real(Float64), dimension(1) :: randomu

    ngamma = 1
    if(particles%axisym.or.(inputs%dist_type.eq.2)) then
        ngamma = ceiling(dble(inputs%n_pfida)/particles%nparticle)
    endif

    if(inputs%verbose.ge.1) then
        write(*,'(T6,"# of markers: ",i10)') int(particles%nparticle*ngamma,Int64)
    endif

    !$OMP PARALLEL DO schedule(dynamic,1) private(iion,igamma,is,fast_ion,vi,ri,phi,tracks,s,c,&
    !$OMP& randomu,plasma,fields,ntrack,jj,rates,rates_is,denn,los_intersect,states,photons,xyz_vi)
    loop_over_fast_ions: do iion=istart,particles%nparticle,istep
        fast_ion = particles%fast_ion(iion)
        if(fast_ion%vabs.eq.0) cycle loop_over_fast_ions
        gamma_loop: do igamma=1,ngamma
            if(particles%axisym) then
                !! Pick random toroidal angle
                call randu(randomu)
                phi = pass_grid%phi(1)+pass_grid%nphi*pass_grid%dphi*randomu(1)
            else
                phi = fast_ion%phi
            endif
            s = sin(phi)
            c = cos(phi)

            !! Calculate position in machine coordinates
            ri(1) = fast_ion%r*c
            ri(2) = fast_ion%r*s
            ri(3) = fast_ion%z

            if(inputs%dist_type.eq.2) then
                !! Get electomagnetic fields
                call get_fields(fields, pos=ri, input_coords=1, output_coords=1)

                !! Correct for gyro motion and get particle position and velocity
                call gyro_correction(fields, fast_ion%energy, fast_ion%pitch, fast_ion%A, ri, vi)
            else !! Full Orbit
                !! Calculate velocity vector
                vi(1) = c*fast_ion%vr - s*fast_ion%vt
                vi(2) = s*fast_ion%vr + c*fast_ion%vt
                vi(3) = fast_ion%vz
            endif
            xyz_vi = matmul(beam_grid%inv_basis,vi)

            !! Track particle through grid
            call track_cylindrical(ri, vi, tracks, ntrack, los_intersect)
            if(.not.los_intersect) cycle gamma_loop
            if(ntrack.eq.0) cycle gamma_loop

            !! Calculate CX probability
            call get_plasma(plasma, pos=ri, input_coords=1)
            rates = 0.d0
            do is=1,n_thermal
                call bt_cx_rates(plasma, plasma%denn(:,is), thermal_mass(is), xyz_vi, rates_is)
                rates = rates + rates_is
            enddo
            if(sum(rates).le.0.) cycle gamma_loop

            !! Weight CX rates by ion source density
            if(particles%axisym) then
                states=rates*fast_ion%weight*(pass_grid%nphi*pass_grid%dphi/(2*pi))  &
                       /(fast_ion%r*pass_grid%dv)/ngamma

            else
                states=rates*fast_ion%weight/(fast_ion%r*pass_grid%dv)/ngamma
            endif

            !! Calculate the spectra produced in each cell along the path
            loop_along_track: do jj=1,ntrack
                call get_plasma(plasma,pos=tracks(jj)%pos,input_coords=1)

                call colrad(plasma, fast_ion%A, xyz_vi, tracks(jj)%time, states, denn, photons)

                call store_fida_photons(tracks(jj)%pos, xyz_vi, beam_lambda0, photons, fast_ion%class,passive=.True.)
                if(inputs%calc_res.ge.1) call store_photon_birth(tracks(1)%pos, photons, spatres%pfida, passive = .True.)
            enddo loop_along_track
        enddo gamma_loop
    enddo loop_over_fast_ions
    !$OMP END PARALLEL DO

#ifdef _MPI
    call parallel_sum(spec%pfida)
    if(inputs%calc_res.ge.1) then
        do jj=1,spec_chords%nchan
            call parallel_merge_reservoirs(spatres%pfida(jj))
        enddo
    endif
#endif

end subroutine pfida_mc