pnpa_mc Subroutine

  1. fidasim.f90
  2. libfida
  3. pnpa_mc

public subroutine pnpa_mc()

Calculate Passive NPA flux using a Monte Carlo fast-ion distribution

Arguments

None

Calls

proc~~pnpa_mc~~CallsGraph proc~pnpa_mc pnpa_mc proc~gyro_surface gyro_surface proc~pnpa_mc->proc~gyro_surface proc~bt_cx_rates bt_cx_rates proc~pnpa_mc->proc~bt_cx_rates proc~hit_npa_detector hit_npa_detector proc~pnpa_mc->proc~hit_npa_detector proc~store_npa store_npa proc~pnpa_mc->proc~store_npa proc~npa_gyro_range npa_gyro_range proc~pnpa_mc->proc~npa_gyro_range proc~get_plasma get_plasma proc~pnpa_mc->proc~get_plasma interface~randu randu proc~pnpa_mc->interface~randu proc~gyro_trajectory gyro_trajectory proc~pnpa_mc->proc~gyro_trajectory proc~get_fields get_fields proc~pnpa_mc->proc~get_fields proc~attenuate attenuate proc~pnpa_mc->proc~attenuate proc~uvw_to_xyz uvw_to_xyz proc~pnpa_mc->proc~uvw_to_xyz interface~parallel_sum parallel_sum proc~pnpa_mc->interface~parallel_sum interface~interpol_coeff interpol_coeff proc~bt_cx_rates->interface~interpol_coeff proc~in_boundary in_boundary proc~hit_npa_detector->proc~in_boundary proc~line_plane_intersect line_plane_intersect proc~hit_npa_detector->proc~line_plane_intersect proc~store_npa->proc~get_fields proc~xyz_to_uvw xyz_to_uvw proc~store_npa->proc~xyz_to_uvw proc~approx_eq approx_eq proc~npa_gyro_range->proc~approx_eq proc~gyro_range gyro_range proc~npa_gyro_range->proc~gyro_range proc~get_plasma->proc~uvw_to_xyz proc~get_position get_position proc~get_plasma->proc~get_position proc~in_plasma in_plasma proc~get_plasma->proc~in_plasma proc~get_plasma->proc~xyz_to_uvw proc~get_fields->proc~uvw_to_xyz proc~get_fields->proc~get_position proc~get_fields->proc~in_plasma proc~calc_perp_vectors calc_perp_vectors proc~get_fields->proc~calc_perp_vectors proc~get_fields->proc~xyz_to_uvw proc~attenuate->proc~get_plasma proc~colrad colrad proc~attenuate->proc~colrad proc~parallel_sum_d1 parallel_sum_d1 interface~parallel_sum->proc~parallel_sum_d1 proc~parallel_sum_i2 parallel_sum_i2 interface~parallel_sum->proc~parallel_sum_i2 proc~parallel_sum_i3 parallel_sum_i3 interface~parallel_sum->proc~parallel_sum_i3 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~parallel_sum_d4 parallel_sum_d4 interface~parallel_sum->proc~parallel_sum_d4 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_d0 parallel_sum_d0 interface~parallel_sum->proc~parallel_sum_d0 proc~cyl_to_xyz cyl_to_xyz proc~get_position->proc~cyl_to_xyz proc~parallel_sum_d1->proc~parallel_sum_d1 mpi_allreduce mpi_allreduce proc~parallel_sum_d1->mpi_allreduce proc~parallel_sum_i2->proc~parallel_sum_i1 proc~parallel_sum_i2->mpi_allreduce proc~parallel_sum_i3->proc~parallel_sum_i1 proc~parallel_sum_i3->mpi_allreduce proc~parallel_sum_d3->proc~parallel_sum_d1 proc~parallel_sum_d3->mpi_allreduce proc~in_plasma->proc~xyz_to_uvw proc~in_plasma->interface~interpol_coeff proc~cyl_to_uvw cyl_to_uvw proc~in_plasma->proc~cyl_to_uvw 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_d2->proc~parallel_sum_d1 proc~parallel_sum_d2->mpi_allreduce proc~parallel_sum_d4->proc~parallel_sum_d1 proc~parallel_sum_d4->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~parallel_sum_d0->mpi_allreduce proc~gyro_range->proc~in_boundary proc~gyro_range->proc~line_plane_intersect proc~gyro_surface_coordinates gyro_surface_coordinates proc~gyro_range->proc~gyro_surface_coordinates proc~line_gyro_surface_intersect line_gyro_surface_intersect proc~gyro_range->proc~line_gyro_surface_intersect proc~boundary_edge boundary_edge proc~gyro_range->proc~boundary_edge proc~in_gyro_surface in_gyro_surface proc~gyro_range->proc~in_gyro_surface proc~get_rate_matrix->interface~interpol_coeff proc~elmhes elmhes proc~eigen->proc~elmhes proc~hqr2 hqr2 proc~eigen->proc~hqr2 proc~balance balance proc~eigen->proc~balance proc~balback balback proc~eigen->proc~balback proc~elmtrans elmtrans proc~eigen->proc~elmtrans proc~cyl_to_xyz->proc~uvw_to_xyz proc~cyl_to_xyz->proc~cyl_to_uvw dgetrs dgetrs proc~linsolve->dgetrs proc~matinv matinv proc~linsolve->proc~matinv dgetrf dgetrf proc~linsolve->dgetrf proc~rswap RSWAP proc~elmhes->proc~rswap proc~hqrvec hqrvec proc~hqr2->proc~hqrvec proc~balance->proc~rswap proc~ludcmp ludcmp proc~matinv->proc~ludcmp proc~lubksb lubksb proc~matinv->proc~lubksb proc~balback->proc~rswap proc~swap swap proc~ludcmp->proc~swap proc~outerprod outerprod proc~ludcmp->proc~outerprod proc~comdiv Comdiv proc~hqrvec->proc~comdiv

Called by

proc~~pnpa_mc~~CalledByGraph proc~pnpa_mc pnpa_mc program~fidasim fidasim program~fidasim->proc~pnpa_mc

Contents

Source Code

pnpa_mc

Source Code

subroutine pnpa_mc
    !+ Calculate Passive NPA flux using a Monte Carlo fast-ion distribution
    integer :: iion,igamma,ngamma,npart
    type(FastIon) :: fast_ion
    real(Float64) :: phi,theta,dtheta
    real(Float64), dimension(3) :: ri, rf, rg, vi
    integer :: det,j,ichan,ir,nrange,it,is
    type(LocalEMFields) :: fields
    type(LocalProfiles) :: plasma
    type(GyroSurface) :: gs
    real(Float64), dimension(nlevs) :: rates, rates_is
    real(Float64), dimension(nlevs) :: states
    real(Float64) :: flux
    integer, dimension(3) :: ind
    real(Float64), dimension(3) :: uvw, uvw_vi
    real(Float64), dimension(2,4) :: gyrange
    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_pnpa)/particles%nparticle)
    endif

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

    !$OMP PARALLEL DO schedule(guided) private(iion,igamma,ind,fast_ion,vi,ri,rf,phi,s,c,ir,it,plasma, &
    !$OMP& randomu,rg,fields,uvw,uvw_vi,rates,rates_is,is,states,flux,det,ichan,gs,nrange,gyrange,theta,dtheta)
    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
            uvw(1) = fast_ion%r*c
            uvw(2) = fast_ion%r*s
            uvw(3) = fast_ion%z

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

                !! Correct for gyro motion and get position and velocity
                call gyro_surface(fields, fast_ion%energy, fast_ion%pitch, fast_ion%A, gs)

                detector_loop: do ichan=1,npa_chords%nchan
                    call npa_gyro_range(ichan, gs, gyrange, nrange)
                    if(nrange.eq.0) cycle detector_loop

                    gyro_range_loop: do ir=1,nrange
                        dtheta = gyrange(2,ir)
                        theta = gyrange(1,ir) + 0.5*dtheta
                        call gyro_trajectory(gs, theta, ri, vi)

                        !! Check if particle hits a NPA detector
                        call hit_npa_detector(ri, vi ,det, rf, det=ichan)
                        if(det.ne.ichan) then
                            if (inputs%verbose.ge.0)then
                                write(*,*) "PNPA_MC: Missed Detector ",ichan
                            endif
                            cycle gyro_range_loop
                        endif

                        !! Calculate CX probability with beam and halo neutrals
                        call get_plasma(plasma, pos=ri)
                        rates = 0.d0
                        do is=1,n_thermal
                            call bt_cx_rates(plasma, plasma%denn(:,is), thermal_mass(is), vi, rates_is)
                            rates = rates + rates_is
                        enddo
                        if(sum(rates).le.0.) cycle gyro_range_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

                        !! Attenuate states as the particle move through plasma
                        call attenuate(ri,rf,vi,states)

                        !! Store NPA Flux
                        flux = (dtheta/(2*pi))*sum(states)*(fast_ion%r*pass_grid%dv)
                        spread_loop: do it=1,25
                            theta = gyrange(1,ir) + (it-0.5)*dtheta/25
                            call gyro_trajectory(gs, theta, ri, vi)
                            call hit_npa_detector(ri, vi ,det, rf, det=ichan)
                            if(det.ne.ichan) then
                                if (inputs%verbose.ge.0)then
                                    write(*,*) "PNPA_MC: Missed Detector ",ichan
                                endif
                                cycle spread_loop
                            endif
                            call store_npa(det,ri,rf,vi,flux/25,fast_ion%class, passive=.True.)
                        enddo spread_loop
                    enddo gyro_range_loop
                enddo detector_loop
            else !! Full Orbit
                !! Convert to beam grid coordinates
                call uvw_to_xyz(uvw, ri)

                !! Calculate velocity vector
                uvw_vi(1) = c*fast_ion%vr - s*fast_ion%vt
                uvw_vi(2) = s*fast_ion%vr + c*fast_ion%vt
                uvw_vi(3) = fast_ion%vz
                vi = matmul(beam_grid%inv_basis,uvw_vi)

                !! Check if particle hits a NPA detector
                call hit_npa_detector(ri, vi ,det, rf)
                if(det.eq.0) cycle gamma_loop

                !! Calculate CX probability with edge neutrals
                call get_plasma(plasma, pos=ri)
                rates = 0.d0
                do is=1, n_thermal
                    call bt_cx_rates(plasma, plasma%denn(:,is), thermal_mass(is), 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

                !! Attenuate states as the particle moves though plasma
                call attenuate(ri,rf,vi,states)

                !! Store NPA Flux
                flux = sum(states)*(fast_ion%r*pass_grid%dv)
                call store_npa(det,ri,rf,vi,flux,fast_ion%class,passive=.True.)
            endif
        enddo gamma_loop
    enddo loop_over_fast_ions
    !$OMP END PARALLEL DO

    npart = pnpa%npart
#ifdef _MPI
    call parallel_sum(npart)
    call parallel_sum(pnpa%flux)
#endif

    if(inputs%verbose.ge.1) then
        write(*,'(T4,"Number of Passive NPA particles that hit a detector: ",i8)') npart
    endif

end subroutine pnpa_mc