ndmc Subroutine

  1. fidasim.f90
  2. libfida
  3. ndmc

public subroutine ndmc()

Calculates neutral beam deposition and spectra

Arguments

None

Calls

proc~~ndmc~~CallsGraph proc~ndmc ndmc proc~init_neutral_population init_neutral_population proc~ndmc->proc~init_neutral_population proc~mc_nbi mc_nbi proc~ndmc->proc~mc_nbi proc~get_plasma get_plasma proc~ndmc->proc~get_plasma proc~store_births store_births proc~ndmc->proc~store_births proc~store_neutrals store_neutrals proc~ndmc->proc~store_neutrals proc~store_nbi_photons store_nbi_photons proc~ndmc->proc~store_nbi_photons interface~randind randind proc~ndmc->interface~randind proc~gyro_step gyro_step proc~ndmc->proc~gyro_step proc~track track proc~ndmc->proc~track proc~colrad colrad proc~ndmc->proc~colrad interface~randu randu proc~ndmc->interface~randu proc~get_fields get_fields proc~ndmc->proc~get_fields proc~parallel_merge_populations parallel_merge_populations proc~ndmc->proc~parallel_merge_populations interface~parallel_sum parallel_sum proc~ndmc->interface~parallel_sum proc~mc_nbi->interface~randu interface~randn randn proc~mc_nbi->interface~randn proc~grid_intersect grid_intersect proc~mc_nbi->proc~grid_intersect proc~in_plasma in_plasma proc~mc_nbi->proc~in_plasma 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~get_plasma->proc~in_plasma proc~get_position get_position proc~get_plasma->proc~get_position proc~update_neutrals update_neutrals proc~store_neutrals->proc~update_neutrals proc~store_photons store_photons proc~store_nbi_photons->proc~store_photons proc~cross_product cross_product proc~gyro_step->proc~cross_product proc~track->proc~get_plasma proc~track->proc~get_fields proc~doppler_stark doppler_stark proc~track->proc~doppler_stark proc~track->proc~in_plasma proc~get_indices get_indices proc~track->proc~get_indices 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~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~parallel_merge_populations->interface~parallel_sum proc~parallel_merge_reservoirs parallel_merge_reservoirs proc~parallel_merge_populations->proc~parallel_merge_reservoirs proc~ind2sub ind2sub proc~parallel_merge_populations->proc~ind2sub 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~parallel_sum_i2->proc~parallel_sum_i1 proc~parallel_sum_i2->mpi_allreduce interface~interpol_coeff interpol_coeff 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~parallel_merge_reservoirs->interface~parallel_sum proc~my_rank my_rank proc~parallel_merge_reservoirs->proc~my_rank proc~num_ranks num_ranks proc~parallel_merge_reservoirs->proc~num_ranks proc~merge_reservoirs merge_reservoirs proc~parallel_merge_reservoirs->proc~merge_reservoirs proc~init_reservoir init_reservoir proc~parallel_merge_reservoirs->proc~init_reservoir proc~in_passive_grid in_passive_grid proc~grid_intersect->proc~in_passive_grid proc~in_plasma->proc~xyz_to_uvw proc~cyl_to_uvw cyl_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~parallel_sum_d3->proc~parallel_sum_d1 proc~parallel_sum_d3->mpi_allreduce dgetrs dgetrs proc~linsolve->dgetrs proc~matinv matinv proc~linsolve->proc~matinv dgetrf dgetrf proc~linsolve->dgetrf proc~cyl_to_xyz cyl_to_xyz proc~get_position->proc~cyl_to_xyz proc~store_photons->proc~get_fields proc~store_photons->proc~uvw_to_xyz proc~store_photons->proc~get_indices proc~get_passive_grid_indices get_passive_grid_indices proc~store_photons->proc~get_passive_grid_indices proc~spectrum spectrum proc~store_photons->proc~spectrum proc~update_reservoir update_reservoir proc~update_neutrals->proc~update_reservoir proc~parallel_sum_d2->proc~parallel_sum_d1 proc~parallel_sum_d2->mpi_allreduce proc~rswap RSWAP proc~elmhes->proc~rswap proc~merge_reservoirs->interface~randu proc~merge_reservoirs->proc~init_reservoir proc~merge_reservoirs->proc~update_reservoir float float proc~merge_reservoirs->float proc~init_reservoir->interface~randu proc~xyz_to_cyl xyz_to_cyl proc~get_passive_grid_indices->proc~xyz_to_cyl proc~hqrvec hqrvec proc~hqr2->proc~hqrvec proc~balance->proc~rswap proc~lubksb lubksb proc~matinv->proc~lubksb proc~ludcmp ludcmp proc~matinv->proc~ludcmp proc~update_reservoir->interface~randind proc~update_reservoir->interface~randu proc~update_reservoir->proc~init_reservoir proc~balback->proc~rswap proc~approx_le approx_le proc~in_passive_grid->proc~approx_le proc~approx_ge approx_ge proc~in_passive_grid->proc~approx_ge proc~uvw_to_cyl uvw_to_cyl proc~in_passive_grid->proc~uvw_to_cyl proc~cyl_to_xyz->proc~uvw_to_xyz proc~cyl_to_xyz->proc~cyl_to_uvw proc~xyz_to_cyl->proc~xyz_to_uvw proc~xyz_to_cyl->proc~uvw_to_cyl proc~comdiv Comdiv proc~hqrvec->proc~comdiv proc~approx_eq approx_eq proc~approx_le->proc~approx_eq proc~approx_ge->proc~approx_eq proc~outerprod outerprod proc~ludcmp->proc~outerprod proc~swap swap proc~ludcmp->proc~swap

Called by

proc~~ndmc~~CalledByGraph proc~ndmc ndmc program~fidasim fidasim program~fidasim->proc~ndmc

Contents

Source Code

ndmc

Source Code

subroutine ndmc
    !+ Calculates neutral beam deposition and spectra
    integer :: neut_type !! full half third energy
    real(Float64)  :: nlaunch   !! nr. of markers
    real(Float64)  :: nneutrals !! # NBI particles
    real(Float64), dimension(3) :: vnbi !! velocities(full..)
    real(Float64), dimension(3) :: rnbi !! initial position

    integer(Int64) :: jj, ii, kk, cnt
    integer :: ntrack, ic
    type(ParticleTrack), dimension(beam_grid%ntrack) :: tracks
    type(LocalProfiles) :: plasma
    type(LocalEMFields) :: fields
    real(Float64), dimension(nlevs) :: states, dens
    real(Float64) :: photons, iflux,flux_tot
    integer(Int32), dimension(3) :: ind
    real(Float64), dimension(3) :: ri,ri_gc,r_gyro
    real(Float64), dimension(1) :: randomu
    integer, dimension(1) :: randi
    logical :: err

    !! Initialize Neutral Population
    call init_neutral_population(neut%full)
    call init_neutral_population(neut%half)
    call init_neutral_population(neut%third)

    if(inputs%verbose.ge.1) then
        write(*,'(T6,"# of markers: ",i10)') inputs%n_nbi
        if(inputs%calc_birth.ge.1) then
            write(*,'(T6,"# of birth markers: 3 x",i10)') int(inputs%n_nbi*inputs%n_birth)
        endif
    endif

    !! # of injected neutrals = NBI power/energy_per_particle
    nneutrals=1.d6*nbi%pinj/ (1.d3*nbi%einj*e0 &
         *( nbi%current_fractions(1)      &
         +  nbi%current_fractions(2)/2.d0 &
         +  nbi%current_fractions(3)/3.d0 ) )

    nlaunch=real(inputs%n_nbi)
    !$OMP PARALLEL DO schedule(guided) &
    !$OMP& private(vnbi,rnbi,tracks,ntrack,plasma,fields,randi,flux_tot, &
    !$OMP& states,dens,iflux,photons,neut_type,jj,ii,kk,ind,err,ri,ri_gc,r_gyro)
    loop_over_markers: do ii=istart,inputs%n_nbi,istep
        energy_fractions: do neut_type=1,3
            !! (type = 1: full energy, =2: half energy, =3: third energy
            if(nbi%current_fractions(neut_type).eq.0.d0) cycle energy_fractions

            call mc_nbi(vnbi,neut_type,rnbi,err)
            if(err) cycle energy_fractions

            call track(rnbi,vnbi,tracks,ntrack)
            if(ntrack.eq.0) cycle energy_fractions

            !! Solve collisional radiative model along track
            flux_tot = 0.d0
            states=0.d0
            states(1)=nneutrals*nbi%current_fractions(neut_type)/beam_grid%dv
            loop_along_track: do jj=1,ntrack
                iflux = sum(states)
                ind = tracks(jj)%ind
                call get_plasma(plasma,pos=tracks(jj)%pos)

                call colrad(plasma,beam_mass,vnbi,tracks(jj)%time,states,dens,photons)
                dens = dens/nlaunch
                call store_neutrals(ind,tracks(jj)%pos,vnbi,neut_type,dens)
                tracks(jj)%flux = (iflux - sum(states))/nlaunch
                flux_tot = flux_tot + tracks(jj)%flux*beam_grid%dv

                if(inputs%calc_birth.ge.1) then
                    call store_births(ind,neut_type,tracks(jj)%flux)
                endif

                if((photons.gt.0.d0).and.(inputs%calc_bes.ge.1)) then
                    call store_nbi_photons(tracks(jj)%pos,vnbi,beam_lambda0, photons/nlaunch,neut_type)
                endif
            enddo loop_along_track
            if((inputs%calc_birth.ge.1).and.(flux_tot.gt.0.d0)) then
                !! Sample according to deposited flux along neutral trajectory
                !$OMP CRITICAL(ndmc_birth)
                do kk=1,inputs%n_birth
                    call randind(tracks(1:ntrack)%flux,randi)
                    call randu(randomu)
                    birth%part(birth%cnt)%neut_type = neut_type
                    birth%part(birth%cnt)%energy = nbi%einj/real(neut_type)
                    birth%part(birth%cnt)%weight = flux_tot/inputs%n_birth
                    birth%part(birth%cnt)%ind = tracks(randi(1))%ind
                    birth%part(birth%cnt)%vi = vnbi
                    ri = tracks(randi(1))%pos + vnbi*(tracks(randi(1))%time*(randomu(1)-0.5))
                    birth%part(birth%cnt)%ri = ri

                    call get_fields(fields,pos=ri)
                    birth%part(birth%cnt)%pitch = dot_product(fields%b_norm,vnbi/norm2(vnbi))
                    call gyro_step(vnbi,fields,beam_mass,r_gyro)
                    birth%part(birth%cnt)%ri_gc = ri + r_gyro
                    birth%cnt = birth%cnt + 1
                enddo
                !$OMP END CRITICAL(ndmc_birth)
            endif
        enddo energy_fractions
    enddo loop_over_markers
    !$OMP END PARALLEL DO

#ifdef _MPI
    !! Combine beam neutrals
    call parallel_merge_populations(neut%full)
    call parallel_merge_populations(neut%half)
    call parallel_merge_populations(neut%third)

    call parallel_sum(nbi_outside)
    if(inputs%calc_birth.ge.1) then
        call parallel_sum(birth%dens)
    endif
    !! Combine spectra
    if(inputs%calc_bes.ge.1) then
        call parallel_sum(spec%full)
        call parallel_sum(spec%half)
        call parallel_sum(spec%third)
    endif
#endif

    if(nbi_outside.gt.0)then
        if(inputs%verbose.ge.1) then
            write(*,'(T4,a, f6.2,a)') 'Percent of markers outside the grid: ', &
                                  100.*nbi_outside/(3.*inputs%n_nbi),'%'
        endif
        if(sum(neut%full%dens).eq.0) stop 'Beam does not intersect the grid!'
    endif

end subroutine ndmc