diff --git a/src/constants.F90 b/src/constants.F90 index 853a40c8d4..936c84499a 100644 --- a/src/constants.F90 +++ b/src/constants.F90 @@ -308,7 +308,7 @@ module constants ! Tally score type -- if you change these, make sure you also update the ! _SCORES dictionary in openmc/capi/tally.py - integer, parameter :: N_SCORE_TYPES = 24 + integer, parameter :: N_SCORE_TYPES = 30 integer, parameter :: & SCORE_FLUX = -1, & ! flux SCORE_TOTAL = -2, & ! total reaction rate @@ -333,8 +333,13 @@ module constants SCORE_INVERSE_VELOCITY = -21, & ! flux-weighted inverse velocity SCORE_FISS_Q_PROMPT = -22, & ! prompt fission Q-value SCORE_FISS_Q_RECOV = -23, & ! recoverable fission Q-value - SCORE_DECAY_RATE = -24 ! delayed neutron precursor decay rate - + SCORE_DECAY_RATE = -24, & ! delayed neutron precursor decay rate + SCORE_N2N = -25, & + SCORE_NGAMMA = -26, & + SCORE_N3N = -27, & + SCORE_N4N = -28, & + SCORE_NP = -29, & + SCORE_NALPHA = -30 ! Maximum scattering order supported integer, parameter :: MAX_ANG_ORDER = 10 diff --git a/src/cross_section.F90 b/src/cross_section.F90 index 46cbb93abc..2484be51b9 100644 --- a/src/cross_section.F90 +++ b/src/cross_section.F90 @@ -2,6 +2,7 @@ module cross_section use algorithm, only: binary_search use constants + use endf use error, only: fatal_error use list_header, only: ListElemInt use material_header, only: Material, materials @@ -45,7 +46,12 @@ contains material_xs % absorption = ZERO material_xs % fission = ZERO material_xs % nu_fission = ZERO - + material_xs % n2n = ZERO + material_xs % n3n = ZERO + material_xs % n4n = ZERO + material_xs % ngamma = ZERO + material_xs % np = ZERO + material_xs % nalpha = ZERO ! Exit subroutine if material is void if (p % material == MATERIAL_VOID) return @@ -129,6 +135,35 @@ contains ! Add contributions to material macroscopic nu-fission cross section material_xs % nu_fission = material_xs % nu_fission + & atom_density * micro_xs(i_nuclide) % nu_fission + + + ! Add contributions to material macroscopic n2n cross section + material_xs % n2n = material_xs % n2n + & + atom_density * micro_xs(i_nuclide) % n2n + + + ! Add contributions to material macroscopic ngamma cross section + material_xs % ngamma = material_xs % ngamma + & + atom_density * micro_xs(i_nuclide) % ngamma + + ! Add contributions to material macroscopic n3n cross section + material_xs % n3n = material_xs % n3n + & + atom_density * micro_xs(i_nuclide) % n3n + + ! Add contributions to material macroscopic n4n cross section + material_xs % n4n = material_xs % n4n + & + atom_density * micro_xs(i_nuclide) % n4n + + ! Add contributions to material macroscopic np cross section + material_xs % np = material_xs % np + & + atom_density * micro_xs(i_nuclide) % np + + + ! Add contributions to material macroscopic nalpha cross section + material_xs % nalpha = material_xs % nalpha + & + atom_density * micro_xs(i_nuclide) % nalpha + + end do end associate @@ -157,7 +192,9 @@ contains real(8) :: f ! interp factor on nuclide energy grid real(8) :: kT ! temperature in eV real(8) :: sigT, sigA, sigF ! Intermediate multipole variables - + integer :: rxn + type(TallyObject) :: t + character(15) :: threshold_rxn_type associate (nuc => nuclides(i_nuclide)) ! Check to see if there is multipole data present at this energy use_mp = .false. @@ -202,6 +239,8 @@ contains else ! Find the appropriate temperature index. kT = sqrtkT**2 + + !rxn = nuc % rxn_index_MT(p % event_MT) select case (temperature_method) case (TEMPERATURE_NEAREST) i_temp = minloc(abs(nuclides(i_nuclide) % kTs - kT), dim=1) @@ -247,12 +286,19 @@ contains (grid % energy(i_grid + 1) - grid % energy(i_grid)) micro_xs(i_nuclide) % index_temp = i_temp + micro_xs(i_nuclide) % index_grid = i_grid + micro_xs(i_nuclide) % interp_factor = f ! Initialize nuclide cross-sections to zero micro_xs(i_nuclide) % fission = ZERO micro_xs(i_nuclide) % nu_fission = ZERO + micro_xs(i_nuclide) % n2n = ZERO + micro_xs(i_nuclide) % n3n = ZERO + micro_xs(i_nuclide) % n4n = ZERO + micro_xs(i_nuclide) % np = ZERO + micro_xs(i_nuclide) % nalpha = ZERO micro_xs(i_nuclide) % thermal = ZERO micro_xs(i_nuclide) % thermal_elastic = ZERO @@ -274,11 +320,66 @@ contains + f * xs % fission(i_grid + 1) ! Calculate microscopic nuclide nu-fission cross section - micro_xs(i_nuclide) % nu_fission = (ONE - f) * xs % nu_fission( & + micro_xs(i_nuclide) % nu_fission = (ONE - f) * xs % nu_fission(& i_grid) + f * xs % nu_fission(i_grid + 1) + end if - end associate - end if + + micro_xs(i_nuclide) % ngamma = (micro_xs(i_nuclide) % absorption - & + micro_xs(i_nuclide) % fission) + if (t % has_threshold_rxn) then + !if (nuclides(i_nuclide)%reaction_index%has_key(p % event_MT)) the + rxn = nuclides(i_nuclide) % rxn_index_MT(p % event_MT) + associate(xs_threshold => nuc % reactions(rxn) % xs(i_temp)) + !end if + threshold_rxn_type = reaction_name(p % event_MT) + if (threshold_rxn_type == 'n2n') then + if (i_grid >= xs_threshold % threshold) then + micro_xs(i_nuclide) % n2n = (ONE - f) * xs_threshold % value(i_grid - xs_threshold & + % threshold + 1) + f * xs_threshold % value(i_grid - xs_threshold & + % threshold + 2) + end if + end if + if (threshold_rxn_type == 'n3n') then + if (i_grid >= xs_threshold % threshold) then + micro_xs(i_nuclide) % n3n = (ONE - f) * xs_threshold % value(i_grid - xs_threshold & + % threshold + 1) + f * xs_threshold % value(i_grid - xs_threshold & + % threshold + 2) + end if + end if + + + if (threshold_rxn_type == 'n4n') then + if (i_grid >= xs_threshold % threshold) then + micro_xs(i_nuclide) % n4n = (ONE - f) * xs_threshold % value(i_grid - xs_threshold & + % threshold + 1) + f * xs_threshold % value(i_grid - xs_threshold & + % threshold + 2) + end if + end if + + if (threshold_rxn_type == 'np') then + if (i_grid >= xs_threshold % threshold) then + micro_xs(i_nuclide) % np = (ONE - f) * xs_threshold % value(i_grid - xs_threshold & + % threshold + 1) + f * xs_threshold % value(i_grid - xs_threshold & + % threshold + 2) + end if + end if + + if (threshold_rxn_type == 'nalpha') then + if (i_grid >= xs_threshold % threshold) then + micro_xs(i_nuclide) % nalpha = (ONE - f) * xs_threshold % value(i_grid - xs_threshold & + % threshold + 1) + f * xs_threshold % value(i_grid - xs_threshold & + % threshold + 2) + end if + end if + + + + end associate + + end if + end associate + end if ! Initialize sab treatment to false micro_xs(i_nuclide) % index_sab = NONE diff --git a/src/endf.F90 b/src/endf.F90 index 0ba2db388b..9aace42d7b 100644 --- a/src/endf.F90 +++ b/src/endf.F90 @@ -66,6 +66,18 @@ contains string = "fission-q-prompt" case (SCORE_FISS_Q_RECOV) string = "fission-q-recoverable" + case (SCORE_N2N) + string = "n2n" + case (SCORE_N3N) + string = "n3n" + case (SCORE_N4N) + string = "n4n" + case (SCORE_NGAMMA) + string = "ngamma" + case (SCORE_NP) + string = "np" + case (SCORE_NALPHA) + string = "nalpha" ! Normal ENDF-based reactions case (TOTAL_XS) @@ -76,12 +88,12 @@ contains string = '(n,level)' case (N_2ND) string = '(n,2nd)' - case (N_2N) - string = '(n,2n)' - case (N_3N) - string = '(n,3n)' - case (N_FISSION) - string = '(n,fission)' + !case (N_2N) + !string = '(n,2n)' + !case (N_3N) + !string = '(n,3n)' + !case (N_FISSION) + ! string = '(n,fission)' case (N_F) string = '(n,f)' case (N_NF) @@ -130,18 +142,18 @@ contains string = '(n,nc)' case (N_DISAPPEAR) string = '(n,disappear)' - case (N_GAMMA) - string = '(n,gamma)' - case (N_P) - string = '(n,p)' + !case (N_GAMMA) + ! string = '(n,gamma)' + !case (N_P) + !string = '(n,p)' case (N_D) string = '(n,d)' case (N_T) string = '(n,t)' case (N_3HE) string = '(n,3He)' - case (N_A) - string = '(n,a)' + !case (N_A) + !string = '(n,a)' case (N_2A) string = '(n,2a)' case (N_3A) diff --git a/src/input_xml.F90 b/src/input_xml.F90 index 91cd57ae29..4587072f5d 100644 --- a/src/input_xml.F90 +++ b/src/input_xml.F90 @@ -2450,6 +2450,17 @@ contains n_nuclides = index_nuclide n_sab_tables = index_sab + + do i=1, n_materials + mat => materials(i) + allocate(mat % mat_nuclide_list(n_nuclides_total)) + mat % mat_nuclide_list(:) = 0 + do j=1, mat % n_nuclides + mat % mat_nuclide_list(mat % nuclide(j))=j + end do + end do + + ! Close materials XML file call doc % clear() @@ -3047,14 +3058,14 @@ contains call fatal_error("n1n score no longer supported for tallies, & &please remove") case ('n2n', '(n,2n)') - t % score_bins(j) = N_2N - + t % score_bins(j) = SCORE_N2N + t % has_threshold_rxn = .true. case ('n3n', '(n,3n)') - t % score_bins(j) = N_3N - + t % score_bins(j) = SCORE_N3N + t % has_threshold_rxn = .true. case ('n4n', '(n,4n)') - t % score_bins(j) = N_4N - + t % score_bins(j) = SCORE_N4N + t % has_threshold_rxn = .true. case ('absorption') t % score_bins(j) = SCORE_ABSORPTION if (t % find_filter(FILTER_ENERGYOUT) > 0) then @@ -3238,9 +3249,12 @@ contains case ('(n,nc)') t % score_bins(j) = N_NC case ('(n,gamma)') - t % score_bins(j) = N_GAMMA + t % score_bins(j) = SCORE_NGAMMA case ('(n,p)') - t % score_bins(j) = N_P + t % score_bins(j) = SCORE_NP + + t % has_threshold_rxn = .true. + case ('(n,d)') t % score_bins(j) = N_D case ('(n,t)') @@ -3248,7 +3262,8 @@ contains case ('(n,3He)') t % score_bins(j) = N_3HE case ('(n,a)') - t % score_bins(j) = N_A + t % score_bins(j) = SCORE_NALPHA + t % has_threshold_rxn = .true. case ('(n,2a)') t % score_bins(j) = N_2A case ('(n,3a)') diff --git a/src/nuclide_header.F90 b/src/nuclide_header.F90 index 2c5de1b140..1f6d45165d 100644 --- a/src/nuclide_header.F90 +++ b/src/nuclide_header.F90 @@ -43,6 +43,14 @@ module nuclide_header real(8), allocatable :: nu_fission(:) ! neutron production real(8), allocatable :: absorption(:) ! absorption (MT > 100) real(8), allocatable :: heating(:) ! heating + real(8), allocatable :: ngamma(:) + real(8), allocatable :: n2n(:) + real(8), allocatable :: n3n(:) + real(8), allocatable :: n4n(:) + + real(8), allocatable :: np(:) + real(8), allocatable :: nalpha(:) + end type SumXS type :: Nuclide @@ -118,6 +126,12 @@ module nuclide_header real(8) :: absorption real(8) :: fission real(8) :: nu_fission + real(8) :: n2n + real(8) :: ngamma + real(8) :: n3n + real(8) :: n4n + real(8) :: np + real(8) :: nalpha real(8) :: thermal ! Bound thermal elastic & inelastic scattering real(8) :: thermal_elastic ! Bound thermal elastic scattering @@ -148,6 +162,13 @@ module nuclide_header real(8) :: absorption ! macroscopic absorption xs real(8) :: fission ! macroscopic fission xs real(8) :: nu_fission ! macroscopic production xs + real(8) :: n2n + real(8) :: n3n + real(8) :: n4n + real(8) :: ngamma + real(8) :: np + real(8) :: nalpha + end type MaterialMacroXS !=============================================================================== diff --git a/src/tallies/tally.F90 b/src/tallies/tally.F90 index c47f94aa76..811ead4f8b 100644 --- a/src/tallies/tally.F90 +++ b/src/tallies/tally.F90 @@ -1131,6 +1131,234 @@ contains end if end if + + case (SCORE_N2N) + + score = ZERO + + if (i_nuclide > 0) then + !if (nuclides(i_nuclide)%reaction_index%has_key(score_bin)) then + !m = nuclides(i_nuclide)%reaction_index%get_key(score_bin) + + ! Retrieve temperature and energy grid index and interpolation + ! factor + !i_temp = micro_xs(i_nuclide) % index_temp + !if (i_temp > 0) then + !i_energy = micro_xs(i_nuclide) % index_grid + !f = micro_xs(i_nuclide) % interp_factor + + !associate (xs => nuclides(i_nuclide) % reactions(m) % xs(i_temp)) + !if (i_energy >= xs % threshold) then + score = micro_xs(i_nuclide) % n2n * atom_density * flux + + + !else + ! This block is reached if multipole is turned on and we're in + ! the resolved range. For (n,gamma), use absorption - + ! fission. For everything else, assume it's zero. + + ! score = ZERO + + else + do l = 1, materials(p % material) % n_nuclides + ! Get atom density + atom_density_ = materials(p % material) % atom_density(l) + + ! Get index in nuclides array + i_nuc = materials(p % material) % nuclide(l) + + score = score + micro_xs(i_nuc) % n2n * atom_density_ * flux + end do + end if + + + case (SCORE_N3N) + + score = ZERO + + if (i_nuclide > 0) then + !if (nuclides(i_nuclide)%reaction_index%has_key(score_bin)) then + !m = nuclides(i_nuclide)%reaction_index%get_key(score_bin) + + ! Retrieve temperature and energy grid index and interpolation + ! factor + !i_temp = micro_xs(i_nuclide) % index_temp + !if (i_temp > 0) then + !i_energy = micro_xs(i_nuclide) % index_grid + !f = micro_xs(i_nuclide) % interp_factor + + !associate (xs => nuclides(i_nuclide) % reactions(m) % + !xs(i_temp)) + !if (i_energy >= xs % threshold) then + score = micro_xs(i_nuclide) % n3n * atom_density * flux + + + !else + ! This block is reached if multipole is turned on and we're in + ! the resolved range. For (n,gamma), use absorption - + ! fission. For everything else, assume it's zero. + + ! score = ZERO + + else + do l = 1, materials(p % material) % n_nuclides + ! Get atom density + atom_density_ = materials(p % material) % atom_density(l) + + ! Get index in nuclides array + i_nuc = materials(p % material) % nuclide(l) + + score = score + micro_xs(i_nuc) % n3n * atom_density_ * flux + end do + end if + + + case (SCORE_N4N) + + score = ZERO + + if (i_nuclide > 0) then + !if (nuclides(i_nuclide)%reaction_index%has_key(score_bin)) then + !m = nuclides(i_nuclide)%reaction_index%get_key(score_bin) + + ! Retrieve temperature and energy grid index and interpolation + ! factor + !i_temp = micro_xs(i_nuclide) % index_temp + !if (i_temp > 0) then + !i_energy = micro_xs(i_nuclide) % index_grid + !f = micro_xs(i_nuclide) % interp_factor + + !associate (xs => nuclides(i_nuclide) % reactions(m) % + !xs(i_temp)) + !if (i_energy >= xs % threshold) then + score = micro_xs(i_nuclide) % n4n * atom_density * flux + + + !else + ! This block is reached if multipole is turned on and we're in + ! the resolved range. For (n,gamma), use absorption - + ! fission. For everything else, assume it's zero. + + ! score = ZERO + + else + do l = 1, materials(p % material) % n_nuclides + ! Get atom density + atom_density_ = materials(p % material) % atom_density(l) + + ! Get index in nuclides array + i_nuc = materials(p % material) % nuclide(l) + + score = score + micro_xs(i_nuc) % n4n * atom_density_ * flux + end do + end if + + case (SCORE_NP) + + score = ZERO + + if (i_nuclide > 0) then + !if (nuclides(i_nuclide)%reaction_index%has_key(score_bin)) then + !m = nuclides(i_nuclide)%reaction_index%get_key(score_bin) + + ! Retrieve temperature and energy grid index and interpolation + ! factor + !i_temp = micro_xs(i_nuclide) % index_temp + !if (i_temp > 0) then + !i_energy = micro_xs(i_nuclide) % index_grid + !f = micro_xs(i_nuclide) % interp_factor + + !associate (xs => nuclides(i_nuclide) % reactions(m) % + !xs(i_temp)) + !if (i_energy >= xs % threshold) then + score = micro_xs(i_nuclide) % np * atom_density * flux + + + !else + ! This block is reached if multipole is turned on and we're in + ! the resolved range. For (n,gamma), use absorption - + ! fission. For everything else, assume it's zero. + + ! score = ZERO + + else + do l = 1, materials(p % material) % n_nuclides + ! Get atom density + atom_density_ = materials(p % material) % atom_density(l) + + ! Get index in nuclides array + i_nuc = materials(p % material) % nuclide(l) + + score = score + micro_xs(i_nuc) % np * atom_density_ * flux + end do + end if + + case (SCORE_NALPHA) + + score = ZERO + + if (i_nuclide > 0) then + !if (nuclides(i_nuclide)%reaction_index%has_key(score_bin)) then + !m = nuclides(i_nuclide)%reaction_index%get_key(score_bin) + + ! Retrieve temperature and energy grid index and interpolation + ! factor + !i_temp = micro_xs(i_nuclide) % index_temp + !if (i_temp > 0) then + !i_energy = micro_xs(i_nuclide) % index_grid + !f = micro_xs(i_nuclide) % interp_factor + + !associate (xs => nuclides(i_nuclide) % reactions(m) % + !xs(i_temp)) + !if (i_energy >= xs % threshold) then + score = micro_xs(i_nuclide) % nalpha * atom_density * flux + + + !else + ! This block is reached if multipole is turned on and we're in + ! the resolved range. For (n,gamma), use absorption - + ! fission. For everything else, assume it's zero. + + ! score = ZERO + + else + do l = 1, materials(p % material) % n_nuclides + ! Get atom density + atom_density_ = materials(p % material) % atom_density(l) + + ! Get index in nuclides array + i_nuc = materials(p % material) % nuclide(l) + + score = score + micro_xs(i_nuc) % nalpha * atom_density_ * flux + end do + end if + + + + case (SCORE_NGAMMA) + if (t % estimator == ESTIMATOR_ANALOG) then + if (survival_biasing) then + ! No absorption events actually occur if survival biasing is on -- + ! just use weight absorbed in survival biasing + score = p % absorb_wgt * flux + else + ! Skip any event where the particle wasn't absorbed + if (p % event == EVENT_SCATTER) cycle SCORE_LOOP + ! All fission and absorption events will contribute here, so we + ! can just use the particle's weight entering the collision + score = p % last_wgt * flux + end if + + else + if (i_nuclide > 0) then + score = micro_xs(i_nuclide) % ngamma * atom_density * flux + else + score = material_xs % ngamma * flux + end if + end if + + + case default if (t % estimator == ESTIMATOR_ANALOG) then ! Any other score is assumed to be a MT number. Thus, we just need @@ -2721,7 +2949,7 @@ contains type(Particle), intent(in) :: p real(8), intent(in) :: distance - + integer :: check_value integer :: i integer :: i_tally integer :: i_filt @@ -2802,6 +3030,7 @@ contains if (p % material /= MATERIAL_VOID) then ! Get pointer to current material mat => materials(p % material) + check_value= mat % mat_nuclide_list(i_nuclide) ! Determine index of nuclide in Material % atom_density array j = mat % mat_nuclide_index(i_nuclide) diff --git a/src/tallies/tally_header.F90 b/src/tallies/tally_header.F90 index 1778474e04..c283c3f16c 100644 --- a/src/tallies/tally_header.F90 +++ b/src/tallies/tally_header.F90 @@ -88,6 +88,7 @@ module tally_header ! reset property - allows a tally to be reset after every batch logical :: reset = .false. + logical :: has_threshold_rxn = .false. ! Number of realizations of tally random variables integer :: n_realizations = 0