fida_weights_mc Subroutine

public subroutine fida_weights_mc()

Calculates FIDA weights

Arguments

None

Calls

proc~~fida_weights_mc~~CallsGraph proc~fida_weights_mc fida_weights_mc proc~get_fields get_fields proc~fida_weights_mc->proc~get_fields proc~randu randu proc~fida_weights_mc->proc~randu proc~track track proc~fida_weights_mc->proc~track proc~get_nlaunch get_nlaunch proc~fida_weights_mc->proc~get_nlaunch proc~get_plasma get_plasma proc~fida_weights_mc->proc~get_plasma interface~randind randind proc~fida_weights_mc->interface~randind proc~write_fida_weights write_fida_weights proc~fida_weights_mc->proc~write_fida_weights proc~gyro_correction gyro_correction proc~fida_weights_mc->proc~gyro_correction proc~get_ep_denf get_ep_denf proc~fida_weights_mc->proc~get_ep_denf proc~in_plasma in_plasma proc~get_fields->proc~in_plasma proc~calc_perp_vectors calc_perp_vectors proc~get_fields->proc~calc_perp_vectors omp_get_thread_num omp_get_thread_num proc~randu->omp_get_thread_num proc~rng_uniform rng_uniform proc~randu->proc~rng_uniform proc~get_indices get_indices proc~track->proc~get_indices proc~track->proc~in_plasma proc~get_plasma->proc~in_plasma h5fcreate_f h5fcreate_f proc~write_fida_weights->h5fcreate_f h5open_f h5open_f proc~write_fida_weights->h5open_f h5ltset_attribute_string_f h5ltset_attribute_string_f proc~write_fida_weights->h5ltset_attribute_string_f h5ltmake_dataset_int_f h5ltmake_dataset_int_f proc~write_fida_weights->h5ltmake_dataset_int_f h5close_f h5close_f proc~write_fida_weights->h5close_f interface~h5ltmake_compressed_dataset_double_f h5ltmake_compressed_dataset_double_f proc~write_fida_weights->interface~h5ltmake_compressed_dataset_double_f h5fclose_f h5fclose_f proc~write_fida_weights->h5fclose_f proc~gyro_correction->proc~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 interface~interpol interpol proc~get_ep_denf->interface~interpol proc~xyz_to_uvw xyz_to_uvw proc~get_ep_denf->proc~xyz_to_uvw proc~cross_product cross_product proc~gyro_step->proc~cross_product proc~in_plasma->proc~xyz_to_uvw interface~interpol_coeff interpol_coeff proc~in_plasma->interface~interpol_coeff proc~interpol2d_arr interpol2D_arr interface~interpol->proc~interpol2d_arr proc~interpol2d_2d_arr interpol2D_2D_arr interface~interpol->proc~interpol2d_2d_arr proc~interpol1d_arr interpol1D_arr interface~interpol->proc~interpol1d_arr proc~h5ltmake_compressed_dataset_double_f_7 h5ltmake_compressed_dataset_double_f_7 interface~h5ltmake_compressed_dataset_double_f->proc~h5ltmake_compressed_dataset_double_f_7 proc~h5ltmake_compressed_dataset_double_f_6 h5ltmake_compressed_dataset_double_f_6 interface~h5ltmake_compressed_dataset_double_f->proc~h5ltmake_compressed_dataset_double_f_6 proc~h5ltmake_compressed_dataset_double_f_5 h5ltmake_compressed_dataset_double_f_5 interface~h5ltmake_compressed_dataset_double_f->proc~h5ltmake_compressed_dataset_double_f_5 proc~h5ltmake_compressed_dataset_double_f_4 h5ltmake_compressed_dataset_double_f_4 interface~h5ltmake_compressed_dataset_double_f->proc~h5ltmake_compressed_dataset_double_f_4 proc~h5ltmake_compressed_dataset_double_f_3 h5ltmake_compressed_dataset_double_f_3 interface~h5ltmake_compressed_dataset_double_f->proc~h5ltmake_compressed_dataset_double_f_3 proc~h5ltmake_compressed_dataset_double_f_2 h5ltmake_compressed_dataset_double_f_2 interface~h5ltmake_compressed_dataset_double_f->proc~h5ltmake_compressed_dataset_double_f_2 proc~h5ltmake_compressed_dataset_double_f_1 h5ltmake_compressed_dataset_double_f_1 interface~h5ltmake_compressed_dataset_double_f->proc~h5ltmake_compressed_dataset_double_f_1 h5dwrite_f h5dwrite_f proc~h5ltmake_compressed_dataset_double_f_7->h5dwrite_f h5dclose_f h5dclose_f proc~h5ltmake_compressed_dataset_double_f_7->h5dclose_f h5sclose_f h5sclose_f proc~h5ltmake_compressed_dataset_double_f_7->h5sclose_f h5ltmake_dataset_double_f h5ltmake_dataset_double_f proc~h5ltmake_compressed_dataset_double_f_7->h5ltmake_dataset_double_f h5pset_chunk_f h5pset_chunk_f proc~h5ltmake_compressed_dataset_double_f_7->h5pset_chunk_f h5pcreate_f h5pcreate_f proc~h5ltmake_compressed_dataset_double_f_7->h5pcreate_f proc~chunk_size chunk_size proc~h5ltmake_compressed_dataset_double_f_7->proc~chunk_size h5pset_deflate_f h5pset_deflate_f proc~h5ltmake_compressed_dataset_double_f_7->h5pset_deflate_f h5pset_shuffle_f h5pset_shuffle_f proc~h5ltmake_compressed_dataset_double_f_7->h5pset_shuffle_f h5screate_simple_f h5screate_simple_f proc~h5ltmake_compressed_dataset_double_f_7->h5screate_simple_f h5pclose_f h5pclose_f proc~h5ltmake_compressed_dataset_double_f_7->h5pclose_f h5dcreate_f h5dcreate_f proc~h5ltmake_compressed_dataset_double_f_7->h5dcreate_f proc~interpol2d_arr->interface~interpol_coeff proc~interpol2d_2d_arr->interface~interpol_coeff proc~h5ltmake_compressed_dataset_double_f_6->h5dwrite_f proc~h5ltmake_compressed_dataset_double_f_6->h5dclose_f proc~h5ltmake_compressed_dataset_double_f_6->h5sclose_f proc~h5ltmake_compressed_dataset_double_f_6->h5ltmake_dataset_double_f proc~h5ltmake_compressed_dataset_double_f_6->h5pset_chunk_f proc~h5ltmake_compressed_dataset_double_f_6->h5pcreate_f proc~h5ltmake_compressed_dataset_double_f_6->proc~chunk_size proc~h5ltmake_compressed_dataset_double_f_6->h5pset_deflate_f proc~h5ltmake_compressed_dataset_double_f_6->h5pset_shuffle_f proc~h5ltmake_compressed_dataset_double_f_6->h5screate_simple_f proc~h5ltmake_compressed_dataset_double_f_6->h5pclose_f proc~h5ltmake_compressed_dataset_double_f_6->h5dcreate_f proc~h5ltmake_compressed_dataset_double_f_5->h5dwrite_f proc~h5ltmake_compressed_dataset_double_f_5->h5dclose_f proc~h5ltmake_compressed_dataset_double_f_5->h5sclose_f proc~h5ltmake_compressed_dataset_double_f_5->h5ltmake_dataset_double_f proc~h5ltmake_compressed_dataset_double_f_5->h5pset_chunk_f proc~h5ltmake_compressed_dataset_double_f_5->h5pcreate_f proc~h5ltmake_compressed_dataset_double_f_5->proc~chunk_size proc~h5ltmake_compressed_dataset_double_f_5->h5pset_deflate_f proc~h5ltmake_compressed_dataset_double_f_5->h5pset_shuffle_f proc~h5ltmake_compressed_dataset_double_f_5->h5screate_simple_f proc~h5ltmake_compressed_dataset_double_f_5->h5pclose_f proc~h5ltmake_compressed_dataset_double_f_5->h5dcreate_f proc~h5ltmake_compressed_dataset_double_f_4->h5dwrite_f proc~h5ltmake_compressed_dataset_double_f_4->h5dclose_f proc~h5ltmake_compressed_dataset_double_f_4->h5sclose_f proc~h5ltmake_compressed_dataset_double_f_4->h5ltmake_dataset_double_f proc~h5ltmake_compressed_dataset_double_f_4->h5pset_chunk_f proc~h5ltmake_compressed_dataset_double_f_4->h5pcreate_f proc~h5ltmake_compressed_dataset_double_f_4->proc~chunk_size proc~h5ltmake_compressed_dataset_double_f_4->h5pset_deflate_f proc~h5ltmake_compressed_dataset_double_f_4->h5pset_shuffle_f proc~h5ltmake_compressed_dataset_double_f_4->h5screate_simple_f proc~h5ltmake_compressed_dataset_double_f_4->h5pclose_f proc~h5ltmake_compressed_dataset_double_f_4->h5dcreate_f proc~h5ltmake_compressed_dataset_double_f_3->h5dwrite_f proc~h5ltmake_compressed_dataset_double_f_3->h5dclose_f proc~h5ltmake_compressed_dataset_double_f_3->h5sclose_f proc~h5ltmake_compressed_dataset_double_f_3->h5ltmake_dataset_double_f proc~h5ltmake_compressed_dataset_double_f_3->h5pset_chunk_f proc~h5ltmake_compressed_dataset_double_f_3->h5pcreate_f proc~h5ltmake_compressed_dataset_double_f_3->proc~chunk_size proc~h5ltmake_compressed_dataset_double_f_3->h5pset_deflate_f proc~h5ltmake_compressed_dataset_double_f_3->h5pset_shuffle_f proc~h5ltmake_compressed_dataset_double_f_3->h5screate_simple_f proc~h5ltmake_compressed_dataset_double_f_3->h5pclose_f proc~h5ltmake_compressed_dataset_double_f_3->h5dcreate_f proc~h5ltmake_compressed_dataset_double_f_2->h5dwrite_f proc~h5ltmake_compressed_dataset_double_f_2->h5dclose_f proc~h5ltmake_compressed_dataset_double_f_2->h5sclose_f proc~h5ltmake_compressed_dataset_double_f_2->h5ltmake_dataset_double_f proc~h5ltmake_compressed_dataset_double_f_2->h5pset_chunk_f proc~h5ltmake_compressed_dataset_double_f_2->h5pcreate_f proc~h5ltmake_compressed_dataset_double_f_2->proc~chunk_size proc~h5ltmake_compressed_dataset_double_f_2->h5pset_deflate_f proc~h5ltmake_compressed_dataset_double_f_2->h5pset_shuffle_f proc~h5ltmake_compressed_dataset_double_f_2->h5screate_simple_f proc~h5ltmake_compressed_dataset_double_f_2->h5pclose_f proc~h5ltmake_compressed_dataset_double_f_2->h5dcreate_f proc~h5ltmake_compressed_dataset_double_f_1->h5dwrite_f proc~h5ltmake_compressed_dataset_double_f_1->h5dclose_f proc~h5ltmake_compressed_dataset_double_f_1->h5sclose_f proc~h5ltmake_compressed_dataset_double_f_1->h5ltmake_dataset_double_f proc~h5ltmake_compressed_dataset_double_f_1->h5pset_chunk_f proc~h5ltmake_compressed_dataset_double_f_1->h5pcreate_f proc~h5ltmake_compressed_dataset_double_f_1->proc~chunk_size proc~h5ltmake_compressed_dataset_double_f_1->h5pset_deflate_f proc~h5ltmake_compressed_dataset_double_f_1->h5pset_shuffle_f proc~h5ltmake_compressed_dataset_double_f_1->h5screate_simple_f proc~h5ltmake_compressed_dataset_double_f_1->h5pclose_f proc~h5ltmake_compressed_dataset_double_f_1->h5dcreate_f proc~interpol2d_coeff interpol2D_coeff interface~interpol_coeff->proc~interpol2d_coeff proc~interpol1d_coeff interpol1D_coeff interface~interpol_coeff->proc~interpol1d_coeff proc~interpol2d_coeff_arr interpol2D_coeff_arr interface~interpol_coeff->proc~interpol2d_coeff_arr proc~interpol1d_coeff_arr interpol1D_coeff_arr interface~interpol_coeff->proc~interpol1d_coeff_arr proc~interpol1d_arr->interface~interpol_coeff proc~interpol2d_coeff_arr->proc~interpol2d_coeff proc~interpol1d_coeff_arr->proc~interpol1d_coeff

Called by

proc~~fida_weights_mc~~CalledByGraph proc~fida_weights_mc fida_weights_mc program~fidasim fidasim program~fidasim->proc~fida_weights_mc

Contents

Source Code


Source Code

subroutine fida_weights_mc
    !+ Calculates FIDA weights
    integer :: i,j,k !! indices  x,y,z of cells
    integer(Int64) :: iion,ip
    real(Float64), dimension(3) :: ri,rg      !! start position
    real(Float64), dimension(3) :: vi     !! velocity of fast ions
    integer,dimension(3) :: ind      !! new actual cell
    integer,dimension(4) :: neut_types=[1,2,3,4]
    logical :: los_intersect

    !! Determination of the CX rates
    type(LocalProfiles) :: plasma
    type(LocalEMFields) :: fields
    real(Float64), dimension(nlevs) :: rates !! CX rates

    !! Collisiional radiative model along track
    integer :: ncell
    integer :: jj      !! counter along track
    type(ParticleTrack),dimension(beam_grid%ntrack) :: tracks

    real(Float64) :: photons !! photon flux
    real(Float64), dimension(nlevs) :: states  !! Density of n-states
    real(Float64), dimension(nlevs) :: denn

    integer :: nwav
    real(Float64) :: etov2, energy, pitch
    real(Float64) :: dE, dP, dEdP
    real(Float64), dimension(:), allocatable :: ebarr, ptcharr
    integer, dimension(1) :: ienergy, ipitch
    real(Float64), dimension(3) :: randomu3

    !! Number of particles to launch
    integer(kind=8) :: pcnt
    real(Float64) :: papprox_tot,inv_maxcnt,cnt,fbm_denf,phase_area
    integer,dimension(3,beam_grid%ngrid) :: pcell
    real(Float64), dimension(beam_grid%nx,beam_grid%ny,beam_grid%nz) :: papprox,nlaunch !! approx. density

    nwav = inputs%nlambda_wght

    !! define arrays
    !! define energy - array
    allocate(ebarr(inputs%ne_wght))
    do i=1,inputs%ne_wght
        ebarr(i)=real(i-0.5)*inputs%emax_wght/real(inputs%ne_wght)
    enddo
    dE = abs(ebarr(2)-ebarr(1))

    !! define pitch - array
    allocate(ptcharr(inputs%np_wght))
    do i=1,inputs%np_wght
        ptcharr(i)=real(i-0.5)*2./real(inputs%np_wght)-1.
    enddo
    dP = abs(ptcharr(2)-ptcharr(1))
    dEdP = dE*dP
    phase_area = dEdP*real(inputs%np_wght)*real(inputs%ne_wght)

    !! allocate storage arrays
    allocate(fweight%weight(nwav,inputs%ne_wght,inputs%np_wght,spec_chords%nchan))
    allocate(fweight%mean_f(inputs%ne_wght,inputs%np_wght,spec_chords%nchan))

    if(inputs%verbose.ge.1) then
        write(*,'(T2,"Number of Channels: ",i5)') spec_chords%nchan
        write(*,'(T2,"Nlambda: ",i4)') nwav
        write(*,'(T2,"Nenergy: ",i3)') inputs%ne_wght
        write(*,'(T2,"Maximum Energy: ",f7.2)') inputs%emax_wght
        write(*,'(T2,"LOS averaged: ",a)') "False"
    endif

    !! zero out arrays
    fweight%weight = 0.d0
    fweight%mean_f = 0.d0

    etov2 = 1.d0/(v2_to_E_per_amu*inputs%ab)

    !! Estimate how many particles to launch in each cell
    papprox=0.d0
    papprox_tot=0.d0
    pcnt=1
    do k=1,beam_grid%nz
        do j=1,beam_grid%ny
            do i=1,beam_grid%nx
                ind =[i,j,k]
                call get_plasma(plasma,ind=ind)
                papprox(i,j,k)=(sum(neut%dens(:,nbif_type,i,j,k)) + &
                                sum(neut%dens(:,nbih_type,i,j,k)) + &
                                sum(neut%dens(:,nbit_type,i,j,k)) + &
                                sum(neut%dens(:,halo_type,i,j,k)))
                if(papprox(i,j,k).gt.0) then
                    pcell(:,pcnt)= ind
                    pcnt=pcnt+1
                endif
                if(plasma%in_plasma) papprox_tot=papprox_tot+papprox(i,j,k)
            enddo
        enddo
    enddo
    pcnt=pcnt-1
    inv_maxcnt=100.0/real(pcnt)
    call get_nlaunch(10*inputs%n_fida,papprox,papprox_tot,nlaunch)
    if(inputs%verbose.ge.1) then
        write(*,'(T6,"# of markers: ",i9)') int(sum(nlaunch),Int64)
    endif

    !! Loop over all cells that have neutrals
    cnt=0.d0
    loop_over_cells: do ip = 1, int(pcnt)
        i = pcell(1,ip)
        j = pcell(2,ip)
        k = pcell(3,ip)
        ind = [i, j, k]
        !$OMP PARALLEL DO schedule(guided) private(iion,vi,ri,rg,ienergy,ipitch, &
        !$OMP tracks,ncell,jj,plasma,fields,rates,denn,states,photons,energy,pitch, &
        !$OMP los_intersect,randomu3,fbm_denf)
        loop_over_fast_ions: do iion=1,int(nlaunch(i, j, k),Int64)
            !! Sample fast ion distribution uniformally
            call randind(inputs%ne_wght, ienergy)
            call randind(inputs%np_wght, ipitch)
            call randu(randomu3)
            energy = ebarr(ienergy(1)) + dE*(randomu3(1)-0.5)
            pitch = ptcharr(ipitch(1)) + dP*(randomu3(2)-0.5)
            if(energy.le.0) cycle loop_over_fast_ions

            call randu(randomu3)
            rg = [beam_grid%xc(i),beam_grid%yc(j),beam_grid%zc(k)] + beam_grid%dr*(randomu3-0.5)

            !! Get velocity
            call get_fields(fields,pos=rg)
            if(.not.fields%in_plasma) cycle loop_over_fast_ions
            call gyro_correction(fields,energy,pitch,ri,vi)

            fbm_denf = 0.0
            if (inputs%dist_type.eq.1) then
                call get_ep_denf(energy,pitch,fbm_denf,coeffs=fields%c)
            endif

            !! Find the particles path through the beam grid
            call track(ri, vi, tracks, ncell, los_intersect)
            if(.not.los_intersect) cycle loop_over_fast_ions
            if(ncell.eq.0) cycle loop_over_fast_ions

            !! Calculate CX probability with beam and halo neutrals
            call get_beam_cx_rate(tracks(1)%ind, ri, vi, beam_ion, neut_types, rates)
            if(sum(rates).le.0.) cycle loop_over_fast_ions
            states=rates*1.d20

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

                call colrad(plasma,beam_ion, vi, tracks(jj)%time, states, denn, photons)

                call store_fw_photons(ienergy(1), ipitch(1), &
                     tracks(jj)%pos, vi, fbm_denf, photons/nlaunch(i,j,k))
            enddo loop_along_track
        enddo loop_over_fast_ions
        !$OMP END PARALLEL DO
        cnt=cnt+1

        if(inputs%verbose.ge.2) then
            WRITE(*,'(f7.2,"% completed",a,$)') cnt*inv_maxcnt,char(13)
        endif
    enddo loop_over_cells

    fweight%weight = ((1.d-20)*phase_area/dEdP)*fweight%weight
    fweight%mean_f = ((1.d-20)*phase_area/dEdP)*fweight%mean_f
    call write_fida_weights()

end subroutine fida_weights_mc