halo Subroutine

public subroutine halo()

Calculates halo neutral density and spectra

Arguments

None

Calls

proc~~halo~~CallsGraph proc~halo halo proc~mc_beam_grid mc_beam_grid proc~halo->proc~mc_beam_grid proc~mc_halo mc_halo proc~halo->proc~mc_halo proc~free_neutral_population free_neutral_population proc~halo->proc~free_neutral_population proc~init_neutral_population init_neutral_population proc~halo->proc~init_neutral_population proc~ind2sub ind2sub proc~halo->proc~ind2sub proc~get_nlaunch get_nlaunch proc~halo->proc~get_nlaunch proc~get_plasma get_plasma proc~halo->proc~get_plasma proc~store_photons store_photons proc~halo->proc~store_photons proc~merge_neutral_populations merge_neutral_populations proc~halo->proc~merge_neutral_populations proc~parallel_merge_reservoirs parallel_merge_reservoirs proc~halo->proc~parallel_merge_reservoirs proc~store_photon_birth store_photon_birth proc~halo->proc~store_photon_birth proc~track track proc~halo->proc~track proc~colrad colrad proc~halo->proc~colrad proc~update_neutrals update_neutrals proc~halo->proc~update_neutrals proc~parallel_merge_populations parallel_merge_populations proc~halo->proc~parallel_merge_populations interface~parallel_sum parallel_sum proc~halo->interface~parallel_sum proc~neutral_cx_rate neutral_cx_rate proc~halo->proc~neutral_cx_rate interface~randu randu proc~mc_beam_grid->interface~randu interface~randn randn proc~mc_halo->interface~randn proc~get_nlaunch->proc~ind2sub proc~rng_init rng_init proc~get_nlaunch->proc~rng_init interface~randind_cdf randind_cdf proc~get_nlaunch->interface~randind_cdf proc~cumsum cumsum proc~get_nlaunch->proc~cumsum 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 proc~get_passive_grid_indices get_passive_grid_indices proc~store_photons->proc~get_passive_grid_indices proc~store_photons->proc~uvw_to_xyz proc~spectrum spectrum proc~store_photons->proc~spectrum proc~get_indices get_indices proc~store_photons->proc~get_indices proc~get_fields get_fields proc~store_photons->proc~get_fields proc~merge_neutral_populations->proc~ind2sub proc~merge_reservoirs merge_reservoirs proc~merge_neutral_populations->proc~merge_reservoirs proc~parallel_merge_reservoirs->interface~parallel_sum proc~init_reservoir init_reservoir proc~parallel_merge_reservoirs->proc~init_reservoir proc~my_rank my_rank proc~parallel_merge_reservoirs->proc~my_rank proc~parallel_merge_reservoirs->proc~merge_reservoirs proc~num_ranks num_ranks proc~parallel_merge_reservoirs->proc~num_ranks proc~store_photon_birth->proc~get_passive_grid_indices proc~update_reservoir update_reservoir proc~store_photon_birth->proc~update_reservoir proc~store_photon_birth->proc~get_indices proc~track->proc~get_plasma proc~doppler_stark doppler_stark proc~track->proc~doppler_stark proc~track->proc~in_plasma proc~track->proc~get_indices proc~track->proc~get_fields 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~update_neutrals->proc~update_reservoir proc~parallel_merge_populations->proc~ind2sub proc~parallel_merge_populations->proc~parallel_merge_reservoirs proc~parallel_merge_populations->interface~parallel_sum 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~bb_cx_rates bb_cx_rates proc~neutral_cx_rate->proc~bb_cx_rates proc~xyz_to_cyl xyz_to_cyl proc~get_passive_grid_indices->proc~xyz_to_cyl proc~parallel_sum_i3->proc~parallel_sum_i1 mpi_allreduce mpi_allreduce proc~parallel_sum_i3->mpi_allreduce proc~rng_init->proc~my_rank proc~rng_seed rng_seed proc~rng_init->proc~rng_seed proc~init_reservoir->interface~randu 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~merge_reservoirs->proc~init_reservoir proc~merge_reservoirs->proc~update_reservoir proc~merge_reservoirs->interface~randu float float proc~merge_reservoirs->float proc~parallel_sum_i2->proc~parallel_sum_i1 proc~parallel_sum_i2->mpi_allreduce proc~update_reservoir->proc~init_reservoir proc~update_reservoir->interface~randu interface~randind randind proc~update_reservoir->interface~randind interface~interpol_coeff interpol_coeff proc~get_rate_matrix->interface~interpol_coeff proc~elmhes elmhes proc~eigen->proc~elmhes proc~balance balance proc~eigen->proc~balance proc~hqr2 hqr2 proc~eigen->proc~hqr2 proc~balback balback proc~eigen->proc~balback proc~elmtrans elmtrans proc~eigen->proc~elmtrans 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 dgetrs dgetrs proc~linsolve->dgetrs proc~matinv matinv proc~linsolve->proc~matinv dgetrf dgetrf proc~linsolve->dgetrf 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~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~calc_perp_vectors calc_perp_vectors proc~get_fields->proc~calc_perp_vectors proc~parallel_sum_d2->proc~parallel_sum_d1 proc~parallel_sum_d2->mpi_allreduce proc~bb_cx_rates->interface~interpol_coeff proc~rswap RSWAP proc~elmhes->proc~rswap proc~xyz_to_cyl->proc~xyz_to_uvw proc~uvw_to_cyl uvw_to_cyl proc~xyz_to_cyl->proc~uvw_to_cyl proc~balance->proc~rswap proc~hqrvec hqrvec proc~hqr2->proc~hqrvec proc~lubksb lubksb proc~matinv->proc~lubksb proc~ludcmp ludcmp proc~matinv->proc~ludcmp proc~balback->proc~rswap proc~cyl_to_xyz->proc~uvw_to_xyz proc~cyl_to_xyz->proc~cyl_to_uvw proc~comdiv Comdiv proc~hqrvec->proc~comdiv proc~outerprod outerprod proc~ludcmp->proc~outerprod proc~swap swap proc~ludcmp->proc~swap

Called by

proc~~halo~~CalledByGraph proc~halo halo program~fidasim fidasim program~fidasim->proc~halo

Contents

Source Code


Source Code

subroutine halo
    !+ Calculates halo neutral density and spectra
    integer :: ic,i,j,k,ncell
    integer(Int64) :: ihalo !! counter
    real(Float64), dimension(3) :: ri    !! start position
    real(Float64), dimension(3) :: vihalo!! velocity bulk plasma ion
    integer,dimension(3) :: ind, tind    !! actual cell
    !! Determination of the CX probability
    type(LocalProfiles) :: plasma
    real(Float64), dimension(nlevs) :: denn    !!  neutral dens (n=1-4)
    real(Float64), dimension(nlevs) :: rates    !! CX rates
    !! Collisiional radiative model along track
    real(Float64), dimension(nlevs) :: states  ! Density of n-states
    integer :: ntrack
    type(ParticleTrack), dimension(beam_grid%ntrack) :: tracks  !! Particle Tracks
    integer :: ii,jj,kk,it,is
    real(Float64) :: tot_denn, photons  !! photon flux
    integer, dimension(beam_grid%ngrid) :: cell_ind
    real(Float64), dimension(beam_grid%nx,beam_grid%ny,beam_grid%nz) :: papprox
    integer(Int32), dimension(beam_grid%nx,beam_grid%ny,beam_grid%nz) :: nlaunch
    real(Float64) :: local_iter_dens

    type(NeutralPopulation) :: cur_pop, prev_pop
    !! Halo iteration
    integer(Int64) :: hh, n_halo !! counters
    real(Float64) :: dcx_dens, halo_iter_dens(2), seed_dcx
    integer :: prev_type = 1  ! previous iteration
    integer :: cur_type = 2 ! current iteration
    real(Float64) :: fi_correction

    !! Initialize Neutral Population
    call init_neutral_population(neut%halo)

    dcx_dens = sum(neut%dcx%dens)

    halo_iter_dens(prev_type) = dcx_dens
    if(dcx_dens.eq.0) then
        if(inputs%verbose.ge.0) then
            write(*,'(a)') 'HALO: Density of DCX-neutrals is zero'
        endif
        stop
    endif
    n_halo = inputs%n_halo

    !! Allocate previous neutral populations
    call init_neutral_population(prev_pop)
    prev_pop = neut%dcx

    seed_dcx = 1.0
    iterations: do hh=1,200
        !! Allocate/Reallocate current population
        call init_neutral_population(cur_pop)
        halo_iter_dens(cur_type) = 0.d0

        !! Calculate how many mc markers to launch in each cell
        papprox = 0.d0
        cell_ind = 0
        ncell = 0
        do ic=1,beam_grid%ngrid
            call ind2sub(beam_grid%dims,ic,ind)
            i = ind(1) ; j = ind(2) ; k = ind(3)

            call get_plasma(plasma,ind=ind)
            if(.not.plasma%in_plasma) cycle
            tot_denn = sum(prev_pop%dens(:,i,j,k))
            papprox(i,j,k)= tot_denn*sum(plasma%deni)
            if(papprox(i,j,k).gt.0.0) then
                ncell = ncell + 1
                cell_ind(ncell) = ic
            endif
        enddo
        call get_nlaunch(n_halo, papprox, nlaunch)

        if(inputs%verbose.ge.1) then
            write(*,'(T6,"# of markers: ",i10," --- Seed/DCX: ",f5.3)') sum(nlaunch), seed_dcx
        endif

        !$OMP PARALLEL DO schedule(dynamic,1) private(i,j,k,ic,ihalo,ii,jj,kk,it,is,ind,vihalo, &
        !$OMP& ri,tracks,ntrack,rates,denn,states,photons,plasma,tind,fi_correction)
        loop_over_cells: do ic=istart,ncell,istep
            call ind2sub(beam_grid%dims,cell_ind(ic),ind)
            i = ind(1) ; j = ind(2) ; k = ind(3)
            loop_over_species: do is=1,n_thermal !This loop has to come first
                !! Loop over the markers
                loop_over_halos: do ihalo=1, nlaunch(i,j,k)
                    !! Calculate ri,vhalo and track
                    call mc_beam_grid(ind, ri)
                    call get_plasma(plasma, pos=ri)
                    call mc_halo(plasma, thermal_mass(is), vihalo)
                    call track(ri,vihalo,tracks,ntrack)
                    if(ntrack.eq.0)cycle loop_over_halos

                    !! Get plasma parameters at particle location
                    call get_plasma(plasma, pos=ri)

                    !! Calculate CX probability
                    tind = tracks(1)%ind
                    ii = tind(1); jj = tind(2); kk = tind(3)
                    call neutral_cx_rate(prev_pop%dens(:,ii,jj,kk), prev_pop%res(ii,jj,kk), vihalo, rates)
                    if(sum(rates).le.0.)cycle loop_over_halos

                    !! Get plasma parameters at mean point in cell
                    call get_plasma(plasma, pos=tracks(1)%pos)

                    !! Weight CX rates by ion source density
                    if(beam_mass.eq.thermal_mass(is)) then
                        states = rates*(plasma%deni(is) + plasma%denf)
                        if(sum(states).eq.0) cycle loop_over_halos
                        fi_correction = max(plasma%deni(is)/(plasma%deni(is)+plasma%denf),0.d0)
                    else
                        states = rates*plasma%deni(is)
                        if(sum(states).eq.0) cycle loop_over_halos
                        fi_correction = 1.d0
                    endif

                    loop_along_track: do it=1,ntrack
                        call get_plasma(plasma,pos=tracks(it)%pos)
                        if(.not.plasma%in_plasma) exit loop_along_track
                        call colrad(plasma,thermal_mass(is),vihalo,tracks(it)%time,states,denn,photons)

                        !! Store Neutrals
                        call update_neutrals(cur_pop, tracks(it)%ind, vihalo, denn/nlaunch(i,j,k))

                        photons = fi_correction*photons !! Correct for including fast-ions in states

                        if((photons.gt.0.d0).and.(inputs%calc_halo.ge.1)) then
                            call store_photons(tracks(it)%pos,vihalo,thermal_lambda0(is), &
                                               photons/nlaunch(i,j,k),spec%halo(:,:,:,is),spec%halostokes(:,:,:,:,is))
                        endif

                        if((photons.gt.0.d0).and.(inputs%calc_res.ge.1)) then
                            call store_photon_birth(tracks(1)%pos, photons/nlaunch(i,j,k), spatres%halo)
                        endif
                    enddo loop_along_track
                enddo loop_over_halos
            enddo loop_over_species
        enddo loop_over_cells
        !$OMP END PARALLEL DO

#ifdef _MPI
        !! Combine densities
        call parallel_merge_populations(cur_pop)
#endif
        halo_iter_dens(cur_type) = sum(cur_pop%dens)

        if(halo_iter_dens(cur_type)/halo_iter_dens(prev_type).gt.1.0) then
            write(*,'(a)') "HALO: Halo generation density exceeded seed density. This shouldn't happen."
            exit iterations
        endif

        !! Set current generation to previous generation
        halo_iter_dens(prev_type) = halo_iter_dens(cur_type)
        prev_pop = cur_pop

        !! merge current population with halo population
        call merge_neutral_populations(neut%halo, cur_pop)

        seed_dcx = halo_iter_dens(cur_type)/dcx_dens
        if(isnan(seed_dcx)) then
            write(*,'(a)') "NAN (Halo generation)/DCX. Exiting..."
            stop
        endif
        n_halo=int(inputs%n_halo*seed_dcx,Int64)

        if(seed_dcx.lt.0.01) exit iterations
    enddo iterations

    call free_neutral_population(cur_pop)
    call free_neutral_population(prev_pop)

#ifdef _MPI
    !! Combine Spectra
    if(inputs%calc_halo.ge.1) then
        call parallel_sum(spec%halo)
    endif
    if(inputs%calc_res.ge.1) then
        do jj=1,spec_chords%nchan
            call parallel_merge_reservoirs(spatres%halo(jj)
        enddo
    endif
#endif

end subroutine halo