diff --git a/docs/source/usersguide/input.rst b/docs/source/usersguide/input.rst index 87b613a45..7372f9e84 100644 --- a/docs/source/usersguide/input.rst +++ b/docs/source/usersguide/input.rst @@ -1233,8 +1233,8 @@ The ```` element accepts the following sub-elements: :type: The type of the filter. Accepted options are "cell", "cellborn", - "material", "universe", "energy", "energyout", "mesh", and - "distribcell". + "material", "universe", "energy", "energyout", "mesh", "distribcell", + and "delayedgroup". :bins: For each filter type, the corresponding ``bins`` entry is given as @@ -1352,6 +1352,15 @@ The ```` element accepts the following sub-elements: not accept more than one cell ID. It is not recommended to combine this filter with a cell or mesh filter. + :delayedgroup: + A list of delayed neutron precursor groups for which the tally should + be accumulated. For instance, to tally to all 6 delayed groups in the + ENDF/B-VII.1 library the filter is specified as: + + .. code-block:: xml + + + :nuclides: If specified, the scores listed will be for particular nuclides, not the summation of reactions from all nuclides. The format for nuclides should be @@ -1381,10 +1390,10 @@ The ```` element accepts the following sub-elements: :scores: A space-separated list of the desired responses to be accumulated. Accepted options are "flux", "total", "scatter", "absorption", "fission", - "nu-fission", "kappa-fission", "nu-scatter", "scatter-N", "scatter-PN", - "scatter-YN", "nu-scatter-N", "nu-scatter-PN", "nu-scatter-YN", "flux-YN", - "total-YN", "current", and "events". These corresponding to the following - physical quantities: + "nu-fission", "delayed-nu-fission", "kappa-fission", "nu-scatter", + "scatter-N", "scatter-PN", "scatter-YN", "nu-scatter-N", "nu-scatter-PN", + "nu-scatter-YN", "flux-YN", "total-YN", "current", and "events". These + correspond to the following physical quantities: :flux: Total flux in particle-cm per source particle. Note: The ``analog`` @@ -1409,6 +1418,10 @@ The ```` element accepts the following sub-elements: Total production of neutrons due to fission. Units are neutrons produced per source neutron. + :delayed-nu-fission: + Total production of delayed neutrons due to fission. Units are neutrons produced + per source neutron. + :kappa-fission: The recoverable energy production rate due to fission. The recoverable energy is defined as the fission product kinetic energy, prompt and diff --git a/openmc/filter.py b/openmc/filter.py index 8c1b6500b..51c60caf8 100644 --- a/openmc/filter.py +++ b/openmc/filter.py @@ -16,7 +16,7 @@ if sys.version_info[0] >= 3: _FILTER_TYPES = ['universe', 'material', 'cell', 'cellborn', 'surface', 'mesh', 'energy', 'energyout', 'mu', 'polar', 'azimuthal', - 'distribcell'] + 'distribcell', 'delayedgroup'] class Filter(object): """A filter used to constrain a tally to a specific criterion, e.g. only @@ -169,7 +169,7 @@ class Filter(object): bins = list(bins) if self.type in ['cell', 'cellborn', 'surface', 'material', - 'universe', 'distribcell']: + 'universe', 'distribcell', 'delayedgroup']: cv.check_iterable_type('filter bins', bins, Integral) for edge in bins: cv.check_greater_than('filter bin', edge, 0, equality=True) @@ -302,7 +302,11 @@ class Filter(object): merged_filter = copy.deepcopy(self) # Merge unique filter bins +<<<<<<< HEAD merged_bins = list(set(list(self.bins) + list(filter.bins))) +======= + merged_bins = list(set(np.concatenate((self.bins, filter.bins)))) +>>>>>>> upstream/develop merged_filter.bins = merged_bins merged_filter.num_bins = len(merged_bins) diff --git a/openmc/statepoint.py b/openmc/statepoint.py index b93624493..cc9da7562 100644 --- a/openmc/statepoint.py +++ b/openmc/statepoint.py @@ -1,5 +1,5 @@ import sys - +import re import numpy as np import openmc @@ -411,13 +411,10 @@ class StatePoint(object): # Add the scores to the Tally for j, score in enumerate(scores): score = score.decode() - # If this is a scattering moment, insert the scattering order - if '-n' in score: - score = score.replace('-n', '-' + moments[j].decode()) - elif '-pn' in score: - score = score.replace('-pn', '-' + moments[j].decode()) - elif '-yn' in score: - score = score.replace('-yn', '-' + moments[j].decode()) + + # If this is a moment, use generic moment order + pattern = r'-n$|-pn$|-yn$' + score = re.sub(pattern, '-' + moments[j].decode(), score) tally.add_score(score) diff --git a/openmc/tallies.py b/openmc/tallies.py index 6cc647911..b2ae4b000 100644 --- a/openmc/tallies.py +++ b/openmc/tallies.py @@ -586,6 +586,21 @@ class Tally(object): if len(self.filters) != len(tally.filters): return False + # Check if only one tally contains a delayed group filter + tally1_dg = False + for filter1 in self.filters: + if filter1.type == 'delayedgroup': + tally1_dg = True + + tally2_dg = False + for filter2 in tally.filters: + if filter2.type == 'delayedgroup': + tally2_dg = True + + # Return False if only one tally has a delayed group filter + if (tally1_dg or tally2_dg) and not (tally1_dg and tally2_dg): + return False + # Look to see if all filters are the same, or one or more can be merged for filter1 in self.filters: mergeable_filter = False diff --git a/src/ace.F90 b/src/ace.F90 index cccc57092..22fcc6b3d 100644 --- a/src/ace.F90 +++ b/src/ace.F90 @@ -215,6 +215,16 @@ contains end do MATERIAL_LOOP3 + ! Show which nuclide results in lowest energy for neutron transport + do i = 1, n_nuclides_total + if (nuclides(i)%energy(nuclides(i)%n_grid) == energy_max_neutron) then + call write_message("Maximum neutron transport energy: " // & + trim(to_str(energy_max_neutron)) // " MeV for " // & + trim(adjustl(nuclides(i)%name)), 6) + exit + end if + end do + end subroutine read_xs !=============================================================================== @@ -482,6 +492,10 @@ contains ! Continue reading elastic scattering and heating nuc % elastic = get_real(NE) + ! Determine if minimum/maximum energy for this nuclide is greater/less + ! than the previous + energy_min_neutron = max(energy_min_neutron, nuc%energy(1)) + energy_max_neutron = min(energy_max_neutron, nuc%energy(NE)) end if end subroutine read_esz @@ -635,6 +649,15 @@ contains ! Allocate space for secondary energy distribution NPCR = NXS(8) + + ! Check to make sure nuclide does not have more than the maximum number + ! of delayed groups + if (NPCR > MAX_DELAYED_GROUPS) then + call fatal_error("Encountered nuclide with " // trim(to_str(NPCR)) & + // " delayed groups while the maximum number of delayed groups & + &set in constants.F90 is " // trim(to_str(MAX_DELAYED_GROUPS))) + end if + nuc % n_precursor = NPCR allocate(nuc % nu_d_edist(NPCR)) @@ -672,6 +695,7 @@ contains else nuc % nu_d_type = NU_NONE + nuc % n_precursor = 0 end if end subroutine read_nu_data diff --git a/src/bank_header.F90 b/src/bank_header.F90 index 499358120..0cb49af35 100644 --- a/src/bank_header.F90 +++ b/src/bank_header.F90 @@ -11,10 +11,11 @@ module bank_header !=============================================================================== type, bind(C) :: Bank - real(C_DOUBLE) :: wgt ! weight of bank site - real(C_DOUBLE) :: xyz(3) ! location of bank particle - real(C_DOUBLE) :: uvw(3) ! diretional cosines - real(C_DOUBLE) :: E ! energy + real(C_DOUBLE) :: wgt ! weight of bank site + real(C_DOUBLE) :: xyz(3) ! location of bank particle + real(C_DOUBLE) :: uvw(3) ! diretional cosines + real(C_DOUBLE) :: E ! energy + integer(C_INT) :: delayed_group ! delayed group end type Bank end module bank_header diff --git a/src/constants.F90 b/src/constants.F90 index 9237dc2d0..9e973f24a 100644 --- a/src/constants.F90 +++ b/src/constants.F90 @@ -257,28 +257,29 @@ module constants EVENT_ABSORB = 2 ! Tally score type - integer, parameter :: N_SCORE_TYPES = 20 + integer, parameter :: N_SCORE_TYPES = 21 integer, parameter :: & - SCORE_FLUX = -1, & ! flux - SCORE_TOTAL = -2, & ! total reaction rate - SCORE_SCATTER = -3, & ! scattering rate - SCORE_NU_SCATTER = -4, & ! scattering production rate - SCORE_SCATTER_N = -5, & ! arbitrary scattering moment - SCORE_SCATTER_PN = -6, & ! system for scoring 0th through nth moment - SCORE_NU_SCATTER_N = -7, & ! arbitrary nu-scattering moment - SCORE_NU_SCATTER_PN = -8, & ! system for scoring 0th through nth nu-scatter moment - SCORE_TRANSPORT = -9, & ! transport reaction rate - SCORE_N_1N = -10, & ! (n,1n) rate - SCORE_ABSORPTION = -11, & ! absorption rate - SCORE_FISSION = -12, & ! fission rate - SCORE_NU_FISSION = -13, & ! neutron production rate - SCORE_KAPPA_FISSION = -14, & ! fission energy production rate - SCORE_CURRENT = -15, & ! partial current - SCORE_FLUX_YN = -16, & ! angular moment of flux - SCORE_TOTAL_YN = -17, & ! angular moment of total reaction rate - SCORE_SCATTER_YN = -18, & ! angular flux-weighted scattering moment (0:N) - SCORE_NU_SCATTER_YN = -19, & ! angular flux-weighted nu-scattering moment (0:N) - SCORE_EVENTS = -20 ! number of events + SCORE_FLUX = -1, & ! flux + SCORE_TOTAL = -2, & ! total reaction rate + SCORE_SCATTER = -3, & ! scattering rate + SCORE_NU_SCATTER = -4, & ! scattering production rate + SCORE_SCATTER_N = -5, & ! arbitrary scattering moment + SCORE_SCATTER_PN = -6, & ! system for scoring 0th through nth moment + SCORE_NU_SCATTER_N = -7, & ! arbitrary nu-scattering moment + SCORE_NU_SCATTER_PN = -8, & ! system for scoring 0th through nth nu-scatter moment + SCORE_TRANSPORT = -9, & ! transport reaction rate + SCORE_N_1N = -10, & ! (n,1n) rate + SCORE_ABSORPTION = -11, & ! absorption rate + SCORE_FISSION = -12, & ! fission rate + SCORE_NU_FISSION = -13, & ! neutron production rate + SCORE_KAPPA_FISSION = -14, & ! fission energy production rate + SCORE_CURRENT = -15, & ! partial current + SCORE_FLUX_YN = -16, & ! angular moment of flux + SCORE_TOTAL_YN = -17, & ! angular moment of total reaction rate + SCORE_SCATTER_YN = -18, & ! angular flux-weighted scattering moment (0:N) + SCORE_NU_SCATTER_YN = -19, & ! angular flux-weighted nu-scattering moment (0:N) + SCORE_EVENTS = -20, & ! number of events + SCORE_DELAYED_NU_FISSION = -21 ! delayed neutron production rate ! Maximum scattering order supported integer, parameter :: MAX_ANG_ORDER = 10 @@ -301,20 +302,21 @@ module constants integer, parameter :: NO_BIN_FOUND = -1 ! Tally filter and map types - integer, parameter :: N_FILTER_TYPES = 12 + integer, parameter :: N_FILTER_TYPES = 13 integer, parameter :: & - FILTER_UNIVERSE = 1, & - FILTER_MATERIAL = 2, & - FILTER_CELL = 3, & - FILTER_CELLBORN = 4, & - FILTER_SURFACE = 5, & - FILTER_MESH = 6, & - FILTER_ENERGYIN = 7, & - FILTER_ENERGYOUT = 8, & - FILTER_DISTRIBCELL = 9, & - FILTER_MU = 10, & - FILTER_POLAR = 11, & - FILTER_AZIMUTHAL = 12 + FILTER_UNIVERSE = 1, & + FILTER_MATERIAL = 2, & + FILTER_CELL = 3, & + FILTER_CELLBORN = 4, & + FILTER_SURFACE = 5, & + FILTER_MESH = 6, & + FILTER_ENERGYIN = 7, & + FILTER_ENERGYOUT = 8, & + FILTER_DISTRIBCELL = 9, & + FILTER_MU = 10, & + FILTER_POLAR = 11, & + FILTER_AZIMUTHAL = 12, & + FILTER_DELAYEDGROUP = 13 ! Mesh types integer, parameter :: & @@ -409,4 +411,14 @@ module constants ! constant for writing out no residual real(8), parameter :: CMFD_NORES = 99999.0_8 + !============================================================================= + ! DELAYED NEUTRON PRECURSOR CONSTANTS + + ! Since cross section libraries come with different numbers of delayed groups + ! (e.g. ENDF/B-VII.1 has 6 and JEFF 3.1.1 has 8 delayed groups) and we don't + ! yet know what cross section library is being used when the tallies.xml file + ! is read in, we want to have an upper bound on the size of the array we + ! use to store the bins for delayed group tallies. + integer, parameter :: MAX_DELAYED_GROUPS = 8 + end module constants diff --git a/src/cross_section.F90 b/src/cross_section.F90 index 4d8fb2f0f..1c56d961e 100644 --- a/src/cross_section.F90 +++ b/src/cross_section.F90 @@ -56,7 +56,7 @@ contains if (grid_method == GRID_MAT_UNION) then call find_energy_index(p % E, p % material) else if (grid_method == GRID_LOGARITHM) then - u = int(log(p % E/1.0e-11_8)/log_spacing) + u = int(log(p % E/energy_min_neutron)/log_spacing) end if ! Determine if this material has S(a,b) tables diff --git a/src/energy_grid.F90 b/src/energy_grid.F90 index 1473e5210..66419f83c 100644 --- a/src/energy_grid.F90 +++ b/src/energy_grid.F90 @@ -73,8 +73,8 @@ contains type(Nuclide), pointer :: nuc ! Set minimum/maximum energies - E_max = 20.0_8 - E_min = 1.0e-11_8 + E_max = energy_max_neutron + E_min = energy_min_neutron ! Determine equal-logarithmic energy spacing M = n_log_bins diff --git a/src/fission.F90 b/src/fission.F90 index 27143bf38..7b2911997 100644 --- a/src/fission.F90 +++ b/src/fission.F90 @@ -63,7 +63,7 @@ contains ! since no prompt or delayed data is present, this means all neutron ! emission is prompt -- WARNING: This currently returns zero. The calling ! routine needs to know this situation is occurring since we don't want - ! to call nu_total unnecessarily if it's already been called + ! to call nu_total unnecessarily if it has already been called. nu = ZERO elseif (nuc % nu_p_type == NU_POLYNOMIAL) then ! determine number of coefficients @@ -89,11 +89,15 @@ contains function nu_delayed(nuc, E) result(nu) - type(Nuclide), pointer :: nuc ! nuclide from which to find nu - real(8), intent(in) :: E ! energy of incoming neutron - real(8) :: nu ! number of delayed neutrons emitted per fission + type(Nuclide), intent(in) :: nuc ! nuclide from which to find nu + real(8), intent(in) :: E ! energy of incoming neutron + real(8) :: nu ! number of delayed neutrons emitted per fission if (nuc % nu_d_type == NU_NONE) then + ! since no prompt or delayed data is present, this means all neutron + ! emission is prompt -- WARNING: This currently returns zero. The calling + ! routine needs to know this situation is occurring since we don't want + ! to call nu_delayed unnecessarily if it has already been called. nu = ZERO elseif (nuc % nu_d_type == NU_TABULAR) then ! use ENDF interpolation laws to determine nu @@ -102,4 +106,60 @@ contains end function nu_delayed +!=============================================================================== +! YIELD_DELAYED calculates the fractional yield of delayed neutrons emitted for +! a given nuclide and incoming neutron energy in a given delayed group. +!=============================================================================== + + function yield_delayed(nuc, E, g) result(yield) + + type(Nuclide), intent(in) :: nuc ! nuclide from which to find nu + real(8), intent(in) :: E ! energy of incoming neutron + real(8) :: yield ! delayed neutron precursor yield + integer, intent(in) :: g ! the delayed neutron precursor group + integer :: d ! precursor group + integer :: lc ! index before start of energies/nu values + integer :: NR ! number of interpolation regions + integer :: NE ! number of energies tabulated + + yield = ZERO + + if (g > nuc % n_precursor .or. g < 1) then + ! if the precursor group is outside the range of precursor groups for + ! the input nuclide, return ZERO. + yield = ZERO + else if (nuc % nu_d_type == NU_NONE) then + ! since no prompt or delayed data is present, this means all neutron + ! emission is prompt -- WARNING: This currently returns zero. The calling + ! routine needs to know this situation is occurring since we don't want + ! to call yield_delayed unnecessarily if it has already been called. + yield = ZERO + else if (nuc % nu_d_type == NU_TABULAR) then + + lc = 1 + + ! loop over delayed groups and determine the yield for the desired group + do d = 1, nuc % n_precursor + + ! determine number of interpolation regions and energies + NR = int(nuc % nu_d_precursor_data(lc + 1)) + NE = int(nuc % nu_d_precursor_data(lc + 2 + 2*NR)) + + ! check if this is the desired group + if (d == g) then + + ! determine delayed neutron precursor yield for group g + yield = interpolate_tab1(nuc % nu_d_precursor_data( & + lc+1:lc+2+2*NR+2*NE), E) + + exit + end if + + ! advance pointer + lc = lc + 2 + 2*NR + 2*NE + 1 + end do + end if + + end function yield_delayed + end module fission diff --git a/src/global.F90 b/src/global.F90 index 87d727829..a4c80daa7 100644 --- a/src/global.F90 +++ b/src/global.F90 @@ -74,6 +74,10 @@ module global integer :: n_sab_tables ! Number of S(a,b) thermal scattering tables integer :: n_listings ! Number of listings in cross_sections.xml + ! Minimum/maximum energies + real(8) :: energy_min_neutron = ZERO + real(8) :: energy_max_neutron = INFINITY + ! Dictionaries to look up cross sections and listings type(DictCharInt) :: nuclide_dict type(DictCharInt) :: sab_dict diff --git a/src/initialize.F90 b/src/initialize.F90 index ca92595c2..9d63eb4e9 100644 --- a/src/initialize.F90 +++ b/src/initialize.F90 @@ -183,22 +183,22 @@ contains subroutine initialize_mpi() - integer :: bank_blocks(4) ! Count for each datatype + integer :: bank_blocks(5) ! Count for each datatype #ifdef MPIF08 - type(MPI_Datatype) :: bank_types(4) + type(MPI_Datatype) :: bank_types(5) type(MPI_Datatype) :: result_types(1) type(MPI_Datatype) :: temp_type #else - integer :: bank_types(4) ! Datatypes + integer :: bank_types(5) ! Datatypes integer :: result_types(1) ! Datatypes - integer :: temp_type ! temporary derived type + integer :: temp_type ! temporary derived type #endif - integer(MPI_ADDRESS_KIND) :: bank_disp(4) ! Displacements + integer(MPI_ADDRESS_KIND) :: bank_disp(5) ! Displacements integer :: result_blocks(1) ! Count for each datatype integer(MPI_ADDRESS_KIND) :: result_disp(1) ! Displacements integer(MPI_ADDRESS_KIND) :: result_base_disp ! Base displacement - integer(MPI_ADDRESS_KIND) :: lower_bound ! Lower bound for TallyResult - integer(MPI_ADDRESS_KIND) :: extent ! Extent for TallyResult + integer(MPI_ADDRESS_KIND) :: lower_bound ! Lower bound for TallyResult + integer(MPI_ADDRESS_KIND) :: extent ! Extent for TallyResult type(Bank) :: b type(TallyResult) :: tr @@ -223,18 +223,19 @@ contains ! CREATE MPI_BANK TYPE ! Determine displacements for MPI_BANK type - call MPI_GET_ADDRESS(b%wgt, bank_disp(1), mpi_err) - call MPI_GET_ADDRESS(b%xyz, bank_disp(2), mpi_err) - call MPI_GET_ADDRESS(b%uvw, bank_disp(3), mpi_err) - call MPI_GET_ADDRESS(b%E, bank_disp(4), mpi_err) + call MPI_GET_ADDRESS(b % wgt, bank_disp(1), mpi_err) + call MPI_GET_ADDRESS(b % xyz, bank_disp(2), mpi_err) + call MPI_GET_ADDRESS(b % uvw, bank_disp(3), mpi_err) + call MPI_GET_ADDRESS(b % E, bank_disp(4), mpi_err) + call MPI_GET_ADDRESS(b % delayed_group, bank_disp(5), mpi_err) ! Adjust displacements bank_disp = bank_disp - bank_disp(1) ! Define MPI_BANK for fission sites - bank_blocks = (/ 1, 3, 3, 1 /) - bank_types = (/ MPI_REAL8, MPI_REAL8, MPI_REAL8, MPI_REAL8 /) - call MPI_TYPE_CREATE_STRUCT(4, bank_blocks, bank_disp, & + bank_blocks = (/ 1, 3, 3, 1, 1 /) + bank_types = (/ MPI_REAL8, MPI_REAL8, MPI_REAL8, MPI_REAL8, MPI_INTEGER /) + call MPI_TYPE_CREATE_STRUCT(5, bank_blocks, bank_disp, & bank_types, MPI_BANK, mpi_err) call MPI_TYPE_COMMIT(MPI_BANK, mpi_err) @@ -306,6 +307,8 @@ contains c_loc(tmpb(1)%uvw)), coordinates_t, hdf5_err) call h5tinsert_f(hdf5_bank_t, "E", h5offsetof(c_loc(tmpb(1)), & c_loc(tmpb(1)%E)), H5T_NATIVE_DOUBLE, hdf5_err) + call h5tinsert_f(hdf5_bank_t, "delayed_group", h5offsetof(c_loc(tmpb(1)), & + c_loc(tmpb(1)%delayed_group)), H5T_NATIVE_INTEGER, hdf5_err) ! Determine type for integer(8) hdf5_integer8_t = h5kind_to_type(8, H5_INTEGER_KIND) diff --git a/src/input_xml.F90 b/src/input_xml.F90 index 8c11fbcf7..627889db5 100644 --- a/src/input_xml.F90 +++ b/src/input_xml.F90 @@ -2155,6 +2155,7 @@ contains subroutine read_tallies_xml() + integer :: d ! delayed group index integer :: i ! loop over user-specified tallies integer :: j ! loop over words integer :: k ! another loop index @@ -2586,6 +2587,28 @@ contains ! Set to analog estimator t % estimator = ESTIMATOR_ANALOG + case ('delayedgroup') + ! Set type of filter + t % filters(j) % type = FILTER_DELAYEDGROUP + + ! Set number of bins + t % filters(j) % n_bins = n_words + + ! Allocate and store bins + allocate(t % filters(j) % int_bins(n_words)) + call get_node_array(node_filt, "bins", t % filters(j) % int_bins) + + ! Check bins to make sure all are between 1 and MAX_DELAYED_GROUPS + do d = 1, n_words + if (t % filters(j) % int_bins(d) < 1 .or. & + t % filters(j) % int_bins(d) > MAX_DELAYED_GROUPS) then + call fatal_error("Encountered delayedgroup bin with index " & + // trim(to_str(t % filters(j) % int_bins(d))) // " that is& + & outside the range of 1 to MAX_DELAYED_GROUPS ( " & + // trim(to_str(MAX_DELAYED_GROUPS)) // ")") + end if + end do + case ('mu') ! Set type of filter t % filters(j) % type = FILTER_MU @@ -2741,6 +2764,21 @@ contains ! Check if total material was specified if (trim(sarray(j)) == 'total') then + + ! Check if a delayedgroup filter is present for this tally + do l = 1, t % n_filters + if (t % filters(l) % type == FILTER_DELAYEDGROUP) then + call warning("A delayedgroup filter was used on a total & + &nuclide tally. Cross section libraries are not & + &guaranteed to have the same delayed group structure & + &across all isotopes. In particular, ENDF/B-VII.1 does & + ¬ have a consistent delayed group structure across & + &all isotopes while the JEFF 3.1.1 library has the same & + &delayed group structure across all isotopes. Use with & + &caution!") + end if + end do + t % nuclide_bins(j) = -1 cycle end if @@ -2790,6 +2828,19 @@ contains allocate(t % nuclide_bins(1)) t % nuclide_bins(1) = -1 t % n_nuclide_bins = 1 + + ! Check if a delayedgroup filter is present for this tally + do l = 1, t % n_filters + if (t % filters(l) % type == FILTER_DELAYEDGROUP) then + call warning("A delayedgroup filter was used on a total nuclide & + &tally. Cross section libraries are not guaranteed to have the& + & same delayed group structure across all isotopes. In & + &particular, ENDF/B-VII.1 does not have a consistent delayed & + &group structure across all isotopes while the JEFF 3.1.1 & + &library has the same delayed group structure across all & + &isotopes. Use with caution!") + end if + end do end if ! ======================================================================= @@ -2902,6 +2953,14 @@ contains end do end if + ! Check if delayed group filter is used with any score besides + ! delayed-nu-fission + if (score_name /= 'delayed-nu-fission' .and. & + t % find_filter(FILTER_DELAYEDGROUP) > 0) then + call fatal_error("Cannot tally " // trim(score_name) // " with a & + &delayedgroup filter.") + end if + ! Check to see if the mu filter is applied and if that makes sense. if ((.not. starts_with(score_name,'scatter')) .and. & (.not. starts_with(score_name,'nu-scatter'))) then @@ -3061,6 +3120,12 @@ contains ! Set tally estimator to analog t % estimator = ESTIMATOR_ANALOG end if + case ('delayed-nu-fission') + t % score_bins(j) = SCORE_DELAYED_NU_FISSION + if (t % find_filter(FILTER_ENERGYOUT) > 0) then + ! Set tally estimator to analog + t % estimator = ESTIMATOR_ANALOG + end if case ('kappa-fission') t % score_bins(j) = SCORE_KAPPA_FISSION case ('current') diff --git a/src/output.F90 b/src/output.F90 index 84146f3d3..3463244e3 100644 --- a/src/output.F90 +++ b/src/output.F90 @@ -309,6 +309,7 @@ contains ! Display weight, energy, grid index, and interpolation factor write(ou,*) ' Weight = ' // to_str(p % wgt) write(ou,*) ' Energy = ' // to_str(p % E) + write(ou,*) ' Delayed Group = ' // to_str(p % delayed_group) write(ou,*) end subroutine print_particle @@ -949,39 +950,41 @@ contains if (n_tallies == 0) return ! Initialize names for tally filter types - filter_name(FILTER_UNIVERSE) = "Universe" - filter_name(FILTER_MATERIAL) = "Material" - filter_name(FILTER_DISTRIBCELL) = "Distributed Cell" - filter_name(FILTER_CELL) = "Cell" - filter_name(FILTER_CELLBORN) = "Birth Cell" - filter_name(FILTER_SURFACE) = "Surface" - filter_name(FILTER_MESH) = "Mesh" - filter_name(FILTER_ENERGYIN) = "Incoming Energy" - filter_name(FILTER_ENERGYOUT) = "Outgoing Energy" - filter_name(FILTER_MU) = "Change-in-Angle" - filter_name(FILTER_POLAR) = "Polar Angle" - filter_name(FILTER_AZIMUTHAL) = "Azimuthal Angle" + filter_name(FILTER_UNIVERSE) = "Universe" + filter_name(FILTER_MATERIAL) = "Material" + filter_name(FILTER_DISTRIBCELL) = "Distributed Cell" + filter_name(FILTER_CELL) = "Cell" + filter_name(FILTER_CELLBORN) = "Birth Cell" + filter_name(FILTER_SURFACE) = "Surface" + filter_name(FILTER_MESH) = "Mesh" + filter_name(FILTER_ENERGYIN) = "Incoming Energy" + filter_name(FILTER_ENERGYOUT) = "Outgoing Energy" + filter_name(FILTER_MU) = "Change-in-Angle" + filter_name(FILTER_POLAR) = "Polar Angle" + filter_name(FILTER_AZIMUTHAL) = "Azimuthal Angle" + filter_name(FILTER_DELAYEDGROUP) = "Delayed Group" ! Initialize names for scores - score_names(abs(SCORE_FLUX)) = "Flux" - score_names(abs(SCORE_TOTAL)) = "Total Reaction Rate" - score_names(abs(SCORE_SCATTER)) = "Scattering Rate" - score_names(abs(SCORE_NU_SCATTER)) = "Scattering Production Rate" - score_names(abs(SCORE_TRANSPORT)) = "Transport Rate" - score_names(abs(SCORE_N_1N)) = "(n,1n) Rate" - score_names(abs(SCORE_ABSORPTION)) = "Absorption Rate" - score_names(abs(SCORE_FISSION)) = "Fission Rate" - score_names(abs(SCORE_NU_FISSION)) = "Nu-Fission Rate" - score_names(abs(SCORE_KAPPA_FISSION)) = "Kappa-Fission Rate" - score_names(abs(SCORE_EVENTS)) = "Events" - score_names(abs(SCORE_FLUX_YN)) = "Flux Moment" - score_names(abs(SCORE_TOTAL_YN)) = "Total Reaction Rate Moment" - score_names(abs(SCORE_SCATTER_N)) = "Scattering Rate Moment" - score_names(abs(SCORE_SCATTER_PN)) = "Scattering Rate Moment" - score_names(abs(SCORE_SCATTER_YN)) = "Scattering Rate Moment" - score_names(abs(SCORE_NU_SCATTER_N)) = "Scattering Prod. Rate Moment" - score_names(abs(SCORE_NU_SCATTER_PN)) = "Scattering Prod. Rate Moment" - score_names(abs(SCORE_NU_SCATTER_YN)) = "Scattering Prod. Rate Moment" + score_names(abs(SCORE_FLUX)) = "Flux" + score_names(abs(SCORE_TOTAL)) = "Total Reaction Rate" + score_names(abs(SCORE_SCATTER)) = "Scattering Rate" + score_names(abs(SCORE_NU_SCATTER)) = "Scattering Production Rate" + score_names(abs(SCORE_TRANSPORT)) = "Transport Rate" + score_names(abs(SCORE_N_1N)) = "(n,1n) Rate" + score_names(abs(SCORE_ABSORPTION)) = "Absorption Rate" + score_names(abs(SCORE_FISSION)) = "Fission Rate" + score_names(abs(SCORE_NU_FISSION)) = "Nu-Fission Rate" + score_names(abs(SCORE_KAPPA_FISSION)) = "Kappa-Fission Rate" + score_names(abs(SCORE_EVENTS)) = "Events" + score_names(abs(SCORE_FLUX_YN)) = "Flux Moment" + score_names(abs(SCORE_TOTAL_YN)) = "Total Reaction Rate Moment" + score_names(abs(SCORE_SCATTER_N)) = "Scattering Rate Moment" + score_names(abs(SCORE_SCATTER_PN)) = "Scattering Rate Moment" + score_names(abs(SCORE_SCATTER_YN)) = "Scattering Rate Moment" + score_names(abs(SCORE_NU_SCATTER_N)) = "Scattering Prod. Rate Moment" + score_names(abs(SCORE_NU_SCATTER_PN)) = "Scattering Prod. Rate Moment" + score_names(abs(SCORE_NU_SCATTER_YN)) = "Scattering Prod. Rate Moment" + score_names(abs(SCORE_DELAYED_NU_FISSION)) = "Delayed-Nu-Fission Rate" ! Create filename for tally output filename = trim(path_output) // "tallies.out" @@ -1403,6 +1406,9 @@ contains E0 = t % filters(i_filter) % real_bins(bin) E1 = t % filters(i_filter) % real_bins(bin + 1) label = "[" // trim(to_str(E0)) // ", " // trim(to_str(E1)) // ")" + case (FILTER_DELAYEDGROUP) + i = t % filters(i_filter) % int_bins(bin) + label = to_str(i) end select end function get_label diff --git a/src/particle_header.F90 b/src/particle_header.F90 index 66df6b6e5..0b9b251ee 100644 --- a/src/particle_header.F90 +++ b/src/particle_header.F90 @@ -1,7 +1,8 @@ module particle_header use bank_header, only: Bank - use constants, only: NEUTRON, ONE, NONE, ZERO, MAX_SECONDARY + use constants, only: NEUTRON, ONE, NONE, ZERO, MAX_SECONDARY, & + MAX_DELAYED_GROUPS use error, only: fatal_error use geometry_header, only: BASE_UNIVERSE @@ -64,10 +65,13 @@ module particle_header integer :: event ! scatter, absorption integer :: event_nuclide ! index in nuclides array integer :: event_MT ! reaction MT + integer :: delayed_group ! delayed group ! Post-collision physical data integer :: n_bank ! number of fission sites banked real(8) :: wgt_bank ! weight of fission sites banked + integer :: n_delayed_bank(MAX_DELAYED_GROUPS) ! number of delayed fission + ! sites banked ! Indices for various arrays integer :: surface ! index for surface particle is on @@ -111,17 +115,19 @@ contains this % alive = .true. ! clear attributes - this % surface = NONE - this % cell_born = NONE - this % material = NONE - this % last_material = NONE - this % wgt = ONE - this % last_wgt = ONE - this % absorb_wgt = ZERO - this % n_bank = 0 - this % wgt_bank = ZERO - this % n_collision = 0 - this % fission = .false. + this % surface = NONE + this % cell_born = NONE + this % material = NONE + this % last_material = NONE + this % wgt = ONE + this % last_wgt = ONE + this % absorb_wgt = ZERO + this % n_bank = 0 + this % wgt_bank = ZERO + this % n_collision = 0 + this % fission = .false. + this % delayed_group = 0 + this % n_delayed_bank(:) = 0 ! Set up base level coordinates this % coord(1) % universe = BASE_UNIVERSE diff --git a/src/physics.F90 b/src/physics.F90 index 889126d7d..581cb04c9 100644 --- a/src/physics.F90 +++ b/src/physics.F90 @@ -1063,14 +1063,15 @@ contains integer, intent(in) :: i_nuclide integer, intent(in) :: i_reaction - integer :: i ! loop index - integer :: nu ! actual number of neutrons produced - integer :: ijk(3) ! indices in ufs mesh - real(8) :: nu_t ! total nu - real(8) :: mu ! fission neutron angular cosine - real(8) :: phi ! fission neutron azimuthal angle - real(8) :: weight ! weight adjustment for ufs method - logical :: in_mesh ! source site in ufs mesh? + integer :: nu_d(MAX_DELAYED_GROUPS) ! number of delayed neutrons born + integer :: i ! loop index + integer :: nu ! actual number of neutrons produced + integer :: ijk(3) ! indices in ufs mesh + real(8) :: nu_t ! total nu + real(8) :: mu ! fission neutron angular cosine + real(8) :: phi ! fission neutron azimuthal angle + real(8) :: weight ! weight adjustment for ufs method + logical :: in_mesh ! source site in ufs mesh? type(Nuclide), pointer :: nuc type(Reaction), pointer :: rxn @@ -1120,6 +1121,11 @@ contains ! Bank source neutrons if (nu == 0 .or. n_bank == size(fission_bank)) return + + ! Initialize counter of delayed neutrons encountered for each delayed group + ! to zero. + nu_d(:) = 0 + p % fission = .true. ! Fission neutrons will be banked do i = int(n_bank,4) + 1, int(min(n_bank + nu, int(size(fission_bank),8)),4) ! Bank source neutrons by copying particle data @@ -1142,15 +1148,24 @@ contains ! Sample secondary energy distribution for fission reaction and set energy ! in fission bank - fission_bank(i) % E = sample_fission_energy(nuc, rxn, p % E) + fission_bank(i) % E = sample_fission_energy(nuc, rxn, p) + + ! Set the delayed group of the neutron + fission_bank(i) % delayed_group = p % delayed_group + + ! Increment the number of neutrons born delayed + if (p % delayed_group > 0) then + nu_d(p % delayed_group) = nu_d(p % delayed_group) + 1 + end if end do ! increment number of bank sites n_bank = min(n_bank + nu, int(size(fission_bank),8)) - ! Store total weight banked for analog fission tallies + ! Store total and delayed weight banked for analog fission tallies p % n_bank = nu p % wgt_bank = nu/weight + p % n_delayed_bank(:) = nu_d(:) end subroutine create_fission_sites @@ -1158,12 +1173,12 @@ contains ! SAMPLE_FISSION_ENERGY !=============================================================================== - function sample_fission_energy(nuc, rxn, E) result(E_out) + function sample_fission_energy(nuc, rxn, p) result(E_out) - type(Nuclide), pointer :: nuc - type(Reaction), pointer :: rxn - real(8), intent(in) :: E ! incoming energy of neutron - real(8) :: E_out ! outgoing energy of fission neutron + type(Nuclide), pointer :: nuc + type(Reaction), pointer :: rxn + type(Particle), intent(inout) :: p ! Particle causing fission + real(8) :: E_out ! outgoing energy of fission neutron integer :: j ! index on nu energy grid / precursor group integer :: lc ! index before start of energies/nu values @@ -1181,10 +1196,10 @@ contains type(DistEnergy), pointer :: edist ! Determine total nu - nu_t = nu_total(nuc, E) + nu_t = nu_total(nuc, p % E) ! Determine delayed nu - nu_d = nu_delayed(nuc, E) + nu_d = nu_delayed(nuc, p % E) ! Determine delayed neutron fraction beta = nu_d / nu_t @@ -1204,7 +1219,7 @@ contains ! determine delayed neutron precursor yield for group j yield = interpolate_tab1(nuc % nu_d_precursor_data( & - lc+1:lc+2+2*NR+2*NE), E) + lc+1:lc+2+2*NR+2*NE), p % E) ! Check if this group is sampled prob = prob + yield @@ -1219,6 +1234,9 @@ contains ! n_precursor -- check for this condition j = min(j, nuc % n_precursor) + ! set the delayed group for the particle born from fission + p % delayed_group = j + ! select energy distribution for group j law = nuc % nu_d_edist(j) % law edist => nuc % nu_d_edist(j) @@ -1227,13 +1245,13 @@ contains n_sample = 0 do if (law == 44 .or. law == 61) then - call sample_energy(edist, E, E_out, mu) + call sample_energy(edist, p % E, E_out, mu) else - call sample_energy(edist, E, E_out) + call sample_energy(edist, p % E, E_out) end if - ! resample if energy is >= 20 MeV - if (E_out < 20) exit + ! resample if energy is greater than maximum neutron energy + if (E_out < energy_max_neutron) exit ! check for large number of resamples n_sample = n_sample + 1 @@ -1248,18 +1266,21 @@ contains ! ==================================================================== ! PROMPT NEUTRON SAMPLED + ! set the delayed group for the particle born from fission to 0 + p % delayed_group = 0 + ! sample from prompt neutron energy distribution law = rxn % edist % law n_sample = 0 do if (law == 44 .or. law == 61) then - call sample_energy(rxn%edist, E, E_out, prob) + call sample_energy(rxn%edist, p % E, E_out, prob) else - call sample_energy(rxn%edist, E, E_out) + call sample_energy(rxn%edist, p % E, E_out) end if - ! resample if energy is >= 20 MeV - if (E_out < 20) exit + ! resample if energy is greater than maximum neutron energy + if (E_out < energy_max_neutron) exit ! check for large number of resamples n_sample = n_sample + 1 diff --git a/src/relaxng/tallies.rnc b/src/relaxng/tallies.rnc index a19121845..75afb0f23 100644 --- a/src/relaxng/tallies.rnc +++ b/src/relaxng/tallies.rnc @@ -24,10 +24,10 @@ element tallies { element filter { (element type { ( "cell" | "cellborn" | "material" | "universe" | "surface" | "distribcell" | "mesh" | "energy" | "energyout" | "mu" | - "polar" | "azimuthal") } | + "polar" | "azimuthal" | "delayedgroup") } | attribute type { ( "cell" | "cellborn" | "material" | "universe" | "surface" | "distribcell" | "mesh" | "energy" | "energyout" | "mu" | - "polar" | "azimuthal") }) & + "polar" | "azimuthal" | "delayedgroup") }) & (element bins { list { xsd:double+ } } | attribute bins { list { xsd:double+ } }) }* & diff --git a/src/relaxng/tallies.rng b/src/relaxng/tallies.rng index cb952b930..d35c17004 100644 --- a/src/relaxng/tallies.rng +++ b/src/relaxng/tallies.rng @@ -145,6 +145,7 @@ mesh energy energyout + delayedgroup mu polar azimuthal @@ -161,6 +162,7 @@ mesh energy energyout + delayedgroup mu polar azimuthal diff --git a/src/source.F90 b/src/source.F90 index c46174947..6226517f3 100644 --- a/src/source.F90 +++ b/src/source.F90 @@ -211,8 +211,9 @@ contains case (SRC_ENERGY_MONO) ! Monoenergtic source site%E = external_source%params_energy(1) - if (site%E >= 20) then - call fatal_error("Source energies above 20 MeV not allowed.") + if (site%E >= energy_max_neutron) then + call fatal_error("Source energy above range of energies of at least & + &one cross section table") end if case (SRC_ENERGY_MAXWELL) @@ -221,8 +222,8 @@ contains ! Sample Maxwellian fission spectrum site%E = maxwell_spectrum(a) - ! resample if energy is >= 20 MeV - if (site%E < 20) exit + ! resample if energy is greater than maximum neutron energy + if (site%E < energy_max_neutron) exit end do case (SRC_ENERGY_WATT) @@ -232,8 +233,8 @@ contains ! Sample Watt fission spectrum site%E = watt_spectrum(a, b) - ! resample if energy is >= 20 MeV - if (site%E < 20) exit + ! resample if energy is greater than maximum neutron energy + if (site%E < energy_max_neutron) exit end do case default diff --git a/src/state_point.F90 b/src/state_point.F90 index c958cbc88..6d1c6f312 100644 --- a/src/state_point.F90 +++ b/src/state_point.F90 @@ -264,6 +264,8 @@ contains call write_dataset(filter_group, "type", "azimuthal") case(FILTER_DISTRIBCELL) call write_dataset(filter_group, "type", "distribcell") + case(FILTER_DELAYEDGROUP) + call write_dataset(filter_group, "type", "delayedgroup") end select call write_dataset(filter_group, "offset", tally%filters(j)%offset) @@ -335,6 +337,8 @@ contains str_array(j) = "fission" case (SCORE_NU_FISSION) str_array(j) = "nu-fission" + case (SCORE_DELAYED_NU_FISSION) + str_array(j) = "delayed-nu-fission" case (SCORE_KAPPA_FISSION) str_array(j) = "kappa-fission" case (SCORE_CURRENT) diff --git a/src/summary.F90 b/src/summary.F90 index d8ac1e3fa..bfa5108a0 100644 --- a/src/summary.F90 +++ b/src/summary.F90 @@ -580,6 +580,8 @@ contains call write_dataset(filter_group, "type", "polar") case(FILTER_AZIMUTHAL) call write_dataset(filter_group, "type", "azimuthal") + case(FILTER_DELAYEDGROUP) + call write_dataset(filter_group, "type", "delayedgroup") end select call close_group(filter_group) @@ -636,6 +638,8 @@ contains str_array(j) = "fission" case (SCORE_NU_FISSION) str_array(j) = "nu-fission" + case (SCORE_DELAYED_NU_FISSION) + str_array(j) = "delayed-nu-fission" case (SCORE_KAPPA_FISSION) str_array(j) = "kappa-fission" case (SCORE_CURRENT) diff --git a/src/tally.F90 b/src/tally.F90 index b1c9b7b34..4d4d6bf93 100644 --- a/src/tally.F90 +++ b/src/tally.F90 @@ -15,6 +15,8 @@ module tally use search, only: binary_search use string, only: to_str use tally_header, only: TallyResult, TallyMapItem, TallyMapElement + use fission, only: nu_total, nu_delayed, yield_delayed + use interpolation, only: interpolate_tab1 #ifdef MPI use message_passing @@ -53,6 +55,10 @@ contains integer :: i_energy ! index in nuclide energy grid integer :: score_bin ! scoring bin, e.g. SCORE_FLUX integer :: score_index ! scoring bin index + integer :: d ! delayed neutron index + integer :: d_bin ! delayed group bin index + integer :: dg_filter ! index of delayed group filter + real(8) :: yield ! delayed neutron yield real(8) :: atom_density_ ! atom/b-cm real(8) :: f ! interpolation factor real(8) :: score ! analog tally score @@ -61,6 +67,7 @@ contains real(8) :: uvw(3) ! particle direction type(Material), pointer :: mat type(Reaction), pointer :: rxn + type(Nuclide), pointer :: nuc i = 0 SCORE_LOOP: do q = 1, t % n_user_score_bins @@ -395,6 +402,186 @@ contains end if + case (SCORE_DELAYED_NU_FISSION) + + ! Set the delayedgroup filter index and the number of delayed group bins + dg_filter = t % find_filter(FILTER_DELAYEDGROUP) + + if (t % estimator == ESTIMATOR_ANALOG) then + if (survival_biasing .or. p % fission) then + if (t % find_filter(FILTER_ENERGYOUT) > 0) then + ! Normally, we only need to make contributions to one scoring + ! bin. However, in the case of fission, since multiple fission + ! neutrons were emitted with different energies, multiple + ! outgoing energy bins may have been scored to. The following + ! logic treats this special case and results to multiple bins + call score_fission_delayed_eout(p, t, score_index) + cycle SCORE_LOOP + end if + end if + if (survival_biasing) then + ! No fission events occur if survival biasing is on -- need to + ! calculate fraction of absorptions that would have resulted in + ! delayed-nu-fission + if (micro_xs(p % event_nuclide) % absorption > ZERO) then + + ! Get the event nuclide + nuc => nuclides(p % event_nuclide) + + ! Check if the delayed group filter is present + if (dg_filter > 0) then + + ! Loop over all delayed group bins and tally to them + ! individually + do d_bin = 1, t % filters(dg_filter) % n_bins + + ! Get the delayed group for this bin + d = t % filters(dg_filter) % int_bins(d_bin) + + ! Compute the yield for this delayed group + yield = yield_delayed(nuc, p % E, d) + + ! Compute the score and tally to bin + score = p % absorb_wgt * yield * micro_xs(p % event_nuclide) & + % fission * nu_delayed(nuc, p % E) / & + micro_xs(p % event_nuclide) % absorption + call score_fission_delayed_dg(t, d_bin, score, score_index) + end do + cycle SCORE_LOOP + else + ! If the delayed group filter is not present, compute the score + ! by multiplying the absorbed weight by the fraction of the + ! delayed-nu-fission xs to the absorption xs + score = p % absorb_wgt * micro_xs(p % event_nuclide) & + % fission * nu_delayed(nuc, p % E) / & + micro_xs(p % event_nuclide) % absorption + end if + end if + else + ! Skip any non-fission events + if (.not. p % fission) cycle SCORE_LOOP + ! If there is no outgoing energy filter, than we only need to + ! score to one bin. For the score to be 'analog', we need to + ! score the number of particles that were banked in the fission + ! bank. Since this was weighted by 1/keff, we multiply by keff + ! to get the proper score. Loop over the neutrons produced from + ! fission and check which ones are delayed. If a delayed neutron is + ! encountered, add its contribution to the fission bank to the + ! score. + + ! Check if the delayed group filter is present + if (dg_filter > 0) then + + ! Loop over all delayed group bins and tally to them individually + do d_bin = 1, t % filters(dg_filter) % n_bins + + ! Get the delayed group for this bin + d = t % filters(dg_filter) % int_bins(d_bin) + + ! Compute the score and tally to bin + score = keff * p % wgt_bank / p % n_bank * p % n_delayed_bank(d) + call score_fission_delayed_dg(t, d_bin, score, score_index) + end do + cycle SCORE_LOOP + else + + ! Add the contribution from all delayed groups + score = keff * p % wgt_bank / p % n_bank * sum(p % n_delayed_bank) + end if + end if + else + + ! Check if tally is on a single nuclide + if (i_nuclide > 0) then + + ! Get the nuclide of interest + nuc => nuclides(i_nuclide) + + ! Check if the delayed group filter is present + if (dg_filter > 0) then + + ! Loop over all delayed group bins and tally to them individually + do d_bin = 1, t % filters(dg_filter) % n_bins + + ! Get the delayed group for this bin + d = t % filters(dg_filter) % int_bins(d_bin) + + ! Compute the yield for this delayed group + yield = yield_delayed(nuc, p % E, d) + + ! Compute the score and tally to bin + score = micro_xs(i_nuclide) % fission * yield & + * nu_delayed(nuc, p % E) * atom_density * flux + call score_fission_delayed_dg(t, d_bin, score, score_index) + end do + cycle SCORE_LOOP + else + + ! If the delayed group filter is not present, compute the score + ! by multiplying the delayed-nu-fission macro xs by the flux + score = micro_xs(i_nuclide) % fission * nu_delayed(nuc, p % E)& + * atom_density * flux + end if + + ! Tally is on total nuclides + else + + ! Get pointer to current material + mat => materials(p % material) + + ! Check if the delayed group filter is present + if (dg_filter > 0) then + + ! Loop over all nuclides in the current material + do l = 1, mat % n_nuclides + + ! Get atom density + atom_density_ = mat % atom_density(l) + + ! Get index in nuclides array + i_nuc = mat % nuclide(l) + + ! Loop over all delayed group bins and tally to them individually + do d_bin = 1, t % filters(dg_filter) % n_bins + + ! Get the delayed group for this bin + d = t % filters(dg_filter) % int_bins(d_bin) + + ! Get the current nuclide + nuc => nuclides(i_nuc) + + ! Get the yield for the desired nuclide and delayed group + yield = yield_delayed(nuc, p % E, d) + + ! Compute the score and tally to bin + score = micro_xs(i_nuc) % fission * yield & + * nu_delayed(nuc, p % E) * atom_density_ * flux + call score_fission_delayed_dg(t, d_bin, score, score_index) + end do + end do + cycle SCORE_LOOP + else + + score = ZERO + + ! Loop over all nuclides in the current material + do l = 1, mat % n_nuclides + + ! Get atom density + atom_density_ = mat % atom_density(l) + + ! Get index in nuclides array + i_nuc = mat % nuclide(l) + + ! Accumulate the contribution from each nuclide + score = score + micro_xs(i_nuc) % fission & + * nu_delayed(nuclides(i_nuc), p % E) * atom_density_ * flux + end do + end if + end if + end if + + case (SCORE_KAPPA_FISSION) if (t % estimator == ESTIMATOR_ANALOG) then if (survival_biasing) then @@ -844,6 +1031,132 @@ contains end subroutine score_fission_eout +!=============================================================================== +! SCORE_FISSION_DELAYED_EOUT handles a special case where we need to store +! delayed neutron production rate with an outgoing energy filter (think of a +! fission matrix). In this case, we may need to score to multiple bins if there +! were multiple neutrons produced with different energies. +!=============================================================================== + + subroutine score_fission_delayed_eout(p, t, i_score) + + type(Particle), intent(in) :: p + type(TallyObject), intent(inout) :: t + integer, intent(in) :: i_score ! index for score + + integer :: i ! index of outgoing energy filter + integer :: j ! index of delayedgroup filter + integer :: d ! delayed group + integer :: g ! another delayed group + integer :: d_bin ! delayed group bin index + integer :: n ! number of energies on filter + integer :: k ! loop index for bank sites + integer :: bin_energyout ! original outgoing energy bin + integer :: i_filter ! index for matching filter bin combination + real(8) :: score ! actual score + real(8) :: E_out ! energy of fission bank site + + ! Save original outgoing energy bin + i = t % find_filter(FILTER_ENERGYOUT) + bin_energyout = matching_bins(i) + + ! Get the index of delayed group filter + j = t % find_filter(FILTER_DELAYEDGROUP) + + ! Get number of energies on filter + n = size(t % filters(i) % real_bins) + + ! Since the creation of fission sites is weighted such that it is + ! expected to create n_particles sites, we need to multiply the + ! score by keff to get the true delayed-nu-fission rate. + + ! loop over number of particles banked + do k = 1, p % n_bank + + ! get the delayed group + g = fission_bank(n_bank - p % n_bank + k) % delayed_group + + ! check if the particle was born delayed + if (g /= 0) then + + ! determine score based on bank site weight and keff + score = keff * fission_bank(n_bank - p % n_bank + k) % wgt + + ! determine outgoing energy from fission bank + E_out = fission_bank(n_bank - p % n_bank + k) % E + + ! check if outgoing energy is within specified range on filter + if (E_out < t % filters(i) % real_bins(1) .or. & + E_out > t % filters(i) % real_bins(n)) cycle + + ! change outgoing energy bin + matching_bins(i) = binary_search(t % filters(i) % real_bins, n, E_out) + + ! if the delayed group filter is present, tally to corresponding + ! delayed group bin if it exists + if (j > 0) then + + ! loop over delayed group bins until the corresponding bin is found + do d_bin = 1, t % filters(j) % n_bins + d = t % filters(j) % int_bins(d_bin) + + ! check whether the delayed group of the particle is equal to the + ! delayed group of this bin + if (d == g) then + call score_fission_delayed_dg(t, d_bin, score, i_score) + end if + end do + + ! if the delayed group filter is not present, add score to tally + else + + ! determine scoring index + i_filter = sum((matching_bins(1:t%n_filters) - 1) * t % stride) + 1 + + ! Add score to tally +!$omp atomic + t % results(i_score, i_filter) % value = & + t % results(i_score, i_filter) % value + score + end if + end if + end do + + ! reset outgoing energy bin + matching_bins(i) = bin_energyout + + end subroutine score_fission_delayed_eout + +!=============================================================================== +! SCORE_FISSION_DELAYED_DG helper function used to increment the tally when a +! delayed group filter is present. +!=============================================================================== + + subroutine score_fission_delayed_dg(t, d_bin, score, score_index) + + type(TallyObject), intent(inout) :: t + integer, intent(in) :: score_index ! index for score + integer, intent(in) :: d_bin ! delayed group bin index + + integer :: bin_original ! original bin index + integer :: filter_index ! index for matching filter bin combination + real(8) :: score ! actual score + + ! save original delayed group bin + bin_original = matching_bins(t % find_filter(FILTER_DELAYEDGROUP)) + matching_bins(t % find_filter(FILTER_DELAYEDGROUP)) = d_bin + + ! Compute the filter index based on the modified matching_bins + filter_index = sum((matching_bins(1:t%n_filters) - 1) * t % stride) + 1 + +!$omp atomic + t % results(score_index, filter_index) % value = & + t % results(score_index, filter_index) % value + score + + ! reset original delayed group bin + matching_bins(t % find_filter(FILTER_DELAYEDGROUP)) = bin_original + + end subroutine score_fission_delayed_dg + !=============================================================================== ! SCORE_TRACKLENGTH_TALLY calculates fluxes and reaction rates based on the ! track-length estimate of the flux. This is triggered at every event (surface @@ -1480,6 +1793,20 @@ contains n + 1, p % E) end if + case (FILTER_DELAYEDGROUP) + + if (survival_biasing .and. t % find_filter(FILTER_ENERGYOUT) <= 0) then + matching_bins(i) = 1 + elseif (active_tracklength_tallies % size() > 0) then + matching_bins(i) = 1 + else + if (p % delayed_group == 0) then + matching_bins = NO_BIN_FOUND + else + matching_bins(i) = p % delayed_group + end if + end if + case (FILTER_MU) ! determine mu bin n = t % filters(i) % n_bins diff --git a/src/tracking.F90 b/src/tracking.F90 index 81f01b411..2e4503e13 100644 --- a/src/tracking.F90 +++ b/src/tracking.F90 @@ -163,6 +163,7 @@ contains ! Reset banked weight during collision p % n_bank = 0 p % wgt_bank = ZERO + p % n_delayed_bank(:) = 0 ! Reset fission logical p % fission = .false. diff --git a/tests/test_filter_delayedgroup/inputs_true.dat b/tests/test_filter_delayedgroup/inputs_true.dat new file mode 100644 index 000000000..21bce00fd --- /dev/null +++ b/tests/test_filter_delayedgroup/inputs_true.dat @@ -0,0 +1 @@ +e771470681d3b4af57a70d148f5eb728df57f1fd7bdb5d6f76ac556b69f10bf0cdd5645fe488f1e17f07cbb72e9fb34af2fd7ede95c8e39c5ffa6ea9b6c5810e \ No newline at end of file diff --git a/tests/test_filter_delayedgroup/results_true.dat b/tests/test_filter_delayedgroup/results_true.dat new file mode 100644 index 000000000..9db6de256 --- /dev/null +++ b/tests/test_filter_delayedgroup/results_true.dat @@ -0,0 +1,15 @@ +k-combined: +9.903196E-01 4.279617E-02 +tally 1: +8.141852E-04 +1.337187E-07 +4.849156E-03 +4.744020E-06 +4.460252E-03 +4.015453E-06 +1.028479E-02 +2.136252E-05 +5.002274E-03 +5.056965E-06 +1.974747E-03 +7.882970E-07 diff --git a/tests/test_filter_delayedgroup/test_filter_delayedgroup.py b/tests/test_filter_delayedgroup/test_filter_delayedgroup.py new file mode 100644 index 000000000..bdad3cc95 --- /dev/null +++ b/tests/test_filter_delayedgroup/test_filter_delayedgroup.py @@ -0,0 +1,30 @@ +#!/usr/bin/env python + +import os +import sys +sys.path.insert(0, os.pardir) +from testing_harness import TestHarness, PyAPITestHarness +import openmc + + +class FilterDelayedgroupTestHarness(PyAPITestHarness): + def _build_inputs(self): + filt = openmc.Filter(type='delayedgroup', + bins=(1, 2, 3, 4, 5, 6)) + tally = openmc.Tally(tally_id=1) + tally.add_filter(filt) + tally.add_score('delayed-nu-fission') + self._input_set.tallies = openmc.TalliesFile() + self._input_set.tallies.add_tally(tally) + + super(FilterDelayedgroupTestHarness, self)._build_inputs() + + def _cleanup(self): + super(FilterDelayedgroupTestHarness, self)._cleanup() + f = os.path.join(os.getcwd(), 'tallies.xml') + if os.path.exists(f): os.remove(f) + + +if __name__ == '__main__': + harness = FilterDelayedgroupTestHarness('statepoint.10.*', True) + harness.main() diff --git a/tests/test_many_scores/results_true.dat b/tests/test_many_scores/results_true.dat index 6c1f62732..9309d0964 100644 --- a/tests/test_many_scores/results_true.dat +++ b/tests/test_many_scores/results_true.dat @@ -109,3 +109,5 @@ tally 1: 7.673412E-04 1.014000E+01 3.427342E+01 +7.652723E-03 +3.578992E-05 diff --git a/tests/test_many_scores/tallies.xml b/tests/test_many_scores/tallies.xml index 1832acf67..2df5597d0 100644 --- a/tests/test_many_scores/tallies.xml +++ b/tests/test_many_scores/tallies.xml @@ -6,7 +6,7 @@ flux total scatter nu-scatter scatter-2 scatter-p2 nu-scatter-2 nu-scatter-p2 transport n1n absorption nu-fission kappa-fission - flux-y2 total-y2 scatter-y2 nu-scatter-y2 events + flux-y2 total-y2 scatter-y2 nu-scatter-y2 events delayed-nu-fission diff --git a/tests/test_score_delayed_nufission/inputs_true.dat b/tests/test_score_delayed_nufission/inputs_true.dat new file mode 100644 index 000000000..5a0ec8a21 --- /dev/null +++ b/tests/test_score_delayed_nufission/inputs_true.dat @@ -0,0 +1 @@ +f3c246a1c83b1283163b22069f78231e3f6623fa3c2eb664ce400d0fff07105b6106d0fa8c6dfd49512a7eb676c80e4ee602e19063a9fc453394fe07a94a1bdd \ No newline at end of file diff --git a/tests/test_score_delayed_nufission/results_true.dat b/tests/test_score_delayed_nufission/results_true.dat new file mode 100644 index 000000000..d797baf18 --- /dev/null +++ b/tests/test_score_delayed_nufission/results_true.dat @@ -0,0 +1,29 @@ +k-combined: +9.903196E-01 4.279617E-02 +tally 1: +1.711611E-02 +5.967549E-05 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +1.026930E-02 +2.198717E-05 +tally 2: +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +1.976462E-02 +1.953328E-04 +tally 3: +1.687894E-02 +5.776176E-05 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +1.061803E-02 +2.308636E-05 diff --git a/tests/test_score_delayed_nufission/test_score_delayed_nufission.py b/tests/test_score_delayed_nufission/test_score_delayed_nufission.py new file mode 100644 index 000000000..0738520e2 --- /dev/null +++ b/tests/test_score_delayed_nufission/test_score_delayed_nufission.py @@ -0,0 +1,32 @@ +#!/usr/bin/env python + +import os +import sys +sys.path.insert(0, os.pardir) +from testing_harness import TestHarness, PyAPITestHarness +import openmc + + +class ScoreDelayedNuFissionTestHarness(PyAPITestHarness): + def _build_inputs(self): + filt = openmc.Filter(type='cell', bins=(21, 22, 23, 27)) + tallies = [openmc.Tally(tally_id=i) for i in range(1, 4)] + [t.add_filter(filt) for t in tallies] + [t.add_score('delayed-nu-fission') for t in tallies] + tallies[0].estimator = 'tracklength' + tallies[1].estimator = 'analog' + tallies[2].estimator = 'collision' + self._input_set.tallies = openmc.TalliesFile() + [self._input_set.tallies.add_tally(t) for t in tallies] + + super(ScoreDelayedNuFissionTestHarness, self)._build_inputs() + + def _cleanup(self): + super(ScoreDelayedNuFissionTestHarness, self)._cleanup() + f = os.path.join(os.getcwd(), 'tallies.xml') + if os.path.exists(f): os.remove(f) + + +if __name__ == '__main__': + harness = ScoreDelayedNuFissionTestHarness('statepoint.10.*', True) + harness.main()