fusion_sources_module.f90

Go to the documentation of this file.

module fusion_sources_module

  use amns_types
  integer, parameter :: nreac=4, nspec=6
  type (amns_handle_type), save :: amns                         ! AMNS global handle
  type (amns_handle_rx_type), save  :: nuclear(nreac)           ! AMNS table handle

contains

  subroutine fusion_sources(coreprof, coresource)

    use itm_types
    use amns_module
    use copy_structures
    use ets_version
    use itm_constants
!    use coresource_types, only: t_fusion=>fusion, get_type_name, get_type_description__ind    ! itmconstants/1.0
    use coresource_identifier, only: t_fusion=>fusion, get_type_name, get_type_description__ind    ! itmconstants/1.3

    implicit none

    type (type_coreprof), pointer :: coreprof(:)
    type (type_coresource), pointer :: coresource(:)

    type (amns_reaction_type)  :: xx_rx
    type (amns_query_type)     :: query
    type (amns_answer_type)    :: answer
    type (amns_set_type)       :: set
    type (amns_reactants_type) :: species
    type (amns_reactant_type)  :: nuclei(nspec)
    real (kind=R8), allocatable :: energy(:), density(:,:), source(:,:), heating(:), rate(:)
    integer :: i, j, nr, ns, nucindex, im, iz
    integer, save :: ions_index(nspec)  ! p, D, T, 3He, 4He, n
    real (kind=R8), save :: mass(nspec) =  &
         (/ itm_mass_h_1, itm_mass_h_2, itm_mass_h_3, itm_mass_he_3, itm_mass_he_4, itm_mn /)
    integer (kind=R8), save :: rm(4,nreac)
    real (kind=R8), save :: fusion_energy(nreac)
    real (kind=R8) :: fraction
    logical, save :: first=.true.
    logical :: found

    if(first) then
       call itm_amns_setup(amns)                                  ! set up the AMNS system
       query%string  = 'version'
       call itm_amns_query(amns,query,answer)
       write(*,*)  'fusion: amns data base version = ',trim(answer%string)

       nuclei(1) = amns_reactant_type(1, 0, 1, 0)
       nuclei(2) = amns_reactant_type(1, 0, 2, 0)
       nuclei(3) = amns_reactant_type(1, 0, 3, 0)
       nuclei(4) = amns_reactant_type(2, 0, 3, 0)
       nuclei(5) = amns_reactant_type(2, 0, 4, 0)
       nuclei(6) = amns_reactant_type(0, 0, 1, 0)

       rm(:,1) = (/ 2, 2, 1, 3 /)   ! D_D_p_T
       rm(:,2) = (/ 2, 2, 6, 4 /)   ! D_D_n_He
       rm(:,3) = (/ 2, 3, 6, 5 /)   ! D_T_n_He
       rm(:,4) = (/ 2, 4, 1, 5 /)   ! D_He_p_He

! sanity check for the reactions
       do nr = 1, nreac
          im = nint(nuclei(rm(1,nr))%mi) + nint(nuclei(rm(2,nr))%mi) - &
               nint(nuclei(rm(3,nr))%mi) - nint(nuclei(rm(4,nr))%mi)
          if(im.ne.0) then
             write(*,*) 'fusion: baryons not conserved for reaction ', nr
             stop 'Error'
          endif
          iz = nint(nuclei(rm(1,nr))%zn) + nint(nuclei(rm(2,nr))%zn) - &
               nint(nuclei(rm(3,nr))%zn) - nint(nuclei(rm(4,nr))%zn)
          if(iz.ne.0) then
             write(*,*) 'fusion: protons not conserved for reaction ', iz
             stop 'Error'
          endif
       enddo
 
! calculate the information we will need later for the reactions
       allocate(species%components(4))                            ! set up reactants
       xx_rx%string='NUC_TT'                                      ! set up reaction
       xx_rx%isotope_resolved='1'

       do nr = 1, nreac
          species%components = (/ &
               nuclei(rm(1,nr)), &
               nuclei(rm(2,nr)), &
               nuclei(rm(3,nr)), &
               nuclei(rm(4,nr)) /)
          species%components(3)%LR=1                              ! product
          species%components(4)%LR=1                              ! product
          if(nuclei(rm(3,nr))%zn .eq.0.0_r8) then
             fraction = mass(rm(3,nr)) / (mass(rm(3,nr))+mass(rm(4,nr)))
          elseif(nuclei(rm(4,nr))%zn .eq.0.0_r8) then
             fraction = mass(rm(4,nr)) / (mass(rm(3,nr))+mass(rm(4,nr)))
          else
             fraction = 1.0_r8
          endif
          fusion_energy(nr) = &
               (mass(rm(1,nr)) + mass(rm(2,nr)) - mass(rm(3,nr)) - mass(rm(4,nr))) * itm_c**2 * &
               fraction
          call itm_amns_setup_table(amns, xx_rx, species, nuclear(nr))   ! set up table
          query%string = 'reactants'
          call itm_amns_query_table(nuclear(nr), query, answer)
          write(*,*) 'fusion: ', nr, trim(answer%string), fusion_energy(nr)/itm_qe/1.0e6
       enddo

       deallocate(species%components)

! calculate the mapping vector from external to internal species
       found = .false.
       ions_index = 0
       do i = 1, size(coreprof(1)%ni%value, dim=2)
          nucindex = coreprof(1)%compositions%ions(i)%nucindex
          do ns = 1, nspec
             found = &
                  (nint(coreprof(1)%compositions%nuclei(nucindex)%zn) .eq. nuclei(ns)%zn) .and. &
                  (nint(coreprof(1)%compositions%nuclei(nucindex)%amn) .eq. nuclei(ns)%mi)
             if(found) then
                ions_index(ns) = i
             endif
          enddo
       enddo
       write(*,*) 'fusion: ions_index = ', ions_index

       first = .false.

    endif

! arrays we will need
    allocate(energy(size(coreprof(1)%ni%value, dim=1)))
    allocate(density(size(energy),nspec))
    allocate(source(size(energy),nspec))
    allocate(heating(size(energy)))
    allocate(rate(size(energy)))

! energy of the reactants
    energy = sum(coreprof(1)%ti%value * coreprof(1)%ni%value, 2) / sum(coreprof(1)%ni%value,2)

! densities of reactants
    density = 0
    do ns = 1, nspec
       if(ions_index(ns) .ne. 0) then
          density(:,ns) = coreprof(1)%ni%value(:,ions_index(ns))
       endif
    enddo

! calculate sources
    source = 0
    heating = 0
    do nr = 1, nreac
       call itm_amns_rx(nuclear(nr), rate, energy)              ! get results
       rate(:) = rate(:) * density(:,rm(1,nr)) * density(:,rm(2,nr))
       if(rm(1,nr) .eq. rm(2,nr)) rate = rate * 0.5             ! this was missed earlier
       source(:,rm(1,nr)) = source(:,rm(1,nr)) - rate(:)
       source(:,rm(2,nr)) = source(:,rm(2,nr)) - rate(:)
       source(:,rm(3,nr)) = source(:,rm(3,nr)) + rate(:)
       source(:,rm(4,nr)) = source(:,rm(4,nr)) + rate(:)
       heating(:) = heating + rate(:) * fusion_energy(nr)
    enddo

    write(*,*) 'fusion: density = ', density(1,:)
    write(*,*) 'fusion: source = ', source(1,:)
    write(*,*) 'fusion: heating = ', heating(1)

! prepare the output CPO
    allocate(coresource(1))
    allocate(coresource(1)%codeparam%codename(1))
    allocate(coresource(1)%codeparam%codeversion(1))
    allocate(coresource(1)%codeparam%output_diag(1))
    coresource(1)%codeparam%codename = 'fusion'
    coresource(1)%codeparam%codeversion = version
    coresource(1)%codeparam%output_diag = 'OK'
    coresource(1)%codeparam%output_flag = 0
    coresource(1)%time = coreprof(1)%time
    call copy_cpo(coreprof(1)%compositions, coresource(1)%compositions)
    allocate(coresource(1)%values(1))
    allocate(coresource(1)%values(1)%sourceid%id(1))
    allocate(coresource(1)%values(1)%sourceid%description(1))
    coresource(1)%values(1)%sourceid%flag = t_fusion
    coresource(1)%values(1)%sourceid%id = get_type_name(t_fusion)
    coresource(1)%values(1)%sourceid%description = get_type_description__ind(t_fusion)
    call copy_cpo(coreprof(1)%rho_tor, coresource(1)%values(1)%rho_tor)
    allocate(coresource(1)%values(1)%si%exp(size(coreprof(1)%ni%value,dim=1),size(coreprof(1)%ni%value,dim=2)))
    coresource(1)%values(1)%si%exp=0
    do ns = 1, nspec
       if(ions_index(ns) .ne. 0) then
          coresource(1)%values(1)%si%exp(:,ions_index(ns)) = source(:,ns)
       endif
    enddo
    allocate(coresource(1)%values(1)%qe%exp(size(coreprof(1)%ni%value,dim=1)))
    coresource(1)%values(1)%qe%exp(:)=heating(:)

! deallocate no longer needed internal arrays
    deallocate(energy, density, source, heating, rate)

    return

  end subroutine fusion_sources

  subroutine fusion_finalize()

    use amns_module

    implicit none

    integer :: nr

    do nr = 1, nreac
       call itm_amns_finish_table(nuclear(nr))                 ! finish with table
    enddo
    call itm_amns_finish(amns)                                 ! finish with amns

  end subroutine fusion_finalize

end module fusion_sources_module