neutron_mc Subroutine

  1. fidasim.f90
  2. libfida
  3. neutron_mc

public subroutine neutron_mc()

Calculate neutron flux using a Monte Carlo Fast-ion distribution

Arguments

None

Calls

proc~~neutron_mc~~CallsGraph proc~neutron_mc neutron_mc proc~time_string time_string proc~neutron_mc->proc~time_string proc~my_rank my_rank proc~neutron_mc->proc~my_rank proc~store_neutrons store_neutrons proc~neutron_mc->proc~store_neutrons proc~get_plasma get_plasma proc~neutron_mc->proc~get_plasma proc~write_neutrons write_neutrons proc~neutron_mc->proc~write_neutrons proc~get_fields get_fields proc~neutron_mc->proc~get_fields proc~get_dd_rate get_dd_rate proc~neutron_mc->proc~get_dd_rate interface~parallel_sum parallel_sum proc~neutron_mc->interface~parallel_sum proc~gyro_correction gyro_correction proc~neutron_mc->proc~gyro_correction 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 h5close_f h5close_f proc~write_neutrons->h5close_f h5open_f h5open_f proc~write_neutrons->h5open_f h5fcreate_f h5fcreate_f proc~write_neutrons->h5fcreate_f h5ltmake_dataset_double_f h5ltmake_dataset_double_f proc~write_neutrons->h5ltmake_dataset_double_f interface~h5ltmake_compressed_dataset_double_f h5ltmake_compressed_dataset_double_f proc~write_neutrons->interface~h5ltmake_compressed_dataset_double_f h5ltmake_dataset_int_f h5ltmake_dataset_int_f proc~write_neutrons->h5ltmake_dataset_int_f h5fclose_f h5fclose_f proc~write_neutrons->h5fclose_f h5ltset_attribute_string_f h5ltset_attribute_string_f proc~write_neutrons->h5ltset_attribute_string_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 interface~interpol_coeff interpol_coeff 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~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~randu randu proc~gyro_correction->interface~randu 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~cross_product cross_product proc~gyro_step->proc~cross_product 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~h5ltmake_compressed_dataset_double_f_4->h5ltmake_dataset_double_f 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 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->h5ltmake_dataset_double_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->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->h5ltmake_dataset_double_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->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->h5ltmake_dataset_double_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->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->h5ltmake_dataset_double_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->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->h5ltmake_dataset_double_f 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->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->h5ltmake_dataset_double_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->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~~neutron_mc~~CalledByGraph proc~neutron_mc neutron_mc program~fidasim fidasim program~fidasim->proc~neutron_mc

Contents

Source Code

neutron_mc

Source Code

subroutine neutron_mc
    !+ Calculate neutron flux using a Monte Carlo Fast-ion distribution
    integer :: iion, ngamma, igamma
    type(FastIon) :: fast_ion
    type(LocalProfiles) :: plasma
    type(LocalEMFields) :: fields
    real(Float64) :: eb, rate
    real(Float64), dimension(3) :: ri
    real(Float64), dimension(3) :: vi
    real(Float64), dimension(3) :: uvw, uvw_vi
    real(Float64)  :: vnet_square
    real(Float64)  :: phi, s, c, factor, delta_phi

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


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

    !! Correct neutron rate when equilibrium is 3D and MC distribution is 4D
    if(particles%axisym.and.(inter_grid%nphi.gt.1)) then
        delta_phi = inter_grid%phi(inter_grid%nphi)-inter_grid%phi(1)
        delta_phi = delta_phi + delta_phi/(inter_grid%nphi-1)/2 !Add half a cell
        factor = delta_phi/(2*pi) * 2 !Riemann sum below assumes coord's are at midpoint of cell
    else
        factor = 1
    endif

    rate=0.0
    ngamma = 20
    !$OMP PARALLEL DO schedule(guided) private(iion,fast_ion,vi,ri,s,c, &
    !$OMP& plasma,fields,uvw,uvw_vi,vnet_square,rate,eb,igamma,phi)
    loop_over_fast_ions: do iion=istart,particles%nparticle,istep
        fast_ion = particles%fast_ion(iion)
        if(fast_ion%vabs.eq.0.d0) cycle loop_over_fast_ions

        !! Calculate position in machine coordinates
        if(particles%axisym) then
            s = 0.d0
            c = 1.d0
        else
            phi = fast_ion%phi
            s = sin(phi)
            c = cos(phi)
        endif

        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)
            if(.not.fields%in_plasma) cycle loop_over_fast_ions

            gyro_loop: do igamma=1,ngamma
                !! Correct for Gyro-motion
                call gyro_correction(fields, fast_ion%energy, fast_ion%pitch, fast_ion%A, ri, vi)

                !! Get plasma parameters
                call get_plasma(plasma,pos=ri)
                if(.not.plasma%in_plasma) cycle gyro_loop

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

                !! Get neutron production rate
                call get_dd_rate(plasma, eb, rate, branch=2)
                rate = rate*fast_ion%weight/ngamma*factor

                !! Store neutrons
                call store_neutrons(rate, fast_ion%class)
            enddo gyro_loop
        else
            !! Get plasma parameters
            call get_plasma(plasma,pos=uvw,input_coords=1)
            if(.not.plasma%in_plasma) cycle loop_over_fast_ions

            !! Calculate effective beam energy
            uvw_vi(1) = fast_ion%vr
            uvw_vi(2) = fast_ion%vt
            uvw_vi(3) = fast_ion%vz
            vi = matmul(beam_grid%inv_basis,uvw_vi)
            vnet_square=dot_product(vi-plasma%vrot,vi-plasma%vrot)  ![cm/s]
            eb = v2_to_E_per_amu*fast_ion%A*vnet_square ![kev]

            !! Get neutron production rate
            call get_dd_rate(plasma, eb, rate, branch=2)
            rate = rate*fast_ion%weight*factor

            !! Store neutrons
            call store_neutrons(rate, fast_ion%class)
        endif
    enddo loop_over_fast_ions
    !$OMP END PARALLEL DO

#ifdef _MPI
    call parallel_sum(neutron%rate)
#endif

    if(inputs%verbose.ge.1) then
        write(*,'(T4,A,ES14.5," [neutrons/s]")') 'Rate:   ',sum(neutron%rate)
        write(*,'(30X,a)') ''
        write(*,*) 'write neutrons:    ' , time_string(time_start)
    endif

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

end subroutine neutron_mc