cfpd_f Subroutine

  1. fidasim.f90
  2. libfida
  3. cfpd_f

public subroutine cfpd_f()

Calculate charged fusion product count rate and weight function using a fast-ion distribution function F(E,p,r,z)

Arguments

None

Calls

proc~~cfpd_f~~CallsGraph proc~cfpd_f cfpd_f proc~convert_sightline_to_xyz convert_sightline_to_xyz proc~cfpd_f->proc~convert_sightline_to_xyz proc~my_rank my_rank proc~cfpd_f->proc~my_rank proc~time_string time_string proc~cfpd_f->proc~time_string proc~get_ep_denf get_ep_denf proc~cfpd_f->proc~get_ep_denf proc~get_pgyro get_pgyro proc~cfpd_f->proc~get_pgyro proc~get_plasma get_plasma proc~cfpd_f->proc~get_plasma proc~pitch_to_vec pitch_to_vec proc~cfpd_f->proc~pitch_to_vec proc~gyro_step gyro_step proc~cfpd_f->proc~gyro_step proc~get_ddpt_anisotropy get_ddpt_anisotropy proc~cfpd_f->proc~get_ddpt_anisotropy proc~write_cfpd_weights write_cfpd_weights proc~cfpd_f->proc~write_cfpd_weights proc~get_fields get_fields proc~cfpd_f->proc~get_fields proc~get_dd_rate get_dd_rate proc~cfpd_f->proc~get_dd_rate interface~parallel_sum parallel_sum proc~cfpd_f->interface~parallel_sum proc~uvw_to_xyz uvw_to_xyz proc~convert_sightline_to_xyz->proc~uvw_to_xyz proc~xyz_to_uvw xyz_to_uvw proc~get_ep_denf->proc~xyz_to_uvw interface~interpol interpol proc~get_ep_denf->interface~interpol proc~get_position get_position proc~get_ep_denf->proc~get_position proc~cross_product cross_product proc~get_pgyro->proc~cross_product proc~get_plasma->proc~xyz_to_uvw proc~get_plasma->proc~uvw_to_xyz proc~in_plasma in_plasma proc~get_plasma->proc~in_plasma proc~get_plasma->proc~get_position proc~gyro_step->proc~cross_product interface~interpol_coeff interpol_coeff proc~get_ddpt_anisotropy->interface~interpol_coeff h5close_f h5close_f proc~write_cfpd_weights->h5close_f h5ltset_attribute_string_f h5ltset_attribute_string_f proc~write_cfpd_weights->h5ltset_attribute_string_f h5open_f h5open_f proc~write_cfpd_weights->h5open_f h5fcreate_f h5fcreate_f proc~write_cfpd_weights->h5fcreate_f interface~h5ltmake_compressed_dataset_double_f h5ltmake_compressed_dataset_double_f proc~write_cfpd_weights->interface~h5ltmake_compressed_dataset_double_f h5ltmake_dataset_int_f h5ltmake_dataset_int_f proc~write_cfpd_weights->h5ltmake_dataset_int_f h5fclose_f h5fclose_f proc~write_cfpd_weights->h5fclose_f 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_dd_rate->interface~interpol_coeff 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~parallel_sum_i3->proc~parallel_sum_i1 mpi_allreduce mpi_allreduce proc~parallel_sum_i3->mpi_allreduce 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~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_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 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_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_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_7 h5ltmake_compressed_dataset_double_f_7 interface~h5ltmake_compressed_dataset_double_f->proc~h5ltmake_compressed_dataset_double_f_7 proc~parallel_sum_i2->proc~parallel_sum_i1 proc~parallel_sum_i2->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~parallel_sum_d4->proc~parallel_sum_d1 proc~parallel_sum_d4->mpi_allreduce 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~parallel_sum_d2->proc~parallel_sum_d1 proc~parallel_sum_d2->mpi_allreduce proc~chunk_size chunk_size proc~h5ltmake_compressed_dataset_double_f_4->proc~chunk_size h5sclose_f h5sclose_f proc~h5ltmake_compressed_dataset_double_f_4->h5sclose_f h5ltmake_dataset_double_f h5ltmake_dataset_double_f proc~h5ltmake_compressed_dataset_double_f_4->h5ltmake_dataset_double_f h5pset_chunk_f h5pset_chunk_f proc~h5ltmake_compressed_dataset_double_f_4->h5pset_chunk_f h5pcreate_f h5pcreate_f proc~h5ltmake_compressed_dataset_double_f_4->h5pcreate_f h5pset_shuffle_f h5pset_shuffle_f proc~h5ltmake_compressed_dataset_double_f_4->h5pset_shuffle_f h5pset_deflate_f h5pset_deflate_f proc~h5ltmake_compressed_dataset_double_f_4->h5pset_deflate_f h5dcreate_f h5dcreate_f proc~h5ltmake_compressed_dataset_double_f_4->h5dcreate_f h5dwrite_f h5dwrite_f proc~h5ltmake_compressed_dataset_double_f_4->h5dwrite_f h5dclose_f h5dclose_f proc~h5ltmake_compressed_dataset_double_f_4->h5dclose_f h5pclose_f h5pclose_f proc~h5ltmake_compressed_dataset_double_f_4->h5pclose_f h5screate_simple_f h5screate_simple_f proc~h5ltmake_compressed_dataset_double_f_4->h5screate_simple_f proc~h5ltmake_compressed_dataset_double_f_2->proc~chunk_size 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->h5pset_shuffle_f proc~h5ltmake_compressed_dataset_double_f_2->h5pset_deflate_f proc~h5ltmake_compressed_dataset_double_f_2->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->h5pclose_f proc~h5ltmake_compressed_dataset_double_f_2->h5screate_simple_f proc~h5ltmake_compressed_dataset_double_f_1->proc~chunk_size 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->h5pset_shuffle_f proc~h5ltmake_compressed_dataset_double_f_1->h5pset_deflate_f proc~h5ltmake_compressed_dataset_double_f_1->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->h5pclose_f proc~h5ltmake_compressed_dataset_double_f_1->h5screate_simple_f proc~h5ltmake_compressed_dataset_double_f_5->proc~chunk_size 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->h5pset_shuffle_f proc~h5ltmake_compressed_dataset_double_f_5->h5pset_deflate_f proc~h5ltmake_compressed_dataset_double_f_5->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->h5pclose_f proc~h5ltmake_compressed_dataset_double_f_5->h5screate_simple_f proc~h5ltmake_compressed_dataset_double_f_3->proc~chunk_size 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->h5pset_shuffle_f proc~h5ltmake_compressed_dataset_double_f_3->h5pset_deflate_f proc~h5ltmake_compressed_dataset_double_f_3->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->h5pclose_f proc~h5ltmake_compressed_dataset_double_f_3->h5screate_simple_f proc~cyl_to_xyz->proc~uvw_to_xyz proc~cyl_to_xyz->proc~cyl_to_uvw proc~h5ltmake_compressed_dataset_double_f_6->proc~chunk_size 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->h5pset_shuffle_f proc~h5ltmake_compressed_dataset_double_f_6->h5pset_deflate_f proc~h5ltmake_compressed_dataset_double_f_6->h5dcreate_f 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->h5pclose_f proc~h5ltmake_compressed_dataset_double_f_6->h5screate_simple_f proc~h5ltmake_compressed_dataset_double_f_7->proc~chunk_size proc~h5ltmake_compressed_dataset_double_f_7->h5sclose_f proc~h5ltmake_compressed_dataset_double_f_7->h5ltmake_dataset_double_f proc~h5ltmake_compressed_dataset_double_f_7->h5pset_chunk_f proc~h5ltmake_compressed_dataset_double_f_7->h5pcreate_f proc~h5ltmake_compressed_dataset_double_f_7->h5pset_shuffle_f proc~h5ltmake_compressed_dataset_double_f_7->h5pset_deflate_f proc~h5ltmake_compressed_dataset_double_f_7->h5dcreate_f proc~h5ltmake_compressed_dataset_double_f_7->h5dwrite_f proc~h5ltmake_compressed_dataset_double_f_7->h5dclose_f proc~h5ltmake_compressed_dataset_double_f_7->h5pclose_f proc~h5ltmake_compressed_dataset_double_f_7->h5screate_simple_f

Called by

proc~~cfpd_f~~CalledByGraph proc~cfpd_f cfpd_f program~fidasim fidasim program~fidasim->proc~cfpd_f

Contents

Source Code

cfpd_f

Source Code

subroutine cfpd_f
    !+ Calculate charged fusion product count rate and weight function using a fast-ion distribution function F(E,p,r,z)
    real(Float64), dimension(3) :: vi, vi_norm, v3_xyz, xyz, r_gyro
    type(LocalProfiles) :: plasma
    type(LocalEMFields) :: fields
    real(Float64) :: pgyro, vnet_square, factor, vabs
    real(Float64) :: eb, pitch, erel, rate, kappa, gyro, fbm_denf
    integer :: ie, ip, ich, ie3, iray, ist, cnt

    if(.not.any(thermal_mass.eq.H2_amu)) then
        write(*,'(T2,a)') 'CFPD_F: Thermal Deuterium is not present in plasma'
        return
    endif
    if(any(thermal_mass.eq.H3_amu)) then
        write(*,'(T2,a)') 'CFPD_F: D-T cfpd production is not implemented'
    endif
    if(beam_mass.ne.H2_amu) then
        write(*,'(T2,a)') 'CFPD_F: Fast-ion species is not Deuterium'
        return
    endif

    allocate(cfpd%flux(ctable%nenergy, ctable%nchan))
    allocate(cfpd%prob(ctable%nenergy, ctable%nchan))
    allocate(cfpd%gam(ctable%nenergy, ctable%nchan))
    allocate(cfpd%weight(ctable%nenergy, ctable%nchan, fbm%nenergy, fbm%npitch))
    cfpd%flux = 0.d0
    cfpd%prob = 0.d0
    cfpd%gam = 0.d0
    cfpd%weight = 0.d0

    rate = 0.d0
    factor = 0.5d0*fbm%dE*fbm%dp*ctable%dl !0.5 for TRANSP-pitch (E,p) space factor
    !$OMP PARALLEL DO schedule(guided) private(vi,vi_norm,v3_xyz,xyz,r_gyro,plasma,fields,pgyro,&
    !$OMP& vnet_square,vabs,eb,pitch,erel,rate,kappa,gyro,fbm_denf,ie,ip,ich,ie3,iray,ist,cnt)
    channel_loop: do ich=1, ctable%nchan
        E3_loop: do ie3=1, ctable%nenergy
            cnt = 0
            ray_loop: do iray=1, ctable%nrays
                step_loop: do ist=1, ctable%nsteps
                    if (ist.gt.ctable%nactual(ie3,iray,ich)) cycle ray_loop

                    !! Calculate position and velocity in beam coordinates
                    call convert_sightline_to_xyz(ie3, ist, iray, ich, xyz, v3_xyz)

                    !! Get fields at sightline position
                    call get_fields(fields, pos=xyz)
                    if(.not.fields%in_plasma) cycle step_loop

                    !! Get plasma parameters at sightline position
                    call get_plasma(plasma, pos=xyz)
                    if(.not.plasma%in_plasma) cycle step_loop

                    !! Loop over energy/pitch/gamma
                    pitch_loop: do ip = 1, fbm%npitch
                        pitch = fbm%pitch(ip)
                        energy_loop: do ie =1, fbm%nenergy
                            eb = fbm%energy(ie)

                            !! Get the probability factor
                            call get_pgyro(fields,ctable%earray(ie3),eb,pitch,plasma,v3_xyz,pgyro,gyro)
                            if (pgyro.le.0.d0) cycle energy_loop
                            cnt = cnt + 1

                            !! Calculate fast-ion velocity
                            call pitch_to_vec(pitch, gyro, fields, vi_norm)
                            vabs = sqrt(eb/(v2_to_E_per_amu*fbm%A))
                            vi = vi_norm*vabs
                            !!Correct for gyro orbit
                            call gyro_step(vi,fields,fbm%A,r_gyro)

                            fbm_denf=0
                            if (inputs%dist_type.eq.1) then
                                !get F at guiding center position
                                call get_ep_denf(eb,pitch,fbm_denf,pos=(xyz+r_gyro))
                            endif
                            if (fbm_denf.ne.fbm_denf) cycle energy_loop

                            !! Calculate effective beam energy
                            vnet_square=dot_product(vi-plasma%vrot,vi-plasma%vrot)  ![cm/s]
                            erel = v2_to_E_per_amu*fbm%A*vnet_square ![kev]

                            !! Get the cfpd production rate and anisotropy term
                            call get_dd_rate(plasma, erel, rate, branch=1)
                            call get_ddpt_anisotropy(plasma, vi, v3_xyz, kappa)

                            !$OMP CRITICAL(cfpd_weight)
                            cfpd%weight(ie3,ich,ie,ip) = cfpd%weight(ie3,ich,ie,ip) &
                                                            + rate * kappa * pgyro &
                                                            * ctable%daomega(ie3,iray,ich) &
                                                            * factor / (fbm%dE*fbm%dp)
                            cfpd%prob(ie3,ich) = cfpd%prob(ie3,ich) + pgyro
                            cfpd%gam(ie3,ich) = cfpd%gam(ie3,ich) + gyro
                            cfpd%flux(ie3,ich) = cfpd%flux(ie3,ich) + rate * kappa * pgyro &
                                                                        * ctable%daomega(ie3,iray,ich) &
                                                                        * fbm_denf * factor
                            !$OMP END CRITICAL(cfpd_weight)

                        enddo energy_loop
                    enddo pitch_loop
                enddo step_loop
            enddo ray_loop
            !$OMP CRITICAL(pweight_cnt)
            cfpd%prob(ie3,ich) = cfpd%prob(ie3,ich) / cnt
            cfpd%gam(ie3,ich) = cfpd%gam(ie3,ich) / cnt
            !$OMP END CRITICAL(pweight_cnt)
        enddo E3_loop
    enddo channel_loop
    !$OMP END PARALLEL DO

#ifdef _MPI
    call parallel_sum(cfpd%flux)
    call parallel_sum(cfpd%weight)
    call parallel_sum(cfpd%prob)
    call parallel_sum(cfpd%gam)
#endif

    if(inputs%verbose.ge.1) then
        write(*,'(30X,a)') ''
        write(*,*) 'write charged fusion products:    ' , time_string(time_start)
    endif

#ifdef _MPI
    if(my_rank().eq.0) call write_cfpd_weights()
#else
    call write_cfpd_weights()
#endif

end subroutine cfpd_f