From ec788085b3baf9c514c8fe73354550dcfe4badcb Mon Sep 17 00:00:00 2001 From: Adam Nelson Date: Tue, 25 Aug 2015 19:42:19 -0400 Subject: [PATCH] Added support for change-in-angle (mu) filters for scattering scores --- openmc/constants.py | 3 ++- openmc/statepoint.py | 3 +-- openmc/tallies.py | 20 +++++++++++++-- src/constants.F90 | 5 ++-- src/hdf5_summary.F90 | 7 +++++- src/input_xml.F90 | 60 ++++++++++++++++++++++++++++++++++++++++++-- src/output.F90 | 14 ++++++++++- src/state_point.F90 | 6 +++-- src/tally.F90 | 8 +++++- 9 files changed, 112 insertions(+), 14 deletions(-) diff --git a/openmc/constants.py b/openmc/constants.py index a6b535e6d1..1423a50029 100644 --- a/openmc/constants.py +++ b/openmc/constants.py @@ -35,7 +35,8 @@ FILTER_TYPES = {1: 'universe', 6: 'mesh', 7: 'energy', 8: 'energyout', - 9: 'distribcell'} + 9: 'distribcell', + 10:'mu'} SCORE_TYPES = {-1: 'flux', -2: 'total', diff --git a/openmc/statepoint.py b/openmc/statepoint.py index 6a4713e91d..8df9369bf7 100644 --- a/openmc/statepoint.py +++ b/openmc/statepoint.py @@ -377,7 +377,7 @@ class StatePoint(object): raise ValueError(msg) # Read the bin values - if FILTER_TYPES[filter_type] in ['energy', 'energyout']: + if FILTER_TYPES[filter_type] in ['energy', 'energyout', 'mu']: bins = self._get_double( n_bins+1, path='{0}{1}/bins'.format(subbase, j)) @@ -417,7 +417,6 @@ class StatePoint(object): path='{0}{1}/n_score_bins'.format(base, tally_key))[0] tally.num_score_bins = n_score_bins - scores = [SCORE_TYPES[j] for j in self._get_int( n_score_bins, path='{0}{1}/score_bins'.format(base, tally_key))] n_user_scores = self._get_int( diff --git a/openmc/tallies.py b/openmc/tallies.py index 6e418a2be4..5bf064f23d 100644 --- a/openmc/tallies.py +++ b/openmc/tallies.py @@ -373,7 +373,7 @@ class Tally(object): if score in self.scores: return else: - self._scores.append(score) + self._scores.append(score.strip()) @num_score_bins.setter def num_score_bins(self, num_score_bins): @@ -896,7 +896,7 @@ class Tally(object): bins = list(itertools.product(*xyz)) # Create list of 2-tuples for energy boundary bins - elif filter.type in ['energy', 'energyout']: + elif filter.type in ['energy', 'energyout', 'mu']: bins = [] for k in range(filter.num_bins): bins.append((filter.bins[k], filter.bins[k+1])) @@ -1251,6 +1251,22 @@ class Tally(object): filter_bins = np.tile(filter_bins, tile_factor) df[filter.type + ' [MeV]'] = filter_bins + elif 'mu' is filter.type: + bins = filter.bins + num_bins = filter.num_bins + + # Create strings for + template = '{0:.2f} - {1:.2f}' + filter_bins = [] + for i in range(num_bins): + filter_bins.append(template.format(bins[i], bins[i+1])) + + # Tile the mu bins into a DataFrame column + filter_bins = np.repeat(filter_bins, filter.stride) + tile_factor = data_size / len(filter_bins) + filter_bins = np.tile(filter_bins, tile_factor) + df[filter.type] = filter_bins + # universe, material, surface, cell, and cellborn filters else: filter_bins = np.repeat(filter.bins, filter.stride) diff --git a/src/constants.F90 b/src/constants.F90 index 3aaf08f0eb..e42848acfb 100644 --- a/src/constants.F90 +++ b/src/constants.F90 @@ -300,7 +300,7 @@ module constants integer, parameter :: NO_BIN_FOUND = -1 ! Tally filter and map types - integer, parameter :: N_FILTER_TYPES = 9 + integer, parameter :: N_FILTER_TYPES = 10 integer, parameter :: & FILTER_UNIVERSE = 1, & FILTER_MATERIAL = 2, & @@ -310,7 +310,8 @@ module constants FILTER_MESH = 6, & FILTER_ENERGYIN = 7, & FILTER_ENERGYOUT = 8, & - FILTER_DISTRIBCELL = 9 + FILTER_DISTRIBCELL = 9, & + FILTER_MU = 10 ! Tally surface current directions integer, parameter :: & diff --git a/src/hdf5_summary.F90 b/src/hdf5_summary.F90 index 375ab74601..4033019b1c 100644 --- a/src/hdf5_summary.F90 +++ b/src/hdf5_summary.F90 @@ -607,7 +607,8 @@ contains ! Write filter bins if (t % filters(j) % type == FILTER_ENERGYIN .or. & - t % filters(j) % type == FILTER_ENERGYOUT) then + t % filters(j) % type == FILTER_ENERGYOUT .or. & + t % filters(j) % type == FILTER_MU) then call su % write_data(t % filters(j) % real_bins, "bins", & length=size(t % filters(j) % real_bins), & group="tallies/tally " // trim(to_str(t % id)) & @@ -653,6 +654,10 @@ contains call su % write_data("energyout", "type_name", & group="tallies/tally " // trim(to_str(t % id)) & // "/filter " // trim(to_str(j))) + case(FILTER_MU) + call su % write_data("mu", "type_name", & + group="tallies/tally " // trim(to_str(t % id)) & + // "/filter " // trim(to_str(j))) end select end do FILTER_LOOP diff --git a/src/input_xml.F90 b/src/input_xml.F90 index 5c939a394f..22f991ab9e 100644 --- a/src/input_xml.F90 +++ b/src/input_xml.F90 @@ -2085,6 +2085,9 @@ contains integer :: imomstr ! Index of MOMENT_STRS & MOMENT_N_STRS logical :: file_exists ! does tallies.xml file exist? real(8) :: rarray3(3) ! temporary double prec. array + integer :: Nmu ! Number of angular bins + real(8) :: dmu ! Mu spacing if using automatic allocation + integer :: imu ! Loop counter for building mu filter bins character(MAX_LINE_LEN) :: filename character(MAX_WORD_LEN) :: word character(MAX_WORD_LEN) :: score_name @@ -2350,8 +2353,9 @@ contains ! Determine number of bins if (check_for_node(node_filt, "bins")) then - if (trim(temp_str) == 'energy' .or. & - trim(temp_str) == 'energyout') then + if ((trim(temp_str) == 'energy' .or. & + trim(temp_str) == 'energyout') .or. & + (trim(temp_str) == 'mu')) then n_words = get_arraysize_double(node_filt, "bins") else n_words = get_arraysize_integer(node_filt, "bins") @@ -2487,6 +2491,40 @@ contains ! Set to analog estimator t % estimator = ESTIMATOR_ANALOG + case ('mu') + ! Set type of filter + t % filters(j) % type = FILTER_MU + + ! Set number of bins + t % filters(j) % n_bins = n_words - 1 + + ! Allocate and store bins + allocate(t % filters(j) % real_bins(n_words)) + call get_node_array(node_filt, "bins", t % filters(j) % real_bins) + + ! Easter egg! Allow a user to input a negative number, if it is only one + ! and that will mean you subivide [-1,1] evenly with the input + ! being the number of bins + if (n_words == 1) then + Nmu = abs(int(t % filters(j) % real_bins(1))) + if (Nmu > 1) then + t % filters(j) % n_bins = Nmu - 1 + dmu = TWO / (real(Nmu,8) - ONE) + deallocate(t % filters(j) % real_bins) + allocate(t % filters(j) % real_bins(Nmu)) + do imu = 1, Nmu + t % filters(j) % real_bins(imu) = -ONE + (imu - 1) * dmu + end do + else + call fatal_error("Must have more than one bin for mu filter & + & on tally " // trim(to_str(t % id)) // ".") + end if + + end if + + ! Set to analog estimator + t % estimator = ESTIMATOR_ANALOG + case default ! Specified tally filter is invalid, raise error call fatal_error("Unknown filter type '" & @@ -2656,6 +2694,7 @@ contains ! scores then strip off the n and store it as an integer to be used ! later. Then perform the select case on this modified (number ! removed) string + n_order = -1 score_name = sarray(l) do imomstr = 1, size(MOMENT_STRS) if (starts_with(score_name,trim(MOMENT_STRS(imomstr)))) then @@ -2701,6 +2740,23 @@ contains end do 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 + if (t % find_filter(FILTER_MU) > 0) then + call fatal_error("Cannot tally " // trim(score_name) //" with a & + &change of angle (mu) filter.") + end if + ! Also check to see if this is a legendre expansion or not. + ! If so, we can accept this score and filter combo for p0, but not + ! elsewhere. + else if (n_order > 0) then + if (t % find_filter(FILTER_MU) > 0) then + call fatal_error("Cannot tally " // trim(score_name) //" with a & + &change of angle (mu) filter unless order is 0.") + end if + end if + select case (trim(score_name)) case ('flux') ! Prohibit user from tallying flux for an individual nuclide diff --git a/src/output.F90 b/src/output.F90 index 939451c6fb..9225e7a586 100644 --- a/src/output.F90 +++ b/src/output.F90 @@ -853,6 +853,17 @@ contains write(unit_,*) ' Outgoing Energy Bins:' // trim(string) end if + ! Write any change-in-angle bins if present + j = t % find_filter(FILTER_MU) + if (j > 0) then + string = "" + do i = 1, t % filters(j) % n_bins + 1 + string = trim(string) // ' ' // trim(to_str(& + t % filters(j) % real_bins(i))) + end do + write(unit_,*) ' Change-in-Angle Bins:' // trim(string) + end if + ! Write nuclides bins write(unit_,fmt='(1X,A)',advance='no') ' Nuclide Bins:' do i = 1, t % n_nuclide_bins @@ -1734,6 +1745,7 @@ contains filter_name(FILTER_MESH) = "Mesh" filter_name(FILTER_ENERGYIN) = "Incoming Energy" filter_name(FILTER_ENERGYOUT) = "Outgoing Energy" + filter_name(FILTER_MU) = "Change-in-Angle" ! Initialize names for scores score_names(abs(SCORE_FLUX)) = "Flux" @@ -2171,7 +2183,7 @@ contains label = "Index (" // trim(to_str(ijk(1))) // ", " // & trim(to_str(ijk(2))) // ", " // trim(to_str(ijk(3))) // ")" end if - case (FILTER_ENERGYIN, FILTER_ENERGYOUT) + case (FILTER_ENERGYIN, FILTER_ENERGYOUT, FILTER_MU) 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)) // ")" diff --git a/src/state_point.F90 b/src/state_point.F90 index 6b983231f6..23c46117dd 100644 --- a/src/state_point.F90 +++ b/src/state_point.F90 @@ -256,7 +256,8 @@ contains group="tallies/tally " // trim(to_str(tally % id)) // & "/filter " // to_str(j)) if (tally % filters(j) % type == FILTER_ENERGYIN .or. & - tally % filters(j) % type == FILTER_ENERGYOUT) then + tally % filters(j) % type == FILTER_ENERGYOUT .or. & + tally % filters(j) % type == FILTER_MU) then call sp % write_data(tally % filters(j) % real_bins, "bins", & group="tallies/tally " // trim(to_str(tally % id)) // & "/filter " // to_str(j), & @@ -843,7 +844,8 @@ contains group="tallies/tally " // trim(to_str(curr_key)) // & "/filter " // to_str(j)) if (tally % filters(j) % type == FILTER_ENERGYIN .or. & - tally % filters(j) % type == FILTER_ENERGYOUT) then + tally % filters(j) % type == FILTER_ENERGYOUT .or. & + tally % filters(j) % type == FILTER_MU) then call sp % read_data(tally % filters(j) % real_bins, "bins", & group="tallies/tally " // trim(to_str(curr_key)) // & "/filter " // to_str(j), & diff --git a/src/tally.F90 b/src/tally.F90 index e78e58f284..d3fce52905 100644 --- a/src/tally.F90 +++ b/src/tally.F90 @@ -990,7 +990,6 @@ contains matching_bins(i) = binary_search(t % filters(i) % real_bins, & k + 1, p % E) end if - end select ! Check if no matching bin was found @@ -1245,6 +1244,13 @@ contains n + 1, p % E) end if + case (FILTER_MU) + ! determine mu bin + n = t % filters(i) % n_bins + ! search to find incoming energy bin + matching_bins(i) = binary_search(t % filters(i) % real_bins, & + n + 1, p % mu) + end select ! If the current filter didn't match, exit this subroutine