diff --git a/openmc/settings.py b/openmc/settings.py index 3efdf1ab7a..658dd8642c 100644 --- a/openmc/settings.py +++ b/openmc/settings.py @@ -577,16 +577,15 @@ class Settings(object): raise ValueError(msg) for key in cutoff: if key == 'weight': - cv.check_type('weight cutoff', cutoff['weight'], Real) - cv.check_greater_than('weight cutoff', cutoff['weight'], 0.0) + cv.check_type('weight cutoff', cutoff[key], Real) + cv.check_greater_than('weight cutoff', cutoff[key], 0.0) elif key == 'weight_avg': - cv.check_type('average survival weight', cutoff['weight_avg'], - Real) + cv.check_type('average survival weight', cutoff[key], Real) cv.check_greater_than('average survival weight', - cutoff['weight_avg'], 0.0) - elif key == 'energy': - cv.check_type('energy cutoff', cutoff['energy'], Real) - cv.check_greater_than('energy cutoff', cutoff['energy'], 0.0) + cutoff[key], 0.0) + elif key in ['energy', 'energy_photon']: + cv.check_type('energy cutoff', cutoff[key], Real) + cv.check_greater_than('energy cutoff', cutoff[key], 0.0) else: msg = 'Unable to set cutoff to "{0}" which is unsupported by '\ 'OpenMC'.format(key) @@ -946,17 +945,9 @@ class Settings(object): def _create_cutoff_subelement(self, root): if self._cutoff is not None: element = ET.SubElement(root, "cutoff") - if 'weight' in self._cutoff: - subelement = ET.SubElement(element, "weight") - subelement.text = str(self._cutoff['weight']) - - if 'weight_avg' in self._cutoff: - subelement = ET.SubElement(element, "weight_avg") - subelement.text = str(self._cutoff['weight_avg']) - - if 'energy' in self._cutoff: - subelement = ET.SubElement(element, "energy") - subelement.text = str(self._cutoff['energy']) + for key, value in self.items(): + subelement = ET.SubElement(element, key) + subelement.text = str(value) def _create_entropy_subelement(self, root): if self._entropy_mesh is not None: diff --git a/src/constants.F90 b/src/constants.F90 index 88bf70b591..53253e23db 100644 --- a/src/constants.F90 +++ b/src/constants.F90 @@ -226,7 +226,9 @@ module constants N_3HEA = 193, N_4N2P = 194, N_4N2A = 195, N_4NPA = 196, N_3P = 197, & N_N3P = 198, N_3N2PA = 199, N_5N2P = 200, N_P0 = 600, N_PC = 649, & N_D0 = 650, N_DC = 699, N_T0 = 700, N_TC = 749, N_3HE0 = 750, & - N_3HEC = 799, N_A0 = 800, N_AC = 849, N_2N0 = 875, N_2NC = 891 + N_3HEC = 799, N_A0 = 800, N_AC = 849, N_2N0 = 875, N_2NC = 891, & + COHERENT = 502, INCOHERENT = 504, PHOTOELECTRIC = 522, & + PAIR_PROD_ELEC = 515, PAIR_PROD = 516, PAIR_PROD_NUC = 517 ! ACE table types integer, parameter :: & diff --git a/src/endf.F90 b/src/endf.F90 index 0ba2db388b..f3c6b70899 100644 --- a/src/endf.F90 +++ b/src/endf.F90 @@ -176,6 +176,20 @@ contains string = '(n,Xa)' case (444) string = '(damage)' + case (COHERENT) + string = 'coherent scatter' + case (INCOHERENT) + string = 'incoherent scatter' + case (PAIR_PROD_ELEC) + string = 'pair production, electron' + case (PAIR_PROD) + string = 'pair production' + case (PAIR_PROD_NUC) + string = 'pair production, nuclear' + case (PHOTOELECTRIC) + string = 'photoelectric' + case (534 : 572) + string = 'photoelectric, ' // trim(SUBSHELLS(MT - 533)) // ' subshell' case (600 : 648) string = '(n,p' // trim(to_str(MT-600)) // ')' case (649) diff --git a/src/global.F90 b/src/global.F90 index 6a37f2d748..d9da251bad 100644 --- a/src/global.F90 +++ b/src/global.F90 @@ -313,7 +313,7 @@ module global logical :: survival_biasing = .false. real(8) :: weight_cutoff = 0.25_8 - real(8) :: energy_cutoff = ZERO + real(8) :: energy_cutoff(3) = [ZERO, 1000.0_8, ZERO] real(8) :: weight_survive = ONE ! ============================================================================ diff --git a/src/input_xml.F90 b/src/input_xml.F90 index c0a91aecbd..0dc1a90a83 100644 --- a/src/input_xml.F90 +++ b/src/input_xml.F90 @@ -599,7 +599,10 @@ contains call get_node_value(node_cutoff, "weight_avg", weight_survive) end if if (check_for_node(node_cutoff, "energy")) then - call get_node_value(node_cutoff, "energy", energy_cutoff) + call get_node_value(node_cutoff, "energy", energy_cutoff(1)) + end if + if (check_for_node(node_cutoff, "energy_photon")) then + call get_node_value(node_cutoff, "energy_photon", energy_cutoff(2)) end if end if diff --git a/src/photon_header.F90 b/src/photon_header.F90 index dfedd5b0b9..a5cdd9791b 100644 --- a/src/photon_header.F90 +++ b/src/photon_header.F90 @@ -1,6 +1,6 @@ module photon_header - use hdf5, only: HID_T, HSIZE_T + use hdf5, only: HID_T, HSIZE_T, SIZE_T use constants, only: ZERO, HALF, SUBSHELLS use dict_header, only: DictIntInt @@ -10,7 +10,7 @@ module photon_header real(8), allocatable :: compton_profile_pz(:) type ElectronSubshell - character(3) :: label + integer :: index_subshell ! index in SUBSHELLS integer :: threshold real(8) :: n_electrons real(8) :: binding_energy @@ -24,7 +24,8 @@ module photon_header end type ElectronSubshell type PhotonInteraction - integer :: Z ! atomic number + character(3) :: name ! atomic symbol, e.g. 'Zr' + integer :: Z ! atomic number ! Microscopic cross sections real(8), allocatable :: energy(:) @@ -82,15 +83,22 @@ contains integer(HID_T) :: rgroup, tgroup integer(HID_T) :: dset_id integer(HSIZE_T) :: dims(1), dims2(2) + integer(SIZE_T) :: name_len integer :: n_energy integer :: n_shell integer :: n_profile integer :: n_transition - ! integer, allocatable :: designators(:) character(3), allocatable :: designators(:) real(8) :: c real(8), allocatable :: matrix(:,:) + ! Get name of nuclide from group + name_len = len(this % name) + this % name = get_name(group_id, name_len) + + ! Get rid of leading '/' + this % name = trim(this % name(2:)) + ! Get atomic number call read_attribute(this % Z, group_id, 'Z') @@ -160,6 +168,7 @@ contains do j = 1, size(SUBSHELLS) if (designators(i) == SUBSHELLS(j)) then call this % shell_dict % add_key(j, i) + this % shells(i) % index_subshell = j exit end if end do diff --git a/src/physics.F90 b/src/physics.F90 index d6a4c0ac52..7582190834 100644 --- a/src/physics.F90 +++ b/src/physics.F90 @@ -45,42 +45,42 @@ contains ! Sample reaction for the material the particle is in if (p % type == NEUTRON) then - call sample_reaction(p) - elseif (p % type == PHOTON) then + call sample_neutron_reaction(p) + else call sample_photon_reaction(p) end if + ! Kill particle if energy falls below cutoff + if (p % E < energy_cutoff(p % type)) then + p % alive = .false. + p % wgt = ZERO + p % last_wgt = ZERO + end if + ! Display information about collision if (verbosity >= 10 .or. trace) then - call write_message(" " // trim(reaction_name(p % event_MT)) & - &// " with " // trim(adjustl(nuclides(p % event_nuclide) % name)) & - &// ". Energy = " // trim(to_str(p % E)) // " eV.") + if (p % type == NEUTRON) then + call write_message(" " // trim(reaction_name(p % event_MT)) & + &// " with " // trim(adjustl(nuclides(p % event_nuclide) % name)) & + &// ". Energy = " // trim(to_str(p % E)) // " eV.") + else + call write_message(" " // trim(reaction_name(p % event_MT)) & + &// " with " // trim(adjustl(elements(p % event_nuclide) % name)) & + &// ". Energy = " // trim(to_str(p % E)) // " eV.") + end if end if - ! check for very low energy - if (p % E < 1.0e-100_8) then - p % alive = .false. - if (master) call warning("Killing neutron with extremely low energy") - end if - - ! Advance URR seed stream 'N' times after energy changes - if (p % E /= p % last_E) then - call prn_set_stream(STREAM_URR_PTABLE) - call advance_prn_seed(size(nuclides, kind=8)) - call prn_set_stream(STREAM_TRACKING) - endif - end subroutine collision !=============================================================================== -! SAMPLE_REACTION samples a nuclide based on the macroscopic cross sections for -! each nuclide within a material and then samples a reaction for that nuclide -! and calls the appropriate routine to process the physics. Note that there is -! special logic when suvival biasing is turned on since fission and -! disappearance are treated implicitly. +! SAMPLE_NEUTRON_REACTION samples a nuclide based on the macroscopic cross +! sections for each nuclide within a material and then samples a reaction for +! that nuclide and calls the appropriate routine to process the physics. Note +! that there is special logic when suvival biasing is turned on since fission +! and disappearance are treated implicitly. !=============================================================================== - subroutine sample_reaction(p) + subroutine sample_neutron_reaction(p) type(Particle), intent(inout) :: p @@ -128,20 +128,19 @@ contains call scatter(p, i_nuclide, i_nuc_mat) ! Play russian roulette if survival biasing is turned on - if (survival_biasing) then call russian_roulette(p) if (.not. p % alive) return end if - ! Kill neutron under certain energy - if (p % E < energy_cutoff) then - p % alive = .false. - p % wgt = ZERO - p % last_wgt = ZERO + ! Advance URR seed stream 'N' times after energy changes + if (p % E /= p % last_E) then + call prn_set_stream(STREAM_URR_PTABLE) + call advance_prn_seed(size(nuclides, kind=8)) + call prn_set_stream(STREAM_TRACKING) end if - end subroutine sample_reaction + end subroutine sample_neutron_reaction !=============================================================================== ! SAMPLE_PHOTON_REACTION samples an element based on the macroscopic cross @@ -184,6 +183,7 @@ contains if (prob > cutoff) then call rayleigh_scatter(elm, alpha, mu) p % coord(1) % uvw = rotate_angle(p % coord(1) % uvw, mu) + p % event_MT = COHERENT return end if @@ -193,6 +193,7 @@ contains call compton_scatter(elm, alpha, alpha_out, mu, .true.) p % E = alpha_out*MASS_ELECTRON p % coord(1) % uvw = rotate_angle(p % coord(1) % uvw, mu) + p % event_MT = INCOHERENT return end if @@ -219,6 +220,7 @@ contains ! E_electron = p % E - elm % shells(i_shell) % binding_energy call atomic_relaxation(p, elm, i_shell) + p % event_MT = 533 + elm % shells(i_shell) % index_subshell p % alive = .false. return end if @@ -239,6 +241,7 @@ contains ! Set energy p % E = MASS_ELECTRON + p % event_MT = PAIR_PROD ! Create photon in opposite direction call p % create_secondary(-p % coord(1) % uvw, MASS_ELECTRON, &